CN106934133B - 基于弹塑性分解的非线性有限元刚度矩阵更新方法 - Google Patents

Abstract

基于弹塑性分解的非线性有限元刚度矩阵更新方法，属于结构分析领域。该方法通过对非线性材料的应变进行分解，使用插值形式描述单元分解出的非线性应变场，并在数值计算过程中仅考虑进入非线性变形状态的局部区域，实现了整体弹性结构与局部非线性结构的状态分离，从而将传统方法中大规模刚度矩阵的更新过程转变为对初始弹性刚度矩阵的低秩摄动过程。本发明具有如下特点：1)非线性变形状态分离，本发明在非线性求解过程中使整体刚度矩阵保持弹性状态不变，而使用较小规模的矩阵代表进入非线性状态的局部结构，迭代时仅需对与局部非线性结构有关的小规模矩阵进行更新；2)适用性广，本发明不局限于特定的单元类型和材料属性，具有较强的通用性。

Description

基于弹塑性分解的非线性有限元刚度矩阵更新方法

技术领域

本发明属于结构分析领域，具体涉及基于弹塑性分解的非线性有限元刚度矩阵更新方法。

背景技术

在结构分析领域，有限元方法作为一种应用广泛的数值计算方法，由于其通用性和有效性，目前已成为计算机辅助设计及数值仿真的重要实现手段。对于极端环境荷载作用下的结构响应计算，需考虑材料的非线性属性，由于此时求解方程中刚度矩阵需根据当前的结构非线性状态实时变化，因此通常采用步进式的增量格式及离散的线性化方法进行求解，这一过程中迭代计算不可避免，且需不断进行整体刚度矩阵的更新和分解运算。当计算模型具有较大规模或较为精细的单元划分时，大规模刚度矩阵的实时更新和分解过程将极为耗时，成为制约非线性有限元技术应用的主要因素。

发明内容

为了克服上述现有方法的不足，本发明通过对材料应变进行分解，提出了一种新的用于材料非线性问题的刚度矩阵更新方法。本发明可将进入非线性变形状态的局部结构区域与整体弹性结构相互分离，求解时仅需对代表局部非线性变形的小规模矩阵进行更新，从而将传统方法中由局部非线性变形引起的大规模刚度矩阵更新过程转变为对初始弹性刚度矩阵的低秩摄动过程。本发明具有较为广泛的适用性，可用于众多领域的非线性结构分析中。

本发明的技术方案：

一种基于弹塑性分解的非线性有限元刚度矩阵更新方法，步骤如下：

(1)应变分解

对于具有非线性本构关系的材料，在任意状态下其应力和应变分别用σ和ε表示，将应变ε分解为线弹性应变ε′与非线性应变ε"之和的形式，即：

ε＝ε′+ε″ (1)

式中，线弹性应变ε′定义为按照材料的初始弹性属性加载至应力σ时对应的应变，非线性应变ε"定义为总应变ε与线弹性应变ε′的差，应力、线弹性应变和非线性应变之间符合如下关系：

σ＝D_eε′＝D_e(ε-ε″) (2)

式中，D_e代表材料的初始弹性本构矩阵，对于单轴材料，其退化为材料初始弹性模量；

(2)使用插值方法建立单元非线性变形机制

对于待分析结构，建立对应的有限元数值模型，并在单元内设置若干非线性变形插值结点，从而单元域内任意一点的非线性应变可通过非线性应变场的插值表达式获得，此时，单元非线性应变场表示为如下插值形式：

ε″＝Cε″_p (3)

式中，C为非线性插值函数矩阵；向量ε″_p集中了单元内所有非线性变形插值点的非线性应变；

(3)在任意增量求解步中，仅考虑模型结构中进入非线性变形状态的局部区域，并使用增量格式，建立分块形式的结构控制方程：

式中，K_e为模型结构的初始弹性刚度矩阵；向量ΔF、ΔX分别代表结构的荷载增量和位移增量；向量ΔΕ″_pr为模型中进入非线性变形状态的局部子结构非线性应变增量；K″_pr为进入非线性变形状态的局部子结构刚度矩阵；K′_r为整体弹性结构与局部非线性结构之间的交互矩阵；

(4)建立模型结构的非线性求解方程，即：

步骤1中，应变分解过程可适用于任意形式的非线性材料本构关系，当卸载刚度与初始弹性刚度相同时，步骤1中定义的线弹性应变与非线性应变将分别等于弹性应变与塑性应变。

