4. Materials, Components, and Property Functions

This chapter specifies the constitutive layer of Pyrolysis.jl: how a pyrolyzing substance is decomposed into discrete phase-typed components, how each component carries its thermophysical property functions, how those functions are evaluated and differentiated, and how components and reactions are assembled into the static Material object that every physics kernel reads. The governing PDE system (Chapter 3, Governing Equations) closes only once these laws are supplied; the homogenization of per-component properties into mixture (effective) properties is the subject of Chapter 5, Effective Properties and Mixing Rules, to which this chapter forwards repeatedly. Throughout, the spatial coordinate is $z$ (axial, through-thickness, $z=0$ substrate / bottom, $z=L$ exposed surface), the primary species state variable is the mass concentration $\xi_j$ [kg·m⁻³], and the heat-of-reaction sign convention is $h>0$ endothermic with the volumetric source $Q_\mathrm{rxn} = -\sum_r h_r r_r$ (overload note H1 of Chapter 2). In the kinetics passages of this chapter the index $i$ denotes a reaction and $j$ a component (overload note I1); elsewhere $i$ is a cell index.

4.1 The phase model: solids, liquids, and gases

4.1.1 Why phase is a type, not a field

A pyrolyzing material is a heterogeneous mixture of condensed-phase species (virgin polymer, char, moisture) and gaseous species (pyrolyzate, water vapor, oxidizer). The three phases obey structurally different conservation laws (Chapter 3): condensed species are Lagrangian — they do not transport relative to the matrix, so their concentration changes only through reaction — whereas gaseous species are Eulerian — they flow through the porous matrix and therefore carry a flux divergence in addition to a reaction source. The gas-vs-condensed dispatch appears in the species residual:

\[\frac{\partial \xi_j}{\partial t} = \begin{cases} -\,\nabla\!\cdot\mathbf{J}_j \;+\; \left.\dfrac{\partial \xi_j}{\partial t}\right|_\mathrm{rxn} & j \text{ is gas},\\[2.2ex] \left.\dfrac{\partial \xi_j}{\partial t}\right|_\mathrm{rxn} & j \text{ solid or liquid.} \end{cases} \tag{4.1}\]

Because this branch is taken inside the innermost residual loop, over every cell and every component, it must not incur a runtime type check. Pyrolysis.jl encodes the phase at the type level: there are three concrete subtypes of AbstractComponent, and the phase predicate is a method that dispatches on the singleton type, so the compiler eliminates the branch entirely at code generation. The runtime Phase enum (SOLID=1, LIQUID=2, GAS=3) exists only for printing, logging, and comparison; it is not consulted in hot paths.

@enum Phase begin SOLID = 1; LIQUID = 2; GAS = 3 end

Implementation note. Phase types and predicates are defined in SolidComponent, LiquidComponent, GaseousComponent, and is_solid / is_liquid / is_gas / phase in src/Materials/components.jl. The compile-time-unrolled application of Eq. (4.1) is apply_species_residual_unrolled!, src/Materials/materials.jl, which emits one statically dispatched is_gas(material.components[j]) branch per component — folded to a constant at codegen.

4.1.2 The three component types

SolidComponent is the matrix-forming condensed phase. It carries seven property-function fields, an optional constant surface absorptivity, and two scalar (Float64) fields:

struct SolidComponent{Pρ,Pc,Pk,Pλ,Pκ,Pε,Pa,Pρi} <: AbstractComponent
    name::Symbol
    density::Pρ                # bulk density ρ_j           [kg/m³]
    heat_capacity::Pc          # c_{p,j}(T)                 [J/(kg·K)]
    thermal_conductivity::Pk   # k_j(T)                     [W/(m·K)]
    gas_transfer_coeff::Pλ     # λ_j(T)                     [m²/s]
    permeability::Pκ           # κ_j(T)                     [m²]
    emissivity::Pε             # ε_j(T)                     [-]
    surface_absorptivity::Pa   # a_j, Nothing = Kirchhoff   [-]
    absorption_coeff::Float64  # α_j (mass basis)           [m²/kg]
    swelling_factor::Float64   # γ_j                        [-]
    intrinsic_density::Pρi     # skeletal density ρ_{i,j}   [kg/m³]
end

The seven property fields are each parameterized by their own type (Pρ, Pc, …), so a component built with a constant conductivity and a polynomial heat capacity is a distinct concrete type from one built with two constants. This preserves type stability and allocation-free evaluation through the heterogeneous component tuple. The emissivity ε_j(T) is a property function like k_j(T) (the ε keyword accepts the same forms) and is evaluated at the cell temperature by the mixing rule (§5.8). surface_absorptivity is the optional constant absorptivity a_j for incident (heater-spectrum) radiation, Pa <: Union{Nothing,Float64}: Nothing selects Kirchhoff's law a_j = ε_j(T) at the type level, so the absorptivity mixture (§5.8.2) is a static sum. The two scalar fields (α_j, γ_j) are stored as plain Float64 because they are not given temperature dependence in the current model.

LiquidComponent has the identical field layout to SolidComponent (bulk and intrinsic density, the five property functions, the optional surface absorptivity, the two scalars). It is the matrix-forming condensed phase used for moisture / liquid water in hygroscopic materials. Its convenience constructor differs only in defaults (see §4.1.5): $\lambda$ defaults to $2\times10^{-5}$ m²·s⁻¹, $\varepsilon$ to 0.95 (condensed liquids are strong emitters — water ≈ 0.96), $\gamma$ to 0 (volumeless bound moisture, ThermaKin L SWELLING: 0), and $\kappa$ to 0 (a liquid normally fills pore space and blocks gas flow).

GaseousComponent represents species that flow through the matrix:

struct GaseousComponent{Pρ,Pc,Pk,Pλ,Pε} <: AbstractComponent
    name::Symbol
    molar_mass::Float64        # M_j                        [kg/mol]
    density::Pρ                # ρ_g(T), default ideal gas  [kg/m³]
    heat_capacity::Pc
    thermal_conductivity::Pk
    gas_transfer_coeff::Pλ
    emissivity::Pε             # ε_j(T)
    absorption_coeff::Float64
    swelling_factor::Float64
end

A gas has no permeability field, no intrinsic density, and no surface absorptivity (get_absorptivity always returns its ε_j(T)) — it is not part of the matrix skeleton. It carries instead a molar_mass $M_j$, used by the ideal- gas density law (§4.3), the ideal-gas pressure closure (Chapter 9), and the Chapman–Enskog diffusivity estimate (§5.10).

4.1.3 Compile-time phase accessors

All cross-phase property access goes through dispatched accessors that return phase-appropriate defaults, so a kernel may ask any component for any quantity without a runtime field-existence test:

Each is @inline and resolves at compile time on the singleton type; none branch on a stored field. The gas's get_permeability and get_intrinsic_density return exactly 0.0, which is the algebraically correct value for the matrix mixing rules (a gas contributes nothing to permeability or to the solid skeleton).

Implementation note. Accessors get_intrinsic_density, get_permeability, has_permeability, get_emissivity, get_absorptivity, is_solid/is_liquid/is_gas, phase, molar_mass, all in src/Materials/components.jl. Kernels read the remaining property fields (density, heat_capacity, …) directly.

4.1.4 Bulk versus intrinsic (skeletal) density

The single most important constitutive distinction at the component level is between two densities (overload note A7 of Chapter 2):

• Bulk (pure-phase) density $\rho_j$ [kg·m⁻³] — the mass of pure phase $j$ per unit of the bulk volume it would occupy in the mixture. This is the density used by every mixing rule, by the mixture-density formula, and by the volume- change kernels (Chapter 10). Field density.

• Intrinsic / skeletal density $\rho_{i,j}$ [kg·m⁻³] — the density of the cell-wall material itself, ignoring any pore volume internal to the component. This is used by one formula only, the porosity (Chapter 5, §5.3):

  \[\phi = \mathrm{clamp}\!\left(1 - \sum_{j\in\{S,L\}} \frac{\xi_j}{\rho_{i,j}},\; 0,\; 1\right). \tag{4.2}\]

  Field intrinsic_density.

When a component is constructed without an explicit intrinsic density, the field defaults to the bulk density, $\rho_{i,j}=\rho_j$: the pure phase is assumed to have no internal pores. The resolution logic enforces that a caller pass either ρ_intrinsic or a pure-phase porosity $\phi$ (from which $\rho_{i,j}=\rho_j/(1-\phi)$ is computed, requiring a constant bulk $\rho_j$), never both:

1. both ρ_intrinsic and φ given        → ArgumentError
2. ρ_intrinsic given                    → to_property(ρ_intrinsic)
3. φ given (requires ρ::Real, 0≤φ<1)    → ConstantProperty(ρ / (1 − φ))
4. neither                              → fall back to bulk ρ

Implementation note. _resolve_intrinsic_density, src/Materials/components.jl. The φ path requires raw_ρ isa Real and 0 ≤ φ < 1, throwing an ArgumentError otherwise. Gas components do not store an intrinsic density; get_intrinsic_density(::GaseousComponent, ::Real) = 0.0, which is consistent with Eq. (4.2) summing only over solid/liquid.

4.1.5 The swelling factor

The swelling factor $\gamma_j\in[0,1]$ controls whether changes in component $j$'s concentration change the bulk volume of the cell. Physically it answers: when this component is created or destroyed, does the material expand/shrink, or does the mass simply enter/leave the existing pore space?

• $\gamma_j = 1$ (default for solids): full contribution. Depletion of the solid shrinks the bulk volume; this is the charring/recession behavior.
• $\gamma_j = 0$ (default for gases): no contribution. Gas fills existing pores and escapes without expanding the bulk.
• $\gamma_j \in (0,1)$: partial contribution, the intumescent case, where evolved gas creates new bulk volume.

$\gamma_j$ enters two places: the mixture-density formula (§4.2.2, where a gas with $\gamma_j=0$ contributes zero specific volume) and the volumetric strain rate $\theta = \sum_j \gamma_j \rho_j^{-1}\,\mathrm{d}\xi_j/\mathrm{d}t$ that drives the ALE mesh motion (Chapter 10). It also weights the matrix volume fractions used for conductivity/emissivity (Chapter 5, §5.1): every component carries weight $\gamma_j\,\xi_j/\rho_j$ there, so a $\gamma_j = 0$ component (the canonical gas) carries no weight in the matrix radiative surface.

Why the cap $\gamma_j \le 1$ is a model decision, not a missing feature. Construction rejects $\gamma_j > 1$ (_validate_swelling_factor) because the same $\gamma$ weight feeds three mixing rules with volume-fraction semantics — the strain kernel, the mixture density, and the $k/\varepsilon$ matrix volume fractions — and a weight above 1 corrupts all three (volume fractions would not sum consistently). Expansion beyond the reactant's volume is already expressible the physical way: give the product a low bulk density ($V_f/V_0 = \sum_j y_j\gamma_j\rho_A/\rho_j$ — the V-CPL-03 P1 verification permutation demonstrates ×2 growth), with $\gamma \in (0,1]$ tuning how much of that specific volume the matrix keeps. A decoupled intumescent expansion factor (strain contribution independent of the mixing weight) would be a separate modeling design; no current use case demands it.

