Formulation and implementation of strain rate‐dependent fracture toughness in context of the phase‐field method

The phase‐field approach is a promising technique for the realistic simulation of brittle fracture processes, both in quasi‐static and transient analysis. Considering fast loading, experimental evidence indicates a strong relationship between the rate of strain and the material's resistance against fracture, which can be considered by a dynamic increase factor for the strength of the material. The paper at hand presents a novel approach within the framework of phase‐field models for brittle fracture. A rate‐dependent fracture toughness is formulated as a function of the rate of crack driving strain components, which results in higher strength for faster loading. Beside the increased amount of energy necessary to evolve a crack at a high strain rate loading situation, the model incorporates quasi‐viscous stress‐type quantities that are not directly related to the formation of the crack and exist only in the phase‐field transition zone between broken and sound material. The governing strong form equations for a transient simulation are derived and the relevant information for an implementation of the model into a finite element code is outlined in detail. The performance of the model is demonstrated for static and dynamic benchmark simulations and for a comparison to experimental findings.


INTRODUCTION
Fracture is one of the most critical failure mechanisms in industrial materials and structural components. The reliable and realistic prediction of crack initiation and propagation is of great importance. Brittle crack formation is usually considered to be the most dangerous mechanism in structural engineering, since it usually occurs at a relatively small deformation in an abrupt and sudden manner. The underlying physical concept is outlined first by Griffith 1 that the separation of the material requires a specific amount of energy, which is irreversibly dissipated within the whole solid. As the source of this energy, Griffith proposed the strain energy due to its deformation. In the sequel, the formation of cracks is the result of the energy transferring from the strain potential energy to the crack surface energy. The amount of energy per unit crack surface necessary is postulated to be a constant, material-dependent parameter, often referred to as the critical energy release rate or fracture toughness  c .
A well-established concept is the rate-dependent strength, also known as the strain rate effect, 2 which is commonly used for the design of civil engineering structures. 3 A representative example of rate-dependent fracture is the findings of Fineberg et al. 4 Based on a microscopical review of the fracture surface of a PMMA specimen after fast crack propagation, they conclude a relationship between the energy dissipation per unit area of the crack and the crack tip propagation speed. In this case, it is rather a rate-dependent characteristic of the crack evolution, than one of the material's viscoelastic properties. Considering macroscopic cracks at brittle fracture, the amount of dissipated energy depends on the size of the crack surface and the fracture toughness. The introduction of rate-dependent fracture toughness  c is the natural choice to account for the phenomenon described.
In principle, the numerical approximation of the topology of cracks by the finite element method (FEM) is divided into two basic concepts, which are discrete and continuous approaches. The discrete approach renders each crack surface as an explicit edge of the finite element (FE) discretization. As the crack and the corresponding jump in the displacement field are actually no longer part of the continuum meshed by finite elements, they are approximated realistically in a straight forward manner. That is true for both the crack's topology and the crack's opening at tensile loads. Additional effort is necessary to prevent interpenetration of the crack surfaces at compressive loading, eg, the introduction of interface elements along the crack surface. 5,6 Furthermore, the fracture evolution requires a method to govern the crack's propagation, eg, determine the propagating direction and branching. A representative approach is the material force method, implemented and studied, eg, by Gurtin, 7 Kienzler and Herrmann, 8 and Maugin. 9,10 Detailed information on the implementation and application within the framework of the FEM are available. [11][12][13] A promising approach of the second type of concepts, a continuous approximation of a crack, is the phase-field method. An additional field parameter is introduced to describe the state of the material, ie, sound or broken. Beside the regularized and mesh independent description of the crack's topology, the phase-field method involves the governing equations for the fracture evolution as well as a continuous approximation of the jump in the displacement field. Initially introduced by Francfort and Marigo 14 as a variational formulation for brittle fracture, the regularization of the approach is outlined by Bourdin et al 15 and analyzed for quasi-static fracture by Hakim and Karma. 16 The regularization length l is the key parameter, as it governs the width of the transition zone between the sound and the cracked phase of the material. In the theory of Γ-convergence, the quality of the crack's approximation by a phase-field crack density functional is evaluated. [17][18][19][20][21] A mathematical proof of Γ-convergence regarding the classical second-order phase-field model is mentioned by Bourdin et al. 15 A similar proof is not available yet for a higher-order phase-field crack density functional. 22 Another key feature of the phase-field method for a realistic crack evolution defines the crack driving force. One of the most common approaches is inspired by Amor et al, 23 where the strain tensor is decomposed into a volumetric and a deviatoric component. Then, the crack driving force is defined depending on the volumetric expansion or compression. Another approach is proposed by Miehe et al, 24 where the strain tensor is spectrally decomposed and the tensile energy is identified for the driving force. A recent approach considers the crack's orientation to define the driving stress and the corresponding driving strain potential. 25 Very few phase-field approaches to account for rate effects are available in the literature. A common rate-dependent phase-field model considers a mobility parameter to govern the speed of the phase-field evolution. This approach is applied as a numerical stabilization technique to the static simulation of brittle fracture at small strains, 26,27 and it is extended to dynamic fracture 28 as well as finite strains. 29,30 Besides, fracture evolution within viscoelastic solids is studied. 31 Another type of history-dependent phase-field models is extended to ductile fracture, [32][33][34][35] where the crack evolution is influenced by a plastic zone. In the contribution at hand, a fundamentally new rate-dependent formulation is introduced. In contrast to existing approaches, a nonconstant fracture toughness is considered. Instead, the fracture toughness is specified as a function of the local strain rate. At small rates of strain, the initial, quasi-static fracture toughness is recovered. For increasing strain rates, the amount of energy necessary for crack formation is increased. The model is motivated by experimental findings that indicate a higher crack resistance for increased loading speed. 2,4 The consistent derivation of the approach and its application to numerical simulations are presented.
The structure of this paper is outlined as follows. In Section 2, a brief overview of the approximation of a crack topology by the phase-field model is given and the governing equations of the standard brittle phase-field model in a transient simulation are outlined. Subsequently, an alternative degradation function is proposed and the rate-dependent phase-field model is introduced by reformulating the fracture toughness  c . Finally, the numerical implementation within the FEM framework is presented in detail. Section 3 contains three numerical benchmark simulations in order to demonstrate the characteristics of the proposed model. The last example particularly investigates the delay of the branching and evaluates the crack propagation speed in a PMMA specimen properly. Finally, Section 4 summarizes the findings and closes the contribution.

