 |
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
|
BELFEM 0.9.0
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-08-25 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 |
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), log-scaled values; T and B are clamped to the table range, the angle is wrapped with π periodicity), 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 E(T) carries the ΔE dip across the Curie point as two Wachtman branches joined by a Bézier bridge. 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.
E(T) = E₀ − b·T·exp(−T₀/T), fitted per metal (mostly against Blanke, Thermophysikalische Stoffgrößen, 1989; Chromium against Armstrong & Brown 1964; Indium against Kim & Ledbetter 1998). ν(T) is not tabulated: Metal::create_mech( E₀, b, T₀, T₂, ν₂ ) takes one anchor (T₂, ν₂), forms K_T = E/(3(1 − 2ν)), the adiabatic K_S = K_T / (1 − T α_V² K_T /(ρ c_p)) and the Grüneisen parameter γ = α_V K_S / (ρ c_p), holds γ constant, and recovers ν(T) = ½ − E/(6K_T(T)) along the spline grid with K_S(T) = γ ρ c_p / α_V. The 0 K value is extrapolated with zero slope; the result is checked to stay inside (−1, ½). Nickel supplies its own E(T) (two Wachtman branches and a Bézier bridge across the Curie point) and only the ν construction.
Isothermal versus adiabatic: the tabulated moduli are dynamic (adiabatic); the difference on E is 0.3–0.5 % at room temperature and vanishes at cryogenic temperatures, and is not corrected.
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