APPENDIX A. FORMULATION OF LAYERED ELASTIC ANALYSIS EQUATIONS

The compatibility equation for linear elastic behavior of a continuous
body in cylindrical coordinates is:

$\left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial}{\partial z^{2}} \right)\left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial}{\partial z^{2}} \right)\phi(r,z) = 0$
(A-1)

$\text{or, }\nabla^{2}\nabla^{2}\varphi(r,z) = 0$ (A-2)

where:

> *z* = vertical axis, positive downwards
>
> *r* = radius from the vertical axis in the horizontal plane
>
>     (stresses and strains are constant in the horizontal plane for any
> *r* and *z*)

$f(r,z)$is a stress function satisfying the equilibrium equations such
that:

$\sigma_{z} = \frac{\partial}{\partial z}\left( (2 - \upsilon)\nabla^{2}\varphi(r,z) - \frac{\partial ^{2}\varphi(r,z)}{\partial z^{2}} \right)$
(A-3)

$\sigma_{r} = \frac{\partial}{\partial z}\left( \upsilon\nabla^{2}\varphi(r,z) - \frac{\partial ^{2}\varphi(r,z)}{\partial r^{2}} \right)$
(A-4)

$\sigma_{t} = \frac{\partial}{\partial z}\left( \upsilon\nabla^{2}\varphi(r,z) - \frac{1}{r}\frac{\partial\varphi(r,z)}{\partial r} \right)$
(A-5)

$\tau_{\text{rz}} = \frac{\partial}{\partial r}\left( (1 - \upsilon)\nabla^{2}\varphi(r,z) - \frac{\partial ^{2}\varphi(r,z)}{\partial z^{2}} \right)$
(A-6)

$w = \frac{1 + \upsilon}{E}\left( 2(1 - \upsilon)\nabla^{2}\varphi(r,z) - \frac{\partial ^{2}\varphi(r,z)}{\partial z^{2}} \right)$
(A-7)

$u = - \ \left( \frac{1 + \upsilon}{E} \right)\ \frac{\partial ^{2}\varphi(r,z)}{\partial r\ \partial z}$
(A-8)

where:

> *σ~z~* = vertical stress
>
> *σ~r~* = radial stress in the horizontal plane
>
> *σ~t~* = tangential stress in the horizontal plane (perpendicular to
> *σ~r~*)
>
> *τ~rz~* = shear stress
>
> *w* = vertical displacement
>
> *u* = radial displacement
>
> *E* = Young's modulus
>
> *υ* = Poisson's ratio

When the stress function satisfies the compatibility equation the
following conditions hold:

-   displacement throughout the body is continuous;

-   each elemental part of the body is in equilibrium; and

-   stress and strain are related by Hooke's law.

The equation given is correct for constant internal body forces, such as
gravity, but not for cases where the body forces vary with *r* or *z*.

The problem now is to find a stress function which satisfies the
compatibility equation. After the stress function has been found, it can
be substituted into the stress and displacement equations to give
numerical solutions.

To find the stress function, start by taking the zero order Hankel
transform of the compatibility equation. This reduces the compatibility
equation to an equation algebraic in the transform variable (*α*) and
with complete differentials in *z*. The zero-order Hankel transform pair
is:

$\begin{array}{r}
f_{H}(\alpha\ ) = \int_{}^{}{F(r)\text{rJ}_{0}(\alpha r)\ \text{dr}} \\
F(r) = \int_{}^{}{f_{H}(\alpha\ ){\alpha J}_{0}(\alpha r)\ d\alpha}
\end{array}$ (A-9)

where $J_{0}$ is the zero order Bessel function.

The Hankel transform of the compatibility equation is:

$\begin{array}{r}
\left( - \alpha^{2} + \frac{d^{2}}{\text{dz}^{2}} \right)\left( - \alpha^{2} + \frac{d^{2}}{\text{dz}^{2}} \right)\ v_{H}(\alpha,z) = 0 \\
\left( - \alpha^{2} + \frac{d^{2}}{\text{dz}^{2}} \right)\left( - \alpha^{2}v_{H}(\alpha,z) + \frac{d^{2}v_{H}(\alpha,z)}{\text{dz}^{2}} \right) = 0 \\
\left( \alpha^{4}v_{H}(\alpha,z) - 2\alpha^{2}\frac{d^{2}v_{H}(\alpha,z)}{\text{dz}^{2}} + \frac{d^{4}v_{H}(\alpha,z)}{\text{dz}^{4}} \right) = 0
\end{array}$ (A-10)

