Adhesive Boundary Element Method Using Virtual Crack Closure Technique

In this study, a new adhesive contact model is built upon a boundary element method (BEM) model developed by Pohrt and Popov (2015). The strain energy release rate (SERR) on the edge of the bonding interface is evaluated using Virtual Crack Closure Technique (VCCT) which is shown to have better accuracy and weaker mesh-size dependency than the closed-form SERR formula derived by Pohrt and Popov. A composite delamination criterion is proposed for crack nucleation and propagation. Numerical results predicted by the present model are in good agreement with the analytical solutions of two classic problems, namely, the axisymmetric parabolic contact and the sinusoidal waviness contact in the plane strain condition. The model of Pohrt and Popov can achieve a similar accuracy for the axisymmetric parabolic contact where the mesh grid is non-conforming to the crack front. Once the conforming mesh grid is used, the accuracy of their model is significantly deteriorated, especially at high work of adhesion and high mesh density. In both BEM models, however, the crack nucleation is found to be mesh-dependent which may be solved by introducing an upper limit for the tensile normal traction.


INTRODUCTION
Two mating surfaces bond to each other through intermolecular attractions, even under the tensile loading. In the idealized environment where the effect of humidity, surface charges, and contaminants can be neglected, the short-range van der Waals force is dominant in the intermolecular attractions. Besides the classic thermal dynamic approach used in the pioneering work of Johnson, Kendall, and Roberts (Johnson et al., 1971) (JKR theory), fracture-adhesion analogy plays a dominant role in finding the analytical solutions in the theoretical study of adhesive contact (Maugis and Barquins, 1978;Greenwood and Johnson, 1981;Maugis, 1992;Johnson, 1995;Persson, 2002;Carbone and Mangialardi, 2008;Xu et al., 2014;Ciavarella, 2015;Menga et al., 2016;Ciavarella et al., 2018;Ciavarella et al., 2019;Jin and Yue, 2020). If the intermolecular attraction outside the contact area is neglected (i.e., the JKR limit), the adhesive interface is equivalent to a brittle crack whose static equilibrium is governed by the Griffith's criterion (Maugis and Barquins, 1978;Greenwood and Johnson, 1981;Johnson, 1995). The intermolecular attractions outside the contact area can be included through the cohesive zone modeling, and the closed-form solution may be obtained if simplified cohesive laws are used [e.g., the Dugdale model (Maugis, 1992;Ciavarella et al., 2019;Jin and Yue, 2020) or the double Hertzian/Westergaard model (Greenwood and Johnson, 1998;Jin et al., 2016)].
The theoretical models can only be applied to few problems with ideal contact geometries, e.g., parabolic (Maugis, 1992) and sinusoidal waviness surfaces (Johnson, 1995). Better accuracy can be achieved using numerical models [e.g., BEM (Greenwood, 1997;Feng, 2000;Rey et al., 2017;Bugnicourt et al., 2018;Ghanbarzadeh et al., 2020) and Green's function molecular dynamics (GFMD) (Pastewka and Robbins, 2014;Müser, 2016;Khajeh Salehani et al., 2020;Wang et al., 2021)] where intermolecular attractions are explicitly quantified using various potentials and phenomenological cohesive laws. In other numerical models, the intermolecular attractions are introduced through the fracture-adhesion analogy where the contact edges are assumed as a Griffith crack. Carbone et al. (Carbone and Mangialardi, 2008;Carbone et al., 2009;Carbone et al., 2015) solved the adhesive contact problem between a thin elastic layer and a rigid rough profile using BEM where the locations of contact edges can be determined based on Griffith's criterion. Pohrt and Popov (Pohrt and Popov, 2015) developed a BEM adhesive contact model which is validated by various adhesion tests (Popov et al., 2017;Popov et al., 2021). The SERR in their model has a closed-form, which was derived based on the Boussinesq solutions. By guessing an overestimated real contact area, the adhesive contact problem in their model becomes an interfacial crack problem where the detachment is governed by a normal traction-based meshdependent criterion. The solution is converged until the crack fronts (i.e., contact edges) are rested.
In this study, a new adhesive contact model based on the work of Pohrt and Popov (Pohrt and Popov, 2015) is developed. The SERR is evaluated using VCCT which is a universal technique for evaluating the SERR in finite element method (Krueger, 2004). New algorithms are proposed for the delamination of the contact edges and the crack nucleation inside the contact area. The rest of this paper is structured as follows. In Section 2, the work of Pohrt and Popov is briefly introduced, including the main BEM algorithm, the mesh-dependent SERR formula, and the delamination criterion. In Section 3, the accuracy of the mesh-dependent SERR formula is explored by revisiting two classic mode-I crack problems. In Section 4, the main algorithm of the present model is given in detail. In Section 5, three classic adhesive contact problems are revisited by the present and previous BEM models.

THE PREVIOUS MODEL
Consider two linear elastic half-spaces in the purely normal contact under a normal load F acting at the far end. The Young's modulus and Poisson's ratio of two bodies are, E i and ] i , where i 1, 2. The geometries of two contact surfaces can be represented by (x, y, h i (x, y)). As proved by Barber (Barber, 2003), this problem is equivalent to a rigid flat in contact with an elastic half-space where the interfacial geometry has a composite form h (x, y) h 1 (x, y) + h 2 (x, y), see Figure 1. The surface deflection linearly depends on the plane strain modulus E* The geometrical relation at the contact interface has the following form where g is the interfacial gap; Ω is the computational domain with a finite size over z 0 plane; u z is the normal surface deflection which can be determined by the convolution between the Boussinesq solution and contact pressure distribution p (x, y) (Johnson, 1987); d is a load-dependent distance. The solutions at the contact interface, p (x, y) and g (x, y), must satisfy the following boundary conditions where Ω c and Ω nc are the contact and out-of-contact regions, respectively. Since g is zero within the contact area, d can be determined by d 1 h(x, y) − u z (x, y) dxdy where real area of contact A r is the size of the contact region, Ω c . Finally, the load equilibrium in the normal direction must be maintained F Ωc p(x, y) dxdy.
The original work of Pohrt and Popov in (Pohrt and Popov, 2015) is slightly different from the one given above due to its indentation-driven condition. The above formulation, in contrast, adopts the normal load-driven condition. Their model can be divided into two parts, a normal contact BEM solver and a mesh-dependent delamination criterion. The former part solves Eqs. 1-4 iteratively with a given contact region Ω c using the conjugate gradient (CG) method (Polonsky and Keer, 1999;Pohrt and Li, 2014). The continuous convolution-fast Fourier transform (CC-FFT) and discrete convolution-fast Fourier Transform (DC-FFT) (Liu et al., 2000;Wang et al., 2020) can be utilized to accelerate the calculation of the normal deflection, u z , for periodic and non-periodic domains, respectively. The resultant p (x, y) is a mixture of tensile and compressive tractions. The later part delaminates the elements where the normal tractions violate the Griffith criterion. Pohrt and Popov (Pohrt and Popov, 2015) creatively developed an analytical solution of SERR (G) using the Boussinesq solutions, and their derivation results in a meshdependent form (Pohrt and Popov, 2015) G x and Δ y are mesh sizes in the x and y directions, respectively; Δ Δ 2 x + Δ 2 y . The analytical solution in Eq. 5 quantifies the released strain energy due to the delamination of a rectangular element per unit area. The Griffith's criterion states that crack propagates when G > w where w is the work of adhesion. Substituting Eq. 5 into G > w, the energy-based delamination criterion can be rewritten as a normal traction-based one, i.e., p < − Σ where For a given Ω c , the corresponding normal traction, p (x, y), is solved using Algorithm 4 in Appendix A. In each iteration, the potentially detached elements are identified using the normal traction-based criterion in Eq. 6, and are considered as part of Ω nc . Implementation of the above delamination criterion is given in Algorithm 5. This process iterates until the convergence criterion is satisfied. To initiate this iterative process, an initial guess of the contact region Ω max c is used of which the true contact area is a subset. Clearly, this BEM model can only be applied during the unloading stage. It is impossible to capture the snap-in instability. In a recent work of Popov et al. (Popov et al., 2021), the model of Pohrt and Popov (Pohrt and Popov, 2015) has been extended to cover the loading stage with the capability of capturing the snap-in instability. The detailed implementation of the iterative process can be found in Algorithm 3. Neglecting the energy dissipation due to the unstable snap-in/pull-off at asperity level (Greenwood, 2017;Popov et al., 2021), this BEM model is valid for adhesive contact in both loading and unloading stages (i.e., loading and unloading stages are assumed to be reversible).

