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 ecient theoretical and compu- tational framework, whose implementation is detailed for 2D axisymmetric and 3D geometries, to simulate indentation-induced cracking phenomena caused by non-conforming contacts with indenter pro les of arbi- trary shape. The formulation hinges on the coupling of the MPJR (eMbedded Pro le for Joint Roughness) interface nite elements which embed the indenter pro le to solve the contact problem between non-planar bodies eciently and the phase- eld 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 rst time in the literature to assess the e ect 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. The nal version of this article has been published on the Journal of the Mechanics and Physics of Solids athttps://doi.org/10.1016/j.jmps.2023.105345 Keywords: spherical indentation, fracture mechanics, contact mechanics, MPJR interface nite elements, phase- eld 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, beingHthe material resistance to plastic deformation, de ned 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: marco.paggi@imtlucca.it (M. Paggi) Preprint submitted to Journal of the Mechanics and Physics of Solids June 23, 2023 Di erent 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 di erent classes of materials, from ceramics and metals, to composites and bio-materials, and at di erent 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 eld 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 cross- section 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 at punch, see Fig. 2. In this case, the theoretical analysis of the at 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 at 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 ratior0=ac varies from 1 :1 to 1:4 for the spherical indentation tests, and it is even higher for the at punch. Moreover, the ratior0=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 e ect 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 at-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 at 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 aw placed on the surface and to assess the evolution of the crack inclination based on the updated stress eld during crack growth, which is clearly di erent from the estimate based on the uncracked con guration 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 nite 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 eld 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 nite element method (XFEM) implemented in Abaqus in (Marimuthu et al., 2017) to improve remeshing opera- tion 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 ful llment of a strength and a fracture criterion that in- volve two well-de ned material parameters, i.e., the (tensile) strength cand the speci c (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=aand 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, di erent 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 Grith'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 (Grith, 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 Grith 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 eld. 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 di erent 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 appro- priate phase- eld formulation to simulate cracking of glass substrates due to at-ended cylindrical punches, proposing a speci c form of the tensile-compressive strain energy speci c 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- eld cohesive zone model (PF-CZM), where the geometric crack function and the energetic degrada- tion 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) di erent 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- eld 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 Pro le for Joint Roughness) interface nite 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 at, but preserving the actual geometry of the indenter by embedding its geometry into an internal variable of the interface nite 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 ne-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 nite elements are augmented to be consis- tent with the degrees of freedom of the nite elements used to discretize the continuum with the PF, which involve the nonlocal damage variable in addition to the displacement eld components. This methodology also opens new perspectives in terms of modeling constitutive coupling e ects 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 eciently 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 rst 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 in uence 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 hydro uoric 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 a ect the critical load. Hence, we expect that mechanical properties and surface roughness resulting from di erent surface treat- ments might in uence 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 formula- tion 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- eld formulation for fracture in Sec.2.2. The axisymmetric MPJR method, presented for the rst time in this work, is applied in a benchmark test to prove its eciency 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 ( iwithi= 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  f2 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 nite 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 nite element separating the in- denter 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 i2R3with (i= 1;2) in the undeformed con guration de ned 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  c2R2 de ned 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' bound- aries such that each boundary can be split into three parts: ( i) a portion where displacements are imposed, @ iu; (ii) a portion where tractions ^tare speci ed, @ it; (iii) the contact interface  iwhere contact trac- tions 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 pro le) 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 pro le of the boundary is embedded point-wise in its exact form into the interface nite element formulation. The geometrical di erence between  2and  2is mathematically described by a functionz(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 at-to- at interface  c. It is solved by intro- ducing the displacement eld of the solids u(i)= (u(i);v(i);w(i))Tsuch that the con guration of the system at the contact interface is described by the gap eld g= (gn;gt1;gt2)Tde ned 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 Hertz- Signorini-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 ifg n0 kpg nifg n<0(5) wherekpis the penalty sti ness. Equation (5) gives a nonzero contact pressure for all the points within the nominally at-to- at 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 nite elements, the surface geometry is embedded directly inside the derivation of the system's sti ness matrix without the need for an explicit discretization of its geometry, which leads to a signi cant 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 e ective method. In such a case, the displacement eld vector of the i-th solid is reduced to two terms, u(i)= (u(i) r;v(i)), whereurandv 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 rst 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 nite elements in 3D and with 4-nodes interface nite 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 xed pairings (Paggi and Wriggers, 2016). Given the nodal displacement vector u= (u1;v1;w1;:::;u 8;v8;w8)Tfor a 3D problem, or u= (ur1;v1;:::;ur4;v4)T for the 2D axisymmetric one, the normal gap inside the interface nite 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 Rde ned by the unit vectors n,t1andt2in 3D, or just 7 Figure 6: 2D and 3D interface nite elements. Figure 7: Element-wise pro le discretization. nandt1in 2D, both related to the local reference system of the interface nite element. In formulae, the discretized gap eld can be written as: g=RNLu (6) The deviation from planarity of the shape of the indenter pro le 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 de ne the pro le shape, then the correction z(x) is computed by introducing the coordinates of the interface nite element nodes. Otherwise, if the surface/pro le data are provided as a discrete set of elevations, as from data acquired from a pro lometer 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 nite element discretization at the interface will be related to the spatial spacing of the external data. If the phase- eld problem requires a discretization of the continuum ner than the sampling spacing of the surface data eld 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 eld is performed to locally re ne 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 nite 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 rdr for the 2D axisymmetric case. The equation leads to the following expressions for the element residual vector Rcand the element contact sti ness matrixKcassociated with the mechanical eld to be used in a Newton-Raphson incremental-iterative solution scheme: Rc=Z cLTNTRTtdA;Kc=Z cLTNTRTCRNL 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- eld interface nite element framework has been herein implemented using a mono- lithic fully implicit solution strategy in the nite element software FEAP8.6 (Zienkiewicz et al., 2013), fol- lowing the methodology discussed in (Reinoso et al., 2017). To make compatible the interface nite element with the phase- eld nite 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 sti ness matrix contributions are expanded by mapping the corresponding terms in reference to the augmented vector of degrees of free- dom ue= (u1;v1;w1; 1;:::;u 8;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, P2DandP3D: Re c=PRc; (9a) Ke c=PKcPT: (9b) 2.2. phase- eld 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- eld 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- eld approach, tends to Grith 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R3withi= 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- eld nonlocal regularization in (b). It is also assumed that the arbitrary solid is subjected to speci c 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 infbudV +Z @ i;ttudA (10) and the free energy functional of the systems is: (u;f) =Z inf () dV +Z fGcdA (11) Where () is the elastic energy density of the body written in terms of the strain eld . The integralR fGcdA identi es the energy dissipation due to fracture events at the crack set  f, whileGcis the critical 9 energy release rate or fracture toughness of the bulk material. This integral cannot be directly evaluated because the crack set  fis apriori unknown. According to the variational approach to fracture proposed in Francfort and Marigo (1998), the problem is solved by substituting the sharp crack with a transition region from undamaged to broken material. The sharp crack is therefore approximated as a band of nite width l0characterized by a crack phase- eld parameter 2[0;1] such that = 0 denotes the intact material and = 1 represents the cracked one. The crack approximation converges to the sharp crack when such a band thickness approaches zero. The energy contribution due to fracture is obtained through its smeared nonlocal approximation: Z fGcdAZ iGc ( ;r ) dV (12) Accordingly, the free energy functional expressed in Eq. (11) becomes (u; ) =Z ig( ) () dV +Z iGc ( ;r ) dV (13) whereg( ) is called energetic degradation function and acts to reduce the elastic sti ness of the material. There are di erent choices for the function g( ). In this work, we used the model introduced by Bourdin et al. in Bourdin et al. (2000): g( ) = (1 )2+kres (14) Where a small positive parameter kresis used to avoid numerical instabilities at the fully cracked state. In Eq. (13), ( ;r ) is the crack surface density function that assumes this generic form for the AT2 approach (Miehe et al., 2010): ( ;r ) =1 2 2 l0+l0jr j2 (15) wherel0is the length scale that de nes the width of the di usive crack band as shown in Fig. 8b. The reader is referred to Amor et al. (2009), Tann e et al. (2018), Strobl and Seelig (2020), and Cavuoto et al. (2022) for a speci c discussion on how to identify the value of l0from experimental data. According to the above state-of-the-art literature, the following variational formulation with respect to the primary eld uand can be derived: u=Z ig( )(u) :(u) dV (16a)  =Z idg( ) d () dVZ Gc l0 +l0r r  dV (16b) whereuand stand for the virtual variation of the displacements and the phase- eld, (u) is the Cauchy stress tensor for the undamaged con guration de ned as =C0:withC0the standard fourth- order sti ness tensor of an isotropic linear elastic material. The strain and stress tensors in Voigt notation for the 3D and the 2D axisymmetric models read, respectively: = [x;y;z; xy; yz; xz]T; = [r;z;; rz]T(17) = [x;y;z; xy; yz; xz]T; = [r;z;; rz]T(18) In Eq. (16b), the elastic strain energy has to be split into its positive and negative counterparts to avoid crack propagation in compression: () = +() + (). Considering only the positive contribution, the same equation becomes:  =Z idg( ) d +() dVZ Gc l0 +l0r r  dV (19) 10 The energy driving the crack growth +has been derived as in Strobl and Seelig (2020): +:=1 + 2E h1i2 ++h2i2 ++h3i2 +  2Ehtr()i2 + (20) wherehiis the Macaulay bracket operator hi+:= (jj+)=2, andiare the eigenvalues of the stress tensor for the intact material. The irreversibility of the phase- eld evolution has been enforced by following the approach in Miehe et al. (2010) by using a history- eld variable of the maximum positive energy contribution: H(u;t) = max t2[0;T] +() (21) wheretis the current pseudo-time step in a quasi-static simulation. With this assumption, Eq. (19) becomes:  =Z idg( ) d H dVZ Gc l0 +l0r r  dV (22) Isoparametric nite elements with standard bilinear shape functions NI() have been used for the spatial discretization of the domain (see Msekh et al. (2015) for more details). The approximated displacement eld, the phase- eld, and their variations read: u=ndX I=1NIuI;  u=ndX I=1NIuI; =ndX I=1NI I;  =ndX I=1NI I (23) wherendstands for the number of nodes for each nite element, uIand Idenote the nodal values of the displacement and phase- eld, respectively, which are collected in the corresponding vectors uand . The strain eld is interpolated through the displacement-strain operator Bu, while the gradient of the phase- eld via B : =Buu;rx =B (24) where Buneeds to be specialised for the 3D and the 2D axisymmetric cases: Bu=2 666666666666666664@NI @x0 0 0@NI @y0 0 0@NI @z @NI @y@NI @x0 0@NI @y@NI @z @NI @z0@NI @x3 777777777777777775; Bu=2 666666664@NI @r0 0@NI @z NI r0 @NI @z@NI @r3 777777775(25) while B for the 3D and the 2D axisymmetric cases read: B =2 666664@NI @x @NI @y @NI @z3 777775; B =2 64@NI @r @NI @z3 75 (26) 11 With the previous interpolation schemes, the residual vectors RuandR associated with the displace- ment and the phase- eld, respectively read: Ru=Z i (1 )2+kres BT udV (27a) R =Z i2(1 )NTHdV +Z iGcl0 BT rx +1 l2 0NT  dV (27b) where the volume integration dVis equal to 2 rdrdz for the axisymmetric setting. The expressions of the element sti ness matrices necessary for linearizing the resulting nonlinear system can be found in Msekh et al. (2015). 3. Benchmark tests: MPJR method for axisymmetric contact problems The MPJR method described in Sec. 2.1 embeds the spherical (in 3D) or circular (in 2D axisymmetric problems) indenter geometry into the interface elements, in order to retrieve the exact solution of the actual problem without the need of explicitly discretizing its shape. For a benchmark test, the reader is referred to (Bonari et al., 2022) for the Hertzian contact problem between a cylinder and a plane, with friction, under plane strain conditions. In the present work, the methodology has been extended to 2D axisymmetric contact problems with spherical indenters, with the aim of studying the problem of indentation-induced cracks in the substrate. The contact solution obtained with the MPJR method, for a radius of the sphere Rs= 10 mm, is compared with: (i) the Hertzian analytical solution for a half-space; ( ii) the numerical solution obtained by explicitly discretizing the spherical indenter following exactly its shape and solving the contact problem with the standard node-to-segment contact algorithm, with the same value of the penalty parameter. The FEM models used for the standard spherical contact problem and the MPJR method are presented in Fig. 9. Figure 9: Finite element model with the indenter explicitly discretized for the standard node-to-segment algorithm on the left and the model with the MPJR interface nite elements on the right. The same substrate mesh has been used in both models. The simulation has been carried out under displacement control, with a maximum displacement perpen- dicular to the substrate of 0 :01 mm, achieved in 100 steps. The substrate size (4 mm wide and 10 mm deep) is large enough with respect to the contact radius according to Hertz's theory. Regarding the mechanical properties, the Young's modulus of the substrate is 63 :40 GPa, and the Poisson ratio is 0 :2, while the in- denter Young's modulus has been set 100 times higher than that of the substrate to simulate a rigid one. 12 No damage has been considered for the solids in this benchmark test. The penalty parameter has been set equal to 1109N/mm. Results of the comparison shown in Figs. 10 and 11 highlight the excellent performance of the MPJR method, which also outperforms the standard node-to-segment contact algorithm, which presents some larger error for the rst pseudo-time step, due to the contact search algorithm which is avoided in the MPJR formulation. The solution obtained with the MPJR method for two uniform FE discretizations and di erent mesh sizes, he, shows that the accuracy can be improved by mesh re nement, as shown in Fig. 10 where the solutions for elements size he= 0:01 mm and he= 0:005 mm have been compared. Figure 10: Comparison between the Hertz equation (black line) and the result of the axisymmetric simulations obtained with a standard node-to-segment contact algorithm and the MPJR interface nite elements and two di erent mesh size, he. Figure 11: Reaction force error with respect to Hertz solution obtained using a standard FEM axisymmetric model with a node-to-segment contact algorithm and the MPJR method with the same mesh size. 4. Simulation of Hertzian cone cracks for smooth or rough spheres The Hertzian indentation test has been simulated by exploiting the axial symmetry of the model in a quasi-static framework. The model geometry and boundary conditions are sketched in Fig. 12, which 13 compares the real geometry of the test with the model geometry which employs the MPJR interface nite elements embedding the actual spherical shape in the at indenter. The model of the substrate consists of a 2510 mm rectangular domain as in the experimental studies in (Conrad et al., 1979; Jyh-Woei et al., 1993), which exploits symmetry conditions along the vertical axis on the left. Symmetry a (a) Real geometry. Interface with embedded sphere pro ileNominally lat indenter (b) Model geometry. Figure 12: Real geometry (a) and equivalent model geometry based on the MPJR interface nite elements embedding the actual shape of the spherical indenter (b). Thanks to the MPJR methodology, the indenter has been modeled as a nominally at punch using a layer of interface nite elements embedding the exact spherical pro le as an analytical function of the interface coordinates: z(x) =Rsp R2skxk, whereRs= 1 mm. The simulation is carried out under displacement control up to a value leading to a maximum contact radius asmaller than the radius of the sphere. Therefore, for this model, the portion of the interface embedding the spherical pro le has been set equal to 0:5 mm, corresponding to the size of the re ned FE discretization in Fig. 13, which is also much larger than the maximum contact radius reached in the simulation. Figure 13: Discretization of the indentation test with the magni cation of the re ned area having a characteristic elements size he= 0:0025 mm. The material properties of the soda-lime glass substrate have been taken from (Conrad et al., 1979) and reported in Tab. 1. The fracture energy of the glass has been estimated to 0 :009 N=mm in (Mouginot and Maugis, 1985). A discussion on the appropriate value for l0and the e ect of the choice on the crack pattern can be found in Strobl and Seelig (2019), (Wu et al., 2022), where di erent values of l0and glass tensile strength ranging between 50 MPa and 150 MPa have been tested in the simulations and compared with experimental data for the case of at-ended cylindrical indenters. The authors showed that the critical displacement required for the crack nucleation increases with c, while the ring radius r0decreases when the tensile strength increase. 14 Moreover, according to Strobl and Seelig (2020), the length scale parameter has to satisfy another criterion: it has to be chosen small enough to capture the extension of a spontaneous ring crack on the surface given in ( l00:04 mm). Considering both aspects, the length scale parameter has been set equal to 0 :01 mm; it corresponds to a tensile strength of the glass of 77 MPa according to the well-known formula for the AT2 phase- eld approach to correlate the strength to the other model parameters: c=q 27 256EGc l0. As shown in the following paragraph, the value l0= 0:01 mm allows a good reproduction of the experimental trends, and the phase- eld discretization of the model has been re ned where cracks are expected to nucleate (element size hel0=4 = 0:0025 mm). Young's Modulus Poisson Ratio Fracture energy Tensile strength Length scale Glass 63:40 GPa 0:20 0:009 N=mm 77 MPa 0:01 mm Table 1: Mechanical and fracture properties of the substrate taken from Conrad et al. (1979); Mouginot and Maugis (1985). As already stated, we remark here that, even though the model in Fig. 18 shows a at indenter, the interface nite elements exactly embed the spherical pro le shown in Fig. 14. 0 0.1 0.2 0.3 0.4 0.500.050.10.15 Figure 14: Pro le of the spherical indenter with radius 1 mm. The vertical displacement eld plotted at di erent pseudo-time steps is shown in Fig. 15, with the expected Hertzian distribution correctly reproduced (see also Bonari and Paggi (2020) for more details on the analysis of the stress eld for this benchmark solution). The crack develops with its typical conical shape, as shown in Fig. 16. The contour plots presented in this section show only a portion of the entire domain to evaluate the crack pattern better. In order to clarify this aspect, the proportion between the dimension of the entire domain and the detail of the crack pattern is shown in Fig. 17. The fracture evolution at di erent pseudo-time steps is shown in Fig. 18 for di erent values of the imposed far- eld displacement  u. 15 (a) u= 5:2103mm (b) u= 5:4103mm (c) u= 6:5103mm (d) u= 8:4103mm Figure 15: Snapshots of the contour plots of the vertical displacement component for the indentation test with a sphere of radiusRs= 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 nite elements with =1 have been removed, for ease of visualization, at a far eld 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 magni ed view. (a) u= 5:2103mm (b) u= 5:4103mm (c) u= 6:5103mm (d) u= 1:102mm 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- eld variable evolution along the interface is shown in Fig. 19(b). The contact pressure and the phase- eld 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 de ned 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 rst point of the interface where the phase- eld 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- eld variable in (b) along the interface at di erent 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 e ect of the sphere radius on the critical load causing crack nucleation has been investigated by considering ve cases with a radius Rsvarying from 1 to 7 :5 mm. The introduction of the interface nite elements allows changing the spherical indenter radius without changing the model geometry since the modi ed spherical geometry is embedded into the nite elements and analytically de ned 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 ratior0=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 con rm 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 at-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 di erent sizes (50 mm  50 mm25: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 con rm 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 di erent types of indenter, see Fig. 22b. It has to be noted that in the case of a at-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) in uences 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 experi- mental data in (Jyh-Woei et al., 1993) on as-received glass substrates. (b) Indentation of at-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 vsthe indenter radius. The experimental data in (b) concern indentation tests on glass specimens with di erent mechanical and geometrical properties. (a) Results of the numerical simulations compared with the experimental data in Conrad et al. (1979). (b) Indentation of at-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 vsthe inverse of the contact radius. To further exploit the advantage of the present approach, which enables to embed any indenter pro le along a nominally at interface and eciently solve the contact problem, we now consider the e ect of surface roughness on cone indentation fracture, using the same axial symmetric nite element model for the smooth sphere, but modifying the embedded indenter pro le. Generally, a real rough surface does not present axial symmetry; however, a 3D indentation model requires a high computational e ort compared to the axisymmetric con guration. 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 con guration has been chosen for the simulations in this work. To the best of the authors' knowledge, even the 2D axisymmetric solutions represent the rst attempt to investigate the e ect 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 pro le has been extracted and superimposed to the spherical shape of the smooth indenter (see the pro le 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 pro le chosen for the simulation has the statistical parameter Rzequal to 1:3m, which measures the average peak-to-valley distance of the pro le. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5-6-4-2024610-4 (a) Pro le obtained with RMD, with Rz= 1:3m and sampling 2 :5m. (b) Final pro le of the indenter embedded into the interface. Figure 23: Pro le of the rough spherical indenter (b) embedded into the MPJR interface nite elements, obtained by superim- posing the rough pro le (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 di erent points on the glass surface before the main crack propagates. (a) u= 9:2103mm (b) u= 9:3103mm (c) u= 9:6103mm (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- eld 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 di erent time steps until the phase- eld variable is 1 at the interface coordinate x= 0:175 mm, see the blue curve corresponding to the imposed displacement  u= 9:2103mm. The main crack does not propagate from that point but from the last peak of the phase- eld 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 radiusRs= 5 mm, in terms of contact pressure and phase- eld 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 pro les. The rough pro le presented in Fig. 23 has been rescaled to simulate three di erent maximum peak-to- valley distances Rz= 1:3m,Rz= 2:6m, andRz= 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 pro les are shown. For each case, the contour plots are taken at the time step corresponding to the crack propagation. The presence of roughness also a ects 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- eld variable for di erent pseudo-time steps up to the occurrence of the rst point along the interface with 1 (blue curve, corresponding to a displacement u= 9:2103mm). 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 pro le with Rz= 5:2m. 23 (a) (b) Figure 26: Rough spherical indenter with Rz= 1:3m: (a) contact pressure; (b) phase- eld variable for di erent pseudo-time steps from the occurrence of the rst point along the interface with 1 (blue curve), up to the situation corresponding to crack propagation for a higher displacement equal to  u= 9:3103mm (red curve). 24 (a) (b) Figure 27: Rough spherical indenter with Rz= 1:3m: (a) contact pressure; (b) phase- eld variable for di erent pseudo-time steps from the occurrence of crack propagation (red curve) onward, up to an imposed displacement equal to  u= 0:01 mm. 25 (a) Contact Pressure comparison. (b) Phase- eld variable comparison. Figure 28: Comparison between smooth and rough spherical indenters ( Rs= 5 mm) at the onset of crack nucleation, respectively at the imposed far- eld displacements of  u= 0:92 mm and  u= 0:93 mm in the two cases. 26 (a) Smooth spherical indenter,  u= 0:0089 mm. (b) Rough spherical indenter with Rz= 1:3m, u= 0:0093 mm. (c) Rough spherical indenter with Rz= 2:6m, u= 0:012 mm. (d) Rough spherical indenter with Rz= 5:2m, u= 0:022 mm. Figure 29: Radial position of the cone-shaped crack ( r0) on the surface of the glass substrate due to smooth and rough indenters with radius Rs= 5 mm and di erent amplitude of roughness, Rz. Figure 30: Reaction force vs. far eld imposed displacement for the smooth and rough indenters (radius of the sphere Rs= 5 mm). 27 5. Conclusion In this work, the simulation of cone-shaped cracks generated by spherical indentation tests has been successfully addressed by combining the MPJR interface nite elements to solve the non-conforming contact problem and the PF nite elements to predict nonlocal damage and fracture in the substrate. Coupling these two sources of nonlinearities, namely the contact mechanics problem and the fracture mechanics one, was one of the major challenges and contributions of the present work from the methodological point of view. In this regard, the MPJR interface nite elements have been demonstrated to be particularly ecient in treating complex interface pro les, as in the case of the simulated rough spherical indenters, where no previous solutions were obtained in the related literature. This computational framework allowed us to assess the indenter radius in uence on the indentation cracking. The comparison with the benchmark experimental results from Conrad et al. (1979) demonstrated that the current variational formalism enables reproducing the critical load for crack initiation and the major features of the physical problem. The in uence of roughness at the contact interface has been investigated by considering three di erent spherical rough pro les with increasing roughness amplitude. The comparison with the smooth case showed an increase in the critical load for crack propagation and in the ring crack radius by increasing roughness, in agreement with the experimental trends available in the literature. As an additional main insight into the problem, we obtained that the presence of roughness can induce higher stresses in the contact zone and therefore localized damage approaching = 1 in that area, which does not occur in case of smooth contacts. However, the origin of the propagating crack leading to cone fracture was always found outside the contact area. In terms of perspectives for future research, the present cases exploited the axial symmetry of the spherical indenter, which is computationally convenient. The framework is already suitable for addressing 3D models where the entire rough surface can be incorporated into the FE simulations, even though it requires a high computational e ort. A current line of research aims at developing a further HPC-enhanced framework to speed up the simulations. Moreover, the developed framework opens new research directions in the eld of indentation-induced fracture, also in the case of coated, FGM, or elasto-plastic substrates. CRediT authorship contribution statement M. R. Marulli : Methodology, Software, Investigation, Validation, Writing - Original Draft, Writing- Reviewing and Editing. J. Bonari : Software, Visualization, Writing- Reviewing and Editing. J. Reinoso : Conceptualization, Supervision, Writing- Reviewing and Editing. M. Paggi : Conceptualization, Supervi- sion, Original Draft, Writing- Reviewing and Editing. Acknowledgements The authors acknowledge the funding received from the European Union's Horizon 2020 research and in- novation program under the Marie Sk lodowska-Curie grant agreement No. 101086342 { Project DIAGONAL (Ductility and fracture toughness analysis of functionally graded materials; HORIZON-MSCA-2021-SE-01 action). Appendix To map the 3D nodal displacement vector u= (u1;v1;w1; 1;:::;u 8;v8;w8; 8)Tinto a global eld vector also including the nodal values of the phase- eld damage variable, ue= (u1;v1;w1; 1;:::;u 8;v8;w8; 8)T, the following matrix operator P3D(3224) is introduced: ue=P3Du; (28) 28 whose expression reads: P3D=2 666666666666666666666666664I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O13 I33O33O33O33O33O33O33O33 O13O13O13O13O13O13O13O133 777777777777777777777777775; (29) where I33is a (33) identify matrix, O33is a (33) null matrix, and O13is a (13) null matrix. For the 2D axisymmetric problem, the displacement eld u= (u1;v1;:::;u 4;v4)Tis mapped onto the vector ue= (u1;v1; 1;:::;u 4;v4; 4)Twith the following matrix operator P2D: P2D=2 66666666666666666641 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 03 7777777777777777775; (30) The above operators have the fundamental property PTP=Iand the following inverse relation holds: u=PTue.