Non-adhesive Contacts With Different Surface Tension Inside and Outside the Contact Area

In the past decade, the influence of surface tension on contact properties has attracted much attention, especially in the context of contact of very soft materials (such as gels) or contacts at the nanoscale. However, in the most current studies it is assumed that the tension of the surface inside and outside the contact area is the same. In practical terms, this means that the object considered is an elastic body “coated” with a tensed membrane. In real contacts, there is no reason why the surface tensions of the “free interface” and the “contact interface” should be equal. On the contrary, especially in contacts of soft bodies with hard solid indenters, one can anticipate that they are completely different. In the present article, we consider an elastic contact taking into account different surface tensions inside and outside the contact area. However, the considered contacts are still “non-adhesive.” This means that the three surface energies in play (two surface energies of both bodies outside the contact and the interface energy in the contact region) fulfill the criterion that the work of separation vanishes. A numerical model based on the Fast Fourier transform–assisted boundary element method is implemented and is illustrated with a few examples.


INTRODUCTION
It is known that surface tension governs the contact behavior of liquids (Brown, 1974) but it is hardly seen in contacts of solid bodies at the macro-and mesoscales. It is clear that it should start to play an essential role if we consider a continuous transition from an elastic solid to a very soft body such as a gel, soft rubber, and biological tissue (Style et al., 2017). In biological studies, it is found that surface tension may determine the tissue growth and the kinetics of the regeneration of organs (Ehrig et al., 2019). Surface tension may become important even in stiff contacts if their size is very small. In particular, it may influence the function of microelectromechanical systems (MEMS) (Syms et al., 2003).
There exists a class of contact problems in which the surface energy has been considered already for almost 50 years; these are adhesive contacts. While non-adhesive contacts of isotropic elastic bodies are completely characterized by the elastic modulus E and the Poisson's ratio ν of bulk material, in adhesive contacts the surface energy plays a key role, as has been shown in the classic theory by Johnson, Kendall, and Roberts (JKR) in 1971 (Johnson et al., 1971). However, when talking about the "surface energy, " we have always to specify what surface energy is meant. At the contact boundary, three interfaces are coming together. The "surface energy" used in the JKR theory is better characterized as "specific work of separation, " which we will also call "work of adhesion, " w.
This quantity was introduced already by Dupre in 1869 (Dupré, 1869); it can be expressed as w = γ 1 + γ 2 − γ 12 (1) where γ 1 and γ 2 are specific free surface energies of two bodies and γ 12 interfacial surface energy inside the contact area (Figure 1). From a general theoretical point of view, the JKR theory is incomplete: It makes a step by including surface energy into consideration but at the same time assumes that the surface energies of bodies outside the contact area are zero. However, there are no theoretical reasons why the surface energy of the interface should be so dominant compared to the surface energies of the surfaces outside the contact. In a general case, all three surface energies are different and may have a comparable order of magnitude.
As academic theoretical models, one can consider various relations of the three surface energies and corresponding limiting cases.
(I) If an elastic body is in contact with a rigid body and the specific surface energy of the elastic body outside the contact area can be neglected compared with the interface energy in the contact area, then we have an adhesive contact with specific work of separation w = γ 2 − γ 12 . It is this case that is considered in the well-known work by Johnson et al. (1971). In reality, this case can be realized only in a contact of a high surface energy solid with low surface energy polymer but is hardly realizable in a contact of two materials with comparable surface energies. (II) If the surface energy of the elastic body outside the contact area is finite, γ 1 = 0, but the work of adhesion is zero, then we have a non-adhesive contact with surface tension. (IIa) The latter case usually has been further simplified with the assumption that the surface energies inside and outside contact are equal: which automatically leads to the conclusion that γ 2 = 0. The corresponding theory was first developed for elastic foundations (Filonenko-Borodich, 1940;Kerr, 1964) and generalized to three-dimensional half-space contacts in Hajji (1978). It attracted much attention in recent years. However, this case cannot be realized physically, as there exist no substances with zero surface energy. In particular, in the cases that motivated consideration of this casecontacts of a hard indenter with a soft solid-one can assume that γ 2 not only does not vanish but is also the largest of the three relevant surface energies [ (Popov, 2017), Chapter 3]. This popular case can practically be realized only as a composite consisting of an elastic body coated with a tensed membrane. (IIb) If the condition of the absence of adhesion is fulfilled, but the validity of the usual (and physical completely unrealistic) condition is not assumed, then we have a "general non-adhesive contact with surface tension, " which is characterized by two independent surface energies. It is this case that will be considered in the present article. (III) Finally, in the general case both work of adhesion and surface tension of the free surface are finite. This leads to a general adhesive contact with surface tension. The latter attracted much interest in the past two decades in the context of indentation of soft matter (e.g., gels) (Carrillo and Dobrynin, 2012;Style et al., 2013;Cao et al., 2014).
It is interesting to note that it took almost 50 years to make the next necessary (and, as a matter of fact, obvious) step after the JKR theory. On the other hand, the appearance of the JKR theory, which is nothing but a direct implementation of the Griffith crack theory, also took 50 years (Popova and Popov, 2018). For the sake of historical truth, let us, however, mention that a solution equivalent to that of JKR was obtained by Sperling already in 1964 in his unpublished doctoral thesis (Sperling, 1964;Borodich, 2014).
In the present article, we will briefly recapitulate the studies of case IIa and then consider a more general (and more complicated) case IIb. As already stated above, case IIa can be physically realized only as an elastic body coated with a stressed membrane. The fundamental solution determining vertical displacement of the surface under the action of a concentrated force was given for this case by Hajji (1978). In this contact problem a characteristic length appears, γ 12 /E * , which is called "elastocapillary length, " determines the influence of surface tension: the surface tension effect is significant if the value of γ 12 /E * is comparable with the contact radius (Long and Chen, 2017). This leads to a "size-dependent" behavior in indentation testing (Style et al., 2013).
In the past few years, adhesive and non-adhesive contacts with surface tension have been studied intensively, for example, contact of homogeneous elastic half space including spherical contact (Long and Wang, 2013;Hui et al., 2015), conical contact (Long and Chen, 2017), and two-dimensional cylindrical contact , all based on the Hajji's fundamental solution, as well as complicated inhomogeneous coatings based on the semianalytical modeling of the surface effect (Zhang et al., 2018).
The present article is structured as follows: We start with a consideration and discussion of the boundary conditions in the general case of arbitrary values of specific surface energies of all three interfaces. We then recapitulate the results of the simplest "membrane model" using a Fast Fourier transform (FFT)-based boundary element method (BEM) implementation. Finally, we generalize the FFT-BEM to the case of non-adhesive contact with different surface energies inside and outside the contact and illustrate this case with a simulation example.