Two Classic Crack Problems
The mesh-dependent SERR formula is essential to the previous BEM model, which is responsible for an accurate delamination criterion. Even though their model has been validated in various numerical and experimental studies (Popov et al., 2017;Popov et al., 2021), the accuracy of Eq. 5 has never been proved. In this section, two classic mode-I crack problems are revisited, namely.
• Problem 1: the mode-I periodic collinear cracks of semiwidth a with period λ embedded in an infinite space, which are subjected to a tensile normal traction p at far end, see Figure 2A. This classic problem was solved analytically by Koiter (Koiter, 1959), and the SERR has the following form (Johnson, 1995;Tada et al., 2000) G analytical 1 p 2 λ 2E * cos πa λ (7) • Problem 2: a penny-shaped mode-I external crack of radius a embedded in an infinite space, which is subjected to a tensile force F at far end, see Figure 2B. The SERR has the following form (Maugis, 1992;Tada et al., 2000) G analytical 2 Two problems are all symmetric about its crack plane (see the darker plane in Figure 2). Therefore, two mode-I crack problems are equivalent to the following adhesive contact problems: • problem 1: a rigid flat and a half-space where half periodic collinear cracks are embedded at the bonding interface; • problem 2: a rigid flat and a half-space where a half external penny-shaped crack is embedded at the bonding interface.
Because the bonded area (contact area) is known in advance, the CG method given in Algorithm 4 can be used to solve the normal traction, p (x, y), within the bonded area. The SERR is evaluated over the elements at the contact edge using Eq. 5, as well as the VCCT. A short introduction of VCCT is given in Section 3.2.

