---

Mechpy Tutorials

a mechanical engineering toolbox

source code - https://github.com/nagordon/mechpy
documentation - https://nagordon.github.io/mechpy/web/

---

Neal Gordon
2017-02-20

---

Composite Plate Mechanics with Python

reference: hyer page 584. 617

The motivation behind this talk is to explore the capability of python as a scientific computation tool as well as solve a typical calcuation that could either be done by hand, or coded. I find coding to be a convient way to learn challenging mathmatics because I can easily replicate my work in the future when I can't remember the details of the calcuation or, if there are any errors, they can typically be easily fixed and the other calcuations re-ran without a lot of duplcation of effort.

Composite mechanics can be very iterative by nature and is easiest to employ linear algebra to find displacements, strains and stresses of composites. Coding solutions is also handy when visualizations are required.

For this example, we are interested in calcuating the stress critical ply in a simple asymteric composite plate with a pressure load applied. We can chooose a variety of boundary conditions of our plate, but this solution is limited to 2 dimensional displacements, x and z. If we are interested in 3 dimensional displacements, the problem becomes much more challenging as partial differentiation of the governing equations gives us a PDE, which is more challenging to solve.

The steps to solving are

• Identify governing and equilibrium equations
• import python required libraries
• declare symbolic variables
• declare numeric variables, including material properties, plate dimensions, and plate pressure
• solve 4th order differntial equation with 7 constants
• apply plate boundary conditions and acquire u(x) and w(x) displacement functions
• acquire strain equations from displacement
• acquire stress equations from strain
• determine critical ply from highest ply stress ratio

# Import Python modules and 
import numpy as np
from sympy import *
from pprint import pprint

# printing and plotting settings 
init_printing(use_latex='mathjax')
get_ipython().magic('matplotlib inline') # inline plotting

x,y,q = symbols('x,y,q')

As mentioned before, if we want to perform a 3 dimensional displacement model of the composite plate, we would have 6 reaction forces that are a function of x and y. Those 6 reaction forces are related by 3 equalibrium equations

# # hyer page 584
# # Equations of equilibrium
# Nxf = Function('N_x')(x,y)
# Nyf = Function('N_y')(x,y)
# Nxyf = Function('N_xy')(x,y)
# Mxf = Function('M_x')(x,y)
# Myf = Function('M_y')(x,y)
# Mxyf = Function('M_xy')(x,y)

# symbols for force and moments
Nx,Ny,Nxy,Mx,My,Mxy = symbols('N_x,N_y,N_xy,M_x,M_y,M_xy')
Nxf,Nyf,Nxyf,Mxf,Myf,Mxyf = symbols('Nxf,Nyf,Nxyf,Mxf,Myf,Mxyf')

Eq(0,diff(Nx(x,y), x)+diff(Nxy(x,y),y))

$$0 = \frac{\partial}{\partial x} \operatorname{N_{x}}{\left (x,y \right )} + \frac{\partial}{\partial y} \operatorname{N_{xy}}{\left (x,y \right )}$$

Eq(0,diff(Nxy(x,y), x)+diff(Ny(x,y),y))

$$0 = \frac{\partial}{\partial x} \operatorname{N_{xy}}{\left (x,y \right )} + \frac{\partial}{\partial y} \operatorname{N_{y}}{\left (x,y \right )}$$

Eq(0,  diff(Mx(x,y),x,2) + 2*diff(Mxy(x,y),x,y) + diff(My(x,y) ,y,2)+ q  )

$$0 = q + \frac{\partial^{2}}{\partial x^{2}} \operatorname{M_{x}}{\left (x,y \right )} + 2 \frac{\partial^{2}}{\partial x\partial y} \operatorname{M_{xy}}{\left (x,y \right )} + \frac{\partial^{2}}{\partial y^{2}} \operatorname{M_{y}}{\left (x,y \right )}$$

What makes composite plates special is the fact that they typically not isotropic. This is handled by the 6x6 ABD matrix that defines the composites properties axially, in bending, and the coupling between the two.

