KWN (precipitation kinetics)¶

Overview¶

The KWN module is a collection of Model primitives for composing Kampmann–Wagner Numerical (KWN) models of precipitation, growth, and coarsening in a multi-component matrix. KWN models track a discretized size distribution of precipitates; nucleation, growth, and coarsening (Ostwald ripening) all emerge from the same population balance equation applied to that distribution. The notation here follows Ury et al., 2023, and the SFFK driving-force formulation follows Svoboda et al., 2004.

A complete KWN model is assembled by composing the precipitation primitives in this module with the Finite volume module, which provides the spatial discretization of the radius-space population balance equation. The KWN primitives compute the physics — volume fractions, matrix concentrations, growth rates, nucleation barriers and prefactors — and hand the resulting cell velocity and source term off to the generic finite-volume advection machinery.

Math¶

The size distribution of precipitate \((j)\) is discretized into finite-volume bins on the radius axis. Subscripts \(i\) index the radius bin (with center \(R^{(j)}_i\) and width \(\Delta R^{(j)}_i = R^{(j)}_{i+1/2} - R^{(j)}_{i-1/2}\)); parenthetical superscripts \((j)\) index the precipitate species. Matrix solute concentrations are stored as mole fractions \(x_k\) over the union \(\mathcal{S} = \bigcup_j \mathcal{S}_j\) of species that participate in any precipitate reaction.

Mass conservation¶

The volume fraction of precipitate \((j)\) is the third moment of its number-density distribution,

\[ f^{(j)} = \sum_i \tfrac{4}{3}\pi \left(R_i^{(j)}\right)^3 n_i^{(j)}, \]

and the current matrix concentration of species \(k\) follows from the closure that everything not locked up in precipitates remains in solution,

\[ x^\infty_k = \frac{x_{0,k} - \sum_j f^{(j)} x^{(j)}_k} {1 - \sum_j f^{(j)}}, \]

where \(x_{0,k}\) is the initial matrix concentration of species \(k\), \(x^{(j)}_k\) is the concentration of species \(k\) in precipitate \((j)\) (taken as zero for \(k \notin \mathcal{S}_j\)), and the initial precipitate volume fraction is assumed to be zero.

Population balance¶

Each precipitate population evolves according to a 1D advection equation in radius space with a nucleation source on the right,

\[ \frac{\partial n^{(j)}}{\partial t} + \frac{\partial \bigl(n^{(j)} \dot{R}^{(j)}\bigr)}{\partial R^{(j)}} = J^{(j)}_{\text{nuc}}. \]

The KWN module supplies \(\dot{R}^{(j)}\) at cell centers and \(J^{(j)}_{\text{nuc}}\); the Finite volume module interpolates the velocity to cell edges, computes the upwinded advective flux, applies boundary conditions, and takes the divergence to recover \(\partial n^{(j)}/\partial t\).

Growth rate¶

Two growth-rate forms are provided.

The simpler rate-limited form solves for a single rate-limiting species \(\star \in \mathcal{S}_j\),

\[ \dot{R}^{(j)} = \frac{1}{R^{(j)}} \, D_\star \, \frac{x^\infty_\star - x^{*,(j)}_\star} {x^{(j)}_\star - x^{*,(j)}_\star}, \]

with \(x^{*,(j)}_\star\) the matrix equilibrium concentration of the rate-limiting species in equilibrium with precipitate \((j)\). The growth rate vanishes as \(x^\infty_\star \rightarrow x^{*,(j)}_\star\).

The more general SFFK form (Svoboda–Fischer–Fratzl–Kozeschnik) handles multi-component diffusion in the dilute, diffusion-limited regime,

(6)¶\[\begin{align} \dot{R}^{(j)} &= \frac{1}{R^{(j)}} \, \frac{\Delta G^{(j)}_{\text{chem}} + \Delta G^{(j)}_{\text{surf}} + \Delta G^{(j)}_{\text{el}}} {R_g T \, S^{(j)}}, \\ S^{(j)} &= \sum_{k \in \mathcal{S}_j} \frac{\bigl(\Delta x^{(j)}_k\bigr)^2}{D_k \, x^\infty_k}, \end{align}\]