Virtual Crack Closure Technique
Virtual crack closure technique (VCCT) is widely used in 2D and 3D finite element model for determining SERR (Krueger, 2004). Consider a simple case where a flexible surface is peeled off from a FIGURE 2 | Schematics of (A) periodic collinear cracks of period λ with semi-width a subjected to a uniform tensile normal traction p at far end, and (B) a penny-shaped external crack of radius a embedded in an infinite space subjected to a normal load F at far end. Only a finite portion of the infinite space is illustrated.
Frontiers in Mechanical Engineering | www.frontiersin.org October 2021 | Volume 7 | Article 754782 rigid flat. The normal traction, p, and crack opening displacement, g, over the bonded and debonded surface elements, respectively, are shown in Figure 3. The discontinuous crack opening displacement and contact pressure distribution are due to the constant elements used in the BEM model (uniform normal traction and interfacial gap within each element). As the crack front propagates from A (at load step k) to A′ (at load step k + 1) by a sufficiently fine element size Δ x , the following self-similarity can be approximately satisfied: Figure 3. This self-similarity enables the determining of SERR within one load step. At step k, the amount of strain energy released due where the superscript (k) is dropped for the sake of simplicity.
For an arbitrary crack front that cannot be perfectly parallel with the element edges (see Figure 4A for an example), the crack opening displacement in Eq. 9 may have more than one candidate to choose among the closest neighbors of element (i, j). Some methods are proposed in the finite element literature to solve this

Normal traction
Interfacial gap

Normal traction Interfacial gap
Bonded area Debonded area ambiguity for non-constant (linear, quadratic, etc) elements, and a comprehensive summary of those methods was given by Wu et al. (Wu et al., 2021). For a constant element, we propose the following simple formula where max (g i j±1 , g i±1 j ) indicates the maximum crack opening displacement at four closest neighbors of element (i, j). The reason for picking the maximum interfacial gap is because that the edge associated with the maximum interfacial gap is likely the weakest edge for the potential crack to propagate through. Unlike the mesh-dependent SERR formula in Eq. 5 which is based on the Boussinesq solutions (Pohrt and Popov, 2015), VCCT does not rely on any specific analytical solutions so that it is a universal method for evaluating SERR in mode I/II/III, regardless of the material properties (elastic, viscoelastic or elastoplastic) and types of the domain (periodic or non-periodic).

