TMM4175 Polymer Composites

Home About Python Links Next Previous Table of Contents

Computational procedures for laminates

Implementing the theory into Python code, step-by-step:

  • 2D compliance matrix and 2D stiffness matrix
  • 2D transformation matrices
  • Transformation of the 2D stiffness matrix
  • 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

The material for subsequent codes and examples:

In [2]:
m1={'E1':40000, 'E2':10000, 'v12':0.3, 'G12':5000, 
    'XT':1200, 'XC':800, 'YT':50, 'YC': 120, 'S12': 60}

The 2D compliance matrix was defined 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)

Example:

In [4]:
temp = S2D(m1)
print(temp)

[[ 2.5e-05 -7.5e-06  0.0e+00]
 [-7.5e-06  1.0e-04  0.0e+00]
 [ 0.0e+00  0.0e+00  2.0e-04]]

The 2D stiffness matrix of the material is simply the inverse of the compliance matrix:

$$ \mathbf{Q} = \mathbf{S}^{-1} $$

where the components are given 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 [5]:
def Q2D(m):
    S=S2D(m)
    return np.linalg.inv(S)

Example:

In [6]:
temp = Q2D(m1)
print(temp)

[[40920.71611253  3069.05370844     0.        ]
 [ 3069.05370844 10230.17902813     0.        ]
 [    0.             0.          5000.        ]]

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 a as the only argument:

In [7]:
def T2Ds(a):
    c,s = np.cos(np.radians(a)), np.sin(np.radians(a))
    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(a):
    c,s = np.cos(np.radians(a)), np.sin(np.radians(a))
    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 Q2Dtransform which takes the material stiffness matrix (Q) and the orientation angle (a) as arguments:

In [8]:
def Q2Dtransform(Q,a):
    return np.dot(np.linalg.inv( T2Ds(a) ) , np.dot(Q,T2De(a)) )

Example:

In [9]:
tempQ=Q2D(m1)
tempQt=Q2Dtransform(tempQ,30) # orientation angle = 30
print('Stiffness matrix in material coordinate system:\n')
print(tempQ)
print('\nTransformed stiffness matrix (30 degrees):\n')
print(tempQt)

Stiffness matrix in material coordinate system:

[[40920.71611253  3069.05370844     0.        ]
 [ 3069.05370844 10230.17902813     0.        ]
 [    0.             0.          5000.        ]]

Transformed stiffness matrix (30 degrees):

[[28558.18414322  7758.95140665  9352.40989125]
 [ 7758.95140665 13212.91560102  3936.98249419]
 [ 9352.40989125  3936.98249419  9689.89769821]]

Defining a layup:

Layups will be defined as dictionaries as described in Definitions and notation:

In [10]:
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}  ]

Example: the orientation of layer number 3 (index 2) is found by:

In [11]:
layup1[2]['ori']
Out[11]:
45

Example: the shear modulus $G_{12}$ for the material used in the first layer:

In [12]:
layup1[0]['mat']['G12']
Out[12]:
5000

Example: number of layers is

In [13]:
n = len(layup1)
print(n)

4

The total thickness::

In [14]:
def laminateThickness(layup):
    return sum([layer['thi'] for layer in layup])

laminateThickness(layup1)
Out[14]:
2.0

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:

In [15]:
temp_A = np.zeros((3,3), float)
print(temp_A)

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

Then add contribution from the individual layers, i.e. $$ \mathbf{A}=\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k-h_{k-1}) $$

In [16]:
# The very bottom is also to bottom of the first layer:
hbot = -laminateThickness(layup1) / 2

# A loop over the layers:
for layer in layup1:
    # material:
    m = layer['mat']
    # stiffness matrix:
    Q = Q2D(m)
    # orientation angle:
    a = layer['ori']
    #transformed stiffness materix:
    Qt = Q2Dtransform(Q,a)
    # compute the top coordinate of current layer:
    htop = hbot + layer['thi']
    # for the next layer, the hbot is the current htop:
    hbot=htop
    # finally, add current layer contribution to the A matrix (temp_A):
    temp_A = temp_A + temp
    
