a mechanical engineering toolbox
source code - https://github.com/nagordon/mechpy
documentation - https://nagordon.github.io/mechpy/
Neal Gordon
2017-02-20
A quick study to show the quasi-isotropic properties of a lamina / laminate for quasi-isotropicity. In plan
Intution says that there should be a discernible "bump" in Ex as the "quasi" laminate is rotated, when in reality, it is identically ZERO. INTERESTING!
material properties are from Engineering Mechanics of Composite Materials, Daniel
In [5]:
from numpy import pi, zeros, ones, linspace, arange, array, sin, cos, sqrt, pi
from numpy.linalg import solve, inv
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot,figure,xlim,ylim,title,legend, \
grid, show, xlabel,ylabel, tight_layout, savefig
import pandas as pd
In [6]:
%matplotlib inline
In [39]:
def T1(th):
'''Stress Transform for Plane Stress
th=ply angle in degrees
voight notation for stress tranform. sigma1 = T1 @ sigmax
recall T1(th)**-1 == T1(-th)'''
n = sin(th*pi/180)
m = cos(th*pi/180)
T1 = array( [[m**2, n**2, 2*m*n],
[n**2, m**2,-2*m*n],
[-m*n, m*n,(m**2-n**2)]])
return T1
def T2(th):
'''Strain Transform for Plane Stress
th=ply angle in degrees
voight notation for strain transform. epsilon1 = T2 @ epsilonx'''
n = sin(th*pi/180)
m = cos(th*pi/180)
T2 = array( [[m**2, n**2, m*n],
[n**2, m**2,-m*n],
[-2*m*n, 2*m*n, (m**2-n**2)]])
return T2
def Qf(E1,E2,nu12,G12):
'''transversly isptropic compliance matrix. pg 58 herakovich
G12 = E1/(2*(1+nu12)) if isotropic'''
nu21 = E2*nu12/E1
Q = array([[E1/(1-nu12*nu21), E2*nu12/(1-nu12*nu21), 0],
[ E2*nu12/(1-nu12*nu21), E2/(1-nu12*nu21), 0],
[0, 0, G12]])
return Q
def plot_properties(layupname, Ex,Ey,nuxy,Gxy,H):
plt.close('all')
Q = Qf(Ex,Ey,nuxy,Gxy)
plyangle = arange(-45, 45.1, 1)
Exbar = zeros(len(plyangle))
Eybar = zeros(len(plyangle))
Gxybar = zeros(len(plyangle))
Qbar = zeros((len(plyangle),3,3))
for i,th in enumerate(plyangle):
Qbar[i] = solve(T1(th), Q) @ T2(th)
#Qbar = [solve(T1(th),Q) @ T2(th) for th in plyangle]
Aij = Qbar*H
# laminate Stiffness
# | Exbar Eybar Gxybar |
# A = | vxybar vyxbar etasxbar |
# | etaxsbar etaysbar etasybar |
# laminate Comnpliance
aij = zeros((len(plyangle),3,3))
for i, _Aij in enumerate(Aij):
aij[i] = inv(_Aij)
# material properties for whole laminate (Daniel, pg183)
Exbar = [1/(H*_aij[0,0]) for _aij in aij]
Eybar = [1/(H*_aij[1,1]) for _aij in aij]
Gxybar = [1/(H*_aij[2,2]) for _aij in aij]
df = pd.DataFrame({'plyangle':plyangle, 'Exbar':Exbar, 'Eybar':Eybar,'Gxybar':Gxybar})
#print(df)
#df.to_csv('Laminate Properties varying angle for E-Glass Epoxy fabric M10E-3783 {}.csv'.format(layupname))
plt.figure(figsize=(12,8))
plot(plyangle, Exbar, label = r"Modulus: $E_x$")
plot(plyangle, Eybar, label = r"Modulus: $E_y$")
plot(plyangle, Gxybar, label = r"Modulus: $G_{xy}$")
title("Laminate Properties varying angle for E-Glass Unidirectional Tape {}".format(layupname))
xlabel("$\Theta$")
ylim([0,6.5e6])
ylabel("modulus, psi")
legend(loc='best')
grid()
#savefig('Laminate Properties varying angle for E-Glass Epoxy fabric M10E-3783 {}.png'.format(layupname))
In [40]:
layupname = '[0]'
Ex= 6000000.00
Ey= 1500000.00
nuxy= 0.28
Gxy= 620000.00
H= 0.1250
plot_properties(layupname, Ex,Ey,nuxy,Gxy,H)
In [41]:
layupname = '[45]'
Ex= 1700027.42
Ey= 1700027.42
nuxy= 0.37
Gxy= 1079136.69
H = 0.1250
plot_properties(layupname, Ex,Ey,nuxy,Gxy,H)
In [42]:
layupname = '[0, 45, -45, 90, 90, -45, 45, 0]'
Ex= 3000925.78
Ey= 3000925.78
nuxy= 0.29
Gxy= 1159143.21
H = 1.0000
plot_properties(layupname, Ex,Ey,nuxy,Gxy,H)
In [43]:
layupname = '[0, 60, -60, -60, 60, 0]'
Ex= 3000925.78
Ey= 3000925.78
nuxy= 0.29
Gxy= 1159143.21
H = 0.7500
plot_properties(layupname, Ex,Ey,nuxy,Gxy,H)
In [ ]:
In [ ]: