Quasi ‐ Static and Dynamic Crack Propagation by Phase Field Modeling: Comparison with Previous Results and Experimental Validation

: In this paper, experimental tensile tests for pre ‐ cracked high Carbon steel ‘C90’ speci ‐ mens were performed for quasi ‐ static and dynamic loading. High loading velocity a ﬀ ects the crack pa tt erns by preventing de ﬂ ection. On the other hand, an e ﬃ cient numerical tool based on the phase ﬁ eld model was developed and validated to predict bri tt le fracture trajectories. A staggered nu ‐ merical scheme was adopted to solve the displacement and damage ﬁ elds separately. Implemen ‐ tation e ﬃ ciency in initiating and propagating cracks, even from an undamaged microstructure, was proved. The e ﬀ ect of the critical fracture energy density G c on the crack path was tested; with smaller G c , the crack pa tt erns become more complex. In addition, the impact of loading velocities was examined, and earlier and faster crack formation and greater crack branching is observed with higher impact velocity. In this study, bidimensional plane stress cases were treated. The phase ﬁ eld model with hybrid formulation was able to predict crack pa tt ern and especially crack arrest and branching found in the literature. The developed model accurately determined the transition zone of the crack path topology that has been observed experimentally.


Introduction
Crack propagation in solids, whether brittle or ductile, can have devastating financial effects and put people's lives in peril.Consequently, both scientists and engineers are interested in looking for predictive description tools regarding this phenomenon.Griffith [1] proposed one of the first successful descriptions which, although functioning, is too cumbersome to measure and implement in engineering cases.Dynamic crack initiation and propagation cannot be predicted using this theory; however, this remains an active topic of research due to its complexity.
Practical engineering cases, with the exception of a few theoretical cases, are generally too complicated to be solved using analytical expressions.Consequently, numerical approaches are extremely important in fracture analysis.Both discrete and continuous methods can be used to model cracks.Although molecular dynamics [2,3] provides a powerful overview of the physics of fracture initiation, the length and time scales are very small for application in practical cases.Discrete element simulations [4] overcome the problem of scale.However, despite the advantages of discrete techniques, continuous description is the most frequently utilized approach.Using various enrichments, the finite elements approach was proven capable of modeling dynamic fracture propagation.The cohesive zone method [5][6][7] is focused on the separation of predetermined surfaces, supposing a finite treatment zone at the fracture's front.A finite displacement step is introduced at the crack front in 'the strong discontinuity approach' [8].The eXtended Finite Element Method (XFEM) [9] has been successfully used to study dynamic fracture propagation in elastic brittle [10] and ductile [11] materials.Gomez et al. [12] used the boundary element method (BEM) to predict crack propagation for fatigue fracture in mixed modes (I + II).Fractures are treated as discrete discontinuities in the methods mentioned above.Consequently, treating crack initiation or branching is complicated and frequently requires special criteria.
Another approach is provided by non-local damage models [13][14][15][16].The thick level set approach calculates the fracture trajectory on the basis of the damage progression at the crack point, and the discontinuities are explicitly taken into account.However, discontinuities are not directly employed in other models, such as the phase field method.Rather, the crack effect is defined by a smooth transition.The completely damaged regions, as well as the undamaged material, are linked by a unique damage variable.The non-local damage models applied to Griffith's theory are based on the variational formulation developed by Francfort and Marigo [17,18] and implemented by Bourdin et al. [19].Complex fracture processes in two-and three-dimensional [20] solids, such as crack initiation, propagation, and branching, can be processed automatically without any additional ad hoc strategies.
The monolithic scheme [21][22][23][24][25][26] and the staggered scheme [19,[27][28][29][30] are the two most important numerical methods adopted to simulate the mechanical problem and determine the fracture pattern.When the fracture propagates closely to an initial crack, monolithic algorithms eventually converge.However, these algorithms struggle to converge in the case of unstable crack initiation and propagation.The non-convexity of global energy and the jumps in the solution are the causes of this problem [31].To avoid this, a variety of modifications [24,31,32] have been proposed.The staggered scheme is derived from the alternate minimization algorithm developed by Bourdin et al. [19].Two well-posed convex problems are solved separately.Even in unstable cases, this iterative method was shown to be very robust [29].Therefore, staggered implementation was considered as the right approach to use in this paper.The only drawback of this technique is that it necessitates a relatively small time step [30]; thus, a high computational time is needed.To overcome this issue, several research works were adopted, such as detecting unstable propagation and reducing the external load variation [22,23,30,[33][34][35], deactivating damage variables until an equivalent stress reaches a critical value σc [36], optimizing the assembly matrices [37], etc.
Decomposing the stored elastic energy into two parts, one associated with a compressive strain and the other with a tensile one, is a common approach in solving the problem of crack propagation.The decomposition of the elastic energy density in volumetric-deviatoric parts [28] and the spectral decomposition of the strain tensor [21] are two examples of decomposition methods that have been investigated.Hybrid formulation [38][39][40] is another technique applied to deal with asymmetric tensile and compressive behavior.Although it is variationally inconsistent, the resulting model is substantially more efficient computationally than the spectral decomposition models [41].
The remainder of this paper is organized as follows: In Section 2, experimental tensile tests on pre-cracked plates with a center hole were performed for different distances between the pre-crack and the hole axis to show their influence on the crack patterns for quasi-static and dynamic loading.Section 3 details the governing equations of the phase field model for brittle fracture based on the variational approach.Section 4 introduces four sets of illustrative numerical examples to prove the capability of the suggested model in depicting the quasi-static and dynamic fracture propagation.Then, the effectiveness of the hybrid formulation in predicting the crack patterns achieved experimentally is proven.The most relevant conclusions are drawn in Section 5 before suggesting some future potential perspectives.

Quasi-Static and Dynamic Tensile Experiments
The available experimental results to confirm the crack trajectory predicted by the phase field formulations are actually scarce in the literature.In this section, we experimentally investigated the crack path for a notched plate under quasi-static and dynamic tension.In the experimental tests, we used a brittle material, i.e., non-alloy steel with a high Carbon percentage C90 in the as-quenched state [37].

Specimen Description
Since there are no norms for testing brittle fracture propagation, we selected a simple and particular specimen geometry conceived to obtain various crack trajectories during the experiments.The simple geometry is set as a rectangular plate with a central hole and an initial pre-crack.Different distances between the pre-crack and the central hole axis, noted as h, were considered (Figure 1).To avoid slipping in the contact regions of the specimen jaw during the tensile tests, holes were machined.

Experimental Design
In these tests an axial displacement is applied, with a tensile machine for quasi-static cases and fast jack for dynamic ones (INSTRON VHS 65/20), with a maximum displacement velocity of 20 m/s for the cross bar.During tensile tests, some specimens were painted with a black and white speckle pattern to allow for a better visualization of crack propagation.A high-speed Manta camera G-1236 from Allied Vision Company-Germany (12.4 million pixels, full resolution: 4112  3008, lens Mount: C) was used to film the tests.Vimba X 6.1 software was used for data acquisition and image processing.The tests consist of applying an imposed displacement load to the top face of the specimen.The crack trajectories observed for the test specimens were obtained after three repeatable tests for each case.

Experimental Results and Analysis
In the quasi-static tensile tests, the controlled displacement load was imposed at a velocity of 0.01 mm/min.Figure 2 shows the experimental crack trajectories.For h = 5 and 6 mm, the crack starts to propagate from the pre-crack to the hole, then a new crack was initiated and propagated from the hole to the right edge, and for h = 9 mm the crack propagated horizontally from the pre-crack to the right edge.
In the dynamic tests, the loading velocity is 1 m/s as shown in Figure 3, and for h = 5 and 9 mm the crack trajectories are the same as those achieved in the quasi-static case.However, for h = 6 mm, the crack propagated horizontally without passing through the hole.
In this test, we can see the effect of loading velocity on the crack trajectory.In fact, the inertial effect can change the crack path.The crack has no time to deviate to the hole and propagate directly toward the opposite edge.

Brittle Fracture Phase Field Model
Let us consider an arbitrary solid Ω, with ∂Ω the external boundary and Γ the internal discontinuity boundary, which represents an initial fracture.The displacement field in the solid is designated by u.A body force acting on the solid b and a surface force f on the boundary ∂ΩF are considered.
According to Griffith's theory, Bourdin et al. [34] proposed a variational approach for fracture problems.This indicates that the critical fracture energy density Gc, also known as the critical energy release rate, is equivalent to the energy needed to generate a fracture surface per unit area.The total potential energy Etot(u, Γ) is formulated as follows: where, ( )    u is the linear strain tensor and dS is the surface element.
The elastic energy density ψ0(ε) in the isotropic linear elasticity is provided by: where µ and λ are Lamé constants.Furthermore, according to the variational approach [34], the crack initiates, propagates and branches when the minimum potential is reached and the irreversible condition is fulfilled.The irreversible condition means that a material cannot regain its initial state after a crack formation.

Diffusive Crack Energy Approximation Using Phase Filed Method
The brittle fracture variational approach was applied with success by Miehe et al. [21,27], Bourdin et al. [34] and Borden et al. [42].To represent the fracture surface Γ (Figure 4a), a phase field α(x, t)  [0,1] is defined (see also Figure 4b).α(x, t) represents a diffusive crack topology.α = 1 indicates the cracked zone and α = 0 indicates the uncracked zone.Therefore, the solid crack surface density γ is defined by Miehe et al. [27] as: where η is the length scale parameter that regulates the width of the transition region of the phase field (Figure 4b).The crack region width increases with the scale parameter η (see Figure 4b).The fracture energy is calculated using Equation (3):

Governing Equations
According to Francfort and Marigo [17], the fracture surface energy is obtained from the elastic energy, which determines the damage variable evolution.The elastic energy is decomposed into tensile and compressive components to ensure that only the tensile component of such density is governing the phase field variable evolution.This assumption cannot predict the crack initiation in the case of uniaxial specimen compression as a concrete compression test.It also discards the influence of the crack roughness [12] on the propagation phenomena.Here, the decomposition method described by Miehe et al. [27] was adopted and the decomposition of the strain tensor is as follows: where ɛ− and ɛ+ are the compressive and tensile strain tensors, respectively.The operators . and . are the positive and the negative parts, respectively.εa is the principal strain and na is its direction vector.
Using the spectral decomposition of the strain tensor, the elastic energy density is illustrated by the following equation: To model the stiffness reduction due to damage phase propagation, the following equation was used [42]: where k is a numerical parameter that prevents numerical singularities if the phase field variable α tends towards 1.Furthermore, this parameter has to be positive and very small compared to 1.
In dynamic fracture propagation, the kinetic energy of the body is added to the energy balance.The kinetic energy of a solid is given by the following equation: where is the velocity vector and  is the mass density of the material.
The total Lagrange energy functional is given by the following equation: The functional L's variation may be calculated, and its first variation must be 0, which leads us to the governing Equation (10) given by Ambati et al. [38]: where is the acceleration vector, σ the stress tensor and I is the unit tensor R d×d .
The law of damaged material behavior given by the thermodynamic of continuums using Equations ( 6) and ( 7) is: To reduce the computational time, we adopted the hybrid formulation of stress tensor proposed by Ambati et al. [38] and given by the following equation: During compression or unloading, the irreversibility condition is required to ensure a monotonically increasing phase field variable.Consequently, a strain-history field H(x, t) [20,27], defined as shown below, was introduced: To solve this problem, the finite element method in the case of plane problems was used.The adopted elements are a three-node triangular model with a linear interpolation.This program was implemented in the Matlab 2015 software, using the Flowcharts of Figures 5 and 6

Numerical Examples
This section deals with quasi-static and dynamic computational examples to prove the proposed model's ability to reproduce the experimental crack paths.

Quasi-Static Tensile Test of Notched Specimens
This example was considered to demonstrate the capability of the staggered scheme with the hybrid formulation to predict the experimental crack trajectory obtained in the previous section for different 'h' values between the pre-crack and the hole under quasi-static loading.The boundary conditions are a fixed edge and an imposed axial dis-placement on the opposite edge (Figure 7).The domain was meshed with triangular elements with a characteristic size s = 0.1 mm.The material properties are E = 210 GPa, Gc = 0.03 kN/mm, and ν = 0.3.The length scale parameter was chosen as η = 2  s = 0.2 mm.
The crack paths given by the numerical simulations are illustrated in Figure 8.The obtained crack paths are in agreement with those found experimentally (see Figure 2).In fact, in the case of h = 5 and 6 mm, the cracks changed direction and pass through the hole and, in the case of h = 9 mm, the cracks propagated directly to the opposite edge.

Compression Test of a Regularly Dispersed Pores Microstructure
In this example, a microstructure built of plaster with holes spaced at regular intervals was considered.Nguyen et al. and Sammis and Ashby [43,44] have investigated this example numerically and experimentally.Figure 9 shows the compression test geometry of a plate with 23 holes of diameter d = 0.2 mm.The specimen's dimensions are 2 mm and 3 mm.
Experimentally, Sammis and Ashby [44] found that the cracks initiate and propagate vertically from the holes.Nguyen et al. [43] obtained, in addition to the vertical cracks observed in the tests, an initiation of cracks from the left and right sides of the holes, which is unphysical, due to high compression stress in these regions.Figure 10 shows the crack path simulation given by our models.This demonstrates the ability of the devel-oped model to initiate cracks from a microstructure that is not damaged, with a correctly predicted crack trajectory after nucleation.

Dynamic Tensile Test of Notched Specimens
This example was dealt with to validate the crack paths obtained experimentally under dynamic loading with the same geometry proposed in Section 2.1, while the elements' size and material properties are as in the example studied in Section 4.1.1.
The numerically obtained crack paths are illustrated in Figure 11.The transition of the crack pattern observed experimentally for h = 6 mm was correctly achieved with our model.The obtained crack trajectories are in agreement with those found experimentally (see Figure 3).Compared to the quasi-static case, a change in the crack path was observed in the case of 6 mm distance between the pre-crack and the central hole axis.

Dynamic Shear Loading of Kalthoff Experiment
The treated example for a dynamic fracture is related to a doubly notched specimen subjected to an impact load, commonly referred to as the Kalthoff-Winkler experiment [45], which has already been investigated by many other researchers [46][47][48][49][50][51][52][53][54].The specimen geometry is shown in Figure 12, and the impact loading is generated by a projectile.The symmetry conditions are adopted to reduce the computational cost.
The following velocity was applied with t0 = 1 µs and v0 = 16.5 m/s.
The material parameters were as stated in [42]: E = 190 GPa,  = 8 × 10 −6 kg/mm 3 , ν = 0.3 and Gc = 22.2 × 10 −3 kN/mm.The length scale parameter η was fixed as 0.4 mm and the time step t was defined as 0.04 µs.The fracture began to propagate at 24 µs.The angle between the horizontal axis was 68°, which corresponded well with the experimental results of Kalthof-Winkler [48] and the numerical results of Zhou et al. [54].
Then, the effect of Gc, on the crack path, was tested.The phase field variable distributions at 90 µs for Gc = 5 × 10 −3 , Gc = 10 × 10 −3 , Gc = 22.2 × 10 −3 and Gc = 30 × 10 −3 kN/mm are displayed in Figure 13.For smaller Gc, more complicated crack schemes were found.When Gc = 5 × 10 −3 and 10 × 10 −3 kN/mm, crack branching occurs.The reflected wave caused secondary cracks in the lower right-hand corner of the specimen, as shown by Zhou et al. [54].However, only one crack was observed for a greater Gc.The principal crack pattern, which began from the point of the pre-existing crack, was unaffected by changes in Gc.
Eventually, higher impact velocities of v0 = 25, 40 and 60 m/s were examined.The crack paths are given in Figure 14.As the impact velocity increased, earlier and faster crack formation and greater crack branching were observed.
As a result, the developed model can easily capture multiple dynamic crack branching with appropriately refined meshes.When several small fractures exist along a major dominant one, this result is intriguing for the study of dynamic crack instability problems [55].It is worth mentioning that the tension-compression slit model [56] affects the development of the crack patterns.In this part, we introduced a third example of the numerical simulations of a dynamic fracture test performed by Grégoire [57].The specimen was directly impacted by a fast jack with an impact velocity of 10 m/s and the compressive loading was transformed by the circular hole into tensile loading in the vicinity of the initial notch.The specimen geometry is shown in Figure 15.The geometry was discretized by 60,000 triangular elements whose characteristic size is s = 0.2 mm.The material properties were those adopted in [57]: E = 2.4 GPa,  = 1.18 × 10 −6 kg/mm 3 , ν = 0.42 and Gc = 7 × 10 −4 kN/mm.The length scale parameter was η = 0.4 mm and the time step was t = 0.5 × 10 −4 ms. Figure 16 illustrates the dynamic crack propagation obtained by simulation.
The crack propagation began from the initial crack with an angle of −76°, which corresponded well to the experimental result of Grégoire [57].In Figure 17, we superposed the experimental result and the numerical result obtained with phase field and XFEM method treated by Grégoire [57].Clearly, the phase field method has a greater ability to reproduce the experimental observations.In fact, the crack path obtained by the developed program and the experimental are practically similar.

Crack Arrest Due to Compressive Stress
The last example considered in this sub-section was the 'one crack two holes' studied by Grégoire [57] and Li et al. [56].The aim was to better understand the effect of holes presence at the cracks; initiation and propagation, as treated by Fageehi and Alshoaibi [58].The specimen was directly impacted by a fast jack with an impact velocity of 10 m/s.The specimen geometry is shown in Figure 18.The material properties and the mesh size are the same as in the preceding example.Figure 19 shows the dynamic crack propagation obtained for different time steps.
The crack propagated from the initial crack and its arrest occurs at 0.13 ms in the vicinity of the high-compression area located under the right circular hole.The predicted crack arrested in the same way as the experimental one.Then, at the boundary of the right circular hole, a secondary crack nucleation was observed at 0.135 ms under high tension.This result is very similar to that obtained by Li et al. [56].However, this phenomenon was not observed experimentally by Grégoire [57].
For length scale η = 0.2 mm, no nucleation of secondary crack was detected at t = 0.135 ms (see Figure 20).This result highlights, again, the importance of internal length as a material parameter.

Conclusions
Several 2D plane stress benchmark examples for quasi-static and dynamic crack propagations were investigated.The effects of critical energy release rate, velocity and length scale parameter on the results were also checked.It can therefore be concluded from this study that, with the hybrid formulation, the phase field method provides a hopeful tool for managing fracture problems and complicated crack path topologies with correct crack patterns and satisfactory accuracy.We conducted tensile tests on pre-cracked specimens with a central hole made in quenched steel C90 due to the limited number of experimental results that demonstrated the dynamic effect on the crack path in brittle materials.The trajectory of the crack is a function of the distance between the pre-crack and the hole axis as well as the loading types: quasi-static or dynamic tests.The numerical analysis revealed the same phenomena.In 2D problems, the implementation of this model is simpler and requires less time to calculate than the case of anisotropic formulation.The results obtained with the phase field method using a hybrid formulation are extremely useful.In fact, this method is capable of capturing various dynamic fracture phenomena including fragmentation, crack arrest, crack branching, and multiple branching.
As future perspectives, we will first focus on the enhancement of the energy expression in order to be able to predict the crack initiation in the case of a uniaxial compression test and will take into account the crack roughness on the propagation phenomena.We will also extend this work to study three-dimensional brittle crack propagation and ductile fracture, taking into account the plastic behavior of the material.
The drawback of this technique is that it requires a small mesh size to ensure correct crack path estimation.Thus, a great deal of computing time and a very powerful computer to carry out the simulations in the case of real engineering structures are needed.The hybrid formulation should be enhanced to predict the crack initiation in the case of a uniaxial compression test and to take into account the effect of crack roughness on propagation phenomena.

Figure 1 .
Figure 1.Geometry of the pre-cracked specimen.
detailed in Appendixes A and B.

Figure 5 .
Figure 5. Flowchart of the staggered scheme with hybrid formulation for dynamic fractures.

Figure 6 .
Figure 6.Flowchart for the phase field solution for a given increment using a staggered scheme at each integration point.

Figure 7 .Figure 8 .
Figure 7. Geometry of the specimen and adopted boundary conditions: tensile test.

Figure 9 .Figure 10 .
Figure 9. Geometry of the specimen and adopted boundary conditions: compression test.

Figure 15 .
Figure 15.Geometry of the specimen and adopted boundary conditions: dynamic compression test.

Figure 17 .
Figure 17.Superposition of the experimental and numerical crack paths.

FF
are the inertial, internal, and external forces for the displacement field and int  F and ext  are the internal and external force terms of the phase field.C is the elasticity tensor.Additionally, the mass and stiffness matrices can be written as