Many physiological systems can be modeled as the layer-foundation structure, and generally, the mechanical properties of surface layer differ from those of bulk interior. Early in 1978, Hajji (1978) investigated the indentation of lung under uniform pressure by considering the sample as an isotropic elastic half space ideally adhered with a tensed membrane. It is assumed that the surface membrane thickness is ultimately small so that its bending rigidity can be neglected. Besides, for small deformation, the membrane tension is assumed keeping constant. With the same assumptions, Kim and Gouldstone (2008) addressed the axisymmetric spherical indentation of elastic solid with strain-independent membrane tension. Moreover, the surface layer of some soft tissues and single cells can be regarded as pre-tensed “plate” or “shell” with finite thickness, which can also resist bending deformation (Zamir and Taber, 2004; Zhang and Zhang, 2009). Recently, Argatov and Sabina (2012) treated the surface layer as a reinforced membrane under generalized plane stress state. They hypothesized that no pre-tension exists in the membrane, and the flexural stiffness is negligible compared with the in-plane tensile stiffness. Such modeling method was employed in the deformation analysis of anisotropic articular cartilage (Argatov and Mishuris, 2016).

It is interesting to note that the above models originally developed for the biological materials are analogous to the well-known surface elasticity theory of Gurtin and Murdoch (1975), which characterizes the material surface with surface tension and surface elasticity. As Hajji stated (Hajji, 1978), the mathematical expressions are basically consistent, although they have completely different physical natures. The constant membrane tension in Hajji (1978) and Kim and Gouldstone (2008) is corresponding to the residual surface tension, and the superficial layer modeled in Argatov and Sabina (2012) and Argatov and Mishuris (2016) can be considered as a solid surface with surface elasticity. Recently, problems of surface loadings and contacts have been widely studied based on the surface elasticity theory (He and Lim, 2006; Wang and Feng, 2007; Zhou and Gao, 2013; Long et al., 2017; Li et al., 2020). It is found that the deformation behavior would be distinctly different from classical models when the size of loading is comparable or even smaller than a critical length, usually at nano/microscale. However, these analytical works are usually limited to simple cases under a given surface pressure or for the symmetric situations.

Physically, the surface elasticity theory is based on the concept of solid surface energy (Gurtin and Murdoch, 1975). For the general contact between two solid bodies, the surface energy of each free surface (γ₁, γ₂) and contacting interface (γ₁₂) should be taken into account simultaneously. Two simplified cases can be recognized according to values of these surface energies. When the energy difference w = γ₁+γ₂-γ₁₂ (called work of adhesion) is appreciable and the effect of surface energy outside the contact can be neglected, adhesive contact model with specific work of adhesion would be appropriate (Johnson et al., 1971). On the other hand, if the energy difference equals to zero (w = 0) but the surface energy of one contact body is much larger than the other (i.e., γ₁ >> γ₂ and γ₁ ≈ γ₁₂), the contact problem can be addressed in the framework of surface elasticity theory (He and Lim, 2006; Wang and Feng, 2007; Zhou and Gao, 2013; Long et al., 2017; Li et al., 2020). In this work, we consider the latter circumstance with constant surface tension or, equivalently, the non-adhesive contact of solids with tensed surface membrane.

The Fast Fourier Transform (FFT) based methods have been widely employed to solve the elastic contact problems for its great advantage on reducing the computation cost. For the efficient Green's function molecular dynamics (GFMD) (Campañá and Müser, 2006; Prodanov et al., 2014), the FFT technique plays an important role in the determination of substrate elastic response to external surface forces. In addition, the boundary element method (BEM) that exploits FFT to accelerate the calculation of displacements induced by given surface forces has been also proved fast and robust for the analysis of various contact problems (Nogi and Kato, 1997; Liu et al., 2000; Polonsky and Keer, 2000; Pohrt and Popov, 2012, 2015; Pohrt and Li, 2014; Rey et al., 2017; Bugnicourt et al., 2018). The basic principle of the FFT-based BEM is to evaluate the linear convolution of fundamental solution and pressure distribution based on the convolution theorem. For the classical elastic problem with half space approximation, the simple Boussinesq's fundamental solution or its form in Fourier space is commonly used in the BEM for normal contact problems. However, when the surface of elastic half space is equi-biaxially tensed, things will be different.

This paper aims to extend the mature FFT-based BEM for the calculation of normal contact of material with membrane/surface tension. Several calculation examples are implemented to demonstrate the validity and usefulness of this method. This method provides a potential numerical approach for the study of contact behavior of some soft structures, which is of particular interest for the measurement of mechanical properties of biological tissues (Zhang et al., 2014).

We consider an ideally smooth elastic substrate. The surface is tensed by a constant, strain-independent tension τ₀, and the bulk material is homogeneous having elastic modulus E and Poisson's ratio ν. The origin of the Cartesian coordinate system (o-xyz) is located at the surface of half space, and the z-axis is perpendicular to the surface, as shown in Figure 1A.

With some specified external pressure applied on the surface, the boundary value problem can be solved in the framework of classical linear elastic theory, provided that the unconventional boundary condition associated with surface tension is employed (Hajji, 1978; He and Lim, 2006; Wang and Feng, 2007; Zhou and Gao, 2013). For the axisymmetric situation under a uniform normal pressure p₀ within a circular region of radius a, the boundary condition reads,

where u(r) is the vertical displacement on the surface, σ_zz and σ_zr are the bulk normal stress and shear stress at z = 0, respectively, and H is the Heaviside step function.

