Analysis of tangential contact boundary value problems using potential functions

This paper presents an analysis technique of high-order contact potential problems and its application to an elastic settlement analysis of a shallow foundation system subjected to a combined traction boundary condition. Closed-form solutions of potential functions are derived for an elastic half-space subjected to bilinear tangential traction boundary conditions over rectangular surface regions. Using the principle of superposition, the present solutions provide a means to form an approximate and continuous solution of elastic contact problems with higher-order tangential boundary conditions. As an application example, an elastic settlement analysis of a rigid footing founded on a dense granular soil is performed under a tangential traction boundary condition prescribed in an analogy with the stress equilibrium states of static sandpiles. A generalized solution approach to combined normal and tangential traction boundary value problems is discussed in the context of foundation engineering.


Introduction
The problem of the equilibrium state of an isotropic elastic halfspace with given conditions on its planar boundary has been known historically as 'the problem of Boussinesq and Cerruti' [1], after the mathematicians who originally developed solutions to the problem in the late nineteenth century. The solutions to the governing equilibrium equations were given by these authors in terms of potential functions that correspond to specified boundary conditions. The problem where normal tractions are prescribed is often specifically associated with the name of Boussinesq. The tangential traction boundary value problem is then referred to as the Cerruti problem [2,3]. Since these scholars presented both solutions nearly simultaneously [4,5], the distinction appears to be mostly one of convenience. contact pressures develop at the interface between the structural footing and a supporting soil body against the external loads transmitted to the soil through the structure. Normal contact pressures are typically prescribed as the traction boundary condition for elastic settlement analysis of shallow foundations subjected to vertically applied concentric loads. The resultants of the contact stresses and the actual distribution in a specific problem are not known a priori, and are the final product of the soil-structure interaction, which is statically indeterminate. Under the requirements of static equilibrium, only the average value of the contact pressure is statically determined simply by dividing the averaged resultant load by the area of contact of the (relatively) rigid footing with the soil. In turn, the spatial distribution of the normal traction boundary condition is often simplified to be uniform, and it is assumed that any shear stresses vanish at the contact interface between the footing and a supporting soil medium [6,7]. However, numerous empirical measurements of stresses beneath piles of sand [8][9][10][11] and laboratory-scale foundations [12][13][14][15] reveal that both normal and shear stresses coexist in a high-order spatial variation that is believed to be dependent on construction and loading histories of the supporting soil and the magnitudes of vertically applied concentric loads. For these reasons, the authors have addressed the problem of calculating the response from highly variable empirical normal traction fields in a past manuscript [16]. Towards the second assumption, it is reasonable that some frictional forces develop in the contact plane where internal resistances of the soil medium may cause confining effects against horizontal expansion when compressed. Thus, it is hypothesized that the tangential traction components may develop even when the footing is subjected to concentric vertical loads. Although there have been few empirical measurements for the shearing stresses beneath vertically loaded shallow foundations (e.g. [13]), a striking similarity between the normal contact pressure distributions of shallow foundations and sandpiles further elaborates the conceptualization of a combined normal and tangential traction boundary condition for shallow foundations. This rationale serves as a motivation for the present study.
The form of the elastic solutions for tangential traction boundary value problems in terms of potential functions is well understood [1]. Yet, closed-form expressions for the potential functions have been calculated only for tractions prescribed by low-order polynomial laws [17,18]. When a traction boundary condition is high order, the calculation of closed-form solutions becomes intractable. It is practical to consider an approximate solution by superimposing known solutions of lower order as per discretization of higher order traction boundary conditions. The conceptual framework of the approximate solution approach [19,20] has shortcomings such that discontinuities in traction fields along the local (discretized) boundaries may lead to singularities in the calculation of stresses when tractions are prescribed over rectangular subdomains in constant and/or unidimensional linear terms only [6]. Recently, progress has been made to overcome these computational difficulties [16,21 -24] in obtaining converged solutions for the elastic contact problems.
In an earlier work [16], the authors recommended an approximate yet continuous solution to the problem of high-order normal traction boundary conditions. The present study aims to provide corresponding solutions to tangential traction boundary value problems as a companion to [16]. A complete set of closed-form solutions of potential functions for bilinear tractions and their derivatives is derived through methods of substitution and direct calculation of area integrals. By the interpolation and superposition of these solutions, the displacement and stress fields of the half-space are approximated under combined normal and tangential traction boundary conditions to a desired degree of accuracy. The method is verified through a study of convergence to the classical Hertz-Mindlin contact solution [25,26]. To further the application to foundation design and analysis, a hypothesis is presented to qualitatively prescribe a plausible tangential contact traction field between a rigid square footing and a supporting granular soil. Parameters of empirical expressions are determined in analogy to the stress-equilibrium states of sandpiles [8 -11], in particular with respect to the continuum-based fixed principle axis (FPA) model [27,28]. The present method of solution is applied to solving a combined boundary value problem for normal and tangential tractions at the base of a square shallow foundation in a Cartesian coordinate system.

Description
We define a right-handed set of Cartesian coordinates such that the positive z-axis points downwards into the body, as in figure 1. Let the displacement at point (x, y, z) be denoted by (u, v, w). The full details of the boundary value problem are available in many texts on the mathematical theory of elasticity [1,29]. In the Neumann problem, surface traction, given by Cartesian components (q x , q y , p), is prescribed on the boundary plane. These expressions for the Cartesian traction boundary conditions are related to the inplane surface stresses and the partial derivatives of the surface displacements by Hooke's Law þ @w @y z¼0 and p ¼ Às zz j z¼0 ¼ À(l þ 2m) @w @z z¼0 À l @u @x z¼0 þ @v @y z¼0 , where l and m are Lamé's constants of an isotropic elastic material. It is useful to specify some subset of the boundary R, within which some given boundary values are finite and outside of which they are zero. This stems from the multitude of applications of this theory to contact problems, within which surface tractions and/or displacements are employed to model a region of contact between two bodies. In general, R can be any subset of the boundary plane. It may even represent a union of connected or disjoint subregions. We take (x 0 , y 0 ) to be the values of the in-plane Cartesian coordinates within the region R. The new closed-form solutions found in this work (see appendix A) will deal specifically with general rectangular regions on the surface of the half-space, which is of the form R ¼ f(x 0 , y 0 )ja 2 x 0 a 1 , b 2 y 0 b 1 g, as depicted in figure 1. The principle of superposition simplifies a solution procedure of the combined traction boundary conditions of equation (2.1). A set of solutions specific to each of the three traction fields can be obtained separately. Each of the three sets of solutions can then be summed as a solution to the combined traction boundary value problem. The solution for non-zero normal traction is detailed in [16]. Here we provide the details of the solution of the case where q x = 0 and q y ¼ p ¼ 0.
The equilibrium solutions of a homogeneous, isotropic, linear-elastic half-space when the traction q x (x 0 , y 0 ) acts in the x-direction over the region R of the planar boundary can be given in terms of three potential functions 2Þ Figure 1. An elastic half-space bounded at z ¼ 0 described by a Cartesian coordinate system. This figure is retained from the authors' previous manuscript [16].

