一种用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法

摘要

本发明涉及一种用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法，首先建立了平面单元特征长度的计算模型，在此基础上建立了三维平行六面体单元的特征长度计算方法，再通过对一般形状单元和楔形单元进行等体积转换，将其特征长度计算问题转化为平行六面体单元的特征长度计算问题。本发明用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法考虑了复合材料结构有限元模型中常用的典型单元构型，包括平行四边形单元，三角形单元，平行六面体单元，一般形状单元及楔形单元，建立了单元特征长度计算模型，提供了准确的单元特征长度计算公式，能够有效地计算单元特征长度，从而计算单元损伤的耗散能量，提高了复合材料结构失效有限元模拟的客观性和准确性。

说明书

技术领域

本发明涉及复合材料结构强度和失效的有限元分析和模拟，适用于工程领域中各种复合 材料结构有限元损伤模拟模型。

背景技术

随着现代航空航天飞行器对材料性能要求的提高，复合材料在航空航天飞行器结构中的 用量越来越高。因此，发展复合材料结构分析技术日益重要。随着计算机和有限元技术的迅 猛发展，以及复合材料在主结构上的逐步应用，复合材料结构的应力分析和强度预测越来越 多地依赖于先进的有限元和数值分析技术。目前广泛应用的方法有裂纹闭合技术（VCCT）， 渐进损伤分析方法（PDM）及内聚力模型（CZM）等方法，这些方法或基于断裂力学、或 基于损伤力学，或基于断裂-损伤力学，并结合先进的有限元分析技术，为深入理解复合材 料结构的失效机理和有效预测复合材料结构破坏提供了有效的途径。

基于损伤力学方法应用有限元分析技术时，面临着损伤局部化问题，若处理不当，将导 致预测结果失真。损伤局部化是指由于连续损伤力学模型中应力软化段的存在，在数值模拟 中损伤会局限在一条窄带内，而由于损伤的数值模拟结果与网格的精细度密切相关，随着有 限元模型网格尺寸的减小，数值计算的结构损伤耗散能量下降，模型预测的结构强度会下降。 为克服损伤局部化，提出了裂纹带模型，假设材料损伤形成的宏观裂纹可以用一层单 元的损伤区域来代替，将损伤裂纹的断裂韧度在整个单元体积内均匀化。裂纹带模型的本质 是将材料损伤过程中耗散的总能量在整个单元体积（微元体）内均匀化，使得有限元模拟过 程中材料损伤耗散的能量恒定。而单元特征长度是裂纹带模型中的重要参数，它不仅与单元 的几何形状有关，而且与材料的属性和裂纹的方向有关，准确地计算单元特征长度是保证准 确地计算耗散能量的前提条件。因此，有必要建立一般单元形状的特征长度的计算方法。

发明内容

本发明要解决的技术问题为：克服现有技术的不足，考虑了复合材料结构有限元模型中 常用的典型单元构型，包括平行四边形单元，三角形单元，平行六面体单元，一般形状单元 及楔形单元，建立了计算其特征长度的方法，该方法能够有效地计算单元特征长度，为裂纹 带模型提供准确的单元特征长度值。

本发明解决上述技术问题采用的技术方案为：一种用于复合材料结构失效有限元模拟中 单元损伤耗散能量的估计方法，该估计方法的步骤如下：

步骤一、裂纹带模型的本质是将材料损伤过程中耗散的总能量在整个单元体积即微元体 内均匀化，使得有限元模拟过程中材料损伤耗散的能量恒定，而单元特征长度是裂纹带模型 中的重要参数，它不仅与单元的几何形状有关，而且与材料的属性和裂纹的方向有关，准确 地计算单元特征长度是保证准确地计算耗散能量的前提条件，复合材料结构失效有限元模拟 中单元特征长度计算方法分为平面单元特征长度计算和三维单元特征长度计算，其中：

平面单元特征长度计算步骤如下：

设裂纹与单元一边的夹角为θ^crc，裂纹长度为l_crc，单元相邻两边的夹角为边长分别 为a和b，与AB边对应的距离为

步骤A1，首先假设一条裂纹穿过一层由n个平行四边形单元或三角形单元组成的区域；

