5. Boundary and Initial Conditions

Boundary conditions (BCs) tell the solver how heat and gas enter and leave the 1-D column; initial conditions set the temperature and composition at t = 0. This chapter shows how to build a BoundaryConditionSet, the complete catalogue of thermal, mass, pressure, and mesh-motion BCs with every field and default, how the surface energy balance is solved, and how to specify initial temperature and concentration profiles. The physics and derivations behind each BC live in the Technical Reference §12 (and the radiation models in §8); this chapter is about how to specify them correctly.

Coordinate and sign reminders (see Technical Reference §2). The axial coordinate is z. z = 0 is the bottom / substrate (tag :bottom, boundary id 2); z = L is the top / exposed surface (tag :top, boundary id 1). Heat normally enters from the top. For all BCs in this package the sign convention is: positive flux = into the domain (heating, or mass entering). A gas escaping out of the top therefore has a negative boundary mass flux in the BC's own convention.

5.1 The BoundaryConditionSet

All boundary conditions are bundled in a single tag-keyed container, BoundaryConditionSet. It holds exactly two kinds of BC — thermal and mass — each as a NamedTuple keyed by boundary tag (:top, :bottom in 1-D):

using Pyrolysis

bc_set = BoundaryConditionSet(
    thermal = (top    = HeatFluxBC(50e3) + RadiativeBC(ε = 0.9, T_inf = 300.0),
               bottom = AdiabaticBC()),
    mass    = (top    = DiffusiveMassBC(0.01),
               bottom = ImpermeableBC()),
)

Both keyword arguments default to fully-insulated, impermeable boundaries:

# Equivalent to the explicit defaults below
bc_set = BoundaryConditionSet()
#   thermal = (top = AdiabaticBC(),   bottom = AdiabaticBC())
#   mass    = (top = ImpermeableBC(), bottom = ImpermeableBC())

Key facts about the constructor:

• Partial specification is allowed. Whatever you supply is merged over the defaults, so BoundaryConditionSet(thermal = (top = HeatFluxBC(50e3),)) leaves the bottom adiabatic and both mass faces impermeable. (Note the trailing comma — a one-element NamedTuple needs it.)
• Type checking is enforced at construction. Passing a mass BC into a thermal slot (or vice-versa) throws an ArgumentError immediately.
• It is type-stable. Each tag compiles to a direct field read via Val{tag} dispatch, so the residual specializes statically with zero runtime overhead.

The set is then handed to PyrolysisProblem via the bc_set keyword:

problem = PyrolysisProblem(
    mesh = mesh,
    material = material,
    bc_set = bc_set,
    T_initial = z -> 300.0,
    ξ_initial = [z -> 1190.0, z -> 0.0],
    tspan = (0.0, 100.0),
)

Note on mesh motion and radiation BCs. The BoundaryConditionSet has only thermal and mass fields. ALE mesh motion needs no BC at all — it is computed automatically from the physics (§5.7) — and the P1 radiation BCs (§5.6) are not stored in this container: they live in a separate RadiationBCSet inside the Experimental module. External incident radiation (cone heater, panel) is specified as a thermal BC (RadiativeFluxBC, §5.2.6), not as a radiation BC.

5.2 Thermal boundary conditions

Thermal BCs are subtypes of AbstractThermalBC. They are the heat-side input to the surface energy balance that the solver closes by a Newton iteration for the surface temperature T_s (§5.4). Every flux is reported with the positive = into domain convention.

All exported by using Pyrolysis.

5.2.1 AdiabaticBC — perfect insulation

bc = AdiabaticBC()

Zero heat flux. This is the default for both faces and the usual choice for the substrate (back face) in a cone-calorimeter idealization where the sample is mounted on insulation.

5.2.2 HeatFluxBC — prescribed flux

bc = HeatFluxBC(50e3)                            # constant 50 kW/m² into the surface
bc = HeatFluxBC(t -> 50e3 * (1 - exp(-t / 10)))  # ramp-up over ~10 s

The single field, flux, is the heat flux q_BC in W/m², positive into the domain. It may be a constant (Real) or a function of time t. Use this for direct-contact heating or for an opaque surface where you want to inject a flux directly (rather than via the radiation models, see RadiativeFluxBC). dq_BC/dT_s = 0, so it does not couple to the surface temperature.

5.2.3 ConvectiveBC — Newton cooling

bc = ConvectiveBC(h = 10.0, T_inf = 300.0)
bc = ConvectiveBC(h = 10.0, T_inf = t -> 300 + 0.1 * t)   # warming ambient

Newton's cooling law q_BC = h·(T_∞ − T_s). Fields:

• h — convective heat-transfer coefficient h_conv in W·m⁻²·K⁻¹. Constant or h(t). Validated to be non-negative.
• T_inf — ambient temperature T_∞ in K. Constant or T_∞(t). Validated to be positive.

The flux is positive (heating) when T_∞ > T_s and negative (cooling) when the surface is hotter than the surroundings — the usual re-radiation/convection loss direction during pyrolysis. Its temperature derivative dq_BC/dT_s = −h makes the Newton solve converge quickly.

Blowing. Strong outgassing thickens the boundary layer and suppresses convective exchange, and a gas denser than the free stream can pool on the face. How a ConvectiveBC responds is its blowing model, NoBlowing() by default (the coefficient stays h whatever leaves the face). The surface balance hands the model the outward mass flux of every gas component at the face at each iterate, so the back face is unaffected automatically (no outflow). Four models ship:

ConvectiveBC(h = 15.0, T_inf = 300.0, blowing = CouetteFilm())        # h·B/(e^B − 1), B = ṁ″·c_p,g/h
ConvectiveBC(h = 11.2, T_inf = 360.0, blowing = CouetteFilm(0.31))    # film law with a fitted constant
ConvectiveBC(h = 11.2, T_inf = 360.0, blowing = PropertyRatioFilm(NITROGEN))
ConvectiveBC(h = 11.2, T_inf = 360.0, blowing = PoolingFilm(NITROGEN; volatiles = pyrolysate))
ConvectiveBC(h = 11.2, T_inf = 360.0, blowing = FunctionBlowing((face, h, T_inf) -> 0.8))

