Computational procedures¶
Implementing the theory into Python codes, step-by-step:
- 2D compliance matrix and 2D stiffness matrix for layer materials
- 2D transformation matrices
- Transformation of the 2D stiffness matrix of a layer
- Defining laminates using Python dictionaries
- Building the laminate stiffness matrix
- Solving the load-deformation relations
- Adding thermo-elastic relations
In [1]:
import numpy as np
Materials¶
A material for subsequent codes and examples as a dictionary:
In [2]:
m1={'E1':40000, 'E2':10000, 'v12':0.3, 'G12':5000,
'XT':1200, 'XC':800, 'YT':50, 'YC': 120, 'S12': 60}
Compliance and stiffness matrices for layer materials¶
The 2D compliance matrix was derived in Plane stress
$$ \mathbf{S}=\begin{bmatrix} 1/E_1 & -v_{12}/E_1 & 0 \\ -v_{12}/E_1 & 1/E_{2} & 0 \\ 0 & 0 & 1/G_{12} \end{bmatrix} $$The function takes a material (m) as the only argument:
In [3]:
def S2D(m):
return np.array([[ 1/m['E1'], -m['v12']/m['E1'], 0],
[-m['v12']/m['E1'], 1/m['E2'], 0],
[ 0, 0, 1/m['G12']]], float)
The 2D stiffness matrix of the material is simply the inverse of the compliance matrix:
$$ \mathbf{Q} = \mathbf{S}^{-1} $$Alternatively, the components of the matrix can be computed by
$$ Q_{11}=\frac{E_1}{1-v_{12}v_{21}}, \quad Q_{22}=\frac{E_2}{1-v_{12}v_{21}}, \quad Q_{12}=\frac{v_{12}E_2}{1-v_{12}v_{21}}, \quad Q_{66}=G_{12} $$In the following implementation we call the previous function for the compliance matrix S with a material (m) as argument. Then, the inverse of the compliance matrix S is returned:
In [4]:
def Q2D(m):
S=S2D(m)
return np.linalg.inv(S)
Example:
In [5]:
Q1 = Q2D(m1)
print(Q1)
Transformations¶
The 2D transformation matrices are
$$ \mathbf{T}_{\sigma}=\begin{bmatrix} c^2 & s^2 & 2cs \\ s^2 & c^2 & -2cs \\ -cs & cs & c^2-s^2 \end{bmatrix} $$and $$ \mathbf{T}_{\epsilon}=\begin{bmatrix} c^2 & s^2 & cs \\ s^2 & c^2 & -cs \\ -2cs & 2cs & c^2-s^2 \end{bmatrix} $$
Both functions takes the orientation angle rot as the only argument:
In [6]:
def T2Ds(rot):
c,s = np.cos(np.radians(rot)), np.sin(np.radians(rot))
return np.array([[ c*c , s*s , 2*c*s],
[ s*s , c*c , -2*c*s],
[-c*s, c*s , c*c-s*s]], float)
def T2De(rot):
c,s = np.cos(np.radians(rot)), np.sin(np.radians(rot))
return np.array([[ c*c, s*s, c*s ],
[ s*s, c*c, -c*s ],
[-2*c*s, 2*c*s, c*c-s*s ]], float)
Transformation of the 2D stiffness matrix is given by: $$ \mathbf{Q{'}} = \mathbf{T}_{\sigma}^{-1}\mathbf{Q} \mathbf{T}_{\epsilon} $$
The function Q2Drotate() takes the material stiffness matrix Q and the orientation angle rot as arguments:
In [7]:
def Q2Drotate(Q,rot):
return np.dot(np.linalg.inv( T2Ds(rot) ) , np.dot(Q,T2De(rot)) )
Example:
In [8]:
tempQ=Q2D(m1)
tempQt=Q2Drotate(tempQ,30) # orientation angle = 30
print('Stiffness matrix in the material coordinate system:\n')
print(tempQ)
print('\nTransformed stiffness matrix (30 degrees rotation):\n')
print(tempQt)
The previously definded function Q2d(m) can me modified to include optional rotation as:
In [9]:
def Q2D(m,rot=0):
S=S2D(m)
Q=np.linalg.inv(S)
if rot==0:
return Q
else:
return np.dot(np.linalg.inv( T2Ds(rot) ) , np.dot(Q,T2De(rot)) )
In the first example that follows, no argument is given for the rotation, while an optional rotation is provided in the second example:
In [10]:
Q1 = Q2D(m1)
print(Q1)
print()
Q2 = Q2D(m1,30)
print(Q2)
Layup¶
A layup will be defined as a list of dictionaries as described in Laminates, defining a layup:
In [11]:
layup1 = [ {'mat':m1 , 'ori': 0 , 'thi':0.5} ,
{'mat':m1 , 'ori': 90 , 'thi':0.5} ,
{'mat':m1 , 'ori': 45 , 'thi':0.5} ,
{'mat':m1 , 'ori':-45 , 'thi':0.5} ]
A useful function for computing the total thickness of the laminate:
In [12]:
def laminateThickness(layup):
return sum([layer['thi'] for layer in layup])
tot = laminateThickness(layup1)
print(tot)
Building the laminate stiffness matrix¶
$$ \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{D} \end{bmatrix} = \begin{bmatrix} A_{xx} & A_{xy} & A_{xs} & B_{xx} & B_{xy} & B_{xs} \\ A_{xy} & A_{yy} & A_{ys} & B_{xy} & B_{yy} & B_{ys} \\ A_{xs} & A_{ys} & A_{ss} & B_{xs} & B_{ys} & B_{ss} \\ B_{xx} & B_{xy} & B_{xs} & D_{xx} & D_{xy} & D_{xs} \\ B_{xy} & B_{yy} & B_{ys} & D_{xy} & D_{yy} & D_{ys} \\ B_{xs} & B_{ys} & B_{ss} & D_{xs} & D_{ys} & D_{ss} \end{bmatrix} $$where
$$ \mathbf{A}=\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k-h_{k-1}) $$$$ \mathbf{B}=\frac{1}{2}\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k^2-h_{k-1}^2) $$$$ \mathbf{D}=\frac{1}{3}\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k^3-h_{k-1}^3) $$Lets start with the A sub-matrix: We initiate a 3x3 numpy array of zeros and proceed by adding contributions from the individual layers, i.e.
$$ \mathbf{A}=\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k-h_{k-1}) $$In [13]:
def computeA(layup):
A=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
A = A + Qt*(htop-hbot) # adding the contribution
hbot=htop # bottom for the next layer
return A
# Test:
A = computeA(layup1)
print(A)
Similar functions for the B and D matrices:
In [14]:
def computeB(layup):
B=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
B = B + (Qt*(htop**2-hbot**2))/2 # adding the contribution
hbot=htop # bottom for the next layer
return B
def computeD(layup):
D=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
D = D + (Qt*(htop**3-hbot**3))/3 # adding the contribution
hbot=htop # bottom for the next layer
return D
Alternatively, a function returning the whole 6x6 laminate stiffness matrix:
In [15]:
def laminateStiffnessMatrix(layup):
ABD=np.zeros((6,6),float)
hbot = -laminateThickness(layup)/2
for layer in layup:
htop = hbot + layer['thi']
Qt = Q2D( layer['mat'], layer['ori'])
ABD[0:3,0:3] = ABD[0:3,0:3] + Qt*(htop -hbot )
ABD[0:3,3:6] = ABD[0:3,3:6] + (1/2)*Qt*(htop**2-hbot**2)
ABD[3:6,0:3] = ABD[3:6,0:3] + (1/2)*Qt*(htop**2-hbot**2)
ABD[3:6,3:6] = ABD[3:6,3:6] + (1/3)*Qt*(htop**3-hbot**3)
hbot=htop
return ABD
In [16]:
layup2 = [ {'mat':m1 , 'ori': 0 , 'thi':0.1} ,
{'mat':m1 , 'ori': 45 , 'thi':0.1} ,
{'mat':m1 , 'ori':-45 , 'thi':0.1} ,
{'mat':m1 , 'ori': 0 , 'thi':0.1} ]
myABD = laminateStiffnessMatrix(layup2)
print(np.round(myABD))
Solving load-deformation relations¶
$$ \begin{bmatrix} \mathbf{N} \\ \mathbf{M} \end{bmatrix} = \begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{D} \end{bmatrix} \begin{bmatrix} \boldsymbol{\varepsilon}^0 \\ \boldsymbol{\kappa} \end{bmatrix} $$Or in details; $$ \begin{bmatrix} N_x \\ N_y \\ N_{xy} \\ M_x \\ M_y \\ M_{xy} \end{bmatrix} = \begin{bmatrix} A_{xx} & A_{xy} & A_{xs} & B_{xx} & B_{xy} & B_{xs} \\ A_{xy} & A_{yy} & A_{ys} & B_{xy} & B_{yy} & B_{ys} \\ A_{xs} & A_{ys} & A_{ss} & B_{xs} & B_{ys} & B_{ss} \\ B_{xx} & B_{xy} & B_{xs} & D_{xx} & D_{xy} & D_{xs} \\ B_{xy} & B_{yy} & B_{ys} & D_{xy} & D_{yy} & D_{ys} \\ B_{xs} & B_{ys} & B_{ss} & D_{xs} & D_{ys} & D_{ss} \end{bmatrix} \begin{bmatrix} \varepsilon_x^0 \\ \varepsilon_y^0 \\ \gamma_{xy}^0 \\ \kappa_x \\ \kappa_y \\ \kappa_{xy} \end{bmatrix} $$
Material and layup for the proceeding examples:
In [17]:
m1={'E1':140000, 'E2':10000, 'v12':0.3, 'G12':5000}
layup = [ {'mat':m1 , 'ori': 0 , 'thi':1} ,
{'mat':m1 , 'ori': 90 , 'thi':1} ,
{'mat':m1 , 'ori': 45 , 'thi':1} ]
Case 1: deformations are known:¶
In this case the procedure is pretty straigth forward:
In [18]:
ABD=laminateStiffnessMatrix(layup)
deformations1=(0.001,0,0,0,0,0)
loads1=np.dot(ABD,deformations1)
print(np.round(loads1))
Case 2: loads are known:¶
Compared to the previous case, this is just a case where we solve for the deformations:
In [19]:
ABD=laminateStiffnessMatrix(layup)
loads2=(1000,0,0,0,0,0)
deformations2=np.linalg.solve(ABD,loads2)
print(np.round(deformations2,5))
#alternatively:
deformations3=np.dot( np.linalg.inv(ABD) ,loads2 )
print(np.round(deformations3,5))
Case 3: some loads and some deformations are known:¶
Now it becomes slightly more tricky. Somehow we must rearrange the set of equations such that the known quantities are split from the unknown quantities.
The function that follows makes use of keyword arguments (**kwargs) in the following way:
- the initial assumtion is that all forces and moments are known and equal to zero
- any argument (such as
Nx=100, orKx=0.1) overrides the initial assumption. - if an argument is prescribing a deformation (such as
ey = 0), the force or moment becomes the unknown (Nyin this case)
In [20]:
def solveLaminateLoadCase(ABD,**kwargs):
# ABD: the laminate stiffness matrix
# kwargs: prescribed loads or deformations
loads=np.zeros((6))
defor=np.zeros((6))
loadKeys=('Nx','Ny','Nxy','Mx','My','Mxy') # valid load keys
defoKeys=('ex','ey','exy','Kx','Ky','Kxy') # valid deformation keys
knownKeys=['Nx','Ny','Nxy','Mx','My','Mxy'] # initially assume known loads
knownVals=np.array([0,0,0,0,0,0],float)
for key in kwargs:
if key in loadKeys:
i=loadKeys.index(key)
knownKeys[i]=key
knownVals[i]=kwargs[key]
if key in defoKeys:
i=defoKeys.index(key)
knownKeys[i]=key
knownVals[i]=kwargs[key]
M1=-np.copy(ABD)
M2= np.copy(ABD)
ID= np.identity(6)
for k in range(0,6):
if knownKeys[k] in loadKeys:
M2[:,k]=-ID[:,k]
else:
M1[:,k]=ID[:,k]
sol=np.dot( np.linalg.inv(M1), np.dot(M2,knownVals))
for i in range(0,6):
if knownKeys[i] == loadKeys[i]:
loads[i]=knownVals[i]
defor[i]=sol[i]
if knownKeys[i] == defoKeys[i]:
defor[i]=knownVals[i]
loads[i]=sol[i]
return loads,defor
Example:
In [21]:
loads, defs = solveLaminateLoadCase(ABD=ABD, Nx=100., ey=0.001, exy=0.005, Mx=50, My=25,Kxy=0.002 )
loadKeys=('Nx','Ny','Nxy','Mx','My','Mxy')
defoKeys=('ex','ey','exy','Kx','Ky','Kxy')
for i in range(0,6):
print(loadKeys[i],'=',loads[i])
for i in range(0,6):
print(defoKeys[i],'=',defs[i])
Adding thermal loading¶
A function for computing the thermal loads:
$$ \mathbf{N}_{th}= \Big[\Delta T \sum_{k=1}^{n} \mathbf{Q}'_k \boldsymbol{\alpha}'_k (h_k-h_{k-1}) \Big] $$$$ \mathbf{M}_{th}=\Big[\Delta T \frac{1}{2} \sum_{k=1}^{n} \mathbf{Q}'_k \boldsymbol{\alpha}'_k (h_k^2-h_{k-1}^2) \Big] $$where
$$ \boldsymbol{\alpha}'= \mathbf{T}_{\varepsilon}^{-1} \boldsymbol{\alpha} $$In [22]:
def thermalLoad(layup,dT):
NTMT=np.zeros(6,float) #thermal loads
hbot = -laminateThickness(layup)/2
for layer in layup:
htop = hbot + layer['thi']
Qt = Q2D( layer['mat'], layer['ori'])
a123=[layer['mat']['a1'], layer['mat']['a2'], 0]
aXYZ=np.dot( T2De(-layer['ori']) , a123 )
NTMT[0:3] = NTMT[0:3] + dT*( np.dot(Qt,aXYZ) )*(htop -hbot )
NTMT[3:6] = NTMT[3:6] + 0.5*dT*( np.dot(Qt,aXYZ) )*(htop**2-hbot**2)
hbot=htop
return NTMT
Example:
In [23]:
cfrp = {'E1':140000, 'E2':8000, 'v12':0.28, 'G12':4000, 'a1':-0.5E-6, 'a2': 25E-6}
layup1 =[ {'mat':cfrp , 'ori': 0 , 'thi':0.5} ,
{'mat':cfrp , 'ori': 90 , 'thi':0.5} ]
thermLoad=thermalLoad(layup1,-100)
ABD=laminateStiffnessMatrix(layup1)
defs= np.dot(np.linalg.inv(ABD) , thermLoad)
print(defs)
In [24]:
from plotlib import illustrateCurvatures
%matplotlib inline
illustrateCurvatures(Kx=defs[3], Ky=defs[4], Kxy=defs[5])
Summary: functions¶
In [25]:
def S2D(m):
return np.array([[ 1/m['E1'], -m['v12']/m['E1'], 0],
[-m['v12']/m['E1'], 1/m['E2'], 0],
[ 0, 0, 1/m['G12']]], float)
def T2Ds(rot):
c,s = np.cos(np.radians(rot)), np.sin(np.radians(rot))
return np.array([[ c*c , s*s , 2*c*s],
[ s*s , c*c , -2*c*s],
[-c*s, c*s , c*c-s*s]], float)
def T2De(rot):
c,s = np.cos(np.radians(rot)), np.sin(np.radians(rot))
return np.array([[ c*c, s*s, c*s ],
[ s*s, c*c, -c*s ],
[-2*c*s, 2*c*s, c*c-s*s ]], float)
def Q2D(m,rot=0):
S=S2D(m)
Q=np.linalg.inv(S)
if rot==0:
return Q
else:
return np.dot(np.linalg.inv( T2Ds(rot) ) , np.dot(Q,T2De(rot)) )
def Q2Drotate(Q,rot):
return np.dot(np.linalg.inv( T2Ds(rot) ) , np.dot(Q,T2De(rot)) )
def computeA(layup):
A=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
A = A + Qt*(htop-hbot) # adding the contribution
hbot=htop # bottom for the next layer
return A
def computeB(layup):
B=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
B = B + (Qt*(htop**2-hbot**2))/2 # adding the contribution
hbot=htop # bottom for the next layer
return B
def computeD(layup):
D=np.zeros((3,3),float)
hbot = -laminateThickness(layup)/2 # bottom of first layer
for layer in layup: # loop, all layers
htop = hbot + layer['thi'] # top of current layer
Qt = Q2D( layer['mat'], layer['ori']) # transformed stiffness matrix
D = D + (Qt*(htop**3-hbot**3))/3 # adding the contribution
hbot=htop # bottom for the next layer
return D
def laminateStiffnessMatrix(layup):
ABD=np.zeros((6,6),float)
hbot = -laminateThickness(layup)/2
for layer in layup:
htop = hbot + layer['thi']
Qt = Q2D( layer['mat'], layer['ori'])
ABD[0:3,0:3] = ABD[0:3,0:3] + Qt*(htop -hbot )
ABD[0:3,3:6] = ABD[0:3,3:6] + (1/2)*Qt*(htop**2-hbot**2)
ABD[3:6,0:3] = ABD[3:6,0:3] + (1/2)*Qt*(htop**2-hbot**2)
ABD[3:6,3:6] = ABD[3:6,3:6] + (1/3)*Qt*(htop**3-hbot**3)
hbot=htop
return ABD
def solveLaminateLoadCase(ABD,**kwargs):
# ABD: the laminate stiffness matrix
# kwargs: prescribed loads or deformations
loads=np.zeros((6))
defor=np.zeros((6))
loadKeys=('Nx','Ny','Nxy','Mx','My','Mxy') # valid load keys
defoKeys=('ex','ey','exy','Kx','Ky','Kxy') # valid deformation keys
knownKeys=['Nx','Ny','Nxy','Mx','My','Mxy'] # initially assume known loads
knownVals=np.array([0,0,0,0,0,0],float)
for key in kwargs:
if key in loadKeys:
i=loadKeys.index(key)
knownKeys[i]=key
knownVals[i]=kwargs[key]
if key in defoKeys:
i=defoKeys.index(key)
knownKeys[i]=key
knownVals[i]=kwargs[key]
M1=-np.copy(ABD)
M2= np.copy(ABD)
ID= np.identity(6)
for k in range(0,6):
if knownKeys[k] in loadKeys:
M2[:,k]=-ID[:,k]
else:
M1[:,k]=ID[:,k]
sol=np.dot( np.linalg.inv(M1), np.dot(M2,knownVals))
for i in range(0,6):
if knownKeys[i] == loadKeys[i]:
loads[i]=knownVals[i]
defor[i]=sol[i]
if knownKeys[i] == defoKeys[i]:
defor[i]=knownVals[i]
loads[i]=sol[i]
return loads,defor
def thermalLoad(layup,dT):
NTMT=np.zeros(6,float) #thermal loads
hbot = -laminateThickness(layup)/2
for layer in layup:
htop = hbot + layer['thi']
Qt = Q2D( layer['mat'], layer['ori'])
a123=[layer['mat']['a1'], layer['mat']['a2'], 0]
aXYZ=np.dot( T2De(-layer['ori']) , a123 )
NTMT[0:3] = NTMT[0:3] + dT*( np.dot(Qt,aXYZ) )*(htop -hbot )
NTMT[3:6] = NTMT[3:6] + 0.5*dT*( np.dot(Qt,aXYZ) )*(htop**2-hbot**2)
hbot=htop
return NTMT