步骤B1，根据单元特征长度的定义，可得：

$l^{*}=\frac{\mathrm{nah}}{l_{\mathrm{crc}}}=\frac{\mathrm{nah}}{h/\sin \left({\theta }^{\mathrm{crc}}\right)}=\mathrm{na}\sin \left({\theta }^{\mathrm{crc}}\right);$

步骤C1，当时，近似地认为n个平行四边形单元的面积为直角 三角形AFE的2倍，则上式可简化为：

步骤D1，时，平行四边形单元或三角形单元的特征长度为

三维单元特征长度计算步骤如下：

步骤A2，首先将六面体与纤维平行的平面按照平行四边形特征长度等效转化方法转换 为两边与纤维平行的矩形E＇F＇G＇H＇，然后将六面体的侧面四条边在平面E＇F＇G＇H＇ 的垂线上进行投影，得到一立方体；

步骤B2，上述立方体A＂B＂C＂D＂E＇F＇G＇H＇各边的长度为：

$\overline{E^{\prime }{H}^{\prime }}=\left(\begin{array}{ccc}\overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})\cos \left({\theta }^f\right)& \mathrm{if}& \tan (\angle \mathrm{HEF})\le \overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})/\overline{\mathrm{EF}}\\ \overline{\mathrm{EF}}\sin (\angle \mathrm{HEF})\sin \left({\theta }^f\right)& \mathrm{if}& \tan (\angle \mathrm{HEF})>\overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})/\overline{\mathrm{EF}}\end{array}\right)$

$\overline{\mathrm{EF}}=S_{\mathrm{EFGH}}/\overline{E^{\prime }{H}^{\prime }}$

$\overline{A^{\prime \prime }{E}^{\prime }}=|\overrightarrow{\mathrm{AE}}·{\overrightarrow{e}}_{(E\bar{F}\times E\bar{H})}|$

式中，分别代表相应线段的长度，∠HEF、∠HEF代表角度，代表 用AE表示的向量，代表单位向量；

步骤C2，由立方体A＂B＂C＂D＂E＇F＇G＇H＇易知，对应于纤维“裂纹”的单元特 征长度为：

$l_{\mathrm{fbr}}^{*}=\overline{E^{\prime }{F}^{\prime }};$

步骤D2，设基体裂纹面与等效立方体单元的E＇F＇G＇H＇面所成夹角为则对应 于基体裂纹的单元特征长度为：

$l_{\mathrm{mtr}}^{*}=\left(\begin{array}{ccc}\overline{A^{\prime }{E}^{\prime }}\cos \left({\theta }_m^{\mathrm{crc}}\right)& \mathrm{if}& \tan \left({\theta }_m^{\mathrm{crc}}\right)\le \overline{A^{\prime \prime }{E}^{\prime }}/\overline{E^{\prime }{H}^{\prime }}\\ \overline{E^{\prime }{H}^{\prime }}\sin \left({\theta }_m^{\mathrm{crc}}\right)& \mathrm{if}& \tan \left({\theta }_m^{\mathrm{crc}}\right)\overline{A^{\prime \prime }{E}^{\prime }}/\overline{E^{\prime }{H}^{\prime }}\end{array}\right)$

其中基体裂纹面与立方体单元边的相对位置可以通过裂纹面在模型整体坐标系中的方 向坐标和平面E′F′H′G′的法向坐标求得；

步骤E2，对于形状更为一般的六面体单元，可以首先将其等体积转化为平行六面体单 元，其中平行六面体单元的底面由四边形EFGH等面积转化来确定，进而采用平行六面体的 计算公式求得单元的特征长度；

步骤F2，对于楔形单元，可以借鉴平面三角形单元特征长度的求解方法，将两个相同 楔形单元组合、扩展为平行六面体后得到其单元特征长度。

步骤二、复合材料结构失效有限元模拟中单元耗散能量计算步骤如下：

步骤A3，在损伤生成阶段，微裂纹分散在整个单元体积内，随着载荷的增大，微裂纹 生长、贯通，最终形成宏观裂纹。因此，单元在损伤（微裂纹）形成、扩展过程中耗散的总 能量W_D就是细观裂纹形成消耗的断裂能W_C，即：

W_D＝W_C；

步骤B3，假设材料的临界应变能释放率G_C恒定，则在单元内材料生成面积为A的裂纹 消耗的断裂能W_C为：

W_C＝G_C*A；

步骤C3，基于上述单元特征长度计算方法得到的特征长度，由W_D＝W_C，得单元的损伤 耗散能密度为：

$\Gamma =\frac{G_c}{V/A}=\frac{G_c}{l^{*}}$

其中，G_C为材料的临界应变能释放率，V为单元的体积，A为与宏观裂纹平行的截面积， l^*为单元的特征长度，由上述单元特征长度计算方法得到；

步骤D3，基于上面得到的单元损伤耗散能密度，最终得到从微裂纹产生、扩展到宏观 裂纹形成过程中，单元损伤所耗散的能量W_D为：

W_D＝Γ*V。

本发明与现有技术相比的优点在于：

（1）、本发明用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法，首先 建立了平面单元特征长度计算模型，进而得到其计算公式，在此基础上建立了三维平行六面 体单元的特征长度近似计算方法，再通过对一般形状单元和楔形单元进行等体积转换，将其 特征长度计算问题转化为平行六面体单元的特征长度计算问题。

（2）、本发明用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法，考虑 了复合材料结构有限元模型中常用的典型单元构型，包括平行四边形单元，三角形单元，平 行六面体单元，一般形状单元及楔形单元，能够有效地计算单元特征长度，为裂纹带模型提 供准确的单元特征长度值。

附图说明

图1是复合材料典型三维单元构型，其中图1（a）和（b）是六面体单元及材料方向， 构型（b）更具有一般性，图1（c）和（d）是楔形单元及材料方向；

图2是单元特征长度的计算原理图；

图3是平面单元与裂纹示意图，其中图3（a）为平行四边形单元与裂纹示意图，图3（b） 为三角形单元与裂纹示意图；

图4是平行六面体等体积转换示意图；

图5是一般形状六面体转换示意图；

图6是楔形单元扩展为六面体示意图。

具体实施方式

本发明用于复合材料结构失效有限元模拟中单元损伤耗散能量的估计方法，该估计方法 的具体步骤如下：

步骤一、裂纹带模型的本质是将材料损伤过程中耗散的总能量在整个单元体积即微元体 内均匀化，使得有限元模拟过程中材料损伤耗散的能量恒定，而单元特征长度是裂纹带模型 中的重要参数，它不仅与单元的几何形状有关，而且与材料的属性和裂纹的方向有关，准确 地计算单元特征长度是保证准确地计算耗散能量的前提条件，复合材料结构失效有限元模拟 中单元特征长度计算方法分为平面单元特征长度计算和三维单元特征长度计算，其中：

平面单元特征长度计算步骤如下：

设裂纹与单元一边的夹角为θ^crc，裂纹长度为l_crc，单元相邻两边的夹角为边长分别 为a和b，与AB边对应的距离为首先假设一条裂纹穿过一层由n个平行四边形 单元或三角形单元组成的区域，根据单元特征长度的定义，可得：

$l^{*}=\frac{\mathrm{nah}}{l_{\mathrm{crc}}}=\frac{\mathrm{nah}}{h/\sin \left({\theta }^{\mathrm{crc}}\right)}=\mathrm{na}\sin \left({\theta }^{\mathrm{crc}}\right)$

当时，近似地认为n个平行四边形单元的面积为直角三角形AFE的2 倍，则上式可简化为：

当时，平行四边形单元或三角形单元的特征长度为：

三维单元特征长度计算步骤如下：

首先将六面体与纤维平行的平面按照平行四边形特征长度等效转化方法转换为两边与 纤维平行的矩形E＇F＇G＇H＇，然后将六面体的侧面四条边在平面E＇F＇G＇H＇的垂线上 进行投影，得到如图4所示的立方体A＂B＂C＂D＂E＇F＇G＇H＇，其各边的长度为：