• CouetteFilm(λ = 0.5): the classical film law h → h·φ/(e^φ − 1), φ = 2λ·ṁ″·c_p,g/h with the heat capacity of the gas leaving the face. λ = 0.5 is the Couette film (φ = B). At a wood-like peak ṁ″ ≈ 10 g·m⁻²·s⁻¹ against h ≈ 10 this is a 35–50 % reduction of h; at PMMA-like peaks (B ≈ 4) convective cooling is essentially extinguished. A smaller λ fitted to a resolved flow holds for that gas only.
• PropertyRatioFilm(free_stream; n_Γ = 0.6): the film law for any light gas, φ = (ṁ″·c_p,∞/h)·(Γ_∞/Γ_g)^n_Γ with Γ = k/c_p of the free stream and of the leaving mixture (heat capacity by mass fraction, conductivity by Wilke's rule). No constant is refitted for a new gas.
• PoolingFilm(free_stream; n_Γ, n_D, ψ0, Λc, Λw): the same with a pooling factor ψ, h → h·φ/(e^{φ/ψ} − 1), that falls from 1 to the floor ψ0 when the gas at the face is denser than the free stream. The face composition comes from a species-by-species film balance with the binary diffusivities into the free stream. The three switch constants belong to the apparatus flow (the defaults are the CAPA chamber's); the gas enters through its properties. Describe the volatiles either through the material's gas components, with Fuller diffusion volumes diffusion_volume = (name = V, ...), or as one lumped GasSpecies (volatiles = ...).
• FunctionBlowing(f): your own rule χ = f(face, h, T_inf) from the BlowingFace (outward species fluxes face.ṁ, surface temperature face.T_s, face.material).

The free stream is a GasSpecies (molar mass, c_p(T), k(T), Fuller diffusion volume, optional viscosity); NITROGEN ships, and nasa7_heat_capacity builds a heat capacity from NASA-7 coefficients. Each convective member of a combined BC applies its own model. The solution records the applied ratio and, for the pooling model, ψ, the face fraction and the face buoyancy in extras.surface_blowing_ratio, surface_pooling_factor, surface_face_fraction and surface_face_buoyancy. The shipped examples were calibrated with the bare coefficient (and without flame feedback — the two partially cancel; see Technical Reference §1.8.1 before recalibrating).

5.2.4 RadiativeBC — surface re-radiation

bc = RadiativeBC(T_inf = 300.0)                  # uses material's effective emissivity
bc = RadiativeBC(ε = 0.9, T_inf = 300.0)         # fixed emissivity
bc = RadiativeBC(T_inf = 300.0, F = 0.7)         # partially enclosed surface
bc = RadiativeBC(ε = T -> 0.85 + 0.1*(T-300)/700, T_inf = 300.0)  # T-dependent ε

Gray-body radiative exchange with the surroundings, q_BC = F·ε_eff·σ·(T_∞⁴ − T_s⁴), where σ = 5.670374419×10⁻⁸ W·m⁻²·K⁻⁴. Fields:

• ε — surface emissivity. Default nothing, which means use the effective emissivity ε_eff = Σ_j v_j ε_j(T) (volume-weighted, each component emissivity evaluated at the boundary-cell temperature) computed from the live composition. This is the recommended setting: it automatically tracks the change from virgin polymer to char as the material decomposes, and any temperature dependence you gave the component ε. You can also pass a constant in [0, 1] or a function ε(T).
• T_inf — surroundings temperature T_∞ in K (constant or T_∞(t)).
• F — view factor in [0, 1], default 1.0 (surface fully sees the surroundings). Use F < 1 for partially enclosed geometries.

Its temperature derivative is dq_BC/dT_s = −4·F·ε_eff·σ·T_s³. When ε is a function of temperature, the solver neglects the dε/dT term in this derivative; this slightly slows Newton convergence but does not change the converged root (see Technical Reference §12.2).

5.2.5 DirichletThermalBC — prescribed temperature

bc = DirichletThermalBC(500.0)            # hold this face at 500 K
bc = DirichletThermalBC(t -> 300 + 5*t)   # ramped wall temperature

Forces the surface to a prescribed temperature T = T_s. This BC bypasses the ordinary surface-balance Newton iteration and is instead enforced with a large penalty conductance (h_large = 1e8 W·m⁻²·K⁻¹) plus a special case in the RHS assembly. Use it for verification problems or when an experiment fixes the surface temperature directly. The single field T may be constant or T(t).

5.2.6 RadiativeFluxBC — external incident radiation

This is the BC for a cone heater, furnace, or radiant panel — an external radiation source heating the surface. It is distinct from RadiativeBC (which is the surface's own re-radiation loss) and from HeatFluxBC (which injects a bare flux with no absorptivity or in-depth treatment).

bc = RadiativeFluxBC(flux = 25e3)                          # incident 25 kW/m², α = a_eff
bc = RadiativeFluxBC(flux = 50e3, absorptivity = 0.9)      # fixed absorptivity
bc = RadiativeFluxBC(flux = t -> 25e3*(1 - exp(-t/10)))    # time-varying flux

Fields:

• flux — incident radiative flux Φ_incident in W/m² (constant or Φ(t)).
• absorptivity — surface absorptivity α in [0, 1]. Default nothing, which uses the material's effective surface absorptivity a_eff = Σ_j v_j a_j: each component's own absorptivity when set, else its emissivity ε_j(T) (Kirchhoff's law). With no component absorptivity set this is exactly ε_eff; setting one lets the absorbed heater flux differ from the re-radiation without splitting the BC. Either way it tracks composition and temperature automatically. You can override with a constant or α(t).

Its behavior depends on the radiation_model keyword passed to solve() (see §7.4 and Technical Reference §8):

The boundary flux from this BC has dq_BC/dT_s = 0 (it is the incoming source and does not depend on the surface temperature).

Do not double-count radiation. Do not combine RadiativeFluxBC and HeatFluxBC for the same radiation source — that adds the energy twice. RadiativeFluxBC is the correct, absorptivity-aware choice for a radiant source.

5.2.7 CombinedThermalBC — adding BCs with +

A real exposed surface usually has several simultaneous heat exchanges: incoming radiation, re-radiation loss, and convective loss. Add them with +:

top_thermal =
    RadiativeFluxBC(flux = 50e3, absorptivity = 0.96) +  # cone heater (incoming)
    RadiativeBC(ε = 0.96, T_inf = 330.0) +               # re-radiation loss
    ConvectiveBC(h = 10.0, T_inf = 330.0)                # convective loss

The combined flux is the sum of the component fluxes, q_BC = Σ_i q_BC,i(T_s, t), and the derivative is summed likewise — assembled by a generated function for type stability. You can chain any number of thermal BCs. Each term keeps its own T_inf (§5.8).

A canonical cone-calorimeter top face:

bc_set = BoundaryConditionSet(
    thermal = (
        top = RadiativeFluxBC(flux = 50e3, absorptivity = 0.96) +
              RadiativeBC(ε = 0.96, T_inf = 330.0) +
              ConvectiveBC(h = 10.0, T_inf = 330.0),
        bottom = AdiabaticBC(),
    ),
    mass = (
        top    = ConvectiveMassBC(h_m = 0.05, ξ_inf = 0.0),
        bottom = ImpermeableBC(),
    ),
)
# Pair with:  solve(problem; radiation_model = SURFACE_ABSORPTION)

5.3 Mass boundary conditions

Mass BCs (subtypes of AbstractMassBC) control how gas escapes (or enters) through a boundary. Solid and liquid components never transport across the boundary; only gaseous components do. The sign convention is positive = mass into the domain, so gas escaping at the top has a negative flux in the BC's own convention (the solver applies the correct face orientation internally).

Pressure BCs (§5.5) are also AbstractMassBC subtypes but only act in Darcy mode. The first four BC types above (and the pressure BCs) are exported by using Pyrolysis; CombinedMassBC is produced by the + operator rather than called directly.

5.3.1 ImpermeableBC — no gas escape

bc = ImpermeableBC()

Zero mass flux for every component. This is the default for both faces and the usual choice for the substrate.

5.3.2 MassFluxBC — prescribed mass flux

bc = MassFluxBC(0.0)             # constant prescribed flux [kg/(m²·s)]
bc = MassFluxBC(t -> 1e-3)       # time-varying

Prescribes the boundary mass flux J directly (positive = into domain), as a constant or J(t). The same value applies to all gas components.

5.3.3 ConvectiveMassBC — external convective transfer

bc = ConvectiveMassBC(h_m = 0.05, ξ_inf = 0.0)   # gas swept away to clean atmosphere

External convective mass transfer, J_j = h_m·(ξ_{j,∞} − ξ_{j,s}). Fields:

• h_m — external (film) mass-transfer coefficient h_m in m/s (velocity-based). Constant or h_m(t); validated non-negative. Typical values: natural convection 0.01–0.05 m/s, forced convection 0.05–0.5 m/s. As a rule of thumb from the heat–mass analogy, h_m ≈ h_conv / (ρ_air · c_{p,air}); for h_conv ≈ 10 W·m⁻²·K⁻¹ this gives h_m ≈ 0.008 m/s.
• ξ_inf — ambient mass concentration ξ_∞ in kg/m³, default 0.0 (gas escapes into clean surroundings). Accepts a scalar (broadcast to every gas species), a function of time t -> ξ_∞(t), or an NTuple with one entry per component (ξ_inf = (0.0, 0.0, 0.27) — e.g. ambient air ingress for a background-air species without injecting pyrolyzate). Entries are validated non-negative; a tuple's length must equal the material's component count.

5.3.4 DiffusiveMassBC — internal diffusion with optional external resistance

bc = DiffusiveMassBC()            # pure internal diffusion (h_m = Inf, default)
bc = DiffusiveMassBC(0.01)        # add a finite external film resistance, h_m = 0.01 m/s

Couples the internal Fickian diffusion gradient in the surface cell to an optional external film resistance. The single field is:

• h_m — external mass-transfer coefficient in m/s. Default Inf (the no-argument constructor DiffusiveMassBC() sets h_m = Inf), meaning no external resistance — the flux is the pure diffusive limit J_j = −λ_j ξ_{j,s}/d, where d is the distance from the surface cell center to the boundary face. A finite h_m adds a series resistance: J_j = −λ_j h_m ξ_{j,s}/(λ_j + h_m d).

This is the most physically consistent escape BC because it ties the surface loss to the same internal transport law used everywhere else (it also keeps the species balance consistent with ALE surface recession). Use it (with Inf or a finite h_m) when you want gas escape governed by the material's own gas transfer coefficient rather than an imposed convective coefficient.

5.3.5 CombinedMassBC — adding mass BCs with +

bc = DiffusiveMassBC(Inf) + ConvectiveMassBC(h_m = 0.01, ξ_inf = 0.0)

Sums component mass fluxes, J_j = Σ_i J_{j,i}, via the + operator (generated function, type-stable). For most pyrolysis runs a single mass BC suffices; combined mass BCs are for modeling competing transport mechanisms or research/validation studies.

5.4 The surface energy balance (how thermal BCs are applied)

Thermal BCs are not applied as a raw flux in the cell; instead the solver finds the surface (face) temperature T_s that makes the boundary energy balance close, then conducts the resulting flux into the adjacent cell. The balance is

q_BC(T_s) = k_eff · (T_s − T_cell) / Δz_{1/2}

— the applied BC flux equals the heat conducted from the surface face into the first cell center, with k_eff the effective conductivity, T_cell the adjacent cell temperature, and Δz_{1/2} the face-to-center distance. This is solved by a damped Newton iteration on T_s (see Technical Reference §12.2). You do not call this yourself — it runs inside the RHS assembly — but it is worth knowing because it explains a few user-visible behaviors:

• Convergence is robust for physical BCs: convective flux is linear in T_s, radiative is T_s⁴; typical solves converge in 2–5 iterations (default tolerance 1e-10, up to 50 iterations).
• Surface temperature is clamped to [200, 3000] K during the iteration to prevent unphysical excursions.
• The solver is AD-aware: it carries ForwardDiff.Dual types through the Newton iteration so the Jacobian's derivatives with respect to the state flow correctly through T_s.
• The converged surface temperature is available afterwards as solution.surface_T, and the flux decomposition (absorbed, re-radiation, convection, conducted) is reported in solution.extras (§9.3).

The balance has no surface gas-enthalpy term. The enthalpy carried by the escaping gas is already in the energy equation as the volumetric advection source S_conv = −Σ_g c_{p,g} J_g ∂T/∂z; a surface term Σ_g ṁ_g Δh_g(T_s, T_∞) would overlap it by exactly ∫S_conv and count that enthalpy a second time. The outward gas flux enters the boundary only through the blowing model of a ConvectiveBC (§5.2.3).

5.5 Pressure boundary conditions (Darcy mode only)

When the material has a permeability (so the solver runs in Darcy–Fick transport mode, see Technical Reference §9), gas transport is pressure-driven and you can prescribe pressure at a boundary. Pressure BCs are AbstractMassBC subtypes, so they go in the mass slot of the BoundaryConditionSet. They do not contribute a diffusive mass flux; instead they set the Darcy gas velocity at the boundary face (via apply_pressure_bc_to_darcy_velocity!). In Fickian (non-Darcy) mode they have no effect.

All three are exported by using Pyrolysis.

5.5.1 DirichletPressureBC — fixed pressure

bc = DirichletPressureBC()           # default: atmospheric, 101325 Pa
bc = DirichletPressureBC(101325.0)   # explicit atmospheric

Holds the boundary at pressure P (Pa, validated positive). The Darcy velocity is computed from the gradient between the boundary and the interior cell, v_g = −(κ/μ)·(P_bc − P_interior)/d. This is the natural exposed-surface condition: the outer face sits at atmospheric pressure and gas flows out down the pressure gradient built up by pyrolysis.

5.5.2 NeumannPressureBC — prescribed gradient

bc = NeumannPressureBC()        # default: zero gradient (impermeable to gas flow)
bc = NeumannPressureBC(0.0)     # same as above
bc = NeumannPressureBC(1000.0)  # prescribed gradient [Pa/m] forcing an outflow

Sets the normal pressure gradient ∂P/∂n (Pa/m), default 0.0. A zero gradient is a no-flow (impermeable) back face for the pressure field — the usual choice for a sealed substrate in Darcy mode.

5.5.3 ConvectivePressureBC — Robin (resistive) outflow

bc = ConvectivePressureBC(1e-6, 101325.0)   # finite surface resistance, atmospheric exterior

A Robin-type condition, v_g·n = h_P·(P_internal − P_external), for a surface with finite resistance to gas flow (a char skin, a coating). Fields:

• h_P — pressure-transfer coefficient h_P in m·Pa⁻¹·s⁻¹.
• P_ext — external pressure P_external in Pa.

Limits: h_P → ∞ approaches DirichletPressureBC; h_P → 0 approaches an impermeable (zero-flow) boundary.

5.6 Radiation boundary conditions (P1 model, Experimental)

The four radiation models are selected by the radiation_model keyword of solve() (default NO_RADIATION), not by a BC in the set. For NO_RADIATION, SURFACE_ABSORPTION, and BEER_LAMBERT, the incoming radiation is specified by a RadiativeFluxBC in the thermal slot (§5.2.6) — there is nothing else to set.

The P1 quasi-steady model (P1_QUASI_STEADY) is different: it solves a volumetric radiation field and needs its own boundary treatment (Marshak, reflective, or Dirichlet on the intensity G). These P1 BCs are experimental and live in Pyrolysis.Experimental, not in the BoundaryConditionSet:

import Pyrolysis.Experimental: MarshakRadiationBC, ReflectiveRadiationBC,
                               DirichletRadiationBC, RadiationBCSet, P1_QUASI_STEADY

The standard P1 configuration is a Marshak condition on the exposed top (carrying the external flux and ambient field) and a reflective condition on the insulated bottom. The Marshak flux is q_BC = c·(G_ext + 4·I_ext − G_surf) with c = ε/(2(2−ε)) and G_ext = 4σT_∞⁴ — the external flux carries the half-range weight 4c = 2ε/(2−ε), the same weight as the ambient term (see Technical Reference §8.6 for the derivation).

Status caveat. P1 radiation is in Pyrolysis.Experimental and runs through a separate residual path. For most cone-calorimeter and TGA work, SURFACE_ABSORPTION (opaque) or BEER_LAMBERT (semi-transparent) with a plain RadiativeFluxBC is the supported, production route. Reach for P1 only when you specifically need in-depth re-radiation in an optically thick sample.

5.7 Mesh motion in ALE simulations (no boundary conditions needed)

ALE mesh motion is not configured through boundary conditions — there are no mesh-motion BC types, and BoundaryConditionSet carries only thermal and mass slots. The mesh-velocity boundary treatment is fixed by the formulation: the bottom node is pinned (w_1 = 0, hard-coded in the residual's mesh-velocity integration) and every other node's velocity follows cumulatively from the volumetric strain as the material shrinks or swells (see Technical Reference §11.3), so the top surface automatically follows the material (recedes or expands). You enable the moving mesh simply with solve(problem; use_ale = true).

5.8 Choosing ambient temperatures

Two thermal BCs carry an ambient temperature T_∞: ConvectiveBC and RadiativeBC each have their own T_inf, which enters only that BC's own flux. BCs without an ambient (HeatFluxBC, AdiabaticBC, RadiativeFluxBC, DirichletThermalBC) carry none. In a CombinedThermalBC each term keeps its own T_inf: a top BC of RadiativeBC(T_inf = 330) plus ConvectiveBC(T_inf = 300) re-radiates to 330 K and convects to 300 K. The surface energy balance (§5.4) has no ambient reference temperature of its own.

Practical guidance: set T_inf consistently on the loss terms (RadiativeBC, ConvectiveBC) to the experimental surroundings temperature. In a cone calorimeter the cone shroud and chamber are warmer than room temperature; values in the 300–360 K range are common and used in the shipped examples.

5.9 Time-dependent parameters

Any BC parameter shown above as "constant or function of t" accepts a single-argument function. The solver evaluates it at the current time. Examples:

# Heat-flux ramp
HeatFluxBC(t -> 50e3 * min(t / 5.0, 1.0))           # linear ramp to 50 kW/m² over 5 s

# Time-varying incident radiation
RadiativeFluxBC(flux = t -> 25e3 * (1 - exp(-t / 10)))

# Warming ambient in a convective loss term
ConvectiveBC(h = 10.0, T_inf = t -> 300 + 0.5 * t)

# Pulsed mass extraction
MassFluxBC(t -> t < 50 ? 0.0 : -1e-3)

Note on time-varying radiation flux. When radiation_model is SURFACE_ABSORPTION or BEER_LAMBERT, solve() re-samples the incident flux of the top RadiativeFluxBC at every RHS evaluation — constant and time-varying fluxes alike, with no up-front constancy detection. A ramped or modulated flux is therefore always tracked exactly. You do not need to do anything special — just supply flux = t -> ….

5.10 Initial conditions

Initial conditions are supplied to PyrolysisProblem as functions of the spatial coordinate z (in meters, with z = 0 at the bottom, z = L at the top). Two inputs:

• T_initial — a function z -> T returning the initial temperature in K. Default: z -> 300.0.
• ξ_initial — the initial mass concentrations ξ_j(z) in kg/m³, one entry per component. Accepts a Vector of functions (one per component), a single function (broadcast to all components), or you may omit it for the auto-default.

problem = PyrolysisProblem(
    mesh = mesh,
    material = material,            # NC = 2 here: virgin solid + gas
    bc_set = bc_set,
    T_initial = z -> 300.0,                  # uniform 300 K
    ξ_initial = [z -> 1190.0, z -> 0.0],     # full virgin density, no gas
    tspan = (0.0, 100.0),
)

5.10.1 Auto-defaults

If you omit ξ_initial, the constructor builds a sensible default: the first component, if it is a solid, is placed at its reference density (its density function evaluated at 300 K) and all other components are set to zero. This is correct for the common "single virgin solid that decomposes into gas(es) and char" setup, where component 1 is the virgin solid. Always double-check it matches your component ordering; if not, pass ξ_initial explicitly.

T_initial defaults to a uniform 300 K.

5.10.2 Uniform vs. depth-varying profiles

A uniform value is just a constant-returning closure:

T_initial = z -> 300.0
ξ_initial = [z -> 1190.0, z -> 0.0]   # uniform composition

For a non-uniform start (e.g. a pre-charred surface layer, or a moisture gradient) make the closures depend on z:

L = 0.02  # domain thickness [m]
# Warmer at the exposed top, cooler at the substrate
T_initial = z -> 300.0 + 50.0 * (z / L)

# A thin char layer near the surface (component 2 = char), virgin below
ξ_initial = [
    z -> z > 0.9 * L ? 200.0 : 1190.0,   # virgin: depleted near top
    z -> z > 0.9 * L ? 150.0 : 0.0,      # char: present near top
    z -> 0.0,                            # gas: none initially
]

The initial conditions are evaluated at each cell center and stored in the mesh. They also populate the per-cell dry-solid mass denominator used by the lateral-shrinkage (χ) dynamics, computed as m_{dry,i} = V_i · Σ_{solid} ξ_{j,i}(0) (see Technical Reference §12.8 and §10.3).

5.10.3 Profile helper functions

Convenience profile builders are available in the Problem submodule (not top-level exported). Import them explicitly when you want them:

import Pyrolysis.Problem: step_profile, linear_profile, gaussian_profile, tanh_profile

z_int = 0.9 * L
# Sharp step: char above the interface, virgin below
ξ_char_initial = z -> step_profile(z; z_interface = z_int,
                                       value_left = 0.0, value_right = 150.0)

# Linear temperature ramp through the thickness
T_initial = z -> linear_profile(z; z_min = 0.0, z_max = L,
                                    value_min = 300.0, value_max = 350.0)

# Smooth (tanh) interface instead of a hard step
ξ_char_initial = z -> tanh_profile(z; z_center = z_int, width = 0.001,
                                      value_left = 0.0, value_right = 150.0)

# Gaussian bump (e.g. a localized moisture pocket)
ξ_moist_initial = z -> gaussian_profile(z; z_center = 0.5 * L, width = 0.002,
                                           amplitude = 100.0, baseline = 0.0)

Their signatures:

A smooth (tanh) interface is generally preferable to a hard step_profile for the solver: sharp discontinuities at t = 0 can provoke small startup transients.

5.11 Common configurations (worked recipes)

5.11.1 Cone calorimeter (opaque charring solid)

External cone flux absorbed at the surface, with re-radiation and convective losses; gas escapes by convective mass transfer; insulated, impermeable back face:

bc_set = BoundaryConditionSet(
    thermal = (
        top = RadiativeFluxBC(flux = 50e3, absorptivity = 0.96) +
              RadiativeBC(ε = 0.96, T_inf = 330.0) +
              ConvectiveBC(h = 10.0, T_inf = 330.0),
        bottom = AdiabaticBC(),
    ),
    mass = (
        top    = ConvectiveMassBC(h_m = 0.05, ξ_inf = 0.0),
        bottom = ImpermeableBC(),
    ),
)
solution = solve(problem; radiation_model = SURFACE_ABSORPTION)

5.11.2 Semi-transparent polymer (in-depth absorption)

Same incoming radiation, but absorbed in depth via Beer–Lambert (requires absorption coefficients on the components):

bc_set = BoundaryConditionSet(
    thermal = (
        top = RadiativeFluxBC(flux = 50e3) +          # α defaults to a_eff (= ε_eff by Kirchhoff)
              RadiativeBC(T_inf = 300.0) +
              ConvectiveBC(h = 10.0, T_inf = 300.0),
        bottom = AdiabaticBC(),
    ),
    mass = (top = DiffusiveMassBC(), bottom = ImpermeableBC()),
)
solution = solve(problem; radiation_model = BEER_LAMBERT)

5.11.3 Furnace temperature program on the exposed face

A prescribed heating program is most directly modeled by holding the surface temperature with a time-dependent DirichletThermalBC:

β = 10.0 / 60.0   # 10 K/min in K/s
bc_set = BoundaryConditionSet(
    thermal = (top = DirichletThermalBC(t -> 300.0 + β * t),
               bottom = AdiabaticBC()),
    mass    = (top = DiffusiveMassBC(), bottom = ImpermeableBC()),
)

5.11.4 Pressure-driven char (Darcy mode)

With a permeable material, set atmospheric pressure on the exposed face and a sealed (zero-gradient) back face. The transport mode switches to Darcy–Fick automatically when the material has permeability:

bc_set = BoundaryConditionSet(
    thermal = (top = HeatFluxBC(50e3) + RadiativeBC(ε = 0.9, T_inf = 300.0),
               bottom = AdiabaticBC()),
    mass    = (top = DirichletPressureBC(101325.0),    # atmospheric exposed face
               bottom = NeumannPressureBC(0.0)),       # sealed substrate
)

5.12 Caveats and pitfalls

• The escaping gas carries no surface enthalpy term. The volumetric gas-advection source S_conv = −Σ_g c_{p,g} J_g ∂T/∂z is always active and already carries the gas enthalpy (Technical Reference §7.4.3, §12.2); a surface term Σ_g ṁ_g Δh_g(T_s, T_∞) on top of it would double-count that enthalpy, a spurious ~5–10 kW/m² sink that under-predicts surface temperature and mass-loss rate. The S_conv-only convention matches ThermaKin, Gpyro, and FDS. Parameter sets calibrated with a code that adds a surface gas-enthalpy term may need re-calibration.

• radiation_model defaults to NO_RADIATION. A RadiativeFluxBC only does something if you pass an appropriate radiation_model to solve(). With the default NO_RADIATION, an incident RadiativeFluxBC is still absorbed at the surface (the absorbed flux is applied as a boundary heat flux), but choose SURFACE_ABSORPTION (opaque) or BEER_LAMBERT (semi-transparent) explicitly to make the intent clear and to enable the in-depth distribution.

• Only thermal and mass BCs go in the set. BoundaryConditionSet accepts only AbstractThermalBC and AbstractMassBC payloads — a wrong-kind BC throws at construction. ALE mesh motion needs no BC at all (it is automatic with use_ale = true), and P1 radiation BCs (MarshakRadiationBC etc.) go through the Experimental P1 residual path, not the set.

• Pressure BCs need Darcy mode. DirichletPressureBC/NeumannPressureBC/ ConvectivePressureBC only have an effect when the material has permeability (Darcy–Fick transport). In Fickian mode they contribute zero flux and are silently inert — use ImpermeableBC/DiffusiveMassBC/ConvectiveMassBC for diffusion-only gas escape.

• Don't double-specify a radiation source. Use RadiativeFluxBC (not also HeatFluxBC) for the cone/panel; use RadiativeBC for the surface's own re-radiation loss. They are different physical terms.

• Match ξ_initial to your component order. The auto-default only fills component 1 (if solid). For multi-solid or moisture-bearing materials, supply ξ_initial explicitly so virgin solid, char, moisture, and gases start with the right concentrations.

For deeper background on the boundary physics, the surface Newton solve, and the radiation models, see Technical Reference §12 (boundary and initial conditions), §8 (thermal radiation), §7 (energy assembly and the S_conv convention), and §9 (gas transport and pressure).