# Second-order phase-field formulations for anisotropic brittle fracture

We address brittle fracture in anisotropic materials featuring two-fold and four-fold symmetric fracture toughness. For these two classes, we develop two variational phase-field models based on the family of regularizations proposed by Focardi (Focardi, M. On the variational approximation of free-discontinuity problems in the vectorial case. Math. Models Methods App. Sci., 11:663684, 2001), for which Gamma-convergence results hold. Since both models are of second order, as opposed to the previously available fourth-order models for four-fold symmetric fracture toughness, they do not require basis functions of C1-continuity nor mixed variational principles for finite element discretization. For the four-fold symmetric formulation we show that the standard quadratic degradation function is unsuitable and devise a procedure to derive a suitable one. The performance of the new models is assessed via several numerical examples that simulate anisotropic fracture under anti-plane shear loading. For both formulations at fixed displacements (i.e. within an alternate minimization procedure), we also provide some existence and uniqueness results for the phase-field solution.

## Authors

• 2 publications
• 11 publications
03/07/2022

### A DG/CR discretization for the variational phase-field approach to fracture

Variational phase-field models of fracture are widely used to simulate n...
05/13/2022

### A fourth-order phase-field fracture model: Formulation and numerical solution using a continuous/discontinuous Galerkin method

Modeling crack initiation and propagation in brittle materials is of gre...
03/01/2021

### Mixed variational formulations for structural topology optimization based on the phase-field approach

We propose a variational principle combining a phase-field functional fo...
10/13/2019

### Variational Phase Field Formulations of Polarization and Phase Transition in Ferroelectric Thin Films

Electric field plays an important role in ferroelectric phase transition...
12/23/2020

### A general theory for anisotropic Kirchhoff-Love shells with embedded fibers and in-plane bending

In this work we present a generalized Kirchhoff-Love shell theory that c...
12/01/2021

### Unconditional well-posedness and IMEX improvement of a family of predictor-corrector methods in micromagnetics

