## Abstract

New solutions of potential functions for the bilinear vertical traction boundary condition are derived and presented. The discretization and interpolation of higher-order tractions and the superposition of the bilinear solutions provide a method of forming approximate and continuous solutions for the equilibrium state of a homogeneous and isotropic elastic half-space subjected to arbitrary normal surface tractions. Past experimental measurements of contact pressure distributions in granular media are reviewed in conjunction with the application of the proposed solution method to analysis of elastic settlement in shallow foundations. A numerical example is presented for an empirical ‘saddle-shaped’ traction distribution at the contact interface between a rigid square footing and a supporting soil medium. Non-dimensional soil resistance is computed as the reciprocal of normalized surface displacements under this empirical traction boundary condition, and the resulting internal stresses are compared to classical solutions to uniform traction boundary conditions.

## 1. Introduction

The classical solution to the problem of a semi-infinite homogeneous, elastic body subjected to vertical loads on its boundary surface was first developed by Boussinesq [1] and discussed in more depth by Love [2]. The solution is derived from the Green's function of the Laplace equation, which was used to determine the stress equilibrium within an elastic half-space. Other solutions can be also obtained using a separate approach involving Bessel functions [3,4]. The technique involving potentials was concurrently applied in line with the theory of elastic contact proposed by Hertz [5]. Further, Newmark [6] provided a concise closed-form solution to the distribution of vertical stress at the corner of a uniform rectangular load. A general solution to this type of boundary-value problem was presented earlier by Love [7], who developed expressions for displacement and stress under any integrable distributions of vertical tractions over an arbitrary domain. Closed-form solutions to uniform tractions over a rectangular area were explicitly presented in his paper (with the exception of those involving vertical displacement). Only partial solutions were provided for surface tractions that linearly and bilinearly vary in a rectilinear domain.

Numerous attempts have been made to complete the integral calculations presented in Love's original work and to further expand the application of his closed-form solutions. Notably, Ahlvin & Ulery [8] developed tables for stress, strain and displacement under a uniform circular load from Love's equations, but they did not produce new closed-form solutions of potential functions. Schmertmann [9] proposed a semi-empirical strain influence method that has been widely used for elastic settlement analysis of shallow foundations based on the assumption that traction fields are uniform at the contact interface between a rigid footing and supporting soils. Recently developed closed-form solutions of the potential functions include those reported by Becker & Bevis [10], who completed Love's discussion of displacement under a uniform rectangular load. Dydo & Busby [11] further discussed linear and bilinear variations in vertical traction fields over a rectangular contact domain and provided one of the most comprehensive sets of closed-form solutions of the potential functions to date. However, the derivatives required to calculate stress, strain and displacement were omitted from their results. Li & Berger [12] developed a corresponding set of closed-form solutions for constant, linear and bilinear tractions over triangular domains. Most recently, Marmo & Rosati [13] suggested a general solution to a problem of polynomial surface-traction conditions over areas defined by arbitrary polygons. Importantly, Kunert [14] developed closed-form solutions for Hertzian-like contact pressures varying over a rectangular area.

Of all the historical contributions to this problem, Boussinesq's point-load and Newmark's uniform surface-load solutions stand out as engineering solutions to foundation design problems. These stress-influence methods are widely accepted in foundation design practice, but Love's closed-form solutions have rarely been viewed as practical design tools. It is laborious to solve the integral equations of the potentials in closed form even for the simplest of polynomial traction boundary conditions. Further, contact traction fields relevant to modern engineering applications may not be of such a low order. As it stands, a new set of closed-form solutions must be developed for higher-order boundary conditions, for which in most cases the calculations are intractable. Alternatively, a numerical model could be designed to obtain an approximate solution to a desired degree of accuracy specific to a prescribed boundary condition, but such a modelling effort can be time-consuming and very costly in iterative design processes.

It is, however, possible to balance mathematical rigorousness and computational efficiency in an alternative solution. A high-order surface traction field can be discretized in the loaded domain using piecewise approximation of lower-order polynomials. Very recently, work has been done to solve elastic contact problems by the superposition of solutions to linear tractions over triangular regions [15]. This technique is analogous to the collocation method used in the indirect boundary element method [16,17]. The lower-order solutions for each discretized subdomain are superimposed to approximate the displacement and stress fields in the elastic body (i.e. the principle of superposition). However, in the present study, it is proposed to calculate the ‘boundary integrals’ analytically for rectangular regions. This way, the superposition of low-order closed-form solutions replaces costly numerical quadrature. This solution approach is not entirely new [18,19]. It has been used to model the micro-contact of rough surfaces using Love's constant load solutions [20]. However, superimposing solutions to only the constant and unidimensional linear boundary conditions cannot lead to a continuous expression across the boundaries of rectangular subdomains. Consequentially, singularities or discontinuities will appear in the resulting stress fields within the loaded domain. The missing component for a robust solution approach to a higher-order traction boundary-value problem using rectangular regions appears to be the closed-form solutions of the bilinear (hyperbolic–paraboloidal) potentials. The authors have not found a complete set of solutions of the bilinear potentials in the literature. The goals of this study are to close this gap in the literature and to provide an approximate yet general solution approach to the Boussinesq–Love problem for higher-order traction boundary conditions.

In this paper, a finite number of closed-form solutions for lower-order boundary conditions are superimposed to accurately approximate a solution for high-order surface traction fields acting on a rectangular surface area of a homogeneous, elastic half-space. We develop a complete set of closed-form solutions for the bilinear hyperbolic–paraboloidal potential functions and their required derivatives. The potentials corresponding to the bilinear interpolants fitting the values of the four corners of each of a given number of discretized subdomains are superimposed, providing an approximation of the potential for the prescribed boundary condition and ensuring continuity across local boundaries. For foundation analysis applications, a surface traction function is empirically formulated by curve fitting pointwise contact-pressure data from past foundation experiments in the literature. This is then prescribed as a boundary condition for the half-space problem. We solve for select displacements, stresses and strains that occur in the body. The predicted distribution of non-dimensionalized vertical resistances (i.e. contact stiffness) is calculated as the reciprocal of normalized vertical displacements in the contact plane. A unique resistance distribution results directly from the corresponding non-uniform vertical boundary traction, which appears to evolve with increased applied load [21,22].

## 2. Governing equations and boundary conditions

A set of right-handed Cartesian coordinates is defined so as to describe an elastic body as half-space bounded by a plane at *z*=0, where the positive *z*-axis points downwards into the body. Let (*x*,*y*,*z*) be an arbitrary point within the body or on the boundary, while (*x*′,*y*′) designates a location within the region of contact *R* as shown in figure 1.

The stress equilibrium of a homogeneous, isotropic and elastic body can be stated as a system of partial differential equations in terms of the displacements (*u*,*v*,*w*) in the directions (*x*,*y*,*z*) as follows [2]:
*D*=∂*u*/∂*x*+∂*v*/∂*y*+∂*w*/∂*z* is the strain dilation; *λ* and *μ* are the Lame's constants; and Δ is the Laplacian differential operator with respect to the spatial coordinates.

This problem reduces to a problem of potential theory when either surface displacements or tractions are given. This is achieved by means of Green's theorem and Betti's reciprocal theorem (see Love's treatise [2], Ch. X for details). The equilibrium state of the half-space is then given in terms of functions which satisfy the Laplace equation:

The boundary conditions are defined as the values of surface displacements or tractions, corresponding to Dirichlet and Neumann conditions for this solution, respectively. It may be reasonable to model the mutual interaction of the foundation and continuum as a mixed boundary value problem, as is required for rigid punch problems [24,25] and Hertz–Mindlin contact theories [5,26]. Solutions to problems regarding this class of boundary conditions are mathematically cumbersome, although work has been done to develop a general solution [27]. If the footing is assumed to be flat, rigid and symmetrically loaded normal to the contact plane, vertical surface displacement *w*(*x*′,*y*′) can be assumed to be uniform within the area of contact (i.e. a constant displacement (Dirichlet) condition on *R*). Subsequently, the lack of contact outside of *R* corresponds to a zero vertical traction, *p*, a Neumann condition elsewhere on the boundary. The boundary conditions for equation (2.2) under these assumptions are written as follows:
*R* [24,25]. The singular values in the stress distribution are non-physical, and the material cannot be in equilibrium.

By contrast, the simplicity and convenience of the stress-influence methods stem from an assumption of a unique distribution of vertical traction prescribed over the area of contact. The resulting Neumann boundary condition can be expressed as
*p*(*x*′,*y*′) are bounded and uniformly equal to zero on the boundary of *R* [7]. Our following discussion will be focused on the application of potentials to solving this particular class of boundary value problems.

## 3. Potential functions for arbitrary contact pressure distributions

Numerous empirical measurements have shown that the distribution of stress between a rigid structure and granular soil has relatively high-order spatial variation. Examples include Terzaghi's model [28], which is associated with general shear failure modes of shallow foundations in which a roughly parabolic distribution of contact pressure may develop. The laboratory test results reported by Bauer *et al*. [29] also show approximate parabolic pressure distributions on a scaled rectangular footing under an assumption of plane strain conditions. Murzenko [21] presents a set of contact pressure distributions that appear to be saddle-shaped (or shaped like the back of a two-humped camel) and correspond to pressure peaks occurring at a certain distance from the centre of the contact plane, a dip at the centre, and zero values along the edges, as shown in figure 2. Further, the contact pressure peaks observed in his experiments tend to move inwards towards the centre of the footing as applied load increases. This phenomenon has been reported in the results of a number of analytical models. Smoltczyk [22] presented an analytical expression for pressure boundary conditions mimicking this behaviour via statistical analysis, while Kerr [30] produced the same expression by introducing a shear membrane and another layer of springs into a Winkler-type [31] model. Furthermore, a number of pressure distributions with these attributes have been derived using elastoplasticity [32–34]. Interestingly, these phenomenological observations and analytical predictions show remarkable resemblance to the normal stress distributions measured empirically in sandpile models [35–38] and corresponding analytical models [39–42].

Considering these various contact pressure distributions in conjunction with the Neumann problem in equation (2.4), we formulate a general method of obtaining approximate solutions for any given surface traction while retaining the continuity of the displacement and stress fields within the body and upon its boundary. The simplest method of continuously approximating an arbitrary surface traction over a rectangular area is bilinear interpolation. Just as a two-dimensional boundary condition can be approximated to any degree of accuracy by a piecewise bilinear function, a solution specific to equation (2.4) can be accurately approximated by the superposition of each solution of the subdomain to a discretized Neumann boundary condition. These in turn are constructed as linear combinations of closed-form solutions for constant, linear and bilinear tractions, which are presented in appendix B.

Consider a given normal pressure *p*(*x*′,*y*′) on a rectangular region *R*={(*x*′,*y*′,0)|*a*_{2}≤*x*′≤*a*_{1},*b*_{2}≤*y*′≤*b*_{1}} on the surface of a half-space. The distribution of pressure can be defined piecewise by a combination of any arbitrary functions. Using the coordinate system defined in figure 1, the distance between an arbitrary point in the body and a point on the surface within the loaded region is given by
*R*:
*V* =∂*χ*/∂*z*, and both satisfy the Laplace equation (i.e. Δ*V* =Δ*χ*=0). The complete three-dimensional displacement, stress and strain fields are determined by these two potential functions and their derivatives, as shown in appendix A.

### 3.1. Superposition of potentials

Let us divide the loaded area *R* into *M* and *N* uniform intervals in the directions of *x*′ and *y*′ , respectively, which creates *M*×*N* disjoint rectangular subdomains:
*i*=1,2,…,*M*, and *j*=1,2,…,*N*, such that *x*=|*a*_{1}−*a*_{2}|/*M* and width of Δ*y*=|*b*_{1}−*b*_{2}|/*N*; the grid points at the corners of the subdomains are given by *x*_{k}=*a*_{2}+(*k*−1)Δ*x*, and *y*_{l}=*b*_{2}+(*l*−1)Δ*y* for *k*=1,2,…,*M*+1 and *l*=1,2,…,*N*+1.

Bilinear interpolation of the values of *p*(*x*′,*y*′) is performed at the four corners of each *R*_{ij}, with the following system of equations:
*m*=0,1 for *x*′ and *n*=0,1 for *y*′, we have the following:
*p*(*x*′,*y*′) is then approximated over *R* via superposition of the interpolated loads over each subdomain:
*R*_{ij} is defined as follows:
*x*′)^{m}(*y*′)^{n} over *R*_{ij} are defined as follows:
*b*_{ij}(*x*′,*y*′)
*R* as follows:

### 3.2. Example calculation for bilinear boundary conditions

Next, let us consider a brief example of the calculation strategy that led to these expressions (in contrast with those reported in [7,11]). Consider the following function:
*x*-directional strain and stress, respectively. It is convenient to pair the spatial and integration coordinates as they appear in the distance function in equation (3.1), as substitution allows one to cancel the integrals and derivatives in the calculations with a change of sign. As (*x*′−*x*)(*y*′−*y*)=*x*′*y*′−*yx*′−*xy*′+*xy*, we see that we may write equation (3.13) as

## 4. Convergence and error assessment

Convergence of the proposed discretization scheme is investigated for interpolation of an arbitrary pressure field and associated error. We define an error function across the discretized loaded region as the differences between an exact boundary condition and the approximation in equation (3.5):
*R*_{s}={(*x*′,*y*′,0)|−*a*≤*x*′≤*a* ,−*a*≤*y*′≤*a*}. The constant coefficient *b*(*x*′,*y*′) be the piecewise bilinear interpolant of *p*(*x*′,*y*′) over *M*^{2} subdomains, as defined by equations (3.5) and (3.6). Figure 3 shows the convergence rate of the bilinear interpolant to the boundary conditions given by equation (4.2) with respect to the total number of subdomains and their area versus the maximum error values in equation (4.1). Similarly, we find that the solution errors relative to the exact solution for the vertical surface displacements, (i.e. *e*_{w}(*x*,*y*)=|*w*_{exact}(*x*,*y*)−*w*_{approx}(*x*,*y*)|) show the same convergence rate of the boundary conditions in equation (4.1).

## 5. Numerical example

Having established a procedure for determining stress and displacement in an elastic half-space under arbitrary pressure boundary conditions, we present a numerical example of its application. When external forces act on a rigid footing with a rough surface resting on a granular material, an apparent area of contact forms along the geometry of the foundation. The resulting equilibrium state, and ultimately the foundation settlement, are dependent upon the system variables (e.g. geometry and loading history [35]) in relation to contact phenomena that occur at the interfaces of the footing and supporting granular soils. To produce a mathematical description within the proposed scope of the present study, we deliberately reduce the multitude of soil–structure interaction phenomena to a single quasi-static boundary-value problem referring only to phenomenological observations at the foundation scale. We assume a stress-free reference as per the theory of elasticity; the theory applied here assumes no body forces and is therefore incapable of modelling geostatic stresses. In the following, we carefully develop a mathematical expression for vertical traction boundary conditions based on laboratory-scale experimental data available from the literature.

Recall the discussion of boundary conditions in §3. Of particular interest are Murzenko's results [21], which exhibit a dip in pressure at the centre of the loaded domain, as shown in figure 2. In the following, we empirically describe a two-dimensional vertical traction distribution that exactly fits these point-value measurements over a square contact area. We also impose zero values of normal traction along the edge of the square region. This enforcement both satisfies zero shear resistances along the edges of the surface footing and maintains continuity in the traction conditions across the entirety of the boundary plane. Otherwise, the discontinuities across the boundary of the loaded region lead to singularities in the resulting stress field [7].

### 5.1. The generation of an empirical pressure surface

Considering a square loaded region *R*_{s} of half-width *a*, we begin to empirically prescribe traction boundary conditions with a parametrized function of a single variable *ρ*:
*et al*. [38] to interpolate stress data measured in sandpiles. As already stated, there is a marked resemblance between Murzenko's data and these sandpile stresses with respect to shape and variation. Based on the results of this review, it was determined that a simple curve fitting of Murzenko's data by equation (5.1) provides a rudimentary approximation of the boundary conditions. The values of the free parameters *A*, *B*, *ω* and *σ* are calculated to fit equation (5.1) for empirical data regarding contact pressure distribution. The four data points (including the zero edge values) across the centre line or diagonal of a square domain are paired with these unknown parameters, which can be solved with a nonlinear equation solver (e.g. the Matlab ‘fsolve’ function).

With reference to the curve fitting across the 0° (centre) and 45° (diagonal) angles of the domain, we define two sets of parameters, i.e. *A*_{0}, *B*_{0}, *ω*_{0}, *σ*_{0}, *A*_{45}, *B*_{45}, *ω*_{45} and *σ*_{45}. We select the case which Murzenko reported as having an average pressure 1 kgf cm^{−2}, depicted in figure 2, for curve-fitting using equation (5.1). We report that there are a number of combinations of parameter values for equation (5.1) which fit the values with slight variations in the shape of the curve; in other words, the curve-fit is not unique. The selected distributions are shown in figure 4, and the relevant parameters are listed in table 1. Although this parametric fitting provides us some understanding of the contact phenomena, it is unknown exactly how the contact stresses vary between these two lines. In an attempt to continuously describe variations between the two lines in a Cartesian domain over a square area, we use coordinate transformation from a mapping function:
*f*(*θ*) is the distance from the centre of a unit square to its boundary at angle *θ*. Taking the variable *ρ* to be radial, and letting *α*(*θ*)=*a*⋅*f*(*θ*), we map the line load to a square domain by letting *a*→*α*(*θ*) in equation (5.1). In turn, functions of *θ* can be defined for all the curve-fitting parameters given in equation (5.1) as follows:

These functions vary continuously from a given curve-fit parameter at 0° to that at 45°. Transforming these auxiliary functions to a Cartesian coordinate system (i.e. letting *P*_{M} denotes the applied load as per the average pressure, 1 kgf cm^{−2}, reported in [21]. There are a number of possible reasons for this discrepancy, one being the extrapolation of the boundary traction between the centre and diagonal lines. However, Murzenko recognized similar discrepancies in total force between his curve-fitted surface tractions and the measurements of applied load, although his method of obtaining a continuous distribution of contact pressures is unknown.

### 5.2. Calculations for displacement, strain and stress fields

We note that the closed form solutions in appendix B may appear to be undefined at the boundary *z*=0. However, given the boundary conditions presented here, these expressions all tend to finite limits on the boundary. The values of the limits at these points are applied so that the displacement, stress and strain fields are all continuous at all points within the body and upon its boundary. Specific to the traction boundary condition described by equation (5.3), the vertical surface displacements are solved using equation (A 3):
*E* and *ν* are the Young's modulus and Poisson's ratio of the elastic material, respectively. The potential *V* is given by equation (3.2). After normalization with respect to the elastic constants and width (2*a*) of the square footing, we have the following:
*R*_{s} of equation (5.4) is shown in figure 7.

Recalling the discussion in §2, a symmetrically loaded rigid plate resting on a homogeneous material will experience uniform surface displacement, while the Neumann problem outlined in equation (2.4) will generally yield non-uniform surface displacements for a given traction distribution. Considering the compatibility requirement, we interpret the results of equation (5.4) as the distribution of non-dimensionalized resistance:

The vertical stresses in the elastic body are independent of the elastic constants from equation (A 12). Notably, *σ*_{zz} at *z*=0 is exactly equal to, but opposite in sign of, the described vertical tractions. Figure 9 compares vertical stresses along depths beneath the loaded area with the classical solution obtained from a uniform pressure boundary condition [6,7]. The non-uniform pressure boundary condition yields higher stresses near the centre and lower stresses towards the edge of the loaded area. In addition, the vertical strains produced within the body from the applied tractions can be obtained from equation (A 6). For illustration purposes only, vertical strain distributions along depth are plotted at the locations of Murzenko's pressure measurements in figure 10. The values of Young's modulus and Poisson's ratio are arbitrarily selected (e.g. values of tri-axial compression tests on very dense sand found in soil mechanics textbooks, [43,44]). Similar to strain-influence methods [9], the definite integration of each curve from the surface to an infinite depth yields displacement from equation (A 3) at Murzenko's sampling locations. The applied surface tractions yield an interesting trend in vertical strain fields; tensile strains are pronounced beneath the corners mainly due to the Poisson effect and the zero traction (resistance) boundary condition imposed on the edges of the loaded domain.

The results for displacement, stress and strain, and in particular the resulting non-uniform distribution of resistance/spring stiffness, correspond uniquely to the traction boundary condition, which is largely dependent upon the nature of the granular soil. The proposed solution could be further refined to incorporate mobilized shear resistance of the particles, such as the discrete models presented by Pasternak [45] and Kerr [30]. In particular, geostatic stress states may introduce another set of initial and/or boundary conditions to foundation systems.

It is important to note that this particular distribution of contact forces *p*(*x*′,*y*′) should not be considered conclusively supportive of a particular model of stress propagation in granular materials, either hyperbolic or elliptic [46]. It is also true that describing the supporting soil medium as an isotropic, homogeneous elastic material with a stress-free reference state is unlikely to produce contact stress distributions like those measured in Murzenko's experiment and/or the sandpile tests. In fact, there is no comprehensive theory that can predict contact force distribution. The stress distribution is probably the phenomenological outcome of numerous multiscale parameters. The presence and degree of inter-granular friction bonding is dependent upon grain-scale geometry, relative density state, dilation and past loading history, which in turn determine the bulk constitutive phenomena of the granular media. For these reasons it was intended only to phenomenologically prescribe a Neumann boundary condition to exemplify the proposed solution procedures. The theoretical question of whether these physical factors may be all together incorporated in the calculations remains unanswered. With more empirical data at various interrelated scales, more detailed analytical models must be developed to describe the system-scale behaviour. The characteristics of granular materials can be further contextualized by many granular physics studies [47–49].

## 6. Concluding remarks

This study has aimed to provide a general solution approach to the classical Boussinesq–Love problem for arbitrary loads over a rectangular surface of an elastic half-space, described as Neumann boundary conditions. It is very challenging to adequately describe contact phenomena in foundation systems because contact stress distributions within the loaded area of even purely elastic bodies is neither uniform nor linear [5,24,25]. The solution procedure presented here, along with the closed-form solutions of the potential functions, are readily applied to arbitrary traction boundary conditions in a Cartesian coordinate system. We attempt to apply these potentials to foundation engineering situations, specifically to elastic settlement analysis of shallow foundations. However, we hope that the ubiquitous nature of the Laplace equation allows the present solution to be feasibly applied to other fields of study. Further, the proposed method of solution can be applied to tangential traction boundary-value problems [50–52]. A set of bilinear solutions to this problem for rectangular regions is currently under development by the present authors.

## Ethics

The authors testify that this material has not been published in whole or in part elsewhere; the manuscript is not currently being considered for publication in another journal; the authors have been personally and actively involved in substantive work leading to the manuscript, and will hold themselves jointly and individually responsible for its content.

## Data accessibility

All data and source code required to generate the results presented here have been included as the electronic supplementary material.

## Authors' 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.

## Acknowledgements

The corresponding author expresses gratitude for support from the Bridge Software Institute at the University of Florida.

## Appendix A. Formula for the displacement vector, stress and strain tensors

The displacements (*u*,*v*,*w*) in the (*x*,*y*,*z*) directions are given by
*λ* and *μ* are the Lame's constants. The strains are derived simply by differentiation as follows:

## Appendix B. Closed-form solutions of potential functions

We here list the closed-form solutions of the potentials for arbitrary rectangular region *R* under the simple polynomial loads employed throughout this work. We define
*F*(*x*′,*y*′) is a function of two variables, the operator acts upon it as follows:

**B.1. Constant load solutions**

**B.2. Linear load solutions**

The solutions for the potentials A^{01}, B^{01} for a linear load in the *y*′ direction can be obtained directly from these solutions by permutation of the terms related to the *x* and *y* variables in the above expressions.

**B.3. Bilinear load solutions**

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4085885.

- Received February 6, 2018.
- Accepted April 3, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.