CN106503292B - 预测低速冲击下复合材料层合板渐进失效的有限元方法 - Google Patents
预测低速冲击下复合材料层合板渐进失效的有限元方法 Download PDFInfo
- Publication number
- CN106503292B CN106503292B CN201610833225.4A CN201610833225A CN106503292B CN 106503292 B CN106503292 B CN 106503292B CN 201610833225 A CN201610833225 A CN 201610833225A CN 106503292 B CN106503292 B CN 106503292B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msubsup
- mover
- damage
- refer
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
本发明涉及复合材料损伤预测，旨在提供预测低速冲击下复合材料层合板渐进失效的有限元方法。该预测低速冲击下复合材料层合板渐进失效的有限元方法包括过程：建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型；建立复合材料弹塑性损伤本构模型；基于ABAQUS‑VUMAT用户动态材料子程序模块，运用后向欧拉算法实现提出的弹塑性损伤本构模型；对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度。本发明利用ABAQUS‑VUMAT用户子程序来数值实现所建立的将塑性和损伤联合的弹塑性损伤本构模型，该模型同时考虑塑性和材料性能退化的影响，能准确预测含塑性特征的复合材料在低速冲击下的渐进损伤失效。
Description
技术领域
本发明是关于复合材料损伤预测领域，特别涉及预测低速冲击下复合材料层合板渐进失效的有限元方法。
背景技术
当前，复合材料正广泛应用于航空航天、风力发电、压力容器、汽车等高新技术领域。但是低速冲击损伤对复合材料的强度、刚度及使用寿命都有较大的影响，因此必须清楚的了解在低速冲击下复合材料的渐进损伤演化过程。
一些复合材料如T300/914，AS4/PEEK等，在横向和剪切加载下具有显著的非线性和塑性，复合材料的不可逆变形是由各种不同的失效机制导致的，比如纤维基体的损伤，纤维/基体的界面分离，塑性变形的累积等，此时弹性损伤本构模型预测这类复合材料的低速冲击响应并不准确，因此提出一种可以准确预测有塑性行为的复合材料层合板的低速冲击响应的弹塑性损伤本构模型尤为必要。近年来针对复合材料层合板提出的弹塑性损伤本构模型未能合理考虑面外应力的影响和合适的失效准则及损伤演化，且基本集中在复合材料的静态损伤研究，而将弹塑性损伤本构模型用来准确预测复合材料层合板在低速冲击下的损伤研究很少。
ABAQUS软件可以对基于二维HASHIN失效准则对弹性复合材料进行低速冲击下的渐进失效研究，但无法直接运用较有优势的PUCK失效准则，更无法直接对有塑性特征的复合材料层合板进行低速冲击研究。
发明内容
本发明的主要目的在于克服现有技术中的不足，提供一种能准确预测含塑性特征的复合材料在低速冲击下的渐进损伤失效的方法。为解决上述技术问题，本发明的解决方案是：
提供预测低速冲击下复合材料层合板渐进失效的有限元方法，包括下述过程：
一、建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型；
二、建立复合材料弹塑性损伤本构模型；
三、基于(使用FORTRAN语言编写的)ABAQUS-VUMAT用户动态材料子程序模块，运用后向欧拉算法实现提出的弹塑性损伤本构模型，求解应力、应变和损伤；
四、对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；
所述过程一中，建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型：基于ABAQUS建立冲锤、复合材料层合板、支撑板部件，分别设置材料属性和划分网格，再进行ASSEMBLY对其组装之后设置分析步和通用接触属性；
其中，复合材料层合板每层根据铺层角度进行铺层，不同铺层之间设置ABAQUS自带的双线性内聚力单元；
所述过程二具体包括下述步骤：
步骤(1)：建立含损伤复合材料层合板层内本构关系；
复合材料应力-应变本构方程：
其中，为未损伤材料的有效应力；S为损伤材料的名义应力；Ep为格林-拉格朗日应变张量的塑性部分，Ee＝E-Ep为格林-拉格朗日应变张量的弹性部分；Cd(C，d1，d2)为含损伤材料的四阶弹性张量，C为未损伤材料的四阶弹性张量，d1和d2是分别对应纤维断裂和基体开裂的损伤变量；
步骤(2)：建立塑性模型，具体建立方式为：
(a)考虑面外应力的复合材料塑性流动准则为：
其中，F为屈服方程；a66为描述各向异性材料与塑性相关的常数，由偏轴拉伸测试确定(其数值约为1.25左右)；f为塑性势函数；k和分别为塑性硬化应力和等效塑性应变；所述是指未损伤材料的有效应力，其中i，j用来确定应力方向；所述指层合板面内垂直于纤维方向的有效应力；所述指面外垂直于纤维的有效应力；所述指面外剪切有效应力；所述指面内剪切有效应力；
(b)塑性硬化应力k为：
其中，β和n是为满足实验硬化曲线的常数(取值分别为600MPa和0.272左右)；塑性变形假定发生在损伤材料的未损伤区域，塑性流动准则和硬化准则表达在有效应力空间中；
(c)等效塑性应变率和塑形应变率分别为：
其中，为一致塑性因子；
步骤(3)：建立基于应变描述的PUCK失效初始判据和损伤演化准则，具体建立方式为：
(d)对于纤维拉伸和压缩，损伤初始判据为：
其中，和分别为纤维拉伸和压缩的初始失效应变；所述T*，C分别指拉伸和压缩；所述E11是指纤维方向应变；所述分别指纤维拉伸和压缩失效判断因子；
纤维拉伸和压缩的损伤演化准则为：
其中，所述是指纤维拉伸和压缩损伤变量；所述E11是指纤维方向应变；所述是指纤维损伤变量达到1的纤维临界拉伸和压缩失效应变；所述是指纤维损伤变量为零的纤维初始拉伸和压缩失效应变；
(e)对于基体拉伸损伤失效初始判据为：
其中，所述指基体拉伸失效判断因子；所述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度；
基体压缩损伤演化准则为：
其中，所述γγ是指断裂面联合剪切应变；和是联合剪切应变的初始和最大应变，γ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，塑性加载出现，根据后向欧拉隐式算法(采用Newton-Raphson迭代)，实现试应力到屈服面的最近点返回，将试应力赋予步骤(4)中的迭代初始条件；变量是塑性一致因子的增量Δλn+1的函数，运用牛顿-拉夫森算法求解再更新应变和有效应力，直到Fn+1≤0，结束迭代，得到n+1增量步状态变量
步骤(6)：应变和有效应力得到更新后，根据有效应力和应变代入到步骤(3)的PUCK失效准则判断是否出现损伤，如有损伤再通过损伤演化公式获得损伤变量，再根据步骤(1)通过有效应力和损伤变量计算名义应力Sn+1；
所述过程四具体为：将过程一建立的模型主文件和过程三建立的ABAQUS-VUMAT用户子程序联合，使用ABAQUS/EXPLICT方法对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；即完成低速冲击载荷下弹塑性复合材料层合板渐进失效特性的预测。
与现有技术相比，本发明的有益效果是：
本发明利用ABAQUS-VUMAT用户子程序来数值实现所建立的将塑性和损伤联合的弹塑性损伤本构模型，该模型同时考虑塑性和材料性能退化的影响，能准确预测含塑性特征的复合材料在低速冲击下的渐进损伤失效。
对于ABAQUS内嵌的基于HASHIN失效准则的弹性损伤本构，本发明基于更准确的PUCK失效准则的弹塑性损伤本构对含塑性的复合材料层合板的低速冲击载荷下的渐进失效过程的预测更加准确。
附图说明
图1为本发明实施例复合材料层合板在低速冲击下有限元模型图。
图2为本发明对所提出的弹塑性损伤本构模型的VUMAT数值实现流程图。
图3为实施例中复合材料层合板在5J能量低速冲击下冲击力-时间数值模拟结果与实验结果对比示意图。
图4为实施例中复合材料层合板在5J能量低速冲击载荷下中心位移-时间数值模拟结果与实验结果对比示意图。
图5为实施例中复合材料层合板在10J能量低速冲击载荷下冲击力-时间数值模拟结果与实验结果对比示意图。
图6为实施例中复合材料层合板在10J能量低速冲击载荷下中心位移-时间数值模拟结果与实验结果对比示意图。
图7为实施例中复合材料层合板在5J能量低速冲击结束后基体拉伸，基体压缩，层间分层损伤云图。
图8为实施例中复合材料层合板在10J能量低速冲击结束后基体拉伸，基体压缩，层间分层损伤云图。
具体实施方式
首先需要说明的是，本发明是计算机技术在复合材料损伤预测领域的一种应用。在本发明的实现过程中，会涉及到多个软件功能模块的应用。申请人认为，如在仔细阅读申请文件、准确理解本发明的实现原理和发明目的以后，在结合现有公知技术的情况下，本领域技术人员完全可以运用其掌握的软件编程技能实现本发明。凡本发明申请文件提及的均属此范畴，申请人不再一一列举。
下面结合附图与具体实施方式对本发明作进一步详细描述：
在ABAQUS/CAE中建立包括石墨/环氧树脂复合材料层合板，冲锤以及支撑板的低速冲击有限元模型，如图1所示，平板的铺层顺序为平板大小为100×100×2.1mm，密度为ρ＝1600kg/m3，总共分为20层，每层均用减缩积分三维实体单元C3D8R来仿真，每两个纤维方向相邻的材料层考虑为一个子材料层(0.21mm)，每个子材料层在厚度方向铺设一个单元，内聚力单元层铺设在纤维方向不同的材料层之间研究分层损伤，冲锤前端建成半球形，冲锤的质量为6.5kg，直径为1.27cm。冲锤的冲击能量分别为5J，10J。
层合板上下各有一个正方形支撑框架模型，框架的每边宽度为10mm，厚度为2mm，冲锤和支撑框架的刚度和泊松比分别为207.6GPa和0.26，密度为ρ＝7830kg/m3下框架的下表面固定，上支撑框架的上表面施加350kPa的压力。
利用ABAQUS/EXPLICT计算模拟冲锤冲击层合板的过程，如图2利用用户子程序VUMAT首先判别材料点是否进入塑性，得到有效应力，根据有效应力以及应变判别单元是否进入损伤，进入损伤阶段则按损伤演化准则计算损伤变量，从而得到名义应力。
图3和图4分别为在5J冲击能量下冲击力-时间和中心位移-时间的曲线图，图5和图6分别为在10J冲击能量下的冲击力-时间和中心位移-时间的曲线图，均与试验值十分准确的吻合。在冲击力-时间曲线中冲击力为0时对应的中心位移为永久中心位移，在中心位移-时间的最低点为冲击过程中的最大位移，可见，层合板在冲击过程中的最大位移和永久位移均与试验值吻合的很好。所以本发明提出的弹塑性损伤本构模型可以准确的捕捉有塑性特征的复合材料层合板在低速冲击下的渐进损伤失效特征。
本发明在ABAQUS软件的基础上进行用户子程序的开发，提出的弹塑性损伤本构模型同时考虑塑性和材料性能退化的影响，能准确预测具有塑性特征的复合材料在低速冲击下的渐进损伤失效，为深入阐明复合材料结构的损伤失效特性、提升轻量化强度设计水平提供了技术支撑。
最后，需要注意的是，以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有很多变形。本领域的普通技术人员能从本发明公开的内容中直接导出或联想到的所有变形，均应认为是本发明的保护范围。
Claims (1)
1.预测低速冲击下复合材料层合板渐进失效的有限元方法，其特征在于，包括下述过程：
一、建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型；
二、建立复合材料弹塑性损伤本构模型；
三、基于ABAQUS-VUMAT用户动态材料子程序模块，运用后向欧拉算法实现提出的弹塑性损伤本构模型，求解应力、应变和损伤；
四、对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；
所述过程一中，建立含冲锤、复合材料层合板以及支撑板的低速冲击有限元模型：基于ABAQUS建立冲锤、复合材料层合板、支撑板部件，分别设置材料属性和划分网格，再进行ASSEMBLY对其组装之后设置分析步和通用接触属性；
其中，复合材料层合板每层根据铺层角度进行铺层，不同铺层之间设置ABAQUS自带的双线性内聚力单元；
所述过程二具体包括下述步骤：
步骤(1)：建立含损伤复合材料层合板层内本构关系；
复合材料应力-应变本构方程：S＝Cd：Ee；
其中，为未损伤材料的有效应力；S为损伤材料的名义应力；Ep为格林-拉格朗日应变张量的塑性部分，Ee＝E-Ep为格林-拉格朗日应变张量的弹性部分；Cd(C，d1，d2)为含损伤材料的四阶弹性张量，C为未损伤材料的四阶弹性张量，d1和d2是分别对应纤维断裂和基体开裂的损伤变量；
步骤(2)：建立塑性模型，具体建立方式为：
(a)考虑面外应力的复合材料塑性流动准则为：
<mrow>
<mi>F</mi>
<mrow>
<mo>(</mo>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mo>,</mo>
<msup>
<mover>
<mi>E</mi>
<mo>~</mo>
</mover>
<mi>p</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msqrt>
<mrow>
<mn>3</mn>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
</msqrt>
<mo>-</mo>
<mi>k</mi>
<mrow>
<mo>(</mo>
<msup>
<mover>
<mi>E</mi>
<mo>~</mo>
</mover>
<mi>p</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>0</mn>
<mo>;</mo>
</mrow>
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mrow>
<mi>i</mi>
<mi>j</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mo>&lsqb;</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mn>22</mn>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mn>33</mn>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<mn>4</mn>
<msubsup>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mn>23</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<mn>2</mn>
<msub>
<mi>a</mi>
<mn>66</mn>
</msub>
<mrow>
<mo>(</mo>
<msubsup>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mn>13</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mn>12</mn>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>;</mo>
</mrow>
其中，F为屈服方程；a66为描述各向异性材料与塑性相关的常数，由偏轴拉伸测试确定；f为塑性势函数；k和分别为塑性硬化应力和等效塑性应变；所述是指未损伤材料的有效应力，其中i，j用来确定应力方向；所述指层合板面内垂直于纤维方向的有效应力；所述指面外垂直于纤维的有效应力；所述指面外剪切有效应力；所述指面内剪切有效应力；
(b)塑性硬化应力k为：
<mrow>
<mi>k</mi>
<mo>=</mo>
<mi>&beta;</mi>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mover>
<mi>E</mi>
<mo>~</mo>
</mover>
<mi>p</mi>
</msup>
<mo>)</mo>
</mrow>
<mi>n</mi>
</msup>
<mo>;</mo>
</mrow>
其中，β和n是为满足实验硬化曲线的常数；塑性变形假定发生在损伤材料的未损伤区域，塑性流动准则和硬化准则表达在有效应力空间中；
(c)等效塑性应变率和塑形应变率分别为：
<mrow>
<msup>
<mover>
<mover>
<mi>E</mi>
<mo>~</mo>
</mover>
<mo>&CenterDot;</mo>
</mover>
<mi>p</mi>
</msup>
<mo>=</mo>
<mover>
<mi>&lambda;</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>;</mo>
</mrow>
<mrow>
<msup>
<mover>
<mi>E</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>p</mi>
</msup>
<mo>=</mo>
<mover>
<mi>&lambda;</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>&part;</mo>
<mi>F</mi>
<mo>/</mo>
<mo>&part;</mo>
<mover>
<mi>S</mi>
<mo>&OverBar;</mo>
</mover>
<mo>;</mo>
</mrow>
其中，为一致塑性因子；
步骤(3)：建立基于应变描述的PUCK失效初始判据和损伤演化准则，具体建立方式为：
(d)对于纤维拉伸和压缩，损伤初始判据为：
<mrow>
<msubsup>
<mi>F</mi>
<mn>11</mn>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<msub>
<mi>E</mi>
<mn>11</mn>
</msub>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mn>1</mn>
<mo>&GreaterEqual;</mo>
<mn>0</mn>
<mo>,</mo>
<msubsup>
<mi>F</mi>
<mn>11</mn>
<mi>C</mi>
</msubsup>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<msub>
<mi>E</mi>
<mn>11</mn>
</msub>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mi>C</mi>
</msubsup>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mn>1</mn>
<mo>&GreaterEqual;</mo>
<mn>0</mn>
<mo>;</mo>
</mrow>
其中，和分别为纤维拉伸和压缩的初始失效应变；所述T*，C分别指拉伸和压缩；所述E11是指纤维方向应变；所述分别指纤维拉伸和压缩失效判断因子；
纤维拉伸和压缩的损伤演化准则为：
<mrow>
<msubsup>
<mi>d</mi>
<mn>11</mn>
<mrow>
<mi>T</mi>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mi>C</mi>
<mo>)</mo>
</mrow>
</mrow>
</msubsup>
<mo>=</mo>
<mfrac>
<msubsup>
<mi>E</mi>
<mrow>
<mi>f</mi>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mi>C</mi>
<mo>)</mo>
</mrow>
</mrow>
</msubsup>
<mrow>
<msubsup>
<mi>E</mi>
<mrow>
<mi>f</mi>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mi>C</mi>
<mo>)</mo>
</mrow>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mi>C</mi>
<mo>)</mo>
</mrow>
</mrow>
</msubsup>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
<mrow>
<mo>(</mo>
<mi>C</mi>
<mo>)</mo>
</mrow>
</mrow>
</msubsup>
<msub>
<mi>E</mi>
<mn>11</mn>
</msub>
</mfrac>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
其中，所述是指纤维拉伸和压缩损伤变量；所述E11是指纤维方向应变；所述是指纤维损伤变量达到1的纤维临界拉伸和压缩失效应变；所述是指纤维损伤变量为零的纤维初始拉伸和压缩失效应变；
(e)对于基体拉伸损伤失效初始判据为：
<mrow>
<msubsup>
<mi>F</mi>
<mn>22</mn>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<msub>
<mi>E</mi>
<mn>22</mn>
</msub>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>2</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mn>1</mn>
<mo>&GreaterEqual;</mo>
<mn>0</mn>
<mo>;</mo>
</mrow>
其中，所述指基体拉伸失效判断因子；所述E22是指基体方向应变；所述是指基体损伤变量为零的基体初始拉伸失效应变；
基体拉伸损伤演化准则为：
<mrow>
<msubsup>
<mi>d</mi>
<mn>22</mn>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<mo>=</mo>
<mfrac>
<msubsup>
<mi>E</mi>
<mrow>
<mi>f</mi>
<mo>,</mo>
<mn>2</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<mrow>
<msubsup>
<mi>E</mi>
<mrow>
<mi>f</mi>
<mo>,</mo>
<mn>2</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>2</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<msubsup>
<mi>E</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mn>2</mn>
</mrow>
<mrow>
<mi>T</mi>
<mo>*</mo>
</mrow>
</msubsup>
<msub>
<mi>E</mi>
<mn>22</mn>
</msub>
</mfrac>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
其中，为基体损伤变量达到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度；
基体压缩损伤演化准则为：
<mrow>
<msubsup>
<mi>d</mi>
<mn>22</mn>
<mi>C</mi>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mi>&gamma;</mi>
<mrow>
<mi>N</mi>
<mi>T</mi>
</mrow>
</msub>
<mo>,</mo>
<msub>
<mi>&gamma;</mi>
<mrow>
<mi>N</mi>
<mi>L</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<msubsup>
<mi>&gamma;</mi>
<mi>&gamma;</mi>
<mi>max</mi>
</msubsup>
<mrow>
<msubsup>
<mi>&gamma;</mi>
<mi>&gamma;</mi>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&gamma;</mi>
<mi>&gamma;</mi>
<mi>f</mi>
</msubsup>
</mrow>
</mfrac>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<msup>
<mi>&gamma;</mi>
<mi>f</mi>
</msup>
<msub>
<mi>&gamma;</mi>
<mi>&gamma;</mi>
</msub>
</mfrac>
<mo>)</mo>
</mrow>
<mo>,</mo>
<msub>
<mi>&gamma;</mi>
<mi>&gamma;</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<msubsup>
<mi>&gamma;</mi>
<mrow>
<mi>N</mi>
<mi>T</mi>
</mrow>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&gamma;</mi>
<mrow>
<mi>N</mi>
<mi>L</mi>
</mrow>
<mn>2</mn>
</msubsup>
</mrow>
</msqrt>
<mo>.</mo>
<mo>;</mo>
</mrow>
其中，所述γγ是指断裂面联合剪切应变；和是联合剪切应变的初始和最大应变，γ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方法对低速冲击进行计算，进一步获得冲击力、位移、速度和加速度；即完成低速冲击载荷下弹塑性复合材料层合板渐进失效特性的预测。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610833225.4A CN106503292B (zh) | 2016-09-20 | 2016-09-20 | 预测低速冲击下复合材料层合板渐进失效的有限元方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610833225.4A CN106503292B (zh) | 2016-09-20 | 2016-09-20 | 预测低速冲击下复合材料层合板渐进失效的有限元方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106503292A CN106503292A (zh) | 2017-03-15 |
| CN106503292B true CN106503292B (zh) | 2018-04-24 |
Family
ID=58290840
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610833225.4A Active CN106503292B (zh) | 2016-09-20 | 2016-09-20 | 预测低速冲击下复合材料层合板渐进失效的有限元方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106503292B (zh) |
Families Citing this family (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107491600B (zh) * | 2017-08-04 | 2021-01-19 | 华北理工大学 | 一种优化冲裁工艺参数的方法 |
| CN107832560B (zh) * | 2017-11-29 | 2021-03-02 | 北京航空航天大学 | 一种全SiC复合材料多钉连接结构失效分析方法 |
| CN108427014B (zh) * | 2018-01-02 | 2020-08-28 | 中国商用飞机有限责任公司北京民用飞机技术研究中心 | 一种对复合材料层合板的撞击位置识别方法 |
| CN108416084B (zh) * | 2018-01-23 | 2022-02-18 | 南京理工大学 | 考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元方法 |
| CN108763635B (zh) * | 2018-04-19 | 2019-03-19 | 中国水产科学研究院渔业工程研究所 | 一种基于文本挖掘的三维网格分离仿真系统及方法 |
| CN109558692B (zh) * | 2018-12-25 | 2019-07-12 | 中国石油大学(华东) | 预测微粒冲击金属件残余应力和马氏体相变的有限元方法 |
| CN109946006B (zh) * | 2019-01-17 | 2019-09-24 | 中国石油大学(华东) | 基于混合硬化模型的微粒流冲击金属材料力学行为预测方法 |
| CN109708969B (zh) * | 2019-02-27 | 2020-05-01 | 西北工业大学 | 一种确定金属材料各向异性和拉压非对称性特征的方法 |
| CN110175396B (zh) * | 2019-05-24 | 2020-10-13 | 吉林大学 | 一种基于统一应力准则的粘接结构断裂失效分析方法 |
| CN110308201B (zh) * | 2019-07-22 | 2023-08-18 | 西安工程大学 | 一种基于磁性的叠层复合材料的损伤检测方法 |
| CN111368389B (zh) * | 2019-10-11 | 2023-11-07 | 暨南大学 | 一种预测复合材料层合板失效强度的方法 |
| CN110909498B (zh) * | 2019-11-15 | 2024-04-05 | 上海交通大学 | 带孔复合材料层合板分层损伤与力学行为的精确预测方法 |
| CN111353228A (zh) * | 2020-02-28 | 2020-06-30 | 山东大学 | 一种复合材料层合板冲击响应建模方法 |
| CN111651869B (zh) * | 2020-05-15 | 2022-03-08 | 西北工业大学 | 一种面向高速切削加工的复合材料塑性本构建模方法 |
| CN111709174B (zh) * | 2020-06-18 | 2024-04-09 | 江南大学 | 一种基于失效面理论的复合材料层合板强度分析方法 |
| CN112329205B (zh) * | 2020-10-12 | 2022-04-29 | 湖北航天技术研究院总体设计所 | 一种复合材料结构低速冲击损伤确定方法及装置 |
| CN112906263B (zh) * | 2021-01-28 | 2022-06-28 | 天津大学 | 含制孔分层损伤的复合材料层合板强度预测方法 |
| CN114486518A (zh) * | 2021-12-31 | 2022-05-13 | 中国航空工业集团公司西安飞机设计研究所 | 一种结构复合材料选用效果评估方法 |
| CN115906570B (zh) * | 2022-11-22 | 2023-07-11 | 大连理工大学 | 一种基于数据驱动的复合材料开孔层合板渐进损伤预测方法 |
| CN117131612B (zh) * | 2023-10-26 | 2024-01-16 | 浙江大学 | 车辆倾覆预测方法、系统和电子设备 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005078556A (ja) * | 2003-09-03 | 2005-03-24 | Yokohama Rubber Co Ltd:The | 構造体の有限要素モデル作成方法および構造体のシミュレーション方法 |
| JP2011053807A (ja) * | 2009-08-31 | 2011-03-17 | Fuji Heavy Ind Ltd | 衝突解析装置及び衝突解析方法 |
| CN105701312A (zh) * | 2015-12-17 | 2016-06-22 | 南京航空航天大学 | 复杂编织结构陶瓷基复合材料疲劳迟滞行为预测方法 |
| CN105740566A (zh) * | 2016-02-18 | 2016-07-06 | 浙江大学 | 一种预测层状复合材料层内损伤和层间分层的有限元方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2009190427A (ja) * | 2008-02-12 | 2009-08-27 | Yokohama Rubber Co Ltd:The | タイヤのシミュレーション方法 |
-
2016
- 2016-09-20 CN CN201610833225.4A patent/CN106503292B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005078556A (ja) * | 2003-09-03 | 2005-03-24 | Yokohama Rubber Co Ltd:The | 構造体の有限要素モデル作成方法および構造体のシミュレーション方法 |
| JP2011053807A (ja) * | 2009-08-31 | 2011-03-17 | Fuji Heavy Ind Ltd | 衝突解析装置及び衝突解析方法 |
| CN105701312A (zh) * | 2015-12-17 | 2016-06-22 | 南京航空航天大学 | 复杂编织结构陶瓷基复合材料疲劳迟滞行为预测方法 |
| CN105740566A (zh) * | 2016-02-18 | 2016-07-06 | 浙江大学 | 一种预测层状复合材料层内损伤和层间分层的有限元方法 |
Non-Patent Citations (2)
| Title |
|---|
| SiC纤维增强复合材料界面破坏与失效机理的研究;刘鹏飞;《中国优秀博硕士学位论文全文数据库 (博士) 工程科技Ⅰ辑》;20070115(第1期);第B020-18页 * |
| 复合材料氢气瓶快充过程温升控制方法研究;刘延雷 等;《太阳能学报》;20120930;第33卷(第9期);第1621-1627页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106503292A (zh) | 2017-03-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106503292B (zh) | 预测低速冲击下复合材料层合板渐进失效的有限元方法 | |
| CN106777769B (zh) | 预测低速冲击下复合材料多层厚板渐进失效的有限元方法 | |
| Faggiani et al. | Predicting low-velocity impact damage on a stiffened composite panel | |
| Kim et al. | Composite damage model based on continuum damage mechanics and low velocity impact analysis of composite plates | |
| Xiao et al. | Progressive damage and delamination in plain weave S-2 glass/SC-15 composites under quasi-static punch-shear loading | |
| Li et al. | Repeated low-velocity impact response and damage mechanism of glass fiber aluminium laminates | |
| Miao et al. | Mechanical behaviors and equivalent configuration of a polyurea under wide strain rate range | |
| Johnson et al. | Influence of delamination on impact damage in composite structures | |
| Cheng et al. | Progressive damage behaviors of woven composite laminates subjected to LVI, TAI and CAI | |
| CN113420376B (zh) | 基于多尺度的碳纤维复合材料抗冲击力学性能仿真方法 | |
| Singh et al. | Reduced order multiscale modeling of fiber reinforced polymer composites including plasticity and damage | |
| Kam et al. | Quasi-static buckling and first-ply failure loads of shear web reinforced glass-fabric composite wind blades | |
| Waas et al. | Prediction of low-velocity face-on impact response and compression after impact (cai) of composite laminates using est and cohesive modeling (dczm) | |
| Minak et al. | Low-velocity impact on laminates | |
| Han et al. | Behaviors of composite laminates under low-energy impact using a novel analytical framework | |
| Pickett et al. | Test and modelling of impact on pre-loaded composite panels | |
| Joglekar et al. | Validation of an efficient finite element analysis approach for simulation of low velocity impact and compression strength after impact response | |
| Yu et al. | The effect of matrix shear strength on the out-of-plane compressive strength of CFRP cross-ply laminates | |
| Li et al. | Characterization of strength of carbon fiber reinforced polymer composite based on micromechanics | |
| CN112926244A (zh) | 一种复合材料层合板开孔件极限载荷确定方法 | |
| CN116882232A (zh) | 低速冲击下的碳纤维复合材料层合板损伤状态的预测方法 | |
| Sun et al. | Study on impact resistance and parameter optimization of patch‐repaired plain woven composite based on multi‐scale analysis | |
| Ji et al. | Study of in-plane fatigue failure and life prediction of weave composites under constant and variable amplitude loading | |
| Sajedi et al. | Numerical investigation on the effect of concrete-FRP bond on the flexural behavior of RC beams | |
| Zhuang et al. | A progressive damage model for carbon fiber-reinforced polymer laminates subjected to fastener pull-through failure |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |