A Unified Numerical Framework for Turbulent Convection and Phase-Change Dynamics in Coupled Fluid-Porous Systems

This work presents a unified numerical framework for simulating incompressible flows within the coupled fluid-porous-medium system and involving heat and solute transport and phase-changing process. A complete set of governing equations is established based on the Darcy-Brinkman equation, the advection-diffusion equations for heat and solute, and a phase field equation describing the evolution of porous medium. Phase-changing process and relevant influences are incorporated as corresponding source terms. A numerical method is then developed to solve the governing equations. Several different types of model problems are simulated with the numerical method. For the incompressible flows inside a coupled fluid-porous-medium system, the channel turbulence over a porous substrate and the thermal convection in a two-layer system are simulated. For the phase-changing flows, the one-dimensional Stefan problem and the two-dimensional flow of pure water freezing are tested. The results agree with the existing simulations. Finally, the full solver is used to simulate the growth of mushy ice during seawater freezing, which can successfully reproduce the experimental results at the exactly same conditions. Therefore, the developed framework provides a versatile and reliable tool for studying complex multiphase, multi-component transport phenomena in fluid-porous-medium systems involving solid-liquid phase change.

💡 Research Summary

**
This paper introduces a comprehensive numerical framework capable of simultaneously handling turbulent convection, heat and solute transport, and solid‑liquid phase change in systems where a clear fluid region is coupled with a porous medium. The authors adopt a single‑domain approach based on the Darcy‑Brinkman formulation, which embeds both viscous and Darcy resistance terms within a unified momentum equation. By introducing a porosity field ϕ that varies smoothly from 0 (pure solid) to 1 (pure liquid), the model automatically guarantees continuity of velocity and pressure across fluid‑porous interfaces without the need for ad‑hoc interface conditions.

The governing equations consist of:

1. Modified phase‑field equation – a second‑order advection‑diffusion‑reaction form, ∂tϕ + U·∇ϕ = κϕ∇²ϕ + G(T − Tϕ), where Tϕ = T₀ − λS represents the local liquidus temperature that depends on solute concentration S. The double‑well potential and fourth‑order term typical of Cahn‑Hilliard models are omitted, allowing larger time steps while still capturing the essential kinetics of melting and solidification.

2. Unified momentum equation – ∂tU = −U·∇(U/ϕ) − (1/ρ₀)∇P + ν∇²U − ϕK(ϕ)U − gϕρ′/ρ₀ e_z. Here K(ϕ) is a permeability function of porosity, ν the kinematic viscosity, and ρ′ the buoyancy‑induced density anomaly. When ϕ = 1 the equation reduces to the Navier‑Stokes form; when ϕ ≈ 0 the Darcy law dominates.

3. Energy equation – expressed as ∂tT = −(c_f/c_m)U·∇T + (1/c_m)∇·(k_e∇T) − (L/c_m)∂tϕ. The effective heat capacity c_m = ϕc_l + (1−ϕ)c_s and effective conductivity k_e = ϕk_l + (1−ϕ)k_s + ε_disp(ϕ)c_m|U| combine solid and liquid contributions and include a velocity‑dependent thermal dispersion term ε_disp(ϕ) that vanishes in pure fluid regions.

4. Solute transport – ∂tS = −(1/ϕ)U·∇S + (1/ϕ)∇·(κ_eS∇S) − S∂tϕ, with κ_eS = κ_S ϕ^τ, τ = ϕ − γ. The authors adopt γ = 1, yielding a second‑order power‑law (Archie‑type) relationship that strongly suppresses diffusion in low‑porosity mushy zones, reflecting the tortuous pathways typical of sea‑ice mushy layers.

All equations are nondimensionalized, introducing Rayleigh, Prandtl, Schmidt, Darcy, and Stefan numbers, among others, to capture the relative importance of buoyancy, viscous, diffusive, and phase‑change effects.

Numerical methodology – Temporal integration uses a second‑order Runge‑Kutta scheme combined with operator splitting, allowing the advection‑diffusion, phase‑field, and momentum sub‑problems to be solved sequentially. A projection (fractional‑step) method enforces incompressibility (∇·U = 0). Spatial discretization is performed on structured grids with central differences for diffusion and a staggered arrangement for velocity–pressure coupling. The framework accommodates Dirichlet, Neumann, and periodic boundary conditions for all fields, and the porosity field itself dictates the transition between fluid and porous physics, eliminating the need for explicit interface treatment.

Validation and application cases – The authors present a hierarchy of benchmark problems:

• Channel turbulence over a porous substrate – Mean velocity profiles, Reynolds stresses, and wall shear stress match direct‑numerical‑simulation (DNS) data, confirming that the Darcy‑Brinkman formulation correctly captures the interaction between turbulent eddies and the underlying porous layer.

• Two‑layer thermal convection – Simulations of a fluid layer atop a porous slab reproduce the expected Nusselt‑Rayleigh scaling and show smooth temperature and velocity continuity across the interface, validating the effective conductivity model that includes thermal dispersion.

• One‑dimensional Stefan problem – Both fixed‑boundary and moving‑boundary versions are solved; interface positions and temperature fields converge with second‑order accuracy and agree with analytical solutions and high‑resolution finite‑difference benchmarks.

• Two‑dimensional pure‑water freezing – The evolution of the solid–liquid front from an initial circular melt to a planar ice front is captured with quantitative agreement to laboratory visualizations, demonstrating the capability of the simplified phase‑field to handle realistic interface kinetics without spurious oscillations.

• Mushy‑ice growth during seawater freezing – This most demanding test couples temperature, salinity, and porosity. The model reproduces experimental measurements of mushy‑layer thickness, porosity distribution, and brine rejection rates under the same cooling conditions. The results confirm that the combined effects of solute‑dependent liquidus, Archie’s diffusion suppression, and velocity‑dependent thermal dispersion are essential for accurate sea‑ice modeling.

Discussion – The unified framework offers several notable advantages: (i) it eliminates the need for separate fluid and porous solvers and the associated interface‑matching procedures; (ii) the second‑order phase‑field reduces computational cost while preserving essential physics of melting/solidification; (iii) the inclusion of thermal dispersion and tortuosity‑based solute diffusion extends applicability to high‑porosity, high‑velocity regimes such as magma chambers or industrial casting. Limitations are acknowledged: the current implementation relies on the Boussinesq approximation and linear equations of state for density, which may be insufficient for strongly stratified or high‑pressure environments. Moreover, three‑dimensional complex geometries and unstructured meshes have not yet been demonstrated, and the model assumes local thermal equilibrium between solid and liquid phases.

Conclusion – By integrating Darcy‑Brinkman momentum, a simplified phase‑field, and comprehensive heat‑solute transport models, the authors deliver a versatile, validated numerical tool for studying coupled turbulent convection and phase‑change phenomena in fluid‑porous systems. The framework successfully bridges the gap between detailed pore‑scale physics and macroscopic engineering applications, and it paves the way for future extensions such as non‑Boussinesq density formulations, adaptive mesh refinement, and high‑performance parallel implementations targeting geophysical, metallurgical, and cryospheric problems.