The energy release rate for non-penetrating crack in poroelastic body by fluid-driven fracture

A new class of constrained variational problems, which describe fluid-driven cracks (that are pressurized fractures created by pumping fracturing fluids), is considered within the nonlinear theory of coupled poroelastic models stated in the incremental form. The two-phase medium is constituted by solid particles and fluid-saturated pores; it contains a crack subjected to non-penetration condition between the opposite crack faces. The inequality-constrained optimization is expressed as a saddle-point problem with respect to the unknown solid phase displacement, pore pressure, and contact force. Applying the Lagrange multiplier approach and the Delfour–Zolésio theorem, the shape derivative for the corresponding Lagrangian function is derived using rigorous asymptotic methods. The resulting formula describes the energy release rate under irreversible crack perturbations, which is useful for application of the Griffith criterion of quasi-static fracture.


Introduction to poroelastic modeling
In the paper, we proceed the development of constrained optimization theory for a new class of variational models arising in poroelasticity and motivated by hydrofracking. A two-phase poroelastic body consisting of solid phase and pores saturated with a Newtonian fluid is considered. We suggest that the body contains a fluid-driven crack (called fractures) since formed by the pressure of a pumped fluid. For physical consistency, the crack is subjected to a non-penetration inequality between opposite faces (the fracture walls). This description allows a compressive pressure at which the crack might close. Here, it would be worthwhile to comment that mutual contact of adjacent crack faces admits the phenomenon of mechanically closed, but hydraulically open cracks (which could arise, e.g., through the presence of debris in an otherwise fluid-conducting crack).
The poroelastic model is described by governing equations stated in incremental form with respect to unknown solid phase displacement, pore pressure, and contact force. The system is endowed with the fluid pressure, which is prescribed inhomogeneous and different on the fracture walls. In the multi-scale formulation, the pressurized fracture equations are coupled with governing equations for the fluid pressure to a single model. Typically, fluid flow in the fracture is governed by the Reynolds lubrication equation, which assumes a local cubic law (see Baykin and Golovin [1]). Modeling of the fluid pressure using a linear diffraction equation was suggested in Mikelic´et al. [2]. In our work, we account for the channelized fluid flow as a prescribed boundary condition. In its turn, the boundary data can be achieved by flow modeling as well as directly from geomechanical data.
The nonlinear theory of solids with non-penetrating cracks and their quasi-static propagation was developed in the variational framework by Khludnev and co-authors [3,4]. For dynamic modeling of cracks, we cite the monograph of Bratov et al. [5]. The non-penetration approach was continued for frictional contact phenomena at the crack in Itou et al. [6] and the limiting small strain in the proceeding works [7,8]. We cite the study [9] for Timoshenko plates with cracks, and the study [10] addressing optimal control problems. Also anti-cracks, rigid, and soft inclusions were incorporated in the theory (see Khludnev et al. [11]). For suitable numerical methods, see Hintermu¨ller et al. [12]. Recently, in Kovtunenko [13], we derived non-penetration conditions at the fluid-driven crack in two-phase poroelastic medium.
Alternatively to the sharp-interface approach, in a brittle zone, the crack surface can be approximated by a phase-field function as described in Mikelic´et al. [14] which may be beneficial for numerical reasoning. Then, the crack and its propagation are determined based on the energy minimization approach to brittle and quasi-brittle fracture (see Kovtunenko [15]). The readers may find helpful the discrete perturbation of global potentials due to crack extension in the vein of variational eigen-erosion methods from Schmidt et al. [16].
The concept of soil and poromechanics was established well by Biot and Terzaghi [17,18] and further developed by Barenblatt et al. [19] and Meirmanov [20] and others. We cite Fellner and Kovtunenko [21] and Kovtunenko and Zubkova [22] for homogenization of a two-phase medium consisted of solid phase and pores, and Sazhenkov et al. [23] for the related multi-scale analysis. In our modeling, we follow the hydraulic fracturing formulation given by Golovin and Baykin [24] and Skopintsev et al. [25] with co-authors as presented next.
For a linear elastic solid phase, the second-order symmetric tensors of linearized strain e and Cauchy stress s are connected by Hooke's law with the help of the fourth-order symmetric tensor of elastic coefficients A, which assumed to be elliptic, and subjected to a prestress t 0 . The prestress admits mechanical stresses of geological layers in reservoir in their natural state as well by fracking (see the influence of the prestress on the failure zone development in Valov et al. [26]). Accounting for the pore pressure p, the effective stress is introduced as where a 2 (0, 1 is the Biot coefficient, and I is the identity tensor. Omitting inertia terms in equations of motion and keeping the minus sign, the quasi-static equilibrium equation reads After substitution into equation (3) of equations (1) and (2) and the symmetric gradient of the displacement vector u where T stands for the transposition; it implies the elliptic equation with respect to unknown u The fluid content in pores is constituted by where S . 0 is the storativity, and tr e implies dilatation according to equation (4). In the mass balance the flow velocity vector q is assumed given by the Darcy flow where k = k r =h r is determined by the permeability k r . 0 and the effective viscosity h r . 0. Inserting equations (6) and (8) into equation (7) results in the parabolic equation with respect to ∂p=∂t and ∂u=∂t From the mathematical point of view, the fully coupled poroelastic equations (5) and (9) present a degenerate elliptic-parabolic system; thus, standard existence theorems are not applicable here. After differentiation of the elliptic equation (5) with respect to time, the system turns into a pure parabolic problem. Its solvability was established by applying the theory of implicit evolution equations (see Showalter [27]). However, the parabolic problem does not conform to the unilateral conditions. On the other side, the governing equations formally coincide with thermoelastic equations when replacing the pore pressure p for temperature. From the literature on thermoelasticity, existence results utilizing the pseudomonotonone theory were known (see Khludnev and Kovtunenko [3], section 3.3), however, restricted to small coupling coefficients a. Avoiding these restrictive assumptions, in Kovtunenko [13] we proved the well-posedness based on Rothe's semi-discretization in time of parabolic equation (9), which reduces it to the elliptic equation with respect to unknown p at fixed time t . d . 0 for given p tÀd := p(t À d) and u tÀd := u(t À d), then passing the time step d to zero.
In the current contribution, we investigate shape differentiability of the poroelastic problem with nonpenetrating crack under irreversible shape perturbations. For this task, we consider the problem in the incremental form (5) and (10), and endow it with a saddle-point formulation. Based on the Lagrange multiplier approach, we apply the formalism of directional differentiability for Lagrangians (see Delfour and Zole´sio [28]) and use rigorous asymptotic methods (see Gonza´lez et al. [29]) to derive a shape derivative for the underlying Lagrangian function implying a free energy. The resulting formula describes the energy release rate under irreversible crack perturbations, which is useful for the Griffith criterion of quasi-static crack evolution (see Charlotte et al. [30]). Other shape derivatives were derived in a series of works for non-penetrating cracks and inclusions in linear elastic bodies by Khludnev and his colleagues [3,4] in Khludnev and Shcherbakov [31] within the Euler-Bernoulli beam theory, in Rudoy and Shcherbakov [32] for Kirchhoff-Love plates, in Lazarev [33], Lazarev and Rudoy [34] for Timoshenko plates, and so on.
The structure of the paper is the following one. In section 2, we state the poroelastic problem with non-penetrating crack in the incremental form. In section 3, variational principle is given for a Lagrangian function, and well-posedness of the corresponding saddle-point problem is established. Also, we formulate the Griffith fracture criterion for the crack quasi-static evolution. In section 4, the shape differentiability of the Lagrangian is proven using the asymptotic methods of analysis based on regular perturbations, thus providing us with a semi-analytic formula for the energy release rate. Special cases of the formula are discussed in section 5; its relation to well-known path-independent integrals (called Cherepanov-Rice, or J-integrals) is presented.

Formulation of poroelastic problem with non-penetrating crack
We start with a description of geometry as illustrated for 2D in Figure 1.
Let O be a reference domain in the Euclidean space of points We assume the Lipschitz continuous boundary ∂O with outward normal vector n = (n 1 , . . . , n d ) T , and the disjoint union ∂O = G D _ S G N with G D 6 ¼ [. Let an oriented manifold of co-dimension one S split O into two sub-domains O 6 with Lipschitz continuous boundaries ∂O 6 such that For a time parameter t 2 (0, T , T . 0, we look for a crack evolution along the interface which is assumed to be C 1, 1 -smooth and irreversible such that We distinguish the crack faces G 6 t & S 6 and chose the normal vector n t at G t outward to O À , thus inward to O + . Physically, G t represents fractures, whereas the complement implies a reservoir. For every fixed t 2 (0, T and x 2 O t in the time-dependent domain from equations (11) to (14), the poroelastic medium is described by the pore pressure p(t, x) and the solid displacement u = (u 1 , . . . , u d ) T (t, x). The latter is involved in the strain e(u) = fe ij (u)g d i, j = 1 according to equation (4) with the entries u i, j := ∂u i =∂x j of the gradient ru = fu i, j g d i, j = 1 . The stress s = fs ij g d i, j = 1 (t, x) and the effective stress t = ft ij g d i, j = 1 (t, x) are introduced in equations (1) and (2), respectively, and the prestress is given by the symmetric tensor For the reason of analysis, we do not consider here the so-called 2.5D models when physical strain and stress are 3 × 3-tensors defined over a 2D-domain O. The system is governed by equations (5) and (10), where componentwisely (div t) i := P d j = 1 t ij, j for i = 1, . . . , d, and the trace tr e(u) = div u :=

The elastic coefficients in equation (5)
for all u, v 2 H 1 (O t ) d , which is uniformly elliptic and bounded: there exists 0\a ł a such that holds for all t 2 ½0, T according to Korn and Poincare´inequalities if u = 0 on G D . In equation (10), the transport coefficient k 2 W 1, ' (O) is assumed uniformly positive and bounded and the time-delayed data in O tÀd for t 2 (d, T ) are given by functions such that the irreversibility of crack evolutions (13) provides the inclusion We decompose the displacement u and the stress tn := ( P d j = 1 t 1j n j , . . . , P d j = 1 t dj n j ) T at the boundary into its normal components implying the vector-vector and matrix-vector multiplications, and tangential components as follows Let the following time-dependent data are prescribed in the reservoir for all t 2 (0, T ) With its help, we state mixed inhomogeneous boundary conditions on the outer boundary The assumed regularity of the data will be used further for asymptotic expansions in sections 4 and 5.
Across the crack G t , functions defined in O t allow discontinuity by the mean of jump We suggest no tangential effective stress at the crack faces and continuity of the fluid pressure over the fracture walls The fluid pressure p 6 re prescribed in equation (19) is different on G 6 t , coincide at the crack-tip, respectively, crack-front in 3D, and can be determined from the lubrication equations in fractures (see Golovin and Baykin [24]).
Assuming at G 6 t the standard boundary condition in the normal direction would lead to interpenetration between the opposite crack faces under compressive stress. For the physical consistency, non-penetration at the crack is suggested see Figure 1. The inequality constraint (24) leads to complementary conditions Conditions (25) imply that equality (23) holds at those points where the crack is open, i.e., n T t ½½u . 0. Otherwise, the closed crack n T t ½½u = 0 in equation (25) has the compressive stress n T t tn t + p re ł 0.

Variational principle for the crack problem
In the domain with crack defined in equation (14), we have the following generalized Green's formula (see Khludnev and Kovtunenko [3], section 1.4) for the elasticity operator Here, the boundary stresses tn on G N and tn t on G 6 t are distributions defined in a generalized sense by duality mappings hÁ , Ái G N and hÁ , Ái G 6 t , which turn into usual integrals for functions. For the stationary transport operator, Green's formula holds for all functions (26) the equilibrium equation (5), using the Neumann condition tn = g on G N from equation (20), and tn t = (n T t tn t )n t at the crack due to the zero tangential stress in equation (21), we obtain

Inserting into equation
By the virtue of ½½n T t tn t + p re = 0 in equation (25), adding and subtracting p re follows the variational equation with respect to t ð in the Lions-Magenes space of functions, which continuation by zero in S belongs to H 1=2 (S). Its counter-part in the duality hÁ , Ái G t is determined in the adjoint space of linear continuous functionals Then, the complementarity conditions (24) and (25) take the weak form and newly introduced variable l in equation (29) implies the contact force. Inserting the transport equation (10) into Green's formula (27) and using equation (18) results in the variational equation with respect to p ð for all test functions q 2 H 1 0 (O t ). Gathering the weak variational formulation (28)-(31) and recalling t = Ae(u) + t 0 À apI, for the triple (u, p, l), we define a Lagrange function L : accounting for the identity apI: e(u) = atr e(u)p and multiplying the quadratic terms by 1/2. With its help existence of a weak solution to the problem is established in the next.
Theorem 1 (solution existence). There exists a triple (u t , p t À p re , l t ) 2 K(O t ) in the feasible set solving uniquely the saddle-point problem for all test functions (v, q, m) 2 K(O t ). Then, it solves the poroelastic problem with non-penetrating crack stated in the weak form of equations (28)- (31), and vice versa.
Proof. With respect to the primal variable u7 !L (u, p, l; G t ), the Lagrangian in equation (32) builds a quadratic bilinear form, which is bounded and positive definite due to the estimates (16) of A. With respect to the dual variable, the quadratic bilinear form p7 !L (u, p, l; G t ) is bounded and negative definite because of estimates (17) of k. The mapping l7 !L (u, p, l; G t ) is linear. Therefore, the unique saddle-point in equation (33) exists by the virtue of minimax theorems. Based on the optimality condition for equation (33), we calculate the Gateaux derivative of the Lagrangian and get the variational equation (28) for u = u t , the stress t t := Ae(u t ) + t 0 À ap t I, and the contact force n T t t t n t + p re = l t according to equation (29). Conversely, from equation (28), it follows by convexity the minimum in equation (33) L u t , p t , l t ; G t ð Þł L v, p t , l t ; G t ð Þ : Similarly, computing the limit results in equation (31) for p = p t and u = u t . The converse assertion that equation (31) implies the maximum is true by the concavity of p7 !L (u, p, l; G t ). The maximum in equation (33) with respect to m taken at q = p t implies the dual complementarity conditions which are equivalent to equation (30) for l = l t and u = u t . The proof is complete. For a perturbation parameter s 2 (0, T À t), we consider an irreversible crack perturbation G t + s satisfying equation (13) (see illustration in Figure 1) and the perturbed domain with crack according to equation (14). Let space points y = (y 1 , . . . , y d ) T be related to the perturbed geometry O t + s . The perturbed Lagrangian L : The perturbed saddle-point problem (33) reads for all test functions (v, q, m) 2 K(O t + s ) in the perturbed feasible set According to Theorem 1, there exists the unique solution (u t + s , p t + s À p re j t + s , l t + s ) 2 K(O t + s ) to equation (37). It is also the solution to the perturbed poroelastic problem with non-penetrating crack from equations (28)- (31) ð for all test functions v 2 H 1 (O t + s ) d with v = 0 on G D , where t t + s = Ae(u t + s ) + t 0 À ap t + s I and l t + s = n T t + s t t + s n t + s + p re ; the perturbed complementarity conditions and the perturbed stationary transport equation for all test functions q 2 H 1 0 (O t + s ). With the help of reduced Lagrangian L in equation (32) calculated on the saddle-point from equation (33), and its perturbation in equation (36) calculated on the saddle-point from equation (37), we define a directional derivative (called the shape derivative) as the one-sided limit Physically, equation (41) implies the energy release rate by extension of the crack. For a constant surface energy density g . 0, let us denote the increase in surface energy due to creation of the new crack by Based on equations (41) and (42), Griffith's fracture criterion can be stated as the following condition Together with irreversibility (13), the inequalities in equation (43) imply the two cases: if ∂=∂tL (u t , p t , l t ; G t ) + G t \0, then G t + s = G t and crack does not grow; if ∂=∂tL (u t , p t , l t ; G t ) + G t ø 0, then jG t + s j . jG t j and crack will begin to grow.
For the reason of fracture criterion (43), the main aim of our further consideration will be to provide a formula for calculating the shape derivative ∂L =∂t in equation (41) (also the limit G t in equation (42)).

Energy release rate by fluid-driven fracture
The crack perturbation can be carried out either in explicit or implicit form. In the explicit case, given a kinematic flow associates a coordinate transformation y = f s (x) and its inverse We suppose that it builds a diffeomorphism of the cracked domains in equations (14) and (35) From equation (44), a time-dependent kinematic velocity is defined as Lj t + s := ½df s =ds8f À1 s . In the implicit case, let a vector of kinematic velocity be given such that preserving the outer boundary and irreversible cracks in equation (13). This determines the flow in equation (45) by means of solutions to the Cauchy problem for non-autonomous and nonlinear ordinary differential equation (ODE) system and to the initial problem for a linear transport equation where the gradient r y f À1 s = (∂(f À1 s ) i =∂y j ) d i, j = 1 , and Lj t + s (y) = L(t + s, y). We assume the both equations (44) and (46) hold.
The following Traits 1-4 are needed to prove the shape differentiability of the Lagrangian.

Trait 1 (bijection of feasible sets). The function composition with f s forms is a bijective map between the feasible sets
Indeed, equation (50) follows straightforwardly from the diffeomorphism in equation (45). Trait 1 allows us to transform one-to-one the perturbed Lagrangian L from equation (36) to the reference geometry by settingL for all (v, q, m) 2 K(O t + s ). Applying equation (51) to the perturbed saddle-point problem (37), we havẽ L s,ũ t + s ,q,m; G t ð Þ łL s,ũ t + s ,p t + s ,l t + s ; G t À Á łL s,ṽ,p t + s ,l t + s ; for all test functions (ṽ,q,m) 2 K(O t ), and (ũ t + s ,p t + s À p re ,l t + s ) 2 K(O t ) is the unique solution to equation (52). Thus, the next trait holds.
Trait 2 (existence of saddle point). The set of saddle-points (ũ t + s ,p t + s À p re ,l t + s ) in equation (52) is a singleton for every s 2 ½0, T À t. We write the s-dependent LagrangianL defined in equation (51) in the explicit form following from equation (36) Here, we have used the chain rule r y v = rf ÀT s rṽ, the notation of d-by-d symmetric tensor such that E(I,ṽ) = e(ṽ) according to equation (4), the Jacobian determinant and the fact that f s is the identity transformation at G N .
Trait 3 (asymptotic expansion). The asymptotic expansion ofL from equation (53) as s ! 0 + holds The partial derivative ∂L =∂s in equation (56) is a continuous function of the first argument given by the explicit representation for t 2 ½0, T À t) where the tangential divergence div G t L := divL À n T t Ln t at G t .
Proof. As s ! 0, the following asymptotic expansion of the terms entering equations (53)-(55) takes place (see, e.g., Sokolowski and Zolesio [35], Chapter 2) Inserting representations (58) into the LagrangianL given by equation (53), we derive its expansions (56) in the first argument. Since Lj t + t is a continuous function of the argument t + t, then the partial derivative t7 !∂L =∂s(t, Á ) in equation (57) is continuous and implies Lj t + t = L at t = 0. This finishes the proof.
The last trait is rather involved and proven in Appendix 1.
Trait 4 (strong convergence). There exists a subsequence of saddle-points (ũ t + s k ,p t + s k À p re ,l t + s k ) in equation (52) converging for s k ! 0 as k ! ' to the saddle-point (u t , p t À p re , l t ) in equation (33) Traits 1-4 satisfy all assumptions in Delfour and Zole´sio [28] (Chapter 10, Theorem 5.1), thus provide the following theorem (see the detailed proof in Gonza´lez et al. [29]). Theorem 2 (shape differentiability of Lagrangian). The shape derivative from equation (41) exists expressed by where (u t , p t , l t ) is the solution to the poroelastic problem with non-penetrating crack (28)- (31), and formula for the partial derivative ∂=∂sL is given in equation (57).
In the following, we specify our main result stated in Theorem 2 with respect to the so-called J-integrals well-known in brittle fracture for linear elastic bodies with cracks.

Representation of the energy release rate as J-integral
For the kinematic velocity L from equations (46) Typically, O surrounds crack-tip, crack-front, kinks, and other singular points, where singular solutions are locally admissible. Inside O, we assume the velocity constant, e.g., equal to one, such that We denote for short O t := OnG t . Based on properties (61) and (62), in the following, we will integrate by parts the expression in equation (57).
where the integrals are the effective stress t t = Ae(u t ) + t 0 À ap t I. In equation (68) t t is a tangential vector at ∂G t positive oriented to n t in 2D, and b t = t t × n t is a binomial vector within the moving frame at ∂G t in 3D.
Proof. We rearrange the terms in formula (57) on the solution (u t , p t , l t ) in the sum The terms not including rL are gathered in Since the assumptions of regularity (61) and (62), the complementarity conditions in equation (25) hold pointwise at G t nO, therefore In I 2 , . . . , I 6 , we integrating by parts with respect to rL. Using t t = Ae(u t ) + t 0 À ap t I, we calculate where the last term L T re(u t ): t t will be shortened when adding and div t t = 0 due to the equilibrium equation (3), then which due to equation (10) constitutes the zeroth term L T rp t S p t À p tÀd ð Þ+ atr e u t À u tÀd ð ÞÀddiv krp t ð Þ ½ = 0, together with I 2 and For I 6 , differentiating along the crack gets where in 2D a tangential vector t t at ∂G t is positive oriented to n t , and in 3D a binomial vector b t = t t × n t builds the moving frame at ∂G t . The integrals from I 3 and I 6 over G t nO, where the stress is determined pointwisely, can be combined using the boundary conditions (21), (23), definition of l t in equation (29), and the following calculation ð by decomposing the vectors into normal and tangential components and using L T n t = 0 at the crack G t . The pointwise product in equation (70) is zero due to the complementarity conditions (25) after application of the tangential differentiation r G t := r À n t (n T t r) to n T t ½½u t = 0. Finally, collecting the integral terms I 1 -I 6 in equation (69) together with equation (70) provides the expressions (63)-(68).
We present several concluding remarks as follows: With the help of tangential gradient and L T n t = 0 at G t , formula (66) can be expressed equivalently Applying the coordinate transformation f s , using expansion of v s in equation (58), and differentiating along the crack, we calculate the limit in equation (42) G t = 2g lim ( It is worth noting that J ∂OnG t from equation (66) is related to J-integrals known for linear elastic bodies. We clarify the relation in the following corollary.
Corollary 1 (non-penetrating crack in linear elastic body under fluid-driven fracture). If the factors S = a = k = 0, then the Lagrangian L in equation (32) implies the strain energy L u, p, l; G t ð Þ= The strain energy release rate is given according to equation (57) by Let the elasticity coefficients A and prestress t 0 be constant. Then, under assumptions (61) and (62), the shape derivative is expressed equivalently by the integrals from equations (63)-(68) If p re is constant and the crack G t is plane in O, then rp re = rn t = 0 at G t \ O such that J G t \O = 0, and the strain energy release rate in equation (72) implies the path-independent sum We finish the paper by emphasizing that formulas of the shape derivative obtained in the paper are of practical use to predict when fluid-driven fractures will begin to grow.