2.1 DPN flow model

In DPN, the spatial arrangements of pores and throats are obtained directly from segmented images. In this work, this was achieved using the maximal ball method (Dong & Blunt Reference Dong and Blunt2009), which can perform direct mapping of the actual sample to an irregular lattice. We calculated the pore body volumes and throat conductances based on the void space geometry. The pores are commonly represented in cubic or spherical form. For pore throats, however, a variety of geometric shapes with various cross-sections such as stars, circles, triangles, squares, etc., could be used. In this work, we use spherical pore bodies and cylindrical throats. Figures 1 and 2 show a typical arrangement of a network, with pore bodies (nodes) indicated by indices $i$ and $j$, and throats $l_{ij}$ connecting both nodes $i$ and $j$. Depending on the driving mechanism, PN can be implemented either quasi-statically or dynamically. In quasi-static PN modelling, invasion of the pores is based solely on capillary entry pressure. Capillary pressure is imposed through boundary conditions on the network lattice. The position of fluid–fluid interfaces is then determined at each stage of the fluid displacement process independent of the time of the simulation, then the simulation runs until an equilibrium is attained between the prescribed capillary pressure and that of the system. If the interface is unstable, it moves through the network until a stable position is found or until it reaches the outlet. No time dependence is included in the calculation and the interface simply moves from one equilibrium position to another (Celia, Reeves & Ferrand Reference Celia, Reeves and Ferrand1995). Quasi-static models can be considered extensions of simple percolation models, with drainage floods being modelled through invasion percolation and imbibition through adapted bond percolation processes (Li, Mcdougall & Sorbie Reference Li, Mcdougall and Sorbie2017).

Figure 1.

A local PN structure.

Figure 2.

A subsample or elementary volume used for coupling simulations.

Quasi-static systems are not time dependent and are, consequently, unsuitable for the problem at hand due to the transient nature of the coupled system whose fluid and structural properties must evolve with time. Therefore, we employ a dynamic approach. In DPN systems, the model solves the pressure evolution problem based on capillary entry pressure and introduces a time dependent pore-filling mechanism. The emergent rate-dependent flow regimes are determined by a competition between capillary and viscous forces. The DPN model calculates the pressure forces by solving pressure evolution based on minimum filling time for each of the pores, which are eventually imparted onto the solid phase. If the pores are located in domain $\unicode[STIX]{x1D6FA}_{f}$, and fluid flow in the pore throats is laminar, given Poiseuille’s law, each fluid phase obeys the following volume conservation law:

(2.1)$$\begin{eqnarray}V_{i}\frac{\unicode[STIX]{x2202}S_{i}}{\unicode[STIX]{x2202}t}+\mathop{\sum }_{i\in N}Q_{ij}=0,\quad \text{in }\unicode[STIX]{x1D6FA}_{f}\end{eqnarray}$$

where $N$ is the number of pore bodies, $V_{i}$ is the volume of pore body $i$, $S_{i}$ represents the local saturation of the reference phase and $Q_{ij}$ is the total volumetric flux from pore body $i$ to pore body $j$. If a channel or throat is completely filled with a single fluid, the flow rate along the throat from pore body $i$ to pore body $j$ can be calculated by

(2.2)$$\begin{eqnarray}Q_{ij}=G_{ij}(\unicode[STIX]{x1D6F7}_{i}-\unicode[STIX]{x1D6F7}_{j}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}_{k}$ is the hydraulic potential in pore body $k=\{i,j\}$ defined as $P_{k}+\unicode[STIX]{x1D70C}\boldsymbol{g}h_{k}$, where $h_{k}$ is the coordinate of the pore $k$ in the direction of the gravitational field, $P_{k}$ is the fluid pressure in pore $k$, $\boldsymbol{g}$ is the gravitational force and $G_{ij}$ is the throat conductance, defined as

(2.3)$$\begin{eqnarray}G_{ij}=\frac{S_{\unicode[STIX]{x1D6FC}}r_{ij}^{4}}{l_{ij}\unicode[STIX]{x1D702}},\end{eqnarray}$$

