A Weakly Compressible SPH Method for the Brinkman Equation and Its Application to Estimating the Bulk Permeability of Digital Cores

1. Introduction

In-depth understanding of fluid flow in multiscale porous media is always required in engineering including but not limited to oil and gas recovery [39,40], carbon sequestration [32,33,34,37,38] and geothermal production [31,35,36]. Porous media in applications of these fields are considered as multiscale because the pores, throats and fractures in them usually span of a wide range in scale, from nanometers to micrometers and from several micrometers to millimeters. Numerical simulations are powerful tools to investigate the characteristic of the flows in these multiscale porous media. However, the classical pore scale simulations by solving the Navier-Stokes equation need too much detailed geometry of pores in all scales of the porous media of interest. Therefore, they require tremendous computational resources to implement the simulation and is impractical in application. A more realistic method is to model the porous media of interest with a properly chosen scale, which depicts the main structure of the porous media with resolved big pores and solid walls on one hand, and hides minor pore-solid structures below that scale in unresolved pixel or voxel on the other. Meanwhile, the sub-grid geometry and flow features in the unresolved pixels or voxels are extracted in terms of the local porosities and the local permeabilities. This procedure is accomplished via industrial microcomputed tomography ($\mu \mathrm{C}\mathrm{T}$). Under a chosen scan resolution, the $\mu \mathrm{C}\mathrm{T}$ converted the porous structure of a piece of mineral core into a digital core of 3D CT values. The digital core is further segmented into ternary components, consisting of the pure pore domains for big pores, the unresolved domains for minor pore-solid structures, and the pure solid domains for solid walls [25]. Fluids which flow through the pure pore and the unresolved domains, bounded by the pure solid domains, are controlled by the Brinkman equation [1,20]. This equation is a generalization of the Navier-Stokes equation by adding a linear dragging force term with respect to the velocity. The dragging term depicts the friction between the fluid and the solid within the unresolved domains, and vanishes in the pure pore domains where the Brinkman equation reduce to the Navier-Stokes equation. In the unresolved domains, the strength of dragging is controlled by the permeability to porosity ratio. When the ratio is high, the dragging term is weak and the Brinkman equation falls to the Navier-Stokes limit; when the ratio is low, the flow is dragging force dominant and the Brinkman equation approaches to the Darcy flow limit. Therefore, the Brinkman equation is a perfect model to characterize multiscale porous flows with automatic scaling behavior. One thing is noted that near the Darcy flow limit, both the porosity and permeability are very low, the porous media is pretty like a solid wall and the dragging force becomes so strong to freeze any movement immediately. In this case, numerical simulation becomes difficult. Because the strong force needs to be resolved by a fine enough spatial solution, otherwise it becomes a stiff force numerically. The stiffness of the dragging force is controlled by the local Darcy number, which we define as the quotient of the permeability to porosity ratio over the resolution squared. As the local Darcy number decreases, the stiffness of the dragging force appears increasingly severe. This usually causes that the numerical solutions strongly deviate from correct solutions or even diverge. We will get back to the local Darcy number later when discussing integration scheme.

Numerical methods for the Brinkman equation are either explicit or implicit. The most frequently applied implicit solver for the Brinkman equation is the PIMPLE solver [24], which combines the PISO (Pressure Implicit Splitting Operator) algorithm [22] and the SIMPLE (Semi-Implicit Method of Pressure-Linked Equations) algorithm [21,23] to simulate the transient velocity field of the Brinkman equation. The PIMPLE solver is a finite volume method (FVM) type of algorithm designed for incompressible flows, in which the divergency of the velocity field must be zero. So, solving for the demanding Poisson pressure equation is required at each time step. This means the PIMPLE solver is both computationally expensive and memory demanding as the number of control volume is growing huge in 3D cases. Alternatively, the explicit methods normally close the momentum equation by an equation of state (EOS). This means solving the Poisson pressure equation is unnecessary and the solver is weakly compressible but fast in time integration. The lattice Boltzmann method (LBM) [16] is one of the typical weakly compressible methods, which has two models especially designed for porous flow, namely, the Gray model LBM [13,14] and the Brinkman force LBM [15,17,18]. The Gray model adds a “partial bounce-back” step in the unresolved domains after the “streaming” and “collision” steps of the classical LBM for modeling the resistance from the sub-grid solid to the flow. This model is subject to the fixed permeability-solid fraction relationship, from which the local permeability is binding to the local porosity. The Brinkman force model directly incorporates the linear dragging term of the Brinkman equation into the classical LBM body forcing. Thus, this model solves the Brinkman equation in the same way as solving the Navier-Stokes equation with a specific body force. Though, both of the LBM models show very high efficiency, their shortcomings are apparent as well. Firstly, the fixed permeability-solid fraction relationship is inconsistent to the real digital core. Moreover, the Gray model is subject to numerical inconsistency due to spurious flow oscillation and viscosity dependent permeability. The Brinkman force model, though suffers from less severe inconsistency, is incorrect when the local Darcy number mentioned above is low. This is because the forcing in LBM uses central scheme, which is only valid when the varying of the force is moderate during one time step. But the typical case of a digital core is that the unresolve domain has very low local Darcy number, which results in a stiff rapidly decaying dragging force. In this case, the central scheme forcing in LBM is unacceptable.

The smoothed particle hydrodynamics (SPH) method [2] is another computational fluid dynamics (CFD) method quite different from FVM and LBM. Essentially, it is a many-body particle dynamics simulator under the Lagrangian coordinate system for a broad range of physical phenomena, including the astrophysics [3], the fluid dynamics [27] and the solid mechanics [4]. In fluid dynamics, SPH uses a bunch of regularly distributed particles to represent the fluid and the solid boundary. The field qualities, such as the velocity, the density and the pressure of the fluid, at any specified position, are interpolated through integration of the quantity with a pre-defined kernel function. During simulation, the velocity and the position of each fluid particle are updated with a time integration scheme. This scheme is either explicit or implicit, depending on the choice of the method to close the momentum equation. The WCSPH [5] uses an EOS and is an explicit method. The incompressible SPH (ISPH) [6,7] requires the velocity divergence equal to zero and is a semi-implicit method. At early stage, both WCSPH and ISPH were blamed for its low accuracy due to tensile instability [4,5], which results in particle clumping and forming void regions. This was largely mitigated by using particle shift technique [7,8] or transport velocity correction scheme [9,10,11]. Meanwhile, as the incorporation of the low dissipation Riemann solver [12] into the community of SPH, the accuracy of both WCSPH and ISPH are apparently improved. But SPH is still rarely utilized in solving the Brinkman equation for evaluating the permeability of digital cores. The reason is that on one hand, as an implicit method, ISPH has no advantage over PIMPLE; on the other, as an explicit method, the conventional WCSPH is still highly inefficient, since the time step is extremely small when the local Darcy number tends to zero in the unresolved domains, where the local porosity and permeability are both very low. Recently, we discovered that the tiny time step issue of WCSPH method for the Brinkman flow can be tackled through a modified scheme. This will encourage the application of the powerful WCSPH method in broader applications.

This paper first expands the variable-porosity viscous term of the Brinkman equation. The velocity-linear contribution is combined with the Darcy drag, and a non-negative limiter is introduced when the combined coefficient becomes negative. This formulation stabilized viscous term when the porosity is violently changing. We then address the small timestep restriction of conventional WCSPH. The Strang splitting method is used along with an exact integration for the drag subproblem when the drag stiffness is moderate. The new integration scheme enables the usage of a considerably larger constrained only by the CFL number, the viscosity and the global maximum acceleration. In strongly drag-dominated regions, the viscous term is neglected and an exact integration scheme is used with the larger timestep. A local criterion based on the local Darcy number selects the integration scheme for each particle. The numerical study evaluates constant-porosity accuracy and demonstrates the behavior of the formulation in variable-porosity and digital-core configurations.

The remainder of this paper is organized as follows. Section 2 presents the reorganized Brinkman equation, the corresponding WCSPH spatial discretization, and the details of the two stiff-drag integration strategies. Section 3 reports channel-flow tests, periodic porous-obstacle examples, and 2D and 3D digital-core demonstrations. Section 4 discusses the scope, limitations, and further validation required for quantitative permeability prediction.

2. Materials and Methods

2.1. The Brinkman Equation

The Brinkman equation in the Eulerian form is written as Eq. (1) [1,41].

