Elastic Solutions to 2D Plane Strain Problems: Nonlinear Contact and Settlement Analysis for Shallow Foundations

: The classical Neumann boundary value problem of an isotropic, homogeneous elastic half-plane under plane strain conditions is readdressed as the limiting case of the fully three-dimensional problem. Analytical solutions of the stress and strain tensors are obtained by taking the limit from known three-dimensional solutions. It is shown that the displacement ﬁelds for the plane strain problem are not well deﬁned. A small number of simple expressions are developed, which provide a general solution for linearly-varying traction over arbitrary regions on the boundary. A simple, efﬁcient, and rapidly convergent algorithm is developed which uses these solutions as analytic elements and provides a solution approach to the general boundary value problem. The method is veriﬁed against known solutions for Hertzian contact between parallel cylinders. Two numerical examples are presented for the analysis of shallow foundation systems. In the ﬁrst, the boundary conditions are informed by analytical elastoplastic calculations and a strain inﬂuence analysis is performed and compared with the Schmertmann method. Subsequently, empirical laboratory contact traction distributions measured by Bauer et al., in both the normal and tangential directions are employed as boundary conditions for an analysis of the underlying stress ﬁeld.


Introduction
The problem of an isotropic and homogeneous elastic half-plane under plane strain with given conditions on its linear boundary has been well-studied for over a century. Various equivalent open-form integral solutions have been obtained for the Neumann problem, where the distributions of boundary traction are prescribed. The problem was first solved by Flamant [1] for concentrated point forces; boundary integral representations of the stress and strain fields under distributed traction can be obtained by making use of the classical Flamant solution as a fundamental solution (e.g., [2]). Alternatively, the solutions can be obtained directly from the governing equations using Fourier transforms [3] or by manipulating potentials in the complex domain [4,5]. Analytical closed-form solutions of these fields, which are easy to evaluate at any point of interest in the material domain, are available only for a select few, simple cases of interest [2,4].
The plane strain problem in linear elasticity has a history of application in the analysis of soil-structure interaction problems. The assumptions that the out of plane strain is zero is suitable for long strip foundations, allowing for analysis of only the center line of contact with the underlying soil and reducing the problem to from three to two dimensions. There are clear limitations to representing a granular soil body as an isotropic and homogeneous linear-elastic solid. The approach represents perhaps the most extreme simplification, in which the complex interaction of discrete grains are homogenized by a representative elastic continuum. However, the ability to solve for closed-form expressions has made elastic solutions very attractive for quick approximations of the stress, strain, and settlement of foundation systems. Boussinesq, who pioneered the earliest closed-form solutions for boundary value problems [6], also wrote explicitly about the mechanics of granular soils decades before soil mechanics were properly formulated in the early 20th century [7]. Later, when Love calculated some of the first closed-form solutions for distributed normal tractions, he framed his work as "an attempt to throw some light on the important technical question of the safety of foundations [8]". These fundamental elastic solutions became the foundation of many methods assessing the equilibrium and settlement of shallow foundation systems. Newmark used these solutions as the basis of his popular stress influence method, providing graphical charts to aid in engineering calculations [9]. The famous strain influence method of Schmertmann [10,11] is similarly based on tabulated values of the strain fields calculated in an elastic half-space. More recent works continue to use fundamental calculations from linear elasticity to assess the settlement of foundation systems [12][13][14][15][16][17]. This list is far from exhaustive.
Recent work [18] has suggested that using physically-accurate boundary conditions (that represent the stress states between footings and real granular materials) can import some of the nonlinear behavior due to the discrete granular assembly into the linear elastic analysis. This approach may lead to a better representation of the stress and strain fields within the loaded elastic body representing an underlying soil. In order to conveniently study the behavior of these fields and their dependence on the boundary data for the plane strain case, mathematical tools are required to represent them (at least approximately) within the half-plane. This could be achieved through either direct analytical or numerical integration of the known open-form integral solutions. The former approach is tedious or intractable for boundary traction fields of arbitrary complexity. On the other hand, pointwise numerical integration of the boundary integrals can be computationally expensive and prohibits an in-depth study of the behavior of stress and strain fields under evolving boundary conditions.
The present paper provides simple analytical tools for the representation of the tensor field solutions for the Neumann plane strain problem. New simplified closed-form solutions are presented for low-order monomial boundary conditions; these are obtained by taking analytic limits from known closed-form solutions of the 3D problem [16,17]. Surprisingly, a relatively small number of simple analytical expressions are found to completely determine the problem. These expressions are then applied as analytic elements to form continuous approximate solutions to higher-order and arbitrary boundary conditions. Applications to geotechnical problems of soil-structure interaction are then presented, which involve analytical and empirical boundary data.

Plane Strain as a Limit of the Fully Three-Dimensional Half-Space Problem
Recent work has been done to provide analytical closed-form solutions for the elastic stress, strain, and displacement fields on a fully three-dimensional isotropic and homogeneous elastic half-space under polynomial Neumann boundary conditions [16,17,19,20]. Consider the half-space H 3 = {(x, y, z) ∈ R 3 |z > 0} as described by the three-dimensional Cartesian coordinate system shown in Figure 1. For static conditions without body forces, the solid is ultimately governed by the system of elliptic partial differentials equations [21]: Here λ, µ are the Lamé elastic constants, (u, v, w) are the elastic displacements in the Cartesian (x, y, z) directions, D = ∂u ∂x + ∂v ∂y + ∂w ∂z is the strain dilation, and ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 + ∂ 2 ∂z 2 is the Laplacian differential operator. The problem is defined when boundary traction functions q x , q y , p are prescribed (acting in the x, y, and z directions, respectively) over a subset R ⊂ ∂H 3 . The half-space boundary ∂H 3 consists of the points in the plane z = 0, and R can in general be the union of any disjoint subsets, outside of which the boundary conditions are assumed to be zero. It is also assumed that the displacements, and consequently the strains and stresses, go to zero as distance from the source region R becomes large. Based originally on the work of Boussinesq [6] and Cerruti [22], the displacement, stress, and strain fields under monomial traction distributions of order m, n can be constructed from three potential functions from H 3 → R: The superscripts mn on the three-dimensional potentials (2)-(4) correspond to the order of the boundary conditions given spatially by monomial terms (x ) m (y ) n . These functions satisfy the three-dimensional Laplace equation, i.e., ∆A mn = ∆B mn = ∆Γ mn = 0. It is also relevant to note that the three functions A mn , B mn , and Γ mn are not fully independent but are related by the following differential operations:

Calculation Approach
The potentials corresponding to the plane strain problem can be obtained from Equations (2)-(4) by taking limits as one dimension of the boundary region R becomes infinite. Let R = {(x, y) ∈ ∂H 3 |a 2 ≤ x ≤ a 1 , −b ≤ y ≤ b}, such that the integrals in Equations (2)-(4) are taken over a rectangle centered at the origin in the y-direction but which lies on an arbitrary cut of the x-axis. Geometrical intuition assures us that as b becomes arbitrarily large, the potentials will exhibit less dependence on the spatial variable y; as b becomes infinite, every y value, having infinite length of rectangle to either side of it, becomes equivalent. This also corresponds to the intuitive geometrical definition of plane strain problems. The potentials which solve the elastic plane strain problem in the xz-plane for monomial boundary traction fields varying with x can formally be defined as: Here, j, k are simply the (integer) orders of the x, z, partial derivatives, respectively. As will be shown below, these expressions do not converge to well-defined limit functions for some values (e.g., j = k = 0). However, the derivatives required to define the stress and strain tensors (see [16,17]) are all well-defined. All partial derivatives with respect to the y variable tend to zero in the limit, while the convergent functions (5)-(7) are independent of y. The analytical properties discussed earlier for Equations (2)-(4) (namely that they satisfy the Laplace equation and are interrelated via differentiation) all hold in the limit for Equations (5)- (7). The following identities can be derived: The result is that the stress and strain tensors for the plane strain problem can be completely determined by only four potential functions for each case, expressible in terms of Equation (5) alone. The stress component fields when the normal traction is p(x ) = (x ) m and the tangential traction q(x ) = 0 are obtained from those for three-dimensional problem given in [16] and are listed as follows: The superscript m again indicates the order of the monomial boundary traction, while the subscript p indicates the direction of the nonzero traction component for Equations (8) Thus, the stress tensor is dependent only on the boundary conditions and the Poisson ratio ν. The stress components for the case where tangential boundary traction is q(x ) = (x ) m and the normal traction p(x ) = 0 are given as follows: The components of Equations (12)- (15) can be superimposed with those of Equations (8)- (11) to obtain the stress fields for plane strain boundary value problems with nonzero p(x ) and q(x ). The required partial derivatives of Equations (2)-(4) are known in closed-form for constant and linear traction fields; the potentials required to satisfy (8)- (15) can be obtained by direct calculation from these closed-form solutions by taking a limit as in Equation (5). For example, it is shown in [16] that for the rectangular region R defined here, By directly calculating the limit of this expression as b → ∞, it can be shown that The solutions for various m can be combined in linear combination to produce solutions for more general polynomial loads.

A Note on the Displacement Field
The discerning reader may question why the solutions here in Equations (8)-(15) are given in terms of the stress components, while the 3D governing Equation (1) are formulated in terms of displacements. In fact, there is some ambiguity in the literature surrounding the specification of displacement fields corresponding to the plane strain Neumann boundary value problem discussed herein. Using the notation of the prior section, it is clear that in the limit associated with the plane strain conditions, v = 0 as a consequence of yy = 0. The governing Equations (1) become By the same argument, the expected the solutions for the remaining nonzero displacement components can be written in terms of Equations (5)-(7) for monomially-distributed traction components p and q: However, a potential defined from Equation (5), namely A m , does not tend to well-defined limits as b → ∞. This apparent lack of well-defined displacement fields for the plane strain elasticity problem appears in the literature and is dealt with in different ways. Johnson [2] deals with the issue by defining the surface displacements in terms of an indefinite integral and presenting results in terms of an arbitrary constant, which he suggests can be chosen to make the displacement field equal zero at some distance from the load. However, the displacement fields defined in this fashion grow without bound as x → ∞, contrary to the assumption that all stress, strain, and displacement fields go to zero as the distance from the source becomes large. Johnson also suggests the study of the rate of change of the surface vertical displacement in the horizontal direction (i.e., ∂w ∂x | z=0 in the present notation), which gives an idea of the relative displacements without the choice of an arbitrary datum. The lack of well-defined boundary displacements is also pointed out by Hemsley [3] in his treatment of the problem involving Fourier transforms and the biharmonic equation. A recent paper calls attention to the contradiction in the context of the problem of indentation of a half-space by a long elastic cylinder [23], in which the normal displacement is illdefined. Those authors highlight a logarithmic dependency on cylinder radius in the solution of the problem of parallel contact of two cylinders, of which the former problem is obtained from a limit as a cylinder radius becomes large. The analytic solutions of Equation (2) [8,16] have a similar logarithmic dependency on b, which makes the taking of the required limit for A m impossible.
It is common practice in the geotechnical community to approximate surface displacements by the integration of strain curves up to an assumed influence depth [10]. Analytical expressions for the internal strain fields under normal load can be obtained from the stresses in Equations (8)-(11) via Hooke's law: for a chosen depth d, the surface displacement w 0 could therefore be approximated as However, this also involves a more-or-less arbitrary datum, as this integral grows without bound as d becomes large. This is equivalent to stating that Equation (17) is ill-defined. In general, it appears that a suitable displacement field, i.e., one that satisfies the plane strain governing equations of elasticity and tends to zero with distance from a source load, simply does not exist. The calculation approach employed here, namely the calculation of closed-form plane strain potentials from fully three-dimensional potentials through the limit procedure, may provide some intuitive insight as to why. As the loaded strip outlined in Figure 1 becomes of infinite length and fixed width with constant force per area, the total force applied to the system becomes infinite in the limit. While the stress and strain fields tend to finite and well-defined limit functions, the displacements grow without bound, as illustrated for a constant traction in Figure 2. Therefore, this lack of definition for the displacement fields appears to be a fundamental feature of the plane strain Neumann boundary value problem, rather than the result of some failure of our analysis. A demonstration that the magnitude of the three-dimensional potential A 00 (Equation (2)) grows without bound as the dimension b of the loaded rectangular region R becomes large. Thus Equations (17) and (18) are not well defined for constant surface traction.

