Spatiotemporal Coupling in High-Res Ghost Fluid Method

• The paper presents a novel Ghost Fluid Method integrating nonlinear geometrical optics with Lax–Wendroff/C–K expansion to accurately simulate multi-material compressible flows.
• It achieves high resolution by reconstructing ghost states as linearly varying functions, preserving local thermodynamic gradients and multidimensional flux effects.
• Numerical experiments show enhanced accuracy in capturing shock interactions and interface dynamics compared to classical Riemann-based solvers.

The spatiotemporal coupling high-resolution Ghost Fluid Method (GFM) represents an advanced computational approach for simulating two-dimensional compressible multi-medium flows with source terms. Distinct from classical Riemann-problem-based solvers, this framework systematically integrates nonlinear geometrical optics for thermodynamic entropy evolution and a Lax–Wendroff/Cauchy–Kowalevski (C–K) spatiotemporal expansion, explicitly accounting for tangential fluxes and source term effects. This methodology achieves higher thermodynamic compatibility, multidimensional fidelity, and numerical accuracy, addressing critical limitations of conventional ghost fluid schemes in multidimensional, multi-material, and source-driven contexts (Huo, 23 Jan 2026).

1. Governing Equations for Compressible Multi-Medium Flow

The method considers the conservative formulation for two-dimensional compressible flow involving $K$ materials, potentially including axisymmetric source terms: $$\partial_t u + \partial_x F(u) + \partial_y G(u) = S(x,y,u),\quad (x,y)\in\Omega,\, t>0$$ with

$$u(x,y,0) = \begin{cases} u_1(x,y), & (x,y)\in\Omega_1^0 \ u_2(x,y), & (x,y)\in\Omega_2^0 \end{cases}$$

for initial partitioning, and closure via a material-dependent equation of state (EOS). The conserved variables and fluxes are

$u=[\rho,\, \rho u_x,\, \rho u_y,\, E]^T$

$$F(u)= [\rho u_x,\, \rho u_x^2 + p,\, \rho u_x u_y,\, u_x(E + p)]^T$$

$$G(u)= [\rho u_y,\, \rho u_x u_y,\, \rho u_y^2 + p,\, u_y(E + p)]^T$$

$E = e + 0.5(u_x^2 + u_y^2)$

with $e = e_k(\rho, p)$ in each region $k$. Axisymmetric flows adopt the general stiffened gas EOS

$$e_k = \frac{p_k + \gamma_k p_\infty^k}{(\gamma_k-1)\rho_k},\quad k=1,2$$

The source term $S(x,y,u)$ can encode geometric and/or physical reaction contributions.

2. Nonlinear Geometrical Optics and Thermodynamic Consistency

Classical Riemann solvers operate under the piecewise-constant assumption for data on either side of discontinuities. Consequently, simulated rarefactions remain strictly isentropic, with entropy variations only permitted as discrete jumps across shocks—violating continuous thermodynamic compatibility in compressible flows.

The Generalized Riemann Problem (GRP) overcomes this by employing piecewise-linear initial data: $u(x,0) = u_k + \sigma_k(x-x_0),\quad x\gtrless x_0$ This supports curved rarefactions carrying smooth, continuous entropy gradients. Thermodynamic consistency is enforced through the Gibbs relation,

$T\, d s = d e - p\, d (1/\rho)$

and explicit tracking of entropy transport is achieved by transforming to characteristic coordinates $(\xi, \tau)$ via nonlinear geometrical optics. For genuinely nonlinear waves,

$$\frac{D u_x}{Dt} + \frac{D p}{Dt} = \frac{D \psi}{Dt}$$

where the total derivative $D/Dt = \partial_t + u_x \partial_x$ and $\psi$ encodes the Riemann invariants. This formulation enables physically correct entropy propagation within the numerical solver (Huo, 23 Jan 2026).

3. Lax–Wendroff/Cauchy–Kowalevski Spatiotemporal Flux Construction

High temporal resolution and correct multidimensionality are achieved through a one-step Lax–Wendroff (or C–K) procedure in the finite-volume update. For the numerical flux in the $x$-direction at face $(i+\tfrac{1}{2},j)$: $$u_{i+\frac{1}{2},j}^{n+\frac{1}{2}} = u_{i+\frac{1}{2},j}^n + \frac{\Delta t}{2}\left(\partial_t u\right)_{i+\frac{1}{2},j}^n + O(\Delta t^2)$$ Here,

$$\partial_t u = -A(u)\,\partial_x u - B(u)\,\partial_y u + S$$

where $A = \partial F/\partial u$, $B = \partial G/\partial u$. In the linearized GRP framework, $A$ is decomposed into $A = R \Lambda R^{-1}$, and

$$(\partial_t u)_{i+\frac{1}{2},j}^{n,x} = -R \Lambda^+ R^{-1} (\partial_x u)|_{i+\frac{1}{2}^-} - R \Lambda^- R^{-1} (\partial_x u)|_{i+\frac{1}{2}^+}$$

Tangential flux $G$ and source term $S$ are embedded via

$$H(x_{i+\frac{1}{2}},y_j,t_n) = -R I^+ R^{-1}[(\partial_y G-S)]|_{i+\frac{1}{2}^-} - R I^- R^{-1}[(\partial_y G-S)]|_{i+\frac{1}{2}^+}$$

with $I^{\pm} = (1 \pm \operatorname{sign}\Lambda)/2$, so the full time-derivative is

$$(\partial_t u)_{i+\frac{1}{2},j}^n = (\partial_t u)_{i+\frac{1}{2},j}^{n,x} + H$$

The time-averaged flux is then

$$\frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} F(u(x_{i+\tfrac{1}{2}},y_j,t))\, dt \approx F(u_{i+\tfrac{1}{2},j}^{n+\tfrac{1}{2}})$$

This procedure markedly improves accuracy in time, with error scaling as $O(\Delta t^2)-O(\Delta t^3)$ versus $O(|\Delta u|)$ for classical Riemann-fluxes.

4. Ghost Fluid Method with Spatiotemporal Coupling

The core of the GFM framework relies on using a level-set $\phi(x,y,t)$ to represent the interface $\Gamma = \{\phi=0\}$, defining narrow “ghost” regions $\Omega_k^{g}$ adjacent to each real fluid domain. The critical methodological distinction is that, whereas classical Riemann-based GFMs reconstruct ghost states as constants, the GRP-based GFM reconstructs ghost states as linearly varying functions, thus preserving local gradients.

Standard Riemann Problem-Based GFM

• Solve a local 1D two-medium Riemann problem normal to the interface, imposing continuity of pressure $[p] = 0$ and normal velocity $[u\cdot n] = 0$ across $\xi=0$, returning post-wave state $u_k^*$.
• Construct a constant ghost state in Cartesian coordinates.

GRP-Based GFM (Spatiotemporal Coupling)

• Solve a local two-medium GRP with piecewise-linear data along the interface normal $\xi$:

$$\tilde{u}(\xi,0) = \tilde{u}_k + \xi (\partial_\xi\tilde{u})_k$$

• Incorporate tangential $\partial_\eta G$ and source term $-\partial_\eta G + S$ into the GRP.
• The GRP solver yields both post-wave states $\tilde{u}_k^*$ and their temporal derivatives $(\partial_t\tilde{u})_k^*$, enabling recovery of spatial derivatives through the local PDE.
• The ghost fluid is reconstructed as

$$u_k^g(x,y) = u_k^* + (x-x_\Gamma)(\partial_x u)_k^* + (y-y_\Gamma)(\partial_y u)_k^*$$

which preserves both thermodynamic and multidimensional effects.

Algorithmic Outline (per time step)

1. Evolve the level set: $\phi^n \rightarrow \phi^{n+1}$ via $\partial_t\phi + u\cdot\nabla\phi=0$.
2. Identify real-fluid and ghost regions.
3. For each interfacial cell:
   • Compute normal and tangential directions,
   • Perform a local GRP along the normal with least-squares slope data,
   • Extract post-wave and slope derivatives,
   • Construct linear ghost states using the above relation.
4. Reconstruct cell-centered slopes.
5. Apply the Lax–Wendroff/C–K flux across all faces.
6. Update finite-volume averages.
7. Advance $t^n \to t^{n+1}$ subject to a CFL stability condition.

5. Validation and Performance in Numerical Experiments

The method has been validated on classical two-medium flow test cases.

Spherical Helium-Air Bubble Interaction (Haas–Sturtevant 1987)

• Domain: $[0,\,0.55]\times[-0.0445,\,0.0445]$, grid: $550\times90$.
• Initial: Helium bubble radius 0.025 at (0.425, 0), Mach 1.25 shock.
• EOS: $\gamma_{He}=1.648$, $\gamma_{air}=1.4$, $p_\infty=0$.
• Boundaries: wall (left/top/bottom), piston (right).
• Results: The GRP-based method reproduces fine-scale instabilities (rollups), sharp interface curvature, and accurate shock/bubble positioning, exceeding both RP-based schemes and experiment in resolution.

Spherical Air Bubble Collapse in Water (Bourne–Field 1992)

• Domain: $[-6,6]\times[0,15]$, grid: $120\times150$.
• Initial: Air bubble at $z=9$ cm, water shock at $z=12$ cm.
• EOS: stiffened gas for water, ideal gas for air.
• Results: The GRP-based solver sharply resolves jet formation and bubble deformation; conventional RP-based methods display excessive numerical diffusion.

No formal $L_2$ error norms or grid-convergence rates are reported, but temporal flux error is $O(\Delta t^2)$–$O(\Delta t^3)$ for the GRP method, compared to $O(|\Delta u|)$ for the classical RP method.

6. Properties, Limitations, and Prospects

Key advantages include:

• Thermodynamic Consistency: Continuous entropy transport is enforced via nonlinear geometrical optics.
• Multidimensional Fidelity: Tangential fluxes and source terms are embedded in the numerical fluxes through Lax–Wendroff/C–K construction.
• High Temporal-Spatial Accuracy: Second to third-order accuracy in time and robust interface capturing through linear ghost states, suppressing spurious oscillations and interface smearing.

Limitations and future considerations comprise:

• Computational Overhead: Each interface requires a local GRP solve, increasing computational demand relative to simpler methods.
• Extension to Reactive Flows: Incorporating chemical source-term coupling into the GRP.
• Generalization: Extension to three dimensions and unstructured mesh geometries remains to be addressed.

The framework thus establishes a robust, multidimensional, thermodynamically compatible, high-resolution GFM for compressible multi-material flows with source terms, leveraging the synergy of advanced nonlinear wave analysis and high-order temporal discretization (Huo, 23 Jan 2026).