 |
BELFEM 0.9.3
Berkeley Lab Finite Element Framework
|
|
BELFEM 0.9.3
Berkeley Lab Finite Element Framework
|
Module: src/physics/materials Purpose: How to create materials, query their properties, extend the roster, and understand the physical models behind every curve Date: 2026-09-15 Revision: 2.0
| Date | Version | Changes |
|---|---|---|
| 2026-01-16 | 1.0 | Initial documentation (reorganized from README.md) |
| 2026-07-14 | 1.1 | Assembly contract: full-signature policy for jc/n dependency routing; undeformed-mesh density rule |
| 2026-08-25 | 2.0 | Rewritten against the current source: real roster and API, RRR-at-construction and the cached lookup tables, cryogenic α, Wachtman moduli, per-metal Kohler curves, formula alloys; the invented Material_Metal sketch, populate_rho_lambda_databases, debye_cp and the textbook-only model formulas removed |
| 2026-09-15 | 2.1 | Elastic moduli: quasi-harmonic K and G on the thermal strain replace the Wachtman E and the constant-Grüneisen ν (§9.7); served moduli are isothermal; Nickel's ΔE dip retired; polycrystalline basis stated |
Read this section first.
set_jc_function(), set_n_function() and load_bh_curve() / set_bh_curve() all transfer ownership. The base Material destructor deletes the J_c and n functions; Metal's destructor deletes the B-H curve. The factory keeps no ownership of anything it returns: the caller deletes the material.
Because every argument is a real, swapped arguments still compile and run. The full overload is rho_powerlaw( normJ, T, normB, angleNxB ); see §7.3 for the assembly contract.
§7.4 has the reasoning.
set_RRR() fixes the residual resistivity ρ₀. Every material constant except mu (μ₀), q (1) and density_correction (1) starts as NaN, so a metal without an RRR has no ρ(T) and no λ(T). Metal::set_RRR also builds the field-dependent lookup table when tables are enabled (the factory's third argument, default true). A formula alloy takes its RRR through the same factory argument; it does not override set_RRR, so the inherited Material::set_RRR aborts with BELFEM_ERROR.
Temperature comes first. For a Metal, these evaluate Kohler's rule directly when no table exists; for an Alloy, they read the table and assert its presence.
MaterialDependency::angleBxJ is the field-current angle used by the metals' Kohler channel. MaterialDependency::angleNxB is the tape-normal-field angle used by the HTS J_c channel, unfolded to [0, π] since 2026-08-16. The FEM calculator routes on depends(), so declare the one you actually consume.
Highest power first (fn_polyval.hpp); the overloads take Cell< real > or std::vector< real >.
A user-defined material has only the properties its init function set. Querying an unset property dispatches through a null routing pointer.
buffer names the thin-shell φ-formulation buffer layer; the factory maps it to Magnesia and rejects a user-defined (usermat) material registered under that name; a built-in section with that label resolves to Magnesia and is not checked.
Supporting classes: MaterialFactory; BhCurve (interface) and BhSplineCurve (HDF5-backed); JcFunction with JcFunctionModifiedKim, JcFunctionDatabase, JcFunctionUserDefined; Abundance (molar masses, isotope mass variance).
Every material property — E, nu, cp, lambda, rho, alpha, density, debye, mu, … — is served through a function pointer. During construction, the material points it at one of three sources: a constant (set_constant), a spline sampled on a uniform grid (create_spline / set_spline), or a custom analytic routine (set_custom, dispatching to the virtual *_custom methods). Derivatives (dcpdT, drhodT, dlambdadT, …) follow the same routing. have( property ) reports whether any source was set; depends( property, dependency ) reports what the source varies with.
UserDefined, Ferro, HTS, PureMetal, LookupAlloy, CompositeAlloy, NonMetal. In practice: the nine metals and the formula Alloy are PureMetal; HastelloyC276 is LookupAlloy; YBCO is HTS; Magnesia is NonMetal. Ferro is not used by any built-in class — Iron and Nickel are PureMetal through Ferromagnetic. The FEM side interrogates depends() rather than the type wherever it can.
with MatFunc1 = real ( const Material *, real ) and the two- and three-argument analogs.
Labels are case-insensitive. Anything that is not a built-in label but parses as <Element><integer>… becomes an Alloy (§4.6). File arguments go through material::data_file(), which searches the run directory first and $BELFEM_DATA/material second.
An interface — B( H ), H( B ), mu( H ), nu( B ) — implemented by BhSplineCurve, which the factory builds from an HDF5 B-H database.
eval( B, angle ) and eval( B, angle, T ), each function overriding the arity that matches the dependencies it declares through JcParameter; both base versions raise an error. Three implementations: JcFunctionModifiedKim (analytic, §9.6), JcFunctionDatabase (tensor-mesh table over (T, log₁₀B, angle) with log-scaled values; T and B above the table are clamped to the table range, the angle is wrapped with π periodicity, and below the table's lowest field, n and a jc table with a measured self-field curve are continued to B = 0 with a C1 join in B and no increase with |B|; see share/material/README.md), JcFunctionUserDefined (plugin).
Molar masses and the isotope/alloy mass-variance parameter Γ from a stoichiometry, used by the metals (M) and by YBCO's Callaway model.
Every metal constructor follows a fixed sequence of create_*() calls; Copper is the reference implementation. The order matters because later steps read earlier curves:
All nine metals provide c_p, α, E, ν, ρ(T), λ(T), θ_D(T), density, and a Kohler curve.
Add to Metal: a reduced-magnetization curve (Crangle & Goodman 1971), the magnon term of the resistivity below the Curie point (rho_mag), a Debye curve inverted from measured resistivity with that term removed, and a Bloch–Grüneisen exponent of 4.5 instead of 5. Nickel's elastic moduli are those of the magnetically saturated state; the ΔE dip of demagnetized nickel is not represented (the class header of cl_Material_Nickel.hpp explains why). Neither class supplies a B-H curve of its own; assign one with load_bh_curve() if the magnetic solver needs it.
A Metal of type LookupAlloy with its own fitted c_p, α, E, ν, ρ(T), λ(T) and magnetic susceptibility χ(T). No field dependence: it declares neither Kohler dependencies nor a table, and only the ( T ) accessors apply.
Normal-state c_p, α, E, ν, ρ(T) and a λ(T) from the Callaway phonon model plus a Wiedemann–Franz electronic term (see callaway_thermal_conductivity.md — and its §8 before trusting the phonon parameters: the electronic term dominates the fit), with T_crit = 92.5 K. The superconducting resistivity is the E–J power law with J_c and n supplied as constants or JcFunction objects:
c_p, α, E, ν and λ(T); no electrical resistivity. This is the material behind the reserved label buffer. Its α, like Hastelloy's and YBCO's, follows c_p below room temperature through the anchored Grüneisen branch (thermal_expansion_from_heat_capacity.md, §7b).
Sn60Pb40, Fe71Cr19Ni10, … in integer mass percent, homogenized from the metal roster with a caller-supplied RRR. What is mixed, how, and — above all — for which alloys the result is trustworthy is the subject of alloy_homogenization.md. Read its §3 before using one for anything cryogenic.
See §6.
material copper -r 100 -b 5 -a 90 -t 4 300 4 prints ρ, λ, c_p, E, ν, α over the range as far as each property exists. Without -c the executable constructs the material with tables disabled, so a Metal still evaluates Kohler's rule directly and an Alloy cannot be asked for field dependence.
A user material is a shared library that exports one C-linkage function named <Label>_init. That function receives the Material and sets its properties. example_user_material.cpp and UserMaterialTemplate.cmake are the templates.
Build with the CMake template (C++17, the project standard), then
The label must match the function prefix; the factory resolves <Label>_init with dlsym and fails loudly otherwise. Two- and three-argument functions declare their dependencies in order (normB, angleNxB, T for a J_c, for instance) and receive the arguments in that order. Superconducting user materials supply J_c and n through JcFunction objects, not through jc_custom (§7.3).
Rules: use extern "C"; make the first parameter of every function const Material *; set only the properties you have — everything else stays have() == false; declare dependencies truthfully because the FEM dispatch trusts them.
depends() is declarative: it tells the caller which overload carries information, but it does not validate the call. cl_FEM_Calculator uses this flag to decide whether a block's resistivity takes the ( T, B, beta ) path.
T, normB (in-plane field magnitude), angleBxJ (field–current angle, Kohler), angleNxB (normal–field angle, HTS), normH (field strength, μ of ferromagnets), normJ, Jc, rho.
angleNxB is unfolded to [0, π] since 2026-08-16: β = acos(n·b/|b|), where +n is the master-side facet normal (the layer-stack direction). β < π/2 means the field has a component along +n. Measured jc(β) tables consume the full range; Kim-type analytic laws are even in β and unaffected. Despite the historical name, this is NOT the angle between current and field — that is the metals' angleBxJ Kohler channel.
The FEM assembly path always calls the full T-bearing HTS overloads —
— regardless of what a given material actually depends on. The material, not the caller, decides which arguments matter: internally, Material::jc_eval() and Material::n_eval() route on the declared dependencies of the attached JcFunction / n-function (depends_on( JcParameter::T ) selects the eval( B, angle, T ) or eval( B, angle ) override) and fall back to the plain constants (MaterialProperty::jc / MaterialProperty::n) when no function is attached.
Consequences for material implementers:
Real density changes with temperature — but all BELFEM computations run on the undeformed mesh, and the transport properties are already corrected for thermal expansion. The thermal mass matrix must therefore use the density at which the mesh is not deformed, which is by default room temperature:
calculator::MaxwellData caches this value once at construction (density() accessor); do not re-derive it per integration point. density( T ) itself is ρ_ref / l(T)³ with l the relative length from the α spline's integral.
A conductor with magnetoresistance
HTS with an analytic J_c
HTS with tabulated J_c and n
A ferromagnet with a B-H curve
A solder from its composition
A plugin material
Each model is described as the code implements it; the theory documents contain the derivations.
ρ(T) = ρ₀ + ρ_i(T), ρ_i(T) = A · J_n(θ(T)/T) / (θ(T)/T)ⁿ
with J_n the Debye integral (tabulated once per exponent by debye.f90), n = 5 for the simple metals and 4.5 for Iron and Nickel (n_bloch_gruen), and θ(T) the effective Debye curve of the class — derived from c_p, from measured ρ(T), or fitted — so that the one-parameter formula reproduces the measured phonon resistivity. set_rho_i_ref( T_ref, ρ_i,ref [, θ] ) fixes A at a literature point. ρ₀ follows from the RRR by the regula falsi in Metal::set_RRR:
RRR = ( ρ_i(T_ref) + ρ₀ ) / ρ₀ at T_ref_rho_i (273.15 K for Copper).
Ferromagnets add a magnon term rho_mag(T) below the Curie point.
λ(T) = 1 / ( w₀ + w_i + w_i0 ), w₀ = ρ₀ / (L₀ T), w_i = hust(P, T), w_i0 = C·w_i·w₀ / (w_i + w₀)
with the eight Hust coefficients per metal (set_lambda_coefficients, NBS Special Publication 260-90). In field, λ(T, B, β) = λ(T) · ρ(T) / ρ(T, B, β) — a Wiedemann–Franz closure with L = L₀.
Sommerfeld–Debye cubic γT + βT³ below 0.02·θ_D, a fifth-order beam polynomial to the first control point, two or three cubic Béziers in ln c_p over ln T, and a linear tail in T. The representation and its constraints are in thermal_expansion_from_heat_capacity.md, §3.
Above a split temperature T* = min(0.618·θ_D, 273.15 K) the tangent coefficient α = (1/L) dL/dT comes from a Bézier fitted to ΔL/L referred to 293.15 K. Below it, α = C(T)·c_p(T) with ln C a cubic matched at T* — the Grüneisen equation of state, which gives the correct T³ approach to zero that a ΔL/L fit cannot. thermal_expansion_from_heat_capacity.md is the full account, including the two guards and the Grüneisen diagnostic that reads O(2) for every metal.
Δρ/ρ_0T = kohler( B, S, β ) with the similarity parameter S = ρ_ref / ρ_0T, where ρ_0T is the zero-field resistivity at the current temperature and ρ_ref the one the curve was fitted at. Each metal carries its own longitudinal and transverse curves in ln(B·S) — Bézier or polynomial segments with a low-field Hermite cubic, a mid-field fit and a saturated or linear tail — and combines them by Pippard's angular interpolation A_∥ cos²β + A_⊥ sin²β. Data sources are cited per class (Lüthi 1960, Fickett 1972, Klaffky & Coleman 1974, Kozlova & Kondorskii 1963, Arajs & Dunmyre 1965, …). The table <label>_RRR<n>.hdf5 stores ln ρ over (T, log₁₀B, β) on a tensor mesh; Metal evaluates the curves directly when no table exists, Alloy requires the table.
E = E_c (J/J_c)ⁿ → ρ_powerlaw = (E_c / J_c) · (J/J_c)^(n−1), with E_c = 10⁻⁴ V/m by default.
JcFunctionModifiedKim( Jc0, B0, k², α ):
J_c(B, θ) = Jc0 / [ 1 + √(k² sin²θ + cos²θ) · B/B0 ]^α
— the anisotropic Kim–Anderson form; the factory's four arguments map to (Jc0, B0, k², α) in that order. Tabulated J_c and n come from JcFunctionDatabase.
The nine pure metals share one closure, Metal::create_mech( E₂, ν₂, T₂, δ_K, δ_G ). It uses a quasi-harmonic picture: an elastic modulus depends on the volume of the lattice, not on the temperature directly, and its logarithmic volume derivative is a material constant. Garai & Laugier 2007 (10.1063/1.2424535, Eq. 8–10) derive this form for the isothermal bulk modulus with a constant Anderson–Grüneisen parameter; BELFEM applies the same form to the shear modulus, which is a modelling choice, not a result of that paper. Using the logarithmic volumetric strain that the expansion curve already integrates for the density,
ε(T) = 3 ∫₀ᵀ α dT = 3 ln( l(T)/l(0) ), K(T) = K₀ exp(−δ_K ε), G(T) = G₀ exp(−δ_G ε), E = 9KG/(3K + G), ν = (3K − 2G)/(6K + 2G).
The two softening constants and the anchor (T₂, E₂, ν₂) are fitted offline per metal to measured elastic constants over 0–300 K (the data span 4 K to room temperature; single-crystal sets are Hill-averaged to the isotropic polycrystal; the sources with DOIs are listed in material_property_sources.md). The form determines three properties independently of the fit: E and ν have zero slope at 0 K because α does; ν rises with temperature exactly when the shear modulus softens faster than the bulk modulus (δ_G > δ_K), which measurement shows for every metal in the roster except iron, whose isothermal ν is flat; and K and G stay positive, so ν stays inside (−1, ½). The constructor samples the served curves over the whole temperature range and rejects constants that make ν fall by more than 1e-3.
Isothermal. The served moduli are the isothermal ones, which a quasi-static stress analysis needs. Ultrasonic data are adiabatic; the shear modulus is unchanged between the two descriptions, the bulk modulus is converted with the exact identity K_T = K_S/(1 + α_V² T K_S/(ρ c_p)) before the fit, using the material's own α, ρ and c_p. For copper the difference is 3 % on K and 0.004 on ν at room temperature and vanishes at low temperature.
Polycrystalline basis. Every metal in the roster represents an isotropic polycrystal. Where the source is a single-crystal measurement, the constants are Hill-averaged (Voigt–Reuss mean of the shear modulus; the bulk modulus of a cubic crystal is unique). For soft, strongly anisotropic metals the Hill average is far above the static polycrystalline value, because the shear modulus relaxes at low frequency; lead is therefore served with static-level moduli (crystal bulk modulus, Blanke's E), see its class header.
Two metals require special treatment. Nickel is fitted to constants measured in a saturating field, because the deep minimum of demagnetized nickel's Young's modulus below the Curie point is a domain-wall effect that vanishes at saturation, and BELFEM's nickel sits in tesla-level fields; the class header of cl_Material_Nickel.hpp documents this decision. Chromium's spin-density-wave anomalies collapse its bulk modulus by 20 % toward the Néel point; no smooth monotone closure follows that, so chromium is served with δ_G = δ_K, a constant ν = 0.2371 taken from its cryogenic plateau, and the anomaly is documented as unrepresented in cl_Material_Chromium.hpp.
Sources and validity ceilings. Measured sources with data points (Ledbetter, Neighbours & Alers, Kamm & Alers, Rayne & Chandrasekhar, Alers et al., Palmer & Lee, Kim & Ledbetter) set the levels and the softening constants. Blanke's compilation cites no sources and shows no data points; it is used only to bound the range: where the quasi-harmonic E(T) departs from Blanke's curve by more than 5 % in shape (both normalized at room temperature), the material's T_max is lowered to that temperature — copper 1000 K, silver 900 K, chromium 570 K, white tin 400 K. Iron (855 K) and aluminum (never beyond 10 %) keep their ceilings. All served moduli are dynamic moduli made isothermal; static moduli, which an elastostatic calculation would prefer, are lower for metals with anelastic relaxation at low frequency; lead, where the difference is 50 %, is served at the static level; for the others it is a few percent.
The Grüneisen parameter is no longer part of the closure. It is printed by the material tool as a diagnostic, γ = α_V K_S/(ρ c_p) at 298.15 K, the temperature of the printed density, and retains its role for the expansion curve: an independent cross-check that reads about 2 for a metal whose α, c_p, ρ and K are mutually consistent. Evaluated along the temperature axis from the served moduli it rises gently from its low-temperature value to the room-temperature one (copper 1.88 → 2.03, silver 2.30 → 2.43, lead 2.42 → 2.74, iron 1.46 → 1.74 between 20 and 293 K) and ends inside the literature range for every metal in the roster.
Used by YBCO only; callaway_thermal_conductivity.md.
This table is generated from cl_Material.hpp. The index is the storage slot, and UNDEFINED sizes the arrays. Constants and non-constant properties share the enumeration; everything from T_crit upward is a constant.
| Property | Index | Unit | Meaning |
|---|---|---|---|
| density | 0 | kg/m³ | density |
| E | 1 | Pa | Young's modulus |
| nu | 2 | - | Poisson's ratio |
| cp | 3 | J/(kg·K) | specific heat capacity |
| lambda | 4 | W/(m·K) | thermal conductivity |
| mu | 5 | H/m | magnetic permeability, defined as ∂B/∂H |
| rho | 6 | Ω·m | electric resistivity |
| alpha | 7 | 1/K | thermal expansion coefficient, defined as (1/l)·∂l/∂T, not (1/l)·Δl/ΔT! |
| Rp02 | 8 | Pa | yield stress |
| debye | 9 | K | Debye temperature |
| rho_i | 10 | Ω·m | inner resistivity of noble metal |
| kohler_trans | 11 | - | Kohler parameter for transverse magnetoresistance |
| kohler_long | 12 | - | Kohler parameter for longitudinal magnetoresistance |
| T_crit | 13 | K | critical temperature |
| M | 14 | kg/mol | molar mass |
| Gamma | 15 | - | impurity parameter |
| R | 16 | J/(kg·K) | specific gas constant |
| T_max | 17 | K | maximum temperature |
| ref_density | 18 | kg/m³ | reference density |
| T_ref_density | 19 | K | temperature at reference density |
| gamma | 20 | J/(kg·K²) | linear Debye parameter, cv = γ·T + β·T³ = ∂cp/∂T at T=0 |
| beta | 21 | J/(kg·K⁴) | cubic Debye parameter, cv = γ·T + β·T³ |
| q | 22 | - | number of atoms per molecule (default: 1) |
| debye0K | 23 | K | Debye temperature at 0 K |
| rho_i_ref | 24 | Ω·m | reference inner resistivity |
| T_ref_rho_i | 25 | K | temperature at reference inner resistivity |
| A_bloch_gruen | 26 | - | A-parameter for Bloch-Grüneisen law |
| n_bloch_gruen | 27 | - | n-parameter for Bloch-Grüneisen law |
| RRR | 28 | - | residual resistivity ratio |
| rho_0 | 29 | Ω·m | residual resistivity |
| layer_thickness | 30 | m | characteristic length, e.g., layer thickness |
| ec | 31 | V/m | critical electric field |
| jc | 32 | A/m² | critical current density |
| n | 33 | - | exponent for power law |
| Tcurie | 34 | K | Curie temperature |
| A_electron_magnon | 35 | Ω · m / K² | electron-magnon scattering parameter |
| A_spin_disorder | 36 | Ω · m | spin disorder amplitude |
| grueneisen | 37 | — | Grüneisen parameter |
| density_correction | 38 | - | density scaling factor, to correct solder thickness |
| rho0_pure | 39 | — | resistivity from TPCR for pure metals |
Units and meanings are the header's own comments. Two of them are imprecise as of 2026-08-25: A_bloch_gruen is marked dimensionless but carries Ω·m, and rho0_pure has no unit and is a resistivity [Ω·m].
The mu row's "defined as ∂B/∂H" is imprecise, and worth unpicking because the Newton path depends on it. mu is index 5, i.e. inside the non-constant range — a material may hold a constant permeability, or a B-H curve, or a custom mu( H ), and loading a curve repoints mFunctionMu at mu_bhcurve (cl_Material.cpp:958-960).
If you need ∂B/∂H — the tangent a Newton iteration wants — build it from both. For B = μ(H)·H,
BELFEM parallelizes with MPI and is deliberately not thread-safe internally (doc/coding_philosophy.md). Within that policy the materials module behaves as follows:
No wall-clock figures are recorded in the tree; measure before relying on any.
1.18.0