CN110245376A - 一种复合材料层裂二维计算方法 - Google Patents

Abstract

本发明公开了一种复合材料层裂二维计算方法，该方法包括以下步骤：S1，采用扩展有限元法求解节点位移；S2，基于虚拟裂纹闭合技术，获得层裂裂纹尖端能量释放率；S3，判断层裂是否扩展；S4，达到预定载荷或结构失效，结束计算。该发明的优点在于：本发明将虚拟裂纹闭合技术和扩展有限元法结合起来，由于采用扩展有限元表征层裂裂纹，网格划分时单元边界无需与裂纹形貌相吻合，模拟层裂扩展时也无需重新划分网格；降低了层裂模拟对单元尺寸的要求，计算结果相对粘聚单元法更加准确；进行层裂扩展分析时，无需迭代，提高了计算效率。

Description

一种复合材料层裂二维计算方法

技术领域

本发明涉及复合材料损伤模拟领域，尤其是一种复合材料层裂二维计算方法。

背景技术

复合材料由于具有高比模量、高比强度、性能易设计等优点，在航空航天、汽车、体育用品、压力容器等领域有着越来越广泛的应用。但复合材料在使用过程中发现了多种损伤模式，包括纤维断裂、基体开裂和层裂等形式。理论和工程经验表明，层裂是复合材料最常见的损伤模式，层裂发生会严重降低复合材料结构强度、刚度和完整性。因此，开展层裂数值分析对复合材料结构设计和强度分析具有重要意义。

常用复合材料层裂数值模拟方法包括内聚力单元、虚拟裂纹闭合技术、扩展有限元等方法。内聚力单元是使用最广泛的方法，但内聚力单元法对网格尺寸的高要求限制了其在大型复合材料结构中的应用；虚拟裂纹闭合技术虽然对网格尺寸要求低，但是要求单元边界与层裂扩展路径吻合，给网格划分带来极大麻烦；扩展有限元中，虽然单元边界无需与层裂扩展路径吻合，但需要结合层裂扩展准则使用。

发明内容

为了克服上述现有技术的不足，为此，本发明提供一种复合材料层裂二维计算方法。

为实现上述目的，本发明采用以下技术方案：

一种复合材料层裂二维计算方法，层裂产生的裂纹包括裂纹面和裂纹尖端，该方法包括以下步骤：

S1，采用扩展有限元法求解节点位移；

S2，基于虚拟裂纹扩展技术，获得层裂裂纹尖端能量释放率；

S3，判断层裂是否扩展；

S4，达到预定载荷或结构失效，结束计算。

详细地说，所述步骤S1包括以下步骤：

S11，对复合材料的几何区域进行网格划分；

S12，采用扩展有限元法表征层裂裂纹，分别表征裂纹面和裂纹尖端，扩展有限元的位移场可由下式表示：

x为点坐标，x_j和x_k为节点坐标，式中u(x)为x点的位移，i,j,k,α分别为求和公式中的循环指数，H(x)为跳跃函数，若x位于裂纹面上方，H(x)＝1；若x位于裂纹面下方，H(x)＝-1，φ_α(x)为奇异增强函数，定义如下：

式中(r,θ)为裂纹尖端局部极坐标；u_i为有限元的节点位移，a_j、分别为跳跃函数、奇异函数增强节点自由度，N_i(x)、N_j(x)、N_k(x)为单元的形函数，N_H、N_α分别为单元内的跳跃函数增强节点数和奇异函数增强节点数，H(x_j)、φ_α(x_k)分别为节点x_j、x_k的跳跃函数值和奇异函数值；

S13，结合扩展有限元表征的裂纹形貌，在网格中确定奇异函数增强节点和跳跃函数增强节点，并对裂纹尖端和裂纹面贯穿的单元进行区域划分；

S14，获得每个单元刚度矩阵，组装形成总体刚度矩阵；

S15，施加载荷和边界条件，求解总体刚度矩阵方程获得节点位移。

详细地说，所述步骤S2包括以下步骤：

S21，虚拟裂纹闭合技术认为裂纹扩展所释放能量等于闭合裂纹时所需的能量，根据裂纹扩展相似性假设，三种模态能量释放率可由以下三式给出，

式中G_I、G_II、G_III分别为I型、II型、III型模态能量释放率，x为虚拟裂纹扩展路径上的点坐标，x-Δa分别为x点沿着裂纹面方向平移Δa后的对应点坐标，δ₁(x-Δa)、δ₂(x-Δa)、δ₃(x-Δa)裂纹面上下对应点张开位移的分量，可由扩展有限元位移u(x)获得，σ₂₂(x)、σ₁₂(x)、σ₁₃(x)分别为x点处应力σ的分量，应力可由下式获得：

