under construction ...

Elastic Layer Analysis (ELA)

Introduction

The following is an example of SymPy in action. This example demonstrates the use of various SymPy functions, specifically the use of matrix algebra and integration techniques to solve an engineering problem i.e. the analysis of the stress, strain and displacement response of an elastic layered structure to imposed surface loads. Typically these problems are solved using finite element methods (FEM) but a simple numerical solution based on the Hankel transformation of a stress function has been in use since it was introduced by Burmister¹.

Elastic layer analysis (ELA) is routinely used in pavement engineering for the design and evaluation of road and airport structures. As with any structure, roads are designed and constructed to carry surface loads - the heavier the load, the stronger the structure needed. Concrete pavements on interstate highways, for example, are durable structures that can withstand very high compressive loads but a thin layer of concrete has a relatively low flexural strength and must therefore be supported to prevent excessive bending. This is typically achieved by paving the concrete layer of top of a cemented or stabilized base layer. The base layer in turn may be placed on top of a strong granular subbase layer, which in turn covers the foundation or subgrade, which is usually the weakest layer in the structure. Building stronger layers on top of the subgrade serves to protect it from excessive deformation. Figure 1 illustrates a typical pavement structure.

Figure 1. Layered pavement structure

For the elastic analysis of the pavement structure shown in Figure 1, the different layers are assumed to be homogeneous and isotropic. The materials used in each layer are defined by elastic parameters i.e. Young's modulus (( E )) and Poisson's ratio (( nu )). Each layer has a finite thickness except the lowest subgrade layer, which is assumed to have an infinite thickness. The nature of the interface between successive layers dictates how stresses and displacements are transferred from one layer to another. For the purpose of this example, full-friction between the layers is assumed, therefore shear stresses and lateral displacements at the lower interface of any layer are equal to those at the top interface of the underlying layer.

Problem Formulation

Stresses and displacements written to show the Hankel form are: [sigma{z}=int{0}^{infty}xi^{2}left(left(-Axi+Cleft(1-2nu-xi zright)right)e^{xi z}+left(Bxi+Dleft(1-2nu+xi zright)right)e^{-xi z}right)\mathrm{J}_0left(rxiright)xi,dxi] [tau{rz}=int{0}^{infty}xi^{2}frac{\mathrm{J}_1left(rxiright)}{\mathrm{J}_0left(rxiright)}left(left(Axi+Cleft(2nu+xi zright)right)e^{xi z}+left(Bxi+Dleft(1-2nu+xi zright)right)e^{-xi z}right)\mathrm{J}_0left(rxiright)xi,dxi] [w=frac{1+nu}{E}int{0}^{infty}xileft(left(-Axi+Cleft(2-4nu-xi zright)right)e^{xi z}-left(Bxi+Dleft(2-4nu+xi zright)right)e^{- xi z}right)\mathrm{J}_0left(rxiright)xi,dxi] [u=frac{1+nu}{E}int{0}^{infty}xifrac{\mathrm{J}_1left(rxiright)}{\mathrm{J}_0left(rxiright)}left(left(Axi+Cleft(1+xi zright)right)e^{xi z}-left(Bxi-Dleft(1-xi zright)right)e^{-xi z}right)\mathrm{J}_0left(rxiright)xi,dxi]

(sigma{z}), (tau{rz}): vertical compressive and shear stress
(u), (w): lateral and vertical displacement
(\mathrm{J}_0, \mathrm{J}_1): Bessel functions of the first kind of orders 0 and 1
(xi), (r), (z): Hankel parameter, radial position in cylindrical coordinates and vertical position
(A,B,C,D): integration constants

The unknowns in the stress and displacement equations are the four integration constants. These may be determined by applying boundary conditions.

Boundary Conditions

The boundaries of the elastic layer system include the surface, interfaces between the layers and the infinite lowest layer.

Surface. For the purpose of this example it is assumed that the loads imposed on the system are circular surface loads defined by a radius (a) and a uniform pressure (p), the circular imprint simulates the impression of vehicle tyres as shown in Figure 2.

Figure 2. Load imprint

If only a vertical uniformly distributed circular load is considered, at a radial position (r) on the surface: [ rleq a left{ begin{array}{l l} sigma{z}left(r,0right)=-p\tau{rz}=0\end{array}right. ] [ rgt a left{ begin{array}{l l} sigma{z}left(r,0right)=0\tau{rz}=0\end{array}right. ]

In the Hankel domain, the load may be determined as: [q=int{0}^{infty}rp\mathrm{J}_0left(xi rright),dr=frac{pa}{xi}\mathrm{J}_1left(xi aright)]

Interlayers. As indicated previously, with full friction, to ensure equilibrium, stresses and displacements at the interlayers are equal: [sigma{z}^{i}left(r,h_{i}right)=sigma{z}^{i+1}left(r,0right)] [tau{rz}^{i}left(r,h_{i}right)=tau{rz}^{i+1}left(r,0right)] [u{r}^{i}left(r,h_{i}right)=u_{r}^{i+1}left(r,0right)] [u{z}^{i}left(r,h_{i}right)=u_{z}^{i+1}left(r,0right)]

Infinite lower layer. All stresses and displacements approach zero at infinite depth ((zrightarrowinfty)). It can be shown that this implies that the integration constants A and C for this layer are zero.

Applying the boundary conditions allows the integration constants to be solved. Depending on the number of layers (n) in the system, the number of equations to solve is 4n-2. The solved constants may then be substituted into the stress and displacement equations, which are integrated to determine response to surface loads at specific evaluation positions r and z. To demonstrate the process, consider a single layer subjected to a single circular load.

One layer solution

Applying the problem formulation and boundary conditions outlined above, the following code listing shows how SymPy can be used to determine the variation in normal stress with depth for a single layer under a circular load.

from sympy import Symbol,Matrix,exp,lambdify
from sympy.mpmath import sqrt,quad,quadosc,inf,pi
from time import time

"""
one.py

Calculate the stress variation with depth beneath
a single layer subjected to a circular load.

"""

xi=Symbol('xi')     # Hankel parameter
j0=Symbol('j0')     # Bessel function J0
j1=Symbol('j1')     # Bessel function J1
nu0=Symbol('nu0')   # Poisson's ratio of layer
a=Symbol('a')       # Load radius
z=Symbol('z')       # Vertical evaluation position

q=j1(xi*a)          # Hankel transformed load

def hsz(z,nu):
    """
    Constants for Hankel transformed normal stress
    
    """ 
    a1=exp(xi*z)
    b1=-exp(-xi*z)
    c1=-(1-2*nu-xi*z)*exp(xi*z)
    d1=-(1-2*nu+xi*z)*exp(-xi*z)
    return a1,b1,c1,d1

def htrz(z,nu):
    """
    Constants for Hankel transformed shear stress
    
    """ 
    a2=exp(xi*z)
    b2=exp(-xi*z)
    c2=(2*nu+xi*z)*exp(xi*z)
    d2=-(2*nu-xi*z)*exp(-xi*z)
    return a2,b2,c2,d2

def s_z(r,z,nu,A,B,C,D):
    """
    Normal stress function
    
    """
    if r==0:
        return (A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z)
    else:
        return ((A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z))*j0(r*xi)

# Apply boundary condition for surface load (z=0)
a1,b1,c1,d1=hsz(0,nu0)
a2,b2,c2,d2=htrz(0,nu0)

# For this single infinite layer the integration constants A and C are zero.
P=Matrix([[b1,d1],[b2,d2]])
R=Matrix(2,1,[1,0])

# 2 equations, 2 unknowns - solve the integration constants
B0,D0=P.LUsolve(R)

# Substitute constants into stress equation (s_z)
# Determine the stress as a function of z beneath the load (r=0). Note: A=C=0
sz=s_z(0,z,nu0,0,B0,0,D0)*q

Note the application of SymPy symbols as well as the SymPy symbolic functions j0 and j1. The function sz in the above listing is the integrand of the normal stress with depth (z). Integrating this function provides the general solution: [ sigma{z}=int{0}^{infty}left(1+xi zright)e^{-xi z}\mathrm{J}_1left(a xiright),dxi ]

Solving the integration constants

For the infinite layer, the integration constants A and C are zero. To solve for the other constants, B and D, the boundary conditions on the surface are applied i.e. Matrix P. To simplify the solution of the integration constants, the response matrix, R, is set to Matrix(2,1,[1,0]) instead of Matrix(2,1,[q,1]) and the transformed load is applied just before integration by multiplying sz by q. Thus we have two equations and two unknowns. Use is made of the LUsolve function to solve for B and D: [ begin{bmatrix} -1 & -1+2nu \ 1 & -2nu end{bmatrix} begin{bmatrix} B \ D end{bmatrix} = begin{bmatrix} 1 \0 end{bmatrix} ]

Instead of using the LUsolve function, an alternative approach would be to multiply the response matrix R by the inverse of the global matrix P:

B0,D0=P.inv()*R

Integrating the stress function

Sympy does not currently include Bessel functions of the first kind, therefore to solve the stress function it is necessary to firstly lambdify [provide help link???] the integrand (see the code segment below). This converts all Sympy functions to corresponding mpmath functions. This was the reason for defining the SymPy functions j0 and j1, which when lambdified are converted into the corresponding mpmath functions j0 and j1.

The integrand sz is a function of the Bessel of order 1 and an exponential function. On the surface (z=0), the exponential function is eliminated, and the integrand is simply a function of the Bessel, which is very highly oscillatory. Figure 3 shows the behaviour of sz at the surface and at a depth of 100 mm. Clearly, the exponential serves to dampen the Bessel oscillations.

Figure 3. Oscillating and damped stress function

Figure 3 emphasises the problem of determining a solution for the stress function at (or near) the surface of the structure. Integration techniques such as the mpmath quad will fail and it becomes necessary to apply integration methods that account for the oscillatory behaviour. One option is to use the mpmath integration function quadosc. For application of quadosc it is necessary to estimate the period of the integrand. This is demonstrated by tagging the following code segment onto that above and running through python to integrate sz. Note how the integrand is multiplied by p and a following the integration to determine the vertical normal stress at the respective depths.

p=0.69              # Load pressure (MPa)
F=20000.            # Load force (Newton)
aa=sqrt(F/(p*pi))   # Load radius (mm)

print "\nIntegration method 1 (z=0)"
t = time()
f=sz
f=f.subs(z,0)
f=f.subs(a,aa)
f=lambdify(xi,f,"mpmath")
szz=-p*aa*quad(f,[0,inf],verbose=True)
print "szz: ",szz
print "Time: ",time()-t

print "\nIntegration method 2 (z=0)"
t = time()
f=sz
f=f.subs(z,0)
f=f.subs(a,aa)
f=lambdify(xi,f,"mpmath")
szz=-p*aa*quadosc(f,[0,inf],period=pi)
print "szz: ",szz
print "Time: ", time()-t

print "\nIntegration method 3 (z=0)"
t = time()
f=sz
f=f.subs(z,0)
f=f.subs(a,aa)
f=lambdify(xi,f,"mpmath")
szz=-p*aa*quadosc(f,[0,inf],period=2*pi/aa)
print "szz: ",szz
print "Time: ", time()-t

print "\nIntegration method 1 (z=100)"
t = time()
f=sz
f=f.subs(z,100)
f=f.subs(a,aa)
f=lambdify(xi,f,"mpmath")
szz=-p*aa*quad(f,[0,inf],verbose=True)
print "szz: ",szz
print "Time: ",time()-t

print "\nIntegration method 3 (z=100)"
t = time()
f=sz
f=f.subs(z,0)
f=f.subs(a,aa)
f=lambdify(xi,f,"mpmath")
szz=-p*aa*quadosc(f,[0,inf],period=2*pi/aa)
print "szz: ",szz
print "Time: ", time()-t

The output:

Integration method 1 (z=0)
Integrating from 0 to +inf (degree 1 of 6)
Integrating from 0 to +inf (degree 2 of 6)
Estimated error: 4.31832e+8
Integrating from 0 to +inf (degree 3 of 6)
Estimated error: 1.0
Integrating from 0 to +inf (degree 4 of 6)
Estimated error: 1.0
Integrating from 0 to +inf (degree 5 of 6)
Estimated error: 1.0
Integrating from 0 to +inf (degree 6 of 6)
Estimated error: 1.0
Failed to reach full accuracy. Estimated error: 1.0
szz:  8058863524.695
Time:  1.09746789932

Integration method 2 (z=0)
szz:  -0.689999990998652
Time:  289.689537048

Integration method 3 (z=0)
szz:  -0.69
Time:  2.73873400688

Integration method 1 (z=100)
Integrating from 0 to +inf (degree 1 of 6)
Integrating from 0 to +inf (degree 2 of 6)
Estimated error: 0.00217811
Integrating from 0 to +inf (degree 3 of 6)
Estimated error: 0.001
Integrating from 0 to +inf (degree 4 of 6)
Estimated error: 1.0e-9
Integrating from 0 to +inf (degree 5 of 6)
szz:  -0.431176903899328
Time:  0.585129022598

Integration method 3 (z=100)
szz:  -0.431176903899328
Time:  0.899243116379

From the output it can be seen that quad fails to integrate the sz function at the surface but quadosc with the period correctly set is very efficient. With the period not set correctly, quadosc provides an accurate solution, but after a long wait. At depths greater than 0, quad provides a rapid and accurate solution of the stress function. In fact, at greater depths, quad will provide a more accurate solution than quadosc. Note the use of the SymPy subs function to substitute numerical values for the symbolic expressions.

To test the solution, the results may be compared to the Boussinesq case² for an infinite half space with a circular load at the origin. The solution for vertical stress beneath the load (r=0): [ sigma{z}=pleft(frac{z^{3}}{left(z^{2}+a^{2}right)^{3/2}}-1right) ]

Two layer solution

With two layers (n=2) the number of integration constants is 4n-2 = 6 i.e. A0,B0,C0,D0 for the first layer and B1,D1 for the infinitely thick second layer with A1 = C1 = 0. In addition to applying the boundary conditions at the surface, the equality of stresses and displacements at the interface between the two layers provides a solution for the integration constants.

On the surface the normal stress is the applied load and the shear stress is zero: [ begin{bmatrix} 1 & -1 & -1 + 2 nu{0} & -1 + 2 nu{0}\1 & 1 & 2 nu{0} & - 2 nu{0} end{bmatrix} begin{bmatrix} A0 \ B0 \C0 \D0 \end{bmatrix} = begin{bmatrix} 1 \0 end{bmatrix} ]

At the interface (z=(h{0})): [ begin{bmatrix}e^{h{0} xi} & - e^{- h{0} xi} & - left(1 - 2 nu{0} - h{0} xiright) e^{h{0} xi} & - left(1 - 2 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & e^{- h{0} xi} & left(2 nu{0} + h{0} xiright) e^{h{0} xi} & left(- 2 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & e^{- h{0} xi} & - left(2 - 4 nu{0} - h{0} xiright) e^{h{0} xi} & left(2 - 4 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & - e^{- h{0} xi} & left(1 + h{0} xiright) e^{h{0} xi} & left(1 - h{0} xiright) e^{- h{0} xi}end{bmatrix}begin{bmatrix} A0 \ B0 \C0 \D0 \end{bmatrix} = begin{bmatrix}-1 & -1 + 2 nu{1}\-1 & -1 + 2 nu{1}\frac{mu{1}}{mu{0}} & frac{mu{1} left(2 - 4 nu{1}right)}{mu{0}}\- frac{mu{1}}{mu{0}} & frac{mu{1}}{mu{0}}end{bmatrix}begin{bmatrix} B1 \D1 \end{bmatrix} ]

Bringing these together into a single matrix: [ begin{bmatrix}e^{h{0} xi} & - e^{- h{0} xi} & - left(1 - 2 nu{0} - h{0} xiright) e^{h{0} xi} & - left(1 - 2 nu{0} + h{0} xiright) e^{- h{0} xi} & 1 & 1 - 2 nu{1}\e^{h{0} xi} & e^{- h{0} xi} & left(2 nu{0} + h{0} xiright) e^{h{0} xi} & left(- 2 nu{0} + h{0} xiright) e^{- h{0} xi} & 1 & 1 - 2 nu{1}\e^{h{0} xi} & e^{- h{0} xi} & - left(2 - 4 nu{0} - h{0} xiright) e^{h{0} xi} & left(2 - 4 nu{0} + h{0} xiright) e^{- h{0} xi} & - frac{mu{1}}{mu{0}} & - frac{mu{1} left(2 - 4 nu{1}right)}{mu{0}}\e^{h{0} xi} & - e^{- h{0} xi} & left(1 + h{0} xiright) e^{h{0} xi} & left(1 - h{0} xiright) e^{- h{0} xi} & frac{mu{1}}{mu{0}} & - frac{mu{1}}{mu{0}}\1 & -1 & -1 + 2 nu{0} & -1 + 2 nu{0} & 0 & 0\1 & 1 & 2 nu{0} & - 2 nu{0} & 0 & 0end{bmatrix} begin{bmatrix} A0 \ B0 \C0 \D0 \B1 \D1 end{bmatrix} = begin{bmatrix} 0 \0 \0 \0 \1 \0 end{bmatrix} ]

The code listing below shows the derivation of these matrices and the solution for the 6 integration constants:

from sympy import Symbol,Matrix,exp,lambdify
from sympy.mpmath import sqrt,quad,quadosc,inf,pi
from time import time

"""
two_a.py

Calculate the stress variation with depth beneath
a two layered system subjected to a circular load.

"""

xi=Symbol('xi')     # Hankel parameter
j0=Symbol('j0')     # Bessel function J0
j1=Symbol('j1')     # Bessel function J1
nu0=Symbol('nu0')   # Poisson's ratio of layer 1
nu1=Symbol('nu1')   # Poisson's ratio of layer 2
h0=Symbol('h0')     # Thickness of layer 1
mu0=Symbol('mu0')   # Lame's parameter for layer 1
mu1=Symbol('mu1')   # Lame's parameter for layer 2
a=Symbol('a')       # Load radius
z=Symbol('z')       # Vertical evaluation position

q=j1(xi*a)          # Hankel transformed load

def hsz(z,nu):
    """
    Constants for Hankel transformed normal stress
    
    """ 
    a1=exp(xi*z)
    b1=-exp(-xi*z)
    c1=-(1-2*nu-xi*z)*exp(xi*z)
    d1=-(1-2*nu+xi*z)*exp(-xi*z)
    return a1,b1,c1,d1

def htrz(z,nu):
    """
    Constants for Hankel transformed shear stress
    
    """ 
    a2=exp(xi*z)
    b2=exp(-xi*z)
    c2=(2*nu+xi*z)*exp(xi*z)
    d2=-(2*nu-xi*z)*exp(-xi*z)
    return a2,b2,c2,d2

def hur(z,nu,R):
    """
    Constants for Hankel transformed radial displacement
    
    """ 
    a3=exp(xi*z)*R
    b3=exp(-xi*z)*R
    c3=-(2-4*nu-xi*z)*exp(xi*z)*R
    d3=(2-4*nu+xi*z)*exp(-xi*z)*R
    return a3,b3,c3,d3

def huz(z,nu,X):
    """
    Constants for Hankel transformed vertical displacement
    
    """ 
    a4=exp(xi*z)*X
    b4=-exp(-xi*z)*X
    c4=(1+xi*z)*exp(xi*z)*X
    d4=(1-xi*z)*exp(-xi*z)*X
    return a4,b4,c4,d4

def s_z(r,z,nu,A,B,C,D):
    """
    Normal stress function
    
    """
    if r==0:
        return (A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z)
    else:
        return ((A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z))*j0(r*xi)

# Apply boundary condition for surface load (z=0)
a00,b00,c00,d00=hsz(0,nu0)
a01,b01,c01,d01=htrz(0,nu0)

# Boundary condition at upper interface
a02,b02,c02,d02=hsz(h0,nu0)
a03,b03,c03,d03=htrz(h0,nu0)
a04,b04,c04,d04=hur(h0,nu0,1)
a05,b05,c05,d05=huz(h0,nu0,1)

# Boundary condition at lower interface
X=mu1/mu0
a06,b06,c06,d06=hsz(0,nu1)
a07,b07,c07,d07=htrz(0,nu1)
a08,b08,c08,d08=hur(0,nu1,X)
a09,b09,c09,d09=huz(0,nu1,X)

P=Matrix([[a02,b02,c02,d02,-b06,-d06],
          [a03,b03,c03,d03,-b06,-d06],
          [a04,b04,c04,d04,-b08,-d08],
          [a05,b05,c05,d05,-b09,-d09],
          [a00,b00,c00,d00,   0,   0],
          [a01,b01,c01,d01,   0,   0]])

R=Matrix(6,1,[0,0,0,0,1,0])

# 6 equations, 6 unknowns - solve the integration constants
A0,B0,C0,D0,B1,D1=P.LUsolve(R)

In the above listing, two elastic parameters (mu_0) and (mu_1) are defined and the ratio of these is evident in the matrices when equating the lateral and vertical displacements at the interlayer. Save the code listing above and import it into isympy:

In [1]: from two_a import *

In [2]: A0
*Too long to display*   

In [3]: len(str(A0))
101784

The integration constant expressions are huge - and this just for a two layer system. Imagine the size of these for 10 or more layers! Of course these constants can be substituted into the stress function as is and integrated to obtain a solution for the two layer system, but given the size of the constants, it will take forever to integrate. It becomes necessary, therefore, to reduce the complexity and size of the integration constants for a faster solution time. This simplication can be achieved through manipulation of the equilibrium equations.

Simplifying the solution

SymPy includes a number of options to simplify expressions. Given the size of the integration constants, however, simplication of these may be very slow. A better approach is to rewrite the equilibrium equations. For example, multiplying each element in a single column in a matrix by a given value divides the value of the corresponding coefficient (A, B, C, or D) by the same amount. This provides a means to eliminate some of the exponential expressions in the global matrix i.e. multiply columns 1 and 3 by (e^{-xi h_0}) and columns 2 and 4 by (e^{xi h_0}) since (e^{xi h_0}e^{-xi h_0}=1). Reconsider the matrix resulting from the boundary conditions at the interface:

[ begin{bmatrix}e^{h{0} xi} & - e^{- h{0} xi} & - left(1 - 2 nu{0} - h{0} xiright) e^{h{0} xi} & - left(1 - 2 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & e^{- h{0} xi} & left(2 nu{0} + h{0} xiright) e^{h{0} xi} & left(- 2 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & e^{- h{0} xi} & - left(2 - 4 nu{0} - h{0} xiright) e^{h{0} xi} & left(2 - 4 nu{0} + h{0} xiright) e^{- h{0} xi}\e^{h{0} xi} & - e^{- h{0} xi} & left(1 + h{0} xiright) e^{h{0} xi} & left(1 - h{0} xiright) e^{- h{0} xi}end{bmatrix}begin{bmatrix} A0 \ B0 \C0 \D0 \end{bmatrix} = begin{bmatrix}-1 & -1 + 2 nu{1}\-1 & -1 + 2 nu{1}\frac{mu{1}}{mu{0}} & frac{mu{1} left(2 - 4 nu{1}right)}{mu{0}}\- frac{mu{1}}{mu{0}} & frac{mu{1}}{mu{0}}end{bmatrix}begin{bmatrix} B1 \D1 \end{bmatrix} ]

Eliminating the exponentials:

[ begin{bmatrix}1 & -1 & -1 + 2 nu{0} + h{0} xi & -1 + 2 nu{0} - h{0} xi\1 & 1 & 2 nu{0} + h{0} xi & - 2 nu{0} + h{0} xi\1 & 1 & -2 + 4 nu{0} + h{0} xi & 2 - 4 nu{0} + h{0} xi\1 & -1 & 1 + h{0} xi & 1 - h{0} xiend{bmatrix}begin{bmatrix} A'0 \ B'0 \C'0 \D'0 end{bmatrix} = begin{bmatrix}-1 & -1 + 2 nu{1}\-1 & -1 + 2 nu{1}\frac{mu{1}}{mu{0}} & frac{mu{1} left(2 - 4 nu{1}right)}{mu{0}}\- frac{mu{1}}{mu{0}} & frac{mu{1}}{mu{0}}end{bmatrix}begin{bmatrix} B1 \D1 \end{bmatrix} ]

where (A'0=A0/e^{xi h_0}), etc.

The inverse of the matrix on the left multiplied by the matrix on the right hand side provides a simple solution for the new integration constants A'0, B'0, C'0 and D'0:

[ begin{bmatrix} 1 & 0 & 0 & 0\0 & 1 & 0 & 0\0 & 0 & 1 & 0\0 & 0 & 0 & 1 end{bmatrix}begin{bmatrix} A'0 \ B'0 \C'0 \D'0 end{bmatrix} = left[LHSright]^{-1}left[RHSright] ]

The following code listing applies these manipulations to the above code:

from sympy import Symbol,Matrix,exp,lambdify
from sympy.mpmath import sqrt,quad,quadosc,inf,pi
from time import time

"""
two_b.py

Calculate the stress variation with depth beneath
a two layered system subjected to a circular load.

"""

xi=Symbol('xi')     # Hankel parameter
j0=Symbol('j0')     # Bessel function J0
j1=Symbol('j1')     # Bessel function J1
nu0=Symbol('nu0')   # Poisson's ratio of layer 1
nu1=Symbol('nu1')   # Poisson's ratio of layer 2
h0=Symbol('h0')     # Thickness of layer 1
mu0=Symbol('mu0')   # Lame's parameter for layer 1
mu1=Symbol('mu1')   # Lame's parameter for layer 2
a=Symbol('a')       # Load radius
z=Symbol('z')       # Vertical evaluation position

q=j1(xi*a)          # Hankel transformed load

def hsz(z,nu):
    """
    Constants for Hankel transformed normal stress
    
    """ 
    a1=exp(xi*z)
    b1=-exp(-xi*z)
    c1=-(1-2*nu-xi*z)*exp(xi*z)
    d1=-(1-2*nu+xi*z)*exp(-xi*z)
    return a1,b1,c1,d1

def htrz(z,nu):
    """
    Constants for Hankel transformed shear stress
    
    """ 
    a2=exp(xi*z)
    b2=exp(-xi*z)
    c2=(2*nu+xi*z)*exp(xi*z)
    d2=-(2*nu-xi*z)*exp(-xi*z)
    return a2,b2,c2,d2

def hur(z,nu,R):
    """
    Constants for Hankel transformed radial displacement
    
    """ 
    a3=exp(xi*z)*R
    b3=exp(-xi*z)*R
    c3=-(2-4*nu-xi*z)*exp(xi*z)*R
    d3=(2-4*nu+xi*z)*exp(-xi*z)*R
    return a3,b3,c3,d3

def huz(z,nu,X):
    """
    Constants for Hankel transformed vertical displacement
    
    """ 
    a4=exp(xi*z)*X
    b4=-exp(-xi*z)*X
    c4=(1+xi*z)*exp(xi*z)*X
    d4=(1-xi*z)*exp(-xi*z)*X
    return a4,b4,c4,d4

def s_z(r,z,nu,A,B,C,D):
    """
    Normal stress function
    
    """
    if r==0:
        return (A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z)
    else:
        return ((A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z))*j0(r*xi)

# Apply boundary condition for surface load (z=0)
a00,b00,c00,d00=hsz(0,nu0)
a01,b01,c01,d01=htrz(0,nu0)

# Boundary condition at upper interface
a02,b02,c02,d02=hsz(h0,nu0)
a03,b03,c03,d03=htrz(h0,nu0)
a04,b04,c04,d04=hur(h0,nu0,1)
a05,b05,c05,d05=huz(h0,nu0,1)

# Boundary condition at lower interface
X=mu1/mu0
a06,b06,c06,d06=hsz(0,nu1)
a07,b07,c07,d07=htrz(0,nu1)
a08,b08,c08,d08=hur(0,nu1,X)
a09,b09,c09,d09=huz(0,nu1,X)

M11=Matrix([[a02,b02,c02,d02],
            [a03,b03,c03,d03],
            [a04,b04,c04,d04],
            [a05,b05,c05,d05]])

M11.col(0,lambda i,j: i*exp(-xi*h0))
M11.col(1,lambda i,j: i*exp(xi*h0))
M11.col(2,lambda i,j: i*exp(-xi*h0))
M11.col(3,lambda i,j: i*exp(xi*h0))

M12=Matrix([[a06,b06,c06,d06],
            [a07,b07,c07,d07],
            [a08,b08,c08,d08],
            [a09,b09,c09,d09]])

N12=M11.LUsolve(M12)
N12.col_del(0)
N12.col_del(1)

M21=Matrix([[a00,b00,c00,d00],
            [a01,b01,c01,d01]])

M21.col(0,lambda i,j: i*exp(-xi*h0))
M21.col(1,lambda i,j: i*exp(xi*h0))
M21.col(2,lambda i,j: i*exp(-xi*h0))
M21.col(3,lambda i,j: i*exp(xi*h0))

P=Matrix([[1,0,0,0,-N12[0,0],-N12[0,1]],
          [0,1,0,0,-N12[1,0],-N12[1,1]],
          [0,0,1,0,-N12[2,0],-N12[2,1]],
          [0,0,0,1,-N12[3,0],-N12[3,1]],
          [M21[0,0],M21[0,1],M21[0,2],M21[0,3],0,0],
          [M21[1,0],M21[1,1],M21[1,2],M21[1,3],0,0]])

R=Matrix(6,1,[0,0,0,0,1,0])

A0,B0,C0,D0,B1,D1=P.LUsolve(R)

A0=A0/exp(xi*h0)
B0=B0/exp(-xi*h0)
C0=C0/exp(xi*h0)
D0=D0/exp(-xi*h0)

Save the code listing above and import it into isympy:

In [1]: from two_b import *

In [2]: A0
*Too long to display*   

In [3]: len(str(A0))
6305

This has significantly reduced the complexity of the integration constants. Let's apply these simplications to a three layer system with two interfaces.

Three layer solution

With three layers there are two interfaces and 10 equations to solve for the integration constants. Using the approaches as outlined above, the code listing that follows shows how to solve for stresses beneath the first layer.

from numpy import arange
import matplotlib.pyplot as plt
from sympy import Symbol,Matrix,lambdify,N,exp
from mpmath import sqrt,quad,quadosc,inf,pi,j1
from time import time

"""
three_a.py

Calculate the stress variation with depth beneath
a three layered system subjected to a circular load.

"""

xi=Symbol('xi')     # Hankel parameter
j0=Symbol('j0')     # Bessel function J0
j1=Symbol('j1')     # Bessel function J1

p=0.69              # Load pressure (MPa)
F=20000.            # Load force (Newton)
a=sqrt(F/(p*pi))    # Load radius (mm)

q=j1(xi*a)          # Hankel transformed load

h=[]
h.append(50.)
h.append(150.)
nu=[]
nu.append(0.45)
nu.append(0.35)
nu.append(0.25)
E=[]
E.append(2000.)
E.append(1000.)
E.append(500.)
mu=[]
for i,j in zip(nu,E):
    mu.append((1+i)/j) 

def hsz(z,nu):
    """
    Constants for Hankel transformed normal stress
    
    """ 
    a1=exp(xi*z)
    b1=-exp(-xi*z)
    c1=-(1-2*nu-xi*z)*exp(xi*z)
    d1=-(1-2*nu+xi*z)*exp(-xi*z)
    return a1,b1,c1,d1

def htrz(z,nu):
    """
    Constants for Hankel transformed shear stress
    
    """ 
    a2=exp(xi*z)
    b2=exp(-xi*z)
    c2=(2*nu+xi*z)*exp(xi*z)
    d2=-(2*nu-xi*z)*exp(-xi*z)
    return a2,b2,c2,d2

def hur(z,nu,X):
    """
    Constants for Hankel transformed radial displacement
    
    """ 
    a3=exp(xi*z)*X
    b3=exp(-xi*z)*X
    c3=-(2-4*nu-xi*z)*exp(xi*z)*X
    d3=(2-4*nu+xi*z)*exp(-xi*z)*X
    return a3,b3,c3,d3

def huz(z,nu,X):
    """
    Constants for Hankel transformed vertical displacement
    
    """ 
    a4=exp(xi*z)*X
    b4=-exp(-xi*z)*X
    c4=(1+xi*z)*exp(xi*z)*X
    d4=(1-xi*z)*exp(-xi*z)*X
    return a4,b4,c4,d4

def s_z(r,z,nu,A,B,C,D):
    """
    Normal stress function
    
    """
    if r==0:
        return (A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z)
    else:
        return ((A-C*(1-2*nu-xi*z))*exp(xi*z)-(B+D*(1-2*nu+xi*z))*exp(-xi*z))*j0(r*xi)


print "Assemble global matrix"
# Apply boundary condition for surface load (z=0)
a00,b00,c00,d00=hsz(0,nu[0])
a01,b01,c01,d01=htrz(0,nu[0])

# Boundary condition at first upper interface
a02,b02,c02,d02=hsz(h[0],nu[0])
a03,b03,c03,d03=htrz(h[0],nu[0])
a04,b04,c04,d04=hur(h[0],nu[0],1)
a05,b05,c05,d05=huz(h[0],nu[0],1)

# Boundary condition at first lower interface
X=mu[1]/mu[0]
a06,b06,c06,d06=hsz(0,nu[1])
a07,b07,c07,d07=htrz(0,nu[1])
a08,b08,c08,d08=hur(0,nu[1],X)
a09,b09,c09,d09=huz(0,nu[1],X)

# Boundary condition at second upper interface
a10,b10,c10,d10=hsz(h[1],nu[1])
a11,b11,c11,d11=htrz(h[1],nu[1])
a12,b12,c12,d12=hur(h[1],nu[1],1)
a13,b13,c13,d13=huz(h[1],nu[1],1)

# Boundary condition at second lower interface
X=mu[2]/mu[1]
a14,b14,c14,d14=hsz(0,nu[2])
a15,b15,c15,d15=htrz(0,nu[2])
a16,b16,c16,d16=hur(0,nu[2],X)
a17,b17,c17,d17=huz(0,nu[2],X)

# Manipulate first layer
M11=Matrix([[a02,b02,c02,d02],
            [a03,b03,c03,d03],
            [a04,b04,c04,d04],
            [a05,b05,c05,d05]])

M11.col(0,lambda i,j: i*exp(-xi*h[0]))
M11.col(1,lambda i,j: i*exp(xi*h[0]))
M11.col(2,lambda i,j: i*exp(-xi*h[0]))
M11.col(3,lambda i,j: i*exp(xi*h[0]))

M12=Matrix([[a06,b06,c06,d06],
            [a07,b07,c07,d07],
            [a08,b08,c08,d08],
            [a09,b09,c09,d09]])

M12.col(0,lambda i,j: i*exp(-xi*h[1]))
M12.col(1,lambda i,j: i*exp(xi*h[1]))
M12.col(2,lambda i,j: i*exp(-xi*h[1]))
M12.col(3,lambda i,j: i*exp(xi*h[1]))

N1=M11.LUsolve(M12)

# Manipulate second layer
M21=Matrix([[a10,b10,c10,d10],
            [a11,b11,c11,d11],
            [a12,b12,c12,d12],
            [a13,b13,c13,d13]])

M21.col(0,lambda i,j: i*exp(-xi*h[1]))
M21.col(1,lambda i,j: i*exp(xi*h[1]))
M21.col(2,lambda i,j: i*exp(-xi*h[1]))
M21.col(3,lambda i,j: i*exp(xi*h[1]))

M22=Matrix([[a14,b14,c14,d14],
            [a15,b15,c15,d15],
            [a16,b16,c16,d16],
            [a17,b17,c17,d17]])

# Last layer An=Cn=0
N2=M21.LUsolve(M22)
N2.col_del(0)
N2.col_del(1)

M31=Matrix([[a00,b00,c00,d00],
            [a01,b01,c01,d01]])

M31.col(0,lambda i,j: i*exp(-xi*h[0]))
M31.col(1,lambda i,j: i*exp(xi*h[0]))
M31.col(2,lambda i,j: i*exp(-xi*h[0]))
M31.col(3,lambda i,j: i*exp(xi*h[0]))

P=Matrix([[       1,       0,       0,       0,-N1[0,0],-N1[0,1],-N1[0,2],-N1[0,3],       0,       0],
          [       0,       1,       0,       0,-N1[1,0],-N1[1,1],-N1[1,2],-N1[1,3],       0,       0],
          [       0,       0,       1,       0,-N1[2,0],-N1[2,1],-N1[2,2],-N1[2,3],       0,       0],
          [       0,       0,       0,       1,-N1[3,0],-N1[3,1],-N1[3,2],-N1[3,3],       0,       0],
          [       0,       0,       0,       0,       1,       0,       0,       0,-N2[0,0],-N2[0,1]],
          [       0,       0,       0,       0,       0,       1,       0,       0,-N2[1,0],-N2[1,1]],
          [       0,       0,       0,       0,       0,       0,       1,       0,-N2[2,0],-N2[2,1]],
          [       0,       0,       0,       0,       0,       0,       0,       1,-N2[3,0],-N2[3,1]],
          [M31[0,0],M31[0,1],M31[0,2],M31[0,3],       0,       0,       0,       0,       0,       0],
          [M31[1,0],M31[1,1],M31[1,2],M31[1,3],       0,       0,       0,       0,       0,       0]])

R=Matrix(10,1,[0,0,0,0,0,0,0,0,1,0])

A0,B0,C0,D0,A1,B1,C1,D1,B2,D2=P.LUsolve(R)

A0=A0/exp(xi*h[0])
B0=B0/exp(-xi*h[0])
C0=C0/exp(xi*h[0])
D0=D0/exp(-xi*h[0])

A1=A1/exp(xi*h[1])
B1=B1/exp(-xi*h[1])
C1=C1/exp(xi*h[1])
D1=D1/exp(-xi*h[1])

# Substitute constants into stress equation (s_z)
# Determine the stress as a function of z beneath the load (r=0)

sz=s_z(0,50,nu[0],A0,B0,C0,D0)*q
print "\nIntegration (z=50)"
t = time()
f=sz
f=lambdify(xi,f,"mpmath")
szz=-p*a*quad(f,[0,inf],verbose=True)
print "szz: ",szz
print "Time: ",time()-t

The solution:

szz:  -0.573809269227193
Time:  376.523229837

Using the SciPy integration functions provides a faster but less precise result:

from scipy.integrate import quad,Inf
..
szz=-p*a*quad(f,0,Inf)[0]

The solution using SciPy's quad:

szz:  -0.573809269152458
Time:  132.026168823

Another option is to calculate the integration constants at every integration step and pass these into the stress function. The code listing below makes use of the SymPy N function which numerically evaluates these constants. This may be the preferred approach when solving for many layers, which would result in very large integration constant expressions. Note that use is made of the chop option, which basically replaces very small numbers with exact zeros. An option that must still be explored is reduction of the integration constants by applying SymPy's simplification functions. This will involve identifying common terms in these expressions that may be factored out or cancelled.

from scipy.integrate import quad,Inf
from numpy import arange
from sympy import Symbol,Matrix,lambdify,N,exp
from mpmath import sqrt,pi,j1
from time import time

"""
three_b.py

Calculate the stress variation with depth beneath
a three layered system subjected to a circular load.

"""

xi=Symbol('xi')     # Hankel parameter
j0=Symbol('j0')     # Bessel function J0

p=0.69              # Load pressure (MPa)
F=20000.            # Load force (Newton)
a=sqrt(F/(p*pi))    # Load radius (mm)

h=[]
h.append(50.)
h.append(150.)
nu=[]
nu.append(0.45)
nu.append(0.35)
nu.append(0.25)
E=[]
E.append(2000.)
E.append(1000.)
E.append(500.)
mu=[]
for i,j in zip(nu,E):
    mu.append((1+i)/j) 

def hsz(z,nu):
    """
    Constants for Hankel transformed normal stress
    
    """ 
    a1=exp(xi*z)
    b1=-exp(-xi*z)
    c1=-(1-2*nu-xi*z)*exp(xi*z)
    d1=-(1-2*nu+xi*z)*exp(-xi*z)
    return a1,b1,c1,d1

def htrz(z,nu):
    """
    Constants for Hankel transformed shear stress
    
    """ 
    a2=exp(xi*z)
    b2=exp(-xi*z)
    c2=(2*nu+xi*z)*exp(xi*z)
    d2=-(2*nu-xi*z)*exp(-xi*z)
    return a2,b2,c2,d2

def hur(z,nu,X):
    """
    Constants for Hankel transformed radial displacement
    
    """ 
    a3=exp(xi*z)*X
    b3=exp(-xi*z)*X
    c3=-(2-4*nu-xi*z)*exp(xi*z)*X
    d3=(2-4*nu+xi*z)*exp(-xi*z)*X
    return a3,b3,c3,d3

def huz(z,nu,X):
    """
    Constants for Hankel transformed vertical displacement
    
    """ 
    a4=exp(xi*z)*X
    b4=-exp(-xi*z)*X
    c4=(1+xi*z)*exp(xi*z)*X
    d4=(1-xi*z)*exp(-xi*z)*X
    return a4,b4,c4,d4

print "Assemble global matrix"
# Apply boundary condition for surface load (z=0)
a00,b00,c00,d00=hsz(0,nu[0])
a01,b01,c01,d01=htrz(0,nu[0])

# Boundary condition at first upper interface
a02,b02,c02,d02=hsz(h[0],nu[0])
a03,b03,c03,d03=htrz(h[0],nu[0])
a04,b04,c04,d04=hur(h[0],nu[0],1)
a05,b05,c05,d05=huz(h[0],nu[0],1)

# Boundary condition at first lower interface
X=mu[1]/mu[0]
a06,b06,c06,d06=hsz(0,nu[1])
a07,b07,c07,d07=htrz(0,nu[1])
a08,b08,c08,d08=hur(0,nu[1],X)
a09,b09,c09,d09=huz(0,nu[1],X)

# Boundary condition at second upper interface
a10,b10,c10,d10=hsz(h[1],nu[1])
a11,b11,c11,d11=htrz(h[1],nu[1])
a12,b12,c12,d12=hur(h[1],nu[1],1)
a13,b13,c13,d13=huz(h[1],nu[1],1)

# Boundary condition at second lower interface
X=mu[2]/mu[1]
a14,b14,c14,d14=hsz(0,nu[2])
a15,b15,c15,d15=htrz(0,nu[2])
a16,b16,c16,d16=hur(0,nu[2],X)
a17,b17,c17,d17=huz(0,nu[2],X)

# Manipulate first layer
M11=Matrix([[a02,b02,c02,d02],
            [a03,b03,c03,d03],
            [a04,b04,c04,d04],
            [a05,b05,c05,d05]])

M11.col(0,lambda i,j: i*exp(-xi*h[0]))
M11.col(1,lambda i,j: i*exp(xi*h[0]))
M11.col(2,lambda i,j: i*exp(-xi*h[0]))
M11.col(3,lambda i,j: i*exp(xi*h[0]))

M12=Matrix([[a06,b06,c06,d06],
            [a07,b07,c07,d07],
            [a08,b08,c08,d08],
            [a09,b09,c09,d09]])

M12.col(0,lambda i,j: i*exp(-xi*h[1]))
M12.col(1,lambda i,j: i*exp(xi*h[1]))
M12.col(2,lambda i,j: i*exp(-xi*h[1]))
M12.col(3,lambda i,j: i*exp(xi*h[1]))

N1=M11.LUsolve(M12)

# Manipulate second layer
M21=Matrix([[a10,b10,c10,d10],
            [a11,b11,c11,d11],
            [a12,b12,c12,d12],
            [a13,b13,c13,d13]])

M21.col(0,lambda i,j: i*exp(-xi*h[1]))
M21.col(1,lambda i,j: i*exp(xi*h[1]))
M21.col(2,lambda i,j: i*exp(-xi*h[1]))
M21.col(3,lambda i,j: i*exp(xi*h[1]))

M22=Matrix([[a14,b14,c14,d14],
            [a15,b15,c15,d15],
            [a16,b16,c16,d16],
            [a17,b17,c17,d17]])

# Last layer An=Cn=0
N2=M21.LUsolve(M22)
N2.col_del(0)
N2.col_del(1)

M31=Matrix([[a00,b00,c00,d00],
            [a01,b01,c01,d01]])

M31.col(0,lambda i,j: i*exp(-xi*h[0]))
M31.col(1,lambda i,j: i*exp(xi*h[0]))
M31.col(2,lambda i,j: i*exp(-xi*h[0]))
M31.col(3,lambda i,j: i*exp(xi*h[0]))

P=Matrix([[       1,       0,       0,       0,-N1[0,0],-N1[0,1],-N1[0,2],-N1[0,3],       0,       0],
          [       0,       1,       0,       0,-N1[1,0],-N1[1,1],-N1[1,2],-N1[1,3],       0,       0],
          [       0,       0,       1,       0,-N1[2,0],-N1[2,1],-N1[2,2],-N1[2,3],       0,       0],
          [       0,       0,       0,       1,-N1[3,0],-N1[3,1],-N1[3,2],-N1[3,3],       0,       0],
          [       0,       0,       0,       0,       1,       0,       0,       0,-N2[0,0],-N2[0,1]],
          [       0,       0,       0,       0,       0,       1,       0,       0,-N2[1,0],-N2[1,1]],
          [       0,       0,       0,       0,       0,       0,       1,       0,-N2[2,0],-N2[2,1]],
          [       0,       0,       0,       0,       0,       0,       0,       1,-N2[3,0],-N2[3,1]],
          [M31[0,0],M31[0,1],M31[0,2],M31[0,3],       0,       0,       0,       0,       0,       0],
          [M31[1,0],M31[1,1],M31[1,2],M31[1,3],       0,       0,       0,       0,       0,       0]])

R=Matrix(10,1,[0,0,0,0,0,0,0,0,1,0])

A0,B0,C0,D0,A1,B1,C1,D1,B2,D2=P.LUsolve(R)

A0=A0/exp(xi*h[0])
B0=B0/exp(-xi*h[0])
C0=C0/exp(xi*h[0])
D0=D0/exp(-xi*h[0])

A1=A1/exp(xi*h[1])
B1=B1/exp(-xi*h[1])
C1=C1/exp(xi*h[1])
D1=D1/exp(-xi*h[1])

# Substitute constants into stress equation (s_z)
# Determine the stress as a function of z beneath the load (r=0)

def coef(x):
    A=N(A0.subs(xi,x),chop=True)
    B=N(B0.subs(xi,x),chop=True)
    C=N(C0.subs(xi,x),chop=True)
    D=N(D0.subs(xi,x),chop=True)
    return A,B,C,D

def szzz(x):
    A,B,C,D=coef(x)
    f=(A-C*(1-2*nu[0]-x*50))*exp(x*50)-(B+D*(1-2*nu[0]+x*50))*exp(-x*50)
    f=f*j1(x*a)
    return N(f,chop=True)

t = time()
print "Integrate ..."
szz=-p*a*quad(szzz,0,Inf)[0]
print "szz: ",szz
print "Time: ",time()-t

The solution:

szz:  -0.573808894095838
Time:  279.138468981

General solution

Still to do.

Burmister, D.M (1945), The General Theory of Stresses and Displacements in Layered Systems, Journal of Applied Physics, Vol. 16, p.89-94↩
Timoshenko, S., and Goodier, J.N. (1951), Theory of Elasticity, Second Edition, McGraw-Hill, p.368↩

Elastic Layer Analysis (ELA) - sympy/sympy GitHub Wiki

Elastic Layer Analysis (ELA)

Introduction

Problem Formulation

Boundary Conditions

One layer solution

Solving the integration constants

Integrating the stress function

Two layer solution

Simplifying the solution

Three layer solution

General solution

⚠️ GitHub.com Fallback ⚠️

Elastic Layer Analysis (ELA) - sympy/sympy GitHub Wiki

Elastic Layer Analysis (ELA)

Introduction

Problem Formulation

Boundary Conditions

One layer solution

Solving the integration constants

Integrating the stress function

Two layer solution

Simplifying the solution

Three layer solution

General solution

⚠️ **GitHub.com Fallback** ⚠️

⚠️ GitHub.com Fallback ⚠️