# composite properties 
A11,A22,A66,A12,A16,A26,A66 = symbols('A11,A22,A66,A12,A16,A26,A66')
B11,B22,B66,B12,B16,B26,B66 = symbols('B11,B22,B66,B12,B16,B26,B66')
D11,D22,D66,D12,D16,D26,D66 = symbols('D11,D22,D66,D12,D16,D26,D66')

## constants of integration when solving differential equation
C1,C2,C3,C4,C5,C6 = symbols('C1,C2,C3,C4,C5,C6')

# plate and composite parameters
th,a,b = symbols('th,a,b')

# displacement functions
u0 = Function('u0')(x,y)
v0 = Function('v0')(x,y)
w0 = Function('w0')(x,y)

Let's compute our 6 displacement conditions which is where our PDE's show up

Nxf = A11*diff(u0,x) + A12*diff(v0,y) + A16*(diff(u0,y) + diff(v0,x)) - B11*diff(w0,x,2) - B12*diff(w0,y,2) - 2*B16*diff(w0,x,y)
Eq(Nx, Nxf)

$$N_{x} = A_{11} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + A_{12} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + A_{16} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - B_{11} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - B_{12} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{16} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Nyf = A12*diff(u0,x) + A22*diff(v0,y) + A26*(diff(u0,y) + diff(v0,x)) - B12*diff(w0,x,2) - B22*diff(w0,y,2) - 2*B26*diff(w0,x,y)
Eq(Ny,Nyf)

$$N_{y} = A_{12} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + A_{22} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + A_{26} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - B_{12} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - B_{22} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{26} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Nxyf = A16*diff(u0,x) + A26*diff(v0,y) + A66*(diff(u0,y) + diff(v0,x)) - B16*diff(w0,x,2) - B26*diff(w0,y,2) - 2*B66*diff(w0,x,y) 
Eq(Nxy,Nxyf)

$$N_{xy} = A_{16} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + A_{26} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + A_{66} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - B_{16} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - B_{26} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{66} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Mxf = B11*diff(u0,x) + B12*diff(v0,y) + B16*(diff(u0,y) + diff(v0,x)) - D11*diff(w0,x,2) - D12*diff(w0,y,2) - 2*D16*diff(w0,x,y)
Eq(Mx,Mxf)

$$M_{x} = B_{11} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + B_{12} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + B_{16} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - D_{11} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - D_{12} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 D_{16} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Myf = B12*diff(u0,x) + B22*diff(v0,y) + B26*(diff(u0,y) + diff(v0,x)) - D12*diff(w0,x,2) - D22*diff(w0,y,2) - 2*D26*diff(w0,x,y)
Eq(My,Myf)

$$M_{y} = B_{12} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + B_{22} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + B_{26} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - D_{12} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - D_{22} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 D_{26} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Mxyf = B16*diff(u0,x) + B26*diff(v0,y) + B66*(diff(u0,y) + diff(v0,x)) - D16*diff(w0,x,2) - D26*diff(w0,y,2) - 2*D66*diff(w0,x,y)
Eq(Mxy,Mxyf)

$$M_{xy} = B_{16} \frac{\partial}{\partial x} \operatorname{u_{0}}{\left (x,y \right )} + B_{26} \frac{\partial}{\partial y} \operatorname{v_{0}}{\left (x,y \right )} + B_{66} \left(\frac{\partial}{\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{0}}{\left (x,y \right )}\right) - D_{16} \frac{\partial^{2}}{\partial x^{2}} \operatorname{w_{0}}{\left (x,y \right )} - D_{26} \frac{\partial^{2}}{\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 D_{66} \frac{\partial^{2}}{\partial x\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

Now, combine our 6 displacement conditions with our 3 equalibrium equations to get three goverening equations

eq1 = diff(Nxf,x) + diff(Nxf,y)
eq1

$$A_{11} \frac{\partial^{2}}{\partial x^{2}} \operatorname{u_{0}}{\left (x,y \right )} + A_{11} \frac{\partial^{2}}{\partial x\partial y} \operatorname{u_{0}}{\left (x,y \right )} + A_{12} \frac{\partial^{2}}{\partial x\partial y} \operatorname{v_{0}}{\left (x,y \right )} + A_{12} \frac{\partial^{2}}{\partial y^{2}} \operatorname{v_{0}}{\left (x,y \right )} + A_{16} \left(\frac{\partial^{2}}{\partial x\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{2}}{\partial x^{2}} \operatorname{v_{0}}{\left (x,y \right )}\right) + A_{16} \left(\frac{\partial^{2}}{\partial y^{2}} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{2}}{\partial x\partial y} \operatorname{v_{0}}{\left (x,y \right )}\right) - B_{11} \frac{\partial^{3}}{\partial x^{3}} \operatorname{w_{0}}{\left (x,y \right )} - B_{11} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{w_{0}}{\left (x,y \right )} - B_{12} \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - B_{12} \frac{\partial^{3}}{\partial y^{3}} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{16} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{16} \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )}$$

eq2 = diff(Nxyf,x) + diff(Nyf,y)
eq2

$$A_{12} \frac{\partial^{2}}{\partial x\partial y} \operatorname{u_{0}}{\left (x,y \right )} + A_{16} \frac{\partial^{2}}{\partial x^{2}} \operatorname{u_{0}}{\left (x,y \right )} + A_{22} \frac{\partial^{2}}{\partial y^{2}} \operatorname{v_{0}}{\left (x,y \right )} + A_{26} \left(\frac{\partial^{2}}{\partial y^{2}} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{2}}{\partial x\partial y} \operatorname{v_{0}}{\left (x,y \right )}\right) + A_{26} \frac{\partial^{2}}{\partial x\partial y} \operatorname{v_{0}}{\left (x,y \right )} + A_{66} \left(\frac{\partial^{2}}{\partial x\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{2}}{\partial x^{2}} \operatorname{v_{0}}{\left (x,y \right )}\right) - B_{12} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{w_{0}}{\left (x,y \right )} - B_{16} \frac{\partial^{3}}{\partial x^{3}} \operatorname{w_{0}}{\left (x,y \right )} - B_{22} \frac{\partial^{3}}{\partial y^{3}} \operatorname{w_{0}}{\left (x,y \right )} - 3 B_{26} \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 2 B_{66} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{w_{0}}{\left (x,y \right )}$$

eq3 = diff(Mxf,x,2) + 2*diff(Mxyf,x,y) + diff(Myf,y,2) + q
eq3

$$B_{11} \frac{\partial^{3}}{\partial x^{3}} \operatorname{u_{0}}{\left (x,y \right )} + B_{12} \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{u_{0}}{\left (x,y \right )} + B_{12} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{v_{0}}{\left (x,y \right )} + B_{16} \left(\frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{3}}{\partial x^{3}} \operatorname{v_{0}}{\left (x,y \right )}\right) + 2 B_{16} \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{u_{0}}{\left (x,y \right )} + B_{22} \frac{\partial^{3}}{\partial y^{3}} \operatorname{v_{0}}{\left (x,y \right )} + B_{26} \left(\frac{\partial^{3}}{\partial y^{3}} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{v_{0}}{\left (x,y \right )}\right) + 2 B_{26} \frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{v_{0}}{\left (x,y \right )} + 2 B_{66} \left(\frac{\partial^{3}}{\partial x\partial y^{2}} \operatorname{u_{0}}{\left (x,y \right )} + \frac{\partial^{3}}{\partial x^{2}\partial y} \operatorname{v_{0}}{\left (x,y \right )}\right) - D_{11} \frac{\partial^{4}}{\partial x^{4}} \operatorname{w_{0}}{\left (x,y \right )} - 2 D_{12} \frac{\partial^{4}}{\partial x^{2}\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} - 4 D_{16} \frac{\partial^{4}}{\partial x^{3}\partial y} \operatorname{w_{0}}{\left (x,y \right )} - D_{22} \frac{\partial^{4}}{\partial y^{4}} \operatorname{w_{0}}{\left (x,y \right )} - 4 D_{26} \frac{\partial^{4}}{\partial x\partial y^{3}} \operatorname{w_{0}}{\left (x,y \right )} - 4 D_{66} \frac{\partial^{4}}{\partial x^{2}\partial y^{2}} \operatorname{w_{0}}{\left (x,y \right )} + q$$

