5.  LEAF - COMPUTER PROGRAM FOR LAYERED ELASTIC ANALYSIS

# 5.1  Introduction.

A complete derivation of the basic layered elastic equations is given in
Appendix A and reference \[Huang, Y.H., Pavement Analysis and Design,
Englewood Cliffs, New Jersey, Prentice-Hall, 1993.\] and only a brief
introduction is given here. Cylindrical coordinates are used with no
response variation in the circumferential direction. The coordinates are
*z* downwards, *r* radially outward from the vertical *z* axis, and *t*
tangent to a circumference at radius *r*. Stresses, strains, and
displacements within a given layer are computed by integrals of the
form:

$\sigma_{z} = - \text{qa}\int_{}^{}{J_{0}(\alpha r)J_{1}(\alpha a)}\begin{array}{r}
\left( \left\{ A(\alpha) - C(\alpha)(1 - 2\upsilon - \alpha z) \right\} e^{\alpha z} \middle| \right) \\

\end{array}$ (5-1)

Vertical stress, *σ~z~*, is shown for illustration. Load is applied at
the top surface of the top layer in a multi-layer system. The load is a
uniform compressive traction, of pressure *q*, applied over a circular
area of radius *a*. The stress computation point is at a distance *z*
from the top of the layer and a distance *r* from the center of the
vertical axis passing through the center of the load. *J*~0~ and *J*~1~
are Bessel functions and *α* is the integration variable (strictly the
Hankel domain variable). In the following, the integrands are referred
to as "elemental" stresses and displacements to signify that the
complete response is composed of the summation of a large number of
small, elemental, components.

FIGURE 5-1.  NOMENCLATURE AND COORDINATES

The coefficients *A*(*α*), *B*(*α*), *C*(*α*), and *D*(*α*) are
determined by the boundary conditions for the layer (usually specified
at the top and bottom surfaces of the layer). They are functions of *α*,
as indicated by the parentheses. In the following the *α* and the
parentheses are dropped, but it should be remembered that the
coefficients are functions of *α* and must be recomputed whenever *α*
changes. Also, because the coefficients are functions of *α*, the terms
in braces in equation 5-1 can be multiplied by any function of *α*,
provided only that the multiplier is consistently applied to all
integrands used to set up the boundary conditions.

Derivation of the equilibrium equations for a multi-layered system with
fully bonded interfaces is described first, followed by derivation of
the equilibrium equations for a system with partially bonded interfaces.
Computer implementation of the integrations to compute system responses
is then described. Methods for solving the equations are described last.
Numerical test cases and comparisons with other methods complete the
chapter.

# 5.2  Derivation of the Interface Equilibrium Equations

The $\text{qaJ}_{1}(\alpha a)$ term represents the applied circular load
at the top surface and always cancels when formulating the equilibrium
equations. Therefore, from Appendix A, and ignoring the load term, the
complete set of elemental stresses and displacements is:

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

Vertical, radial, and tangential stresses are denoted by *σ~z~*, *σ~r~*,
and *σ~t~* respectively.

Shear stress in the vertical-radial plane is denoted by *τ~zr~*.

Vertical and radial displacements are denoted by *w* and *u*
respectively.

Layers are numbered downward starting at 1 for the top layer. Interface
*i* is the interface at the bottom of layer *i*. See figure 5-2.

FIGURE 5-2.  LAYER AND INTERFACE NUMBERING

For a circular load and zero shear stress at the surface, the
equilibrium equations for the surface plane (interface 0) are:

$$\begin{matrix}
 \\
\sigma_{z}(\text{vertical stress}) \\
\tau_{\text{zr}}(\text{shear stress}) \\
\begin{pmatrix}
 - 1 & 1 & (1 - 2\upsilon_{1}) & (1 - 2\upsilon_{1}) \\
1 & 1 & 2\upsilon_{1}\  & - 2\upsilon_{1}
\end{pmatrix}\begin{pmatrix}
A_{1} \\
B_{1} \\
C_{1} \\
D_{1}
\end{pmatrix} = \begin{pmatrix}
 - 1 \\
0
\end{pmatrix}
\end{matrix}$$

$\begin{pmatrix}
1 & - 1 & - (1 - 2\upsilon_{1}) & - (1 - 2\upsilon_{1}) \\
1 & 1 & 2\upsilon_{1}\  & - 2\upsilon_{1}
\end{pmatrix}\begin{pmatrix}
A_{1} \\
B_{1} \\
C_{1} \\
D_{1}
\end{pmatrix} = \begin{pmatrix}
1 \\
0
\end{pmatrix}$ (5-3)

The minus one (--1) term in the right hand side vector in the first set
of equations represents the compressive stress applied by the circular
load (the stress intensity *qa* appears in the response equations). The
signs are changed to make the load vector positive just for convenience.
The left hand side of the equations is found by setting *z* in layer 1
(the top layer) equal to zero for the vertical and shear stress
elemental equations.

For two layers fully bonded at a common interface, the boundary
conditions are determined by equilibrium and continuity at the
interface. Therefore, at the interface, the stresses and displacements
are equal in the two layers. The equilibrium equations are set up by
equating, in turn, vertical stress, shear stress, vertical displacement,
and radial displacement for the two layers bounding the interface. The
appropriate values of vertical distance, *z*, are substituted into the
elemental equations for two adjacent layers *i* and *i*+1. Interface *i*
is at the bottom of layer *i* and at the top of layer *i*+1. The
distance, *z*, from the top of layer *i* to the interface is the
thickness of layer *i*, *h~i~*. Therefore, set *z* = *h~i~* in the
elemental equations to obtain the left hand side of the equation. The
distance, *z*, from the top of layer *i*+1 to the interface is zero.
Therefore, set *z* = 0 in the elemental equations to obtain the right
hand side of the equation. Also set:

  ------------ -----------------------------------------------------------------
  *E~i~*       = Young's modulus of layer *i*

  *E~i~*~+1~   = Young's modulus of layer *i*+1

  *ν~i~*       = Poisson's ratio of layer *i*

  *ν~i~*~+1~   = Poisson's ratio of layer *i*+1
  ------------ -----------------------------------------------------------------

$$R_{i} = \frac{1 + \upsilon_{i + 1}}{1 + \upsilon_{i}}\frac{E_{i}}{E_{i + 1}}$$

Then:

$\begin{array}{r}
\begin{matrix}
\sigma_{z}(\text{vertical stress}) \\
\tau_{\text{zr}}(\text{shear stress}) \\
w(\text{vertical defl}\text{.}) \\
u(\text{radial defl}\text{.})
\end{matrix}\begin{pmatrix}
e^{{\alpha h}_{i}} & - e^{- {\alpha h}_{i}} & - (1 - 2\upsilon_{i} - {\alpha h}_{i})e^{{\alpha h}_{i}} & - (1 - 2\upsilon_{i} + {\alpha h}_{i})e^{- {\alpha h}_{i}} \\
e^{{\alpha h}_{i}} & e^{- {\alpha h}_{i}} & (2\upsilon_{i}\  + \alpha\ h_{i})e^{{\alpha h}_{i}} & - (2\upsilon\ _{i} - \alpha\ h_{i})e^{- {\alpha h}_{i}} \\
e^{{\alpha h}_{i}} & e^{- {\alpha h}_{i}} & - (2 - 4\upsilon_{i} - {\alpha h}_{i})e^{{\alpha h}_{i}} & (2 - 4\upsilon_{i} + {\alpha h}_{i})e^{- {\alpha h}_{i}} \\
e^{{\alpha h}_{i}} & - e^{- {\alpha h}_{i}} & (1 + {\alpha h}_{i})e^{{\alpha h}_{i}} & (1 - {\alpha h}_{i})e^{- {\alpha h}_{i}}
\end{pmatrix}\begin{pmatrix}
A_{i} \\
B_{i} \\
C_{i} \\
D_{i}
\end{pmatrix} = \\
\begin{pmatrix}
1 & - 1 & - (1 - 2\upsilon_{i + 1}) & - (1 - 2\upsilon_{i + 1}) \\
1 & 1 & 2\upsilon_{i + 1} & - 2\upsilon_{i + 1} \\
R_{i} & R_{i} & - (2 - 4\upsilon_{i + 1})R_{i} & (2 - 4\upsilon_{i + 1})R_{i} \\
R_{i} & - R_{i} & R_{i} & R_{i}
\end{pmatrix}\begin{pmatrix}
A_{i + 1} \\
B_{i + 1} \\
C_{i + 1} \\
D_{i + 1}
\end{pmatrix}
\end{array}$ (5-4)

## 5.2.1  Equilibrium Equations with Arbitrary Interface Bonding.

In order to allow relative horizontal movement between two layers at an
interface, a uniformly distributed shear spring is assumed to connect
the layers. The spring acts in the radial direction and has the
following law:

$\begin{array}{r}
\tau_{i} = k_{i}(u_{i} - u_{i + 1})\text{, or} \\
\tau_{i + 1} = k_{i}(u_{i} - u_{i + 1})\text{ since }\tau_{i} = \tau_{i + 1}
\end{array}$ (5-5)

where:

  ------------- --------------------------------------------------------------
  *τ~i~*        = radial shear stress at the interface between layers *i* and
                *i+*1.

  (*u~i~* --    = relative radial displacement across the interface.
  *u~i~*~+1~)   

  *k~i~*        = interface spring stiffness, with units of lb / inch relative
                radial displacement / inch along the radius / inch along the
                circumference at radius *r* (lb/in^3^). That is, a radial
                spring connects elemental areas either side of the interface.
                The spring resists relative radial displacement across the
                interface.
  ------------- --------------------------------------------------------------

To reduce numerical complications, the computer implementation can use
the relationship:

$\begin{array}{r}
k_{i} = \frac{l_{i}}{1 - l_{i}}\text{, giving} \\
(1 - l_{i})\tau_{i} = l_{i}(u_{i} - u_{i + 1})
\end{array}$ (5-6)

This change of variable is used in JULEA, with a further logarithmic
transformation in the input data parameter. LEAF uses the same variable,
*l*, both within the computational routines and as the input data
parameter.

For fully bonded layers,
$k_{i} = \infty\text{,  }l_{i} = 1\text{, and }u_{i} = u_{i + 1}$.

For fully unbonded layers,
$k_{i} = 0\text{,  }l_{i} = 0\text{, and }\tau_{i} = 0$.

For reference, the transformation in JULEA is as follows:

$m_{i} =$input parameter for interface *i*.

$E_{2} =$modulus of layer 2 (directly below the surface layer).

$$\text{If}\ m_{i} \geq \text{100},\text{000}\ \text{Then}\ l_{i} = 0\ \text{Else}\ l_{i} = \text{10}^{- \ \frac{m_{i}}{E_{2}}}$$

Note that *l* is a parameter. It does not have units and it has no
physical significance other than the manner in which it varies the
spring stiffness. It has a value from 0 to 1 (fully unbonded to fully
bonded).

Writing out the elemental equations for shear stress and radial
deflection in full, a new equilibrium equation for arbitrary interface
bonding is derived as follows.

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

$$\begin{array}{r}
\text{using }(1\text{-l}_{i})\tau_{i} = l_{i}(u_{i} - u_{i + 1}), \\
(1\text{-l}_{i})\tau_{i} - l_{i}u_{i} = - l_{i}u_{i + 1} \\
\text{and setting }z = h_{i}\text{ for layer }i \\
z = 0\text{ for layer }i + 1
\end{array}$$