print('A=')
print(temp_A)
    

A=
[[163682.86445013  12276.21483376      0.        ]
 [ 12276.21483376  40920.71611253      0.        ]
 [     0.              0.          20000.        ]]

Now, lets put all the stuff into a function for the A-matrix:

In [17]:
def computeA(layup):
    A=np.zeros((3,3),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q,a)
        htop = hbot + layer['thi']   # top of current layer
        A = A + Qt*(htop-hbot)
        hbot=htop                    # for the next layer
    return A

Example:

In [18]:
computeA(layup1)
Out[18]:
array([[4.48976982e+04, 1.23913043e+04, 0.00000000e+00],
       [1.23913043e+04, 4.48976982e+04, 9.09494702e-13],
       [0.00000000e+00, 9.09494702e-13, 1.62531969e+04]])

Continue with the B-matrix: $$ \mathbf{B}=\frac{1}{2}\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k^2-h_{k-1}^2) $$

In [19]:
def computeB(layup):
    B=np.zeros((3,3),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q,a)
        htop = hbot + layer['thi']   # top of current layer
        B = B + (1/2)*Qt*(htop**2-hbot**2)
        hbot=htop                    # for the next layer
    return B

And then, the D-matrix: $$ \mathbf{D}=\frac{1}{3}\sum_{k=1}^n\mathbf{Q{'}}_{k}(h_k^3-h_{k-1}^3) $$

In [20]:
def computeD(layup):
    D=np.zeros((3,3),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q,a)
        htop = hbot + layer['thi']   # top of current layer
        D = D + (1/3)*Qt*(htop**3-hbot**3)
        hbot=htop                    # for the next layer
    return D

Example: for the laminate with layup

In [21]:
layup2 = [ {'mat':m1 , 'ori':  0  , 'thi':1} , 
           {'mat':m1 , 'ori': 90  , 'thi':1} ,
           {'mat':m1 , 'ori':  0  , 'thi':1}  ]

we find the sub-matrices as

In [22]:
computeA(layup2)
Out[22]:
array([[9.20716113e+04, 9.20716113e+03, 1.73830940e-13],
       [9.20716113e+03, 6.13810742e+04, 1.70542246e-12],
       [1.73830940e-13, 1.70542246e-12, 1.50000000e+04]])
In [23]:
computeB(layup2)
Out[23]:
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
In [24]:
computeD(layup2)
Out[24]:
array([[8.95140665e+04, 6.90537084e+03, 1.44859116e-14],
       [6.90537084e+03, 2.55754476e+04, 1.42118538e-13],
       [1.44859116e-14, 1.42118538e-13, 1.12500000e+04]])

For efficiency we may improve the code by merging the computeA, computeB and computeD into one function that return A, B and D as individual variables simultanously:

In [25]:
def computeABD(layup):
    A=np.zeros((3,3),float)
    B=np.zeros((3,3),float)
    D=np.zeros((3,3),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q,a)
        htop = hbot + layer['thi']   # top of current layer
        A = A +       Qt*(htop   -hbot   )
        B = B + (1/2)*Qt*(htop**2-hbot**2)
        D = D + (1/3)*Qt*(htop**3-hbot**3)
        hbot=htop                    # for the next layer
    return A,B,D

Example:

In [26]:
myA, myB, myD = computeABD(layup2)

print(myA)
print('---------------------------------------------')
print(myB)
print('---------------------------------------------')
print(myD)

[[9.20716113e+04 9.20716113e+03 1.73830940e-13]
 [9.20716113e+03 6.13810742e+04 1.70542246e-12]
 [1.73830940e-13 1.70542246e-12 1.50000000e+04]]
---------------------------------------------
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
---------------------------------------------
[[8.95140665e+04 6.90537084e+03 1.44859116e-14]
 [6.90537084e+03 2.55754476e+04 1.42118538e-13]
 [1.44859116e-14 1.42118538e-13 1.12500000e+04]]

The A, B and D can be concatenated into one 6 x 6 matrix as shown in this example:

In [27]:
myABBD=np.concatenate(( np.concatenate((myA,myB),axis=1) , np.concatenate((myB,myD),axis=1)),axis=0)
print(np.round(myABBD))

[[92072.  9207.     0.     0.     0.     0.]
 [ 9207. 61381.     0.     0.     0.     0.]
 [    0.     0. 15000.     0.     0.     0.]
 [    0.     0.     0. 89514.  6905.     0.]
 [    0.     0.     0.  6905. 25575.     0.]
 [    0.     0.     0.     0.     0. 11250.]]

And finally, another function that returns the whole 6 x 6 laminate stiffness matrix more efficiently than the previous approches. The local intermediate variables are also eliminated:

In [28]:
def laminateStiffnessMatrix(layup):
    ABD=np.zeros((6,6),float)
    hbot = -laminateThickness(layup)/2        # bottom of first layer
    for layer in layup:
        Qt = Q2Dtransform(Q2D(layer['mat']), layer['ori'])
        htop = hbot + layer['thi']   # top of current layer
        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                    # for the next layer
    return ABD

Example:

In [29]:
myABBD = laminateStiffnessMatrix(layup2)
print(np.round(myABBD))

[[92072.  9207.     0.     0.     0.     0.]
 [ 9207. 61381.     0.     0.     0.     0.]
 [    0.     0. 15000.     0.     0.     0.]
 [    0.     0.     0. 89514.  6905.     0.]
 [    0.     0.     0.  6905. 25575.     0.]
 [    0.     0.     0.     0.     0. 11250.]]

Solving laminate load-deformation

$$ \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 [30]:
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 [31]:
ABD=laminateStiffnessMatrix(layup)
deformations1=(0.001,0,0,0,0,0)
loads1=np.dot(ABD,deformations1)
print(np.round(loads1))

[195.  40.  33. -97.  31.  33.]

Case 2: loads are known:

Compared to the previous case, this is just a case where we solve for the deformations:

In [32]:
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))