步骤2中，由插值方法建立的单元非线性变形机制不局限于特定的单元类型，但非线性变形插值点的个数和位置影响单元计算精度。

步骤3和步骤4中，需首先判断结构中所有单元非线性变形插值点的非线性状态，并在构造矩阵K′_r和K″_pr的过程中使其仅包含处于非线性变形状态的单元和插值点信息。

步骤3和步骤4中，矩阵K_e代表整体弹性结构，在分析过程中始终不变；向量ΔΕ″_pr由所有位于局部非线性变形区域内的非线性插值点的非线性应变集成得到，其中不包括处于线弹性状态的插值点信息；K″_pr代表进入非线性变形状态的局部子结构，包含位于其中的非线性插值点材料非线性本构信息，需在分析时根据材料非线性状态实时更新；矩阵K′_r与数值模型中产生非线性变形局部区域的位置有关，而与材料的非线性程度无关，其中的每个元素仅需计算一次，无需在非线性求解过程中重复计算。

步骤4中，求解方程式(5)通过对式(4)进行凝聚获得，等号左边括号中的矩阵表达式为初始弹性刚度矩阵的低秩摄动，与整体结构的切向刚度矩阵等效。

本发明的有益效果：

1.非线性变形状态的分离。对于非线性仅产生于结构局部的有限元计算问题，本发明在控制方程中隔离出了代表非线性变形状态的小规模矩阵，进而在非线性求解时仅需对这一小规模矩阵进行更新，避免了对大规模整体刚度矩阵的更新。

2.广泛的适用范围。本发明所提方法在结构非线性分析时不局限于特定的材料本构及单元类型，可应用到众多学科及研究领域，例如土木、水利、机械、力学等，因而具有广泛的适用性。

附图说明

图1为实施例一的应变分解示意图。

图2为实施例二的有限元模型。

图3为实施例二的单元模型。

图4为实施例二的单元内任一点的应力与应变。

具体实施方式

本发明提供了一种基于应变分解思想的刚度矩阵更新方法，以下结合附图实例和技术方案，进一步说明本发明的具体实施方式。附图和实施例仅为说明本发明的实施方式，不构成对本发明的任何限制。

实施例一：材料应变分解实施方式

以图1所示单轴非线性材料的应力-应变关系为例说明本发明中材料应变分解的实施方式。图中C点代表当前的应力应变状态，其应力和应变分别为σ和ε，E_e代表材料初始弹性模量，延长直线OA交应力σ于B点，定义B点横坐标为线弹性应变ε′，定义ε与线弹性应变ε′的差为非线性应变ε"，即：

ε″＝ε-ε′

实施例二：非线性问题求解

以图2所示某个由三结点平面应力单元组成的有限元模型为例，说明本发明所提出方法在非线性有限元计算中的实施方式。对于平面应力单元，单元内任意一点的应变可由3个分量表示，即ε＝(ε_xx ε_yy γ_xy)^T，其中ε_xx、ε_yy为正应变，γ_xy为剪应变，相应的应力分量为σ＝(σ_xx σ_yy τ_xy)^T，其中σ_xx、σ_yy为正应力，τ_xy为剪应力。此时材料的应变分解表达式为：

式中，ε′＝(ε′_xx ε′_yy γ′_xy)^T代表材料线弹性应变；ε″＝(ε″_xx ε″_yy γ″_xy)^T代表材料非线性应变。

在每个单元中设置一个非线性变形插值点，如图3所示，则单元非线性应变场可表示为如下插值形式：

ε″＝Cε″_p (2)

式中，C为非线性插值函数矩阵；ε″_p为非线性变形插值点处的非线性应变。本例中，由于仅有一个插值点，非线性应变在单元内按常数分布，C为单位矩阵。

非线性求解采用Newton-Raphson迭代方法，在非线性求解前，首先计算各个单元弹性刚度矩阵k_e及交互矩阵k′，对于任意单元，相应的k_e和k′的计算公式为：

式中，B为单元应变矩阵，代表的单元应变分布和结点位移的关系；D_e为材料初始弹性本构矩阵。其次，对单元矩阵k_e进行集成，获得结构整体弹性刚度矩阵K_e。

在非线性求解过程中，任一增量迭代步中的求解方程更新过程可按如下步骤进行：

