An anisotropic phase-field model based on the equivalent crack surface energy density at finite strain

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 369 (2020) 113202 www.elsevier.com/locate/cma

An anisotropic phase-field model based on the equivalent crack surface energy density at finite strain Bo Yin, Michael Kaliske ∗ Institute for Structural Analysis, TU Dresden, 01187, Germany Received 29 October 2019; received in revised form 29 May 2020; accepted 2 June 2020 Available online xxxx

Abstract The promising phase-field method has been intensively studied and demonstrated to provide reliable predictions of structural failure. Motivated by experimental findings, the fracture patterns within several engineering materials are characterized as anisotropic, and are governed by the elastic and the fracture properties simultaneously. In this work, a physically anisotropic phase-field model is developed at finite strain in the framework of the second order phase-field theory. By using only one phasefield variable to evaluate the fracture status of the homogenized matrix and fiber materials simultaneously, an equivalent crack surface energy density function is established. Furthermore, the constitutive relations of the coupled problem are consistently derived by a straightforward variational principle, and are subsequently implemented into the context of the Finite Element Method. To demonstrate the capabilities of this approach, representative examples are studied and discussed. Consequently, corresponding findings and potential perspectives close this paper. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Anisotropic fracture; Phase-field modeling; Composite material; Finite deformation

1. Introduction Fracture plays a crucial role in failure mechanisms of materials and the reliable prediction of fracture evolution is of great importance and necessity. Numerous engineering materials, e.g. fiber reinforced materials, are designed to exhibit strong direction-dependent elastic response, which means they show anisotropic behavior. A number of pioneering works have studied the materials’ anisotropic behavior and developed several classical numerical models depending upon experimental observations, see [1–4]. Incorporating material failure mechanisms, the anisotropic damage evolution regarding engineering application is generally summarized in [5]. In particular, a representative study of anisotropic damage is theoretically developed within the framework of finite deformation in [6] and a subsequent extension to a gradient-enhanced model is formulated in [7]. Meanwhile, a couple of contributions studied anisotropic failure in the field of biomechanics, e.g. see [8–10]. In additional to anisotropic continuum damage modeling, brittle fracture theories have also been developed for anisotropic media. Motivated by experimental investigation, G RIFFITH [11] outlines the classical brittle fracture theory that crack formation irreversibly dissipates a specific amount of elastic strain energy. The amount of energy ∗ Corresponding author.

E-mail address: [email protected] (M. Kaliske). https://doi.org/10.1016/j.cma.2020.113202 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

consumed to generate a unit crack surface is defined as the fracture toughness Gc , also known as the critical energy release rate. As a material constant, the fracture toughness Gc is supposed to be direction-independent for isotropic solids. Nevertheless, this quantity is also identified to be directional-dependent as well in several anisotropic materials. Incorporating classical G RIFFITH theory, the smeared representative of a crack by means of phase-field modeling is intensively studied and demonstrated to be able to appropriately approximate crack initiation, propagation as well as branching with complex patterns. The fundamentals of phase-field fracture are initially developed in [12,13] as a straightforward variational model by minimizing the total energy, and it is subsequently regularized by [14–17] to name a few, for brittle fracture simulation. Meanwhile, several approaches, see [18–21], are developed to obtain a reliable fracture prediction, including the trustworthy crack evolution and physical crack kinematics. In the sequel, brittle phase-field modeling has been extended to several specific problems, e.g. [22–25] for dynamic fracture studies, [26–29] for rate-dependent fracture analysis, [30–42] for ductile fracture simulation and [43–45] for fatigue failure prediction. As one of the important features for composite materials, the anisotropic fracture characteristics have also been incorporated into various phase-field models as well. The work in [46] considers soft biological tissues as elastically anisotropic, which is coupled to an isotropic phase-field evolution. Similarly, in the framework of an isotropic phase-field evolution, [47] defines the driving force consisting of the contributions of both the matrix material and the fibers, which eventually yields anisotropic crack patterns as well. Due to experimental findings in [48] that certain crystalline materials show orientation-dependent fracture toughness but an isotropic elastic response, the anisotropic phase-field evolution is also employed to simulate the failure of polycrystals, see e.g. [49,50]. In the sequel, both the anisotropic elastic response and the anisotropic phase-field evolution are coupled simultaneously, see [51,52]. Furthermore, the anisotropic phase-field model with multiple phase-field variables is also developed, where [53] degrades the stress components independently by different phase-field variables along the spatial directions and [54] includes an equivalent phase-field variable by applying an algorithm to different phase-field components. The aforementioned anisotropic phase-field models are conceptually based upon the framework of weak anisotropy, which is outlined in [15,55]. Generally, the crack surface energy density function is restricted to involve the first order gradient of the phase-field, which is known as the second order phase-field theory, see [56]. Nevertheless, the concept of strong anisotropic fracture is also accordingly proposed in [57]. By evaluating the experimental results in [58], a strong anisotropic phase-field model is formulated in [59] by extending the second order phase-field theory to the fourth order phase-field approach. Meanwhile, a cubic symmetric fourth order structural tensor is introduced to govern the direction-dependent fracture resistance, also see [52]. Based on this formulation, the crack energy density function consists of both the first order and the second order gradient of the phase-field. Due to the existence of higher order gradient terms in the governing equations obtained by a variational principle, see [52,60], the continuity requirement of the element formulation in the FE application has to increase accordingly, which naturally increases the implementation complexity. In the paper at hand, a physically motivated anisotropic phase-field formulation, which is conceptually based on the second order theory, is postulated to simulate structural failure in anisotropic media. Both the elastic response and the fracture features are characterized by anisotropic properties simultaneously. The phase-field constitutive law involves the fracture evolution with respect to both the matrix material and the fiber components using one homogenized phase-field parameter, and consequently ends up with an equivalent crack surface energy density function. In published models, see [15,49–52] just to name a few, similarities exist with respect to the crack surface energy density function, which consists of a second order structural tensor to govern the anisotropic fracture property. Nevertheless, the remaining phase-field term is different. This slight dissimilarity effectively resolves the issue of the β-dependent crack profile, especially for large value of the scaling factor β. Furthermore, by applying a straightforward and consistent variational algorithm, the constitutive governing equations for the coupled problem are obtained and eventually implemented into the FE framework. The framework of this paper is outlined as follows. In Section 2, the fundamental constitutive formulation of the composite material is introduced along with a derivation of the stress tensor and the consistent tangent. In Section 3, an equivalent crack energy density function for anisotropic phase-field modeling is physically postulated. The following section develops the governing equations for the coupled problem. To validate and evaluate this present phase-field model, several representative and demonstrative numerical simulations are outlined in Section 5. Subsequently, Section 6 summarizes the findings and closes the paper with potential perspectives.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

3

2. Constitutive law of fiber reinforced material 2.1. Kinematics Before covering the concept of anisotropic material description, the fundamental kinematic relations need to be introduced. Within a continuous solid Ω0 , let the position vectors X and x stand for the material point in the reference and current configuration, respectively. The motion during the deformation is denoted by ϕ, which is defined as x = ϕ (X, t). The deformation gradient F = ∇ X ϕ (X, t) is one of the basic quantities for the kinematic description, where ∇ X stands for the gradient operator with respect to the reference configuration. To describe the distinct volumetric and isochoric material behavior of e.g. elastomers, the deformation gradient F is further split into 1 ¯ F = J 3 F,

(1)

¯ describes the isochoric response and the JACOBIAN J is the determinant of F. As common strain measures where F for the elastic response, the right and left C AUCHY–G REEN tensors and their isochoric quantities are defined as C = F T · F,

b = F · FT

and

¯ T · F, ¯ C¯ = F

¯ ·F ¯ T, b¯ = F

(2)