[ 0.01397 -0.00108 -0.00766  0.01048 -0.00326 -0.00972]
[ 0.01397 -0.00108 -0.00766  0.01048 -0.00326 -0.00972]

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, or Kx=0.1) overrides the initial assumption.
  • if an argument is prescribing a deformation (such as ey = 0), the force or moment becomes the unknown (Ny in this case)
In [33]:
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 [34]:
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])

Nx = 100.0
Ny = 267.96584767608005
Nxy = 179.77587666830732
Mx = 50.0
My = 25.0
Mxy = 123.05527887701813
ex = -0.0010226711477403423
ey = 0.001
exy = 0.005
Kx = -0.0010821318614371504
Ky = -0.002375894915094696
Kxy = 0.002

Adding thermal loading

A function to compute 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 [35]:
def thermalLoad(layup,dT):
    NTMT=np.zeros(6,float) #thermal loads
    hbot = -laminateThickness(layup)/2
    for layer in layup:
        m = layer['mat']
        Q = Q2D(m)
        a = layer['ori']
        Qt = Q2Dtransform(Q,a)
        a123=[m['a1'], m['a2'], 0]
        aXYZ=np.dot( T2De(-layer['ori']) , a123 )
        htop = hbot + layer['thi'] 
        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 [36]:
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)

[-5.40034911e-04 -5.40034911e-04  2.37179524e-19 -1.93504616e-03
  1.93504616e-03  2.36974809e-19]
In [37]:
from plotlib import illustrateCurvatures
%matplotlib inline
illustrateCurvatures(Kx=defs[3], Ky=defs[4], Kxy=defs[5])

Next Previous

Disclaimer:This site is about polymer composites, designed for educational purposes. Consumption and use of any sort & kind is solely at your own risk.
Fair use: I spent some time making all the pages, and even the figures and illustrations are my own creations. Obviously, you may steal whatever you find useful here, but please show decency and give some acknowledgment if or when copying. Thanks! Contact me: nils.p.vedvik@ntnu.no www.ntnu.edu/employees/nils.p.vedvik

Copyright 2021, All right reserved, I guess.