where $S_{\unicode[STIX]{x1D6FC}}$ is the shape factor, $r_{ij}$ is the radius of the connecting throat, $l_{ij}$ is the length of the throat taken from the centroids of pore body $i$ and pore body $j$, and $\unicode[STIX]{x1D702}$ is the fluid viscosity. The conductance is calculated at every throat location and is dependent on its radius and length. By applying a pressure gradient across the network (via boundary conditions), the pressure field inside the network can eventually be obtained by applying mass conservation at each interior node $i$,

(2.4)$$\begin{eqnarray}\mathop{\sum }_{i\in N_{i}}Q_{ij}=0.\end{eqnarray}$$

For more than one fluid occupying the throat, the effects of capillary pressure must be accounted for across the fluid–fluid interface. The DPN model can be further simplified by using a circular throat cross-section, which prevents corner flow. For the case of drainage, if the displacement occurs in a piston-like manner and the effects of gravity are neglected, equation (2.2) is rewritten as

(2.5)$$\begin{eqnarray}Q_{ij}=G_{ij}(P_{i}-P_{j}-\unicode[STIX]{x0394}P_{ij}^{c}),\end{eqnarray}$$

where

(2.6)$$\begin{eqnarray}G_{ij}=S_{\unicode[STIX]{x1D6FC}}r_{ij}^{4}\left(\frac{1}{l_{ij}^{(i)}\unicode[STIX]{x1D702}_{i}}+\frac{1}{l_{ij}^{(j)}\unicode[STIX]{x1D702}_{j}}\right).\end{eqnarray}$$

In equation (2.6), $\unicode[STIX]{x1D702}_{i}$ and $\unicode[STIX]{x1D702}_{j}$ are the viscosities of the fluids in segments of the pore throats adjacent to pore $i$ and pore $j,$ while $l_{ij}^{(i)}$ and $l_{ij}^{(j)}$ are the corresponding segment lengths. The capillary pressure jump between the fluid–fluid interface is expressed by the Young–Laplace equation,

(2.7)$$\begin{eqnarray}\unicode[STIX]{x0394}P_{ij}^{c}=\unicode[STIX]{x1D6FE}_{sf}\cos \unicode[STIX]{x1D703}\left(\frac{1}{R_{1}}+\frac{1}{R_{2}}\right),\end{eqnarray}$$

where $R_{1}$ and $R_{2}$ are the radii of curvature of the interface in any two orthogonal planes, $\unicode[STIX]{x1D703}$ is the contact angle and $\unicode[STIX]{x1D6FE}_{sf}$ is the surface tension. For a given distribution of fluids, the phase conductances are calculated from (2.6). The pressure fields can be computed by substituting (2.5) into (2.4), which yields the following linear system of equations:

(2.8)$$\begin{eqnarray}[\unicode[STIX]{x1D642}]\{\boldsymbol{P}\}=[\boldsymbol{B}],\end{eqnarray}$$

where $\unicode[STIX]{x1D642}$ is the matrix of the conductivities, $\boldsymbol{P}$ is the vector of unknown nodal pressures and $\boldsymbol{B}$ is a vector formed from boundary pressures and capillary pressures. This equation can be solved using the biconjugate gradient solution method. The updated pressure field is then used in (2.5) to compute fluxes through the pore throats. Simulations can also be carried out with constant injection rate boundary conditions (Aker et al. Reference Aker, Jørgen Måløy, Hansen and Batrouni1998), which we used in this work. After solving for the volumetric fluxes of all capillary elements, the current invading phase saturation in each pore body is used to determine the minimum filling time for each of the pores; this is set as the time step size $\unicode[STIX]{x0394}t$. This time step corresponds to the shortest time required to fill one non-wetting phase-filled pore. Then, equation (2.1) is used to update the saturation in each pore body.

2.2 Finite element method

In this paper, the solid response is quantified using finite element formulation. We assume that the solid undergoes large deformations and, as such, we model the solid domain as nonlinear elastic. For a spatial solid domain $\unicode[STIX]{x1D6FA}_{s}^{0}\subset \mathbb{R}^{d}$ bounded by $\unicode[STIX]{x1D6E4}$, the continuum equations can be described in total Lagrangian form as

(2.9a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}\ddot{\boldsymbol{u}}=\unicode[STIX]{x1D735}_{0}\boldsymbol{\cdot }(\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{F}^{\text{T}})+\unicode[STIX]{x1D70C}_{0}\boldsymbol{b}\quad \text{in }\unicode[STIX]{x1D6FA}_{s}^{0}, & \displaystyle\end{eqnarray}$$

(2.9b)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\overline{\overline{\boldsymbol{u}}}\quad \text{on }\unicode[STIX]{x1D6E4}_{D}^{0}, & \displaystyle\end{eqnarray}$$

(2.9c)$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{F}^{\text{T}})\boldsymbol{\cdot }\boldsymbol{n}=\overline{\overline{\boldsymbol{T}}}\quad \text{on }\unicode[STIX]{x1D6E4}_{N}^{0}, & \displaystyle\end{eqnarray}$$

$0$

$\boldsymbol{u}$

$\boldsymbol{F}$

$\boldsymbol{F}=\unicode[STIX]{x1D735}_{0}\boldsymbol{u}$

$\overset{}{\boldsymbol{u}}$

$\overline{\overline{\boldsymbol{T}}}$

$\boldsymbol{n}$

$\unicode[STIX]{x1D70C}_{0}$

$\boldsymbol{b}$

$\unicode[STIX]{x1D64E}$

$\unicode[STIX]{x1D640}$

$\unicode[STIX]{x1D64E}$

$\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{F}^{\text{T}}$

$B_{oh}=\bigcup _{e=1}^{E}\unicode[STIX]{x1D6FA}_{s}^{e,\ast }$

$\unicode[STIX]{x1D6FA}_{0}^{e}$

$\unicode[STIX]{x1D6E4}_{0}^{e}=\unicode[STIX]{x1D6E4}_{D}^{0}\cup \unicode[STIX]{x1D6E4}_{N}^{0}$

$\boldsymbol{u}_{h}$

$\unicode[STIX]{x1D64B}_{h}$

(2.10)$$\begin{eqnarray}\displaystyle & \displaystyle X_{h}^{k}=\{\boldsymbol{u}_{h}\in L^{2}(B_{oh})|[\boldsymbol{u}_{h}|_{\unicode[STIX]{x1D6FA}_{0}^{e}}\in \mathbb{M}^{k}(\unicode[STIX]{x1D6FA}_{0}^{e})\forall \unicode[STIX]{x1D6FA}_{0}^{e}\in B_{oh}]\}\subset X^{f}(B_{oh})=\mathop{\prod }_{e}(H^{1}(\unicode[STIX]{x1D6FA}_{0}^{e})),\qquad & \displaystyle\end{eqnarray}$$

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle S_{h}^{k}=\{\unicode[STIX]{x1D64B}_{h}\in [L^{2}(B_{oh})]^{2}|[\unicode[STIX]{x1D64B}_{h}|_{\unicode[STIX]{x1D6FA}_{0}^{e}}\in \mathbb{M}^{k}(\unicode[STIX]{x1D6FA}_{0}^{e})^{2}\forall \unicode[STIX]{x1D6FA}_{0}^{e}\in B_{oh}]\}\subset S^{f}(B_{oh})=\mathop{\prod }_{e}(H^{1}(\unicode[STIX]{x1D6FA}_{0}^{e}))^{2},\qquad & \displaystyle\end{eqnarray}$$

where $\mathbb{M}^{k}$ is the set of polynomial functions up to degree $k\geqslant 1$. Let $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}$ be an arbitrary test function defined in the space

(2.12)$$\begin{eqnarray}\displaystyle X_{hc}^{k} & = & \displaystyle \{\!\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}\in X_{h}^{k}|[\!\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}=0\forall \boldsymbol{X}\in \unicode[STIX]{x2202}B_{oh},\forall t\in T\text{ and }\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}(t_{0})=0\forall \boldsymbol{X}\in B_{oh}\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}(t_{f})\nonumber\\ \displaystyle & & \displaystyle \qquad =0\forall \boldsymbol{X}\in B_{oh}\! ]\!\}.\end{eqnarray}$$