$${\displaystyle \frac{\rho }{\varphi }}\left[{\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\left(\mathit{u}·\nabla \right){\displaystyle \frac{\mathit{u}}{\varphi }}\right]=-\nabla p+\nabla ·\left({\displaystyle \frac{\mu }{\varphi }}\nabla \mathit{u}\right)-{\displaystyle \frac{\mu }{\kappa }}\mathit{u}+\rho \mathit{g},$$

(1)

where $\mathit{u}$is the flow velocity, $\mathit{g}$is the external acceleration, $\rho$ is the density of the fluid, $p$ is the pressure of the flow, $\mu$ is the dynamic viscosity of the fluid. $\varphi$ and $\kappa$ are two position-based functions of the porous media, representing the local porosity and the local permeability, respectively.

In designed testing porous flow problems, the image of the porous media of interest is well segmented into three categories, representing the pure pores, the unresolved regions and the pure solid, where values of the porosity and the permeability are imposed for each pixel (in 2D) or each voxel (in 3D). However, in the real porous flow problems for a digital core, the CT values (or reduced CT values) of the porous media along with the resolution $\delta x$, are all we have at hand. Thus, it is necessary to discuss a little bit about how we obtain the local porosity and the local permeability from the CT values.

In order to obtain the local porosity, we need to know the bulk porosity of the core and this is down by experiment at large. Then, assuming the density and the effective atomic number of each mineral in the unresolved domain are similar, we use a truncated linear function to model the CT value to pixel (in 2D) or voxel (in 3D) porosity relationship, in Eq. (2) [25]:

$$\varphi \left(c\right)=\left\{\begin{array}{cc}1,& c\le C_0,\\ {\displaystyle \frac{C_1-c}{C_1-C_0}},& C_0<c<C_1,\\ 0,& c\ge C_1,\end{array}\right.$$

(2)

where $C_0$ and $C_1$ are the threshold CT values for segmenting pure pore domains and pure solid domains, respectively. We can set and adjust the two threshold values with the help of any visualization software, by ensuring the averaged porosity equal to the bulk porosity of the core. For instance, Figure 1 shows that the blue region with low CT value represents the pure pore domain, the red region with high CT value represents the pure solid domain, and the region with gradually changing colors from blue to red represents the unresolved domain with sub-grid porous structure.

In order to obtain the local permeability, we select a few positions of the core with different porosities, drill the core at these positions for little cuts, scan the cuts with a much finer resolution by a CT device, segment the scanned CT imaging for binary structure containing only pores and solid, compute the porosity of each cut. Then, we have two methods for the permeability of each cut. The slower method is to use the pore-scale CFD simulation to solve the Navier-Stokes equation for the averaged velocity in the cut and then compute the permeability with Darcy’s law. The faster method is to use the pore-morphology method to compute the pore size distribution of the cut, compute the weighted average pore diameter $d$of the cut with the pore size distribution, and then use Eq. (3) [25,26]:

$$\kappa ={\displaystyle \frac{d^2\varphi }{12{\left(1-0.49\mathrm{l}\mathrm{n}\varphi \right)}^2\min \left\{400{\left(1-\varphi \right)}^2,1\right\}}},$$

(3)

to approximate the permeability of the cut. The limiter in Eq. (3) ensures $\kappa \to \infty ,$ when $\varphi \to 1$. Once the permeability of the cuts with different porosities are available, we can interpolate or fit a curve for the local porosity-permeability relation, then locate permeability value for each voxel is obtained.

Back to the Brinkman equation, the Eulerian form of Eq. (1) can be rewritten in the Lagrangian form as Eq. (4):

$${\displaystyle \frac{D\mathit{v}}{Dt}}=-{\displaystyle \frac{\nabla p}{\rho }}+{\displaystyle \frac{\mu }{\rho }}\nabla ·\left[{\displaystyle \frac{1}{\varphi }}\nabla \left(\mathsf{\varphi }\mathit{v}\right)\right]-{\displaystyle \frac{\varphi \mu }{\kappa \rho }}\mathit{v}+\mathit{g},$$

(4)

where $\mathit{v}=\mathit{u}/\varphi$, $D/Dt=\partial /\partial t+(\mathit{v}·\nabla )$. The viscous term in Eq. (4) can be split into three terms as Eq. (5):

$${\displaystyle \frac{\mu }{\rho }}\nabla ·\left[{\displaystyle \frac{1}{\varphi }}\nabla \left(\mathsf{\varphi }\mathit{v}\right)\right]={\displaystyle \frac{\mu }{\rho }}\left\{{\nabla }^2\mathit{v}+\left[\nabla \left(\mathrm{l}\mathrm{n}\varphi \right)·\nabla \right]\mathit{v}+\mathit{v}{\nabla }^2\left(\mathrm{l}\mathrm{n}\varphi \right)\right\}.$$

(5)

This formulation has two meanings in physics. The first is to provide viscous dissipation to the flow. The second is to make disturbance when the fluid flows through the porous region with varying porosities. When $\varphi$ is a spatial constant, the formulation reduces to the normal viscous term as in the Navier-Stokes equation. When $\varphi$ is spatially varying, the second and the third term in Eq. (5) start to disturb the flow field. It is noticed that the second term contains derivatives of velocity, which is similar to the first term. Whereas the third term is linear in velocity, which is in accordance with the dragging force term. Therefore, we reorganize the Brinkman equation by absorbing the third term of Eq. (5) into the dragging force as Eq. (6):

$$\begin{array}{c}{\displaystyle \frac{D\mathit{v}}{Dt}}=-{\displaystyle \frac{\nabla p}{\rho }}+{\displaystyle \frac{\mu }{\rho }}\left\{{\nabla }^2\mathit{v}+\left[\nabla \left(\mathrm{l}\mathrm{n}\varphi \right)·\nabla \right]\mathit{v}\right\}-{\displaystyle \frac{\mu }{\rho }}C\left(\varphi ,\kappa \right)\mathit{v}+\mathit{g},\\ C\left(\varphi ,\kappa \right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[{\displaystyle \frac{\varphi }{\kappa }}-{\nabla }^2\left(\mathrm{l}\mathrm{n}\varphi \right),0\right].\end{array}$$

(6)

where, $C\left(\varphi ,\kappa \right)$ is the coefficient of the dragging force. Though $\varphi /\kappa$ is dominant in the coefficient for most cases, the use of limiter ensures that the dragging term is absolutely dissipative. We use the artificial equation of state in Eq. (7) to close Eq. (6).

$$p=c^2\left(\rho -{\rho }_0\right),$$

(7)

where $c$ is the sound speed and ${\rho }_0$ is the reference density. We set $c\ge 10V_{max}$, where $V_{max}$ is the maximum expected speed of the flow, to suppress the pressure oscillation down to below 1%. Since, the continuity equation is not directly solved in the chosen WCSPH formulation, where density is evaluated by interpolation with the kernel function and the mass conservation is strictly guaranteed, Eq. (6) and Eq. (7) are the equations we are solving for the rest of the paper. In order to estimate the bulk permeability of the digital cores, we use the no-slip boundary condition at the solid walls and periodic boundary condition at the physical boundary of the digital cores.

2.2. The Classical WCSPH of the Brinkman Equation

2.2.1. The WCSPH Formulation

For a digital core with pure pore domains, unresolved domains and pure solid domains, each pixel (in 2D) or voxel (in 3D) is initially assigned a particle at the center position of its cell. The particles in pure solid domains are labeled as fixed solid particles, and the rest particles are labeled as fluid particles. Each fluid particle $i$interacts with its neighboring fluid particles and solid particles within the range of the support of its kernel function $W\left({|\mathit{r}-\mathit{r}}_i\right|)$, where ${\mathit{r}}_i$ is the position of particle $i$. One can be chosen any positive function with compact support as the kernel function $W\left(x\right)$. Here, we use the Wendland C2 function [28,29] in Eq. (8) as the kernel function:

$$W\left(x\right)=\alpha {\left(2-x\right)}^4\left(2x+1\right),$$

(8)

where $x={|\mathit{r}-\mathit{r}}_i|/h$, $h$ is the smoothing length deciding the influence range of the kernel function. Here, 2 times the CT imaging resolution is used as the smoothing length, namely, $h=2\delta x$. The coefficient $\alpha ={\displaystyle \frac{7}{64\pi h^2}}$ in 2D, and $\alpha ={\displaystyle \frac{21}{256\pi h^3}}$ in 3D, to ensure the exact integration of $W\left(x\right)$ is unity. The density of particle $i$ is evaluated by the interpolation in Eq. (9):

$${\rho }_i=m{\sum }_j{W}_{ij},$$

(9)

where $m$ is the mass of the fluid particle, $W_{ij}=W\left(r_{ij}\right)$, where ${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$ and ${\overline{r}}_{ij}=\left|{\mathit{r}}_{ij}\right|$, particle $j$ is located in the support of $W\left({|\mathit{r}-\mathit{r}}_i\right|)$. The density of all solid particles is ${\rho }_0$. The SPH discretization of Eq. (6) is Eq. (10):

$$\begin{array}{c}\begin{array}{c}{\displaystyle \frac{D{\mathit{v}}_i}{Dt}}={\mathit{a}}_i^p+{\mathit{a}}_i^{\mu }+{\mathit{a}}_i^d+\mathit{g},\\ {\mathit{a}}_i^p=-{\displaystyle \frac{2}{{\rho }_i}}{\sum }_j{P}_{ij}^{*}{V}_j{\nabla }_i{W}_{ij},\\ {\mathit{a}}_i^{\mu }={\displaystyle \frac{\mu }{{\rho }_i}}{\sum }_j{\displaystyle \frac{{\beta }_{ij}{V}_j{\mathit{v}}_{ij}}{{\overline{r}}_{ij}}}{\displaystyle \frac{\partial W_{ij}}{\partial {\overline{r}}_{ij}}},\end{array}\\ \begin{array}{c}{\beta }_{ij}=\left\{\begin{array}{cc}2-{\mathit{r}}_{ij}·\nabla \left[\mathrm{l}\mathrm{n}\varphi \left({\mathit{x}}_i\right)\right],& j\notin \mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d},\\ 2,& j\in \mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d},\end{array}\right.\\ {\mathit{a}}_i^d=-{\displaystyle \frac{\mu }{{\rho }_i}}C\left(\varphi ,\kappa ;{\mathit{x}}_i\right){\mathit{v}}_i,\\ C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{{\displaystyle \frac{\varphi \left({\mathit{x}}_i\right)}{\kappa \left({\mathit{x}}_i\right)}}-{\nabla }^2\left[\mathrm{l}\mathrm{n}\varphi \left({\mathit{x}}_i\right)\right],0\right\},\end{array}\end{array}$$

(10)

where ${\mathit{a}}_i^p, {\mathit{a}}_i^{\mu }$ and ${\mathit{a}}_i^d$are the acceleration of particle $i$ due to the pressure force, the viscous force and the dragging force, respectively, ${\mathit{v}}_{ij}={\mathit{v}}_i-{\mathit{v}}_j$ is the velocity difference of two fluid particles, but when $j$ is a solid particle, in order to implement the no-slip boundary condition at the solid wall, ${\mathit{v}}_{ij}={\mathit{v}}_i-(2{\mathit{v}}_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}-{\mathit{v}}_i)=2{\mathit{v}}_i$ in this paper for stationary walls. $V_j=m/{\rho }_j$ is the volume of particle $j$. ${\mathit{x}}_i$ is the coordinate of the pixel or voxel associated with ${\mathit{r}}_i$. Supposing the data for the structure of a 3D porous media has the voxel dimensions of $L_1, L_2$and $L_3$ in the length (X axis), the width (Y axis) and the depth (Z axis), ${\mathit{x}}_i$is computed in Eq. (11) as:

$$\begin{array}{c}{X}_i=⌊r_{i,X}/\delta x⌋ \mathrm{m}\mathrm{o}\mathrm{d} L_1,\\ Y_i=⌊r_{i,Y}/\delta x⌋ \mathrm{m}\mathrm{o}\mathrm{d} L_2,\\ Z_i=⌊r_{i,Z}/\delta x⌋ \mathrm{m}\mathrm{o}\mathrm{d} L_3,\end{array}$$

(11)

where $X_i, Y_i$ and $Z_i$ are the three components of ${\mathit{x}}_i$in dimension of the length, the width and the depth of the 3D data. Similarly, $r_{i,X}, r_{i,Y}$ and $r_{i,Z}$ are the three components of ${\mathit{r}}_i$. The operator $⌊·⌋$ is to take the closest integer no greater than the target value. The operator mod is short for modulo. For 2D case, ${\mathit{x}}_i$is computed likewise. Since the values of porosity $\varphi$and permeability $\kappa$ depend only on the specific pixel (in 2D) or voxel (in 3D), the values of $\varphi$and $\kappa$, as well as the derivatives of $\mathrm{l}\mathrm{n}\varphi$, for each pixel or voxel, are computed by using Eq. (2)-(3) or other given information and are saved before the simulation. The derivatives of $\mathrm{l}\mathrm{n}\varphi$ are calculated by central difference schemes with the periodic extension of the CT imaging to account for boundary effects. It is noticed that $\mathrm{l}\mathrm{n}\varphi$ is singular at the solid pixels or voxels, where $\varphi =0$. To avoid this singularity, the central difference scheme reduces to one-sided forward or backward difference scheme adjacent to the solid boundary, where $\mathrm{l}\mathrm{n}\varphi$ at the singularity is not used. In the pressure acceleration term, $P_{ij}^{*}$ is the mean pressure of particle $i$ and particle $j$, with a Riemann solver correction, proposed by C. Zhang, X.Y. Hu and N. A. Adams [12] for more accurate pressure fields. The scheme of $P_{ij}^{*}$ is different depending on whether $j$ is a fluid particle or a solid particle. When particle $j$is fluid, the scheme is written in Eq. (12.1):

$$P_{ij}^{*}={\displaystyle \frac{c^2\left({\rho }_i+{\rho }_j-2{\rho }_0\right)}{2}}-{\displaystyle \frac{{\beta }_1{\rho }_0}{2}}{\displaystyle \frac{{\mathit{v}}_{ij}·{\mathit{r}}_{ij}}{{\overline{r}}_{ij}}},$$

(12.1)

where ${\beta }_1$ is a limiter defined in Eq. (12.2):

$${\beta }_1=\min \left\{{\eta }_1\max \left(-{\displaystyle \frac{{\mathit{v}}_{ij}·{\mathit{r}}_{ij}}{{\overline{r}}_{ij}}},0\right),c\right\},$$

(12.2)

where the coefficient is suggested to take ${\eta }_1=3$ and is adopted in this paper. When particle $j$is solid, the scheme of $P_{ij}^{*}$ is written in Eq. (12.3):

$$P_{ij}^{*}=c^2\left({\rho }_i-{\rho }_0\right)-{\displaystyle \frac{{\rho }_0}{2}}\left({\displaystyle \frac{\mathit{g}·{\mathit{r}}_{ij}}{{\overline{r}}_{ij}}}+2{\beta }_2{\mathit{v}}_i·{\mathit{n}}_j\right)$$

(12.3)

where${\mathit{n}}_j$ is the unit normal vector of the boundary pointing from the solid wall to the fluid domain. The unit normal vector is computed by ${\mathit{n}}_j=\mathsf{\Phi }\left({\mathit{r}}_j\right)/\left|\mathsf{\Phi }\right({\mathit{r}}_j\left)\right|$, where the vector $\mathsf{\Phi }\left({\mathit{r}}_j\right)$ is interpolated by the solid neighbor particles of particle $j$ in Eq. (12.4):

$$\mathsf{\Phi }\left({\mathit{r}}_j\right)=-{\sum }_{k\in \mathrm{S}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}{\nabla }_j{W}_{jk}.$$

(12.4)

The limiter ${\beta }_2$ is defined in Eq. (12.5):

$${\beta }_2=\min \left\{{\eta }_2\max \left(-{\mathit{v}}_i·{\mathit{n}}_j,0\right),c\right\},$$

(12.5)

where ${\eta }_2$is also set to 3. Since the solid particles are stationary, the unit normal vector field ${\mathit{n}}_j$of the solid wall is computed once before the simulation.

2.2.2. The Particle Shift Technique

The particle shift technique proposed by C. Huang, D.H. Zhang, Y.X. Shi, Y.L. Si and B. Huang [8] is applied to mitigate the tensile instability. The particle shifting displacement and corresponding change of velocity are computed in Eq. (13.1) and Eq. (13.2):

$$\delta {\mathit{r}}_i=-{\displaystyle \frac{U_{\mathrm{m}\mathrm{a}\mathrm{x}}h\delta t}{2}}{\sum }_j\left\{1+0.24{\left[{\displaystyle \frac{W_{ij}}{W\left(\delta x\right)}}\right]}^4\right\}{V}_j{\nabla }_i{W}_{ij},$$

(13.1)

$$\delta {\mathit{v}}_i=\delta {\mathit{r}}_i·\nabla {\mathit{v}}_i,$$

(13.2)

where $\delta t$ and $U_{\mathrm{m}\mathrm{a}\mathrm{x}}$are the current time step and maximum speed, respectively. We apply this particle shift technique at the end of each round of time integration.

2.2.3. The Classical Time Integration Scheme

The position-based Verlet scheme [30] is applied for the above WCSPH formulation. Position-based Verlet scheme is a symplectic integrator, which comprises of a half step of pre-update for position, a full step of update for velocity and a half step of post-update for position, sequentially. Denoting the current timestamp by the superscript $n$, the position-based Verlet time integration scheme is written in Eq. (14.1) – Eq. (14.3):

$${\mathit{r}}_i^{n+1/2}={\mathit{r}}_i^n+{\mathit{v}}_i^n\delta t_1^{n-1}/2,$$

(14.1)

$${\mathit{v}}_i^{n+1}={\mathit{v}}_i^n+{\displaystyle \frac{D\mathit{v}}{Dt}}\left({\mathit{r}}_i^{n+1/2},{\mathit{v}}_i^n,{\rho }_i^{n+1/2},{\mathit{x}}_i\right)\delta t_1^n,$$

(14.2)

$${\mathit{r}}_i^{n+1}={\mathit{r}}_i^{n+1/2}+{\mathit{v}}_i^{n+1}\delta t_1^n/2,$$

(14.3)

where $\delta t_1^{n-1}$ is the previous time step. After obtaining the intermediate position from Eq. (14.1), we compute the intermediate density and the coordinate ${\mathit{x}}_i$of the pixel (in 2D) or voxel (in 3D) where ${\mathit{r}}_i^{n+1/2}$ is associated with. Then, we compute the acceleration D$\mathit{v}/Dt$ by Eq. (10) with the intermediate position ${\mathit{r}}_i^{n+1/2}$, the current velocity ${\mathit{v}}_i^n$, the intermediate density ${\rho }_i^{n+1/2}$ and the pixel’s or voxel’s coordinate ${\mathit{x}}_i$. The time step $\delta t_1^n=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\delta t^c,\delta t^{\mu },\delta t^a,\delta t^d\right\}$is applied, according to the restriction of the CFL number, the viscosity, the acceleration and the dragging, respectively, in Eq. (15.1) – Eq. (15.4):

$$\delta t^c={\displaystyle \frac{1}{4}}{\displaystyle \frac{h}{U_{\mathrm{m}\mathrm{a}\mathrm{x}}+c}},$$

(15.1)

$$\delta t^{\mu }={\displaystyle \frac{1}{8}}{\displaystyle \frac{h^2{\rho }_0}{\mu }},$$

(15.2)

$$\delta t^a={\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{h}{\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|D\mathit{v}/Dt\right|\right\}}}}$$

(15.3)

$$\delta t^d=\mathrm{m}\mathrm{i}\mathrm{n}\left\{{\displaystyle \frac{\left|{\nabla }^2\mathit{v}\right|}{C\left(\varphi ,\kappa \right)\left|\mathit{v}\right|}}\right\}\delta t^{\mu }.$$

(15.4)

The detailed implementation of the classical WCSPH for the Brinkman equation is listed in Algorithm 1.

Algorithm 1: Algorithm of the classical WCSPH for the Brinkman equation

1. $\mathrm{Read} \mathrm{the} \mathrm{image} \mathrm{of} \mathrm{the} 2\mathrm{D} \mathrm{or} 3\mathrm{D} \mathrm{porous} \mathrm{media} \mathrm{of} \mathrm{interest} \mathrm{from} \mathrm{an} \mathrm{imaging} \mathrm{file}. \mathrm{Decide} \mathrm{the} \mathrm{resolution} \delta x$ $\mathrm{of} \mathrm{the} \mathrm{image} \mathrm{and} \mathrm{obtain} \mathrm{the} \mathrm{local} \mathrm{porosity} \varphi$$\mathrm{and} \mathrm{the} \mathrm{local} \mathrm{permeability} \kappa$$\mathrm{for} \mathrm{each} \mathrm{pixel} \mathrm{or} \mathrm{voxel}, \mathrm{either} \mathrm{from} \mathrm{Eq}. \left(2\right) \mathrm{and} \mathrm{Eq}. \left(3\right), \mathrm{or} \mathrm{from} \mathrm{other} \mathrm{given} \mathrm{information}. \mathrm{Compute} \mathrm{the} \mathrm{dragging} \mathrm{coefficient} C\left(\varphi ,\kappa \right)$;