4.2 Component property fields and their roles

4.2.1 Inventory

Each property field of a component, its symbol, its SI units, and where it is consumed downstream:

Note that $\alpha_j$ is a mass-basis absorption coefficient [m²·kg⁻¹]; it is multiplied by concentration to give the volumetric (extinction) absorption coefficient $\alpha_\mathrm{eff}=\sum_j\xi_j\alpha_j$ [m⁻¹] used in Beer–Lambert and P1 radiation (Chapter 8). A solid usually carries a nonzero $\lambda_j$ in practice (e.g. $2\times10^{-5}$ m²·s⁻¹; the SolidComponent default is 0, and validate_problem warns when a gas-producing Fickian material leaves it there, §5.6): although solids do not themselves diffuse, their $\lambda_j$ contributes to the mixture gas-transfer coefficient through the volume-weighted mixing rule, matching the ThermaKin convention.

4.2.2 Mixture density with per-component swelling

The mixture density combines per-component bulk densities through specific volumes (volume per unit mass), every component weighted by its own swelling factor (the universal γ semantics shared with θ and the mixing-rule volume fractions, §5.1.3):

\[\frac{1}{\rho} = \sum_{j} \gamma_j\,\frac{Y_j}{\rho_j}, \qquad Y_j = \frac{\xi_j}{\rho_\mathrm{tot}},\quad \rho_\mathrm{tot}=\sum_k \xi_k. \tag{4.3}\]

Here $Y_j$ is the mass fraction. A non-swelling gas ($\gamma_j=0$) contributes zero specific volume — its mass is "free" in the pore space and does not change the bulk volume — whereas a swelling gas ($\gamma_j=1$) contributes its full specific volume like a solid. This is the discrete realization of the swelling semantics of §4.1.5.

Important — correction to a common gloss. The swelling weight in Eq. (4.3) is applied per component, with each component's own swelling_factor — not an averaged $\gamma_\mathrm{avg}$ multiplying the gas sum, and not gases only. The source evaluates γ_j · Y_j/ρ_j term by term for every component (a γ = 0 term contributes nothing).

Implementation note. mixture_density and the inline, allocation-free kernel _specific_vol_with_swelling_direct, src/Properties/mixture_rules.jl. The inline form computes $Y_j=\xi_j\,/\rho_\mathrm{tot}$ on the fly to avoid allocating an intermediate mass-fraction tuple. Degeneracy guards: if $\rho_\mathrm{tot}<10^{-20}$ or the accumulated specific volume $<10^{-20}$, the function returns 0.0.

4.2.3 Porosity and its single use

Porosity $\phi$ (Eq. 4.2) is the void fraction of the bulk not occupied by the solid/liquid skeletal frame. It is computed from intrinsic densities and is the canonical definition $\phi = \mathrm{clamp}(1-\sum_{j\in\{S,L\}}\xi_j/\rho_{i,j},0,1)$. Its derivatives are exactly (in the unclamped interior; at an active clamp the derivative is 0, and the AD path used in production returns exactly that)

\[\frac{\partial \phi}{\partial \xi_j} = \begin{cases} -\,1/\rho_{i,j} & j\in\{S,L\},\\ 0 & j\in G. \end{cases} \tag{4.4}\]

A gas does not enter $\phi$ because it occupies the void that the solid frame leaves; adding gas mass does not displace skeletal material. Porosity feeds (i) the ideal-gas pressure closure $P=N(\xi,T)/\phi$ (Chapter 9) and (ii) the Carman–Kozeny permeability closure (Chapter 5, §5.7). The all-phase gas volume fraction $\sum_{j\in G} v_j$ is a different quantity; the intrinsic-density form Eq. (4.2) is canonical and is what porosity computes and update_darcy_properties! uses (overload note A8 of Chapter 2).

Implementation note. Canonical porosity: porosity → _porosity_void_frac, in src/Properties/mixture_rules.jl (derivatives are taken by forward-mode AD through porosity itself; §5.3.2). In _porosity_void_frac a component contributes $\xi_j/\rho_{i,j}$ only when is_solid(c) || is_liquid(c) and $\rho_{i,j}>0$. Full development is in Chapter 5.

4.3 Ideal-gas density and the gaseous-component default

By default a GaseousComponent constructed via the convenience constructor is given a density obeying the ideal-gas law at the reference pressure:

\[\rho_g(T) = \frac{P_\mathrm{ref}\,M_j}{R_g\,T}, \tag{4.5}\]

with $P_\mathrm{ref}=101325$ Pa, $R_g=8.31446261815324$ J·mol⁻¹·K⁻¹, and $M_j$ the gas molar mass. The constructor wraps this in a CallableProperty:

ρ_func = CallableProperty(T -> P_REF * M_float / (R_GAS * T))

Two modeling points deserve emphasis:

1. The density used for mixing/volume is at $P_\mathrm{ref}$, not the local pressure $P$. The component's density(T) is consulted for volume fractions, specific volumes, and the mixture density. The local pressure-dependent gas density (which matters for the pore-pressure closure $P=N/\phi$) is handled separately by the ideal-gas formula in the Darcy machinery (Chapter 9), not by the component's density field. The component-level density at $P_\mathrm{ref}$ is what enters the specific volume in mixing rules. In the swelling-weighted bases (mixture density, conductivity/emissivity weights) the gas's contribution is usually gated to zero by $\gamma_j=0$ (Eq. 4.3); in the all-phase basis used for λ mixing (§5.6) the gas volume is present regardless of γ and is over-weighted by ~$P/P_\mathrm{ref}$ in a pressurized cell — acceptable at the ≤ atm-scale overpressures of the validated cases. The pressure closure $P=N/\phi$ and the porosity are unaffected (they never use this density).

2. The user may override the gas density by passing any property function (constant, polynomial, table, callable) — Eq. (4.5) is merely the default produced when one constructs a gas with GaseousComponent(:name; M=…, c=…, k=…, λ=…).

Implementation note. GaseousComponent(name; M, c, k, λ, ε=0.0, α=0.0, γ=0.0) constructor, src/Materials/components.jl. Defaults for a gas are $\varepsilon=0$, $\alpha=0$, $\gamma=0$ (a gas escapes without swelling and is not part of the radiating surface).

4.4 Property-function types

A property function is a callable object $p$ such that $p(T)$ (or, for state-dependent properties, $p(T,\xi)$) returns the property value at the given state. All eight concrete types subtype AbstractPropertyFunction and implement both the (p)(T) and (p)(T, ξ) call signatures (the latter ignoring $\xi$ for all but StateDependentProperty). This uniformity lets a single component field hold any property law, and lets the kinetics evaluator pass $\xi$ to a heat-of- reaction property without knowing whether it is composition-dependent.

The property bracket. Every temperature-dependent law evaluates on the clamped argument $\tilde T = \operatorname{clamp}(T, T_{\min}, T_{\max})$ with $[T_{\min}, T_{\max}] = [200, 3000]$ K (T_PROPERTY_MIN, T_PROPERTY_MAX, src/Materials/constants.jl; the same bracket clamps the surface-temperature Newton iterate, §12). Inside the bracket the clamp is the identity with unit derivative, so nothing changes for any physical state; outside it the law sees the nearest bracket end. A Newton trial state with a nonphysical (even negative) temperature therefore produces a failed Newton step, not a DomainError out of a fractional power or a fitted polynomial's far extrapolation. No derivative is hand-coded: the Jacobian and sensitivity passes evaluate each law on ForwardDiff.Dual arguments (§4.5), so the derivatives quoted below are what automatic differentiation of the law returns — and, through the clamp, zero outside the bracket.

The component constructors sample every law on the bracket at 1 K steps and reject a component whose bulk or intrinsic density, heat capacity or thermal conductivity is not positive there (a gas may have $c_p = 0$), whose gas transfer coefficient or permeability is negative, or whose emissivity leaves $[0, 1]$. The error names the first offending temperature: a linear fit $k = 0.3 - 4.2\times10^{-4}\,T$ is rejected at 715 K because it turns a cell into a perfect insulator once it reaches that temperature. Hold a fit constant outside its measured range (T -> law(clamp(T, T_lo, T_hi))) or pass validate_properties = false for laws built programmatically from search parameters, where an unphysical value should reach the objective function instead of throwing.

4.4.1 ConstantProperty

\[p(T) = \texttt{value}, \qquad \frac{\mathrm{d}p}{\mathrm{d}T}=0.\]

The workhorse for properties given a single number. It has a fast path in the parameter-evaluation interface (§4.7) that bypasses the dynamic-dispatch hop.

4.4.2 LinearProperty

\[p(T) = a + bT, \qquad \frac{\mathrm{d}p}{\mathrm{d}T}=b.\]

Fields a (intercept at $T=0$) and b (slope). Constructed automatically from a two-tuple via to_property((a,b)).

4.4.3 PolynomialProperty

\[p(T) = \sum_{m=0}^{N-1} c_m\, T^{m} = c_0 + c_1 T + c_2 T^2 + \cdots,\]