Solution by superposition and example calculations
Solutions for high-order contact traction boundary conditions may be obtained from the superposition of discretized solutions corresponding to lower-order interpolation functions. A tangential traction field q x (x 0 , y 0 ) of arbitrary complexity can be approximated over a discretized loaded domain using the piecewise interpolation of bilinear polynomials. Subsequently, the lower-order solutions for each discretized subdomain are superimposed to approximate the displacement and stress fields in the elastic body. The discretization scheme is detailed in the authors' previous work [16], and thus, is omitted in the present paper. The procedure of the present solution is analogous to the collocation techniques used in the indirect boundary element method [30,31]. However, the present method incorporates the analytic calculation of the boundary integrals and superimposes the resulting functions over the whole domain, replacing a point-wise numerical quadrature scheme which would be computationally more expensive. The proposed solution method is analogous to the analytic element method, originally proposed by O.D.L. Strack to study groundwater flow [32,33]. Under this interpretation, the closed-form solutions presented here can be viewed as new 'analytic elements' derived for the solution of Neumann boundary value problems in an elastic half-space.
In the present work, the solution for an arbitrary distribution of tangential traction is approximated as a linear combination of potentials for constant, linear and bilinear tractions. If we divide the loaded area R into M and N uniform intervals in the directions of x 0 and y 0 , respectively, then the discretization creates M Â N disjoint rectangular subdomains: We define expressions for the potentials under the simplest polynomial traction laws over subdomains R ij as such that A 00 ij is the potential of equation (2.4) when q x ¼ 1, G 01 ij is that of equation (2.2) when q x ¼ y 0 , etc. The present authors have calculated closed-form solutions of these functions and all derivatives present in equations (2.5)-(2.7) for constant, linear and bilinear loads, i.e. for cases m, n ¼ 0, 1. These are listed in appendix A.
The potential defined by equation (2.2), which Boussinesq called the 'second logarithmic potential function', does not appear in the solution of the problem when only normal traction components have non-zero boundary values. The closed-form solutions for the corresponding equation (2.10) were therefore omitted from [16]. However, the second logarithmic potential function is one of the major difficulties in obtaining a closed-form solution to the tangential traction problem. We therefore provide some example calculations which focus on this expression.
The function @ 3 c x /@x 3 required for the solution of the stress component s xx is approximated as where c mn ij are interpolation constants corresponding to the order of the polynomial components of the boundary tractions (x 0 ) m ( y 0 ) n . For example, the last term in the bracket on the right-hand side is defined as Owing to the linearity of the integral, equation (2.11) can be rewritten as The open-form integro-differential expression on the right-hand side can be found to be xi y jþ1 y j through direct calculation by means of substitutions and reversing the orders of integrations and differentiations. This way, the function of equation (2.11) can be written in its precise closed form and evaluated over the half-space using any mathematical software package. The remaining terms on the right-hand side of equation (2.11) can be calculated in a similar manner (see appendix A). As a second example, we consider another potential that contains partial derivatives taken with respect to both x and y: The closed-form solution of equation (2.12), which is required for the calculation of s yy , can be obtained in a straightforward manner. The two differentiations with respect to x and y, cancel out the integrations with respect to x 0 and y 0 , as per the calculation procedure of the first example: @ 2 G 00 @x@y À y @ 2 G 00 @y 2 À 2 @G 00 @y : Using the property of the harmonic function @G ij 11/@ x (i.e. it satisfies the Laplace equation D(@G ij 11/@ x) ¼ (@ 2 /@x 2 þ @ 2 /@y 2 þ @ 2 /@z 2 )(@G ij 11/@ x) ¼ 0 everywhere in the half-space), we can deduce: closed-form solutions for the hyperbolic-paraboloidal (bilinear) tangential tractions that are listed in appendix A. Note that the variables associated with the x-and y-directions are interchangeable. Thus, the solution for the case where q y = 0 and q x ¼ p ¼ 0 can be easily obtained. For the solution of the problem corresponding to the normal loading case, that is where p = 0 and q x ¼ q y ¼ 0, along with the complete details of the computational procedure applied here, refer to [16].