Recently, Kim Wilkening (Convergence of a mass-lumped finite element...
05/20/2021

### A cahn-Hilliard multiphase system with mobilities for wetting simulation

This paper tackles the simulation of the wetting phenomenon using a phas...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Predicting the fracture behavior in brittle materials with anisotropic fracture toughness is still a challenge. The extension to the anisotropic case of selection criteria for the crack direction which have proved successful in fracture mechanics of isotropic materials, such as the principle of local symmetry [22, 12] and the maximum energy release rate criterion [28], is by no means trivial [34]. The generalized energy release rate criterion [37, 24, 34], where the crack propagation direction is chosen by maximizing the ratio of the energy release rate to the fracture toughness with respect to the propagation angle, was shown to correlate successfully with experiments on materials exhibiting both weak [29] and strong [39] anisotropy in the fracture toughness, see also the detailed review and discussion in [33].

The variational phase-field approach to fracture, pioneered by Bourdin et al. [5] and first proposed as the regularization of Francfort and Marigo’s variational fracture formulation [17], was recently extended to materials with anisotropic fracture toughness. Hakim and Karma [26, 27] proposed a phase-field model for materials featuring two-fold weak anisotropy and elucidated several features of crack propagation in these materials. This model was extended to the three-dimensional setting including geometric nonlinearity and nonlinear elasticity [11] and to the consideration of multiple cleavage planes using multiple phase-field variables [36]. With the purpose of tackling phenomena such as sawtooth crack patterns and forbidden crack directions, which are directly related to four-fold strongly anisotropic surface energies, Li et al. [32] formulated a fourth-order anisotropic variational fracture model inspired by phase-field models of crystal growth. Similarly, Teichtmeister et al. [40]

formulated second- and fourth-order phase-field models based on structural tensors for two-fold and four-fold symmetry, respectively, and extended the framework to large deformations. Li and Maurini

[33] improved and simplified the phase-field model in [32] in order to enable an analytical solution for the optimal crack profile, while keeping the fundamental properties of the phase-field approximation. Once again, four-fold symmetry was addressed with a fourth-order model.

In this paper, we address brittle fracture in anisotropic materials featuring two-fold and four-fold symmetric fracture toughness and develop two variational phase-field models based on the family of regularizations proposed by Focardi [16], who also proved their -convergence. The main novelty lies in the model for four-fold symmetric anisotropy: as opposed to all previously available models for the same type of anisotropy, our proposed model is of second order, hence it does not require higher-continuity basis functions nor mixed variational principles for finite element discretization. For the four-fold symmetric formulation we show that the standard quadratic degradation function is unsuitable and propose a procedure to derive a suitable one. The performance of the new models is assessed via several numerical examples that simulate anisotropic fracture under anti-plane shear loading. For both formulations at fixed displacements, we also provide some existence and uniqueness results for the phase-field solution.

The remainder of this manuscript is organized as follows. Section 2 introduces basic concepts on anisotropy of fracture properties that serve as basis for the following investigation. In Section 3, we recall the phase-field formulation of (isotropic) brittle fracture based on the Ambrosio-Tortorelli regularization and introduce the main results in [16], which are our point of departure for the new developments in the subsequent sections. We show that the models based on the Focardi regularization combined with the Euclidean norm, which we term isotropic Foc- models, can be viewed as a family of models which includes the well-known AT-2 model as a special case (corresponding to the isotropic Foc- model). In Section 4, we focus on the isotropic Foc- model. We show that the standard quadratic degradation function is unsuitable for this model, we derive a new suitable one, and we test the ensuing formulation on an anti-plane shear numerical setup. In Section 5, we develop the anisotropic counterparts of the isotropic Foc-2 and Foc-4 formulations, which model an anisotropic fracture behavior with two- and four-fold symmetric fracture toughness, respectively. Section 6 illustrates the numerical tests on the two anisotropic formulations for different anti-plane shear setups. We draw the main conclusions in Section 7. Appendix A reports some existence and uniqueness results on both anisotropic formulations.

## 2 Basic concepts on anisotropic fracture toughness

Materials with anisotropic fracture behavior are characterized by an orientation-dependent fracture toughness or critical energy release rate . One of the examples is sketched in Figure 1a, where it is assumed that there exists a plane of anisotropy within the body such that for any point and with the polar coordinate system within associated to , is viewed as a function of the polar angle , that is, . This is in contrast to the simple isotropic case when . In this paper, we consider the two-dimensional setting by taking , which implies that the anisotropy plane and coincide. For modeling purposes, it is furthermore assumed that

 Gc(θ)=G0γ(θ), (1)

where is a positive constant and is a periodic trigonometric function. A suitable form of this function reads

 γ(θ)=γk(θ):=1+τcos(k(θ−ω)), (2)

with , , and . The parameter defines the direction corresponding to the largest fracture toughness and hence it is referred to as the strongest or principal material direction. The parameter defines the magnitude of the largest fracture toughness and is called the anisotropy strength. A useful induced quantity for characterization of the anisotropy strength is the ratio . Finally, is the periodicity parameter such that the function is -periodic. The above gives rise to the notion of a -fold symmetric fracture toughness .

In the present work, we are interested in the cases and in (2), i.e. in the so-called two- and four-fold anisotropy cases. We give the corresponding explicit expressions of for future reference:

 γ2(θ)= 1+τcos(2(θ−ω)) = 1+τ[cos(2ω)(cos2(θ)−sin2(θ))+2sin(2ω)sin(θ)cos(θ)], (3)

with , , and

 γ4(θ)=1+ τcos(4(θ−ω)) =1+ τ[cos(4ω)(cos4(θ)−6sin2(θ)cos2(θ)+sin4(θ)) +4sin(4ω)sin(θ)cos(θ)(cos2(θ)−sin2(θ))], (4)

with , . The polar plots of both functions are continuous closed curves surrounding the origin of the polar coordinate system. According to standard terminology, the anisotropy is said to be weak if is convex and strong if is non-convex. The convexity of the function is known to be equivalent to the convexity of the region enclosed by the polar plot of , which is guaranteed if

 γ(θ)+γ′′(θ)≥0∀θ∈[0,2π) (5)

With the expressions of and in (3) and (4), it is straightforward to verify that

 γ2(θ)convex⇔τ∈[0,13] (6)

and

 γ4(θ)convex⇔τ∈[0,115] (7)

regardless of the value of . The polar plots of and , as well as and for a fixed and different magnitudes of are depicted in Figure 1b. With the increase of , the elliptic shape of and turns into a peanut-like shape, whereas for and one has a transformation from super-ellipse to ’butterfly’.

In the following, when considering the cases of isotropic and anisotropic fracture toughness we refer to them as ’isotropic fracture’ and ’anisotropic fracture’, respectively. In the former case, we assume such that, owing to (1) where , it is simply . In the latter case, the term ’anisotropic fracture’ is accompanied by the specification ’with -fold symmetric ’.

## 3 Variational formulation of brittle fracture and phase-field regularization

In this section, we briefly recall the phase-field formulation of (isotropic) brittle fracture from [5, 6, 7, 8], based on the so-called Ambrosio-Tortorelli regularization [2, 1]. This is followed by the presentation of the main results by Focardi [16], which are our point of departure for the new developments in the subsequent sections. Finally, we briefly outline the incremental minimization problem which is solved for computing a quasi-static phase-field evolution.

Let , or be an open and bounded domain representing the configuration of a -dimensional body, and let and be the (non-overlapping) portions of the boundary of which are fixed and subjected to the imposed displacement , respectively (Figure 2a). Let also represent the crack set, that is, the set of points where the displacement field is discontinuous.

### 3.1 Ambrosio-Tortorelli regularization

The variational approach to (isotropic) brittle fracture by Francfort and Marigo [17] relies on the energy functional

 E(u,Γc)=Eel.(u,Γc)+ES(Γc)=∫Ω∖ΓcΨ(ε(u))dx+G0∫ΓcdHd−1, (8)

and the related incremental minimization problem. In (8), is the elastic energy density, a function of the infinitesimal strain , where is the displacement field and is the gradient operator, and is the Hausdorff measure. The terms and respectively represent the elastic energy stored in the cracked body and the surface energy associated to the crack. Upon integration, the latter reads where, in simple terms, and are the length and the surface area of when and , respectively. Thus the expression of the surface energy associated to the crack is consistent with Griffith’s theory [23] with a fracture toughness .

The regularization of (8) à la Ambrosio-Tortorelli [2] developed by Bourdin-Francfort-Marigo [5, 6, 7, 8], which forms the basis for a variety of (isotropic) fracture phase-field formulations, reads as follows:

 E(u,α)=Eel.(u,α)+ES(α)=∫Ωg(α)Ψ(ε(u))dx+G0cw∫Ω(1ℓw(α)+ℓ|∇α|2)dx, (9)

with as the phase-field variable, which takes the value on , decays smoothly to in a subset of and vanishes in the rest of the domain, as sketched in Figure 2b. With this definition, the limits and represent the fully broken and the intact (undamaged) material phases, respectively, whereas the intermediate range mimics the transition zone between them. The function , often denoted as the degradation or the coupling function, is responsible for the material stiffness degradation, the function , often termed the local dissipation function [38], defines the (decaying) profile of , whereas the parameter controls the thickness of the localization zone of , i.e. of the transition zone between the two material states.

The specific choice of the functions and in (9) establishes the rigorous link between (8) and (9) when via the notion of -convergence, see e.g. [9, 10, 1], also giving a meaning to the induced constant . Thus, is a continuous monotonic function that fulfils the properties: , , and for , see e.g. [38]. The quadratic polynomial is the simplest choice. The function is continuous and monotonic such that , and for . The constant is a normalization constant in the sense of -convergence. The two suitable candidates for which are widely adopted read , such that , respectively.

The two instances of formulation (9) containing the aforementioned choices for and are typically termed AT-1 and AT-2 models, see Table 1. AT stands for Ambrosio-Tortorelli type of regularization, see [2, 1].

The main difference between the two models is that AT-1 leads to the existence of an elastic stage before the onset of fracture, whereas using AT-2 the phase-field variable starts to evolve as soon as the material is loaded, see e.g. [38, 3, 35] for a more detailed explanation, as well as the discussion in Section 4.1.

### 3.2 Focardi regularization

In Focardi [16], the Ambrosio-Tortorelli regularization of the free-discontinuity problem is generalized so as to include an orientation-dependent surface energy, and the related -convergence result is obtained. Using the notation in (8) and (9), the unregularized energy functional considered by Focardi and its regularized (phase-field) counterpart read

 E(u,Γc):=∫Ω∖ΓcΨ(ε(u))dx+G0∫Γcφ(n)dHd−1, (10)

and

 E(u,α):=∫Ωg(α)Ψ(ε(u))dx+G0∫Ω(1ℓqw(α)+ℓp−1pφp(∇α))dx, (11)

respectively. In the above,

is the normal unit vector to

, is a norm, and ; with . The coupling function is supposed to fulfill the standard properties listed in Section 3.1. Theorem 3.1 in [16] establishes the -convergence result between and under no specific assumptions for the norm .

For any choice of , one may say that (11) provides a -parametric family of regularizations to the formulation in (10). We denote this family as the Foc- models.

Suppose now that is a Euclidean norm. In this case, since is a unit vector, we simply arrive at in (10) — this then becomes identical to (8) — and we also obtain in (11). In the following, we term the -parametric family of formulations in (11) with given by the Euclidean norm the isotropic Foc- models (when the reference to isotropy is clear from the context, the term isotropic is omitted).

The isotropic Foc- model with can be viewed as an alternative to the Ambrosio-Tortorelli approximations AT-1 and AT-2, and, in particular, a generalization of the AT-2 case. In Table 2 the isotropic Foc-2 and Foc-4 models are reported. One immediately notices that the Foc-2 formulation is identical to the AT-2 model, provided that the degradation function in (11) is taken as . The Foc-4 formulation turns out to be a more interesting case and is thoroughly investigated in Section 4.

However, the generality of (10) and (11) is significantly broader, since does not have to be a Euclidean norm. In Section 5, we explain how the norm in (10) — hence in (11) — as well as the parameter in (11) can be chosen based on the anisotropy functions , from (2), and thus used to model anisotropic fracture (we will find out that the most natural choice is ). This results in what we term the anisotropic Foc- models.

Table 3 summarizes the terminology introduced so far.

### 3.3 Incremental variational problem

We consider a quasi-static loading process in which the monotonically increasing imposed displacement , with denoting the discrete pseudo-time step parameter, is prescribed on . It is assumed that, during this process, the crack surface evolves incrementally and in an irreversible manner, that is, .

With the phase-field energy functional defined by (9) or (11), the state of the system at a given loading step is given by

 (u,α)=argmin{E(v,β):v∈V¯un,β∈Dαn−1}. (12)

Here

 V¯un:={u∈W1,2(Ω;Rd):u=0onΓD,0,u=¯unonΓD,1}

is the space of kinematically admissible displacements at step , with as the usual Sobolev space of functions with values in , and

 Dαn−1:={α∈W1,p(Ω;R),p>1:α≥αn−1inΩ}

is the admissible space for at step , with as the usual Sobolev space of functions with values in and known from the previous step. The condition in is used to enforce the irreversibility of the crack phase-field evolution. It is the backward difference quotient form of in . Also, in , we take for the AT-1 and AT-2 models in (9), and for the Foc- model in (11).

 (13)

where and denote the directional derivatives of with respect to and , respectively. The displacement test space in (13) is

 V0:={u∈W1,2(Ω;Rd):u=0onΓD,0∪ΓD,1}.

Conditions (13) characterize, in general, a local minimum of .

### 3.4 Handling of irreversibility

In what follows, we incorporate the irreversibility constraint via simple penalization111 A detailed discussion about various available options to handle the crack phase-field irreversibility conditions, their equivalence and related numerical comparisons can be found in our earlier publications [20] and [15]. , thus arriving at

 F(u,α):=E(u,α)+^λ2∫Ω⟨α−αn−1⟩2−dx, (14)

where is given by (9) or (11), and is the penalty parameter. In [20], an analytical procedure for the ’reasonable’ choice of a lower bound for that guarantees a sufficiently accurate enforcement of the constraint is devised. The optimal is given by

 (15)

where is a user-prescribed irreversibility tolerance threshold (in [20], a practical value for is suggested as ).

With defined by (14), the sought solution at a given loading step is given by

 (u,α)=argmin{F(v,β):v∈V¯un,β∈W1,p(Ω;R)}. (16)

The optimality condition in this case simplifies to

 (17)

where , and

 Fα=Eα+^λ∫Ω⟨α−αn−1⟩−βdx.

Notice that in (16) and (17) we again take for the AT-1 and AT-2 models, and for the Foc- model.

## 4 The isotropic Foc-4 model

In Section 3.2 we have seen that the isotropic Foc- models with can be viewed as a family of models which includes AT-2 as a special case. In view of the future extension to anisotropic fracture (to be addressed in Section 5), and in particular to the two- and four-fold anisotropy cases, we are especially interested in the isotropic Foc- and Foc- models. Since the former (with the standard quadratic degradation function) is identical to the AT-2 model, whose properties are already well known, in this section we concentrate on the isotropic Foc- model (which we simply denote as Foc- model throughout this section). The corresponding energy functional is obtained from (11) by setting and , which yields

 E(u,α)=∫Ωg(α)Ψ(ε(u))dx+G04∫Ω(3bwα4ℓ+ℓ3|∇α|4)dx, (18)

where .

As follows, we first carry out the analysis of a one-dimensional bar under homogeneous conditions. Through this analysis we show that the standard quadratic degradation function is not suitable for the Foc-4 formulation, and we then derive analytically a new suitable expression for . Subsequently, we study the one-dimensional localization behavior and construct the so-called optimal profile of induced by the fracture energy term in (18). We compare it with the optimal profiles of the AT-1 and AT-

2 models. The knowledge of this profile also allows us to estimate the magnitude of the penalty parameter

in (14) to be used for the enforcement of irreversibility in computations. Finally, the new formulation obtained combining the Foc-4 model and the new degradation function is tested on a simple anti-plane shear numerical setup.

### 4.1 One-dimensional homogeneous behavior

The main tool for assessing the fundamental properties of a fracture phase-field formulation — in particular, the impact of ingredients such as — is the analysis of a one-dimensional elastic bar under tension [38, 35, 4, 31]. In the following, we recall the corresponding results for the AT-1 and AT-2 models, and apply the same analysis to understand the role of the degradation function and to devise a suitable one for the Foc-4 model.

Following [38], we consider a one-dimensional elastic bar of length , clamped at one end () and subjected to an imposed displacement at the other end (). The elastic energy density function reads with as the Young’s modulus. In this pure tensile loading case, the phase-field irreversibility constraint is not necessary, hence we can set . The strong form of the equilibrium equations in reads

 ¯σ′(x)=0with¯σ=g(α)¯ϵ, (19)

(regardless of the model), whereas the damage criterion reads

 12g′(α)¯ϵ2ℓ+1cw(w′(α)ℓ−2ℓα′′(x))=0, (20)

for the AT-1 and AT-2 models, and

 12g′(α)¯ϵ2ℓ+3(1bwα3ℓ−ℓ3[α′(x)]2α′′(x))=0, (21)

for the Foc-4 model. Here, we use the re-scaled quantities and .

Let us focus on the study of the homogeneous solution, where the damage and the stress and strain fields are uniform along the bar [38, 35]. In this case, due to the disappearance of the spatial derivatives of , equations (19)-(21) provide the explicit - and - relations:

• for the AT-1 model,

 ¯ϵ=√3811−α,¯σ=√38(1−α)32, (22)
• for the AT-2 model,

 ¯ϵ=√α1−α,¯σ=α12(1−α)32, (23)
• for the Foc-4 model,

 ¯ϵ=√−6bwα3g′(α),¯σ=g(α)√−6bwα3g′(α). (24)

In the first two cases, we already substituted for the standard quadratic expression, whereas we are considering a generic for the third case. Figure 3 depicts the - and - plots for the AT-1 and AT-2 models, as well as for the Foc-4 model with , .

The curves obtained with the AT-1 and AT-2 models are well known. From them it can be inferred that the AT-1 model features an initial linearly-elastic phase of behavior, after which damage starts to occur. The peak value of the stress is reached for a value of given by . Conversely, the AT-2 model has no purely elastic phase, as the evolution of starts from the onset of loading and the peak value of the stress is reached for .

The curves corresponding to the Foc-4 model are qualitatively similar to those of the AT-2 model. However, the non-linearity prior to the peak stress is even more pronounced and the value of is even larger. Indeed, for it is . If the polynomial degree of the degradation function is increased, is reduced, however reaches at increasingly large values of .

We now aim at devising a degradation function for the Foc-4 model that, in addition to fulfilling the basic properties of degradation functions (see Section 3.1), leads to a - behavior similar to the one of the AT-1 model. We thus enforce the following property222Clearly, the representation on the right-hand side of (25) is not unique. Even if we keep the polynomial representation at the denominator, another option that preserves the same behavior of the right-hand side is e.g. . However, it can be shown that the corresponding -parametric family of does not fulfil . As a result, ceases to be automatically bounded from above by 1, and incorporation of the constraint in within the minimization problem for in (18) becomes necessary.

 −α3g′(α)!=s1−αm,α∈(0,1), (25)

where and , . Solution of the differential equation in (25), with the boundary conditions and , reads

 g(α)=1−(1+4m)α4+4mαm+4,m≥1. (26)

Secondly, for given by (24) with as in (26), we seek to yield in . Figure 4 illustrates as a function of . It can be clearly seen that the desired property for this quantity is fulfilled for , yielding

 g(α)=(1−α4)2. (27)

Figure 5 shows the -, - and - relations for the Foc-4 model with the degradation function (27), also in comparison with the AT-1 and AT-2 cases. It is clear that the constructed eighth-order polynomial function in (27) leads to the desired behavior.

### 4.2 One-dimensional localization behavior

If is a fully developed crack, the phase-field profile that represents the regularized counterpart of can be obtained by minimizing the fracture surface energy functional on a suitable admissible set, namely

 α=argmin{ES(β),β∈W1,p(Ω;R),p>1:β(Γc)=1,β≥0}, (28)

see Figure 6a for an illustrative sketch in two dimensions. Using a slicing argument, one may resort to a one-dimensional setting in (28): we assume with such that is represented by the point , and we also denote . In this setting, the minimizer can be constructed explicitly.

The strong formulation of the boundary value problem to compute for the AT-1 and AT-2 models is given by

 (29)

Note that in (29) the inequalities hold for AT-1 and the corresponding equalities for AT-2. The solutions representing the optimal phase-field profiles read

 (30)

for AT-1, and

 ˆα(t)≈exp(−|t|ℓ),t∈(−L,L), (31)

for AT-2, see e.g. [20] for details.

For the Foc- model, the optimal profile must be computed from

 (32)

and is straightforwardly obtained as 333 Assuming a general solution of the form , with as an unknown parameter and , and substituting it into the differential equation, one obtains . Using the boundary conditions, and are derived. With the assumption , the simplified approximate result in (33) follows immediately.

 ˆα(t)≈exp(−14√bw|t|ℓ),t∈(−L,L). (33)

The plots of in (30), (31) and (33) are shown in Figure 6b. It is evident that the AT- and the Foc- models have very similar optimal profiles. A practical outcome of this observation is that, when using the penalty method to enforce irreversibility, we can choose the penalty parameter for the Foc- formulation using the same estimate valid for the AT- formulation, see eq. (15).

### 4.3 Anti-plane shear test results

We now test the formulation consisting of the Foc- model with the new degradation function in (27) on an anti-plane shear loading numerical setup, see Figure 7. A two-dimensional rectangular domain containing a slit along is subjected to incremental anti-plane displacements , applied on and , respectively.

The loading conditions imply with , such that and . As a result, with as the shear modulus (Lamé’s second parameter). In the computations, we furthermore set , , , . As already mentioned, the degradation function is taken as . Also, we set with , see Remark 3.4. The applied incremental displacement is given by , , with as the loading increment. We employ the numerical package FreeFem++ [25]. Both the displacement field and the crack phase-field are approximated using -triangles. The finite element mesh is a pre-adapted one, which is refined in the region of where crack propagation is expected. The minimum and maximum characteristic mesh sizes are given by , where and stand for the mesh size inside and outside of the refined region, respectively.

Figure 8a-b visualizes the fields computed with the AT-2 and Foc-4 models with quadratic degradation function at the same loading step after the peak load. It is evident that, in both cases, does not vanish outside of the localization zone, which is a direct consequence of the accumulation of damage since the onset of loading discussed in Section 4.1. For the given setup at the chosen loading step, it is for AT-2 and for Foc-4. Interestingly, both values are significantly lower than the value of obtained for the one-dimensional homogenous responses of the respective models, see Section 4.1. Finally, Figure 8c shows results for the Foc-4 model with the degradation function in (27), where no spurious damage accumulation outside of the localization region is obtained.

## 5 The anisotropic Foc-2 and Foc-4 models

In this section, we develop the anisotropic counterparts of the isotropic Foc-2 and Foc-4 formulations. To this end, we suitably choose the norm in (10) — and, as a result, in (11) — based on the knowledge of the anisotropy functions , in (3) and (4). In particular, is used to induce the norm in the Foc-2 formulation, whereas is employed for defining in the Foc-4 formulation. We term the resulting formulations anisotropic Foc- and Foc- models, see Table 3.

### 5.1 Derivation of the anisotropic models

Since in (10) is assumed to be a function of the unit normal , whereas the anisotropy functions take as argument the polar angle , the first straightforward step is to express via the components of . Let represent a pre-existing crack, point be the coordinate of its tip, and be a possible extension of starting from . Also, let us define the Cartesian and polar coordinate systems associated to , as well as the unit normal to at . The angle is defined as the angle between the positive direction of the -axis and the crack extension , see Figure 9. With this definition, is related to the components of by and . Substituting in (3) and (4), we obtain

 γ2(n)=n2x+n2y+τ[cos(2ω)(n2y−n2x)−2sin(2ω)nxny], (34)

and

 γ4(n)=(n2x+n2y)2+τ[cos(4ω)(n4y−6n2xn2y+n4x)−4sin(4ω)nxny(n2y−n2x)]. (35)

where we also exploited the following representations of the unity

 (36)

We can now introduce the mappings for two- and four-fold anisotropy, which respectively read:

 φ(a,b):={a2+b2+τ[cos(2ω)(b2−a2)−2sin(2ω)ab]}12, (37)

and

 φ(a,b):={(a2+b2)2+τ[cos(4ω)(b4−6a2b2+a4)−4sin(4ω)ab(b2−a2)]}14. (38)

We call in (37) and (38) - and -induced norms, respectively, to be used in the Focardi discrete variational formulation in (10). Noting now that the above norms to be used in (10) have the form , , and that the consequent terms , to appear in (11) read , the most natural choice is to set . This yields

 φ2(∇α)=|∇α|2+τ[cos(2ω)(α2,y−α2,x)−2sin(2ω)α,xα,y], (39)

as the ingredient of (11) when and is the -induced norm, and

 φ4(∇α)=|∇α|4+τ[cos(4ω)(α4,y−6α2,xα2,y+α4,x)−4sin(4ω)α,xα,y(α2,y−α2,x)], (40)

as the ingredient of (11) when and is the -induced norm. We denoted as the components of .

The norm in (39) can be also written as

 φ2(∇α)=B∇α⋅∇α=Bijα,iα,j (41)

with . Here

 B=I−τDwithD=[cos(2ω)sin(2ω)sin(2ω)−cos(2ω)] (42)

is a symmetric positive definite second-order tensor, and is the two-dimensional second-order unit tensor. The norm in (40) can be written as

 φ4(∇α)=B∇α⊗∇α⋅∇α⊗∇α=Bijklα,iα,jα,kα,l (43)

with . Here

 B=I⊗I−τD (44)

is a positive definite fourth-order tensor exhibiting minor and major symmetries, with components

 D1111=D2222 = −cos(4ω) (45) D1122=D2211=D1212=D2121=D1221=D2112 = cos(4ω) (46) D1112=D1121=D1211=D2111 = −sin(4ω) (47) D2221=D2212=D2122=D1222 = sin(4ω) (48)
###### Remark 1

The representation of in terms of the components of the unit vector is not unique. However, the choice given by (36) — and resulting in the appearance of the terms and in (37) and (38), respectively — is the only appropriate one in this case, as it preserves an identical dimension for all terms entering the corresponding .

###### Remark 2

The choice of the exponents and in (37) and (38) is a natural one, since the limiting case of , with , turns the corresponding into the Euclidean norm.

###### Remark 3

The mappings in in (37) and (38) indeed define two norms, since owing to

 (49)

and

 (50)

the following norm-equivalence result holds:

 (1−τ)|(a,b)|≤φ(a,b)≤(1+τ)|(a,b)|,

for all and belonging to the appropriate range.

We can now finalize the proposed phase-field models for anisotropic fracture with two- and four-fold symmetric fracture toughness. For two-fold symmetry, i.e. for , we propose the following energy functional:

 E(u,α):=∫Ω(1−α)2Ψ(ε(u))dx+G02∫Ω{α2ℓ+ℓφ2(∇α)}dx, (51)

with given by (39) or equivalently by (41), whereas for four-fold symmetry, i.e. for , we propose

 E(u,α):=∫Ω(1−α4)2Ψ(ε(u))dx+G04∫Ω{3bwα4ℓ+ℓ3φ4(∇α)}dx. (52)

with given by (40) or equivalently by (43), and where also . We term the formulations (51) and (52) the anisotropic Foc- and Foc- models, respectively.

The questions of existence and uniqueness of the solution to the minimization problem for in (51) and (52) is a delicate task. Indeed, it is very well known that for the isotropic case, may admit infinitely many solutions, or none. In the light of the so-called alternate minimization scheme, see e.g. [5, 6, 7, 8], applied to resolve numerically, one typically restricts the discussion to studying the well-posedness of the problems when is fixed, and when is fixed. The former case is simple: the problem is well-posed due to the convexity of with respect to . Well-posedness of the latter minimization problem with given by (51) and (52) is discussed in Appendices A.1 and A.2, respectively. In the appendices, this is formulated as an -parametric result.

### 5.2 Relation to existing phase-field models for anisotropic fracture

It is straightforward to recognize that the surface energy integral in (51) is a specific instance of the class of integrals proposed in the literature for two-fold anisotropy [26, 27, 33]

 ES(α)=G0cw∫Ω(w(α)ℓ+ℓB∇α⋅∇α) (53)

where and are specialized to the AT-2 model. It is more interesting to compare the surface energy integral in (52) with the integral proposed in [33] for four-fold anisotropy (a simplified version of the model in [32]):

 ES(α)=G0cw∫Ω(w(α)ℓ+ℓ3B∇2α⋅∇2α)