A coupled approach to predict cone-cracks in spherical indentation tests with smooth or rough indenters

A coupled approach to predict cone-cracks in spherical indentation tests with smooth or rough indenters M.R. Marullia, J. Bonaria, J. Reinosob, M. Paggia,∗ aIMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy bGrupo de Elasticidad y Resistencia de Materiales, Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de Los Descubrimientos s/n, 41092 Sevilla, Spain Abstract Indentation tests are largely exploited in experiments to characterize the mechanical and fracture properties of the materials from the resulting crack patterns. This work proposes an efficient theoretical and computational framework, whose implementation is detailed for 2D axisymmetric and 3D geometries, to simulate indentation-induced cracking phenomena caused by non-conforming contacts with indenter profiles of arbitrary shape. The formulation hinges on the coupling of the MPJR (eMbedded Profile for Joint Roughness) interface finite elements which embed the indenter profile to solve the contact problem between non-planar bodies efficiently and the phase-field for brittle fracture to simulate crack evolution and nonlocal damage in the substrate. The novel framework is applied to predict cone-crack formation in the case of indentation tests with smooth spherical indenters, with validation against experimental data. Then, the methodology is employed for the very first time in the literature to assess the effect of surface roughness superimposed on the shape of the smooth spherical indenter. In terms of physical insights, numerical predictions quantify the dependencies of the critical load for crack nucleation and the crack radius on the amplitude of roughness in comparison with the behavior of smooth indenters. Again, the consistency with available experimental trends is noticed. This article has been published on the Journal of the Mechanics and Physics of Solids at https://doi. org/10.1016/j.jmps.2023.105345 Keywords: spherical indentation, fracture mechanics, contact mechanics, MPJR interface finite elements, phase-field fracture, roughness. 1. Introduction An indentation test consists of a sharp or spherical hard tip pressed onto a sample to be mechanically characterized. When the tip is removed after indentation, an impression is left on the target material, and the information on the resulting inelastic deformation can be exploited to assess the mechanical properties of the substrate, see (Hutchings, 2009) for a review of the development of the method before and after 1950. For instance, the method has been largely applied to estimate the hardness Hof an elasto-plastic material, being Hthe material resistance to plastic deformation, defined as the indentation load divided by the area of impression: H=Pm/AC. In the above equation, ACis the projected contact area between the indenter tip and the surface evaluated at the indentation load Pm. In case of brittle and quasi-brittle materials, like glass, silicon, ceramics, etc., cracks develop during indentation. The dimension of these cracks and the indentation loads are related to other material properties of the specimen, and in particular, to the fracture toughness (Lawn, 1998). ∗Corresponding author Email address: [email protected] (M. Paggi) Preprint submitted to Journal of the Mechanics and Physics of Solids June 26, 2023 arXiv:2306.13539v1 [cond-mat.mtrl-sci] 23 Jun 2023 Different types of indentation tests can be performed according to the tip’s dimension and shape. Sharp pyramidal indenters are used in Vickers and Berkovich indentation tests, while the Brinell indentation exploits a spherical tip. For this reason, indentation tests can be applied to different classes of materials, from ceramics and metals, to composites and bio-materials, and at different scales. In this regard, the crack morphology generally depends on the indenter shape: median, radial, and lateral cracks can be recognized in the post-indentation samples (Johanns, 2014; Schneider et al., 2012). Spherical indentation tests on quasi-brittle materials (see the sketch in Fig. 1) present the complexity related to the fact that the contact is non-conformal, i.e., the contact area between the sphere and the substrate grows with the applied load. Fracture is characterized by the occurrence of cone-shaped cracks. Moreover, by further increasing the load level, the expanding contact circle may reach and overcome the position of the surface ring crack, resulting in the generation of secondary ring cracks (Lawn, 1998). The formation of such cracks has been pioneeringly studied by Hertz (1882) at the end of the 19th Century. By assuming frictionless contact between two elastic solids, the Hertzian description of the stress field under a spherical indenter was used to interpret the fracture nucleation and propagation stages. The fracture pattern develops under the form of a surface ring crack close to the contact boundaries, which then propagates as a cone-crack in the substrate by further increasing the applied load, see the sketch in Fig. 1. Most of the early experiments on cone fracture were conducted on glass, notably soda–lime glass, or PMMA, under static loading or by considering a free-fall of a spherical indenter onto the specimen. The transparency of these brittle materials allows easy visualization of the crack growth (Lawn, 1998; Lee et al., 2012; Puttick, 1978; Ritter et al., 1988; Schneider et al., 2012). (a) (b) Figure 1: Hertzian fracture topology: (a) sketch of a spherical indentation test and cone-shaped crack; (b) top view and crosssection of the Hertzian indentation crack in a Silicon Nitride specimen caused by a Tungsten Carbide sphere with Rs= 1.98 mm, adapted from (Lawn, 1998). A cone-crack topology can also be obtained by indenting a substrate with a cylindrical flat punch, see Fig. 2. In this case, the theoretical analysis of the flat cylinder problem is easier than the spherical one since the type of contact is conformal, i.e., the extension of the contact area does not change with the applied load, and, in this case, it is given by the punch base area. An experimental comparison between the flat punch and the spherical indenter geometries has been conducted in (Mouginot and Maugis, 1985). For both indenting geometries, the experimental tests showed that the ring-crack radius r0is always larger than the contact radius at the crack onset ac: the ratio r0/ac varies from 1.1 to 1.4 for the spherical indentation tests, and it is even higher for the flat punch. Moreover, the ratio r0/acwas found to be a decreasing function of the indenter radius Rs(Conrad et al., 1979; Mouginot and Maugis, 1985). After Hertz, Auerbach (1893) extensively studied the effect of the indenter size of spherical indenters on cone-crack initiation and gave an empirical linear relationship between the critical indentation load Pcand 2 (a) (b) Figure 2: Indentation fracture caused by a flat-ended cylinder: (a) sketch of the indentation test and crack topology; (b) cone-shaped crack obtained with a cylindrical punch in soda-lime glass, adapted from (Lawn, 1998). the sphere radius (Pc∝Rs, known as Auerbach’s law). Later, it was noticed that the linear relationship holds only for small indenters when ais comparable with Rs. The linear dependency has also been questioned since experimental observations in (Tillett, 1956) showed that Pcasymptotically grows with R2 sfor larger radii. As an advancement from stress theories that have been applied to assess the critical load for damage nucleation and to predict the direction of crack growth based on the principle stress direction, see e.g. (Mouginot and Maugis, 1985) for flat punches, the application of linear elastic fracture mechanics (LEFM) has been fostered to study more in detail the problem of crack growth. In (Kocer and Collins, 1998; Roesler, 1956; Tillett, 1956), LEFM has been used to simulate the propagation of a crack from an existing flaw placed on the surface and to assess the evolution of the crack inclination based on the updated stress field during crack growth, which is clearly different from the estimate based on the uncracked configuration on the base of stress theories. For more details on the experimental and theoretical analyses of the indentation-induced fracture, see, e.g., the review articles in (Guin and Gueguen, 2019; Kocer, 2003). As far as numerical approaches are concerned, an incremental finite element model has been employed in (Kocer, 2003) to determine the trajectories of the Hertzian cone cracks for various initiation radii by evaluating the change of the stress field caused by crack propagation. The preferred direction of crack growth for its incremental extension was determined using the criterion of maximum strain energy release rate (Sun and Jin, 2012). Cone cracking resulting from spherical indentation has also been simulated using the extended finite element method (XFEM) implemented in Abaqus in (Marimuthu et al., 2017) to improve remeshing operation due to crack growth. A more recent analysis can be found in (Strobl et al., 2017), where the Hertzian cone crack obtained with a cylindrical punch has been analyzed in the context of Finite Fracture Mechanics (FFM). This approach requires the simultaneous fulfillment of a strength and a fracture criterion that involve two well-defined material parameters, i.e., the (tensile) strength σcand the specific (i.e. per unit area) fracture energy Gc(Doitrand et al., 2022; Leguillon, 2002; Leguillon et al., 2018). The FFM approach in (Strobl et al., 2017) reproduced the decrease of the ratio r0/a and the increase of Pc/a3/2with the increase of the indenter radius, consistently with experiments. The FFM coupled criterion has also been applied in (Hahn and Becker, 2021) where the authors well reproduced the experimental results in (Mouginot and Maugis, 1985) using, however, different Gcand σccompared to the values available in the literature. Cone-shaped fracture has also been recently studied using the Phase-Field (PF) approach. This method relies on the seminal work by Francfort and Marigo (1998), see also the subsequent developments in (Amor et al., 2009; Bourdin et al., 2000; Miehe et al., 2010). The PF method is based on Griffith’s idea of the competition between the elastic and fracture energy counterparts: a pre-existing crack propagates if the increase in the surface energy necessary to create a new crack front is balanced by the reduction of the elastic strain energy stored in the body (Griffith, 1921). In the PF variational approach, the fracture phenomenon is treated as nonlocal damage, and it is solved through an energy minimization which, in the 3 Γ-convergence limit, consistently reproduces the Griffith theory of fracture. Therefore, a damage variable at the material point level is introduced as an additional primary unknown of the problem in addition to the displacement field. A recent review of the method can be found in (Wu et al., 2020). The approach has been proved to successfully capture crack nucleation and propagation not only in isotropic linear elastic materials but also in the case of plasticity (Ambati et al., 2015; Miehe et al., 2016), hyperelasticity (Mandal et al., 2020; Marulli et al., 2022; Miehe and Sch¨anzel, 2014; Russ et al., 2020) and heterogeneous composites (Carollo et al., 2018; Dean et al., 2020; Guill´en-Hern´andez et al., 2020). Experimental validation of two different variants of the PF method, with a rigorous convergence check of the staggered formulation, has been discussed in (Cavuoto et al., 2022) in relation to PMMA specimen geometries with circular holes and V-notches, leading to complex patterns featuring crack nucleation, propagation, and even secondary propagation events. Concerning indentation problems, the PF approach has been applied to cone-shaped fracture in (Kumar et al., 2022; Strobl and Seelig, 2019, 2020; Wu et al., 2022), where the authors discussed the most appropriate phase-field formulation to simulate cracking of glass substrates due to flat-ended cylindrical punches, proposing a specific form of the tensile-compressive strain energy specific for the indentation tests (Strobl and Seelig, 2019, 2020). An alternative formulation has been given in (Wu et al., 2022) with the so-called phase-field cohesive zone model (PF-CZM), where the geometric crack function and the energetic degradation function have been adapted to reproduce a brittle fracture law with a linear softening curve. A revised PF approach can also be found in the work of Kumar et al. (2022), where the authors introduced a strength of the glass boundary layer on the surface (60 MPa) different from that in the bulk (150 MPa), to account for microscopical defects located on the material surface, frequent in glass. According to the above state of the art, it has to be noticed that less attention has been so far devoted to the case of spherical indentation, which has been discussed numerically using the PF approach only in (Phase Field Modeling of Hertzian Cone Cracks Under Spherical Indentation, 2020) for the case of a sphere with radius 1 mm. The problem of spherical indentation requires, in fact, the solution of the contact problem in addition to the simulation of fracture in the bulk, which implies the solution of a strong nonlinear mechanical problem with two forms of nonlinearities. In (Phase Field Modeling of Hertzian Cone Cracks Under Spherical Indentation, 2020), the continuous update of the contact area was addressed by solving the contact problem with a penalty method, coupling the mechanical and the phase-field evolution equations using a staggered solution scheme. To make further progress on the prediction of cone cracks resulting from non-conforming contact, we propose here an innovative computational approach where we simulate the contact problem with the MPJR (eMbedded Profile for Joint Roughness) interface finite element formulation proposed in (Reinoso and Paggi, 2014), and extended in (Bonari and Paggi, 2020; Bonari et al., 2022) for frictional and adhesive contact problems. The method allows considering the contact surfaces as nominally flat, but preserving the actual geometry of the indenter by embedding its geometry into an internal variable of the interface finite element to correct the normal gap function. As proved in the previous publications, the method can be successfully applied to 2D and 3D contact problems with a smooth shape of the indenter (spherical, wavy) and also to the more challenging contact problem involving a random fine-scale roughness. Therefore, this predictive methodology precludes the need to explicitly discretize the actual boundary geometry, which is, however, fully taken into account as an exact correction to the gap function. In this study, for the simulation of crack growth, the state-of-the-art version of the PF approach to fracture in (Wu et al., 2022) is also implemented for the continuum substrate. To make both formulations compatible, the nodal degrees of freedom of the MPJR interface finite elements are augmented to be consistent with the degrees of freedom of the finite elements used to discretize the continuum with the PF, which involve the nonlocal damage variable in addition to the displacement field components. This methodology also opens new perspectives in terms of modeling constitutive coupling effects between damage in the bulk surrounding the surface and the tribological properties, such as friction, wear, and adhesion. Moreover, the 3D formulation has also been detailed for 2D axisymmetric problems since they are particularly relevant to efficiently simulate spherical indentation fracture reducing the computation cost as compared to full-scale 3D simulations. The resulting overall computational formulation, equipped with both nonlinearities due to contact and 4 fracture, is herein tested for the first time in relation to spherical indentation with smooth or rough spheres. It is, in fact, known that another critical aspect often neglected in indentation-induced fracture simulations is the influence of surface morphology on the test results. In this concern, experimental results provide non-consistent trends depending on the indenter material. In the case of steel indenters, the critical load and the ring crack radius in the abraded cases were greater than on the as-received surfaces, as reported in (Conrad et al., 1979; Jyh-Woei et al., 1993; Lu et al., 1995; Mouginot and Maugis, 1985). Similarly, etching treatments in hydrofluoric acid of the glass surface increased the critical load Pcwith respect to the simply polished surface in (Hamilton and Rawson, 1970). Similarly, steel ball indentation on silicon wafer substrate in (Jyh-Woei et al., 1993) showed that the critical load was higher for the polished silicon surface than for the abraded one, while the crack ring radius increased. On the other hand, in the case of a glass indenter on a glass substrate (Johnson et al., 1973), surface abrasion was reported not to affect the critical load. Hence, we expect that mechanical properties and surface roughness resulting from different surface treatments might influence crack formations, and a computational framework accounting for surface roughness is desirable to gain insight into these issues. Therefore, the present article aims at developing a new simulation framework with strong mechanical foundations that allows testing the potential of the PF approach coupled with the MPJR interface formulation for non-conformal indentation cracking problems. The article is structured as follows: Sec. 2 provides the overall coupled problem formulation, and it includes the treatment of the contact problem in Sec. 2.1 and the phase-field formulation for fracture in Sec.2.2. The axisymmetric MPJR method, presented for the first time in this work, is applied in a benchmark test to prove its efficiency against a standard node-to-segment approach in Section 3. Section 4 shows the application of the computational method to smooth and rough spherical indenters, with novel insights into the physical problem. 2. Governing equations of the coupled contact and fracture problem The present section describes the general variational framework of the system which includes: (i) the contact interaction between two solid domains (Ωiwith i= 1,2) at the contact interface Γc; and (ii) the treatment of the brittle fracture problem. A crack Γfis assumed to nucleate and propagate in the substrate, as shown in Fig. 3. Figure 3: Non-conformal contact between two solids Ω1and Ω2that will lead to a contact domain Γc, function of the applied load or displacement, and an evolving crack Γf∈Ω1. The free energy functional of the system reads: Π(u,Γf,Γc)=ΠΩ1+ ΠΩ2+ ΠΓf+ ΠΓc(1) 5 where ΠΩ1and ΠΩ2denote the total potential energy of the solids, ΠΓfis the energy dissipated due to fracture, and ΠΓcis the contribution due to the contact interactions. Both contributions are discussed in the following subsections. 2.1. Contact contribution to the weak form This section focuses on the treatment of the contact problem between the indenter and the substrate. The MPJR interface finite element formulation introduced in (Paggi and Reinoso, 2018) is herein adopted to discretize the interface. The formulation is herein specialized to frictionless and adhesiveless indentation problems, although the method also applies to such scenarios as shown in (Bonari et al., 2022). The proposed numerical method consists of a zero-thickness interface finite element separating the indenter and the substrate. Both materials can be modeled with their actual elastic properties, as in (Reinoso and Paggi, 2014). It can also be applied to the case of a rigid indenter in contact with an elastic substrate as a special case. Let the two solids occupy the 3D domains Ωi∈R3with (i= 1,2) in the undeformed configuration defined by the reference system Oe1e2e3, as shown in Fig. 4. The position of a point in the body is given by the vector of its Cartesian coordinates x. The bodies are separated by an interface Γc∈R2 defined by the two opposite sides Γ1and Γ2. Figure 4: Solid domains Ω1and Ω2interacting through the contact interface Γc. Kinematic and traction boundaries conditions can be prescribed on disjointed parts of the solids’ boundaries such that each boundary can be split into three parts: (i) a portion where displacements are imposed, ∂Ωiu; (ii) a portion where tractions ˆ tare specified, ∂Ωit; (iii) the contact interface Γiwhere contact tractions are exchanged. Let the body Ω1be the substrate with a smooth contact interface, while Ω2is the indenter whose contact interface Γ2has an arbitrary shape that can be described by an analytical function (e.g., a parabolic or harmonic profile) or by a set of discrete data related to a more complex rough topology. The core of the approach consists in simplifying the original boundary of the indenter Γ2into a smooth surface Γ∗ 2, while the actual profile of the boundary is embedded point-wise in its exact form into the interface finite element formulation. The geometrical difference between Γ2and Γ∗ 2is mathematically described by a function z(x) where x= (x, y)Tis the coordinate position vector on the surface, see Fig. 5. This method reduces the contact problem to a nominally flat-to-flat interface Γ∗ c. It is solved by introducing the displacement field of the solids u(i)= (u(i), v(i), w(i))Tsuch that the configuration of the system at the contact interface is described by the gap field g= (gn, gt1, gt2)Tdefined as the projection of the relative displacement u(1) −u(2) onto the normal and tangential directions of the interface: gn=n·(u(1) −u(2)), gt1=t1·(u(1) −u(2)), gt2=t2·(u(1) −u(2)) (2) The normal gap g∗ nis then corrected to restore the exact shape of the indenter as follows: g∗ n=gn+z(x) (3) 6 Figure 5: From original boundary to embedded geometry In the hypothesis of a frictionless and adhesiveless contact problem, the contact traction p= (pn, τ1, τ2)T reduces to the normal component pn. The associated contact conditions expressed using the standard HertzSignorini-Moreau inequalities read: g∗ n≥0, pn≤0, g∗ npn= 0 on Γi(4) and are treated with a penalty approach (Wriggers, 2006)[p. 118]: pn=(0 if g∗ n≥0 kpg∗ nif g∗ n<0(5) where kpis the penalty stiffness. Equation (5) gives a nonzero contact pressure for all the points within the nominally flat-to-flat interface where the corrected normal gap is negative valued. Since the exact geometrical corrective term z(x) is provided at each integration point of the interface finite elements, the surface geometry is embedded directly inside the derivation of the system’s stiffness matrix without the need for an explicit discretization of its geometry, which leads to a significant advantage in terms of modeling and discretization of complex surface topologies. In the case of spherical indentation problems, it is convenient to specialize the above 3D formulation by introducing the 2D axisymmetric model to devise a computationally effective method. In such a case, the displacement field vector of the i-th solid is reduced to two terms, u(i)= (u(i) r, v(i)), where urand v denote, respectively, the displacement in the radial direction and on the vertical direction. The gap vector also reduces to two terms, g= (gn, gt), computed according to the first two equations (2). The corrective term of the normal gap vector is still evaluated according to Eq.(3), where the position vector xcoincides with the radial coordinate r, such that z(x) = z(r). The contact constraint conditions in Eq.(5) still apply, provided that the contact traction vector p= (pn, τ)Twith only two components is considered. The numerical treatment of the contact problem requires the discretization of the interface, which is conducted with 8-nodes MPJR interface finite elements in 3D and with 4-nodes interface finite elements in 2D, see Fig. 6. A conformal mesh discretization is adopted for the continuum at the interface between the two solids. The solution method has, therefore, the same features as a segment-to-segment contact algorithm with fixed pairings (Paggi and Wriggers, 2016). Given the nodal displacement vector ¯ u= (u1, v1, w1, ..., u8, v8, w8)Tfor a 3D problem, or ¯ u= (ur1, v1, ..., ur4, v4)T for the 2D axisymmetric one, the normal gap inside the interface finite element can be derived by computing the relative displacement ∆ubetween the couples of nodes from the opposite sides of the interface using a matrix operator L. The result is then interpolated at any point inside the interface element through standard linear shape functions collected in the matrix operator N. Finally, the normal and tangential components of the gap are determined using a rotation matrix Rdefined by the unit vectors n,t1and t2in 3D, or just 7 Figure 6: 2D and 3D interface finite elements. Figure 7: Element-wise profile discretization. nand t1in 2D, both related to the local reference system of the interface finite element. In formulae, the discretized gap field can be written as: g=RNLu (6) The deviation from planarity of the shape of the indenter profile can be taken into account by computing the corrected gap vector g∗at each interface integration point, according to Eq. (3). If an analytical function is used to define the profile shape, then the correction z(x) is computed by introducing the coordinates of the interface finite element nodes. Otherwise, if the surface/profile data are provided as a discrete set of elevations, as from data acquired from a profilometer or AFM, then those data are provided in input to the software. Such input data are stored in a history variable inside the user element routine only once, at the initialization of the problem. A mapping routine connects the external data to the proper node. The finite element discretization at the interface will be related to the spatial spacing of the external data. If the phase-field problem requires a discretization of the continuum finer than the sampling spacing of the surface data field at the interface, which can be due to the internal length-scale constraint of the method (see Sec. 2.2), then a linear interpolation of the input heights field is performed to locally refine the data assigned to the nodes of the conformal interface discretization. Such an issue is not at stake in case of pure contact problems without fracture. The contribution of a single interface finite element to the variational formulation of the system is given by: δΠ = ZΓ∗ c δg∗ npndA,(7) where dA denotes the area of the surface patch in 3D, while it is equal to 2πr dr for the 2D axisymmetric case. The equation leads to the following expressions for the element residual vector Rcand the element contact stiffness matrix Kcassociated with the mechanical field to be used in a Newton-Raphson incremental-iterative solution scheme: Rc=ZΓ∗ c LTNTRTtdA,Kc=ZΓ∗ c LTNTRTCRNL dA (8) 8 where Cis the linearized interface constitutive matrix having only one nonzero component C11 =∂pn/∂gn= kpfor the contact points. The combined phase-field interface finite element framework has been herein implemented using a monolithic fully implicit solution strategy in the finite element software FEAP8.6 (Zienkiewicz et al., 2013), following the methodology discussed in (Reinoso et al., 2017). To make compatible the interface finite element with the phase-field finite element for the bulk which presents an additional degree of freedom associated to the damage variable in each node, the element residual vector and the element stiffness matrix contributions are expanded by mapping the corresponding terms in reference to the augmented vector of degrees of freedom ¯ ue= (u1, v1, w1, ϕ1, ..., u8, v8, w8, ϕ8)Tfor a 3D problem, or ¯ ue= (ur1, v1, ϕ1, ..., ur4, v4, ϕ4)Tfor the 2D axisymmetric case. Operatively, this is done by using a matrix operator Pwhose expression is collected in the Appendix for the 3D and the 2D cases, P2D and P3D: Re c=PRc,(9a) Ke c=PKcPT.(9b) 2.2. phase-field contribution to the weak form The present subsection deals with the governing equations of the continuum solids in contact, which are treated according to the phase-field approach for fracture. The spherical indentation problem does not lead to micro-cracks, phase transformation in the materials, or plasticity, in contrast to sharp indentation. Hence, the tested materials will be considered linear elastic with nonlocal damage that, in the phase-field approach, tends to Griffith fracture in the limit of a vanishing internal length scale governing the nonlocality of the method. The variational formulation is developed considering a standard solid domain Ωi⊂R3with i= 1,2 having boundaries denoted by ∂Ωi, see Fig. 8. (a) (b) Figure 8: A cracking solid Ωiwith a sharp crack Γfin (a) and its phase-field nonlocal regularization in (b). It is also assumed that the arbitrary solid is subjected to specific body forces ¯ b, imposed displacements ¯ uon Ωi,t, and boundary tractions ¯ ton Ωi,t, possibly caused by contact. The external potential energy of the system reads: P=ZΩi\Γf ¯ b·udV + Z∂Ωi,t ¯ t·udA (10) and the free energy functional of the systems is: Π(u,Γf) = ZΩi\Γf ψ(ϵ) dV + ZΓf GcdA (11) Where ψ(ϵ) is the elastic energy density of the body written in terms of the strain field ϵ. The integral RΓfGcdA identifies the energy dissipation due to fracture events at the crack set Γf, while Gcis the critical 9 (a) ¯u= 5.2×10−3mm (b) ¯u= 5.4×10−3mm (c) ¯u= 6.5×10−3mm (d) ¯u= 8.4×10−3mm Figure 15: Snapshots of the contour plots of the vertical displacement component for the indentation test with a sphere of radius Rs= 1 mm. In (d), the displacement discontinuity due to fracture is evident from the jump from blue to red colors. Figure 16: 3D post-process view of the 2D axisymmetric crack pattern where fractured finite elements with ϕ∼ =1 have been removed, for ease of visualization, at a far field vertical displacement of ¯u= 0.01 mm. 16 Figure 17: Crack pattern due to the indentation with a sphere having R= 1 mm shown in the entire domain and in a magnified view. (a) ¯u= 5.2×10−3mm (b) ¯u= 5.4×10−3mm (c) ¯u= 6.5×10−3mm (d) ¯u= 1.×10−2mm Figure 18: Snapshots of the fracture evolution for the indentation test with a sphere radius 1 mm. 17 A quantitative analysis of the results shows that the interface contact pressure increases according to the Hertzian pressure distribution as shown in Fig. 19(a), where the red curve corresponds to the crack onset. For the same pseudo-time steps, the phase-field variable evolution along the interface is shown in Fig. 19(b). The contact pressure and the phase-field plots are used to evaluate, respectively, the contact radius at crack nucleation, ac, and the ring crack radius, r0. The contact radius at crack nucleation is defined by the coordinate where the contact pressure becomes vanishing after having been negative valued, for the red curve in Fig. 19(a) that corresponds to the first point of the interface where the phase-field variable reaches unity. The coordinate of that point, which can be assessed from Fig. 19(b), gives the ring crack radius, r0. Hence, the propagation of the crack outside the contact area can be clearly noted in Fig. 19 by comparing the two red curves which gives r0/ac= 1.28. (a) (b) Figure 19: Contact pressure in (a) and phase-field variable in (b) along the interface at different pseudo-time steps for the indentation test with Rs= 1 mm. The red curves correspond to the variables at the point of crack nucleation. The effect of the sphere radius on the critical load causing crack nucleation has been investigated by considering five cases with a radius Rsvarying from 1 to 7.5 mm. The introduction of the interface finite elements allows changing the spherical indenter radius without changing the model geometry since the modified spherical geometry is embedded into the finite elements and analytically defined as a function of its radius. The results of the simulations are collected in Fig. 20 where current simulations are compared with the values obtained in the experimental campaign in (Conrad et al., 1979) for smooth spheres. This graph shows a very satisfactory agreement between the experimental and numerical data. 18 Figure 20: Critical load variation with respect to the sphere indenter radius and comparison with experimental data from Conrad et al. (1979) for the case of as-received (smooth) specimens. The ring crack radius r0is always greater than the critical contact radius ac, consistently with results in (Conrad et al., 1979; Jyh-Woei et al., 1993; Mouginot and Maugis, 1985), see the plots in Fig. 21 where the ratio r0/acis related to the indenter radius, Rs. Since no experimental data are available in (Conrad et al., 1979), the simulations results have been compared in Fig. 21a against the indentation tests presented in (Jyh-Woei et al., 1993) where the authors used the same type of glass and test geometry for the specimens. In this case, the scatter in the experimental data does not allow to confirm the obtained numerical trend. However, although a one-to-one comparison is not possible, the predicted trend is consistent with that found in the experiments in (Mouginot and Maugis, 1985) for flat-ended cylindrical punch (Fig. 21b, upper panel) and for spherical indenters on glass specimens abraded with 1000 grit silicon carbide paper (Fig. 21b, lower panel). Although those experimental trends have been obtained with glass specimens of different sizes (50 mm× 50 mm ×25.4 mm) and material properties (borosilicate glass with E= 80 GPa, ν= 0.22 ) than those used in the simulations in Fig. 21a, the trends are fully consistent and confirm that the ratio r0/acdecreases by increasing the indenter radius. The variation of the ratio r0/achas been related to the inverse of the contact radius 1/acin Fig. 22a and compared with the experimental data in (Conrad et al., 1979) which show again a very high scatter. For this reason, the trend resulting from the numerical model has been compared also with that resulting from the experimental tests in (Mouginot and Maugis, 1985) for different types of indenter, see Fig. 22b. It has to be noted that in the case of a flat-end cylindrical punch, the radius of the indenter Rsalways coincides with the contact radius, which does not happen in the case of non-conformal contacts, as explained in the introduction. The data on the spherical indenters in Fig. 21b, as well as the investigations in (Conrad et al., 1979; Jyh-Woei et al., 1993), show that the surface treatment (abrasion with 1000 grit SiC paper, or with 7 µm diamond paste) influences the indentation tests and increases the r0/acratio. This latter aspect will be treated in the following paragraphs dealing with the case of rough indenters. 19 (a) Simulations results compared with the experimental data in (Jyh-Woei et al., 1993) on as-received glass substrates. (b) Indentation of flat-ended cylinders and spheres on abraded glass substrates with 1000 grit silicon carbide paper in (Mouginot and Maugis, 1985). Figure 21: Ratio between the crack radius and the critical contact radius vs the indenter radius. The experimental data in (b) concern indentation tests on glass specimens with different mechanical and geometrical properties. (a) Results of the numerical simulations compared with the experimental data in Conrad et al. (1979). (b) Indentation of flat-ended cylinders and spheres on glass abraded with 1000 grit SiC paper or 7 µm diamond paper in Mouginot and Maugis (1985) Figure 22: Ratio between the crack radius and the critical contact radius vs the inverse of the contact radius. To further exploit the advantage of the present approach, which enables to embed any indenter profile along a nominally flat interface and efficiently solve the contact problem, we now consider the effect of surface roughness on cone indentation fracture, using the same axial symmetric finite element model for the smooth sphere, but modifying the embedded indenter profile. Generally, a real rough surface does not present axial symmetry; however, a 3D indentation model requires a high computational effort compared to the axisymmetric configuration. For example, the 2D axisymmetric model of the smooth spherical indentation test presented at the beginning of this section requires a computational time of around 30 min to solve 100 pseudo-time steps up to an imposed displacement 20 of 0.01 mm. In the case of a 3D simulation of the same problem, because of the drastic increase of the degrees of freedom, parallel computing facilities are required: with 20 cores, the CPU time for a 3D FE simulation becomes comparable with the CPU time for a single CPU calculation of the 2D axisymmetric problem. Moreover, for the 3D simulations of a rough sphere, the CPU time increases up to 5 h, even with parallel computing, because of the more complex damage pattern, as shown later. For this reason, the 2D axisymmetric configuration has been chosen for the simulations in this work. To the best of the authors’ knowledge, even the 2D axisymmetric solutions represent the first attempt to investigate the effect of roughness for this kind of complex nonlinear coupled problem involving contact and fracture. For this purpose, one rough surface has been generated using the Random Midpoint Displacement (RMD) algorithm already widely exploited in the contact mechanics literature to simulate realistic surface roughness (Paggi and Ciavarella, 2010), with a spatial resolution of 2.5µm, and fractal dimension D= 2.1. From the numerically generated rough surface, one profile has been extracted and superimposed to the spherical shape of the smooth indenter (see the profile in Fig. 23(a) and the resulting shape of the indenter in Fig. 23(b). In the experimental data in (Conrad et al., 1979; Jyh-Woei et al., 1993; Mouginot and Maugis, 1985), the authors did not directly measure the specimens’ roughness due to the abrasion process; however, they reported the presence of defects on the surface of a few microns after the treatment with grit papers or diamond paste. The rough profile chosen for the simulation has the statistical parameter Rzequal to 1.3µm, which measures the average peak-to-valley distance of the profile. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 -6 -4 -2 0 2 4 610-4 (a) Profile obtained with RMD, with Rz= 1.3µm and sampling 2.5µm. (b) Final profile of the indenter embedded into the interface. Figure 23: Profile of the rough spherical indenter (b) embedded into the MPJR interface finite elements, obtained by superimposing the rough profile (a) to a sphere of radius Rs= 5 mm. The fracture evolution for the case with rough indenter radius Rs= 5 mm is detailed in Fig. 24. The crack propagates with a conical shape, as in the case of a smooth spherical indenter. However, the presence 21 of roughness causes local damage at different points on the glass surface before the main crack propagates. (a) ¯u= 9.2×10−3mm (b) ¯u= 9.3×10−3mm (c) ¯u= 9.6×10−3mm (d) ¯u= 0.01 mm Figure 24: Evolution of damage and the main crack for the indentation of glass with a rough spherical indenter with Rz= 1.3µm. The local damage can also be seen by looking at the phase-field variable plotted along the interface in Fig.s 25, 26 and 27 together with the contact pressure at the interface. In Fig. 25, the two variables have been plotted at different time steps until the phase-field variable is ≈1 at the interface coordinate x= 0.175 mm, see the blue curve corresponding to the imposed displacement ¯u= 9.2×10−3mm. The main crack does not propagate from that point but from the last peak of the phase-field red curve in Fig. 26 at x= 0.285 mm, representing the ring radius r0. This point is again outside the contact area, as in the case of the smooth indenter, since the last point in contact has a radial coordinate x= 0.2625 mm. No further damage along the interface can be seen after the main crack propagates, as shown in Fig. 27. The case of the rough indenter has been compared with the smooth case, considering the same spherical radius Rs= 5 mm, in terms of contact pressure and phase-field along the interface at the onset of main crack propagation, which happens at an imposed displacement equal to ¯u= 0.92 mm for the smooth sphere and ¯u= 0.93 mm for the rough indenter. The comparison shown in Fig. 28 highlights that roughness induces stress concentrations and damage at the points in contact, which do not occur for smooth profiles. The rough profile presented in Fig. 23 has been rescaled to simulate three different maximum peak-tovalley distances Rz= 1.3µm, Rz= 2.6µm, and Rz= 5.2µm. The comparison with the corresponding smooth spherical indenter shows that the ring crack radius increases by amplifying roughness, see Fig. 29 where the predictions for the smooth case and the three rough profiles are shown. For each case, the contour plots are taken at the time step corresponding to the crack propagation. The presence of roughness also affects the critical load for the crack initiation, as highlighted in Fig. 30. The result is in line with the experimental data in (Conrad et al., 1979), where an increase of the critical load 22 (a) (b) Figure 25: Rough spherical indenter with Rz= 1.3µm: (a) contact pressure; (b) phase-field variable for different pseudo-time steps up to the occurrence of the first point along the interface with ϕ≈1 (blue curve, corresponding to a displacement ¯u= 9.2×10−3mm). in the case of abrasion of the glass substrate was reported. In particular, for the indenter radius Rs= 5 mm analyzed here, (Conrad et al., 1979) reported an increase of the critical load from 210N to 300N in case of abrasion with 600 SiC paper, while the simulations show an increase from 173N to 346N for the rough profile with Rz= 5.2µm. 23 (a) (b) Figure 26: Rough spherical indenter with Rz= 1.3µm: (a) contact pressure; (b) phase-field variable for different pseudo-time steps from the occurrence of the first point along the interface with ϕ≈1 (blue curve), up to the situation corresponding to crack propagation for a higher displacement equal to ¯u= 9.3×10−3mm (red curve). 24 (a) (b) Figure 27: Rough spherical indenter with Rz= 1.3µm: (a) contact pressure; (b) phase-field variable for different pseudo-time steps from the occurrence of crack propagation (red curve) onward, up to an imposed displacement equal to ¯u= 0.01 mm. 25