Analysis of Hertzian indentation fracture using a phase field approach

The formation and further evolution of cracks caused by the compression of a stiff indenter onto the surface of an initially defect‐free brittle solid is a fascinating problem of fracture mechanics. Its prediction, however, is still a challenging task since crack nucleation is caused by a rather weak stress concentration in the contact near‐field. The present contribution focuses on phase field simulation of indentation fracture, including crack formation at some a priori unknown location outside of the contact region and the subsequent formation of a cone crack. While the phase field method, at first glance, appears to be a promising tool to simulate the current problem we elaborate critical issues and discuss essential modifications. Finally, the indentation fracture process is simulated showing the effect of varying indenter radii on crack initiation and the influence of Poisson's ratio on the angle of cone crack propagation in good agreement with experimental findings and other theoretical studies.


Introduction
Since the early investigations by H. Hertz at the end of the nineteenth century, indentation fracture, i.e. the formation and further evolution of cracks caused by the compression of a stiff indenter onto the initially defect-free surface of a brittle solid, has been subject of numerous experimental and theoretical studies, e.g. [1][2][3][4][5][6]. In the present work, we regard a flat and initially defect-free surface of an isotropic linear elastic body (E, ν) subjected to indentation loading by a rigid cylindrical punch of radius a (Fig. 1). The body is taken large enough relative to a (e.g. half-space) that axisymmetry can be assumed. Loading is imposed through a prescribed, monotonically increasing vertical displacement d in the contact region r ≤ a, which is assumed frictionless. Hertzian fracture proceeds essentially in two stages: First, a ring crack of radius r 0 > a and small extension z 0 from the free surface spontaneously forms somewhat outside the contact region at a critical load. Afterwards and with increasing load, it turns into a cone-shaped crack that grows in a stable manner at an approximately constant angle ϕ which depends on Poisson's ratio ν of the material (Fig. 1). A particular challenge in the theoretical analysis of indentation fracture is the initial ring crack formation from the defect-free surface which cannot be described by classical fracture mechanics concepts. By applying the concept of finite fracture mechanics (FFM) in [6], it was possible to predict the critical load, the radius and the length (extension from free surface) of the initially forming ring crack from the analysis without any additional assumptions. Through considering strength σ t and toughness G c as independent material parameters in a hybrid fracture criterion, FFM shows some analogy to phase field models of fracture. The present work aims at a holistic description of the entire fracture process including fracture initiation and subsequent cone crack growth. Phase field modeling has in recent years become a popular computational technique and appears to be a promising tool to analyze the current fracture problem. For comparison with the experimental study by Mouginot and Maugis [1], we consider material parameters of borosilicate glass with E = 80 GPa, ν = 0.22, G c = 9 J/m 2 in the following. Since the tensile strength of glass is a rather unclear material parameter we investigate strengths in the range σ t = 115 . . . 150 MPa. For indenter radii in the range of a = 0.5 . . . 2.5 mm as considered in [1] a length (depth from free surface, see Fig. 1) of the spontaneously formed initial ring crack of z 0 ≈ 0.04 mm was estimated from the FFM study [6].

Suitable phase field approaches
In phase field approaches to fracture a sharp crack Γ is approximated in a smooth manner by an additional scalar field, the phase field S(x, t). The phase field describes the current state of the material and varies continuously between S = 1 for intact and S = 0 for broken material. Thus, the surface energy of a discrete crack is represented by a volume integral containing the surface density γ ℓ = (1 − S) 2 /(2ℓ) + ℓ |∇S| 2 /2, where the internal length parameter ℓ characterizes the transition zone between the intact and broken state. The coupled set of governing equations are the momentum balance equation ρü − div σ = 0 with u =ū on ∂Ω u and σ · n =t on ∂Ω t (1) and the phase field evolution equation Since the above equations do not take the irreversible character of fracture into account, an additional irreversibility con-straintṠ ≤ 0 is applied where the phase field becomes smaller than a critical value S ≤ S 0 ≈ 0 in order to describe a "crack-like" behavior only when the phase field is close to the broken state.

Critical issues with regard to simulation of indentation fracture
A severe problem of the classical phase field approaches is the unphysical occurrence of fracture under compression and the inadequate treatment of the non-interpenetration condition at crack closure [7]. Since in the present problem large regions, especially below the indenter, are loaded by compression this becomes particularly relevant here. Widespread isotropic decompositions of the elastic energy density like the volumetric-deviatoric split [8] or the spectral decomposition [9] do not account for the crack orientation. Thus, such approaches generally violate crack boundary conditions and lead to a spurious phase field evolution [10]. They also completely fail to reproduce indentation fracture. While the energy consumption of a smeared crack is approximated quite well for a small, but non-vanishing internal length ℓ the crack tip consumes additional energy which is not accounted for in the method [11]. Thus, the deviation of the approximated surface energy from the theoretical value becomes significant unless the crack is many times longer than ℓ, as illustrated in [12]. Obviously, this requirement cannot be fulfilled during crack formation in absence of a pre-crack. Regarding the indentation fracture problem, the surface energy of the initiation process should be properly approximated. Therefore, the internal length ℓ has to be chosen sufficiently small compared to the length z 0 of the "spontaneously" formed ring crack (cf. Fig. 1). This imposes restrictions on the degradation function g(S) which couples the field equations (1) and (2) by describing the change of elastic properties, e.g. σ := g(S) C 0 : ε ε ε act + C 0 : ε ε ε pas , and the energy transfer from the elastic field to the phase field, e.g. ψ drive := ∂g(S)/∂S ψ el (ε ε ε act ). The internal length ℓ is directly related to the tensile strength σ t in conjunction with the degradation function g(S). For given material parameters (E, G c , σ t ) the internal length ℓ can only be reduced by choosing an appropriate degradation function.

Constitutive choice of stiffness degradation and phase field driving energies
As a pragmatic approach to avoid unphysical phase field behavior under compression we chose the degradation of stiffness and the elastic contribution to phase field evolution independently according to σ := g(S) C 0 : ε ε ε and ψ drive := ∂g(S) ∂S Since in the present situation of indentation fracture under monotonic loading cracks once they have formed remain open, i.e. traction-free, it is feasible here to omit any split in the calculation of stress (1) and consider a full stiffness degradation (3) 1 . By contrast, a tension-compression split is required in the phase field evolution equation (2) to avoid the phase field to be driven by elastic energy caused by compression. Since brittle materials tend to fail normal to the maximum principal stress, the phase field S should be driven by positive principal stresses rather than strains. In the present case, it is assumed that all positive effective stresses contribute to phase field evolution (3) 2 , withσ i being the eigenvalues of the effective stress tensorσ = C 0 : ε ε ε. It should be noted that, despite the separate choice of stiffness degradation and crack driving energy, the structure of the governing equations still corresponds to the variational approach [7] by using the derivative ∂g/∂S as scaling factor in (3) 2 to keep the energy transfer from the elastic field to the phase field consistent. As pointed out in Sect. 2.1, the internal length ℓ depends on the choice of the degradation function g(S) and has to be sufficiently small compared to z 0 . Referring to the discussion in [12], the degradation function g = (4 − a s ) S 3 + (a s − 3) S 4 of fourth-order with a s = 1 × 10 −4 is an optimal compromise to reduce the internal length sufficiently, but also guarantees that the phase field approaches the broken state. For fixed E and G c the internal length thus varies monotonically as a function of the tensile strength as worked out in [13] with values in the range According to the present boundary value problem of indentation fracture an axisymmetric FE formulation with bilinear shape functions is employed. Since the physically motivated separate choice of stiffness degradation (3) 1 and crack driving energy (3) 2 would cause unsymmetric entries in the tangent matrix of the fully coupled set of equations, we use an alternate minimization procedure where displacements u and phase field S are computed separately in an alternating manner as proposed by [7]. To achieve reliable solutions small time/load steps are necessary and, in addition, multiple iterations between the solutions of both equation systems are needed in each load step until the change of the phase field between consecutive iterations is below a small tolerance. Numerical studies are carried out on a quadratic domain of size 100a×100a for different indenter radii a = 0.5 . . . 2.5 mm. Boundary conditions of zero normal displacement and zero shear stress are applied at the bottom and on the symmetry axis (r = 0), whereas the outer boundary (r = 100a) is traction-free. Loading is specified in terms of a monotonically increasing indenter displacement d. Corresponding load increments are adaptively controlled to accurately capture fracture initiation. Since the FFM study [6] indicates that the quasi-static energy balance is satisfied in the mean during formation of the finite ring crack and the subsequent cone crack propagation proceeds in stable manner dynamic terms are neglected. The element size is chosen to be h e ≈ ℓ/4 in the zone of interest. For any finite discretization the mechanical field inevitably influences the phase field solution at a crack and thus leads to an overestimation of the approximated surface energy [12,14]. This effect is almost eliminated by a simultaneous re-adjustment of the fracture toughness and the internal length parameter according to The phase field method with the extensive modifications discussed in Section 2.2 is in the following applied to the simulation of the entire process of indentation fracture. This includes crack formation at some a priori unknown location outside the contact region as well as the subsequent formation of a cone crack as discussed in Section 1.