σ＝Dε

式中ε为应变张量，可由ε＝Bu给出，其中B为位移-应变关系矩阵，可由单元形状函数N_i(x)，跳跃增强函数H(x)、奇异增强函数φ_α(x)及它们的导数给出。

B_l＝N_l(x)、N_l(x)(H(x)-H(x_l))或N_l(x)(φ_α(x)-φ_α(x_l))，l为矩阵组装指数，N_node为节点数目，包括普通节点和跳跃函数和奇异函数增强节点。u为单元自由度位移向量，包括u_i、a_j和D为材料的刚度矩阵，对于平面应力问题，D可表示为：

式中E为杨氏模量，υ为泊松比，

S22,假设层裂向前扩展Δa，根据虚拟扩展裂纹与网格单元之间的位置关系，对虚拟扩展裂纹Δa进行分段积分，若扩展裂纹Δa位于单个单元内，则无需进行分段；若扩展裂纹Δa经过多个单元，则利用所经过单元对扩展裂纹Δa进行自然分段；同时将分段积分点沿着裂纹面方向平移Δa，获得对应的张开位移分量；

S23，在S21和S22的基础上，对虚拟扩展裂纹各分段分别进行高斯积分并求和，获得层裂裂纹尖端的能量释放率，其数值积分由如下三式给出：

式中N_seg为虚拟扩展裂纹分段数目，N_gau为高斯积分点数，x_mn为第m分段第n高斯积分点坐标，可由裂纹分段两端坐标及高斯积分点位置确定，w_n为第n高斯积分点的权重。

详细地说，所述步骤S3包括以下步骤：

S31,根据算得的能量释放率，判断是否满足层裂扩展条件；

式中G_Ic、G_IIc、G_IIIc分别为I型断裂韧度、II型断裂韧度、III型断裂韧度；

S32，若满足上式，则层裂向前扩展Δa长度，更新裂纹几何尺寸，回到S12步，继续计算；若不满足上式，则增加载荷，继续计算。

本发明的优点在于：

1)本发明将虚拟裂纹闭合技术和扩展有限元法结合起来，由于采用扩展有限元表征层裂裂纹，网格划分时单元边界无需与裂纹形貌相吻合，同时层裂扩展时也无需重新划分网格。

2)相对于虚拟裂纹闭合技术，裂纹扩展长度受网格尺寸限制，即每步裂纹扩展长度为单元长度，本发明由于采用增强函数表征裂纹面和裂纹尖端，所以裂纹扩展长度不受网格尺寸限制。

3)本发明提出了通过分段积分的方法，该方法解决了扩展有限元裂纹尖端能量释放率计算困难的问题。

4)本发明在模拟层裂扩展过程中不涉及材料性能衰退,这一材料非线性因素，因此避免了求解过程中由于材料非线性带来的反复迭代，提高了计算效率。

5)本发明基于虚拟裂纹闭合技术，相较于内聚力单元法，该方法在模拟层裂扩展时，对网格尺寸要求较低。

附图说明

图1为本发明的流程示意图；

图2为本发明中节点分类、裂纹面贯穿单元和裂纹尖端单元区域划分示意图；

图3为扩展有限元裂纹尖端能量释放率计算模型示意图；

图4为实施例中模型示意图；

图5为裂纹扩展方向单元尺寸为2mm的网格示意图；

图6为裂纹扩展方向单元尺寸为1mm的网格示意图；

图7为本方法使用的模型和内聚力单元模型的计算结果对比图；

图8为I型模态、II模态、III型模态的裂纹示意图。

具体实施方式

本发明所采用的实施例为双悬臂梁模型，几何尺寸如图4所示，该模型左端固定，右端上下顶点各施加大小相等、方向相反一对平衡力，大小为P，模型包括两层复合材料层和中间层，中间层为层裂发生层，在层裂发生层预制30mm长的层裂裂纹，根据图8所示，可以判断该裂纹为I型裂纹。复合材料层的材料常数为：纤维方向杨氏模量E₁₁＝150GPa，横向杨氏模量E₂₂＝11GPa，剪切模量G₁₂＝60GPa，泊松比υ₁₂＝0.21。中间层为各向同性材料，其杨氏模量为E＝10GPa,泊松比υ＝0.3，断裂韧性G_IC＝0.026Nmm^-1。

对该实施例采用基于虚拟裂纹闭合技术的扩展有限元方法进行计算，如图1所示，主要步骤如下：

S1，采用扩展有限元法求解节点位移

具体步骤如下：