2. $\mathrm{Initiate} \mathrm{a} \mathrm{particle} \mathrm{at} \mathrm{the} \mathrm{center} \mathrm{of} \mathrm{each} \mathrm{pixel} \mathrm{or} \mathrm{voxel} \mathrm{and} \mathrm{classify} \mathrm{all} \mathrm{the} \mathrm{particles} \mathrm{as} \mathrm{the} \mathrm{set} \mathrm{of} \mathrm{fluid} \mathrm{and} \mathrm{the} \mathrm{set} \mathrm{of} \mathrm{solid} \mathrm{according} \mathrm{to} \mathrm{the} \mathrm{value} \mathrm{of} \mathrm{porosity} \mathrm{of} \mathrm{the} \mathrm{pixel} \mathrm{or} \mathrm{voxel}. \mathrm{Build} \mathrm{up} \mathrm{the} \mathrm{Cell}-\mathrm{Linked} \mathrm{lists} \mathrm{with} \mathrm{the} \mathrm{length} \mathrm{of} \mathrm{the} \mathrm{cells} \mathrm{being} 4\delta x$$, \mathrm{for} \mathrm{the} \mathrm{fluid} \mathrm{particles} \mathrm{and} \mathrm{the} \mathrm{solid} \mathrm{particles}, \mathrm{respectively}. \mathrm{Setup} \mathrm{up} \mathrm{the} \mathrm{hydrodynamic} \mathrm{parameters}, \mathrm{as} \mathrm{the} \mathrm{reference} \mathrm{density} {\rho }_0$$\mathrm{and} \mathrm{the} \mathrm{dynamic} \mathrm{viscosity} \mu$$. \mathrm{Set} \mathrm{the} \mathrm{value} \mathrm{for} \mathrm{the} \mathrm{external} \mathrm{acceleration} \mathit{g}$;

3.  $\mathrm{Set} \mathrm{the} \mathrm{total} \mathrm{simulation} \mathrm{time} T$$. \mathrm{Set} \mathrm{all} \mathrm{fluid} \mathrm{particles} \mathrm{to} \mathrm{stationary} \mathrm{state}, \mathrm{and} \mathrm{compute} \delta t_1^0$$\mathrm{from} \mathrm{Eq}. \left(15.1\right) \text{–} \mathrm{Eq}. \left(15.4\right). \mathrm{Set} t=0;$

4. $\mathrm{While} t<T$, do

5.    $\mathrm{For} n=\mathrm{0,1},\dots ,4$

6.      $\mathrm{Update} \mathrm{the} \mathrm{position} {\mathit{r}}_i^{n+1/2}$$\mathrm{of} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$ with Eq. (14.1);

7.      $\mathrm{Compute} \mathrm{the} \mathrm{intermediate} \mathrm{density} {\rho }_i^{n+1/2}$ with Eq. (9);

8.      $\mathrm{Compute} \mathrm{the} \mathrm{pixel}\text{’}\mathrm{s} \mathrm{or} \mathrm{voxel}\text{’}\mathrm{s} \mathrm{coordinate} {\mathit{x}}_i$$\mathrm{associated} \mathrm{with} {\mathit{r}}_i^{n+1/2}$;

9.      $\mathrm{Compute} \mathrm{the} \mathrm{net} \mathrm{acceleration} {\displaystyle \frac{D\mathit{v}}{Dt}}\left({\mathit{r}}_i^{n+1/2},{\mathit{v}}_i^n,{\rho }_i^{n+1/2},{\mathit{x}}_i\right)$ with Eq (10);

10.      $\mathrm{Compute} \mathrm{the} \mathrm{maximum} \mathrm{speed} U_{\mathrm{m}\mathrm{a}\mathrm{x}}$$\mathrm{and} \mathrm{the} \mathrm{maximum} \mathrm{acceleration}, \mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|D\mathit{v}/Dt\right|\right\},$$\mathrm{then} \mathrm{update} \mathrm{the} \mathrm{time} \mathrm{step} \delta t_1^n$with Eq. (15.1) – Eq. (15.4);

11.      $\mathrm{Update} \mathrm{the} \mathrm{velocity} {\mathit{v}}_i^{n+1}$$\mathrm{of} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$ with Eq. (14.2);

12.     $\mathrm{Update} \mathrm{the} \mathrm{position} {\mathit{r}}_i^{n+1}$$\mathrm{of} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$ with Eq. (14.3);

12.      $\mathrm{Compute} \mathrm{the} \mathrm{density} {\rho }_i^{n+1}$ with Eq. (9);

13.      $\mathrm{Shift} \mathrm{the} \mathrm{position} \mathrm{and} \mathrm{velocity} \mathrm{of} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$ with Eq. (13.1) and Eq. (13.2);

14.      $t:=t+(\delta t_1^{n-1}+\delta t_1^n)/2;$

15.   End for loop

16.   Update the Cell-Link list for all fluid particles;

17. End while loop

18. Terminate the simulation.

2.3. Time Integration for Stiff Drag

2.3.1. The Strang Splitting Method

From Eq. (15.4), we know that when the coefficient $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)$ of the dragging force in Eq. (10) is huge, the time step size $\delta t^d$ due to the dragging is likely to become very small, which causes a tiny time step size $\delta t$for the time integration. This is a common situation for the flows through a digital core but a disaster to the explicit time integration for the classical WCSPH simulation. However, we can circumvent this obstacle by treating Eq. (10) with the Strang splitting method [27]. Namely, we rewrite Eq. (10) as a splitting form in Eq. (16.1) and Eq. (16.2):

$${\displaystyle \frac{D{\mathit{v}}_i^{*}}{Dt}}={\mathit{a}}_i^p+{\mathit{a}}_i^{\mu }+\mathit{g},$$

(16.1)

$${\displaystyle \frac{D{\mathit{v}}_i^{**}}{Dt}}={\mathit{a}}_i^d.$$

(16.2)

During the time integration, the Strang splitting method is used to update the velocity of fluid particles, like Eq. (14.2) in the position-based Verlet scheme. The Strang splitting method consists of three steps.

Step 1: we first compute the maximum speed $U_{\mathrm{m}\mathrm{a}\mathrm{x}}$and the maximum acceleration, $\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|D{\mathit{v}}^{*}/Dt\right|\right\}$from Eq. (16.1), then update the time step by $\delta t_2^n=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\delta t^c,\delta t^{\mu },\delta t^a\right\}$with Eq. (15.1) – Eq. (15.3). Then, we integrate the velocity for a half time step with Eq. (16.1) to obtain the intermediate velocity in Eq. (17.1):

$${\mathit{v}}_i^{\left(1\right)}={\mathit{v}}_i^n+{\displaystyle \frac{1}{2}}{\displaystyle \frac{D{\mathit{v}}_i^{*}}{Dt}}\delta t_2^n.$$

(17.1)

Step 2: We solve Eq. (16.2) exactly, with the intermediate velocity ${\mathit{v}}_i^{\left(1\right)}$ as initial state. The solution with one time step is written in Eq. (17.2):

$${\mathit{v}}_i^{\left(2\right)}={\mathit{v}}_i^{\left(1\right)}\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\displaystyle \frac{\mu }{{\rho }_i}}C(\varphi ,\kappa ;{\mathit{x}}_i)\delta t_2^n\right].$$

(17.2)

Step 3: Taking ${\mathit{v}}_i^{\left(2\right)}$as the initial state, we integrate the velocity for another half time step with Eq. (16.1) to obtain the final velocity at the next time step in Eq. (17.3):