$\overline{E^{\prime }{H}^{\prime }}=\left(\begin{array}{ccc}\overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})\cos \left({\theta }^f\right)& \mathrm{if}& \tan (\angle \mathrm{HEF})\le \overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})/\overline{\mathrm{EF}}\\ \overline{\mathrm{EF}}\sin (\angle \mathrm{HEF})\sin \left({\theta }^f\right)& \mathrm{if}& \tan (\angle \mathrm{HEF})>\overline{\mathrm{EH}}\sin (\angle \mathrm{HEF})/\overline{\mathrm{EF}}\end{array}\right)$

$\overline{\mathrm{EF}}=S_{\mathrm{EFGH}}/\overline{E^{\prime }{H}^{\prime }}$

$\overline{A^{\prime \prime }{E}^{\prime }}=|\overrightarrow{\mathrm{AE}}·{\overrightarrow{e}}_{(E\bar{F}\times E\bar{H})}|$

式中，分别代表相应线段的长度，∠HEF、∠HEF代表角度，代表用AE 表示的向量，代表单位向量。由立方体A＂B＂C＂D＂E＇F＇G＇H＇易知，对应于纤维 “裂纹”的单元特征长度为：

$l_{\mathrm{fbr}}^{*}=\overline{E^{\prime }{F}^{\prime }}$

设基体裂纹面与等效立方体单元的E＇F＇G＇H＇面所成夹角为则对应于基体裂 纹的单元特征长度为：

$l_{\mathrm{mtr}}^{*}=\left(\begin{array}{ccc}\overline{A^{\prime }{E}^{\prime }}\cos \left({\theta }_m^{\mathrm{crc}}\right)& \mathrm{if}& \tan \left({\theta }_m^{\mathrm{crc}}\right)\le \overline{A^{\prime \prime }{E}^{\prime }}/\overline{E^{\prime }{H}^{\prime }}\\ \overline{E^{\prime }{H}^{\prime }}\sin \left({\theta }_m^{\mathrm{crc}}\right)& \mathrm{if}& \tan \left({\theta }_m^{\mathrm{crc}}\right)\overline{A^{\prime \prime }{E}^{\prime }}/\overline{E^{\prime }{H}^{\prime }}\end{array}\right)$

其中基体裂纹面与立方体单元边的相对位置可以通过裂纹面在模型整体坐标系中的方向坐 标和平面E′F′H′G′的法向坐标求得。对于形状更为一般的六面体单元，可以首先将其等体积 转化为平行六面体单元，其中平行六面体单元的底面由四边形EFGH等面积转化来确定，进 而采用平行六面体的计算公式求得单元的特征长度，其等体积转化过程如图5所示；对于楔 形单元，可以借鉴平面三角形单元特征长度的求解方法，将两个相同楔形单元组合、扩展为 平行六面体后得到其单元特征长度，楔形单元转化为平行六面体的过程如图6所示。

步骤二、复合材料结构失效有限元模拟中单元耗散能量计算步骤如下：

在损伤生成阶段，微裂纹分散在整个单元体积内，随着载荷的增大，微裂纹生长、贯通， 最终形成宏观裂纹。因此，单元在损伤（微裂纹）形成、扩展过程中耗散的总能量WD就是 细观裂纹形成消耗的断裂能W_C，即：

W_D＝W_C

假设材料的临界应变能释放率G_C恒定，则在单元内材料生成面积为A的裂纹消耗的断裂 能W_C为：

W_C＝G_C*A

基于上述单元特征长度计算方法得到的特征长度，由W_D＝W_C，得单元的损伤耗散能密 度为：

$\Gamma =\frac{G_c}{V/A}=\frac{G_c}{l^{*}}$

其中，G_C为材料的临界应变能释放率，V为单元的体积，A为与宏观裂纹平行的截面积，l^*为单元的特征长度，由上述单元特征长度计算方法得到。基于单元损伤耗散能密度，最终得 到从微裂纹产生、扩展到宏观裂纹形成过程中，单元损伤所耗散的能量W_D为：

W_D＝Γ*V。

本发明未详细公开的部分属于本领域的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理 解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员 来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显 而易见的，一切利用本发明构思的发明创造均在保护之列。