Two-scale approaches for fracture in fluid-saturated porous media

. A derivation is given of two-scale models that are able to describe deformation and flow in a fluid-saturated and progressively fracturing porous medium. From the micromechanics of the flow in the cavity, identities are derived that couple the local momentum and the mass balances to the governing equations for a fluid-saturated porous medium, which are assumed to hold on the macroscopic scale. By exploiting the partition-of-unity property of the finite element shape functions, the position and direction of the fractures are independent from the underlying discretization. The finite element equations are derived for this two-scale approach and integrated over time. The resulting discrete equations are nonlinear due to the cohesive crack model and the nonlinearity of the coupling terms. A consistent linearization is given for use within a Newton-Raphson iterative procedure. Finally, examples are given to show the versatility and the efficiency of the approach.


Introduction
Since the pioneering work of Terzaghi (1943) and Biot (1965) the flow of fluids in deforming porous media has received considerable attention.Recently, Lewis and Schrefler (1998) have given an account of this topic which is crucial for understanding and predicting the physical behaviour of many systems of interest, for example, in geotechnical and petroleum engineering, but also for soft tissues.Because of the complicated structure and functioning of human tissues, the classical two-phase theory has been extended to three and four-phase media, taking into account ion transport and electrical charges (Huyghe and Janssen 1997, Snijders et al. 1997, Van Loon et al. 2003).A general framework for accommodating multi-field problems has been presented by Jouanna and Abellan (1995).
In spite of the importance of the subject, flow in damaged porous media has received little attention.Yet, the presence of damage, such as cracks, faults, and shear bands, can markedly change the physical behaviour.Furthermore, the fluid can transport contaminants which can dramatically reduce the strength of the solid skeleton.To account for such phenomena, the fluid flow must be studied also in the presence of discontinuities in the solid phase.The physics of the flow within such discontinuities can be very different from that of the interstitial fluid in the deforming bulk material.These differences affect the flow pattern and therefore also the deformations in the vicinity of the discontinuity.As we will show at the end of the paper, the local differences in flow characteristics can even influence the flow and deformations in the entire body of interest.
In this contribution, we will describe a general numerical methodology to capture deformation and flow in progressively fracturing porous media, summarising and unifying the recent work reported in de Borst et al. (2006), Réthoré et al. (2007a, 2007b, 2007c).The model allows for flow inside an evolving crack to be in the tangential direction.This is achieved by a priori adopting a two-scale approach.At the fine scale the flow in the cavity created by the (possibly cohesive) crack can be modelled in various ways, e.g., as a Stokes flow in an open cavity, or using a Darcy relation for a damaged porous material.Since the crosssectional dimension of the cavity is small compared to its length, the flow equations can be averaged over the width of the cavity.The resulting equations provide the momentum and mass couplings to the standard equations for a porous medium, which are assumed to hold on the macroscopic scale.
Numerically, the two-scale model which ensues, imposes some requirements on the interpolation of the displacement and pressure fields near the discontinuity.The displacement field must be discontinuous across the cavity.Furthermore, the micromechanics of the flow within the cavity require that the flow normal to the cavity is discontinuous, and in conformity with Darcy's relation which, at the macroscopic scale, is assumed to hold for the surrounding porous medium, the normal derivative of the fluid pressure field must also be discontinuous from one face of the cavity to the other.For arbitrary discretisations, these requirements can be satisfied by exploiting the partition-ofunity property of finite element shape functions (Babuska & Melenk 1997), as has been done successfully in applications to cracking in single-phase media (Belytschko & Black 1999, Moës et al. 1999, Wells & Sluys 2001, Wells et al. 2002a, 2002b, Remmers et al. 2003, Samaniego & Belytschko 2005, Réthoré et al. 2005a, 2005b, Areias & Belytschko 2006, de Borst et al. 2006).
To provide a proper setting, we will first briefly recapitulate the governing equations for a deforming porous medium under quasi-static loading conditions.The strong as well as the weak formulations will be considered, since the latter formulation is crucial for incorporating the micromechanical flow model properly.This micromechanical flow model is discussed in the next section, where it will be demonstrated how the momentum and mass couplings of the micromechanical flow model to the surrounding porous medium can be accomplished in the weak formulation.Time integration and a consistent linearisation of the resulting equations, which are nonlinear due to the coupling terms and the cohesive crack model complete the numerical model.Finally, example calculations are given of a body with stationary and propagating cracks.

Balance equations
We consider a two-phase medium subject to the restriction of small displacement gradients and small variations in the concentrations (Jouanna & Abellan 1995).Furthermore, the assumptions are made that there is no mass transfer between the constituents and that the processes which we consider, occur isothermally.With these assumptions, the balances of linear momentum for the solid and the fluid phases read: (1) with σ π the stress tensor, ρ π the apparent mass density, and v π the absolute velocity of constituent π.As in the remainder of this paper, π = s, f, with s and f denoting the solid and fluid phases, respectively.Further, g is the gravity acceleration and is the source of momentum for constituent π from the other constituent, which takes into account the possible local drag interaction between the solid and the fluid.Evidently, the latter source terms must satisfy the momentum production constraint: (2) We now neglect convective, gravity and acceleration terms, so that the momentum balances reduce to: (3) Adding both momentum balances, and taking into account Eq. ( 2), one obtains the momentum balance for the mixture: (4) where the stress is, as usual, composed of a solid and a fluid part, σ = σ s + σ f .
In a similar fashion as for the balances of momentum, one can write the balance of mass for each phase as: (5) Again neglecting convective terms, the mass balances can be simplified to give: We multiply the mass balance for each constituent π by its volumetric ratio n π , add them and utilise the constraint (7) to give: The change in the mass density of the solid material is related to its volume change by: -------= with K s the bulk modulus of the solid material and K t the overall bulk modulus of the porous medium.Using the Biot coefficient, α =1− K t / K s (Lewis & Schrefler 1998), this equation can be rewritten as (10) For the fluid phase, a phenomenological relation is assumed between the incremental changes of the apparent fluid mass density and the fluid pressure p (Lewis & Schrefler 1998): (11) with the overall compressibility, or Biot modulus (12) where K f is the bulk modulus of the fluid.Inserting relations (10) and ( 11) into the balance of mass of the total medium, Eq. ( 8), gives: (13) The field equations, i.e., the balance of momentum of the saturated medium, Eq. ( 4), and the balance of mass, Eq. ( 13), are complemented by the boundary conditions n Γ • σ = t p , v s = v p (14) which hold on complementary parts of the boundary Ω t and Ω v , with Γ = Ω = Ω t ∪Ω v , Ω t Ω v =, t p being the prescribed external traction and v p the prescribed velocity, and which hold on complementary parts of the boundary Ω q and Ω p , with Γ = Ω = Ω q ∪Ω p and Ω q Ω p =, q p and p p being the prescribed outflow of pore fluid and the prescribed pressure, respectively.The initial conditions specify the displacements u π , the velocities v π , and the pressure field at t =0, u π (x, 0) = , v π (x, 0) = , p (x, 0) = p 0 (16) 3 Constitutive equations The effective stress increment in the solid skeleton, dσ ' s is related to the strain increment dε s by an incrementally linear stress-strain relation for the solid skeleton, dσ ' s = tan :dε s (17) where tan is the fourth-order tangent stiffness tensor of the solid material and the d-symbol denotes a small increment.Since the effective stress in the solid skeleton is related to the partial stress by σ ' s = σ s / n s , the above relation can be replaced by α 1 - where the notation D tan = n s tan has been used.In the examples, a linear-elastic behaviour of the bulk material has been assumed, and we have set D tan = D, the linear-elastic stiffness tensor.
At the discontinuity Γ d a discrete relation holds between the interface tractions t d and the relative displacements δ : ) with κ a history parameter.After linearisation, necessary to use a tangential stiffness matrix in an incremental-iterative solution procedure, one obtains: with T the material tangent stiffness matrix of the discrete traction-separation law: A key element is the presence of a work of separation or fracture energy, G c , which governs crack growth and enters the interface constitutive relation ( 19) in addition to the tensile strength f t .It is defined as the work needed to create a unit area of fully developed crack: (22) with σ the stress across the fracture process zone.It thus equals the area under the decohesion curves as shown in Fig. 1.Evidently, cohesive-zone models as defined above are equipped with an internal length scale, since the quotient G c / E, with E a stiffness modulus for the surrounding continuum, has the dimension of length.
In a standard manner we adopt Darcy's relation for the flow in the bulk of the porous medium, (23) with k f the permeability coefficient.The equations for the flow in the cavity, which close the initial value problem, will be detailed in a subsequent section.

Weak form of the balance equations
To arrive at the weak form of the balance equations, we multiply the momentum balance (4) and the mass balance (13) by kinetically admissible test functions for the displacements of the skeleton, η, and for the pressure, ζ.After substitution of Darcy's relation for the flow in the porous medium, Eq. ( 23), integrating over the domain Ω and using the divergence theorem then leads to the corresponding weak forms: (24) and ( 25) Because of the presence of a discontinuity inside the domain Ω, the power of the external tractions on Γ d and the normal flux through the faces of the discontinuity are essential features of the weak formulation.Indeed, these terms enable the momentum and mass couplings between the discontinuitythe microscopic scale-and the surrounding porous medium -the macroscopic scale.
The momentum coupling stems from the tractions across the faces of the discontinuity and the pressure applied by the fluid in the discontinuity onto the faces of the discontinuity.We assume stress continuity from the cavity to the bulk, so that we have: with t d the cohesive tractions, which are given by Eq. ( 19).Therefore, the weak form of the balance of momentum becomes: Since the tractions have a unique value across the discontinuity, the pressure p must have the same value at both faces of the discontinuity, and, consequently, this must also hold for the test function for the pressure, ζ.Accordingly, the mass transfer coupling term for the water can be rewritten as follows: (28) where q d = n f [v f − v s ] represents the difference in the fluid fluxes through both faces of the discontinuity.
The above identity for the coupling of the mass transfer can be interpreted as follows.Part of the fluid that enters the cavity through one of its faces flows away tangentially, that is in the cavity.Therefore, the fluid flow normal to the cavity is discontinuous.Because the fluid flow between the cavity and the surrounding porous medium has to be continuous at each of the faces of the discontinuity, and because the fluid velocity is related to the pressure gradient via Darcy's law, the gradient of the pressure normal to the discontinuity must be discontinuous across the crack.

Micro-macro coupling
To quantify the influence of the 'micro'-flow inside the discontinuity on the 'macro'-scale, the balances of mass and momentum at the microscale must be considered.Different assumptions can now be made for the cavity.First, we will assume an open cavity which is totally filled with a Newtonian fluid.Then, the balance of mass for the 'micro'-flow in the cavity reads: subject to the assumptions of small changes in the concentrations and that convective terms can be neglected, cf., Eq. ( 6).We assume that the first term can be neglected because the problem is monophasic in the cavity and the velocities are therefore much higher than in the porous medium.With this assumption, and focusing on a two-dimensional configuration, the mass balance inside the cavity simplifies to: and w = v • the tangential and normal components of the fluid velocity in the discontinuity, respectively, see Fig. 2. Accordingly, the difference in the fluid velocity components that are normal to both crack faces is given by: (30) To proceed, the velocity profile of the fluid flow inside the discontinuity must be known.From the balance of momentum for the fluid in the cavity and the assumption of a Newtonian fluid, the following velocity profile results: where an integration has been carried out from y = −h to y = h and µ the viscosity of the fluid.The essential boundary v = v f has been applied at both faces of the cavity, and stems from the relative fluid velocity in the porous medium at y = ± h: 31) into Eq.( 30) and again integrating with respect to y then leads to: (32) This equation gives the amount of fluid attracted in the tangential fluid flow.It can be included in the weak form of the mass balance of the 'macro'-flow to ensure the coupling between the 'micro'flow and the 'macro'-flow.Since the difference in the normal velocity of both crack faces is given by the mass coupling term becomes: Another possibility for closure of the initial value problem is to assume that the cavity is partially filled with solid material.The mass balance for the fluid inside the cavity then reads: Because the width of the cavity 2h is negligible compared to its length, the mass balance is again enforced in an average sense over the cross section.For the first term, we obtain: (34) where v s and w s are the component of the solid velocity tangential and normal to the crack, respectively.Because α depends on the capillarity pressure only, it can be assumed as constant over the cross section.Assuming furthermore that v s varies linearly with y, and defining <v s > =1/2(v s (h)+ v s (−h)), the integral can be solved analytically: (35) Applying the same operations to the second term of Eq. ( 30), the following expression ensues: − (36) We now introduce Darcy's relation projected onto the axis tangential to the crack, (37) with k d the permeability of the damaged, porous material inside the cavity.The dependence of the permeability inside the cavity, k d , on y will be negligible compared to that on x and accordingly, k d can be assumed not to depend on y.However, the decohesion inside the cavity will affect the permeability, and therefore k d = k d (h).Upon substitution of Eq. ( 37) into Eq.( 36), the following relation is obtained: At the faces of the cavity the permeability equals that of the bulk, and for continuity reasons = , so that the second term becomes: (39) For the third term, neglecting variations of the pressure over the height of the cavity, one obtains: (40) Recalling that the term n f [w f − w s ] can be identified as the coupling term • q d of the weak form of the mass balance, we obtain: (41)

Discontinuities in a two-phase medium
A finite element method that can accommodate the propagation of discontinuities through 3 Body composed of continuous displacement fields at each side of the discontinuity Γ d elements was proposed by Belytschko and co-workers (Belytschko & Black 1999, Möes et al. 1999), exploiting the partition-of-unity property of finite element shape functions (Babuska & Melenk 1997).Since finite element shape functions ϕ j form partitions of unity, ϕ j =1 with n the number of nodal points, the components v i of a velocity field v can be interpolated as ( 42) with a j the 'regular' nodal degrees-of-freedom for the displacements, ψ k the enhanced basis terms, and the additional displacement degrees-of-freedom at node j which represent the amplitude of the kth enhanced basis term ψ k .Next, we consider a domain Ω that is crossed by a single discontinuity at Γ d (see Fig. 3).The velocity field v s can then be written as the sum of two continuous velocity fields s and : where H Γ d is the Heaviside step function centred at the discontinuity.The decomposition in eq. ( 43) has a structure similar to the interpolation in Eq. ( 42).Accordingly, the partition-of-unity property of finite element shape functions enables the direct incorporation of discontinuities, including cracks and shear bands, in finite element models such that the discontinuous character of cracks and shear bands is preserved.With the standard small-strain assumption that the strain-rate field of the solid, s , is derived from the symmetric part of the gradient of the velocity field, we obtain: away from the discontinuity Γ d , with the superscript s denoting the symmetric part of the gradient operator.
With respect to the interpolation of the pressure p we note that the fluid flow normal to the discontinuity can be discontinuous.Since the fluid velocity is related to the pressure gradient via Darcy's relation, the gradient of the pressure normal to the discontinuity can therefore also be discontinuous across the discontinuity.Accordingly, the enrichment of the interpolation of the pressure must be such that the pressure itself is continuous, but has a discontinuous normal first spatial derivative.The distance function D Γ d defined as 45) satisfies this requirement, and accordingly, the interpolation of p will be such that: We now discretise the fields v s and p and the test functions η and ζ for the velocity and the pressure, respectively, in a Bubnov-Galerkin sense: (47) where the matrix N contains the shape functions N i used as partition of unity for the interpolation of the velocity field v s , and H contains the shape functions H i used as partition of unity for the interpolation of the pressure field p. and are the nodal arrays assembling the amplitudes that correspond to the standard and enhanced interpolations of the velocity field, and and assemble the amplitudes that correspond to the standard and enhanced interpolations of the pressure field.The choice for the interpolants in N and H is driven by modelling requirements.Indeed, the modelling of the fluid flow inside the cavity requires the computation of second derivatives of the pressure, see Eqs. ( 33) or ( 41).Hence, the order of the finite element shape functions H i has to be sufficiently high, otherwise the coupling between the fluid flow in the cavity and the bulk will not be achieved.
Further, the order of the finite element shape functions N i must be greater than or equal to that of H i for consistency in the discrete momentum balance.Inserting Eqs. ( 47) into Eqs.( 27) and ( 28) and requiring that the result holds for all admissible , , and gives: (48) with the external force vectors: (49) and the internal force vectors: (50) where, for two dimensions, m T = [1, 1, 0].
() ∆t ----------== respectively.Substitution of this identity into the set of semi-discrete equations (48) results in: (52) where while the other internal force vectors remain unchanged except that they are now evaluated at t + ∆t.
For use in a Newton-Raphson solution method, the set (52) must be linearised.To this end, the stress and the pressure are decomposed as σ j = σ j − 1 +dσ, p j = p j − 1 +dp (53) with the subscripts j − 1 and j signifying the iteration numbers.The linearisation of the set ( 52  The matrices , and result from differentiating the mass coupling term with respect to , and , respectively.The resulting expressions depend on the subscale model that has been chosen for the microflow, and are generally very complicated.In any case, these coupling terms and the coupling term cause the tangent stiffness matrix to become unsymmetric.To restore symmetry, the non-symmetric contributions that arise from the mass and momentum couplings have been removed from the tangent stiffness matrix, and the iterations in the example calculations in the next section have been carried out with the resulting symmetric stiffness matrix.

Examples: stationary and propagating cracks
First, an open, stationary crack is considered in a linear-elastic medium and no tractions are transferred between both sides of the crack.Hence, the first subscale model for the fluid flow in the cavity is used and singularity functions that stem from linear-elastic fracture mechanics are added as enrichment functions at the crack tips (Réthoré et al. 2007b).The specimen is a square-shaped fractured block under plane-strain conditions.A normal fluid flux q o =10 −4 ms −1 starting at t =0 s is imposed at the bottom, 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 in order to remove rigid body motions.The block is 10 m × 10 m and consists of a porous material with a fluid volume fraction n f = 0.3.The absolute mass densities are ρ' s = ρ s / n s =2000kg/m 3 for the solid phase and ρ' f = ρ f / n f = 1000 kg/m 3 for the fluid phase.The solid constituent is assumed to behave in a linear elastic manner with a Young's modulus E = 9 GPa and a Poisson's ratio ν = 0.4.The Biot coefficient α has been set equal to 1, and the Biot modulus has been assigned a value Q =10 18 GPa so as to simulate a quasiincompressible fluid.This is not a limitation of the model, but more clearly brings out the influence of a fault.The bulk material has a permeability k f =10 −9 m 3 /Ns.Initially, one fault is considered.This fault, at the centre of the specimen, is 2 m long, and is inclined with the horizontal axis.A reference simulation has been run for a period of 10 s using 75 time steps.For this reference simulation, the coarse mesh shown in Fig. 4 has been used with quadrilateral elements equipped with quadratic shape functions.The 'macro'-flow is oriented from the bottom to the top of the specimen.The 'micro'-flow in the cavity is oriented by the direction of the fault.In the cavity, the velocity of the fluid is very high compared with the velocity in the bulk because there is no resisting solid skeleton.
Next, a set of ten stationary faults is randomly generated.The length of the faults is between 1 and 3 m, while the fault angles vary from 10 o to +30 o .Fig. 8 shows the influence of the faults on the norm of the pressure gradient.The global fluid flow is strongly affected by the 'micro'-flows inside the faults.From Figure it is observed that the main effect is due to the two longest faults.Evidently, the effect of a fault on the 'macro'-flow increases with its length because more fluid can flow inside the cavity.
Finally, a cohesive crack is simulated which propagates in a plate under plane-strain conditions.Now, the second model for fluid flow inside the cavity has been used, in which progressive decohesion within the crack is assumed (Réthoré et al. 2007c).The plate has sides of 0.25 m and an initial notch which is located at the symmetry axis and is 0.05 m deep.A fixed vertical velocity v = 2.35 × 102 µm/s is prescribed in a opposite direction at the bottom and at the top of the plate (tensile loading).All other boundaries of the plate are assumed to be impervious.The mesh consists of 20 × 20 quadrilateral elements with bilinear shape functions for the pressure and the displacement fields.The time step size is 1 s and the analysis is continued until the crack reaches the right side of the plate.
Obviously, the amount of fluid that flows into the cavity is closely related to the crack opening displacement and vanishes at the crack tip.As a consequence the water pressure decreases in the vicinity of the crack.Fig. 9 gives the pressure fields for a simulation with a full coupling pling at the interface as derived in the previous section and for a simulation without a mass transfer coupling term (and also without pressure enrichment).In the latter case, the crack is not seen as a discontinuity in the pressure field and the fluid phase flows through the crack as it does in the bulk.Accordingly, the pressure and pressure gradient are continuous at the interface.Indeed, the results presented in the two graphs of Fig. 9 are very different.Because of the mass transfer coupling, the water is sucked into the crack and high negative pressures occur around the crack.As a consequence, the water saturation decreases in this zone and intense cavitation occurs.Because of the negative values of the water pressure, sucking tractions modify the global response of the plate.As shown in Fig. 10, the load-displacement curve obtained with the coupling term results in a higher load-carrying capacity.Indeed, the effect of the mass transfer coupling is strong and changes the fluid flow in the entire domain.Moreover, the fluid velocity is increased by an order of magnitude.Finally, the profile of the tangential component of the Darcian velocity is shown in Fig. 11.The graphs are plotted such that the velocity is positive when it is oriented from the actual crack tip to the initial pre-notch.When the coupling is activated, the water flows from the actual tip to the initial notch with a value depending on the crack opening displacement whereas the orientation of the flow in the uncoupled case depends on the position on the interface.

Concluding remarks
A methodology has been proposed to insert discontinuities such as cracks, faults, or shear bands, in a porous medium.The discontinuities can be located arbitrarily, not related to the underlying discretisation.For the fluid flow in the fractured porous medium a two-scale approach has been chosen, where the flow of the fluid inside the discontinuity (the 'micro'-scale) is modelled independently from the flow of the pore fluid in the surrounding porous medium (the 'macro'scale).The mechanical and the mass transfer couplings between the two scales are obtained by inserting the homogenised 'constitutive' relations of the 'micro'-flow into the weak form of the balance equations of the bulk.The assumptions made for the fluid flow in and near the discontinuity require the addition of special enrichment functions for the displacement and the pressure fields.These conditions are satisfied by exploiting the partition of unity property of the finite element polynomial shape functions.

Fig. 2
Fig. 2 Geometry and local coordinate system in cavity

Fig. 4
Fig. 4 Fine and coarse mesh composed of quadrilateral elements

Fig. 6
Fig. 6 Normal derivative of the pressure for a crack angle of 30 o .The left picture, where the data have been plotted as a piecewise-constant field, shows the results for the coarse mesh.Right, the results for the fine mesh have been smoothed and are depicted with filled iso-values n Γ d ∇p ⋅

Fig. 8
Fig. 8 L 2 -Norm of the pressure gradient

Fig. 9
Fig. 9 Pressure field in Pa.Left: The case with full coupling; Right: The case without coupling