with \(\Delta x^{(j)}_k = x^{(j)}_k - x^{*,(j)}_k\) and \(D_k\) the matrix diffusivity of species \(k\). The driving force is split into chemical, surface, and elastic contributions:

(7)¶\[\begin{align} \Delta G^{(j)}_{\text{chem}} &= \sum_{k \in \mathcal{S}_j} \Delta x^{(j)}_k \bigl[\mu^{\text{matrix}}_k(x^\infty) - \mu^{\text{equil},(j)}_k\bigr], \\ \Delta G^{(j)}_{\text{surf}} &= -\frac{2 \gamma^{(j)} V_m^{(j)}}{R^{(j)}}, \\ \Delta G^{(j)}_{\text{el}} &= -3 V_m^{(j)} \varepsilon_0^{(j)} \sigma_h + \frac{6 V_m^{(j)} \mu K}{3K + 4\mu} \bigl(\varepsilon_0^{(j)}\bigr)^2. \end{align}\]

Here \(V_m^{(j)}\) is the molar volume of the precipitate, \(\gamma^{(j)}\) its surface energy, \(\varepsilon_0^{(j)}\) the misfit strain (coherent, isotropic spherical inclusion in an isotropic matrix with shear modulus \(\mu\) and bulk modulus \(K\)), and \(\sigma_h\) the matrix hydrostatic stress.

Nucleation flux¶

Nucleation deposits new precipitates at the critical radius,

\[ J^{(j)}_{\text{nuc}} = Z^{(j)} \beta^{(j)} N_0^{(j)} \exp\!\left(-\frac{\Delta G^{*,(j)}}{k_B T}\right) \delta\!\left(R^{(j)} - R^{(j)}_{\text{crit}}\right), \]

with the classical-nucleation-theory expressions

(8)¶\[\begin{align} R^{(j)}_{\text{crit}} &= \frac{2 \gamma^{(j)} V_m^{(j)}}{\Delta g^{(j)}_v}, \\ \Delta G^{*,(j)} &= \frac{16 \pi}{3} \, \frac{\bigl(\gamma^{(j)}\bigr)^3 \bigl(V_m^{(j)}\bigr)^2} {\bigl(\Delta g^{(j)}_v\bigr)^2}, \\ Z^{(j)} &= \frac{V_m^{(j)}}{2\pi N_a \bigl(R^{(j)}_{\text{crit}}\bigr)^2} \sqrt{\frac{\gamma^{(j)}}{k_B T}}, \\ \beta^{(j)} &= \frac{4 \pi \bigl(R^{(j)}_{\text{crit}}\bigr)^2 N_a^{4/3}} {\bigl(V_m^{(j)}\bigr)^{4/3}} \, \frac{1}{S^{(j)}}. \end{align}\]

\(\Delta g^{(j)}_v\) is the volumetric driving force for nucleation of precipitate \((j)\) (in J/mol — the per-mole convention used throughout NucleationBarrierAndCriticalRadius), \(N_0^{(j)}\) is the nucleation site density, \(N_a\) is Avogadro’s number, \(k_B\) is the Boltzmann constant, and \(S^{(j)}\) is the same projected diffusivity sum that appears in the SFFK growth rate.

The volumetric driving force \(\Delta g^{(j)}_v\) can either be supplied as a CALPHAD-derived tabulation against composition and temperature (the pattern used by the Al–Cu example below) or assembled from the participating species’ matrix and equilibrium concentrations with IdealSolutionVolumetricDrivingForce, which evaluates the Hu–Cocks ideal-solution form \(\Delta g_v = R_g T \sum_k w_k \ln(c_k / c_k^{\text{eq}})\). Setting every weight to 1 recovers the “product of all components” convention used for compound precipitates.

Note

For most systems \(R_{\text{crit}}\) is much smaller than the smallest cell size used in a practical finite-volume discretization. The Dirac delta at \(R^{(j)}_{\text{crit}}\) is therefore approximated discretely by DumpInSmallestBin from the Finite volume module, which deposits the entire flux magnitude into the smallest bin.

Catalog¶

The KWN module exposes the following primitives, each implementing one of the equations above. They are designed to be composed inside a ComposedModel (typically driven by an ImplicitUpdate) rather than used in isolation.

Type	Role
PrecipitateVolumeFraction	Sums \(\tfrac{4}{3}\pi (R_i^{(j)})^3 n_i^{(j)}\) over bins to obtain \(f^{(j)}\)
CurrentConcentration	Mass-balance closure for \(x^\infty_k\) given the precipitate populations
ProjectedDiffusivitySum	Computes the multi-species sum \(S^{(j)}\) shared by SFFK and nucleation
ChemicalGibbsFreeEnergyDifference	Assembles \(\Delta G^{(j)}_{\text{chem}}\) from per-species potentials
IdealSolutionVolumetricDrivingForce	Assembles the molar driving force \(R_g T \sum_k w_k \ln(c_k / c_k^{\text{eq}})\)
RateLimitedPrecipitateGrowthRate	Single-species, equilibrium-driven growth rate \(\dot{R}^{(j)}\)
SFFKPrecipitationGrowthRate	SFFK multi-component growth rate \(\dot{R}^{(j)}\)
NucleationBarrierAndCriticalRadius	Computes \(R^{(j)}_{\text{crit}}\) and \(\Delta G^{*,(j)}\)
ZeldovichFactor	Computes the Zeldovich factor \(Z^{(j)}\)
KineticFactor	Computes the attachment frequency \(\beta^{(j)}\)
NucleationFluxMagnitude	Assembles \(Z^{(j)} \beta^{(j)} N_0^{(j)} \exp(-\Delta G^{*,(j)} / k_B T)\)

Example: growth-only Al–Cu¶

The input below is the growth-only Al–Cu regression case tests/regression/kwn/growth-only-scaled/model.i. A uniform grid in a scaled radius coordinate \(\xi \in [0, 1]\) is mapped to a semi-infinite physical radius \(R = s\,\xi / (1 - \xi)\) via the precomputed true_centers, center_jacobian, and center_inverse_jacobian tensors. The matrix concentration and growth rate are recomputed from the current size distribution each step, and the population balance is advanced in time by an ImplicitUpdate Newton solve.

# neml2
# Units: length = microns, time = hours

# Initial setup: Al with 4 wt% Cu
# That gives 0.01738 mole fraction of Cu
x0_Cu = 0.01738

# Let's assume a temperature of 300 K, I can run a pycalphad simulation to predict the equilbrium Cu concentration in the matrix, precipitate, and the difference
T = 300.0
xp_Cu = 3.333291e-01
diff_FCC_Cu = 3.333270e-01

# Diffusivity ~ 0.150 exp[-30200/(RT)] cm^2/s which gives 297794 microns^2/hour
D = 297794

[Tensors]
  [edges]
    type = Python
    expr = 'Scalar.linspace(0.0, 1.0, 101).sub_batch.retag(1)'
  []
  [centers]
    type = Python
    expr = 'Scalar(0.5 * (edges.data[..., :-1] + edges.data[..., 1:]), sub_batch_ndim=1)'
  []
  [dx_centers]
    type = Python
    expr = 'Scalar(centers.data[..., 1:] - centers.data[..., :-1], sub_batch_ndim=1)'
  []
  [dx]
    type = Python
    expr = 'Scalar(edges.data[..., 1:] - edges.data[..., :-1], sub_batch_ndim=1)'
  []

  [scale_factor]
    type = Python
    expr = 'Scalar(5.0)'
  []

  [true_centers]
    type = Python
    expr = 'Scalar(scale_factor.data * centers.data / (1.0 - centers.data), sub_batch_ndim=1)'
  []

  [center_jacobian]
    type = Python
    expr = 'Scalar(scale_factor.data / (1.0 - centers.data) ** 2, sub_batch_ndim=1)'
  []
  [center_inverse_jacobian]
    type = Python
    expr = 'Scalar((1.0 - centers.data) ** 2 / scale_factor.data, sub_batch_ndim=1)'
  []

  [unscaled_ic]
    type = Python
    expr = 'Scalar(1e-12 * torch.ones(100, dtype=torch.float64), sub_batch_ndim=1)'
  []

  [ic]
    type = Python
    expr = 'Scalar(unscaled_ic.data * center_jacobian.data, sub_batch_ndim=1)'
  []

  [time]
    type = Python
    expr = 'Scalar.linspace(0.0, 100.0, 500)'
  []

  [x0_Cu]
    type = Python
    expr = 'Scalar(torch.tensor(${x0_Cu}, dtype=torch.float64))'
  []
  [xp_Cu]
    type = Python
    expr = 'Scalar(torch.tensor(${xp_Cu}, dtype=torch.float64))'
  []

  [X_Cu_vary]
    type = Python
    expr = 'Scalar(torch.tensor([1e-10, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05, 0.055, 0.06, 0.065, 0.07, 0.075, 0.08, 0.085, 0.09, 0.095], dtype=torch.float64), sub_batch_ndim=1)'
  []
  [chem_diff]
    type = Python
    expr = 'Scalar(torch.tensor([-8296.144765949808, 4829.582140453559, 4829.582140453546, 4829.582140453556, 4829.582140453537, 4829.582140453549, 4829.5821404535445, 4829.5821404535445, 4829.582140453554, 4829.582140453549, 4829.582140453542, 4829.582140453542, 4829.582140453532, 4829.582140453542, 4829.582140453554, 4829.5821404535445, 4829.582140453552, 4829.582140453537, 4829.5821404535445, 4829.582140453542], dtype=torch.float64), sub_batch_ndim=1)'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'time'
    ic_Scalar_names = 'number_density'
    ic_Scalar_values = 'ic'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'implicit_rate'
    unknowns = 'number_density'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [input_temperature]
    type = ScalarConstantParameter
    value = ${T}
    parameter = 'forces/T'
  []
  [volume_fraction]
    type = PrecipitateVolumeFraction
    radius = 'true_centers'
    number_density = 'true_number_density'
    volume_fraction = 'vf'
  []
  [x_Cu]
    type = CurrentConcentration
    initial_concentration = 'x0_Cu'
    precipitate_volume_fractions = 'vf'
    precipitate_concentrations = 'xp_Cu'
    current_concentration = 'x_Cu'
  []

  [chemical_potential_difference]
    type = ScalarLinearInterpolation
    argument = 'x_Cu'
    abscissa = 'X_Cu_vary'
    ordinate = 'chem_diff'
  []

  [diffusivity_sum]
    type = ProjectedDiffusivitySum
    concentration_differences = ${diff_FCC_Cu}
    diffusivities = ${D}
    far_field_concentrations = 'x_Cu'
    projected_diffusivity_sum = 'diff_sum'
  []
  [growth_rate]
    type = SFFKPrecipitationGrowthRate
    radius = 'true_centers'
    projected_diffusivity_sum = 'diff_sum'
    gibbs_free_energy_difference = chemical_potential_difference
    temperature = 'forces/T'
    gas_constant = 8.314
    growth_rate = 'growth_rate'
  []

  [scaled_cell_velocity]
    type = ScalarMultiplication
    from = 'growth_rate'
    scaling = 'center_inverse_jacobian'
    to = 'internal/scaled_cell_velocity'
  []
  [advection_velocity]
    type = LinearlyInterpolateToCellEdges
    cell_values = 'internal/scaled_cell_velocity'
    cell_centers = 'centers'
    cell_edges = 'edges'
    edge_values = 'v_edge'
  []
  [advective_flux]
    type = FiniteVolumeUpwindedAdvectiveFlux
    u = 'number_density'
    v_edge = 'v_edge'
    flux = 'J'
  []
  [left_bc]
    type = FiniteVolumeAppendBoundaryCondition
    input = 'J'
    bc_value = 0.0
    side = 'left'
  []
  [right_bc]
    type = FiniteVolumeAppendBoundaryCondition
    input = 'J_with_bc_left'
    bc_value = 0.0
    side = 'right'
  []
  [flux_divergence]
    type = FiniteVolumeGradient
    u = 'J_with_bc_left_with_bc_right'
    dx = 'dx'
    grad_u = 'number_density_rate'
  []
  [integrate_u]
    type = ScalarBackwardEulerTimeIntegration
    variable = 'number_density'
  []
  [implicit_rate]
    type = ComposedModel
    models = 'input_temperature growth_rate scaled_cell_velocity advection_velocity advective_flux left_bc right_bc flux_divergence integrate_u unscale volume_fraction x_Cu diffusivity_sum'
  []
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_Scalar = 'number_density'
  []
  [model_scaled]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [unscale]
    type = ScalarMultiplication
    from = 'number_density'
    scaling = 'center_inverse_jacobian'
    to = 'true_number_density'
  []
  [model]
    type = ComposedModel
    models = 'model_scaled unscale volume_fraction x_Cu'
    additional_outputs = 'number_density true_number_density vf x_Cu'
  []
[]

Walkthrough¶

PrecipitateVolumeFraction (volume_fraction) integrates the physical number density true_number_density over the physical bin centers true_centers to give the precipitate volume fraction vf. The unscale step below converts the solution variable number_density (which lives on the scaled grid) into true_number_density by multiplying by the inverse Jacobian.
CurrentConcentration (x_Cu) closes the matrix composition with \(x^\infty_{\text{Cu}} = (x_{0,\text{Cu}} - f\, x^{(j)}_{\text{Cu}}) / (1 - f)\), given the initial Cu mole fraction x0_Cu and the precipitate Cu fraction xp_Cu.
A ScalarLinearInterpolation (chemical_potential_difference) evaluates the tabulated chemical-potential difference \(\mu^{\text{matrix}}_{\text{Cu}}(x^\infty_{\text{Cu}}) - \mu^{\text{equil}}_{\text{Cu}}\) at the current matrix concentration. This is the externally supplied \(\Delta G^{(j)}_{\text{chem}}\) contribution; in a CALPHAD-coupled workflow it would come from a pycalphad tabulation against composition and temperature.
ProjectedDiffusivitySum (diffusivity_sum) computes \(S^{(j)} = (\Delta x^{(j)}_{\text{Cu}})^2 / (D_{\text{Cu}} x^\infty_{\text{Cu}})\) with a single contributing species, since the only diffusing solute in this example is Cu.
SFFKPrecipitationGrowthRate (growth_rate) divides the chemical driving force by \(R_g T S^{(j)} R^{(j)}\) to give \(\dot{R}^{(j)}\) at each cell center. (Surface and elastic terms are omitted in this single-species, growth-only example; both are optional inputs on the SFFK model.)
scaled_cell_velocity, advection_velocity, advective_flux, left_bc, right_bc, and flux_divergence are generic finite-volume primitives from the Finite volume module. They map the cell-center growth rate to the scaled grid, interpolate it to cell edges, take the upwinded advective flux, append zero-flux Dirichlet boundary conditions at both ends, and take the divergence to recover \(\partial n^{(j)}/\partial t\).
integrate_u adds the backward-Euler residual; implicit_rate bundles every model in the residual graph; and model_scaled drives the Newton solve through ImplicitUpdate. The outer model then unscales the solution and exposes the diagnostic outputs (number_density, true_number_density, vf, x_Cu) for post-processing.

To add nucleation, the same skeleton is augmented with NucleationBarrierAndCriticalRadius, ZeldovichFactor, KineticFactor, and NucleationFluxMagnitude; the source term is then dumped into the smallest bin and added to the flux divergence via a ScalarLinearCombination. The full composition lives in tests/regression/kwn/growth-nucleation-scaled/model.i.

Examples¶

End-to-end examples live in tests/regression/kwn/:

growth-only-scaled/ — single-species growth with a tabulated chemical potential, semi-infinite scaled radius grid.
growth-nucleation-scaled/ — full nucleation + growth on a semi-infinite scaled radius grid.
growth-nucleation-unscaled/ — same physics on a fixed-extent radius grid (useful for isolating the effect of the semi-infinite scaling).