The weak formulation of the expression is carried out by integration in the reference coordinate system by finding $\boldsymbol{u}_{h}\in X_{h}^{k}$ and $\unicode[STIX]{x1D64B}_{h}\in S_{h}^{k}$, such that

(2.13)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{e}\int _{\unicode[STIX]{x1D6FA}_{0}^{e}}(\unicode[STIX]{x1D70C}_{0}\ddot{\boldsymbol{u}}_{h}-\unicode[STIX]{x1D735}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{h})\cdot \unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}\,\text{d}V=\mathop{\sum }_{e}\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{0}^{e}}\overline{\overline{\boldsymbol{T}}}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}\,\text{d}S+\mathop{\sum }_{e}\int _{\unicode[STIX]{x1D6FA}_{0}^{e}}\unicode[STIX]{x1D70C}_{0}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}\,\text{d}V\forall \unicode[STIX]{x1D6FF}\boldsymbol{u}_{h}\nonumber\\ \displaystyle & & \displaystyle \qquad \in X_{hc}^{k},\forall t\in T.\end{eqnarray}$$

Based on the nonlinear dynamic behaviour of the solid, the following discrete nonlinear system of linear equations can be obtained:

(2.14)$$\begin{eqnarray}\unicode[STIX]{x1D648}\ddot{\boldsymbol{u}}(t)+\unicode[STIX]{x1D63E}\dot{\boldsymbol{u}}(t)+\boldsymbol{f}_{int}(t)=\boldsymbol{f}_{ext}(t),\end{eqnarray}$$

where $\unicode[STIX]{x1D648}$ is the $3N\times 3N$ mass matrix, $\unicode[STIX]{x1D63E}$ is the damping matrix, $\boldsymbol{f}_{int}(t)$ and $\boldsymbol{f}_{ext}(t)$ are the internal force vector and the external force vector from the fluid in $3N$ dimensions, respectively. Since the solution is solved using Lagrangian meshes, $\unicode[STIX]{x1D648}$ is time invariant. Vibration effects are not considered in this work, so the damping coefficient is zero. The mass matrix, internal and external forces follow as

(2.15)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}_{IJ}=\int _{\unicode[STIX]{x1D6FA}_{s}}\unicode[STIX]{x1D70C}_{0}N_{I}N_{J}\,\text{d}V, & \displaystyle\end{eqnarray}$$

(2.16)$$\begin{eqnarray}\displaystyle & \displaystyle (\boldsymbol{f}_{int,I})^{\text{T}}=\int _{\unicode[STIX]{x1D6FA}_{s}}\frac{\unicode[STIX]{x2202}N_{I}}{\unicode[STIX]{x2202}\boldsymbol{X}}\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{F}^{\text{T}}\,\text{d}V, & \displaystyle\end{eqnarray}$$

(2.17)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{ext,I}=\int _{\unicode[STIX]{x1D6FA}_{s}}\unicode[STIX]{x1D70C}_{0}\boldsymbol{b}N_{I}\,\text{d}V+\int _{\unicode[STIX]{x1D6E4}_{N}^{0}}\overline{\overline{\boldsymbol{T}}}\boldsymbol{\cdot }N_{I}\,\text{d}S, & \displaystyle\end{eqnarray}$$

where $N_{I}$ is the shape function corresponding to node $I$. The internal and external nodal forces are functions of nodal displacement and time. The external loads are prescribed as a function of time incrementally. Therefore, assembling the fluid load requires the computation of the time-varying force exerted by the encompassing fluid on the solid across the interface.

2.2.1 Time discontinuous Galerkin method

To solve for the solid response, we apply the time-discontinuous Galerkin method (TDG) (Mancuso & Ubertini Reference Mancuso and Ubertini2003; Noels & Radovitzky Reference Noels and Radovitzky2008). A similar procedure has been used by De Rosis et al. (Reference De Rosis, Falcucci, Porfiri, Ubertini and Ubertini2014) for finite element method (FEM) coupling with the lattice Boltzmann method (LBM). However, their coupling methodology was applied only for linear dynamics. The TDG they applied can be easily adapted for nonlinear analysis, such as for Eulerian formulation (for example, see Denoël & Detournay (Reference Denoël and Detournay2011)). The TDG method considered in this work uses piecewise linear time interpolation. It is based on prediction and correction schemes, which can provide the requisite accuracy in a few iterations. The implementation is based on the Nørsett algorithm, which is optimal from the computational standpoint among higher-order unconditionally stable methods (Mancuso & Ubertini Reference Mancuso and Ubertini2002).

