Analysis of progressive fracture in fluid‐saturated porous medium using splines

Abstract Powell‐Sabin B‐splines are employed to model progressive fracturing in a fluid‐saturated porous medium. These splines are defined on triangles and are 𝒞1‐continuous throughout the domain, including the crack tips, so that crack initiation can be evaluated directly at the tip. On one hand, the method captures stresses and fluid fluxes more accurately than when using standard Lagrange elements, enabling a direct assessment of the fracture criterion at the crack tip and ensuring local mass conservation. On the other hand, the method avoids limitations for discrete crack analysis which adhere to isogeometric analysis. A crack is introduced directly in the physical domain. Due to the use of triangles, remeshing and crack path tracking are straightforward. During remeshing transfer of state vectors (displacement, fluid pressure) is required from the old to the new mesh. The transfer is done using a new approach which exploits a least‐squares fit with the energy balance and conservation of mass as constraints. The versatility and accuracy to simulate free crack propagation are assessed for mode‐I and mixed‐mode fracture problems.

(A) (B) F I G U R E 1 (A) A solid body with an internal discontinuity Γ c . Γ c is an interface boundary with positive and negative sides, Γ + c and Γ − c , respectively. Boundary Γ u is prescribed with a displacementū; Γ t with a prescribed tractiont; Γ p with a prescribed fluid pressure p; Γ q with a prescribed inflowq; (B) pressure around the internal discontinuity Γ c infinitesimal strains and linear elastic material behavior are assumed. There is no mass transfer or chemical interaction between the solid and the fluid. In most poroelastic systems, the deformation of the solid occur fast in comparison with the pressure change of the interstitial fluid. Thus, the solid deformation could be assumed as a quasi-static process. Not considering gravity, convective or inertia terms, the strong form of the hydro-static momentum equations and the boundary conditions reads: 32 in which n refers to the normal vector at the boundaries.ū andt denote prescribed displacements and tractions, respectively. u 0 andu 0 represent initial displacements and velocities separately. ⋅ □ denotes the time derivative: and denotes the total stress defined as where is the Biot coefficient, p represents the apparent fluid pressure and I denotes the unit tensor. s is the stress inside the solid material, which relates to the infinitesimal strain via with D the fourth-order elastic stiffness tensor. The tractions t c act on the crack interface Γ c , which is defined as t c = t (⟦u⟧) − p d n, (5) in which p d is the pressure inside the crack Γ c , see Figure 1B. t (⟦u⟧) are tractions due to the influence of the crack interface.
In this contribution, a cohesive-zone model is used to model tractions across the crack. The model relates the tractions t (⟦u⟧) on Γ c to the displacement jump across it. In practice, to obtain t (⟦u⟧), we first compute the cohesive tractions t d in the local coordinate system (s, n), see Figure 1A.
with ⟦v⟧ being the displacement jump across Γ c , given in the local coordinate system (s, n). The tractions t (⟦u⟧) and the displacement jump [[u]] in the global coordinate system (x 1 , x 2 ) are then obtained via a standard transformation: with R the rotation matrix.
In the current study, an exponential traction-opening law is used to obtain the tractions on Γ c : where t u is the fracture strength,  c denotes the fracture energy, d int represents the initial crack shear stiffness (when In the case of unloading (f < 0), the tractions are obtained from a secant relation. To avoid interpenetration, a penalty stiffness k p is specified in the normal direction. Darcy's law is used to describe the flow of Newtonian fluids in an isotropic porous medium with v andu the velocity of the fluid and solid separately; n f the porosity of the medium and k f the effective permeability coefficient of the porous medium, k f = k∕ . k and are the intrinsic permeability of the porous medium and the viscosity of the fluid separately. Using Darcy's law, the mass balance equations of the fluid and solid phases, and boundary conditions read: 17 in which p andq represent the prescribed fluid pressure and inflow, respectively. p 0 denotes the initial fluid pressure. q c ⋅ n is the inflow from the fracture, q is the fluid flux, q = −k f ∇p and M represents the Biot modulus: with K s and K f the solid and fluid bulk moduli, respectively.

POWELL-SABIN FINITE ELEMENT DISCRETIZATION OF GOVERNING EQUATIONS
In this section, the governing equations for the porous medium are solved by Powell-Sabin finite elements. We use Powell-Sabin B-splines for the trial functions in the solution space, and also for the parametrization of the geometry. Poromechanical interface elements are employed at the crack.

Weak form of the governing equations
The weak form of the balance equations (1) and (11) is obtained through multiplication by the test functions and for the solid skeleton and the interstitial pressure, respectively. Considering the internal boundaries Γ + c and Γ − c as well as the conditions at the external boundaries Γ u , Γ t , Γ p , and Γ q , and using the divergence theorem and integration by parts leads to the weak form: The interface at the crack, Γ c , introduces two terms in the weak form: cohesive tractions and the normal fluid flux through the interface faces. Taking force equilibrium conditions at both faces of the crack interface, Equation (5), we have: with n = n − = −n + . We reformulate Equation (13a) using Equation (14) ∫ If the pressure is discontinuous at the crack interface Γ c , and using the equation n = n − = −n + , the weak form of the mass balance can be rewritten as:

Powell-Sabin finite elements
To discretize Equations (15) and (16), Powell-Sabin B-splines are employed. Powell-Sabin B-splines are defined on triangles, holding  1 -continuity throughout the entire domain, even at crack tips. Other  1 -continuous triangular elements have been constructed in the past, like the Argyris element, 33 the Hsieh-Clough-Tocher (HCT) element, 34 and natural elements. 35 Powell-Sabin B-splines describe the geometry and interpolate the displacement field u and the fluid pressure p in an isoparametric sense: where X j k represent the coordinates of the corners Q j k of the Powell-Sabin triangles, see Figure 2B. U j k and p j k denote the displacement and pressure degrees of freedom (PDOFs), respectively, at Q j k , and N v is the total number of vertices. The indices j = 1, 2, 3 imply that three Powell-Sabin B-splines are defined on each vertex k, Figure 2B.
We now give a concise description of Powell-Sabin B-splines. 36 We consider a cracked domain Ω ( Figure 2A) with a triangulation  , Figure 2B.  can be generated by any package for standard triangular elements, such as Gmsh. 37 e = 1, 2, … , E triangles and N v vertices are defined on  , which is denoted by thick black lines in Figure 2B. To construct Powell-Sabin B-splines, each triangle e is split into six (n = 1, 2, … , 6) mini-triangles, compare Figure 2B. This leads to the Powell-Sabin refinement  * and Powell-Sabin points (plotted in green in Figure 2B). 36 A Powell-Sabin triangle, drawn in red, is then defined for each vertex k. 38 We constrain the Powell-Sabin triangles on the boundary as follows: (i) for an angle < 180 • , two sides of the Powell-Sabin triangle must be aligned with the two boundary edges, (ii) for an angle = 180 • , one side of the Powell-Sabin triangle must be aligned with the boundary edge. Powell-Sabin triangles are not restricted by these constrains along an internal discontinuity (crack interface), see Figure 3B. Each mini-triangle n of element e has a barycentric coordinate system = [ . 38 The Powell-Sabin B-splines N j k over each mini-triangle, Figure 2D, can be obtained using the Bézier ordinates b r,s,t : 38 where B 2 r,s,t ( ) represent the Bernstein polynomials. The Bézier ordinates b r,s,t are obtained by considering the properties of Powell-Sabin B-splines at each vertex k. 36 We can implement the Powell-Sabin B-splines in existing finite element codes by Bézier extraction with C e n a matrix filled by Bézier ordinates, N e n Powell-Sabin B-splines associated with each mini-triangle n in element e, and B Bernstein polynomials.
Powell-Sabin B-splines do not satisfy the Kronecker-delta property and are non-interpolatory at the vertex. 39 Thus, imposing essential boundary conditions on Γ u and Γ p is not as trivial as in standard finite elements. In this contribution, we will employ Lagrange multipliers to weakly impose essential boundary conditions. The Lagrange multiplier method introduces a new unknown vector . We discretize the Lagrange multipliers and the field variables , that is, the displacement u and the fluid pressure p on essential boundaries by using Powell-Sabin B-splines: with N bv total number of vertices on the essential boundary Γ u . N e b denotes the element boundary shape function matrix. e b and e b denote degrees of freedom associated with element boundaries. We take the enforcement of the displacement boundary condition u =ū as an illustration. Considering = u in Equation (20), and subsequently applying the Lagrange multiplier method, Equation (15) can be reformulated as 39 Lagrange multipliers introduce an additional unknown on the boundary Γ u , which can be interpreted as a reactive traction, that is, = (u) ⋅ n on Γ u . Linearizing Equation (21) results in the corresponding system of equations.

Poromechanical interface elements in the Powell-Sabin finite element scheme
The fluid inside the crack Γ c will induce pressures on the crack faces. To consider this effect we enhance the interface elements with porosity. The standard interface element is augmented with PDOFs. Since Powell-Sabin B-splines do not satisfy the Kronecker-delta property and are non-interpolatory at the vertex, the augmentation with pressure DOFs is not as standard as when using Lagrange basis functions. We will consider the model of two pressure degrees of freedom (2PDOF) as an example. For a 2PDOF model 2PDOF are added, one on each side of the crack, which allows a discontinuous pressure across the crack Γ c . In the 2PDOF model the pressures at each side of the crack are independent, which allows for a discontinuous pressure across the interface element, see Figure 3C. Analogous to Darcy's law, the fluid transport across the crack interface is then given as 32 with k i the interface permeability, H p is the pressure jump matrix, ; h + p and h − p are the shape function matrices associated with the crack interfaces Γ + c and Γ − c , respectively. The arrays p + and p − contain the PDOFs at both sides of the interface, andp = In 2PDOF model there is no independent fluid pressure within the crack. Consequently, the pressure term vanishes in the traction continuity condition Equation (5), and Equation (23) is rewritten as: Substituting the interface term, Equation (22), into the weak form of the mass balance, Equation (16), leads to: The mass conservation Equation (24) contains time derivatives, which are discretized using the Backward Euler scheme: In combination with the Powell-Sabin approximation Equation (17) the weak form, Equations (23) and (24), yields: where m = [1 1 0] T . Matrices N, B, and H are related to the displacement degree of freedom U, and contain shape functions, their derivatives, and relative displacements, respectively. 38 N p and B p are shape functions and their derivatives related to the pressure degree of freedom p. Linearization of Equation (26) leads to equations for the Newton-Raphson iterative scheme: with F ext and F int being external and internal force vectors, which are obtained from Equation (26). The stiffness matrices are defined as: with R the rotation matrix and T d the tangent stiffness of traction-opening law at the interface Γ c . 38

ADAPTIVE ANALYSIS FOR CRACK GROWTH
We now present the algorithm to insert a new crack segment, including remeshing after a crack insertion. Then, we will discuss how to map state vectors, displacement and fluid pressure, onto the new mesh. The mapping is performed under constraints from the energy balance and mass conservation.

Crack insertion and domain remeshing
Considering the  1 -continuity of Powell-Sabin B-splines at the crack tip, point A in Figure 4A, the crack initiation can be evaluated directly at this point by comparing the major principal stress 1 and the tensile strength t u . When 1 ⩾ t u , a crack is inserted through the entire element e 0 in front of the crack tip, see Figure 4B. The new crack tip is now at point C. Due to lack of information about the possible curvature of the crack, it is introduced as a straight line within the element. 38 The normal vector n 1 of the new crack, line AC in Figure 4B, can be obtained directly from the stress tensor at the crack tip due to the  1 -nature of Powell-Sabin B-splines. However, in this study, to further improve the quality of the prediction of the crack propagation direction, we employ an average stress tensor over a small finite domain to compute n 1 . The averaging procedure is performed using a Gaussian weight function: 40 with r the distance to the crack tip, and l a smoothing length, typically taken around three times a representative element size. The crack insertion strategy bears resemblance to the "rotation of edges" strategy, 41 but the remeshing now ensures a proper aspect ratio of the elements. After inserting a new crack segment, element e 0 is split into two triangles e 1 and e 2 , Figure 4B. Consequently, the element next to the new crack tip, e 3 in Figure 4B, has four vertices, which is not possible for the definition of Powell-Sabin B-splines. Thus we need to remesh the domain Ω c with the new crack tip, see Figure 4C. Vertices on and outside the domain Ω c and crack tips will not move. In determining the domain Ω c , we stand at the element with the newly inserted crack segment, shaded grey in Figure 4A. Then, a radial marching is done until three elements have been crossed in all directions, see Figure 4. The elements along one side of the crack interface are excluded, which avoids updating the field variables along the crack interface.
Remeshing then proceeds by first adjusting the elements in order to avoid triangles with four vertices and a bad aspect ratio, subsequently by fixing a polygon Ω c with previously adjusted elements, and finally by solving an optimization problem of the interior angles of each triangle: max j min (j = 1, 2, … ) subject to: where i k is the ith interior angle (i = 1, 2, 3) of triangle k. The remeshing procedure proceeds sequentially: we first obtain the minimum interior angle 1 min by solving Equation (29), then we further remesh the domain by using Equation (29) to maximize the second minimum interior angle 2 min of all triangles inside Ω c . We repeat the procedure until all interior angles have attained a maximum value.

State vector update after crack insertion
New elements and vertices are introduced after the insertion of a new crack segment. Remeshing Ω c is necessary to ensure elements with a suitable aspect ratio, which yields a modification of the mesh. 38 Consequently, Powell-Sabin B-spline functions must be computed on new triangles. Here, the mesh before remeshing is denoted as Ω b c , while the mesh after remeshing is represented as Ω r c . For non-linear problems, remeshing also requires a transfer or mapping of state vector, displacement and fluid pressure, from the old mesh Ω b c to the new mesh Ω r c at time step t and t + Δt, respectively. Next, we take the transfer of state vector at time step t + Δt as an illustrative case.
We first map the displacement t+Δt U from the old mesh Ω b c onto the new mesh Ω r c . We define as the displacement (fluid pressure) shape function matrix associated with the old mesh Ω b c and the new mesh Ω r c , respectively. Furthermore, U b (u b ) and U r (u r ) are displacements related to Ω b c and Ω r c , respectively. Due to the non-interpolatory property of Powell-Sabin B-splines a least-square fit is employed to carry out the mapping: In general least-square fit according to Equation (30) does not guarantee conservation of energy during the mapping from mesh Ω b c to mesh Ω r c . 42 To minimize the energy difference between Ω b c and Ω r c , Equation (30) is solved by matching the energies related to the meshes Ω b c and Ω r c , which is given as where only the internal work W int and the work W coh related to the cohesive traction on the crack surface are considered, due to their links with the displacement u.
For fracture propagation in fluid-saturated media, we also need to transfer the fluid pressure from the old mesh Ω b c to the new mesh Ω r c . This is done in a similar way to that in Equation (30): Equation (32) does not guarantee conservation of mass between the mesh Ω b c and Ω r c . To minimize the difference, Equation (32) is optimized by conserving the mass related to Ω b c and Ω r c : where only mass terms related to the pressure p are considered. In solving Equations (30) and (32) we have to consider Dirichlet boundary conditions: where Γ u c and Γ p c are boundaries with prescribed displacement and fluid pressure, respectively. Here we fix the degree of freedom on the red polygonal boundary, see Figure 4C for instance, and along the crack path Γ c .
In sum, to map the displacement and the fluid pressure from the old mesh Ω b c to the new mesh Ω r c we must solve the following optimization problems with boundary condition Equation (34): and Herein, the MATLAB function fmincon has been used to find the optimum in Equation (35). In general, the constraint equations from the energy balance and mass conservation reduce the error of the state vector update. An error analysis of the state vector update with constraint equations has been given in Reference 42.
The computational efficiency of the proposed method is lower than that in standard FEM. In the evaluation of Equation (35) we need to find the state vector of Gauss points on the refined mesh from the old mesh. 43 For each triangular element there are six mini-triangles used to perform the integration. Thus, the number of triangles used in the integration is N e × 6, and the number of Gauss points on the refined mesh is N g × N e × 6, where N e denotes the total number of triangular elements and N g is the number of Gauss integration points inside each mini-triangle. In standard FEM the number of Gauss points is N g × N e , which is smaller than that in the proposed method. Thus, the computation time in the proposed method is increased compared to standard FEM.

NUMERICAL EXAMPLES
Below we will consider three examples. The first example deals with a pre-fractured specimen (stationary crack), in order to benchmark the proposed method. The last two examples feature crack propagation under pure mode-I and mixed-mode loading conditions, demonstrating the ability of the method to analyze the propagation of cracks.

Stationary crack: Square plate with a center crack
First, a square plate with an inclined center crack is considered, see Figure 5. The tilt angle of the crack is chosen as = ∕6. A constant fluxq = 10 −4 m∕s is imposed at the bottom of the plate. All boundaries are impermeable except for the top, where the fluid is allowed to flow freely. Furthermore, the bottom edge can move freely in the horizontal direction, and the right and left edges can move freely in the vertical direction. A plane-strain condition is assumed and the crack faces are assumed to be traction-free. The material properties of the solid are given as: Young's modulus E = 9 GPa, Poisson's ratio = 0.4, Biot modulus M = 10 18 MPa, Biot coefficient = 1, porosity n f = 0.3 and the intrinsic permeability of the porous medium k = 10 −12 m 2 . The interface is assumed to be impermeable, rendering the interface permeability k i = 0. The fluid viscosity is taken as = 10 −9 MPa s. The plate has been discretized as shown in Figure 5B. A time step size Δt = 1 s is used until a steady state is reached at t = 40 s. Dirichlet boundary conditions are enforced by Lagrange multipliers as described in Section 3.2. To avoid interpenetration, a penalty stiffness k p = 10 10 MPa/mm is specified in the normal direction of the crack opening. Figure 6 shows the pressure, the flux, the displacement and the stress distribution in the steady state (t = 40 s). In Figure 6A, the fluid pressure is discontinuous across the crack. There is no flow across the interface. The fluid flows from one side of the interface to another side through the crack tips, illustrated in Figure 6B, which shows that the crack behaves like a barrier for the flow. Due to the impermeable crack interface, the displacement and the stress are also discontinuous across the crack. The pressure and displacement profile show a good agreement with the results reported by an interface element approach using Lagrange interpolations. 17

Progressive fracturing in a single-edge notched plate
A square plate with dimensions 250 mm × 250 mm with an edge crack is now considered, see Figure 7. The length of the initial crack is a = 50 mm. We set the Young's modulus 25.85 × 10 3 MPa, the Poisson's ratio = 0.18, the porosity n f = 0.2,  (8) is used with the tensile strength t u = 2.7 MPa and fracture energy  c = 0.095 N/mm. We only consider mode-I fracture, that is, d int = 0 in Equation (8). To avoid interpenetration, a penalty stiffness k p = 10 10 MPa/mm is specified in the normal direction of the crack opening. Vertical velocitieṡū = 2.35 × 10 −2 mm/s are applied at the top and the bottom sides of the plate. All boundaries are assumed to be impermeable, except for the crack interface for which the permeability, k i , is set equal to the effective permeability coefficient of the porous medium k f . The plate has been discretized by the triangulation shown in Figure 7B. The analysis is carried out with a time increment Δt = 0.01 s. Displacement control is employed to apply the velocitẏū in the simulation, Δu =̇ū × Δt = 2.35 × 10 −4 mm. The displacement is imposed by the Lagrange multiplier method, compare Equation (21).
The computed load-displacement curve is shown in Figure 8A with different number of elements. A good agreement is attained with the extended isogeometric analysis (XIGA), 8 in which only the results, before the crack along the middle plane is fully opened, are shown. The curve exhibits structural softening due to progressive crack propagation. In the figure, a close agreement is observed between the two discretizations. 8 Figure 8B presents the pressure p + along the top crack interface Γ + c after the crack is fully developed along the middle plane of the plate. In addition, a comparison between the two discretizations is given for the fluid pressure p + here. A slight difference is observed in the figure due to the influence of mesh size. Obviously, the pressure acts as a tensile stress on the crack faces and shows a nonlinear distribution along the crack path, due to the nonlinear opening of the crack, see Figure 9B. The displacement of the plate increases with time, resulting in an increased crack opening, see Figure 9B. In the process of crack propagation, the pressure will concentrate around the crack tip ( Figure 9A), due to the flow bypassing the crack tip. When the crack has fully opened, the fluid will gradually flow to the porous medium, progressively reducing the fluid pressure concentration, see the last picture in Figure 9A. Consequently, the pressure along the crack path gradually reduces to zero.

Arbitrary propagation: A plate with two propagating cracks
This example serves to demonstrate the ability of the proposed method for mixed-mode crack problems. The setup of the problem is similar to the Nooru-Mohamed tension-shear test. 44 The Nooru-Mohamed tension-shear test has been carried out on a double-edge notched plane concrete specimen with a thickness of 50 mm. 44 Figure 10A shows the  (8). 45 Plane-stress conditions are assumed and the loading condition is set up as: Step 1 Displacement control is considered to fully track the load-displacement path with steps of Δu x along the upper left and bottom right edges. The top and bottom edges are fixed in the y-direction at this stage, see Figure 10A. The total time is T = 0.6 s.
Step 2 Displacement control is employed to fully track the load-displacement path with steps of Δu y along the top and bottom edges, while keeping û x = 0.012 mm constant along the upper left and bottom right edges, see Figure 10A.
The response curve is given in terms of the vertical resultant force F y versus the vertical displacement u y on the top edge at the loading Step 2, see Figure 11A. The parameter settings k i = k f and k i = 0 give almost identical results, with a maximum difference of 1.67%, which indicates the fluid flow within the interface has a negligible influence on the response. In addition, we present the results in the case of dry material (no fluid) with identical set-ups. Almost the same results are obtained for the cases of fluid-saturated medium and dry material medium. This is also validated for the prediction of the crack path, given in Figure 11B. Obviously, the responses are dominated by the fracturing of the solid medium and the cohesive zone model. In the fracturing process, the crack opening gradually increases, which causes the decrease of the fluid pressure along the interface, as validated in Figure 8B. Consequently, the fluid pressure is considerably smaller than cohesive tractions, rendering a negligible influence on the responses in Figure 11. Figure 12 presents the fluid pressure, the flux, and the displacement distribution for the interface permeability k i = k f and k i = 0, respectively. Obviously, different patterns are observed in the fluid pressure and flux profile. The fluid pressure difference p + − p − and the flux difference q + − q − along the crack interface, in the case of k i = k f , are smaller than those in the case of k i = 0, as illustrated in Figure 12A,B. This is caused by the impermeable interface set-up in the case of k i = 0. For the displacement profile, both settings of the interface permeability k i induce almost identical patterns, as presented in Figure 12C.

CONCLUSIONS
A Powell-Sabin finite element scheme has been developed for progressive fracturing in a fluid-saturated medium. A cohesive zone law governs the interface behavior of the solid part, while a two-pressure degrees of freedom (2PDOF) model manages the fluid flow within the crack.  1 -continuous Powell-Sabin B-splines, which are based on triangles, are employed for the geometry, the displacement and the fluid pressure approximation. Due to the  1 continuity property, one can directly assess the crack initiation at the crack tip.
In the process of crack propagation, the crack is directly introduced in the physical domain, different from the discrete modeling in the isogeometric analysis framework. Due to the use of triangular elements, remeshing is straightforward. After remeshing new Powell-Sabin B-spline functions are introduced, which requires a mapping of the state vector (displacement and fluid pressure) from the old to the new mesh. A novel least-squares fit methodology has been proposed to carry out the mapping. To preserve the energy and mass in the transfer, the energy balance and mass conservation are chosen as constraints.
Numerical examples are presented to validate the proposed method and to show the refinement ability of the Powell-Sabin B-splines. For the chosen examples, the interface permeability plays a minor influence on the response, except for the fluid pressure and the flux. This is due to the fact that the fluid pressure decreases upon an increase of the crack opening. Consequently, the fracture behavior of the solid and the cohesive zone law control the crack propagation.
Due to certain constraints with neighboring tetrahedrons, the extension of Powell-Sabin B-splines to three-dimensional objects is non-trivial, 46 but one can construct prisms as a tensor product of two-dimensional Powell-Sabin B-splines and NURBS in the third dimension.
In the analysis, we employed the same basis functions for the interpolation of the displacement field and the pressure field, which does not comply with the LBB or inf-sup condition. 47 Nevertheless no oscillations in the pressure field were observed, which may be due to the higher-order continuity of the Powell-Sabin B-splines.