步骤一：进行单元状态的确定，用于计算在当前位移状态下各单元的实际非线性状态，包括计算非线性插值点处的切线材料本构矩阵D_t、判断单元是否进入非线性变形状态。若进入非线性变形状态，即D_t≠D_e，则计算单元的非线性刚度矩阵k″_p，本例中，采用高斯积分算法并令非线性变形插值点与高斯积分点重合，则k″_p可表示为

k″_p＝AtD_e(D_e-D_t)^-1D_e (5)

其中，A为单元面积；t为单元厚度。上式中D_e与D_t均为非线性变形插值点处的材料属性。若塑性插值点未进入塑性变形状态，即D_t＝D_e，则不计算矩阵k″_p。

步骤二：根据由步骤一所得的各单元的非线性状态，更新矩阵K″_pr和K′_r。其中矩阵K″_pr由所有进入非线性变形状态的单元矩阵k″_p集成得到，矩阵K′_r由所有进入非线性变形状态的单元矩阵k′集成得到。

步骤三：更新非线性求解方程并完成求解：

Claims (1)

1.一种基于弹塑性分解的非线性有限元刚度矩阵更新方法，其特征在于，步骤如下：

(1)应变分解

对于具有非线性本构关系的材料，在任意状态下其应力和应变分别用σ和ε表示，将应变ε分解为线弹性应变ε′与非线性应变ε"之和的形式，即：

ε＝ε′+ε″ (1)

式中，线弹性应变ε′定义为按照材料的初始弹性属性加载至应力σ时对应的应变，非线性应变ε"定义为总应变ε与线弹性应变ε′的差，应力、线弹性应变和非线性应变之间符合如下关系：

σ＝D_eε′＝D_e(ε-ε″) (2)

式中，D_e代表材料的初始弹性本构矩阵，对于单轴材料，其退化为材料初始弹性模量；

(2)使用插值方法建立单元非线性变形机制

对于待分析结构，建立对应的有限元数值模型，并在单元内设置若干非线性变形插值结点，从而单元域内任意一点的非线性应变通过非线性应变场的插值表达式获得，此时，单元非线性应变场表示为如下插值形式：

ε″＝Cε″_p (3)

式中，C为非线性插值函数矩阵；向量ε″_p集中了单元内所有非线性变形插值点的非线性应变；

非线性求解采用Newton-Raphson迭代方法，在非线性求解前，首先计算各个单元弹性刚度矩阵k_e及交互矩阵k′，对于任意单元，相应的k_e和k′的计算公式为：

式(4)中，B为单元应变矩阵，代表的单元应变分布和结点位移的关系；D_e为材料初始弹性本构矩阵；其次，对单元矩阵k_e进行集成，获得结构整体弹性刚度矩阵K_e；

在非线性求解过程中，任意增量迭代步中的求解方程更新过程按如下步骤进行：

步骤一：进行单元状态的确定，用于计算在当前位移状态下各单元的实际非线性状态，包括计算非线性插值点处的切线材料本构矩阵D_t、判断单元是否进入非线性变形状态；若进入非线性变形状态，即D_t≠D_e，则计算单元的非线性刚度矩阵k″_p，如公式(5)所示：

k″_p＝AtD_e(D_e-D_t)^-1D_e (5)

其中，A为单元面积；t为单元厚度；式(5)中D_e与D_t均为非线性变形插值点处的材料属性；若塑性插值点未进入塑性变形状态，即D_t＝D_e，则不计算矩阵k″_p；

步骤二：根据步骤一所得的各单元的非线性状态，更新矩阵K″_pr和K′_r；其中矩阵K″_pr由所有进入非线性变形状态的单元矩阵k″_p集成得到，矩阵K′_r由所有进入非线性变形状态的单元矩阵k′集成得到；

步骤三：更新非线性求解方程并完成求解：

(3)在任意增量求解步中，仅考虑模型结构中进入非线性变形状态的局部区域，并使用增量格式，建立分块形式的结构控制方程：

式中，K_e为模型结构的初始弹性刚度矩阵；向量ΔF、ΔX分别代表结构的荷载增量和位移增量；向量ΔΕ″_pr为模型中进入非线性变形状态的局部子结构非线性应变增量；K″_pr为进入非线性变形状态的局部子结构刚度矩阵；K′_r为整体弹性结构与局部非线性结构之间的交互矩阵；

(4)建立模型结构的非线性求解方程，即：