- - - -
# 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




```python
# 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


```python
# # 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')
```


```python
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 )}$$




```python
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 )}$$




```python
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.


```python
# 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


```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
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


```python
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 )}$$




```python
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 )}$$




```python
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.


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


```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
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 )}$$




```python
Mxyf = 0
Eq(Mxy,Mxyf)
```




$$M_{xy} = 0$$



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


```python
dsolve(diff(Nx(x)))
```




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




```python
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


```python
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}$$




```python
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}$$




```python
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)$$




```python
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}}$$




```python
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)$$




```python
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)$$



- - - - 

22d478c344305fe81892507b412ae64daddc2a41
ipynb
Jupyter Notebook
tutorials/Composite_Plate_Mechanics_with_Python_Theory.ipynb
nagordon/mechpy
aae2315b883f6af7cd90a8451d170744bbf1053a
MIT
2017-01-27T04:40:30.000Z
2021-12-03T03:46:07.000Z
tutorials/Composite_Plate_Mechanics_with_Python_Theory.ipynb
Lunreth/mechpy
aae2315b883f6af7cd90a8451d170744bbf1053a
MIT
2016-03-01T00:42:38.000Z
2020-03-04T15:45:39.000Z
tutorials/Composite_Plate_Mechanics_with_Python_Theory.ipynb
Lunreth/mechpy
aae2315b883f6af7cd90a8451d170744bbf1053a
MIT
2016-04-25T14:12:34.000Z
2021-07-07T17:46:35.000Z
Qwen/Qwen-72B
1. YES
2. YES