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

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.

Figure 1.

A description of the coordinate system and geometry of a rectangular loaded domain on the planar surface of a semi-infinite half-space.

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

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.

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)|a2≤x′≤a1,b2≤y′≤b1} 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 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 . Each individual subdomain has a length of Δx=|a1−a2|/M and width of Δy=|b1−b2|/N; the grid points at the corners of the subdomains are given by x=a2+(k−1)Δx, and y=b2+(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, with the following system of equations: Solving for the coefficients of Lagrange polynomials of first order per superscript indices m=0,1 for x′ and n=0,1 for y′, we have the following: 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 is defined as follows: In turn, the potentials corresponding to the polynomial component (x′)(y′) over R 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 , and both satisfy the Laplace equation within the whole of the half-space (i.e. ). The potentials are then written for the interpolated loads b(x′,y′) and By superposition, equations (3.2) and (3.3) are approximated over the entire domain R as follows: and 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 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 , 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 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={(x′,y′,0)|−a≤x′≤a ,−a≤y′≤a}. The constant coefficient 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 M2 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(x,y)=|wexact(x,y)−wapprox(x,y)|) show the same convergence rate of the boundary conditions in equation (4.1).

Figure 3.

Convergence of maximum error of the bilinear interpolant to the biquadratic function in equation (4.2), with respect to (a) the total number of subdomains, and (b) the area of subdomains. Results are for discretization schemes ranging from M=N=4 to M=N=50.

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

The generation of an empirical pressure surface

Considering a square loaded region R 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. A0, B0, ω0, σ0, A45, B45, ω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 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:

Figure 4.

Curve-fit representations of equation (5.1) with respect to Murzenko's experimental data [21] for the case where he reported an average pressure of approximately 1 kgf cm−2, (a) across the centre line, and (b) across the diagonal line of the square footing.

Table 1.

Input parameters for equation (5.3) fitting the Murzenko data in the case with reported average pressure 1 kgf cm−2.

A0	A45	B0	B45	ω0	ω45	σ0	σ45
1.7939	1.8292	0.2939	0.3292	0.4119	0.8114	0.2715	0.2853

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 , ), 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: where is the total force of the continuous surface traction of equation (5.3), and P 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.

Figure 5.

(a) A surface plot of a pressure distribution with pronounced pressure dip, derived to fit the data from Murzenko [21] (for the case with reported average pressure 1 kgf cm−2) by way of equation (5.3). (b) A contour plot of the same (note that the average pressure under this surface is calculated to be ≈1.095 kgf cm−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): 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. ). 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 of equation (5.4) is shown in figure 7.

Figure 6.

Convergence of the bilinear interpolant to the empirical boundary condition shown in figure 5, with respect to (a) the total number of subdomains and (b) the area of subdomains. Results are for discretization schemes ranging from M=N=4 to M=N=100.

Figure 7.

(a) Cross-sectional and (b) contour plots of the normalized displacement field of equation (5.4) specific to the vertical traction boundary condition of figure 5.

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

Figure 8.

(a) Cross-sectional and (b) contour plots of the resistance field of equation (5.5) specific to the vertical traction boundary condition of figure 5.

The vertical stresses in the elastic body are independent of the elastic constants from equation (A 12). Notably, σ 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.

Figure 9.

Variation of vertical stress at different depths under the semi-empirical boundary conditions (a) across the centre, and (b) across the diagonal of the loaded area. Compare with the results for a uniform average pressure on the boundary, for (c) the centre and (d) the diagonal. The results are normalized by the average pressure.

Figure 10.

Vertical strains with depth under points across (a) the centre line, and (b) the diagonal of the loaded region. We choose example elastic constants E=1000 kgf cm−2, ν=0.4 for this calculation.

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

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.