$$\begin{array}{r}
\ (1 - l_{i})J_{1}(\alpha r)\left\{ \ \left\lbrack A_{i} + C_{i}(2\upsilon_{i}\  + \alpha\ h_{i}) \right\rbrack\ e^{{\alpha h}_{i}} + \left\lbrack B_{i} - D_{i}(2\upsilon_{i}\  - \alpha\ h_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 - \frac{(1 + \upsilon_{i})l_{i}}{E_{i}}\frac{J_{1}(\alpha r)}{\alpha}\left\{ \ \left\lbrack A_{i} + C_{i}(1 + {\alpha h}_{i}) \right\rbrack\ e^{{\alpha h}_{i}} - \left\lbrack \ B_{i} - D_{i}(1 - {\alpha h}_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 =  - \frac{(1 + \upsilon_{i + 1})l_{i}}{E_{i + 1}}\frac{J_{1}(\alpha r)}{\alpha}\left\{ \ A_{i + 1} + C_{i + 1} - B_{i + 1} + D_{i + 1}\  \right\}
\end{array}$$

$$\text{Multiply}\ \text{through by }\frac{{\alpha E}_{i}}{(1 + \upsilon_{i})}\text{ and cancel }J_{1}(\alpha r)$$

$$\begin{array}{r}
\ \frac{{\alpha E}_{i}(1 - l_{i})}{(1 + \upsilon_{i})}\left\{ \ \left\lbrack A_{i} + C_{i}(2\upsilon_{i}\  + \alpha\ h_{i}) \right\rbrack\ e^{{\alpha h}_{i}} + \left\lbrack \ B_{i} - D_{i}(2\upsilon_{i}\  - \alpha\ h_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 - l_{i}\left\{ \ \left\lbrack A_{i} + C_{i}(1 + {\alpha h}_{i}) \right\rbrack\ e^{{\alpha h}_{i}} - \left\lbrack \ B_{i} - D_{i}(1 - {\alpha h}_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 =  - \frac{(1 + \upsilon_{i + 1})E_{i}l_{i}}{(1 + \upsilon_{i})E_{i + 1}}\left\{ \ A_{i + 1} + C_{i + 1} - B_{i + 1} + D_{i + 1}\  \right\}
\end{array}$$

Setting
$G_{i} = \frac{E_{i}(1 - l_{i})}{(1 + \upsilon_{i})}\text{, and }R_{i} = \frac{(1 + \upsilon_{i + 1})E_{i}}{(1 + \upsilon_{i})E_{i + 1}}$
(as before), and rearranging the equation:

$\begin{array}{r}
{\alpha G}_{i}\left\{ \ \left\lbrack A_{i} + C_{i}(2\upsilon_{i}\  + \alpha\ h_{i}) \right\rbrack\ e^{{\alpha h}_{i}} + \left\lbrack B_{i} - D_{i}(2\upsilon_{i}\  - \alpha\ h_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 - l_{i}\left\{ \ \left\lbrack A_{i} + C_{i}(1 + {\alpha h}_{i}) \right\rbrack\ e^{{\alpha h}_{i}} - \ \left\lbrack B_{i} - D_{i}(1 - {\alpha h}_{i}) \right\rbrack\ e^{- {\alpha h}_{i}} \right\} \\
 =  - R_{i}l_{i}\left\{ A_{i + 1} + C_{i + 1} - B_{i + 1} + D_{i + 1} \right\}
\end{array}$ (5-7)

Setting *l~i~* = 1, *G~i~* = 0, reduces the equation to the radial
deflection continuity equation derived for a fully bonded interface.
Setting *l~i~* to zero decouples the coefficients in layer *i* from
those in layer *i*+1. This causes numerical difficulties with a fully
unbonded interface, as discussed later.

Replacing the radial deflection continuity equation with the new spring
equilibrium equation (relating shear stress to radial deflection).

   $\begin{array}{r}
\begin{matrix}
\sigma_{z} \\
\tau_{\text{zr}} \\
w \\
\tau_{\text{zr}}\text{ vs}\text{.}u
\end{matrix} \\

\end{array}$ (5-8)

For a subgrade of infinite depth, the coefficients *A~N~* and *C~N~* are
zero because the responses must decay to zero at infinity and the
$e^{\alpha z}$ terms increase unbounded as z increases (see, also,
Appendix A). The equation for the interface above the subgrade is then:

$\begin{array}{r}
 \\

\end{array}$ (5-9)

where:

+-------+--------------------------------------------------------------+
| ###   | = the subgrade (last) layer, and                             |
| # *N* |                                                              |
+-------+--------------------------------------------------------------+
| *i*   | = *N*-1                                                      |
+-------+--------------------------------------------------------------+

Assembling the equilibrium equations 5-3, 5-4, and 5-9, for a sample
structure having three layers with fully bonded interfaces gives the
matrix equation 5-10. The same banded structure is repeated as more
layers are added.

$\begin{pmatrix}
1 & - 1 & - (1 - 2\upsilon_{1}) & - (1 - 2\upsilon_{1}) & & & & & & \\
1 & 1 & 2\upsilon_{1} & - 2\upsilon_{1} & & & & & & \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{3,3}}e^{{\alpha h}_{1}} & a_{\text{3,4}}e^{- {\alpha h}_{1}} & - 1 & 1 & (1 - 2\upsilon_{2}) & (1 - 2\upsilon_{2}) & & \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{4,3}}e^{{\alpha h}_{1}} & a_{\text{4,4}}e^{- {\alpha h}_{1}} & - 1 & - 1 & - 2\upsilon_{2} & 2\upsilon_{2} & & \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{5,3}}e^{{\alpha h}_{1}} & a_{\text{5,4}}e^{- {\alpha h}_{1}} & - R_{1} & - R_{1} & (2 - 4\upsilon_{2})R_{1} & - (2 - 4\upsilon_{2})R_{1} & & \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{6,3}}e^{{\alpha h}_{1}} & a_{\text{6,4}}e^{- {\alpha h}_{1}} & - R_{1} & R_{1} & - R_{1} & - R_{1} & & \\
 & & & & e^{{\alpha h}_{2}} & - e^{- {\alpha h}_{2}} & a_{\text{7,7}}e^{{\alpha h}_{2}} & a_{\text{7,8}}e^{- {\alpha h}_{2}} & 1 & (1 - 2\upsilon_{3}) \\
 & & & & e^{{\alpha h}_{2}} & e^{- {\alpha h}_{2}} & a_{\text{8,7}}e^{{\alpha h}_{2}} & a_{\text{8,8}}e^{- {\alpha h}_{2}} & - 1 & 2\upsilon_{3} \\
 & & & & e^{{\alpha h}_{2}} & e^{- {\alpha h}_{2}} & a_{\text{9,7}}e^{{\alpha h}_{2}} & a_{\text{9,8}}e^{- {\alpha h}_{2}} & - R_{2} & - (2 - 4\upsilon_{2})R_{2} \\
 & & & & e^{{\alpha h}_{2}} & - e^{- {\alpha h}_{2}} & a_{\text{10,7}}e^{{\alpha h}_{2}} & a_{\text{10,8}}e^{- {\alpha h}_{2}} & - R_{2} & R_{2}
\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}
1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix}$ (5-10)

The presence of exponentials in the equations can lead to very large and
very small numbers and consequent numerical difficulties due to
underflow and overflow. Introducing additional exponential functions can
reduce such difficulties. For example, multiplying each element in a
single column in the matrix by a given value divides the value of the
corresponding coefficient (*A*, *B*, *C*, or *D*) by the same amount.
Some control is thus given over the magnitudes of the elements of the
matrix and the magnitudes of the coefficients *A*, *B*, *C*, and *D*,
provided a corresponding adjustment is made in the elemental equations.
The procedure used in LEAF is to multiply the first and second columns
of the left hand side of each equilibrium equation by
$e^{- {\alpha h}_{\text{Oi},1}}$and multiply the second and fourth
columns of the left hand side by $e^{- {\alpha h}_{\text{Oi},2}}$.
Columns in the right hand sides are multiplied by
$e^{- {\alpha h}_{\text{Oi} + 1,1}}$ and
$e^{- {\alpha h}_{\text{Oi} + 1,2}}$. That is, the multiplying factor in
the right hand side is that for the left hand side of the next lower
layer. These operations are equivalent to shifting the origin of the *z*
coordinate, as discussed later. The procedure for selecting the values
of the exponents is given after describing how the equations are solved
for the coefficients *A*, *B*, *C*, and *D*.

The corresponding change in the response equation for *σ~z~*, for
example, is:

$\begin{array}{r}
\sigma_{z} = - \text{qa}\int_{}^{}{J_{0}(\alpha r)J_{1}(\alpha a)(A_{i}^{'} - C_{i}^{'}(1 - 2\upsilon_{i} - \alpha z))e^{- \alpha(h_{\text{Oi},1} - z)}d\alpha} \\
 + \text{qa}\int_{}^{}{J_{0}(\alpha r)J_{1}(\alpha a)(B_{i}^{'} + D_{i}^{'}(1 - 2\upsilon_{i} + \alpha z))e^{- \alpha(h_{\text{Oi},2} + z)}d\alpha}
\end{array}$ (5-11)

where:

  -------- --------------------------------------------------------------
  *σ~z~*   = vertical stress in layer *i*.

  *z*      = distance from the top of layer *i* and the primes denote
           that the numerical values of the *A*s, *B*s, *C*s, and *D*s
           are different after the multiplications.
  -------- --------------------------------------------------------------

Making the indicated changes in the equilibrium equations 5-8 gives.

$\begin{array}{r}
 \\

\end{array}$Equation (12)

Now assembling the previously used three-layer system with the
multiplying factors above the columns they are multiplied into shows how
the left and right sides of the equilibrium equations are associated
with the current and next lower layers, respectively. (The equation
should be read as if the multiplying factors are within the elements of
the full matrix. The matrix was not written out fully to keep it at a
reasonable size for display on a single page.)

   $\begin{array}{r}
\begin{pmatrix}
e^{- {\alpha h}_{O1,1}} & e^{- {\alpha h}_{O1,2}} & e^{- {\alpha h}_{O1,1}} & e^{- {\alpha h}_{O1,2}} & e^{- {\alpha h}_{O2,1}} & e^{- {\alpha h}_{O2,2}} & e^{- {\alpha h}_{O2,1}} & \ e^{- {\alpha h}_{O2,2}} & \ e^{- {\alpha h}_{O3,1}} & e^{- {\alpha h}_{O3,2}}
\end{pmatrix} \\
\begin{pmatrix}
1 & - 1 & a_{\text{1,3}} & a_{\text{1,4}} & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & a_{\text{2,3}} & a_{\text{2,4}} & 0 & 0 & 0 & 0 & 0 & 0 \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{3,3}}e^{{\alpha h}_{1}} & a_{\text{3,4}}e^{- {\alpha h}_{1}} & - 1 & \ 1 & a_{\text{3,7}} & a_{\text{3,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{4,3}}e^{{\alpha h}_{1}} & a_{\text{4,4}}e^{- {\alpha h}_{1}} & - 1 & \  - 1 & a_{\text{4,7}} & a_{\text{4,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{5,3}}e^{{\alpha h}_{1}} & a_{\text{5,4}}e^{- {\alpha h}_{1}} & - R_{1} & \  - R_{1} & a_{\text{5,7}} & a_{\text{5,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{6,3}}e^{{\alpha h}_{1}} & a_{\text{6,4}}e^{- {\alpha h}_{1}} & - R_{1} & \ R_{1} & - R_{1} & - R_{1} & 0 & 0 \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \  - e^{- {\alpha h}_{2}} & a_{\text{7,7}}e^{{\alpha h}_{2}} & a_{\text{7,8}}e^{- {\alpha h}_{2}} & 1 & \ a_{\text{7,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \ e^{- {\alpha h}_{2}} & a_{\text{8,7}}e^{{\alpha h}_{2}} & a_{\text{8,8}}e^{- {\alpha h}_{2}} & - 1 & \ a_{\text{8,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \ e^{- {\alpha h}_{2}} & a_{\text{9,7}}e^{{\alpha h}_{2}} & a_{\text{9,8}}e^{- {\alpha h}_{2}} & - R_{2} & \ a_{\text{9,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \  - e^{- {\alpha h}_{2}} & a_{\text{10,7}}e^{{\alpha h}_{2}} & a_{\text{10,8}}e^{- {\alpha h}_{2}} & R_{2} & - R_{2}
\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}
1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix}
\end{array}$ (5-13)

Multiplying the *A*s, *B*s, *C*s, and *D*s by the inverse of the
multiplying factors gives

$\begin{array}{r}
\begin{pmatrix}
e^{- {\alpha h}_{O1,1}} & e^{- {\alpha h}_{O1,2}} & e^{- {\alpha h}_{O1,1}} & e^{- {\alpha h}_{O1,2}} & e^{- {\alpha h}_{O2,1}} & e^{- {\alpha h}_{O2,2}} & e^{- {\alpha h}_{O2,1}} & \ e^{- {\alpha h}_{O2,2}} & \ e^{- {\alpha h}_{O3,1}} & e^{- {\alpha h}_{O3,2}}
\end{pmatrix} \\
\begin{pmatrix}
1 & - 1 & a_{\text{1,3}} & a_{\text{1,4}} & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & a_{\text{2,3}} & a_{\text{2,4}} & 0 & 0 & 0 & 0 & 0 & 0 \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{3,3}}e^{{\alpha h}_{1}} & a_{\text{3,4}}e^{- {\alpha h}_{1}} & - 1 & \ 1 & a_{\text{3,7}} & a_{\text{3,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{4,3}}e^{{\alpha h}_{1}} & a_{\text{4,4}}e^{- {\alpha h}_{1}} & - 1 & \  - 1 & a_{\text{4,7}} & a_{\text{4,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & e^{- {\alpha h}_{1}} & a_{\text{5,3}}e^{{\alpha h}_{1}} & a_{\text{5,4}}e^{- {\alpha h}_{1}} & - R_{1} & \  - R_{1} & a_{\text{5,7}} & a_{\text{5,8}} & 0 & 0 \\
e^{{\alpha h}_{1}} & - e^{- {\alpha h}_{1}} & a_{\text{6,3}}e^{{\alpha h}_{1}} & a_{\text{6,4}}e^{- {\alpha h}_{1}} & - R_{1} & \ R_{1} & - R_{1} & - R_{1} & 0 & 0 \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \  - e^{- {\alpha h}_{2}} & a_{\text{7,7}}e^{{\alpha h}_{2}} & a_{\text{7,8}}e^{- {\alpha h}_{2}} & 1 & \ a_{\text{7,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \ e^{- {\alpha h}_{2}} & a_{\text{8,7}}e^{{\alpha h}_{2}} & a_{\text{8,8}}e^{- {\alpha h}_{2}} & - 1 & \ a_{\text{8,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \ e^{- {\alpha h}_{2}} & a_{\text{9,7}}e^{{\alpha h}_{2}} & a_{\text{9,8}}e^{- {\alpha h}_{2}} & - R_{2} & \ a_{\text{9,10}} \\
0 & 0 & 0 & 0 & e^{{\alpha h}_{2}} & \  - e^{- {\alpha h}_{2}} & a_{\text{10,7}}e^{{\alpha h}_{2}} & a_{\text{10,8}}e^{- {\alpha h}_{2}} & R_{2} & - R_{2}
\end{pmatrix}\ \begin{pmatrix}
A_{1}e^{{\alpha h}_{O1,1}} \\
B_{1}e^{{\alpha h}_{O1,2}} \\
C_{1}e^{{\alpha h}_{O1,1}} \\
D_{1}e^{{\alpha h}_{O1,2}} \\
A_{2}e^{{\alpha h}_{O2,1}} \\
B_{2}e^{{\alpha h}_{O2,2}} \\
C_{2}e^{{\alpha h}_{O2,1}} \\
D_{2}e^{{\alpha h}_{O2,2}} \\
B_{3}e^{{\alpha h}_{O3,1}} \\
D_{3}e^{{\alpha h}_{O3,2}}
\end{pmatrix} = \begin{pmatrix}
1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix}
\end{array}$ (5-14)

and shows that
${A_{i}^{'} = A_{i}e^{{\alpha h}_{\text{Oi},1}}\text{,   \{}B}_{i}^{'} = B_{i}e^{{\alpha h}_{\text{Oi},2}}{\text{,   \{}C}_{i}^{'} = C_{i}e^{{\alpha h}_{\text{Oi},1}}{\text{, and  \{}D}_{i}^{'} = D_{i}e^{{\alpha h}_{\text{Oi},2}}$.

An alternative derivation is

$\begin{array}{r}
y_{1} = A_{i}e^{\alpha z} = A_{i}^{'}e^{- \alpha(h_{\text{Oi},1} - z)} = A_{i}^{'}e^{- {\alpha h}_{\text{Oi},1}}e^{\alpha z} \\
\therefore A_{i}^{'} = A_{i}e^{{\alpha h}_{\text{Oi},1}}\text{, and } \\
y_{2} = B_{i}e^{- \alpha z} = B_{i}^{'}e^{- \alpha(h_{\text{Oi},2} + z)} = B_{i}^{'}e^{- {\alpha h}_{\text{Oi},2}}e^{- \alpha z} \\
\therefore B_{i}^{'} = B_{i}e^{{\alpha h}_{\text{Oi},2}}
\end{array}$ (5-15)

Applying the multiplying factors therefore corresponds to changes in the
independent variable of

$\begin{array}{r}
z_{1} = (h_{\text{Oi},1} - z)\text{  and} \\
z_{2} = (h_{\text{Oi},2} + z)
\end{array}$ (5-16)

with origin shifts of *h~Oi,~*~1~ and --*h~Oi,~*~2~, respectively. The
same relationships hold for *C~i~* and *D~i~*. Figures 5-3 and 5-4
illustrate the relationships between the original and transformed
variables.

FIGURE 5-3.  ORIGIN SHIFT FOR POSITIVE EXPONENTS,
$y_{1} = \text{Ae}^{\alpha z} = A^{'}e^{- {\alpha z}_{1}}$

FIGURE 5-4.  ORIGIN SHIFT FOR NEGATIVE EXPONENTS,
$y_{2} = \text{Be}^{- \alpha z} = B^{'}e^{- {\alpha z}_{2}}$

## 5.2.2 Final Manipulation of the Equilibrium Equations.

The left hand side of the equilibrium equation 5-12 has a simple
analytic inverse. The complete equation is rewritten in equation 5-17
with the left and right hand sides premultiplied by the inverse of the
left hand side. (The first matrix on the right hand side of 5-17 is the
inverse of the left hand side of 5-12.)

$\begin{array}{r}
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
A_{i}^{'} \\
B_{i}^{'} \\
C_{i}^{'} \\
D_{i}^{'}
\end{pmatrix} = \frac{1}{4l_{i}(1 - \upsilon_{i})} \\
 \\

\end{array}$Equation (17)

Multiplying the two 4x4 matrices gives terms in the elements of the
product matrix of:

$\begin{array}{r}
e^{\alpha(h_{\text{Oi},1} - h_{i})} \times e^{- \alpha(h_{\text{Oi} + 1,1})} = e^{\alpha(h_{\text{Oi},1} - h_{\text{Oi} + 1,1} - h_{i})} \\
e^{\alpha(h_{\text{Oi},2} + h_{i})} \times e^{- \alpha(h_{\text{Oi} + 1,1})} = e^{\alpha(h_{\text{Oi},2} - h_{\text{Oi} + 1,1} + h_{i})} \\
e^{\alpha(h_{\text{Oi},1} - h_{i})} \times e^{- \alpha(h_{\text{Oi} + 1,2})} = e^{\alpha(h_{\text{Oi},1} - h_{\text{Oi} + 1,2} - h_{i})} \\
e^{\alpha(h_{\text{Oi},2} + h_{i})} \times e^{- \alpha(h_{\text{Oi} + 1,2})} = e^{\alpha(h_{\text{Oi},2} - h_{\text{Oi} + 1,2} + h_{i})}
\end{array}$ (5-18)

To stop these terms from increasing with *α*, the following conditions
must hold:

$\begin{array}{r}
h_{\text{Oi} + 1,1} \geq (h_{\text{Oi},2} + h_{i})\text{, and} \\
h_{\text{Oi} + 1,2} \geq (h_{\text{Oi},2} + h_{i})
\end{array}$ (5-19)

The current condition used in LEAF is that equality be satisfied for the
above two equations. However, the conditions required to minimize
numerical difficulties are very poorly understood. Suitable conditions
probably vary with different configurations of the structure (number and
thickness of layers, material properties, and interface parameter) and
different coordinate values for the evaluation point, particularly when
numerical stability has to be traded off against solution time.

# 5.3 Solution of the Interface Equilibrium Equations.

Because of the simple banded nature of the equations for a complete
multi-layered system and, secondarily, the existence of a simple
analytic inverse for the left hand side of the equilibrium equations,
the equations can be solved more efficiently than by the use of a
typical linear equation solver. Using a three-layer system for
illustration, assembly of the complete system of equations gives:

$\begin{pmatrix}
1 & 0 & 0 & 0 & a_{1,5} & a_{1,6} & a_{1,7} & a_{1,8} & 0 & 0 \\
0 & 1 & 0 & 0 & a_{2,5} & a_{2,6} & a_{2,7} & a_{2,8} & 0 & 0 \\
0 & 0 & 1 & 0 & a_{3,5} & a_{3,6} & a_{3,7} & a_{3,8} & 0 & 0 \\
0 & 0 & 0 & 1 & a_{4,5} & a_{4,6} & a_{4,7} & a_{4,8} & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & a_{5,9} & a_{5,\text{10}} \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & a_{6,9} & a_{6,\text{10}} \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & a_{7,9} & a_{7,\text{10}} \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & a_{8,9} & a_{8,\text{10}} \\
a_{9,1} & a_{9,2} & a_{9,3} & a_{9,4} & 0 & 0 & 0 & 0 & 0 & 0 \\
a_{\text{10},1} & a_{\text{10},2} & a_{\text{10},3} & a_{\text{10},4} & 0 & 0 & 0 & 0 & 0 & 0
\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}
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
1 \\
0
\end{pmatrix}$ (5-20)

The vertical and shear stress equations for the top surface have been
moved to the bottom of the coefficients matrix. The matrix can now be
partitioned for part inversion \[[^1]\].

$$\begin{pmatrix}
A_{\text{11}} & A_{\text{12}} \\
A_{\text{21}} & A_{\text{22}}
\end{pmatrix}\begin{pmatrix}
x \\
y
\end{pmatrix} = \begin{pmatrix}
X \\
Y
\end{pmatrix}$$

where *A*~11~ is the top left 8×8 partition, *A*~22~ is the bottom right
2×2 partition (all zeros), and the remaining two partitions (2x8 and 8x2
respectively) are the rectangular matrices representing the top surface
and the subgrade layer. Also:

$$x = \begin{pmatrix}
A_{1}^{'} \\
B_{1}^{'} \\
C_{1}^{'} \\
D_{1}^{'} \\
A_{2}^{'} \\
B_{2}^{'} \\
C_{2}^{'} \\
D_{2}^{'}
\end{pmatrix};\ X = \begin{pmatrix}
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0 \\
0
\end{pmatrix};\ y = \begin{pmatrix}
B_{3}^{'} \\
D_{3}^{'}
\end{pmatrix};\ Y = \begin{pmatrix}
1 \\
0
\end{pmatrix}$$

Expanding the matrix equation:

$$\begin{array}{r}
A_{\text{11}}x + A_{\text{12}}y = X = 0 \\
A_{\text{21}}x + A_{\text{22}}y = Y
\end{array}$$

and part inverting:

$$\begin{array}{r}
x = - A_{\text{11}}^{- 1}A_{\text{12}}y \\
A_{\text{21}}x = Y \\
 - A_{\text{21}}A_{\text{11}}^{- 1}A_{\text{12}}y = Y
\end{array}$$

The product $A_{\text{11}}^{- 1}A_{\text{12}}$ is an 8×2 matrix which
can be formed by back substitution of the columns of $A_{\text{12}}$
into $A_{\text{11}}$.

$$\begin{array}{r}
 - A_{\text{11}}(A_{\text{11}}^{- 1}A_{\text{12}})y + A_{\text{12}}y = 0 \\
\therefore A_{\text{11}}(A_{\text{11}}^{- 1}A_{\text{12}}) = A_{\text{12}}
\end{array}$$

The structure of the equation
$A_{\text{11}}(A_{\text{11}}^{- 1}A_{\text{12}}) = A_{\text{12}}$for a
three layer system is:

$\begin{pmatrix}
1 & 0 & 0 & 0 & a_{1,5} & a_{1,6} & a_{1,7} & a_{1,8} \\
0 & 1 & 0 & 0 & a_{2,5} & a_{2,6} & a_{2,7} & a_{2,8} \\
0 & 0 & 1 & 0 & a_{3,5} & a_{3,6} & a_{3,7} & a_{3,8} \\
0 & 0 & 0 & 1 & a_{4,5} & a_{4,6} & a_{4,7} & a_{4,8} \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
b_{1} & c_{1} \\
b_{2} & c_{2} \\
b_{3} & c_{3} \\
b_{4} & c_{4} \\
b_{5} & c_{5} \\
b_{6} & c_{6} \\
b_{7} & c_{7} \\
b_{8} & c_{8}
\end{pmatrix} = \begin{pmatrix}
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
a_{5,9} & a_{5,\text{10}} \\
a_{6,9} & a_{6,\text{10}} \\
a_{7,9} & a_{7,\text{10}} \\
a_{8,9} & a_{8,\text{10}}
\end{pmatrix}$ (5-21)

Where the right hand side is $A_{\text{12}}$ and *b* and *c* are the
columns of $A_{\text{11}}^{- 1}A_{\text{12}}$.

Full inversion of *A*~11~, or, equivalently, equation solution, is not
necessary because the lower triangle is already zero. (Numerical
efficiency can be further improved because the diagonal elements are
unity and the upper triangle contains zeros.) Multiplication by *A*~21~,
a 2×8 matrix, gives a 2×2 matrix, $A^{}$, and a set of 2 equations in 2
unknowns:

$\begin{array}{r}
 - A_{\text{21}}A_{\text{11}}^{- 1}A_{\text{12}}y = Y \\
A_{\text{21}}A_{\text{11}}^{- 1}A_{\text{12}}y = - Y \\
\begin{pmatrix}
A_{\text{11}}^{} & A_{\text{12}}^{} \\
A_{\text{21}}^{} & A_{\text{22}}^{}
\end{pmatrix}\begin{pmatrix}
y_{1} \\
y_{2}
\end{pmatrix} = \begin{pmatrix}
 - Y_{1} \\
 - Y_{2}
\end{pmatrix} \\
\begin{pmatrix}
A_{\text{11}}^{} & A_{\text{12}}^{} \\
A_{\text{21}}^{} & A_{\text{22}}^{}
\end{pmatrix}\begin{pmatrix}
B_{3}^{'} \\
D_{3}^{'}
\end{pmatrix} = \begin{pmatrix}
 - 1 \\
0
\end{pmatrix}
\end{array}$ (5-22)

with solution

$\begin{array}{r}
B_{3}^{'} = \frac{- 1}{A_{\text{11}}^{} - \frac{A_{\text{12}}^{}A_{\text{21}}^{}}{A_{\text{22}}^{}}} \\
D_{3}^{'} = - \frac{A_{\text{21}}^{}B_{3}^{'}}{A_{\text{22}}^{}}
\end{array}$ (5-23)

The rest of the constants for solving the layered elastic equations are
found directly from $x = - A_{\text{11}}^{- 1}A_{\text{12}}y$. The
product $A_{\text{11}}^{- 1}A_{\text{12}}$ is stored in a separate array
and the computation of each constant requires 2 multiplications and 2
additions.

$x = - A_{\text{11}}^{- 1}A_{\text{12}}y;\ \text{or }\begin{pmatrix}
A_{1}^{'} \\
B_{1}^{'} \\
C_{1}^{'} \\
D_{1}^{'} \\
A_{2}^{'} \\
B_{2}^{'} \\
C_{2}^{'} \\
D_{2}^{'}
\end{pmatrix} = \begin{pmatrix}
b_{1} & c_{1} \\
b_{2} & c_{2} \\
b_{3} & c_{3} \\
b_{4} & c_{4} \\
b_{5} & c_{5} \\
b_{6} & c_{6} \\
b_{7} & c_{7} \\
b_{8} & c_{8}
\end{pmatrix}\begin{pmatrix}
 - B_{3}^{'} \\
 - D_{3}^{'}
\end{pmatrix}$ (5-24)

# 5.4 Integration of the Elemental Response Equations.

After solving the boundary condition equations for the constants *A*,
*B*, *C*, and *D*, the elemental response equations are integrated. (The
primes are dropped in the following although it should still be assumed
that the origin shift is applied as necessary). Analytic integrals do
not exist for the general form of the response equations and numerical
integration must be used. Gauss-Laguerre integration was selected for
LEAF. Using, as before, the response equation for vertical stress as an
example:

$\begin{array}{r}
\sigma_{z} = - \text{qa}\int_{}^{}{J_{0}(\alpha r)J_{1}(\alpha a)(A_{i} - C_{i}(1 - 2\upsilon_{i} - \alpha z))e^{- \alpha(h_{\text{Oi},1} - z)}d\alpha} \\
 + \text{qa}\int_{}^{}{J_{0}(\alpha r)J_{1}(\alpha a)(B_{i} + D_{i}(1 - 2\upsilon_{i} + \alpha z))e^{- \alpha(h_{\text{Oi},2} + z)}d\alpha}
\end{array}$ (5-24)

Defining new variables (see also equations 5-16):

$\begin{array}{r}
\alpha_{1} = \alpha(h_{\text{Oi},1} - z) = {\alpha z}_{1} \\
d\alpha = \frac{{d\alpha}_{1}}{z_{1}}\text{, and} \\
\alpha_{2} = \alpha(h_{\text{Oi},2} + z) = {\alpha z}_{2} \\
d\alpha = \frac{{d\alpha}_{2}}{z_{2}}
\end{array}$ (5-25)

and substituting into the full integral gives:

$\begin{array}{r}
\sigma_{z} = - \text{qa}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 2\upsilon_{i} - \frac{{z\alpha}_{1}}{z_{1}} \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
 + \text{qa}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 1 - 2\upsilon_{i} + \frac{{z\alpha}_{2}}{z_{2}} \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}}
\end{array}$ (5-26)

The Gauss-Laguerre method replaces an integral over an infinite range
with a summation over a finite range (see, for example, Numerical
Methods). By this method, the integral is evaluated with the expression:

$\begin{array}{r}
\sigma_{z} = - \text{qa}\sum_{j = 1}^{N}{J_{0}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1j}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 2\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}} \\
 + \text{qa}\sum_{j = 1}^{N}{J_{0}\left( \frac{{r\alpha}_{2j}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2j}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 1 - 2\upsilon_{i} + \frac{{z\alpha}_{2j}}{z_{2}} \right) \right\}\frac{w_{j}}{z_{2}}}
\end{array}$ (5-27)

Where *α*~1*j*~ and *α*~2*j*~ are the abscissa values for the summation
of *N* points and *w~j~* are weights corresponding to the *α*~1*j*~ and
*α*~2*j*~ values.

One of the conditions for proper convergence of the integrals is that,
since the Bessel functions are periodic, the increment in the *α~j~*
values must be small enough to put more than two samples in each period
to avoid aliasing. Figures 5-5 through 5-7 illustrate the problem.
Figure 5-5 shows plots of *J*~0~ and *J*~1~. Figures 5-6 and 5-7 show
the components of the integral of vertical strain for an extreme set of
structure and load conditions. The evaluation point is at the top of the
subgrade and a radius of 150 inches (3.81 m) from the center of the
load. Load radius is 9.4 inches (239 mm), load pressure is 180 psi
(1,240 kPa), and all interfaces are fully bonded. The material
properties for the layers are shown in table 5-1.

TABLE 5-1.  STRUCTURE LAYER PROPERTIES USED TO ILLUSTRATE ALIASING IN
INTEGRATION OF THE BESSEL FUNCTIONS.

  ------------------ ---------------- ---------------- -------------------
  Layer Number       Thickness, in    Modulus, psi     Poisson's Ratio

  1                  5                200,000          0.35

  2                  6.25             73,690           0.35

  3                  6.25             47,220           0.35

  4                  1                25,000           0.35

  5                  Infinite         22,500           0.35
  ------------------ ---------------- ---------------- -------------------

![](media/image1.wmf)

FIGURE 5-5.  PLOTS OF ZERO AND FIRST ORDER BESSEL FUNCTIONS

![](media/image2.wmf)

FIGURE 5-6.  PLOTS OF *w~j~*, (*B*+*D*)*/z*~2~, LOAD FUNCTION
*J*~1~(*αa/z*~2~), AND *J*~0~(*αr/z*~2~) FOR *r* = 150 INCHES AND
*N *= 2,000

![](media/image3.wmf)

FIGURE 5-7.  PLOTS OF *w~j~*, (*B*+*D*)*/z*~2~, LOAD FUNCTION
*J*~1~(*αa/z*~2~), AND *J*~0~(*αr/z*~2~) FOR *r* = 150 INCHES AND
*N* = 200

Each marker in the plots in figures 5-6 and 5-7 corresponds to one
sample value *α~j~*. In figure 5-6, with *N* = 2,000, it can be seen
that up to *α* = 0.5 there are two or more sample points in each period
of *J*~0~(*α~2~r/z*~2~). Thereafter, the increments in *α* increase in
magnitude and aliasing occurs. The aliasing would clearly cause a
significant error in the summation, but the (*B*+*D*)*/z*~2~ term decays
rapidly past *α* = 0.5. At *α* = 0.48, (*B*+*D*)*/z*~2~ = ‑0.018, and at
*α* = 0.95, (*B*+*D*)*/z*~2~ = ‑0.00057. In figure 5-7, with N = 200,
aliasing occurs as early as the second period and the error in the full
summation is unacceptably large, despite the rapid decay in the
(*B*+*D*)*/z*~2~ term. (Values for *A* and *C* are not given because
they are zero when the evaluation point is located in the subgrade.)

Instead of increasing the value of *N* when required to achieve the
desired accuracy, the strategy used in LEAF is to change the values of
*z*~1~ and *z*~2~ by changing the values of *h~O~*~,1~ and *h~O~*~,2~
for the layer in which the evaluation point is located (see equations
5-16). The condition is that the number of sample points in one period
should always be greater than 3 to 4. At the current writing, the
settings in LEAF are:

$\begin{array}{r}
N = \text{500} \\
h_{\text{Ol},1} = \frac{\ r_{\text{max}}}{4} + z \\
h_{\text{Ol},2} = \frac{\ r_{\text{max}}}{4} + z
\end{array}$ (5-28)

where:

  ---------- --------------------------------------------------------------
  *l*        refers to the evaluation layer.

  *z*        = the distance from the top of the evaluation layer to the
             evaluation point.

  *r~max~*   = the greater of the load radius, *a*, and the radius to the
             evaluation point, *r*.
  ---------- --------------------------------------------------------------

The value of *h~Ol~*~,1~ is changed only if the new value is greater
than the existing value found from equations 5-19. If *h~Ol~*~,1~ is
changed, then the values of *h~Oi~*~,1~ and *h~Oi~*~,2~ for all other
layers are changed to satisfy equations 5-19 in relation to the new
values in the evaluation layer.

# 5.5 Computer program implementation of the integrations.

Computer program implementation of Gauss-Laguerre integration is
basically the same for stresses, strains, and displacements. The only
difference is that the equations evaluated in the program loops are
different. A primary objective was to reduce the run-time for the
pavement thickness design procedures implemented in the computer program
LEDFAA. These implementations require the computation of one, or, at
most, two responses for a single pavement structure of rigid or flexible
type with or without an overlay and with many landing gear
configurations in the design aircraft mix. The strategy adopted was
based on the following.

1.  Allow for independent computation of only those responses required
    for a particular type of pavement design. For example, new flexible
    design requires computation of only vertical strain at the top of
    the subgrade in order to iterate to the design thickness.

2.  Arrange the program loops so that individual elements of the
    integrand are computed the minimum number of times.

3.  Assume that wheel loads and tire pressures are the same for each
    wheel in an individual landing gear wheel-group.

An outline of the program code for vertical strain computation is given
below to illustrate the implementation of the integration subroutines.

An elemental equation has not so far been derived for vertical strain.
The required equation is derived from the set of equations 5-2 in either
of two ways: $\varepsilon_{z} = \frac{\partial w}{\partial z}$ or
$\varepsilon_{z} = \frac{1}{E}(\sigma_{z} - \upsilon(\sigma_{r} + \sigma_{t}))$.
The resulting integral for computing vertical strain is:

$\begin{array}{r}
\varepsilon_{z} = - \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1}}{z_{1}} \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
 + \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2}}{z_{2}} \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}}
\end{array}$ (5-29)

and the Gauss-Laguerre summation is:

$\varepsilon_{z} = \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\begin{array}{r}
\left\lbrack - \sum_{j = 1}^{N}{J_{0}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1j}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}} \middle| \right\rbrack \\

\end{array}$ (5-30)

Data is passed to the subroutine for a single structure and a complete
traffic mix. The traffic mix has NAC aircraft. An "aircraft" is a group
of NTires wheels, usually representing a landing gear, but not
necessarily. For example, all sixteen main landing gear wheels on a
B-747 could be defined as a wheel-group. Each wheel in an aircraft
wheel-group has the same wheel load and tire pressure. Each aircraft
wheel-group has NevalPoints evaluation points at a single depth measured
from the top of the pavement structure. Vertical strain is to be
calculated at each evaluation point for each aircraft.

The expression to be summed can be split into four levels. The levels
are, in decreasing order in the hierarchy:

1.  Under the assumption that all wheels on a landing gear have the same
    load and the same tire pressure, the term
    $\frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}$ only needs to be
    computed once. If pressure and/or load are different on one or more
    of the wheels the *qa* multiplier would have to be moved to level 3.

2.  Terms like
    $\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}$
    only need to be computed once every integration point. That is, they
    are common to every evaluation point on every wheel. These terms are
    also the most expensive to compute because the coefficients *A*,
    *B*, *C*, and *D*, must be computed from the interface equations
    each time. Special cases can also apply to these terms. For example,
    in the computation of vertical strain at the top of the subgrade
    only the term containing *B* and *D* needs to be computed because
    *A* and *C* are zero in the subgrade layer.

3.  Terms like$J_{1}\left( \frac{{a\alpha}_{1j}}{z_{1}} \right)$,
    containing *aα*, where *a* is the tire radius, only need to be
    computed once for each wheel. That is, they are common to every
    evaluation point on each wheel taken separately. However, if both
    load and pressure are the same on all of the wheels on a gear then
    the *J*~1~ terms can be shifted up to level 2.

4.  Terms like $J_{0}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)$ and
    $J_{1}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)$must be calculated
    for every evaluation point on every wheel.

Computation loops in the program are therefore ordered as shown in the
Visual Basic code fragment below.

' Depth from top of layer to evaluation point.

ZLayer = ZEval - ZInterface(I - 1)

' RMax = maximum of evaluation point radius or load radius.

Call SetOShifts(OShift(), I, ZLayer, RMax)

Z1 = OShift(I, 1) - ZLayer

Z2 = OShift(I, 2) + ZLayer

ZLayer1 = ZLayer / Z1

ZLayer2 = ZLayer / Z2

Poisx2 = Poissons(EvalLayer) \* 2

IConstants = (EvalLayer - 1) \* 4

' GLAlpha() = Gauss-Laguerre abscissa values.

' GLWeight() = Gauss-Laguerre weight values.

For IG = 1 To GLNGauss ' Gauss-Laguerre integration points.

StrainWIGforConverge = 0 ' Initialize convergence criterion.

If EvalLayer \< NLayers Then ' A and C are zero for subgrade.

' Put constants for all layers in B().

Call FindConstants(GLAlpha(IG) / Z1, B(), OShift(), BadCondition)

AK = B(IConstants + 1) ' A for evaluation layer.

AlphaZ = GLAlpha(IG) \* ZLayer1

CK = B(IConstants + 3) \* (1 - Poisx2 \* 2 - AlphaZ)

StrainWIG1 = -(AK - CK) \* GLWeight(IG) / Z1

=$- \left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}$

StrainWIGforConverge = (Abs(AK) + Abs(CK)) \* GLWeight(IG) / Z1

End If

' Do the same for constants B and D.

Call FindConstants(GLAlpha(IG) / Z2, B(), OShift(), BadCondition)

BK = B(IConstants + 2)

AlphaZ = GLAlpha(IG) \* ZLayer2

DK = B(IConstants + 4) \* (1 - Poisx2 \* 2 + AlphaZ)

StrainWIG2 = (BK + DK) \* GLWeight(IG) / Z2

=
$\left\{ B_{i} + D_{i}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2j}}{z_{2}} \right) \right\}\frac{w_{j}}{z_{2}}$

StrainWIGforConverge = StrainWIGforConverge \_

\+ (Abs(BK) + Abs(DK)) \* GLWeight(IG) / Z2

' Loop over all aircraft.

For IAC = 1 To NAC

If Not ACConverged(IAC) Then

' Skip if all evaluation points for aircraft IAC have converged.

A2 = TireRadius(IAC, 1) \' Move down 1 loop for varying pressures.

A1 = A2 / Z1

A2 = A2 / Z2

If EvalLayer \< NLayers Then

' First line below is load function for uniform pressure.

' Second line is for selectable load functions. Parabolic is currently

' available as an alternative to uniform pressure.

' J1AlphaA1 = bessj1(GLAlpha(IG) \* A1) \* StrainWIG1

J1AlphaA1 = LoadFunction(GLAlpha(IG) \* A1) \* StrainWIG1

=$- J_{1}\left( \frac{{a\alpha}_{1j}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}$

End If

J1AlphaA2 = LoadFunction(GLAlpha(IG) \* A2) \* StrainWIG2

=$J_{1}\left( \frac{{a\alpha}_{2j}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2j}}{z_{2}} \right) \right\}\frac{w_{j}}{z_{2}}$

' Loop over all tires on each gear (aircraft).

For ITire = 1 To NTires(IAC)

' Loop over all evaluation points on each gear (aircraft).

For IEval = 1 To NEvalPoints(IAC)

' R2 is Horizontal distance between evaluation point and center of

' tire. Previously computed.

R2 = Radius(IAC, ITire, IEval)

R1 = R2 / Z1

R2 = R2 / Z2

If EvalLayer \< NLayers Then

J0AlphaR1 = bessj0(GLAlpha(IG) \* R1)

=$J_{0}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)$

End If

J0AlphaR2 = bessj0(GLAlpha(IG) \* R2)

=$J_{0}\left( \frac{{r\alpha}_{2j}}{z_{2}} \right)$

' StrainWI = elemental strain for current integration point.

StrainWI = J0AlphaR1 \* J1AlphaA1 + J0AlphaR2 \* J1AlphaA2

' Accumulate sum at each evaluation point for each gear.

StrainW(IAC, IEval) = StrainW(IAC, IEval) + StrainWI

' Convergence criterion always greater than zero and conservative.

If Abs(StrainWIGforConverge / StrainW(IAC, IEval)) \< 0.00001 Then

If Not Converged(IAC, ITire, IEval) Then

' Keep track of the number of evaluation points which have

' converged.

NEvalsConverged(IAC, ITire) = NEvalsConverged(IAC, ITire) + 1

If NEvalsConverged(IAC, ITire) = NEvalPoints(IAC) Then

' All evaluation points for tire Itire on aircraft IAC have

' converged.

NtiresConverged(IAC) = NtiresConverged(IAC) + 1

IterationstoConverge = IG

End If

Converged(IAC, Itire, IEval) = True

If NtiresConverged(IAC) = Ntires(IAC) Then

ACConverged(IAC) = True

NACConverged = NACConverged + 1

End If

End If

End If

Next IEval

Next ITire

End If

Next IAC

If NACConverged = NAC Then

' All evaluation points on all aircraft have converged.

NConverge = IG

Exit For ' Quit integration.

Next IG

For IAC = 1 To NAC

Factor = TirePress(IAC, 1) \* TireRadius(IAC, 1)

Factor = Factor \* (1 + Poissons(EvalLayer)) / Youngs(EvalLayer)

=$\frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}$

For IEval = 1 To NEvalPoints(IAC)

' Apply final factor to compute strains at all evaluation points.

StrainW(IAC, IEval) = StrainW(IAC, IEval) \* Factor

If StrainW(IAC, IEval) \> StrainWmax(IAC) Then

' Maximum strain for each gear (aircraft).

StrainWmax(IAC) = StrainW(IAC, IEval)

End If

Next Ieval

Next IAC

The second IAC loop is more complicated when computing horizontal
responses because the responses at a single evaluation point from each
tire must be transformed into rectangular coordinates before they can be
summed. Vertical responses can be summed directly in the first loop as
shown.

For the vertical strain integral, convergence is determined from the
expression:

$\frac{\left\{ \left| A_{\text{iIG}} \right| + \left| C_{\text{iIG}}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1\text{IG}}}{z_{1}} \right) \right| \right\}\frac{w_{\text{IG}}}{z_{1}} + \left\{ \left| B_{\text{iIG}} \right| + \left| D_{\text{iIG}}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2\text{IG}}}{z_{2}} \right) \right| \right\}\frac{w_{\text{IG}}}{z_{2}}}{\begin{array}{r}
\left| \sum_{j = 1}^{\text{IG}}{J_{0}\left( \frac{{r\alpha}_{1j}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1j}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}\frac{w_{j}}{z_{1}}} \middle| \right| \\

\end{array}}$ (5-31)

Where IG is the current integration step number. Convergence is assumed
when the inequality is satisfied. Similar expressions are derived for
the other integrals.

The convergence criterion is based on the following observations:

The Bessel functions are oscillatory functions and including them in the
convergence criterion can lead to a zero of a Bessel function
prematurely terminating the integration unless some kind of running
total is used. All Bessel functions have a decaying envelope which is
always between +/- 1. An asymptotic approximation is
$J_{n}(\alpha) = \sqrt{\frac{2}{\text{πα}}}\text{cos}\left( \alpha - \frac{n\pi}{2} - \frac{\pi}{4} \right)$for
α \> 25. The multiplying factor $\sqrt{\frac{2}{\text{πα}}}$ has values
of 0.16 and 0.08 for α = 25 and 100 respectively. The envelope therefore
decays slowly and has very little influence on convergence of the
integrals except when the evaluation depth is close to the top of the
top layer. Since the absolute value of the magnitude of the bessel
functions is always less than 1, ignoring the Bessel functions in
convergence is always conservative.

It is likely that the factors
$\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1j}}{z_{1}} \right) \right\}$
and
$\left\{ B_{i} + D_{i}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2j}}{z_{2}} \right) \right\}$
are never zero, but this has not been proved. Taking the sum of the
absolute values will, under this assumption, therefore guarantee that
false termination of the integration will not occur, and also gives a
conservative convergence criterion.

# 5.6 Fully Unbonded Interfaces.

With a fully unbonded interface, the interface stiffness is zero.
Radial/vertical shear stress is therefore zero and the contacting layers
are uncoupled in the boundary condition equations. Physically, the
layers are completely free to slide across each other and one of the
conditions for response in one of the layers must be independent of the
corresponding response condition in the other layer.

Solution of the boundary condition equations using an equation solver
with full pivoting works, but is slow. For an efficient solution, the
offending equations must be eliminated from the boundary condition
equations for each interface. That is, rewrite the equations as three
equations in six unknowns and reformulate for solution. But this would
give a special condition in the implementation and consequent
complication in the computer program.

Setting the interface stiffness to a large enough value in the equations
as formulated above gives a numerically stable solution and is also more
reasonable physically (concrete pavement slabs do not freely slide on
the supporting layers except in exceptional circumstances). However,
selecting a standard interface stiffness for "fully unbonded" is
difficult because physical conditions are generally unknown for this
construction method and the spring model is probably a poor
representation of the physical conditions. Tests were run with the
layered elastic formulation given above for new rigid and unbonded
overlaid rigid pavement structures. Results for a typical new rigid case
(table 5-2), with a single load of 50,000 pounds (222.4 kN) (radius of
9.4 inches (239 mm) and pressure of 180 psi (1,240 kPa )) are shown in
table 5-3.

TABLE 5-2.  RIGID PAVEMENT STRUCTURE LAYER PROPERTIES USED TO ILLUSTRATE
EFFECTS OF INTERFACE PROPERTIES WITH ONE UNBONDED INTERFACE

  -------------- ------------- ------------- --------------- ---------------
  Layer Number   Thickness, in Modulus, psi  Poisson's Ratio Interface

  1              14            4,000,000     0.15            fully unbonded

  2              6             700,000       0.20            fully bonded

  3              6             35,430        0.35            fully bonded

  4              Infinite      15,000        0.35            
  -------------- ------------- ------------- --------------- ---------------

TABLE 5-3.  HORIZONTAL AND SHEAR STRESSES WITH DIFFERENT PARAMETERS AT
INTERFACE 1 FOR TWO DIFFERENT EVALUATION POINT RADII

+-----------+-----------+-----------+-----------+-----------+-----------+
| Interface | Interface | *X~E~*\*  |           | *X~E~* =  |           |
|           | P         | = 0,      |           | 0, *Y~E~* |           |
| S         | arameter, | *Y~E~* =  |           | = 100,    |           |
| tiffness, | *l*       | 10,       |           | *Z~E~* =  |           |
| *k*, psi  |           | *Z~E~* =  |           | 14\*      |           |
|           |           | 14\*\*    |           |           |           |
+-----------+-----------+-----------+-----------+-----------+-----------+
|           |           | H         | Shear     | H         | Shear     |
|           |           | orizontal | Stress,   | orizontal | Stress,   |
|           |           | Stress,   | *τ~yz~*,  | Stress,   | *τ~yz~*,  |
|           |           | *σ~x~*,   | psi       | *σ~x~*,   | psi       |
|           |           | psi       |           | psi       |           |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-12^   | 10^-12^   | 6.59      | 8         | 10.04     | 1         |
|           |           |           | .45×10^0^ |           | .61×10^0^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-11^   | 10^-11^   | 41.47     | 4         | 53.44     | -1        |
|           |           |           | .89×10^2^ |           | .11×10^2^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-10^   | 10^-10^   | 280.09    | -4        | 24.21     | -1        |
|           |           |           | .84×10^0^ |           | .19×10^0^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-9^    | 10^-9^    | 260.56    | 1.        | 22.80     | -9.       |
|           |           |           | 06×10^-1^ |           | 77×10^-3^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-8^    | 10^-8^    | 258.18    | -9.       | 23.02     | -1.       |
|           |           |           | 55×10^-3^ |           | 81×10^-2^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-7^    | 10^-7^    | 258.02    | 1.        | 23.01     | 8.        |
|           |           |           | 95×10^-3^ |           | 79×10^-4^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-6^    | 10^-6^    | 258.01    | -1.       | 23.01     | -1.       |
|           |           |           | 04×10^-4^ |           | 91×10^-5^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-5^    | 10^-5^    | 258.01    | -2.       | 23.01     | -2.       |
|           |           |           | 81×10^-4^ |           | 93×10^-6^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-4^    | 10^-4^    | 258.01    | -2.       | 23.01     | -6.       |
|           |           |           | 86×10^-4^ |           | 58×10^-6^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.0010    | 10^-3^    | 258.01    | -2.       | 23.01     | -5.       |
|           |           |           | 86×10^-4^ |           | 07×10^-6^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.0101    | 10^-2^    | 258.02    | -2.       | 23.01     | 6.        |
|           |           |           | 78×10^-4^ |           | 34×10^-6^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.1111    | 10^-1^    | 258.02    | -1.       | 23.02     | 1.        |
|           |           |           | 86×10^-4^ |           | 32×10^-4^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 1.0       | 0.50      | 258.04    | 6.        | 23.03     | 1.        |
|           |           |           | 17×10^-4^ |           | 24×10^-3^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 3.0       | 0.75      | 258.16    | 2.        | 23.12     | 3.        |
|           |           |           | 43×10^-3^ |           | 77×10^-3^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| infinity  | 1.00      | 166.21    | -1        | 16.15     | 8.        |
|           |           |           | .49×10^1^ |           | 74×10^-1^ |
+-----------+-----------+-----------+-----------+-----------+-----------+

\*\* *X~E~*, *Y~E~*, and *Z~E~* are the coordinates of the evaluation
point

\* Nominal value. Value used = 13.999986

Good numerical stability is shown for all values of the interface
parameter in the range 10^-1^ to 10^-7^. A default value of *l* = 0.01
is currently used as the default for new rigid pavements, which have one
unbonded interface under the concrete slab (layer 1).

The value of the evaluation depth, *Z~E~*, determines which layer will
be selected for computing the responses. For the example given, input
values less than 14 inches (355.6 mm) would be expected to select layer
1 and input values greater than 14 inches (355.6 mm) would be expected
to select layer 2. However, due to numerical round off effects, this is
not necessarily so. Table 5-4 shows the layer selected and the shear
stress at the evaluation depth for evaluation depths close to 14 inches
(355.6 mm) and with the interface parameter set to a value of 0.01.

TABLE 5-4.  COMPARISON OF INTERFACE SHEAR STRESSES FOR FULLY SLIDING AND
FULLY BONDED INTERFACE PARAMETER VALUES

  ------------------ ------------------ ----------------- ----------------
  Evaluation Depth   Evaluation Layer   *τ~yz~* with *l*  *τ~yz~* with *l*
                                        = 0.01            = 1.0

  13.999900          1                  -2.04×10^-3^      -14.8810

  13.999986          1                  -2.78×10^-4^      -14.8798

  13.999987          2                  2.59×10^-5^       -14.8798

  14.000072          2                  8.26×10^-5^       -14.8795
  ------------------ ------------------ ----------------- ----------------

Shear stress is also shown with the interface parameter set to a value
of one. Good continuity of shear stress is shown at the transition from
layer 1 to layer 2.

Rigid pavements with an unbonded overlay are represented in the current
LEDFAA design procedure as structures with two unbonded interfaces; one
under the new slab and one under the existing concrete layer. With this
condition in the layered elastic equations, stiffnesses of close to 1.0
are required at both interfaces in order to obtain reasonable numerical
stability combined with sufficiently accurate stress computation (the
interface parameter has a value of 0.5 when stiffness has a value of
1.0). Shear stresses at the interfaces are of the order of 10^-3^ psi
(6.9 Pa) with *k* = 1.0. The current default value in LEDFAA of *k* for
rigid-over-rigid overlay pavements is 1.0 lb/in^3^ (271 kN/m^3^).

Results for a typical overlay structure (table 5-5) are shown in table
5-6.

TABLE 5-5.  RIGID PAVEMENT STRUCTURE LAYER PROPERTIES USED TO ILLUSTRATE
EFFECTS OF INTERFACE PROPERTIES WITH TWO UNBONDED INTERFACES

  -------------- ------------- ------------- --------------- ---------------
  Layer Number   Thickness, in Modulus, psi  Poisson's Ratio Interface

  1              12            4,000,000     0.15            fully unbonded

  2              14            1,500,000     0.15            fully unbonded

  3              6             250,000       0.35            fully bonded

  4              6             35,430        0.35            fully bonded

  5              Infinite      6,000         0.35            
  -------------- ------------- ------------- --------------- ---------------

TABLE 5-6.  HORIZONTAL AND SHEAR STRESSES WITH DIFFERENT PARAMETERS AT
INTERFACES 1 AND 2

+-----------+-----------+-----------+-----------+-----------+-----------+
| Interface | Interface | *X~E~* =  |           | *X~E~* =  |           |
|           | P         | 0, *Y~E~* |           | 0, *Y~E~* |           |
| S         | arameter, | = 10,     |           | = 10,     |           |
| tiffness, | *l*       | *Z~E~* =  |           | *Z~E~* =  |           |
| *k*,      |           | 12\*      |           | 26\*\*    |           |
| lb/in^3^  |           |           |           |           |           |
+-----------+-----------+-----------+-----------+-----------+-----------+
|           |           | H         | Shear     | H         | Shear     |
|           |           | orizontal | Stress,   | orizontal | Stress,   |
|           |           | Stress,   | psi       | Stress,   | psi       |
|           |           | psi       |           | psi       |           |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-5^    | 10^-5^    | 442.57    | 6         | -0.17     | 1.        |
|           |           |           | .67×10^1^ |           | 50×10^-8^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 10^-4^    | 10^-4^    | -89.62    | 2         | 425.63    | 1.        |
|           |           |           | .61×10^1^ |           | 15×10^-5^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.0010    | 10^-3^    | 253.43    | 2         | 83.94     | -6.       |
|           |           |           | .18×10^1^ |           | 32×10^-5^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.0101    | 10^-2^    | 261.55    | 7.        | 103.92    | 1.        |
|           |           |           | 65×10^-1^ |           | 35×10^-4^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 0.1111    | 10^-1^    | 268.36    | 3.        | 104.06    | -1.       |
|           |           |           | 91×10^-3^ |           | 94×10^-5^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 1.0       | 0.50      | 269.04    | 5.        | 103.04    | 8.        |
|           |           |           | 41×10^-3^ |           | 60×10^-4^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| 3.0       | 0.75      | 269.30    | 3.        | 102.97    | 2.        |
|           |           |           | 50×10^-3^ |           | 85×10^-3^ |
+-----------+-----------+-----------+-----------+-----------+-----------+
| infinity  | 1.00      | 468.92    | -2.       | 61.56     | -2        |
|           |           |           | 77×10^-3^ |           | .53×10^0^ |
+-----------+-----------+-----------+-----------+-----------+-----------+

\* Nominal value. Value used = 11.999988

\*\* Nominal value. Value used = 25.999774

# 5.7 Computation of Responses.

## 5.7.1 LEDFAA and Backcalculation Responses.

LEDFAA currently requires the computation of three individual responses
in the thickness design procedures.

1.  Vertical strain at the top of the subgrade as the criterion of shear
    failure in the subgrade of a flexible pavement.

2.  Horizontal strain at the bottom of the surface asphalt layer as the
    criterion of fatigue failure in the surface asphalt layer of a
    flexible pavement.

3.  Horizontal stress at the bottom of the surface PCC layer as the
    criterion of fatigue failure in the surface PCC layer of a new rigid
    pavement and horizontal stress at the bottom of the existing PCC
    layer as the criterion of fatigue deterioration in the existing PCC
    layer of an overlaid rigid pavement.

However, elemental equations for responses of the axi-symmetric layered
elastic system include only vertical stress, *σ~z~*, radial stress,
*σ~r~*, tangential stress, *σ~t~*, vertical-radial shear stress,
*τ~zr~*, vertical deflection, *w*, and radial deflection, *u*. The
desired responses have to be computed indirectly from the available
response equations.

The derivation of the integral equation for vertical strain has already
been given:

$$\begin{array}{r}
\varepsilon_{z} = - \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 1 - 4\upsilon_{i} - \frac{{z\alpha}_{1}}{z_{1}} \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
 + \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 1 - 4\upsilon_{i} + \frac{{z\alpha}_{2}}{z_{2}} \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}}
\end{array}$$

and a subroutine listed for independent computation of vertical strain.

Horizontal stress is computed from the *X* and *Y* components of the
radial and tangential stresses, *σ~r~* and *σ~t~*.

$\begin{array}{r}
\sigma_{r} = \text{qa}\int_{}^{}{J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ 2\upsilon_{i}J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)C_{i} + \left( J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right) - \frac{J_{1}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)}{\left( \frac{{r\alpha}_{1}}{z_{1}} \right)} \right)\left( A_{i} + C_{i}\left( 1 + \frac{{z\alpha}_{1}}{z_{1}} \right) \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
 + \text{qa}\int_{}^{}{J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ 2\upsilon_{i}J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)D_{i} - \left( J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right) - \frac{J_{1}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)}{\left( \frac{{r\alpha}_{2}}{z_{2}} \right)} \right)\left( B_{i} - D_{i}\left( 1 - \frac{{z\alpha}_{2}}{z_{2}} \right) \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}}
\end{array}$ (5-32)

$\begin{array}{r}
\sigma_{t} = \text{qa}\int_{}^{}{J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ 2\upsilon_{i}J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)C_{i} + \frac{J_{1}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)}{\left( \frac{{r\alpha}_{1}}{z_{1}} \right)}\left( A_{i} + C_{i}\left( 1 + \frac{{z\alpha}_{1}}{z_{1}} \right) \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
 + \text{qa}\int_{}^{}{J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ 2\upsilon_{i}J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)D_{i} - \frac{J_{1}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)}{\left( \frac{{r\alpha}_{2}}{z_{2}} \right)}\left( B_{i} - D_{i}\left( 1 - \frac{{z\alpha}_{2}}{z_{2}} \right) \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}}
\end{array}$ (5-33)

If

  -------- --------------------------------------------------------------
  *x~E~* , = *X* and *Y* coordinates of the evaluation point
  *y~E~*   

  *x~T~* , = *X* and *Y* coordinates of the center of the tire load
  *y~T~*   

  *R*      = ((*x~E~* - *x~T~*)^2^ + (*y~E~* - *y~T~*)^2^)^1/2^

           = radius from the center of the tire load to the evaluation
           point

  cos *θ*  = (*x~E~* - *x~T~*) / *R*

  sin *θ*  = (*y~E~* - *y~T~*) / *R*
  -------- --------------------------------------------------------------

Then

$$\sigma_{x} = \sigma_{r}\text{cos}^{2}\theta + \text{sin}^{2}\theta$$

$$\sigma_{y} = \sigma_{r}\text{sin}^{2}\theta + \text{cos}^{2}\theta$$

$$\tau_{\text{xy}} = \text{sin}\theta\text{cos}\theta(\sigma_{r} - \sigma_{t})$$

$$\tau_{\text{xz}} = \tau_{\text{rz}}\text{cos}\theta$$

$$\tau_{\text{yz}} = \tau_{\text{rz}}\text{sin}\theta$$

$\sigma_{H\text{max}} = \frac{(\sigma_{x} + \sigma_{y})}{2} + \left( \left( \frac{\sigma_{x} + \sigma_{y}}{2} \right)^{2} + \tau_{\text{xy}}^{2} \right)^{\frac{1}{2}}$=
maximum horizontal stress

$\sigma_{H\text{min}} = \frac{(\sigma_{x} + \sigma_{y})}{2} - \left( \left( \frac{\sigma_{x} + \sigma_{y}}{2} \right)^{2} + \tau_{\text{xy}}^{2} \right)^{\frac{1}{2}}$=
minimum horizontal stress

The subroutine for horizontal stress sums the *σ~x~*, *σ~y~*, and
*τ~xy~* components at each evaluation point for all of the tires in a
wheel group. The total sum for each evaluation point is then converted
to the maximum horizontal stress at that evaluation point.

Horizontal strain is, currently, only required once per design as a
check for fatigue failure of the asphalt surface layer. This response is
therefore computed from the components of "all responses" as described
below (see equation 5-34 below).

Computation of vertical deflection is also required in a separate
backcalculation computer program. This is directly computed from the
integration of the elemental equation for vertical deflection.

## 5.7.2 All Responses.

The final step in any design in LEDFAA is to compute all structural
responses at each evaluation point for each aircraft wheel group and to
write these responses to a text file for reference by the designer. LEAF
therefore has an option to compute, and return, all stress, strain, and
deflection components in rectangular coordinates and all principal
stresses and strains.

First, *σ~x~*, *σ~y~*, *τ~xy~*, *τ~xz~*, *τ~yz~*, and *ε~z~* are
computed from the dedicated subroutines. The remaining stress and strain
components are then computed from the following expressions.

$\sigma_{z} = {E\varepsilon}_{z} + \upsilon(\sigma_{x} + \sigma_{y})$

$\varepsilon_{x} = \frac{1}{E}\left( \sigma_{x} - \upsilon(\sigma_{y} + \sigma_{z}) \right)$

$\varepsilon_{y} = \frac{1}{E}\left( \sigma_{y} - \upsilon(\sigma_{x} + \sigma_{z}) \right)$

$G = \frac{E}{2(1 + \upsilon)} = \text{Modulus of Rigidity}$

$$\varepsilon_{\text{xy}} = \frac{\tau_{\text{xy}}}{G}$$

$\varepsilon_{\text{xz}} = \frac{\tau_{\text{xz}}}{G}$

$\varepsilon_{\text{yz}} = \frac{\tau_{\text{yz}}}{G}$

Principal stresses and strains are computed using the procedures given
in Irons \[Error: Reference source not found\].

Vertical deflection, *δ~z~*, is computed using the dedicated subroutine.

$\begin{array}{r}
w = - \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{\frac{z_{1}}{\alpha_{1}}J_{0}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ A_{i} - C_{i}\left( 2 - 4\upsilon_{i} - \frac{{z\alpha}_{1}}{z_{1}} \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
\  - \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{\frac{z_{2}}{\alpha_{2}}J_{0}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ B_{i} + D_{i}\left( 2 - 4\upsilon_{i} + \frac{{z\alpha}_{2}}{z_{2}} \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}} \\
\text{and }\delta_{z} = w
\end{array}$ (5-34)

Horizontal deflection is computed with another dedicated subroutine
(finding horizontal deflection as a component in all of the responses is
its only use). Radial deflection, *u*, is computed first and then
transformed to rectangular coordinates.

$\begin{array}{r}
u = \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{\frac{z_{1}}{\alpha_{1}}J_{1}\left( \frac{{r\alpha}_{1}}{z_{1}} \right)J_{1}\left( \frac{{a\alpha}_{1}}{z_{1}} \right)\left\{ A_{i} + C_{i}\left( 1 + \frac{{z\alpha}_{1}}{z_{1}} \right) \right\} e^{- \alpha_{1}}\frac{{d\alpha}_{1}}{z_{1}}} \\
\  - \frac{\text{qa}(1 + \upsilon_{i})}{E_{i}}\int_{}^{}{\frac{z_{2}}{\alpha_{2}}J_{1}\left( \frac{{r\alpha}_{2}}{z_{2}} \right)J_{1}\left( \frac{{a\alpha}_{2}}{z_{2}} \right)\left\{ B_{i} - D_{i}\left( 1 - \frac{{z\alpha}_{2}}{z_{2}} \right) \right\} e^{- \alpha_{2}}\frac{{d\alpha}_{2}}{z_{2}}} \\
\text{and }\delta_{x} = u\text{cos}\theta \\
\delta_{y} = u\text{sin}\theta
\end{array}$ (5-35)

The equations for *w* and *u* both have *α* as a divisor inside the
integral sign. The equations for stress and strain do not. This probably
accounts for the better convergence of the displacement equations for
small values of *z~eval~* as demonstrated below.

# 5.8 Testing Against Other Methods.

Considering the wide range of load and structural conditions encountered
in airport pavement structures, it is very difficult to comprehensively
test a layered elastic computer program for accuracy and stability.
Tests are described in this section covering performance when simulating
a uniform half-space, with and without interface bonding, and when
simulating a new rigid pavement structure. A run time comparison with
JULEA is also described.

## 5.8.1 Response Comparisons.

The first concern is with convergence when the evaluation point is close
to the surface. The uniform pressure circular load is a discontinuous
function applied to the surface and causes slow decays in the transform
domain. (Rapid changes in a function with respect to radius have slow
changes with respect to the Hankel domain variable *α*. The same kind of
reciprocal behavior is also true of Laplace and Fourier transforms.)
JULEA uses an extrapolation technique to sidestep the poor convergence
when the evaluation point is close to the surface. However, the offsets
added in LEAF improve convergence under these circumstances and the
inaccuracies in response computation are not serious, particularly for
the computations made within the top layer for use in LEDFAA or for FWD
backcalculation. A contrived structure (table 5-7) representing a
uniform half-space is used to compare LEAF computations with Boussinesq
analytic solutions for responses under the load center.

TABLE 5-7.  STRUCTURE LAYER PROPERTIES USED TO ILLUSTRATE COMPARISON
BETWEEN BOUSSINESQ HALF-SPACE RESPONSES AND LEAF RESPONSES

  --------------- ------------- -------------- --------------- --------------
  Layer Number    Modulus, psi  Thickness, in  Poisson's Ratio Interface

  1               2,000          5             0.25            fully-bonded

  2               2,000         18             0.25            fully-bonded

  3               2,000          1             0.25            fully-bonded

  4               2,000         infinite       0.25            
  --------------- ------------- -------------- --------------- --------------

The load is 5 psi (34.5 kPa) on a radius of 10 inches (254 mm). The
results are shown in tables 5-8 through 5-10.

TABLE 5-8.  COMPARISON BETWEEN BOUSSINESQ AND LEAF -- VERTICAL
DEFLECTION

  ------------------ ----------------- ----------------- -----------------
  *z~eval~*          Boussinesq        LEAF              \% Difference

  0.0                0.0468750         0.0468960         0.045

  0.01               0.0468593         0.0468791         0.042

  0.1                0.0467190         0.0467262         0.018

  0.2                0.0465594         0.0465564         -0.0064

  0.5                0.0460743         0.0460744         0.00022

  1.0                0.0452353         0.0452355         0.00044
  ------------------ ----------------- ----------------- -----------------

TABLE 5-9.  COMPARISON BETWEEN BOUSSINESQ AND LEAF -- VERTICAL STRESS

  ------------------ ----------------- ----------------- -----------------
  *z~eval~*          Boussinesq        LEAF              \% Difference

  0.0                -5.00000          -5.38574          7.715

  0.01               -4.99999          -5.38065          7.613

  0.1                -4.99999          -5.18782          3.757

  0.2                -4.99996          -4.92312          -1.537

  0.5                -4.99938          -4.99854          -0.096

  1.0                -4.99507          -4.99507          0.0
  ------------------ ----------------- ----------------- -----------------

TABLE 5-10.  COMPARISON BETWEEN BOUSSINESQ AND LEAF -- HORIZONTAL STRESS

  ------------------ ----------------- ----------------- -----------------
  *z~eval~*          Boussinesq        LEAF              \% Difference

  0.0                -3.75000          -4.03931          7.715

  0.01               -3.74375          -3.95939          5.760

  0.1                -3.74375          -3.67988          -0.207

  0.2                -3.62504          -3.64225          0.475

  0.5                -3.43820          -3.43851          0.009

  1.0                -3.13056          -3.13056          0.0
  ------------------ ----------------- ----------------- -----------------

The accuracy of vertical deflection at the surface is very good and can
be used with confidence in FWD backcalculations. Comparison with JULEA
and BISAR also shows that vertical deflections computed at the surface
by LEAF are accurate to the same degree for arbitrary horizontal offsets
of the evaluation point.

Reasonably large errors are present for stresses computed at the
surface. Although these errors are less than 10 percent, it is not known
how the errors may change for different configurations of the structure.
The error is very small one half of an inch (13 mm) below the surface.
LEDFAA does not require the computation of any stress or strain response
in the upper half of the top layer and the inaccuracy is not of
significance for the LEDFAA application. FWD backcalculation requires
only vertical deflection at the surface, and this response is accurate
with LEAF.

Tables 5-11 through 5-13 show additional comparisons with Boussinesq and
JULEA for evaluation points deeper in the structure. Particular
attention is paid to performance close to the interfaces. The structure
and load are the same as above.

TABLE 5-11.  COMPARISON BETWEEN BOUSSINESQ, LEAF, AND JULEA, FULLY
BONDED INTERFACES, VERTICAL DISPLACEMENT, INCHES

  ------------------ ----------------- ----------------- -----------------
  *z~eval~*          Boussinesq        LEAF              JULEA

  0.00               0.0468750         0.0468960         0.0470

  4.99               0.0376274         0.0376276         0.0376

  5.01               0.0375878         0.0375880         0.0376

  10.00              0.0285692         0.0285693         0.0186

  22.99              0.0157157         0.0157159         0.0157

  23.01              0.0157041         0.0157042         0.0157

  23.99              0.0151497         0.0151497         0.0151

  24.01              0.0151387         0.0151388         0.0151
  ------------------ ----------------- ----------------- -----------------

TABLE 5-12.  COMPARISON BETWEEN BOUSSINESQ, LEAF, AND JULEA, FULLY
BONDED INTERFACES, HORIZONTAL STRESS, PSI

  ------------------ ------------------ ----------------- -----------------
  *z~eval~*          Boussinesq         LEAF              JULEA

  0.00               -3.750000          -4.03931          -3.75

  4.99               -1.18192           -1.18192          -1.18

  5.01               -1.17513           -1.17513          -1.16

  10.00              -0.214466          -0.214466         -0.215

  22.99              0.0535116          0.0535116         0.0534

  23.01              0.0535043          0.0535043         0.0534

  23.99              0.0529210          0.0529210         0.0527

  24.01              0.0529050          0.0529050         0.0527
  ------------------ ------------------ ----------------- -----------------

TABLE 5-13.  COMPARISON WITH BOUSSINESQ AND JULEA, FULLY BONDED
INTERFACES, VERTICAL STRAIN, MILLI-STRAIN

  ------------------ ------------------ ----------------- -----------------
  *z~eval~*          Boussinesq         LEAF              JULEA

  0.00               -1.56250           -1.68304          -1.57

  4.99               -1.98198           -1.98198          -1.98

  5.01               -1.98154           -1.98154          -2.00

  10.00              -1.56250           -1.56250          -1.56

  22.99              -0.585598          -0.585598         -0.585

  23.01              -0.584797          -0.584796         -0.584

  23.99              -0.547276          -0.547276         -0.546

  24.01              -0.546545          -0.546545         -0.546
  ------------------ ------------------ ----------------- -----------------

LEAF shows good performance for this test case.

Tables 5-14 through 5-16 show comparisons with JULEA for the first
interface (between the first and second layers) fully unbonded.
Otherwise the test case is the same as above.

TABLE 5-14.  COMPARISON WITH JULEA, FIRST INTERFACE FULLY UNBONDED,
VERTICAL DISPLACEMENT, INCHES

  ------------------------ ----------------------- -----------------------
  *z~eval~*                LEAF                    JULEA

  0.00                     0.0571609               0.0572

  4.99                     0.0452916               0.0453

  5.01                     0.0452456               0.0453

  10.00                    0.0359154               0.0359

  22.99                    0.0188570               0.0189

  23.01                    0.0188411               0.0188

  23.99                    0.0180932               0.0181

  24.01                    0.0180784               0.0181
  ------------------------ ----------------------- -----------------------

TABLE 5-15.  COMPARISON WITH JULEA, FIRST INTERFACE FULLY UNBONDED,
HORIZONTAL STRESS, PSI

  ------------------------ ----------------------- -----------------------
  *z~eval~*                LEAF                    JULEA

  0.00                     -3.28587                -3.02

  4.99                     1.64421                 1.63

  5.01                     -3.83583                -3.86

  10.00                    -1.06495                -1.07

  22.99                    0.0253571               0.0261

  23.01                    0.0255122               0.0262

  23.99                    0.0319912               0.0324

  24.01                    0.0321025               0.0325
  ------------------------ ----------------------- -----------------------

TABLE 5-16.  COMPARISON BETWEEN LEAF AND JULEA, FIRST INTERFACE FULLY
UNBONDED, VERTICAL STRAIN, MILLI-STRAIN

  ------------------------ ----------------------- -----------------------
  *z~eval~*                LEAF                    JULEA

  0.00                     -1.87157                -1.75

  4.99                     -2.97356                -2.96

  5.01                     -1.60355                -1.62

  10.00                    -1.90294                -1.90

  22.99                    -0.79165                -0.792

  23.01                    -0.79051                -0.791

  23.99                    -0.73696                -0.737

  24.01                    -0.73591                -0.736
  ------------------------ ----------------------- -----------------------

The comparison with JULEA is good except for stress and strain at the
surface. This error is ignored as discussed above.

The structure used in the above comparisons is representative of a
uniform half-space. Rigid pavement structures represent the extreme in
non-uniform layered systems. Table 5-17 gives the parameter values for a
typical rigid pavement structure. A 50,000 lb (222 kN) load with 180 psi
(1,240 kPa) contact pressure was applied to the structure. Vertical and
horizontal responses in the top (PCC) layer are shown in figures 5-8
through 5-11 for LEAF, JULEA, and Bisar.

TABLE 5-17.  TYPICAL RIGID PAVEMENT STRUCTURE USED TO ILLUSTRATE TOP
LAYER STRESS DISTRIBUTION

  -------------- ------------- ------------- -------------- ----------------
  Layer Number   Modulus, psi  Thickness, in Poisson's      Interface
                                             Ratio          

  1              4,000,000     14            0.15           fully-unbonded

  2              700,000        6            0.20           fully-bonded

  3              35,430         6            0.35           fully-bonded

  4              15,000        infinite      0.35           
  -------------- ------------- ------------- -------------- ----------------

![](media/image4.wmf)

FIGURE 5-8.  VERTICAL STRESS IN THE TOP LAYER OF A RIGID PAVEMENT
STRUCTURE (see figure 5-9 for legend)

![](media/image5.wmf)

FIGURE 5-9.  FIGURE 5-8 EXPANDED OVER THE TOP 2 INCHES

![](media/image6.wmf)

FIGURE 5-10.  HORIZONTAL STRESS IN THE TOP LAYER OF A RIGID PAVEMENT
STRUCTURE(see figure 5-11 for legend)

![](media/image7.wmf)

FIGURE 5-11.  FIGURE 5-10 EXPANDED OVER THE TOP 0.5 INCHES

## 5.8.2 Run time Comparison.

A run time comparison was made against JULEA using a modified version of
LEDFAA which has JULEA and leaf selectable as the layered elastic
computational routine. The version of leaf used in the comparison did
not have the additional offsets and run times were probably less than
they will be with the version of leaf described in this document.
However, the trends shown in the results are still valid. New runtime
tests will be made when the new version of leaf has been fully
integrated into LEDFAA.

Figure 5-12 shows complete run times for execution of the "Life" button
in LEDFAA for a set of test cases with 20 B-777 aircraft in the traffic
mix. The test cases included up to 16 layers in a flexible pavement
structure with all interfaces fully bonded. The version of JULEA used in
the comparison was in Fortran compiled as a Windows dynamic link
library. LEAF was fully integrated into the Visual Basic code for LEDFAA
and the complete program compiled into native code under VB 5.0. Times
shown are for total runtime, including all LEDFAA housekeeping and CDF
computation, etc.

![](media/image8.wmf)

FIGURE 5-12.  RUNTIME COMPARISON BETWEEN JULEA AND LEAF -- TIMES FOR
LEDFAA SOLUTIONS

General observations are:

1.  LEAF is more efficient than JULEA for all configurations tested
    (although it is not known in detail what influence the Fortran
    versus Visual Basic implementations had on runtimes).

2.  Runtimes for JULEA increase significantly as the number of layers
    increase.

3.  Part inversion gives significant benefit only with large numbers of
    layers.

[^1]: . Irons, Bruce, and Ahmad, Sohrab, Techniques of Finite Elements,
    Prentice Hall, 1981.