Loaded crack surfaces in two and three dimensions with XFEM

Applied Mathematical Modelling 78 (2020) 863–885

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Loaded crack surfaces in two and three dimensions with XFEM Markus Schätzer∗, Thomas-Peter Fries Institute of Structural Analysis, Graz University of Technology, Lessingstr. 25/II, Graz 8010, Austria

a r t i c l e

i n f o

Article history: Received 24 September 2018 Revised 4 October 2019 Accepted 8 October 2019 Available online 17 October 2019 Keywords: Crack propagation Critical pressure Superposition Data transfer

a b s t r a c t A numerical framework in two and three dimensions is proposed for crack propagation induced by loaded crack surfaces. The method is based on the XFEM with a hybrid explicitimplicit crack description. The load on the crack surface is given on the explicit surface mesh and transferred to the zero-level set of the implicit description. Therefore, integration points must be identified on the zero-level set which are also used to determine stress intensity factors by crack opening displacements. A crack propagation model is developed where the load factor of the crack surface loading is determined such that the criterion for crack propagation is met. Numerical results show the versatility and success of the approach. © 2019 Elsevier Inc. All rights reserved.

1. Introduction In the analysis of fractures one often assumes stress-free crack surfaces [1–3]. However, there are many applications where the crack surface is loaded, in fact, this load may be the dominant driving force for the potential propagation of the crack. Examples are corrosion and freezing processes in fractured structures. Another example is hydraulic fracturing (HF) for reservoir stimulation where a fracking fluid is pumped into some rock formation to induce fractures and increase the permeability and, thereby, the productivity [4–6]. Loaded crack surfaces similar to HF also occur, e.g., in dam break scenarios or in magma-driven dykes [7]. Hence, the investigation and numerical analysis of cracks under loading conditions which apply to the crack surface itself is important. For a general assessment, it is very useful to assume sharp crack surfaces where the crack surface is identified by a surface mesh or level-set function and, consequently, the locations where the stresses are acting are clearly identified. This is in contrast to smeared crack representations as in phase-field methods [8–12] where the consideration of general loading conditions within cracks is more difficult. Additionally, these methods are often mesh-dependent and require extremely fine meshes to resolve the sharp discontinuity across the crack surfaces [13]. The extended finite element method (XFEM) has developed to be a standard tool in fracture mechanics with a sharp crack representation [14–17] and is also used herein. Many research activities have been done to improve the accuracy and efficiency of the XFEM in fracture mechanics. Adaptive local mesh refinement XFEM techniques ensure an accurate modeling of the crack front, see, e.g., [18–22]. In [23], an XFEM approach has been developed for modelling multiple branched cracks, voids and cracks emanating from holes. A junction enrichment was introduced in [24] to describe intersecting fractures in deformable porous media. Intersections between hydraulic and natural fractures are investigated in the context of XFEM in [25]. An XFEM formulation with higher-order elements for curved cracks was introduced, e.g., in [26,27], however, the ∗

Corresponding author. E-mail addresses: [email protected] (M. Schätzer), [email protected] (T.-P. Fries).

https://doi.org/10.1016/j.apm.2019.10.020 0307-904X/© 2019 Elsevier Inc. All rights reserved.

864

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

higher-order accuracy may only hardly be extended to propagating and three-dimensional cracks, wherefore we restrict to linear elements herein. Kang [28–30] applied a consecutive interpolation procedure in XFEM which provides smoother stress distributions also in the vicinity of the crack tip leading to improved stress intensity factors (SIFs). The modeling of cohesive cracks with the XFEM was investigated, e.g., in [31–34]. For further literature on the XFEM, we refer the interested reader to the overviews given in [35,36] and text books [37,38]. We also note that there is a wide range of other approaches which may be used in the context of fracture mechanics with sharp crack representations. For example, meshfree methods [39–42], boundary element methods (BEM) [43,44], scaled boundary finite element methods (SBFEM) [45–47]. One may also find combinations of the isogeometric analysis (IGA) and the (i) BEM [48–51], the (ii) XFEM [52,53], the (iii) SBFEM [54] and (iv) meshfree methods [55]. We adopt the XFEM with a hybrid explicit-implicit crack description from [56,57]. This enables a very general setup for loaded crack surfaces: On the one hand, an explicit crack representation based on surface meshes is used which simplifies the definition of the crack update during propagation. It also enables a very general framework to consider for general loadings on the crack surfaces: The stress fields may be defined on the surface mesh or even result from solving additional model equations on the surface mesh. On the other hand, also an implicit crack representation based on zero level-sets is generated from the explicit one which is useful for addressing major XFEM-related issues, most importantly (i) the definition of enriched nodes, (ii) the enrichment functions based on customary coordinate systems, and (iii) the numerical integration. The present work focuses on providing a general numerical framework for crack propagation in two or three dimensions with arbitrarily curved and loaded crack surfaces. The focus is on interpolating and integrating the stress fields for which data are transferred from the explicit to the implicit crack description. Furthermore, a very general procedure for the extraction of SIFs is presented which is based on fitting the approximated crack opening displacements (CODs) to the expected displacement fields of the individual crack modes. This is not trivial because the displacement fields are determined based on the implicit description, that is, for the fitting, it is required to explicitly locate points on the zero-level sets. Finally, we present a new approach to extract a crack propagation criterion for which the overall loading is seperated into two states: one is the load on the crack surface, the other the remaining load including traction and displacement boundary conditions. It is then determined by which factor the crack loading may be increased until the crack propagates. Because the most important field of application of loaded cracks seems to be in HF, we relate to this topic frequently, although the context herein is more general and the numerical results are more academic than this concrete field of application. For an overview of XFEM-based approaches in the context of HF, see, e.g., [21,58–60], where the focus is often more on the modeling and application side than on the numerical treatment. We note that HF is typically modeled as a coupled multi-physics problem with (at least) three different fields: (i) a model which represents the deformation of the rock, typically considered as a linear-elastic body, (ii) a model for crack propagation, and (iii) a model for the fracking fluid inside the fracture, typically represented by the Reynolds equation. A coupled formulation of the problem generally leads to a non-linear and time-dependent system of equations, see, for example, in [59,61–63] where the problem is discussed in the context of porous media. It is obvious that HF (i) involves a large number of physical processes (opening of existing fracture networks versus creation of new crack surfaces, fluid lag and leak-off, energy release by toughness and viscosity), (ii) takes place on various time and length scales which are not easily bridged/upscaled and (iii) the uncertainties are immense with existing experimental and field data being sparse. This makes the modeling and simulation an extremely challenging task [64]. A current trend is to resolving these processes and influencing factors individually to judge on the importance in practical applications of HF. For example, the influences of different asymptotic behaviours in the crack tip region are investigated in [60,65] with respect to different propagation regimes. However, in the light of the comments made above, in particular the lack of data for validation, one may also think of coarser models which summarize physical phenomena in a largely simplified manner. These models could provide rough relations between the structural analysis (deformations, crack openings and propagation) and the exerted loading on the crack surface so that no explicit fluid model is, in fact, needed. It is noted that even such models fall into the context of the proposed numerical framework and, in fact, we suggest how such a largely simplified model may be set up. A major advantage of the proposed method is the modularity which allows to improve individual aspects of the model and provides interfaces to other disciplines. For example, the use of arbitrary tractions on the crack surfaces allows the use of advanced fluid models in HF. However, we emphasize that this work is not to be seen as a major contribution to the modeling of HF but to the numerical treatment of arbitrarily loaded crack surfaces in general. The original contributions are mostly (i) considerable improvements of the hybrid explicit-implicit XFEM concerning the definitions of level-set functions and resulting coordinate systems, (ii) the two-way data transfer between the explicit crack surface mesh and the implicit level set functions (iii), the extraction of SIFs in 2D and 3D based on CODs in the XFEM, and (iv) the crack propagation criterion for loaded crack surfaces based on the load factor method and superposition of loading scenarios. We find applications in HF are representative for situations where very general and complex pressure profiles are present and, therefore, frequently relate to this application. Many other applications with loaded crack surfaces feature simpler pressure distributions such as linear or constant functions. An outline of the paper is as follows: Section 2 briefly repeats the XFEM approximation of a linear elastic fracture mechanics (LEFM) problem. Herein, the definition of the enrichment functions based on customary coordinate systems is discussed in the context of a hybrid explicit-implicit crack description. Section 3 deals with loads applied on the crack surfaces which may be defined on an explicitly defined surface mesh. A focus is placed on the surface integration of the fluid pressure on a purely implicitly defined crack geometry which also demands a discussion of data transfer between two non-matching

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

865

meshes. In addition to the general discussion of this problem, investigations of pressure distributions in the context of HF due to different propagation regimes are done, yet for rather academic setups. In Section 4, the propagation of crack surfaces is discussed. Herein, the evaluation of a critical pressure magnitude and its corresponding time is investigated based on the superposition principle of LEFM. Another chapter of this section is the computation of SIFs. Verifications and numerical results of the proposed method can be found in Section 5. Finally, Section 6 concludes this paper. 2. XFEM with a hybrid explicit-implicit crack description In this section, the solution of a linear elastic fracture problem with pressurized crack surfaces with the XFEM and a hybrid explicit-implicit crack description is described based on [56,57]. The involved coordinate systems are discussed which simplify the detection of the enriched nodes and are basis for the definition of the enrichment functions. It is shown how the enrichment functions are defined based on these coordinate systems without any inconsistencies and discontinuities between the cut element edges. Then the computation of SIFs for an arbitrary three-dimensional mixed-mode loaded crack configuration is discussed, which works analogously for all three modes. 2.1. XFEM approximation The standard finite element method (FEM) offers a powerful numerical approach for solving PDEs, however, it requires smooth solution features for an optimal convergence rate. Singular stresses as they appear at the crack front in LEFM, however, hinder optimal convergence. The FEM also requires a suitable mesh that accounts for the independent apertures of √ the upper and lower crack faces and the r-behaviour of the displacement field in the vicinity of the crack front. That is, local mesh-refinements are required at the crack front which, in the context of crack propagation, usually demands frequent remeshings. The XFEM is an extension of the FEM where additional enrichment functions enable the approximation of inner-element discontinuities and singularities with optimal accuracy only based on uniform, regular background meshes [14,15,35]. For LEFM problems the following enriched approximation has proven useful [15]