$${\mathit{v}}_i^n={\mathit{v}}_i^{\left(2\right)}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{D{\mathit{v}}_i^{*}}{Dt}}\delta t_2^n.$$

(17.3)

It is noticed that step 1 and step 3 in the Strang splitting method are independent with the dragging force. So, we are safe to use the explicit time integration scheme with the time step $\delta t_2^n,$which is only restricted by the CFL number, the viscosity, the acceleration. The exact solution ensures the accuracy of step 2 and the Strang splitting method has second order of precision locally in time. During the time integration of the modified method, Eq. (17.1) – Eq. (17.3) replace Eq. (14.2) and the original time step $\delta t_1^n$ in Eq. (14.1) and Eq. (14.3) is replaced by $\delta t_2^n.$

2.3.2. The Exact Integration Scheme Neglecting Viscous Contribution

When the resolved viscous acceleration is apparently small relative to the drag acceleration, the viscous contribution is neglected over one step and Eq. (10) is reduced to Eq. (18). This approximation is intended for the drag-dominated regime where the stiffness of the drag is severe.

$${\displaystyle \frac{D{\mathit{v}}_i}{Dt}}=-{\displaystyle \frac{2}{{\rho }_i}}{\sum }_j{P}_{ij}^{*}{V}_j{\nabla }_i{W}_{ij}-{\displaystyle \frac{\mu }{{\rho }_i}}C(\varphi ,\kappa ;{\mathit{x}}_i){\mathit{v}}_i+\mathit{g}.$$

(18)

Considering as an ODE with respect to the velocity ${\mathit{v}}_i$, Eq. (18) can be solved analytically with the solution in Eq. (19.1) - Eq. (19.3):

$${\mathit{v}}_i\left(t\right)=\left[{\mathit{v}}_i\left(t_0\right)-{\mathit{v}}_i\left(\infty \right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left[-\left(t-t_0\right)/t_i^{\infty }\right]+{\mathit{v}}_i\left(\infty \right),$$

(19.1)

$$t_i^{\infty }={\displaystyle \frac{{\rho }_i}{\mu C(\varphi ,\kappa ;{\mathit{x}}_i)}},$$

(19.2)

$${\mathit{v}}_i\left(\infty \right)=\left(\mathit{g}-{\displaystyle \frac{2}{{\rho }_i}}{\sum }_j{P}_{ij}^{*}{V}_j{\nabla }_i{W}_{ij}\right)t_i^{\infty }.$$

(19.3)

Since the velocity solution is exact, we can apply any time step with it. To be accordance with the Strang Splitting method above, the time step $\delta t_2^n$is used. The velocity update scheme with time step $\delta t_2^n$is written in Eq. (20):

$${\mathit{v}}_i^{n+1}=\left[{\mathit{v}}_i^n-{\mathit{v}}_i\left(\infty \right)\right]\mathrm{e}\mathrm{x}\mathrm{p}\left(-\delta t_2^n/t_i^{\infty }\right)+{\mathit{v}}_i\left(\infty \right).$$

(20)

Meanwhile, the position of fluid particle $i$ at timestamp $t$can be exactly updated in Eq. (21):

$${\mathit{r}}_i\left(t\right)={\mathit{r}}_i\left(t_0\right)-t_i^{\infty }\left[{\mathit{v}}_i\left(t_0\right)-{\mathit{v}}_i\left(\infty \right)\right]\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{t-t_0}{2t_i^{\infty }}}\right)-1\right]+{\mathit{v}}_i\left(\infty \right)\left(t-t_0\right).$$

(21)

The half time step and the complete time step position update with time step $\delta t_2^n$ are written in Eq. (22.1) and Eq. (22.2), respectively:

$${\mathit{r}}_i^{n+1}={\mathit{r}}_i^{n+1/2}-t_i^{\infty }[{\mathit{v}}_i^n-{\mathit{v}}_i\left(\infty \right)]\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\delta t_2^n}{2t_i^{\infty }}}\right)-1\right]+{\mathit{v}}_i\left(\infty \right)\delta t_2^n/2,$$

(22.1)

$${\mathit{r}}_i^{(n+1)+1/2}={\mathit{r}}_i^{n+1/2}-t_i^{\infty }\left[{\mathit{v}}_i^n-{\mathit{v}}_i\left(\infty \right)\right]\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{\delta t_2^n}{t_i^{\infty }}}\right)-1\right]+{\mathit{v}}_i\left(\infty \right)\delta t_2^n.$$

(22.2)

2.3.3. The Integration Selector

The stiffness of the linear drag is related to the local Darcy number ${\mathrm{D}\mathrm{a}}_l={\displaystyle \frac{\kappa }{{\varphi \delta x}^2}}.$ Because $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)$ is approximately ${\displaystyle \frac{\varphi }{\kappa }}$ away from strong porosity curvature, $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)h^2$ scales with $1/{\mathrm{D}\mathrm{a}}_l$ when $h/\delta x$ is fixed. We therefore use $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)h^2$ as a practical indicator of drag stiffness and set a threshold $\xi$ to distinguish between the severely and the moderately stiff situations of the drag, by using Inequality (23):

$$C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)h^2>\xi ,$$

(23)

When Inequality (23) is satisfied, the exact integration scheme by neglecting the viscous contribution is applied; otherwise, the Strang splitting method retaining viscosity is used. The value $\xi =3$ is calibrated from the constant-body-force channel tests reported below because the two updates have comparable errors near this value. It should be noted that the selector should be interpreted as an empirical resolution-dependent criterion rather than a universal physical boundary between viscous and drag-dominated flow.

The integration selector is used by each fluid particle before step 9 of Algorithm 1. The decision is made by checking Inequality (23). If it is satisfied, we skip step 9 and the exact integration is applied; otherwise, we continue with step 9 and follow the position-based Verlet scheme with Strang splitting method. Besides, it is noticed that, in step 10 of Algorithm 1, the maximum acceleration is obtained only from fluid particles with the Strang splitting method.

The detailed implementation of the modified WCSPH for the Brinkman equation is listed in Algorithm 2.

Algorithm 2: Algorithm of the modified WCSPH for the Brinkman equation

The parameters setup, the particles and Cell-Linked lists initiation, and other preparation steps before the simulation are the same as step 1 – step 3 in Algorithm 1 and we don’t restate here;

2. $\mathrm{Set} \mathrm{a} \mathrm{label} \mathrm{variable} l$$\mathrm{for} \mathrm{all} \mathrm{fluid} \mathrm{particles} \mathrm{and} \mathrm{initialize} \mathrm{it} \mathrm{as} l\equiv 0$$. \mathrm{Set} \mathrm{an} \mathrm{auxiliary} \mathrm{position} \mathrm{variable} {\mathit{r}}^{*}$$\mathrm{for} \mathrm{all} \mathrm{fluid} \mathrm{particles}. \mathrm{Compute} \delta t_2^0$$\mathrm{from} \mathrm{Eq}. \left(15.1\right) \text{–} \mathrm{Eq}. \left(15.3\right). \mathrm{Set} t=0$;

3. $\mathrm{While} t<T$, do

4.    $\mathrm{For} n=\mathrm{0,1},\dots ,4$

5.       $\mathrm{If} l_i=0$$, \mathrm{update} \mathrm{the} \mathrm{position} {\mathit{r}}_i^{n+1/2}$$\mathrm{of} \mathrm{fluid} \mathrm{particle} i$$\mathrm{with} \mathrm{Eq}. \left(14.1\right), \mathrm{otherwise} \mathrm{set} {\mathit{r}}_i^{n+1/2}={\mathit{r}}_i^{*}$;

6.       $\mathrm{Compute} \mathrm{the} \mathrm{intermediate} \mathrm{density} {\rho }_i^{n+1/2}$ with Eq. (9);

7.       $\mathrm{Shift} \mathrm{the} \mathrm{position} \mathrm{and} \mathrm{velocity} \mathrm{of} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$ with Eq. (13.1) and Eq. (13.2);

8.       $\mathrm{Compute} \mathrm{the} \mathrm{intermediate} \mathrm{density} {\rho }_i^{n+1/2}$ with Eq. (9);

9.       $\mathrm{Compute} \mathrm{the} \mathrm{pixel}\text{’}\mathrm{s} \mathrm{or} \mathrm{voxel}\text{’}\mathrm{s} \mathrm{coordinate} {\mathit{x}}_i$$\mathrm{associated} \mathrm{with} {\mathit{r}}_i^{n+1/2}$;

10.     $\mathrm{For} \mathrm{each} \mathrm{fluid} \mathrm{particle} i$$, \mathrm{check} \mathrm{whether} \mathrm{Inequality} \left(23\right) \mathrm{is} \mathrm{satisfied} \mathrm{or} \mathrm{not}. \mathrm{If} \mathrm{it} \mathrm{is} \mathrm{satisfied}, \mathrm{set} l_i=1$$; \mathrm{otherwise}, \mathrm{compute} \mathrm{the} \mathrm{net} \mathrm{acceleration} {\displaystyle \frac{D\mathit{v}}{Dt}}\left({\mathit{r}}_i^{n+1/2},{\mathit{v}}_i^n,{\rho }_i^{n+1/2},{\mathit{x}}_i\right)$$\mathrm{with} \mathrm{Eq} \left(10\right) \mathrm{and} \mathrm{set} l_i=0$;

11.       $\mathrm{Compute} \mathrm{the} \mathrm{maximum} \mathrm{speed} U_{\mathrm{m}\mathrm{a}\mathrm{x}}$$\mathrm{and} \mathrm{the} \mathrm{maximum} \mathrm{acceleration}, \mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|D\mathit{v}/Dt\right|\right\}$$\mathrm{for} \mathrm{fluid} \mathrm{particles} \mathrm{with} l_i=0$$. \mathrm{Then} \mathrm{update} \mathrm{the} \mathrm{time} \mathrm{step} \delta t_2^n=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\delta t^c,\delta t^{\mu },\delta t^a\right\}$with Eq. (15.1) – Eq. (15.3);

12.       $\mathrm{Update} \mathrm{the} \mathrm{velocity} {\mathit{v}}_i^{n+1}$$\mathrm{of} \mathrm{fluid} \mathrm{particle} i$$\mathrm{with} \mathrm{Eq}. \left(17.1\right) \text{–} \mathrm{Eq}. \left(17.3\right), \mathrm{if} l_i=0$$, \mathrm{or} \mathrm{with} \mathrm{Eq}. \left(20\right), \mathrm{if} l_i=1$;

13.       $\mathrm{Update} \mathrm{the} \mathrm{position} {\mathit{r}}_i^{n+1}$$\mathrm{of} \mathrm{fluid} \mathrm{particle} i$$\mathrm{with} \mathrm{Eq}. \left(14.3\right), \mathrm{if} l_i=0$$; \mathrm{otherwise}, \mathrm{update} \mathrm{the} \mathrm{position} {\mathit{r}}_i^{n+1}$$\mathrm{with} \mathrm{Eq}. \left(22.1\right) \mathrm{and} \mathrm{compute} \mathrm{the} \mathrm{auxiliary} \mathrm{position} {\mathit{r}}_i^{*}$ with Eq. (22.2);

14.       $t:=t+(\delta t_2^{n-1}+\delta t_2^n)/2;$

15.   End for loop

16.   Update the Cell-Link list for all fluid particles;

17. End while loop

18. Terminate the simulation.

3. Numerical Results

In this section, 2D channel flows with constant and varying porosity are used to validate the proposed method. Then, the method is applied to 2D porous flows over a periodic lattice of cylinders and flows through 2D and 3D digital cores. For all test cases, the fluid starts at rest with zero initial pressure; flows are driven by a constant acceleration; periodic boundary conditions are applied in all axial directions; no-slip boundary condition is applied at the interface between the fluid and pure solid regions; water is used as the fluid, with density ${\rho }_0={10}^3\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ and dynamic viscosity $\mu =1\mathrm{m}\mathrm{P}\mathrm{a}·\mathrm{s}$.

3.1. 2D Channel Flows with Constant Porosity and Permeability

As shown in Figure 2, we study the 2D flows through a channel along the X axis confined by two parallel plates with separation distance $H$. Each of the plates consists of four layers of pure solid particles with spacing $\delta x$. The channel is filled with unresolved porous media possessing constant porosity $\varphi$ and constant permeability $\kappa$. The coefficient $C\left(\varphi ,\kappa \right)$ of the dragging force reduces to $\varphi /\kappa$. Two dimensionless numbers to characterize this flow are the Reynolds number $\mathrm{R}\mathrm{e}=U_{\mathrm{m}\mathrm{a}\mathrm{x}}H{\rho }_0/\mu$ and the Darcy number $\mathrm{D}\mathrm{a}=\kappa /H^2$, where $U_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum speed of the flow. Since our interest to the flow focuses on the low Reynolds number and low Darcy number regime, we can adjust the values of $\kappa , H$ and the magnitude of the acceleration, i.e., $\left|\mathit{g}\right|=g_X$, to control the values of Re and Da. With low Reynolds number, the well-developed flow is steady flow with an exact solution $\mathit{v}\left(Y\right)=<v_X\left(Y\right), 0>$, where the origin of the coordinate is located at the left side of the channel center and the X-component of the velocity is given in Eq. (24):

$$v_X\left(Y\right)=\left\{1-{\displaystyle \frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left[\sqrt{\frac{\varphi }{\kappa }}\left(Y-\frac{H}{2}\right)\right]}{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\sqrt{\frac{\varphi }{\kappa }}\frac{H}{2}\right)}}\right\}{\displaystyle \frac{{\rho }_0{g}_x\kappa }{\mu \varphi }},$$

(24)

From Eq. (24) the maximum flow speed is $U_{\mathrm{m}\mathrm{a}\mathrm{x}}={\varphi v}_X\left(H/2\right)$.

In the following, we set the channel separation $H=N\delta x$ and implement serval of simulations for the channel flows with five empirical porosity-permeability pairs $\left(\varphi ,\kappa \right)=\left(\mathrm{0.009,0.1}\mathrm{m}\mathrm{D}\right), \left(\mathrm{0.028,1.0}\mathrm{m}\mathrm{D}\right), \left(\mathrm{0.070,10}\mathrm{m}\mathrm{D}\right), \left(0.150,{10}^2\mathrm{m}\mathrm{D}\right), \left(0.302,{10}^3\mathrm{m}\mathrm{D}\right)$ from former digital core analysis. In all simulations, we set the magnitude of acceleration $g_X=$5m/s². The numerical solutions of the flows are obtained along the vertical centerline of the channel via interpolation, and are compared with the exact solutions. The accuracy of the numerical solutions is measured by the root mean square error percentage (RMSEP) [19] $\epsilon$given by Eq. (25):

$$\epsilon =\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{n=1}^N{\left[{\displaystyle \frac{{\tilde{v}}_X\left(Y_n\right)-v_X\left(Y_n\right)}{v_X\left(Y_n\right)}}\right]}^2},$$

(25)

where $v_X\left(Y_n\right)$ and ${\tilde{v}}_X\left(Y_n\right)$ are the exact solution and numerical solution, respectively. Here, we set $N=200.$

In the first group of simulations, we show that $\xi =3$ is a proper threshold for Inequality (23). First, we solve the equation $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)h^2=3$ for $\delta x$ with the five porosity-permeability pairs, respectively. Then, we implement the simulations with the solved resolutions and corresponding porosity-permeability pairs. Actually, for each parameter set, two simulations are conducted, one for the Strang splitting method and the other for the exact integration without viscous force. Thus, during the time integration step of the two simulations, we deliberately switch off the integration selector and impose one of the integration schemes, respectively. Figure 3 compares the two numerical solutions with the corresponding exact solution for the five parameter sets.

The solved critical resolution $\delta x$ when $\xi =3,$ the RMSEP ${\epsilon }_1$ for Strang splitting method and the RMSEP ${\epsilon }_2$ for the exact integration scheme neglecting viscous force, corresponding to the porosity-permeability pairs, are listed in Table 1.

It is observed that for each porosity-permeability pair, the corresponding RMSEPs of the two integration schemes are very close. And all RMSEPs values are below 5%. This demonstrates that the gap between severe and moderate stiff regions of the dragging force can be bridged smoothly by setting the threshold $\xi =3$. Meanwhile, near this threshold, numerical solutions all present satisfactory precision.

Then, we fix the resolution and channel separation to test the accuracy of the proposed method when simulating flows with different Reynolds numbers and Darcy numbers. In Figure 4, comparisons of the velocities are made for $\delta x=5\times {10}^{-8}\mathrm{m}$ and $H={10}^{-5}\mathrm{m}$.

The porosity-permeability pairs of the channel, along with corresponding Reynolds number and Darcy number of the flows, and the RMSEP of the numerical solutions when the channel separation $H={10}^{-5}\mathrm{m}$ are listed in Table 2.

It is noted that during this group of simulations, the integration selectors all point to the Strang splitting method. This shows that for the five numerical cases, the resolution $\delta x=5\times {10}^{-8}\mathrm{m}$ is fine enough to keep the stiffness of the dragging force at a low level. The RMSEP of the five cases are all below 5%, showing good accuracy of the Strang splitting method when the stiffness of the dragging force is weak.

The porosity-permeability pairs of the channel, along with corresponding Reynolds number and Darcy number of the flows, and the RMSEP of the numerical solutions when the channel separation $H={10}^{-4}\mathrm{m}$ is listed in Table 3.

Figure 5 compares the velocities for $\delta x=5\times {10}^{-7}\mathrm{m}$ and $H={10}^{-4}\mathrm{m}$.

During this group of simulations, the integration selectors point to the Strang splitting method when the porosity-permeability pairs $\left(\varphi ,\kappa \right)= \left(0.150,{10}^2\mathrm{m}\mathrm{D}\right),$ $\left(0.302,{10}^3\mathrm{m}\mathrm{D}\right)$; and point to the exact integration scheme without viscous force when the porosity-permeability pairs$\left(\varphi ,\kappa \right)=\left(\mathrm{0.009,0.1}\mathrm{m}\mathrm{D}\right), \left(\mathrm{0.028,1.0}\mathrm{m}\mathrm{D}\right),\left(\mathrm{0.070,10}\mathrm{m}\mathrm{D}\right)$. The RMSEP of the former two cases are between 5%-8%, showing the accuracy of the Strang splitting method decreases a little as the stiffness of the dragging force enhances. When the stiffness of the dragging force surpasses the threshold, the RMSEP of the latter three cases are all below 4%, showing good accuracy of the exact integration scheme without viscous force when the dragging force is severely stiff.