BOUNDARY CONDITIONS TAKING INTO ACCOUNT SURFACE TENSION
Analysis of a contact between a rigid body and an elastic half space taking into account surface tension, carried out in Karpitschka et al. (2016) and Popov (2020), shows that the Young's law determining the contact angle of liquids is also valid for solids. This leads to the following boundary conditions: p(x, y) = γ 1 u(x, y), outside the contact area (5) γ 1 cos θ = γ 2 − γ 12 , at the contact boundary (6) where f (x, y) is the profile of rigid indenter and u(x, y) the normal displacement of the surface of elastic half space, d is indentation depth, θ is contact angle, is Laplace operator, and p(x, y) is the normal pressure at the surface immediately below the (infinitely thin) surface layer. Note that for definitions of the surface profile we used other direction of the z-axis (out of the elastic halfspace) than for displacement and other quantities (into the halfspace). It is also important to note that Equation (5) is written in approximation of small slopes of the surface of elastic body, which, however, is the necessary condition also for application of the superposition principle used throughout the paper. Equation (4) states that the surfaces of the indenter and the elastic halfspace coincide in the contact area. Equation (5) determines the elastic stress at the surface under the tensioned surface layer. Equation (6) is the equilibrium condition of the contact boundary under the action of three surface tensions. Similarly to the contact angle in liquid contacts, θ is a thermodynamic property of the system. It does not depend on the shape of the body and deformation of the surface.
Using Equation (1) for the work of adhesion, w, Equation (6) can be rewritten as This equation is known as Young-Dupre equation (Young, 1805;Dupré, 1869). For non-adhesive contacts with surface tension it has w = 0. From Equation (7), it follows that cos θ = −1 and This equation means that at the boundary the slope of the surface profile is equal to the slope of the elastic half-space. Thus, in this specific case the condition (6) can be reformulated as a continuity of the surface slope at the contact boundary ( Figure 1B). In this article, we focus on this non-adhesive contact with surface tension and its numerical simulation.

