Hydromechanical model for hydraulic fractures using XFEM

Abstract

In this study, a hydromechanical model for fluid flow in fractured porous media is presented. We assume viscous fluids and the coupling equations are derived from the mass and momentum balance equations for saturated porous media. The fluid flow through discrete cracks will be modelled by the extended finite element method and an implicit time integration scheme. We also present a consistent linearization of the underlying non-linear discrete equations. They are solved by the Newton-Raphson iteration procedure in combination with a line search. Furthermore, the model is extended to includes crack propagation. Finally, examples are presented to demonstrate the versatility and efficiency of this two-scale hydromechanical model. The results suggest that the presence of the fracture in a deforming, porous media has great impact on the fluid flow and deformation patterns.

Keywords

multi-phase medium / porous / fracture / multi-scale method

Introduction

In order to gain a better understanding of hydraulic fracturing, it is of utmost importance to consider the fluid flow through the propagating fracture network [¹]. In particular in the past two decades, there has been a tremendous improvement of computational methods for fracture. One of the most popular method includes the extended finite element method [²,³] which has been developed in 1999 in order to model crack growth without remeshing. It has been also applied to fluid-structure interaction [⁴,⁵] and also to fluid flow through cracks [⁶,⁷]. An interesting alternative to extended finite element method (XFEM) is the so-called phantom node method [⁸]. The key idea has been proposed by Hansbo and Hansbo [⁹]. Instead of enrichment functions, the phantom node method exploits overlapping elements which significantly facilitates its implementation [¹⁰,¹¹]. While the original phantom node method requires the crack to cross the entire element, special crack tip elements have been proposed improving computational efficiency and allowing for more complex crack patterns [¹²,¹³]. Instead of a nodal enrichment the embedded finite element method [¹⁴] enriches the element and exploits static condensation removing the additional degrees of freedom at element level. We also would like to mention some recent very efficient remeshing techniques [¹⁵–¹⁹] which are strong competitors to enriched method; these methods do not require modification of very efficient elements and exploit the advantages of vast improvements in remeshing done the past years. Remeshing has also been used in cohesive element approaches [²⁰,²¹] which are popular methods for dynamic fracture. An alternative to finite element methods are meshfree methods [²²–²⁵]. While early approaches capture discrete cracks with techniques such as the visibility criterion [²⁶], diffraction criterion [²⁷] and the transparency criterion [²⁷], newer approaches exploit also the concept of enrichment [²⁸–³²]. Those methods have also been extended to very complex three-dimensional fracture [³³–³⁵] including to fracture due to fluid-structure interaction [⁶].

There are other efficient methods for fracture such as phase field approaches [³⁶,³⁷] or the cracking particles method [³⁸–⁴⁰]. However, they smear the crack over a certain width and cannot model fluid flow through discrete cracks. The cracking particles method is a discrete crack approach and does not smear the crack over a certain region. However, the crack is modelled as a set of cracked nodes and therefore it cannot easily model fluid flow through these crack segments neither. A very interesting approach to fracture is peridynamics (PD) [⁴¹]. PD is a non-local theory which elegantly unifies continuous and discontinuous deformations and hence fracture is part of the solution and not part of the problem. However, to model fluid flow through discrete cracks is quite challenging in PD. Other approaches to fracture include the discrete element method (DEM) [⁴²–⁴⁴], the particle flow code (PFC) [⁴⁵–⁴⁷] or the discontinuous deformation analysis (DDA) [⁴⁸–⁵¹]. However, calibration of the ’microscopic’ material parameters for those approaches in order to predict “macroscopic” properties for complex material behavior is complicated.

In this study, we present an extended finite element formulation for fluid flow in porous media (Fig. 1). Our approach is based on the method proposed in Ref. [⁷].

In the next section, we present the governing equations in strong and weak form. Sections 3 and 4 describe the XFEM formulation and discrete system of equations. Subsequently, we present computational results before we conclude the manuscript in Sections 5 and 6.

Governing equations for the bulk

Strong form

Balance of momentum

We assume small displacement theory, no mass transfer between the constituents and isothermal conditions [⁵²]. Hence, the balance equation of momentum for a two phase material reads

(1)\nabla ·{{\displaystyle \mathit{\sigma }}}_{\pi }+{{\displaystyle \widehat{\mathit{p}}}}_{\pi }+{{\displaystyle \rho }}_{\pi }\mathit{g}={\displaystyle \frac{\partial \left({{\displaystyle \rho }}_{\pi }{{\displaystyle \mathit{v}}}_{\pi }\right)}{\partial t}}+\nabla ·\left({{\displaystyle \rho }}_{\pi }{{\displaystyle \mathit{v}}}_{\pi }\otimes {{\displaystyle \mathit{v}}}_{\pi }\right),

where the subscripts \pi =s,f, denote the solid and fluid phase, respectively, \rho represents the densities, v the absolute velocities and \mathit{\sigma } the stress of the constituents. Furthermore, {{\displaystyle \widehat{\mathit{p}}}}_{\pi } represents the source of the momentum from the other constituent, which accounts for the the load drag between the solid and fluid phase. For a closed system, the source of the momentum must fulfill:

(2){{\displaystyle \widehat{\mathit{p}}}}_s+{{\displaystyle \widehat{\mathit{p}}}}_f=\mathrm{0.}

Neglecting body forces and the inertia terms, Eq. (1) is reduced to

(3)\nabla ·{{\displaystyle \mathit{\sigma }}}_{\pi }+{{\displaystyle \widehat{\mathit{p}}}}_{\pi }=\mathrm{0.}

Balance of mass

For a two phase medium, the continuity equation is given by Ref. [⁵³]

(4){\displaystyle \frac{\partial {{\displaystyle \rho }}_{\pi }}{\partial t}}+\nabla ·\left({{\displaystyle \rho }}_{\pi }{{\displaystyle \mathit{v}}}_{\pi }\right)-\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{{\displaystyle \rho }}_{\pi }·{{\displaystyle \mathit{v}}}_{\pi }=\mathrm{0},

where \rho is the mass density and v the absolute velocity. The volume ratio of the solid and fluid phase n_s and n_f, respectively, have to fulfill the following equation:

(5){{\displaystyle n}}_s+{{\displaystyle n}}_f=\mathrm{1.}

The apparent mass density for each constituent is obtained by

(6){{\displaystyle \rho }}_{\pi }={{\displaystyle n}}_{\pi }\times {{\displaystyle \overline{\rho }}}_{\pi },

where \overline{\rho } is the absolute mass density. Substituting Eq. (6) into Eq. (4) gives:

(7){{\displaystyle n}}_s\nabla ·{{\displaystyle \mathit{v}}}_s+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \overline{\rho }}}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}=\mathrm{0},

(8){{\displaystyle n}}_f\nabla ·{{\displaystyle \mathit{v}}}_f+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \overline{\rho }}}_f}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_f}{\partial t}}=\mathrm{0.}

Combining Eq. (5) with Eq. (7) and Eq. (8) gives:

(9)\nabla ·{{\displaystyle \mathit{v}}}_s+{{\displaystyle n}}_f\nabla ·\left({{\displaystyle \mathit{v}}}_f-{{\displaystyle \mathit{v}}}_s\right)+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho }}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho }}_f}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_f}{\partial t}}=\mathrm{0.}

For a compressible solid, the time derivative of the density of the solid phase is obtained from the mass conservative equation [⁵³] by

(10){\displaystyle \frac{\partial \left({{\displaystyle \rho }}_s{{\displaystyle \mathit{v}}}_s\right)}{\partial t}}=\mathrm{0.}

From the entropy inequality [⁵⁴] for unsaturated flow accounting for interfaces, the pressure in the solid phase is:

(11){{\displaystyle p}}_s={{\displaystyle p}}_f\times {{\displaystyle n}}_f.

Assuming the solid density is a function of pressure p_s and temperature, yields

(12){\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho }}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \mathrm{v}}}_s}}{\displaystyle \frac{\partial {{\displaystyle \mathrm{v}}}_s}{\partial t}}={\displaystyle \frac{\mathrm{1}}{{{\displaystyle K}}_s}}{\displaystyle \frac{\partial {{\displaystyle p}}_s}{\partial t}}+{{\displaystyle \beta }}_s{\displaystyle \frac{\partial {{\displaystyle p}}_s}{\partial t}}.

where K_s designates the bulk modulus of solid phase,

{{\displaystyle \beta }}_s

the thermal expansion coefficient and T_s the temperature. Since the whole process is under isothermal, the last item can be omitted:

(13){\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho }}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \mathrm{v}}}_s}}{\displaystyle \frac{\partial {{\displaystyle \mathrm{v}}}_s}{\partial t}}={\displaystyle \frac{\mathrm{1}}{{{\displaystyle K}}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}.

Let us define the Biot’ constant [⁵³] as

(14)\mathrm{1}-\alpha ={\displaystyle \frac{{{\displaystyle K}}_T}{{{\displaystyle K}}_s}},

with K_T being the overall bulk modulus of the two phase medium and the Biot coefficient. The change of the solid mass density is related to its volume change by

(15)\nabla ·{{\displaystyle \mathit{v}}}_s=-{\displaystyle \frac{{{\displaystyle K}}_s}{{{\displaystyle K}}_T}}{\displaystyle \frac{{{\displaystyle n}}_s}{{{\displaystyle \rho }}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}.

Substituting now Eq. (14) into Eq. (15) leads to

(16)\left(\mathrm{1}-\alpha \right)\nabla ·{{\displaystyle \mathit{v}}}_s={\displaystyle \frac{{{\displaystyle n}}_s}{{{\displaystyle \rho }}_s}}{\displaystyle \frac{\partial {{\displaystyle \rho }}_s}{\partial t}}.

For the fluid phase, the relationship between the incremental change of fluid density and of the fluid pressure reads:

(17){\displaystyle \frac{\mathrm{1}}{Q}}\mathrm{d}p={\displaystyle \frac{{{\displaystyle n}}_f}{{{\displaystyle \rho }}_f}}\mathrm{d}{{\displaystyle p}}_f,

with the compressibility modulus:

(18){\displaystyle \frac{\mathrm{1}}{Q}}={\displaystyle \frac{\alpha -{{\displaystyle n}}_f}{{{\displaystyle K}}_s}}+{\displaystyle \frac{{{\displaystyle n}}_f}{{{\displaystyle K}}_f}},

where K_f is the bulk modulus of the fluid phase. Inserting Eq. (17) and Eq. (18) into Eq. (9) the balance equation of mass:

(19)\alpha \nabla ·{{\displaystyle \mathit{v}}}_s+{{\displaystyle n}}_f\nabla ·\left({{\displaystyle \mathit{v}}}_f-{{\displaystyle \mathit{v}}}_s\right)+{\displaystyle \frac{\mathrm{1}}{Q}}{\displaystyle \frac{\partial p}{\partial t}}=\mathrm{0.}

Kinematic relation

Assuming a linear elastic solid, the kinematic relation for small strain theory reads:

(20){{\displaystyle \mathit{\epsilon }}}_s={\nabla }^s{{\displaystyle \mathit{u}}}_s,

with {{\displaystyle \mathit{\epsilon }}}_s and u_s being the linear strain tensor and displacement field of the solid phase.

Constitutive relation

The effective stress increment d{{\displaystyle {\mathit{\sigma }}^{\prime }}}_s in the solid reads:

(21)\mathrm{d}{{\displaystyle {\mathit{\sigma }}^{\prime }}}_s={\displaystyle \frac{\mathrm{d}{\mathit{\sigma }}_s}{{{\displaystyle n}}_s}},

the incremental stress-strain relationship for the solid media:

(22)\mathrm{d}{{\displaystyle \mathit{\sigma }}}_s=\mathit{D}:\mathrm{d}{{\displaystyle \mathit{\epsilon }}}_s,

with D being the fourth-order elasticity tensor.

Boundary conditions

The boundary conditions for the two phase media (see Fig. 2) are given by

(23){{\displaystyle \mathit{n}}}_{\Gamma }·\mathit{\sigma }={{\displaystyle \mathit{t}}}_p\mathrm{on}\partial {{\displaystyle \mathrm{\Omega }}}_{\mathrm{t}}\mathrm{;}\mathit{u}={{\displaystyle \mathit{u}}}_p\mathrm{on}\partial {{\displaystyle \mathrm{\Omega }}}_u,

with the von Neumann and Dirichlet boundaries \partial {{\displaystyle \mathrm{\Omega }}}_{\mathit{t}} and \partial {{\displaystyle \mathrm{\Omega }}}_{\mathit{u}} and \partial \mathrm{\Omega }=\partial {{\displaystyle \mathrm{\Omega }}}_{\mathit{t}}\cup \partial {{\displaystyle \mathrm{\Omega }}}_{\mathit{u}}, t_p and u_p describe the traction and displacement, respectively. According to Darcy’s law [⁵⁵] for an isotropic media, the fluid velocity (v) is related to the Darcy flux (q) by the porosity (\phi).

(24)\phi \mathit{v}=\mathit{q}.

On complementary part of the boundary \partial {{\displaystyle \mathrm{\Omega }}}_q and \partial {{\displaystyle \mathrm{\Omega }}}_p, with \partial \Omega =\partial {{\displaystyle \mathrm{\Omega }}}_q\cup \partial {{\displaystyle \mathrm{\Omega }}}_p hold:

(25){{\displaystyle n}}_f\left({{\displaystyle \mathit{v}}}_f-{{\displaystyle \mathit{v}}}_s\right)·{{\displaystyle n}}_{\Gamma }={{\displaystyle \mathit{q}}}_p·{{\displaystyle n}}_{\Gamma }\mathrm{on}\partial {{\displaystyle \mathrm{\Omega }}}_q\mathrm{;}p={{\displaystyle p}}_p\mathrm{on}\partial {{\displaystyle \mathrm{\Omega }}}_p,

with {{\displaystyle \mathit{q}}}_p·{{\displaystyle n}}_{\Gamma } and p are describing the outflow of pore fluid and pressure, respectively.

Weak form and coupling

To derive the weak form from the strong form, we multiply the momentum balance Eq. (3) and the mass balance Eq. (19) with kinematically admissible test functions for the displacements and pressure. Integrating by parts, Gauss divergence theorem and applying Darcy’s law finally leads to the well-known weak form:

(26){\displaystyle {\int }_{\Omega }\left(\nabla ·\mathit{\eta }\right)}·\mathit{\sigma }\mathrm{d}\mathrm{\Omega }+\begin{array}{|l|}\hline {\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}〚\mathit{\eta }〛·\mathit{\sigma }·{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}}\mathrm{d}\mathrm{\Omega }\hfill \\ \hline\end{array}={\int }_{{\displaystyle \Gamma }}\mathit{\eta }·{\mathit{t}}_{\mathit{p}}\mathrm{d}\mathrm{\Omega },

and

(27)-{\displaystyle {\int }_{\mathrm{\Omega }}\alpha \mathit{\xi }\nabla }·{{\displaystyle \mathit{v}}}_s\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}{{\displaystyle k}}_f}\nabla \mathit{\xi }·\nabla p\mathrm{d}\mathrm{\Omega }-{\displaystyle {\int }_{\mathrm{\Omega }}\mathit{\xi }{Q}^{-\mathrm{1}}\dot{p}\mathrm{d}\mathrm{\Omega }}

where \dot{p}={\displaystyle \frac{\partial p}{\partial t}} denotes the time derivatives. The traction force on {{\displaystyle \Gamma }}_d is induced by the flow pressure inside cavity. Due to the presence of the discontinuity inside the domain, the traction force \mathit{\sigma }·{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d} on {{\displaystyle \Gamma }}_d and fluid flux 〚\mathit{\xi }{{\displaystyle n}}_f\left({{\displaystyle \mathit{v}}}_f-{{\displaystyle v}}_s\right)〛 through the discontinuity face are essential parts of the weak form. Since the length to width ratio of the cavity is considerable large, one could assume that the traction force on each side of the cavity is equal. Because of the continuity from the bulk to cavity, this gives:

(28)\mathit{\sigma }·{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}=-p{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}\mathrm{on}{{\displaystyle \Gamma }}_d.

Substituting Eq. (28) into Eq. (26), the final weak form of the balance equation of momentum gives

(29){\displaystyle {\int }_{\mathrm{\Omega }}\left(\nabla ·\mathit{\eta }\right)}·\mathit{\sigma }\mathrm{d}\mathrm{\Omega }-{\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}〚\eta 〛}·p{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}\mathrm{d}\mathrm{\Omega }={\displaystyle {\int }_{\Gamma }\mathit{\eta }·{{\displaystyle \mathit{t}}}_p}\mathrm{d}\mathrm{\Omega }.

Darcy’s law is expressed by

(30)\mathit{v}={\displaystyle \frac{\mathit{q}}{\phi }},

where \phi is the porosity of the bulk and v=v_f−v_s the velocity. The pressure values for both faces of the cavity are identical leads to the following coupling term of the weak form for the mass balance equation:

(31){\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}·〚\xi {{\displaystyle n}}_f\left({{\displaystyle \mathit{v}}}_f-{{\displaystyle \mathit{v}}}_s\right)〛}\mathrm{d}\Gamma ={\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}\xi }{{\displaystyle n}}_f{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}·〚{{\displaystyle \mathit{v}}}_f-{{\displaystyle \mathit{v}}}_s〛\mathrm{d}\Gamma ={\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}\xi }{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}·{{\displaystyle \mathit{q}}}_d\mathrm{d}\Gamma ,

with q_d being the flow flux through the discontinuity.

Fluid flow inside the cavity

Assuming the fluid flow inside the cavity is a Newtonian fluid, the general balance equation of momentum reads [⁵⁶]

(32)\underset{\mathrm{V}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{\frac{\partial \mathbf{u}}{\partial t}}}{\displaystyle }+\underset{\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{\left(\mathit{u}·\nabla \right)\mathit{u}}}-\underset{\mathrm{D}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{\mu {\nabla }^{\mathrm{2}}\mathit{u}}}=\underset{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}{\underbrace{-\nabla w}}+\underset{\mathrm{E}\mathrm{x}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}{\underbrace{g}},

where \mathit{\mu } being the viscosity of the fluid, u being the velocity of the fluid.

Assuming small displacements, no mass exchange and neglecting body forces and inertia forces, the balance equation of momentum is simplified to:

(33)\underset{\mathrm{D}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{u}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{-\mu {\nabla }^{\mathrm{2}}\mathit{u}}}=\underset{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}{\underbrace{-\nabla w}}.

In the two phase medium, the interface pressure on the two sides of the cavity serves as internal source of momentum, which gives the final balance equation of momentum:

(34)\mu {\nabla }^{\mathrm{2}}{{\displaystyle \mathit{v}}}_{\overline{f}}=\nabla p,

with {{\displaystyle \mathit{v}}}_{\overline{f}} denotes the fluid flow velocity, the subscript \overline{f} denotes fluid inside cavity.

Due to the high length to cross-section ratio (Fig. 3), the fluid flow inside cavity can be considered as quasi one dimensional flow given by

(35){\displaystyle \frac{\partial {{\displaystyle p}}_{\overline{f}}}{\partial \overline{y}}}=\mathrm{0},

in the normal direction {{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d} and in the tangential direction of {{\displaystyle \mathit{t}}}_{{{\displaystyle \Gamma }}_d}:

(36){\displaystyle \frac{\partial {{\displaystyle p}}_{\overline{f}}}{\partial \overline{x}}}=\mu {\displaystyle \frac{{\partial }^{\mathrm{2}}{{\displaystyle v}}_{\overline{f}}}{\partial {{\displaystyle \overline{y}}}^{\mathrm{2}}}},

with \overline{x} and \overline{y} coordinations respect to the normal and tangent directions of the cavity, respectively. In the normal direction {{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d} of the cavity, the pressure {{\displaystyle p}}_{\overline{f}} is constant through the cross section of the cavity, for continuity restriction, {{\displaystyle p}}_{\overline{f}} must be equal to p.

(37){{\displaystyle p}}_{\overline{f}}=p.

In order to derive the axial velocity, Eq. (36) is integrate twice respecting the y coordinate:

(38){\displaystyle {\int }_{-h}^h{\displaystyle {\int }_{-h}^h{\displaystyle \frac{\partial p}{\partial \overline{x}}}}}={\displaystyle {\int }_{-h}^h{\displaystyle {\int }_{-h}^h\mu {\displaystyle \frac{{\partial }^{\mathrm{2}}{{\displaystyle v}}_{\overline{f}}}{\partial {{\displaystyle \overline{y}}}^{\mathrm{2}}}}}},

yielding

(39)v\left(y\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mu }}{\displaystyle \frac{\partial p}{\partial \overline{x}}}\left({{\displaystyle \overline{y}}}^{\mathrm{2}}-h^{\mathrm{2}}\right)+{{\displaystyle v}}_f,

where {{\displaystyle v}}_f={{\displaystyle \mathit{v}}}_f·{{\displaystyle \mathit{t}}}_{{{\displaystyle \Gamma }}_d} serves as the essential boundary conditions on both sides of the cavity.

Under the assumption of smaller changes in concentrations, the balance of mass for fluid inside cavity reads:

(40){\displaystyle \frac{\partial \rho }{\partial t}}+\rho \nabla ·{{\displaystyle \mathit{u}}}_{\overline{f}}=\mathrm{0.}

Assuming the fluid inside cavity as mono-phase (no mass transfer between cavity and bulk), the balance equation of mass simplifies to

(41)\nabla ·{{\displaystyle \mathit{u}}}_{\overline{f}}=\mathrm{0.}

Obviously, the flow velocity inside the cavity is much higher than in the bulk. The mass balance equation can be rewritten:

(42){\displaystyle \frac{\partial v}{\partial \overline{x}}}+{\displaystyle \frac{\partial w}{\partial \overline{y}}}=\mathrm{0},

where w={{\displaystyle \mathit{v}}}_{\overline{f}}·{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d} denotes the normal velocity of the fluid flow inside the cavity

The mass balance equation is averaged over the width of the cavity:

(43){\displaystyle {\int }_{-h}^h{\displaystyle \frac{\partial v}{\partial \overline{x}}}}\mathrm{d}\overline{y}+{\displaystyle {\int }_{-h}^h{\displaystyle \frac{\partial w}{\partial \overline{y}}}}\mathrm{d}\overline{y}=\mathrm{0.}

The difference of the fluid flow between two sides of the cavity is given by

(44)〚{{\displaystyle w}}_{\overline{f}}〛=w\left(h\right)-w\left(-h\right)=-{\displaystyle {\int }_{-h}^h{\displaystyle \frac{\partial v}{\partial \overline{x}}}}\mathrm{d}\overline{y}.

Substituting Eq. (39) into Eq. (44) gives:

(45)〚{{\displaystyle w}}_{\overline{f}}〛={\displaystyle \frac{\mathrm{2}}{\mathrm{3}\mu }}{\displaystyle \frac{\partial }{\partial \overline{x}}}\left({\displaystyle \frac{\partial p}{\partial \overline{x}}}{h}^{\mathrm{3}}\right)-\mathrm{2}h{\displaystyle \frac{\partial {{\displaystyle v}}_{\overline{f}}}{\partial \overline{x}}}.

This equation describes the total amount of fluid attracted in the tangential flow. It can be included in the weak form for the mass coupling, which ensures the coupling between fluid inside the cavity and fluid in the bulk. Indeed, the coupling term {{\displaystyle n}}_f{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}·〚{{\displaystyle v}}_f-{{\displaystyle v}}_s〛 can be written as:

(46){{\displaystyle n}}_f{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}·〚{{\displaystyle v}}_f-{{\displaystyle v}}_s〛={{\displaystyle n}}_f〚{{\displaystyle w}}_{\overline{f}}-{{\displaystyle w}}_s〛={{\displaystyle n}}_f\left(\frac{\mathrm{2}}{\mathrm{3}\mu }{\displaystyle \frac{\partial }{\partial \overline{x}}}\left({\displaystyle \frac{\partial p}{\partial \overline{x}}}{{\displaystyle h}}^{\mathrm{3}}\right)-\mathrm{2}h{\displaystyle \frac{\partial {{\displaystyle v}}_f}{\partial \overline{x}}}\mathrm{2}\frac{\partial h}{\partial t}\right),

with W_s the normal velocity of the solid skeleton, and the difference between the two sides of the cavity gives:

(47)〚{{\displaystyle W}}_s〛=\mathrm{2}{\displaystyle \frac{\partial h}{\partial t}}\mathrm{.}

Following the Darcy’s law, the tangential velocity reads:

(48){{\displaystyle v}}_{\overline{f}}=\left({{\displaystyle \mathit{v}}}_s-{\displaystyle \frac{{{\displaystyle k}}_f}{{{\displaystyle n}}_f}}\nabla p\right)·{{\displaystyle \mathit{t}}}_{{{\displaystyle \Gamma }}_d}.

Discretization

null

The crack (or cavity) leads to a discontinuous displacement field while the pressure field across the cavity is continuous. Moreover, the spatial derivatives of the pressure orthogonal to the cracks is also discontinuous. Hence, the discretization of the displacement field is given by

{{\displaystyle \mathit{u}}}_s\left(x\right)={\displaystyle {\sum }_{i\in N}{{\displaystyle N}}_i}\left(x\right){{\displaystyle \overline{\mathit{u}}}}_i+{\displaystyle {\sum }_{i\in {{\displaystyle N}}_{cut}}{{\displaystyle N}}_i}\left(x\right){{\displaystyle \mathbf{\aleph }}}_{{{\displaystyle \Gamma }}_d}\left(x\right){{\displaystyle \widehat{\mathit{u}}}}_i

(49)+{\displaystyle {\sum }_{i\in {{\displaystyle N}}_{tip}}{\displaystyle {\sum }_{j\in \left[\mathrm{1},\mathrm{4}\right]}{{\displaystyle N}}_i}}\left(x\right){{\displaystyle \Psi }}_j\left(x\right){{\displaystyle \tilde{u}}}_{i,j},

where N_i are the standard finite element shape functions, {{\displaystyle \overline{u}}}_i, {{\displaystyle \widehat{u}}}_i and {{\displaystyle \tilde{u}}}_{\mathit{i，}j} are nodal parameters and the Heaviside function is defined by

(50){{\displaystyle \mathbf{\aleph }}}_{{{\displaystyle \Gamma }}_d}\left(x\right)=\{\begin{array}{l}\mathrm{0}\mathrm{if}\phi \left(x\right)\le \mathrm{0}\hfill \\ \mathrm{1}\mathrm{i}\mathrm{f}\phi \left(x\right)>\mathrm{0}\hfill \end{array},

where \phi \left(x\right) being the level- set function.

For the node set N_tip, we include the well known crack tip enrichment functions in {{\displaystyle \Psi }}_j:

(51){{\displaystyle \Psi }}_{\mathrm{1}}\left(\mathit{x}\right)=\sqrt{r}\sin {\displaystyle \frac{\theta }{\mathrm{2}}},

(52){{\displaystyle \Psi }}_{\mathrm{2}}\left(\mathit{x}\right)=\sqrt{r}\sin {\displaystyle \frac{\theta }{\mathrm{2}}}\sin \theta ,

(53){{\displaystyle \Psi }}_{\mathrm{3}}\left(\mathit{x}\right)=\sqrt{r}\cos {\displaystyle \frac{\theta }{\mathrm{2}}},

(54){{\displaystyle \Psi }}_{\mathrm{4}}\left(\mathit{x}\right)=\sqrt{r}\cos {\displaystyle \frac{\theta }{\mathrm{2}}}\sin \theta .

These functions depend on a local coordinate system (r, \theta) shown in Fig. 4. Eq. (49) can be rewritten as

(55){{\displaystyle \mathit{u}}}_s=\mathit{N}\mathit{U},

where the matrix N contains the standard and enriched shape functions. Accordingly, the array U includes the displacement for the standard and enriched degrees of freedom {{\displaystyle \overline{u}}}_i,{{\displaystyle \widehat{u}}}_i,{\tilde{u}}_{i,j}.

For the discretization of the pressure, the node set N_press is enriched with the signed distance function {{\displaystyle \mathit{D}}}_{{{\displaystyle \Gamma }}_d}. The enriched distance function {{\displaystyle \mathit{D}}}_{{{\displaystyle \Gamma }}_d} is continues to the discontinues, but its normal derivative is discontinues.

(56){{\displaystyle \mathit{D}}}_{{{\displaystyle \mathrm{\Gamma }}}_d}\left(x\right)=\{\begin{array}{l}-d\mathrm{i}\mathrm{f}\phi \left(x\right)\le \mathrm{0}\hfill \\ d\mathrm{i}\mathrm{f}\phi \left(x\right)>\mathrm{0}\hfill \end{array},

with d the absolute distance to the discontinuity.

The node set N_press includes the nodes affected by the crack. Hence, the discretization of the pressure field can be expressed as

(57)p\left(\mathit{x}\right)={\displaystyle {\sum }_{i\in N}{{\displaystyle H}}_i\left(\mathit{x}\right)}{{\displaystyle \overline{p}}}_i+{\displaystyle {\sum }_{{i\in N}_{press}}{{\displaystyle H}}_i\left(\mathit{x}\right)}{{\displaystyle \mathit{D}}}_{{\Gamma }_d}\left(\mathit{x}\right){{\displaystyle \widehat{p}}}_i,

where H_i is the standard FE shape function for the pressure. We can also rewrite this expression in matrix-vector form by

(58)p=\mathit{H}\mathit{P},

where H contains the standard as well as the enriched shape functions and P contains the degrees of freedom for the pressure {{\displaystyle \overline{p}}}_i,{{\displaystyle \widehat{p}}}_i. The order of the shape function N_i and H_i should be adequate to fulfill the modeling requirements. For the consideration of the consistency of the momentum balance equation, the order of the displacement shape function N_i should be greater than or equal to the order of pressure shape function H_i. This study used the quadrilateral elements equipped with linear shape functions.

Discrete equations

The vector of external force F_ext and external fluid flux Q_ext are given by:

(59){{\displaystyle \mathit{F}}}_{ext}={\displaystyle {\int }_{\Gamma }{\mathit{N}}^{\mathrm{T}}}{{\displaystyle \mathbf{t}}}_p\mathrm{d}\Gamma ,

(60){{\displaystyle \mathit{Q}}}_{ext}={\displaystyle {\int }_{\Gamma }{\mathit{H}}^{\mathrm{T}}}{{\displaystyle \mathit{n}}}_t{{\displaystyle \mathit{q}}}_p\mathrm{d}\Gamma .

Using backward finite difference approximation.

(61)\left({\displaystyle \frac{\partial ·}{\partial t}}\right)i={\displaystyle \frac{{{\displaystyle i^{-}}}_{i-\mathrm{1}}^{\mathrm{.}}}{\Delta t}},

with \Delta x the time increment, {}^{·}i and {}^{·}i-\mathrm{1} denoting the unknown at time step i and i−1, respectively

The coupling force vector F_coupling (on the crack boundary) derives from Eq. (29):

(62){{\displaystyle \mathit{F}}}_{coupling}=-\left({\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}{〚\eta 〛}^{\mathrm{T}}{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}\mathit{H}}\right)\mathit{P}.

Integrating Eq. (31) along {{\displaystyle \Gamma }}_d gives the internal fluid flux Q_coupling:

(63){{\displaystyle Q}}_{coupling}={\displaystyle {\int }_{{{\displaystyle \Gamma }}_d}{\mathit{H}}^{\mathrm{T}}{{\displaystyle \mathit{n}}}_{{{\displaystyle \Gamma }}_d}^{\mathrm{T}}{{\displaystyle \mathit{q}}}_d\mathrm{d}{{\displaystyle \Gamma }}_d}.

The formation of the internal fluid flux Q_coupling is nonlinear, hence, a iteration procedure must be conducted at each time step \Delta t in order to control the accuracy of the resolution. The iteration residual vector Rⁱ at i^th iteration is defined:

（64）{\mathit{R}}^i=\left[\begin{array}{ll}\mathrm{0}\hfill & \mathrm{0}\hfill \\ {{\displaystyle K}}_{up}^T\hfill & {{\displaystyle K}}_{pp}^{\left(\mathrm{1}\right)}\hfill \end{array}\right]{{\displaystyle \left(\begin{array}{l}\Delta U\hfill \\ \Delta P\hfill \end{array}\right)}}^i+\left[\begin{array}{ll}{{\displaystyle K}}_{uu}\hfill & {{\displaystyle K}}_{up}\hfill \\ \mathrm{0}\hfill & \Delta t{{\displaystyle K}}_{pp}^{\left(\mathrm{2}\right)}\hfill \end{array}\right]{{\displaystyle \left(\begin{array}{l}U\hfill \\ P\hfill \end{array}\right)}}^i+{{\displaystyle \left(\begin{array}{l}{{\displaystyle F}}_{coupling}\hfill \\ {{\displaystyle Q}}_{coupling}\hfill \end{array}\right)}}^i-\left(\begin{array}{l}{{\displaystyle F}}_{ext}\hfill \\ {{\displaystyle Q}}_{ext}\hfill \end{array}\right),

with stiffness matrix:

(65){{\displaystyle \mathit{K}}}_{uu}={\displaystyle {\int }_{\mathrm{\Omega }}{\mathit{B}}^{\mathrm{T}}}{\mathit{D}}^{\mathit{tan}}\mathit{B}\mathrm{d}\mathrm{\Omega },

(66){{\displaystyle \mathit{K}}}_{up}=-{\displaystyle {\int }_{\mathrm{\Omega }}{\alpha \mathit{B}}^{\mathrm{T}}}\mathit{m}\mathit{H}\mathrm{d}\mathrm{\Omega },

(67){{\displaystyle \mathit{K}}}_{pp}^{\left(\mathrm{1}\right)}=-{\displaystyle {\int }_{\mathrm{\Omega }}{Q^{-\mathrm{1}}\mathit{H}}^{\mathrm{T}}}\mathit{H}\mathrm{d}\mathrm{\Omega },

(68){{\displaystyle \mathit{K}}}_{pp}^{\left(\mathrm{2}\right)}=-{\displaystyle {\int }_{\mathrm{\Omega }}{{{\displaystyle k}}_f\nabla \mathit{H}}^{\mathrm{T}}\nabla }\mathit{H}\mathrm{d}\mathrm{\Omega },

where m=[1,1,0] for two dimensional problem.

In the Newtown-Raphson algorithm, the iteration matrix Kⁱ derived from:

(69)K^i={\displaystyle \frac{\partial f}{\partial ·}},

with f and · denote the residual function and unknowns, respectively. In this study, the iteration matrix Kⁱ has the form:

(70)K^i=\left[\begin{array}{ll}{{\displaystyle K}}_{uu}\hfill & {{\displaystyle K}}_{up}+{\displaystyle \frac{\partial {{\displaystyle \mathit{F}}}_{coupling}}{\partial P}}\hfill \\ {{\displaystyle K}}_{up}^T+\Delta t{\displaystyle \frac{\partial {{\displaystyle \mathit{Q}}}_{coupling}}{\partial U}}\hfill & {{\displaystyle K}}_{pp}^{\left(\mathrm{1}\right)}+\Delta t{{\displaystyle K}}_{pp}^{\left(\mathrm{2}\right)}+{\displaystyle \frac{\partial {{\displaystyle \mathit{Q}}}_{coupling}}{\partial P}}\hfill \end{array}\right],

with all items evaluated at iteration i.

The coupling term F_coupling and Q_coupling cause the Jacobian matrix of the residual Rⁱ to become asymmetric, to regain the symmetric the coupling terms are omitted in the Jacobian matrix. This may decrease the convergence rate of the Newtown-Raphson algorithm. Nevertheless, the symmetric matrix allows flexible implementation as well as better condition of the matrix structure. The simplified Jacobian matrix reads:

(71)K^i=\left[\begin{array}{ll}{{\displaystyle K}}_{uu}\hfill & {{\displaystyle K}}_{up}+{\displaystyle \frac{\partial {{\displaystyle \mathit{F}}}_{coupling}}{\partial P}}\hfill \\ {{\displaystyle K}}_{up}^T\hfill & {{\displaystyle K}}_{pp}^{\left(\mathrm{1}\right)}+\Delta t{{\displaystyle K}}_{pp}^{\left(\mathrm{2}\right)}\hfill \end{array}\right].

Example calculation

Consider a two-dimensional specimen under plane-strain conditions as depicted in Fig. 5. A normal fluid flux q₀=10⁻⁴ m·s⁻¹ starting at t = 0 s is imposed at the bottom face while the top face is assigned a drained condition with zero pressure. Both left and right faces have undrained boundary conditions. No mechanical load is applied, but essential boundary conditions have been applied. The block is 10 m × 10 m and consists of a porous material with a fluid volume fraction of n_f= 0.3. The absolute mass densities are {{\displaystyle \rho }}_s^{\text{'}}={{\displaystyle \rho }}_c/{{\displaystyle n}}_c = 2000 kg/m³ for the solid phase and {{\displaystyle \rho }}_f^{\text{'}}={{\displaystyle \rho }}_f/{{\displaystyle n}}_f= 1000 kg/ m³ for the fluid phase. A linear elastic isotropic solid is assumed with a Young’s modulus of E = 9 GPa and a Poisson ratio of n = 0.4. The Biot coefficient a has been set equal to 1, and the Biot modulus has been assigned a value Q = 10¹⁸ GPa so as to simulate a quasi-incompressible fluid. This is not a limitation of the model, but shows the influence of cracks more clearly. The bulk material has a permeability k_f = 10⁻⁹ m³/(N·s) while the viscosity of the fluid is m = 10⁻³ N/ (m²·s). One cavity is included in the model, the length of the cavity is 3 m, the inclination angle regarding to the x direction is 30°. The total time duration of the process is 10 s. The whole fluid flow process in the specimen can be summarized as following: Due to the applied fluid flow flux on the bottom edge of the domain, the fluid starts to flow into the bulk. The flow direction is normal to the bulk. After the fluid enters the cavity, it flows inside the cavity tangentially to the crack until it reaches the crack tips and reenters the bulk. In the cavity, the velocity of the fluid is much higher compared to the velocity in the bulk because there is no resisting solid skeleton. The pressure increases in the bulk and the cavity in time. Fig. 6 shows the fluid pressure and stress profile in x direction at 10 s.

As shown in Fig. 7, an edge crack with length of 0.05 m is located in the middle of the plate. Tensile loading is imposed by applying opposite velocities on the top and bottom edges. The boundary of plate is assumed to be impervious. Figure 8 shows the pressure distribution and displacement field.

Conclusions

We have presented an extended finite element method (XFEM) in porous media for fluid flow through discrete cracks. An incompressible Newtonian fluid is assumed and the material in the bulk is assumed to be isotropic and linear elastic. The crack is modeled by the Heaviside and crack tip enrichment for linear elastic fracture mechanics (LEFM) in the context of partition of unity enrichment. An abs-enrichment is chosen for the pressure field which is continuous across the crack but its gradient orthogonal to the crack is discontinuous. A monolithic approach has been proposed where non-linearities are resolved by the Newton method with line search. In the future, we aim to extend this approach to more complex material behavior, such as: fracture propagation under the influence of the pressurized voids and three-dimensions model with the final goal to model hydraulic fracturing for shale gas development.

References

[1] Bažant Z P. Size effect in blunt fracture: concrete, rock, metal. Journal of Engineering Mechanics, 1984, 110(4): 518–535

[2] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131–150

[3] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[4] Wang H, Belytschko T. Fluid–structure interaction by the discontinuous-galerkin method for large deforma-tions. International Journal for Numerical Methods in Engineering, 2009, 77(1): 30–49

[5] Gerstenberger A, Wall W A. An extended finite element method/lagrange multiplier based approach for fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering, 2008, 197(19‒20): 1699–1714

[6] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 22: 48

[7] RéthoréJ, Borst R, Abellan M A. A two-scale approach for fluid flow in fractured porous media. International Journal for Numerical Methods in Engineering, 2007, 71(7): 780–800

[8] Song J H, Areias P, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[9] Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(33‒35): 3523–3540

[10] Rabczuk T, Areias P, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[11] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee C K, Zi G, Zhuang X, Liu G R, Rabczuk T. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. Journal of Applied Mathematics, 2013, 2013: 978026

[12] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92–93: 242–256

[13] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[14] Sukumar N, Chopp D L, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite-element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190(46‒47): 6183–6200

[15] Areias P, Rabczuk T, Camanho P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[16] Areias P, Dias-da Costa D, Sargado J, Rabczuk T. Element-wise algorithm for modeling ductile fracture with the rousselier yield function. Computational Mechanics, 2013, 52(6): 1429–1443

[17] Areias P, Rabczuk T, Dias-da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[18] Areias P, Rabczuk T, Camanho P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[19] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[20] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal for Numerical Methods in Engineering, 1999, 44(9): 1267–1282

[21] Remmers J, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. Computational Mechanics, 2003, 31(1‒2): 69–77

[22] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1–4): 3–47

[23] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[24] Dolbow J, Belytschko T. Numerical integration of the galerkin weak form in meshfree methods. Computational Mechanics, 1999, 23(3): 219–230

[25] Belytschko T, Guo Y, Liu W K, Xiao S P. A unified stability analysis of meshless particle methods. International Journal for Numerical Methods in Engineering, 2000, 48(9): 1359–1400

[26] Belytschko T, Tabbara M. Dynamic fracture using element-free galerkin methods. International Journal for Numerical Methods in Engineering, 1996, 39(6): 923–938

[27] Organ D, Fleming M, Terry T, Belytschko T. Continuous meshless approximations for nonconvex bodies by diffraction and transparency. Computational Mechanics, 1996, 18(3): 225–235

[28] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[29] Rabczuk T, Areias P. A meshfree thin shell for arbitrary evolving cracks based on an extrinsic basis. Computer Modeling in Engineering & Sciences, 2006, 16: 115–130

[30] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367–382

[31] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[32] Liew K M, Zhao X, Ferreira A J M. A review of meshless methods for laminated and functionally graded plates and shells.  Composite Structures, 2011,  93(8): 2031–2041

[33] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[34] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75: 943–960

[35] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[36] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273–1311

[37] Areias P, Msekh M, Rabczuk T. Damage and fracture algorithm using the screened poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[38] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[39] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29‒30): 2777–2799