Effect of indenter radius on crack initiation
After the onset of loading, we first recognize a small amount of diffuse phase field evolution before the phase field abruptly localizes well outside the contact region and features a ring crack, which extends almost normal from the free surface. The contour plots in Fig. 2 show ring cracks caused by a small indenter of a = 0.5 mm and a large indenter of a = 2.5 mm for varying material strengths somewhat after initiation. They show that (for a constant material strength) a larger indenter provokes ring crack formation at a relatively smaller radial position r 0 /a. With increasing material strength σ t the radial crack position r 0 decreases. Furthermore, both figures indicate a decreasing crack length z 0 with increasing material strength σ t , which is in general agreement with FFM [6]. The "spontaneous" character of this ring crack formation is clearly indicated by the abrupt increase of the fracture surface energy ∆E s in Fig. 3. While the radius of crack initiation r 0 can be easily identified (vertical line in Fig. 2), the vertical extension z 0 of the spontaneously formed ring crack is much less clear. However, by assuming a cylindrical crack surface of radius r 0 its axial extension is estimated from the crack surface energy as z 0 ≈ ∆E s /(2π r 0 G c ) ≈ 2.51 × 10 −2 mm (for a = 0.5 mm and σ t = 115 MPa), a value not too much different from the length estimated using FFM [6].  radius a and the tensile strength σ t . Radial crack positions obtained by the phase field as well as corresponding results obtained by finite fracture mechanics [6] are somewhat smaller compared to the experimental findings [1]. However, the trend regarding relatively smaller radial positions for larger indenter radii a is displayed by both methods. In addition, the influence of the material strength σ t on the initial crack position r 0 , which has been observed in the FFM investigation [6], is also well reproduced by the phase field simulations, see Fig. 4a). Quantitatively, regarding varying indenter radii and tensile strengths, www.gamm-proceedings.com FFM and phase field results for the ring crack position r 0 are in the same order of magnitude, especially for indenter radii above a = 1.5 mm. Figure 4b) shows the critical indenter displacement d c /a to initiate fracture which is relatively larger for smaller indenters. This result is also in accordance with experimental findings [1]. Moreover, in line with FFM results [6], larger tensile strengths σ t require higher loading d c to initiate the ring crack.

Effect of Poisson's ratio on cone crack propagation
As discussed in Sect. 1, once the initial ring crack is established subsequent crack growth forms a cone. The angle of this cone depends on Poisson's ratio as reported by, e.g., [4]. In the following, special focus is given to the reproduction of this relation with the phase field method. A key assumption is that the cone angle depends only on ν. This allows comparison of different materials without taking different E, G c and σ t into account. Figure 5 shows the formation of cone cracks for different Poisson's ratios in the range ν = 0.1 . . . 0.3 (with the same E, G c and σ t ). The larger is Poisson's ratio the smaller is the cone crack angle (dϕ/dν < 0), which is in agreement with studies from literature. In Figure 6 cone crack angles ϕ obtained by the present simulations are compared with experimental [1,4] and other numerical results [4,5]. Cone angles obtained with the phase field method vary almost linearly between ϕ = 38.1 • for ν = 0 and ϕ = 19.3 • for ν = 0.3 and are somewhat larger compared to the other results. This might be explained by the fact that the cone angle slightly decreases with further crack progress. Cone angles from experiments [4] are measured at a larger radius and therefore are apparently smaller.   Fig. 6: Cone crack angle ϕ depending on Poisson's ratio ν. Results of phase field simulations in comparison with experimental findings reported in [1] and [4] and numerical results from FEM [4] and X-FEM [5].