NUMERICAL SIMULATION OF NON-ADHESIVE CONTACT WITH SURFACE TENSION
A numerical solution to non-adhesive contact with surface tension for the case IIa (γ 1 = γ 12 : = γ ) is based on Hajji's equation for the vertical surface displacement u at the position x, y caused by the pressure distribution p c x ′ , y ′ : with where H 0 and Y 0 are Struve and Bessel functions of the second kind of zero order (Hajji, 1978). It should be stressed that the pressure p c (x, y) is the pressure acting on the surface of the upper tensed layer and is not equal to the pressure p(x, y) in Equation (5), which is the pressure under the tensed layer. Qualitatively speaking, p c (x, y) is the "external pressure" acting on the elastic body coated with membrane, while p(x, y) is the pressure in the elastic body immediately under the membrane. For vanishing surface tension the first term approximations to (9) and (10) become the classic Boussinesq's equation In this case, p = p c . As already stated, according to the theory of surface elasticity (Gurtin et al., 1998), the contact with surface tension can be considered as a contact with an elastic body covered by an (infinitely) thin membrane. The equilibrium and constitutive equations in the bulk of the elastic body still follow the theory of elasticity. Therefore, we can find a solution to contact with surface tension by use of the fundamental solution for elastic halfspace without surface tension. Using the definitions of pressure p(x, y) and p c (x, y), as shown in Figure 1B, the equilibrium condition of the interface in contact the contact region is p c (x, y) = p(x, y) − γ 12 u x, y , inside the contact. (12) The other two boundary conditions (4) and (5) should of course also be fulfilled. The deformation u x, y is generated by pressure p(x, y), which can be obtained from the Boussinesq's equation.

Realization in the Boundary Element Method
The BEM was developed in the past 30 years and represents an effective numerical tool for solving various contact problems including homogeneous material and layered systems, adhesive and non-adhesive contact, rough contact, and thermal contact (Nogi and Kato, 1997;Liu et al., 2000;Campañá and Müser, 2006). Some BEM formulations are based on integral equations of the form of (9) or (11). In a discretized form they can be written as where u ij is the displacement of surface element at position i, j in two-dimensional discretization, p i ′ j ′ c is stress acting on an element located on i ′ , j ′ , andK iji ′ j ′ is the influence coefficient and its value is calculated analytically or numerically depending on the form of the fundamental solution. The matrix of influence coefficients can be found, e.g., in Pohrt and Li (2014). The convolution nature of integral equations of the form (9) allows formulating and solving a contact problem using the FFT. Using this technique increases computing efficiency by replacing the integration by multiplication and therefore reducing the complexity from N 4 to N 2 log N 2 (for a system with discretization of N × N elements). For example, the algorithm of evaluating Equation (13) can be represented as where p c is the pressure vector matrix, u is the matrix of discretized displacements, and K is the matrix of the influence coefficient. Note that for the non-periodic contact in a finite domain, the techniques of zero padding and wraparound order of matrix of pressure and influence coefficient in a doubled domain should be used to execute cyclic convolution. Thus, the p c and K should be expanded to dimension 2N × 2N when FFT is performed. The displacement u is then extracted from the obtained displacement with the same dimension 2N × 2N after the IFFT. The details of theories of the linear and cyclic convolutions and their applications to periodic and non-periodic contacts as well as numerical procedures can be found in Liu et al. (2000) and Ju and Farris (1996). In Equation (14) the fundamental solution (given by the matrix K) should be first Fourier transformed. Alternatively, one can directly use the fundamental equation in Fourier space, which in most cases can be calculated very easily. For example, for an elastic half-space it is equal to where q = q 2 x + q 2 y . This gives the relation of the Fourier components of pressure and displacement u q x , q y = 2 E * p q x , q y q .
Displacement in the coordinate space in then calculated as The general form of the fundamental equation in Fourier space can be written as where the Fourier transform of the kernel of convolution, C q , depends on a particular system. The only general prerequisite for this simple form is the lateral homogeneity of the system. In the direction perpendicular to the surface, it can be heterogeneous (e.g., a layered system or functionally graded material). C q is of course also a function of the effective elastic modulus in the case of a homogeneous material and coating system. Equations (12) and (5) for non-adhesive contact with surface tension can be written in Fourier space as follows: p c (q) = p(q) + γ 12 · q 2 · u q , in contact, p(q) + γ 1 · q 2 · u q = 0, out of contact.
We first consider the simple case with equal surface energy inside and outside the contact area, γ 1 = γ 12 = γ . Substitution of (18) into (19) gives the relation between pressure p c and deformation u u q x , q y = p c (q x , q y ) For the displacement field in coordinate space, it follows Thus, any available BEM program can be used for simulating this system covered with a tensed membrane just by adding the term "γ q 2 " in the Green's function. Note that this equation connects the surface displacement with the external pressure. Therefore, the boundary condition and numerical procedure can be dealt in the same way as in all existing BEM: the pressure outside the contact is zero: p c x, y = 0, outside of the contact, while the surface displacement inside the contact area meets the geometry condition. Of course, the zero padding and wraparound order are still needed for the non-periodic contacts as discussed in Liu et al. (2000). For an isotropic homogeneous material, Equation (22) becomes

Case Study 1: Isotropic Homogeneous Elastic Half-Space Coated With a Tensed Membrane
Let us discuss numerical simulation of a "Hertzian" contact with surface tension: a contact of a rigid parabolic indenter with elastic half-space taking into account surface tension, which we first consider equal inside and outside the contact. The result obtained using the Hajji's fundamental solution was given in Long and Wang (2013). Here we would like to compare this solution with the above-described solution based on the fundamental solution (1/C(q) + γ · q 2 ) −1 in the Fourier space. Parameters used in the simulation are as follows: E * = 0.011 MPa; sphere radius R = 1 mm; indentation depth d = 0.01 mm; and two different specific surface energies, γ = 0.2 N/m and γ = 0.4 N/m, were considered, corresponding to the values γ /(E * a 0 ) = 0.18 and 0.36 of the parameter γ /(E * a 0 ) (ratio of the elastocapillary length and the Hertzian contact radius a 0 ). Figure 2 shows numerical solution obtained with discretization 512 × 512 points, where the triangles in Figure 2A are pressure distributions calculated based on the Hajji's fundamental equation, and stars are results based on the alternative approach using the fundamental solution (1/C(q) + γ · q 2 ) −1 . The calculating time by using this fundamental solution is 12 times smaller than that with the Hajji's fundamental solution. The contact behavior found in Long and Wang (2013) practically coincides with that found with the above-described procedure. In particular, the pressure at the contact boundary has a "jump" in the case with surface tension (Figure 2A), and the surface tension leads to a reduction of contact area (Figures 2A,B).

Case Study 2: Layered System
As explained above, any existing BEM formulation for a contact without surface tension can be trivially extended to include surface tension by a small change in the corresponding fundamental solution in Fourier space, according to Equation (21) or (24). Let us illustrate this with an example of an elastic layer with surface tension. We proceed from the BEM formulation without surface tension described in Li et al. (2019). Adding the term "γ q 2 " in the fundamental solution, we include the surface tension. Here we show a case of soft layer bounded on the elastic half space. The layer had elastic modulus E * 1 = 0.011 MPa, specific surface energy γ = 0.2 N/m, and thickness h 0 = 0.03 mm. The elastic modulus of the underlying halfspace is E * 2 = 10E * 1 . The profile of indenter is a sphere with radius R = 1 mm superposed with a two-dimensional waviness with wavelength λ = 0.03 mm and a small amplitude of h = 0.0005 mm. The indentation depth was d = 0.01 mm. The simulation results without and with surface tension are shown in Figure 3. The color map presents the pressure distribution in the contact region. Two cross sections of this map are selected to show the details. Under the above conditions, the contact area in the system without surface tension is compact (left figure), while "switching on" the surface tension makes it much more heterogeneous and even not simply connected (right figure). This leads to much more intensive oscillations of pressure. Note that the pressure shown in Figure 3 is the pressure immediately under the surface of the indenter. This pressure is relevant for estimating the possible damage of the indenter.

Different Surface Energies Inside and Outside of the Contact Area
When surface energies inside and outside the contact are different, γ 1 = γ 12 , the FFT of Equations (19) and (20) cannot be directly carried out for the whole region. However, we can rewrite Equation (12) in the following form: p 1 (x, y) = p(x, y) − γ 1 u x, y , in the whole area, (25) p(x, y) = p 1 (x, y) + (γ 1 − γ 12 ) u x, y , in the contact area, where p 1 (x, y) is an auxiliary function. In Fourier space, Equation (25) takes the form p 1 q = p q + γ 1 · q 2 · u q = 1/C q + γ · q 2 · u q .
The computation algorithm is the following. Initially we assume that the contact area is the geometrical intersection A and the displacement of the surface of the elastic half-space in this area coincides with the profile of indenter u A (x, y); the corresponding FFT (u A ) is calculated. Using this FFT (u A ) one can obtain the auxiliary stress function p 1 q as well as p 1 (x, y) appearing in (25). After that, the deformation u x, y can be calculated within the contact area. With p 1 (x, y) and u x, y as well as the curvature of the surface, u x, y , one gets the pressure distribution p(x, y) from Equation (26). After that, the usual correction and iteration procedure starts: the elements that have negative pressure or geometrical penetration out of contact are marked as "detached" so that a new contact area, A ′ , is generated. With this new contact area, the above iteration is repeated until both pressure and geometry conditions for all elements in contact area are met. Example 1 We give an example with the same parameters as studied in the case of Figure 2, with the only difference that the surface energy outside contact area, γ 1 , is not equal to that inside, γ 12 . Surface energy in the noncontact region was equal to γ 1 = 0.2 N/m. Two cases were studied: one with smaller surface energy inside contact γ 12 = 0.1 N/m and the other one with larger γ 12 = 0.4 N/m. The obtained pressure distribution and surface displacement are shown in Figure 4. The results with equal surface energy are also shown with stars for a comparison. One can see that when the surface tension inside the contact is larger, the loading needed for generating the same indentation depth is also higher due to the surface tension effect (curve with triangles). Interestingly the contact area and the surface displacement are the same in all three cases. It is clear that the surface displacement in the contact coincides with the profile of the indenter shifted by the indentation depth. Displacement outside the contact is determined by the pressure p and surface tension γ 1 , where p is numerically calculated by Boussinesq's equation with the given surface deformation inside the contact.
So surface displacement will be the same for unequal surface energy inside the contact if the contact area is unchanged. The difference of surface energy changes only the pressure p c with a constant value in the case of sphere contact because of the constant mean curvature. Figure 4 shows that the shape of the surface of elastic body does not depend on the surface tension inside the contact area. A posteriori, this seems to be a trivial conclusion. However, it is trivial only under the assumption of non-adhesive contact, which guarantees the unchanged condition at the boundary of the contact area, which does not depend on the interface surface energy density. The only quantity that is influenced by the   specific surface energy of the contact interface is the pressure distribution. This means that the computation procedure in the case of two different surface energies can be simplified as shown in Figure 5. 1. In the first step, the contact problem is solved with a constant specific surface energy equal to the specific surface energy outside the contact area. 2. In the second step, the pressure distribution inside the contact area is corrected by the term δp(x, y) = (γ 1 − γ 12 ) u x, y .
Example 2 Consider a wavy surface similar to that shown in Figure 3 with the radius of curvature, R = 1 mm; wavelength, λ = 0.02 mm; and amplitude of waviness, h = 0.0001 mm. The elastic modulus of the half space is assumed E * 1 = 0.013 MPa and the indentation depth d = 0.012 mm. In the following two cases are considered: (1) equal specific surface energies γ 1 = γ 12 = 0.1 N/m and (2) two different specific surface energies, γ 1 = 0.1 N/m and γ 12 = 0.2 N/m, inside and outside the contact correspondingly. Numerical results are shown in Figure 6. The contact area and surface displacement are the same in both cases.

CONCLUSION
Non-adhesive contact with surface tension is usually solved using the fundamental solution of Hajji. This fundamental solution has a complicated form in the coordinate space but can be derived in an extremely simple way in the Fourier space. More than that, this derivation has a universal character and can be applied for more complicated situations, as, e.g., for layered systems or functionally graded materials. With an example of Hertzian contact this method was validated by comparison with the fundamental solution in the coordinate space, and it was then applied to the rough contact of a coated system.
We argue that the approach based on the fundamental solution can be used only in the case when the surface tension inside and outside of the contact area are equal. There are no physical reasons why it should be the case. Therefore, in the general case, the approach based on the use of the fundamental solution and superposition principle does not work. We show how this problem can be reduced to the simpler contact problem with constant surface tension inside and outside the contact area.
Of course, the case of non-adhesive contact with different surface energies inside and outside the contact is a little bit an academic exercise. In a general "dry" contact of two solids, the condition of non-adhesive contact normally will not be fulfilled. However, the contact can be made non-adhesive by introducing an intermediate liquid with dielectric constant equal to that of one of the bodies. According to the theory of van der Waals' interactions by Dzyaloshinskii et al. (1961), this leads to suppression of van der Waals force (and thus of the separation energy), while the surface tensions of both bodies remain generally non-zero.