CN108416084A - 考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法 - Google Patents

Abstract

本发明公开了一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，属于复合材料力学性能有限元分析领域；以基于后退欧拉积分法则及牛顿迭代求解复合材料弹塑性损伤的模型为基础，通过在初始求解试用应力以及塑性屈服迭代中求解弹性应力时均考虑前一增量步中损伤状态量对当前应力的软化影响，从而在弹塑性求解过程中考虑损伤退化作用，更真实地反映树脂基复合材料在复杂加载条件下塑性行为与损伤之间的相互作用关系。本发明可有效地描述树脂基复合材料弹塑性与损伤的耦合作用关系，可应用于有限元软件中复合材料弹塑性损伤本构模型的开发、新型树脂基复合材料力学性能研究以及相关的复合材料工程结构精细化建模分析等科研及工程技术领域。

Description

考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法

技术领域

本发明涉及复合材料力学性能有限元分析技术，尤其涉及一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法。

背景技术

目前，复合材料已被广泛应用于军事装备、航空航天、民用工业等诸多领域，而复合材料在受载时表现出一定的非线性力学行为，并伴随着材料力学性能、结构承载能力和使用寿命的退化。因此，揭示复合材料在受载时的非线性力学行为并准确地预测分析该力学特性，对复合材料结构承载性能分析、结构力学性能优化设计、疲劳寿命预测具有重要的科研价值及工程意义。

关于复合材料受载时的非线性力学行为，可归结为非线性的损伤退化行为以及复合材料基体塑性变形的累积。这两种非线性行为在加载作用下存在着相互耦合作用的关系。目前，已有关于复合材料的弹塑性损伤有限元算法被提出，并被应用于科研及工程技术领域中。而关于复合材料弹塑性损伤行为预测的有限元算法中，主要分别孤立地考虑弹塑性阶段及损伤退化阶段，忽略了塑性屈服阶段中损伤与塑性之间的耦合作用关系，其主要表现在算法实现过程中采用前一增量步的状态量求解当前增量步的塑性屈服状态时，并未考虑到前一增量步中材料的损伤状态对当前判断屈服状态以及进入塑性屈服阶段的影响，导致最终无法描述树脂基复合材料在加载过程中非线性损伤与塑性之间的相互作用关系。

并且，现有的商用大型有限元软件中，未见存在能够直接提供给用户对复合材料弹塑性连续损伤行为进行分析的模型，均需要用户自主编写子程序予以实现，而合适的可用于有限元中预测复合材料弹塑性连续损伤行为的算法是紧缺的，进一步增加了工程技术人员对复合材料结构的弹塑性损伤特性进行有限元分析的难度。

发明内容

本发明的目的在于提供一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法。

实现本发明目的的技术方案为：一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，包括以下步骤：

步骤1，开始当前增量步，读取前一时刻收敛状态量及当前增量步中应变增量；

步骤2，假定当前步开始时无塑性应变增量，则当前步试用弹性应变由前一时刻收敛的弹性应变量与当前增量步应变增量线性叠加获取，通过弹性本构关系根据当前步的弹性试用应变获取当前的试用弹性应力，考虑前一增量步损伤状态对当前试用弹性应力的软化影响，得出当前试用弹性应力状态量；

步骤3，判断材料屈服状态；若未进入屈服阶段，计算当前增量步状态量；若进入屈服阶段，对含损伤的塑性屈服阶段通过牛顿迭代法求解，迭代过程中每一迭代步的试用弹性应力计算均考虑前一增量步损伤状态的影响；

步骤4，根据损伤力学模型，求解当前增量步中的损伤指标，更新损伤状态量；

步骤5，求解并输出含塑性及损伤退化作用的状态量，予以保存，开始下一增量步。

与现有的复合材料弹塑性损伤有限元算法相比，本发明的有益效果为：本发明通过在求解试用应力以及塑性屈服阶段迭代求解中考虑前一增量步损伤状态量的退化作用，能够描述复合材料受载出现非线性行为中损伤与塑性同时出现时的相互作用关系，形成一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，从而能有效地预测复合材料的弹塑性损伤耦合行为。

附图说明

图1为本发明考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法流程图。

具体实施方式

结合图1，一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，包括以下步骤：

步骤1，开始当前增量步，读取前一时刻收敛状态量及当前增量步中应变增量；

步骤2，假定当前步开始时无塑性应变增量，则当前步试用弹性应变由前一时刻收敛的弹性应变量与当前增量步应变增量线性叠加获取，通过弹性本构关系根据当前步的弹性试用应变获取当前的试用弹性应力，考虑前一增量步损伤状态对当前试用弹性应力的软化影响，得出当前试用弹性应力状态量；

步骤3，判断材料屈服状态；若未进入屈服阶段，计算当前增量步状态量；若进入屈服阶段，对含损伤的塑性屈服阶段通过牛顿迭代法求解，迭代过程中每一迭代步的试用弹性应力计算均考虑前一增量步损伤状态的影响；

步骤4，根据损伤力学模型，求解当前增量步中的损伤指标，更新损伤状态量；

步骤5，求解并输出含塑性及损伤退化作用的状态量，予以保存，开始下一增量步。

步骤3在塑性求解迭代过程中计算弹性应力时考虑前一增量步的损伤状态量的退化作用，考虑弹塑性与损伤同时发生的耦合作用。

进一步的，步骤1中，读取前一增量步即第n步结束时材料单元网格积分点上的有效应力应变ε_n、塑性应变等效塑性应变值损伤状态量D_n以及当前增量步即第n+1增量步的应变增量△ε_n+1。

进一步的，步骤2具体为：

假定当前无塑性增量以得出上标含trail表示试用状态量，上标含p表示塑性变量，上标含e表示弹性变量，为第n+1增量步的试用塑性增量，为第n+1增量步的试用弹性增量，为第n+1增量步的试用等效塑性应变，采用弹性本构模型并考虑前一增量步损伤状态量的退化作用得出当前试用应力状态量其中C为材料弹性刚度矩阵。

进一步的，步骤3具体为：

选取屈服准则，通过给定的关于复合材料的屈服函数判断是否进入塑性屈服状态，若屈服函数成立，则无塑性变形，若不成立则产生塑性屈服；

未进入塑性阶段表示无塑性增量产生，则当前的塑性应变和等效塑性应变为当前应力和总应变为

进入塑性屈服阶段后，关于塑性加载准则对当前的塑性增量进行迭代求解，根据求解塑性乘子△λ_n+1，并且△λ_n+1>0，该求解过程基于后向欧拉积分法和牛顿迭代法，求解过程如下：

(1)迭代初始步k＝0中变量初始化 为第n+1增量步的试用损伤状态量，分别为第n+1增量步初始迭代步的塑性乘子、塑性应变量、等效塑性应变量和弹性应力；

(2)牛顿迭代法开始迭代，第k+1步获取则当前增量步n+1步中的第k+1迭代步的塑性应变为等效塑性应变为并考虑前一增量步损伤状态量对应力的影响，则应力为

(3)判断当前迭代步屈服函数是否满足数值偏差e为数值偏差，若满足则塑性求解收敛，输出当前迭代步所得出的塑性应变、等效塑性应变及考虑前一步损伤影响的应力；若不满足则返回步骤(2)中继续进行迭代求解，直至结果收敛；

(4)输出迭代求解的收敛值至此塑性阶段求解结束。

进一步的，步骤4具体为：

根据塑性阶段求解结果，选择损伤判断准则，获取当前材料单元积分点处的损伤指标，并选择损伤演化准则，求解得出当前的损伤状态量，更新损伤变量矩阵D_n+1以及损伤变量矩阵的增量△D_n+1。

损伤判断准则包括Hashin损伤判断准则、Puck损伤判断准则，损伤演化准则包括基于应变的线性损伤变量演化准则、基于断裂能的非线性损伤变量演化准则、基于断裂密度的非线性损伤变量演化准则。

进一步的，步骤5具体为：

基于当前增量步中的损伤变量矩阵及其增量求解损伤后的应力状态将当前增量步结束时的有效应力应变ε_n+1、塑性应变等效塑性应变值损伤状态量D_n+1予以保存，用以传递到下一增量步的计算中。

下面结合具体实施例对本发明进行详细说明。

实施例

结合图1，一种考虑了纤维增强树脂基复合材料弹塑性与损伤耦合作用的弹塑性损伤有限元子程序算法，包括以下步骤：

步骤1，开始当前增量步，读取前一时刻收敛状态量及当前增量步中应变增量；具体为：

读取前一增量步即t_n时刻材料单元积分点上的有效应力应变ε_n、塑性应变等效塑性应变值损伤状态量D_n以及当前增量步即第n+1增量步t_n～t_n+1时刻的应变增量

有效应力应变塑性应变损伤状态量D_n为

其中，d₁,d₂分别是表征纤维损伤和基体损伤的损伤变量，d₁₂,d₂₃,d₁₃为表征剪切变形的损伤变量。

步骤2，假定初始试用状态量，考虑损伤状态的软化影响；具体为：

假定当前无塑性增量，得出当前步开始时的试用状态量，则

采用弹性本构模型并考虑前一增量步损伤状态量的退化作用得出当前试用应力状态量

所述的试用状态变量的求解公式中，I为6×6的单位矩阵，C为材料的刚度矩阵，表示如下:

所述的材料刚度矩阵表达式中E₁,E₂,E₃分别是沿着纤维方向、面内横向和面外横向的杨氏模量，G₁₂,G₂₃,G₃₁是剪切模量，v_ij(i≠j且i,j＝1,2,3)是泊松比，Ω表达式如下

Ω＝1-v₁₂v₂₁-v₂₃v₃₂-v₁₃v₃₁-2v₁₂v₃₁v₂₃

步骤3，判断材料屈服状态；

若未进入屈服阶段，计算当前步状态量；

若进入屈服阶段，对含损伤的塑性屈服阶段通过牛顿迭代法求解；具体为：

选取复合材料常用的屈服函数，通过屈服函数判断是否进入塑性屈服状态式中a₆₆,β,α均为材料常数，对应于不同纤维复合材料可从相关论文中获取，若关于当前步的试用状态量屈服函数成立，则无塑性变形，若不成立则产生塑性屈服；

若未进入塑性阶段则表示无塑性增量产生，则总的应变ε_n+1＝ε_n+△ε_n+1，当前的塑性应变等效塑性应变值弹性应变当前应力和总应变为

进入塑性屈服阶段后，关于塑性加载准则对当前的塑性增量进行迭代求解，即求解塑性乘子△λ_n+1>0满足式中e为数值小量可取值为10^-6，该求解过程主要基于后向欧拉积分法和牛顿迭代法，求解过程如下：

a.迭代初始步k＝0中变量初始化

b.牛顿迭代法开始迭代，第k+1步获取则当前迭代步的塑性应变为等效塑性应变为并考虑前一增量步损伤状态量对应力的影响，则应力为

c.判断当前迭代步屈服函数是否满足数值偏差e为数值偏差，可取值为1E-6，若满足则塑性求解收敛，输出当前迭代步所得出的塑性应变、等效塑性应变及考虑前一步损伤影响的应力；若不满足则返回步骤b中继续进行迭代求解，直至结果收敛；

d.输出迭代求解的收敛值塑性应变等效塑性应变

总应变ε_n+1＝ε_n+△ε_n+1，弹性应变弹性应力至此塑性阶段求解结束。

步骤4，根据损伤力学模型，求解当前增量步中的损伤指标，更新损伤状态量；具体为：

选择基于应变Hashin损伤判断准则，获取当前步中纤维的拉伸和压缩损伤指标f_FT,f_FC，以及基体的拉伸和压缩损伤指标f_MT,f_MC，该损伤判断准则具体数学形式参见Feng和Aymerich于2014年发表在Composite Structure上的文章“Finite element modeling ofdamage induced by low-velocity impact on composite laminates”中的3.1.1节和3.1.2节内容；

选择基于应变的线性损伤变量演化准则获取纤维拉伸损伤变量和纤维压缩损伤变量基体拉伸损伤变量和压缩损伤变量该基于应变的线性损伤变量演化准则具体数学形式同样可参见Feng和Aymerich于2014年发表在Composite Structure上的文章“Finite element modeling of damage induced by low-velocity impact oncomposite laminates”中的3.1.1节和3.1.2节内容，则纤维损伤变量基体损伤变量剪切损伤变量 则根据所得出的当前步的损伤变量更新损伤变量矩阵D_n+1，以及损伤变量矩阵的增量△D_n+1＝D_n+1-D_n。

步骤5，求解并输出含塑性及损伤退化作用的状态量，予以保存，开始下一增量步。具体为：

基于当前增量步中的损伤变量矩阵及其增量求解损伤后的应力状态将当前增量步结束时的有效应力应变ε_n+1、塑性应变等效塑性应变值损伤状态量D_n+1予以保存，用以传递到下一增量步的计算中，此时当前步的计算过程结束，下一增量步的计算按步骤1至步骤5进行计算，如此往复即可。

Claims (6)

1.一种考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，包括以下步骤：

步骤1，开始当前增量步，读取前一时刻收敛状态量及当前增量步中应变增量；

步骤2，假定当前步开始时无塑性应变增量，则当前步试用弹性应变由前一时刻收敛的弹性应变量与当前增量步应变增量线性叠加获取，通过弹性本构关系根据当前步的弹性试用应变获取当前的试用弹性应力，考虑前一增量步损伤状态对当前试用弹性应力的软化影响，得出当前试用弹性应力状态量；

步骤3，判断材料屈服状态；若未进入屈服阶段，计算当前增量步状态量；若进入屈服阶段，对含损伤的塑性屈服阶段通过牛顿迭代法求解，迭代过程中每一迭代步的试用弹性应力计算均考虑前一增量步损伤状态的影响；

步骤4，根据损伤力学模型，求解当前增量步中的损伤指标，更新损伤状态量；

步骤5，求解并输出含塑性及损伤退化作用的状态量，予以保存，开始下一增量步。

2.根据权利要求1所述的考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，步骤1中，读取前一增量步即第n步结束时材料单元网格积分点上的有效应力应变ε_n、塑性应变等效塑性应变值损伤状态量D_n以及当前增量步即第n+1增量步的应变增量△ε_n+1。

3.根据权利要求2所述的考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，步骤2具体为：

假定当前无塑性增量以得出上标含trail表示试用状态量，上标含p表示塑性变量，上标含e表示弹性变量，为第n+1增量步的试用塑性增量，为第n+1增量步的试用弹性增量，为第n+1增量步的试用等效塑性应变；采用弹性本构模型并考虑前一增量步损伤状态量的退化作用得出当前试用应力状态量其中C为材料弹性刚度矩阵。

4.根据权利要求3所述的考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，步骤3具体为：

选取屈服准则，通过给定的关于复合材料的屈服函数判断是否进入塑性屈服状态，若屈服函数成立，则无塑性变形，若不成立则产生塑性屈服；

未进入塑性阶段表示无塑性增量产生，则当前的塑性应变和等效塑性应变为当前应力和总应变为

进入塑性屈服阶段后，关于塑性加载准则对当前的塑性增量进行迭代求解，根据求解塑性乘子△λ_n+1，并且△λ_n+1>0，该求解过程基于后向欧拉积分法和牛顿迭代法，求解过程如下：

(1)迭代初始步k＝0中变量初始化 为第n+1增量步的试用损伤状态量，分别为第n+1增量步初始迭代步的塑性乘子、塑性应变量、等效塑性应变量和弹性应力；

(2)牛顿迭代法开始迭代，第k+1步获取则当前增量步n+1步中的第k+1迭代步的塑性应变为等效塑性应变为并考虑前一增量步损伤状态量对应力的影响，则应力为

(3)判断当前迭代步屈服函数是否满足数值偏差e为数值偏差，若满足则塑性求解收敛，输出当前迭代步所得出的塑性应变、等效塑性应变及考虑前一步损伤影响的应力；若不满足则返回步骤(2)中继续进行迭代求解，直至结果收敛；

(4)输出迭代求解的收敛值至此塑性阶段求解结束。

5.根据权利要求4所述的考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，步骤4具体为：

根据塑性阶段求解结果，选择损伤判断准则，获取当前材料单元积分点处的损伤指标，并选择损伤演化准则，求解得出当前的损伤状态量，更新损伤变量矩阵D_n+1以及损伤变量矩阵的增量△D_n+1。

6.根据权利要求5所述的考虑复合材料弹塑性与损伤耦合的弹塑性损伤有限元算法，其特征在于，步骤5具体为：

基于当前增量步中的损伤变量矩阵及其增量求解损伤后的应力状态将当前增量步结束时的有效应力应变ε_n+1、塑性应变等效塑性应变值损伤状态量D_n+1予以保存，用以传递到下一增量步的计算中。