evaluated with evalpoly (Horner's scheme) for numerical stability and to avoid recomputing powers. Forward-mode AD through Horner's scheme returns the exact derivative

\[\frac{\mathrm{d}p}{\mathrm{d}T} = \sum_{m=1}^{N-1} m\,c_m\,T^{m-1} = c_1 + 2c_2 T + 3c_3 T^2 + \cdots . \tag{4.6}\]

4.4.4 ArrheniusProperty

\[p(T) = A\,\exp\!\left(-\frac{E}{R_g\,T}\right), \qquad \frac{\mathrm{d}p}{\mathrm{d}T} = p(T)\,\frac{E}{R_g\,T^2}. \tag{4.7}\]

Fields A (pre-exponential), E (activation energy [J·mol⁻¹]), and T_ref (a reference temperature, defaulting to $T_\mathrm{ref}=298.15$ K, carried for normalization purposes but not used in the value itself). The law shares its functional form with the reaction-rate constant (§4.5.3, Chapter 6).

4.4.5 TableProperty

Tabulated $(T_m, v_m)$ pairs with linear interpolation and flat extrapolation. The constructor validates that the temperatures are strictly increasing (error otherwise). For a query $T_q$:

\[p(T_q) = \begin{cases} v_1 & T_q \le T_1,\\[1ex] v_m + \dfrac{T_q - T_m}{T_{m+1}-T_m}\,(v_{m+1}-v_m) & T_m \le T_q < T_{m+1},\\[2ex] v_N & T_q \ge T_N. \end{cases} \tag{4.8}\]

The interval search is linear for $N\le 8$ (fewer comparisons than a binary search on small tables) and binary for $N>8$ (logarithmic). Under AD the derivative is the slope of the segment containing $T_q$ (the right-hand segment at a knot) and zero in the flat extrapolation; it jumps at every knot.

4.4.6 CallableProperty

Wraps an arbitrary callable $f(T)$. Its value is $p(T)=f(\tilde T)$ on the clamped argument of the property bracket; its derivative is whatever AD propagates through $f$, so $f$ must accept ForwardDiff.Dual arguments (no Float64 annotations or conversions inside). This is the type produced by to_property(f) for a function f and is how the ideal-gas density default (Eq. 4.5) is stored.

4.4.7 SutherlandProperty

Sutherland's law for a gas transport property, carrying its own three constants:

\[p(T) = \mu_{\text{ref}} \left(\frac{T}{T_{\text{ref}}}\right)^{3/2} \frac{T_{\text{ref}}+S}{T+S}.\]

SutherlandProperty() is the canonical air triple ($\mu_{\text{ref}}=1.716\times10^{-5}$ Pa·s, $T_{\text{ref}}=273.15$ K, $S=110.4$ K) and is the default Material.gas_viscosity (§5.9), so a material that says nothing about its pore gas uses air. Its derivative, as AD returns it, is

\[\frac{\mathrm{d}p}{\mathrm{d}T} = p(T)\left(\frac{3}{2T} - \frac{1}{T+S}\right),\]

(the AD-cleanliness of the type is tested in test/materials/test_gas_viscosity_property.jl).

Like every property law it evaluates on the clamped argument $\tilde T$ of the property bracket. The positivity floor MU_FLOOR is applied once, centrally, in effective_gas_viscosity(T, material) (§5.9), so that every curve a material may carry is protected rather than just this one.

4.4.8 StateDependentProperty

A property that depends on both temperature and the full concentration tuple, $p(T,\xi)=f(T,\xi)$. Calling it with $T$ alone raises an error (state is mandatory). Its canonical use is a moisture-dependent heat of sorption: the heat released/absorbed when water desorbs depends on the current moisture content, so

\[h(T,\xi) = f\!\left(T,\; \mathrm{MC}(\xi)\right),\]

where $\mathrm{MC}$ is, e.g., the moisture content reconstructed from the concentration tuple. The kinetics evaluator always passes $\xi$ to the heat property (§4.5, Chapter 6), so a StateDependentProperty heat is evaluated with the live composition.

Its derivatives in $T$ and in each $\xi_j$ are propagated by AD through $f$ (§4.5), so $f$ must accept a Dual temperature and a tuple of Dual concentrations.

4.4.9 to_property conversions

The convenience converter normalizes user inputs to property functions, so that component constructors and reaction constructors accept either property objects or bare values:

Implementation note. All property types, their call methods, linear_interpolate, and to_property are in src/Materials/property_functions.jl.

4.5 Property derivatives and their role in the Jacobian

Implicit time integration (Chapter 15) requires the Jacobian $\partial(\text{RHS})/\partial u$. Because effective properties depend on $T$ and $\xi$, and the RHS depends on effective properties, the per-component derivatives $\mathrm{d}p_j/\mathrm{d}T$ and $\partial p_j/\partial \xi_k$ propagate through the chain rule into the structured Jacobian (Chapter 5, §5.15 for the mixing-rule chain; Chapter 15 for assembly). The property types carry no derivative methods. Each Jacobian operator either runs ForwardDiff over a stencil-local kernel that evaluates the property laws on Dual arguments, or implements an analytical override tested against a ForwardDiff reference (Chapter 15). Either way the property derivative is the derivative of the law as written.

4.5.1 Temperature derivative $\mathrm{d}p/\mathrm{d}T$

What forward-mode AD returns for each type, inside the property bracket:

These values are exact derivatives of the evaluated law, with no truncation error. A TableProperty has a discontinuous derivative at every knot, which the Newton iteration sees as a kink; a smooth law (polynomial, Arrhenius, Sutherland) avoids it.

4.5.2 Concentration derivative $\partial p/\partial \xi_j$

For every property type except StateDependentProperty, the two-argument call p(T, ξ) ignores $\xi$, so $\partial p/\partial \xi_j = 0$ exactly. For a StateDependentProperty the derivative is $\partial f/\partial \xi_j$, propagated by AD through f(T, ξ) with the concentration tuple carrying Dual entries. The kinetics evaluator passes the live $\xi$ to every heat-of-reaction property, so a moisture-dependent heat of sorption contributes its composition derivative to the reaction-source block of the Jacobian.

Implementation note. The call methods (p)(T) and (p)(T, ξ) are in src/Materials/property_functions.jl; the per-cell reaction-source Jacobian that differentiates the heats is ReactionSourceOperator (src/Jacobian/operators/reaction_source.jl).

4.5.3 The reaction-rate connection

The Arrhenius derivative (Eq. 4.7) is identical in form to the temperature-derivative of the reaction-rate constant $k(T)=A\exp(-E/(R_g T))$, since $\mathrm{d}k/\mathrm{d}T = k\,E/(R_g T^2)$; the reaction-source Jacobian obtains it by ForwardDiff over the per-cell kernel (Chapter 15). The full reaction-rate derivative including the depletion limiter and temperature gates is developed in Chapter 6.

4.6 Material assembly

4.6.1 The Material type

A Material is the static, fully type-specialized aggregate that every physics kernel reads. It is immutable and heterogeneously typed so the compiler can specialize every loop over its components and reactions:

struct Material{NC,NR,C<:Tuple,RT<:Tuple,LSL,Pμ} <: AbstractMaterial
    name::Symbol
    components::C          # heterogeneous component tuple
    reactions::RT          # heterogeneous reaction tuple
    mixing::NamedTuple{(:k,:λ,:κ),Tuple{MixingSpec,MixingSpec,MixingSpec}}
    lateral_shrinkage_law::LSL   # A(t)/A₀ = law(χ̄), or nothing
    depletion_limiter::DepletionLimiter   # tanh(ξ/threshold) rate roll-off (§6.3)
    gas_viscosity::Pμ      # pore-gas μ(T), default SutherlandProperty() (§5.9)
end

The type parameters carry compile-time information:

• NC — number of components, a static integer driving every @generated unroll (Val(NC)).
• NR — number of reactions, likewise static.
• C, RT — the exact concrete tuple types of the components and reactions. Because these are full tuple types (not Vector{AbstractComponent}), iterating over them is type-stable and allocation-free; the cost is that each distinct material shape is a distinct compiled type.
• LSL — the lateral-shrinkage-law type, a concrete Function subtype or Nothing. Parameterizing it here lets calls through material.lateral_shrinkage_law specialize statically rather than dispatch through a Union{Function,Nothing} (Chapter 10). The inner constructor enforces LSL === Nothing || LSL <: Function.
• Pμ — the type of the pore-gas viscosity property function, parameterized for the same reason so material.gas_viscosity(T) stays monomorphic. The inner constructor probes the curve at 200–1800 K and rejects one that is not positive and finite there (§5.9).

The mixing field bundles three MixingSpecs — one each for thermal conductivity ($k$), gas transfer ($\lambda$), and permeability ($\kappa$) — in a fixed-order NamedTuple so that material.mixing.k compiles to a direct field read. Each MixingSpec carries (rule, β, ℓ): the mixing rule, the WEIGHTED blend weight $\beta$ (default 0.5), and the Carman–Kozeny characteristic length $\ell$ (NaN unless the rule is CARMAN_KOZENY). The default channel rules are

\[k:\ \text{WEIGHTED},\ \beta=0.5;\qquad \lambda:\ \text{PARALLEL};\qquad \kappa:\ \text{WEIGHTED},\ \beta=0.5,\]

built by MaterialMixing(). The MixingSpec constructor validates $0\le\beta\le1$ and, for CARMAN_KOZENY, requires a finite positive length_scale, failing fast otherwise. The semantics of each rule (PARALLEL, SERIES, WEIGHTED, BRUGGEMAN, CARMAN_KOZENY), the BRUGGEMAN rejection for $\lambda$ and $\kappa$, and the Newton solve for BRUGGEMAN, are all developed in Chapter 5.

The depletion_limiter field carries the material's own DepletionLimiter (default threshold = 1.0 kg·m⁻³, enabled = true), set via the constructor's depletion_limiter keyword and read by evaluate_reactions! for every reaction-rate evaluation (§6.3).

Implementation note. Material, its inner/outer constructors, _construct_material, MixingRule, MixingSpec, MaterialMixing, in src/Materials/materials.jl. Counts via n_components, n_reactions.

4.6.2 PendingReaction and symbolic resolution

Reactions may be written against component symbols rather than integer indices, which is far less error-prone for the user:

Reaction(:decomposition, :virgin => :MMA;          A=8.5e12, E=188e3, h=870e3)
Reaction(:charring,      :virgin => (:char, :gas); A=1e15,  E=2e5, h=5e5, yields=(0.2,0.8))

A symbolic constructor produces a PendingReaction, which stores the symbol payloads plus the kinetic parameters. When Material(...) is called, it builds a name → index map from the components tuple (erroring on a duplicate component name) and resolves every PendingReaction into a numeric Reaction. The resolution path is taken only if at least one reaction is pending — the common all-numeric case skips the Dict allocation entirely:

needs_resolution = any(rxn isa PendingReaction for rxn in reactions)
needs_resolution || return reactions          # fast path, no Dict

Numeric reactions (built against integer indices directly) bypass resolution.

Implementation note. PendingReaction, the symbolic Reaction constructors, resolve_reaction, and _resolve_reactions (called from the Material constructor) in src/Materials/reactions.jl and src/Materials/materials.jl.

4.6.3 The single-reactant mass-balance convention

The Reaction type is n-ary in both reactants and products in principle, but the current model supports only single-reactant reactions ($N_\mathrm{in}=1$). The mass-balance convention (developed fully in Chapter 6) is:

• The lone reactant is consumed with implicit mass stoichiometry $-1$ per unit reaction extent.
• Each product $j$ is produced with an explicit mass yield $\nu_j > 0$.
• Mass conservation requires $\sum_{\text{products}} \nu_j = 1$.

The Reaction inner constructor enforces this at build time — every construction path, including the bare positional constructor, runs the validators. Reaction orders must be $n \ge 0$ (_validate_reaction_order); every product yield must be $\ge 0$ (_validate_yield, enforced unconditionally — a negative yield consumes a species the rate law does not see and drives its concentration negative without bound; note $(1.5, -0.5)$ sums to 1 and would pass a sum-only check); and for a single-reactant reaction the product yields must sum to 1 (_validate_stoichiometry, relative tolerance $\max(10^{-10}\cdot(1+\Sigma\nu), 10^{-14})$), raising an error on violation. Passing validate_mass=false disables only the sum check (not recommended); the non-negativity and order checks always run.

Implementation note. _validate_stoichiometry, verify_stoichiometry (tol $10^{-10}$, returns false + @warn on failure), and mass_balance_error in src/Materials/reactions.jl.

4.6.4 Dry-pyrolysis-gas identification

For variable-cross-section (lateral shrinkage, Chapter 10) the pyrolysis-progress state $\chi$ counts only dry pyrolysis gas — water vapor is excluded because moisture evaporation should not drive radial shrinkage. The set of dry-gas component indices is computed at compile time by a recursive tuple walk that includes a component index iff it is a gas and its name is not a water-vapor name. The water-vapor test is a substring/equality heuristic on the component name ("H2O", "water"/"Water", or exactly "h2o" case-insensitively):

is_gas(c) && !_is_water_vapor_name(c.name)  ⇒ include index

This name-based identification is a documented heuristic (open issue §4.9); a gas whose name does not contain a recognized water token is treated as dry pyrolyzate regardless of its actual chemistry.

Implementation note. dry_pyrolysis_gas_indices → _collect_dry_pyrolysis_gas_indices and the heuristic _is_water_vapor_name, src/Materials/materials.jl. The result is a statically known tuple of indices, constant-folded by the compiler.

4.6.5 Type-stable iteration over heterogeneous tuples

The defining performance pattern of the Materials subsystem is the use of @generated functions and recursive tuple decomposition to unroll every loop over the (heterogeneous) component and reaction tuples at compile time. A naive for c in material.components over a heterogeneous tuple would box each element and allocate; instead, e.g., the mixture heat capacity is

@generated function _sum_heat_capacity(ξ, T, comps::Tuple, ::Val{NC}) where {NC}
    expr = quote
        result = zero(_prop_type(ξ, T))
    end
    for i = 1:NC
        push!(expr.args, :(result += ξ[$i] * comps[$i].heat_capacity(T)))
    end
    push!(expr.args, :(return result))
    return expr
end

which emits $N_C$ statically dispatched terms with no runtime loop. The same idiom appears in evaluate_reactions!, the species-source and gas-flux kernels, the mixing-rule sums (Chapter 5), and the porosity sums. The net effect is that the entire property-evaluation loop in update_properties! is allocation-free once the cache is sized (Chapter 5, §5.13): all composition-dependent dispatch is generated, all intermediate tuples are statically sized, and there are no dictionary lookups in hot paths.

4.7 Time- and temperature-parameter evaluation

Several kernels accept "flexible" parameters that may be a constant, a function of time, or a property function — for example a boundary heat flux $q(t)$, a ramped convective coefficient, or a temperature-dependent property given as a property object. Two thin, AD-preserving dispatch helpers normalize these.

Time-dependent parameters (evaluate_param(param, t)):

\[\texttt{evaluate\_param}(p, t) = \begin{cases} \texttt{convert}(\mathrm{typeof}(t),\,p) & p::\text{Real},\\ p(t) & p::\text{Function}. \end{cases}\]

The convert(typeof(t), p) for a constant is the crucial detail: when $t$ is a ForwardDiff.Dual (during a sensitivity solve, Chapter 17), a bare Float64 constant is promoted to the Dual type so the computation stays in the AD number type end to end. Returning the raw Float64 would break the Dual propagation.

Temperature-dependent parameters (evaluate_param_T(param, T)):

\[\texttt{evaluate\_param\_T}(p, T) = \begin{cases} \texttt{convert}(\mathrm{typeof}(T),\,p) & p::\text{Real},\\ p(T) & p::\text{Function},\\ p(T) & p::\text{AbstractPropertyFunction}. \end{cases}\]

with the same Dual-preserving convert for constants. A dedicated fast-path method short-circuits the most common case:

@inline evaluate_param_T(param::ConstantProperty, T) = param.value

This bypasses the param(T) dynamic-dispatch hop and reads the value field directly — measurable in tight loops where the same constant property is queried per cell per step.

Implementation note. evaluate_param, evaluate_param_T, and the ConstantProperty fast path in src/Materials/property_functions.jl.

4.8 Comparison to ThermaKin2Ds, Gpyro, and FDS

Pyrolysis.jl's material model is a deliberate synthesis of the three reference condensed-phase pyrolysis codes. We map their notation to ours on first use.

4.8.1 Component / species representation

ThermaKin2Ds represents a material as a set of components with mass concentrations as the primary state — the convention Pyrolysis.jl adopts wholesale ($\xi_j$ in kg·m⁻³ is the species state, not mass fraction). ThermaKin distinguishes condensed from gaseous components and applies the solid-volume-fraction convention to conductivity and emissivity — exactly Pyrolysis.jl's solid_volume_fractions (Chapter 5, §5.1), which excludes non-swelling gases (ThermaKin: "when $\gamma=0$, gaseous components do not contribute to the thermal conductivity"). ThermaKin's stoichiometric coefficient $\theta_i^{\,j}$ maps to our $\nu_{i,j}$ (overload note A2).

Gpyro uses a control-volume formulation with $\bar\rho\,\Delta z$ (mass per unit area) and mass fractions $Y_i$ of condensed-phase species; its gas species are tracked separately. Gpyro's volume-fraction property averaging $X_\alpha$ maps to our $v_j$, and its power-law temperature dependence of properties is representable here as a PolynomialProperty or CallableProperty. Gpyro's solid-fraction $\mathrm{SF}_k$ / conversion $\chi_k$ mechanism for selecting non-charring / charring / intumescent / shrinkage behavior corresponds to our swelling-factor $\gamma_j$ plus the lateral-shrinkage law (Chapter 10).

FDS (solid phase) uses a layered material model with per-material density, $c_p$, $k$, and emissivity, plus Arrhenius reaction parameters; its pre-exponential $Z$ maps to our $A_i$ (overload note A3) and its solid absorption coefficient $\kappa_s$ maps to our volumetric $\alpha$ (overload note A5). FDS's material density-ratio swell/shrink rule corresponds to our $\gamma_j$-weighted volume change (Chapter 10).

4.8.2 Heat capacity and the gas-storage convention

ThermaKin's energy equation (its Eq. 17) sums the volumetric heat capacity over all components, gas included: $\rho c_p^{\,\text{ThermaKin}} = \sum_{j\in\text{all}} \xi_j c_{p,j}$. Pyrolysis.jl uses the same sum in the temperature-storage term (Chapter 3):

\[\rho c_p^\mathrm{eff} = \sum_{j=1}^{N_C} \xi_j\, c_{p,j}(T), \tag{4.9}\]

because the pore gas shares the matrix temperature and its concentration is an integrated state, so its sensible storage is available at no cost (Chapter 5, §5.4 and Chapter 7). Gpyro's default thermal-equilibrium form drops the gas term (restored by its FULL_QSG option) and FDS has no in-depth gas; both are approximations of (4.9) that Pyrolysis.jl does not make. The gas flux relative to the skeleton enters the energy equation only through the advective source $S_\mathrm{conv}$.

Implementation note. effective_heat_capacity → _sum_heat_capacity (every component, no phase branch), src/Properties/effective_properties.jl.

4.8.3 Reaction kinetics representation

ThermaKin supports bimolecular rates $\propto \xi_{\text{COMP1}}\xi_{\text{COMP2}}$; Gpyro uses an nth-order-in-conversion rate with optional O₂ sensitivity $g(Y_{O_2})$ and a critical-temperature Heaviside; FDS uses $\rho^n\cdot Y_{O_2}^{n}\cdot$ a temperature factor with kinetic-parameter recovery from a TGA reference temperature and peak rate. Pyrolysis.jl restricts to single-reactant mass-action kinetics in mass-concentration basis with smooth tanh temperature gates and a smooth depletion limiter (Chapter 6). The representation of heats of reaction as property functions — including the state-dependent moisture-sorption heat via StateDependentProperty — generalizes all three codes' scalar heats. The full kinetic comparison is in Chapter 6.

4.8.4 Property functions and temperature dependence

ThermaKin and Gpyro both permit temperature-dependent properties (tabulated or power-law). Pyrolysis.jl's property-function abstraction subsumes both: a Gpyro power law $p(T)=p_0(T/T_0)^m$ is a CallableProperty; a ThermaKin table is a TableProperty; FDS ramp-style piecewise data is a TableProperty with linear interpolation and flat extrapolation. The novel addition is the uniform $(T,\xi)$ call signature, which allows composition-dependent properties (heat of sorption) without a special-case code path.

4.9 Limitations and open issues

The following are documented limitations of the constitutive layer; several are revisited where the relevant physics is developed.

1. Single-reactant reactions only. $N_\mathrm{in}=1$ is enforced. Multi-step schemes must be modeled as sequences of single-reactant reactions. Multi-reactant (e.g. genuine bimolecular char oxidation) would require extending the reactant tuple to 3-ary (idx, order, mass_stoich) and a revised mass-balance convention (Chapter 6).

2. Intrinsic-density fallback may be unphysical. When $\rho_{i,j}$ is not supplied it defaults to the bulk $\rho_j$, asserting the pure phase has no internal pores. For genuinely microporous solids this under-predicts porosity and over-predicts the matrix volume.

3. Name-based dry-gas identification. Water vapor is excluded from the pyrolysis-progress source by a substring match on the component name (_is_water_vapor_name). A water species named outside the recognized tokens would be (incorrectly) counted as dry pyrolyzate, and a non-water species named with a water token would be (incorrectly) excluded.

4. Kinked and AD-dependent property laws. A TableProperty is differentiated segment by segment, so its derivative jumps at every knot and the Newton iteration sees a kink there. CallableProperty and StateDependentProperty are differentiated by propagating ForwardDiff.Dual numbers through the user's function, which must therefore be generic in its argument types (§4.5).

5. Heterogeneous-tuple compilation cost. Preserving the exact component and reaction tuple types yields allocation-free, type-stable hot paths, but at the cost of recompiling the residual for each distinct material shape. Materials that differ only in numerical parameters share a compiled type; materials that differ in the types of their property functions do not.

6. No temperature-dependent pre-exponential or pressure-dependent kinetics. The Arrhenius property and Reaction carry constant $A$ and $E$ (no Kissinger/ Ozawa $T$-dependent activation energy, no distributed-activation-energy tooling, no explicit pressure dependence). See Chapter 6.

7. Heat-of-reaction sign mapping to reference codes. Internally $h>0$ is endothermic ($Q_\mathrm{rxn}=-\sum_r h_r r_r$), whereas ThermaKin/Gpyro publish $h>0$ as exothermic. Users porting parameters from those codes must flip the sign (Chapter 2, overload note H1; Chapter 6).