[发明专利]预测低速冲击下复合材料层合板渐进失效的有限元方法

摘要

本发明涉及复合材料损伤预测，旨在提供预测低速冲击下复合材料层合板渐进失效的有限元方法。该预测低速冲击下复合材料层合板渐进失效的有限元方法包括过程建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型；建立复合材料弹塑性损伤本构模型；基于ABAQUS‑VUMAT用户动态材料子程序模块，运用后向欧拉算法实现提出的弹塑性损伤本构模型；对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度。本发明利用ABAQUS‑VUMAT用户子程序来数值实现所建立的将塑性和损伤联合的弹塑性损伤本构模型，该模型同时考虑塑性和材料性能退化的影响，能准确预测含塑性特征的复合材料在低速冲击下的渐进损伤失效。

主权项

预测低速冲击下复合材料层合板渐进失效的有限元方法，其特征在于，包括下述过程：一、建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型；二、建立复合材料弹塑性损伤本构模型；三、基于ABAQUS‑VUMAT用户动态材料子程序模块，运用后向欧拉算法实现提出的弹塑性损伤本构模型，求解应力、应变和损伤；四、对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；所述过程一中，建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型：基于ABAQUS建立冲锤、复合材料层合板、支撑板部件，分别设置材料属性和划分网格，再进行ASSEMBLY对其组装之后设置分析步和通用接触属性；其中，复合材料层合板每层根据铺层角度进行铺层，不同铺层之间设置ABAQUS自带的双线性内聚力单元；所述过程二具体包括下述步骤：步骤(1)：建立含损伤复合材料层合板层内本构关系；复合材料应力‑应变本构方程：S＝Cd：Ee；其中，为未损伤材料的有效应力；S为损伤材料的名义应力；Ep为格林‑拉格朗日应变张量的塑性部分，Ee＝E‑Ep为格林‑拉格朗日应变张量的弹性部分；Cd(C，d1，d2)为含损伤材料的四阶弹性张量，C为未损伤材料的四阶弹性张量，d1和d2是分别对应纤维断裂和基体开裂的损伤变量；步骤(2)：建立塑性模型，具体建立方式为：(a)考虑面外应力的复合材料塑性流动准则为：F(S‾,E~p)=3f(S‾)-k(E~p)=0;]]>f(S‾ij)=12[(S‾22-S‾33)2+4S‾232+2a66(S‾132+S‾122)];]]>其中，F为屈服方程；a66为描述各向异性材料与塑性相关的常数，由偏轴拉伸测试确定；f为塑性势函数；k和分别为塑性硬化应力和等效塑性应变；所述是指未损伤材料的有效应力，其中i，j用来确定应力方向；所述指层合板面内垂直于纤维方向的有效应力；所述指面外垂直于纤维的有效应力；所述指面外剪切有效应力；所述指面内剪切有效应力；(b)塑性硬化应力k为：k=β(E~p)n;]]>其中，β和n是为满足实验硬化曲线的常数；塑性变形假定发生在损伤材料的未损伤区域，塑性流动准则和硬化准则表达在有效应力空间中；(c)等效塑性应变率和塑形应变率分别为：E~·p=λ·;]]>E·p=λ·∂F/∂S‾;]]>其中，为一致塑性因子；步骤(3)：建立基于应变描述的PUCK失效初始判据和损伤演化准则，具体建立方式为：(d)对于纤维拉伸和压缩，损伤初始判据为：F11T*=(E11E0,1T*)2-1≥0,F11C=(E11E0,1C)2-1≥0;]]>其中，和分别为纤维拉伸和压缩的初始失效应变；所述T*，C分别指拉伸和压缩；所述E11是指纤维方向应变；所述分别指纤维拉伸和压缩失效判断因子；纤维拉伸和压缩的损伤演化准则为：d11T*(C)=Ef,1T*(C)Ef,1T*(C)-E0,1T*(C)(1-E0,1T*(C)E11);]]>其中，所述是指纤维拉伸和压缩损伤变量；所述E11是指纤维方向应变；所述是指纤维损伤变量达到1的纤维临界拉伸和压缩失效应变；所述是指纤维损伤变量为零的纤维初始拉伸和压缩失效应变；(e)对于基体拉伸损伤失效初始判据为：F22T*=(E22E0,2T*)2-1≥0;]]>其中，所述指基体拉伸失效判断因子；所述E22是指基体方向应变；所述是指基体损伤变量为零的基体初始拉伸失效应变；基体拉伸损伤演化准则为：d22T*=Ef,2T*Ef,2T*-E0,2T*(1-E0,2T*E22);]]>其中，为基体损伤变量达到1时基体临界拉伸失效应变；所述是指基体拉伸损伤变量；(f)对于基体压缩损伤初始判据为：其中，N是关于失效断裂面的法向方向，T和L是关于失效断裂面的切向方向；Yc是横向压缩强度，断裂平面上的应力Sij(i，j＝L，T，N)由笛卡氏坐标系下的Piola‑Kirchhoff应力Sij(i，j＝1，2，3)通过旋转矩阵T(α)旋转获得，T(α)为笛卡尔坐标系到断裂面坐标系的旋转矩阵；所述SNN是指断裂面的法向应力；SNT，SNL是指断裂面的切向应力，μNL，μTN为断裂面面内两个切向方向摩擦系数，θf为断裂面的断裂角；所述指基体压缩失效判断因子；所述是指笛卡尔坐标系下面内剪切强度；所述是指在断裂平面内的横向剪切强度；所述S123是指在笛卡尔坐标系下的六个Piola‑Kirchhoff应力Sij(i，j＝1，2，3)；所述SLTN是指在断裂面坐标系下的六个Piola‑Kirchhoff应力Sij(i，j＝L，T，N)；所述T(α)T是指T(α)的转置矩阵；所述90°是指采用角度制计量的90度；基体压缩损伤演化准则为：d22C(γNT,γNL)=γγmaxγγmax-γγf(1-γfγγ),γγ=γNT2+γNL2.;]]>其中，所述γγ是指断裂面联合剪切应变；和是联合剪切应变的初始和最大应变，γNT和γNL是断裂面的剪切应变；是指基体压缩损伤变量；所述过程三具体包括下述步骤：步骤(4)：通过用户子程序VUMAT的SDV定义第n+1增量步的开始时的初始状态变量值，同时也是第n增量步结束时的状态变量值在第n+1增量步开始时，VUMAT读入；其中，所述n是指第n增量步，所述En是指第n增量步结束时的格林‑拉格朗日总应变张量，所述是指第n增量步结束时的格林‑拉格朗日塑性应变张量，所述是指第n增量步结束时的等效塑性应变，所述是指第n增量步结束时的未损伤材料的有效应力，所述kn是指第n增量步结束时的塑性硬化应力，所述dij，n是指第n增量步结束时的损伤变量；步骤(5)：VUMAT由应变增量驱动，计算试应力，将试应力代入到步骤(2)屈服方程中；该公式是步骤(2)中的公式在第n+1增量步的特定计算；如果Fn+1≤0，则处于弹性阶段，将所有试应力和应变更新为n+1增量步状态变量如果Fn+1＞0，塑性加载出现，根据后向欧拉隐式算法，实现试应力到屈服面的最近点返回，将试应力赋予步骤(4)中的迭代初始条件；变量是塑性一致因子的增量Δλn+1的函数，运用牛顿‑拉夫森算法求解再更新应变和有效应力，直到Fn+1≤0，结束迭代，得到n+1增量步状态变量步骤(6)：应变和有效应力得到更新后，根据有效应力和应变代入到步骤(3)的PUCK失效准则判断是否出现损伤，如有损伤再通过损伤演化公式获得损伤变量，再根据步骤(1)通过有效应力和损伤变量计算名义应力Sn+1；所述过程四具体为：将过程一建立的模型主文件和过程三建立的ABAQUS‑VUMAT用户子程序联合，使用ABAQUS/EXPLICT方法对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；即完成低速冲击载荷下弹塑性复合材料层合板渐进失效特性的预测。