A phase-field approach to pneumatic fracture with anisotropic crack resistance

Phase-field models of fracture allow the prediction of crack propagation and crack patterns. In this contribution, externally driven fracture processes in linear and finite elasticity are investigated. Different approaches to consider pneumatic pressure and materials with non-isotropic crack resistance are studied, combined, and examined in detail. The versatility of the proposed models is proven by a series of numerical simulations in two and three dimensions.

among many others. Starting from the original phasefield fracture model of brittle elastic solids, the method has been extended to ductile fracture (Ambati and De Lorenzis 2016;Kuhn et al. 2016), to different coupled field problems like temperature or electric fields (Abdollahi and Arias 2012;Miehe et al. 2015;Zhou et al. 2018), to heterogenous solids (Teichtmeister et al. 2017;Bleyer and Alessi 2018;Li and Maurini 2019), to hydraulic fracture (Bourdin et al. 2012;Heider and Markert 2017;Wilson and Landis 2016;Zhou et al. 2018;Mikelić et al. 2019) and so on. In this respect, the present work contributes to the modeling of driven crack growth in solids with anisotropic behavior-a constellation which occurs in geological formations, in toughened glass or in strained biological tissues and often comes together with large deformations.
Various studies have been made on external load driven fracture, in particular, on hydraulic fracturing where a fluid under high pressure is injected into compressed soil to open cracks. Several phase-field approaches to hydraulic fracturing have been suggested, e.g. by Miehe et al. (2015Miehe et al. ( , 2016, Mikelic et al. (2015Mikelic et al. ( , 2015 and others Luo and Ehlers (2015), Heider and Markert (2017) and Wilson and Landis (2016), usually assuming a homogeneous or porous bulk material. Pneumatic fracture can be interpreted as a simplification of hydraulic fracture in the sense that only the crack opening pressure is relevant and the coupling with fluid flow is negligible. Related loading regimes can be found in subsurface volcanic activities or pressurized vessel and comparable situations may also arise in the context of bio-materials. For example, a high blood pressure vasculature may lead to the rupture of vascular walls evoking major complications, cf. Yamashima and Friede (1984) and Ghilardi et al. (1993), Fig. 1.
By now most phase-field fracture simulations refer to homogeneous isotropic material. For the mentioned driven fracture situations this is not a realistic assumption. For instance, the typical shale formations of hydraulic fracturing have a heterogeneously layered structure. Such layers may have interfaces and directions with reduced crack resistance but also the stone itself may have inhomogeneities or be direction dependent. There are different ways to consider such fracture anisotropy. An explicit modeling of interfaces and its interaction with crack propagation is presented in Hansen-Dörr et al. (2017) and Bilgen et al. (2018) in the phase-field fracture context. Ghamgosar et. al. (2015) suggest to introduce various critical fracture energy densities G I c , G I I c , G I I I c for each rock direction. This technique was used in the early phase-field models of Adda-Bedia et al. (1999) to capture the preferred direction of crack propagation. Such adapted solution techniques facilitate the description of crack propagation in specific situations like pure mode I or mode II loading. In the last years it has been generalized to a Griffith energy tensor or, alternatively, a structural tensor included in the regularized variational phase-field formulation, cf. Clayton and Knap (2015), Teichtmeister et al. (2017), Liu and Juhre (2018) and Li et al. (2015). Such a tensor weights the partial derivatives in different directions separately.
In this work we present a model for pneumatic fracture in solids with non-isotropic fracture resistance. The model is novel, allows for linearized and finite deformations, and is based on a strict variational derivation. The direction dependent fracture resistance is mapped by a stochastic and by a deterministic approach. The paper is structured as follows: In Sect. 2 the phase-field fracture approach is introduced shortly. The equations of pneumatic fracture are specified in Sect. 3 where the crack opening pressure is incorporated as a moving boundary condition. Next, in Sect. 4 the non-isotropic fracture resistance is realized by a random field, a nonisotropic crack density function and a combination of both. All models are evaluated by series of numerical examples. Numerical simulations for pressure induced crack growth in materials with layer-wise fracture resistance are presented in Sect. 5. A short summary in Sect. 6 concludes the paper.

The phase-field approach to fracture
In phase-field fracture the material state is characterized by a continuous field z(x, t) ∈ [0, 1] -the phasefield. The limit z = 1 indicates the crack and z = 0 marks the initial, unbroken state.
In the following we consider a solid with domain Ω ⊂ R 3 and boundary ∂Ω = Γ ⊂ R 2 deforming in a time interval t ∈ [0, t tot ] ⊂ R + . Its total potential energy shall be given by an elastic bulk energy with a Helmholtz free energy density Ψ e and additional surface contributions due to cracking: The surface integral in (1) results from the fact that crack growth involves the creation of new surfaces Γ (t); the specific energy G c corresponds to Griffith's critical energy release rate in brittle fracture. In phasefield fracture the sharp crack surface is approximated by a crack density function γ = γ (z, ∇z). In the present work we apply the second-order function where the length parameter l c weights the influence of the gradient term and determines the width of the diffuse zone. This form of the crack density function has commonly been used in the past years, cf. Bourdin et al. (2008), Miehe et al. (2010), Kuhn and Müller (2010) and Borden et al. (2012). The regularized total potential energy now has the form of an Ambrosio-Tortorelli functional Ambrosio and Tortorelli (1990) and states a non-local fracture theory here. The phase field z follows from potential energy optimization but is commonly subjected to a viscous regularization. We write here in general forṁ with a mobility parameter M [m 2 /Ns] and a local force symbolized by Y . The latter is usually the variation of There are different ways to formulate the crack driving force Y e . In most theoretical works it is simply Y e = δ z Ψ e . Practical prerequisite, though, is a decomposition of the energy density which takes the asymmetry of fracture into account, i. e. only tensile states lead to crack propagation, Y e = δ z Ψ e+ . However, there is no unique way to formulate such a decomposition and other physically motivated crack driving forces may be applied; we refer to Bilgen and Weinberg (2019) for a discussion. As an example, we consider here the local tensile stress state (Rankine stress) to be responsible for crack growth, i. e. the crack driving force is active when the maximum principal stress σ I = max(σ a ), a ∈ {1, 2, 3} exceeds a critical value σ c :

Finite elasticity
The phase-field model is adapted to finite deformations with a deformation mapping χ(X, t) : Ω × [0, t tot ] → R 3 and the deformation gradient where the quantities in capitals refer to the initial configuration. The determinant J = det F represents the mapping of volume elements, J = dv/ dV . The elastic strain energy density is defined as isotropic tensor function of F and can be formulated with its invariants. Specifically, we set here with modified invariants based on the isochoric split of the deformation gradient, F = J 1/3F .
For later reference, we also define the right Cauchy-Green deformation tensor C = F T F and the cofactor tensor H = cof(F) which represents the mapping of surface elements in a deformation. With the latter, the second invariant can be formulated as As proposed in Hesch et al. (2017) and analyzed in Thomas et al. (2020) we consider the tensioncompression asymmetry of fracture by a decomposition of the invariants into tensile (+) and compressive (−) components This allows to formulate the elastic strain energy density in a decomposed way, The degradation function g(z) is introduced to account for the material loss; we set here g(z) = (1 − z) 2 . In the following we apply a Neohookean material model with strain energy density where K is the bulk modulus and μ the initial shear modulus. Naturally, the simple volumetric term in (12) can easily be replaced by a convex expression, which then renders the energy density to be polyconvex. For detailed information we refer to Hesch et al. (2017) and Weinberg and Hesch (2015).
The decomposition of the invariants into tensile and compressive components gives Then, the volumetric part of the energy Ψ vol = 1 2 K (J − 1) 2 corresponds to Ψ vol = 1 2 K (tr ) 2 in linear elasticity. Accordingly, the decomposition of the energy conforms with the split of the linear elastic energy proposed by Amor et al. (2009) and we conclude that the proposed decomposition of the strain energy density (11) can be adapted to small deformations as well.
For later reference we state here another common energy split, namely the λ − μ split, with Lámeconstants λ = K − 2 3 μ and μ, 3 Pneumatic fracture In this section we consider pneumatic fracture, i.e., crack growth in a compressed medium driven by a pressurized state within the crack. The applied pressure is modeled as a moving boundary condition. The model is similar to the approach of Bourdin et al. (2012), who considered small deformations only, and derived here in a variational form for small and finite deformations.

Pressure driven crack growth
We start with the potential energy (1) of the body in its regularized form (3) and focus on the strain energy density (12), its small strain counterpart or the corresponding formulation (16). The pressurized gas injected for pneumatic cracking is associated with another energy contribution and so we add a pneumatic term which has the form or, for linearized strains, Herebyp(z; t) is the applied pressure. This pressure field is coupled with the phase field z and may vary in time. Within a crack at position x it is modeled as: The pressure is active when the phase-field is = 0 and acts on the crack then. At the crack flanks the static boundary condition σ · n = 0 typically holds. For a crack with pressurep it changes to σ · n =p or σ · n −p = 0. However, in the diffuse crack zone of the phase-field model the normal to a crack flank n cannot be computed directly and so the additional pressure component has to be added in the variational approach.
The specific pressure field (20) is a matter of modeling; the linear form corresponds to the model of Bourdin (2012) for small strains, see Remark 3.6. The first Piola-Kirchhoff stress tensor is derived from the strain energy density as and relates to the Cauchy stress tensor σ by σ = J −1 P F T , cf. Holzapfel (2006). With (17) and (18) it follows where we denote the Neohookean part of the stresses by P e and H = cof(F) as above.
Alternatively we derive the linear-elastic stresses as the work conjugate to the linearized strains, which gives with σ e = ∂Ψ e /∂ the modified stresses Please note that this is a (diffuse) moving boundary condition since the loaded boundary grows with crack propagation.

Proof of concept
To demonstrate applicability we study a two dimensional notched plate of size 100 × 100 mm 2 . The mesh consists of 256×256 quadratic B-spline finite elements.
Here we use the spline functions for convenience only; they provide here a similar approximation as classical Lagrangian elements of quadratic order. While the boundaries are constrained in all directions, a notch of length 16 mm in the middle of the domain is prescribed by fixing the phase-field parameter to z = 1, see left plot of Fig. 2. The numerical simulation uses the Neohookean model with the decomposed invariants (10). The material parameters correspond to concrete and are given in Table 1; the mobility of Eq. (4) is M = 1mm 2 /Ns whereby the time is normalized, t = t n · pΔt, Δt = 1s. The nonlocal length parameter is l c = 4.5 · 10 −4 mm. For the applied pressure the values p 0 = 10 N/mm 2 and p 1 = 10 4 · t n N/mm 2 are chosen, with t n ∈ {1, . . . , 75}.
The growing crack is shown in Fig. 2  We remark that the irreversibility of crack growth in (17) is enforced through an inequality constraintṡ ≤ 0. In Wheeler et al. (2014) an augmented Lagrangian algorithm for energy optimization is proposed instead.

Validation
In order to validate our solution we make use of analytical expressions for the two-dimensional Griffith crack, which go back to Sneddon & Lowengrub (1969), Sects. 2.1 and 2.4. The critical energy release rate of a pressurized crack of length 2a is G c = π p 2 a/E with the elastic parameter E = E/(1 − ν 2 ). Consequently, the critical pressure of crack opening follows as Accordingly, the crack opening displacement can be integrated to give the 2D volume of the crack, V =  (27) 2π pa 2 /E , which has a critical value at pressure p c , i.e.
Injection of a volume less then V c gives no crack growth but a linear pressure increase in the elastically deforming body. After V c is reached, the crack starts propagating. Eliminating a from p c and V c , we obtain the corresponding pressure-volume relation in closed form: In Fig. 3 the result of our numerical simulation is compared with the analytical solution (27). Hereby the pressure is prescribed and the resulting 2D volume is calculated from the displacement of the crack flanks. We see the linear pressure increase in the body until the critical pressure is reached; then the crack grows and the pressure decays. Both curves have the same shape but the pressure during crack growth is smaller in our simulation than in the analytical solution. We remark that we prescribe the crack pressure but do not control the volume. Here and below, we define a crack boundary at phase-field values of z ≥ 0.9; other values may result in different volumes. The same holds for further phasefield model details like the degeneration function g(z) and length parameter l c . For related studies we refer to Bourdin et al. (2012) and Mikelic et al. (2015). Additionally we compare the crack opening displacements. According to Sneddon & Mott (1946), the opening of a fully pressurized crack follows the equation with c = 2a/E and radius r ≡ x for the twodimensional model, and with c = 4a/(π E ) for a threedimensional penny-shaped crack.
In the left plot of Fig. 4, analytical and numerical results for the two-dimensional model are plotted. Shown are the displacements at two points, namely at position x 1 = 0.4a and x 2 = 0.8a. The analytical solution is determined according to Eq. (28) for each time step with the respective applied pressure p(t) for values p 0 , p 1 as above and time step Δt = 0.01 s. The simulations show a good agreement with the analytical calculations for the first stages (in the elastic range) but crack growth changes the overall situation. The analytical solution considers the displacements at a given (initial) crack length, whereas the computations now refer to larger cracks and so a relaxation and an overall smaller displacement follows. To illustrate the situation we plotted the displacement at position x = 1.2a, i.e., on a point which only opens during crack growth.
In the right plot of Fig. 4, the corresponding results of a three-dimensional computation are displayed, see Sect. 5.2 for the details. Here the displacements are measured at material points with local radius r 1 = 0.4a and r 2 = 0.8a. In both cases, we see again that the crack expansion is smaller than the analytical solution suggests, and due to the linear pressure increase, the crack width also increases linearly. The computation is stopped here after the crack has reached approximately 1.5 of the initial crack length a because the configuration then deviates from the assumptions of the analytical solution considerably.

Two notch crack propagation
In another two-dimensional simulation we keep mesh, boundary conditions, material and pressure values as above but adapt the phase field in such a way that now two notches are prescribed and loaded with the applied pressure. The resulting phase-field snapshots are presented in Fig. 5. It can be seen that the cracks grow towards the boundaries until the horizontal notch reaches the vertical one. A corresponding crack path can be found in Miehe and Mauthe (2016).

Crack propagation with a Rankine stress model
To illustrate the effect of the alternative crack driving force (5) in the context of pressure driven crack growth, the simulation of Sect. 3.2 is repeated. The setting and also the pressure remain the same and the related phasefield snapshots are shown in Fig. 6. Compared to the results from the energetic crack driving force, it can be seen that the crack branches later such that the boundaries influence the crack tip in their vicinity only. This effect stems from the fact that only tensile stress states contribute to crack growth here whereby in the energy driven model also shearing has a certain contribution, cf. Bilgen and Weinberg (2019).

Remarks
In our approach the applied pressure results in an additional energy contribution. Bourdin et al. (2012) presented a similar way to incorporate the pressure as an additional crack driving force, preserving Γconvergence of the regularized formulation. They extended the total potential energy by the work of a pressure forcep such that: For small strains and one pressurized crack this formulation is basically identical to ours in Eq. (17). To show this we set p 0 = 0 and p 1 =p in Eq. (19), and separate the pneumatic energy contribution of (17) to state whereby the last step presumes neglected boundary effects. Then the result equals the pressure related term of Bourdin et al. (2012), eq. (29). Clearly, the model of Bourdin et al. (2012) is restricted to linear elasticity whereas we state it in a general form here. Additionally we remark that the proposed type of pressure boundary condition can be used in the same way for two and three dimensional simulations, cf. Sect. 5.2.

Phase-field model with non-isotropic crack resistance
Geological and also organic materials are commonly not smoothly homogenous but have inclusions, softer and harder layers, and different strength in different directions. That influences also their resistance to fracture and therefore we will consider such properties in the following. We remark that we do not aim at modeling the material's anisotropy because there exist established strategies for linear and finite elasticity Marsden et al. (1994) and Holzapfel (2006). Here we focus on non-isotropic crack resistance.

Random variations in tensile resistance
At first we model variations of the crack resistance due to inhomogeneities in the solid. This can be done by , or using a stochastically 'smeared' approach. Here, the basic idea is to vary the critical parameters of the phase-field model, in particular the critical energy release rate G c and/or the critical stress σ c , by means of a superposed function ξ(x) that influences the crack resistance Y f of Eq. (4). The total potential energy (3) then has the form: One way to incorporate local variation is to model ξ(x) stochastically distributed over the whole domain.
To do so, we apply an operator scaling stable random field (OSSRF). Such a OSSRF maps the deviations of the material from its ideal state and has been used for phase-field modeling before, cf. Anders et al. (2011). The OSSRF is based on time and space continuous Lévy walks which are, figuratively speaking, a mix of long trajectories and short random movements. A stable Lévy motion is a self-similar process with a scaling parameter H > 0, the Hurst index. Operatorscaling self-similar processes allow the Hurst index to vary with the coordinates. Thus we choose for (x, t) ∈ Ω × [0, t tot ] and parameters 0 < H 1 , H 2 , H 3 < 1 a random field of the form where The variables θ 1 , θ 2 are normalized eigenvectors of a scaling matrix with given eigenvalues H −1 1 , H −1 2 , respectively. The field (32) is anisotropic if the spatial scaling relation has different Hurst indices in the -not necessarily orthogonal -directions of the eigenvectors θ 1 , θ 2 . The measure P(d y × ds) indicates a complex-valued independently scattered isotropic random Gaussian measure. The Hurst index H 3 refers to time; however, we set ξ(x, t) = ξ(x) and fix the temporal coordinate here. For mass conservation the condition has to be fulfilled. In Figure 7 the orientations and the random field are demonstrated exemplarily. With this condition, the OSSRF does not affect the physics of the model, see Anders et al. (2011) for a discussion. In contrast to Gaussian white noise the OSSRF allows to model long range dependencies in space and can map anisotropic behavior.

Proof of concept
To demonstrate the influence of the OSSRF ξ(x) we perform a mode-I tension test. A square domain of size 100 × 100 mm 2 with a prescribed notch is loaded by tension, see left plot of Fig. 8. The mesh consists of 200 × 200 quadratic B-spline elements; the lengthscale parameter is l c = 1 mm. The inhomogeneity of the material with parameters of Table 1 is modeled with field (32) with H 1 = 0.9, H 2 = 1 and t = t 0 , where t 0 is fixed. The results for two realizations are demonstrated in Fig. 8. In both cases the crack path is influenced by the random field. We remark that these investigations are only exemplarily. The OSSRF enables layer-like properties within a material, e.g. in shale formations or in plagued arteries. However, it may be difficult to adapt it to a given material behavior with preferred directions. Therefore, we suggest a deterministic tensor formulation in the following.

Deterministic conductance tensor
The elastic strain energy density (7) is per definition an isotropic tensor function whose arguments are based on the right Cauchy-Green deformation tensor C. It is invariant against rigid body movements, i. e. it holds Ψ e (C) = Ψ e QC Q T = Ψ e Q F T F Q T for arbitrary rotations Q ∈ SO(3) Holzapfel (2006); SO (3) is the special orthogonal group of three dimensional rotations. The formulation in invariants (7) describes material with isotropic properties, i. e.
for all Q ∈ SO(3). For materials with a preferred direction, for instance as a consequence of fibers, the energy formulation needs to be modified. The simplest way to account for a reinforced direction a is to introduce an additional term, a structural tensor a ⊗ a, in the energy form. One can easily verify that this also gives an isotropic tensor function, now in two arguments, Ψ e (C, a ⊗ a) = Ψ e QC Q T , Qa ⊗ a Q T . Additional directions may be considered by additional structural tensors. We adapt this approach here and consider material with isotropic elastic properties but a direction dependent resistance to cracking. As above we assume a modified resistance in a direction a = αa 0 whereby a 0 is a normalized vector and α its relative length. Each direction defines a conductance tensor α 2 a 0 ⊗ a 0 which now enters the crack resistance part of energy density (3) as an additional term. Formulating the crack density function (2) for two (or more) preferred directions it reads Fig. 7 Random field ξ with θ 1 = e i·0.3π , θ 2 = e i·0.2π , H 1 = H 2 = 0.5; the colors describe the intensity/magnitude of ξ (at a fixed time) Fig. 8 Mode-I crack propagation with a crack resistance controlled the random field ξ(x, t 0 ) γ (z, a 1 , a 2 , . . . ) = 1 2l c z 2 + l c 2 ∇z · ∇z + l c 2 α 2 1 ∇z · a 1 ⊗ a 1 · ∇z + l c 2 α 2 2 ∇z · a 2 ⊗ a 2 · ∇z + . . .
Similar to (34) it holds for an isotropic crack density function the equation for all Q ∈ SO(3); in the case of anisotropy the condition is valid only for all Q ∈ B ⊂ SO(3) with an appropriate subset B of the orthogonal group defined by the conductance vectors a 1 , a 2 …of preferred crack directions. There may be several directions of modified crack resistance which result in additional terms in (35) but in the following we focus on a unidirectional conductance for simplicity. Then, the crack density function can be written in the form where the relative vector length α plays the role of a weight, α 2 . Choosing α = 0 leads to isotropic material behavior as Eq. (37) reduces to Eq. (2). High values of α weight the conductance tensor and pronounce the preferred crack direction; a negative sign −α 2 would result in a higher crack resistance of that direction.

Mode-I tension test with different weights
To illustrate the influence of the weight α 2 in the conductance tensor, the mode-I tension test of Sect. 4.1 is repeated. The mesh consists here of 200×200 quadratic B-spline elements; the length-scale parameter arises Fig. 9 Phase-field crack of the mode-I tension test for a material with direction vector a 0 = (1/ √ 2, 1/ √ 2) and weights α = 1, α = √ 5, α = 5 (from left to right); crack driven by the λ-μ energy split as l c = 1 mm. The preferred direction of the material is declined by an angle of 45 • which corresponds to a direction vector a 0 = (1/ √ 2, 1/ √ 2). Varying the parameter α with α ∈ {1, √ 5, 5} gives different weights to the inclination; the related crack paths are shown in Fig. 9. It is noticeable that with increasing α the crack deviates stronger in the preferred direction; for α → ∞ the crack would always go in direction of vector a 0 . This result is confirmed by other authors who presented similar approaches, cf. Teichtmeister et al. (2017) and Liu and Juhre (2018).

Mode-II shear test with different weights and crack driving forces
To take both tensile and shear components into account a mode-II shear test is simulated in the following. While the mesh and the material parameters remain the same, the boundary conditions on the upper side are adapted, prescribing a displacement in x-direction. The direction vector is given by a 0 = (1, 0) to simulate a preference in horizontal direction. We perform the simulation with two different crack driving forces here, namely the energetic force with energy split (16) and the Rankine stress model (5). Comparing both crack paths we see a similar effect of the conductance tensor, i.e., the crack path is influenced and changes in the preferred direction. For bigger weights α the crack propagates stronger in the horizontal direction, see Figs. 10 and 11.
The main difference is that the crack computed with the λ − μ split is bent, likely because of the coupled influence of tension and shear, whereas the Rankine stress model results in a straight crack. The results coin-cide with the results of the tension test and confirm the effect of weight α. The related load-deflection curves are presented in Fig. 12. In both cases it is observed that with growing weight the applied force and/or the prescribed displacement at cracking increase. We remark that this model has been used with the λ-μ energy split of linear elasticity in Clayton and Knap (2015), Teichtmeister et al. (2017) and Liu and Juhre (2018) for polycrystals but it works for deformations of rubbery material as well.

Mode-I tension test with different directions
In a next step we perform the mode-I tension test with different direction vectors. Based on the previous results the weight α is fixed at 5 √ 2, α 2 = 50. Two direction vectors are chosen, a 01 = (1/ √ 2, 1/ √ 2) and a 02 = (1/ √ 5, −2/ √ 5). With both crack driving forces, energy split (16) and Rankine stress (5), we get coinciding results, exemplarily displayed in Fig. 13. The left plot shows the crack path for an isotropic material model, middle and right plot show the crack growing in the direction of a 01 and a 02 . It can be observed that the cracks immediately grow in the preferred direction.

Randomly weighted direction dependence
Here we now combine both approaches to anisotropic crack growth, the stochastic field (32) and the conductance tensor (37). The crack resistance of the material is assumed to change locally, prescribed by a direction vector a = αa 0 . Because a random change in direction may easily result in a nearly isotropic behavior, Phase-field crack of the mode-I tension test for an isotropic material behavior and for an anisotropic crack resistance with different direction vectors a 01 = (1/ √ 2, 1/ √ 2) and a 02 = (1/ √ 5, −2/ √ 5) and length α = 5 √ 2 (from left to right); crack driven by the Rankine stress we define here a conductance direction a 0 and change only its influence, namely the vector length. It shall be defined by the OSSRF, i. e α 2 = ξ(x).
The effect is demonstrated by means of the mode-II shear test of Sect. 4.2.2. The direction vector is a 0 = (1, 0) again, define a compliance in horizontal direction. The results are displayed in Fig 14. For the random field the eigenvectors and the Hurst indices are fixed but the amplitude of ξ(x) varies between [0, 10] in the first case and between [0,25] in the second case. It can be observed that the influence of anisotropy is differently strong, resulting in meaningful looking crack paths.

Numerical Examples
In this section we provide some sample calculations for pneumatic fracture in a material with layer-like anisotropic crack resistance. Such materials are, for example, shale formations, laminated safety glass, multiply rubber or vascular vessels. The crack growth is driven by the moving boundary model of Sect. 3.1 and includes now the direction dependent crack density function (37). In the following we will consider linear material behavior, but rather the theory can be also adapted to finite deformations.

Pressure driven crack growth in plane stress
We start with the example in Sect. 3.1 where a quadratic plate with a prescribed notch in its center is considered, see left plot of Fig. 15. The material parameters and boundary conditions remain unchanged. For crack anisotropy three conductance vectors are defined, all with relative length α = 5 √ 2. The directions are chosen as a 01 = (1/ √ 2, 1/ √ 2), a 02 = (2/ √ 5, 1/ √ 5) and a 03 = (1/ √ 5, −2/ √ 5). The initial pressure is given by The notch is initially horizontal and so the crack would expand as displayed in Fig. 2. Because of the preferred direction the crack here declines and propagates immediately in the preferred direction as shown in the three snapshots of Fig. 15. In the displayed anisotropic computations the crack driving force is based on the maximum of principal stress (5) and as soon as the critical stress is exceeded the crack grows in the preferred direction. The crack patterns for the energetic crack driving force are identical and therefore not displayed. In both cases the anisotropic crack resistance is modeled reasonably.
Next a layered structure is investigated, see Fig. 16. The mesh consists of 200 × 200 quadratic B-spline elements, all boundaries are fixed in each direction and the material parameters remain as previously introduced. The pressure is given by p 0 = 10 N/mm 2 as well. In each layer another preferred direction occurs. The direction vectors are chosen alternatingly to be a 01 = (1/ √ 2, 1/ √ 2) and a 02 = (−1/ √ 2, 1/ √ 2). In the center of the plate a small notch is prescribed by fixing z = 1. The propagating phase-field cracks are shown in Fig. 16 at different instants of time. It can be seen how the crack kinks at the layers' interfaces. The crack zigzags until it reaches the boundaries and so follows the prescribed anisotropic crack resistivity.

Pressure driven crack growth in a cylinder with centered crack
Although our examples have been plane, the proposed formulations are not restricted to two dimensions. In this subsection we examine the well known pennyshaped crack problem in three dimensions, studied for example in Lee et al. (2016), Miehe and Mauthe (2016) and Zhou et al. (2018). We consider a cylinder with a diameter of 50 mm and a height of 50 mm. The mesh consists of 115 000 linear brick elements (Lagrange ansatz functions). The material parameters remain as before, are based on linear material behavior and correspond to concrete, Table 1. The outer boundaries of the structure are constrained in all directions. In the center of the specimen an initial crack is modeled by a squared surface where the phase-field parameter is set z = 1. The mesh and the initial configuration are demonstrated in Fig. 17. The pressure is induced by the moving boundary condition suggested of Sect. 3.1 with an amount of p 0 = 10 N/mm 2 and an anisotropic crack resistance is assumed. In one example the direction vector is chosen to be a 01 = (1/ √ 3, 1/ √ 3, 1/ √ 3) so that the preferred direction is a spacial diagonal. In another example the direction vector is set to a 02 = 1 √ 27 (1, 5, 1) which means that the y-direction is strongly preferred. The computed crack paths with the related isosurfaces are presented for z = 0.9 in Fig. 18. We use the Rankine stress model (5) here, i.e., the crack surface grows as 2), α = 5 √ 2, and phase-field cracks at different time steps, which t n ∈ {1, 20, 32} soon as the critical stress is exceeded. It is observed that a crack surface arises which spreads through the whole cylinder and catches the anisotropic crack resistance.

Concluding remarks
In this contribution we present an extended phase-field model for pneumatic anisotropic fracture. The pressure driven crack growth is modeled with a moving bound-ary condition introduced in such a way that the pressure acts as a stress in the phase-field crack zone. The suitability of this approach is demonstrated with several examples. For modeling non-isotropic crack resistance two different approaches are discussed, namely a random field function in the Griffith energy and a structural conductance tensor in the crack density function. Both methods give the expected results but the conductance tensor approach is found to be more convenient for modeling layer-like materials. Eventually, regard- (1, 5, 1) (lower row) ing the applications in geological and biological materials both methods are combined.
To show the capabilities of the proposed strategy, two and three dimensional pneumatic fracture simulations have been performed with different phase-field crack driving forces. Hereby we found that all formulations can be used within our approach to pressuredriven crack growth leading to similar results and catching the direction dependent crack resistance properly.