Yikes, I do not want to solve that (at least right now). If we make the assumption that the plate has equal displacement of y in the x and y direction, then we can simply things ALOT! These simplifications are valid for cross ply unsymmetric laminates plate, Hyer pg 616. This is applied by setting some of our material properties to zero. $ A16=A26=D16=D26=B16=B26=B12=B66=0 $ Almost like magic, we now have some equations that aren't so scary.

u0 = Function('u0')(x)
v0 = Function('v0')(x)
w0 = Function('w0')(x)

Nxf = A11*diff(u0,x) + A12*diff(v0,y) - B11*diff(w0,x,2)
Eq(Nx, Nxf)

$$N_{x} = A_{11} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )} - B_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )}$$

Nyf = A12*diff(u0,x) + A22*diff(v0,y) - B22*diff(w0,y,2)
Eq(Ny,Nyf)

$$N_{y} = A_{12} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )}$$

Nxyf = A66*(diff(u0,y) + diff(v0,x))
Eq(Nxy,Nxyf)

$$N_{xy} = A_{66} \frac{d}{d x} \operatorname{v_{0}}{\left (x \right )}$$

Mxf = B11*diff(u0,x) - D11*diff(w0,x,2) - D12*diff(w0,y,2)
Eq(Mx,Mxf)

$$M_{x} = B_{11} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )} - D_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )}$$

Myf = B22*diff(v0,y) - D12*diff(w0,x,2) - D22*diff(w0,y,2)
Eq(My,Myf)

$$M_{y} = - D_{12} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )}$$

Mxyf = 0
Eq(Mxy,Mxyf)

$$M_{xy} = 0$$

Now we are getting somewhere. Finally we can solve the differential equations

dsolve(diff(Nx(x)))

$$\operatorname{N_{x}}{\left (x \right )} = C_{1}$$

dsolve(diff(Mx(x),x,2)+q)

$$\operatorname{M_{x}}{\left (x \right )} = C_{1} + C_{2} x - \frac{q x^{2}}{2}$$

Now solve for u0 and w0 with some pixie dust

eq4 = (Nxf-C1)
eq4

$$A_{11} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )} - B_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )} - C_{1}$$

eq5 = Mxf -( -q*x**2 + C2*x + C3 )
eq5

$$B_{11} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )} - C_{2} x - C_{3} - D_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )} + q x^{2}$$

eq6 = Eq(solve(eq4,diff(u0,x))[0] , solve(eq5, diff(u0,x))[0])
eq6

$$\frac{1}{A_{11}} \left(B_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )} + C_{1}\right) = \frac{1}{B_{11}} \left(C_{2} x + C_{3} + D_{11} \frac{d^{2}}{d x^{2}} \operatorname{w_{0}}{\left (x \right )} - q x^{2}\right)$$

w0f = dsolve(eq6, w0)
w0f

$$\operatorname{w_{0}}{\left (x \right )} = - \frac{A_{11} C_{2} x^{3}}{6 A_{11} D_{11} - 6 B_{11}^{2}} + \frac{A_{11} q x^{4}}{12 A_{11} D_{11} - 12 B_{11}^{2}} + C_{1} + C_{5} x + \frac{x^{2} \left(- A_{11} C_{3} + B_{11} C_{4}\right)}{2 A_{11} D_{11} - 2 B_{11}^{2}}$$

eq7 = Eq(solve(eq6, diff(w0,x,2))[0] , solve(eq4,diff(w0,x,2))[0])
eq7

$$\frac{1}{A_{11} D_{11} - B_{11}^{2}} \left(- A_{11} C_{2} x - A_{11} C_{3} + A_{11} q x^{2} + B_{11} C_{1}\right) = \frac{1}{B_{11}} \left(A_{11} \frac{d}{d x} \operatorname{u_{0}}{\left (x \right )} - C_{1}\right)$$

u0f = dsolve(eq7)
u0f

$$\operatorname{u_{0}}{\left (x \right )} = \frac{1}{A_{11} D_{11} - B_{11}^{2}} \left(- \frac{B_{11} C_{1}}{2} x^{2} - B_{11} C_{3} x - B_{11} C_{4} D_{11} + \frac{B_{11} q}{3} x^{3} + C_{2} D_{11} x + \frac{B_{11}^{3} C_{4}}{A_{11}}\right)$$

---