python计算复合材料层合板ABD刚度矩阵、预测层合板强度

鄙人不才，在学校的时候没有学python，复合材料力学也是一知半解，后来工作的时候遇到了需要计算复合材料层合板ABD刚度矩阵的内容，然后恰好在学习python，于是花时间编写了下预测这方面的内容，然后后期还编写了通过Tsai-Wu，Tsai-Hill失效准则预测层合板强度的模块，不过代码编的比较烂，后期会学习python项目管理，然后封装好。


备注：由于本人代码水平有限，现有代码还需要更改，主要是模块间的数据读入和数据传递比较混乱，但模块的计算能力没有问题，各位观众老爷在学习的时候可以参考《复合材料力学》相关书籍的经典层合板理论对照阅读。


@Time : 2018-12-28
@Author : Allen Pen
@E-mail : [e-mail]

# 文件读取 这里需要进行数据读入操作，可以通过GUI填入，或者是读取.txt文档。或者是通过细观力学计算得来。

理论部分

待更新。

复合材料基础理论

待更新。

单层板的宏观力学分析

层合板的宏观力学分析

单层板的细观力学分析

待更新。

计算逻辑与实现逻辑

计算逻辑

# 计算单层板弹性常数

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018-12-28
# @Author  : Allen Pen
# @E-mail  : [e-mail]

import numpy as np
import numpy.linalg as lg
from CLPTC.calculator.input_file_reader import ply_properties,mat_properties

class Lamina():
    '''
    注意：
    正轴    -- 与铺层纤维方向相同的方向为正轴向，以1,2,3方向表示
    偏轴    -- 一般为层合板自然轴定义的坐标系，以x,y,z表示，相对于铺层的坐标系
    铺层角  -- 偏轴x转至正轴1的夹角，逆时针转向为正
    泊松比  -- 泊松比使用国外的书和软件的定义方式，即ν_12=-ε_2/ε_1 ，ν_21=-ε_1/ε_2
              使用nu12,nu12（OptiStruct） = nuxy（ANSYS） = ν_12(国外教材)=ν_1（复合材料力学）=ν_21（复合材料力学）=-ε_2/ε_1

    matid      material id
    e1         Young Modulus in direction 1
    e2         Young Modulus in direction 2
    g12        in-plane shear modulus

    nu21       Poisson's ratio 21,nu21=-ε2/ε1
    nu12       Poisson's ratio 12: use formula nu21/e1 = nu12/e2


    st1,st2    allowable tensile stresses for directions 1 and 2
    sc1,sc2    allowable compressive stresses for directions 1 and 2
    ss12       allowable in-plane stress for shear
    strn       allowable strain for direction 1


    plyid      id of the composite lamina
    t           ply thickness
    theta       ply angle

    S           （单层的）正轴柔量矩阵
    Q           （单层的）正轴模量矩阵

    T_stress    应力转换矩阵
    T_strain    应变转换矩阵

    S_offaxis   （单层的）偏轴柔量矩阵
    Q_offaxis   （单层的）偏轴模量矩阵

    Sij	       compliance component 	柔量分量
    Qij	       modulus component    	模量分量

    '''
    def __init__(self):
        self.matid    = None
        self.plyid    = None
        self.t        = None
        self.theta    = None


        self.e1   = None
        self.e2   = None
        self.g12  = None
        self.nu21 = None
        self.nu12 = None
        self.st1  = None
        self.st2  = None
        self.sc1  = None
        self.sc2  = None
        self.ss12 = None
        self.strn = None



        self.S        = None
        self.Q        = None
        self.T_stress = None
        self.T_strain = None
        self.S_offaxis= None
        self.Q_offaxis= None

        self.lamina_Ex  = None
        self.lamina_Ey  = None
        self.lamina_Gxy = None
        self.lamina_nuxy= None



        self.laminates= []
        self.cos      = None
        self.cos2t    = None
        self.sin      = None
        self.sin2t    = None

    def calc_SQ(self):

        '''
        计算单层的正轴刚度
        :return: 只返回 S、Q，中间的s11等变量无法访问
        '''
        e1 = self.e1
        e2 = self.e2
        nu21 = self.nu21
        nu12 = self.nu12
        g12 = self.g12

        s11 = 1 / e1
        s22 = 1 / e2
        s66 = 1 / g12
        s12 = -nu12 / e2
        s21 = -nu21 / e1
        s16 = s61 = s26 = s62 = 0

        qm  = (1 - nu21*nu12)
        q11 = e1/qm
        q22 = e2/qm
        q66 = g12
        q12 = nu12 * e1/qm
        q21 = nu21 * e2/qm
        q16 = q61 = q26 = q62 = 0

        self.S = np.array(
            [[s11,s12,s16],
             [s21,s22,s26],
             [s61,s62,s66]],dtype = float)

        self.Q = np.array(
            [[q11,q12,q16],
             [q21,q22,q26],
             [q61,q62,q66]],dtype = float)

        # 计算各单层的偏轴刚度
        theta = self.theta

        self.cos   = np.cos( np.deg2rad(   theta ) )
        self.sin   = np.sin( np.deg2rad(   theta ) )
        # self.cos2t = np.cos( np.deg2rad( 2*self.theta ) )
        # self.sin2t = np.sin( np.deg2rad( 2*self.theta ) )
        cos    = self.cos
        sin    = self.sin
        cos2   = cos**2
        sin2   = sin**2
        sincos = sin*cos
        # cos2t  = self.cos2t
        # sin2t  = self.sin2t

        # 应力转换矩阵，用于将偏轴应力转换至正轴应力
        self.T_stress = np.array(
            [[   cos2,  sin2,  2*sincos],
             [   sin2,  cos2, -2*sincos],
             [-sincos,sincos, cos2-sin2]],dtype=float)

        # 应变转换矩阵，用于将偏轴应变转换至正轴应变
        self.T_strain = np.array(
            [[     cos2,    sin2,    sincos],
             [     sin2,    cos2,   -sincos],
             [-2*sincos,2*sincos, cos2-sin2]],dtype=float)

        # 计算偏轴刚度矩阵
        self.Q_offaxis = np.dot(np.dot(lg.inv(self.T_stress), self.Q), self.T_strain)
        self.S_offaxis = np.dot(np.dot(lg.inv(self.T_strain), self.S), self.T_stress)

        # 铺层等效的工程弹性常数
        self.lamina_Ex  = 1/self.S_offaxis[0,0]
        self.lamina_Ey  = 1/self.S_offaxis[1,1]
        self.lamina_Gxy = 1/self.S_offaxis[2,2]
        self.lamina_nuxy= -self.lamina_Ex*self.S_offaxis[0,1]

def Get_lamina_prop(plyid):
    lamina_prop = Lamina()

    #铺层属性
    lamina_prop.matid = ply_properties['材料索引值'][plyid]
    lamina_prop.t     = ply_properties['厚度'][plyid]
    lamina_prop.theta = ply_properties['铺层角'][plyid]

    # 铺层的材料属性
    lamina_prop.e1   = mat_properties['E1'][lamina_prop.matid]
    lamina_prop.e2   = mat_properties['E2'][lamina_prop.matid]
    lamina_prop.nu21 = mat_properties['nu21'][lamina_prop.matid]
    lamina_prop.g12  = mat_properties['G12'][lamina_prop.matid]
    lamina_prop.nu12 = lamina_prop.nu21/lamina_prop.e1*lamina_prop.e2

    lamina_prop.calc_SQ()


    return lamina_prop

层合板ABD刚度矩阵计算

```python
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018-12-28
# @Author  : Allen Pen

import numpy as np
import numpy.linalg as lg
import matplotlib.pyplot as plt
from CLPTC.calculator.lamina import Get_lamina_prop #从
from CLPTC.calculator.input_file_reader import total_ply_num,z_coord_dict,z_total,load_info,ply_properties,mat_properties

class Laminate():

    """

    equiv_e1         equivalent laminate modulus in 1 direction
    equiv_e2         equivalent laminate modulus in 2 direction
    equiv_g12        equivalent laminate shear modulus in 12 direction
    equiv_nu21       equivalent laminate Poisson ratio in 21 direction
    equiv_nu12       equivalent laminate Poisson ratio in 12 direction

    plyid      id of the composite lamina
    plies      list of plies
    plyts      ply thicknesses for this laminate
    t          total thickness of the laminate

    A          面内刚度系数矩阵
    B          耦合刚度系数矩阵
    D          弯曲刚度系数矩阵
    ABD        ABD刚度矩阵

    用于计算柔度矩阵的中间变量
    A_temp   = A^-1
    B_temp   = -A^-1*B
    H_temp   = B*A^-1
    D_temp   = -B*A^-1*B+D

    A_inverted          面内柔度系数矩阵
    B_inverted          耦合柔度系数矩阵
    H_inverted          一般情况下 H_inverted = B_inverted
    D_inverted          弯曲柔度系数矩阵
    ABHD_inverted      ABD柔度矩阵

    epsilon_0 板中面应变
    K       板中面翘曲率
    """

    def __init__(self):

        self.plyid= None
        self.plies= None
        self.plyts= None
        self.t    = None

        self.equiv_e1   = None
        self.equiv_e2   = None
        self.equiv_g12  = None
        self.equiv_nu21 = None
        self.equiv_nu12 = None

        self.A = None
        self.B = None
        self.D = None
        self.ABD = None

        self.A_temp = None
        self.B_temp = None
        self.H_temp = None
        self.D_temp = None

        self.A_inverted = None
        self.B_inverted = None
        self.H_inverted = None
        self.D_inverted = None
        self.ABHD_inverted = None

        self.epsilon_K_offaxis = None
        self.epsilon_offaxis_0 = None #板中面应变
        self.K = None    #板中面翘曲率
        self.load_info = None #用于在最大强度比小于1的情况下更新载荷

        self.N_M = None #层合板当前能承受的载荷


        self.count = 0 #用于层合板失效层数的统计

        _ = self.init_build() #在类实例化时，将ABD 刚度阵、柔度阵，以及层的应变计算出来

        # 设置np数组输出的精度，注意，数组依然按float64储存
        np.set_printoptions(formatter={'float': '{: 10.3f}'.format})

    def calc_ABD(self):
        """
        用于计算A,B,D刚度矩阵，以及A',B',D'柔度矩阵
        ref. 沈观林.《复合材料力学》
        :return:
                ABD             A,B,D刚度矩阵
                ABHD_inverted   A',B',D'柔度矩阵
        """
        self.A = np.zeros([3,3], dtype=float)
        self.B = np.zeros([3,3], dtype=float)
        self.D = np.zeros([3,3], dtype=float)

        self.ply_prop ={}
        for plyid in range(1,total_ply_num + 1):
            ply = Get_lamina_prop(plyid)

            #将计算后的偏轴刚度、应变转换矩阵储存起来
            self.ply_prop['Q_'+str(plyid)] = ply.Q
            self.ply_prop['T_strain_' +str(plyid)] = ply.T_strain
            self.ply_prop['T_stress_'+str(plyid)] = ply.T_stress

            zk = z_coord_dict['z' + str(plyid)]
            zk_1 = z_coord_dict['z' + str(plyid-1)]


            for i in range(3):
                for j in range(3):
                    self.A[i][j] +=     ply.Q_offaxis[i][j]*(zk    - zk_1  )
                    self.B[i][j] +=1/2.*ply.Q_offaxis[i][j]*(zk**2 - zk_1**2)
                    self.D[i][j] +=1/3.*ply.Q_offaxis[i][j]*(zk**3 - zk_1**3)

            part1 = np.concatenate([self.A, self.B], axis=1)
            part2 = np.concatenate([self.B, self.D], axis=1)

            self.ABD = np.concatenate([part1, part2], axis=0)


    def calc_ABHD_inverted(self):

        # 计算ABD'柔度矩阵
        self.ABHD_inverted = lg.inv(self.ABD)

    def calc_equivalent_modulus(self):
        """
        用于计算层合板等效的工程弹性常数
        :return:
                equiv_e1     等效E1
                equiv_e2     等效E2
                equiv_nu21   等效nu21
                equiv_nu12   等效nu12
                equiv_g12    等效G12
        """
        AI = lg.inv(self.ABD)
        a11, a12, a22, a33 = AI[0,0], AI[0,1], AI[1,1], AI[2,2]
        self.equiv_e1    = 1./(z_total*a11)
        self.equiv_e2    = 1./(z_total*a22)
        self.equiv_nu21  = - a12 / a11
        self.equiv_nu12  = - a12 / a22
        self.equiv_g12   = 1./(z_total*a33)

层合板强度预测

待更新。

最后的最后

欢迎大家点赞、评论及转载，转载请注明出处！


如果觉得我帮助到了你：
  为我打call，不如为我打款！