Governing Equations

AdditiveFOAM solves incompressible momentum conservation with Boussinesq buoyancy and mushy-zone drag, temperature transport with latent heat, and a tabulated temperature–solid-fraction relation. The formulation includes semi-implicit thermodynamic coupling, phase-dependent thermal properties, and irreversible powder consolidation.

Continuum model

AdditiveFOAM solves incompressible mass and momentum conservation for the mixture velocity \(\mathbf{u}\). In the kinematic form used by the solver,

\[\nabla\cdot\mathbf{u}=0,\] \[\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u} =\nabla\cdot(\nu\nabla\mathbf{u})-\nabla p+\rho_k\mathbf{g}-D_\nu\mathbf{u}.\]

The Boussinesq density factor is

\[\rho_k=1-\beta(T-T_{\mathrm{liq}}),\]

and the mushy-zone resistance is based on a capped Kozeny–Carman form,

\[\widetilde f_s=\min(f_s,0.95),\] \[D_\nu=\nu\frac{180}{\lambda^2} \frac{\widetilde f_s^2}{(1-\widetilde f_s)^3}.\]

Here \(\mathbf{u}\) is mixture velocity, \(t\) is time, \(\nu=\mu/\rho\) is kinematic viscosity, \(\mu\) is dynamic viscosity, \(\rho\) is reference density, \(p\) is kinematic pressure, \(\mathbf{g}\) is gravitational acceleration, and \(\nabla\) is the spatial gradient operator. The temperature is \(T\), \(T_{\mathrm{liq}}\) is the liquidus, \(\beta\) is the thermal-expansion coefficient, and \(\rho_k\) is the dimensionless Boussinesq density factor. The solid fraction is \(f_s\), \(D_\nu\) is the kinematic mushy-zone drag coefficient, and \(\lambda\) is the dendrite-arm spacing. The 0.95 cap keeps the resistance finite in fully solid cells. When a rotating frame is configured, OpenFOAM’s MRF acceleration is added to the left-hand side.

Temperature and phase coupling

The energy equation is written in temperature form:

\[\rho c_p\frac{\partial T}{\partial t} +\rho c_p\nabla\cdot(\mathbf{u}T) =\nabla\cdot(k\nabla T)+\rho L_f\frac{\partial f_s}{\partial t}+\dot{q}.\]

Here \(c_p\) is specific heat capacity, \(k\) is thermal conductivity, \(L_f\) is latent heat of fusion, and \(\dot q\) is volumetric heat input. The thermoPath table defines the solid-fraction function \(f_s=F(T)\) and can represent Gulliver–Scheil, lever-rule, or linear behavior. AdditiveFOAM interpolates piecewise linearly between the table points and applies \(\operatorname{clip}_{[0,1]}\!\left[x\right]=\min[\max(x,0),1]\).

Semi-implicit correction from the tabulated path

Let \(m\) denote the thermodynamic-correction index at the new time level. For a temperature inside the solidification interval, the solver brackets \(T^{(m)}\) by lower and upper table points \((T_l,F_l)\) and \((T_h,F_h)\) and evaluates the local path slope

\[g^{(m)}=\left.\frac{\mathrm{d}F}{\mathrm{d}T}\right|_m =\frac{F_h-F_l}{T_h-T_l}.\]

The solver maps \(f_s^{(m)}\) back onto the local table segment:

\[T_0^{(m)}=T_l+\frac{f_s^{(m)}-F_l}{g^{(m)}}.\]

This places the iterated phase state on the local linear model. The next solid-fraction estimate is

\[f_s^{(m+1)}= \operatorname{clip}_{[0,1]}\! \left[f_s^{(m)}+g^{(m)}\left(T^{(m+1)}-T_0^{(m)}\right)\right].\]

For a linear multistep time derivative, only the new-time value changes during the nonlinear correction. If \(a_0\) is the new-time derivative coefficient and \(\Delta t\) is the current time-step duration, then

\[\left.\frac{\partial f_s}{\partial t}\right|^{(m+1)} \approx \left.\frac{\partial f_s}{\partial t}\right|^{(m)} +\frac{a_0}{\Delta t}g^{(m)} \left(T^{(m+1)}-T_0^{(m)}\right).\]

Substitution into the energy equation gives the temperature solve

\[\mathcal{L}_T\!\left(T^{(m+1)}\right) =\rho L_f\left[ \left.\frac{\partial f_s}{\partial t}\right|^{(m)} +\frac{a_0}{\Delta t}g^{(m)} \left(T^{(m+1)}-T_0^{(m)}\right) \right]+\dot q,\]

where \(\mathcal{L}_T\) contains sensible-energy storage, advection, and diffusion. Moving the term proportional to \(T^{(m+1)}\) to the matrix gives

\[\mathcal{L}_T\!\left(T^{(m+1)}\right) -\rho L_f\frac{a_0}{\Delta t}g^{(m)}T^{(m+1)} =\rho L_f\left[ \left.\frac{\partial f_s}{\partial t}\right|^{(m)} -\frac{a_0}{\Delta t}g^{(m)}T_0^{(m)} \right]+\dot q.\]

Within the usual decreasing solid-fraction path, \(g^{(m)}\leq0\), so the latent-heat Jacobian adds a non-negative contribution to the temperature-equation diagonal. This is the semi-implicit stabilization: the tabulated path remains explicit data, while its local slope enters the linear system implicitly.

The new temperature is used in the clipped phase update above. Iteration stops after at least two corrections when

\[\max_{\text{cells}}\left|f_s^{(m+1)}-f_s^{(m)}\right| <\texttt{thermoTolerance},\]

or when nThermoCorrectors is reached. Outside the tabulated solidification interval, limiting branches drive inconsistent phase states back toward the liquidus or solidus.

Time-integration methods

The new-time coefficient used in the correction is

Scheme	\(a_0\)	Behavior
Euler	\(1\)	First-order implicit Euler
backward	\(1+\dfrac{\Delta t}{\Delta t+\Delta t_0}\)	Second-order backward differencing with variable time steps; \(a_0=3/2\) for equal steps
CrankNicolson	\(1+\psi\)	Crank–Nicolson with off-centering coefficient \(\psi\); \(\psi=1\) is fully centred

Here \(\Delta t_0\) is the preceding time-step duration and \(\psi\) is the OpenFOAM Crank–Nicolson off-centering coefficient.

Select the scheme in system/fvSchemes:

ddtSchemes
{
    default backward;
}

and configure the nonlinear correction in system/fvSolution:

PIMPLE
{
    nThermoCorrectors 20;
    thermoTolerance   1e-8;
    explicitSolve     false;
}

explicitSolve true requests explicit sensible-energy transport when its diffusion-number check passes and forces forward Euler regardless of fvSchemes. Use false to guarantee implicit Euler, backward, or CrankNicolson transport.

Phase properties and powder transition

Let \(f_p\) be the alpha.powder field and \(f_s\) the solid fraction obtained from thermoPath. The liquid fraction is

\[f_l=1-f_s.\]

For either thermal conductivity or heat capacity, denoted generically by \(\psi\in\{k,c_p\}\), AdditiveFOAM evaluates the mixture property as

\[\psi_{\mathrm{eff}} =(1-f_p)\left[f_s\psi_s(T_s)+f_l\psi_l(T_l)\right] +f_p\psi_p(T_s),\]

where \(\psi_s\), \(\psi_l\), and \(\psi_p\) are the configured solid, liquid, and powder property polynomials. Their evaluation temperatures are limited to

\[T_s=\operatorname{clip}_{[300\ \mathrm{K},\,T_{\mathrm{solidus}}]}\!\left[T\right],\] \[T_l=\operatorname{clip}_{[T_{\mathrm{liquidus}},\,2T_{\mathrm{liquidus}}]}\!\left[T\right].\]

Here \(f_p\), \(f_s\), and \(f_l\) are the powder, solid, and liquid fractions; \(T_s\) and \(T_l\) are the temperatures used to evaluate the solid/powder and liquid polynomials, respectively. The subscript on \(\psi_i\) identifies phase \(i\in\{s,l,p\}\).

Powder consolidation is represented by the irreversible update

\[f_p^{\,n+1}= \begin{cases} 0, & f_s^{\,n+1}<1,\\ f_p^{\,n}, & f_s^{\,n+1}=1. \end{cases}\]

Here \(n\) and \(n+1\) denote the old and new time levels.

A cell retains its initialized powder fraction while fully solid. Melting sets the powder fraction to zero permanently, so subsequent solidification uses the consolidated solid–liquid mixture.