$$\text{where }v_{H}(\alpha\text{,z})\text{ is the Hankel transform of }f(r,z)\text{ with respect to }r\text{.}$$

Performing the indicated operations results in an ordinary differential
equation in *z*:

$\frac{d^{4}v_{H}(\alpha,z)}{\text{dz}^{4}} - 2\alpha^{2}\frac{d^{2}v_{H}(\alpha,z)}{\text{dz}^{2}} + \alpha^{4}v_{H}(\alpha,z) = 0$
(A-11)

Now solve using the Laplace transform method:

$\begin{array}{r}
s^{4}v_{H}(\alpha,s) - s^{3}v_{H}(\alpha,0) - s^{2}v_{H}^{'}(\alpha,0) - sv_{H}^{''}(\alpha,0) - v_{H}^{'''}(\alpha,0) \\
 - s^{2}2\alpha^{2}v_{H}(\alpha,s) + s2\alpha^{2}v_{H}(\alpha,0) + 2\alpha^{2}v_{H}^{'}(\alpha,0) \\
 + \alpha^{4}v_{H}(\alpha,s) = 0
\end{array}$ (A-12)

where $v_{H}(\alpha,0),\ v_{H}^{'}(\alpha,0)$, etc., are initial
conditions of
$\ v_{H}(\alpha,z),\ \frac{\text{dv}_{H}(\alpha,z)}{\text{dz}}$, etc.

Therefore:

$(s^{4} - 2\alpha^{2}s^{2} + \alpha^{4})\ v_{H}(\alpha,s) = A^{}s^{3} + B^{}s^{2} + C^{}s + D^{}$
(A-13)

$\begin{array}{r}
v_{H}(\alpha,0) = A^{} \\
v_{H}^{'}(\alpha,0) = B^{} \\
v_{H}^{''}(\alpha,0) = C^{} + 2\alpha^{2}A^{} \\
v_{H}^{'''}(\alpha,0) = D^{} + 2\alpha^{2}B^{}
\end{array}$ (A-14)

and

$\begin{array}{r}
\ v_{H}(\alpha,s)\  = \frac{A^{}s^{3} + B^{}s^{2} + C^{}s + D^{}}{s^{4} - 2\alpha^{2}s^{2} + \alpha^{4}} \\
\frac{A^{}s^{3} + B^{}s^{2} + C^{}s + D^{}}{(s^{2} - \alpha^{2})^{2}}
\end{array}$ (A-15)

Solving for each term on the right hand side, and making use of
mathematical identities (for example):

$\frac{s^{3}}{(s^{2} - \alpha ^{2})^{2}} = \frac{1}{2(s - \alpha )} + \frac{1}{2(s + \alpha )} + \frac{\alpha}{4(s - \alpha )^{2}} - \frac{\alpha}{4(s + \alpha )^{2}}$

and

$\begin{array}{r}
F(\alpha,z) = \frac{1}{2}(e^{\alpha z} + e^{- \alpha z}) + \frac{1}{4}\alpha z(e^{\alpha z} - e^{- \alpha z}) \\
\text{cosh}(\alpha z) + \frac{1}{2}\alpha z\text{sinh}(\alpha z)
\end{array}$

gives:

$\begin{array}{r}
v_{H}(\alpha,z) = \frac{A^{}}{2}(e^{\alpha z} + e^{- \alpha z}) + \frac{A^{}\alpha z}{4}(e^{\alpha z} - e^{- \alpha z}) - \frac{B^{}}{4\alpha }(e^{\alpha z} + e^{- \alpha z}) + \frac{B^{}z}{4}(e^{\alpha z} + e^{- \alpha z}) \\
 + \frac{C^{}z}{4}(e^{\alpha z} - e^{- \alpha z}) - \frac{D^{}}{4\alpha }(e^{\alpha z} - e^{- \alpha z}) + \frac{D^{}z}{4}(e^{\alpha z} - e^{- \alpha z})
\end{array}$ (A-16)

Or, expressed in hyperbolic functions:

$\begin{array}{r}
v_{H}(\alpha,z) = A^{}\text{cosh}(\alpha z) + \frac{A^{}\alpha z}{2}\text{sinh}(\alpha z) \\
 - \frac{B^{}}{2\alpha }\text{cosh}(\alpha z) + \frac{(B^{} + D^{})z}{2}\text{cosh}(\alpha z) \\
 - \frac{D^{}}{2\alpha }\text{sinh}(\alpha z) + \frac{C^{}z}{2}\text{sinh}(\alpha z)
\end{array}$ (A-17)

However, since the starred constants are undetermined functions of *α*,
equation (A-17) can be simplified to:

$v_{H}(\alpha,z) = {Ae}^{\alpha z} + {Be}^{- \alpha z} + \text{Cz}e^{\alpha z} + \text{Dz}e^{- \alpha z}$
(A-18)

The stress function can now be found by inverse Hankel transformation:

$\begin{array}{r}
\varphi\ (r,z) = \int_{}^{}{v_{H}(\alpha,z){\alpha J}_{0}(\alpha\ r)\ d\alpha} \\
\int_{}^{}{({Ae}^{\alpha z} + {Be}^{- \alpha z} + \text{Cz}e^{\alpha z} + \text{Dz}e^{- \alpha z}){\alpha J}_{0}(\alpha\ r)\ d\alpha}
\end{array}$ (A-19)

Equation (A-19) does not have an analytic solution and must be solved
numerically for given values of *α* and *z*. The constants *A*, *B*,
*C*, and *D* are found from the boundary conditions. The constants are
functions of *α* and new values must be found at each step of the
integration, which means that the complete stress function cannot be
used to find the constants. However, the compatibility equation
operations are all performed within the integral sign of the stress
function and each elemental part of the integration is, itself, a stress
function satisfying the compatibility equation. The integrand can
therefore be substituted into the stress and displacement equations to
find the constants at each value of *α*. Substituting the integrand into
the equations for vertical and shear stresses and for vertical and
radial displacements gives the "elemental" stresses and displacements:

$\begin{array}{r}
\sigma_{z}^{} = \alpha^{3}J_{0}(\alpha r)\left\{ \left\lbrack - A\alpha + C(1 - 2\upsilon - \alpha z) \right\rbrack e^{\alpha z} + \left\lbrack B\alpha + D(1 - 2\upsilon + \alpha z) \right\rbrack e^{- \alpha z} \right\} \\
\tau_{\text{rz}}^{} = \alpha^{3}J_{1}(\alpha r)\left\{ \left\lbrack A\alpha + C(2\upsilon\  + \alpha\ z) \right\rbrack e^{\alpha z} + \left\lbrack B\alpha - D(2\upsilon\  - \alpha\ z) \right\rbrack e^{- \alpha z} \right\} \\
w^{} = \frac{1 + \upsilon}{E}\alpha^{2}J_{0}(\alpha r)\left\{ \left\lbrack - A\alpha + C(2 - 4\upsilon - \alpha z) \right\rbrack e^{\alpha z} - \left\lbrack B\alpha + D(2 - 4\upsilon + \alpha z) \right\rbrack e^{- \alpha z} \right\} \\
u^{} = \frac{1 + \upsilon}{E}\alpha^{2}J_{1}(\alpha r)\left\{ \left\lbrack A\alpha + C(1 + \alpha z) \right\rbrack e^{\alpha z} - \left\lbrack B\alpha - D(1 - \alpha z) \right\rbrack e^{- \alpha z} \right\}
\end{array}$ (A-20)

And the complete stresses and displacements, written to show the Hankel
form, are:

$\begin{array}{r}
\sigma_{z} = \int_{}^{}{\alpha^{2}\left\{ \left\lbrack - A\alpha + C(1 - 2\upsilon - \alpha z) \right\rbrack e^{\alpha z} + \left\lbrack B\alpha + D(1 - 2\upsilon + \alpha z) \right\rbrack e^{- \alpha z} \right\}{\alpha J}_{0}(\alpha r)d\alpha} \\
\tau_{\text{rz}} = \int_{}^{}{\frac{\alpha^{2}J_{1}(\alpha r)}{J_{0}(\alpha r)}\left\{ \left\lbrack A\alpha + C(2\upsilon\  + \alpha\ z) \right\rbrack e^{\alpha z} + \left\lbrack B\alpha - D(2\upsilon\  - \alpha\ z) \right\rbrack e^{- \alpha z} \right\}{\alpha J}_{0}(\alpha r)d\alpha} \\
w = \frac{1 + \upsilon}{E}\int_{}^{}{\alpha\left\{ \left\lbrack - A\alpha + C(2 - 4\upsilon - \alpha z) \right\rbrack e^{\alpha z} - \left\lbrack B\alpha + D(2 - 4\upsilon + \alpha z) \right\rbrack e^{- \alpha z} \right\}{\alpha J}_{0}(\alpha r)d\alpha} \\
u = \frac{1 + \upsilon}{E}\int_{}^{}{\frac{{\alpha J}_{1}(\alpha r)}{J_{0}(\alpha r)}\left\{ \left\lbrack A\alpha + C(1 + \alpha z) \right\rbrack e^{\alpha z} - \left\lbrack B\alpha - D(1 - \alpha z) \right\rbrack e^{- \alpha z} \right\}{\alpha J}_{0}(\alpha r)d\alpha}
\end{array}$ (A-21)

Considering the top layer, which is a special case having a prescribed
vertical circular load, *p*(*r*), and zero shear stress at the top
surface, the complete stresses when *z* = 0 are:

$\begin{array}{r}
\int_{}^{}{\alpha^{2}\left\{ - A\alpha + C(1 - 2\upsilon) + B\alpha + D(1 - 2\upsilon) \right\}{\alpha J}_{0}(\alpha r)d\alpha} = p(r) \\
\int_{}^{}{\frac{\alpha^{2}J_{1}(\alpha r)}{J_{0}(\alpha r)}\left\{ A\alpha + 2C\upsilon + B\alpha - 2D\upsilon\  \right\}{\alpha J}_{0}(\alpha r)d\alpha} = 0
\end{array}$ (A-22)

The Hankel transform of a circular load at the origin of intensity *q*
and radius *a* is (for tension positive):

$H(p(r)) = - \frac{\text{qa}}{\alpha}J_{1}(\alpha\ a)$ (A-23)

Consequently, the Hankel transforms of the vertical and shear stress
equations at *z* = 0 are:

$\begin{array}{r}
\alpha^{2}( - A\alpha + C(1 - 2\upsilon) + B\alpha + D(1 - 2\upsilon)) = \frac{- \text{qa}}{\alpha}J_{1}(\alpha\ a) \\
\frac{\alpha^{2}J_{1}(\alpha r)}{J_{0}(\alpha r)}(A\alpha + 2C\upsilon + B\alpha - 2D\upsilon) = 0
\end{array}$ (A-24)

and:

$\begin{array}{r}
 - A\alpha + C(1 - 2\upsilon) + B\alpha + D(1 - 2\upsilon) = \frac{- \text{qa}}{\alpha^{3}}J_{1}(\alpha\ a) \\
A\alpha + 2C\upsilon + B - 2D\text{αυ} = 0
\end{array}$ (A-25)

The constants for the top layer are therefore expressed as functions of
$P(\alpha\ a) = \frac{\text{qa}}{\alpha^{3}}J_{1}(\alpha\ a)$, where two
equations can be written in four unknowns (the constants). Writing out
the stress and displacement equations at either side of each of the
remaining interfaces leads to a series of four equations in eight
unknowns (the constants for the two layers) which satisfy the boundary
conditions. Satisfying the boundary conditions at each interface relates
the constants for each layer to those of the layer above. All constants
are therefore functions of $P(\alpha\ a)$. By writing out the equations
in matrix form, and synthesizing Gaussian elimination, it can be shown
that each constant is a simple multiple of $P(\alpha\ a)$. Therefore,
$P(\alpha\ a)$ can, if desired be moved into the elemental stress
equations and the first element of the right-hand-side vector be set to
a value of one.

Boundary conditions for the bottom of the bottom layer require special
treatment in the same way as the top of the top layer. If the bottom
layer is of infinite extent, two of the constants can be found by
multiplying the Hankel transforms of the vertical and shear stress
equations through by $e^{- \alpha z}$, and rearranging:

$\begin{array}{r}
\frac{H(\sigma_{z})e^{- \alpha z}}{\alpha^{2}} = \left\{ - A\alpha + C(1 - 2\upsilon - \alpha z) + \left\lbrack B\alpha + D(1 - 2\upsilon + \alpha z) \right\rbrack e^{- 2\alpha z} \right\} \\
\frac{H(\tau_{z})J_{0}(\alpha r)e^{- \alpha z}}{\alpha^{2}J_{1}(\alpha r)} = A\alpha + C(2\upsilon\  + \alpha\ z) + \left\lbrack B\alpha - D(2\upsilon\  - \alpha\ z) \right\rbrack e^{- 2\alpha z}
\end{array}$ (A-26)

For *z = infinity*:

$\begin{array}{r}
0 = - A\alpha + C - 2C\upsilon - C\alpha z \\
0 = A\alpha + 2C\upsilon + C\alpha z
\end{array}$ (A-27)

Adding and subtracting the two equations (A-27) gives *C =* 0 and *A =*
0.

Stresses and displacements therefore decay as $e^{- \alpha z}$ and
${ze}^{- \alpha z}$ as *z* increases downward from the subgrade
interface.

For a rigid constraint at the bottom of the structure, the two
displacement equations must be used to find the two equations in four
unknowns.

The form of the boundary condition equations is shown below for a
three-layer structure.

$\begin{pmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & & & & & & & & \\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & & & & & & & & \\
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} & a_{3,6} & a_{3,7} & a_{3,8} & & & & \\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} & a_{4,6} & a_{4,7} & a_{4,8} & & & & \\
a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} & a_{5,6} & a_{5,7} & a_{5,8} & & & & \\
a_{6,1} & a_{6,2} & a_{6,3} & a_{6,4} & a_{6,5} & a_{6,6} & a_{6,7} & a_{6,8} & & & & \\
 & & & & a_{7,5} & a_{7,6} & a_{7,7} & a_{7,8} & a_{7,9} & a_{7,\text{10}} & a_{7,\text{11}} & a_{7,\text{12}} \\
 & & & & a_{8,5} & a_{8,6} & a_{8,7} & a_{8,8} & a_{8,9} & a_{8,\text{10}} & a_{8,\text{11}} & a_{8,\text{12}} \\
 & & & & a_{9,5} & a_{9,6} & a_{9,7} & a_{9,8} & a_{9,9} & a_{9,\text{10}} & a_{9,\text{11}} & a_{9,\text{12}} \\
 & & & & a_{\text{10},5} & a_{\text{10},6} & a_{\text{10},7} & a_{\text{10},8} & a_{\text{10},9} & a_{\text{10},\text{10}} & a_{\text{10},\text{11}} & a_{\text{10},\text{12}} \\
 & & & & & & & & a_{\text{11},9} & a_{\text{11},\text{10}} & a_{\text{11},\text{11}} & a_{\text{11},\text{12}} \\
 & & & & & & & & a_{\text{12},9} & a_{\text{12},\text{10}} & a_{\text{12},\text{11}} & a_{\text{12},\text{12}}
\end{pmatrix}\ \begin{pmatrix}
A_{1} \\
B_{1} \\
C_{1} \\
D_{1} \\
A_{2} \\
B_{2} \\
C_{2} \\
D_{2} \\
A_{3} \\
B_{3} \\
C_{3} \\
D_{3}
\end{pmatrix} = \begin{pmatrix}
P(\alpha\ a) \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix}$

(A-28)

and

$P(\alpha\ a) = \frac{\text{qa}}{\alpha^{3}}J_{1}(\alpha\ a)$ (A-29)

When the bottom layer is of infinite extent, *A~3~* and *C~3~* are zero
and the corresponding rows and columns can be removed from the
coefficient matrix. The number of equations for this case is therefore
4 x 3 - 2 = 10.

$\begin{pmatrix}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & & & & & & \\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & & & & & & \\
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & a_{3,5} & a_{3,6} & a_{3,7} & a_{3,8} & & \\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4} & a_{4,5} & a_{4,6} & a_{4,7} & a_{4,8} & & \\
a_{5,1} & a_{5,2} & a_{5,3} & a_{5,4} & a_{5,5} & a_{5,6} & a_{5,7} & a_{5,8} & & \\
a_{6,1} & a_{6,2} & a_{6,3} & a_{6,4} & a_{6,5} & a_{6,6} & a_{6,7} & a_{6,8} & & \\
 & & & & a_{7,5} & a_{7,6} & a_{7,7} & a_{7,8} & a_{7,9} & a_{7,\text{10}} \\
 & & & & a_{8,5} & a_{8,6} & a_{8,7} & a_{8,8} & a_{8,9} & a_{8,\text{10}} \\
 & & & & a_{9,5} & a_{9,6} & a_{9,7} & a_{9,8} & a_{9,9} & a_{9,\text{10}} \\
 & & & & a_{\text{10},5} & a_{\text{10},6} & a_{\text{10},7} & a_{\text{10},8} & a_{\text{10},9} & a_{\text{10},\text{10}}
\end{pmatrix}\ \begin{pmatrix}
A_{1} \\
B_{1} \\
C_{1} \\
D_{1} \\
A_{2} \\
B_{2} \\
C_{2} \\
D_{2} \\
B_{3} \\
D_{3}
\end{pmatrix} = \begin{pmatrix}
P(\alpha\ a) \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix}$ (A-30)

# The Boussinesq Case.

The Boussinesq case is for an infinite half space with a circular load
at the origin. The solution for vertical stress when radius is zero is
(see reference \[A-1\], page 368):

$\sigma_{z} = q\left( \frac{z^{3}}{\left( z^{2} + a^{2} \right)^{32}} - 1 \right)$
(A-31)

The surface boundary condition for a circular load at the surface is:

$- A + C(1 - 2\upsilon\ ) + B + D(1 - 2\upsilon\ ) = - \frac{\text{qa}}{\alpha^{3}}J_{1}(\alpha\ a)$
(A-32)

The surface boundary condition for zero shear stress at the surface is:

$A + 2C\upsilon + B - 2D\upsilon = 0$ (A-33)

Vertical stress is:

$\sigma_{z} = \int_{}^{}{\alpha^{2}\left\lbrack ( - A + C(1 - 2\upsilon))e^{\alpha z} - {C\alpha ze}^{\alpha z} + (B + D(1 - 2\upsilon))e^{- \alpha z} + {D\alpha ze}^{- \alpha z} \right\rbrack\ }{\alpha J}_{0}(\alpha\ r)d\alpha$
(A-34)

For a solution when *z* = ∞, *A* = *C* = 0, solving the boundary
condition equations for *B* and *D* gives:

$B = - 2\upsilon\frac{\text{qa}}{\alpha^{3}}J_{1}(\alpha a)\text{  and  }D = - \frac{\text{qa}}{\alpha^{3}}J_{1}(\alpha a)$
(A-35)

Substitution into the vertical stress equation (A-34) and setting *r* =
0 gives:

$\sigma_{z} = - \text{qa}\int_{}^{}{e^{- \alpha z}J_{1}(\alpha a)d\alpha} - \text{qa}\int_{}^{}{{\alpha ze}^{- \alpha z}J_{1}(\alpha a)d\alpha}$
(A-36)

Using relations 6.611(1.) and 6.623(1.) from reference \[A-2\]:

$\sigma_{z} = - \text{qa}\left( \frac{\left( z^{2} + a^{2} \right)^{12} - z}{a\left( z^{2} + a^{2} \right)^{12}} + \frac{a}{\left( z^{2} + a^{2} \right)^{32}} \right) = - \text{qa}\left( \frac{\left( z^{2} + a^{2} \right)^{32} - z^{3} - a^{2}z}{a\left( z^{2} + a^{2} \right)^{32}} + \frac{\text{az}}{\left( z^{2} + a^{2} \right)^{32}} \right)$
(A-37)

and, rearranging:

$s_{z}\text{=-}q\left( 1 - \frac{z^{3}}{\left( z^{2} + a^{2} \right)^{32}} \right)$
(A-38)

which agrees with the Boussinesq solution (A-31).

[References]{.underline}.

A-1. Timoshenko, S., and Goodier, J.N., Theory of Elasticity, Second
Edition, McGraw-Hill, 1951.

A-2. Gradshteyn, I.S., and Ryzhik, I.M., Table of Integrals, Series, and
Products, Academic Press, 1980