Figure 6 compares the velocities for $\delta x=5\times {10}^{-6}\mathrm{m}$ and $H={10}^{-3}\mathrm{m}$. The porosity-permeability pairs of the channel, along with corresponding Reynolds number and Darcy number of the flows, and the RMSEP of the numerical solutions when the channel separation $H={10}^{-3}\mathrm{m}$ is listed in Table 4.

In this group of simulations, the integration selectors all point to the exact integration scheme without viscous force. This shows that for the five numerical cases, the resolution $\delta x=5\times {10}^{-6}\mathrm{m}$ is not fine enough to resolve the dragging force. But the RMSEP of the five cases are all below 4%. Moreover, the stiffer the dragging force is, the lower RMSEP appears.

In the last group of simulation with constant porosity and permeability, we verify that the numerical solution of the Strang splitting method will greatly deviate from the exact solution in lack of fine resolution. We keep $\delta x=5\times {10}^{-6}\mathrm{m}$ and $H={10}^{-3}\mathrm{m}$ and redo the previous group of simulation. But this time, we switch off the integration selector and impose Strang splitting method during the time integration step. Figure 7 compares the numerical solutions with the corresponding exact solutions.

It is observed that even if the Strang splitting method is stable for all simulation, the numerical solutions all deviate from respective exact solutions when the spatial resolution is not fine enough. Therefore, when the simulation is in lack of fine resolution, or in another word, the dragging force is severely stiff, we need to use the exact integration scheme without viscous force.

3.2. 2D Channel Flows with Varying Porosity and Permeability

When the channel is filled with unresolved porous media with varying porosity and permeability, the derivatives of $\mathrm{l}\mathrm{n}\varphi$ in the coefficient $C\left(\varphi ,\kappa \right)=C\left(\varphi ,\kappa ;X,Y\right)$ of the dragging force start to impact the fluid flow. In this study, we only set varying porosity and permeability along the streamwise direction. Thus, the same as in the previous cases, the flow is only in the streamwise direction and the spanwise component of velocity is zero. The channel length $L=M\delta x$. We set the porosity changing along X axis with the cosine profile $\varphi \left(X\right)={\varphi }_m+{\varphi }_a\mathrm{c}\mathrm{o}\mathrm{s}\left(2\pi X/L\right)$ in one period, where ${\varphi }_m$ and ${\varphi }_a$ are the mean porosity and porosity amplitude, respectively. This setting ensure that the porosity is smoothly changing in the streamwise direction without any jump. Besides, we build a pseudo relation between the porosity and the local average pore diameter: $d\left(X\right)=3\times {10}^{-5}\varphi \left(X\right)+9.7\times {10}^{-7}\left(m\right)$. Substituting this pseudo relation into Eq. (3), we obtain the position dependent permeability $\kappa \left(X\right)$ and the simulations are ready.

We conduct three numerical tests, in which we fix $\delta x=5\times {10}^{-5}m, M=10$ and $N=100.$ So, the channel length $L=5\times {10}^{-4}m$ and the channel separation $H=5\times {10}^{-3}m$. The magnitude of acceleration is fixed at $g_X=$15m/s². With is acceleration, the fully developed flow is steady. We choose five vertical lines: $X=X_i=\left(i+0.5\right)\delta x, i=\mathrm{0,1},\mathrm{2,3},4,$for observation. At each of the vertical lines, we compute the X-component of the reference velocity with Eq. (26):

$${\tilde{v}}_X\left(Y;X_i\right)=\left\{1-{\displaystyle \frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left[\sqrt{C\left(\varphi ,\kappa ;X_i,Y\right)}\left(Y-\frac{H}{2}\right)\right]}{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\sqrt{C\left(\varphi ,\kappa ;X_i,Y\right)}\frac{H}{2}\right)}}\right\}{\displaystyle \frac{{\rho }_0{g}_x}{\mu C\left(\varphi ,\kappa ;X_i,Y\right)}},$$

(26)

After each simulation, we compare the numerical solution with the reference velocity at all vertical observation lines by computing the RMSEP. It should be noted that the reference velocity is not the exact solution of flow. Since the porosity is periodically varying in the streamwise direction, the fluid in the flow is either dragged or pushed by its upwind or downwind neighbors. Although the dragging forces and the pushing forces negate in general, they do not locally. As a result, the flow velocity is attracted toward its mean value. Therefore, at the observation vertical lines with higher porosity, the computed velocity is supposed to lag behind the reference velocity, and at those with lower porosity, the computed velocity is expected to surpass the reference velocity.

Figure 8 compares the numerical solution with the reference velocities at five observation lines when ${\varphi }_m=0.13$ and ${\varphi }_a=0.05$. It is observed that with this strong porosity amplitude, the numerical solution has apparent difference with the reference velocity at respective observation lines, where the RMSEP varies from 4% to 27%. The computed velocity at $X=X_0$ and $X=X_1$ is lower than the corresponding reference velocity, where the porosities $\varphi \left(X_0\right)=0.178$ and $\varphi \left(X_1\right)=0.159$ are greater than mean porosity ${\varphi }_m$. While the computed velocity at $X=X_2$, $X=X_3$ and $X=X_4$ is higher than the corresponding reference velocity, where the porosities $\varphi \left(X_2\right)=0.130$, $\varphi \left(X_3\right)=0.101$ and $\varphi \left(X_4\right)=0.082$ are less than or equal to ${\varphi }_m$. This agrees well with our expectation.

In the second test, Figure 9 compares the numerical solution with the reference velocities at five observation lines when ${\varphi }_m=0.17$ and ${\varphi }_a=0.01$. With weaker porosity magnitude, the numerical solution has minor difference with the reference velocity at respective observation lines, where the RMSEP only varies from 0.4% to 3%. Likewise, the computed velocity is lower than the reference velocity at position with higher porosity, and is higher than the reference velocity at position with lower porosity.

In the last test, Figure 10 compares the numerical solution with the reference velocities at the observation lines when ${\varphi }_m=0.178$ and ${\varphi }_a=0.002$. With this tiny porosity magnitude, the RMSEP at five observation lines all drop to below 1%. The velocity difference between the computed and the reference maintains the same pattern as that appeared in Figure 8 and Figure 9. The magnitude of all velocities at the five observation lines shrinks to a narrow band between $9.3\times {10}^{-8}\mathrm{m}/\mathrm{s}$ and $9.8\times {10}^{-8}\mathrm{m}/\mathrm{s}$. This shows the numerical solution is asymptotically convergent.

3.3. 2D Porous Flows over a Periodic Lattice of Cylinders

The 2D channel flows validated the accuracy of the proposed method and show that the numerical solution is asymptotically convergent when the local porosity and local permeability is continuously varying. However, these test cases are only one-dimensional or simple two-dimensional flows. In order to test two-dimensional flows with apparent variation both in velocity and pressure field, we study the 2D porous flows over a square lattice of cylinders. As shown in Figure 11, the configuration of this problem is converted to a single cylinder placed at the center of a 2D square lattice. The length of the lattice is $L$ and the radius of the cylinder is $R$. The origin of the coordinate is located at the lower left corner of the lattice. The constant acceleration $\mathit{g}$ driving the flow points to the X-axis direction and the periodic boundary conditions at the four edges of the lattice duplicate an infinite periodic array of the same configuration.

Within the square lattice, the area is divided into two regions, one inside the cylinder and the other outside the cylinder. In the numerical tests, we fill the region outside the cylinder using unresolved porous media with constant porosity and constant permeability. While, the region inside the cylinder is either pure solid or filled with unresolved porous media with varying porosity and permeability. The purpose of designing such porous regions is threefold: 1) testing porous flows over a pure solid cylinder and comparing the impact of different porosities and permeabilities to the flow fields; 2) observing a high porosity and high permeability porous flow over a gradually denser obstacle; 3) observing a porous flow with median porosity and permeability over a gradually looser hole.

In the following, we choose the resolution $\delta x={10}^{-6}\mathrm{m}$ and set the lattice length $L=100\delta x={10}^{-4}\mathrm{m}$, the cylinder radius $R=20\delta x=2\times {10}^{-5}\mathrm{m}$. The magnitude of acceleration is fixed at $g_X=$15m/s². We use the same pseudo relation between the porosity and the local average pore diameter: $d\left(\varphi \right)=3\times {10}^{-5}\varphi +9.7\times {10}^{-7}\left(m\right)$. Along with Eq. (3), we obtain the porosity dependent permeability $\kappa \left(\varphi \right)$. The porosity of the region outside and inside the cylinder are denoted as ${\varphi }_{out}$ and ${\varphi }_{in}$, respectively. All simulations are conducted for 0.01s to ensure the flows are well developed.

In the first two numerical tests, the cylinder is pure solid, ${\varphi }_{out}$ is set as 0.9 and 0.15, respectively. The velocity and pressure contours, along with the velocity vector field and streamlines are shown in Figure 12 and Figure 13, respectively. When ${\varphi }_{out}=0.9$, the corresponding permeability is as high as 53.8Da, resulting in a weak dragging force. The maximum speed of the flow is around ${10}^{-3}\mathrm{m}/\mathrm{s}$ and is located at the midpoint of two vertically adjacent cylinders. From Figure 12, it is observed that all contours and the streamlines are very smooth. This flow is very similar to the classical low Reynolds number flow over a periodic lattice of pure solid cylinders, where no porous media is present.