The integration is accomplished through the incremental solution procedure, in which the time interval of interest ($T$) is discretized into $n^{f}$ time steps such that $T=\bigcup _{n=0}^{n=n^{f}-1}[t^{n},t^{n+1}]$ and $\unicode[STIX]{x0394}t=t^{n+1}-t^{n}$ is the time step. The scheme treats the displacement and velocity vectors as independent variables and uses piecewise linear time interpolants. The iterative scheme is summarized by the following equation:

(2.18)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D648}^{\ast }\overline{\boldsymbol{v}}_{k}^{l+1}(t)=\boldsymbol{F}_{k}(t)-\frac{\unicode[STIX]{x2202}\boldsymbol{f}_{int}}{\unicode[STIX]{x2202}\boldsymbol{u}}\tilde{\boldsymbol{u}}_{k}^{l+1}(t),\\ \displaystyle \overline{\boldsymbol{u}}_{k}^{l+1}(t)=\tilde{\boldsymbol{v}}_{k}^{l+1}(t)+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}t\overline{\boldsymbol{v}}_{k}^{l+1}(t),\end{array},\quad k=0,1;l=0,1\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x2202}\boldsymbol{f}_{int}/\unicode[STIX]{x2202}\boldsymbol{u}=\unicode[STIX]{x1D646}^{int}$ is the Jacobian matrix of the internal forces or the tangent stiffness matrix, and the barred terms represent time derivatives. The predictor displacements $\tilde{\boldsymbol{u}}_{k}$ and velocities $\tilde{\boldsymbol{v}}_{k}$ are given as

(2.19a)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{v}}_{0}^{l+1}=\boldsymbol{v}_{\ast }^{t}+\unicode[STIX]{x1D706}^{\ast }\unicode[STIX]{x0394}t(\overline{\boldsymbol{v}}_{0}^{l}-\overline{\boldsymbol{v}}_{1}^{l}), & \displaystyle\end{eqnarray}$$

(2.19b)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{u}}_{0}^{l+1}=\boldsymbol{u}_{\ast }^{t}+\unicode[STIX]{x1D706}^{\ast }\unicode[STIX]{x0394}t(\overline{\boldsymbol{u}}_{0}^{l}-\overline{\boldsymbol{u}}_{1}^{l})+\unicode[STIX]{x1D6FC}^{\ast }\unicode[STIX]{x0394}t\tilde{\boldsymbol{v}}_{0}^{l+1}, & \displaystyle\end{eqnarray}$$

(2.19c)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{v}}_{1}^{l+1}=\boldsymbol{v}_{\ast }^{t}+(1-\unicode[STIX]{x1D6FC}^{\ast })\unicode[STIX]{x0394}t\overline{\boldsymbol{v}}_{0}^{l+1}, & \displaystyle\end{eqnarray}$$

(2.19d)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{u}}_{1}^{l+1}=\boldsymbol{u}_{\ast }^{t}+(1-\unicode[STIX]{x1D6FC}^{\ast })\unicode[STIX]{x0394}t\overline{\boldsymbol{u}}_{0}^{l+1}+\unicode[STIX]{x1D6FC}^{\ast }\unicode[STIX]{x0394}t\tilde{\boldsymbol{v}}_{1}^{l+1}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}_{\ast }^{t}$ and $\boldsymbol{v}_{\ast }^{t}$ are the vector of unknown displacements and its time derivative at time $t$ of the FE solution. Here, $\unicode[STIX]{x1D648}^{\ast }$ is the effective mass matrix given by

(2.20)$$\begin{eqnarray}\unicode[STIX]{x1D648}^{\ast }=\unicode[STIX]{x1D648}+{\unicode[STIX]{x1D6FC}^{\ast }}^{2}\unicode[STIX]{x0394}t^{2}\frac{\unicode[STIX]{x2202}\boldsymbol{f}_{int}}{\unicode[STIX]{x2202}\boldsymbol{u}}.\end{eqnarray}$$

It should be noted that this effective mass matrix is equally symmetric and positive-definite, just as $\unicode[STIX]{x1D648}$ and $\unicode[STIX]{x1D646}^{int}$. The effective load $\boldsymbol{F}_{k}$ (in which the fluid loads would be transmitted) is obtained within two iterations. Based on the on left Radau abscissas, the effective loads are defined as

(2.21)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{0}=\left(\frac{3\sqrt{2}}{4}-\frac{1}{2}\right)\boldsymbol{f}_{ext,0}+(6-3\sqrt{2})\boldsymbol{f}_{ext,2/3}, & \displaystyle\end{eqnarray}$$

(2.22)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{1}={\textstyle \frac{3}{2}}\,\boldsymbol{f}_{ext,2/3}-{\textstyle \frac{1}{2}}\boldsymbol{f}_{ext,0}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{f}_{ext,0}$ is the value of $\boldsymbol{f}_{ext}$ at time $t$ and $\boldsymbol{f}_{ext,2/3}$ is the value at time $t+2/3\unicode[STIX]{x0394}t$. Coefficients $\unicode[STIX]{x1D706}^{\ast }$ and $\unicode[STIX]{x1D6FC}^{\ast }$ are given as follows:

(2.23a,b)$$\begin{eqnarray}\unicode[STIX]{x1D706}^{\ast }=\sqrt{2}-\frac{4}{3},\quad \unicode[STIX]{x1D6FC}^{\ast }=1-\frac{\sqrt{2}}{2}.\end{eqnarray}$$

The model is iterated twice to achieve third-order accuracy. To complete the TDG algorithm, the following equations are solved at the end of the current time step:

(2.24a,b)$$\begin{eqnarray}\boldsymbol{v}_{\ast }^{n+1}=\overline{\boldsymbol{u}}_{1}^{l},\quad \boldsymbol{u}_{\ast }^{n+1}=\boldsymbol{u}_{\ast }^{n}+\unicode[STIX]{x0394}t[(1-\unicode[STIX]{x1D6FC}^{\ast })\overline{\boldsymbol{u}}_{0}^{l}+\unicode[STIX]{x1D6FC}^{\ast }\overline{\boldsymbol{u}}_{1}^{l}].\end{eqnarray}$$

2.3 Coupled model

As stated previously, this work aims to develop a three-dimensional fluid–solid model incorporating FE and DPN. Both solution domains are illustrated in figure 3. Pore fluids exert forces on the solid walls, resulting in deformation and eventual modification of force balance. Additionally, the deformation of the solid results in the rearrangement of its structure which can alter fluid properties and the network topology. Thus, the need to update the network regularly during the simulation to track the deformation of the solid. To achieve the solution of such a system, the fluidic and structural computations are performed sequentially. Therefore, at each time step, we began by solving the DPN equation, from which we obtain the pressure at each node and throat. The DPN solver then feeds the hydrodynamic load to the solid at each time-exchange step. The FE solver computes the solid displacement $\boldsymbol{u}$, which is then used to determine the new structure of the DPN.

Figure 3.

Schematic of the superimposed domains (fluid and solid). (a) The voxelized solid phase along with the PN and (b) discretized finite element domain and PN.

More precisely, the updated set of displacements are used to compute the new position of the solid boundary. The structural computation updates the position of the structural surface. Thus, an updated fluid domain is needed to accommodate the new interface location. We achieved this by performing a re-extraction process for the network. This provides a new value for the node radius $\{r_{i}^{\prime }\}$ and volume $\{V_{i}^{\prime }\}$, as well as a new throat length $\{l_{ij}^{\prime }\}$ and radius $\{r_{ij}^{\prime }\}$. An iterative process is also required at each time step to ensure that the coupled system converges to a steady solution before advancing to the next time step. This is done to reduce the magnitude of errors that could occur during the network re-extraction process. We examine this further in § 3. Another important aspect of the fluid–solid coupling procedure is the hydrodynamic load calculation. The force on the pore body is calculated with the following equation:

(2.25)$$\begin{eqnarray}\boldsymbol{F}_{i}=\int _{\unicode[STIX]{x1D6E4}_{i,f}}P_{i}\boldsymbol{n}\,\text{d}s,\quad \text{on }\unicode[STIX]{x1D6E4}_{i,f}\end{eqnarray}$$

where $\boldsymbol{F}_{i}$ is the pressure force applied by the $i$th pore with pore pressure $P_{i}$. The boundary of the $i$th pore is $\unicode[STIX]{x1D6E4}_{i,f}$ and $\boldsymbol{n}$ is the unit vector pointing from the centroid of the pore to the nearest solid boundary. Since we are using spherical DPN nodes, we assume the force is the same across the same pore-body node. The fluid load from the pore-body surface along the normal direction is imparted to the solid nodes. One of the central features in this coupling implementation is the exchange of information across the fluid and solid domains. As demonstrated in figure 3(a), the transmission of the fluid load to the solid would require a secondary domain for interfacing both subsystems since the adjoining pore bodies and solid do not match exactly.

2.3.1 Convex hull

In this work, the fluid loads are transmitted to solid nodal points that border an arbitrary volume surrounding the nearest pore body and throat. This initially requires test point values (spatial surface data) generated on the pore network, followed by a sequential neighbourhood search for FE nodal points within a prescribed region of the test point data. Several volumetric shapes could be used to enclose the points (e.g. the smallest cube or sphere). However, we use a convex hull, which is defined as the smallest convex set enclosing the points. The convex hull helps reduce the amount of empty space and saves memory. Calculation of convex hulls is a well-studied problem in computational geometry and this method has diverse applications in other fields, such as cluster analysis, collision detection, image processing, statistics, sphere packing and point location (Barber, Dobkin & Huhdanpaa Reference Barber, Dobkin and Huhdanpaa1996). It has also been applied in the natural element method (Sukumar, Moran & Belytschko Reference Sukumar, Moran and Belytschko1998), which relies on Delaunay triangulations in $\mathbb{R}^{n}$, computed from convex hull in $\mathbb{R}^{n+1}$, to construct its interpolants. Convex hull has also been applied for modelling geophysical phenomena and complex FSI problems (Sukumar et al. Reference Sukumar, Moran and Belytschko1998). We define a convex hull as follows: for a set of randomly generated points $\unicode[STIX]{x1D6F6}^{k}=\{\unicode[STIX]{x1D6FE}_{\mathbf{1}},\unicode[STIX]{x1D6FE}_{\mathbf{2}},\ldots ,\unicode[STIX]{x1D6FE}_{\boldsymbol{N}}\}$ on the $k$th pore in $\mathbb{R}^{n}$, the convex hull $C$ of the points is expressed mathematically in the form

(2.26)$$\begin{eqnarray}C=\left\{\mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D706}_{i}\unicode[STIX]{x1D6FE}_{i}:\unicode[STIX]{x1D706}_{i}\geqslant 0\forall i\in \{1,\ldots ,n\}\wedge \mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D706}_{i}=1\right\}.\end{eqnarray}$$

Equation (2.26) describes the convex combinations of points in set $\unicode[STIX]{x1D6F6}$. For an $n$-dimensional convex hull set of nodal points, we apply the Quickhull (known as Qhull) algorithm developed by Barber et al. (Reference Barber, Dobkin and Huhdanpaa1996) to calculate the convex hull of sets of the multidimensional points of that hull. The hull is created based on random points (test points) generated on the surface of pores. Now for such a set of points, if a distinctive portion of the finite element nodal points ($\boldsymbol{n}=\{n_{1},n_{2},\ldots ,n_{M}\}$) and $\unicode[STIX]{x1D6F6}$ are found within the same space, they are labelled unity, while points outside are labelled zero. This process is graphically shown in figure 4. The forces $\boldsymbol{F}_{i}$ are then transmitted to nodes labelled unity while others outside this space, labelled zero, are left out.

Thus, at each time step, the coupling algorithm proposed in this work has the following sequence of operations: