This paper presents an immersed phase field model designed to predict the fracture-induced flow due to brittle fracture in vuggy porous media. Due to the multiscale nature of pores in the vuggy porous material, crack growth may connect previously isolated pores, which leads to flow conduits. This mechanism has important implications for many applications such as disposal of carbon dioxide and radioactive materials and hydraulic fracture and mining. To understand the detailed microporomechanics that causes the fracture-induced flow, we introduce a new phase field fracture framework where the phase field is not only used as an indicator function for damage of the solid skeleton but also used as an indicator of the pore space. By coupling the Stokes equation that governs the fluid transport in the voids, cavities, and cracks and Darcy’s flow in the deformable porous media, our proposed model enables us to capture the fluid–solid interaction of the pore fluid and solid constituents during crack growth. Numerical experiments are conducted to analyze how the presence of cavities affects the accuracy of predictions based on the homogenized effective medium during crack growth.

Geomaterials such as carbonate rocks, sandstone, or limestone often contain geometrical features such as cracks, joints, vugs, or cavities. When the defects are partially or fully saturated with the pore fluid, the geometry of the features may affect effective stiffness, permeability, water retention characteristics, and drained or undrained shear strength of the material (Juanes et al., 2006; Sun et al., 2011b; Kang et al., 2016; Suh et al., 2017; Selvadurai et al., 2017; Wang and Sun, 2017; and Sun and Wong, 2018). Furthermore, brittle fracture in materials that possess geometrical features may lead to the pore fluid in cavities migrating into the flow channels and causing flow conduits that lead to often undesirable outcomes. Modeling geometrical features in porous media is, thus, highly important and, at the same time, a challenging subject for the hydromechanically coupled analysis in geomechanics problems such as hydrocarbon resource recovery or development of enhanced geothermal energy reservoirs (Paterson and Fermigier, 1997; Class et al., 2002; Rutqvist et al., 2007; Grant, 2013; Wagner et al., 2015; Heider and Markert, 2017; and Suh and Sun, 2020).

One possible modeling choice is to consider a fictitious effective medium at a scale where representative elementary volume exists. In this case, the geometrical features of the material are not explicitly modeled, but the influences of these geometrical features are incorporated in the constitutive relations by treating defects as a different pore system that interacts with the matrix pores (Choo and Borja, 2015; Choo et al., 2016; Liu and Abousleiman, 2017; and Wang and Sun, 2018). The upshot of the multi-porosity and multi-permeability models is mainly the simple numerical treatment, since there is no need for complex meshing techniques or embedded strong discontinuities, and the computational efficiency compared to pore-scale models that require extremely large domains in order to reproduce hydromechanical behavior at large scales (Ghaboussi and Barbosa, 1990; Spaid and Phelan, 1997; Blunt et al., 2002; Tang et al., 2005; Arson and Pereira, 2013; Pereira and Arson, 2013; and Suh and Yun, 2018). However, the drawback of this approach is that the homogenized effective medium may not sufficiently represent the microstructural details. This makes the identification of material parameters more complicated since the effective permeability of multiple interacting systems is not isotropic and the constitutive law for fluid mass exchange inherently depends on the microstructure.

Another common alternative to model the interaction between cavities and crack growth is to conduct simulations via a fracture network model (Ozkan et al., 2010; Leung and Zimmerman, 2012; Fu et al., 2013; and Hyman et al., 2015). However, the obvious drawback is that the fracture in those models must be either a straight line (in the two-dimensional case) or a plane (in the three-dimensional case), and hence, the geometrical effect on the porous media cannot be captured precisely.

In this research, we introduce a phase field framework that allows us to enable a unified treatment to simulate the evolving geometry of cracks and cavities. By introducing the phase field as a unified representation of the void space that is not suitable to be treated as an effective medium, we introduce a framework that enables us to analyze how crack propagation in vuggy porous media may affect the flow mechanism differently than that in porous media with pores well distributed in the host matrix. Our result indicates that the interaction between the propagating cracks of the cavities is important for capturing the hydromechanical responses of the porous media, and that existing effective medium approach that characterizes the pore space with a single hydraulic model such as the cubic law and Kozeny–Carman model may not be sufficient to capture the cavity–crack–host matrix interactions.

As for notations and symbols, bold-faced and blackboard bold-faced letters denote tensors (including vectors that are rank-one tensors); the symbol “·” denotes a single contraction of adjacent indices of two tensors (e.g., a · b = aibi or c · d = cijdjk); the symbol “:” denotes a double contraction of adjacent indices of tensors of rank two or higher (e.g., C:ε = Cijklεkl); and the symbol “⊗” denotes a juxtaposition of two vectors (e.g., ab = aibj) or two symmetric second order tensors [e.g., (αβ)ijkl = αijβkl]. We also define identity tensors as I = δij and I=δikδjl, where δij is the Kronecker delta. As for sign conventions, unless specified, the directions of the tensile stress and dilative pressure are considered positive.

We consider a fully saturated Biot–Stokes system (Fig. 1) that consists of two regions (intact porous matrix BD and cracks or cavities BS) separated by the sharp interface Γ*, where we assume that the solid phase in BD forms a deformable porous matrix while solid particles in BS are in suspension. In this case, both the solid and fluid phases coexist in both regions. By considering our material of interest as a multi-phase continuum, we utilize the effective stress principle for the intact porous matrix where the fluid flow is modeled with Darcy’s law, while the motion of the solid–fluid mixture is modeled by the Stokes equation (Li et al., 2018). Two distinct regions are then coupled by properly imposing three transmissibility conditions at the interface. The model problem with the sharp interface will be later on extended to a diffuse Biot–Stokes model by introducing the phase field in Sec. III.

FIG. 1.

Schematic representation of the Biot–Stokes system that possesses a sharp interface Γ*.

FIG. 1.

Schematic representation of the Biot–Stokes system that possesses a sharp interface Γ*.

Close modal

Although the Biot–Stokes system only contains two immiscible solid and fluid phases, for mathematical convenience, we idealize the material of interest as a three-phase continuum where each constituent [i.e., solid (s), pore fluid (fD), and free fluid (fS)] occupies a fraction of volume at the same material point. By letting dV = dVs + dVf denote the representative elementary volume of the material, we define the volume fractions of the constituents as

ϕαdVαdV  ;  α=s,f,
(1)

where the index s refers to the solid phase and f indicates the fluid phase. Since the sharp interface separates our system of interest into two regions, the volume fraction of the pore and free fluids can be expressed as

ϕf=(1HΓ*)dVfDdVϕfD+HΓ*dVfSdVϕfS,
(2)

where HΓ* is the Heaviside function that satisfies

HΓ*=0inBD1inBS.
(3)

In addition, by letting ρs and ρf denote the intrinsic mass densities of the solid and fluid, respectively, the partial mass densities of each constituent (ρα, where α = s, fD, fS) are given by

ρsϕsρs  ;  ρfDϕfDρf  ;  ρfSϕfSρf  ;  ρρs+ρfD+ρfS,
(4)

where ρ is the mass density of the entire system. In this study, we assume that both the solid and fluid phases are incompressible, so that intrinsic mass densities ρs and ρf are considered constants.

This section briefly reviews the balance principles, constitutive laws in the bulk volume of a porous medium (Sec. II B 1), the region where the solid–fluid mixture flows freely (Sec. II B 2), and the sharp interface between two regions (Sec. II B 3).

1. Conservation laws for an intact porous matrix

For the region where the solid forms an intact porous matrix, we adopt the effective stress principle (Lade and De Boer, 1997; Borja, 2006) so that the external loading imposed on the matrix is assumed to be carried by both the solid skeleton and the pore fluid. In this case, the region BD is governed by the following system of equations (Borja and Alarcón, 1995; White and Borja, 2008; and Sun et al., 2013):

(σBpfDI)+ρg=0   inBD,
(5)
vs+wfD+1Mṗ=0   inBD,
(6)

where σ′ is the effective stress, B = 1 − K/Ks is Biot’s coefficient, M is Biot’s modulus, pfD is the pore pressure, g is the gravitational acceleration, vα is the intrinsic velocity of the constituent α, and wfD=ϕfD(vfDvs) is the Eulerian relative flow vector of the pore fluid (i.e., Darcy’s velocity). Here, we assume that B ≈ 1 and 1/M ≈ 0 to simplify the formulation. Note that Biot’s coefficients of many sandstone and shale specimens are often less than one, whereas it is more reasonable to assume Biot’s coefficient equal to 1 for granite (e.g., Westerly granite) (Berryman and Wang, 1995; Zimmerman, 2000; and Jaeger et al., 2009). In either case, the damage of the solid skeleton may reduce the elastic bulk modulus of the solid skeleton. Therefore, both Biot’s coefficient and modulus may evolve according to the solid deformation. This nonlinear effect is not considered in this study but will be considered in the future. We also assume that the behavior of the intact matrix in BD is linear and isotropic elastic, and hence, only two independent elastic moduli are needed to capture the elastic response. The constitutive relation for the solid skeleton can, therefore, be written as follows:

σ0=λtr(ε)I+2με  inBD,
(7)

where σ0 indicates the effective stress of the undamaged matrix. The actual and undamaged effective stresses are related by a degradation function, which will later be discussed in Sec. III B. Furthermore, ε=(us+usT)/2 is the infinitesimal solid strain tensor that depends on the solid displacement us, and parameters λ and μ are the Lamé constants. For the constitutive equation that describes the laminar pore fluid flow in BD, we use the generalized Darcy’s law that linearly relates the relative velocity wfD and pore pressure gradient pfD, i.e.,

wfD=kμf(pfDρfg)  inBD,
(8)

where μf is the dynamic viscosity of the pure fluid phase and k is the effective permeability of the porous matrix. Additionally, in order to incorporate the effect of deformation of the matrix on the porous medium flow (Mauran et al., 2001; Schutjens et al., 2004), this study adopts the Kozeny–Carman equation to empirically capture the porosity–permeability relation (Chapuis and Aubertin, 2003; Costa, 2006; and Wang and Sun, 2017). Note that the Kozeny–Carman equation is often considered a rough approximation of the porosity–permeability relation. A more precise prediction of permeability may require new geometrical attributes such as tortuosity (Sun et al.,2011b; 2011a), formation factor (Worthington, 1993; Sun and Wong, 2018), and percolation threshold (Mavko and Nur, 1997). This extension is out of the scope of this study but will be considered in the future work.

Recall Sec. II A that ϕs+ϕfD=1 in BD. Then, by letting ϕϕfD denote the porosity of the matrix, the Kozeny–Carman equation reads

k=k0(1ϕ0)2ϕ03ϕ3(1ϕ)2  inBD,
(9)

where k0 and ϕ0 denote the reference permeability and porosity, respectively.

2. Conservation laws for solid–fluid mixture

This study attempts to model suspension flow in BS, where mass and linear momentum balances for both solid and fluid phases should be satisfied. We, therefore, write the governing balance equations for BS as

σfS+ρfSg=0   inBS,
(10)
σs+ρsg=0   inBS,
(11)
vs+wfS=0   inBS,
(12)

where σα is the Cauchy stress tensor of the α constituent, and the relative flow vector of the free fluid can be defined as wfS=ϕfS(vfSvs). By assuming that the free fluid resides in BS with low Reynolds number (i.e., Re ≪ 1), we adopt a simplified version of the Navier–Stokes model, i.e., the Stokes equation. The Stokes model for the steady-state motion of an incompressible fluid yields the following relationship for the free fluid stress tensor σfS:

σfS=pfSI+μeff(vfS+vfST)  inBS,
(13)

where pfS is the free fluid pressure and μeff is the effective viscosity of the solid–fluid mixture (Mooney, 1951; Cheng and Law, 2003), i.e.,

μeff=μfexp2.5c1c/cmax  inBS,
(14)

where c1ϕfS indicates the solid particle concentration and cmax denotes its upper bound. Again, notice that we introduce only one solid constituent for the entire system since it is convenient for us to later on impose interface conditions and further adopt the phase field fracture model that simulates the evolving interface. This approach may not be suitable for modeling complete suspension flow where vs=vfS. However, we assume that solid particles in BS follow the same constitutive relations as the free fluid in order to replicate the suspension flow as close as possible, i.e.,

σs=pfSI+μeff(vs+vsT)  inBS.
(15)

3. Conservation laws for the sharp interface between intact matrix and solid–fluid mixture

In order to properly model the interaction between the porous matrix (BD) and the vugs or cavities (BS), complete mass conservation and force equilibrium for the entire system should be satisfied. Since we have two different constituents for the same type of fluid (fD and fS) while considering only one solid constituent (s), coupling two subsystems, thus, requires the enforcement of fluid transmissibility conditions at the sharp interface Γ* that models the coupled Stokes–Darcy flow (Arbogast and Lehr, 2006; Arbogast and Brunson, 2007; Badia et al., 2009; Wu and Mirbod, 2018; and Bergkamp et al., 2020).

The first interface condition is the fluid continuity that ensures mass conservation. Since we assume that the fluid phase is incompressible, the interfacial fluid fluxes for each subsystem (MfD* and MfS*) can be expressed as follows:

MfD*=Γ*wfDnD*mfD*dΓ  ;  MfS*=Γ*wfSnS*mfS*dΓ,
(16)

where nD* and nS* denote the outward-oriented normal vectors from BD and BS, respectively. From Eq. (16), mass continuity (MfD*+MfS*=mfD*+mfS*=0) yields the following transmissibility condition:

wfDnD*+wfSnS*=(wfSwfD)n*=0  onΓ*,
(17)

where we take n*=nS*=nD* for notational convenience (Fig. 1). Here, Eq. (16) implies that the normal component of the fluid velocities (wfS and wfD) should be identical in order to guarantee that the exchange of fluid mass between BS and BD is conservative.

The second condition is the force equilibrium at the interface Γ*. From each subsystem, total forces acting on the interface (fD* and fS) may be written as

fD*=Γ*pfDn*tfD*dΓ  ;  fS*=Γ*σfSn*tfS*dΓ,
(18)

where tfD* and tfS* indicate the tractions at the interface. The force equilibrium requires fD*+fS*=tfD*+tfS*=0, implying that the normal and shear components should be balanced at the same time. By decomposing the traction vectors as

ti*=(ti*n*)n*+j=12(ti*mj*)mj*  ;  i=fD,fS,
(19)

where m1* and m2* are the interfacial tangent vectors, we get two more transmissibility conditions that describe normal and shear force equilibrium, respectively,

tfS*n*+pfD=0   onΓ*,
(20)
tfS*mj*+μfαSDk(wfSwfD)mj*=0   onΓ*.
(21)

Equation (21) is the Beavers–Joseph–Saffman condition (Beavers and Joseph, 1967; Saffman, 1971; Layton et al., 2002; and Arbogast and Brunson, 2007). This idealized condition relates the slip velocity and the shear stress through the dimensionless slippage coefficient αSD, which depends on the microstructural attributes of the interfaces, such as surface roughness and irregular patterns, as well as the flow velocity (Beavers and Joseph, 1967; Terzis et al., 2019; and Guo et al., 2020). The validity and limitations of the Beavers–Joseph–Saffman condition are documented in a number of literature studies such as Auriault (2010), Mikelic and Jäger (2000), and Monchiet et al. (2019) and will not be repeated here. Possible extensions of the interface conditions to turbulent and multiphase flows are an active research area that is clearly out of the scope of this study but will be considered in the future.

This section introduces the mathematical model that uses smooth implicit function, i.e., the phase field, to approximate evolving sharp interfaces due to damage. We first review the general procedure that employs an implicit function to approximate sharp interfaces (Sec. III A) shown in Fig. 2. Since the phase field is a smooth representation of the Heaviside function, we derive the corresponding mathematical model that approximates interfacial transmissibility conditions suitable for the diffuse representation of the interface. To capture crack growth according to Griffith’s theory, we adopt the classical variational fracture model to allow for crack growth represented by the evolution of the phase field defined over the spatial domain (Sec. III B). These techniques are then applied to the derivation shown in Sec. III C in which a mathematical model to capture the hydromechanical coupling of pore fluid flows in both the host matrix and evolving interfaces in brittle porous media is discussed. The resultant model does not require locally defined enrichment function or remeshing and can be implemented in a standard finite element or finite element/volume solver.

FIG. 2.

Diffuse representation of the interface where an exemplary 1D domain consists of BS in x/L ∈ [0.4, 0.6] sandwiched between undamaged porous matrices BD.

FIG. 2.

Diffuse representation of the interface where an exemplary 1D domain consists of BS in x/L ∈ [0.4, 0.6] sandwiched between undamaged porous matrices BD.

Close modal

This study employs a diffuse approximation for the sharp interface Γ* by introducing a phase field variable d ∈ [0, 1] that varies smoothly from 0 in BD to 1 in BS. Specifically, we approximate the interfacial area AΓ* as AΓd*, which can be expressed in terms of volume integration of the surface density functional Γd*(d,d) over B=BDBS¯ (Miehe et al., 2010; Borden et al., 2012; Suh and Sun, 2019; and Suh et al., 2020),

AΓ*AΓd*=BΓd*(d,d)dV.
(22)

Here, the size of the diffusive zone [i.e., transition zone where d ∈ (0, 1)] is controlled by the regularization length scale parameter l* such that AΓd*Γ-converges to AΓ* (Mumford and Shah, 1989), i.e.,

AΓ*=liml*0AΓd*.
(23)

Based on this approach, the phase field d and its gradient ∇d can be regarded as smooth approximations of the Heaviside function HΓ* and the Dirac delta function δΓ*, respectively (Stoter et al., 2017; Suh and Sun, 2020). Therefore, the volume integrals of an arbitrary function G̃ over BD and BS can, respectively, be approximated as

BDG̃dV=BG̃(1HΓ*)dV=liml*0BG̃(1d)dVBG̃(1d)dV,
(24)
BSG̃dV=BG̃HΓ*dV=liml*0BG̃ddVBG̃ddV.
(25)

Similarly, the surface integral of the function G̃ along the sharp interface Γ* can be approximated as

Γ*G̃dΓ=Γ*G̃δΓ*dV=liml*0BG̃ddVBG̃ddV,
(26)

and we also approximate the normal vector n* as

n*dd.
(27)

For completeness, this section reviews the phase field model for brittle fracture. We consider the following surface density functional, which is widely used in modeling brittle or quasi-brittle fracture (Bourdin et al., 2008; Miehe et al., 2010; Borden et al., 2012; Bryant and Sun, 2018; and Suh et al., 2020) that possesses quadratic local dissipation function:

Γd*(d,d)=d22l*+l*2(dd).
(28)

At this point, we highlight that the evolution of the phase field (i.e., propagation of cavities or cracks) is a mechanical process driven by the effective stress σ′. In other words, we assume that the solid skeleton is completely damaged in the liquefied zone BS, whereas, in BD, the solid skeleton remains undamaged. We, thus, omit the terms that are unrelated to the deformation and fracture in this section. Having critical energy Gc that is required to create new free surfaces, the potential energy density ψ reads

ψ=g(d)ψe+(ε)+ψe(ε)ψbulk(ε,d)+GcΓd*(d,d),
(29)

where ψbulk(ε, d) is the degrading elastic bulk energy and g(d) = (1 − d)2 is the degradation function that induces energy dissipation. Following Amor et al. (2009), we adopt the additive decomposition scheme that splits the elastic energy ψe into compressive (ψe) and tensile and deviatoric (ψe+) modes, where we only degrade ψe+ in order to avoid crack propagation under compression (Na and Sun, 2018; Bilgen and Weinberg, 2019; and Heider and Sun, 2020), i.e.,

ψe+=12Kεvol+2+μ(εdev:εdev),
(30)
ψe=12Kεvol2,
(31)

where K = λ + 2 μ/3 is the bulk modulus of the porous matrix and ±=(±||)/2 is the Macaulay bracket operator. In this case, the effective stress tensor σ′ can also be decomposed as follows:

σ=g(d)σ0++σ0,
(32)

where σ0±=ψe±/ε is the fictitious undamaged effective stress, in which we previously assumed σ0′ to be linear elastic [Eq. (7)].

Based on the fundamental lemma of calculus of variations, the damage evolution equation can be obtained by seeking the stationary point where the functional derivative of Eq. (29) with respect to d vanishes, i.e.,

ψdψd=0,
(33)

where

ψd=g(d)ψe++Gcl*d  ;  ψd=Gcl*2d.
(34)

Here, the superposed prime denotes the derivative with respect to d and ∇2(•) = ∇ · ∇(•) is the Laplacian operator. Furthermore, by following the treatment used by Miehe et al. (2010), we introduce a history function H that is the pseudo-temporal maximum of the positive energy density (ψe+) in order to ensure the crack irreversibility constraint,

H=maxτ[0,t]ψe+.
(35)

By replacing ψe+ in Eq. (34) with H, Eq. (33) finally yields the following phase field equation that governs the evolution of the interface:

g(d)H+Gcl*(dl*22d)=0  inB.
(36)

Note that we can obtain the diffuse representation of the interface by solving Eq. (36), as shown in Fig. 2.

In this study, we leverage the phase field as an indicator function not only for the location of cracks but also for the location of other defects such as cavities or geometrically complicated voids that do not fit for computational homogenization. This approach may efficiently couple the Stokes flow inside the vugs (BS) that interact with the pore fluid in the intact porous matrix (BD) while both regions are evolving due to crack growth. A major advantage of this work is that free flow inside the fracture is explicitly replicated, and hence, there is no need to introduce permeability enhancement models (e.g., cubic law) (Witherspoon et al., 1980; Konzuk and Kueper, 2004; and Jin and Arson, 2020). This explicit treatment enables the simulations to remain physical even in the situations (e.g., high Reynolds number, rough surface, aperture variation) where the validity of the cubic law is questioned (Witherspoon et al., 1980; Miehe and Mauthe, 2016; Heider and Markert, 2017; Wang and Sun, 2017; Sun et al., 2017; Choo and Sun, 2018; Chukwudozie et al., 2019; and Wang and Sun, 2019).

Verification and experimental validation of the phase field fracture model for brittle solids have been well documented in the literature. For brevity, similar studies are not provided in this paper. Interested readers may refer to, for instance, studies by Nguyen et al. (2016) and Pham et al. (2017).

We present an immersed phase field Biot–Stokes model designed to simulate the coupled hydro-mechanical behaviors of flow of vuggy porous media with evolving fractures in the brittle regime. This section omits the gravitational effects for brevity (i.e., g = 0).

The model problem with the sharp interface (Sec. II) is the system that possesses two distinct boundaries BD and BS that can both be decomposed into Dirichlet (BDu, BDp, BSw, and BSp) and Neumann (BDt, BDq, BSt, and BSq) boundaries satisfying

BD=BDuBDt¯=BDpBDq¯;=BDuBDt=BDpBDq,
(37)
BS=BSwBSt¯=BSpBSq¯;=BSwBSt=BSpBSq,
(38)

where the union of BS and BD is B and the boundary domain follows the same treatment. Here, we capture the transition of the constitutive responses of the solid constituent in the intact and liquefied states through a partition of unity argument in the local constitutive responses. As such, we adopt only one solid constituent, and the balance equations of linear momentum in the sub-domains BD and BS [Eqs. (5) and (11)] are combined into one set of equations over the domain B. The governing equations for the model problem are summarized as follows:

(1d)(σpfDI)+dσs=0   inB,
(39)
vs+wfD=0   inBD,
(40)
σfS+ρfSg=0   inBS,
(41)
vs+wfS=0   inBS,
(42)
g(d)H+Gcl*(dl*22d)=0   inB,
(43)
(wfSwfD)n*=0   onΓ*,
(44)
tfS*n*+pfD=0   onΓ*,
(45)
tfS*mj*+μfαSDk(wfSwfD)mj*=0   onΓ*,
(46)

where the natural and essential boundary conditions are not included for brevity. Following the standard weighted-residual procedure, we multiply Eqs. (39)–(43) with proper weight functions (ηs, ξfD, ηfS, ξfS, and ζ) and integrate over their corresponding domain. The resultant weighted-residual statement reads (Badia et al., 2009; Stoter et al., 2017)

Bηs:(σpfDI)(1d)dV+Bηs:σsddVBDtηst^DdΓ=0,
(47)
BDξfD(u̇s)dVBDξfDwfDdVBDqξfDq^DdΓ+Γ*ξfDwfD(n*)=mfD*dΓ=0,
(48)
BSηfS:σfSdVBStηfSt^SdΓΓ*ηfSσfSn*=tfS*dΓ=0,
(49)
BSξfS(u̇s)dV+BSξfS(wfS)dV=0,
(50)
Bζg(d)H+Gcl*ddV+BζGcl*ddV=0,
(51)

where t^D and q^D are the prescribed traction and flux in the porous matrix, respectively, and t^S is the fluid traction. Then, we directly impose the interfacial transmissibility conditions [Eqs. (44)–(46)] into the field equations [Eqs. (48) and (50)]. Due to the fluid mass continuity (i.e., mfD*+mfS*=0), the fourth term on the left-hand side of Eq. (48) becomes

Γ*ξfDwfD(n*)dΓ=Γ*ξfDwfSn*dΓ,
(52)

while normal and shear force equilibrium (i.e., tfD*+tfS*=0) can be imposed at the third term on the left-hand side of Eq. (49), i.e.,

Γ*ηfSσfSn*dΓ=Γ*ηfS(pfDn*)dΓ+j=12Γ*ηfS×μfαSDk(wfSwfD)mj*mj*dΓ.
(53)

Finally, we apply Eqs. (24)–(27) in order to convert subdomain integrals (BD and BS) into an integral over the entire domain (B) and to also transform the interface equations [Eqs. (52) and (53)] into a set of immersed boundary conditions. As a result, we get the weak statements for a phase field Biot–Stokes model, which is to find {us,pfD,wfS,pfS,d} such that, for all {ηs,ξfD,ηfS,ξfS,ζ},

Gu=GDp=GSw=GSp=Gd=0,
(54)

where

Gu=Bηs:(σpfDI)(1d)dV+Bηs:σsddVBDtηst^DdΓ,
(55)
GDp=BξfD(u̇s)(1d)dVBξfDwfD(1d)dV+BξfD(wfSd)dVBDqξfDq^DdΓ,
(56)
GSw=BηfS:σfSddVBηfS(pfDd)dV+j=12BηfSμfαSDk(wfSwfD)mj*mj*ddVBStηfSt^SdΓ,
(57)
GSp=BξfS(u̇s)ddV+BξfS(wfS)ddV,
(58)
Gd=Bζg(d)H+Gcl*ddV+BζGcl*ddV.
(59)

Here, as pointed out by Stoter et al. (2017), the Γ-convergence ensures that the immersed boundary conditions imposed in Eqs. (56) and (57) are consistent with the interface conditions [Eqs. (44)–(46)] if l* → 0, which, in turn, confirms the mass conservation and force equilibrium for the entire system B.

This section highlights the capacities of the immersed phase field model to capture the hydromechanical interactions among the pore fluids in the cavities, the cracks, and the homogenized pore space and the host matrix in two numerical experiments. Our focus is on modeling the problems that involve the mechanically driven pore fluid migration due to deformation and crack growth inside the solid skeleton. The first example simulates the consolidation process of the porous material that contains a semi-circular cavity at the bottom that serves as a pore fluid outlet, while the second problem showcases the fracture-induced Stokes–Darcy flow in a vuggy porous medium.

In order to solve Eqs. (55)–(59) numerically, we adopt the standard finite element method where the solution procedure is based on operator-split (Miehe et al., 2010; Heister et al., 2015; and Suh et al., 2020) that successively updates the field variables. In other words, the phase field d is updated first by solving Gd = 0, while all other field variables are held fixed, and the solver then advances the remaining variables by solving {Gu,GDp,GSw,GSp}T = 0. The implementation of our proposed model including finite element discretization and the solution scheme relies on the finite element package FEniCS (Logg et al., 2012a; 2012b; and Alnæs et al., 2015). It is noted that there exist multiple different strategies to solve the same system of equations. Since the exploration of different solution schemes is out of the scope of this study, we omit the details of the implementation for brevity.

We first simulate a consolidation problem, which has always been one of the key problems in geotechnical engineering. While the classical consolidation problem considers time-dependent water expulsion from the homogeneous porous material, as illustrated in Fig. 3, this numerical example explores the case where the system includes a cavity at the bottom that serves as a pore fluid outlet. This specific setting is designed to simulate the mechanically driven Stokes–Darcy flow without significant changes in microstructural attributes.

FIG. 3.

Schematic of geometry and boundary conditions for the consolidation problem.

FIG. 3.

Schematic of geometry and boundary conditions for the consolidation problem.

Close modal

The problem domain is a water-saturated 1 × 2 m2 sized rectangular porous matrix (BD) that contains a semi-circular cavity (BS) whose diameter is 0.2 m. We prescribe a 1 kPa compressive mechanical traction at the top, while zero pressure boundary is imposed at the bottom of the cavity so that the time-dependent dissipation of pore pressure can be observed. The material parameters for this example are chosen as follows: intrinsic mass densities of the solid and fluid: ρs = 2700 kg/m3 and ρf = 1000 kg/m3; Young’s modulus and Poisson’s ratio of the solid skeleton: E = 100 MPa and ν = 0.25; initial permeability and initial porosity of the matrix: k0 = 1.0 × 10−8 m2 and ϕ0 = 0.4; dynamic viscosity of the fluid phase: μf = 1.0 × 10−3 Pa s; slippage coefficient αSD = 0; and regularization length of the interface l* = 0.002 m. Furthermore, we assume that the solid constituent remains intact in BD throughout the simulation, while the free fluid inside the cavity has zero particle concentration (i.e., c = 0).

Figure 4 shows the spatial distributions of the prime variables at t = 1.0 × 10−3 s. Here, we compute fluid pressure and relative fluid velocity of the entire system as pf=(1d)pfD+dpfS and wf=(1d)wfD+dwfS, respectively, since we have separate degrees of freedom for pore and free fluids residing in regions BD and BS. The results imply that the applied mechanical load t^D at t = 0 builds up the pore pressure, which, in turn, affects the pore fluid to migrate toward the cavity. Furthermore, the free fluid inside the cavity tends to exhibit higher velocity and lower pressure compared to those of the pore fluid because of different constitutive relations (i.e., Stokes equation and Darcy’s law) in each region. As illustrated in Fig. 5, we also investigate the time-dependent response of the system that clearly describes the consolidation process and, at the same time, highlights the continuous pressure and velocity fields along the y-axis (i.e., from the center point of the cavity to the top-central point of the external boundary). At t = 0, the entire load is taken by the incompressible pore water, which triggers the fluid flow inside the medium. This fluid flow is accompanied by a dissipation of pore pressure over time and an increase in the compression of the entire system, which is consistent with previous studies on homogeneous materials (Li et al., 2004; White and Borja, 2008; Kim et al., 2009; and Wang and Sun, 2016). In addition, the continuous pressure and velocity profiles imply that our model is capable of imposing mass continuity and force equilibrium at the interface as a set of immersed boundary conditions, which confirms the validity of the model.

FIG. 4.

Spatial distributions of the (a) phase field d, (b) solid displacement us (m), (c) fluid pressure pf=(1d)pfD+dpfS (Pa), and (d) relative fluid velocity wf=(1d)wfD+dwfS (m/s), at t = 1.0 × 10−3 s.

FIG. 4.

Spatial distributions of the (a) phase field d, (b) solid displacement us (m), (c) fluid pressure pf=(1d)pfD+dpfS (Pa), and (d) relative fluid velocity wf=(1d)wfD+dwfS (m/s), at t = 1.0 × 10−3 s.

Close modal
FIG. 5.

Response of the saturated Biot–Stokes system under 1 kPa consolidation pressure: (a) solid displacement; (b) fluid pressure; (c) fluid velocity.

FIG. 5.

Response of the saturated Biot–Stokes system under 1 kPa consolidation pressure: (a) solid displacement; (b) fluid pressure; (c) fluid velocity.

Close modal

In the second set of experiment, we conduct numerical simulations within two different types of domains that possess horizontal edge crack (Fig. 6): one explicitly captures the geometry of the large cavities in the porous media; another one captures the influence of the cavities by increasing the porosity of the homogenized effective medium. While the former approach adopts a more explicit representation of the pore geometry and, hence, may provide more detailed information on the interactions between the vugs and the propagating cracks, the latter approach could be numerically more efficient. Our objective is to demonstrate, quantitatively, the difference in the two approaches such that a fuller picture on the trade-off between computational efficiency, accuracy, and precision of the predictions can be established.

FIG. 6.

Schematic of geometry and boundary conditions for the fracture problem: (a) a domain with explicitly modeled cavities and (b) its homogenized counterpart.

FIG. 6.

Schematic of geometry and boundary conditions for the fracture problem: (a) a domain with explicitly modeled cavities and (b) its homogenized counterpart.

Close modal

1. Modeling vuggy porous media

As illustrated in Fig. 6(a), we first consider a domain that consists of the porous matrix with explicitly modeled cavities. Our first representation consists of totally nine cavities with different major and minor radii (Table I) that share the same aspect ratio of 2:1 and are tilted by 45°, such that the volume fraction of the cavities θcav is 0.056. Here, we assume that the solid skeleton inside the cavities is completely damaged (i.e., d = 1), while the porous matrix initially remains completely undamaged (i.e., d = 0). The material properties for this case are chosen as follows: ρs = 2700 kg/m3, ρf = 1000 kg/m3, E = 20 GPa, ν = 0.2, k0 = 1.0 × 10−12 m2, μf = 1.0 × 10−3 Pa s, αSD = 0.1, Gc=20 J/m2, and l* = 0.125 × 10−3 m. In addition, the initial particle concentration is chosen as c0 = 0.6 and its upper bound as cmax = 0.7, in order to mimic the mudflow inside the cracks or cavities (O’Brien and Julien, 1988; Iverson, 1997).

TABLE I.

Major and minor radii of the explicitly modeled elliptical cavities in Fig. 6(a).

Index123456789
Major radius ra (mm) 0.400 0.600 0.500 0.500 0.820 0.800 0.580 0.700 0.650 
Minor radius rb (mm) 0.200 0.300 0.250 0.250 0.410 0.400 0.290 0.350 0.325 
Index123456789
Major radius ra (mm) 0.400 0.600 0.500 0.500 0.820 0.800 0.580 0.700 0.650 
Minor radius rb (mm) 0.200 0.300 0.250 0.250 0.410 0.400 0.290 0.350 0.325 

In contrast, our second domain in Fig. 6(b) is a homogenized representation of Fig. 6(a), where all the cavities are considered a part of matrix pores. In this case, the porosity of the homogenized medium is determined as ϕhom = (1 − θcav)ϕ0 + θcav = 0.433. It is noted the correct homogenized effective properties often depend on the geometry of vugs or inclusions, which can be determined from computed tomographic images or directly obtained from the experiment (Sun et al.2011a; 2011b; Kim et al., 2016; and Lee et al., 2017). Since the micro-structural attributes are not always available, this study adopts an alternative approach where the effective material properties are determined by using the equivalent inclusion method (Hashin, 1960; Zimmerman, 1991; and Ramakrishnan and Arunachalam, 1993). Following Ramakrishnan and Arunachalam (1993) and by assuming that the matrix shares the same material properties as those chosen in Fig. 6(a), the effective bulk modulus (Khom) and shear modulus (μhom) for the homogenized representation [Fig. 6(b)] are determined as follows:

Khom=K(1θcav)21+1+ν2(12ν)θcav  ;  μhom=μ(1θcav)21+1119ν4(1+ν)θcav,
(60)

so that the effective Young’s modulus Ehom = 16.50 GPa and Poisson’s ratio νhom = 0.206. Assuming that the inclusion permeability is much higher than the matrix permeability, we approximate the effective permeability (khom) by following Markov et al. (2010), which is obtained based on Maxwell’s formula, i.e.,

khom=k0(1+2θcav)1θcav=1.18×1012[m2].
(61)

In addition, since all the cavities in Fig. 6(a) are completely isolated, we adopt the following effective critical energy Gc,hom proposed by Jelitto and Schneider (2018) for the homogenized representation, which depends on the volume fraction of the cavities:

Gc,hom=Gc(1θcav2/3)=17.07[J/m2].
(62)

2. Mechanically driven fracture-induced flow

As illustrated in Fig. 6, we conduct two different types of simulations within each domain: the tension tests with the prescribed vertical displacement rate of 0.01 × 10−3 m/s and the shear tests with the horizontal displacement rate of 0.01 × 10−3 m/s. In both tension and shear tests, the displacements are prescribed at the upper boundary, whereas the bottom part of the domain is held fixed. We also impose hydraulically insulated boundary conditions for the left and right boundaries, while we permit water intake from the upper and lower boundaries by imposing p^fD=0.

Figure 7 illustrates the evolution of the phase field for both tension and shear tests in a computational domain where the cavities are explicitly modeled, compared with the crack trajectories obtained from the homogenized domain. The domain without cavities exhibits the crack patterns that are similar to the results of previous studies on homogeneous solids (Miehe et al., 2010; Borden et al., 2012; Bryant and Sun, 2018; and Suh et al., 2020), while the domain with explicitly modeled cavities exhibits distinct crack patterns. More importantly, Figs. 8 and 9 reveal that neglecting the interaction between the cavity and the crack in the homogenized model may lead to over-simplified global responses that lack the distinctive characteristics of the cavity–crack coalescence.

FIG. 7.

Evolution of the phase field of the specimens subjected to the numerical tension and shear tests.

FIG. 7.

Evolution of the phase field of the specimens subjected to the numerical tension and shear tests.

Close modal
FIG. 8.

Force–displacement curves obtained from (a) tension and (b) shear tests.

FIG. 8.

Force–displacement curves obtained from (a) tension and (b) shear tests.

Close modal
FIG. 9.

Fluid influx at the top surface over time measured from (a) tension and (b) shear tests.

FIG. 9.

Fluid influx at the top surface over time measured from (a) tension and (b) shear tests.

Close modal

During the numerical experiments, the porous matrix initially undergoes linear elastic deformation until the crack nucleation takes place. At this point, since tensile loading directly influences the volume change of the material, both specimens under tensile load exhibit higher fluid influx at the top, compared to those measured from the shear tests. After the first peaks shown in Figs. 8(a) and 8(b), cracks start to initiate from the tips of the pre-existing flaw since they experience higher stress concentration compared to the matrix–cavity interface.

In both tensile and shear experiments performed on the vuggy specimen, the crack nucleation increases the surface influx rate at the permeable boundaries as the pore fluid starts to leak from the intact matrix to the damaged regions regardless of the spatial homogenization. The two numerical specimens, nevertheless, begin to behave differently when the cracks propagate toward the adjacent vugs and coalesce with each other in both tension and shear tests in the vuggy specimen (Fig. 7) (Qinami et al., 2019; Suh and Sun, 2020). These changes in surface influx cannot be replicated in the homogenized porous specimen as homogenization takes away the possibility of simulating the coalescence between the cavity and the crack.

After the coalescence of the cavity and the crack in the vuggy specimen, the reaction force in both cases increases again with lower influx rates until it reaches the second peak (i.e., where crack nucleation takes place at the matrix–cavity interface), and the crack eventually reaches both ends of the specimens. This result implies that the existence of vugs or cavities has a profound impact on the material behavior that cannot be easily replicated in the homogenized effective medium. Consequently, either a more effective macroscopic theory or a suitable multiscale technique is needed to incorporate the cavity–crack interaction into the predictions.

Figure 10 illustrates the pressure (pf) and x-directional velocity (wfS) fields from a domain with explicitly modeled cavities under tensile and shear loadings, at t = 0.476 s and t = 1.012 s, respectively, where cracks start to propagate from the cavities. Here, the superimposed arrows in Fig. 10 indicate the direction of the velocity vector wf=(1d)wfD+dwfS. In both cases, the leakage of pore fluid takes place toward the interconnected cracks and cavities at the middle, while the free fluid in BS tends to migrate toward the center, the region where a large crack opening displacement occurs. However, it is worthy to note that the fluid flow occurs from the region that has negative pressure to the damaged zone where pfS0. Unlike previous studies that use the cubic law to predict the hydraulic responses of the flow conduit (Mauthe and Miehe, 2017; Heider and Markert, 2017; Wang and Sun, 2017; and Chukwudozie et al., 2019), the pore pressure distribution inside the void space is governed by the Stokes equation directly. This set of numerical experiments again highlight that our proposed model is capable of simulating the fracture–cavity interaction with the evolving interface, which may not be easily captured either by using hydraulic phase field fracture models or by adopting the classical Biot–Stokes model with a sharp interface.

FIG. 10.

Snapshots of the pressure pf=(1d)pfD+dpfS (Pa) and velocity wfS (m/s) fields obtained from the tension (t = 0.476 s) and shear (t = 1.012 s) tests.

FIG. 10.

Snapshots of the pressure pf=(1d)pfD+dpfS (Pa) and velocity wfS (m/s) fields obtained from the tension (t = 0.476 s) and shear (t = 1.012 s) tests.

Close modal

To assess the computational efficiency of the proposed model, we record the CPU time for both simulations. A laptop with an Intel Core i9-9880H Processor CPU with 16 GB memory at 2667 MHz (DDR4) is used to run both simulations on a single core. Both simulations are solved by the same Scalable Nonlinear Equation Solver (SNES) available in FEniCS. In the case where vuggy pores are explicitly modeled, the time taken to assemble the system of equations is 1.13 s and the average time taken to advance one time step with (on average) five Newton–Raphson iterations is 35.69 s. Meanwhile, in the homogenized case, it takes 1.17 s to assemble the system of equations and 33.34 s to advance one time step also with (on average) five Newton–Raphson iterations. In this example, simulations with the explicitly captured vuggy pores require about 7% more CPU time to run the same simulation.

The future work may consider flow with a higher Reynolds number suitable for the Navier–Stokes equation in the fluid domain. Such an extension is, nevertheless, out of the scope of the current study.

This article presents a new immersed phase field model that captures the hydro-mechanical coupling mechanisms in vuggy porous media where brittle cracks filled with water may coalescence with pores that trigger both redistribution of flow and macroscopic softening that cannot be captured without the Stokes–Darcy flow. By generalizing the phase field as an indicator of defects, we introduce a simple and unified treatment to handle the evolving geometries due to crack growth and the resultant changes in constitutive responses without the need of re-meshing or introduction of enrichment functions. By directly simulating the flow inside the cracks, we bypass the need of introducing a phenomenological permeability enhancement model to replicate the flow conduit. This explicit approach can be advantageous over the embedded discontinuity approach when there are a substantial crack opening and a flow near the locations with void–crack interactions where a homogenized pore pressure jump would not be sufficient to capture the pattern of the pore pressure field in the defects. The future work may include the extension of the proposed model to three-dimensional cases, as well as extending the Stokes–Darcy flow model to the generalized Navier–Stokes–Darcy flow for injection and other problems with higher Reynolds numbers.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The suggestions of the two anonymous reviewers are gratefully acknowledged. The authors are thankful for the fruitful discussion with Professor Robert Zimmerman while the second author visited London in 2017. The first author was supported by the Earth Materials and Processes program from the US Army Research Office under Grant No. W911NF-18-2-0306. The second author was supported by the NSF CAREER grant from the Mechanics of Materials and Structures program at National Science Foundation under Grant No. CMMI-1846875 and the Dynamic Materials and Interactions program from the Air Force Office of Scientific Research under Grant No. FA9550-17-1-0169. These supports are gratefully acknowledged. The views and conclusions contained in this article are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the sponsors, including the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation herein.

1.
Alnæs
,
M.
,
Blechta
,
J.
,
Hake
,
J.
,
Johansson
,
A.
,
Kehlet
,
B.
,
Logg
,
A.
,
Richardson
,
C.
,
Ring
,
J.
,
Rognes
,
M. E.
, and
Wells
,
G. N.
, “
The FEniCS project version 1.5
,”
Arch. Numer. Software
3
(
100
),
9
(
2015
).
2.
Amor
,
H.
,
Marigo
,
J.-J.
, and
Maurini
,
C.
, “
Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments
,”
J. Mech. Phys. Solids
57
(
8
),
1209
1229
(
2009
).
3.
Arbogast
,
T.
and
Brunson
,
D. S.
, “
A computational method for approximating a Darcy–Stokes system governing a vuggy porous medium
,”
Comput. Geosci.
11
(
3
),
207
218
(
2007
).
4.
Arbogast
,
T.
and
Lehr
,
H. L.
, “
Homogenization of a Darcy–Stokes system modeling vuggy porous media
,”
Comput. Geosci.
10
(
3
),
291
302
(
2006
).
5.
Arson
,
C.
and
Pereira
,
J.-M.
, “
Influence of damage on pore size distribution and permeability of rocks
,”
Int. J. Numer. Anal. Methods Geomech.
37
(
8
),
810
831
(
2013
).
6.
Auriault
,
J.-L.
, “
About the beavers and Joseph boundary condition
,”
Transp. Porous Media
83
(
2
),
257
266
(
2010
).
7.
Badia
,
S.
,
Quaini
,
A.
, and
Quarteroni
,
A.
, “
Coupling Biot and Navier–Stokes equations for modelling fluid-poroelastic media interaction
,”
J. Comput. Phys.
228
(
21
),
7986
8014
(
2009
).
8.
Beavers
,
G. S.
and
Joseph
,
D. D.
, “
Boundary conditions at a naturally permeable wall
,”
J. Fluid Mech.
30
(
1
),
197
207
(
1967
).
9.
Bergkamp
,
E. A.
,
Verhoosel
,
C. V.
,
Remmers
,
J. J. C.
, and
Smeulders
,
D. M. J.
, “
A staggered finite element procedure for the coupled Stokes–Biot system with fluid entry resistance
,”
Comput. Geosci.
24
,
1497
(
2020
).
10.
Berryman
,
J. G.
and
Wang
,
H. F.
, “
The elastic coefficients of double-porosity models for fluid transport in jointed rock
,”
J. Geophys. Res.: Solid Earth
100
(
B12
),
24611
24627
, (
1995
).
11.
Bilgen
,
C.
and
Weinberg
,
K.
, “
On the crack-driving force of phase-field models in linearized and finite elasticity
,”
Comput. Methods Appl. Mech. Eng.
353
,
348
372
(
2019
).
12.
Blunt
,
M. J.
,
Jackson
,
M. D.
,
Piri
,
M.
, and
Valvatne
,
P. H.
, “
Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow
,”
Adv. Water Resour.
25
(
8-12
),
1069
1089
(
2002
).
13.
Borden
,
M. J.
,
Verhoosel
,
C. V.
,
Scott
,
M. A.
,
Hughes
,
T. J. R.
, and
Landis
,
C. M.
, “
A phase-field description of dynamic brittle fracture
,”
Comput. Methods Appl. Mech. Eng.
217-220
,
77
95
(
2012
).
14.
Borja
,
R. I.
, “
On the mechanical energy and effective stress in saturated and unsaturated porous continua
,”
Int. J. Solids Struct.
43
(
6
),
1764
1786
(
2006
).
15.
Borja
,
R. I.
and
Alarcón
,
E.
, “
A mathematical framework for finite strain elastoplastic consolidation part 1: Balance laws, variational formulation, and linearization
,”
Comput. Methods Appl. Mech. Eng.
122
(
1-2
),
145
171
(
1995
).
16.
Bourdin
,
B.
,
Francfort
,
G. A.
, and
Marigo
,
J.-J.
, “
The variational approach to fracture
,”
J. Elasticity
91
(
1-3
),
5
148
(
2008
).
17.
Bryant
,
E. C.
and
Sun
,
W.
, “
A mixed-mode phase field fracture model in anisotropic rocks with consistent kinematics
,”
Comput. Methods Appl. Mech. Eng.
342
,
561
584
(
2018
).
18.
Chapuis
,
R. P.
and
Aubertin
,
M.
, “
On the use of the Kozeny–Carman equation to predict the hydraulic conductivity of soils
,”
Can. Geotech. J.
40
(
3
),
616
628
(
2003
).
19.
Cheng
,
N.-S.
and
Law
,
A. W.-K.
, “
Exponential formula for computing effective viscosity
,”
Powder Technol.
129
(
1-3
),
156
160
(
2003
).
20.
Choo
,
J.
and
Borja
,
R. I.
, “
Stabilized mixed finite elements for deformable porous media with double porosity
,”
Comput. Methods Appl. Mech. Eng.
293
,
131
154
(
2015
).
21.
Choo
,
J.
and
Sun
,
W.
, “
Cracking and damage from crystallization in pores: Coupled chemo–hydro-mechanics and phase-field modeling
,”
Comput. Methods Appl. Mech. Eng.
335
,
347
379
(
2018
).
22.
Choo
,
J.
,
White
,
J. A.
, and
Borja
,
R. I.
, “
Hydromechanical modeling of unsaturated flow in double porosity media
,”
Int. J. Geomech.
16
(
6
),
D4016002
(
2016
).
23.
Chukwudozie
,
C.
,
Bourdin
,
B.
, and
Yoshioka
,
K.
, “
A variational phase-field model for hydraulic fracturing in porous media
,”
Comput. Methods Appl. Mech. Eng.
347
,
957
982
(
2019
).
24.
Class
,
H.
,
Helmig
,
R.
, and
Bastian
,
P.
, “
Numerical simulation of non-isothermal multiphase multicomponent processes in porous media
,”
Adv. Water Res.
25
(
5
),
533
550
(
2002
).
25.
Costa
,
A.
, “
Permeability-porosity relationship: A reexamination of the Kozeny–Carman equation based on a fractal pore-space geometry assumption
,”
Geophys. Res. Lett.
33
(
2
),
L02318
, (
2006
).
26.
Fu
,
P.
,
Johnson
,
S. M.
, and
Carrigan
,
C. R.
, “
An explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in arbitrary discrete fracture networks
,”
Int. J. Numer. Anal. Methods Geomech.
37
(
14
),
2278
2300
(
2013
).
27.
Ghaboussi
,
J.
and
Barbosa
,
R.
, “
Three-dimensional discrete element method for granular materials
,”
Int. J. Numer. Anal. Methods Geomech.
14
(
7
),
451
472
(
1990
).
28.
Grant
,
M.
,
Geothermal Reservoir Engineering
(
Elsevier
,
2013
).
29.
Guo
,
C.
,
Li
,
Y.
,
Nian
,
X.
,
Xu
,
M.
,
Liu
,
H.
, and
Wang
,
Y.
, “
Experimental study on the slip velocity of turbulent flow over and within porous media
,”
Phys. Fluids
32
(
1
),
015111
(
2020
).
30.
Hashin
,
Z.
, “
The elastic moduli of heterogeneous materials
,” Technical Report No. 9,
Harvard University
,
Cambridge, MA
,
1960
.
31.
Heider
,
Y.
and
Markert
,
B.
, “
A phase-field modeling approach of hydraulic fracture in saturated porous media
,”
Mech. Res. Commun.
80
,
38
46
(
2017
).
32.
Heider
,
Y.
and
Sun
,
W.
, “
A phase field framework for capillary-induced fracture in unsaturated porous media: Drying-induced vs. hydraulic cracking
,”
Comput. Methods Appl. Mech. Eng.
359
,
112647
(
2020
).
33.
Heister
,
T.
,
Wheeler
,
M. F.
, and
Wick
,
T.
, “
A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach
,”
Comput. Methods Appl. Mech. Eng.
290
,
466
495
(
2015
).
34.
Hyman
,
J. D.
,
Karra
,
S.
,
Makedonska
,
N.
,
Gable
,
C. W.
,
Painter
,
S. L.
, and
Viswanathan
,
H. S.
, “
dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport
,”
Comput. Geosci.
84
,
10
19
(
2015
).
35.
Iverson
,
R. M.
, “
The physics of debris flows
,”
Rev. Geophys.
35
(
3
),
245
296
(
1997
).
36.
Jaeger
,
J. C.
,
Cook
,
N. G. W.
, and
Zimmerman
,
R.
,
Fundamentals of Rock Mechanics
(
John Wiley & Sons
,
2009
).
37.
Jelitto
,
H.
and
Schneider
,
G. A.
, “
A geometric model for the fracture toughness of porous materials
,”
Acta Mater.
151
,
443
453
(
2018
).
38.
Jin
,
W.
and
Arson
,
C.
, “
Fluid-driven transition from damage to fracture in anisotropic porous media: A multi-scale XFEM approach
,”
Acta Geotech.
15
(
1
),
113
144
(
2020
).
39.
Juanes
,
R.
,
Spiteri
,
E. J.
,
Orr
, Jr.,
F. M.
, and
Blunt
,
M. J.
, “
Impact of relative permeability hysteresis on geological CO2 storage
,”
Water Resour. Res.
42
(
12
),
W12418
, (
2006
).
40.
Kang
,
D. H.
,
Yun
,
T. S.
,
Kim
,
K. Y.
, and
Jang
,
J.
, “
Effect of hydrate nucleation mechanisms and capillarity on permeability reduction in granular media
,”
Geophys. Res. Lett.
43
(
17
),
9018
9025
, (
2016
).
41et al..
Kim
,
J.
,
Tchelepi
,
H. A.
,
Juanes
,
R.
 et al, “
Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics
,” in
SPE Reservoir Simulation Symposium
(
Society of Petroleum Engineers
,
2009
).
42.
Kim
,
K. Y.
,
Suh
,
H. S.
,
Yun
,
T. S.
,
Moon
,
S.-W.
, and
Seo
,
Y.-S.
, “
Effect of particle shape on the shear strength of fault gouge
,”
Geosci. J.
20
(
3
),
351
359
(
2016
).
43.
Konzuk
,
J. S.
and
Kueper
,
B. H.
, “
Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture
,”
Water Resour. Res.
40
(
2
),
W02402
, (
2004
).
44.
Lade
,
P. V.
and
De Boer
,
R.
, “
The concept of effective stress for soil, concrete and rock
,”
Geotechnique
47
(
1
),
61
78
(
1997
).
45.
Layton
,
W. J.
,
Schieweck
,
F.
, and
Yotov
,
I.
, “
Coupling fluid flow with porous media flow
,”
SIAM J. Numer. Anal.
40
(
6
),
2195
2218
(
2002
).
46.
Lee
,
C.
,
Suh
,
H. S.
,
Yoon
,
B.
, and
Yun
,
T. S.
, “
Particle shape effect on thermal conductivity and shear wave velocity in sands
,”
Acta Geotech.
12
(
3
),
615
625
(
2017
).
47.
Leung
,
C. T. O.
and
Zimmerman
,
R. W.
, “
Estimating the hydraulic conductivity of two-dimensional fracture networks using network geometric properties
,”
Transp. Porous Media
93
(
3
),
777
797
(
2012
).
48.
Li
,
C.
,
Borja
,
R. I.
, and
Regueiro
,
R. A.
, “
Dynamics of porous media at finite strain
,”
Comput. Methods Appl. Mech. Eng.
193
(
36-38
),
3837
3870
(
2004
).
49.
Li
,
Z.
,
Lai
,
M.-C.
,
Peng
,
X.
, and
Zhang
,
Z.
, “
A least squares augmented immersed interface method for solving Navier–Stokes and Darcy coupling equations
,”
Comput. Fluids
167
,
384
399
(
2018
).
50.
Liu
,
C.
and
Abousleiman
,
Y. N.
, “
Shale dual-porosity dual-permeability poromechanical and chemical properties extracted from experimental pressure transmission tests
,”
J. Eng. Mech.
143
(
9
),
04017107
(
2017
).
51.
Logg
,
A.
,
Mardal
,
K.-A.
, and
Wells
,
G.
,
Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book
(
Springer Science & Business Media
,
2012a
), Vol. 84.
52.
Logg
,
A.
,
Wells
,
G. N.
, and
Hake
,
J.
, “
Dolfin: A C++/python finite element library
,” in
Automated Solution of Differential Equations by the Finite Element Method
(
Springer
,
2012b
), pp.
173
225
.
53.
Markov
,
M.
,
Kazatchenko
,
E.
,
Mousatov
,
A.
, and
Pervago
,
E.
, “
Permeability of the fluid-filled inclusions in porous media
,”
Transp. Porous Media
84
(
2
),
307
317
(
2010
).
54.
Mauran
,
S.
,
Rigaud
,
L.
, and
Coudevylle
,
O.
, “
Application of the Carman–Kozeny correlation to a high-porosity and anisotropic consolidated medium: The compressed expanded natural graphite
,”
Transp. Porous Media
43
(
2
),
355
376
(
2001
).
55.
Mauthe
,
S.
and
Miehe
,
C.
, “
Hydraulic fracture in poro–hydro-elastic media
,”
Mech. Res. Commun.
80
,
69
83
(
2017
).
56.
Mavko
,
G.
and
Nur
,
A.
, “
The effect of a percolation threshold in the Kozeny–Carman relation
,”
Geophysics
62
(
5
),
1480
1482
(
1997
).
57.
Miehe
,
C.
,
Hofacker
,
M.
, and
Welschinger
,
F.
, “
A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits
,”
Comput. Methods Appl. Mech. Eng.
199
(
45-48
),
2765
2778
(
2010
).
58.
Miehe
,
C.
and
Mauthe
,
S.
, “
Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in hydro–poro-elasticity and hydraulic fracturing of fluid-saturated porous media
,”
Comput. Methods Appl. Mech. Eng.
304
,
619
655
(
2016
).
59.
Mikelic
,
A.
and
Jäger
,
W.
, “
On the interface boundary condition of Beavers, Joseph, and Saffman
,”
SIAM J. Appl. Math.
60
(
4
),
1111
1127
(
2000
).
60.
Monchiet
,
V.
,
Ly
,
H.-B.
, and
Grande
,
D.
, “
Macroscopic permeability of doubly porous materials with cylindrical and spherical macropores
,”
Meccanica
54
(
10
),
1583
1596
(
2019
).
61.
Mooney
,
M.
, “
The viscosity of a concentrated suspension of spherical particles
,”
J. Colloid Sci.
6
(
2
),
162
170
(
1951
).
62.
Mumford
,
D. B.
and
Shah
,
J.
, “
Optimal approximations by piecewise smooth functions and associated variational problems
,”
Commun. Pure Appl. Math.
39
,
2673
(
1989
).
63.
Na
,
S.
and
Sun
,
W.
, “
Computational thermomechanics of crystalline rock, Part I: A combined multi-phase-field/crystal plasticity approach for single crystal simulations
,”
Comput. Methods Appl. Mech. Eng.
338
,
657
691
(
2018
).
64.
Nguyen
,
T. T.
,
Yvonnet
,
J.
,
Bornert
,
M.
,
Chateau
,
C.
,
Sab
,
K.
,
Romani
,
R.
, and
Le Roy
,
R.
, “
On the choice of parameters in the phase field method for simulating crack initiation with experimental validation
,”
Int. J. Fract.
197
(
2
),
213
226
(
2016
).
65.
O’Brien
,
J. S.
and
Julien
,
P. Y.
, “
Laboratory analysis of mudflow properties
,”
J. Hydraul. Eng.
114
(
8
),
877
887
(
1988
).
66et al..
Ozkan
,
E.
,
Raghavan
,
R. S.
,
Apaydin
,
O. G.
 et al, “
Modeling of fluid transfer from shale matrix to fracture network
,” in
SPE Annual Technical Conference and Exhibition
(
Society of Petroleum Engineers
,
2010
).
67.
Paterson
,
A.
and
Fermigier
,
M.
, “
Wetting of heterogeneous surfaces: Influence of defect interactions
,”
Phys. Fluids
9
(
8
),
2210
2216
(
1997
).
68.
Pereira
,
J.-M.
and
Arson
,
C.
, “
Retention and permeability properties of damaged porous rocks
,”
Comput. Geotech.
48
,
272
282
(
2013
).
69.
Pham
,
K. H.
,
Ravi-Chandar
,
K.
, and
Landis
,
C. M.
, “
Experimental validation of a phase-field model for fracture
,”
Int. J. Fract.
205
(
1
),
83
101
(
2017
).
70.
Qinami
,
A.
,
Bryant
,
E. C.
,
Sun
,
W.C.
, and
Kaliske
,
M.
, “
Circumventing mesh bias by r-and h-adaptive techniques for variational eigenfracture
,”
Int. J. Fract.
220
(
2
),
129
142
(
2019
).
71.
Ramakrishnan
,
N.
and
Arunachalam
,
V. S.
, “
Effective elastic moduli of porous ceramic materials
,”
J. Am. Ceram. Soc.
76
(
11
),
2745
2752
(
1993
).
72.
Rutqvist
,
J.
,
Birkholzer
,
J.
,
Cappa
,
F.
, and
Tsang
,
C.-F.
, “
Estimating maximum sustainable injection pressure during geological sequestration of Co2 using coupled fluid flow and geomechanical fault-slip analysis
,”
Energy Convers. Manage.
48
(
6
),
1798
1807
(
2007
).
73.
Saffman
,
P. G.
, “
On the boundary condition at the surface of a porous medium
,”
Stud. Appl. Math.
50
(
2
),
93
101
(
1971
).
74.
Schutjens
,
P. M. T. M.
,
Hanssen
,
T. H.
,
Hettema
,
M. H. H.
,
Merour
,
J.
,
De Bree
,
P.
,
Coremans
,
J. W. A.
,
Helliesen
,
G.
 et al, “
Compaction-induced porosity/permeability reduction in sandstone reservoirs: Data and model for elasticity-dominated deformation
,”
SPE Reservoir Eval. Eng.
7
(
03
),
202
216
(
2004
).
75.
Selvadurai
,
A. P. S.
,
Couture
,
C.-B.
, and
Rezaei Niya
,
S. M.
, “
Permeability of wormholes created by CO2-acidized water flow through stressed carbonate rocks
,”
Phys. Fluids
29
(
9
),
096604
(
2017
).
76.
Spaid
,
M. A. A.
and
Phelan
,
F. R.
, Jr.
, “
Lattice Boltzmann methods for modeling microscale flow in fibrous porous media
,”
Phys. Fluids
9
(
9
),
2468
2474
(
1997
).
77.
Stoter
,
S. K. F.
,
Müller
,
P.
,
Cicalese
,
L.
,
Tuveri
,
M.
,
Schillinger
,
D.
, and
Hughes
,
T. J. R.
, “
A diffuse interface method for the Navier–Stokes/Darcy equations: Perfusion profile for a patient-specific human liver based on MRI scans
,”
Comput. Methods Appl. Mech. Eng.
321
,
70
102
(
2017
).
78.
Suh
,
H. S.
,
Kang
,
D. H.
,
Jang
,
J.
,
Kim
,
K. Y.
, and
Yun
,
T. S.
, “
Capillary pressure at irregularly shaped pore throats: Implications for water retention characteristics
,”
Adv. Water Resour.
110
,
51
58
(
2017
).
79.
Suh
,
H. S.
,
Sun
,
W.
, and
O’Connor
,
D. T.
, “
A phase field model for cohesive fracture in micropolar continua
,”
Comput. Methods Appl. Mech. Eng.
369
,
113181
(
2020
), http://www.sciencedirect.com/science/article/pii/S0045782520303662.
80.
Suh
,
H. S.
and
Sun
,
W. C.
, “
An open-source fenics implementation of a phase field fracture model for micropolar continua
,”
Int. J. Multiscale Comput. Eng.
17
(
6
),
639
(
2019
).
81.
Suh
,
H. S.
and
Sun
,
W. C.
, “
An immersed phase field fracture model for fluid-infiltrating porous media with evolving Beavers–Joseph–Saffman condition
,”
E3S Web Conf.
205
,
03009
(
2020
).
82.
Suh
,
H. S.
and
Yun
,
T. S.
, “
Modification of capillary pressure by considering pore throat geometry with the effects of particle shape and packing features on water retention curves for uniformly graded sands
,”
Comput. Geotech.
95
,
129
136
(
2018
).
83.
Sun
,
W.
,
Cai
,
Z.
, and
Choo
,
J.
, “
Mixed Arlequin method for multiscale poromechanics problems
,”
Int. J. Numer. Methods Eng.
111
(
7
),
624
659
(
2017
).
84.
Sun
,
W.
,
Ostien
,
J. T.
, and
Salinger
,
A. G.
, “
A stabilized assumed deformation gradient finite element formulation for strongly coupled poromechanical simulations at finite strain
,”
Int. J. Numer. Anal. Methods Geomech.
37
(
16
),
2755
2788
(
2013
).
85.
Sun
,
W.
and
Wong
,
T.-f.
, “
Prediction of permeability and formation factor of sandstone with hybrid lattice Boltzmann/finite element simulation on microtomographic images
,”
Int. J. Rock Mech.Min. Sci.
106
,
269
277
(
2018
).
86.
Sun
,
W. C.
,
Andrade
,
J. E.
, and
Rudnicki
,
J. W.
, “
Multiscale method for characterization of porous microstructures and their impact on macroscopic effective permeability
,”
Int. J. Numer. Methods Eng.
88
(
12
),
1260
1279
(
2011a
).
87.
Sun
,
W. C.
,
Andrade
,
J. E.
,
Rudnicki
,
J. W.
, and
Eichhubl
,
P.
, “
Connecting microstructural attributes and permeability from 3D tomographic images of in situ shear-enhanced compaction bands using multiscale computations
,”
Geophys. Res. Lett.
38
(
10
), (
2011b
).
88.
Tang
,
G. H.
,
Tao
,
W. Q.
, and
He
,
Y. L.
, “
Lattice Boltzmann method for gaseous microflows using kinetic theory boundary conditions
,”
Phys. Fluids
17
(
5
),
058101
(
2005
).
89.
Terzis
,
A.
,
Zarikos
,
I.
,
Weishaupt
,
K.
,
Yang
,
G.
,
Chu
,
X.
,
Helmig
,
R.
, and
Weigand
,
B.
, “
Microscopic velocity field measurements inside a regular porous medium adjacent to a low Reynolds number channel flow
,”
Phys. Fluids
31
(
4
),
042001
(
2019
).
90.
Wagner
,
J. L.
,
Casper
,
K. M.
,
Beresh
,
S. J.
,
Hunter
,
P. S.
,
Spillers
,
R. W.
,
Henfling
,
J. F.
, and
Mayes
,
R. L.
, “
Fluid-structure interactions in compressible cavity flows
,”
Phys. Fluids
27
(
6
),
066102
(
2015
).
91.
Wang
,
K.
and
Sun
,
W.
, “
A semi-implicit discrete-continuum coupling method for porous media based on the effective stress principle at finite strain
,”
Comput. Methods Appl. Mech. Eng.
304
,
546
583
(
2016
).
92.
Wang
,
K.
and
Sun
,
W.
, “
A unified variational Eigen–Erosion framework for interacting brittle fractures and compaction bands in fluid-infiltrating porous media
,”
Comput. Methods Appl. Mech. Eng.
318
,
1
32
(
2017
).
93.
Wang
,
K.
and
Sun
,
W.
, “
A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deep learning
,”
Comput. Methods Appl. Mech. Eng.
334
,
337
380
(
2018
).
94.
Wang
,
K.
and
Sun
,
W.
, “
An updated Lagrangian LBM-DEM-FEM coupling model for dual-permeability fissured porous media with embedded discontinuities
,”
Comput. Methods Appl. Mech. Eng.
344
,
276
305
(
2019
).
95.
White
,
J. A.
and
Borja
,
R. I.
, “
Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients
,”
Comput. Methods Appl. Mech. Eng.
197
(
49-50
),
4353
4366
(
2008
).
96.
Witherspoon
,
P. A.
,
Wang
,
J. S. Y.
,
Iwai
,
K.
, and
Gale
,
J. E.
, “
Validity of cubic law for fluid flow in a deformable rock fracture
,”
Water Resour. Res.
16
(
6
),
1016
1024
, (
1980
).
97.
Worthington
,
P. F.
, “
The uses and abuses of the Archie equations, 1: The formation factor-porosity relationship
,”
J. Appl. Geophys.
30
(
3
),
215
228
(
1993
).
98.
Wu
,
Z.
and
Mirbod
,
P.
, “
Experimental analysis of the flow near the boundary of random porous media
,”
Phys. Fluids
30
(
4
),
047103
(
2018
).
99.
Zimmerman
,
R. W.
, “
Elastic moduli of a solid containing spherical inclusions
,”
Mech. Mater.
12
(
1
),
17
24
(
1991
).
100.
Zimmerman
,
R. W.
, “
Coupling in poroelasticity and thermoelasticity
,”
Int. J. Rock Mech. Min. Sci.
37
(
1-2
),
79
87
(
2000
).