Strain Energy Release Rate Results
In problem 1, the conforming mesh is used where the element edges are parallel with the crack front (see the straight crack front in Figure 3 as an example). Since the crack problem is in the plane strain condition, the variation of SERR along the crack front is negligible. The SERR determined by VCCT, Eq. 10, associated with three different mesh densities all have excellent agreement with the analytic solution G analytic 1 , see Table 1. Equation 5 is derived based on the non-periodic, point load (Boussinesq) solutions under the three-dimensional (3D) elasticity. It is still valid to apply Eq. 5 to evaluate G in problem 1. Firstly, the 3D BEM model is inherently in the plane strain condition due to the usage of two-dimensional FFT in the xy plane so that the crack has an infinite length in the y direction; Secondly, the semi-width of the crack, 0.05 (mm), is negligible compared with the period λ 2.5 (mm), the coupling between neighboring cracks can be neglected. Thus, each collinear crack can be studied individually. Finally, the element sizes, Δ x and Δ y , are extremely small compared with the period. Therefore, it is valid to apply Eq. 5 to any periodic problem as long as the mesh density is sufficiently high. Surprisingly, the SERR determined by the mesh-dependent formula in Eq. 5 is greatly underestimated by a maximum of 47% associated with mesh grid of 1,024 × 1,024, see Table 1.
Additionally, the SERR determined by Eq. 5 decreases by more than 10% with respect to G analytic 1 when the mesh grid is changed from 256 × 256 to 1,024 × 1,024. The variation of the SERR determined by VCCT, in contrast, is only about 1%. This implies that the mesh-dependent formula has a stronger meshdependency compared with VCCT when applied to the conforming mesh. Equation 10 associated with the minimum gap, min (g i j±1 , g i±1 j ), and the average gap, mean (g i j±1 , g i±1 j ), are also tested, and the corresponding mean SERRs are nearly the same.
The mesh grid in problem two is non-conforming, since element edges do not always align perfectly with the crack front (see the non-smoothed crack front highlighted in blue in Figure 4A as an example). The distributions of SERR over the circular crack front G(θ), determined by both VCCT and meshdependent formula, oscillate about mean levels with the amplitude as high as 50% of G analytic 1 , see Figure 4B. The oscillation of G(θ) at the penny-shaped crack front commonly occurs in the finite element model due to the non-smoothed crack front (Wu et al., 2021). As we shown in Section 4, the accuracy of the SERR determined by VCCT, as well as the mesh-dependent formula, is related to an accurate estimation of the strain energy U and the virtually closed area per element. The accuracy of U is strongly related to the choice of crack opening displacement. Due to the fact that the crack is non-smoothed, crack opening displacement is ambiguous, and the virtually closed area per element is hard to determine. As an example in Figure 4A, some elements shown at the crack front are only closed partially. For an adhesive contact where the crack front is not known in advance, it is nearly impossible to correctly determine the virtually closed area per element. Therefore, the oscillation amplitude of G(θ) is high. Many VCCT formulations, similar to Eq. 10, were proposed to improve the accuracy of G(θ) associated with non-conforming mesh (Krueger, 2004;Xie and Biggers, 2006;Wu et al., 2021). The extent of oscillation may be slightly relieved, but the oscillation cannot be completely removed.
The SERR solved in problem two are summarized in Table 2 where the maximum, average, and standard deviation values of G(θ) all slightly reduce with the increase of the mesh density for both VCCT and mesh dependent formula. This is expected since the smoothness of the crack front is improved at a higher mesh density. However, the accuracy improvement is not significant TABLE 1 | The SERR of the crack problem 1 determined by mesh-dependent formula, Eq. 5 and VCCT, Eq. 10. Size of computational domain, Ω, is 2.5 × 2.5 (mm 2 ); tensile normal traction p 1 (MPa); λ 2.5 (mm); a 0.05 (mm); E* 1 (MPa). Three different mesh grids are used to discretize the computational domain, namely, 256 × 256, 512 × 512, and 1,024 × 1,024. Collinear cracks are distributed in the x direction. Only one period is modeled within the computational domain where two semi-cracks are located at the vicinity of x ±1.5 (mm). The periodic distribution of the collinear cracks in the x direction, and the plane strain condition (with infinite crack length in the y direction) are achieved simultaneously by using FFT in the CG method given in Algorithm 4.  of an internal penny-shaped crack determined by VCCT and mesh-dependent formula, i.e., Eq. 10 and Eq. 5. max (), mean () and std () are the functions of evaluating the maximum, average and standard deviation values of G′(θ ∈ [0, 360°)). Key parameters are the same as that in the caption of Figure 4.

Method
Mesh grid max (G9) mean (G9) std (G9) . This implies that the inaccuracy introduced by the non-smoothed crack front may be relieved if the SERR is determined as the mean of all SERR at the crack front. This averaging procedure can be applied to any purely normal adhesive contact problem, because the SERR at any point of the static contact edge is the same as the work of adhesion. The mesh-dependent formula captures a similar oscillation of G(θ) with the mean level about 40% less than G analytic 2 . Surprisingly, the maximum (see Table 2) and local maxima (see Figure 4B) of G(θ) determined by the mesh-dependent formula are approximately the same as G analytic 2 for all three mesh densities. The row "Maximum change" in Table 2 implies that Eq. 5 has a much stronger mesh-dependency than VCCT when applied to the non-conforming mesh. Equation 10 associated with the minimum gap, min (g i j±1 , g i±1 j ), and the average gap, mean (g i j±1 , g i±1 j ), are also tested, and the corresponding mean SERRs are nearly the same.

Remarks
According to Table 1 and 2, we can conclude that the meshdependent SERR formula, Eq. 5, cannot accurately estimate SERR associated with the conforming mesh nor the mean SERR value associated with the non-conforming mesh. Its accuracy becomes worse when a finer mesh is used. Pohrt and Popov assumes that the strain energy U due to the delamination of a rectangular area of size Δ x × Δ y is equivalent to the work done by indenting an elastic half-space within a rectangular area of the same size under the uniform normal traction of σ This may be correct for evaluating U associated with the crack nucleation within a bonded area, since its neighboring elements are all closed. It is not appropriate, however, to evaluate U of the delaminated element at the crack front where some of its neighboring elements have already been opened.
Surprisingly, the BEM model adopts the mesh-dependent SERR formula has been extensively validated by the analytical solutions and empirical data (Pohrt and Popov, 2015;Popov et al., 2017;Popov et al., 2021). The above-contradicted conclusions may be due to the following reasons: 1). As far as the authors know, the previous BEM model (Pohrt and Popov, 2015) has never been applied to the case where the contact edge is perfectly aligned with the element edge, e.g., adhesive contact in the plane strain/stress condition. Approximately 50% underestimation of SERR by Eq. 5 would result in an earlier crack rest, and cause an overestimation of the contact area. This expectation is confirmed by an example of sinusoidal waviness contact in the plane strain condition given in Section 5.2. 2). Instead of Eq. 5, its equivalent form in terms of normal traction, Eq. 6, is used in the previous BEM model (Pohrt and Popov, 2015). Because Σ ∼ G √ is proportional to the square root of SERR, the inaccuracy introduced by SERR is weakened. 3). Only the elements associated with high peak values of SERR are detached. According to Figure 4B and In a summary, the mesh-dependent formula is less accurate to evaluate the SERR of the elements at the crack front.

THE CURRENT MODEL
In this section, a new BEM model is built upon the previous one (Pohrt and Popov, 2015) to fix the inaccurate delamination criterion. Like the previous one discussed in Section 2, the current model can only be applied during the unloading stage which is impossible to capture the snap-in instability. The conjugate gradient method (Algorithm 4) is adopted in the current model to numerically determine the normal traction over a given contact area. The debonding of the contact interface can be divided into two stages, namely, crack nucleation and crack propagation. A composite delamination criterion is developed below to cover both stages.

Crack Nucleation
Crack nucleation occurs inside the bonded area at the vicinity of closed valleys (dimples). However, VCCT cannot be applied due to the absence of local crack fronts. Cohesive zone modeling (CZM) may be used to circumvent this difficulty. A typical failure of a mode-I cohesive crack consists of two stages, i.e., damage initiation and damage evolution. Before the damage initiation, the cohesive (tensile) normal traction p linearly increases with the interfacial separation δ where the proportionality K is the normal stiffness, see Figure 5. Beyond the damage initiation point (p c , δ c ), the cohesive interface enters the damage evolution (softening) stage where the normal stiffness degrades gradually until the complete failure is reached at δ δ f . The normal traction-based delamination criterion, Eq. 6, proposed by Pohrt and Popov (Pohrt and Popov, 2015) is an extrinsic cohesive law (Zhou et al., 2020) where the interfacial separation is strictly zero within the bonded area. Therefore, corresponding p(δ) curve before the damage initiation has an infinite slope K → ∞ (see the red thicker line in Figure 5). The maximum tensile stress, Σ, determined by Eq. 6 could be treated as the critical tensile stress |p c | Σ at the damage initiation point. For a brittle fracture, catastrophic failure occurs immediately after the damage initiation point. Thus the damage evolution stage could be neglected. Equation 6 also removes the need of assigning multiple parameters (e.g., K, δ c and δ f ) in the phenomenological laws. Algorithm 1 below illustrates how the crack nucleation is implemented in the present model.

Crack Propagation
The VCCT is used in the present model to calculate the SERR on the contact edge. Because the accuracy of the SERR is deteriorated by the non-smoothed crack front associated with the non-conforming mesh, it is not appropriate to identify the element delamination based on the local SERR values. Taking the stable penny-shaped external crack in Figure 4A as an example. If the delamination of nonconforming elements at the crack front are checked using G(θ) > w G analytic 1 where G(θ) is shown in Figure 4B, it is clear that this crack front is not stable and nearly half of elements at the crack front should be delaminated. The delamination continues until the max(G(θ)) < G analytic 1 . This contradicts the stable assumption and underestimates the bonded area.
To overcome this contradiction, a new crack propagation criterion is proposed based on the mean SERR: mean(G) > w where mean(G) is obtained by averaging the SERR values of all elements at all crack fronts. Once mean(G) > w, only the elements associated with larger SERR values are delaminated so that the mean value of the SERR of the residual bonded elements at the crack front is approximately the same as the work of adhesion w. This delamination process can guarantee that no crack propagation at any individual elements once mean(G) ≤ w is satisfied. The implementation of this delamination process can be found in Algorithm 2.
The algorithm of the present BEM model is almost the same as that given in Algorithm 3.