S11，对复合材料的几何区域进行网格划分；在本方案中分成两种大小的网格，如图5、6所示，图5在层裂扩展方向上单元尺寸为2mm，图6在层裂扩展方向上单元尺寸为1mm。

S12，如图2、3所示，采用扩展有限元法表征层裂裂纹，分别表征裂纹面和裂纹尖端。

S13，结合扩展有限元表征的裂纹形貌，在网格中确定奇异函数增强节点和跳跃函数增强节点，并对裂纹尖端和裂纹面贯穿的单元进行区域划分。

S14，获得每个单元的刚度后，组成单元刚度矩阵，组装总体形成的刚度矩阵；

S15，施加载荷和边界条件，求解获得节点位移。

S2，基于虚拟裂纹扩展技术，获得层裂裂纹尖端能量释放率。

S21，假设裂纹扩展Δa，根据虚拟扩展裂纹与单元之间的位置关系，确定是否需要对虚拟扩展裂纹长度Δa进行分段。在本实例中，需要将Δa分成两段，分别在裂纹扩展Δa的两分段内进行高斯积分，由于该实例中能量释放率只包含I型模态，故只需下式所示的积分即可，

式中σ₂₂(x_mn)为虚拟扩展裂纹上高斯积分点的应力σ分量，应力可由下式获得：

σ＝Dε

式中ε为应变张量，可由ε＝Bu给出，其中B为位移-应变关系矩阵，可由单元形状函数N_i(x)，跳跃增强函数H(x)、奇异增强函数φ_α(x)及它们的导数给出。

B_l＝N_l(x)、N_l(x)(H(x)-H(x_l))或N_l(x)(φ_α(x)-φ_α(x_l))，l为矩阵组装指数，N_node为节点数目，包括普通节点和跳跃函数和奇异函数增强节点。u为单元自由度位移向量，包括u_i、a_j和D为材料的刚度矩阵，对于平面应力问题，D可表示为：

式中E为杨氏模量，υ为泊松比，

δ₂(x_mn-Δa)为裂纹面上高斯积分点对应位置张开位移在垂直方向上的分量，裂纹面上高斯积分点对应位置是通过将高斯积分点向裂纹面方向平移Δa获得，N_seg为虚拟扩展裂纹分段数目，本实例中为2，N_gau为高斯积分点数，w_n为高斯积分点的权重，在所有分段内将上式进行累积求和，最终将高斯积分结果除以2Δa，获得能量释放率G_I。

S3，判断层裂裂纹是否扩展

S31，根据算得的能量释放率，判断是否满足扩展条件，由于该裂纹为I型裂纹，公式

简化成

S32，若满足该式，则裂纹扩展Δa，并更新裂纹几何形貌，回到步骤(2)继续进行计算；若不满足裂纹扩展条件，则继续增加载荷，并返回步骤(5)，继续进行计算。

S4，达到预定载荷或结构失效，结束计算。

计算结果如图7所示，图中横坐标为裂纹张开位移，纵坐标为所施加载荷大小，图中Analsolu代表该问题的解析解，Cohelem2和Cohelem1代表采用内聚力单元法在单元尺寸为2mm和1mm时的计算结果，XFEM-VCCTelem2和XFEM-VCCTelem1代表采用基于虚拟裂纹闭合技术的扩展有限元方法在单元尺寸为2mm和1mm时的计算结果。从图中可以看出曲线分两个阶段，即上升和下降阶段。在上升阶段，采用基于虚拟裂纹闭合技术的扩展有限元方法的结果曲线梯度要比内聚力单元梯度大，这是由于在内聚力单元法中考虑了材料刚度的下降，而本发明的方法中并未考虑，因此结果曲线梯度较大。在下降阶段，无论是单元尺寸为2mm或1mm，本发明所述方法的计算结果均比内聚力单元法更接近解析解，甚至本方法在2mm时产生的结果要比内聚力单元法在1mm的结果更准确。另外，由于内聚力单元法考虑了材料损伤，力学性能下降，会产生锯齿状计算结果，而本发明所述方法不会产生锯齿状结果。再者，由于内聚力单元涉及到材料非线性问题，计算过程中需要大量迭代求解，耗费了大量的计算时间，而本方法无需考虑材料损伤，则没有此类问题。由此可见，相较于内聚力单元法，本方法可以提供更高精度的解，同时所需时间也会大大降低，极大的提高了计算效率问题，克服了内聚力单元法面临的网格尺寸要求高的问题。

以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。

Claims (4)

1.一种复合材料层裂二维计算方法，其特征在于，层裂产生的裂纹包括裂纹面和裂纹尖端，该方法包括以下步骤：

S1，采用扩展有限元法求解节点位移；

S2，基于虚拟裂纹扩展技术，获得层裂裂纹尖端能量释放率；

S3，判断层裂是否扩展；

S4，达到预定载荷或结构失效，结束计算。

2.根据权利要求1所述的一种复合材料层裂二维计算方法，其特征在于，所述步骤S1包括以下步骤：

S11，对复合材料的几何区域进行网格划分；

S12，采用扩展有限元法表征层裂裂纹，分别表征裂纹面和裂纹尖端，扩展有限元的位移场可由下式表示：

x为点坐标，x_j和x_k为节点坐标，式中u(x)为x点的位移，i,j,k,α分别为求和公式中的循环指数，H(x)为跳跃函数，若x位于裂纹面上方，H(x)＝1；若x位于裂纹面下方，H(x)＝-1，φ_α(x)为奇异增强函数，定义如下：

式中(r,θ)为裂纹尖端局部极坐标；u_i为有限元的节点位移，a_j、分别为跳跃函数、奇异函数增强节点自由度，N_i(x)、N_j(x)、N_k(x)为单元的形函数，N_H、N_α分别为单元内的跳跃函数增强节点数和奇异函数增强节点数，H(x_j)、φ_α(x_k)分别为节点x_j、x_k的跳跃函数值和奇异函数值；

S13，结合扩展有限元表征的裂纹形貌，在网格中确定奇异函数增强节点和跳跃函数增强节点，并对裂纹尖端和裂纹面贯穿的单元进行区域划分；

S14，获得每个单元刚度矩阵，组装形成总体刚度矩阵；

S15，施加载荷和边界条件，求解总体刚度矩阵方程获得节点位移。

3.根据权利要求2所述的一种复合材料层裂二维计算方法，其特征在于，所述步骤S2包括以下步骤：

S21，虚拟裂纹闭合技术认为裂纹扩展所释放能量等于闭合裂纹时所需的能量，根据裂纹扩展相似性假设，三种模态能量释放率可由以下三式给出，

式中G_I、G_II、G_III分别为I型、II型、III型模态能量释放率，x为虚拟裂纹扩展路径上的点坐标，x-Δa分别为x点沿着裂纹面方向平移Δa后的对应点坐标，δ₁(x-Δa)、δ₂(x-Δa)、δ₃(x-Δa)裂纹面上下对应点张开位移的分量，可由扩展有限元位移u(x)获得，σ₂₂(x)、σ₁₂(x)、σ₁₃(x)分别为x点处应力σ的分量，应力可由下式获得：

σ＝Dε

式中ε为应变张量，可由ε＝Bu给出，其中B为位移-应变关系矩阵，可由单元形状函数N_i(x)，跳跃增强函数H(x)、奇异增强函数φ_α(x)及它们的导数给出；

B_l＝N_l(x)、N_l(x)(H(x)-H(x_l))或N_l(x)(φ_α(x)-φ_α(x_l))，l为矩阵组装指数，N_node为节点数目，包括普通节点和跳跃函数和奇异函数增强节点；u为单元自由度位移向量，包括u_i、a_j和D为材料的刚度矩阵，对于平面应力问题，D可表示为：

式中E为杨氏模量，υ为泊松比，

S22,假设层裂向前扩展Δa，根据虚拟扩展裂纹与网格单元之间的位置关系，对虚拟扩展裂纹Δa进行分段积分，若扩展裂纹Δa位于单个单元内，则无需进行分段；若扩展裂纹Δa经过多个单元，则利用所经过单元对扩展裂纹Δa进行自然分段；同时将分段积分点沿着裂纹面方向平移Δa，获得对应的张开位移分量；

S23，在S21和S22的基础上，对虚拟扩展裂纹各分段分别进行高斯积分并求和，获得层裂裂纹尖端的能量释放率，其数值积分由如下三式给出：

式中N_seg为虚拟扩展裂纹分段数目，N_gau为高斯积分点数，x_mn为第m分段第n高斯积分点坐标，可由裂纹分段两端坐标及高斯积分点位置确定，w_n为第n高斯积分点的权重。

4.根据权利要求3所述的一种复合材料层裂二维计算方法，其特征在于，所述步骤S3包括以下步骤：

S31,根据算得的能量释放率，判断是否满足层裂扩展条件；

式中G_Ic、G_IIc、G_IIIc分别为I型断裂韧度、II型断裂韧度、III型断裂韧度；

S32，若满足上式，则层裂向前扩展Δa长度，更新裂纹几何尺寸，回到S12步，继续计算；若不满足上式，则增加载荷，继续计算。