Composite A and D matrix calculation from ply stacking sequence.py

# calculate equivalent G matrix terms

import numpy as np
import matplotlib.pyplot as plt

# define the plate input parameters
class input_params():
    def __init__(self):
        self.fs_length = (362-2*20)*0.001 # m
        self.fs_width = 0.080 # m
        self.fs_divs = 40
        # material properties AS4/3501-6 from NASA report
        self.fs_material_E1 = 138.0E9 # Pa
        self.fs_material_E2 = 9.0E9 # Pa
        self.fs_material_G12 = 6.9E9# Pa
        self.fs_material_nu = 0.30
        # self.fs_material_S1_T = 303*10**3 * 6894.76 # Pa
        # self.fs_material_S1_C = 245*10**3 * 6894.76 # Pa
        # self.fs_material_S2_T = 5.47*10**3 * 6894.76 # Pa
        # self.fs_material_S2_C = self.fs_material_S2_T*5 # Pa - estimate
        # self.fs_material_S12 = 16.5*10**3 * 6894.76 # Pa
#         self.fs_material_eps1_T = 13900 # micro strain
#         self.fs_material_eps1_C = 11700 # micro strain
        self.fs_material_t = 0.0001175 # m
        self.fs_ply_ori = [45,0] # half-layup from outer ply in, orientation of each ply wrt span direction (X), Z-up

def calc_D(INPUT):
    # Analytical: Calculate the plate bending stiffness properties (D matrix)
    # calculate the material invariants
    nu_21 = INPUT.fs_material_nu*INPUT.fs_material_E2/INPUT.fs_material_E1
    Q11 = INPUT.fs_material_E1 / (1 - INPUT.fs_material_nu * nu_21)
    Q22 = INPUT.fs_material_E2 / (1 - INPUT.fs_material_nu * nu_21)
    Q12 = INPUT.fs_material_nu * INPUT.fs_material_E2 / (1 - INPUT.fs_material_nu * nu_21)
    Q66 = INPUT.fs_material_G12
    U1 = 3 / 8 * (Q11 + Q22) + 1 / 4 * Q12 + 1 / 2 * Q66
    U2 = 1 / 2 * (Q11 - Q22)
    U3 = 1 / 8 * (Q11 + Q22) - 1 / 4 * Q12 - 1 / 2 * Q66
    U4 = 1 / 8 * (Q11 + Q22) + 3 / 4 * Q12 - 1 / 2 * Q66
    U5 = 1 / 8 * (Q11 + Q22) - 1 / 4 * Q12 + 1 / 2 * Q66
    #print(Q11, Q22, Q12, Q66)

    # sum over all plies
    nplies = range(len(INPUT.fs_ply_ori) * 2)
    h = np.arange(len(INPUT.fs_ply_ori) + 1) * INPUT.fs_material_t
    h = np.hstack([-h[-1:-len(INPUT.fs_ply_ori) - 1:-1], h])
    ply_ori = np.hstack([INPUT.fs_ply_ori,
                         INPUT.fs_ply_ori[-1:-len(INPUT.fs_ply_ori) - 1:-1]])
    D = np.zeros([3, 3])
    A = np.zeros([3, 3])
    for ply in nplies:
        rot = ply_ori[ply] * np.pi / 180  # in rad
        Q11_bar = U1 + U2 * np.cos(2 * rot) + U3 * np.cos(4 * rot)
        Q22_bar = U1 - U2 * np.cos(2 * rot) + U3 * np.cos(4 * rot)
        Q12_bar = U4 - U3 * np.cos(4 * rot)
        Q66_bar = U5 - U3 * np.cos(4 * rot)
        Q16_bar = 1 / 2 * U2 * np.sin(2 * rot) + U3 * np.sin(4 * rot)
        Q26_bar = 1 / 2 * U2 * np.sin(2 * rot) - U3 * np.sin(4 * rot)
        D += 1 / 3 * np.array([[Q11_bar, Q12_bar, Q16_bar],
                               [Q12_bar, Q22_bar, Q26_bar],
                               [Q16_bar, Q26_bar, Q66_bar]]) * ((h[ply + 1] ** 3) - (h[ply] ** 3))
        A += np.array([[Q11_bar, Q12_bar, Q16_bar],
                       [Q12_bar, Q22_bar, Q26_bar],
                       [Q16_bar, Q26_bar, Q66_bar]]) * (h[ply + 1] - h[ply])
    return A, D, h, ply_ori, Q11, Q12, Q22, Q66


# plate properties
INPUT = input_params()
print(vars(INPUT))

# Analytical: Calculate the plate bending stiffness properties (D matrix)
A, D, h, ply_ori, Q11, Q12, Q22, Q66 = calc_D(INPUT)

pass

note that this calculates the stiffness matrices from the laminate, but assumes the same material for all plies.