respectively. 2.2. H ELMHOLTZ energy density function The typical H ELMHOLTZ strain energy density function for s sets of fiber reinforced materials reads s ( ) ∑ ( ) Φ0 = U (J ) + Φ M C¯ + ΦkF F, m0k ⊗ m0k

(3)

k=1

with respect to the reference configuration, where U (J ) represents the volumetric energy in a common form U (J ) = κ M (J − ln (J ) − 1) .

(4)

The model parameter κ is known as the bulk modulus of the matrix material. The term Φ denotes the isochoric strain energy density of the matrix material, which is characterized as isotropic. Several alternative H ELMHOLTZ energy density functionals are available to approximate the matrix material. In this paper, a simple N EO - HOOKEAN model [4] is chosen for the matrix material, reading M

M

) ( ) µM ( Φ M C¯ = I¯1 − 3 , 2

(5)

¯ where the material constant µ M is known as the shear modulus. The scalar I¯1 is the first invariant of C¯ or b, ( ) ( ) F ¯ ¯ ¯ i.e. I1 = tr C = tr b . The last energetic quantity Φk in Eq. (3) denotes the strain energy potential for the fiber reinforcement. It is worthy to note that in order to simplify the phase-field fracture modeling, debonding between the fiber and the matrix material is excluded from the scope of this work. Nevertheless, to address the fiber debonding issue, several promising approaches based on the interface element framework are introduced to simulate the fiber pull out failure, see e.g. [61–63]. The homogenized fiber energy density ΦkF is characterized to depend not only on the deformation gradient tensor F but also on the second order structural tensor m0k ⊗ m0k . The orientation of the fiber regarding the set k in the reference configuration is expressed by a unit vector m0k , and it can be obtained in the current configuration by a push-forward operation mk = F · m0k . A variety of functional forms can be used to describe the fiber energy density, e.g. see [4]. In this work, a classical form is taken into consideration for the constitutive description of the fiber reinforcement reading λkF ( )2 I4k + µkF I5k , (6) 2 where λkF and µkF are two material constants regarding the fiber properties. The fourth invariant I4k and the fifth invariant I5k are two non-negative scalars, which are defined by ΦkF =

I4k = E : m0k ⊗ m0k

and

I5k = (E · E) : m0k ⊗ m0k ,

(7)

4

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

respectively. These two invariants numerically describe the elastic properties of the fiber families. The G REEN – L AGRANGIAN strain tensor E is obtained by E = 12 (C − 1), where 1 represents the second order identity tensor. 2.3. Material stress and consistent tangent tensors Having the strain energy function of the composite material at hand, the constitutive relations can be derived by a straightforward manner, where a detailed derivation for fiber reinforcement at small and finite strains is found in [64]. In the sequel, the K IRCHHOFF stress tensor and the consistent tangent tensor regarding the current configuration are obtained by the push-forward operation, which read τ vol = J ∂ J U (J ) 1, ( ) 1 ¯ M M ¯ b − I1 1 , τ =µ 3 ( ) ( ) ( ) µkF F F τ k = λk I4k mk ⊗ mk + mk ⊗ (b − 1) · mk + (b − 1) · mk ⊗ mk . 2

(8)

( ) Cvol = J ∂ J U (J ) + J 2 ∂ J2 J U (J ) 1 ⊗ 1 − 2 J ∂ J U (J ) I, ) 2( M ¯ µ I1 P − τ M ⊗ 1 − 1 ⊗ τ M , CM = 3( ) ( ) F F F Ck = λk mk ⊗ mk ⊗ mk ⊗ mk + µk (mk ⊗ mk ) ⊙ b + b ⊙ (mk ⊗ mk ) ,

(9)

and

respectively. With the definition of K RONECKER delta δab , the two frequently used fourth order tensors Iabcd = 1 δ + δad δbc ) and Pabcd = Iabcd − 13 δab δcd are defined. Additionally, another fourth order tensor denoted by 2 (δac bd the symbol ⊙ is given in index notation as ) 1( Aac B bd + Aad B bc . (10) ( A ⊙ B)abcd = 2 Thus, the total stress and consistent tangent are obtained as τ = τ vol + τ

M

+

s ∑

τ kF

and

M

C = Cvol + C +

k=1

s ∑

CkF ,

(11)

k=1

respectively. 3. Anisotropic phase-field modeling 3.1. Fundamental isotropic phase-field topology Within a continuous solid domain, a sharp crack is numerically approximated by a diffusive phase-field distribution, which is illustrated in Fig. 1. A scalar order parameter, namely the phase-field variable d(X, t), is introduced to identify the material state, i.e. the sound material is denoted by d = 0 and the fully cracked state is represented by d = 1. Mathematically motivated by a one-dimensional bar with an infinite length L ∈ [−∞, +∞], which is assumed to be cracked at position X = 0, the closed form solution for a continuous phase-field according to [16] is approximated by an exponential function |X | d(X ) = exp(− ). (12) l This solution is appropriately and naturally bounded by (0 < d ≤ 1). Another important parameter, the length scale l, is employed to govern the width of the transition zone between fractured and sound state. Extending to a twoor three-dimensional framework, the second order functional of the crack surface density is defined as ) 1 ( 2 γl = d + l 2 |∇ X d|2 , (13) 2l

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

5

Fig. 1. Diffusive crack topology.

where ∇ X (∗) denotes the gradient operator with respect to the reference configuration. Highlighted in [56], the realistic crack length Γ is approximately obtained by the phase-field crack length Γl as long as l → 0, i.e. ∫ γl d V, (14) Γ ≈ lim Γl = lim l→0

l→0

Ω0

which is known as the classical Γ -convergence condition for fracture. The defined crack surface density function in Eq. (13) is commonly known as the AT2 model, which yields an exponentially shaped crack profile as described by Eq. (12). Another common alternative is the classical AT1 model, see [17] for a detailed comparison to the AT2 approach. Furthermore, the recent works [65,66] have discussed the accuracy and reliability of phase-field modeling regarding the length scale and mesh sensitivity, nevertheless, they are mainly restricted to small strain problems. Addressing fracture for finite deformation, the paper at hand largely simplifies the problem by using the conventional setup for FE application in [67] to satisfy h e ≤ l/2 by considering the AT2 model. 3.2. Equivalent γlE Q for anisotropic phase-field modeling The isotropic crack surface density functional is outlined above in Eq. (13). According to the isotropic phase-field theory regarding brittle fracture, the total energy dissipated in the solid to generate crack surfaces Φ f rac is expressed as ∫ Φ f rac ≈ Gc γl d V, (15) Ω0

where Gc is known as the critical energy release rate or the fracture toughness. It is defined as the critical value of the strain energy releases within an infinitesimal crack length, see [11]. For most brittle fracture problems, this quantity is characterized as a material constant. Nevertheless, as highlighted in [25] and [42], this quantity depends on for example the loading rate and hardening state during plastic deformation, respectively. In this context, in order to incorporate the effects of fiber fracture properties on the overall fracture pattern, the homogenized fracture toughness Gc is postulated to be decomposed into the matrix material part GcM and the reinforcing fiber part GcFk . Accordingly, the phase-field crack surface energy density function is decomposed into γlM and γlkF as well. It is important to notice that only one phase-field variable is employed to evaluate the fracture status of the homogenized matrix and fiber materials simultaneously, and that the decomposition of Gc and γl does not simply imply a multiplicative relation. However, the relation can be expressed by rewriting Eq. (15) as ∫ ∫ s ∑ GcFk γlkF d V. (16) Φ f rac = GcM γlM d V + Ω0

k=1

Ω0

On the one hand, the matrix material is assumed to behave isotropically with respect to both elastic and fracture responses. Therefore, the energy density function of the phase-field crack surface γlM takes the form of Eq. (13)

6

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

with an alternative representation ( ) 1 d 2 + l 2 1 : ∇X d ⊗ ∇X d . (17) γlM = 2l On the other hand, the fiber reinforcement naturally exhibits a direction dependent response regarding the elastic behavior. To understand the direction dependent fracture of the fiber reinforcement, one may consider a group of un-bonded fibers. Separating these fibers does not consume any external work, since there are no bonding effects between the fibers. However, breaking the fibers requires a certain level of external load along the fiber’s longitudinal direction. Thus, the anisotropic fracture of fiber reinforcement is physically motivated as well. A direction dependent energy density functional of the crack surface is proposed by incorporating a consistent form of Eq. (13) as ( ) 1 γlkF = d 2 + l 2 m0k ⊗ m0k : ∇ X d ⊗ ∇ X d , (18) 2l where the unit vector m0k represents the orientation of the fiber layout in the reference configuration. Hence, the direction-dependent energy density γlkF is a function depending on the phase-field d, the phase-field gradient ∇ X d and the orientation vector m0k of fibers. Having the definition of γlM and γlkF in Eqs. (17) and (18), respectively, Eq. (16) is further extended as (( ) ( ) ) ∫ s s ∑ ∑ GcFk GcFk 0 1 M 2 2 0 Φ f rac = Gc 1+ d +l 1+ m ⊗ mk : ∇ X d ⊗ ∇ X d d V. (19) GM GM k Ω0 2l k=1 c k=1 c By introducing the ratio β k and the structural tensor W as s ∑ GcF and W =1+ β k m0k ⊗ m0k . β k = Mk Gc k=1 respectively, Eq. (19) is consequently simplified as Φ f rac = GcM crack surface energy density function ) ) (( s ∑ 1 EQ k 2 2 γl = β d + l W : ∇X d ⊗ ∇X d . 1+ 2l k=1

(20) ∫ Ω0

γlE Q d V with the definition of the equivalent

(21)

Physically, the ratio β k has the range as [0, +∞). In the specific case, β k = 0 returns Eq. (21) to Eq. (15), which yields an isotropic phase-field evolution. Comprehensively, for any non-zero β k , γlE Q is generally characterized as a function of ( ) γlE Q → F d, ∇ X d, β k , m0k . (22) Furthermore, recalling several existing representative anisotropic phase-field models [15,49–52], which are postulated in the framework of the second order phase-field theory, a general form of the anisotropic surface density function is defined as ( ) 1 0 2 2 γl = d + l W : ∇X d ⊗ ∇X d . (23) 2l ( ) ∑ Compared to the present model, the factor of 1 + sk=1 β k in front of the term d 2 in Eq. (21) is not considered in Eq. (23). In particular, these existing models do not take the ratio of the fracture toughnesses between the matrix and the fiber into account, nevertheless, they introduce another factor ζ k to adjust the length scale parameter along the fiber direction m0k to obtain an anisotropic phase-field distribution. In Eq. (23), a similar form of the second order structural tensor is defined as s ∑ 0 W =1+ ζ k m0k ⊗ m0k . (24) k=1

One the main issues regarding this definition is that the phase-field distribution is investigated to be strongly depending on the scaling factor ζ k . Namely, a large value ζ k leads to a significantly diffusive phase-field crack profile. Nevertheless, this issue does not happen in the present model. To illustrate this fact, a representative demonstration test is performed by considering both the anisotropic surface energy density functions of Eqs. (21) and (23) in Section 5.1.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

7

4. Governing equations for coupled problem The constitutive evolution for phase-field modeling is developed based on the balance of energy. The total internal energy density consists of the effective strain energy and the fracture energy within an anisotropic solid domain Ω0 as ∫ ( ∫ ) Q= γlE Q d V, (25) g (d) Φ0+ + Φ0− d V + GcM Ω0

Ω0

where the function g (d) in Eq. (25) prescribes the degradation of the elastic strain energy to evolve the phase-field and satisfies the mathematical properties g (0) = 1,

∂d g (d) ≤ 0

g (1) = 0,

and

∂d g (1) = 0.

(26)

The first two conditions in Eq. (26) indicate the boundary limitation of the sound and broken states, respectively. The two remainder conditions in Eq. (26) prescribe a monotonic degradation during the fracture evolution procedure. Particularly, at a fully broken state, ∂d g (1) = 0 circumvents the elastic driving force in Eq. (25). Satisfying those properties, an early proposal by B OURDIN et al. in [14] gives a quadratic function g (d) = (1 − d)2 .

(27)

Several alternatives are studied as well by [20,60,68]. In this work, a sinusoidal form function (π ) g (d) = 1 − sin ·d (28) 2 is applied, where a more brittle fracture property is investigated compared to the quadratic function Eq. (27), see [25,26]. Following the proposal of A MOR et al. [18], a volumetric–deviatoric split is employed to decompose the total strain energy into positive and negative parts to identify a proper driving force for the phase-field evolution, where the positive and negative parts are defined as Φ0+

= H (J ) U (J ) + Φ

M

+

s ∑

ΦkF

and Φ0− = (1 − H (J )) U (J ) ,

(29)

k=1

respectively. The operator H (J ) is applied to distinguish the volume compression or expansion state, which reads { 1, if J ≥ 1, H (J ) = (30) 0, if J < 1. In addition to the volumetric–deviatoric split, several other alternatives are also available, see [19–21] to name a few. A variational principle is applied to derive the governing equations for the elastic response and the phase-field evolution. Based on the balance of energy, the relation DK + Pint = Pext Dt exists, where the kinetic energy K and the external power Pext in Eq. (31) are defined as ∫ ∫ ∫ ∫ ( ) 1 ˙ ˙ ˙ ˙ B · U dV + T· U dA − χ d˙ d˙ d V, K= ρ0 U · U d V and Pext = 2 Ω0 Ω0 ∂Ω0 Ω0

(31)

(32)

respectively. The material density, velocity, body force and surface traction in the reference configuration are ˙ B and T, respectively. The last power type quantity in Eq. (32) considers a pseudo viscous expressed by ρ0 , U, resistance against the phase-field evolution. The non-negative parameter χ, which is known as the mobility number, is defined to evaluate the amount of viscous resistance during the phase-field evolution, i.e. brittle fracture is recovered at χ = 0. Usually, by applying a sufficiently small value of χ, the results are obtained as quasi-brittle fracture and the difference compared to brittle fracture at χ = 0 is quite minor and negligible. Meanwhile, the numerical stability is significantly increased. Nevertheless, using a comparatively large χ increases the resistance of phase-field evolution, which leads to delayed fracture compared to quasi-brittle fracture (χ ≈ 0) or brittle fracture (χ = 0). For detailed discussions, it is referred to [69] for quasi-static analysis and to [22] for transient analysis.

8

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

The internal power Pint in Eq. (31) is expanded as DQ Pint = } ∫D t{( ) ( ) ˙ + ∂d g (d) Φ + + GcM γlE Q d˙ + GcM ∂∇ d γlE Q · ∇ X d˙ d V, = g (d) ∂ F Φ0+ + ∂ F Φ0− : F 0 X       Ω0 P+

(33)

P−

∂ F Φ0±

where the quantities P ± = are known as the first P IOLA –K IRCHHOFF stress tensors and the rate of ˙ = ∇ X U. ˙ Inserting Eqs. (32) and (33) into Eq. (31) and by further algorithmic deformation gradient is F manipulations, the governing equations for the coupled problem are consistently obtained. Namely, the evolution equation for the elastic response reads ( ) ( ) ρ0 U¨ − ∇ X · g (d) P + + P − − B = 0 in Ω0 and N · g (d) P + + P − = T at ∂ Ω0 (34) regarding the reference configuration and can also be expressed as ( ) ( ) ρ u¨ − ∇ x · g (d) τ + + τ − − b = 0 in Ω and n · g (d) τ + + τ − = t at ∂ Ω

(35)

with respect to the current configuration. The divergence operator regarding the reference configuration and the ¨ N, B and T current configuration are denoted by ∇ X · (∗) and ∇ x · (∗), respectively. Similarly, the quantities U, in Eq. (34) are replaced by u¨ , n, b and t in Eq. (35) with respect to the current configuration. Particularly, the K IRCHHOFF stress tensors are obtained by τ ± = P ± · F T . Furthermore, the governing equation for the anisotropic phase-field evolution yields (( ) ) s ( ) ∑ GcM + k 2 1+ β d − l ∇X · W · ∇X d + χ d˙ = 0 in Ω0 , ∂d g (d) Φ0 + l (36) k=1 ( ) W · ∇ X d · N = 0 at ∂ Ω0 . In order to avoid healing during phase-field evolution, several options to define the driving force are available. Following M IEHE’s proposal [19], Φ0+ in Eq. (36) is replaced by an internal variable H = Max Φ0+ (F, τ ) , τ ≤tn+1

(37)

where a monotonically increasing phase-field is numerically guaranteed. In this contribution, M IEHE’s approach [19] is employed for anisotropic fracture simulation at finite strains, yielding ) ) (( s ( ) ∑ GcM 2 k ∂d g (d) H + + χ d˙ = 0. (38) 1+ β d − l ∇X · W · ∇X d l k=1 Alternatively, another approach outlined by K UHN et al. [70] imposes D IRICHLET boundary conditions to the phase-field degree of freedom as soon as the material is assumed to be fully cracked, e.g. d ≥ 0.9999, which also yields numerical irreversibility appropriately. The detailed derivation and the numerical implementation into the FE framework are not presented in this context. An inhouse element subroutine is developed, which consists of the residual and the consistent stiffness tangent and is referred to the appendix in [26]. 5. Representative simulations In this section, several representative numerical examples are studied to demonstrate the capability of the present anisotropic phase-field model. In order to investigate the anisotropic fracture evolution more efficiently and to simplify parametric studies, all simulations consider only one set of fiber reinforcement. 5.1. Two-dimensional Mode I test In the first example, a classical Mode I test is performed and the simulation results are compared to those obtained by the existing model according to Eq. (23). A two-dimensional boundary value problem, a specimen with

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

9

Fig. 2. Geometry setup of the classical Mode I test. The initial notch can be represented by (a) a discrete crack and (b) a phase-field crack, alternatively. Comparison of the initial phase-field crack profile by the existing model in Eq. (23) and by the present model in Eq. (21) regarding different scaling factors ζ and β in (c) and (d), respectively.

an initial notch, is considered to be subjected to tension up to the final failure. The geometrical setup is shown in Fig. 2(a), where a discrete crack is prescribed vertically in the middle of the edge. Another efficient way to prescribe an existing crack is defining a phase-field crack by setting the phase-field variables along the crack path as d = 1. Nevertheless, it is important to notice that a line of nodes with d = 1 is insufficient to represent a crack, since a fully evolved crack is approximated by a column of elements with fully degraded stiffness. Therefore, the nodes attached to the non-stiffness elements are given the phase-field initial boundary condition as d = 1. It is referred to [71] for

10

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

a detailed discussion. Thus, an alternative setup of the initial phase-field crack is depicted in Fig. 2(b). Regarding the numerical model, the finite element discretization consists of 15625 uniform 4-node quadratic elements and the element size is h e = 0.8 mm. The matrix material is assumed to behave purely hyperelastically with the material constants κ M = 1.96 MPa and µ M = 0.2 MPa. Additionally, fiber reinforcement is also considered by the model parameters λ F = 2.02 MPa and µ F = 0.18 MPa, where the orientation is depicted in Fig. 2(a) and (b). With respect to the fracture properties, the fracture toughness of the matrix material is GcM = 5 e−3 J/mm2 and the length-scale parameter is given as l = 0.2 mm. A set of factors β = [1, 10, 100] is studied. Before performing the fracture simulation, the profiles of the initial phase-field crack are investigated and compared in Fig. 2(c) and (d) with regard to the surface density functions of Eqs. (21) and (23), respectively. By increasing the factors ζ , the crack profile becomes more diffusive along the prescribed fiber orientation according to the surface density function Eq. (23), which leads to an unexpected modification of the initial damage status of the specimen. Nevertheless, the present phase-field model by Eq. (21) does not display such diffusivity. It successfully avoids the β-dependent phase-field distribution and presents an ideal and physical initial crack even at a large value of β. Nevertheless, applying β → +∞ leads to unphysical negative phase-field values adjacent to the crack path due to numerical reasons, but the neglectable negative values are quite small and do not affect the numerical simulations. To conclude according to this fact, the initial material damage status does not vary by increasing the factor β. After investigating the initial phase-field crack profile, tension failure is evaluated in Fig. 3. For both the existing and the present phase-field models, an inclined crack propagation is obtained in this anisotropic medium. It also can be found that by increasing ζ and β, the crack orientation is closer to the fiber direction. Furthermore, the relations of the resultant load with respect to the displacement are evaluated for both models. Since the existing phase-field model (Eq. (23)) significantly changes the diffusive crack profile of the initially prescribed notch, Fig. 3(c) does not properly capture the realistic relationship between load and displacement. Nevertheless, by applying the present model (Eq. (21)), the simulation results successfully predict a reasonable load–displacement property, see Fig. 3(d). Furthermore, different loading states during the fracture evolution, which is simulated by the present model with β = 100, are shown in Fig. 4. A post-processing technique is applied to blank the phase-field values d ≥ 0.95 to obtain a vivid crack propagation process. 5.2. Three-dimensional out-plane tear test This example studies the out-plane tear with different fiber orientations to demonstrate the capability of the present phase-field evolution. Opposite to [42], which studies an isotropic viscoelastic polymer tear fracture, this example does not take any rate-dependent properties into account but includes anisotropic characteristics. A threedimensional boundary value problem is analyzed and the geometric setup for two different fiber orientations is depicted in Fig. 5. The finite element model is discretized by 16240 uniform 8-node brick elements and the element size is h e = 0.5 mm. The model parameters of the matrix hyperelastic material are given as κ M = 1.96 MPa, µ M = 0.133 MPa and the material constants for the fiber reinforcement are λ F = 0.016 MPa and µ F = 0.0013 MPa. Regarding the fracture properties, the energy release rate of the matrix material is GcM = 0.244 e−3 J/mm2 and the length-scale parameter is l = 1 mm. The ratio along the fiber direction is β = 100. Firstly, the simulation studies the fracture evolution when the fibers are vertically oriented, see Fig. 5(a). As expected, the crack initiates at the notch tip and propagates straight towards the bottom edge, see Fig. 6. During this procedure, the crack path is totally parallel to the fiber orientation and it shows exactly the same crack path compared to an isotropic phase-field modeling, which is not shown here. In the sequel, the fracture evolution is simulated when the fiber orientation is changed to be along the horizontal direction, see Fig. 5(b). Crack initiation still happens at the notch tip, nevertheless, it does not follow the isotropic crack path anymore. Instead, the crack tries to evolve approximately along the horizontal direction, see Fig. 7. The main reason for this is that the ratio factor β is given as a relatively large value, which indicates that the fiber’s strength is quite high. Hence, breaking the fiber requires a much larger deformation than that for breaking the matrix material. Nevertheless, before reaching the state of fiber breaking, the stored strain energy due to large deformation reaches the level that can break the matrix material by shear effects. Therefore, the eventual fracture presents an approximately horizontal crack path. This phenomenon can also be physically understood by the fact that splitting the fibers might be easier than breaking the fibers. These two simulations are showing the fracture approximately parallel to the fiber direction due to a relatively large fracture toughness for fiber material.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

11

Fig. 3. Comparison of the eventual crack profile and orientation by (a) the existing model and (b) the present model with different factors ζ, β = [1, 10, 100]. Meanwhile, the load–displacement relationships are also evaluated and compared in (c) and (d) regarding different anisotropic phase-field characteristics.

5.3. Three-dimensional multi-layer tension test This example further studies a composite material’s fracture with multiple material components. Generally, two types of materials are considered, see Fig. 8(a). The fiber directions are perpendicular to each other, i.e. 45◦ and −45◦ for layer A and B, respectively. In addition to the varying fiber direction, the factor β for these two layers are also different, e.g. βA = 8 and βB = 10. The other fundamental material constants are assumed to be the same for

12

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 4. Fracture evolution of the simulation by the present phase-field model with β = 100.

Fig. 5. Geometric setup of the out-plane tear regarding two fiber orientations.

both layers, which are given as κ M = 2800 MPa, µ M = 1200 MPa, λ F = 30 MPa and µ F = 10 MPa. The matrix material’s fracture toughness is GcM = 5 J/mm2 and the length scale parameter is 2 mm. In order to demonstrate the capabilities of the present anisotropic phase-field modeling applied to multi-layer composite material fracture, three tests with different geometric setup are investigated, which are outlined in Fig. 8(b)–(d). The finite element model is discretized by uniform 8-node brick elements with the edge length of 1 mm, which consists of 100750 elements for tests 1 and 2 and 201500 elements for test 3. Regarding the first geometric setup in Fig. 8(b), the structure is composed of materials A and B alternately along the height direction. The interface debonding between two materials is not considered in this simulation. Thus, the crack is supposed to propagate straight in one material and subsequently to change direction when the crack evolves into another material. Therefore, according to the prescription of the fiber orientation for layers A and B, an expected zig-zag crack is obtained, see Fig. 9. A similar example is performed in [52] as well. The second test studies a complex crack pattern in a two-layer specimen, see Fig. 8(c). Unlike the setup of the first test, two types of materials are bonded along the thickness direction. Hence, when the specimen is subjected to tension, the crack initiates at the discrete crack tip in both materials. In the sequel, the crack tries to evolve along the fiber direction within each material, which should yield a nearly perpendicular crack pattern at the inside and outside surfaces of the specimen. However, since two materials are fully bonded, a portion of resistance exists to restrict the crack developing independently. Therefore, only a small angle, which is less than 90◦ , is initially observed between these two crack patterns at the inside and outside surfaces. Then, the cracks within layers A and B merge at the interface of the two materials and a thorough crack along the thickness is formed, see Fig. 10(a). Due to the different fracture properties (βA ̸= βB ), the thorough crack evolves up to a certain state and, subsequently, the crack at the inside layer starts to govern the overall direction. Then, the crack at the outside layer changes its direction and follows the inside layer crack up to the consequent failure, as shown in Fig. 10(b)–(d).

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

13

Fig. 6. Phase-field crack evolution of the specimen with vertically embedded fibers.

The third test increases the complexity of the second test by taking two more layers into consideration. Thus, a “sandwich” model is presented by the setup in Fig. 8(d), which consists of two layers A and two layers B bonded successively along the thickness direction. Similar to the second test, the crack initiates and propagates with resistance from the neighboring layers, which leads to an S-shaped crack pattern when observing from the bottom viewpoint. Compared to the second test, the final crack path does not incline too much and stays approximately at the center part of the specimen. Nevertheless, the crack surface shows complex concave–convex patterns, see Fig. 11. In realistic engineering applications, a non-smooth crack surface is quite common regarding the fracture of composites, which can be appropriately captured by this present phase-field modeling. 5.4. Failure of a wooden specimen In this example, the fracture of a wood structure conducted by M URATA et al. [72] is studied by the present anisotropic phase-field modeling. The elastic response of wood exhibits orthotropic behavior along the longitudinal, radial and tangential directions. The specimen in [72] made of cherry wood is tested in a tangential–radial setup for the classical Mode I fracture, which does not focus on the influence of the longitudinal property on the fracture evolution. Therefore, the insight is only on the fracture evolution in the tangential–radial plane and the crack surface perpendicular to the longitudinal direction is out of the scope. The experimental setup and the numerical threedimensional boundary value problem are depicted in Fig. 12(a) and (b), respectively. The finite element model is discretized by 10324 8-node linear brick elements with an element size h e = 0.02 mm at the potentially cracked

14

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 7. Phase-field crack evolution of the specimen with horizontally embedded fibers.

region. According to [72], the YOUNG’s modulus along the tangential direction E t is identified to be larger than that along the radial direction E r . Therefore, by incorporating the present model, this anisotropic tangential–radial plane is assumed to be decomposed into an isotropic plane with a portion of resistance along the tangential direction. The material parameters are identified as κ iso = 44.5 MPa, µiso = 33.5 MPa, λt = 36.35 MPa and µt = 25.46 MPa. The fracture toughness for the assumed isotropic plane is Gciso = 0.15 J/mm2 and the length-scale parameter is l = 0.1 mm. In total, eight experiments are conducted in [72]. The experiments show the load–displacement curves, which are different from one another, e.g. either with respect to the peak loads or fracture displacements. One of the important reasons for this phenomenon is that the specimens are obtained from different parts of the wood and the fracture properties may vary from one specimen to another. Therefore, in order to capture the range of the load–displacement relations, several factors β are taken into account, e.g. β = [0, 0.5, 1, 1.5]. The phase-field evolution with β = 1 is shown in Fig. 13, where the crack is represented by an isosurface contour plot of d = 0.95. Furthermore, the experimental investigation in Fig. 14(a) shows that the crack propagates through the whole specimen straight crossing the tangential direction. The simulation result is shown in Fig. 14(b) appropriately capturing this property and the load–displacement relations are compared in Fig. 14(c) with good agreement, which demonstrates the capabilities of the present phase-field model to fit experiments.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

15

Fig. 8. Material properties specification and the geometry setup for 3 different tests.

Fig. 9. Zig-zag crack evolution in the specimen with the setup of Fig. 8(b). The phase-field crack surface is emphasized by plotting the isosurface at d = 0.95 in a transparent geometry.

16

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 10. Curved crack surface and path obtained in the two-layer specimen of Fig. 8(c). The phase-field crack surface is emphasized by plotting the isosurface at d = 0.95 in a transparent geometry. The figures above are viewed from the side and the figures below are from the bottom viewpoint. The crack initiates in (a), changes direction at the outside layer in (b) and follows the crack of the insider layer in (c) up to final failure in (d).

Fig. 11. Curved crack surface and path obtained in the four-layer specimen in Fig. 8(d). The phase-field crack surface is emphasized by plotting the isosurface at d = 0.95 in a transparent geometry. The figures above are viewed from the side and the figures below are from the bottom viewpoint.

Fig. 12. Geometrical setup.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

17

Fig. 13. Phase-field crack evolution regarding the setup of Fig. 12 is emphasized by plotting the isosurface at d = 0.95 in a transparent geometry.

Fig. 14. Comparison of experimental results from [72] and the simulation results for both (a) and (b) the crack path as well as (c) the load–displacement relations.

5.5. Airbag structure opening simulation The last example is performed to simulate the opening process of an airbag structure, which is based on experimental investigations according to our industrial partner B ENECKE -K ALIKO AG. The comprehensive fracture mechanism of the airbag simulation involves not only the elastomer tear fracture but also adhesion failure, which is experimentally observed between the elastomer layer and the fixed rigid bottom. The adhesive damage is achieved by classical cohesive interface elements, which are characterized by a traction–separation law outlined in [61–63]. Three important parameters, the unit separation work φ, the maximum separation Tmax and the contact stiffness K c , are identified to evaluate the constitutive response of the standard cohesive approach. The implementation of the cohesive element is out of the scope of this paper and is described in the Appendix for a general insight. The geometry setup and the boundary conditions of the airbag structure are shown in Fig. 15(a). According to the realistic setup, a predetermined crack is considered along the edge of the lid with a thickness of 0.6 mm. The opening rate of the airbag cover is measured as 20 m/s at the opening edge. In total, 13488 finite elements are utilized for the airbag model. Among them, 8992 8-node brick elements are used to simulate the elastomeric layer and about 2880 rigid elements are employed to represent the stiff board. The rest are cohesive elements and placed under the elastomer layer. Different from standard volume elements, the upper and lower surfaces of the cohesive element share the same coordinates initially, see Fig. 15(b). The cohesive behavior is identified by the parameters of φ = 0.01 J/mm2 , Tmax = 7 mm and K c = 200 MPa. A number of experimental observations indicate that the elastic property of the elastomer layer shows isotropic nonlinear elasticity, while the fracture path does not show isotropic characteristics. Therefore, in this work, the elastic property of the elastomer layer is assumed to be isotropic with the

18

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 15. Geometry setup and finite element discretization strategy.

parameters κ M = 1.96 MPa, µ M = 0.204 MPa and λ F = µ F = 0 MPa. Meanwhile, the parameters for phase-field evolution are given as l = 6 mm, GcM = 0.13 J/mm2 and β = 0 to simulate isotropic fracture as a reference. In the sequel, the anisotropic fracture is performed with the parameters of GcM = 0.013 J/mm2 and β = [5, 10] with the pseudo fiber direction −45◦ from the x axis in the x–y plane. First of all, the reference isotropic fracture simulation shows the crack evolution in Fig. 16. It can be found that the isotropic setup leads to an exactly symmetric crack pattern. Fracture initiates at the left and right corners simultaneously. In the sequel, two crack tips merge at the middle region of the lifting edge of the lid. Meanwhile, the other two crack tips also propagate towards the rotation edge and show a symmetrical phase-field distribution eventually. The following simulations employ anisotropic phase-field evolutions, where the crack propagation processes are shown in Figs. 17 and 18 regarding the factors β = 5 and 10, respectively. As the results indicate, fracture initiates earlier at the left corner than that at the right corner for both cases. The eventual crack at the right side follows the edge of the flapping board. Whereas, the crack at the left side does not propagate along the pre-crack notch, which consequently results in a curved phase-field path. Especially for the case β = 10, at the lifting edge of the flapping board, the crack from the left side changes its direction significantly to merge the crack from the right side, ending up with a sharp corner path. These characteristics obtained by anisotropic phasefield modeling show good agreement compared to the experimental observation from B ENECKE -K ALIKO AG, see Fig. 19. Taking two failure mechanisms into consideration increases the complexity of the numerical simulation to certain extent. Moreover, the phase-field evolution is also driven by different loading modes during the deformation process. Initially, the cracks initiating at the corners and merging along the lifting edge are mainly a result of tension failure (Mode I). In the sequel, the cracks propagation along two sides are governed by a combination of both tension and tearing (Mode I and Mode III). Furthermore, the numerical stability is also influenced by several largely distorted elements, which are mainly located at the crack tip. Therefore, avoiding such element distortion is of

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 16. Fracture evolution of the airbag opening by an isotropic phase-field modeling (β = 0).

Fig. 17. Fracture evolution of the airbag opening by the present phase-field modeling with the factor β = 5.

19

20

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

Fig. 18. Fracture evolution of the airbag opening by the present phase-field modeling with the factor β = 10.

Fig. 19. Comparisons of the eventual phase-field cracks with the experimental observation from B ENECKE -K ALIKO AG.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

21

great importance for the numerical simulation, e.g. properly reducing the load increment for each step or applying a small mobility number to enhance the numerical stability for the phase-field solution. 6. Conclusions Phase-field modeling has been proved to be a promising approach to simulate the G RIFFITH-type brittle fracture. In this paper, a physically motivated anisotropic phase-field evolution is postulated within the framework of a second order phase-field theory. The standard isotropic crack energy density function γl is appropriately extended to an equivalent crack energy density function γlE Q , which depends not only upon the spatial direction of fiber but also on the fracture toughness ratio between the fiber reinforcement and the matrix material. The proposed anisotropic phase-field model for large strain problems is developed by applying a straightforward variational principle, and, in the sequel, is implemented consistently into the framework of the finite element method. The present model accounts for the constitutive elastic anisotropy as well as the fracture toughness anisotropy simultaneously. The classical anisotropic elastic property is achieved by considering a direction-dependent energy quantity to the total strain energy potential function. Meanwhile, the anisotropic phase-field evolution is physically proposed by considering the homogenized fracture of the matrix material and fibers. As a result, the equivalent crack surface γlE Q is formulated, which consists of a structural tensor W multiplying ∇ X d ⊗ ∇ X d and a factor ( ) ∑senergy 1 + k=1 β k multiplying the phase-field term d 2 to govern the directional dependency. The kernel development of this approach is obtaining a β k -independent crack profile, which does not yield a wide and diffusive phase-field distribution with large values of β k . Nevertheless, this issue is not properly resolved by several existing alternatives in [15,49–52] yet. Beyond the theoretical work, several representative numerical examples are studied as well. The first demonstrative example is motivated to examine if the phase-field crack profile is β-dependent within an anisotropic media. It is found that by applying the existing models (Eq. (23)), the diffusivity of the approximated crack largely depends on the factor ζ , nonetheless, this issue does not happen to the present model (Eq. (21)). Meanwhile, the evaluation of the load–displacement curves also show a reasonable characteristics by the present model. Then, two three-dimensional tear examples are simulated to investigate crack propagation. The crack path is observed to approximately follow the fiber orientation when applying a large fiber strength. In the sequel, a multi-layer composite structure is studied by defining different characteristics in each layer. Due to the structural setup, complex crack paths including zigzag and non-smooth S-shape patterns are consequently observed. The subsequent example is motivated by [72] to simulate the fracture property of a wooden structure. According to the experimental evidence, the crack is identified to propagate across the tangential direction instead of the radial direction due to the setup, which is properly predicted by the numerical simulation. And the range of the experimentally obtained load–displacement data is also effectively captured by controlling the factor β. The last example increases the simulation complexity by considering anisotropic fracture and interface failure simultaneously. An airbag opening process is numerically performed and fracture in the elastomer layer is characterized as anisotropic. Compared to the experimental investigation, the simulated phase-field crack pattern shows an amazing agreement. The contribution at hand introduces a modified phase-field model for anisotropic fracture simulation, which shows promising capabilities, as aforementioned. Nevertheless, further parametric studies need to be carried out in order to apply it to more complex engineering tasks. Besides, this work is only formulated in the pure elastic framework and the extension to inelasticity is also the following priority. Furthermore, the challenging higher order phase-field theory is an interesting and open topic that can be further developed. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment The authors would like to acknowledge the financial support of ANSYS I NC ., Canonsburg, USA, the experimental data supplied by B ENECKE -K ALIKO AG, Hannover, Germany, and the technical support of the Centre for Information Services and High Performance Computing of TU Dresden for providing access to the Bull HPC-Cluster.

22

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

¯ and ∂B± Fig. 20. Illustration of cohesive element from the reference to current configurations, ¯t and ∂B± representing the quantities of T 0 in the current configuration, respectively.

Appendix Cohesive zone element fundamentals This section gives a general description on the failure mechanism of the cohesive zone element approach, since the finite element implementation is dissimilar to classical volume elements due to non initial thickness. Meanwhile, the domain integration is manipulated within two crack surfaces instead of that within a continuous volume. Hence, special algorithmic treatments are required for the numerical interpolation by incorporating the conventional isoparametric concept. The three-dimensional setup of the cohesive zone element is illustratively shown in Fig. 20. In order to keep a consistent description, the weak form with respect to the reference configuration neglecting inertia effects reads (∫ ) ∫ ∫ ∫ ∫ + − + + − − ¯ · δ U¯ d A¯ + ¯ · δ U¯ d A¯ T · δUd A − T T B · δUd V − P : δ Fd V − = 0, (39) ∂B0+

∂Ω0

Ω0

Ω0

∂B0−







F

where the first three terms are aforementioned and last two terms are known as the interface force F. Regarding ± ± the field variables and their virtual quantities U¯ and δ U¯ at the upper and lower surfaces (∂B0− and ∂B0+ ) of the interface element, the finite element interpolations read − U¯ =

4 ∑

N I (ξ, η, ζ = −1) U¯   

I

N I (ξ, η, ζ = −1) δ U¯   

I

I =1

−

δ U¯ =

4 ∑ I =1

and

+ U¯ =

I =5

I N−

I N−

8 ∑

and

+

δ U¯ =

I

N I (ξ, η, ζ = +1) U¯ ,   

8 ∑ I =5

I N+

(40) I

N I (ξ, η, ζ = +1) δ U¯ .    I N+

¯ ± with respect to the reference configuration in Eq. (39) exist The first P IOLA –K IRCHHOFF traction vectors T + − ¯ ¯ ¯ −T = T = T. Meanwhile, the areas for the upper and lower crack surfaces are equal to each other, i.e. ∂B0+ = ∂B0− = ∂B0 . As a result, the interface force F is written as ) ( 4 ∫ 8 ∑ ∑ I I ¯ ¯ F= T N −I δ U¯ − N +I δ U¯ d A. (41) ∂B0

I =1

I =5

¯ this work adopts the classical exponential functional due to the To approximate the crack surface traction T, advantages of continuous and differentiable properties, which reads ( )( ) + − ¯ = φ exp − Tmax T U¯ − U¯ . (42) 2 D D

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

23

The model parameter φ represents the work required to completely separate two crack surfaces and Tmax denotes the maximum separation. The characteristic opening length is numerically obtained by D = φ/ (Tmax · exp (1)). Furthermore, if a contact deformation occurs, the traction law in Eq. (42) cannot be used anymore. Instead, the crack surfaces are supposed to transfer the load resistance and the traction law is re-defined as ( )2 P ¯ ¯ T = −K c n, (43) D where K c is a material constant to represent the contact stiffness due to the numerical penetration P. The unit vector n¯ is normal to the lower surface at each integration point in the current configuration, which is obtained ( +by the−cross) − − ¯ product of the tangent vectors ∂ξ x and ∂η x . The contact penetration is approximated by P = U¯ − U¯ · n. Having the traction Eqs. (42) and (43) at hand, in the sequel, the straightforward linearization of Eq. (41) yields the consistent tangent quantity. It is particularly referred to [63] for the detailed derivation and algorithmic implementation. 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]

S.W. Tsai, H.T. Hahn, Introduction to Composite Materials, Technomic Publishing Company, Lancaster, 1980. C.T. Herakovich, Mechanics of Fibrous Composites, John Wiley & Sons, New York, 1998. R.M. Jones, Mechanics of Composite Materials, second ed., Taylor & Francis, Philadelphia, 1999. G.A. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, John Wiley & Sons, Chichester, 2000. J. Lemaitre, R. Desmorat, Engineering Damage Mechanics, Springer, Berlin, 2005. A. Menzel, P. Steinmann, A theoretical and computational framework for anisotropic continuum damage mechanics at large strains, Int. J. Solids Struct. 38 (2001) 9505–9523. T. Waffenschmidt, C. Polindara, A. Menzel, S. Blanco, A gradient-enhanced large-deformation continuum damage model for fibre-reinforced materials, Comput. Methods Appl. Mech. Engrg. 268 (2014) 801–842. B. Calvo, E. E. Peña, M. Marínez, M. Doblaré, An uncoupled directional damage model for fibred biological soft tissues. Formulation and computational aspects, Internat. J. Numer. Methods Engrg. 69 (2007) 2036–2057. V. Alastrué, M. Martínez, M. Doblaré, A. Menzel, Anisotropic micro-sphere-based finite elasticity applied to blood vessel modelling, J. Mech. Phys. Solids 57 (2009) 178–203. D. Balzani, S. Brinkhues, G.A. Holzapfel, Constitutive framework for the modeling of damage in collagenous soft tissues with application to arterial walls, Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 139–151. A.A. Griffith, The phenomena of rupture and flow in solids, Philos. Trans. R. Soc. Lond. Ser. A 221 (1921) 163–198. L. Ambrosio, V.M. Tortorelli, Approximation of functionals depending on jumps by elliptic functionals via convergence, Comm. Pure Appl. Math. 43 (1990) 999–1036. G.A. Francfort, J.J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (1998) 1319–1342. B. Bourdin, G.A. Francfort, J.J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48 (2000) 797–826. V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, J. Mech. Phys. Solids 57 (2009) 342–368. C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, Int. J. Numer. Methods Eng. 83 (2010) 1273–1311. K. Pham, H. Amor, J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, Int. J. Damage Mech. 20 (2011) 618–652. H. Amor, J.J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, J. Mech. Phys. Solids 57 (2009) 1209–1229. 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. Engrg. 199 (2010) 2765–2778. C. Steinke, M. Kaliske, A phase-field crack approximation approach based on directional stress decomposition, Comput. Mech. 63 (2019) 1019–1046. J. Storm, D. Supriatna, M. Kaliske, The concept of Representative Crack Elements (RCE) for phase-field fracture - Anisotropic elasticity and thermo-elasticity, Internat. J. Numer. Methods Engrg. 121 (2020) 779–805. C. Kuhn, R. Müller, A continuum phase field model for fracture, Eng. Fract. Mech. 77 (2010) 3625–3634. A. Schlüter, A. Willenbücher, C. Kuhn, R. Müller, Phase field approximation of dynamic brittle fracture, Comput. Mech. 54 (2014) 1141–1161. C. Steinke, K. Özenç, G. Chinaryan, M. Kaliske, A comparative study of the r-adaptive material force approach and the phase-field method in dynamic fracture, Int. J. Fract. 201 (2016) 97–118. B. Yin, C. Steinke, M. Kaliske, Formulation and implementation of strain rate dependent fracture toughness in context of the phase-field method, Internat. J. Numer. Methods Engrg. 121 (2020) 233–255. B. Yin, M. Kaliske, Fracture simulation of viscoelastic polymers by the phase-field method, Comput. Mech. 65 (2020) 293–309. R. Shen, H. Waisman, L. Guo, Fracture of viscoelastic solids modeled with a modified phase field, Comput. Methods Appl. Mech. Engrg. 346 (2019) 862–890.

24

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

[28] L.M. Schänzel, Phase Field Modeling of Fracture in Rubbery and Glassy Polymers at Finite Thermo-Viscoelastic Deformations (Ph.D. thesis), Universität Stuttgart, 2015. [29] P.J. Loew, B. Peters, L.A.A. Beex, Rate-dependent phase-field damage modeling of rubber and its experimental parameter identification, J. Mech. Phys. Solids 127 (2019) 266–294. [30] F.P. Duda, A. Ciarbonetti, P.J. Sánchez, A.E. Huespe, A phase-field/gradient damage model for brittle fracture in elastic–plastic solids, Int. J. Plast. 65 (2015) 269–296. [31] M. Ambati, T. Gerasimov, L. De Lorenzis, Phase-field modeling of ductile fracture, Comput. Mech. 55 (2015) 1017–1040. [32] C. Miehe, M. Hofacker, L. Schänzel, F. Aldakheel, Phase field modeling of fracture in multi-physics problems. part II. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic–plastic solids, Comput. Methods Appl. Mech. Engrg. 294 (2015) 486–522. [33] C. Kuhn, T. Noll, R. Müller, On phase field modeling of ductile fracture, Surv. Appl. Math. Mech. 39 (2016) 35–54. [34] M.J. Borden, T.J.R. Hughes, C.M. Landis, A. Anvari, I.J. Lee, A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects, Comput. Methods Appl. Mech. Engrg. 312 (2016) 130–166. [35] C. Miehe, S. Teichtmeister, F. Aldakheel, Phase-field modelling of ductile fracture: a variational gradient-extended plasticity-damage theory and its micromorphic regularization, Philos. Trans. R. Soc. A 374 (2016). [36] C. Miehe, F. Aldakheel, A. Raina, Phase field modeling of ductile fracture at finite strains. A variational gradient-extended plasticity-damage theory, Int. J. Plast. 84 (2016) 1–32. [37] C. Miehe, D. Kienle, F. Aldakheel, S. Teichtmeister, Phase field modeling of fracture in porous plasticity: A variational gradient-extended eulerian framework for the macroscopic analysis of ductile failure, Comput. Methods Appl. Mech. Engrg. 312 (2016) 3–50. [38] M. Ambati, R. Kruse, L. De Lorenzis, A phase-field model for ductile fracture at finite strains and its experimental verification, Comput. Mech. 57 (2016) 149–167. [39] R. Alessi, J.J. Marigo, C. Maurini, S. Vidoli, Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples, Int. J. Mech. Sci. 149 (2018) 559–576. [40] F. Aldakheel, P. Wriggers, C. Miehe, A modified Gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling, Comput. Mech. 62 (2018) 815–833. [41] D. Kienle, F. Aldakheel, M.A. Keip, A finite-strain phase-field approach to ductile failure of frictional materials, Int. J. Solids Struct. 172–173 (2019) 147–162. [42] B. Yin, M. Kaliske, A ductile phase-field model based on degrading the fracture toughness: theory and implementation at small strain, Comput. Methods Appl. Mech. Eng. 366 (2020) 113068. [43] R. Alessi, S. Vidoli, L. De Lorenzis, A phenomenological approach to fatigue with a variational phase-field model: The one-dimensional case, Eng. Fract. Mech. 190 (2018) 53–73. [44] P. Carrara, M. Ambati, R. Alessi, L. De Lorenzis, A framework to model the fatigue behavior of brittle materials based on a variational phase-field approach, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112731. [45] M. Seiler, T. Linse, P. Hantschke, M. Kästner, An efficient phase-field model for fatigue fracture in ductile materials, Eng. Fract. Mech. 224 (2020) 106807. [46] A. Raina, C. Miehe, A phase-field model for fracture in biological tissues, Biomech. Model. Mechanobiol. 15 (2015) 1–18. [47] O. Gültekin, H. Dal, G.A. Holzapfel, A phase-field approach to model fracture of arterial walls: theory and finite element analysis, Comput. Methods Appl. Mech. Engrg. 312 (2016) 542–566. [48] M.H.B. Nasseri, B. Mohanty, Fracture toughness anisotropy in granitic rocks, Int. J. Rock Mech. Min. Sci. 45 (2008) 167–193. [49] J.D. Clayton, J. Knap, A geometrically nonlinear phase field theory of brittle fracture, Int. J. Fract. 189 (2014) 139–148. [50] J.D. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169. [51] O. Gültekin, H. Dal, G.A. Holzapfel, Numerical aspects of anisotropic failure in soft biological tissues favor energy-based criteria: A rate-dependent anisotropic crack phase-field model, Comput. Methods Appl. Mech. Engrg. 331 (2018) 23–52. [52] S. Teichtmeister, D. Kienle, F. Aldakheel, M.A. Keip, Phase field modeling of fracture in anisotropic brittle solids, Int. J. Non-Linear Mech. 97 (2017) 1–21. [53] J. Bleyer, R. Alessi, Phase-field modeling of anisotropic brittle fracture including several damage mechanisms, Comput. Methods Appl. Mech. Engrg. 336 (2018) 213–236. [54] T.T. Nguyen, J. Réthoré, M.C. Baietto, Phase field modelling of anisotropic crack propagation, Eur. J. Mech. A 65 (2017) 279–288. [55] V. Hakim, A. Karma, Crack path prediction in anisotropic brittle materials, Phys. Rev. Lett. 95 (2005) 235501. [56] M.J. Borden, T.J.R. Hughes, C.M. Landis, C.V. Verhoosel, A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework, Comput. Methods Appl. Mech. Engrg. 273 (2014) 100–118. [57] S. Torabi, J. Lowengrub, A. Voigt, S. Wise, A new phase-field model for strongly anisotropic systems, Proc. R. Soc. A 465 (2009) 1337–1359. [58] A. Takei, B. Roman, J. Bico, E. Hamm, F. Melo, Forbidden directions for the fracture of thin anisotropic sheets: An analogy with the wulff plot, Phys. Rev. Lett. 110 (2013) 144301. [59] B. Li, C. Peco, D. Millán, I. Arias, M. Arroyo, Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy, Internat. J. Numer. Methods Engrg. 102 (2014) 711–727. [60] M.J. Borden, Isogeometric Analysis of Phase-Field Models for Dynamic Brittle and Ductile Fracture (Ph.D thesis), The University of Texas at Austin, 2012. [61] J.W. Foulk, D.H. Allen, K.L.E. Helms, Formulation of a three-dimensional cohesive zone model for application to a finite element algorithm, Comput. Methods Appl. Mech. Engrg. 183 (2000) 51–66.

B. Yin and M. Kaliske / Computer Methods in Applied Mechanics and Engineering 369 (2020) 113202

25

[62] G. Geissler, M. Kaliske, M. Nase, W. Grellmann, Peel process simulation of sealed polymeric film computational modelling of experimental results, Int. J. Comput.-Aided Eng. Softw. 24 (2007) 586–607. [63] I. Zreid, R. Fleischhauer, M. Kaliske, A thermomechanically coupled viscoelastic cohesive zone model at large deformation, Int. J. Solids Struct. 50 (2013) 4279–4291. [64] M. Kaliske, A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains, Comput. Methods Appl. Mech. Engrg. 185 (2000) 225–243. [65] X. Zhang, C. Vignes, S.W. Sloan, D.C. Sheng, Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale, Comput. Mech. 59 (2017) 737–752. [66] T. Mandal, V.P. Nguyen, J.Y. Wu, Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture, Eng. Fract. Mech. 217 (2019) 106532. [67] M. Hofacker, A Thermodynamically Consistent Phase Field Approach to Fracture (Ph.D. thesis), Universität Stuttgart, 2013. [68] C. Kuhn, A. Schlüter, R. Müller, On degradation functions in phase field fracture models, Comput. Mater. Sci. 108 (2015) 374–384. [69] C. Miehe, L.M. Schänzel, Phase field modeling of fracture in rubbery polymers. Part I: Finite elasticity coupled with brittle failure, J. Mech. Phys. Solids 65 (2014) 93–113. [70] C. Kuhn, Numerical and Analytical Investigation of a Phase Field Model for Fracture (Ph.D. thesis), Technischen Universität Kaiserslautern, 2013. [71] C. Steinke, K. Özenç, G. Chinaryan, M. Kaliske, A comparative study of the r-adaptive material force approach and the phase-field method in dynamic fracture, Int. J. Fract. 201 (2016) 97–118. [72] K. Murata, E.V. Bachtiar, P. Niemz, Determination of mode I and mode II fracture toughness of walnut and cherry in TR and RT crack propagation system by the Arcan test, Holzforschung 71 (2017) 985–990.