uh ( x ) =



Ni (x )ui +



Ni (x )ψstep (x )ai +

I∗

I

4  

j Ni (x )ψtip (r, θ )bij .

(1)

j=1 J ∗

The first term on the right hand side represents the standard finite element approach while the other terms represent the non-smooth features of the solution. Herein, we use meshes composed by linear tensor-product elements but most conclusions also apply to simplex elements. ψ step is the Heaviside function, which is used as the global enrichment function to consider the discontinuous displacements along the crack surfaces. ψtip are the global enrichment functions to consider the singular stresses at the crack front, j



j ψtip (r, θ ) =

√ r sin

θ 2

,

√ r cos

θ 2

,

√ r sin

θ 2

sin θ ,

√

r cos

θ 2



sin θ .

(2)

These functions are often based on a polar coordinate system which aligns with the crack surface and has its origin at j the crack front. Ni are the standard finite element shape functions, ai and bi are additional degrees of freedom at the corresponding enriched nodes I∗ and J∗ . Either enriched nodes J∗ are identified as element nodes whose elements contain the crack-tip/front or as set of element nodes within a fixed enrichment area near the tip/front. These two enrichment schemes differ in terms of achievable convergence rates and the condition numbers of the corresponding system of equations [27,66]. Better convergence properties are achieved with a fixed enrichment area (geometrical enrichment) which is independent of the element sizes. However, such a procedure negatively influences the resulting condition number. Many research activities have been done to stabilize the system of equations, see [67–71]. We prefer a topological enrichment herein where only the crack front elements are enriched with (2). In contrast to the crack tip enrichment functions defined in Eq. (2), Duarte [72] and Malekan [73] proposed the use of enrichment functions which are derived numerically by means of global-local techniques. Instead of the asymptotic neartip functions, discontinuous linear enrichment functions were investigated for three-dimensional fracture problems in [71]. More recently, an improved XFEM (IXFEM) was introduced in [74] which utilizes interpolative least-squares fitting instead of extra degrees of freedom to characterize the singularity at the crack front. As the resulting shape functions of the enriched approximation according to Eq. (1) do not have Kronecker-δ property it j follows for non-zero values of ai or bi that uh (xi ) = ui which complicates the imposition of Dirichlet boundary conditions and the interpretation of the solution. A shifted version of the standard XFEM approach was introduced in [75] which restores Kronecker-δ property with a minimum effort. This shifted version modifies Eq. (1) in order to maintain the desired properties as follows

uh ( x ) =

 I

Ni (x )ui +

 I∗

shift Ni (x )ψstep ( x )ai +

4   j=1

J∗

j,shift Ni (x )ψtip (r, θ )bij ,

(3)

866

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

Fig. 1. (a) Explicitly defined crack and and selected isosurfaces of the three level-set functions where (b) φ1 = 1, (c) φ2 = 1, (d) φ3 = 0; 1; 2.

with shift ψstep = ψstep (x ) − ψstep (xi ),

and

j,shift j j ψtip = ψtip (r, θ ) − ψtip ( r i , θi ) .

(4)

We also note that other versions of the XFEM may useful, e.g., the corrected XFEM [76] or the stable XFEM [77] which avoid issues in blending elements. 2.2. Hybrid explicit-implicit crack description An explicit crack description, for example, by means of straight line-segments for crack paths in 2D or flat triangles for crack surfaces in 3D has the advantage that the update of the crack geometry during the propagation is simple as new segments can be added explicitly after a propagation step. This description also enables a very general framework to consider for general loadings on the crack surfaces: The stress fields may be defined on the surface mesh or even result from solving additional model equations on the surface mesh. However, a purely explicit description complicates the definition of the enrichment functions and the detection of the enriched nodes in the XFEM. On the other hand, an implicit description, by means of level-set functions simplifies the detection of enriched nodes and the definition of the enrichment functions and is the standard in most XFEM-applications. However, a purely implicit description requires an effort for updating the crack geometry during the propagation as, for example, a Hamilton-Jacobi equation has to be solved [78] which may require a sophisticated FEM approach, e.g., due to a necessary stabilization of the governing weak form. A fast marching method was used in [79,80] for updating level-set functions during propagation. The work of Duflot [81] provides an overview of different level-set approaches for the description and update of crack geometries. In [56], a method was introduced which combines the advantages of both descriptions: The crack geometry is firstly determined explicitly, by means of straight linesegments in two dimensions or flat triangles in three dimensions, see Fig. 1(a). Then, level-set functions are derived from this master configuration based on distance computations. That is, all benefits of level-sets in XFEM can be used during the simulation while the crack propagation is considered explicitly. For the description of the crack geometry, three level-set functions are used, see [56] for details: • • •

φ1 (x ) is the (unsigned) distance function to the crack path/surface. φ2 (x ) is the (unsigned) distance function to the crack tip/front. φ3 (x ) is a signed distance function to the extended crack path/surface.

A graphical representation of these functions is given by some example level-sets in Fig. 1(b–d). To minimize the computational effort and required memory space, level-set functions are just computed in the element nodes and are interpolated within the elements by the corresponding shape functions

φ hj (x ) =



Ni φ j (xi ).

(5)

It is noted that the interpolation always occurs in the reference element of the physical element. In general, this would lead to non-planar crack descriptions within the cut tensor-product elements. As a simplification we subdivide each cut element into several simplex elements which leads to piecewise straight or planar crack geometries, see Section 3.1. 2.3. Coordinate systems based on level-set functions Based on the level-set functions φ 1 , φ 2 , and φ 3 a local coordinate system (a, b) is implied which simplifies the detection of the enriched nodes. Herein, the local coordinates are defined by two scalar functions



a (x ) = ±

φ 2 ( x ) 2 − φ3 ( x ) 2

(6)

and

b ( x ) = φ3 ( x ) .

(7)

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

867

Fig. 2. The coordinate r (x ) in the crack tip element, (a) using r = φ2 interpolated from the nodes which is not adequate, (b) evaluating (a, b) at the nodes based on φ i , interpolating within the element and computing r based on Eq. (8) gives a better result when compared to the exact solution shown in (c).

The sign of a(x ) depends on the position of point x. It is negative for points on the crack surface or behind the crack front and positive for points in front of the crack front, see [56] for further information. This coordinate system is also related to the polar coordinate system (r, θ ) which is used for the definition of the crack tip/front enrichment functions of Eq. (3). The relation is given as

r (x ) = and

θ (x ) =



a ( x ) 2 + b( x ) 2

(8)

⎧ −1 b(x)  tan for: a > 0; b = 0 ⎪ a (x ) ⎪ ⎪

 ⎪ b ( x ) −1 ⎪ ⎪ + π for: a < 0; b = 0 ⎨tan a (x ) π

2 ⎪ ⎪ ⎪ 3π ⎪ ⎪ 2 ⎪ ⎩

0

for: a = 0; b > 0

(9)

for: a = 0; b < 0 for: a > 0; b = 0.

Similar to the interpolation of the level-set functions within the elements also the local coordinates (a, b) and polar coordinates (r, θ ) can be interpolated based on the element shape functions. However, in crack tip enriched elements it is recommended to interpolate the polar coordinates based on the local coordinates (a, b) rather than by the level-set functions as, otherwise, the unsigned distance functions (φ 1 , φ 2 ) lead to some issues. For example, in the element containing the crack tip, the coordinate r (x ) should be zero at the crack tip. However, as the level-set values at the element nodes are all further away from the crack tip they are all significantly larger than 0. An interpolation will only yield values inbetween the nodal values, i.e., the desired 0 at the crack tip is not achievable. In contrast, when the local coordinates (a, b) are evaluated at the nodes based on Eqs. (6) and (7) and then interpolated within the elements, it is possible to determine the polar coordinates (r, θ ) based on these interpolated coordinates (a, b) through Eqs. (8) and (9) which will yield the desired r = 0 at the tip. See Fig. 2 for a comparison. The required global derivatives of the polar coordinates can be easily computed with the global derivatives of the shape functions which is a standard procedure in the FEM.

f,ih (x ) =



N j,i f (x j ).

(10)

However, this leads to discontinuous derivatives between the element edges due to the C0 -continuity of the shape functions. Therefore, it is proposed to determine the ‘exact’ derivatives in the element nodes and interpolate them within the element by the corresponding shape functions

f,ih (x ) =



N j f,i (x j ).

(11)

A related setting is used in [82] where an average nodal gradient avoids the discontinuities between the elements. According to the chain rule this interpolation scheme leads to the following formulation of the global derivatives of the polar coordinate r (x ) and θ (x ):

dr dr da dr db = · + · dxi da dx db dx   Nja j Njb j dr = · N j a j,xi + · N j b j,xi dxi r( N j a j , N j b j ) r( N j a j , N j b j ) dθ dθ da dθ db = · + · dxi da dx db dx

(12) (13) (14)

868

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

(a)

(b)

10 -1

10

-2

10

-3

10 -4 10 -5

10

-6

10

-7

10 -8 10 -3

(c)

discontinuous derivatives smooth derivatives

10 -2

10 -1

10 0

(d)

Fig. 3. Interpolated dr/dy based on (a) discontinuous derivatives obtained from Eq. (10), (b) smooth derivatives obtained from Eqs. (13) and (15), (c) correct derivatives, and (d) convergence.

  − Njb j Nja j dθ = 2 · N j a j,xi + 2 · N j b j,xi . dxi r ( Nja j, Njb j ) r ( Nja j, Njb j )

(15)

Fig. 3 (a and b) shows the interpolated derivatives of the polar coordinate r (x ) due to the standard interpolation (Eq. (10)) and the proposed interpolation (Eq. (13)). The ‘correct’ derivatives according to continuous level-set functions φ i (contrary to φih ) are illustrated in Fig. 3(c). The interpolation error is quantified by means of the relative L2 -norm over the complete domain  as follows:

 

 f ref − f h 2 d ε2 =   2 . f ref d

(16)

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

869

Fig. 4. Possible shapes of the zero level-sets within one sub-element in (a) two dimensions, and in three dimensions when (b) one node has a different sign or (c) two nodes have the same sign.

Herein, fref represents the reference results which are achieved when the ‘exact’ level-set values by means of no interpolation within the elements are used in all points. fh represents the approximated results of the two interpolation schemes which are based on Eq. (10) or (13) and (15). The convergence rate for both interpolation schemes are shown in Fig. 3(d), where the dashed line represents the optimal convergence rate. These results show that the proposed interpolation scheme provides a very good accuracy. This, together with the generation of the polar coordinate system as mentioned above, may be seen as an improvement over the standard XFEM with hybrid explicit-implicit crack description from Fries and Baydoun [56]. 3. Loading on the crack surface 3.1. Integration of the load on implicit crack surface In this section, the surface integration of some arbitrary stress fields on the crack surface is described. Later, some examples for realistic, general and complex loadings are discussed in the field of HF. The focus is placed on the surface integration on a purely implicitly described crack geometry. A goal is the detection of the integration points on the implicitly defined crack geometry and the modification of the corresponding integration weights. In general, the interpolation of the level-set functions within the elements with the corresponding shape functions leads to curved zero level-sets in bilinear quadrilateral or trilinear hexahedral elements which complicates the detection. Therefore, a decomposition of the corresponding reference elements into simplex elements is employed. This simplifies the detection of the zero level-sets as they are straight or planar within simplex elements. The decomposition of cut elements into polygonal subcells has become a standard in the XFEM and is often used for the numerical integration of the weak form, e.g., in [15,17,35]. A similar setting was introduced for higher-order elements in [27,83,84]. In this paper, a quadrilateral element is subdivided into two triangles and a hexahedral element is subdivided into six tetrahedrons which lead to three possible shapes of the zero level-sets within each simplex element: a straight line in two dimensions or a flat triangle or quadrilateral in three dimensions, see Fig. 4. That is, the dimension of the zero level-sets is one order lower than the problem itself and can be described with (d − 1 )dimensional standard finite elements which live in Rd . Because these elements describe the crack geometry they are called ‘crack elements’ and build the basis for the surface integration. Then, the surface integration requires four steps which are explained for a linear hexahedral element, representatively. (i) (ii) (iii) (iv)

Define the integration points r i in the reference crack element. Map these points to the reference tetrahedral element. Map the reference tetrahedral element to the corresponding sub-tetrahedron in the reference hexahedral element. Map the reference hexahedral element to the physical element of the mesh.

It is noted that the mapping of the reference hexahedral element to the physical element of the mesh in step (iv) with the trilinear shape functions leads to curved zero level-sets in general. This procedure applies analogously to tetrahedral elements where the situation is simple because the zero level-sets stay naturally flat even in the physical element. The general procedure for generating integration points on the implicit crack surface is illustrated in Fig. 5 for an arbitrary three-dimensional hexahedral element. As the shape functions are discontinuous along the crack surface, the integration points have to be splitted into two opposite but infinitesimally close points being on either side of the crack surface. This is done between steps (iii) and (iv) of the surface integration based on the normalized normal vector nξ and the shifting magnitude ε ξ

ri ± (ξ i ) = r i (ξ i ) ± ε ni .

(17) ξ

ξ

This normal vector is given by the cross product of the two tangential vectors t 1 and t 2 which are extracted from the roots of the zero level-sets at the element edges. Based on these tangent vectors t and the Jacobi-matrices J, the modification of

870

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

Fig. 5. Mapping of the Gauss points for the integration of the loading on the crack surface in three dimensions.

Fig. 6. Data projection from the explicit to the implicit crack geometry.

the integration weights ωi is given by

ωix = ωia ·



det (JTr Jr ) ·

||Jξ t r1 × Jξ t r2 || ||Jxt ξ1 × Jxt ξ2 || · . ||t r1 × t r2 || ||t ξ1 × t ξ2 ||

(18)

It is noted that the first mapping requires the Gram’s determinant as a (d − 1)-dimensional element is mapped to Rd . The other mappings take into account that an area is projected instead of a volume. These integration points are also the basis for the computation of SIFs which is discussed in Section 4.1. A discussion about the behaviour of the pressure along the crack surface is done in Section 3.3 in the context of HF. 3.2. Data transfer between the explicit and implicit crack description As discussed in Section 2.2, a hybrid explicit-implicit crack description has many advantages over a purely implicit or explicit description. However, this combined formulation leads to two different crack geometries which may not coincide exactly, wherefore, data has to be transferred between these descriptions. 3.2.1. Projection from the explicit to the implicit crack geometry Data may be often given on the explicitly defined surface mesh, for example, in the context of HF, as a result from solving flow models on the crack surface. The implicitly defined zero level-sets are basis for the surface and domain integration and therefore also for the positioning of the integration points x i , see Section 3.1. As mentioned above, this implicitly defined crack geometry may not coincide exactly with the explicit description which complicates the interpolation of the required functions. We describe how to project data from the explicit surface mesh to the integration points x i of the implicit crack geometry. In the following, r i represent integration points within the reference element and x i represent the mapped integration points in the global Cartesian coordinates of the mesh. The detection of the zero level-sets and the mapping of r i to x i was already discussed in Section 3.1. The task is now to obtain a proper function value at x i based on given nodal function values at the nodes of the explicit surface mesh. Therefore, the closest point on the surface mesh has to be identified and labelled x˜ i . Such a closest point may be inside a surface element, on the edges of a surface element or at the nodes of the explicit mesh. The corresponding function value at x˜ i has to be interpolated from the nodal values and is used as the representative function value at x i . Such a projection is, e.g., known from fluid-structure interaction as described in [85]. Examples are seen in Fig. 6 where the closest point x˜ 1 is within a triangle, x˜ 2 on an edge, and x˜ 3 at a

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

871

1.5

0.6

1 0.4

0.5 0.2

0 -0.5

0

-1 -0.2

-1.5 -0.4

-2 0

0.2

0.4

0.6

0.8

1

0

(a)

0.2

0.4

0.6

0.8

1

(b)

Fig. 7. First order approximation of the pressure m0 for a viscosity-dominated (a) KGD fracture, and (b) penny-shaped fracture.

vertex, respectively. Since the distance between the explicit and implicit crack description is generally very small and the changes of the explicitly defined functions are often moderate over the crack surface it is seen that this approach provides satisfactory results in general. 3.2.2. Projection from the implicit to the explicit crack geometry A data transfer from the implicit to the explicit crack description may also be desired from the implicit to the explicit description. An example is found in HF where the crack width based on the XFEM-approximation may be needed for the flow model on the crack surface. A very simple and fast procedure is based on a nearest neighbour interpolation [86] where the required information are obtained from the closest point of the other mesh. For example, displacement quantities such as openings are computed at different points on the implicit crack surface in a post-processing of the XFEM-approximation. These points are mapped with their stored information to the global coordinate system. For all required points on the explicit surface mesh the value of the closest point to this point cloud is used. It is noted that the maximum element size of the surface mesh and the background mesh should be in a similar range to be able to represent physically meaningful quantities on both meshes. 3.3. Pressure distributions in the context of HF Because HF is among the most important fields of applications of this work, some typical loading conditions in this case are discussed here. HF processes are usually described as non-linear coupled problems, where the fluid pressure p is related to the crack width w. This relation can be expressed by a non-linear model known as Reynolds equation [87]. It is obtained by combining Poiseulle’s law for fluid flux with the continuity equation. Injection rates are assumed to be small enough that the assumption of a laminar flow subjected to friction between two infinite parallel plates is valid. Based on the Reynolds equation the fluid pressure mainly depends on the injection rate Q0 , the crack width w, the fluid density ρ , the fluid viscosity μ, the fluid velocity v, the leak-off qL , the existence of a fluid lag lf and the time t. For simple crack geometries, resulting pressure distributions are given in the literature according to different propagation regimes, for example, in [88]. It was observed that the assumption of a uniform pressure is suitable for large toughness (zero viscosity) problems [88]. On the other hand, zero toughness (viscosity dominated) problems lead to singular pressures at the crack tip. Dimensionless pressure distributions are given for a viscosity-dominated propagation based on the KGD-model and a penny-shaped crack, for example, in [89,90]. The first order approximations of the dimensionless pressure distributions due to a viscosity-dominated propagation for the KGD-model can be written as [89]

KGD m0 =

√ 1 B(1/2, 2/3 )[ 3 ·2 F1 (−1/6, 1, 1/2, ρ 2 ) 3π −10/7 · 0.156 ·2 F1 (−7/6, 1, 1/2, ρ 2 )] + 0.0663(2 − π ρ ),

(19)

and for a penny shaped crack as [90]

Penny = 0.795 − m0

0.239

( 1 − ρ ) 1/3

− 0.0911 · ln(0.5ρ ).

(20)

Herein, B is the Euler beta function, 2 F1 is the hypergeometric function, and ρ is the scaled coordinate which is 0 at the injection point and 1 at the crack tip/front. Fig. 7 shows a graphical illustration of Eqs. (19) and (20).

872

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

The coupling of the elasticity and fluid problem is often associated with numerical problems often leading to a large number of iterations of giving rise to instabilities [6]. Furthermore, due to the fact that the evaluation of the pressure is computationally expensive and there are only few experimental data available for the validation, there is considerable motivation to simplify this delicate coupled problem. In the light of sparse experimental data this may be even more desirable. Furthermore, in most industrial applications of HF the focus may rather be on the resulting shape of the fracture, pressure and fluid volume over time and less complicated models may serve this purpose as well. In this contribution, we introduce a simplified model which has the advantage that it is very robust and computationally cheap because no iterations between different models, e.g., fluid-structure interactions, are required. It is noted that the proposed model is not restricted to a specific application but can be used for a wide range of applications with loaded crack surfaces. Herein, the simplified model is discussed in the context of HF, respectively. In HF, one may suggest to replace the Reynolds equation by pre-defined pressure distributions similar to those from above but scaled by only few parameters which increases the robustness and decreases the computational efforts. Herein, the assumed distribution is based on the knowledge of the dominated regimes (toughness and viscosity) which can vary during the time. The aim is to develop a simple, robust and computationally cheap simulation tool for quasi-static HF processes which enable a satisfactory prediction of the obtained crack configuration, in particular to the corresponding crack geometry or volume and the corresponding pressure over time t. We assume that all pressure distributions during the propagation are reproduced as a linear combination of the toughness-dominated distribution p( p0 , ρ ) = p0 and the viscosity-dominated case p( p0 , ρ ) = p0 · m0 . The influence of each regime is considered by a time-dependent scaling function 0 ≤ (t) ≤ 1 which represents the dominance of the toughness-dominated regime. That is, the propagation is with (t ) = 1 purely toughnessand with (t ) = 0 purely viscosity-dominated. Based on this scaling function the current pressure distribution at a time t is given by

p( p0 , t ) = p0 [ (t ) + (1 − (t ) )m0 ].

(21)

It is noted that this kind of pressure distribution should only highlight the construction of some general yet physically motivated pressure distributions as possible loadings on crack surfaces. In this work, different distributions are used for the validation of the proposed method and the focus is here rather on the availability of reference solutions than practical examples of industrial relevance. The critical pressure magnitude which leads to a crack propagation, see Section 4, is determined based on the superposition principle of LEFM. It is important that the original Reynolds model may also be used within the proposed approach. However, as this may induce numerical problems we rather prefer to allow for arbitrary yet user-prescribed pressure profiles in this approach as mentioned above. The transition between the different regimes may itself be a (time-dependent) function which is an additional free parameter in the proposed approach. 4. Crack propagation In LEFM, the decision whether a crack propagates is often based on the stress intensity factor KI or the material toughness KIc [2]. Herein, it is assumed that an unstable, planar propagation (without deflection) occurs when the current mode I stress intensity factor KI reaches the material toughness. Threshold definitions such as those used for fatigue fractures, i.e., Kth enable stable crack growths which already occur below the critical state [91]. It is noted that the proposed method is not limited to a specific propagation criterion. For simplicity, KIc is used as a limitation value herein. This criterion takes the form

KI = KIc ,

(22)

and only considers mode I loadings. For an arbitrary three-dimensional crack configuration a more general formulation is based on the energy release rate G which is a measure of the energy available for the creation of new crack faces. Similar to Eq. (22) a propagation occurs when the current energy release rate G reaches a critical value Gc

G = Gc .

(23)

Techniques for the computation of G are, for example, the path-independent J-integral [92] and the interaction integral [93]. In LEFM, the J-integral is equivalent to G [82] which enable a simpler evaluation, as a line-integral has to be solved instead of a contour integral, however, the extension to three dimensions causes some difficulties [93]. An alternative is to use the relation between the energy release rate and the SIFs which is given by



G(Km ) =

1 − ν2 E







KI2 + KII2 +

1 2 K , 2μ I I I

(24)

where E is the Young’s modulus, ν the Poisson’s ratio, μ the Lamé-Constant and Km the stress intensity factor of mode m. Although, all modes are considered in this representation of G it assumes a self-similar crack growth, for example, a planar crack is assumed to remain planar and maintain a constant shape as it grows [1]. For fluid-driven fractures this is often given by the fact that the fluid pressure dominates the propagation against other loadings. The proposed method is not limited to this restricted relation, wherefore, also formulations of infinitesimally kinked crack extensions [1,94], are possible. Nevertheless, we are using Eq. (24) for reasons of simplicity. In three dimensions, the fact that each point along the crack front has its own behaviour complicates the assessment of the propagation as most of the time Gc is only reached at one

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

873

Fig. 8. Coordinate system at the crack front.

point of the crack front. This would mean that only this point propagates which is not very meaningful in simulations. For a more natural propagation we assume that the critical point where Gc is reached grows with an increment length rmax and all other points grow with scaled lengths. These scalings are based on the corresponding ratios of the current energy release rates G and the critical Gc which is also used in [95]. For this it follows that a point i with an energy release rate Gi grows with an increment length ri of

r i =

Gi · rmax . Gc

(25)

For keeping the number of observed points reasonable we restrict the assessment of the propagation to crack front nodes which are the nodes of the explicit crack surface mesh at the front. Their connection forms the three-dimensional C0 continuous polygon which bounds the explicitly defined polyhedral crack geometry of Section 2.2. For the evaluation of the propagation direction it is often useful to know what a planar crack extension or an extension out of a certain plane, for example, due to mode III means. This is often based on a local coordinate system which is defined by three vectors n, t and q, see Fig. 8. These vectors represent the normal vector n of the crack surface, the tangential vector q of the crack front and the co-normal vector t which is normal to these two vectors and can be interpreted as a planar extension of the crack. The definition of these directions within a line-segment of the crack front is quite easy as they only depend on the corresponding triangle element. However, at crack front nodes the definition can not be based on a single triangle element as it would lead to discontinuities. A smooth coordinate system is achieved by averaging the required vectors of the corresponding elements. That procedure is, for example, used in [56] where the normal vector n is averaged through the element areas and the tangential vector q is averaged with the lengths of the neighboring element edges. In the literature, a wide range of reliable and accurate propagation criteria exists for the evaluation of the propagation direction. Common criteria are, for example, the maximum circumferential stress criterion [96], the maximum strain criterion [97], the maximal strain energy release rate criterion [98], and the minimum strain energy density criterion [99]. In this contribution, the evaluation of the propagation direction is based on the maximum circumferential stress criterion [96] where the propagation angle θ c can be evaluated in the n, t -plane based on SIFs [15]. When KII = 0 this criterion leads to following expression



θc

1 = 2 arctan 4

KI − sign(KII ) · KII



 K 2 I

KII



+8 .

(26)

When KII = 0 the crack propagates straight which means that the propagation angle θ c is zero. The propagation out of the n, t -plane due to mode III is neglected in this criterion. 4.1. Computation of SIFs In this section, the computation of SIFs is discussed for an arbitrarily curved crack in three dimensions. Three independent crack modes fully characterize the situation at the crack front and may be used to determine if and how the crack propagates. For example, the SIFs may be used to model the crack growth speed and direction and are basis for the decision whether a crack propagates. A challenging task is the computation of these factors for arbitrary non-planar three-dimensional crack geometries. One approach is based on the interaction integral [82,100,101] which, however, is cumbersome in three dimensions and prone to numerical problems. A more intuitive procedure is based on CODs where the approximated aperture of the crack, which is extracted from the XFEM simulation, is compared with a reference state in which the expected openings for a pure mode I, II and III are known [102–104]. One advantage of the approach over the interaction integral is that only displacements are needed and stresses do not have to be determined. Another advantage of this procedure is that it works for two- and three-dimensional planar or non-planar crack geometries in a similar manner. Due to the fact that this procedure is simple and robust it is used in this work. In [105], it is shown how to use this method in the context of LEFM with the XFEM and a hybrid implicit-explicit crack description. Herein, the integration points of Section 3.1 are used to evaluate the approximated CODs in the vicinity of the crack tip. These openings are compared with

874

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

Γu

Γt

Γu

=

p

Γu

+

Γt

(a)

(b)

p

(c)

Fig. 9. Superposition of an externally and internally loaded crack by means of the principle of superposition. (a) Total load, (b) external load, (c) internal load.

the reference state defined by the local coordinates a and b which were described in Section 2.2. A third direction c is defined, which is normal to a and b by the cross product of ∇ a and ∇ b. It is necessary to describe both, the approximated and the expected CODs in the same coordinate system, wherefore, the approximated CODs which are given in the global coordinate system (x, y, z) are transformed into the local coordinate system (a, b, c). This transformation is based on the Jacobi-matrix J.

⎡

uha = J · uhx

with

a,x

a,y

a,z

⎤

J = ⎣b,x

b,y

b,z ⎦.

c,x

c,y

c,z

(27)

It was shown that KI mainly leads to displacements in direction of b, KII to displacements in direction of a and KIII to displacements in direction of c. Therefore, SIFs are directly computed by using these selected directions with

KI = uhb /uIb ;

KII = uha /uIIa ;

KIII = uhc /uIcII .

(28)

uij

Herein, uhi describes the approximated CODs in direction of i and the expected CODs of mode j in direction of i. The accuracy of the determined SIFs based on the proposed method has already been demonstrated for stress-free and uniform loaded crack surfaces by the authors in [105]. Recently developed improvements or extensions of the standard XFEM such as, e.g., adaptive mesh refinements [18–20] or applied consecutive interpolation procedures [28–30] may further increase the accuracy of the computed SIFs, however, they are not the focus of this contribution and could be considered in later works. The evaluated SIFs are used in Section 4 to assess the crack configuration and determine whether a propagation occurs and in which direction. 4.2. External and internal loading cases Two different loadings are distinguished here. The first is called ’external’ and contains all stresses which come from the structure, such as body forces, prescribed tractions at the boundary, prescribed displacements etc. The second type of forces are stresses on the crack surface, typically pressure distributions, and are called ‘internal’. Fig. 9(a) shows a cracked domain where both loading types are present. This distinction is required as an increase of the internal load does not influence the external loads. f p In the following, Gm or Gm represents the energy release rate due to the external loads f or the internal load p of mode m. Based on the superposition principle of LEFM the total energy release rate Gtotal is given by a combination [106] of both loading types



Gtotal =

GIf +

 2 GIp



+

GIIf +

 2 GIIp



+

GIfII +



2

GIpII

.

(29)

As Eq. (29) shows, a modification of the pressure leads to a non-linear behaviour of the energy release rate. In the following, it is shown how to evaluate the critical pressure magnitude of a pre-defined pressure distribution based on Eq. (29). 4.3. Evaluation of critical pressure values In this section, the evaluation of a critical pressure magnitude is discussed based on the superposition principle of LEFM. Although, the behaviour of a crack is strongly dependent on the crack geometry, e.g., the crack length, for a given crack configuration in a linear elastic material with a neglected plastic zone, there holds for an initial prescribed pressure p0 that their associated SIFs scale proportionally with

Kmp (λ · p0 ) = λ · Kmp ( p0 ).

(30)

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

875

That is, the critical pressure pc for a purely internally loaded crack configuration can be computed by a combination of Eqs. (23), (24) and (30) as



pc = λ · p0

with

λ=

G



Gc Kmp

( p0 )

.

(31)

Eq. (31) only holds when the initial pressure p0 of an assumed pressure distribution p(p0 , t) has no effect on the distribution which means, a linear behaviour of p0 is demanded. When also external loads are present, a separate observation of the crack behaviour caused by the external and internal loadings is necessary. Herein, a simple access is possible by the superposition principle of LEFM where the crack configuration is decoupled into a part with only internal loadings and a part with only external loadings as shown in Fig. 9. The f p total for the combined case. For a associated SIFs are denoted by Km for the external case, Km for the internal case and Km total λ-scaled pressure state the total SIFs Km are given by the sum of the external and the λ-scaled internal case. total Km = Kmf + λKmp .

(32)

By using the crack propagation criterion of Eq. (23), a quadratic equation for the evaluation of the scaling factor λ is obtained, when Eq. (32) is inserted into Eqs. (24) and (29), hence,

0 = aλ2 + bλ + c

with

a = G(Kmp ), b=

 1 2 (1 − ν 2 ) f p KI KI + KIIf KIIp + KIfII KIpII , E μ

c = G(Kmf ) − Gc . p G(Km )

(33)

f G(Km )

Herein, and are the energy release rates due to the internal and external loads, and Gc is the critical energy release rate. This scaling factor λ enables the evaluation of the critical pressure pc based on Eq. (31). The critical pressure fulfills the propagation criterion of Eq. (23) and the crack geometry can be explicitly updated based on Eqs. (25) and (26). The major steps of the proposed procedure are summarized in Fig. 10. 5. Numerical examples This section presents test cases with loaded crack surfaces in two and three dimensions based on LEFM to verify and illustrate the performance of the proposed method. Achieved results are compared with analytical or empirical results when available. For all examples a brittle, isotropic and linearly elastic material is assumed with a Young’s modulus E = 37.8 GPa, a Poisson’s ratio ν = 0.25 and a critical energy release rate Gc = 680 J/m2 . Furthermore, plane stress conditions are assumed in two dimensions. It is noted that our interest in this work is only in the propagation of an arbitrarily, already existing mixedmode loaded macro-fracture within an isotropic material with prescribed pressure distributions. No crack initiations, microfracture fields, dynamic effects, are considered. Instead, the obtained final crack paths or surfaces and their corresponding critical pressure magnitudes are presented. For the examples related to HF, a pseudo-time is used for the propagation which is obtained based on the assumption of an impermeable material without fluid lag. Herein, the injected fluid is equal to the crack volume, wherefore, the time t can be expressed based on a surface integral over the crack surface with



t=

ub dA . Q0

(34)

There, ub are the CODs in direction of the local coordinate b of Section 2.2 and Q0 is the injection rate of the fluid. For all examples an injection rate of Q0 = 1 ml/s is assumed. Again, the HF examples are rather academic with focus on the presence of possibly complex loadings on the crack surfaces. The aim is to show the versatility and flexibility of the proposed numerical treatment which is also suited for other applications with loaded crack surfaces. The first two-dimensional test cases verify the estimation of the critical pressure with the procedure presented in Section 4.3. Herein, static planar crack configurations with well known analytical solutions of different pressure distributions are used. Example 5.4.1 shows the same procedure for a constant pressure distribution in three dimensions. All other examples investigate the propagation of a fracture due to the critical pressure. All systems of equations are solved with the Parallel Sparse Direct Solver PARDISO [107–109]. Domains are meshed with bilinear 4-node quadrilateral or trilinear 8-node hexahedral elements. 5.1. Pressurized crack in plate with finite width The first example investigates the evaluation of SIFs due to different pressure distributions within an edge cracked plate of finite width. The plate is 2m wide and 4m high and exhibits a 1m long initial crack, see Fig. 11(a). The performance of the proposed method is examined based on different mesh refinements, see Fig. 11(b) for an example mesh. Dirichlet boundary conditions are prescribed at the top and bottom nodes of the right side. Only in this example, a geometrical enrichment

876

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

Fig. 10. Algorithm of the proposed procedure.

is used where all element nodes of elements are crack tip enriched where the minimum distance of at least one element node to the crack tip is smaller than rtip = 10 cm. The domain is loaded by five different pressure distributions where the first three represent the simplest cases by polynomial functions, the fourth is a more general sinus-distribution and the last one the most complex case where a singular pressure at the crack tip is applied. A graphical illustration of the different distributions and their corresponding crack openings are shown in Fig. 12. Herein, the red lines represent the openings which are obtained from the configuration shown in Fig. 11(a). The green lines are the openings of the same configuration, however, with fixed displacements on the top and bottom of the plate. Analytical SIFs are given for all cases in [110,111] with

 KI =

0

a

σ (x )m(x, a )dx.

(35)

Herein, a is the crack length, σ (x) the pressure applied within the fracture, and m(x, a) the weight function. In our example, an initial pressure magnitude p0 of 1MPa is assumed. Fig. 13 shows the achieved results where the relative errors of the computed SIFs are plotted over the used number of elements. For this test case, the proposed method converges for finite pressure distributions with a rate of 1 to the analytical solution which is in good agreement with the achieved results

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

877

2 a=1 σ(x) x 2

2 (a)

(b)

Fig. 11. Test case with one edge crack in a finite width plate: (a) dimensions, supports, and loading; and (b) mesh.

Fig. 12. Applied pressure distributions and their corresponding crack openings.

in [28]. The singular pressure distribution leads to a lower convergence rate. This can be explained by the fact that singular pressure values are applied at a location where also singular stress and strain fields are present. 5.2. Mixed-mode pressurized crack under compression In this example, the accuracy of the computed critical pressure magnitudes is investigated based on a 3m long pressurized crack in a domain subjected to far field compressive bi-axial tractions σx = 10 MPa and σy = 7 MPa. The crack is oriented at an angle α with respect to the direction of horizontal traction σ x . A structured background-mesh consisting of bilinear quadrilateral elements is used to describe a 10 m · 10 m domain, see Fig. 14. Critical pressure magnitudes according to a constant prescribed pressure distribution are evaluated for different orientations of the initial crack based on different mesh refinements. The computed magnitudes are compared with an analytical solution which is based on Eq. (24). Herein, the resulting SIFs of a two-dimensional infinite domain are given by Rice [21,112]

KI =

√

 

π a p − σx sin2 α + σy cos2 α

KII = 0.5

√

π a[σx − σy ] sin 2α ,

(36)

where a is the half crack length, α is the orientation of the crack and p is the internal pressure. Achieved results are shown in Fig. 15 where the ratio of the computed and analytical critical pressure magnitude is plotted over the used number of elements.

878

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

Fig. 13. Relative error of KI due to different pressures within a finite width plate.

σy = 7MPa

σx = 10MPa

2a

=

3m

α

10m

σx = 10MPa

10m

σy = 7MPa

(a)

(b)

Fig. 14. Mixed-mode pressurized crack under compression: (a) dimensions and loading, and (b) mesh.

Numerical results are in good agreement with those available in the literature, where most of the meshes achieve pressure magnitudes with an error of less than 3%. In addition, these solutions also show a mesh independency as the same background mesh is used for all orientations of the crack and no significant deviations are observed.

5.3. Propagation of a mixed-mode pressurized crack under compression 5.3.1. Constant far field tractions This example investigates the influence of the external loads on the resulting crack path in a HF-related context with injection. A similar setting is investigated in the works of Dong [113] and Gupta [21]. In this contribution, a 3m long crack is pressurized with a constant pressure within the fracture and is subjected to far field compressive stresses σ x and σ y . The crack is oriented at an angle α = 40◦ with respect to the horizontal traction σ x . The domain is described with 59 × 59 bilinear quadrilateral elements. Variable ratios of the vertical and horizontal tractions are used where the ratio is defined by

η=

σy . σx

(37)

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

879

Fig. 15. Mixed-mode pressurized crack under compression.

4 3.5 3 2.5 2 1.5 1 0.5 0 0

(a)

1000

2000

3000

4000

5000

(b)

Fig. 16. Rotated pressurized crack loaded by constant far field tractions. (a) Resulting crack paths and (b) corresponding pressure magnitudes.

The resulting crack paths are shown in Fig. 16 with their corresponding pressures. Comparable pressure magnitudes are achieved as the resulting external load

R = σx · l y + σy · l x

(38)

is kept constant at a value of 40MN. The dimensions of the domain are lx = 10 m and ly = 10 m. It can be observed that the crack propagates in direction of the maximum load but it should be noted that this is just the case when the magnitudes of external loads and the required pressure are similar. When the internal pressure dominates the propagation, quasi straight crack extensions are expected as internal loads mainly lead to a mode I behaviour. An almost straight propagation is achieved where the horizontal and vertical tractions are equal as this behaviour presents a hydrostatic pressure. 5.3.2. Variable far field tractions The example of Section 5.3.1 is also used for linearly changed far field tractions, see Fig. 17(a). This far field behaviour leads to different loaded crack tips as the internal pressure is still prescribed as constant over the whole crack. For comparable pressure magnitudes the resulting load

R = 0.5(σx,o + σx,u ) · ly + σy · lx

(39)

880

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

9 8 7 6 5 4 3 2 1 0 0

1000

(a)

2000

3000

4000

5000

(b)

Fig. 17. Rotated pressurized crack loaded by linear far field tractions. (a) Resulting crack paths and (b) corresponding pressure magnitudes.

Fig. 18. Pressurized penny-shaped crack.

Fig. 19. Critical pressure for a penny-shaped crack.

is kept constant by a value of 120 MN. Different gradients of the horizontal tractions are investigated and characterized with

η=

σx,u σx,o

(40)

which represents the ratio of the upper and lower pressure magnitude of the horizontal tractions. The resulting crack paths and their corresponding pressure magnitudes are illustrated in Fig. 17. It can be observed that in this case the internal pressure dominates the propagation as the crack extension is almost straight. However, the propagation of the upper crack tip slightly tilts to the vertical direction whereas the lower one slightly

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

881

tilts to the horizontal direction. This can be explained by the fact that the horizontal traction is at the lower crack tip higher than the vertical traction whereas the behaviour is opposite at the upper crack tip. It should be noted that our model does not (yet) consider contact of the crack faces so that the internal pressure should dominate the external loading throughout the simulation.

Fig. 20. Mixed-mode loaded penny-shaped crack. (a) Configuration and final crack surfaces η = (b) 0.0, (c) 0.5, (d) 1.0, (e) 1.5 and (f) 2.0.

882

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

5.4. Penny-shaped crack In three dimensions, numerical results of a penny-shaped crack are generated. The domain which surrounds the embedded fracture is discretized with 39 × 39 × 39 trilinear hexahedral elements. Herein, the dimension of each edge of the domain is 2m and the radius r of the initial crack is set to r = 30 cm. The fracture is described with 96 flat triangles in which 16 edges represent the crack front, see Fig. 18. A critical energy release rate of Gc = 68 J/m2 is assumed. Two different configurations are considered in the following. 5.4.1. Pressurized penny-shaped crack within an infinite domain The first three-dimensional example, is motivated by a horizontal penny-shaped crack within an infinite domain which is pressurized with a constant fluid pressure p. In this example, we observe the development of the pressure during time and compare our results with a reference solution which is given, for example, in [21,114]. Herein, the critical crack radius rc at a certain time t is given as



rc =

9Q02 t 2 E  64π Gc

1 / 5

,

(41)

where Q0 is the fluid injection rate, E the plane-strain elastic modulus and Gc the critical energy release rate. The corresponding critical pressure is given with

 pc =

π Gc E  4rc

1 / 2 .

(42)

For the simulation of the propagation 15 time-steps are considered where it is assumed that the crack propagates with a maximum increment length of r = 1cm. Results are shown in Fig. 19 for a quasi static propagation simulation. In the reference solution we can see when the fluid is injected into the fracture the pressure begins to rise linearly until a critical value is reached. During this time the fluid fills and opens the fracture. After a critical value, a propagation starts and the pressure decreases which is in good agreement with our numerical results. 5.4.2. Mixed-mode loaded penny shaped crack The last example is an extension of the two-dimensional example presented in Section 5.3.1 where we want to show the influences of the external loads on the resulting crack surfaces of an initial penny-shaped crack. The discretization of the domain and the penny-shaped crack is the same as in Section 5.4.1. However, here the fracture is rotated by 40◦ about the x-axis. Variable ratios of the vertical and horizontal tractions are considered while the horizontal traction σ x is kept equal to σ y . The ratio of the vertical to the horizontal traction is defined by

σz σx

or

η=

σz . σy

(43)

3.5

3

2.5

critical pressure

η=

2

1.5

η = 0.0 η = 0.5 η = 1.0 η = 1.5 η = 2.0

1

0.5

0 0

20

40

60

time

80

100

Fig. 21. Critical pressure for a rotated penny-shaped crack.

120

140

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

883

In this test case, 10 time-steps are considered where the maximum crack increment r is set to 5cm. Similar to Eq. 38 the resulting force is kept constant at a value of 12MN. The resulting crack surfaces are shown in Fig. 20(b–f) for different ratios of the tractions. In Fig. 21, the change of the pressure is shown over time. The final crack surfaces and pressure distributions are in good agreement with those of the two-dimensional case which are shown and discussed in Section 5.3.1 and are also in good agreement with those obtained in [21,115]. 6. Conclusion A numerical framework is presented which considers crack propagation induced by loaded crack surfaces which is relevant in a number of applications. Among the most important ones is HF wherefore we link to that context frequently and discuss concrete examples how to generate complex and meaningful pressure fields. The XFEM with a hybrid explicitimplicit crack description is used for the accurate approximation of displacement quantities in the cracked domain. The method is significantly improved for a more accurate representation of the crack front and the associated coordinate systems. A major focus is on the data transfer between the explicit and implicit description and the proper consideration in the XFEM to consider for arbitrary load distributions. Therefore, integration points must be placed properly on the zero-level set of the implicit description. Then, closest point projections are employed to extract data from the explicit description. A crack propagation model is developed where internal and external loadings are distinguished. The computation of a load factor for the internal loading on the crack surface while keeping the external loading constant is outlined. Finally, the computation of SIFs based on CODs is described in the context of implicit crack descriptions in the XFEM. The resulting model and numerical approach is simple and versatile and numerical results show the success of the approach. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, CRC Press, 2005. M. Kanninen, C. Popelar, Advanced Fracture Mechanics, Oxford engineering science series, Oxford University Press, 1985. J. Knott, Fundamentals of Fracture Mechanics, Butterworths, 1973. T. Perkins, L. Kern, et al., Widths of hydraulic fractures, J. Petrol. Tech. 13 (09) (1961) 937–949. S. Khristianovic, Y. Zheltov, Formation of vertical fractures by means of highly viscous fluids, in: Proceedings of the 4th World Petroleum Congress, Rome, 2, 1955, pp. 579–586. A. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci. 44 (5) (2007) 739–757, doi:10. 1016/j.ijrmms.20 06.11.0 06. A. Rubin, Propagation of magma-filled cracks, Annu. Rev. Earth Planet. Sci. 23 (1995) 287–336, doi:10.1146/annurev.ea.23.050195.001443. C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Eng. 199 (45-48) (2010) 2765–2778, doi:10.1016/j.cma.2010.04.011. C. Kuhn, R. Müller, A continuum phase field model for fracture, Eng. Fract. Mech. 77 (18, SI) (2010) 3625–3634, doi:10.1016/j.engfracmech.2010.08.009. N. Moës, C. Stolz, P.E. Bernard, N. Chevaugeon, A level set based model for damage growth: The thick level set approach, Int. J. Numer. Methods Eng. 86 (3) (2011) 358–380, doi:10.1002/nme.3069. M. Ambati, T. Gerasimov, L. De Lorenzis, A review on phase-field models of brittle fracture and a new fast hybrid formulation, Comput. Mech. 55 (2) (2015) 383–405, doi:10.10 07/s0 0466- 014- 1109- y. A. Egger, U. Pillai, K. Agathos, E. Kakouris, E. Chatzi, I.A. Aschroft, S.P. Triantafyllou, Discrete and phase field methods for linear elastic fracture mechanics: a comparative study and state-of-the-art review, Appl. Sci. Basel 9 (12) (2019), doi:10.3390/app9122436. J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (5) (1993) 277–296. T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Methods Eng. 45 (5) (1999) 601–620, doi:10. 1002/(SICI)1097- 0207(19990620)45:5<601::AID- NME598>3.0.CO;2- S. N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Methods Eng. 46 (1) (1999) 131–150, doi:10.1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J. J. Dolbow, N. Moës, T. Belytschko, Discontinuous enrichment in finite elements with a partition of unity method, Finite Elem. Anal. Des. 36 (3) (20 0 0) 235–260. N. Sukumar, N. Moës, B. Moran, T. Belytschko, Extended finite element method for three-dimensional crack modelling, Int. J. Numer. Methods Eng. 48 (11) (20 0 0) 1549–1570, doi:10.1002/1097-0207(20000820)48:11<1549::AID- NME955>3.0.CO;2- A. T. Yu, T.Q. Bui, Numerical simulation of 2-D weak and strong discontinuities by a novel approach based on XFEM with local mesh refinement, Comput. Struct. 196 (2018) 112–133. Z. Wang, T. Yu, T.Q. Bui, S. Tanaka, C. Zhang, S. Hirose, J.L. Curiel-Sosa, 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks, Comput. Methods Appl. Mech. Eng. 313 (2017) 375–405. Y. Jin, O. González-Estrada, O. Pierard, S. Bordas, Error-controlled adaptive extended finite element method for 3d linear elastic crack propagation, Comput. Methods Appl. Mech. Eng. 318 (2017) 319–348. P. Gupta, C.A. Duarte, Simulation of non-planar three-dimensional hydraulic fracture propagation, Int. J. Numer. Anal. Methods Geomech. 38 (13) (2014) 1397–1430, doi:10.1002/nag.2305. T.-P. Fries, A. Byfut, A. Alizada, K.W. Cheng, A. Schröder, Hanging nodes and XFEM, Int. J. Numer. Methods Eng. 86 (4-5) (2011) 404–430. C. Daux, N. Moës, J. Dolbow, N. Sukumar, T. Belytschko, Arbitrary branched and intersecting cracks with the extended finite element method, Int. J. Numer. Methods Eng. 48 (12) (20 0 0) 1741–1760. M. Sheng, G. Li, S. Shah, A.R. Lamb, S.P. Bordas, Enriched finite elements for branching cracks in deformable porous media, Eng. Anal. Boundary Elem. 50 (2015) 435–446. F. Cruz, D. Roehl, E. do Amaral Vargas Jr, An XFEM element to model intersections between hydraulic and natural fractures in porous rocks, Int. J. Rock Mech. Min. Sci. 112 (2018) 385–397. F. Stazi, E. Budyn, J. Chessa, T. Belytschko, An extended finite element method with higher-order elements for curved cracks, Comput. Mech. 31 (1-2) (2003) 38–48. K.W. Cheng, T.-P. Fries, Higher-order XFEM for curved strong and weak discontinuities, Int. J. Numer. Methods Eng. 82 (5) (2010) 564–590. Z. Kang, T.Q. Bui, T. Saitoh, S. Hirose, et al., An extended consecutive-interpolation quadrilateral element (xcq4) applied to linear elastic fracture mechanics, Acta Mech. 226 (12) (2015) 3991–4015.

884

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

[29] Z. Kang, T.Q. Bui, T. Saitoh, S. Hirose, Quasi-static crack propagation simulation by an enhanced nodal gradient finite element with different enrichments, Theor. Appl. Fract. Mech. 87 (2017) 61–77. [30] Z. Kang, T.Q. Bui, S. Hirose, et al., Dynamic stationary crack analysis of isotropic solids and anisotropic composites by enhanced local enriched consecutive-interpolation elements, Compos. Struct. 180 (2017) 221–233. [31] N. Moës, T. Belytschko, Extended finite element method for cohesive crack growth, Eng. Fract. Mech. 69 (2002) 813–833, doi:10.1016/S0013-7944(01) 00128-X. [32] A.R. Khoei, M. Hirmand, M. Vahab, M. Bazargan, An enriched FEM technique for modeling hydraulically driven cohesive fracture propagation in impermeable media with frictional natural faults: Numerical and experimental investigations, Int. J. Numer. Methods Eng. 104 (6) (2015) 439–468, doi:10.1002/nme.4944. [33] J. Jas´ kowiec, F.P. van der Meer, A consistent iterative scheme for 2d and 3d cohesive crack analysis in XFEM, Comput. Struct. 136 (2014) 98–107. [34] X. Zhang, T.Q. Bui, A fictitious crack XFEM with two new solution algorithms for cohesive crack growth modeling in concrete structures, Eng. Comput. 32 (2) (2015) 473–497. [35] T.-P. Fries, T. Belytschko, The extended/generalized finite element method: An overview of the method and its applications, Int. J. Numer. Methods Eng. 84 (3) (2010) 253–304, doi:10.1002/nme.2914. [36] T. Belytschko, R. Gracie, G. Ventura, A review of extended/generalized finite element methods for material modeling, Modell. Simul. Mater. Sci. Eng. 17 (4) (2009), doi:10.1088/0965-0393/17/4/043001. [37] A.R. Khoei, Extended Finite Element Method: Theory and Applications, John Wiley & Sons, 2014. [38] S. Pommier, A. Gravouil, N. Moës, A. Combescure, Extended Finite Element Method for Crack Propagation, John Wiley & Sons, 2013. [39] B. Rao, S. Rahman, An efficient meshless method for fracture analysis of cracks, Comput. Mech. 26 (4) (20 0 0) 398–408. [40] X. Zhuang, C. Augarde, K. Mathisen, Fracture modeling using meshless methods and level sets in 3D: framework and modeling, Int. J. Numer. Methods Eng. 92 (11) (2012) 969–998. [41] T. Rabczuk, T. Belytschko, A three-dimensional large deformation meshfree method for arbitrary evolving cracks, Comput. Methods Appl. Mech. Eng. 196 (29-30) (2007) 2777–2799. [42] M. Duflot, A meshless method with enriched weight functions for three-dimensional crack propagation, Int. J. Numer. Methods Eng. 65 (12) (2006) 1970–2006. [43] G.E. Blandford, A.R. Ingraffea, J.A. Liggett, Two-dimensional stress intensity factor computations using the boundary element method, Int. J. Numer. Methods Eng. 17 (3) (1981) 387–404. [44] Y. Mi, M. Aliabadi, Dual boundary element method for three-dimensional fracture mechanics analysis, Eng. Anal. Boundary Elem. 10 (2) (1992) 161–171. [45] Z. Yang, Fully automatic modelling of mixed-mode crack propagation using scaled boundary finite element method, Eng. Fract. Mech. 73 (12) (2006) 1711–1731. [46] S. Goswami, W. Becker, Computation of 3-D stress singularities for multiple cracks and crack intersections by the scaled boundary finite element method, Int. J. Fract. 175 (1) (2012) 13–25. [47] C. Song, E.T. Ooi, S. Natarajan, A review of the scaled boundary finite element method for two-dimensional linear elastic fracture mechanics, Eng. Fract. Mech. 187 (2018) 45–73. [48] B.H. Nguyen, H.D. Tran, C. Anitescu, X. Zhuang, T. Rabczuk, An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems, Comput. Methods Appl. Mech. Eng. 306 (2016) 252–275, doi:10.1016/j.cma.2016.04.002. [49] X. Peng, E. Atroshchenko, P. Kerfriden, S.P.A. Bordas, Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth, Comput. Methods Appl. Mech. Eng. 316 (SI) (2017) 151–185, doi:10.1016/j.cma.2016.05.038. [50] S.G. Ferreira Cordeiro, E.D. Leonel, Mechanical modelling of three-dimensional cracked structural components using the isogeometric dual boundary element method, Appl. Math. Modell. 63 (2018) 415–444, doi:10.1016/j.apm.2018.06.042. [51] F.L. Sun, C.Y. Dong, H.S. Yang, Isogeometric boundary element method for crack propagation based on Bezier extraction of NURBS, Eng. Anal. Boundary Elem. 99 (2019) 76–88, doi:10.1016/j.enganabound.2018.11.010. [52] T.Q. Bui, Extended isogeometric dynamic and static fracture analysis for cracks in piezoelectric materials using NURBS, Comput. Methods Appl. Mech. Eng. 295 (2015) 470–509, doi:10.1016/j.cma.2015.07.005. [53] S.S. Ghorashi, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, Int. J. Numer. Methods Eng. 89 (9) (2012) 1069–1101, doi:10.1002/nme.3277. [54] S. Natarajan, J. Wang, C. Song, C. Birk, Isogeometric analysis enhanced by the scaled boundary finite element method, Comput. Methods Appl. Mech. Eng. 283 (2015) 733–762, doi:10.1016/j.cma.2014.09.003. [55] N. Nguyen-Thanh, J. Huang, K. Zhou, An isogeometric-meshfree coupling approach for analysis of cracks, Int. J. Numer. Methods Eng. 113 (10) (2018) 1630–1651, doi:10.1002/nme.5713. [56] T. Fries, M. Baydoun, Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description, Int. J. Numer. Methods Eng. 89 (12) (2012) 1527–1558, doi:10.1002/nme.3299. [57] M. Baydoun, T.P. Fries, Crack propagation criteria in three dimensions using the XFEM and an explicit-implicit crack description, Int. J. Fract. 178 (1-2, SI) (2012) 51–70, doi:10.1007/s10704- 012- 9762- 7. [58] S. Salimzadeh, N. Khalili, A three-phase XFEM model for hydraulic fracturing with cohesive crack propagation, Comput. Geotech. 69 (2015) 82–92, doi:10.1016/j.compgeo.2015.05.001. [59] T. Mohammadnejad, A.R. Khoei, An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model, Finite Elem. Anal. Des. 73 (2013) 77–95, doi:10.1016/j.finel.2013.05.005. [60] B. Lecampion, An extended finite element method for hydraulic fracture problems, Commun. Numer. Methods Eng. 25 (2) (2009) 121–133, doi:10. 1002/cnm.1111. [61] J. Réthoré, R.d. Borst, M.-A. Abellan, A two-scale approach for fluid flow in fractured porous media, Int. J. Numer. Methods Eng. 71 (7) (2007) 780–800, doi:10.1002/nme.1962. [62] J. Réthoré, R.d. Borst, M.-A. Abellan, A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks, Comput. Mech. 42 (2) (2008) 227–238, doi:10.1007/s00466-007- 0178- 6. [63] A.R. Khoei, M. Vahab, E. Haghighat, S. Moallemi, A mesh-independent finite element formulation for modeling crack growth in saturated porous media based on an enriched-FEM technique, Int. J. Fract. 188 (1) (2014) 79–108, doi:10.1007/s10704- 014- 9948- 2. [64] E. Detournay, Mechanics of Hydraulic Fractures, in: S.H Davis, P. Moin (Eds.), Annual Review of Fluid Mechanics, 48, 2016, pp. 311–339, doi:10.1146/ annurev- fluid- 010814- 014736. [65] E. Gordeliy, A. Peirce, Implicit level set schemes for modeling hydraulic fractures using the XFEM, Comput. Methods Appl. Mech. Eng. 266 (2013) 125–143, doi:10.1016/j.cma.2013.07.016. [66] P. Laborde, J. Pommier, Y. Renard, M. Salaün, High-order extended finite element method for cracked domains, Int. J. Numer. Methods Eng. 64 (3) (2005) 354–381. [67] E. Béchet, H. Minnebo, N. Moës, B. Burgardt, Improved implementation and robustness study of the X-FEM for stress analysis around cracks, Int. J. Numer. Methods Eng. 64 (8) (2005) 1033–1056. [68] A. Menk, S.P. Bordas, A robust preconditioning technique for the extended finite element method, Int. J. Numer. Methods Eng. 85 (13) (2011) 1609–1632. [69] S. Loehnert, A stabilization technique for the regularization of nearly singular extended finite elements, Comput. Mech. 54 (2) (2014) 523–533.

M. Schätzer and T.-P. Fries / Applied Mathematical Modelling 78 (2020) 863–885

885

[70] V. Gupta, C.A. Duarte, I. Babuska, U. Banerjee, Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics, Comput. Methods Appl. Mech. Eng. 289 (2015) 355–386, doi:10.1016/j.cma.2015.01.014. [71] K. Agathos, G. Ventura, E. Chatzi, S.P. Bordas, Stable 3D XFEM/vector level sets for non-planar 3D crack propagation and comparison of enrichment schemes, Int. J. Numer. Methods Eng. 113 (2) (2018) 252–276. [72] C.A. Duarte, D.-J. Kim, Analysis and applications of a generalized finite element method with global–local enrichment functions, Comput. Methods Appl. Mech. Eng. 197 (6-8) (2008) 487–504. [73] M. Malekan, F.B. Barros, Well-conditioning global–local analysis using stable generalized/extended finite element method for linear elastic fracture mechanics, Comput. Mech. 58 (5) (2016) 819–831. [74] R. Tian, L. Wen, L. Wang, Three-dimensional improved XFEM (XFEM) for static crack problems, Comput. Methods Appl. Mech. Eng. 343 (2019) 339–367. [75] T. Belytschko, N. Moës, S. Usui, C. Parimi, Arbitrary discontinuities in finite elements, Int. J. Numer. Methods Eng. 50 (4) (2001) 993–1013. [76] T. Fries, A corrected XFEM approximation without problems in blending elements, Int. J. Numer. Methods Eng. 75 (5) (2008) 503–532, doi:10.1002/ nme.2259. [77] I. Babuška, U. Banerjee, Stable Generalized Finite Element Method (SGFEM), Comput. Methods Appl. Mech. Eng. 201 (2012) 91–111, doi:10.1016/j.cma. 2011.09.012. [78] N. Gravouil A. and, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets - Part II: Level set update, Int. J. Numer. Methods Eng. 53 (11) (2002) 2569–2586, doi:10.1002/nme.430. [79] D.L. Chopp, N. Sukumar, Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method, Int. J. Eng. Sci. 41 (8) (2003) 845–869. [80] N. Sukumar, D.L. Chopp, B. Moran, Extended finite element method and fast marching method for three-dimensional fatigue crack propagation, Eng. Fract. Mech. 70 (1) (2003) 29–48. [81] M. Duflot, A study of the representation of cracks with level sets, Int. J. Numer. Methods Eng. 70 (11) (2007) 1261–1302, doi:10.1002/nme.1915. [82] N. Moës, A. Gravouil, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets - Part I: Mechanical model, Int. J. Numer. Methods Eng. 53 (11) (2002) 2549–2568, doi:10.1002/nme.429. [83] T.-P. Fries, S. Omerovic´ , Higher-order accurate integration of implicit geometries, Int. J. Numer. Methods Eng. 106 (5) (2016) 323–371. ´ D. Schöllhammer, J. Steidl, Higher-order meshing of implicit geometries-Part I: Integration and interpolation in cut elements, [84] T. Fries, S. Omerovic, Comput. Methods Appl. Mech. Eng. 313 (2017) 759–784, doi:10.1016/j.cma.2016.10.019. [85] A. de Boer, A.H. van Zuijlen, H. Bijl, Review of coupling methods for non-matching meshes, Comput. Methods Appl. Mech. Eng. 196 (8) (2007) 1515– 1525, doi:10.1016/j.cma.2006.03.017. [86] P. Thevenaz, T. Blu, M. Unser, Interpolation revisited, IEEE Trans. Med. Imaging 19 (7) (20 0 0) 739–758, doi:10.1109/42.875199. [87] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 20 0 0. [88] E. Detournay, Propagation regimes of fluid-driven fractures in impermeable rocks, Int. J. Geomech. 4 (1) (2004) 35–45. [89] J. Adachi, E. Detournay, Self-similar solution of a plane-strain fracture driven by a power-law fluid, Int. J. Numer. Anal. Methods Geomech. 26 (6) (2002) 579–604. [90] A. Savitski, E. Detournay, Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions, Int. J. Solids Struct. 39 (26) (2002) 6311–6337. [91] R. Ritchie, Mechanisms of fatigue-crack propagation in ductile and brittle solids, Int. J. Fract. 100 (1) (1999) 55–83, doi:10.1023/A:1018655917051. [92] J.R. Rice, A path independent integral and the approximate analysis of strain concentration by notches and cracks, J. Appl. Mech. 35 (2) (1968) 379–386. [93] J.F. Yau, S.S. Wang, H.T. Corten, A Mixed-Mode Crack Analysis of Isotropic Solids Using Conservation Laws of Elasticity, J. Appl. Mech. 47 (2) (1980) 335–341. [94] H. Awaji, T. Kato, S. Honda, T. Nishikawa, Criterion for combined mode I-II of brittle fracture, J. Ceram. Soc. Jpn. 107 (10) (1999) 918–924. [95] W. Weber, G. Kuhn, An optimized predictor-corrector scheme for fast 3D crack growth simulations, Eng. Fract. Mech. 75 (3-4) (2008) 452–460, doi:10.1016/j.engfracmech.2007.01.005. [96] F. Erdogan, G. Sih, On the crack extension in plates under plane loading and transverse shear, J. Basic. Eng. 85 (4) (1963) 519–525. [97] S.K. Maiti, R. Smith, Comparison of the criteria for mixed mode brittle fracture based on the preinstability stress-strain field, Int. J. Fract. 24 (1) (1984) 5–22. [98] M. Hussain, S. Pu, J. Underwood, Strain energy release rate for a crack under combined mode i and mode ii, in: Proceedings of the Fracture Analysis: National Symposium on Fracture Mechanics, Part II, ASTM International, 1974. [99] G. Sih, B. Macdonald, Fracture mechanics applied to engineering problems-strain energy density fracture criterion, Eng. Fract. Mech. 6 (2) (1974) 361–386. [100] J. Dolbow, M. Gosz, On the computation of mixed-mode stress intensity factors in functionally graded materials, Int. J. Solids Struct. 39 (9) (2002) 2557–2574, doi:10.1016/S0 020-7683(02)0 0114-2. [101] M. Walters, G. Paulino, R. Dodds, Interaction integral procedures for 3-D curved cracks including surface tractions, Eng. Fract. Mech. 72 (11) (2005) 1635–1663, doi:10.1016/j.engfracmech.2005.01.002. [102] S.K. Chan, I.S. Tuba, W.K. Wilson, On the finite element method in linear fracture mechanics, Eng. Fract. Mech. 2 (1970) 1–17. [103] M. Nejati, A. Paluszny, R. Zimmerman, On the use of quarter-point tetrahedral finite elements in linear elastic fracture mechanics, Eng. Fract. Mech. 144 (2015) 194–221, doi:10.1016/j.engfracmech.2015.06.055. [104] N. Muthu, S.K. Maiti, B.G. Falzon, I. Guiamatsia, A comparison of stress intensity factors obtained through crack closure integral and other approaches using eXtended element-free Galerkin method, Comput. Mech. 52 (3) (2013) 587–605, doi:10.10 07/s0 0466- 013- 0834- y. [105] M. Schätzer, T.-P. Fries, Stress intensity factors through crack opening displacements in the XFEM, in: Advances in Discretization Methods, Springer, 2016, pp. 143–164. [106] C. Rans, R. Alderliesten, Formulating an effective strain energy release rate for a linear elastic fracture mechanics description of delamination growth, in: Proceedings of 17th International Conference on Composite Materials (ICCM-17), 2009. [107] O. Schenk, K. Gärtner, Solving unsymmetric sparse systems of linear equations with pardiso, Future Gener. Comput. Syst. 20 (3) (2004) 475–487. [108] O. Schenk, K. Gärtner, On fast factorization pivoting methods for sparse symmetric indefinite systems, Electron. Trans. Numer. Anal. 23 (1) (2006) 158–179. [109] G. Karypis, V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput. 20 (1) (1998) 359–392. Ref: Pardiso [110] G. Glinka, G. Shen, Universal features of weight functions for cracks in mode i, Eng. Fract. Mech. 40 (6) (1991) 1135–1146. [111] A. Kaya, F. Erdogan, Stress intensity factors and cod in an orthotropic strip, Int. J. Fract. 16 (2) (1980) 171–190. [112] J.R. Rice, Mathematical analysis in the mechanics of fracture, Fract. Adv. Treat. 2 (1968) 191–311. [113] C. Dong, C. de Pater, Numerical implementation of displacement discontinuity method and its application in hydraulic fracturing, Comput. Methods Appl. Mech. Eng. 191 (8-10) (2001) 745–760, doi:10.1016/S0045-7825(01)00273-0. [114] B. Bourdin, C.P. Chukwudozie, K. Yoshioka, et al., A variational approach to the numerical simulation of hydraulic fracturing, in: Proceedings of the SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2012. [115] A. Sadeghirad, D.L. Chopp, X. Ren, E. Fang, J. Lua, A novel hybrid approach for level set characterization and tracking of non-planar 3d cracks in the extended finite element method, Eng. Fract. Mech. 160 (2016) 1–14.