PHASE-FIELD MODEL FOR QUASI-VISCOUS FRACTURE
In this section, the variational phase-field model for dynamic brittle fracture is outlined and the governing equations are derived based on the Hamiltonian principle. Then, the formulation of rate-dependent fracture toughness and the corresponding governing equations are pointed out in detail.

Continuous approximation of crack topology
Consider the solid domain Ω with the discrete crack Γ and its regularized approximation Γ l , see Figure 1. The regularization is based on a scalar field parameter d, that is introduced as an order parameter to distinguish the unbroken solid (d = 0) from the broken domain (d = 1). The transition zone between broken and unbroken region is continuous and governed by the regularization length parameter l. A one-dimensional solution of the phase-field at x = 0 is approximated by the function which is outlined by Miehe et al. 27 Its extension to two or three dimensions is based on the isotropic, second-order functional of the crack surface density where the spatial gradient operator is denoted by ∇. Mentioned by Borden et al, 22 Γ-convergence is proved for this functional, ie, the discrete crack surface Γ is obtained for the regularization length parameter l → 0, where Within the framework of FEM application, according to Hofacker, 36 a reasonable accuracy is only obtained as long as l is larger than twice the size of the cracked element h e .

Governing equations for a phase-field model considering transient brittle fracture
The governing equations of transient brittle fracture for a phase-field crack approximation originate from a Hamiltonian principle. The general procedure involves the Lagrangian energy density . 37  specifies the amount of total energy density at a material point for time t n and contains energetic contributions of inertia, strain potential, and fracture energy density. The variation functional generally reads where the unknown field parameters are the displacement vector u and the scalar phase-field d. The vector of surface tractions for the displacements and the phase-field are denoted by t * u and t * d , respectively. It is necessary to point out that the surface traction terms are not considered in order to simplify the derivation work. The Lagrangian energy density is a function of the unknown field variables u and d as well as their derivatives ∇u, ∇ . u, . d, ∇d, and ∇ . d. Accordingly, the strong forms of both equations are given in general forms as respectively. The Lagrangian energy density is defined as which involves three quantities, the kinetic energy density , the effective strain energy densitŷe, and the fracture energy density Γ . The kinetic energy density is related to the velocity . u and the material density , which reads The effective strain energy densitŷe is composed by the fracture driving part + e and the fracture nondriving part − e , which is expressed aŝe The fracture driving part is multiplied by the degradation function g(d), which is discussed in detail in the next section. In order to identify the crack driving force + e , two common approaches are available in the literature. The volumetric-deviatoric split by Amor et al 23 decomposes the strain tensor into volumetric and deviatoric components. Only at expansion state, the volumetric energy is considered to evolve fracture. Another split according to Miehe 24 decomposes the strain tensor in its eigenspace as where i and n i are the eigenvalues and eigenvectors of the strain tensor, respectively. Moreover, the operator ⟨•⟩ ± is defined by Accordingly, the positive and negative strain energy density functions read where and are two material parameters and tr(•) is the trace operator. The positive and negative stress tensors are obtained by ± = ± e . Further constitutive derivations are outlined by Miehe. 38 In this contribution, the spectral split is applied to evaluate the proposed phase-field evolution. The last quantity Γ in Equation (8) is the fracture surface energy, which is computed by Considering the aforementioned definitions of ,̂e, and Γ , the unknowns in the Lagrangian energy density in Equation (5) are reduced to   → fun( which leads to the strong forms in Equation (6) and Equation (7) as and respectively. Δ (•) is the Laplace operator and d denotes the first-order derivative of the degradation function with respect to the phase-field. 39 In order to avoid healing of the phase-field, the driving force + e in Equation (17) is replaced by the internal variable Thus, a monotonically increasing phase-field is numerically described. By this assumption, Equation (17) turns into An alternative approach to achieve an irreversible fracture evolution is outlined by Kuhn et al 40

Alternative degradation function
The degradation function g(d) is used as a scalar multiplier in Equation (10), which specifies the local amount of strain energy degradation due to the value of the phase-field at this point. An increasing of phase-field locally yields both an increasing of crack surface due to Equation (2) and a reduction of the effective strain energy density specified in Equation (10), ie, it mathematically describes the energy transformation from the strain energy potential to the crack surface. To define the g(d), the limit values for the broken and the unbroken state are necessarily considered. At the unbroken state, the material is perfectly sound and no strain energy is degraded yielding Nevertheless, at the broken state, the driving force is fully degraded, ie, Within this range (0 ≤ d ≤ 1), the function g(d) needs to decrease monotonically, which reads Finally, the maximum degradation at d = 1 must describe an energetic minimum to restrict the phase-field from evolving larger than 1, which is obtained by d is a multiplier to the driving force of the phase-field in Equation (19) and d 2 is the second-order derivative of g(d) with respect to the phase-field. Furthermore, the degradation function needs to be continuous and differentiable. The most common choice is the quadratic function which is proposed in the early work of Bourdin et al. 15 The cubic function and the quartic function are investigated by Kuhn et al. 39 They conclude that the use of Equation (25) or Equation (26) results in a more brittle characteristic of the fracture process compared to the use of Equation (24). The concept "more brittle" for linear elastic material is understood as the stress-strain relationship as linear as possible up to the critical stress state. Another common phase-field model, namely, the AT1 formulation, 42 can guarantee an ideally linear elastic response up to the elastic stress limit, which is not discussed in this work. However, using Equation (25) or Equation (26) also requires a perturbation of the phase-field in order to initiate the phase-field evolution when d = 0, since g(d) d | d=0 = 0 leads to the driving force  multiplying 0 in Equation (19) at unbroken state d = 0 initially. This problem does not exist for the quadratic degradation function Equation (24) due to In order to achieve more brittle fracture without the initial perturbation of the phase-field, the exponential function is studied, 25 where the parameter is introduced in order to control the shape of the degradation function. A comparison of the aforementioned degradation functions is shown in Figure 2.
In this contribution, a sinusoidal form of the degradation function is proposed as which does not require any numerical perturbation of the phase-field to initiate its evolution. Furthermore, a comparatively more brittle property can be investigated compared to the quadratic function in Equation (24), which is evaluated numerically in Section 3.1. Comparisons of Equation (28) and Equation (24) as well as their first and second-order derivatives are shown in Figure 3 to Figure 5. According to Equation (19), the first-order derivative d is a multiplier to the driving force . It should be noted that using the sinusoidal degradation function Equation (28) yields a smaller scaling factor for d ≤ 0.3 than the quadratic degradation function Equation (24), see Figure 4. Therefore, the evolution of the phase-field is delayed for d ≤ 0.3, which leads to less driving force dissipated for the phase-field evolution. Furthermore, the linearization of Equation (19) for the FE implementation also involves the second-order derivative 2 g(d) d 2 , see Figure 5.

FIGURE 4
Comparison of the first-order derivative of the quadratic and the sinusoidal degradation function with respect to the phase-field,

Rate-dependent fracture toughness
Experimental observations indicate that the instantaneous fracture toughness is not constant at different strain rates for brittle materials like concrete 2 or PMMA. 4 Another example of rate-dependent fracture is found in polymeric materials. 12 The underlying mechanism is that the strength of the polymer chains and the cross-links is different. When fast loading is applied, the external energy required to break the material is actually higher, because the chemical bonds of the chains themselves have to break. In contrast, for slow loading, the chain segments can rearrange uniformly and the external energy is distributed to the cross-links, which require a lower level of energy in order to be broken. Approaches considering rate-dependent fracture toughness as a function of crack propagation velocity . a are available as Here,  0 c , . a 0 , and are the quasi-static values for the fracture toughness, the reference crack propagation velocity and the amplification factor, respectively. The approach is implemented and discussed for discrete crack approximation. 12 Within the framework of continuous crack approximation by the phase-field method, the definition and measurement of the crack propagation velocity . a is not trivial. The main reason is, that the crack tip is not defined explicitly in the phase-field method. 43 In this study, the basic assumption is an increased fracture toughness depending on the rate of the strain tensor . , ie, a large strain rate results in a large fracture toughness. Furthermore, the rate-dependent fracture toughness additionally depends on the state of the material, ie, The nonnegative scalar is a scaling factor, which is used to specify the amount of rate dependency, eg, = 0 results in a constant fracture toughness, see Figure 6. The multiplication with the degradation function g(d) prevents an increase of the fracture toughness for those regions, where the phase-field has a large value, ie, the crack is already evolved. Here, the opening of the crack may result in high strain rates, which should not be considered to evaluate the fracture toughness at all. Rather, the increased fracture toughness is meaningful only for those regions, where a crack might propagate. The rate of strain tensor is a local quantity, which preserves the local formulation of the phase-field approach.
It is important to note that . evaluates the rate of the strain regardless of the strain state, ie, the change of both principal tensile and principal compressive strain may influence the fracture toughness. Considering the phase-field being only driven by tensile energy, it seems plausible to consider only the tensile strain rates, ie, Here, the rate-dependent fracture toughness depends on the rate of the tensile strain and the rate of the compressive strain is totally neglected. In the following, this will be referred to as a spectral rate-dependent fracture toughness. A detailed discussion of the numerical calculation of . + is part of Section 2.6.

Modified governing equations
The evaluation of the Hamiltonian principle for the mechanical degree of freedom, given in Equation (6), yields the strong form commonly known as the balance of linear momentum, see Equation (16). Here, the divergence of the stress tensor is involved. Considering a rate-dependent fracture toughness, an additional stress-type quantity arises, which requires a modification of the balance of linear momentum where ∇ · (•) is the divergence operator and the total stress is obtained. The energy density  =̂e + Γ (34) involves both the effective strain energy densitŷe and the fracture surface energy density Γ . In the following, the equations and derivations are based on the assumption of the spectral rate-dependent fracture toughness, which is consistent with the phase-field evolution. Based on the definition in Equation (31), a rate-dependent fracture surface energy density is introduced as being a function of d, ∇d, as well as . + . The temporal discretization of the rate of the positive strain tensor yields the where the time step size, the positive strain tensors at the current and the previous steps are denoted by Δt, + t n+1 , and + t n , respectively. Thus, its first-order derivative with respect to the current strain is given by where the fourth-order positive projection tensor ℙ + t n+1 is employed. Detailed information on the derivation and implementation of ℙ + t n+1 is outlined by Miehe. 38 It should be noted that ℙ + t n+1 is based on the principal directions of the strain tensor, which can change with respect to time. The total stress now explicitly reads Beside the effective stress eff , an additional quasi-viscous stress-type quantity h + contributes to the total stress tensor as well. It should be noted, that both + and h + are multiplied by the degradation function, ie, they are completely degraded for d = 1. The rate-dependent fracture toughness according to Equation (31) is a function of the phase-field d. Therefore, its derivative has to be considered within the strong form of the phase-field evolution specified in Equation (7) as well. The modified phase-field evolution equation reads Remark 1. Dissipation with a rate-dependent fracture toughness According to the second law of thermodynamics, the entropy of a closed system is constant or growing. 44,45 An FE simulation is expected to be a closed system, ie, the amount of internal energy at a specific time t n is equal to the amount of external energy, that has been introduced into the system by the boundary conditions from the beginning at t 0 up to time t n . In a transient simulation with the phase-field model, three components of internal energy are specified explicitly, ie, the kinetic energy, the strain energy potential, and the fracture energy, see Equations (9), (10), and (35), respectively. In a phase-field approach with a constant fracture toughness, the fracture energy is directly related to the amount of strain potential energy degraded and may be calculated as a postprocessing step at any time of the simulation. Nevertheless, given the rate-dependent fracture toughness, the fracture energy depends on the crack surface as well as on the current rate of strain during the formation of the crack. Therefore, the total fracture energy is a result of an accumulation over the whole process of crack evolution. In addition, there is a quasi-viscous stress-type quantity, which also affects the strain energy. A significant portion of this energy remains in the transition zone between broken and unbroken material (g(d) ≠ 0), which actually constitutes an additional amount of internal dissipation. The dissipation can be specified by the evaluation of the Clausius-Plank inequality for isothermal processes, ie, where .
 represents the internal dissipation rate.
.  denotes the internal energy rate that can be written as Therefore, the rate of internal dissipation is approximated within an infinitesimal time step by The computation of the accumulated total internal dissipation by is performed numerically as where  t n+1 and  t n are the amount of total internal dissipation at the current and at the previous step, respectively.

Numerical implementation
Within the isoparametric concept of the FEM, N I represents the shape function at node I and its spatial gradient is written as N I x k , where k = 1, 2, 3 indicates the spatial direction for a three-dimensional problem. The residual of the mechanical equilibrium is obtained in an index notation form as The linearization of the total stresses yields the consistent tangent modulus where and Therefore, the mechanical stiffness matrix reads Furthermore, the mass matrix is given by with a Kronecker symbol ij . The residual of the phase-field equilibrium is written as and the tangent for the phase-field reads For the monolithic solution of the problem, the coupling terms read and Special attention should be paid tō+ in Equation (56), as it is not equal to + at all conditions. Considering an unloading process, the current strain energy + is not equal to the maximum value  and  = 0 holds. Therefore,̄+ is defined bȳ+ Additional derivatives are given by 2  c The matrices of stiffness and mass read and respectively, where n represents the number of nodes in element e. Moreover, the total residual is formed by .
By means of assembling those matrices of all elements within the global system, the total stiffness and inertia matrices as well as the residuals are obtained.

Quasi-static Mode I fracture test
In a first example, a Mode I benchmark by a two-dimensional boundary value problem is performed, where the geometric setup is shown in Figure 7. The left and right edges of the specimen are constraint and a monotonically increasing displacement is applied horizontally at the right edge. The numerical model is discretized by 9018 four-node linear elements and a refinement is applied at the potential fracture zone, where the element size is h ≈ 2 mm. The model parameters are = 2.8 GPa, = 1.2 GPa,  0 c = 500 J∕m 2 , and l = 4 mm. The inertia effects are neglected by setting = 0 kg/m 3 . The effect of both the strain rate and the loading rate is studied by this example, investigating different scaling factors and different loading rates . u i . The solution scheme is based on a staggered method. The simulation is started with the simplified staggered solution up to d = 0.3, which only solves the mechanical evolution and the phase-field evolution once at each loading step. As soon as the phase-field exceeds d > 0.3, an enhanced staggered scheme is applied until the end of the simulation, which solves the two evolution equations iteratively up to both equilibriums reach converged states.
It is known that using a finite length scale parameter l for the AT2 phase-field model results in a nonlinear relation between load and displacement before the peak load. The deviation between the linear and the nonlinear behavior may be decreased by reducing the length scale, which also requires decreasing the element size in the fracture zone. At the same length scale, the approximation of brittle fracture is enhanced by the sine degradation function compared to the quadratic degradation function. For the purpose of illustration, a linear elastic reference is given in Figure 8. The quadratic degradation function results in a peak force of f = 888 kN at u = 0.468 mm. The reaction force of the linear elastic reference at the same displacement is f = 1008 kN, ie, 13.5% higher. By using the sine degradation function, fracture occurs at a later stage at a displacement of u = 0.482 mm with the peak force f = 948.3kN. The linear elastic reference at this displacement is f = 1038 kN, ie, the difference is only 9.5%. With respect to this aspect, an enhanced approximation of brittle fracture by the sine degradation function is obtained and it is applied to all simulations in the sequel.
The scaling factor governs the rate-dependent  c . Given a constant loading rate 2.5 mm∕s, different scaling factors = 10 −4 , = 10 −3 , = 10 −2 , and = 10 −1 are analyzed in comparison to a nearly rate-independent fracture with = 10 −10 . In general, increasing the scaling factor results in an increase in both the peak force and the displacement of the total fracture. Both aspects are based on the same two reasons. On the one hand, the increased fracture toughness results in a larger amount of energy required per unit length of the crack length, which eventually leads to a larger amount of loading to evolve cracks. On the other hand, the amount of quasi-viscous stress is increased with increased , too, and this leads to larger internal dissipation. The relation between load and displacement with respect to the scaling factor is shown in Figure 9 and the detailed data are given in Table 1. Another aspect observed is that an increased value of    modifies the profile of the phase-field approximation, see Figure 10, which results in widening of the transition zone. From the energetical point of view, this results in additional energy dissipated into the formation of the fracture surface, both due to the increased amount of crack length and the increased local fracture toughness during the evolution. As already mentioned, this is also a major reason for the increase of both the peak load and total displacement necessary for complete fracture.
The loading rate . u governs the rate-dependent  c as well. Given a constant scaling factor = 10 −3 , the specimen is subjected to different loading rates  Figure 11. In general, increasing the loading rate results in an increase for both the peak reaction force and the displacement for total fracture as well, which can be explained similarly by increasing the scaling factor . In order to demonstrate the impact of different loading rates on the speed of crack propagation, the crack path at a displacement of u = 0.5 mm is visualized in Figure 12 with different loading rates. In general, the faster the loading is applied, the shorter the crack is at a specific displacement, which is in agreement with the previous explanations. .
Furthermore, in order to investigate the mesh-dependent influence of the presented formulation, six different FE discretizations shown in Figure 13 are studied for identical material parameters and loading conditions. The number of elements increases from 4042 to 11 667 and the element size h at the potential cracked region reduces from 6 mm to 1.5 mm, respectively. The length scale parameter, the scaling factor, as well as the loading rate are identically applied to all the simulations, which are l = 12 mm, = 2 · 10 −4 , and . u = 2.5mm∕s, respectively. The relations of load and displacement are compared in Figure 14. Concentrating on the fracture evolution process, Figure 15A gives clear visualization of the numerical convergence. Additionally, Figure 15B compares the displacements at the peak loads regarding different discretized models, which also shows a numerical convergence by increasing the element number.
Some remarks regarding the convergence are discussed. The positive strain rate . + is singular, which leads to ( as a singular quantity as well. Especially in the case of crack forming, the strain rate . + largely increases due to the full degradation and so does the quantity . However, the current fracture toughness  c is returned back to its initial value  0 c , since the degradation function is g(d) = 0 at a cracked state. Therefore, the singularity of  c vanishes with the regularization of the degradation function. Meanwhile, the quasi-viscous stress-type quantity g(d) · h + vanishes as well. As a consequence, the singularity of the strain or strain rate does not really affect the numerical convergence.

Spallation of a bar
The second example is a transient analysis of a three-dimensional setup. The main purpose of the simulation is the investigation of the proposed model with a rate-dependent fracture toughness in combination with wave propagation phenomena. Special attention is paid to the energetical balance during the simulation. The specimen investigated is a bar with a square cross-section, which is shown in Figure 16A. The length of the bar is 20 m and the height and width of the quadratic cross-section are both 0.1 m. It is motivated by the well-known Hopkinson-Bar spallation test that is used in order to examine the fracture toughness under transient tensile loading.  The lateral movement of the cross-section is constraint. The longitudinal direction of the bar is free and the front surface is subjected to a surface load in order to introduce a compressive stress wave into the specimen. The surface load is defined by (t) =̂· (t), where the magnitude is given aŝ= 1 MPa. f(t) is a half-sinusoidal function specified in Equation (64) and visualized in Figure 16B. Therefore, the external energy is introduced only during the initial period 0 ≤ t ≤ 32 μs and the sum of internal energy of the whole system is expected to be constant for t > 32 μs The discretization of the setup is based on eight-node linear element with the size h e = 1 mm. The model parameters are = 8.9 GPa, = 13.3 GPa, = 2300 kg∕m 3 , l = 2 mm, and  0 c = 0.1 J∕m 2 . The Newmark time integration scheme 46 with the constants = 0.25 and = 0.5 is employed, which leads to computations without numerical dissipation and unconditionally stable. The time step Δt = 1.25 μs is employed to obtain a realistic simulation of the wave propagation. 43 A major aspect of the simulations in this example is the investigation of the global energy balance. Here, four contributions to the total energy are relevant, ie, the effective elastic strain energy Φ e , the kinetic energy K, the fracture energy Φ Γ , and the internal dissipation energy D. Each component is a result of the volume integral of its specific density over the whole domain Ω. The sum of K and Φ e is denoted as the elastic internal energy Π e . Γ represents the fracture energy and the additional internal dissipation energy density D is obtained numerically according to Equation (45). These two terms specify the total dissipated energy Π d .

Remark 2. Computation of Γ
The fracture energy density can be computed based on Equation (14) with a constant fracture toughness  c . However, regarding a rate-dependent fracture toughness, Equation (14) cannot be used anymore, since  c is changing locally at each time step. Therefore, the approximation is employed to compute the fracture energy density, where the increment of the fracture surface density from previous to current step is expressed by Δ l,t n+1 = ( l,t n+1 − l,t n ). As a first step, a simulation with no phase-field evolution is considered in order to demonstrate the process of wave propagation. This is achieved numerically by a very large fracture toughness of  0 c = 10 20 J∕m 2 . Therefore, both Γ and  remain zero due to no fracture evolution and only the kinetic and elastic strain energy are transferred into each other, see Figure 17.
The total energy E increases until t = 32 μs, since no further loading is introduced to the system after this moment. Both the kinetic and the strain energy remain at a constant level until t = 53 μs. At this point, the compressive stress wave propagates to the end of the specimen and reaches the rear surface, where a tensile stress wave is reflected. While the compressive wave is propagating (32 μs ≤ t ≤ 53 μs), equal amounts of  and̂e are present. During the reflection of the wave, there is a short moment, where the local deformation of the bar is completely transferred into kinetic energy, ie, K = E and Φ e = 0. When the reflected tensile wave propagates in the opposite direction, the internal energy is split into equal amounts of  and̂e again. For the next example, a fracture evolution with a constant fracture toughness ( = 0) is investigated. The energy balance during the simulation is shown in Figure 18. For the first period, the fracture energy Γ remains zero, because no crack evolves for a compressed state. Once the compressive wave is reflected, the specimen is at tensile state locally and the phase-field starts to evolve at t ≈ 66 μs, where the subsequent fracture pattern is shown in Figures 19E, 19F, and 19G. At the beginning, the wave propagates through the whole bar from the front to the rear side and back again. After the evolution of the phase-field crack, the bar is broken into two separated parts by the crack surface. Within each of these segments, stress waves propagate and are reflected in a similar manner. However, the maximum tensile strain energy within each of these bars is not large enough to evolve other new phase-field cracks. Due to that reason, the fracture energy Γ is constant after a single crack is formed. A small numerical error results in the total energy E increased by approximately 0.5% at t = 200 μs.
A rate-dependent fracture toughness ( > 0) is considered for the third part of this example. The fracture patterns obtained at t = 100 μs for different values of are visualized in Figure 20. As already discussed in the previous Mode I test, the phase-field profile is widened for increased values of . Furthermore, for = 10 −3 , the resistance against fracture is large enough to circumvent the evolution of a phase-field crack completely, ie, d < 0.5 in this case. The energy balance during the simulation (0 < t < 200 μs) for = 10 −5 , = 10 −4 , and = 10 −3 is shown in Figures 21,22, and 23, respectively. Detailed information on the split of the total energy into the four components Φ e , K, Φ Γ , and D at the end of the simulation is summarized in Table 2. Here, also the categories Π e and Π d are explicitly stated.
The simulation results indicate that an increase of the scaling factor strongly affects the distribution between Π e and Π d . Generally, increasing results in a decrease of Π e and an increase of Π d . Another interesting aspect observed in the diagrams, Figure 21 and Figure 22, is that the separation of the bar into two individual segments leads to a dominant movement (high K) and a rather small deformation (low Φ e ). This is a typical feature of comparable experimental spallation tests, where the rear segment is flying away at a specific velocity after being separated from the front part of the structure. Nevertheless, this is not the case for = 10 −3 , because no separation of the bar occurs (see Figures 20D and  22). In this case, D is substantially increasing, since the quasi-viscous stress h + is only present within the transition zone of the phase-field (0 < d < 1).

Branching test
In this example, the branching phenomenon during fracture evolution is investigated. The geometry and boundary condition setup is shown in Figure 24. The width and height of the specimen are 380 mm and 440 mm, respectively. An initial crack is present at the right side with a length of 4 mm. The upper and lower edges are fixed along the vertical direction and a monotonic loading is applied at both edges in order to introduce tensile loading at rate    47 The specimen is discretized by 65 195 four-node linear elements and the mesh is refined with a size of h = 0.5 mm at the region in front of the initial crack tip, see Figure 24. The initial fracture toughness is defined as  0 c = 500 J∕m 2 and a length scale of l = 2 mm is used. Rate-dependent fracture is investigated for different values of . Regarding the phase-field evolution process, three stages are considered. Initially, a static simulation is performed by a simplified staggered solution until the phase-field near the crack tip reaches a value of d = 0.1. At this point, the transient simulation is activated with an enhanced staggered solution until d = 0.8. Up to this point, the scaling factor is not considered at all. As soon as d ≥ 0.8, the rate-dependent fracture toughness is activated. Meanwhile, the transient simulation with a monolithic solution scheme is applied. Two main advantages of this approach are outlined. First of all, it ensures that the crack initiates at the same time for all values of , which makes it possible to compare and analyze the impact of on the branching characteristics in a straightforward manner. Secondly, a nonzero results in a significant delay of the phase-field evolution from d = 0.1 to d = 0.8, especially for larger values of . As the focus of this example is at the propagation and branching of the crack instead of its initiation, the above considerations improve both the computational efficiency and the scientific significance of this study. Figure 25 shows the final crack path in the fine mesh region with respect to the scaling factor . Generally, a single crack is propagating from the initial crack tip in horizontal direction. At a certain point, which depends on , the profile of the phase-field crack starts to widen and sequentially continues up to another point, where two independent branches are established. The simulation with the rate-independent fracture toughness ( = 0) shows an additional feature that some branches stop to propagate shortly after the branching and only a single crack propagates further. As observed in the previous examples, a larger increases the widening effect of the phase-field.
Regarding the branching effect, generally a larger crack-driving force results in a larger phase-field value, which results in a widened profile due to the gradient term within the fracture surface formulation (see Equation (2)). As soon as the phase-field is wide enough to establish two spatially separated crack tips around the widened region, both tips start to propagate independently from each other. The strong influence of the local stress state on the angle between the branches is discussed by Steinke et al 43 in detail. Nevertheless, increasing the fracture toughness by results in a smaller value of the phase-field locally for the same level of crack-driving force. Therefore, an increased value of leads to a longer distance a crack may propagate with a widened profile before the independent crack tips form and separate from each other. Meanwhile, the total number of branches is reduced. In detail, a single straight crack in horizontal direction is observed for a length of 30 mm, ≈ 73 mm, and ≈ 140 mm for = 0, = 2·10 −9 , and = 10 −8 , respectively. The simulation with = 0 results in four major branches with several additional microbranches. For = 2 · 10 −9 , again, four major branches are observed. However, the crack patterns are smooth and no additional microbranches are observed. Finally, the simulation with = 10 −8 only results in two major branches, which propagate very close next to each other until the end of the refined area. Another interesting finding, the angle obtained between these major branches decrease with increasing . For the rate-independent fracture, after the two branches separated, an angle of ≈ 38 • between the branches is observed, which remains constant until the next major branching occurs. In contrast, for the rate-dependent fracture, the branches seem to stick together with a very small branching angle, eg, ≈ 10 • for = 10 −8 .
A major aspect of the evaluation of the experimental data is the record of the crack propagation velocity. 47 However, for the phase-field approach, the definition and tracking of a crack tip is a nontrivial task. The most common approaches are based on the definition of an isocurve at a critical phase-field value, which is evaluated in order to obtain the crack tip and its propagation velocity. In the contribution at hand, the approach presented by Steinke et al 43 is employed and the velocity is part of the postprocessing of the simulation. Here, the velocity evaluated always considers the tip of the upper, farthest left branch.
The approach to obtain the crack propagation speed experimentally is described in detail by Fineberg et al 4 and Sharon et al. 47 Two fundamental findings are noted. First of all, the maximum speed is limited to a value of ≈ 60% of the Rayleigh wave speed in the specific material. This is demonstrated by records of the crack propagation speed with respect

FIGURE 27
Crack velocity with respect to the time to both time and length of the crack. Secondly, a straight macroscopic crack cutting the whole specimen horizontally has been observed. However, a microscopical investigation of the crack surface revealed a typical sequence of surface characteristics. From the onset of crack propagation, the crack tip speed is continuously increasing. At the low level of a crack propagating speed, a plane and flawless crack surface is investigated. At increased crack propagation speed, the crack surface becomes rougher, which actually increases the total crack surface. Thus, although the macroscopic crack pattern is straight, a comparatively large amount of energy is dissipated. In the experimental evaluation, a relation between the crack tip propagation speed and the effective surface energy is established, which can be described by the crack tip speed-dependent fracture toughness given in Equation (29). Figures 26 and 27 show the crack propagation velocity regarding crack length and time, respectively. The experimental results are compared to the simulation results with different values of the scaling factor . It has been already shown, that the maximum speed is captured very well by the rate-independent phase-field model, 43 which results in a crack tip speed jumping from zero to the limiting value suddenly. In contrast, the proposed rate-dependent fracture model achieves the smooth acceleration of the crack tip correctly. Therefore, all fundamental aspects observed experimentally (smooth acceleration of the crack tip, limiting crack propagation speed, larger energy dissipation for fast propagation) can be modeled with the approach presented here.

CONCLUSIONS
The phase-field framework is a powerful approach to model fracture by a continuous crack approximation for static and transient simulations. In this work, a new phase-field approach for rate-dependent fracture is proposed, analyzed, and discussed. In contrast to rate-dependent material behavior, the fracture evolution is characterized as rate-dependency. The fracture toughness is reformulated as a strain rate-dependent quantity, which is governed by a scaling factor . In addition, the degradation function is introduced to regularize the rate-dependent fracture toughness. Due to the definition of  c totally dependent on local quantities, the mathematical derivation and the numerical implementation for a transient simulation within the framework of the FEM is straightforward. Furthermore, an additional rate-dependent stress-type quantity is derived based on the proposed model, which exists only in the transition zone between sound and broken material. The second law of thermodynamics is discussed and an additional amount of energy is dissipated due to the rate-dependent effects.
Three numerical examples are investigated in order to demonstrate the characteristics and performance of the new model. The first example is a static simulation of a Mode I fracture test. First of all, the sinusoidal degradation function is evaluated to yield a more brittle fracture compared to the quadratic degradation function. Furthermore, a study on the impact of the scaling factor and the loading rate for Mode I fracture is performed, which concludes that larger and faster load generally leads to higher peak load and larger displacement for total fracture as well as a widened phase-field profile. Additionally, the mesh sensitivity of this formulation is studied by comparing six different discretizations with the same model parameters and the loading conditions, where a proper numerical convergence is obtained. The second example studies rate-dependent fracture in a transient simulation. A compressive wave is applied to a slender bar initially. After reflection at the rear end, a tensile wave is formed, which leads to the spallation of the bar. First of all, the simulation demonstrates the ability of the phase-field method to model spallation in principle and a balanced energetic evolution is identified. A parametric study on the scaling factor yields the conclusion that the amount of internal dissipation increases for increased values of . Finally, the third example is motivated by experimental observations on the propagation of fast cracks in PMMA. Here, the experimental results for both the crack path and the acceleration of the crack tip are recovered by the simulation results. Nevertheless, by phase-field simulations with rate-independent fracture, the limiting crack propagation velocity is obtained correctly but it cannot recover the smooth acceleration of the crack propagation velocity. Instead, the proposed rate-dependent phase-field evolution can fulfill this property.
The contribution at hand proposes a rate-dependent fracture evolution. From a numerical point of view, the crack propagation is stabilized and experimentally observed characteristics of dynamic crack evolution can be simulated. Further steps are the extension into the framework of finite deformation, the consideration of other approaches to split the strain energy density, and the application of the approach compared to other additional experimental results.