Whereas, when ${\varphi }_{out}=0.15$, the corresponding permeability is 101.8mD, a typical value for a mineral core with median permeability. The stiffness of the dragging force is moderate under current resolution. Comparing with the previous test, the flow is apparently resisted from the porous media outside the cylinder. The maximum speed of the flow is around $1.4\times {10}^{-5}\mathrm{m}/\mathrm{s}$ and is located at all the 1/3 the of lattice region to the right of the cylinder. As viewed from Figure 13, the velocity and pressure contours, as well as the streamlines and the velocity field, show huge differences with the previous case. It is noticed that the streamlines in the middle 1/2 of the lattice region are strongly perturbed by the solid cylinder. The streamlines in front of the cylinder are squeezed and pushed towards the upper and lower side the cylinder in a gathering manner until they pass over it.

In the next two tests, we apply a position dependent porosity ${\varphi }_{in}$ for the region inside the cylinder. The porosity ${\varphi }_{in}$ is changing along the radial direction, given by Eq. (27):

$${\varphi }_{in}\left({\mathit{x}}_i\right)={\displaystyle \frac{{\varphi }_2-{\varphi }_1}{20}}⌊{\displaystyle \frac{|{\mathit{x}}_i-{\mathit{x}}_0|}{\delta x}}⌋+{\varphi }_1,$$

(27)

Where ${\varphi }_1$ and ${\varphi }_2$ are the imposed porosities at the cylinder center and edge, respectively. ${\mathit{x}}_0$ is the coordinate of the cylinder center. The operator $⌊·⌋$ is to take the closest integer no greater than the target value. For the third test, we set ${\varphi }_{out}=0.9, {\varphi }_1=0$ and ${\varphi }_2=0.4$, where the porosity decreases from the edge to the center of the cylinder, modeling a gradually denser obstacle.

From Figure 14, it is observed that the flow pattern is quite similar to that of the first test in this subsection. Again, all contours and the streamlines are very smooth. Differences only appear inside the cylinder, where detailed flow structures in both velocity and pressure are clearly viewed. The streamlines inside the cylinder are only mildly squeezed by the gradually denser porous media.

In the fourth test, we set ${\varphi }_{out}=0.15, {\varphi }_1=0.30$ and ${\varphi }_2=0.15$, where the porosity increases from the edge to the center of the cylinder, modeling a gradually looser hole. It is observed from Figure 15, that this flow has a nearly uniform flow field outside the cylinder, where the velocity aligns with streamwise direction, with the magnitude of around ${10}^{-5}\mathrm{m}/\mathrm{s}$. Inside the cylinder, the velocity field is mildly perturbed by the porous media with varying porosity but reserves nearly symmetric flow structure, where maximum speed of around ${2.1\times 10}^{-5}\mathrm{m}/\mathrm{s}$ appears at the center of the cylinder.

3.4. Flows Through 2D and 3D Digital Cores

The method is next applied to digital-core flows driven by a constant acceleration $\mathit{g}$. The superficial streamwise velocity is evaluated as the bulk-volume average ${\overline{u}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{m}}=\overline{\varphi v}={\displaystyle \frac{1}{v}}{\int }_V\varphi v_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{m}}dV$, with zero contribution from solid pixels or voxels. The apparent bulk permeability ${\kappa }_b$ is then calculated from Darcy’s law in Eq. (28):

$${\kappa }_b={\displaystyle \frac{\mu {\overline{u}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{m}}}{\rho \left|\mathit{g}\right|}}.$$

(28)

A simulation is stopped when the relative change in ${\overline{u}}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{m}}$over 100 timesteps is below ${10}^{-2}$. This operational criterion identifies an approximately steady average flux, but it is not a substitute for a full momentum-residual and mass-balance analysis. The following 2D and 3D cases therefore demonstrate permeability calculation under the adopted model assumptions. Both simulations use$\left|\mathit{g}\right|=10\mathrm{m}/{\mathrm{s}}^2$.

The 2D calculation is a feasibility illustration rather than an independent permeability validation. A core with an experimentally measured bulk porosity of approximately 7% is imaged by $\mu$CT at 15$\mathsf{\mu }\mathrm{m}$ resolution and trimmed to a ${400}^3$ voxel cube. The central slice normal to the Z-axis is used as the 2D domain. Reduced CT values are mapped to local porosity using Eq. (2), with $C_0=16$ and $C_1=148$selected to reproduce the measured average porosity. The same illustrative porosity-pore-diameter relation used in the preceding tests, together with Eq. (3), supplies voxel permeability. A constant acceleration is applied in the X-axis direction.

Figure 16 shows the contour of the velocity component in the streamwise direction, the pressure contour and the porosity contour of the flow through the 2D digital core.

Figure 16 shows that the largest velocities occur along the approximately 0.3-porosity crack, while the background porosity is mostly below 0.1. Under the adopted CT-to-porosity and porosity-to-permeability mappings, Eq. (28) gives an apparent permeability of 39.1mD. Because no measured permeability or independent numerical reference is available for this 2D slice, this value should be interpreted as a model prediction rather than a validated estimate.

The 3D example uses a core with an experimentally measured bulk porosity of approximately 10%, imaged by $\mu$CT at 46$\mathsf{\mu }\mathrm{m}$ resolution and trimmed to 200³ voxels. Equation (2), with $C_0=16$ and $C_1=208$, maps reduced CT values to local porosity. Three selected sub-samples are imaged at 1 μm resolution. Their porosity and pore-size distributions give the pairs$\left(\varphi ,d\right)=\left(0.171, 3.7\mathsf{\mu }\mathrm{m}\right), \left(0.215, 4.7\mathsf{\mu }\mathrm{m}\right)$ and $\left(0.223, 4.7\mathsf{\mu }\mathrm{m}\right)$. The relation $d\left(\varphi \right)=2.0\times {10}^{-5}\varphi +2.2\times {10}^{-7}\left(\mathrm{m}\right)$ and Eq. (3) are then used to assign voxel permeability. Because the three calibration points occupy a narrow porosity interval, use of this relation outside that interval is an extrapolation and is a major source of uncertainty. The body acceleration is applied in the Z-axis direction.

Figure 17 shows preferential flow through higher-porosity voxels. Under the adopted image mappings and convergence criterion, Eq. (28) gives an apparent z-direction permeability of 169.3mD. This prediction has not been compared with laboratory permeability, pore-scale DNS, FVM, or LBM for the same sample. Quantitative validation and calculation of the full permeability tensor remain necessary before the result can be used predictively.

The 200³ voxel calculation was performed on one NVIDIA Titan V GPU and reached the operational flux criterion after 3000 steps, requiring approximately 45 min of wall-clock time and 478 MB of GPU memory.

4. Discussion

This study develops an explicit WCSPH framework for simulating stiff Darcy-Brinkman flow in ternary digital-core domains comprising resolved pores, unresolved porous regions, and solid grains. The variable-porosity viscous operator is decomposed so that its velocity-linear contribution can be combined with the Darcy drag term, yielding a stabilized formulation for strongly heterogeneous porosity fields. To alleviate the severe timestep restriction imposed by stiff drag, two complementary integration strategies are introduced. For moderately stiff drag, Strang splitting isolates the linear drag term, which is integrated analytically within the split update. In strongly stiff drag regions, the viscous contribution is neglected and particle updates are achieved by an exact integration scheme. A local selector based on $C\left(\varphi ,\kappa ;{\mathit{x}}_i\right)h^2$ determines the appropriate update for each particle. These treatments substantially relax the classical drag-induced stability restriction, although the global timestep remains constrained by acoustic propagation, viscous diffusion, particle acceleration.

The constant-porosity channel-flow tests demonstrate that the proposed formulation maintains satisfactory accuracy over the investigated Reynolds and Darcy number ranges, with RMSEP values generally below 8%. The variable-porosity channel tests further show that the numerical response approaches the constant-porosity limit as the amplitude of porosity variation decreases. This behavior supports the internal consistency of the formulation but should not be interpreted as a formal spatial or temporal convergence demonstration because the underlying physical problem changes with the porosity amplitude. The cases of the periodic lattice of cylinders illustrate the capability of the method to represent heterogeneous velocity and pressure fields in the presence of spatially varying porosity and permeability.

The 2D and 3D digital-core examples demonstrate the practical implementation of the method for multiscale porous structures. Under the adopted CT-to-porosity and porosity-to-permeability mappings, the calculated apparent permeabilities are 39.1mD for the 2D slice and 169.3mD in the Z-axis direction for the 3D sample. These values should be regarded as model predictions rather than independently validated permeability measurements. They depend on image segmentation, the limited porosity-pore-size calibration, the stabilized continuum formulation, and the operational steady-state criterion. In particular, the three calibration samples used for the 3D core cover a relatively narrow porosity interval, and extrapolation of the fitted relationship outside this range represents an important source of uncertainty.

Accordingly, the present results establish the proposed WCSPH method as a computationally feasible framework for stiff multiscale digital-core flow, rather than a fully validated permeability predictor. Further work should compare the method with high-resolution FVM, LBM, or pore-scale DNS solutions and, where possible, laboratory permeability measurements for the same samples. Additional investigations should evaluate the sensitivity to the integration threshold, limiter activation, CT segmentation, and local permeability calibration; and calculate the full permeability tensor in the three principal directions under consistent boundary conditions. These developments will be essential for establishing the quantitative accuracy, robustness, and general applicability of the proposed framework.