An integrated numerical model for coupled poro-hydro-mechanics and fracture propagation using embedded meshes
AA N INTEGRATED NUMERICAL MODEL FOR COUPLEDPORO - HYDRO - MECHANICS AND FRACTURE PROPAGATIONUSING EMBEDDED MESHES
A P
REPRINT
Guotong Ren
Department of Petroleum EngineeringUniversity of TulsaTulsa, OK 740104 [email protected]
Rami M. Younis
Department of Petroleum EngineeringUniversity of TulsaTulsa, OK 74104 [email protected]
August 10, 2020 A BSTRACT
Integrated models for fluid-driven fracture propagation and general multiphase flow in porous mediaare valuable to the study and engineering of several systems, including hydraulic fracturing, under-ground disposal of waste, and geohazard mitigation across such applications. This work extendsthe coupled model multiphase flow and poromechanical model of [1] to admit fracture propagation(FP). The coupled XFEM-EDFM scheme utilizes a separate fracture mesh that is embedded on astatic background mesh. The onset and dynamics of fracture propagation (FP) are governed by theequivalent stress intensity factor (SIF) criterion. A domain-integral method (J integral) is appliedto compute this information. An adaptive time-marching scheme is proposed to rapidly restrict andgrow temporal resolution to match the underlying time-scales. The proposed model is verified withanalytical solutions, and shows the capability to accurately and adaptively co-simulate fluid transportand deformation as well as the propagation of multiple fractures. K eywords Fracture propagation · Coupled hydro-mechanics · Porous media · Extended finite element method · Embedded discrete fracture method
Numerical models are an important enabling technology towards the advancement of a number of engineered systemswithin the nexus of energy, water, and the environment (e.g. subsurface energy resource extraction, waste disposalor storage, and geological intermittent energy storage systems). The first-order response of such systems is oftendriven by coupled poro-thermo-hydro-mechanics, including fracturing and multiphase flow. While these underlyingprocesses typically occur with local spatial and temporal support within the extent of the system-scale, their long-range interactions and causal relations cannot be ignored. Subsequently, ubiquitous (across engineering subprocesses)models are necessary in the context of joint assimilation of multiple types of observations (e.g. displacement, fluidflow, and thermodynamic state) acquired over the duration of multiple operational substeps (e.g. hydraulic fracturing,injection, and production). The focus of this work is the development of an efficient, robust, and ubiquitous numericalapproximation for geological systems undergoing concurrent or consecutive types of operations. In particular, thetarget processes are fracture, displacement, and multiphase flow and transport at the system-scale.An important branch of hydromechanical models with fracture seeks to homogenize fracture within the context of acontinuum model. In terms of mechanics, continuum damage models (CDM) (e.g. continuum damage mechanics [2]and crack band theory [3]) can describe the extent of rock failure within a representative element of the continuum by adamage tensor that depends on the underlying stress field. As a continuum model, there is naturally the requirement forsufficient scale-separation allowing such homogenization. In terms of hydrodynamics, continuum flow models in frac-tured rock are also well-studied. Models such as dual-porosity and dual-permeability (e.g. [4]) or multiple interacting a r X i v : . [ c s . C E ] A ug PREPRINT - A
UGUST
10, 2020 (a) (b) (c)
Figure 1: Examples of discrete fracture representations: (a) a fitted triangular mesh; (b) an embedded Cartesian mesh;and (c) two intersecting fractures (dashed lines) represented using reduced-dimension immersed boundaries (solidlines).continua MINC (e.g. [5]) homogenize sub-continuum damage to accommodate flow and transport. Continuum mod-els abstract challenges in accommodating complex fracture intersections (e.g. [6, 7]) and can lead to computationallyefficient simulation models. On the other hand, by their very construction, such models may be inappropriate whenthere is no clear scale-separation in terms of damage or flow and the subscale heterogeneity is severe and perhaps non-isotropic (see for example [8] for flow). Additional considerations such as stress-locking when damage localizes intoa fracture, and mesh-dependency have been reported in applications of CDM ([9]). Recent work ([10, 11]) proposeda CDM model in a multi-scale context whereby it is applied as a standalone propagation criterion, and is combinedwith a discrete fracture method (explicit representation of discontinuity). Additionally hybrid approaches that combineCDM models for diffuse and dense damage with explicit discrete models for large fractures have been proposed in thecontext of coupled flow and poromechanics (e.g. [12, 13]). One challenge ahead of applying such models in dynamiccontexts is in representing the interaction between collocated CDM and explicit processes and in delineating the two.Discrete Fracture and Matrix methods approximate the continuity of mass and momentum directly on the fractureand porous media continua with transmission conditions across them. Within this class of approach, fractures maybe represented using diffusive phase-field indicators (e.g. [14, 15, 16]) or explicitly. In the latter approach, fractureaperture and geometry-dependent friction and stick-slip conditions may be modeled directly. In explicit approaches,fractures may be modeled as lower-dimensional, evolving, sharp-interfaces, or immersed boundaries with an impliedaperture field along the reduced dimension. Within this class of model, fitted methods (e.g. Figure 1a) utilize un-structured meshes to represent the matrix continuum (rock skeleton), and edges, or faces to represent fractures withimplied aperture (e.g. [17, 18, 19]). Using fitted representations, Fracture Propagation (FP) requires that the mesh beadapted dynamically. On the other hand, embedded representations (e.g. Figure 1b) offer a convenience in the choiceof matrix-mesh topology since it, in principle, need not be constrained to the fracture geometry. Subsequently, usingembedded representations, the background mesh can remain static as the fracture propagates locally. Examples ofnumerical models using embedded meshes include the extended finite element method (XFEM) (e.g. [20, 21, 22, 23])and the embedded discrete fracture method (EDFM) (e.g. [24, 25, 26]).While embedded models have been proposed to approximate coupled multiphase flow and poromechanics (e.g. [27, 1,13]), and more recently to also include hydraulic FP (e.g. [28, 29, 30, 31, 32]), several challenges remain. Firstly, thestructure of the resulting algebraic systems involve local and variable degrees-of-freedom that can hinder the efficiencyof preconditioned iterative solution methods (e.g. [1]). This aspect is universal to both monolithic and partitioneddiscretizations, and nonlinear and linear solution methods. A second challenge concerns coupling models betweenmultiphase hydromechanics in the matrix and fractures to the criteria for failure under FP, particularly involvingbranching or intersection. Finally, adaptive methods must be developed to accommodate transitions into, and out ofperiods during which fractures are to propagate. This work develops a mixed XFEM-EDFM embedded model with FP,with a focus on addressing i) solution efficiency and ii) temporal adaptivity to accommodate the onset of propagationat arbitrary times during the simulation process. In this work, we consider that fractures are open to flow as in theconditions believe to prevail during a hydraulic-fracturing operation. That is, we neglect the enforcement of contactconditions as well as associated stick and frictional slip contact models.In Section 2, the initial boundary value problem is formulated along with the fracture propagation constraints. InSections 3 and 4, the mixed discretization scheme and solution methods are developed. Section 5 presents several2
PREPRINT - A
UGUST
10, 2020numerical results that verify correctness and efficiency, accuracy using reference problems, empirical consistency, aswell as computational examples of consecutive and concurrent propagation and hydromechanics.
As illustrated in Figure 1c, we consider a spatial domain Ω ⊂ R with external boundary Γ and its associated outward-oriented unit-normal n Γ . Dirichlet and Neumann boundaries for fluid flow are denoted as Γ p and Γ Q respectively,while Γ u and Γ t are the counterparts for mechanics. The boundary segments are disjoint ( Γ p (cid:84) Γ Q = Γ u (cid:84) Γ t = ∅ ),and Γ p ∪ Γ Q = Γ u ∪ Γ t = Γ .We consider a collection of fractures Ω F = {C c ⊂ Ω , c = 1 , . . . , N F } , each of which is parameterized over a realopen interval I = ( a , a ) by a sufficiently smooth mapping, γ c : I × R + → Ω . Subsequently, a fracture’s tips are D c = { γ c ( η, t ) : η ∈ ∂I, t ≥ } , and their instantaneous velocities are ∂ γ c ∂t ( ∂I, t ) . Geometrically, the instantaneousunit-tangent at η ∈ I is ([33]), t c = 1 (cid:107) ∂ γ c ∂η (cid:107) ∂ γ c ∂η , (1)and the oriented unit-normal to one side of the fracture is subsequently defined as, n + c = 1 (cid:107) ∂ t c ∂η (cid:107) (cid:20) ∂ γ c ∂η − (cid:18) ∂ γ c ∂η · t c (cid:19) t c (cid:21) , (2)and on the other side (orientation) as, n − c = − n + c . (3)A scalar aperture field is defined on fractures ω c ( C c ) and it is assumed that ω c << . With this assumption, it isreasonable to define a matrix domain as Ω M := Ω \ Ω F . State fields defined on matrix or fracture are distinguishedby the subscripts M and F respectively.Fields defined on the matrix, Z M , are assumed to take two limiting values at each fracture location and across itsaperture; these are denoted Z + M and Z − M and may be distinct (resulting in a discontinuous field). The associated jumpacross fractures is denoted as, (cid:74) Z M (cid:75) = Z + M − Z − M . (4)The projection matrices N c := n c ⊗ n c and T c := I − N c are introduced such that the tangential gradient anddivergence operators on fractures are defined as, ∇ c i := T c ∇ i and ∇ c · i := T c : ∇ i respectively. We consider the flow of two immiscible phases; wetting and non-wetting fluid phases denoted by subscript κ ∈ { w, n } respectively. We define the following field variables, supported independently on the matrix and fracture: p κ is thephase pressure; S κ the saturation; ρ κ the mass density; and ˆ Q κ as point-sources of mass. Additionally, in the matrix,we have the Lagrangian porosity field denoted as φ ∗ . The Lagrangian porosity augments the volumetric strain fieldand thereby introduces a coupling to the deformation as is described in a subsequent section. Continuity equations areposed for each of the two phases within Ω M and Ω F . In the matrix and for each phase, we have, (cid:90) Ω M ∂ t ( ρ κ S κ,M φ ∗ ) d x + (cid:90) Γ ρ κ v κ,M · n Γ d x = (cid:90) Ω F (cid:74) ρ κ v κ,M · n c (cid:75) d x + ˆ Q κ,M , (5)and in the fracture, for c = 1 , . . . , N F ,∂ t (cid:90) C c ρ κ S κ,F d x + (cid:90) C c ∇ c · ( ρ κ v κ,F ) d x = − (cid:90) C c (cid:74) ρ κ v κ,M · n c (cid:75) d x + ˆ Q κ,F , (6)where v κ,M/F are the fluid phase velocities. Effectively, these forms of continuity imply full coupling betweenmechanics and flow. In the matrix, the Lagrangian porosity integrates strain and we assume two-phase Darcy flow,whereas in fracture, aperture is dictated by displacement and it effects the fluid velocity. We assume Poiseuille flow infracture so that the fluid occupies the entire fracture volume. The constitutive relations applied in our computational3 PREPRINT - A
UGUST
10, 2020examples are described in a later section. The problem is closed by enforcing pressure continuity across fracture andfixed fluid velocity across the outer-boundary; i.e., p κ,M = p κ,F on C c (7a) v κ,M · n Γ = 0 on Γ Q . (7b) We define the average pore pressure in the matrix or fracture using the saturation weighted average; i.e. p M/F = p w,M/F S w + p n,M/F S n . The independent variable in the geomechanical system is the displacement field u : Ω × R + → R . Under the assumption of infinitesimal deformation ( (cid:107)∇ u (cid:107) (cid:28) ), the quasi-static continuity of momentumrequires that, ∇ · ( σ ( u , p M I )) + ρ b f = 0 , on Ω M , (8)where σ is the total stress; I is the identity matrix; ρ b is the average density of the rock and fluid; and f is thebody force per unit volume. On the outer boundary, Neumann (force) and Dirichlet (displacement) conditions areconsidered, σ · n Γ = t on Γ t (9a) u = ˆ u on Γ u , (9b)while on immersed fracture boundaries, the fluid pressure p F is imposed onto the oriented surfaces of the fracture, σ · n + c = − σ · n − c = p F I · n c on C c , (10a) Adopting Biot’s single phase poroelasticity theory [34] the effective stress law is written as, σ = σ (cid:48) − αp M I , (11)where σ is the total stress tensor; σ (cid:48) is the effective stress tensor; and α ∈ (0 , is a Biot coefficient. The effectivestress σ (cid:48) is modeled using linear elasticity theory, σ (cid:48) = λ ( ∇ · u ) I + 2 G ε ( u ) , (12)where λ, G > are Lam ´ e coefficients evaluated from properties of the matrix skeleton, and the strain ε is a secondorder tensor. Under infinitesimal deformation the strain tensor is a function of displacement as, ε = ∇ sym u := 12 ( ∇ T u + ∇ u ) , (13)and the volumetric strain, (cid:15) is equal to the trace of the strain tensor, Tr ( ε ) . The Lagrangian porosity in the matrixunder linear poroelastic infinitesimal deformation is modeled after [35] as, φ ∗ = φ + α ( (cid:15) − (cid:15) ) + 1 M ( p M − p M ) , (14)where (cid:15) and p M are the reference states of volumetric strain and matrix pressure, respectively, both of which arechosen as in this work; φ is the porosity at the state of (cid:15) and p M ; and M is a Biot coefficient.The multiphase extension to the Darcy velocity in the matrix is modeled as, v κ,M = − k M · k rκ,M µ κ ( ∇ p κ,M − ρ κ g dhdz ) . (15)where k M is the second order permeability tensor; g is the gravitational force; and dhdz is the gradient of elevationwith respect to gravity. The relative permeability for the wetting and non-wetting phases in matrix is calibrated by theCorey relationship, k rn,M = k endrn ( S n − S n,r − S n,r − S w,r ) n n , (16a) k rw,M = k endrw ( S w − S w,r − S n,r − S w,r ) n w . (16b)4 PREPRINT - A
UGUST
10, 2020where Kr endw , Kr endn are the end points of the wetting and non-wetting phase relative permeability, S w,r and S n,r arethe residual saturation of the wetting and non-wetting phases n n and n w are the exponential numbers.In fracture, fluid velocity is modeled according to Poiseuille’s law, v κ,F = − ω c S κ,F µ κ ( ∇ c p κ,F − ρ κ g dhdz ) , (17)where ω c = (cid:74) u (cid:75) · n c is fracture aperture. The macroscopic concepts of Irwin’s law are adopted to formulate the model for fracture growth rate. In particular, weassume the availability of an (empirical) critical fracture toughness, K c , that is independent of fracture growth rate.Given a fracture tip a ( t ) ∈ D c , we affix a local instantaneous frame, and consider the equivalent stress intensity factor K eqI which is defined as a function of its mode I and II counterparts [21], K eqI = 12 cos( θ/ { K I (1 + cos( θ )) − K II sin( θ ) } , (18)and θ is the direction of maximum tensile stress θ = 2 arctan 14 (cid:18) K I /K II − sign ( K II ) (cid:112) ( K I /K II ) + 8 (cid:19) . (19)The relationship between θ and the ratio K I /K II is illustrated in Figure 2, which suggests that the absolute value ofFigure 2: θ vs. K I /K II θ decreases monotonically as K I /K II increases.The J integral adopted here to extract SIF is ([36]), J = lim Γ I → (cid:90) Γ I (cid:20) − σ ij ∂u i ∂x + σ ij ε ij δ j (cid:21) n j d Γ , (20)where Einstein summation rules are adopted; Domain B ρ is confined by a circle around the crack tip of radius ρ , andthen Γ I = ∂ B ρ ∪ ( B ρ ∩ C c ) ; n is the outward unit-normal to the neighborhood.In order to model K I and K II , the following relationship holds ([21]) J = 1 E (cid:48) ( K I + K II ) . (21)5 PREPRINT - A
UGUST
10, 2020where E (cid:48) is equal to E under the plane stress condition and equal to E/ (1 − ˜ ν ) under the plane strain condition. ˜ ν isthe Poisson’s ratio.In this work, the creation of initial fracture or defect is not modeled; rather, it is assumed that initially, there arepreexisting fractures represented by the mappings γ c . Under linear elastic fracture mechanics theory, equilibrium(or static) FP requires that the strain energy release rates are less than or equal to the critical [37], i.e. K eqI ≤ Kc .Subsequently, during water injection, pressure build-up within the fracture leads to an increase in K eqI . If K eqI achievesa critical threshold K c , the tip is assumed to advance, relieving the local build-up in K eqI . In a quasi-static manner,this process may continue while the local relief due to advancement and the build-up due to injection are stableand balanced. This quasi-static process terminates when the build-up due to pressure causes K eqI to drop below thethreshold criterion. Finally, a stable propagation scheme of linear elastic fracture that is regulated by the K eqI reads, (cid:26) Stable Propagation K eqI = K c Static K eqI < K c . (22) We consider the mesh T h as a subdivision of the domain Ω M into disjoint elements Ω e (quadrilaterals) in two dimen-sions. Subsequently, we define ¯Ω = ∪ Ω e ∈T h Ω e as the union set of all disjoint matrix elements. The diameter of eachelement Ω e ∈ T h is denoted by h . We denote partitioned boundaries of the domain ¯Ω as Γ et , Γ eu , Γ ep , Γ eq . Another setof meshes denoted by ˆ T h and corresponding to fracture γ c is attached to the base mesh resulting in a set of disjointelements γ ec . The fracture elements γ ec are simply chosen as the portions that are partitioned by matrix elements Ω e .Thereafter, ¯ γ c = ∪ γ ec ∈ ˆ T h γ ec is defined as the union set of all disjoint fracture elements. Generally, the fracture mesh(choice of segmentation of fracture) can be independent of the choice of matrix mesh. Accommodating this may beachieved by a preprocessing of connections and transmissibility in the context of EDFM, and by a treatment of enrich-ment in the XFEM context. While this feature may be desirable in practice, in this work we assume that the fracturesegments conform to the edges of the matrix mesh in which they are embedded. We also introduce a partition of thetime interval, I nT = ( t n , t n +1 ) with ∆ t = |I nT | , such that ¯ I T = ∪ n I nT for n ∈ , ..., N t .Flow equations are discretized using a finite-volume approximation whereas mechanical equations are approximatedusing a finite element method. The unknowns are staggered; pressure p κ,M/F and saturation S κ,M/F are cell-centeredwithin Ω e and γ ec , and displacement nodes are located at the vertices of grid Ω e . The current setting simultaneouslyensures both local mass conservation and eliminates spurious spatial instability at early times for compressible systems(see for example, [38, 39, 40, 41]). First, we follow the Bubnov–Galerkin approach that the test ( S ) and trial functional space ( S ) share the same func-tional space, S (Ω) = (cid:8) δ u | δ u ∈ H (Ω) : δ u = 0 on Γ u (cid:9) (23) S (Ω) = (cid:8) u | u ∈ H (Ω) : u = ˆ u on Γ u (cid:9) . (24)Furthermore, the test and trial functional space of admissible strain field δ ε and ε shall be denoted by E and E ,respectively. Subsequently, variational forms of the displacement δ u and strain δ ε are decomposed into standard andenhanced parts in order to account for the discontinuity, δ u = δ ¯ u + δ ˜ u , δ u ∈ S δ ε = δ ¯ ε + δ ˜ ε , δ ε = ∇ sym δ u ∈ E , (25)where the δ ¯ u and δ ¯ ε are corresponding standard parts of the displacement and strain fields while δ ˜ u and δ ˜ ε areenhanced parts of the displacement and strain fields. Here, we propose the following sets of weak forms of eq. (8)combining the Biot poroelastic theory: (cid:90) Ω δ ¯ ε : σ (cid:48) d x − (cid:90) Ω δ ¯ ε : αp M I d x − W ext ( δ ¯ u ) = 0 ∀ δ ¯ u ∈ S (cid:90) Ω δ ˜ ε : σ (cid:48) d x − (cid:90) Ω δ ˜ ε : αp M I d x − W ext ( δ ˜ u ) = 0 ∀ δ ˜ u ∈ S , (26)6 PREPRINT - A
UGUST
10, 2020where external virtual work terms can be expressed as, W ext ( δ ¯ u ) = (cid:90) Γ t δ ¯ u · t d x ∀ δ ¯ u ∈ S W ext ( δ ˜ u ) = (cid:90) γ c (cid:74) δ ˜ u (cid:75) · p F n c d x ∀ δ ˜ u ∈ S , (27)Note that p M is applied in the strain energy term δ ε : σ (cid:48) while p F is used in the external virtual work term W ext ( δ ˜ u ) .These contribute to coupling between mechanics and fluid flow. Figure 3: Illustration of enrichment nodes; Black line represents the fracture; All grid vertices are included in the nodeset I ; Squares consist of subset K while circles consist of subset L .Following the approach of [21], crack surfaces are modeled by step functions and tips by asymptotic near-tip fields.Three clusters of nodes shown in Figure 3 are defined: • I : the set of all nodes in the region ¯Ω . • L : the sub set of I that are enriched by the Heaviside function, H γ c , and L are additional nodes that are usedto capture the displacement discontinuity of the fracture body. • K : the sub set of I that are enriched by the asymptotic branch functions, F l = √ r sin θ l = 1 √ r cos θ l = 2 √ r sin θ sin θ l = 3 √ r cos θ sin θ l = 4 (28), and nodes in K are used to capture the stress singularity at the fracture tip.Therefore, the standard and enhanced terms and their variational counterparts are expressed as, ¯ u = (cid:88) i ∈ I N i ¯ u i δ ¯ u = (cid:88) i ∈ I N i δ ¯ u i ˜ u = (cid:88) i ∈ L N i ( H γ c − H i γ c )˜ a i + (cid:88) i ∈ K (cid:88) l =1 N i ( F l − F il )˜ b li δ ˜ u = (cid:88) i ∈ L N i ( H γ c − H i γ c ) δ ˜ a i + (cid:88) i ∈ K (cid:88) l =1 N i ( F l − F il ) δ ˜ b li , (29)7 PREPRINT - A
UGUST
10, 2020 ¯ ε = (cid:88) i ∈ I B ui ¯ u i δ ¯ ε = (cid:88) i ∈ I B ui δ ¯ u i ˜ ε = (cid:88) i ∈ L B ai ˜ a i + (cid:88) i ∈ K (cid:88) l =1 B b l i ˜ b li δ ˜ ε = (cid:88) i ∈ L B ai δ ˜ a i + (cid:88) i ∈ K (cid:88) l =1 B b l i δ ˜ b li , (30)where N i are finite element linear basis functions for quadrilateral elements; Enrichment functions are shifted by thefunction values evaluated at vertices i , H i γ c and F il respectively. ¯ u i is the displacement for standard nodes; ˜ a i is thedisplacement for Heaviside enriched nodes; ˜ b li is the displacement for asymptotic branch function enriched nodes, and B ui , B ai , B b l i can be derived correspondingly by substitution of eq. (29) into eq. (13).The discretized version of terms in eqs. (26) and (27) is expressed in the table 1. The subscripts i and j identifythe nodal number, D is the elastic matrix. The numerical integration in elements cut by fractures is performed usingtriangulation and coordinate transformation of the tip element ([1]).Table 1: Discretized forms of the geomechanics equation (cid:82) Ω δ ¯ ε : σ (cid:48) d x and (cid:82) Ω δ ˜ ε : σ (cid:48) d x (cid:82) Ω e ( B ri ) T D B sj d x ( r, s = u, a, b l ) (cid:82) Ω δ ¯ ε : αp M I d x and (cid:82) Ω δ ˜ ε : αp M I d x (cid:82) Ω e ( B ri ) T αp M (1 , , T d x ( r = u, a, b l ) (cid:82) γ c (cid:74) δ ˜ u (cid:75) · p F n c d x (cid:82) γ ec N i (cid:74) F l (cid:75) p F n c d x and (cid:82) γ ec N i (cid:74) H γ c (cid:75) p F n c d x (cid:82) Γ t δ ¯ u · t d x (cid:82) Γ et N i t d x A first order backward Euler fully implicit scheme is adopted for the time discretization of flow equations. Fluxesbetween different connections, i.e. matrix-matrix (M-M), matrix-fracture (M-F), fracture-fracture (F-F), are approx-imated using a two-point flux approximation (TPFA). A standard finite-volume scheme is employed as the spatialdiscretization for M-M and F-F connections. We follow the strategy proposed by [42] for M-F flux calculation. Thefully-implicit discretization for eq. (5) and eq. (6) are listed in Table 2,Table 2: The discretization of flow equationsContinuous form Discrete form for a single cell Ω e or γ ec (cid:82) Ω M ∂ t (cid:0) ρ κ S κ,M φ ∗ (cid:1) d x / ∆ t (cid:2)(cid:0) ( φ + α(cid:15) + M p M ) ρ κ S κ,M (cid:1) n +1 i − (cid:0) ( φ + α(cid:15) + M p M ) ρ κ S κ,M (cid:1) ni (cid:3) ∂ t (cid:82) γ c ( t ) ρ κ S κ,F d x L i / ∆ t (cid:2)(cid:0) ω c ρ κ S κ,M (cid:1) n +1 i − (cid:0) ω c ρ κ S κ,M (cid:1) ni (cid:3)(cid:82) Γ ρ κ v κ,M · n Γ d x (cid:80) j ∈ adjMM ( i ) T MM ( k rκ ρ κ /µ κ ) n +1( i + j ) / ( p κ,M i − p κ,M j ) n +1 (cid:82) γ c ( t ) ∇ c · ( ρ κ v κ,F ) d x (cid:80) j ∈ adjF F ( i ) T F F ( k rκ ρ κ /µ κ ) n +1( i + j ) / ( p κ,F i − p κ,F j ) n +1 (cid:82) γ c ( t ) (cid:74) ρ κ v κ,M · n c (cid:75) d x β (cid:80) j ∈ adjMF ( i ) T MF ( k rκ ρ κ /µ κ ) n +1( i + j ) / ( p κ,M i − p κ,F j ) n +1 where the superscripts n + 1 and n represent the current and previous time-steps, and L i is the length of fracturesegment i . We also assume that while fracture porosity is one, the volume of the fracture is affected by its aperture ω c . The connectivity routines, adjM M ( i ) and adjF F ( i ) return matrix or fracture cell index j adjacent to matrixor fracture cell index i or given matrix cell index i , adjM F ( i ) returns the fracture cell index j . The definitions oftransmissibility terms T MM , T F F and T MF are as proposed in [8]. First order phase potential upwinding (PPU) isemployed to complete the k rκ on the interface, while ρ κ and µ κ are interpolated from two adjacent cells. Remark 1.
If FP does not occur at a given time-step, then v = and the computational domain remains static.Therefore, ∂ t (cid:90) γ c ( t ) ρ κ S κ,F d x ≈ (cid:90) γ c ( t n +1 ) / ∆ t (cid:2)(cid:0) ρ κ S κ,F (cid:1) n +1 − (cid:0) ρ κ S κ,F (cid:1) n (cid:3) d x . (31)8 PREPRINT - A
UGUST
10, 2020
Under FP however, a new fracture segment is to be added into the system at time t n +1 , and subsequently, ∂ t (cid:90) γ c ( t ) ρ κ S κ,F d x ≈ t (cid:34)(cid:90) γ c ( t n +1 ) (cid:0) ρ κ S κ,F (cid:1) n +1 d x − (cid:90) γ c ( t n ) (cid:0) ρ κ S κ,F (cid:1) n d x (cid:35) . (32)Figure 4: Schematics of fluid interaction between fracture and matrix during FP Remark 2.
Correction for M-F fluid transfer under propagation is necessary. As illustrated in Figure 4, consider anewly established tip element that propagates from length D nc to length D n +1 c over a time-step ∆ t = t n +1 − t n . Thetotal fluid losses from the fracture to matrix is, Q MF = (cid:90) t n +1 t n (cid:90) D n +1 c D nc q MF H ( t − τ ) dxdt, (33) where H ( t − τ ) is the Heaviside function and τ is the arrival time of the fracture tip at position x . Assuming a constantpropagation velocity c and a uniform flux q MF over the time period t n +1 − t n , and utilizing a trapezoidal rule for theintegration, we have, Q MF = (cid:90) t n +1 t n (cid:90) D nc + ct D nc q n +1 MF dxdt = 12 q n +1 MF ( D n +1 c − D nc )( t n +1 − t n ) . (34) This model is applied to M-F source terms, and the β appearing in the Table 2 equals / in this model. After we assemble the Jacobian system resulted from Table 1 and Table 2, the fully coupled system solved by theNewton-Raphson method is, ∂R ∂p n,M ∂R ∂S w,M ∂R ∂p n,F ∂R ∂S w,F ∂R ∂ u ∂R ∂p n,M ∂R ∂S w,M ∂R ∂p n,F ∂R ∂S w,F ∂R ∂ u ∂R ∂p n,M ∂R ∂S w,M ∂R ∂p n,F ∂R ∂S w,F ∂R ∂ u ∂R ∂p n,M ∂R ∂S w,M ∂R ∂p n,F ∂R ∂S w,F ∂R ∂ u ∂R ∂p n,M ∂R ∂S w,M ∂R ∂p n,F ∂R ∂S w,F ∂R ∂ u δp n,M δS w,M δp n,F δS w,F δ u = − R R R R R , (35)where gradients can be derived by using the formulations in Table 1 and 2. A direct local solution for FP speed via eq. (22) is intractable. Rather, at a given time-step, the condition in eq. (22) istested for each tip within the domain. FP tip-advancement algorithms are required to indicate the onset of propagation9
PREPRINT - A
UGUST
10, 2020over a time-step, ∆ t , and to approximate a length, ∆ a with orientation, θ , of the propagation, assuming a lineartrajectory.Broadly, there are two classes of tip-advancement algorithm differing in the assumptions applied to the tip advance-ment ∆ a over time-step ∆ t , and their relation to the propagation criterion as illustrated in Figure 5. (a) Scheme A (b) Scheme B Figure 5: Comparison of two schemes for fracture propagation. Black line: old fracture segment, Blue dashed/solidline: new fracture segment; Red color stands for the parameter that needs to be determined.
Scheme A
As illustrated in Figure 5a, given a fixed target time-step, ∆ t , fracture advancement ∆ a that will satisfythe criterion eq. (22) is determined. One class of algorithm requires frequent reordering of enriched nodes and re-assembly of elastic stiffness matrices associated with enriched nodes until a suitable step-size is found. Example ofmethods in this category include, [43, 44]. The unknown, ∆ a , must satisfy eq. (22), and can be determined by solving, ∆ a = (cid:26) K eqI ( t ) − K c < | K eqI − K c | otherwise . (36) Scheme B
Alternately, as illustrated in Figure 5b, a fixed propagation length ∆ a is specified, and the correspondingtime-step size ∆ t and its associated state variables ( p κ,M/F , S κ,M/F , u ) are determined to satisfy the constraint ineq. (22), e.g. [31]. To obtain such a step, the coupled problem eq. (36) can be solved for the time-step. In thisapproach however, the modeling of multiple, simultaneously propagating fractures is challenging since the various tippropagation speeds may differ widely.An alternative combining ideas from both approaches entails selection of a target time-step size and utilizing anapproximate corresponding advancement step. Since the approximate step may not satisfy the constraint in eq. (36), asequence of subsequent solution substeps are computed such that the terminal substep satisfies the equality constraintin eq. (22). Hence, the time of this final substep can be considered as the correct time-step for the given advancementstep. In [30] and [45] a regularized form of the constraint equation is applied to identify approximate advancementincrements and at each substep in [30], the mesh is adapted and the coupled system is resolved. In that work, the fluidflow and its coupling to mechanics within the matrix are neglected. Therefore, intermittent periods of flow simulationprocesses are not modeled. This is critical for the application of such models in data assimilation of combined flowand fracture systems.In this work, we adapt this general approach to the embedded method under general two-phase flow and mechanics infracture and matrix. We address the following four technical components to achieve this:1. Accurate estimation of the SIF in embedded fractures through a numerical evaluation of J integral,2. Geometric algorithms to update fracture geometry γ c and insertion of new u into numerical discretization,3. Reliable extrapolations of state variables ( p κ , S κ and u ) following propagation in order to improve subse-quent nonlinear algebraic solution,4. Adaptive time-step selection criteria to improve computational performance at the onset of end of FP, and toimprove performance within the internal FP iteration.10 PREPRINT - A
UGUST
10, 2020
The general solution procedure is listed in Algorithm 4, and Algorithms 1-3 are sub-routines to complete the algorithm.The input I m in Algorithm 1 is a Boolean record of whether tips meet the propagation criterion. For an arbitraryelement I im ∈ I m , I im = (cid:26) TRUE < ( K eq,iI − K c ) /K c < (cid:15) SIF
FALSE ( K eq,iI − K c ) /K c < , (37)where (cid:15) SIF is a user specified tolerance. At the start of each time-step solve from line 3 and line 4 in Algorithm 4,we first update the domain based on the I m , then solve for { p κ,M/F , S κ,M/F , u } . In the post processing, from lines 5to 10, we either update I m if FP is stable or restart the solve, otherwise. Sections 4.2 - 4.4 develop the numericalapproximation to the SIF, the proposed geometric updates, and the state-variable initial guess and updates. Algorithm 1:
Domain Update Input I m Output p κ , S κ , u for i ∈ I m do if i = TRUE then Follow the routine in Section 4.3 if ∃ i ∈ I m = TRUE then Algorithm 6: Initialization() in Section 4.4
Algorithm 2:
Check StablePropagation Input K eqI Output
IS STABLE IS STABLE ← TRUE; for k ∈ D c do if ( K eqI,k − K c ) /K c > (cid:15) SIF then IS STABLE ← FALSE Break
An interaction-integral approach is adopted that superimposes the actual field with an auxiliary state within the expres-sion for energy release rate. In this work, the auxiliary is chosen as the asymptotic tip analytical solution of mode I andII fractures. Denoting the solution fields (stress, strain and displacements) with superscript (1) , and the auxiliary fieldsunder modes I and II with superscripts (2) I and (2) II respectively. For the sake of the computational convenience,the contour integral of eq. (20) can be converted into an equivalent area integral as (the derivation is shown in B), Algorithm 3:
Do Flagging Input K eqI Output I m I m ← {} ; for k ∈ D c do if < ( K eqI,k − K c ) /K c < (cid:15) SIF then I m ← I m ∪ TRUE else if K eqI,k − K c < then I m ← I m ∪ FALSE 11
PREPRINT - A
UGUST
10, 2020
Algorithm 4:
FP solver Input ∆ t, I m Output p κ,M/F , S κ,M/F , u Algorithm 1: Domain Update( I m ) Solve the fully coupled system eq. (35) Algorithm 5: K eqI ← SIF Calculation() if Algorithm 2: Check StablePropagation ( K eqI ) = TRUE then Algorithm 3: I m ← Do Flagging( K eqI ) else Apply eq. (48) to update ∆ t go to 𝑝 𝐹 𝒏 𝑐 𝒞 𝑐 (𝛼𝑝 𝑀 𝑰) ⋅ 𝒏 𝑐 𝜕ℬ 𝜌 Ω𝒟 𝑐 Figure 6: J integral graphical illustration. I I,II = (cid:90) B ρ (cid:34) σ (1) ij ∂u (2) I,IIi ∂x + σ (2) I,IIij ∂u (1) i ∂x − σ (1) ij ε (2) I,IIij δ j (cid:35) ∂q∂x j d x + (cid:90) Γ L p ∗ F q ∂u (2) I,II ∂x d Γ , (38)where Γ L = ( γ + c ∪ γ − c ) ∩ B ρ ; p ∗ F := p F − n Tc ( αp M I ) n c = p F − αp M is the net pressure applied on the fracturesurface. The graphical illustration to the line integral of eq. (38) is shown in Figure 6. Along the fracture line, while p F offers tensile forces to open C c , the effective pore pressure αp M applies compression to the surface. The weightfunction q : B ρ → [0 , is q = (cid:26) (cid:13)(cid:13) x − x ∗ tip (cid:13)(cid:13) (cid:54) ρ (cid:13)(cid:13) x − x ∗ tip (cid:13)(cid:13) > ρ . (39)The value of scalar field q is depicted in Figure 7b. In addition to the contribution from the bulk material, the secondterm in eq. (38) considers contribution of fracture pressure to the interaction integral. Note that the tensors in eq. (38)require projection onto the local reference coordinates located at the fracture tip. The transformation matrix T isdefined as, T = (cid:20) cos θ sin θ − sin θ cos θ (cid:21) , (40)and, the coordinate transform from the global frame x to the local frame x ∗ is defined as, x ∗ = T ( x − x Tip ) . (41)12 PREPRINT - A
UGUST
10, 2020Subsequently, the tensor and vector transformation of an arbitrary field Z and X between the two coordinates arecomputed as, Z ( x ∗ ) = T Z ( x ) T T (42) X ( x ∗ ) = T X ( x ) . (43) (a) (b) Figure 7: (a) Fracture tip local coordinate, (b) q value distribution, red and blue colors correspond to value 1 and 0,respectively. Algorithm 5:
SIF Calculation Output K eqI Define domain B ρ and calculate scalar field q and transformation matrix T Define two element sets, Ω B ρ and γ B ρ c that satisfy { Ω B ρ | Ω B ρ ∩ ∂ B ρ (cid:54) = ∅ } and { γ B ρ c | γ B ρ c ∩ B ρ (cid:54) = ∅ } for Ω e ∈ Ω B ρ do for each quadrature point x G with a weight W G in each elements or sub-triangles do compute σ (1) ( x ∗ ) , ∂ u (1) ∂ x ( x ∗ ) and ∂q∂ x ( x ∗ ) in the local coordinate use mode I or II fracture analytical solutions to get σ (2) ij ( x ∗ ) and ε (2) ij ( x ∗ ) integrand is calculated I ← I + (cid:0) σ (1) ij ∂u (2) i ∂x + σ (2) ij ∂u (1) i ∂x − σ (1) ij ε (2) ij δ j (cid:1) ∂q/∂x j W G for γ ec ∈ γ B ρ c do for each quadrature point x G with a weight W G do compute ∂q∂ x ( x ∗ ) in the local coordinates I ← I + (cid:0) q ∂p ∗ F ∂x u (2)2 + 2 p ∗ F ∂q∂x u (2)2 (cid:1) W G Apply eq. (18) to get K eqI Algorithm 5 develops the process applied to perform the numerical approximation. The accuracy of the computedestimate is naturally dependent not only on the quadrature approximation in Algorithm 5, but also on the accuracy ofthe underlying estimates of the state variables.
We consider piece-wise linear approximations to curved fracture paths. A given element may contain multiple frac-ture segments. To accommodate the XFEM approximation, a modification of the signed-distance function is firstintroduced, followed by proposed updates to the EDFM discretization. PREPRINT - A
UGUST
10, 2020 (a) (b)
Figure 8: A kinked fracture in an element, (a) x ∗ is closest point to x on the crack, e , e are unit norm vectors of γ c and γ c , (b) a fracture is described by H ( x ) = 0 , ψ ( x ) > and ψ ( x ) > .In order to constrain the fracture in both perpendicular and lateral directions, functions, H ( x ) , ψ ( x ) and ψ ( x ) , areused here. ψ ( x ) and ψ ( x ) are descriptions of lines in 2-D that are orthogonal to the fracture tip as illustrated inFigure 8b. Subsequent updates of these functions are necessary during FP. In a matrix cell containing a piece-wiselinear fracture C c , we rely on a signed normal distance function as follows. Given an arbitrary point x in the cell, let x ∗ ∈ C c be the closest point x on the fracture. Then the signed normal distance function φ ( x ) is defined as, φ ( x ) = sgn (( x − x ∗ ) · n c ( x ∗ )) (cid:107) x − x ∗ (cid:107) . (44)where (cid:107) (cid:107) calculates the length of a segment. As illustrated in Figure 8a, this leads in an ambiguity whenever x ∗ coincides with a vertex in the piece-wise linear construct. This is characterized by the condition that x lies within thecone formed by the two unit normal vectors. This is illustrated by the shaded region in the figure formed by n c and n c . As a convention we apply the choice n c in the signed distance function.Additionally, the Heaviside function is defined as, H γ c ( x ) := φ ( x ) > φ ( x ) = 0 − otherwise (45) Two possible schemes are proposed here as illustrated in Figure 9. In the scheme 1, the two fracture segments ( γ c , γ c )are treated as a whole, and hence, there is only one set of degrees of freedom ( p κ and S κ ) assigned to the grid. Theintegration region is confined by the dashed lines in the rectangular region, and the average distance between a fractureis computed as, (cid:104) d (cid:105) = (cid:82) V ∗ (cid:107) x − x ∗ (cid:107) dvV ∗ , (46)where V ∗ is the integration region. A constrained triangulation is applied to partition the domain, and numericalintegration is performed on each sub triangle. This scheme facilitates computation of the average distance in caseswith multiple fracture segments in an element while limiting the total number of degree of freedoms introduced in thesystem. However, once a new fracture segment is introduced during FP, the properties of the old segment γ c has to bemapped to the combined segment ( γ c + γ c ).A second approach is to treat γ c and γ c separately with two distinct sets of dofs as illustrated in Figure 9b. Theaverage normal distances are computed using eq. (46) to each segment. All computational examples in this workapply the second approach. As illustrated in Figure 10a, once the stress intensity factor exceeds the FP criterion, fracture tips are to be advanced bya fixed advancement length ∆ a along a direction θ . Pressure and saturation initialization within the newly introduced14 PREPRINT - A
UGUST
10, 2020 (a) Scheme 1 (b) Scheme 2
Figure 9: Two approaches to model multiple fracture segments in an element; Red and blue solid lines representfracture segments, dashed lines confine the integration region for (cid:104) d (cid:105) calculations. (a) (b) Figure 10: (a) Fracture propagation algorithm schematics; (b) Initialization of newly added fracture segmentsfracture segments is necessary for the solution of the resulting nonlinear algebraic system. As illustrated in Figure10b, the pressure and saturation initialization can be performed without any further geometric processing since theconnection-list data structure provides constant time look-up of the face and its two connected cells; subsequently theupdate is, p κ,F j := p κ,F adjFF ( j ) . (47)Note that, while this update serves as an initial guess after propagation, the initial fluid volume prior to the additionof the segment is assumed to be zero. On the mechanics side, a drained condition computation is performed to as-sign displacements of the newly introduced enriched nodes. This step is critical since the fracture aperture is directlyassociated to the displacement of enriched nodes, and a poor initial guess results in severely degraded nonlinear con-vergence behavior. By fixing p κ,M/F and S κ,M/F , the discretized equations of mechanics in table 1 result in a linearproblem. The algorithm for initialization of new fracture segments p κ,M/F , S κ,M/F and u is listed in Algorithm6. Algorithm 6:
Initialization Apply eq. (47) to assign initial value to the newly added p κ,F and S κ,F Solve eqs. (26) and (27) for u by fixing p κ,M/F and S κ,M/F .Once the initialization is complete, the nonlinear system at time-step t n +1 is solved by means of an Inexact Newtonmethod. After primary unknowns are obtained, a time-step controller will decide the next taken time-step. However,if the propagation condition is violated, a linear interpolation is applied for the calculation of the time-step cut, ∆ t := α SIF K c (1 + (cid:15) SIF ) − K eq,oldI K eqI − K eq,oldI ∆ t, (48)15 PREPRINT - A
UGUST
10, 2020where α SIF is a damping factor, K eq,oldI and K eqI are previous and current equivalent SIF respectively, and K c areperturbed by (cid:15) SIF to ensure the positive coefficient. Figure 11 illustrates the chord iterative process of eq. (48).Starting from an initial guess t n +11 , the corresponding SIF falling on the blue curve equals K eq,n +1 I, that is over thecritical, K c . An interpolation is made to find the next time-step t n +12 − t n = α SIF K c (1+ (cid:15) SIF ) − K eq,nI K eq,n +1 I, − K eq,nI ( t n +11 − t n ) . Theprocedure is performed iteratively until the criterion is satisfied.Figure 11: Graphical illustration of the time-step selection strategy. Blue curve: the relationship between time andSIF; Red lines: interpolation lines.The proposed time-step strategy is not only suitable for FP simulation, but also for production processes. For instance,a hydraulic fracturing process followed by long-range fluid flow will lead to a geometric increase in time-step sizeunder the same framework, if ∆ t max is set to the largest time-scale. This occurs since at the transition, the SIF willdecrease. On the other hand, under a transition from fluid production to onset of FP, the proposed time-step eq. (48)will geometrically be rescaled thereby avoiding excessive redundant solves. The proposed methods are implemented within an in-house software framework [8, 25, 46]. Validation results illus-trating accuracy are presented for benchmark models of elastic mechanics with no fluid flow. This is followed byseveral numerical examples of propagation with coupled multiphase flow and poromechanics.
The viscosity- and toughness-dominated propagation regimes of the KGD fracture model are considered. The dimen-sionless toughness variable K m characterizing these regimes is defined as, K m = K (cid:48) (2 E (cid:48) µ (cid:48) q ) / , (49)where K (cid:48) = 8 K C / √ π , µ (cid:48) = 12 µ , E (cid:48) = E/ (1 − v ) is the equivalent Young’s modulus in the plain strain condition,and q is the flow rate into one wing of the fracture. If K m < , the flow lies in the viscosity dominated regime andthe solution is called the M-vertex solution when K m = 0 . Under this scenario, rock is very brittle ( K C → ) so thatenergy dissipation is dominant in viscous flow. On the other hand, if K m > , the process is toughness-dominated,and in the limit that K m → ∞ , the solution is called a K-vertex solution. In that case, energy is mostly utilized tobreak the rock and influences from either small aperture or highly viscous fluid could be neglected. Model instancesof the two regimes are created using the parameters listed in Table 3.16 PREPRINT - A
UGUST
10, 2020Table 3: Rock and fluid properties for KGD model verification E (GPa) ˜ ν K c ( MPa √ m ) (cid:15) SIF µ ( Pa · s ) q ( m / s ) K m toughness dominated 8.3 0.25 0.5 0.001 e − × isapplied unless noted otherwise. Incompressible fluid is injected into the center point of the fracture at a constant flowrate of q (refer to Table 3). In order to model the unbounded domain assumed in the analytical reference solutions,the mid-points on the left and right edges of the domain are restricted from displacement in the y-direction, while thoseon the top and bottom of the domain are restricted in x-direction. The remainder of the boundary is traction free. Theinitial fracture length ( m ) is also small enough relative to the simulation domain size. Analytical solutions for eitherpropagation regime are reviewed in A.With reference to the analytical solutions, the relative error in propagation length L is, ε L = | L ( t ) − L ( t ) | L ( t ) , (50)where L is analytical solutions. Similarly, the relative error in aperture ω c at the wellbore is defined as, ε ω c = | ω c ( t ) − ω c ( t ) | ω c ( t ) , (51)where ω c is the analytical solution.Asymptotic mesh refinement results for the K- and M-vertex cases are presented in Figures 13a and 13b respectively.In these computations, finely resolved temporal and tip-length steps are fixed ( ∆ t = 1 e − s , (cid:15) SIF = 2 e − forthe K vertex case and (cid:15) SIF = 1 e − for the M-vertex case) whereas a refinement path along h is studied. The J-integral approximations apply radii of h , and the tip element enrichment scheme is based on the topological geometrywhere only the element containing the tip is enriched by branch functions. The computed error compares the stablenumerical solution for tip length in each simulation with the corresponding analytical solution. Convergence rates inboth scenarios are approximately 0.5 order. Note that the propagation scheme relies not only on state variables, butalso on the numerical approximation of the equivalent SIF using the converged state variables. Similar convergencerates applied to elastic problems have been reported previously [47]. To improve the accuracy of SIF estimate, thefixed-area tip enrichment scheme proposed in [48] could be applied here. Applying the toughness dominated properties listed in Table 3, Figures 14a and 14b show the time-series of fracturelengths computed (blue) using advancement step-lengths of ∆ a = 0 . and . respectively compared to the evolutionof the analytical solution (red). Note that in these figures, the circular data markers represent solution-steps at whichthe propagation is stable as per Algorithm 2, whereas the stair-casing connectors indicate intermediate sub-steps of thealgorithms. Clearly, selecting a larger advancement step-size for a given propagation velocity leads to a proportional17 PREPRINT - A
UGUST
10, 2020 (a) K vertex (b) M vertex
Figure 13: mesh refinement study of the K and M vertex solutionsincrease in the time-separation between successive stable solutions. As is reflected by Figures 14e and 14f, this comesat no extra cost in terms of computation in this case. Moreover, Figures 14c and 14d present the computed lengthand aperture relative-error time-series respectively using either step-size. These figures also compare the applicationof two different time-step algorithms. The Na¨ıve scheme increases or cuts the time-step size by a constant factor,whereas the adaptive scheme applies the estimate in Equation 48. In Figure 14e, we compare the cumulative time-stepsolves over the course of the simulation for two different schemes. In this case, the proposed estimate can reducecomputational cost by 70% percent using (cid:15)
SIF = 0 . . The computational performance and relative error obtained in the viscosity-dominated scenario are similar to thoseof the toughness-dominated case. They are summarized by Figures 15a through 15d. The pressure profile along thefracture at a fixed time is extracted and plotted in Figure 15e, where the axes are normalized by the value of inletpressure and current length (3.95m). The computed pressure profile matches the analytical solution very well withinthe interior segment of the fracture, away from the regions near the tips. The analytical solution dictates that pressureexhibits a singular behavior near the tip, and obeys an asymptotic function O ((1 − xL ) − / ) in the M vertex solution([49]). In the numerical model, the pressure behavior is not explicitly constrained by this asymptotic function andlinear interpolation is used instead. A model problem is considered where an initial horizontal fracture of length m is centered within a square domainof sides m in length. The matrix permeability is e − m and porosity is 0.1. Two phase fluid flow in the porousmedia and fracture are considered; fluid density and viscosity are assumed to be equal (1000 kg/m ) with isothermalcompressibility of . e Pa − and . Pa · s, respectively. Corey-Brook relative permeability is applied to flowwithin the matrix, and linear models to flow within fracture. Initially, the pressure is p M = 5 MPa; water saturationis S w = 0 . ; the Biot coefficient is α = 0 . ; fracture toughness is K c = 4 e Pa √ m; and the tolerance is set at (cid:15) SIF = . m / s.Two variations of the model problem are considered: mixed-mode type 1 (MMT-1) as depicted in Figure 16a, andmixed-mode type 2 (MMT-2) as depicted in Figure 16b. Under MMT-1 or 2, the shear stress will re-direct the fracturepath along a direction that is perpendicular to the local maximum tensile stress. “Mixed mode” refers to load typesincluding both opening and in-plane shear ([50]). the difference between MMT-1 and MMT-2 is the magnitude offorces applied to the left and right boundaries.In both MMT-1 and MMT-2, the shear force applied to the boundary alters the direction of the principal stresses andhence the FP path does not follow a straight line. Due to the different stress magnitudes acting on the left and rightboundaries in the cases MMT-1 and MMT-2, the curvature of the path manifests differently.18 PREPRINT - A
UGUST
10, 2020 (a) (b)(c) (d)(e) (f)
Figure 14: Simulation results of toughness dominated FP. (a) ∆ a = 0 . : fracture length vs. time; (b) ∆ a = 1 . :fracture length vs. time; (c) relative error of the aperture at the wellbore. Dashed lines indicate intermediate stateswhile circles are accepted solutions that satisfy the FP criterion. (d) relative error of the propagation length; (e) ∆ a = 1 . : cumulative time-step solves vs. time; (f) ∆ a = 1 . : cumulative nonlinear iterations vs. time;19 PREPRINT - A
UGUST
10, 2020 (a) (b)(c) (d)(e)
Figure 15: Simulation results of viscosity dominated FP. (a) fracture length vs. time; (b) relative error of the fracturelength; (c) aperture at the wellbore vs. time; (d) relative error of the aperture at the wellbore; (e) pressure profile at FPlength = 3.95; (f) Number of Nonlinear iterations at first 1 sec period time. Dashed lines indicate intermediate stateswhile circles are accepted solutions that satisfy the FP criterion.20
PREPRINT - A
UGUST
10, 2020 (a) Mixed Mode type 1 (MMT-1) (b) Mixed mode type 2 (MMT-2)
Figure 16: Initial fracture configuration and boundary conditions, y displacement is fixed at the mid point of left andright edges and x displacement is fixed at the mid point of top and bottom edges. The asymptotic convergence of the predicted fracture path is studied with reference to ∆ a and h . In particular, defininga relative error as, ε c = (cid:13)(cid:13) Y ( x ) − Y ( x ) (cid:13)(cid:13) (cid:107) Y ( x ) (cid:107) , (52)where the mapping function Y : R + → R + takes a vector of x coordinates of reference points and yields a vector of y coordinates of the same group of points on the propagation path. The reference state, Y is selected as the predictionusing the finest mesh in the refinement paths listed in Table 4, where the minimum ∆ a chosen on the finest mesh size h is larger or equal to the radius r of B ρ to improve the accuracy of J-integral approximation.Table 4: Sensitivity analysis: parameters of mesh sizes ( h ) and propagation length step ( ∆ a ) in porous mediaSensitivity Tests Meshes ∆ a MMT-1 MMT-2 ∆ a ×
173 4.4 2.2187 ×
173 2.2 1.1187 ×
173 1.1 0.55187 ×
173 0.55 0.275 h ×
37 2.2 2.287 ×
83 2.2 2.2137 ×
127 2.2 2.2187 ×
173 2.2 2.2As we can see from Figures 17a, 17b, 17c and 17d, all scenarios exhibit asymptotic decay in the relative error withrespect to ∆ a and h . Note that asymptotic decay in relative error is not sensitive to the higher curvature of the path inMMT-2 relative to MMT-1 in these cases.Figures 18a and 18c present the paths and fracture length evolutions obtained for cases MMT-1 and MMT-2 respec-tively. Note that in the MMT-2 test, the initial fracture surface is not perpendicular to the minimum horizontal stressthat is 0.05 times of the initial reservoir pressure p i . This is intentionally designed to investigate the sensitivity to ahighly curved fracture path. As expected, the fracture path exhibits a larger curvature than that in MMT-1. Despitedistinct FP path curvature present in Figure 19a, differences in paths reduce with ∆ a . Therefore, in order to track theFP path closely, a finer mesh h and smaller ∆ a are required. The FP speed reflected in Figure 19e is not affected by ∆ a significantly except that the speed associated with ∆ a = 2 . becomes slower at late times. Mesh refinement h does not have an significant impact on the propagation path and its speed shown from Figure 19c and 19e.Figure 20 presents snapshots of the pore pressure and saturation fields under MMT-1. The pressure fields show ahigh pressure band in the matrix cut through by the hydraulic fracture. The band diffuses along the normal direction.Lower fluid pressure happens behind the fracture tip region due to high tensile stresses and the resultant dilation in the21 PREPRINT - A
UGUST
10, 2020 (a) (b)(c) (d)
Figure 17: (a) ∆ a consistency test in MMT-1; (b) h consistency test in mixed MMT-1; (c) ∆ a sensitivity test in mixedMMT-2; (d) h sensitivity test in MMT-2.corresponding matrix rock. Overall, in these cases, only a small mount of water leaks into the rock matrix as a resultof the low matrix permeability as well as the scale contrast between the fracture and matrix. Using the model problem illustrated in Figure 21a, the influences of the matrix permeability and Biot parameters onpredicted simulations are investigated. Initial reservoir pressure in all cases is 5 MPa; water saturation is 0.6; waterinjection rate is . m / s; critical SIF K c is set . e Pa √ m; the tolerance parameter is (cid:15) SIF is ; and ∆ a is . m.By comparing the resultant fracture lengths obtained at a particular time while assuming k M ∈ { e − , e − , e − m } and α = 0 . , it is apparent that lower matrix permeability leads to a higher propagation rate. Indeed, Figure22 shows that a two-order of magnitude increase in matrix permeability can lead to a reduction in propagationspeed. This is associated with the increased fluid leak-off from fracture to matrix. Figure 21b shows the FP direction isslightly altered by the fluid leak-off. K II is a strong function of shear deformation which, in this case, is predominantlydriven by the lateral forces applied to the top and bottom boundaries. The increased local matrix pressure induced bythe fluid leak-off leads to a larger K I /K II ratio. Subsequently, this results in reduced curvature in the propagationpath observed for the high permeability field. The results agree with the observation in [51].Figure 21c shows results obtained while fixing permeability at k M is fixed at e − m , and selecting α as either or . . Figure 21c shows less curvature in the high α scenario. Increasing the coefficient α in this case causes thepore pressure to counter the compressive forces applied at the boundary and subsequently leads to an increase in theratio K I /K II . This in turn leads to a reduction in the curvature of the resulting FP path. The value and direction(white lines) of the maximum principal stress of two different scenarios are compared in Figure 23. We can observedifferences in the direction of the stress field near the fracture tip region that causes the deviation in the FP path. In the22 PREPRINT - A
UGUST
10, 2020
10 15 20 25 30 35 40
X axis(m) Y ax i s ( m ) a = 4.40 a = 2.20 a = 1.10 a = 0.55 (a) (b)
10 15 20 25 30 35 40
X axis(m) Y ax i s ( m )
37 3787 83137 127187 173 (c) (d)(e) (f)
Figure 18: Left column: fracture patterns at t = 50 s; right column: fracture length vs. time. (a), (b) and (e) ∆ a sensitivity tests on MMT-1; (c), (d) and (f) h sensitivity tests on MMT-1. (b) and (d) are zoomed counterparts of (a)and (c), respectively.rest of the domain, as compared to the decoupled case ( α = 0 ), pore pressure also alters the direction of the maximumprincipal stress along the fracture body. 23 PREPRINT - A
UGUST
10, 2020 (a) (b)
10 15 20 25 30 35 40
X axis(m) Y ax i s ( m )
37 3787 83137 127187 173 (c)
14 15 16 17 18 19 20
X axis(m) Y ax i s ( m )
37 3787 83137 127187 173 (d)(e)
Figure 19: Left column: fracture patterns at t = 50 s; right column: fracture length vs. time. (a) and (b) ∆ a sensitivitytests on MMT-2; (c), (d) and (e) h sensitivity tests on MMT-2. (d) is the zoomed counterparts of (c). The test problem is illustrated in Figure 24. A × mesh is applied, and water is injected into the center ofeither fracture at the same rate of . m /s . The critical SIF is set at e Pa √ m with tolerance and ∆ a is chosenas . m.Figure 25 shows snapshots of the fracture path, pore pressure, saturation and displacement distributions obtained forvarious times during the simulation. The two fractures initially propagate symmetrically until the stress field in the24 PREPRINT - A
UGUST
10, 2020 (a) 37 ×
37 (b) 87 ×
83 (c) 137 ×
123 (d) 187 × ×
37 (f) 87 ×
83 (g) 137 ×
123 (h) 187 × Figure 20: pore pressure (top, unit: Pa) and saturation (bottom) distribution at the end of simulation 50 seconds underMMT-1 FP. (a) domain schematics(b) permeability effects (c) Biot coefficient effects
Figure 21: Poroelastic influences on FP in porous mediavicinity of the fracture ends responds to their interaction ([47]). Subsequently, the localization is evident as the twotips overtake each other and arrest as they propagate towards the interior of the domain. This localization effect issimilar to results reported using the phase field and XFEM methods (e.g. [47, 48]).25
PREPRINT - A
UGUST
10, 2020
Time(s) F r ac t u r e l e ng t h ( m ) k M = 5e-13 = 0.8k M = 5e-14 = 0.8k M = 5e-15 = 0.8 Figure 22: fracture length evolution under different k M (a) α = 0 (b) α = 0 . Figure 23: Comparison of maximum principal stress (unit: Pa) and its direction (white lines) for different α A key objective is to co-simulate periods of fluid flow and deformation with no propagation, and the onset and dynam-ics of propagation of one or more fracture. Timely engineering questions relying on such ability pertain to the designof hydraulic fracture treatments and implementation strategies in domains undergoing complex sequences of fracture,infill drilling, and production. We consider a hypothetical scenario to model such a system.The initial fracture pattern is showed in Figure 26. Stimulation-alternating-production is assumed to occur over fouroperational stages: (1) fractures A and B are hydraulically fractured by injection of water into their centers at the rateschedule illustrated in Figure 27a, (2) fractures C and D are hydraulically fractured by injection of water into theircenters at the rate schedule illustrated in Figure 27b, (3) fluid withdrawal (production) is subsequently conducted for30 days by enforcing a constant pressure of . e Pa at the centers of all fractures, and (3) production is halted, andwater is injected into the center of fracture E using the injection rate schedule in Figure 27c. Water injection schedulesfor these three consecutive operations are illustrated in Figure 27. A mesh of × is applied; K c is e Pa √ m; thetolerance (cid:15) SIF is . ; the propagation length step ∆ a is . m; the initial time-step is set at . secs; the maximumstep is set at . e secs (5 days); α SIF is chosen as 1.0; and the Biot coefficient α is set to 0.4.26 PREPRINT - A
UGUST
10, 2020Figure 24: Initial fracture configuration and boundary conditions, where x and y displacement is fixed at four cornersof the rectangular.Figure 25: Snapshots (by column) of the FP path, pore pressure (unit: Pa), water saturation, and the displacement(unit: m) at t = 3 . , . , and . secs.Figure 28, shows snapshots of fracture paths, pore pressure, water saturation, and displacement fields obtained at threetimes during the simulation (the ends of stages 1, 2, and 4). The stress shadow effects and localization instability areevident in the paths obtained following stages 1 and 2. Moreover, Fracture E is observed to propagate in a symmetric,straight line since the stress field following symmetric production is also symmetric. Finally, while a large pressuregradient is observed in the vicinity of the fractures, saturation change due to water leak-off into the matrix is ratherlimited.Figure 29 presents the evolution of the nonlinear iteration count and time-step sizes over the duration of the simula-tion. During the hydraulic fracturing stages, time-step sizes remained below 10 secs, and no more than 10 nonlineariterations were required for convergence. During the transition from hydraulic fracturing to production (HF-Prod)the time-step size was exponentially adapted to reach the specified maximum of . e secs (5 days). At the secondtransition from production to hydraulic fracturing (Prod-HF), the time-step size was rapidly adapted to the stimulationtime scales. This is an indication of the efficacy of eq. (48) to prevent wasted effort due to unconverged time-stepattempts. Figure 30 compares the evolution in cumulative nonlinear iterations computed using a naive time-step selec-tion strategy with that of eq. (48).There are further implications of utilizing the proposed time-step selection process. In particular, Algorithm 4 showsthat in between accepted solutions satisfying both the propagation constraint with equality ( K eqI = K c ) and conser-vation, there may be a number of sub-steps that violate the constraint ( K eqI ≤ K c ) and over which the effective SIFis to build-up. Another possibility is for states that overshoot the constraint ( K eqI ≥ K c ), and lead to a time-step27 PREPRINT - A
UGUST
10, 2020Figure 26: Initial fracture configuration and boundary conditions, x and y displacement is fixed at four corners of therectangular; fracture coordinates are (97, 160) to (103, 160), (97, 140) to (103, 140), (97, 100) to (103, 100), (97, 60)to (103, 60) , (97, 40) to (103, 40) from the top to the bottom. (a) Injection stage 1 (b) Injection stage 2(c) Injection stage 3 Figure 27: Water rate injection schedules for three stimulation periods. Stages 1 and 2 are separated from the stage 3by a 30 day period of fluid withdrawal.cut. Figures 31a and 31a present the numbers of intermediate time-steps (blue) and failed time-steps (red) for twelveaccepted propagation solutions using SIF tolerances of e − and e − respectively. Comparing the na¨ıve and theproposed strategies, both the total number of intermediate time-steps, and failed steps can be reduced significantly. Themore strict tolerance criterion leads to a larger number of intermediate and steps, and show a more drastic differencebetween the approaches. 28 PREPRINT - A
UGUST
10, 2020Figure 28: Snapshots of (by column) the FP path, pore pressure (unit: Pa), water saturation, and displacement fields(unit: m) at the conclusion of stages (by row) 1,2, and 4.Figure 29: Nonlinear iterations and time-stepsFigure 30: Cumulative time-step solves over the course of simulation29
PREPRINT - A
UGUST
10, 2020 (a) (cid:15)
SIF = 0 . (b) (cid:15) SIF = 0 . Figure 31: Number of time-step solves to reach at the FP criterion for the first 12 FP steps.Finally, to emphasize the implications to systems involving intermittent fracturing and withdrawal operations, thesimulation is repeated after slightly modifying the third operational (production) stage. In particular, during thisperiod, production is only allowed to occur through fractures A and B instead of A, B, C, and D. Figure 32 presentsthe pore pressure snapshots at the end of the production period and the final fracture patterns obtained for the twoscenarios. The result demonstrates the utility of the integrated model to enable co-design of such operations.30
PREPRINT - A
UGUST
10, 2020 (a) (b)(c) (d)
Figure 32: Comparison of pore pressure fields (top row, unit: Pa) and fracture paths (bottom row) obtained using thesymmetric production stage (left column) and the asymmetric plan (right column).
Figure 33: Initial domain configuration.In this case, simultaneous and interspersed periods of production and stimulation are studied to demonstrate the sig-nificant influence of poroelasticity on fracture propagation. The scale of the model is m × m with × grid blocks. The preexisting fractures are shown in fig. 33 and are located at wells A and B . The matrix permeabilityis − m ; Porosity φ is . ; Biot poroelastic coefficient α is . ; Fluid viscosity of oil and water is identical, . cp; Compressibility of oil and water is . e − Pa − ; Young’s modulus E is . GPa and Poisson’s ratio ν is . ; K c is31 PREPRINT - A
UGUST
10, 2020 MPa √ m, ∆ a = 3 m, (cid:15) SIF = 0 . . The stress at the left and right boundaries is 28 MPa while at the top and bottomboundaries it is 26 MPa. The well schedules are as follows:Wells A and B are set to produce fluids for 5 years, and the well pressures are fixed at 5 MPa;Well C is not operated for the first three years. Then, an injection process is triggered, and hydraulic fracturingtakes place over 13 minutes. Following this period of injection at constant rate, Well C is turned to producefluids for the remainder of time in the simulation. The injection rate is controlled at 0.001 m /s whileproduction is controlled by setting the well pressure to 5 MPa, which is lower that the ambient pressure in themodel.In order to avoid the closure of fractures due to pressure depletion, proppant is assumed to be present uniformly and ismodeled by application of a uniform (elastic) force to the fracture surfaces as an inner boundary condition. This forceis assumed to be elastic through an inversely linear relationship with aperture. (a) 3 years(b) 5 years Figure 34: Pore pressure distribution of the poroelastic case; Short black lines are the minimum principal stress; Longblack lines are the propagation path.fig. 34a shows the pressure field at three years, prior to the 13 minute fracturing process that is to occur in Well C. Inthe figure, the quivered black lines depict the minimum principal stress field at this time. Moreover, the solid blackline depicts the fracture path that will ultimately be formed at the end of the 13 minute injection period. Due to theproduction and resulting pressure depletion drive by Wells A and B , the direction of the minimum principal stressfield is rotated to 90 degrees in the infill region, and becomes perpendicular to the minimum horizontal stress that isapplied at the boundary. During fracture propagation, the path will follow the direction perpendicular to the minimumprincipal stress and hence, extend vertically through the infill region. fig. 34b presents the pressure field at the end ofthe simulation period (five years).The above simulation is repeated while ignoring the pore pressure effects on effective stress by setting the Biot coeffi-cient to α = 0 and retaining all other conditions. The results corresponding to conditions in fig. 34 but for this case areshown in fig. 35. Without poroelastic effects, the stress reorientation in the infill region is not observed, and the mini-mum principal stress field remains vertical and aligned with the minimum horizontal stress. The fracture propagationpath is horizontal and hits the preexisting fractures at Wells A and B . This effect is significant to the hydrodynamicsas well: fig. 36 shows the cumulative evolution of fluid volumes that are produced by Well C in either case, and wecan observe a difference of 67% over the five period between the elastic and poroelastic models.32 PREPRINT - A
UGUST
10, 2020 (a) 3 years(b) 5 years
Figure 35: Pore pressure distribution of the elastic case; Short black lines are the minimum principal stress; Longblack lines are the propagation path.Figure 36: Comparison of the cumulative oil production of poroelastic and elastic scenarios.The summary of the nonlinear solver performance for the poroelastic case is shown in fig. 37. The average nonlineariterations per time step over the course of the simulation are 4.36 while the stimulation period takes around 46.7% ofthe total nonlinear iterations. 33
PREPRINT - A
UGUST
10, 2020Figure 37: Nonlinear solver performance of the poroelastic scenario.
A coupled multiphase flow and mechanical model is proposed allowing the co-simulation of fluid-driven fracturing aswell as injection and depletion processes in porous media using embedded meshes. The proposed mixed discrete frac-ture and matrix EDFM-XFEM discretization is augmented with an adaptive time-step controller, extended J-integralcomputation for poromechanics, and state initializations for propagated segments. The computational results show:1. The extension of the J integral estimation to poromechanics provides sufficient accuracy and/or leads toconsistent solutions.2. The proposed state-initialization strategy in newly propagated segments leads to improved nonlinear solverperformance.3. The proposed time-step controller is effective in automatically adapting to transitions, leading to computa-tional efficiency.4. Using the proposed unified model to simulate intermittent production and fracture requires computationaleffort that is on par with that of tying separate propagation and hydromechanical models (e.g. [52, 53]). Thisis facilitated by the time-step adaptivity to factor propagation onset.Extension of the proposed methodology to three-dimensions relies primarily on the availability of efficient and stablethree-dimensional computational geometry infrastructure. In particular, this would be required to extract of geo-metric information such as fracture piece-wise planar intersections with the background mesh, fracture leading-edgelocations, and fracture intersection. Nevertheless, the functional forms remain the same; for example the J-integralestimates presented here for SIF computations need to be performed in three-dimensions.
This material is based upon work supported by the U.S. Department of Energy under Award Number DE-FE-0031777.The authors also acknowledge partial funding from the members of the TU Future Reservoir Simulation Systems &Technology (FuRSST) Industry-University Consortium.
References [1] Guotong Ren, Jiamin Jiang, and Rami M Younis. A model for coupled geomechanics and multiphase flow infractured porous media using embedded meshes.
Advances in Water Resources , 122:113–130, 2018.[2] Lazar M Kachanov. Rupture time under creep conditions.
International journal of fracture , 97(1-4):11–18, 1999.[3] Zdenˇek P Baˇzant and Byung H Oh. Crack band theory for fracture of concrete.
Mat´eriaux et construction ,16(3):155–177, 1983. 34
PREPRINT - A
UGUST
10, 2020[4] JE Warren, P Jj Root, et al. The behavior of naturally fractured reservoirs.
Society of Petroleum EngineersJournal , 3(03):245–255, 1963.[5] Yu-Shu Wu, Karsten Pruess, et al. A multiple-porosity method for simulation of naturally fractured petroleumreservoirs.
SPE Reservoir Engineering , 3(01):327–336, 1988.[6] Cunbao Li, Viet T Chau, Heping Xie, and Zdenˇek P Baˇzant. Recent advances in mechanics of fracking and newresults on 2d simulation of crack branching in anisotropic gas or oil shale.
Acta Mechanica , 229(2):975–992,2018.[7] Kumchol Yun, Tae-Jong Kim, Paek San Jang, Zhenqing Wang, and Sakaya Ronald. An improved crack track-ing algorithm with self-correction ability of the crack path and its application in a continuum damage model.
International Journal for Numerical Methods in Engineering , 117(2):249–269, 2019.[8] Jiamin Jiang, Rami M Younis, et al. Hybrid coupled discrete-fracture/matrix and multicontinuum models forunconventional-reservoir simulation.
SPE Journal , 21(03):1–009, 2016.[9] Milan Jir´asek and Thomas Zimmermann. Analysis of rotating crack model.
Journal of engineering mechanics ,124(8):842–851, 1998.[10] Simon-Nicolas Roth, Pierre L´eger, and Azzeddine Soula¨ımani. A combined xfem–damage mechanics approachfor concrete crack propagation.
Computer Methods in Applied Mechanics and Engineering , 283:923–955, 2015.[11] Elena Tamayo-Mas, Jordi Feliu-Fab`a, Montserrat Casado-Antolin, and Antonio Rodr´ıguez-Ferran. Acontinuous-discontinuous model for crack branching.
International Journal for Numerical Methods in Engi-neering , 120(1):86–104, 2019.[12] Guotong Ren, Jiamin Jiang, Rami M Younis, et al. Fully-coupled xfem-edfm hybrid model for geomechanicsand flow in fractured reservoirs. In
SPE Reservoir Simulation Conference . Society of Petroleum Engineers, 2017.[13] Xia Yan, Zhaoqin Huang, Jun Yao, Zhao Zhang, Piyang Liu, Yang Li, and Dongyan Fan. Numerical simula-tion of hydro-mechanical coupling in fractured vuggy porous media using the equivalent continuum model andembedded discrete fracture model.
Advances in Water Resources , 126:137–154, 2019.[14] Timo Heister, Mary F Wheeler, and Thomas Wick. A primal-dual active set method and predictor-correctormesh adaptivity for computing fracture propagation using a phase-field approach.
Computer Methods in AppliedMechanics and Engineering , 290:466–495, 2015.[15] Sanghyun Lee, Mary F Wheeler, and Thomas Wick. Pressure and fluid-driven fracture propagation in porousmedia using an adaptive finite element phase field model.
Computer Methods in Applied Mechanics and Engi-neering , 305:111–132, 2016.[16] Chukwudi Chukwudozie, Blaise Bourdin, and Keita Yoshioka. A variational phase-field model for hydraulicfracturing in porous media.
Computer Methods in Applied Mechanics and Engineering , 2019.[17] Benoit Carrier and Sylvie Granet. Numerical modeling of hydraulic fracture problem in permeable mediumusing cohesive zone model.
Engineering fracture mechanics , 79:312–328, 2012.[18] Randolph R Settgast, Pengcheng Fu, Stuart DC Walsh, Joshua A White, Chandrasekhar Annavarapu, and Fred-erick J Ryerson. A fully coupled method for massively parallel simulation of hydraulically driven fractures in3-dimensions.
International Journal for Numerical and Analytical Methods in Geomechanics , 41(5):627–653,2017.[19] TT Garipov, M Karimi-Fard, and HA Tchelepi. Discrete fracture model for coupled flow and geomechanics.
Computational Geosciences , 20(1):149–160, 2016.[20] Jens M Melenk and Ivo Babuˇska. The partition of unity finite element method: basic theory and applications.
Computer methods in applied mechanics and engineering , 139(1-4):289–314, 1996.[21] Nicolas Mo¨es, John Dolbow, and Ted Belytschko. A finite element method for crack growth without remeshing.
International journal for numerical methods in engineering , 46(1):131–150, 1999.[22] Bernd Flemisch, Inga Berre, Wietse Boon, Alessio Fumagalli, Nicolas Schwenck, Anna Scotti, Ivar Stefans-son, and Alexandru Tatomir. Benchmarks for single-phase flow in fractured porous media.
Advances in WaterResources , 111:239–258, 2018.[23] Lars H Odsæter, Trond Kvamsdal, and Mats G Larson. A simple embedded discrete fracture–matrix model for acoupled flow and transport problem in porous media.
Computer Methods in Applied Mechanics and Engineering ,343:572–601, 2019.[24] Liyong Li, Seong H Lee, et al. Efficient field-scale simulation of black oil in a naturally fractured reser-voir through discrete fracture networks and homogenized media.
SPE Reservoir Evaluation & Engineering ,11(04):750–758, 2008. 35
PREPRINT - A
UGUST
10, 2020[25] Jiamin Jiang and Rami M Younis. An improved projection-based embedded discrete fracture model (pedfm) formultiphase flow in fractured reservoirs.
Advances in water resources , 109:267–289, 2017.[26] Matei T¸ ene, Sebastian BM Bosma, Mohammed Saad Al Kobaisi, and Hadi Hajibeygi. Projection-based embed-ded discrete fracture model (pedfm).
Advances in Water Resources , 105:205–216, 2017.[27] Rajdeep Deb and Patrick Jenny. Modeling of shear failure in fractured reservoirs with a porous matrix.
Compu-tational Geosciences , 21(5-6):1119–1134, 2017.[28] T Mohammadnejad and AR Khoei. An extended finite element method for hydraulic fracture propagation indeformable porous media with the cohesive crack model.
Finite Elements in Analysis and Design , 73:77–95,2013.[29] Saeed Salimzadeh and Nasser Khalili. A three-phase xfem model for hydraulic fracturing with cohesive crackpropagation.
Computers and Geotechnics , 69:82–92, 2015.[30] P Gupta and CA Duarte. Coupled hydromechanical-fracture simulations of nonplanar three-dimensional hy-draulic fracture propagation.
International Journal for Numerical and Analytical Methods in Geomechanics ,42(1):143–180, 2018.[31] Elizaveta Gordeliy and Anthony Peirce. Coupling schemes for modeling hydraulic fracture propagation usingthe xfem.
Computer Methods in Applied Mechanics and Engineering , 253:305–322, 2013.[32] Elizaveta Gordeliy and Anthony Peirce. Implicit level set schemes for modeling hydraulic fractures using thexfem.
Computer Methods in Applied Mechanics and Engineering , 266:125–143, 2013.[33] Heinrich Walter Guggenheimer.
Applicable Geometry: Global and Local Convexity . RE Krieger PublishingCompany, 1977.[34] MA Biot and DG Willis. The elastic coefficients of the theory of consolidation.
J. appl. Mech , 24:594–601,1957.[35] Rick H Dean, Xiuli Gai, Charles M Stone, Susan E Minkoff, et al. A comparison of techniques for couplingporous flow and geomechanics.
Spe Journal , 11(01):132–140, 2006.[36] Matthew C Walters, Glaucio H Paulino, and Robert H Dodds Jr. Interaction integral procedures for 3-d curvedcracks including surface tractions.
Engineering Fracture Mechanics , 72(11):1635–1663, 2005.[37] Zdenˇek P Baˇzant, Marco Salviato, Viet T Chau, Hari Viswanathan, and Aleksander Zubelewicz. Why frackingworks.
Journal of Applied Mechanics , 81(10):101010, 2014.[38] PA Vermeer and A Verruijt. An accuracy condition for consolidation by finite elements.
International Journalfor numerical and analytical methods in geomechanics , 5(1):1–14, 1981.[39] M´arcio A Murad and Abimael FD Loula. On stability and convergence of finite element approximations of biot’sconsolidation problem.
International Journal for Numerical Methods in Engineering , 37(4):645–667, 1994.[40] Birendra Jha and Ruben Juanes. A locally conservative finite element framework for the simulation of coupledflow and reservoir geomechanics.
Acta Geotechnica , 2(3):139–153, 2007.[41] Jihoon Kim, Hamdi A Tchelepi, and Ruben Juanes. Stability and convergence of sequential methods for coupledflow and geomechanics: Drained and undrained splits.
Computer Methods in Applied Mechanics and Engineer-ing , 200(23-24):2094–2116, 2011.[42] Seong H Lee, MF Lough, and CL Jensen. Hierarchical modeling of flow in naturally fractured formations withmultiple length scales.
Water resources research , 37(3):443–455, 2001.[43] Michael J Hunsweck, Yongxing Shen, and Adri´an J Lew. A finite element approach to the simulation of hydraulicfractures with lag.
International Journal for Numerical and Analytical Methods in Geomechanics , 37(9):993–1015, 2013.[44] Qinglei Zeng, Zhanli Liu, Tao Wang, Yue Gao, and Zhuo Zhuang. Fully coupled simulation of multiple hydraulicfractures to propagate simultaneously from a perforated horizontal wellbore.
Computational Mechanics , 61(1-2):137–155, 2018.[45] Nathan Shauer and Carlos Armando Duarte. Improved algorithms for generalized finite element simulations ofthree-dimensional hydraulic fracture propagation.
International Journal for Numerical and Analytical Methodsin Geomechanics , 43(18):2707–2742, 2019.[46] G Ren and RM Younis. A coupled xfem-edfm numerical model for hydraulic fracture propagation. In
ECMORXVI-16th European Conference on the Mathematics of Oil Recovery , 2018.36
PREPRINT - A
UGUST
10, 2020[47] Bertrand Paul, Maxime Faivre, P Massin, Richard Giot, Daniele Colombo, F Golfier, and Alexandre Martin. 3dcoupled hm–xfem modeling with cohesive zone model and applications to non planar hydraulic fracture propa-gation and multiple hydraulic fractures interference.
Computer Methods in Applied Mechanics and Engineering ,342:321–353, 2018.[48] HanYi Wang. Numerical investigation of fracture spacing and sequencing effects on multiple hydraulic fractureinterference and coalescence in brittle and ductile reservoir rocks.
Engineering Fracture Mechanics , 157:107–124, 2016.[49] Jose I. Adachi.
Fluid-driven fracture in permeable rock . PhD thesis, University of Minnesota, 2001.[50] TL Anderson.
Fracture Mechanics: Fundamentals and Applications . Taylor & Francis, 1995.[51] Tomoya Usui, Saeed Salimzadeh, Adriana Paluszny, and Robert W Zimmerman. Effect of poroelasticity onhydraulic fracture interactions. In
Poromechanics VI , pages 2008–2015. 2017.[52] Xuyang Guo, Kan Wu, Cheng An, Jizhou Tang, John Killough, et al. Numerical investigation of effects ofsubsequent parent-well injection on interwell fracturing interference using reservoir-geomechanics-fracturingmodeling.
SPE Journal , 2019.[53] Ali Rezaei, Birol Dindoruk, and Mohamed Y Soliman. On parameters affecting the propagation of hydraulicfractures from infill wells.
Journal of Petroleum Science and Engineering , 182:106255, 2019.[54] CF Shih, B Moran, and T Nakamura. Energy release rate along a three-dimensional crack front in a thermallystressed body.
International Journal of fracture , 30(2):79–102, 1986.
A Analytical solutions for viscosity and toughness dominated solutions
The K vertex solution is well known problem of a uniformly pressurized crack also called as Griffith’s crack. E (cid:48) is theplain strain stress tensor, µ (cid:48) is the µ . Q is the flow rate into two fracture wings. The analytical expressions for thecrack length and fracture aperture at the wellbore with time are, l = (cid:0) E (cid:48) Qt √ πK C (cid:1) / (53) ω c (0 , t ) = 4 K C √ πE (cid:48) l . (54)The M vertex solution assumes a zero toughness fracture and all of the energy overcome the flow in the channel. Weadopt the 10th order solution from [49], l = 0 . (cid:0) E (cid:48) Q t µ (cid:48) (cid:1) / (55) ω c (0 , t ) = 1 . (cid:0) µ (cid:48) E (cid:48) t (cid:1) / (cid:0) E (cid:48) Q t µ (cid:48) (cid:1) / (56)The self similar solutions for aperture and pressure profiles for both M and K vertex solutions are also derived in theliterature and not repeated here. Please refer to the original paper ([49]) for more details. B J integral derivation
The derivation for the eq. (38) is shown in the section. The analytical expression of the J integral is J = lim Γ I → (cid:90) Γ I (cid:20) − σ ij ∂u i ∂x + σ ij ε ij δ j (cid:21) n j d Γ , (57)where the effective stress is applied here and the superscript is neglected for the conciseness. Following the divergencetheorem, a contour integral is casted into a volume and a line integral in 2 dimension ([54]), J = (cid:90) B ρ (cid:18) σ ij ∂u i ∂x − σ ij ε ij δ j (cid:19) ∂q∂x j + ∂ x j (cid:18) σ ij ∂u i ∂x − σ ij ε ij δ j (cid:19) qd x − (cid:90) Γ L t j q ∂u j ∂x d Γ , (58)37 PREPRINT - A
UGUST