|
MFC
Exascale flow solver
|
|
MFC
Exascale flow solver
|
This document catalogs every equation solved by MFC, organized by physical model. Each section notes the input parameter(s) that activate the corresponding physics module and cross-references the relevant source files.
The models and algorithms described here are detailed in Wilfong et al. [54] (MFC 5.0) and Bryngelson et al. [8]. Foundational references for each model are cited inline; see the Bibliography for full details.
For parameter details and allowed values, see Case Files and the Case Parameters reference.
MFC solves the compressible Navier-Stokes equations (or Euler equations when viscosity is off) in a finite volume framework. The general semi-discrete form is:
\[\frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u} = \mathbf{s}(\mathbf{q})\]
where:
The parameter model_eqns (1, 2, 3, or 4) selects the governing equation set.
Key source files: src/simulation/m_rhs.fpp (RHS evaluation), src/common/m_variables_conversion.fpp (EOS and variable conversion).
The main flow solver (Navier-Stokes equations, Riemann solvers, viscous stress, body forces, surface tension, etc.) is unit-agnostic: whatever units the user provides for the initial and boundary conditions, the solver preserves them throughout the computation. If the user inputs SI units, the outputs are in SI units. If the user inputs CGS, the outputs are in CGS. No internal non-dimensionalization is performed by the flow solver.
This means that for simulations without sub-grid bubble models, the user can work in any consistent unit system without additional effort.
Several EOS and transport parameters use transformed stored forms that differ from the standard physical values. This is the most common source of input errors:
| Parameter | Physical quantity | What MFC expects (stored form) |
|---|---|---|
| fluid_pp(i)%gamma | Heat capacity ratio \(\gamma\) | \(\Gamma = \frac{1}{\gamma - 1}\) |
| fluid_pp(i)%pi_inf | Stiffness pressure \(\pi_\infty\) [Pa] | \(\Pi_\infty = \frac{\gamma\,\pi_\infty}{\gamma - 1}\) [Pa] |
| fluid_pp(i)%Re(1) | Dynamic viscosity \(\mu\) | \(1/\mu\) (inverse viscosity) |
| fluid_pp(i)%Re(2) | Bulk viscosity \(\mu_b\) | \(1/\mu_b\) (inverse bulk viscosity) |
These transformations arise because MFC internally solves the energy equation using the transformed variables \(\Gamma\) and \(\Pi_\infty\) (see Section 3.1), and the viscous stress is computed by dividing by Re rather than multiplying by \(\mu\).
Common mistake: setting fluid_pp(1)%gamma = 1.4 for air. The correct value is 1.0 / (1.4 - 1.0) = 2.5. Setting gamma = 1.4 corresponds to a physical \(\gamma \approx 1.71\), which is not a standard gas.
Pre-computed stored-form values for common fluids (SI units):
| Material | \(\gamma\) | \(\pi_\infty\) [Pa] | gamma (stored) | pi_inf (stored) [Pa] |
|---|---|---|---|---|
| Air | 1.4 | 0 | 2.5 | 0 |
| Helium | 5/3 | 0 | 1.5 | 0 |
| Water (Tait) | 4.4 | 6.0e8 | 0.2941 | 7.76e8 |
| Water (Le Metayer et al. [27]) | 6.12 | 3.43e8 | 0.1953 | 4.10e8 |
Example for an air-water simulation:
For viscous cases, provide the reciprocal of the dynamic viscosity:
The solver does not check or convert units. All inputs must use the same consistent unit system (e.g., all SI or all CGS). Mixing units — for example, pressures in atmospheres with densities in kg/m³ — will produce silently incorrect results.
The sub-grid bubble models (bubbles_euler = .true. or bubbles_lagrange = .true.) solve the bubble wall dynamics in non-dimensional form. The bubble wall pressure equation as implemented is:
\[p_{bw} = \left(\text{Ca} + \frac{2}{\text{We}_b\,R_0}\right)\left(\frac{R_0}{R}\right)^{3\gamma} - \text{Ca} - \frac{4\,\text{Re}_{\text{inv}}\,\dot{R}}{R} - \frac{2}{R\,\text{We}_b}\]
Here \(R\) and \(R_0\) are non-dimensional radii (scaled by \(x_0\)), and \(\dot{R}\) is a non-dimensional wall speed (scaled by \(u_0\)); the entire bubble ODE is solved in non-dimensional variables.
The dimensionless groups are:
| Dimensionless group | Definition | Code variable | Computed from |
|---|---|---|---|
| \(\text{Ca}\) (Cavitation number) | \(p_{0,\text{ref}} - p_v\) | Ca | bub_pp%p0ref - bub_pp%pv |
| \(\text{Eu}\) (Euler number) | \(p_{0,\text{ref}}\) | Eu | bub_pp%p0ref |
| \(\text{We}_b\) (bubble Weber number) | \(1/\sigma\) | Web | 1 / bub_pp%ss |
| \(\text{Re}_{\text{inv}}\) (inverse bubble Reynolds number) | \(\mu_l\) | Re_inv | bub_pp%mu_l |
Because the bubble equations use these dimensionless numbers directly, all bub_pp%% inputs are interpreted by the code as already non-dimensional. The code does not non-dimensionalize bubble quantities internally. Therefore, when bubbles are enabled, the simulation must be run in a fully non-dimensional form: all inputs — flow ICs/BCs, EOS parameters, domain lengths, dt, and bub_pp%% values — must be scaled with the same \((x_0, p_0, \rho_0, u_0, t_0, T_0)\) reference quantities, or the coupled solution will be physically incorrect.
When using bubble models, the user must choose reference scales and non-dimensionalize all inputs (flow and bubble) consistently. The standard convention used in the MFC examples is:
| Reference quantity | Symbol | Typical choice |
|---|---|---|
| Length | \(x_0\) | \(R_{0,\text{ref}}\) (reference bubble radius) |
| Pressure | \(p_0\) | \(p_{0,\text{ref}}\) (reference bubble pressure) |
| Density | \(\rho_0\) | \(\rho_{0,\text{ref}}\) (reference liquid density) |
| Velocity | \(u_0\) | \(\sqrt{p_0 / \rho_0}\) (derived) |
| Time | \(t_0\) | \(x_0 / u_0\) (derived) |
| Temperature | \(T_0\) | \(T_{0,\text{ref}}\) (reference temperature) |
The following table lists every bub_pp%% parameter and its required non-dimensionalization:
| Parameter | Physical meaning | Non-dimensional form |
|---|---|---|
| bub_pp%R0ref | Reference bubble radius | \(R_{0,\text{ref}} / x_0\) |
| bub_pp%p0ref | Reference bubble pressure | \(p_{0,\text{ref}} / p_0\) |
| bub_pp%rho0ref | Reference liquid density | \(\rho_{0,\text{ref}} / \rho_0\) |
| bub_pp%T0ref | Reference temperature | \(T_{0,\text{ref}} / T_0\) (typically 1) |
| bub_pp%ss | Surface tension \(\sigma\) | \(\sigma / (\rho_0\,x_0\,u_0^2)\) |
| bub_pp%pv | Vapor pressure | \(p_v / p_0\) |
| bub_pp%mu_l | Liquid dynamic viscosity | \(\mu_l / (\rho_0\,x_0\,u_0)\) |
| bub_pp%mu_v | Vapor dynamic viscosity | \(\mu_v / (\rho_0\,x_0\,u_0)\) |
| bub_pp%mu_g | Gas dynamic viscosity | \(\mu_g / (\rho_0\,x_0\,u_0)\) |
| bub_pp%vd | Vapor diffusivity | \(D / (x_0\,u_0)\) |
| bub_pp%k_v | Vapor thermal conductivity | \(k_v\,T_0 / (x_0\,\rho_0\,u_0^3)\) |
| bub_pp%k_g | Gas thermal conductivity | \(k_g\,T_0 / (x_0\,\rho_0\,u_0^3)\) |
| bub_pp%cp_v | Vapor specific heat | \(c_{p,v}\,T_0 / u_0^2\) |
| bub_pp%cp_g | Gas specific heat | \(c_{p,g}\,T_0 / u_0^2\) |
| bub_pp%R_v | Vapor gas constant | \(R_v\,T_0 / u_0^2\) |
| bub_pp%R_g | Gas gas constant | \(R_g\,T_0 / u_0^2\) |
| bub_pp%gam_v | Vapor heat capacity ratio | Already dimensionless (no scaling) |
| bub_pp%gam_g | Gas heat capacity ratio | Already dimensionless (no scaling) |
| bub_pp%M_v | Vapor molar mass | Consistent units; only ratios are used (no scaling needed) |
| bub_pp%M_g | Gas molar mass | Consistent units; only ratios are used (no scaling needed) |
When the reference scales match the bubble reference values (e.g., \(x_0 = R_{0,\text{ref}}\), \(p_0 = p_{0,\text{ref}}\), \(\rho_0 = \rho_{0,\text{ref}}\)), the reference parameters simplify to unity: bub_pp%R0ref = 1, bub_pp%p0ref = 1, bub_pp%rho0ref = 1.
When bubbles are enabled, the flow-level parameters must also be non-dimensionalized with the same reference scales:
| Parameter | Non-dimensional form |
|---|---|
| x_domain%beg, x_domain%end | Domain bounds divided by \(x_0\) |
| patch_icpp(i)%pres | Pressure divided by \(p_0\) |
| patch_icpp(i)%alpha_rho(j) | Partial density divided by \(\rho_0\) |
| patch_icpp(i)%vel(j) | Velocity divided by \(u_0\) |
| fluid_pp(i)%gamma | \(1/(\gamma_i - 1)\) (dimensionless, same as without bubbles) |
| fluid_pp(i)%pi_inf | \(\gamma_i\,\pi_{\infty,i} / [(\gamma_i - 1)\,p_0]\) (scaled by reference pressure) |
| fluid_pp(i)%Re(1) | \(\rho_0\,x_0\,u_0 / \mu_i\) (Reynolds number, inverse viscosity) |
| dt | Time step divided by \(t_0\) |
MFC has two conceptually distinct viscosity-related parameters that serve different physical roles:
\[\tau_{ij} \propto \frac{\nabla u}{\text{Re}}\]
Stored in the physical_parameters derived type (src/common/m_derived_types.fpp).\[p_{bw} \ni -\frac{4\,\text{Re}_{\text{inv}}\,\dot{R}}{R}\]
Stored in the subgrid_bubble_physical_parameters derived type (src/common/m_derived_types.fpp).These two parameters represent viscous effects at fundamentally different scales — bulk flow dissipation vs. single-bubble-wall damping — and are stored in separate derived types with separate code paths. They are not interchangeable: fluid_pp%Re(1) is an inverse viscosity while bub_pp%mu_l is a viscosity (non-dimensionalized).
A typical bubble case setup in case.py follows this pattern:
Note the inverse relationship: fluid_pp%Re(1) = 1 / bub_pp%mu_l when both use the same reference scales and the same physical viscosity. This is expected — they encode the same physical viscosity but in reciprocal forms for their respective equations.
The primary workhorse model (Allaire et al. [1]; Wilfong et al. [54] Sec. 2.1). The state vector is:
\[\mathbf{q} = \bigl(\alpha_1 \rho_1,\;\alpha_2 \rho_2,\;\ldots,\;\rho u_1,\;\rho u_2,\;\rho u_3,\;\rho E,\;\alpha_1,\;\alpha_2,\;\ldots\bigr)^T\]
Continuity (one per component):
\[\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0\]
Momentum:
\[\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot \bigl(\rho\,\mathbf{u} \otimes \mathbf{u} + p\,\mathbf{I} - \boldsymbol{\tau}^v\bigr) = 0\]
Energy:
\[\frac{\partial (\rho E)}{\partial t} + \nabla \cdot \bigl[(\rho E + p)\,\mathbf{u} - \boldsymbol{\tau}^v \cdot \mathbf{u}\bigr] = 0\]
Volume fraction advection:
\[\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}\]
where the \(K\) term enforces interface conditions via the Wood sound speed:
\[K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\displaystyle\frac{\rho_1 c_1^2}{\alpha_1} + \displaystyle\frac{\rho_2 c_2^2}{\alpha_2}}\]
Setting alt_soundspeed = .true. enables the \(K\) correction (Kapila et al. [25], with Wood sound speed). Setting alt_soundspeed = .false. uses the Allaire variant without the \(K\) correction, which is conservative but does not strictly obey the second law of thermodynamics.
Mixture rules:
\[1 = \sum_i \alpha_i, \qquad \rho = \sum_i \alpha_i \rho_i, \qquad \rho e = \sum_i \alpha_i \rho_i e_i\]
Allows pressure disequilibrium between phases (Saurel et al. [41]; Wilfong et al. [54] Sec. 2.1).
Continuity and momentum: Same as the five-equation model.
Separate phasic internal energy:
\[\frac{\partial (\alpha_i \rho_i e_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i e_i\,\mathbf{u}) + \alpha_i p_i\,\nabla \cdot \mathbf{u} = -\mu\,p_I\,(p_2 - p_1) - \alpha_i\,\boldsymbol{\tau}_i^v : \nabla \mathbf{u}\]
Volume fraction:
\[\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)\]
Interfacial pressure:
\[p_I = \frac{z_2\,p_1 + z_1\,p_2}{z_1 + z_2}, \qquad z_i = \rho_i\,c_i\]
Infinite pressure relaxation is applied at each Runge-Kutta stage to drive toward pressure equilibrium.
Mixture speed of sound:
\[c^2 = \sum_k Y_k\,c_k^2\]
With phase change (relax = .true.), additional source terms appear in the phasic energy and volume fraction equations:
See Section 8 (Phase Change) below for details.
The primary closure for each phase:
\[p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}\]
Equivalently:
\[e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}\]
Total energy relation:
\[\rho E = \Gamma\,p + \Pi_\infty + \frac{1}{2}\rho\,|\mathbf{u}|^2 + q_v\]
where MFC internally tracks the transformed thermodynamic quantities:
\[\Gamma_k = \frac{1}{\gamma_k - 1}, \qquad \Pi_{\infty,k} = \frac{\gamma_k\,\pi_{\infty,k}}{\gamma_k - 1}\]
and the mixture rules are arithmetic averages of these transformed quantities:
\[\Gamma = \sum_i \frac{\alpha_i}{\gamma_i - 1}, \qquad \Pi_\infty = \sum_i \frac{\alpha_i\,\gamma_i\,\pi_{\infty,i}}{\gamma_i - 1}, \qquad q_v = \sum_i \alpha_i\,\rho_i\,q_{v,i}\]
The pressure is recovered from the total energy as:
\[p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}\]
Phasic speed of sound:
\[c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}\]
Wood mixture sound speed:
\[\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}\]
Input parameters per fluid: gamma ( \(\Gamma_k = 1/(\gamma_k - 1)\)), pi_inf ( \(\Pi_{\infty,k} = \gamma_k\,\pi_{\infty,k}/(\gamma_k - 1)\)), cv ( \(c_{v,k}\)), qv ( \(q_{v,k}\)), qvp ( \(q'_{v,k}\)). Note that gamma and pi_inf are stored in transformed form, not as the raw physical values (see Section 1b).
For reacting gas mixtures:
\[p = \frac{\rho\,R_u\,T}{W}, \qquad W = \left(\sum_m \frac{Y_m}{W_m}\right)^{-1}\]
Temperature is obtained from the internal energy by Newton iteration:
\[e_g - \sum_m e_m(T)\,Y_m = 0\]
Species internal energy from enthalpy:
\[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\]
NASA polynomial enthalpies (McBride et al. [31]):
\[\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}\]
Newtonian viscous stress (no bulk viscosity by default):
\[\boldsymbol{\tau}^v = 2\,\eta\left(\mathbf{D} - \frac{1}{3}\,\text{tr}(\mathbf{D})\,\mathbf{I}\right)\]
where the strain rate tensor is:
\[\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)\]
With bulk viscosity:
\[\tau_{ij} = \mu\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) + \left(\zeta - \frac{2\mu}{3}\right)\delta_{ij}\,\frac{\partial u_k}{\partial x_k}\]
Cartesian components:
\[\tau_{xx} = \mu\left(2\,\frac{\partial u}{\partial x} - \frac{2}{3}\nabla\cdot\mathbf{u}\right), \qquad \tau_{xy} = \mu\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\]
and similarly for all other components. Cylindrical coordinate formulations include additional \(1/r\) terms.
Viscosity averaging:
\[\frac{1}{\text{Re}_\text{mix}} = \sum_j \frac{\alpha_j}{\text{Re}_j}\]
Input parameters: Re_inv (shear and volume Reynolds numbers per fluid).
Additional geometric source terms appear with \(1/r\) factors in the continuity, momentum, and energy equations. Key modifications:
Viscous stress: \(\tau_{yy}\) includes \(v/r\) corrections:
\[\tau_{yy} = \mu\left(\frac{4}{3}\frac{\partial v}{\partial r} - \frac{2}{3}\frac{\partial u}{\partial x} - \frac{2}{3}\frac{v}{r}\right)\]
Source: src/simulation/m_bubbles_EE.fpp, src/simulation/m_bubbles.fpp
Modified mixture pressure:
\[p = (1 - \alpha)\,p_l + \alpha\left(\frac{R^3\,p_{bw}}{\bar{R}^3} - \frac{\rho\,R^3\,\dot{R}^2}{\bar{R}^3}\right)\]
Modified stiffened gas for the liquid phase:
\[\Gamma_l\,p_l + \Pi_{\infty,l} = \frac{1}{1 - \alpha}\left(E - \frac{1}{2}\rho\,|\mathbf{u}|^2\right)\]
Bubble wall pressure (polytropic):
\[p_{bw} = \left(p_0 + \frac{2\sigma}{R_0}\right)\left(\frac{R_0}{R}\right)^{3\gamma} - \frac{4\mu\,\dot{R}}{R} - \frac{2\sigma}{R}\]
Void fraction transport:
\[\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = \frac{3\,\alpha\,\bar{R}^2\,\dot{R}}{\bar{R}^3}\]
Number density conservation:
\[\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0\]
where \(n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}\).
Polydispersity (polydisperse = .true.): Log-normal PDF discretized into \(N_\text{bin}\) equilibrium radii with standard deviation poly_sigma, integrated via Simpson's rule.
\[R\,\ddot{R} + \frac{3}{2}\,\dot{R}^2 = \frac{p_{bw} - p_\infty}{\rho_l}\]
\[R\,\ddot{R}\left(1 - \frac{\dot{R}}{c}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3c}\right) = \frac{p_{bw} - p_\infty}{\rho_l}\left(1 + \frac{\dot{R}}{c}\right) + \frac{R\,\dot{p}_{bw}}{\rho_l\,c}\]
Enthalpy-based formulation with compressibility corrections via the Tait EOS:
\[R\,\ddot{R}\left(1 - \frac{\dot{R}}{C}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3C}\right) = H\left(1 + \frac{\dot{R}}{C}\right) + \frac{R\,\dot{H}}{C}\left(1 - \frac{\dot{R}}{C}\right)\]
where the enthalpy difference is:
\[H = \frac{n_\text{tait}(1 + B)}{n_\text{tait} - 1}\left[\left(\frac{p_{bw}}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} - \left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}}\right]\]
and the local liquid sound speed:
\[C = \sqrt{n_\text{tait}(1+B)\left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} + (n_\text{tait} - 1)\,H}\]
Internal bubble pressure ODE:
\[\dot{p}_b = \frac{3\gamma_b}{R}\left(-\dot{R}\,p_b + R_v\,T_{bw}\,\dot{m}_v + \frac{\gamma_b - 1}{\gamma_b}\,k_{bw}\left.\frac{\partial T}{\partial r}\right|_R\right)\]
Vapor mass flux:
\[\dot{m}_v = \frac{D\,\rho_{bw}}{1 - \chi_{vw}}\left.\frac{\partial \chi_v}{\partial r}\right|_R\]
Population balance equation:
\[\frac{\partial f}{\partial t} + \frac{\partial (f\,\dot{R})}{\partial R} + \frac{\partial (f\,\ddot{R})}{\partial \dot{R}} = 0\]
Moment transport:
\[\frac{\partial (n_\text{bub}\,\mu_i)}{\partial t} + \nabla \cdot (n_\text{bub}\,\mu_i\,\mathbf{u}) = n_\text{bub}\,\dot{\mu}_i\]
where moments \(\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}\).
CHyQMOM inversion recovers 4 quadrature nodes \((w_j, R_j, \dot{R}_j)\) from 6 moments via:
\[\bar{u} = \frac{\mu_{10}}{\mu_{00}}, \quad \bar{v} = \frac{\mu_{01}}{\mu_{00}}, \quad c_{20} = \frac{\mu_{20}}{\mu_{00}} - \bar{u}^2, \quad c_{11} = \frac{\mu_{11}}{\mu_{00}} - \bar{u}\bar{v}, \quad c_{02} = \frac{\mu_{02}}{\mu_{00}} - \bar{v}^2\]
Source: src/simulation/m_bubbles_EL.fpp
Volume-averaged carrier flow equations with bubble source terms:
Continuity:
\[\frac{\partial \rho_l}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l) = \frac{\rho_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right]\]
Momentum:
\[\frac{\partial (\rho_l\,\mathbf{u}_l)}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l \otimes \mathbf{u}_l + p\,\mathbf{I} - \boldsymbol{\tau}_l) = \frac{\rho_l\,\mathbf{u}_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{I} - \boldsymbol{\tau}_l)\]
Energy:
\[\frac{\partial E_l}{\partial t} + \nabla \cdot \bigl[(E_l + p)\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l\bigr] = \frac{E_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l)\]
The left-hand side is the standard conservation law for the liquid phase; the right-hand side source terms capture the effect of the bubbles on the host liquid.
Void fraction via regularization kernel:
\[\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)\]
where \(\delta_\sigma\) is a Gaussian kernel:
\[\delta_\sigma(\mathbf{r}) = \frac{1}{(2\pi\sigma^2)^{3/2}}\exp\!\left(-\frac{|\mathbf{r}|^2}{2\sigma^2}\right)\]
with \(\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble})\).
Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration.
Source: src/simulation/m_hypoelastic.fpp
Cauchy stress decomposition:
\[\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}\]
Elastic energy contribution to total energy:
\[E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}\]
Elastic stress evolution:
\[\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e\]
Source term:
\[\mathbf{S}^e = \rho\bigl(\mathbf{l} \cdot \boldsymbol{\tau}^e + \boldsymbol{\tau}^e \cdot \mathbf{l}^T - \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) + 2G\,\mathbf{D}^d\bigr)\]
where \(\mathbf{l} = \nabla \mathbf{u}\) is the velocity gradient and \(\mathbf{D}^d = \mathbf{D} - \frac{1}{3}\text{tr}(\mathbf{D})\,\mathbf{I}\) is the deviatoric strain rate.
Lie objective temporal derivative (Kelvin-Voigt):
\[\hat{\boldsymbol{\tau}}^e = \frac{D\boldsymbol{\tau}^e}{Dt} - \mathbf{l} \cdot \boldsymbol{\tau}^e - \boldsymbol{\tau}^e \cdot \mathbf{l}^T + \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) = 2G\,\mathbf{D}^d\]
This adds 6 additional transport equations in 3D (symmetric stress tensor: \(\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e\)).
Source: src/simulation/m_hyperelastic.fpp
Reference map evolution:
\[\frac{\partial (\rho\,\boldsymbol{\xi})}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\xi} \otimes \mathbf{u}) = 0\]
Deformation gradient from reference map:
\[\mathbf{F} = (\nabla \boldsymbol{\xi})^{-1}\]
Left Cauchy-Green tensor:
\[\mathbf{b} = \mathbf{F}\,\mathbf{F}^T\]
Neo-Hookean Cauchy stress:
\[\boldsymbol{\tau}^e = \frac{G}{J}\left(\mathbf{b} - \frac{\text{tr}(\mathbf{b})}{3}\,\mathbf{I}\right)\]
where \(J = \det(\mathbf{F})\).
Hyperelastic energy:
\[e^e = \frac{G}{2}\bigl(I_{\mathbf{b}} - 3\bigr), \qquad I_{\mathbf{b}} = \text{tr}(\mathbf{b})\]
Source: src/common/m_phase_change.fpp
\(N\)-fluid pressure-temperature equilibrium. The equilibrium condition is:
\[f(p) = \sum_i \alpha_i - 1 = 0\]
Temperature from energy conservation:
\[T = \frac{\rho e + p - \sum_i (\alpha_i \rho_i)\,q_{v,i}}{\sum_i (\alpha_i \rho_i)\,c_{v,i}\,\gamma_i}\]
Newton residual:
\[g(p) = \sum_i \frac{(\gamma_i - 1)\,(\alpha_i \rho_i)\,c_{v,i}}{(p + \pi_{\infty,i})} \cdot \frac{\rho e + p - \sum_j (\alpha_j \rho_j)\,q_{v,j}}{\sum_j (\alpha_j \rho_j)\,c_{v,j}\,\gamma_j}\]
Solved via Newton's method for the equilibrium pressure.
Two coupled equations for \((\alpha_1 \rho_1,\;p)\):
Gibbs free energy equilibrium (Clausius-Clapeyron):
\[F_1 = T\left[(c_{v,l}\gamma_l - c_{v,v}\gamma_v)(1 - \ln T) - (q'_l - q'_v) + c_{v,l}(\gamma_l - 1)\ln(p + \pi_{\infty,l}) - c_{v,v}(\gamma_v - 1)\ln(p + \pi_{\infty,v})\right] + q_{v,l} - q_{v,v} = 0\]
Energy conservation constraint:
\[F_2 = \rho e + p + m_l\,(q_{v,v} - q_{v,l}) - m_T\,q_{v,v} - m_{qD} + \frac{m_l\,(c_{v,v}\,\gamma_v - c_{v,l}\,\gamma_l) - m_T\,c_{v,v}\,\gamma_v - m_{cpD}}{m_l\left(\frac{c_{v,l}\,(\gamma_l - 1)}{p + \pi_{\infty,l}} - \frac{c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}}\right) + \frac{m_T\,c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}} + m_{cvgp}} = 0\]
where \(m_T\) is the total mass, \(m_l = \alpha_l \rho_l\) is the liquid partial density, and \(m_{qD}\), \(m_{cpD}\), \(m_{cvgp}\) are auxiliary thermodynamic sums over additional fluids (beyond the phase-changing pair).
Solved via 2D Newton-Raphson.
Source: src/common/m_chemistry.fpp
Species transport:
\[\frac{\partial (\rho_g\,Y_m)}{\partial t} + \frac{\partial (\rho_g\,u_i\,Y_m)}{\partial x_i} = W_m\,\dot{\omega}_m\]
Net production rate:
\[\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n\]
Reaction rate (law of mass action):
\[\mathcal{R}_n = k_n(T)\left[\prod_j \left(\frac{\rho_g\,Y_j}{W_j}\right)^{\nu'_{jn}} - \frac{1}{K_n}\prod_k \left(\frac{\rho_g\,Y_k}{W_k}\right)^{\nu''_{kn}}\right]\]
Arrhenius rate:
\[k_n(T) = A_n\,T^{b_n}\exp\!\left(-\frac{T_{a,n}}{T}\right)\]
Molecular diffusion (transport_model):
Enthalpy flux with diffusion:
\[q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k\]
Reaction mechanisms are code-generated via Pyrometheus (Cisneros-Garibay et al. [12]), which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms.
Source: src/simulation/m_surface_tension.fpp, src/simulation/include/inline_capillary.fpp
Color function advection:
\[\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0\]
Capillary stress tensor (CSF model):
\[\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)\]
In component form, with \(\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|\):
\[\Omega_{xx} = -\sigma\,(\hat{w}_y^2 + \hat{w}_z^2)\,\|\nabla c\|, \qquad \Omega_{xy} = \sigma\,\hat{w}_x\,\hat{w}_y\,\|\nabla c\|\]
The capillary stress divergence is added to the momentum and energy equations. The total energy equation becomes:
\[\frac{\partial (\rho E + \varepsilon_0)}{\partial t} + \nabla \cdot \bigl[(\rho E + \varepsilon_0 + p)\,\mathbf{u} + (\boldsymbol{\Omega} - \boldsymbol{\tau}^v) \cdot \mathbf{u}\bigr] = 0\]
Capillary mixture energy:
\[\varepsilon_0 = \sigma\,\|\nabla c\|\]
Continuity:
\[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0\]
Momentum:
\[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \left[\rho\,\mathbf{u} \otimes \mathbf{u} + \left(p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{I} - \mathbf{B} \otimes \mathbf{B}\right] = 0\]
Energy:
\[\frac{\partial \mathcal{E}}{\partial t} + \nabla \cdot \left[\left(\mathcal{E} + p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\right] = 0\]
Induction:
\[\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0\]
Total energy:
\[\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}\]
Fast magnetosonic speed:
\[c_f = \sqrt{\frac{1}{2}\left(c_s^2 + v_A^2 + \sqrt{(c_s^2 + v_A^2)^2 - 4\,c_s^2\,v_A^2\cos^2\theta}\right)}\]
Alfven speed:
\[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\]
Uses the HLLD Riemann solver (riemann_solver = 4). Hyperbolic divergence cleaning (hyper_cleaning = .true.) via the GLM method (Dedner et al. [15]).
Conserved variables:
\[\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T\]
where:
\[D = \Gamma\,\rho, \qquad \mathbf{m} = \Gamma^2\rho h\,\mathbf{u} + |\mathbf{B}|^2\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\]
\[\tau = \Gamma^2\rho h - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{u}|^2|\mathbf{B}|^2 - (\mathbf{B} \cdot \mathbf{u})^2}{2} - \Gamma\,\rho\]
Primitive recovery uses Newton-Raphson on the nonlinear conserved-to-primitive relation.
Source: src/simulation/m_igr.fpp
Modified momentum with entropic pressure \(\Sigma\)**:**
\[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \bigl[\rho\,\mathbf{u} \otimes \mathbf{u} + (p + \Sigma)\,\mathbf{I} - \boldsymbol{\tau}\bigr] = 0\]
Elliptic PDE for \(\Sigma\)**:**
\[\alpha\left[\text{tr}(\nabla \mathbf{u})^2 + \text{tr}^2(\nabla \mathbf{u})\right] = \frac{\Sigma}{\rho} - \alpha\,\nabla \cdot \left(\frac{\nabla \Sigma}{\rho}\right)\]
where \(\alpha \sim \Delta x^2\) (regularization strength proportional to mesh spacing squared):
\[\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2\]
RHS strain-rate source (3D):
\[\text{RHS} = \alpha\left[2\left(\frac{\partial u}{\partial y}\frac{\partial v}{\partial x} + \frac{\partial u}{\partial z}\frac{\partial w}{\partial x} + \frac{\partial v}{\partial z}\frac{\partial w}{\partial y}\right) + \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial w}{\partial z}\right)^2 + (\nabla \cdot \mathbf{u})^2\right]\]
Iterative solver: Jacobi (igr_iter_solver = 1) or Gauss-Seidel (igr_iter_solver = 2), up to num_igr_iters iterations (default 5).
Uses Lax-Friedrichs flux (replaces WENO + Riemann solver).
Source: src/simulation/m_body_forces.fpp
Time-dependent acceleration:
\[a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)\]
Momentum source:
\[\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)\]
Energy source:
\[\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)\]
Source: src/simulation/m_acoustic_src.fpp
Source terms added to the RHS of the governing equations.
Formulation follows Maeda and Colonius [29].
Discrete delta function (spatial support):
\[\delta_h(r) = \frac{1}{(2\pi\sigma^2)^{d/2}}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)\]
Forcing form (added to conservative variables):
\[\mathbf{s}_\text{ac} = \Omega_\Gamma\,f(t)\left[\frac{1}{c_0},\;\cos\theta,\;\sin\theta,\;\frac{c_0^2}{\gamma - 1}\right]^T\]
Temporal profiles:
Spatial supports: planar, spherical transducer, cylindrical transducer, transducer array (arcuate, annular, circular).
Source: src/simulation/m_weno.fpp
Weighted sum of candidate polynomials at cell interfaces:
\[f_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}\]
WENO-JS (Jiang and Shu [22], default):
\[\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}\]
where \(d_r\) are ideal weights, \(\beta_r\) are smoothness indicators, and \(\varepsilon\) is a small regularization parameter (weno_eps).
WENO-M (mapped_weno = .true.): Henrick et al. [21] mapped weights for improved accuracy at critical points:
\[\omega_M^{(r)} = \frac{d_r\bigl(1 + d_r - 3\omega_0^{(r)} + (\omega_0^{(r)})^2\bigr)\,\omega_0^{(r)}}{d_r^2 + \omega_0^{(r)}(1 - 2d_r)}, \qquad \omega^{(r)} = \frac{\omega_M^{(r)}}{\sum_s \omega_M^{(s)}}\]
WENO-Z (wenoz = .true.): Borges et al. [6] improved weights with global smoothness measure:
\[\alpha_r = d_r\left(1 + \left(\frac{\tau}{\beta_r + \varepsilon}\right)^q\right), \qquad \tau = |\beta_0 - \beta_{k-1}|\]
The parameter \(q\) controls the convergence rate at critical points (typically \(q = 1\) for fifth-order reconstruction, as used in MFC).
TENO (teno = .true.): Fu et al. [16] targeted ENO with smoothness threshold \(C_T\) (teno_CT):
\[\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}\]
If \(\xi_r < C_T\), set \(\alpha_r = 0\) (stencil excluded).
Primitive variable reconstruction is used to avoid spurious oscillations at interfaces.
Source: src/simulation/m_muscl.fpp
\[q_L(j) = q(j) - \frac{1}{2}\,\phi(r)\,\Delta q, \qquad q_R(j) = q(j) + \frac{1}{2}\,\phi(r)\,\Delta q\]
Five slope limiters (muscl_lim):
where \(a\) and \(b\) are the left and right slope differences.
THINC interface compression (int_comp = .true.): applies hyperbolic tangent compression near material interfaces:
\[q_\text{THINC} = q_\min + \frac{q_\max}{2}\left(1 + \text{sign}(s)\,\frac{\tanh(\beta) + A}{1 + A\,\tanh(\beta)}\right)\]
where \(A = \frac{\exp(\text{sign}(s)\,\beta\,(2C - 1))\,/\,\cosh(\beta) - 1}{\tanh(\beta)}\) and \(\beta\) controls compression steepness.
5th-order or 3rd-order polynomial interpolation without WENO nonlinear weights, using Lax-Friedrichs numerical flux. Stencil coefficients:
5th-order: \(\text{coeff}_L = [-3/60,\;27/60,\;47/60,\;-13/60,\;2/60]\)
3rd-order: \(\text{coeff}_L = [2/6,\;5/6,\;-1/6]\)
Source: src/simulation/m_riemann_solvers.fpp
\[\mathbf{F}^\text{HLL} = \frac{S_R\,\mathbf{F}_L - S_L\,\mathbf{F}_R + S_L\,S_R\,(\mathbf{q}_R - \mathbf{q}_L)}{S_R - S_L}\]
with wave speed estimates \(S_L = \min(u_L - c_L,\;u_R - c_R)\), \(S_R = \max(u_R + c_R,\;u_L + c_L)\).
Four-state solver resolving the contact discontinuity. Star-state satisfies:
\[p_L^* = p_R^* = p^*, \qquad u_L^* = u_R^* = u^*\]
Iterative exact Riemann solver.
Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (Miyoshi and Kusano [35]). The Riemann fan is divided by outer wave speeds \(S_L\), \(S_R\), Alfven speeds \(S_L^*\), \(S_R^*\), and a middle contact \(S_M\):
\[S_M = \frac{(S_R - u_R)\rho_R u_R - (S_L - u_L)\rho_L u_L - p_{T,R} + p_{T,L}}{(S_R - u_R)\rho_R - (S_L - u_L)\rho_L}\]
\[S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}\]
where \(p_T = p + |\mathbf{B}|^2/2\) is the total (thermal + magnetic) pressure. Continuity of normal velocity and total pressure is enforced across the Riemann fan.
Source: src/simulation/m_time_steppers.fpp
RK1 (Forward Euler):
\[\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\]
RK2:
\[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\]
\[\mathbf{q}^{n+1} = \frac{1}{2}\mathbf{q}^n + \frac{1}{2}\mathbf{q}^{(1)} + \frac{1}{2}\Delta t\,\mathbf{L}(\mathbf{q}^{(1)})\]
RK3 (SSP):
\[\mathbf{q}^{(1)} = \mathbf{q}^n + \Delta t\,\mathbf{L}(\mathbf{q}^n)\]
\[\mathbf{q}^{(2)} = \frac{3}{4}\mathbf{q}^n + \frac{1}{4}\mathbf{q}^{(1)} + \frac{1}{4}\Delta t\,\mathbf{L}(\mathbf{q}^{(1)})\]
\[\mathbf{q}^{n+1} = \frac{1}{3}\mathbf{q}^n + \frac{2}{3}\mathbf{q}^{(2)} + \frac{2}{3}\Delta t\,\mathbf{L}(\mathbf{q}^{(2)})\]
Embedded RK pairs for error estimation with Hairer-Norsett-Wanner algorithm for initial step size.
For equations of the form \(\partial \mathbf{q}/\partial t = -\nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{s}(\mathbf{q})\), the Strang splitting scheme integrates three sub-equations per time step:
\[\frac{\partial \mathbf{q}'}{\partial t} = \mathbf{s}(\mathbf{q}'), \quad t \in [0,\,\Delta t/2], \quad \mathbf{q}'(0) = \mathbf{q}^n\]
\[\frac{\partial \mathbf{q}''}{\partial t} = -\nabla \cdot \mathbf{F}(\mathbf{q}''), \quad t \in [0,\,\Delta t], \quad \mathbf{q}''(0) = \mathbf{q}'(\Delta t/2)\]
\[\frac{\partial \mathbf{q}'''}{\partial t} = \mathbf{s}(\mathbf{q}'''), \quad t \in [0,\,\Delta t/2], \quad \mathbf{q}^{n+1} = \mathbf{q}'''(\Delta t/2)\]
The stiff bubble source terms are integrated using an adaptive embedded Runge-Kutta scheme with error control.
Inviscid:
\[\Delta t = \text{CFL} \cdot \min_{i,j,k}\left(\frac{\Delta x_\text{cell}}{|u| + |v| + |w| + c}\right)\]
Viscous:
\[\Delta t_v = \text{CFL} \cdot \min\left(\frac{\Delta x^2}{\nu}\right)\]
where \(c\) is the mixture sound speed and \(\nu\) is the kinematic viscosity.
Source: src/common/m_finite_differences.fpp
Used for viscous fluxes and velocity gradients.
1st-order:
\[f'_i = \frac{f_{i+1} - f_i}{x_{i+1} - x_i}\]
2nd-order centered:
\[f'_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1} - x_{i-1}}\]
4th-order centered:
\[f'_i = \frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{-x_{i+2} + 8x_{i+1} - 8x_{i-1} + x_{i-2}}\]
Boundary-adjusted one-sided stencils (3rd-order):
\[f'_i\big|_\text{left} = \frac{-3f_i + 4f_{i+1} - f_{i+2}}{x_{i+2} - x_i}\]
\[f'_i\big|_\text{right} = \frac{3f_i - 4f_{i-1} + f_{i-2}}{x_i - x_{i-2}}\]
Source: src/simulation/m_cbc.fpp, src/simulation/m_compute_cbc.fpp, src/common/m_constants.fpp
| BC Code | Type |
|---|---|
| -1 | Periodic |
| -2 | Reflective |
| -3 | Ghost cell extrapolation |
| -4 | Riemann extrapolation |
| -14 | Axis symmetry |
| -15 | Slip wall |
| -16 | No-slip wall |
Characteristic decomposition:
\[\frac{\partial \mathbf{q}_c}{\partial t} + \mathbf{R}_x\,\boldsymbol{\Lambda}_x\,\mathbf{L}_x\,\frac{\partial \mathbf{q}_c}{\partial x} = 0\]
Characteristic amplitudes \(\mathcal{L}_1\) through \(\mathcal{L}_5\) define wave interactions at boundaries:
\[\mathcal{L}_1 = (u - c)\left(\frac{\partial p}{\partial x} - \rho\,c\,\frac{\partial u}{\partial x}\right), \qquad \mathcal{L}_2 = u\left(c^2\,\frac{\partial \rho}{\partial x} - \frac{\partial p}{\partial x}\right)\]
\[\mathcal{L}_3 = u\,\frac{\partial v}{\partial x}, \qquad \mathcal{L}_4 = u\,\frac{\partial w}{\partial x}, \qquad \mathcal{L}_5 = (u + c)\left(\frac{\partial p}{\partial x} + \rho\,c\,\frac{\partial u}{\partial x}\right)\]
For non-reflecting boundaries, the incoming wave amplitudes are set to zero.
GRCBC (incoming wave from ghost point):
\[\mathcal{L} = -\boldsymbol{\Lambda}^o\,\mathbf{L}_x\,\frac{\partial \mathbf{q}_c}{\partial x} + n_x\,\lambda^i\,\mathbf{L}_x\,\frac{\mathbf{P}\,(\mathbf{q}_p^{(g)} - \mathbf{q}_p)}{\Delta}\]
8 types: slip wall (-5), non-reflecting buffer (-6), non-reflecting sub. inflow (-7), non-reflecting sub. outflow (-8), force-free sub. outflow (-9), constant-pressure sub. outflow (-10), supersonic inflow (-11), supersonic outflow (-12).
Source: src/simulation/m_ibm.fpp
Ghost-cell IBM with level set field \(\phi\).
Image point:
\[\mathbf{x}_{ip} = \mathbf{x}_{gp} + 2\,\phi(\mathbf{x}_{gp})\,\hat{\mathbf{n}}\]
Velocity boundary conditions:
Neumann BC for pressure, density, and volume fractions.
Supports STL/OBJ geometry import with ray tracing for inside/outside determination.
Chen correction (low_Mach = 1, Chen et al. [11]): anti-dissipation pressure correction (APC) added to the HLLC flux:
\[p_d = \frac{\rho_L\,\rho_R\,(S_R - u_R)(S_L - u_L)}{\rho_R(S_R - u_R) - \rho_L(S_L - u_L)}\]
\[\text{APC}_{p\mathbf{u}} = (z - 1)\,p_d\,\hat{\mathbf{n}}, \qquad \text{APC}_{pE} = (z - 1)\,p_d\,S_*\]
The corrected HLLC flux in the star region becomes \(\mathbf{F}^{*} + \text{APC}\).
Thornber correction (low_Mach = 2, Thornber et al. [49]): modifies the reconstructed velocities at cell interfaces:
\[u'_L = \frac{u_L + u_R}{2} + z\,\frac{u_L - u_R}{2}, \qquad u'_R = \frac{u_L + u_R}{2} - z\,\frac{u_L - u_R}{2}\]
\[z = \min\!\left(\max\!\left(\frac{|u_L|}{c_L},\;\frac{|u_R|}{c_R}\right),\;1\right)\]
This reduces numerical dissipation at low Mach numbers while recovering the standard scheme at high Mach numbers.
Additionally, the artificial Mach number technique (pi_fac) modifies \(\pi_\infty\) to artificially increase the local Mach number.
Volume fraction limiting: enforces \(0 \le \alpha_k \le 1\) and rescales as:
\[\alpha_k \leftarrow \frac{\alpha_k}{\sum_j \alpha_j}\]
Advective flux limiting based on local volume fraction gradient \(\chi\) to prevent spurious oscillations near material interfaces.