[40] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[41] Silling S A, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88(2): 151–184

[42] Munjiza A, Andrews K, White J. Combined single and smeared crack model in combined finite-discrete element analysis. International Journal for Numerical Methods in Engineering, 1999, 44(1): 41–57

[43] Scholtès L, Donzé F V. Modelling progressive failure in fractured rock masses using a 3D discrete element method. International Journal of Rock Mechanics and Mining Sciences, 2012, 52: 18–30

[44] Tan Y, Yang D, Sheng Y. Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline sic. Journal of the European Ceramic Society, 2009, 29(6): 1029–1037

[45] Kulatilake P, Malama B, Wang J. Physical and particle flow modeling of jointed rock block behavior under uniaxial loading. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(5): 641–657

[46] Zeng Q y, Zhou J. Analysis of passive earth pressure due to various wall movement by particle flow code (2D). Rock and Soil Mechanics, 2005, 26: 43–47

[47] Wang C, Tannant D. Rock fracture around a highly stressed tunnel and the impact of a thin tunnel liner for ground control. International Journal of Rock Mechanics and Mining Sciences, 2004, 41: 676–683

[48] Jing L, Ma Y, Fang Z. Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(3): 343–355

[49] Shi G H. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of de-formable block structures. Engineering Computations, 1992, 9(2): 157–168

[50] Pearce C, Thavalingam A, Liao Z, Bićanić N. Computational aspects of the discontinuous deformation analysis framework for modelling concrete fracture. Engineering Fracture Mechanics, 2000, 65(2–3): 283–298

[51] Ning Y, Yang J, An X, Ma G. Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework. Computers and Geotechnics, 2011, 38(1): 40–49

[52] Abellan M A, De Borst R. Wave propagation and localisation in a softening two-phase medium. Computer Methods in Applied Mechanics and Engineering, 2006, 195(37–40): 5011–5019

[53] Lewis R W, Schrefler B A. The finite element method in the deformation and consolidation of porous media. New York: John Wiley and Sons Inc., 1987

[54] Hassanizadeh S M, Gray W G. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Advances in Water Resources, 1990, 13(4): 169–186

[55] Darcy H. The public fountains of the city of Dijon: exhibition and application. Victor Dalmont, 1856 (in French)

[56] Bachelor G. An Introduction to Fluid Mechanics. Cambridge University Press, 1967

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature