Impact statement
Porous media play a vital role in various engineering applications, including combustion systems, heat exchangers, aero-acoustics and filtration. However, interactions between their complex geometries and turbulent flows pose significant challenges for numerical simulation, especially in mesh generation and time-step restrictions. This work introduces a high-fidelity simulation framework that overcomes these challenges by enabling accurate low-Mach-number simulations of turbulent flows through porous structures, while significantly reducing meshing complexity and computational cost. Additionally, the introduction of a synthetic turbulence injection at the inlet provides realistic inflow conditions, enhancing the practical relevance of simulations. These capabilities are particularly relevant for developing ultra-low emission combustion technologies, where turbulence–porous medium interactions critically influence flame stability, efficiency and pollutant formation.
1. Introduction
Porous media play a crucial role in various industries, including energy, chemical processing and environmental engineering, due to their ability to influence the flow dynamics and for example promote heat and mass transfer or optimise fluid mixing. These properties make them essential in applications such as filtration, catalysis (Lucci et al., Reference Lucci, Della Torre, Montenegro, Kaufmann and Eggenschwiler2017) or heat exchangers (Rashidi et al., Reference Rashidi, Kashefi, Kim and Samimi-Abianeh2019). In turbulent combustion systems, porous media contribute to ultra-low emission technologies, such as heterogeneous combustion (Trimis & Durst, Reference Trimis and Durst1996; Wood & Harris, Reference Wood and Harris2008, Boigné et al., Reference Boigné, Muhunthan, Mohaddes, Wang, Sobhani, Hinshaw and Ihme2019), by improving mixing and stabilising flames.
Furthermore, porous media are widely studied for noise attenuation in aeronautical applications, such as liners in engine nacelle (Ma & Su, Reference Ma and Su2020) and trailing edges of airfoils (Teruna et al., Reference Teruna, Manegar, Avallone, Ragni, Casalino and Carolus2020). However, their intricate flow dynamics, particularly at high Reynolds numbers, presents significant experimental and computational challenges (Wood et al., Reference Wood, He and Apte2020; Jin & Kuznetsov, Reference Jin and Kuznetsov2024).
Despite the widespread use of porous media, their geometric complexity and specific physical behaviours in turbulent flows make accurate numerical simulations challenging, even with high performance computing. Experimental studies are all the more difficult due to limited optical access to the pore structures, which hinders direct observations without disturbing the natural flow. As a result, many historical studies have focused on laminar flow regimes and simplified geometries, such as packed spheres or two-dimensional (2-D) porous structures made of circular or squared cylinders, especially through numerical approaches (Jin et al., Reference Jin and Kuznetsov2017; Srikanth et al., Reference Srikanth, Huang, Su and Kuznetsov2021; Rao & Jin, Reference Rao and Jin2022). While several macroscopic models have been proposed, as reviewed by Jin and Kuznetsov (Reference Jin and Kuznetsov2024), present understanding of turbulence in porous media remains limited due to the lack of high-fidelity data, which hinders the development of ad hoc turbulence models.
Modelling in porous media typically relies on macroscopic equations, derived via time, volume or combined averaging of the microscopic governing equations, leading to the RANS (Reynolds-averaged Navier–Stokes (NS)), VANS (volume-averaged NS) and DANS (double-averaged NS) equations. Initial approaches employed traditional RANS models, such as the
$\kappa$
-
$\varepsilon$
(Kuwahara et al., Reference Kuwahara, Kameyama, Yamashita and Nakayama1998) or
$v^2$
-
$f$
models (Kazerooni & Hannani, Reference Kazerooni and Hannani2007), where
$\kappa$
is the turbulent kinetic energy,
$\varepsilon$
is its dissipation rate,
$v^2 \equiv v^{\prime2}$
is the squared wall-normal velocity fluctuation, and
$f$
is the elliptic-relaxation function accounting for pressure-strain redistribution. Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) examined flows in rough-walled channels, where the walls were covered with porous elements, using eight classical RANS models and compared the results with direct numerical simulations (DNSs) data. Their findings revealed errors exceeding 20% in the friction coefficient predictions across all examined models, demonstrating that first- and second-order RANS models could yield incorrect trends.
Consequently, significant efforts have been devoted to developing turbulence models tailored for porous media, particularly in the turbulent regime. The prevailing approach employs the DANS equations with a two-equation
$\kappa$
-
$\varepsilon$
model closure, introducing an effective viscosity through the Boussinesq hypothesis (Masuoka & Takatsu, Reference Masuoka and Takatsu1996; Pedras & de Lemos, Reference Pedras and De Lemos2003; de Lemos, Reference De Lemos2005; Nakayama & Kuwahara, Reference Nakayama and Kuwahara2008). However, these models have two key limitations. First, their coefficients are typically calibrated for specific applications (e.g. periodic cells at moderate Reynolds numbers), restricting their general character; second, they rely on strong closure assumptions, such as isotropy, which can introduce inaccuracies and uncertainties, particularly in complex porous wall systems where anisotropy is dominant (de Lemos, Reference De Lemos2012). Developing a turbulence model requires balancing accuracy and efficiency, which demands a fundamental understanding of turbulence, supported by high-fidelity experiments or numerical simulations like pore-resolved DNS
Recent advances in computational resources and high performance computing have enabled more advanced and detailed numerical studies of simplified porous media, such as periodically arranged cylinders or spheres, in transitional and fully turbulent regimes (see the review by Wood 2020). These advancements provide an opportunity to better understand the physics of the flow inside porous media and refine turbulence models for improved accuracy and broader applicability (Jin & Kuznetsov, Reference Jin and Kuznetsov2024). Simulations such as DNS and large eddy simulation (LES) of porous medium flows are typically performed using the finite volume method (FVM) (Jin et al., Reference Jin, Uth, Kuznetsov and Herwig2015, Gasow et al., Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020, Reference Gasow, Kuznetsov, Avila and Jin2021, Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2022; Wang et al., Reference Wang, Chu, Lozano-Durán, Helmig and Weigand2021, Reference Wang, Lozano-Durán, Helmig and Chu2022; Srikanth et al., Reference Srikanth, Huang, Su and Kuznetsov2021) or the lattice Boltzmann method (LBM) (Suga, Reference Suga2016; Jin et al., Reference Jin and Kuznetsov2017; He et al., Reference He, Apte, Schneider and Kadoch2018, Reference He, Apte, Finn and Wood2019; Liu et al., Reference Liu, Sun, Wu, Wei and Hou2021, Reference Liu, Shi and Liu2023; Diao et al., Reference Diao, Chen, Liu, Wei and Hou2023). Finite difference and spectral methods are less commonly employed in this context. Among these techniques, LBM has gained popularity due to its parallelisability. However, further research is required to enhance LBM’s applicability, particularly regarding numerical stability and its extension to compressible flows (Fattahi et al., Reference Fattahi, Waluga, Wohlmuth, Rüde, Manhart and Helmig2016).
Commonly, DNS and LES employ body-fitted approaches that suffer from strong CPU limitation due to the no-slip boundary conditions (Jin et al., Reference Jin, Uth, Kuznetsov and Herwig2015). A key challenge in simulating porous media is the fine resolution required, which, along with body-fitted meshing constraints, increases computational costs. Mesh-related issues, such as poor-quality cells at contact regions, can be addressed with Cartesian grid-based methods like immersed boundary conditions (IBCs) (Mittal & Iaccarino, Reference Mittal and Iaccarino2005) and LBM.
In this context, a critical gap persists in understanding the interaction between turbulent flows and complex porous geometries, particularly in applications involving combustion. Addressing this challenge requires the development of advanced numerical techniques capable of accurately capturing the dynamics between turbulent flows and porous media. In the present work, we propose a high-fidelity numerical approach to investigate such interactions. A spectral difference (SD) method-based solver for the compressible NS equations, coupled with IBCs and pressure gradient scaling (PGS), is used to perform DNSs. These enhancements significantly improve both the accuracy and computational efficiency of the simulations, thereby enabling a deeper understanding of turbulence–porous medium interactions under realistic engineering conditions.
This study evaluates the robustness and reliability of the proposed numerical framework (Section 2), through benchmark tests described in Section 3. First, the implementation of IBCs within the SDM solver is compared with the classical body-fitted approach (Section 3.1). Second, the influence of the PGS method is assessed (Section 3.2). Afterwards, the full solver is benchmarked against the FVM-DNS data from Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) in Section 3.3, ensuring consistency with high-fidelity reference data and confirming its suitability for simulating turbulent flows in porous geometries. Finally, an application case on a complex porous medium geometry is provided in Section 4 and conclusions are drawn in Section 5.
2. Numerical methodology
This study uses the JAGUAR code (ONERA–CERFACS ©) to perform DNSs, relying on an SD method to solve the unsteady compressible NS equations on unstructured grids (Cassagne et al., Reference Cassagne, Boussuge, Villedieu, Puigt, D’Ast and Genot2015; Vanharen et al., Reference Vanharen, Puigt, Vasseur, Boussuge and Sagaut2017; Veilleux et al., Reference Veilleux, Puigt, Deniau and Daviller2022). The code employs a nodal, polynomial-based spectral representation of order
$p$
, in combination with advanced time integration techniques and Riemann solvers for flux reconstruction (Marchal et al., Reference Marchal, Deniau, Boussuge, Cuenot and Mercier2023). Regarding spatial discretisation, the solver uses Lagrange polynomials of order
$p$
to represent cell solutions and quadrature-based interpolation for flux computations. In this study
$p=4$
is adopted. Solution points are distributed according to Gauss–Chebyshev nodes of the first kind, yielding
$(p+1)^3$
points per cell, while flux points follow Gauss–Legendre quadrature, providing
$(p+2)^3$
points per cell, ensuring high-order accuracy. Diffusion terms are treated following the classical approach of Kopriva (Reference Kopriva1998), without
$p$
-refinement. Inter-element numerical fluxes are computed using the Harten–Lax–van Leer–contact Riemann solver (Toro, Reference Toro2009, pp. 315–344). Temporal discretisation is performed via the explicit five-stage, fourth-order total variation diminishing Runge–Kutta scheme of Spiteri and Ruuth (Reference Spiteri and Ruuth2002), which offers a favourable trade-off between computational efficiency and stability. A Courant–Friedrichs–Lewy (CFL) number of 0.15 is used in all simulations, and adaptive time stepping ensures numerical stability throughout the simulations. The PGS method is available (Section 2.1). Concerning the boundary conditions, IBCs are adopted (Section 2.2) and the possibility of injecting turbulence at the inlet is implemented via the synthetic random Fourier mode (SRFM) method based on a von Kármán-Pao (VKP) energy spectrum (Section 2.3).
2.1. Pressure gradient scaling
The method aims to improve the computational efficiency of low-Mach-number simulations performed with a compressible NS solver using explicit time integration schemes (Ramshaw et al., Reference Ramshaw, O’Rourke and Stein1985). It introduces a scaling coefficient
$\alpha \gt 1$
applied to the pressure
$p$
, such that
$p_{\textrm {PGS}} = p_{\textrm { noPGS}}/\alpha ^2$
, which reduces the effective speed of sound
$c = \sqrt {\gamma p / \rho }$
by a factor
$1/\alpha$
, where
$\gamma = C_{p}/C_{v}$
is the ratio of specific heats, with
$C_{p}$
the specific heat at constant pressure and
$C_{v}$
the specific heat at constant volume, and
$\rho$
is the fluid density. As a result, the CFL constraint on the time step,
$\Delta t = \text{CFL} \cdot \Delta x /(c + |\vec {u}|)$
, where
$\Delta x$
is the minimum mesh size,
$\vec {u}$
is the velocity vector and
$\text{CFL}$
is chosen to ensure stability without preconditioning, is relaxed. Hence, the maximum admissible time step increases by a factor
$\alpha$
, i.e.
$\Delta t_{\textrm {PGS}} = \alpha \, \Delta t_{\textrm {noPGS}}$
. Physically, this corresponds to reducing the sensitivity of the momentum equation to pressure gradients while preserving the velocity field, provided that pressure inhomogeneities remain small. It can also be interpreted as an `artificial increase of compressibility', since the effective Mach number
$Ma = |\vec {u}|/c$
is augmented by a factor
$\alpha$
due to the reduced sound speed.
Mathematically, let us consider a decomposition of the pressure field into a spatially averaged component and a fluctuation term,
$p(\vec {x},t) = \bar {p}(t) + p^\prime (\vec {x},t)$
. If the condition
$p^\prime / \bar {p} \ll 1$
holds, one may choose
$\alpha$
such that
$\alpha ^2 p^\prime /\bar {p} \ll 1$
. Under this assumption, the pressure gradient term in the momentum equation can be approximated as
$\nabla p \simeq \nabla (\hat {p}/\alpha ^2)$
, with
$\hat {p} = \bar {p} + \alpha ^2 p^\prime$
, where
$\hat {p}$
is referred to as the `modified pressure field' in the PGS framework. Therefore, the pressure-gradient-scaled compressible NS equations in conservative form read
\begin{equation} \left \{ \begin{array}{lll} \partial _t \rho + \vec {\nabla } \cdot (\rho \vec {u}) &= 0, \\[5pt] \partial _t (\rho {\vec {u} \,}) + \vec {\nabla } \cdot (\rho \vec {u} \otimes \vec {u}) &= - 1/\alpha ^2 \vec {\nabla } \hat {p} + \vec {\nabla } \cdot \overline {\overline {\tau }} + \vec {f}, \\[5pt] \partial _t (\rho E) + \vec {\nabla } \cdot \left ( \vec {u} (\rho E + \hat {p}) \right ) &= \vec {\nabla } \cdot \left ( {\vec {u} \, } \overline {\overline {\tau }} + \lambda \vec {\nabla } \hat {T} \right ), \end{array} \right . \end{equation}
where
$E = e + \frac {1}{2} |\vec {u}|^2$
is the total energy per unit mass,
$\overline {\overline {\tau }} = \mu (\vec {\nabla } \vec {u} + \vec {\nabla } \vec {u}^T) - \frac {2}{3} \mu \vec {\nabla } \cdot \vec {u} \overline {\overline {I}}$
is the shear stress tensor,
$\mu$
is the dynamic viscosity,
$\lambda = C_p \mu / \textit{Pr}$
is the thermal conductivity,
$C_p$
is the specific heat at constant pressure, Pr is the Prandtl number,
$\vec {f}$
represents the source term accounting for external or artificial forces and
$\hat {T} = f( {\hat {p}}, \rho )$
is the temperature, which follows the state equation based on the modified pressure
$\hat {p}(\mathbf{x}, t)$
. Note that the PGS equations (2.1) maintain the same formal structure as the classical compressible NS equations but with a modified pressure field.
2.2. Immersed boundary conditions
The IBCs are implemented in JAGUAR code (D’Ayer et al., Reference D’Ayer, Colombié, Dounia and Chedevergne2025) via the Brinkman penalisation method (Liu & Vasilyev, Reference Liu and Vasilyev2007), and integrated using the Strang splitting approach, as described by MacNamara and Strang (Reference MacNamara and Strang2016). This approach belongs to the family of volume penalisation (VP) methods (Arquis & Caltagirone, Reference Arquis and Caltagirone1984), and is thus classified as a continuous immersed boundary method (IBM). Hereafter, we will use the term IBC to refer specifically to the IBM implementation employed in this work.
A time-independent mask function
$\chi (\vec {x}) = \{ 1 \text{ if } \vec {x} \in \Omega , \ 0 \text{ if } \vec {x} \in \mathcal{F} \}$
distinguishes the fluid domain
$\mathcal{F}$
from the solid regions
$\Omega$
. A VP source term is added to the NS equations (2.1) everywhere in the domain to enforce the no-slip and no-penetration conditions at the fluid–solid interface
\begin{equation} \mathbf{S}(\vec {x},t) = -\chi (\vec {x})/\eta \text{ } \{\vec {0} ,\ \rho \vec {u_\Omega }(\vec {x},t) - \rho \vec {u}(\vec {x},t), \ \rho E_\Omega (\vec {x},t) - \rho E(\vec {x},t) \}, \end{equation}
where
$\vec {u}_\Omega (\vec {x},t)$
and
$E_\Omega (\vec {x},t)$
are the target velocity and total energy at position
$\vec {x}$
. In practice, the source term is applied only in the solid region (
$\chi = 1$
) and the target velocity is set to zero,
$\vec {u}_\Omega (\vec {x},t) = \vec {0}$
, so that the total energy reduces to the internal energy and a target temperature
$T_\Omega (\vec {x},t)$
can be introduced for a calorically perfect gas:
$E_\Omega = e_\Omega = c_v T_\Omega$
, where
$c_v$
is the specific heat at constant volume of the fluid. In this work, the walls are assumed to be isothermal with the inlet temperature so that
$T_\Omega = T_0 = 300\,\mathrm{K}$
. A crucial role in the continuous VP approach is played by the Brinkman penalisation coefficient
$\eta$
, which acts as a normalised viscous permeability (
$0 \lt \eta \ll 1$
). This coefficient controls the effectiveness of the method at the fluid–solid interface: decreasing
$\eta$
sharpens the diffuse representation of this interface where the velocity is very small but not zero and whose exact position is not known a priori. In fact, the error associated with this approach can be expressed as
$ \|\vec {u} - \vec {u}^{[n]}_\eta \| \le \|\vec {u} - \vec {u}_\eta \| + \|\vec {u}_\eta - \vec {u}^{[n]}_\eta \|,$
where
$\vec {u}$
is the exact solution of the NS equations,
$\vec {u}_\eta$
is the theoretical solution obtained using the VP-NS approach and
$\vec {u}^{[n]}_\eta$
is the corresponding numerical solution. It is known that
$\|\vec {u} - \vec {u}_\eta \| = \mathcal{O}(\eta ^{-1/2})$
decays exponentially with
$\eta$
(D’Ayer et al., 2025); thus, at the wall, we expect
$\vec {u}_\eta \sim \mathcal{O}(\eta ^{-1/2})$
. However, the discretisation error
$\|\vec {u}_\eta - \vec {u}^{[n]}_\eta \|$
may also depend on
$\eta$
. It is evident that small values of
$\eta$
are desirable, as they enforce the boundary conditions more accurately. However, decreasing
$\eta$
increases the system stiffness, and an explicit treatment of the penalisation term imposes a severe time-step restriction,
$\Delta t \le \eta$
, for numerical stability (Kolomenskiy & Schneider, Reference Kolomenskiy and Schneider2009). To circumvent this limitation, the penalisation term is treated separately from the other NS terms using Strang splitting (MacNamara & Strang, Reference MacNamara and Strang2016).
In the present computations,
$\eta$
is set to
$10^{-10}$
, several orders smaller than the time step, ensuring stability and effective boundary-condition enforcement. This yields normalised velocities in the diffuse interface of
$10^{-5}\!-\!10^{-7}$
, consistent with the theoretical prediction
$\vec {u}_\eta \sim \mathcal{O}(\eta ^{-1/2})=\mathcal{O}(10^{-5})$
.
Mesh resolution is also critical for accuracy. While refinement generally helps, the relative position of cells to the boundary matters: partially immersed cells can cause spurious near-wall fluctuations (Runge phenomena) since fixed-degree polynomials poorly approximate discontinuities (D’Ayer et al., Reference D’Ayer, Colombié, Dounia and Chedevergne2025, p. 6). Still, with
$\eta \ll 1$
, a sufficiently fine mesh gives results close to body-fitted simulations, as shown by D’Ayer et al. (Reference D’Ayer, Colombié, Dounia and Chedevergne2025) for cylinder and rough-channel flows.
2.3. Synthetic turbulence injection
The turbulent velocity field at the inlet,
$\vec {u}_{in} = \vec {\bar {u}}_{in} + \vec {u}_{in}^\prime$
, is modelled using Reynolds decomposition, where
$\vec {\bar {u}}_{in}$
is the time-averaged velocity and
$\vec {u}^\prime _{in}$
represents turbulent fluctuations. The mean velocity
$\vec {\bar {u}}_{in}$
is prescribed according to the desired profile, while
$\vec {u}^\prime _{in}$
is injected using the SRFM (Marchal et al., Reference Marchal, Deniau, Boussuge, Cuenot and Mercier2023). The method generates homogeneous isotropic turbulence by summing
$N$
Fourier modes with wavenumbers
$\kappa _n$
, following a user-defined energy spectrum. In this study, we use a VKP energy spectrum, detailed in Appendix C, with
$N = 300$
, turbulence intensity
$\text{Tu} = u^\prime _{rms} / \bar {u}_{in} = 10\%$
, energy-containing length scales
$L_{e,{\text{max}}} = 0.66$
m, and Kolmogorov length scales
$\eta _{k} = (\kappa _{Kol}^{-1}) = 6.33 \times 10^{-3}$
m. It should be emphasised that the applied boundary condition is not inherently periodic and therefore does not preserve perfect inlet periodicity under periodic boundary conditions.
3. Benchmark tests
3.1. The IBC validation
To validate the effectiveness of IBCs in modelling solid walls within complex geometries, we perform DNS of flow through a high-porosity periodic 2-D porous medium. The results obtained using IBCs are compared with those obtained with a traditional body-fitted meshing. The test geometry comprises square-edged obstacles, a particularly challenging case for IBCs due to the sharp fluid–solid interfaces.
Note that the compressible NS equations solved (2.1) include the energy equation. However, no heat transfer tests were performed for the configurations considered, since the present work focuses on velocity and turbulence statistics. Relevant validations of our IBC approach, also for heat transfer, are provided in D’Ayer et al. (Reference D’Ayer, Colombié, Dounia and Chedevergne2025) for a flow past a circular cylinder and a heated channel with academic roughness.
3.1.1. Numerical set-up
A turbulent flow past a squared cylinder (Figure 1 a) at Reynolds number
$\textit{Re}=u_m d/\nu = 1200$
is considered, where the kinematic viscosity is
$\nu = 1.33 \times 10^{-3}$
m2s−1, the square size is
$d=0.05$
m and the mean velocity component along the
$x_1$
-axis is
$u_m=32.0$
ms−1. The flow is driven by a source term
$\vec {f} = \{g_1,0,0\}$
in the NS equations (2.1), with
$g_1=1500$
Pa/m in both simulations, consistent with Section 3.3. The relatively high Reynolds number ensures conservative validation of the IBC approach. Periodic boundary conditions are imposed in all directions, yielding an effective porosity
$\varphi =0.93$
, with a computational domain (
$L_x \times L_y \times L_z$
) = (
$5d \times 3d \times 4d$
). The IBC case uses a Cartesian grid with (
$N_{x} \times N_{y} \times N_{z}$
) =
$(60 \times 36 \times 30)$
hexahedral cells and SD polynomial order
$p=4$
, giving approximately
$8.1$
million degrees of freedom (DoF), consistent with Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) for
$\textit{Re}=1000$
. The body-fitted case employs the same grid with the solid volume removed, yielding
$7.6$
million DoF. Average resolutions are
$(\Delta x_1/d)_{avg} = (\Delta x_2/d)_{avg} = 1.7 \times 10^{-2}$
and
$(\Delta x_3/d)_{avg} = 2.7 \times 10^{-2}$
.
Figure 1. Computational domain configurations and geometries used in the study.
Figure 2. The IBC validation: first-order velocity statistics
$\langle u_i \rangle$
on the averaged longitudinal plane (
$x_1$
,
$x_2$
).
Solver details are provided in Section 2. To isolate boundary-condition effects, no PGS method was used. Simulations ran on the IRENE Skylake Partition at TGCC with 500 processors for approximately
$25\,000$
CPU hours, advanced until statistical values such as
$\langle u \rangle$
and
$\langle uu \rangle$
stabilised. This required
$\sim$
10 wash-out cycles at
$p=4$
after initialising at lower resolution, where a wash-out cycle corresponds to the mean flow traversal time
$T_{\text{char}} = 5d/u_m \approx 1.35 \times 10^{-2}$
s. The time step from the CFL condition was
$\Delta t= 3.8\times 10^{-7}$
s in both cases.
3.1.2. Mean velocity field
In Figure 2, the 2-D distributions of the mean velocity field components
$\langle u_i \rangle$
over the averaged longitudinal plane
$(x_1, x_2)$
are presented. Since the volume average of
$\langle w \rangle$
is theoretically zero due to the 2-D configuration, this component is not shown. The streamwise velocity component,
$\langle u \rangle$
, exhibits excellent agreement, particularly in the critical square edge region (more details in Appendix A), where the boundary layer (BL) originates. The cross-wise component,
$\langle v \rangle$
, shows slight discrepancies, particularly in the impingement and wake regions. The differences in the impingement zone are expected, because the IBCs model the wall as a permeable medium, preventing the velocity from dropping abruptly to zero as it would under a no-slip condition. In the wake region, we attribute the observed discrepancies to the lower dissipation of IBCs compared with body-fitted approaches. Since IBCs represent the solid boundary as a low permeability medium, they capture weaker velocity gradients at the wall, which directly affects the computation of dissipation. As a result, IBCs tend to slightly underestimate the velocity magnitude, as they do not fully account for all dissipation mechanisms at the wall. The majority of the observed error occurs in the wake region due to a combination of slightly different upstream conditions that can amplify the upstream–downstream wake interactions in a periodic configuration.
To quantify these effects, Figure 3 shows velocity profiles along the
$x_1$
and
$x_2$
-directions, confirming nearly perfect agreement for
$\langle u \rangle$
and only minor, localised discrepancies in
$\langle v \rangle$
. Figure 4 presents the Reynolds stress tensor (RST) components
$\langle u^\prime _i u^\prime _i \rangle$
on the
$(x_1,x_2)$
plane, showing excellent agreement in both shape and magnitude, even near walls. Table 1 summarises volume-averaged errors for first- and second-order velocity statistics, with all
$L_2$
and
$L_\infty$
norms below
$1$
% and
$5.5$
%, respectively.
Figure 3. The IBC validation:
$\langle u_i \rangle$
on lines along
$x_1$
and
$x_2$
from the averaged longitudinal plane (
$x_1$
,
$x_2$
). For visualisation, IBC results are shown at reduced resolution; however, both body-fitted results and the error between the methods are presented at full simulation resolution.
Figure 4. The IBC validation: second-order velocity statistics
$\langle u^\prime _i u^\prime _i \rangle$
on the averaged longitudinal plane (
$x_1$
,
$x_2$
).
3.2. The PGS validation
While the efficacy of PGS has been demonstrated by various authors for FVM solvers (Papageorgakis & Assanis, Reference Papageorgakis and Assanis1999), its performance within an SD solver using IBCs for complex geometries has yet to be thoroughly assessed. The objective of this validation is, therefore, to compare two simulations of a porous medium using both PGS and non-PGS approaches. To ensure a fair comparison, both simulations are conducted over the same physical time span. However, the non-PGS case requires significantly more iterations. Since PGS modifies the flow dynamics by artificially adjusting the pressure field, allowing for a relaxed time step via an increased Mach number, the critical concern is whether this method alters the mean flow or fails to capture the full spectrum of temporal frequencies required for statistical analysis.
Table 1. The IBC validation: volume-averaged errors using different norms (
$L_1$
,
$L_2$
and
$L_\infty$
) for first- and second-order velocity statistics. Norms are normalised by the bulk velocity
$u_m$
of the body-fitted case
3.2.1. Numerical set-up
We consider a simplified porous medium configuration as shown in Figure 1b. This set-up models a turbulent flow through a single pore of a Gyroid-type porous topology, belonging to the triply periodic minimal surface (TPMS) family (Wang et al., Reference Wang, Chen, Zeng, Ma, Wang and Cheng2023), whose solid structure is defined by positive values of the function:
$ f(x,y,z) = \sin (k_x x) \cos (k_y y) + \sin (k_y y) \cos (k_z z) + \sin (k_z z) \cos (k_x x) - \sigma _t$
. In this case the structural factor
$\sigma _t = 0.75$
, meaning a porosity
$\varphi = 0.7$
, with wavenumber coefficients
$k_x = k_y = k_z = 6 \pi$
$\mathrm{m}^{-1}$
. The computational domain is a hexahedral box
$(L_x \times L_y \times L_z) = (4s \times s \times s)$
, where a cubic pore cell of dimension
$s=1/3$
m is located within the region
$x \in [s, 2s]$
. Periodic boundary conditions are applied in the directions transverse to the flow. A turbulent flow is injected at the inlet along the flow direction with a unity mean axial velocity
$u_{0}=1$
ms−1, a Reynolds number
$\textit{Re}=(3s) u_{0}/\nu = 2500$
and a turbulent VKP spectrum, described in Section 2.3. A sponge zone is implemented just before the outlet (
$x \in [3s, 4s]$
) to prevent acoustic wave reflections. A pressure outlet condition is imposed with a reference pressure
$p_0 = 10^5$
Pa. Both simulations employ IBCs on a Cartesian grid with
$(N_{x} \times N_{y} \times N_{z}) = (80 \times 20 \times 20)$
hexahedral cells, using a SD polynomial order
$p = 4$
, resulting in approximately 4 million DoFs. This grid resolution is selected to ensure a ratio of
$\eta _k / \Delta x_{\text{avg}} = 1.5$
(Pope, Reference Pope2000, pp. 344–351), where
$\eta _k$
is the injected Kolmogorov length scale.
Two simulations were performed: one with PGS (
$1/\alpha ^2 = 0.05$
) and one without. The average time steps were
$6.2 \times 10^{-6}$
and
$1.4 \times 10^{-6}$
s, respectively. Low-pressure inhomogeneity assumptions were confirmed, with negligible overall pressure gradients. Simulations ran on the Skylake partition of the IRENE supercomputer at TGCC using 512 processors, requiring 30 000 CPU hours with PGS and 125 000 hours without. They continued until statistical quantities, such as time correlations, reached steady state – approximately 20 wash-out cycles, after an initial ramp-up of SD polynomial order and PGS influence. A representative period of over 10 wash-out cycles (
$T_{\text{wash-out}} = (L_x - L_{x,\text{sponge}})/u_0 = 1$
s) was selected for robust statistics. The PGS simulation also served to initialise the no-PGS case, providing approximate initial conditions for high-fidelity computations.
It is important to note that the time step in our simulations is constrained by the combined effect of the CFL condition and the high spatial resolution required for accurately modelling complex geometries, rather than by the physical time scales of turbulence. Specifically, the estimated Kolmogorov time scale
$\tau _\eta$
in homogeneous isotropic turbulence (HIT) conditions (Pope, Reference Pope2000, p. 185) of the injected turbulent spectrum is of the order of
$10^{-2}$
$s$
. In contrast, the simulation time steps are of the order of
$ 1.0 \times 10^{-6}$
s (for a
$ \text{CFL} = 0.15$
), nearly four orders of magnitude smaller. As a result, the simulations resolve extremely high frequencies that are not only unnecessary for our analysis but also computationally expensive. This highlights the importance of validating the PGS method.
3.2.2. Mean velocity field
In Figure 5, the mean velocity components
$\langle u_i (\vec {x}_0, t) \rangle$
in the streamwise longitudinal pore-crossing section
$(x, y)$
are analysed. Good agreement is observed between the PGS and no-PGS approaches in both the shape and magnitude of the velocity field. Local minima and maxima are also reasonably well predicted, and detachment phenomena and wake structures remain consistent. To better quantify the differences, velocity profiles along the
$x$
- and
$y$
-axes are extracted from Figure 5 and shown in Figure 6. The PGS approach effectively captures the mean flow features, with only minor accuracy losses in separation and shear regions, where local relative errors may reach up to
$10\%$
. These discrepancies likely arise from a slight shift in the inlet–outlet pressure drop induced by the PGS, which appears to cause a mild suction effect at the inlet, slightly advancing wake development. Notably, PGS performs relatively well near the walls, suggesting that any potential PGS–IBC interaction does not significantly amplify the errors introduced by IBCs.
Figure 5. The PGS validation: mean velocity components
$\langle u_i \rangle /u_0$
on the longitudinal pore-crossing plane.
Figure 6. The PGS validation:
$\langle u_i \rangle$
on lines along
$x$
and
$y$
extracted from the longitudinal pore-crossing plane. For visualisation, PGS results are shown at reduced resolution; however, both no PGS results and the error between the methods are presented at full simulation resolution.
Table 2 presents volume-averaged errors for both first and second-order velocity statistics using different norms to compare PGS and no-PGS simulations. The mean relative errors in norm
$L_2$
are all below
$4.5\%$
and
$3.5\%$
for the velocity components
$\langle u_i \rangle$
and RST components
$\langle u^\prime _i u^\prime _j \rangle$
, respectively. The largest local errors are concentrated downstream of the pore, in the wake region, as illustrated in Figure 6. This allows us to say that the slight difference in pressure field affects only the wake area, suggesting that the inside pore physics is not contaminated by the PGS method. This is crucial for applications where we aim to simulate multiple pores in the streamwise direction.
Table 2. The PGS validation: volume-averaged errors using different norms (
$L_1$
,
$L_2$
and
$L_\infty$
) for first- and second-order velocity field statistics, averaged over the computational domain (excluding the sponge zone and the very close region to it)
3.2.3. Energy spectra
Since PGS primarily affects time resolution, it is natural to compare temporal signals at the same location with and without PGS. Here, energy spectra
$E_{ii} (\vec {x_0}, f)$
, where
$f$
is the temporal frequency, are analysed at two spatial probes: one upstream (A) of the porous medium in
$\vec {x}_A = (0.1, 0.0, 0.0)$
and the other downstream (B) in
$\vec {x}_B = (0.8, 0.0, 0.0)$
. The energy spectra definition is the one used by Pope (Reference Pope2000, pp. 65–73) as the fast fourier transform (FFT) of the (one-point space) two-time covariance
$R_{ii}(\vec {x}_0,t') = \langle u^\prime _i(\vec {x}_0, t) u^\prime _i(\vec {x}_0, t+t') \rangle _T$
. In Figure 7, energy spectra of the RST normal components
$E_{ii}$
are presented and show strong coherence in energy distribution and overall captured energy.
In the PGS case, the larger time step (lower time resolution) leads to a slight overestimation of energy at the inlet. However, the energy cascade remains consistent. At very high frequencies, numerical noise appears at the inlet, but these non-physical scales are negligible since their energy content is several orders of magnitude lower than the injected Kolmogorov-scale frequency (
$f_\eta = 1/\tau _\eta \simeq 12.5$
Hz). At the upstream probe, a peak at
$f = 3.0$
Hz corresponds to the injected turbulence, where the energy-containing eddies have a characteristic size
$l_0 = 0.33$
m and velocity
$u_0 = 1$
m/s, yielding a characteristic time
$\tau _0 = l_0/u_0 = 0.33$
s and frequency
$f = 1/\tau _0 = 3.0$
Hz. Downstream, two key observations arise. First, the energy spectrum is enriched due to turbulence production by the pore, increasing both kinetic energy and the turbulence-scale range at high frequencies. Second, a peak at
$f = 5.5$
Hz corresponds to vortex shedding behind the pore, behaving as a blunt body. This frequency aligns with a turbulent Strouhal number
$St = fL/U = 0.25$
, considering an approximated obstacle diameter of
$L \simeq s/10$
and a mean axial velocity at the pore outlet
$U \simeq 1.8$
m/s.
Overall, the robustness of the PGS method is confirmed by the minimal differences in captured energy (integral of the spectrum), which remain below
$10\%$
at the downstream probe, the most critical location for this case.
Figure 7. The PGS validation: comparison of the energy spectra
$E_{ii}(\vec {x}, f)$
at the inlet probe location
$\vec {x}_A$
(green) and the outlet probe location
$\vec {x}_B$
(red). Solid coloured lines (—) correspond to simulations without PGS; dashed coloured lines (- - -) correspond to simulations with PGS; dotted black line (
$\cdot \cdot \cdot$
) represents the Kolmogorov
$-5/3$
power law.
3.3. Porous medium reference case
Following the demonstration of the influence of IBCs and PGS method, we now proceed to a comparison of the complete numerical approach, integrating both PGS and IBCs to enhance the SD solver’s performance, with a classical FVM solver, which have been extensively employed in the literature. For this comparative analysis, the reference work of Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015), who conducted high-fidelity DNS using both FVM and LBM on a simplified porous medium consisting of an array of squared cylinders, is reproduced.
As those authors demonstrated comparable performances between their numerical approaches, the present analysis will compare our SD solver to the FVM reference case. This choice is also motivated by the FVM’s widespread adoption within the porous medium turbulence literature and its generally perceived higher reliability for such complex flows.
The following comparison between SD-DNS and FVM-DNS focuses primarily on two-point velocity correlations, as we are interested in turbulent flow applications.
3.3.1. Numerical set-up
The geometry considered is the one used by Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015), a generic porous medium composed of squared rods of size
$d$
arranged periodically to form a porous structure with pore size
$s$
(Figure 1c). This geometry has been also used by other authors in their numerical studies (Kuwahara et al., Reference Kuwahara, Yamane and Nakayama2006; Suga, Reference Suga2016; Chu et al., Reference Chu, Weigand and Vaikuntanathan2018). Specifically, we consider the configuration where
$s/d=2$
(meaning a porosity
$\varphi = 0.88$
) and
$Re=700$
, which corresponds to case E of the high-resolution FVM-DNS simulations presented by Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015).
The size of the computational domain was determined using the concept of REV-T, which defines the representative elementary minimum volume necessary to capture turbulent behaviour without loss of information. Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) conducted a REV-T size study, systematically reducing the number of bars in the computational domain with periodic boundary conditions until a minimum size was identified that preserved large-scale turbulent structures. Based on their findings, the computational domain is a box of size (
$L_{x} \times L_{y} \times L_{z}$
) = (
$12d \times 8d \times 4d$
), for a total of 12 squared cylinders of size
$d \times d \times 4d$
, with
$d=0.05$
m.
Periodic boundary conditions are imposed on all boundaries of the computational domain in all three spatial directions. The turbulent flow analysed is self-generated by the porous medium (Jin et al., Reference Jin, Uth, Kuznetsov and Herwig2015) and is imposed by adding a pressure gradient source term in the
$x_1$
-direction, as detailed previously in Section 3.1. In this way, a Reynolds number
$Re=u_m d/\nu \simeq 700$
is reached, where
$u_m = 18.40$
ms−1 is the bulk velocity along the
$x_1$
-axis. However, we underline that Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) concluded in their study that the Reynolds number (in the considered range of
$500{-}1000$
) has a negligible effect on the dimensionless two-point correlations
$R_{ii}/u_m^2$
for the porosity values being examined.
Moreover, Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) also performed a mesh convergence study by introducing an accuracy parameter based on dissipation calculations and pressure loss imposed through the pressure gradient source term. Following their approach, we adopted the same mesh density. Using IBCs, this resulted in a Cartesian grid with (
$N_{x} \times N_{y} \times N_{z}$
) =
$(144 \times 96 \times 30)$
cells and a SD polynomial order of
$p=4$
, yielding approximately
$51.8$
million DoF.
To enhance simulation efficiency, PGS method with a relatively mild value
$\alpha ^{-2}=0.7$
is used, allowing for a 20% increase in the time step compared with the case without PGS. We exercise caution when applying PGS in this case because the presence of sharp geometric features can locally induce high Mach numbers, and artificially increasing them through PGS could lead to undesirable effects. However, for smoother geometries, as demonstrated in Section 2.1, more aggressive PGS can be safely employed.
The simulations were conducted on the IRENE Skylake Partition at TGCC using 2500 processors. The total CPU time for the simulation was approximately
$150\,000$
h. The simulation was run until statistical values, such as spatial correlations, exhibited clear convergence, with time statistics remaining stable. This process required up to
$200$
wash-out cycles in one pore element. For validation purposes, we do not extend the simulation further to achieve perfect convergence of the two-point correlations, as this process is computationally expensive, especially for the
$x_3$
-direction, where convergence is slow. However, we ensure that the observed convergence was sufficient and that the simulated time span of
$200$
wash-out cycles was comparable to the 250 wash-out cycles performed by Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015), particularly for comparing two-point correlation behaviour. Note that one wash-out cycle is computed considering the characteristic time
$T_{\text{char}} = 4d/u_m \simeq 1.1 \times 10^{-2}$
s. The time steps for
$p=4$
are around
$3.8\times 10^{-7}$
s.
In contrast, Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) employed a FVM solver that directly solves the incompressible NS equations. Time integration was performed using a second-order implicit backward Euler scheme, while spatial discretisation relied on a second-order central difference method. Pressure–velocity coupling was managed through the pressure-implicit with splitting of operators (PISO) algorithm, and periodic boundary conditions were applied. The simulation used approximately
$51.8$
million DoFs, with mesh refinement near the walls, and required approximately
$50,000$
CPU hours (using 480 processors) for the test case considered.
3.3.2. Two-point correlations
One effective approach to detecting turbulent structures and analysing their scales is by determining two-point correlations in the velocity field, as the latter fully characterise an incompressible turbulent flow.
The two-point correlation function between the velocity fluctuations
$u'_i (\vec {x_0},t)$
and
$u'_j (\vec {x_0} + \vec {r},t)$
at a given time
$t$
is referred to as the two-point (one-time) autocovariance by Pope (Reference Pope2000, pp. 77−78):
$R_{ij} (\vec {r}, \vec {x_0}) = \langle u'_i (\vec {x_0}, t) u'_j (\vec {x_0} + \vec {r}, t) \rangle$
. This function characterises how the velocity fluctuations at
$\vec {x_0}$
correlate with those at neighbouring points at a distance
$\vec {r}$
. The peak of this function occurs at
$\vec {r}=\vec {0}$
, indicating that the strongest correlation is self-referential. Usually, two-point correlation functions are adimensionalised by their peak value
$R_{ij}(\vec {r}=\vec {0})$
but here we use instead the bulk velocity squared, as done also by Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) to remove the influence of the Reynolds number.
Figure 8 presents the 2-D distributions of the normalised two-point correlation functions
$R_{11}(\vec {x_0}, \vec {r})/u_m^2$
and
$\hat {R}_{11}(\vec {x_0}, \vec {r})/u_m^2$
in the longitudinal midplane. The results closely match the FVM-DNS data of Jin et al., (Reference Jin, Uth, Kuznetsov and Herwig2015, Figures 4 and 5) in both magnitude and structure. This agreement further confirms the ability of the implemented IBCs to accurately resolve boundary-layer features and the wake dynamics behind the cylinders.
In Figure 9 normalised two-point correlation profiles along the
$x_1$
- and
$x_2$
-directions are extracted from the planar data of Figure 8. The SD solver accurately predicts both the complex functional behaviour and the peak magnitude of the correlation functions. The same agreement is also found for the turbulent correlation
$\hat {R}_{11}/u_m^2$
(Jin et al., Reference Jin, Uth, Kuznetsov and Herwig2015, pp. 84−85), not presented here. This implies that both the turbulent and non-turbulent correlation components are faithfully solved. Minor oscillations indicative of incomplete convergence persist, particularly in the
$x_3$
-direction. Despite these marginal convergence artefacts, the overall agreement is substantial, with the most significant discrepancy being an approximate 5%–10% underestimation of the
$\hat {R}_{11}/u_m^2$
peak. We posit that part of this underestimation would be considerably reduced upon achieving complete convergence of the simulation, as shown by the 90% confidence intervals (CIs). As previously noted, although a trend toward convergence along the
$x_3$
-direction was observed, the simulation was stopped due to the high computational cost and time required to reach full convergence, which was considered beyond the primary scope of this validation.
Figure 8. Two-point correlation
$R_{11}(\vec {x_0}, \vec {r})/u_m^2$
(top) and turbulent two-point correlation
$\hat {R}_{11}(\vec {x_0}, \vec {r})/u_m^2$
(bottom) in the plane
$x_3 = L_3/2$
at the correlation point
$\vec {x_0}=(3s,2s,s)$
, marked with `
$+$
'.
Figure 9. Two-point correlations
$R_{ij}(\vec {x_0}, \vec {r})/u_m^2$
at the correlation point
$\vec {x_0}=(3s,2s,s)$
along
$x_1$
-axis (left) and
$x_2$
-axis (right). — (blue) indicates the SD results with a 90% CI region highlighted; - - - (red) indicates FVM results of Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015).
4. Application case
In this section, the whole numerical approach is applied to a complex porous medium (Figure 1d). A turbulent flow through a TPMS porous geometry of the Gyroid type is considered, the geometry being defined using the equation defined in Section 3.2 with a structural factor
$\sigma _t = 1.20$
, corresponding to a porosity
$\varphi = 0.90$
and
$k_x = k_y = k_z = 6\pi$
$\mathrm{m}^{-1}$
. The computational domain is a box
$(L_x \times L_y \times L_z) = (8L \times L \times L)$
, where
$L$
is unity. This configuration extends the elementary pore simulation analysed in Section 3.2, using the same boundary conditions set-up and the same inlet Reynolds number
$\textit{Re} = Lu_0/\nu = 2500$
. The porous matrix is embedded in the region
$x \in [2L, 6L]$
and it accounts for
$12 \times 3 \times 3$
pores in the
$x$
,
$y$
and
$z$
directions, respectively, with a pore cell size
$s=L/3$
. A sponge zone is implemented in the region
$x \in [7.5L, 8L]$
, sufficiently far from the porous matrix to avoid interaction with the wakes immediately downstream.
The simulation employs IBCs on a Cartesian grid of
$(N_x \times N_y \times N_z) = (480 \times 60 \times 60)$
hexahedral cells with SD polynomial order
$p=4$
, yielding
$\sim$
216 million DoF (
$\sim$
1 million per pore cell). The resolution satisfies
$\eta _\kappa / \Delta x_{\text{avg}} = 1.5$
(Pope, Reference Pope2000), with
$\eta _\kappa = 6.33 \times 10^{-3}$
m. Mesh refinement confirmed convergence, and a domain size study (configurations
$3 \times 3$
vs.
$5 \times 5$
pores) showed no significant differences in two-point correlations, so the
$3 \times 3$
set-up was adopted.
Figure 10. Vorticity magnitude normalised by the local velocity magnitude
$\|\vec {\omega }\|/\|\vec {u}\|$
.
The PGS method uses
$1/\alpha ^2 = 0.05$
, allowing a much
$\alpha \simeq 4.5$
times larger time step than the baseline (
$1.95 \times 10^{-6}$
s). Runs were performed on the IRENE Skylake partition with 3500 processors at a cost of
$\sim$
2 million CPU hours. Statistics stabilised after approximately 60 wash-out cycles (
$T_{\text{wash-out}} = s/u_0 = 0.33$
s). A representative window of
$\sim$
100 wash-out cycles was then used to collect robust turbulence data.
4.1. Preliminary results
Our primary interest lies in the investigation of turbulent flow behaviour in porous media, particularly in the context of turbulent combustion applications. A key focus is understanding the evolution of turbulent length scales as the flow develops within the porous matrix. Figure 10 presents the vorticity magnitude normalised by the local velocity magnitude, visualised on a longitudinal plane. This highlights the characteristic eddy sizes along the streamwise direction. At the inlet, the imposed turbulence contains energy-bearing eddies with characteristic dimensions approximately twice the pore cell size. As the flow enters the porous structure, these eddies undergo a progressive scale reduction, adapting to the pore size over several layers in the longitudinal direction. We refer to this region as the turbulence length scale adaptation region, which is of practical importance. In combustion applications, this region often represents a critical zone where flame stabilisation depends on precise control of turbulence properties.
Figure 11. Averaged mean fields along the longitudinal axis in the
$x_1$
-direction: velocity components
$\langle u_i \rangle$
, normal RST components
$\langle u_i^\prime u_i^\prime \rangle$
and integral time scales
$\tau _i = \int _0^{+\infty } \rho _{ii}(\vec {x_0}, t') \, dt'$
, where
$\rho _{ii}(\vec {x_0}, t')$
is the temporal autocorrelation function of the velocity component
$u_i$
at the spatial location
$\vec {x_0}$
for a time shift
$t'$
(Pope, Reference Pope2000, pp. 65−73). Empty circles (’o’) represent full transversal section averages while lines (’—’) are averages along the centre line of each pore.
In the specific configuration considered here, defined by the pore scale, porosity and inlet flow conditions, the adaptation region spans approximately three pores. Beyond this region, turbulence is fully adapted to the pore scale, in agreement with the pore-scale prevalence hypothesis (Uth et al., Reference Uth, Jin, Kuznetsov and Herwig2016), which posits that, in a porous medium, the size of turbulent structures is constrained by the pore scale. In this fully adapted region, the flow can be considered statistically developed. At the outlet, as the flow exits the porous matrix, the eddies begin to merge once more, in a manner reminiscent of decaying grid-generated turbulence.
An illustrative example of physical outcomes that can be extracted from this simulation is given in Figure 11. However, results analysis and further physical discussions are beyond the scope of this paper.
5. Conclusions
A numerical approach has been developed to study turbulent flows through complex porous media, achieving a good compromise between accuracy and computational cost. The method relies on an SD solver, the JAGUAR code (ONERA–CERFACS ©), combining IBCs with the PGS approach to efficiently tackle complex geometries in low-Mach-number conditions. The approach has been validated using an SD polynomial order
$p=4$
through a series of three benchmark tests on simplified porous configurations.
First, for IBC validation, SD simulations using IBCs were compared with those based on classical body-fitted meshes. The main challenge with IBCs lies in accurately capturing near-wall phenomena such as BLs and stagnation or recirculation zones. Our results show that the average relative error in the
$L_2$
norm for first- and second-order mean velocity fields remains below 1%, with local discrepancies, primarily near walls and in wake regions, limited to 5.5% in the
$L_\infty$
norm. These differences do not alter the overall mean flow behaviour, confirming that the IBC implementation faithfully reproduces the mean dynamics while significantly reducing meshing complexity around intricate geometries.
Second, for PGS validation, SD simulations with and without the PGS method were compared. The selected PGS factor increases the effective time step by a factor of
$\alpha = 4.5$
. While the PGS method introduces slight modifications to the pressure field under incompressibility assumptions and filters out high-frequency content, we observed good agreement in both mean flow statistics and energy spectra. The average
$L_2$
norm errors for first- and second-order velocity moments are below 4.5% and 3.5%, respectively. The energy spectra correctly capture the relevant peaks and magnitudes of dominant frequencies, with very limited relative errors, mainly concentrated in the wake downstream of the porous section. These results confirm that the PGS method preserves the essential flow physics while significantly accelerating simulations.
Third, for the porous medium reference case, the SD-DNS approach, employing both PGS and IBCs, was compared with the classical porous medium simulation by Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015) using FVM-DNS. The SD solver reproduces the normalised turbulent two-point correlations
$R_{ii}/u_m^2$
with high accuracy. A slight underestimation of the peak values is observed, attributed to reduced near-wall dissipation introduced by the IBCs, which leads to marginally lower fluctuation levels.
Finally, we applied the developed approach to a complex porous geometry involving approximately a quarter of a billion degrees of freedom. The configuration consists of a TPMS porous medium of Gyroid type, with porosity
$\varphi = 0.90$
, traversed by an injected turbulent flow at an inlet Reynolds number of
$\textit{Re} = 2500$
, turbulence intensity of
$\text{Tu}=10\%$
, and energy-containing turbulent scales twice the pore cell size. This case introduces the novel consideration of injected turbulence as a boundary condition, representative of realistic engineering applications.
Overall, this work presents an alternative high-fidelity numerical strategy to traditional FVM-DNS and LBM-DNS approaches for the study of complex porous flows. In future work, we aim to apply this methodology to TPMS porous media to investigate how pore shape, porosity and injection parameters influence the interaction between turbulence and porous structures.
Supplementary material
Animations illustrating the application case, described in Section 4, are available online at the following link: https://www.youtube.com/@porous\media\channel. The supplementary material for this article can be found at https://doi.org/10.1017/flo.2025.10036.
Data and coding availability statement
Raw data and post-treatment codes are available from the corresponding author.
Author contributions
Adrian Rusnak: Investigation, Formal analysis, Data curation, Visualisation, Writing – Original draft. François Chedevergne: Conceptualisation, Methodology, Supervision, Writing – Review & Editing. Rémi Roncen: Conceptualisation, Methodology, Supervision, Project administration, Funding acquisition, Writing – Review & Editing.
Funding statement
The authors thank H. Deniau, T. Marchal, E. D’Ayer and A. Colombié for their support with the JAGUAR simulations and post-processing. Co-funded by the European Union (ERC, POROLEAF, 101103502). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Computational resources were provided by GENCI at TGCC under the grants 2024-SS012A15391 and 2025-A0172B15655 on the supercomputer IRENE’s SKL Partition.
Declaration of generative AI and AI-assisted technologies in the writing process
During the writing of the original draft, A. Rusnak used ChatGPT in order to improve language and readability. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.
Declaration of interests
The authors declare no conflict of interest.
Appendix A. The IBC impact: near-wall flow and pressure drop
Using IBCs naturally introduces near-wall errors, which are intrinsic to the artificial modelling. In Section 3.1, we showed that IBC results closely match classical body-fitted simulations for both first- and second-order velocity statistics. Here, we focus on the near-wall region to illustrate the diffused fluid–solid interface and the influence of VP on the pressure gradient.
Figure A1 shows the velocity streamlines at the edge
$(x,y)=(0.1,0.1)$
of the solid bar. The corresponding streamwise velocity profile along
$x=1$
is shown in Figure A2. A thin diffused layer of size
$h/d \sim 10^{-3}$
, characteristic of continuous IBM approaches, appears with normalised BL velocities of
$10^{-5}$
–
$10^{-7}$
, while the velocity inside the solid remains effectively zero. The sharp solid edge is preserved, although a slight curvature proportional to
$h/d$
smooths the boundary locally.
Figure A1. The IBC validation (case with IBCs): streamlines around the upstream corner of the solid. Red lines indicate the exact solid edges.
Figure A2. The IBC validation (case with IBCs): cross-wise profile of the streamwise velocity at the solid edge (x = 1 in Figure A1). The diffused interface layer due to the continuous IBC is highlighted.
Figure A3. The IBC validation: pressure drop measured between upstream and downstream probes (top) and corresponding frequency spectrum (bottom).
Figure A3 reports the pressure drop
$\Delta p = p_{\mathrm{out}} - p_{\mathrm{in}}$
measured by probes upstream and downstream of the obstacle. The relative difference between IBC and body-fitted cases is approximately
$+0.45\%$
for the mean and
$9.0\%$
for the variance (Figure A3, top). The corresponding pressure-drop spectrum (Figure A3, bottom) shows only marginal spectral discrepancies.
Appendix B. The PGS impact: pressure drop
The main requirement of the PGS method is the low-inhomogeneity assumption. In our study, pressure variations remained well within this condition (
$p^\prime / \bar {p} \ll 1\%$
) across most of the domain. Slight deviations up to 2% were observed near the walls due to the IBCs, but these are very localised and negligible in magnitude. Moreover, the averaged pressure drop in the streamwise direction was almost unchanged after introducing the PGS, as illustrated in Figure B1.
Figure B1. The PGS validation: averaged pressure drop (
$p-p_0$
) along the streamwise direction
$x$
, where
$p_0$
is the inlet pressure.
Appendix C. Von Kármán-Pao (VKP) spectrum definition
The VKP energy spectrum
$E(\kappa )$
in wavenumber space
$\kappa$
as defined by Pao (Reference Pao1965) as
$E(\kappa ) = \beta u^{\prime 2}_{rms}/\kappa _{e} \left (\kappa /\kappa _{e}\right )^4 \left (1 + (\kappa /\kappa _{e}\right )^2)^{-17/6} \exp (-2\left (\kappa /\kappa _{\text{Kol}}\right )^2)$
, where
$\beta$
is a scaling coefficient,
$\kappa _{e}$
corresponds to the VKP spectrum energy-containing wavenumber and
$\kappa _{\text{Kol}}$
denotes the Kolmogorov wavenumber. Assuming HIT in the limit of infinite Reynolds number, meaning
$\left (\kappa _{\text{Kol}} / \kappa _{e}\right )^2 \to \infty$
, the integral of the spectrum can be enforced to match the turbulent kinetic energy as
$ \int E(\kappa ) \, d\kappa = \text{TKE} = 3/2 u^{\prime 2}_{rms}$
which leads to the approximate value
$\beta \approx 1.45$
. Note that the maximum of the VKP spectrum occurs at the wavenumber
$\kappa _{e_{max}} = \sqrt {12/5} \, \kappa _{e}$
. Thus the other parameters of the isotropic VKP spectrum are given by
$ L_{e} = \sqrt {12/5} \pi /(2u^{\prime 2}_{rms}) \int E(\kappa )/\kappa \, d\kappa$
and
$ \kappa _{\text{Kol}} = \left (\epsilon /\nu ^3\right )^{1/4} = \left [ (2\nu \int \kappa ^2 E(\kappa ) \, d\kappa )/\nu ^3\right ]^{1/4}$
. Therefore, the spectrum can be fully characterised by
$u^{\prime }_{rms}$
and
$\kappa _{e}$
. However, in practical applications, where the Reynolds number is finite (albeit large), the VKP spectrum limits
$[\kappa _{e},\kappa _{\text{Kol}}]$
have to be explicitly defined as HIT assumptions are no longer valid. The strategy adopted in this case consists in constructing a spectrum that approximates as closely as possible the one under HIT conditions. We prescribe the energy-containing eddy length scale
$L_{e,\text{VKP}}$
, the Kolmogorov length scale
$\kappa _{Kol}$
and the turbulent intensity rate
$Tu = u^{\prime }_{\textit{rms}} / u_0$
, such that a target root-mean-square turbulence intensity
$u^\prime _{rms}=\text{Tu} \cdot u_0$
at the inlet Reynolds number
$Re=u_0L/\nu$
is achieved. The amplitude
$\beta$
of the spectrum is then obtained by solving an optimisation problem of the type
$\min _{\beta } \left ( \epsilon _L + \epsilon _k \right )$
aimed at minimising the relative errors on the TKE
$\epsilon _k=\left | \int E(\kappa ) \, d\kappa - 3/2 u^{\prime 2}_{rms} \right | / u^{\prime 2}_{rms}$
and the integral length scale
$\epsilon _L=\left | L_{e_{\text{max}}} - \pi /(2 u^{\prime 2}_{rms}) \int E(\kappa )/\kappa \, d\kappa \right | / L_{e_{\text{max}}}$
if compared with the HIT solutions.
Open access
