一种用于求解局部材料非线性问题的高效计算方法

摘要

本发明公开了一种高效求解局部材料非线性问题的计算方法，属于结构分析技术领域。该方法将非线性变形转变为作用于弹性结构上的虚拟荷载，从而刚度矩阵在非线性分析中保持不变，仅需进行一次分解；在方程构建过程中仅考虑进入非线性变形状态的局部结构信息，则迭代求解过程中仅需对一个代表局部非线性的小规模矩阵进行分解运算，避免了传统方法中对大规模整体刚度矩阵的实时更新和分解，有效提升了局部非线性问题的计算效率。此外，本发明不局限于特定的材料本构关系和数值模型，可应用到众多涉及结构非线性计算的学科及研究领域，具有较强的通用性和较广泛的适用性。

权利要求

1.一种用于求解局部材料非线性问题的高效计算方法，其特征在于，步骤如下：(1)建立用于结构分析的数值模型，并在结构中预设若干非线性结点；

(2)计算整体结构的初始弹性刚度矩阵Ke并对其进行三角分解，设结构总位移自由度数是n，则矩阵Ke的规模为n×n；

(3)对于任意非线性增量计算步，在增量荷载ΔF作用下，仅考虑进入非线性变形状态局部区域内的非线性结点并建立相应增量形式的非线性变形向量：ΔΕ″p＝[Δε″1 Δε″2 ... Δε″m]T   (1)式中，m为该增量步中进入非线性变形状态局部区域内的非线性结点个数；Δε″i代表其中第i个非线性结点的增量非线性变形，其中1≤i≤m；设每个非线性结点的非线性变形由p个分量组成，则ΔΕ″p为一个pm阶向量，其中每个元素均代表一个非线性自由度；

(4)求解增量荷载ΔF作用下的结构位移增量，即：

式中，ΔX代表位移增量；ΔFp＝K′rΔΕ″p，为代表结构非线性的附加虚拟荷载，其中矩阵K′r为非线性变形向量ΔΕ″p与附加虚拟荷载ΔFp与之间关系的刚度矩阵，其阶数为n×pm；给定荷载增量ΔF作用下的非线性变形向量ΔΕ″p可使用下式求解：式中， 是K′r的转置，为pm×n阶矩阵，代表结点位移与非线性结点处的弹性内力之间关系，弹性内力指具有初始弹性属性的结构发生位移改变时非线性结点处产生的应力；

K″pr代表非线性结点处的非线性变形与弹性内力之间的关系，为pm×pm阶矩阵，需根据非线性结点处的材料本构关系确定。

说明书

一种用于求解局部材料非线性问题的高效计算方法

技术领域

[0001] 本发明属于结构分析技术领域，更具体的说，涉及一种用于求解局部材料非线性问题的高效计算方法。

背景技术

[0002] 利用数值仿真技术计算结构在极端环境荷载作用下的非线性性能已在众多研究领域得到了广泛的应用。由于结构荷载-变形关系的非线性属性，通常需采用步进式的增量格式并结合适当的迭代算法进行求解，这一过程中，结构的整体刚度矩阵需在计算过程中不断进行更新和分解运算，随着计算模型规模和精细化程度的增加，刚度矩阵的规模亦迅速增加，此时，大规模刚度矩阵的实时更新和分解成为制约非线性问题计算效率的关键。对于材料非线性问题，通常仅结构局部发生非线性变形，而其余大部分将始终处于初始弹性变形状态，传统方法计算过程中并不考虑这一特点，因而不能实现弹性和非线性部分的差别化对待，为此，如何充分利用非线性变形的局部特征实现高效数值求解成为这类问题的研究重点。

发明内容

[0003] 为了克服上述现有方法的不足，本发明利用Woodbury公式，通过考虑材料非线性问题的局部特征，提出了一种用于结构局部非线性问题的高效计算方法，计算过程仅需对一个代表局部非线性变形的小规模矩阵进行更新和分解运算，避免了传统方法中对大规模整体刚度矩阵进行的相应的运算。发明具有较为广泛的适用性，可用于众多领域的非线性结构分析中。

[0004] 本发明的技术方案：

[0005] 一种用于求解局部材料非线性问题的高效计算方法，步骤如下：

[0006] (1)建立用于结构分析的数值模型，并在结构中预设若干非线性结点；

[0007] (2)计算整体结构的初始弹性刚度矩阵Ke并对其进行三角分解，设结构总位移自由度数是n，则矩阵Ke的规模为n×n；

[0008] (3)对于任意非线性增量计算步，在增量荷载ΔF作用下，仅考虑进入非线性变形状态局部区域内的非线性结点并建立相应增量形式的非线性变形向量：

[0009] ΔΕ″p＝[Δε″1 Δε″2 ... Δε″m]T        (1)

[0010] 式中，m为该增量步中进入非线性变形状态局部区域内的非线性结点个数；Δε″i代表其中第i个非线性结点的增量非线性变形，其中1≤i≤m；设每个非线性结点的非线性变形由p个分量组成，则ΔΕ″p为一个pm阶向量，其中每个元素均代表一个非线性自由度；

[0011] (4)求解增量荷载ΔF作用下的结构位移增量，即：

[0012]

[0013] 式中，ΔX代表位移增量；ΔFp＝K′rΔΕ″p，为代表结构非线性的附加虚拟荷载，其中矩阵K′r为非线性变形向量ΔΕ″p与附加虚拟荷载ΔFp与之间关系的刚度矩阵，其阶数为n×pm。给定荷载增量ΔF作用下的非线性变形向量ΔΕ″p可使用下式求解：

[0014]

[0015] 式中， 是K′r的转置，为pm×n阶矩阵，代表结点位移与非线性结点处的弹性内力之间关系，弹性内力指具有初始弹性属性的结构发生位移改变时非线性结点处产生的应力；K″pr代表非线性结点处的非线性变形与弹性内力之间的关系，为pm×pm阶矩阵，需根据非线性结点处的材料本构关系确定。

[0016] 步骤1中，结构中的非线性结点位置和数量需根据具体问题的精度要求确定，非线性变形场可通过插值方法获得。

[0017] 步骤3中，任意状态下的非线性变形等于按照材料初始弹性刚度卸载后的残余变形；向量ΔΕ″p中包含当前增量步下结构中所有进入非线性变形状态的非线性结点变形，但不包括处于初始弹性变形状态的非线性结点变形。

[0018] 步骤3中，矩阵K′r及 与数值模型中非线性出现部位有关，而与材料的非线性程度无关，矩阵中的每个元素仅需计算一次，无需在非线性迭代过程中重复计算。

[0019] 步骤3和步骤4中，刚度矩阵Ke的逆矩阵与向量或矩阵的乘法运算可通过对Ke的三角分解和回代计算完成，由于该矩阵在分析过程中保持不变，其三角分解仅需在分析前进行一次。

[0020] 步骤4中，任意非线性增量计算步的位移增量计算方程与Woodbury公式等价，可写为如下形式的数学表达式：

[0021]

[0022] 本发明中，每个计算步对矩阵 分解所需的计算量成为使用求解时的主要计算消耗，该矩阵的规模为pm×pm(与非线性自由度数目相当)，代表进入非线性变形状态的局部区域。当结构仅少部分进入非线性时，有pm<分解所需的计算量，从而可显著提升非线性问题的计算效率。

[0023] 本发明的有益效果：

[0024] 1.高效的非线性计算。对于非线性仅产生于结构局部的材料非线性问题，本发明在求解时仅需对一代表局部非线性的小规模进行分解运算，避免了现有方法中对规模较大的整个刚度矩阵进行的实时更新和分解运算，可在保证必要计算精度的前提下有效提升非线性问题的求解效率。

[0025] 2.广泛的适用范围。本发明不局限于特定的材料本构关系和数值模型，可应用到众多涉及结构非线性计算的学科及研究领域，因而具有广泛的适用性。

附图说明

[0026] 图1为实施例的分析模型图。

[0027] 图2为本发明的计算流程图。具体实施方案

[0028] 本发明提供了一种用于求解局部材料非线性问题的高效数值计算方法，以下结合附图实例和技术方案，进一步说明本发明的具体实施方式。附图和实施例仅为说明本发明的实施方式，不构成对本发明的任何限制。

[0029] 实施例：

[0030] 1.以图1所示某个由二结点桁架单元组成的数值计算为例，说明本发明所提出方法在非线性计算中的实施方式。对于平面应力单元，考虑其截面轴向应变ε与轴力N的非线性关系，在单元内设置一个非线性结点并令其位于单元中间，如图1所示，其截面变形的分解表达式为：

[0031] ε＝ε′+ε″        (1)

[0032] 式中，ε′代表截面线弹性应变，即使用截面初始弹性刚度加载至轴力N时的截面变形；ε″代表截面非线性应变。

[0033] 2.合成整体结构的初始弹性刚度矩阵Ke并对其进行三角分解，即：

[0034] Ke＝LDLT         (2)

[0035] 式中，L为单位下三角矩阵，D为对角矩阵。

[0036] 3.结合图2给出的本发明在非线性求解过程中的迭代计算流程，任一迭代增量步的计算可按如下步骤进行：

[0037] 步骤一：进行单元状态的确定，用于计算在当前位移条件下各单元中非线性结点处的材料状态，包括截面内力N及材料切线模量Et，并进一步判断单元是否进入非线性变形状态，若某单元非线性结点处有Et≠Ee，则判定该单元进入非线性变形状态，其中Ee为该单元非线性结点处材料的初始弹性模量。

[0038] 步骤二：根据由步骤一所得的各单元中非线性结点处的材料状态，更新矩阵K′r和并计算矩阵K″pr，设共有m个单元进入非线性变形状态(本例中，有m≤4)，则矩阵K″pr可表示为

[0039]

[0040] 其中，A为截面面积；l为单元长度；(.)i代表与第i个非线性结点相关的材料属性。

[0041] 步骤三：使用Woodbury公式计算结构增量荷载ΔF作用下的位移增量ΔX，即：

[0042]

[0043] 利用式(2)给出的矩阵Ke的三角分解形式，则上式涉及到的与矩阵 有关的乘法运算均可通过回代计算完成。

[0044] 步骤四：收敛判断并更新用于下一个计算步的荷载增量。此时需检查所得计算结果是否达到规定的误差范围之内，若收敛，则进行下一个增量步的计算，若不收敛，则继续进行迭代。