Application of low-order potential solutions to higher-order vertical traction boundary problems in an elastic half-space

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.


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] 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. 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].

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]: where 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: For a comprehensive description of complete solutions in terms of potentials, see §44 of reference [23]. The form of the potentials used to satisfy equation (2.1) is determined uniquely by the given boundary values.
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: The solution to this boundary value problem yields a vertical stress distribution with an absolute minimum at the centre and infinite stresses along the edges of 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 1 2π The resulting displacement and stress fields satisfy equation (2.1), are continuous within the body and lack singularities, providing that the surface tractions 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.

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][33][34]. Interestingly, these phenomenological observations and analytical predictions show remarkable resemblance to the normal stress distributions measured empirically in sandpile models [35][36][37][38] and corresponding analytical models [39][40][41][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   Figure 2. (Not scaled) Pressure variation for five loading cases extrapolated from data measured by Murzenko [21] at points (a) across the centre of a square footing, and (b) across its diagonal. This behaviour is assumed to be symmetrical across the rest of the footing, outlining a two-dimensional pressure surface.
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 (3.1) As tractions normal to the surface boundary are assumed to be zero outside the region of contact, the potential function can be written as a double integral over the geometry of the region R: This is the general solution to the Laplace equation subjected to the Neumann boundary condition in the half-space. To complete the equilibrium solutions of equation (2.1), it is necessary to introduce Boussinesq's logarithmic potential [1,2]: This function is related to equation (3.2) by 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.

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: for i = 1, 2, . . . , M, and j = 1, 2, . . . , N, such that R = R ij . Each individual subdomain has a length of 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: Solving for the coefficients c mn ij of Lagrange polynomials of first order per superscript indices m = 0, 1 for x and n = 0, 1 for y , we have the following: x y , An arbitrary p(x , y ) is then approximated over R via superposition of the interpolated loads over each subdomain: where the bilinear interpolant over subdomain R ij is defined as follows: In turn, the potentials corresponding to the polynomial component (x ) m (y ) n over R ij are defined as follows: and Equations (3.7) and (3.8) are special cases of (3.2) and (3.3), respectively, with specified traction distributions over a given subdomain. We then have that A mn ij = ∂B mn ij /∂z, and both satisfy the Laplace equation within the whole of the half-space (i.e. A mn ij = B mn ij = 0). The potentials are then written for the interpolated loads b ij (x , y ) and By superposition, equations (3.2) and (3.3) are approximated over the entire domain R as follows: The derivatives of the potentials in equations (3.11) and (3.12) can be constructed from the derivatives of those in equations (3.7) and (3.8), respectively, in the same manner. Appendix B contains the complete set of closed-form solutions to equations (3.7) and (3.8) and the spatial derivatives required to satisfy the stress, strain and displacement formulae presented in appendix A.

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: This function would appear in equations (A 4) and (A 10) for 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 The first term of the right-hand side of the above expression is then calculated as follows: A closed-form solution of the required second partial derivative of B 11 ij is then shown to be: where the derivatives of the potentials with respect to the lower-order pressure fields are already determined, as in appendix B. Owing to the symmetry of B 11 ij , ∂ 2 B 11 ij /∂y 2 can be immediately deduced as follows: A third result can be obtained from the calculation above by considering some facts about the potentials involved. For instance, we know that ∂ 2 B 11 ij /∂z 2 = ∂A 11 ij /∂z by the definition for the potentials in equations (3.7) and (3.8). Owing to the fact that these are harmonic functions, we can therefore write the following: This provides a closed-form solution for the expression appearing in the vertical displacement and strain, as well as all three normal stresses in appendix A. The result can be verified by direct calculation and/or relevant solutions from appendix B. All other calculations for the required closed-form solutions proceeded in analogous ways.

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): We have been unable to derive an analytical bound on the error function given in equation (4.1) or find one in the literature pertaining to bilinear interpolation. Thus, we opt to directly assess errors associated with the proposed interpolation scheme in a convergence study. Consider the following biquadratic pressure distribution: This distribution is defined over a square area, R s = {(x , y , 0)| − a ≤ x ≤ a , − a ≤ y ≤ a}. The constant coefficient 9 4 is conveniently used to normalize the pressure profile by an average pressure over the square region. Let 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

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 pointvalue 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].

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 ρ: This function was originally developed by Ai 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: This function varies from 1 to √ 2 such that 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:     (2) , (2) and σ (θ) = σ 0 f (θ) 2 log(σ 45 /σ 0 )/log(2).
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 θ = tan −1 (y /x ), ρ = x 2 + y 2 ), we can write the vertical traction as follows: The results of equation (5.3) are shown in figure 5. The function exactly predicts all five experimental data points (along with zero values at the boundary) with the parameters in table 1, and reveals a continuous traction field across the entire loaded domain. However, there is a discrepancy in comparison to the total force reported in the literature: a −a p(x , y ) dx dy is the total force of the continuous surface traction of equation (5.3), and 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.

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 where 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 (2a) of the square footing, we have the following: The normalized displacement field can be calculated by the sum of the potential functions of equation (3.9) over 40 × 40 discretized subdomains (i.e. w * (x , The loaded domain is discretized with this number of subdomains by trial and error to satisfy the convergence criterion presented in figure 6. We note that the higher-order boundary condition applied here shows a slower rate of convergence than the biquadratic function in §4. A normalized displacement field over 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 nondimensionalized resistance: The results from equation (5.5) are shown in figure 8. This expression of soil resistance can be viewed as a continuous distribution of elastic spring stiffnesses that is analogous to an extended Winkler foundation model. 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 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 intergranular 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 systemscale behaviour. The characteristics of granular materials can be further contextualized by many granular physics studies [47][48][49].

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][51][52]. A set of bilinear solutions to this problem for rectangular regions is currently under development by the present authors.

B.2. Linear load solutions
z log((y − y) + r) + (y − y) log(z + r) 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.