By using Hankel integral transform and substituting the boundary condition into the general solutions, the vertical displacement caused by the uniform pressure p₀ can be derived (Hajji, 1978). Moreover, by implementing a limiting process in which the radius a decreases to zero with the resultant force held constant unity πa²p₀ = 1, the fundamental solution that is the vertical displacement at the point (x, y) induced by unit concentrate force acting at the coordinate origin can be derived as,

where s = 2τ₀/E^* is the characteristic length, E^* = E/(1–ν²) is the reduced modulus, H_n and Y_n are the Struve function and the Bessel function of second kind of order n, respectively.

The FFT of this fundamental solution, G~(ω), is plotted in Figure 2, where ω=ωx2+ωy2 is the magnitude of wave vector (ω_x is the angular frequency in the x-direction and ω_y is the angular frequency in the y-direction). It is found that the results are consistent with the theoretical prediction (Li and Popov, 2020),

which reduces to the Boussinesq's solution in Fourier space when the characteristic length s decreases to zero, G~(ω)=2/(E*ω ).

For the indentation of a rigid body with known profile f (x, y) into a depth δ, see Figure 1A, the vertical displacement on the surface u(x, y) has to fulfill the following condition:

where Ω is the domain of contact.

Assume that the positive contact pressure between the substrate and the indenter distributes as p(x, y). Then, according to the superposition principle, the vertical displacement generated by p(x, y) can be expressed by

which is exactly a two dimensional continuous convolution of fundamental solution and contact pressure distribution. As a result, substitution of Equation (5) into the displacement boundary condition Equation (4) leads to

Note that the inequality characterizing the non-overlapping condition outside the contact region is always accompanied with Equation (6) no matter which expression form is utilized. The resultant load P acting on the rigid indenter equals to the summation of the contact pressure over entire contact area. For a given indentation depth, the corresponding contact area and contact pressure are to be determined.

Next, this problem will be solved in the light of FFT based boundary element method. In contrast to finite element method, only part of substrate surface needs to be discretized. As shown in Figure 1B, squared elements with side length Δ are used to mesh the surface. Assuming that the pressure within each single element distributes uniformly, Equation (6) can be transformed into linear algebraic equations,

where

Note that the computational domain that is discretized into surface elements should include all potential contact region. If the computational domain contains N×N elements, the influence coefficients K_ijkl have N⁴ values. Using direct multiplication method has the complexity of O(N⁴). Instead, we prefer to apply the FFT technique to evaluate the vertical displacement for a given pressure, i.e., the left side of Equation (7). The complexity is then reduced to O(N²logN). It is worth mentioning that the multilevel multi-integration method suggested by Brandt and Lubrecht (1990) allows for the same computation acceleration.

The basic idea of the FFT based speed-up method is to interpret the left side of Equation (7) as a two-dimensional discrete convolution and applying the circular convolution theorem (Liu et al., 2000; Pohrt and Li, 2014). In this work, the influence coefficient matrix (K_kl)_N×N is first prepared in real space by performing numerical integration of Equation (8) with x_i = y_j = 0. It should be pointed out that appropriate zero padding and wrap-around order operations to the original influence coefficient matrix (K_kl)_N×N and pressure matrix (p_kl)_N×N are necessary for the proper application of circular convolution theorem, which can avoid the periodic error (Liu et al., 2000). Next, FFT is implemented on the two matrices. The results can be easily obtained through element-wise multiplication of the influence coefficient matrix and pressure matrix in Fourier space. Then, performing an inverse FFT (IFFT) yields the vertical displacement caused by the pressure distribution (p_kl)_N×N in real space. As a result, the problem becomes,

where ∘ represents the element-wise product, and h is given by the indenter shape and the indentation depth, h_ij = δ +f (x_i, y_j).

With the FFT based acceleration technique, this inverse problem can be solved using an iteration scheme based on conjugate gradient method (Liu et al., 2000; Polonsky and Keer, 2000; Pohrt and Li, 2014). The final obtained contact pressure (p_kl)_N×N should be positive at all contact elements where the induced vertical displacement matches with the given rigid profile, and equal to zero in the non-contact region where the non-overlapping condition is always satisfied. Finally, the overall contact responses under a specified indentation loading δ are ready to be evaluated. The resultant force P applied on the indenter is obtained by summing up the forces at all discrete contacting elements. The real contact area A is computed by multiplying the number of contacting elements with the area of a single element Δ². It should be pointed out that the contact area obtained by this direct summation method can be only accurate when the contacting surface is discretized by sufficiently refined grids. Area correction to continuum limit as done by Prodanov et al. (2014) is necessary. Inevitably, the BEM calculation with finer surface discretization would require more time cost. In addition, it is worth mentioning that Yastrebov et al. (2017) presented a correction route to compute accurate contact area with relatively coarse discretization for such kind of numerical method.

In summary, the FFT-based boundary element method is extended for normal contact of soft material, of which the surface is tensed by equi-biaxial tension. Three type of indenters compressing on the elastic half space with constant membrane/surface tension are considered. The calculation results show excellent agreement with related available solutions. Typical size-dependent contact behavior is demonstrated. In particular, due to the presence of membrane/surface tension, the eccentricity of contact ellipse will be slightly increased for the contact with a rigid elliptical indenter. For a given substrate indented by a specific rigid indenter, the mean contact pressure defined by the ratio of indentation load to the real contact area can be increased by increasing membrane/surface tension on the substrate surface. This numerical method would be helpful for the deformation analysis of some soft systems with membrane/surface tension.

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

WY and GW conceived the study. WY carried out numerical simulations. All authors drafted and reviewed the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.