Approximate Solutions to Higher-Order Problems Via Superposition
The analytical solutions which provide the plane strain stress and strain fields for loworder (constant and linear) monomial traction boundary conditions in the Appendices A.1 and A.2 can be applied to construct approximate analytical solutions to general boundary value problems. The linearity of the governing partial differential equations assures that the superposition of solutions under varying boundary conditions provides a solution for more complex boundary data. The uniqueness of the solutions then assures that a determined series of solutions which converges to a given expression on the boundary converges to the solution of the corresponding boundary value problem within the problem domain (see [16,17]). The plane strain assumption makes it so this can be achieved through simple linear interpolation of arbitrary p(x ), q(x ). This is shown by analytically approximating functions of the form where f (x ) is an arbitrary, bounded, continuous (C 0 ) function representing p or q.

A Simple Algorithm
Discretize the region a 2 ≤ x ≤ a 1 into N uniform subregions of width dx as shown in Figure 3; the N + 1 nodal interpolation points are then N . The necessary derivatives of Equation (21) can be analytically approximated as a sum of N linear combinations of solutions of those of Equation (5) The relevant interpolation constants c m i can be written as This procedure is valid for all partial derivatives of the required potential functions, and also for arbitrary q(x ). It provides a simple but general means of approximating solutions to arbitrary boundary value problems via the linear combination and superposition of the few and simple expressions in the Appendices A.1 and A.2. The simplicity of the linear interpolation scheme also allows for ease in error estimation. For example, let e p (x) be the error function between a normal traction distribution p and its linear interpolation as shown in Figure 3; it is known (e.g., [24]) that a sharp error bound can be given by The uniqueness of the elastic field guarantees that the internal conditions converge alongside the conditions on the boundary. The present approximation approach can be categorized as a numerical method as an instance of the Analytical Element Method (AEM) in elasticity. This method was originally developed by Strack [25,26] for the analysis of groundwater flow problems. In essence, AEM is equivalent to the Boundary Element Method (BEM), except the boundary integrals, dependent on the boundary conditions, are evaluated analytically. The closedform solutions in the Appendices A.1 and A.2 can therefore be considered analytic elements for plane strain elastic contact problems.

Verification Problem: Elastic Contact of Parallel Cylinders
The plane strain condition serves as an approximation of the problem of two long parallel elastic cylinders in contact. The contact for an arbitrary cross-section is also equivalent to the limit of a classic Hertzian contact [27] as the out-of-plane dimension becomes large. Similarly, the stress and strain fields can be considered as limit of the threedimensional solutions given by Kunert [28] for Hertzian contact traction over a rectangular area. The simple algorithm presented in the last section can be verified in comparison to exact solutions for this classical problem.
Assuming a frictionless contact (i.e., q(x ) = 0), the contact traction field for an area of width 2a, centered at the origin, is given by Equation (26) is shown in Figure 4a alongside its linear intepolation. Closed-form solutions for the internal stress fields beneath the center of the contact area are given in [29]. Figure 4b compares these internal stress fields with those approximated by the present method for N = 20 discretized subregions.   (26), along with a linear interpolation consisting of N = 20 discretized subregions. (b) Comparison of the exact internal stress fields given in [29] with those approximated using the superposition of 20 potentials corresponding to the linear interpolants. The out-of-plane stress σ yy is calculated for a Poisson's ratio ν = 0.3.
The global error bound given by Equation (25) can not be calculated due to the fact that the slope of the tangent line of Equation (26) is infinite at the edge of the contact area; however, the stress fields still show strong convergence to the exact analytical solutions. Figure 5 shows that the the maximum error in the stress approximations of Figure 4 converge to zero with increasing N. The rate of convergence of the stress fields given by the present method is facilitated by the two-dimensional nature of plane strain problems. It is noted that simple Matlab code converges quickly for small N, and is noticeably faster than similar code employing bilinear interpolation for the fully three-dimensional contact problem (see [16,17]).

Applications to Shallow Foundation Analysis
The classical Boussinesq-Cerruti half-space problem has historically been applied to the analysis of soil-structure interaction, in particular the settlement and stress states in foundation systems [8,9,16,17,30]. The plane strain problem has been used as an approximation of strip foundations, wherein the aspect ratio of the footing is often large enough to neglect the length dimension. However, the contact and penetration behavior of granular soils is highly complex, and the observed contact traction and displacement fields generally vary non-linearly with a large range of multi-scale material and loading variables. There are therefore a large range of possibilities when selecting boundary conditions for this type of analysis; the generality of the calculation approach outlined here allows for the calculation and comparison of results for evolving contact traction fields from a wide array of sources. Two examples are chosen to highlight this application. The first involves an evolving normal contact traction field obtained analytically through an elastoplastic analysis [31]. This is used to investigate the strain response due to an increasingly loaded foundation. The second example applies traction fields which have been measured experimentally by Bauer et al. [32], acting in both the normal and tangential directions. These empirical traction fields are used to investigate the static stress behavior beneath a loaded strip foundation.

Example 1: Foundation Contact as an Elastoplastic Rigid Punch Problem
A loaded shallow foundation system is naturally modelled as a rigid punch problem, due to the fact that the footing is often relatively rigid in comparison to the underlying granular soil body. This type of problem is initially formulated using a mixed boundary condition [33][34][35], with displacements prescribed inside the contact region, and zero traction conditions applied outside on the surface of the half-space. Once the resulting surface traction fields are obtained corresponding to the prescribed displacements, equivalent solutions can be obtained by prescribing purely Neumann (traction) boundary conditions.
Boussinesq obtained the first solution [6,35] for a flat punch and showed that for purely elastic conditions, uniform surface displacement leads to singularities (infinite values) in the normal traction field at the edge of contact for any magnitude of indentation. Thus, any real elastic material will yield locally when indented any amount by a flat rigid punch. Schultze [36] initially combined Boussinesq's result for elastic plane strain contact with the classical geotechnical theory of bearing capacity [37]; this led to a contact traction distribution which mimicked Boussinesq's in a central elastic region but varied linearly to zero in an external plastic region. The contact traction therefore evolves with increasing load magnitude. Balakrishna [38] later extended this elastoplastic analysis to the three-dimensional case of an axisymmetric circular punch. Abdullah [31] later augmented Schultze's results so that the traction in the plastic region varies quadratically. This was shown to agree better with numerical results; it also coincides more accurately with Terzaghi's initial result for the traction at ultimate bearing capacity [37]. For these reasons, Abdullah's traction field will be employed in the present example. Evolution of the contact traction from a purely elastic state to a quadratic plastic state at ultimate bearing capacity is shown in Figure 6 for a surface footing on a dry, dense sand.  Figure 6. Example normal contact traction field distributions for an elastoplastic rigid punch foundation loaded on the surface of a dry, cohesionless sand (from Abdullah [31]). The traction field evolves inwards towards the bounding quadratic curve for a load P ult , determined by the classical bearing capacity theory [37]. The traction fields are normalized by the ultimate average pressure p ult = P ult /2a.
The evolving traction fields shown in Figure 6 suggest a nonlinear displacement response of the foundation with increasing load, and hence a nonlinear elastic stiffness path. Recall from Section 2.1.2 that the displacement fields are not uniquely defined for the boundary conditions at infinity. However, the strain fields are well-defined, and the surface displacements can be approximated (at least qualitatively) by Equation (20), integrating a know strain field through the depth. Figure 7a-c depicts the distributions of vertical strain zz with depth beneath three points of the loaded surface region. zz can be calculated from the stress components (Equations (8)-(11)) using Hooke's law. The strain values are dependent on the Lame's parameters λ, µ; these are chosen so that that corresponding Young's modulus is E = 98, 000 kPa, and the Poisson's ratio ν = 0.4. The peak strain location tends to decrease with increasing load at the center (x = 0), while it increases at the midpoint from the edge (x = 0.5a) and at the edge location (x = a). Furthermore, there exists a small region of elastic vertical tension beneath the edge of the load. This agrees with the three-dimensional elastic calculations of [16], and can be viewed as a necessary outcome of Hooke's law and the zero-stress boundary conditions at this point. This phenomenon is therefore an outcome of the continuum-based approximation and must be considered non-physical if dealing with a granular soil.  Figure 7d shows approximate displacement points obtained by the numerical integration of the strain curves of Figure 7a-c. This provides a qualitative understanding of the evolution of stiffness across the width of the loaded footing predicted by the present elastoplastic model. This corresponds predictably with the evolution of contact traction shown in Figure 6. The bulk of the contact traction converges towards the center of the footing with increasing load, which corresponds to a softening of the stiffness curve at x = 0. The contact stiffens at the footing's edge (x = a) as the load becomes less concentrated in the neighborhood of this point, while the stiffness path halfway between the center and edge remains more or less linear. Importantly, non-uniform settlement develops, from which the vertical soil-foundation stiffness may affect bending moment distributions across the width of a strip footing even when relative length b/a > 10.
The strain curves and displacements of Figure 7 are compared with a simplified Schmertmann strain influence analysis [10,11], in particular the type outlined in [39] for plane strain problems. For a critical analysis of Schmertmann's method and variations, see [40]. This method estimates the settlement as: where I z in this case is given as a bilinear curve ranging linearly from 0.2 at z/2a = 0 to 0.6 at z/2a = 0.5, and then linearly to zero at z/2a = 4; p is the average applied overpressure, and E the Young's modulus, assumed constant in this case. The summation is taken over even intervals so that ∆z = 8a/n. Note the difference of the characteristics of Schmertmann's influence factor and the strain curve under load P ult at x = 0 in Figure 7a.
The Schmertmann strain influence curve takes its peak at a deeper location beneath the load, and by assumption varies linearly to zero at z/2a = 4; the calculated strain curves from elasticity have shallower peak locations and still take relatively large nonzero values at z/2a = 4, as they go to zero in the limit at infinity. Note that Schmertmann's method is included in Figure 7 only for the sake of comparison and to highlight these differences.