RESULTS AND DISCUSSION
Two classic adhesive contact problems, namely, axisymmetric parabolic contact and sinusoidal waviness contact in the plane strain condition, are revisited by the previous and present models.

Axisymmetric Parabolic Contact
Consider a rigid parabolic indenter of radius R in a purely normal contact with an elastic half-space under a normal load of F, see the inset in Figure 6B for the schematic. Ignoring the intermolecular attractions outside the contact area, this classic problem can be solved using JKR theory (Johnson et al., 1971;Maugis, 1992). The closed-form solutions of the radius of adhesive contact area, a, and the indentation depth, δ, are For a given normal load F, the corresponding contact pressure distribution p (r < a) inside the contact area is and the interfacial gap distribution g (r > a) outside the contact area is The radius of contact area a(F) and indentation depth δ(F) associated with a sequence of F (increasing from F min to 5|F min |) are solved by two BEM models and JKR theory where F min − 3 2 πwR is the pull-off force associated with the fixed load condition. a and δ are normalized by a 0 and δ 0 which are the corresponding radius a 0 ( 9 2 πwR 2 E * ) 1/3 and indentation depth δ 0 a 2 0 /R − 2πwa 0 /E * at zero normal load, respectively. An initial guess of the contact area completely contains the largest contact area under the maximum normal load. The critical tensile stress p c is not necessary in this problem, since crack front exists in the initial guess of contact area and no crack nucleates inside the bonded area.
Excellent agreement of a(F) and δ(F) predicted by two BEM models and JKR theory can be found in Figure 6A. A similar agreement of contact pressure distribution p (r < a) and interfacial gap distribution g (r > a) under a normal load of F 5|F min | is shown in Figure 6B. In Table 1 and 2, the SERR determined by a mesh-dependent formula adopted in the previous model is found to be strongly influenced by the mesh size. To investigate how this mesh-dependency influences the results of the previous model, a mesh convergence test is conducted and the relative percentage difference of a and δ with respect to the JKR solutions (i.e., a JKR and δ JKR ) are given in Table 3. The relative percentage difference of a and δ are negligible. a and δ predicted by both BEM models monotonically converge to the JKR solutions as the mesh density increases. Numerical data in Table 3 confirms that the error brought by the meshdependent SERR formula can be neglected in the previous model when the non-conforming mesh is used.

Plane Strain Sinusoidal Waviness Contact
Consider a rigid periodic sinusoidal waviness in the adhesive normal contact with an elastic half-space, see Figure 7. The sinusoidal profile z h (x, y) has the following geometry Dimensionless coordinate Dimensionless load FIGURE 6 | (A) Dimensionless contact radius a(F)/a 0 and dimensionless indentation δ(F)/δ 0 ; (B) Contact pressure p (x, 0) and interfacial gap g (x, 0) at F 5|F min | where the contact radius is a max . E* 1 (MPa); L x L y 2.5 a max ; n x n y 128; w 10 × 10 -6 (mJ/mm 2 ); R 10 (mm).
TABLE 3 | Results of mesh convergence test of the present and previous BEM models. All key parameters are the same as that give in the caption of Figure 6. Normal load F 5|F min |. a JKR 0.3478 (mm) and δ JKR 7.4188 × 10 -3 (mm) are solutions of JKR theory.

Mesh
(a − a JKR )/a JKR × 100 (%) ( δ − δ JKR )/δ JKR × 100 (%) Present model Pohrt and Popov (2015) Present model Pohrt and Popov (2015) 128 Ignoring the intermolecular attraction outside the contact area, Johnson (Johnson, 1995) solved this problem analytically by superposing a collinear periodic crack problem on a normal non-adhesive sinusoidal waviness contact problem. The contact width 2a within one period λ can be determined from the following non-linear equation where ψ a πa/λ, p* πE*Δ/λ, and α 2E * w λ 1 p* . Differentiating the left hand side of Eq. 17 with respect to ψ a , we can have a polynomial of tan(ψ a /2) which results in either two distinct real roots (stable and unstable branches), two identical roots (instability point) or no real root. The contact pressure p (x ∈ [− a, a]) within the bonded zone is 1 p(x) 2 p ′ cos(ψ) sin 2 (ψ a ) sin 2 (ψ a ) − sin 2 (ψ) where p ′ p* sin 2 (ψ a ) and p ′′ −p* α tan(ψ a ) . The interfacial gap, g (x′ ∈ [ − l, l]), can be derived in an integral form based on the fracture mechanics approach developed by Xu and Jackson where x′ is the local coordinate centered at the trough of the sinusoidal waviness, i.e., x′ x + λ/2, see Figure 7 and l λ/2 − a. The Green's function G P uy (x ′ , x ′′ ) has the closed-form in Ref. , and the average interfacial gap is g 1 L l −l g(x ′ )dx ′ . According to Johnson's solution, the loading and unloading stages are not reversible. Therefore BEM (VCCT) can only be applied to the unloading stage. As pointed out by Johnson (Johnson, 1995), the sinusoidal waviness bonding interface needs an infinite tensile force to separate if complete contact occurs. Following Johnson's suggestion, we assume there is a flaw (non-contact region) of width 2b f ≈ 0.11λ which is symmetric 1 Several typos are found in Johnson's original derivation (Johnson, 1995). The contact pressure p(x) is decomposed into p′(x) and p″(x) terms. The former term is Westergaard solution, but the square root is missing in Eq. 5 in Ref. (Johnson, 1995). The latter part contains an extra negative sign in front of p ′′ , see Eq. 6 in Ref. (Johnson, 1995 about x′ Nλ where N 0, ±1, ±2, /. The initial interfacial gap is illustrated using blue dashed line in Figure 7 where the compressible air with negligible internal pressure is trapped within the initial non-contact region. Substituting b f into Eq. 17, the critical average contact pressure p c can be derived p c p* cos 2 πb f /λ − α cot πb f /λ As p > p c , the contact width remains the same. The trapped volume, as well as the average interfacial gap, gradually increases as p drops. As p drops below p c , an unstable point is reached where the contact width a instantaneously reduces and quickly rests at a stable value. Similarly, the average interfacial gap g instantaneously increases and stabilized at a larger value, see Figure 8C for a clear illustration. Following the stable branch of the unloading curve, jump-off occurs when p is further reduced to a minimum value where the local slope of a( p) → ∞. Beyond that point, the sinusoidal waviness is completely detached from the rigid flat.
Figures 8A-C indicates a good agreement in the unloading stage between the numerical solutions of contact width, 2a, and the average interfacial gap, g, solved by the present model and Johnson's solutions until jump-off contact occurs for all three values of α. The numerical solutions of the previous model, however, gradually deviate from Johnson's solution as α increases. Contact pressure p(x) and interfacial gap g(x) predicted by the present model at α 0.2 are also in good agreement with Johnson's solution within both the tensile and compressive regime, see Figure 8D. The previous model underestimates the interfacial gap and overestimates the contact pressure, and the accuracy is even worse associated with higher α. The sinusoidal waviness contact problem clearly shows that the previous model is not valid to be applied to the plane strain contact problem where the mesh is conforming.

Discussion
In Section 3.3, the accuracy of the present model is validated by two classic adhesive contact problems. The previous model is found to be sufficiently accurate only when a non-conforming mesh is used. Up till now, we have shown that 1) VCCT is suitable for evaluating the SERR in the purely normal adhesive contact problem; and 2) the composite delamination criterion is correctly implemented for the crack propagation stage. The crack fronts initially exist in the above two classic problems, thus the crack nucleation part of the composite crack delamination criterion has not been tested yet. The following 3D sinusoidal waviness contact problem is chosen to explore the composite delamination criterion at the crack nucleation stage.
Consider a rigid 3D sinusoidal waviness in a purely normal adhesive contact with an elastic half-space. The sinusoidal surface z h (x, y) has the following geometry The contact pair is subjected to an average contact pressure of p. Under the JKR limit, only the asymptotic solutions are available (Johnson, 1995). When the rigid waviness is in complete contact with the elastic half-space, the corresponding normal traction p c (x, y) is (Johnson, 1995;Xu et al., 2015) p c (x, y) p* cos 2π λ x x cos 2π λ y y + p where p* πE * Δ λ −2 x + λ −2 y When the wavy surface is in a complete contact (solid interaction occurs everywhere on z 0 plane) with the elastic half-space, it is subjected to an average contact pressure p ≥ p*. If the wavy surface is unloaded from the complete contact status, crack nucleation occurs at the vicinity of the valley where maximum tensile traction ( p − p* < 0) is located. According to the composite delamination criterion, crack nucleates at the trough of the waviness once the following inequality is satisfied At this instance, crack nucleates and instantaneously propagates until it rests on a stable branch, see Figure 9. The prediction of contact ratio predicted by the two BEM models shows an excellent agreement. As the work of adhesion w increases, the crack nucleation | Contact ratio to average contact pressure relation associated with multiple work of adhesion w. E* 1 (MPa), λ x λ y 10 (mm), Δ 0.1 (mm), n x n y 128, and w 0 ∼ 60 × 10 -6 (mJ/mm 2 ).
TABLE 4 | A mesh convergence study of two BEM models. p/p* 0.5, w 10 × 10 -6 (mJ/mm 2 ), and the other parameters are the same as that given in the caption of Figure 9.