Solution to a Cattaneo -Mindlin boundary value problem
For the validation of the solution procedure, we model two elastic spherical bodies in contact subjected to relative tangential displacements within a Hertzian normal contact interface, and benchmark on an analytical solution given by Cattaneo [34] and Mindlin [26]. A circular region of contact of radius a is formed between two elastically similar spheres under normal force P, and a subsequent tangential force Q is applied in the x-direction relative to the plane of contact. No relative displacement occurs in a 'stick' zone of radius a stick , a, while slip occurs in the outer annulus. The resulting tangential traction q x is assumed to be governed by the Coulomb model of static friction. Given the normal traction p, q x c f p within the extent of the contact region, where c f is the coefficient of friction. The tangential contact traction q x can be written as for r a stick , is a radial polar coordinate within the circular contact region. Mindlin [26] provides an analytical solution of uniform tangential displacement within the stick region r a stick : where n ¼ l/2(l þ m) is Poisson's ratio of the material. Thus, applying the tractions (3.1) in the tangential x-direction as a Neumann boundary condition to a homogeneous, isotropic, elastic half-space should yield the constant x-directional displacements of the stick region given by (3.2). The elastic parameters used in the verification of convergence are listed in table 1. The convergence of the approximation of equation (3.1) by piecewise-bilinear surfaces is estimated in terms of the corresponding total applied tangential force. The boundary of the Hertzian circular contact region is approximated by discretization with subdomains dividing the circumscribed square, as schematically illustrated in figure 2a. The error in the approximation of the tangential traction field is calculated as ð3:3Þ Figure 3 shows the convergence of the displacement ratio within the stick zone using an error estimate Both stress and displacement fields simultaneously satisfy equilibrium and compatibility. The approximation improves as a finer discretization is used to converge toward the exact circular boundary of the loaded domain using the rectilinear subdomains. However, the present discretization approximation scheme would yield a faster convergence rate when the geometry of loaded domain is rectilinear.

Application to settlement analysis of shallow foundations
Having verified the solution method, we undertake an application example of foundation engineering. When contact forces develop in the plane of non-conformal contact surfaces between a rigid footing and a supporting granular medium, a highly variable traction field forms at the contact plane. The corresponding equilibrium state of the supporting granular medium, and ultimately the foundation settlement, are dependent upon the boundary conditions that span multiple scales. These include, but are not limited to, variables at both the soil-structure system scale, for example, geometry and loading history [10], and the individual grain constituent scale, for example, the rheology of a roughsurface contact. In principle, the effects of intergranular forces can be quantified through statistical analysis using discrete element contact models [35 -37]. However, the computational analysis required to conduct a full parametric study is daunting and costly. An alternative for practical engineering analysis is to develop a mathematical model in which these multitudes of particulate interaction phenomena are deliberately reduced to a single quasi-static boundary value problem referring only to phenomenological observations at the foundation (system) scale. Based on a review of experimental [8][9][10][11][12][13][14][15] and analytical studies [27,28,38,39], the authors note a striking resemblance between the stress fields beneath shallow foundations and those of body force distributions in static sandpiles. The system-scale phenomenology of sandpile stresses may be best described by the FPA model [27,28]. In the FPA model, self-equilibrating tangential tractions corresponds to a unique body-force distribution that produces non-uniform normal contact pressures  in the plane of contact with a rigid base plate. Initial stress equilibrium states of the shallow foundation system are then conjectured so that (i) the very existence of internal friction causes the development of self-equilibrating tangential tractions beneath vertically loaded shallow foundations and (ii) these tangential traction fields develop as per the loading history leading to a 'saddle-shaped' distribution of normal contact pressures, in qualitatively similar ways to the stress states observed in the sandpile experiments. The notion of coupled normal and tangential traction boundary conditions of the foundation is borne of the sandpile analogy. The working hypothesis is that tangential tractions can be empirically correlated to measured normal contact pressures, for example, measured in [12]. Thus, self-equilibrating tangential traction components can be introduced as a secondary Neumann boundary condition of the vertically loaded shallow foundation problem.
In an earlier work [16], the authors suggested an empirical function prescribing the normal-traction boundary condition for the foundation problem, which was given as ð4:1Þ Here, A, B, z and v form a set of auxiliary functions which vary continuously between curve-fit parameters for the point-wise contact pressures measured by Murzenko [12]; the function a serves effectively as mapping from a circular to a square region. A semi-empirical expression of tangential tractions can be formulated in a similar manner. We have already discussed the analogy between the normal pressures beneath shallow foundations and those found in static sandpiles, epitomized by the analytical FPA model [27,28]. The lack of empirical data of tangential contact tractions beneath foundations necessitates assumptions on these self-equilibrating tangential contact tractions in the contact plane of the vertically loaded shallow foundation. It is assumed that the location of maximum tangential traction coincides with the occurrence of the peak normal pressure. In turn, the radial location j* of the peak normal pressure varies with respect to the angle of repose f where the ratio of the shear and normal stresses 6 ¼ s rz /s zz can be calculated, as shown in figure 4. In addition, the tangential traction distributions have zero values at both the centre and the edge of the contact plane.
In the following, we will make use of these hypotheses to prescribe a tangential contact traction field in relation to a given normal pressure distribution measured by Murzenko [12]. It is important to note that the present numerical example is intended only to illustrate the application of the present method of solution to foundation engineering problems where traction boundary conditions are assumed to be known a priori.