Example 2: Empirically Measured Contact Traction Fields from Real Granular Materials
A wide variety of experimental literature exists regarding contact traction fields beneath model or real footings founded on granular soils, e.g., [32,36,[41][42][43][44]. In general, measurements of the frictional tangential contact traction fields are more difficult and rarely appear alongside the more easily-obtained normal distributions. Bauer et al. [32] and Muhs and Bub [43] provide measurements of both normal and tangential contact traction for various load magnitudes. Of these, Bauer's experimental apparatus was specifically constructed to simulate plane strain conditions. For this reason, contact traction data from one of Bauer's experimental cases is selected for an example study of the elastostatic stress equilibrium state beneath a strip footing.
The normal and tangential traction fields from the case selected from Bauer's data are shown in Figure 8. It is assumed that both traction components take zero values at the footing's edge; this is physically justified since the edge normal resistance of a surface footing must be zero, while the tangential value is bounded by the normal under a finite coefficient of friction. The assumption of zero edge values also allows for the removal of a jump discontinuity in the boundary conditions which is known to lead to singularities in the resulting stress tensor [8]. The distributions in Figure 8 are each the result of a curve fit of the resulting eight data points with a polynomial of the form The curve-fit coefficients are listed in Table 1. The traction fields are normalized by the average normal pressure p avg , obtained by the integration of the polynomial resulting from the curve-fit of the normal traction data. The lack of perfect symmetry in the traction data may stem from a combination of eccentricity in the load and inhomogeneity in the soil samples in [32].  [32]. Data points from test 1, load #2 of [32] curve-fit using 7th order polynomial functions to provide continuous expressions for boundary conditions. Table 1. Curve-fit polynomial coefficients for Equation (27) which generate the traction distributions in Figure 8.  Figure 9 shows the vertical stress (σ zz ) distributions at various depths beneath the loaded boundary region. Figure 9a,b show the components of σ zz due to the normal and tangential components of the boundary conditions, respectively. Note that the lack of symmetry in the boundary conditions in rapidly smoothed out with depth within the elastic media; this can be viewed as a necessary outcome of Saint Venant's principle [45]. It is clear that in this case σ zz is dominated by the normal traction component, as one might expect. The total combined values of σ zz differ negligibly from those under the influence of p(x ) alone, as shown in Figure 9c. This is in contrast to the results of [17,18], where the assumed radial tangential contact traction was large enough to have a considerable qualitative and quantitative effects on the internal vertical stress distributions.  Figures 10 and 11 show the same distributions for the horizontal and shear stresses σ xx and σ xz , respectively. In both the cases the stress components are more noticeably influenced by the presence of the surface frictional tangential traction q. This component of contact traction will therefore influence the overall behavior of the principle stress fields in the present case. This suggests that the modes of distribution of the shearing traction component may be non-negligible when assessing the equilibrium states of shallow foundation systems. The preceding stress tensor calculations can be combined to form the following two quantities: the hydrostatic pressure P = 1 3 (σ xx + σ yy + σ zz ), (28) and the von Mises stress Recall that the presence of σ yy in quantities (28) and (29) implies a dependence on the Poisson ratio ν. The hydrostatic pressure and von Mises stress are plotted in Figure 12 for a value of ν = 0.4. Finally, Bauer's paper appears unique in that it provides load-settlement data alongside the contact traction measurements. This allows for a calibration of the strain-influence method of settlement estimation used in Example 1 via Equation (20). For reference, the footing in Bauer's experiment had total width of 0.305 m; the traction values for our example shown in Figure 8 have been normalized by a factor of p avg = 119.09 kPa, calculated by integrating the polynomial fit of Bauer's data using Equation (27). The soil sample was a crushed quartz sand with known density (15.7 kN/m 3 ) and internal friction angle (45 degrees). Bauer reported a settlement of approximately 1 cm for this loading case (test 1, load 2). Figure 13 shows the strain under the center of the footing estimated by the present method with a fixed Poisson ratio of 0.4. The strain values are multiplied by the Young's modulus, which is unknown in this case. Note that the presence of the shearing traction has negligible effect on the normal strain (and thus the soil-structure stiffness): this is in good agreement with the results of [17,18]. Integrating the strain curves from the combined loading yields a value of wE = 48.798 kPa-m; thus, the model requires a Young's modulus of around 5000 kPa to match the measured settlement value.

Discussion and Conclusions
Constant and linear analytic elements have been developed for the Neumann boundary value problem of a homogeneous and isotropic elastic half-plane under plane strain conditions. These are derived as the spatial limit of similar three-dimensional solutions over rectangular areas. A simple algorithm for superposition of solutions corresponding to linear interpolants of boundary conditions provides a convenient method for the analysis of general contact problems prescribing arbitrary distributions of normal and tangential surface traction.
Examples related to soil-structure interaction were selected to highlight the application of the proposed method. The elastoplastic normal tractions from Example 1 were selected in part to highlight the fact that contact is almost always a nonlinear phenomenon; evolving boundary conditions will introduce nonlinearity to an otherwise linear problem. It is worth noting that some error in this case is likely, due to the application of boundary conditions derived in part from plasticity in a purely elastic domain. However, based on Figure 7, it is clear that some of the expected nonlinearity in the load-displacement behavior is captured by the present method, while absent from the Schmertmann example. Interestingly, despite the discussed physical differences between the strain curves in the two methods, the present method yields settlement estimates at the center of the load which are strikingly close to the Schmertmann values. The similarity in the estimated settlements suggests that analysis can be improved by employing well-derived or measured boundary conditions, even when using a simple linear elastic constitutive model.
The selection of the empirically-measured normal and tangential tractions in Example 2 were motivated by recent work that suggests that the presence of the tangential contact traction fields may have a non-negligible effect on the analysis [17,18]. In particular, the tangential traction field should be included in a contact problem with friction if accuracy in the local stress tensor is desired. The data from [32] was chosen because it is one of few places where both normal and tangential fields were recorded and published for a problem approximating plane strain conditions. Far more research is needed before the normal and tangential contact traction fields can be accurately predicted a priori for a given soil-structure system, as they are dependent on the load magnitude and geometry and a plethora of uncertain multiscale properties of the granular material. The uncertainty and heterogeneity related to the evolution of contact traction fields generated between structural footings and underlying soils may make the present method attractive to geotechnical engineers and analysts wishing to assess the elastic settlement and equilibrium states of shallow foundation systems under varying assumptions.
One major limitation of the application of this model to soil-structure interaction is the problem of scaling: there is a known scale-dependency of soil-foundation vertical stiffness. The present examples have normalized the boundary conditions by footing size and load magnitude, while some of the real scale-dependency stems from the interactions with a complex granular body. Further work may try to calibrate the model to site-scale settlement measurements (e.g., [46]), but the true utility of this model can only be assessed if load-settlement data is presented alongside the measured contact traction fields. The argument here is that some of the nonlinearity from the scale and load dependence can be captured by accurate boundary conditions.
There is further uncertainty introduced in the present model when assessing the material constants. Clearly, there are natural limitations when using any continuum method to describe the mechanics of granular materials. The uncertainty increases when using the simplest possible constitutive model with two elastic constants. The effects of varying the Poisson ratio in the present model is non-negligible when exact results are required, but small when compared with uncertainty in the boundary conditions, influence depth, and Young's modulus; thus, the value of 0.4 was deemed reasonable for use throughout this entire paper.
The value of Young's modulus used in Example 1 (i.e., 98,000 kPa) was deemed a decent estimate for a dense homogeneous sand by assessing values listed in a classic geotechnical reference [47]. However, the value required to match the settlement measurement of Bauer in Example 2 was an order of magnitude smaller (approximately 5000 kPa). This value also seems reasonable when when compared to those used by Schmertmann when estimating the modulus from cone penetration test data [10]. Even a "homogeneous" soil will have variations in the modulus from this estimate with depth. This further exacerbates the issue of sensitivity to strain influence depth and the problem with scaling. The Schmertmann method is designed to account for this material heterogeneity, while the present solutions are for a truly homogeneous material sans body forces. Applying the strain predictions from this model over layers of varying elastic moduli has not been tested or verified against any known data. One approximate workaround is to do the analysis based on in-situ measurements and then take an average to assess the effective Young's modulus for the homogeneous problem. Ultimately, the "elastic modulus" is not a real attribute of a granular body, but instead a parameter of idealized models which may take different values in different contexts. The value determined from triaxial test data may not be the same as the value correlated from a cone penetration test (or that which is required to match a measured settlement when analyzing the soil-structure contact problem). The analyst must be ever-cognizant of this fact.
The solutions presented here are analytic solutions to a type of boundary value problem that has historically been applied to geotechnical problems, but the assumptions in the model are vast over-simplifications of any real soil-structure interaction problem. The purpose of this work was to provide incremental progress by generalizing the types of boundary conditions that can easily be assessed. However, the solutions and method are applicable to more general non-Hertzian contacts and may also be of use to practitioners in the field of tribology. They can be applied to assess the internal stress and strain field for any contact situation in which the plane strain assumption is valid and when surface traction fields can be measured or deduced a priori.

Conflicts of Interest:
The authors declare no conflict of interest.