Mesh grid
Present model Pohrt and Popov (2015)  occurs at lower p. This is due to the fact that maximum tensile stress Σ given in Eq. 6 monotonically increases with w.
To explore the influence of mesh size on the solutions of BEM models where crack nucleates inside the bonded area, the sinusoidal waviness contact problem is repeatedly solved using three different mesh densities, and the selected results are given in Table 4. When mesh density is low (e.g., 128 × 128), crack nucleation occurs under the average contact pressure of p 0.5p*. As density increases to 256 × 256 and 512 × 512, the interface remains closed. This mesh-dependent crack nucleation phenomenon is due to the fact that Σ is also meshdependent. Consider a special case where Δ x Δ y , the maximum tensile stress Σ in Eq. 6 simplifies to (Pohrt and Popov, 2015): which indicates a monotonic increase of Σ with the reduction of element size Δ x . For a fixed normal load p, the maximum tensile stress at the complete contact (closure) is p − p* < 0. If the mesh size is smaller than a critical value, Σ would exceed | p − p*| and the interface remains closure. This mesh-dependent nucleation paradox may be solved by setting an upper limit for Σ, i.e., max(Σ) Σ c . However, the empirical approach of measuring Σ c is not clear. For the rough surface contact, two BEM models could result in unexpected closures in the vicinity of a great number of valleys at the unloading stage. Commonly, a large number of elements with small sizes are used to represent the fine details of rough surfaces. In those cases, both models tend to overestimate the interfacial toughness, as well as the real area of contact.
In this study, only the purely normal contact problem associated with mode-I crack is modeled. The present model can also be extended to other complex cases in mode-I, II, and III and/or with bi-material interface cracks, e.g., the purely normal contact considering the interfacial friction (Khajeh Salehani et al., 2018) and adhesive friction (Khajeh Salehani et al., 2019). The calculation of tangential surface displacement, as well as new crack propagation and nucleation criterion, must be added to the present model. The former is easily implemented using the Boussinesq-Cerruti solution and FFT. A typical crack propagation criterion in a mixed mode may be written as G f(G I , G II , G III , / ) > w where G is the total SERR which is composed of the SERRs in mode-I, II, and III, as well as other parameters, e.g., the mode mixity. SERRs of all modes are determined locally using VCCT. An averaging process is applied to get the mean value of G based on the local value of each element at the contact edge. For an absence of crack fronts, the critical state of interfacial stress may be obtained following the derivation of the mesh-dependent SERR formula where the strain energy contribution due to the tangential load must be added. For better accuracy, the authors suggest developing an empirical approach to curve-fit a phenomenological law for the critical state of interfacial stress at the instance when the crack is nucleated within the bonding area.

CONCLUSION
In this study, a new BEM model is developed upon the previous one proposed by Port and Popov. VCCT is used to evaluate SERR on the edge of the bonding interface. It is a universal tool independent of the material models and domain types. By revisiting two classic crack problems, it is shown that VCCT has better accuracy and weaker mesh-size dependency than the closed-form SERR formula adopted in the previous model. A composite delamination criterion is proposed for crack nucleation and crack propagation. Two classic adhesive contact problems, namely, the axisymmetric parabolic contact and the sinusoidal waviness contact in the plane strain condition, are revisited by the present and previous BEM models. Numerical results predicted by the present model have good agreement with the analytical solutions of both problems. The previous model can achieve a similar accuracy for the axisymmetric parabolic contact where the mesh grid is non-conforming. Once the mesh grid is conforming to the crack front, the accuracy of the previous model is significantly deteriorated, especially at high work of adhesion and high mesh density. In both BEM models, however, the crack nucleation is found to be mesh-dependent which may be solved by introducing an upper limit for the tensile normal traction.
As far as the authors know, the previous model has never been applied to any plane-strain contact problems where the mesh grid is strictly conforming to the crack front, the comparisons and conclusions made upon the results of the previous model are still valid. The authors recommend replacing the mesh-dependent SERR formula with the VCCT in the adhesive contact models for better accuracy, weaker mesh-dependency, and broader applications.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YX and RZ contributed to the conception and design of the study. YX performed the numerical analysis. YX wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.