Definition of tangential traction boundary conditions
Similar to equation (4.1) for normal contact traction fields, a spatial function is selected to describe the hypothesized characteristics of the tangential traction distribution across the centre and diagonal lines of a contact interface between a vertically loaded rigid footing and supporting dense granular soil. Accordingly, a single-variable function can be defined as where C and D are free parameters for the a is the distance from the centroid to the edge of the loaded region. j is a local coordinate defined along the lines, as outlined in figure 5. Referring to the peak location of the normal traction j* , the free parameters of equation (4.2) are subjected to two conditions. Stated mathematically, the first condition is dq/djj j¼j * ¼ 0. The second condition is that the peak magnitude of the tangential traction, q( j*), is set equal to 6p( j*), where 6 is the proportionality at the peak location. It is physically plausible to let 6 be the stress ratio s rz /s zz of the FPA model. We obtain the peak locations across the centre and diagonal lines for the normal contact pressure data of Murzenko [12] as per equation (4.1) with the parameters given in table 2. These peak locations are then correlated with angles of repose from the FPA model using figure 4a. The angles of repose in the present example are intended as an intermediary curve-fit parameter of the working hypothesis, used to correlate point-wise measurements of normal contact pressure with statically admissible stress states of granular materials from the FPA model. They are not meant to be interpreted as physical material properties of Murzenko's sand. From figure 4b, the values of the peak stress ratio are found to be 6 0 ¼ 0.3429 and 6 45 ¼ 0.2890, where the subscripts refer to the angles across which the fitting occurs, that is the centre and diagonal lines of the loaded square region. In figure 6, the tangential tractions are shown for comparison with the curve-fitted normal contact pressure distributions. All relevant parameters for the generation of these expressions from equations (4.1) and (4.2) are given in table 2.
Next, let r and u be a set of polar coordinates over the loaded square region. Using the same procedure as outlined in [16], one can obtain a continuous expression of tangential radial traction q r (r, u): q r (r,u) ¼ C(u)r(a(u) À r)e ((a(u)Àr)=D(u)) 2 : ð4:3Þ Here, the auxiliary functions are defined as Surface tractions q x and q y , which are obtained from equations (4.3) and (4.4) for the parameters in table 2, are graphically illustrated in surface plots and contours in the Cartesian coordinate system in figure 7. Solutions to both the x-and y-directional surface traction boundary conditions can be superimposed with the solution to the normal traction boundary condition [16] to obtain a new stress equilibrium state.

Numerical results from combined boundary conditions
With a discretization scheme of 40 Â 40 subdomains across the square domain R S of the rigid footing, total z-directional displacement w at the contact boundary plane z ¼ 0 is investigated in relation to the combined influence of the tangential and normal surface traction boundary conditions. These influences are defined separately as follows: Here, the potentials x and V are those defined in equations (2.3) and (2.4); the subscripts (x, y, z) are used here to specify these potentials under the action of the Cartesian surface tractions (q x , q y , p). Subsequently, the total vertical displacement is written as wj z¼0 ¼ w q j z¼0 þ w p j z¼0 . We can define corresponding expressions for these surface displacements normalized by the elastic constants and the width of the loaded area as follows: @x y @y z¼0 and w Ã p ; > > > = > > > ; ð4:6Þ As shown in figure 8, the components of normalized vertical displacement defined in equations (4.6) are given as the cross-sections across the centre and diagonal lines of the loaded region. It is noted that the vertical displacement due to the applied tangential traction is of opposite sign to that resulting from the normal (compressive) traction. Therefore, the inclusion of the tangential tractions has the effect of Although the difference is small in the case of vertical displacements, which are mainly due to normal tractions, the presence of tangential tractions has an interesting effect on the vertical stress field s zz developed in the elastic body ( figure 9). In comparison with the vertical stresses obtained from the normal tractions alone [16], their distributions show more pronounced pressure dips for a larger region with respect to depth. This suggests that the tangential traction boundary condition of the linear elastic foundation system may be calibrated to account for the confining effect of geostatic stresses in the development of vertical resistances.
It is noteworthy that tensile horizontal stresses s xx develop in the centre within a shallow depth beneath the contact plane ( figure 10). This simulated surface tension counteracts the load effect of   applied vertical force in the elastic body. In the context of soil mechanics, the pre-stressed contact surface in tension can be interpreted as the confining effects everywhere within a shallow depth beneath a rigid footing (i.e. mobilized shearing resistance of a granular material against downward vertical displacements). However, the use of tangential traction boundary conditions to mimic the observed quasi-static bulk shear behaviours of granular materials [40] is subject to further validation. A systematic phenomenology is required that correlates the self-affine nature of granular soil at the underlying scales, for example, the topological and metric properties of fractal structures [41], to the system-scale contact phenomena.

Concluding remarks
In the present study, an analysis technique is developed to solve generalized Neumann boundary value problems (i.e. a homogeneous and isotropic elastic half-space subjected to an arbitrary distribution of combined normal and tangential tractions on its surface). The method for the proposed solution uses the interpolation and superposition of the closed-form solutions of potential functions for bilinear boundary conditions over rectangular loaded domains. The motivation has been computationally efficient approximate solutions to foundation engineering contact problems. The simulations presented here predict that the hypothetical presence of unique tangential traction fields impedes vertical displacement and produces pretension zones that lead to the development of shearing resistance. In addition, these zones alter the stress states of the supporting elastic body and tend to attenuate the rise of compressive vertical stress at the centre of the foundation. It is therefore evident that tangentially induced stresses in the contact plane are related not only to externally applied tangential loading, but also to internal friction (e.g. angle of repose [27,28]) and loading history [10,12].    The model presented here offers an improved tool for studying multidimensional contact evolution inside a loaded region subjected to combined loading conditions. For example, a comprehensive experimental study of contact stresses that is beyond the scope of this paper is necessary to characterize the range of physically plausible boundary conditions for the foundation engineering problem. The authors hope that the present solution's method, which is an extension of earlier work [16], will serve as a basis for researchers to explore a wide range of other applications in contact mechanics on various scales [42].
Data accessibility. All data and source code required to generate the results presented here have been included as the electronic supplementary material.
Author's contributions. The authors jointly conceived of the presented ideas. A.G.T. developed the theoretical formalism, performed the analytic calculations, developed the numerical implementation and generated the corresponding figures. J.H.C. verified the analytical methods and provided critical feedback and discussion at every stage. J.H.C. also provided the initial intellectual motivation, encouraged A.G.T. to investigate the problem and supervised the findings of this work. Both authors discussed the results and contributed to the final manuscript.
Competing interests. We declare we have no competing interests. Funding. This study received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

Appendix A. Library of closed-form potentials
The closed-form solutions required to apply our computational procedure to the solutions in this paper are listed here. We give the solutions in general for arbitrary rectangular regions R ¼ f(x 0 , y 0 , 0)ja 2 x 0 a 1 , b 2 y 0 b 1 g. Note that in the following, the notation []j a 1 a 2 j b 1 b 2 is used to compact the multiple repeated terms generated by the double definite integration by the fundamental theorem of calculus. In general, if F(x 0 , y 0 ) is a function of two variables, the operator acts upon it as follows: Ellipses ' . . . ' are used to carry on to the next line when the analytic expressions become too long. For the sake of brevity, we omit a number of expressions (e.g. @G 00 /@y, @ 2 G 01 /@y 2 ), which can be obtained through symmetry from other results given (e.g. @G 00 /@x, @ 2 G 10 /@x 2 ).