CN118171519A - 一种高效求解几何非线性问题的组合法 - Google Patents

Abstract

本发明公开了一种高效求解几何非线性问题的组合法，属于高性能计算仿真技术领域。所述方法首先需要建立待分析的几何模型，其次对待分析的几何模型进行任意单元类型的网格划分，然后设置待分析模型的边界条件及荷载情况。最后将上述条件输入至求解算法中。本发明将更新拉格朗日法、有限元方法、非线性算法、几何非线性理论结合起来，在现有非线性求解算法的理论基础上进行创新，以计算机卓越的数值计算能力为载体来求解几何非线性问题，在保证求解结果数值精度的前提下，提升几何非线性问题求解速度。

Description

一种高效求解几何非线性问题的组合法

技术领域

本发明属于高性能计算仿真技术领域，具体是一种高效求解几何非线性问题的组合法。

背景技术

结构变形问题广泛存在于工程领域，在机械制造、材料冶金、能源动力中常规或者极端环境下关键部件服役的可靠程度均有着一定程度的影响，而对于几何非线性问题的准确求解，可以克服构件大变形的劣势，针对大变形的结果进行改善，继而改进相应的构件性能。

结构变形情况可以分为小变形与大变形两种情况，弹性力学中利用小变形假设则可导出小变形假设下结构的有限元方程，该方程是线性微分方程，通过既定的边界条件则可完成解析解的求解。传统有限元理论的线性分析正是建立在弹性体小变形基础上，即通过已知构型的荷载及位移的线性关系可推算出任意构型下荷载及位移的关系，完成求解。而对于大变形问题，需要不断更新当前构型下构件几何变形对构件刚度的影响。因此几何大变形问题的求解，是多个以直逼曲的解法。通过多次迭代来求解一个超越方程，最终获得其数值解。而对于非线性求解算法，又应该具有较强的自适应性，即在非线性区域强化收敛性，而在非线性较弱的区域加快求解速度。

因此，本发明建立一种高效求解几何非线性问题的组合法，利用发展成熟的有限元数值计算理论、大变形理论、更新拉格朗日法、非线性数值算法求解及控制理论。着力于求解几何非线性问题，在准确求解几何非线性问题的同时，提升求解速度。不仅可以应用于高端装备大区域几何非线性分析，还可以在为梁壳单元非线性问题提供更加精确的求解方法，能够极好的促进材料和结构设计的力学性能优化。

发明内容

本发明公开了一种高效求解几何非线性问题的组合法，该方法将更新拉格朗日法、有限元方法、非线性算法、几何非线性理论结合起来，在现有非线性求解算法的理论基础上进行创新，以计算机卓越的数值计算能力为载体来求解几何非线性问题，在保证求解结果数值精度的前提下，提升了几何非线性问题的求解速度。

本发明解决上述技术问题所提供的技术方案是：一种高效求解几何非线性问题的组合法，包括以下步骤：

S1、建立待定几何非线性分析平面梁的几何模型或者几何装配体模型；

S2、对待分析几何模型进行空间离散化，即网格划分，生成计算所需的网格数据；

S3、对划分的网格模型进行边界条件、载荷类型大小、截面材料设定；

S4、设定组合法初始荷载因子数值的大小λ_initial、迭代容差值大小t、荷载因子符号S、设定广义刚度参数参考值GSP_ref、总荷载因子数值大小λ；

S5、使用组合法针对既定荷载边界条件下的问题求解，需计算当前参考荷载外力向量{F_ref}。进入荷载增量步预测阶段求解。计算当前构型下各单元的切线刚度矩阵[K_T]、弹性刚度矩阵[K_E]、几何刚度矩阵[K_G]、各单元截面面积A、长度L、单元夹角θ、截面惯性矩I并存储、单元内力{f_i}。求解当前增量步下的切线增量。

S6、在当前构型的刚度情况下求解并更新广义刚度参数GSP的大小，并记录当前荷载增量步的GSP大小，并与上一增量步的GSP数值进行比较。并判定当前广义刚度参数数值GSP与广义刚度参数参考值GSP_ref的关系。依据判定结果选择当前构型的校正算法。并更新增量步荷载因子及位移增量的值。

S7、进入校正求解阶段，首先更新当前构型下的外荷载向量{F_e}、内力向量{F_i}、不平衡力向量其次判断当前荷载增量步是否达到收敛标准。然后依据判断结果进入迭代校正求解流程。最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小；

S8、完成几何非线性求解，收集最终位移场的位移向量结果；提取节点位移各增量步向量分量及荷载因子值，绘制荷载位移曲线；记录求解器总耗时。

进一步的技术方案是，所述步骤S1具体实现方法为：

建立宏观尺度任意规模平面梁单元的几何模型或者几何装配体模型，建立一个连续的求解域；

进一步的技术方案是，所述步骤S2具体实现方法为：

S21、离散参数的设定；

结合大规模的几何连续模型的几何特征，设定几何模型的全局网格种子、部分特征几何边的局部种子，在部分几何特征较复杂的特征边处设置局部细化参数；

S22、对几何非线性问题进行数学建模，会产生偏微分方程(PDE)，数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解，其离散方程

K(u)＝F,F∈Rⁿ

式中：n是离散化系统的总自由度，K是在当前边界荷载问题下系统的线性方程组系数矩阵，F为荷载向量，可使用线性代数方法求解未知量u；

S23、建立离散化几何模型；

设定网格划分参数，使用网格划分工具开展几何模型网格划分，可实现不同形状的大规模几何模型的空间离散，将模型离散为梁单元的连续离散体，同时将模型离散后的数据信息写出至本地文件；

进一步的技术方案是，所述步骤S3具体实现方法为：

S31、设置下列材料参数，材料杨氏模量E、单元截面形状类型S_type、单元截面长度b，单元截面高度h，并将求解所需材料参数输出到本地文件；

S32、对离散后的模型施加边界条件及载荷，指定边界约束的几何位置、载荷大小以及载荷施加的几何位置。并将受约束载荷的几何位置集合信息输出到本地文件；

进一步的技术方案是，所述步骤S4具体实现方法为：

S41、改良弧长法求解时需设置初始荷载因子值的大小，荷载因子初始值会影响求解收敛性。迭代容差值t是校正阶段求解时迭代步收敛标准。荷载因子符号S初始状态默认为正号，并在求解过程中依据广义传导参数进行调整。荷载因子最大值λ_max为1，即满载情况—参考荷载向量{F_ref}等于{F_e}外力向量。预设初始荷载因子值为0.1、迭代容差值为10^-3、广义刚度参数参考值为10^-2。

进一步的技术方案是，所述步骤S5具体实现方法为：

S51、进入组合法求解器之前,需要读取模型信息，其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储。首先根据既定载荷以及边界情况，计算当前构型下的参考荷载热流率向量{F_ref}。其次从本地文件中读取既定的初始荷载因子λ_initial、迭代容差t、荷载因子最大值λ_max。读取梁单元截面材料杨氏模量E、单元截面形状类型S_type、单元截面长度高度b、h，单元夹角θ。当前构型下的初始单元内力{f_i}等于上一构型平衡时的内力。

具体计算步骤如下：

首先在第i个荷载增量步下第j个迭代步，读取既定荷载以及边界情况，将受约束自由度的边界情况及载荷大小进行存储；

其次依据有限元理论对荷载数值大小在所约束区域上进行计算。

式中：为第i个荷载增量步下第j个迭代步的整体刚度矩阵；Ω为单元域；N为平面梁单元形函数；Γ为模型离散后荷载加载区域。f为荷载数值。

S52、进入组合法求解器荷载增量步主循环，首先遍历模型单元，对离散后模型每个单元进行节点位移的提取，求解当前构型下的初始几何构型信息即求解各单元当前构型下的单元长度L，截面面积A，单元内力{f_i}。下面给出第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵的计算方法。

由虚功原理可得平面梁单元弹性刚度、几何刚度矩阵在求解域离散后的数值表达式如下：

式中：L为单元长度，A为单元截面面积，E为单元杨氏模量,I为单元截面惯性矩，P为单元轴力大小。W^L为虚功线性项，W^NL为虚功非线性项。

具体计算步骤如下：

首先计算当前平面梁单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s)下梁单元的形函数N；

其次形函数在母单元局部坐标系下对各个坐标分量求偏导，获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量

轴力杆单元的单元刚度矩阵计算方法如下：

对于轴力杆单元其形函数N为：

式中：N_i,j为平面杆单元节点自由度的形函数，r为局部坐标系下单元长度方向数值积分积分点坐标。

则杆单元的单元刚度矩阵[K_r]^e为：

对于弯曲梁单元其形函数N为：

N_i＝1-3r²+2r³,N_j＝(r-2r²+r³)l

N_k＝3r²-2r³,N_l＝(r³-r²)l

式中：N_i,j,k,l为平面梁单元节点自由度的形函数，r为局部坐标系下r坐标轴方向数值积分积分点坐标。

则弯曲梁单元的单元刚度矩阵[K_b]^e为：

则一般平面梁单元的弹性刚度矩阵[K_E]^e为：

则一般平面梁单元的几何刚度矩阵[K_G]^e为：

则第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵为：

另外为避免占用过多计算资源，存储该荷载增量步下计算所得的各单元弹性刚度矩阵方便后续单元内力的更新。

S53、将上述整体坐标系下的单元切线传导矩阵，保存其非0元素的位置索引及其对应位置矩阵元素的值，使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储；完成求解域在第i个荷载增量步下第j个迭代步的整体切线刚度矩阵的组装。

S54、通过置1法对第i个荷载增量步下第j个迭代步的整体切线刚度矩阵外部荷载热流率向量{F_e}进行处理。依据边界条件，在整体切线刚度矩阵及荷载向量中，将受约束自由度的行与列除主对角线置1外其余全为0，其意义为让受约束的自由度在求解时等于设定的边界值。该值可从模型离散后的文本文件将受约束的自由度和边界数值读取可得。

置1法具体形式如下：

S55、配置第三方库Eigen，设定线性方程组求解器为Cholesky分解求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。将第i个荷载增量步下第j个迭代步处理过后的整体切线刚度矩阵热流率向量{F}传入，调用语句compute(F)进行线性方程组的求解，求得预测阶段切线增量/>方程形式为：

S56、判断当前荷载增量步是否为首个荷载增量步，若当前增量步为首个荷载增量步时，则不需要计算弧长大小。若当前增量步不为首个荷载增量步时，需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量/>的二范数的大小，依据判定结果对算法进行加速。其中上标i为当前荷载增量步，下标1代表为预测阶段迭代步。

弧长计算公式为：ΔSⁱ＝{ΔU^i-1}^T{ΔU^i-1}式中上标i代表当前荷载增量步；

式中：{ΔU}^T上标T为转置含义，ΔU为荷载增量步位移增量；上标i代表荷载增量步步数，下标1代表预测阶段迭代步；

S57、若当前荷载增量步为首个荷载增量步，荷载因子大小为初始荷载因子大小λ_initial，记录荷载增量步预测阶段切线增量向量。若当前荷载增量步不为首个荷载增量步荷载因子/>具体的具体计算步骤如下:

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；T为转置符号。

进一步的技术方案是，所述步骤S6具体实现方法为：

S61、首先依据更新的第i个荷载增量步的切线增量计算当前荷载增量步的广义刚度参数GSP，具体计算步骤如下：

式中：为首个荷载增量步预测阶段所求得的切线增量；T为转置符号，为i-1荷载增量步预测阶段所求得的切线增量；

S62、比较相邻两个荷载增量步的广义刚度参数值，若相邻两个荷载增量步的广义刚度参数正负性不同，则更新求解过程中荷载因子符号S。并且对于求解过程中，广义刚度参数GSP在不断逼近荷载极限点时，由于切线斜率不断逼近于水平，相邻增量步切线增量趋近于无穷，GSP逐渐减小。直至相邻增量步切线增量不为锐角，此时GSP正负性发生改变。当越过荷载极限点时，相邻增量步切线增量又恢复为锐角，GSP数值逐渐增加。广义刚度参数数值在荷载极限点处有如下规律：

S63、依据求解过程中越过荷载极限点的个数及相邻增量步广义刚度参数的数值大小与既定的广义刚度参数参考数值大小进行判定，依据判定结果选择合适的迭代算法对其校正阶段求解进行加速。结构初始时刻默认为加载状态，则荷载因子符号S为负时构型为卸载状态，反之结构则为加载状态。则组合法校正阶段算法的加速策略如下所示：

S64、预测阶段荷载因子计算完成，保存当前荷载增量步预测阶段切线增量向量，更新增量步位移向量,以及荷载因子λ大小。结束荷载增量步预测阶段求解。

计算公式为：

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；上标T为转置符号

进一步的技术方案是，所述步骤S7具体实现方法为：

S71、首先依据更新第i个荷载增量步下第j个迭代步的荷载因子更新当前构型下需加载的外部荷载向量/>以及总的荷载因子λ。具体计算步骤如下；

式中：i为荷载增量步；j为迭代步；为荷载因子增量步增量；/>为荷载因子迭代步增量

S72、其次遍历单元，计算第i个荷载增量步下第j个迭代步当前构型下单元的弹性矩阵乘以局部坐标系下的节点相对位移量{u}及旋转变换矩阵[T]得到整体坐标系下的各单元内力向量/>将各分量组装至整体内力向量/>完成当前构型下各单元内力的求解。

具体计算步骤如下：

式中：线性传导矩阵的计算同S52所述，此处不再赘述。下标j-1是由于预测阶段并未进行该内力的计算。

S73、更新当前构型下的残差向量判断残差向量/>是否小于既定残差。此处的残差向量/>要依据边界条件对受约束的自由度位置进行处理，自由度受位移边界的的内力为0；

第i个荷载增量步下第j个迭代步的残差向量计算步骤如下：

若则程序跳出校正阶段循环，若荷载因子λ小于λ_max时，程序执行下一荷载增量步计算。反之，完成几何非线性问题求解。若/> 则进入迭代校正阶段。

S74、迭代校正阶段主要是为了消除系统内外不平衡力。依据增量步判定的算法类型对位移以及残差量采用弧长法或者荷载法的校正算法进行求解。首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是位移校正向量/>以及不平衡力校正向量/>求解，然后依据弧长约束方程或荷载方程更新荷载因子大小，最后更新总的位移向量以及总荷载因子大小。

依据增量步的判定结果可将弧长法校正阶段替换为荷载控制法校正阶段进行计算，加速校正阶段。第i个荷载增量步下第j个校正阶段荷载控制法的计算过程如下：

而对于恒定弧长法其及迭代校正过程的计算步骤如下：

算法中采用的约束方程为柱面约束方程，其方程形式为：

{ΔU^i-1}^T{ΔU^i-1}＝ΔS²

式中：i代表荷载增量步；ΔU^i-1表示上一荷载增量步位移增量；

对于第i个荷载增量步下第j个迭代步有限元方程如下：

而对于第i个荷载增量步下第j个迭代步的外力载荷如下：

而在第i个荷载增量步下第j个迭代步的位移增量如下：

因此在第i个荷载增量步下第j个迭代步方程的最终形式如下：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；代表不平衡力校正向量；

因此校正阶段求解公式为：

单元及整体切线传导矩阵的计算步骤同S52，此处不再赘述。

将参考荷载外力向量{F_ref},更新后的不平衡力向量依据边界条件进行处理，对自由度被约束的位置采用置1法进行处理。处理方法同S54，此处不再赘述。

配置第三方库Eigen，设定线性方程组求解器为Cholesky求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。传入即可求得位移校正向量/>同理可求得不平衡力校正向量/>

弧长法迭代校正阶段荷载因子求解方法如下所示：

在荷载增量步i下第j迭代步的位移增量可表示为如下形式：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；代表不平衡力校正向量；

位移增量代入到柱面约束方程中可得以当前迭代步荷载因子/>为变量的一元二次方程。因此在校正阶段时关于荷载因子的一元二次方程为：

则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为：

式中：i代表荷载增量步步数；j代表迭代步步数；代表预测阶段切线增量；代表位移校正向量；/>代表不平衡力校正向量；

则对于弧长法校正阶段的算法可表示为：

S75、依据S63中的判定结果选择校正阶段计算方式，计算荷载因子更新增量步位移向量以及增量步荷载因子Δλ大小。直至残差满足小于预设迭代容差停止迭代。并判断总荷载因子大小λ是否大于预设的荷载因子最大值λ_max。

进一步的技术方案是，所述步骤S8具体实现方法为：

S81、当总荷载因子数值大小大于预设最大荷载因子大小时，跳出荷载增量步主循环，记录当前构型下总位移向量。求解器求解完成并存储当前模型各个节点的位移向量。

S82、选择梁单元模型中某一点，提取其各增量步内位移分量及荷载因子数值大小，以荷载因子数值为纵坐标，位移分量为横坐标。绘制荷载因子曲线并输出。同时依据初始构型下及最终构型，绘制单元的变形效果图。最终完成该模型的几何非线性问题求解。

本发明的有益效果是：本发明基于自编C++程序完成有限元方法下几何模型的建立，并以梁单元类型对几何模型进行离散。然后将模型信息导入到自主研发的求解器中。算法依据求解中刚度参数数值的变化来甄别当前位置是否处于极限点位置，对于极限点位置采用收敛性强的求解算法，而对于非极限点区域，则采用求解速度较快的求解算法。从而实现在准确求解几何非线性问题的前提下，不使用并行技术对算法求解过程进行加速，缩短几何非线性问题的求解耗时。保证了结果的数值精度，提升了几何非线性问题求解速度。对求解速度有高要求的几何非线性计算场景提供了解决方案。

附图说明

图1是本发明的一种高效求解几何非线性问题的组合法的流程图；

图2是本发明的几何计算模型—Lee Frame示意图；

图3是离散的单元类型示意图；

图4是几何模型离散后的网格图；

图5是组合法求解几何非线性问题的最终变形效果图；

图6是组合法求解几何非线性问题的荷载位移曲线图；

图7是组合法求解几何非线性问题的GSP位移曲线与荷载位移曲线图；

表1是组合法求解几何非线性问题耗时表

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明提供了一种高效求解几何非线性问题的组合法，包括以下步骤：

S1、建立待定几何非线性分析平面梁的几何模型或者几何装配体模型；几何模型示意图如图2所示。

具体的实现步骤为：

建立宏观尺度任意规模平面梁单元的几何模型或者几何装配体模型，建立一个连续的求解域；

S2、对待分析几何模型进行空间离散化，即网格划分，生成计算所需的网格数据；

具体的实现步骤为：

S21、离散参数的设定；

结合大规模的几何连续模型的几何特征，设定几何模型的全局网格种子、部分特征几何边的局部种子，在部分几何特征较复杂的特征边处设置局部细化参数；单元示意图如图3所示，几何模型离散后的网格图如图4所示。

S22、对几何非线性问题进行数学建模，会产生偏微分方程(PDE)，数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解，其离散方程

K(u)＝F,F∈Rⁿ

式中：n是离散化系统的总自由度，K是在当前边界荷载问题下系统的线性方程组系数矩阵，F为荷载向量，可使用线性代数方法求解未知量u；

S23、建立离散化几何模型；

设定网格划分参数，使用网格划分工具开展几何模型网格划分，可实现不同形状的大规模几何模型的空间离散，将模型离散为梁单元的连续离散体，同时将模型离散后的数据信息写出至本地文件；

S3、对划分的网格模型进行边界条件、载荷类型大小、截面材料设定；

具体的实现步骤为：

S31、设置下列材料参数，材料杨氏模量E、单元截面形状类型S_type、单元截面长度b，单元截面高度h，并将求解所需材料参数输出到本地文件；

S32、对离散后的模型施加边界条件及载荷，指定边界约束的几何位置、载荷大小以及载荷施加的几何位置。并将受约束载荷的几何位置集合信息输出到本地文件；

S4、设定组合法初始荷载因子数值的大小λ_initial、迭代容差值大小t、荷载因子符号S、设定广义刚度参数参考值GSP_ref、总荷载因子数值大小λ；

具体的实现步骤为：

S41、改良弧长法求解时需设置初始荷载因子值的大小，荷载因子初始值会影响求解收敛性。迭代容差值t是校正阶段求解时迭代步收敛标准。荷载因子符号S初始状态默认为正号，并在求解过程中依据广义传导参数进行调整。荷载因子最大值λ_max为1，即满载情况—参考荷载向量{F_ref}等于{F_e}外力向量。预设初始荷载因子值为0.1、迭代容差值为10^-3、广义刚度参数参考值为10^-2。

S5、使用组合法针对既定荷载边界条件下的问题求解，需计算当前参考荷载外力向量{F_ref}。进入荷载增量步预测阶段求解。计算当前构型下各单元的切线刚度矩阵[K_T]、弹性刚度矩阵[K_E]、几何刚度矩阵[K_G]、各单元截面面积A、长度L、单元夹角θ、截面惯性矩I并存储、单元内力{f_i}。求解当前增量步下的切线增量。

具体的实现步骤为：

S51、进入组合法求解器之前,需要读取模型信息，其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储。首先根据既定载荷以及边界情况，计算当前构型下的参考荷载热流率向量{F_ref}。其次从本地文件中读取既定的初始荷载因子λ_initial、迭代容差t、荷载因子最大值λ_max。读取梁单元截面材料杨氏模量E、单元截面形状类型S_type、单元截面长度高度b、h，单元夹角θ。当前构型下的初始单元内力{f_i}等于上一构型平衡时的内力。

具体计算步骤如下：

首先在第i个荷载增量步下第j个迭代步，读取既定荷载以及边界情况，将受约束自由度的边界情况及载荷大小进行存储；

其次依据有限元理论对荷载数值大小在所约束区域上进行计算。

式中：为第i个荷载增量步下第j个迭代步的整体刚度矩阵；Ω为单元域；N为平面梁单元形函数；Γ为模型离散后荷载加载区域。f为荷载数值。

S52、进入组合法求解器荷载增量步主循环，首先遍历模型单元，对离散后模型每个单元进行节点位移的提取，求解当前构型下的初始几何构型信息即求解各单元当前构型下的单元长度L，截面面积A，单元内力{f_i}。下面给出第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵的计算方法。

由虚功原理可得平面梁单元弹性刚度、几何刚度矩阵在求解域离散后的数值表达式如下：

/>

式中：L为单元长度，A为单元截面面积，E为单元杨氏模量,I为单元截面惯性矩，P为单元轴力大小。W^L为虚功线性项，W^NL为虚功非线性项。

具体计算步骤如下：

首先计算当前平面梁单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s)下梁单元的形函数N；

其次形函数在母单元局部坐标系下对各个坐标分量求偏导，获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量

轴力杆单元的单元刚度矩阵计算方法如下：

对于轴力杆单元其形函数N为：

式中：N_i,j为平面杆单元节点自由度的形函数，r为局部坐标系下单元长度方向数值积分积分点坐标。

则杆单元的单元刚度矩阵[K_r]^e为：

对于弯曲梁单元其形函数N为：

N_i＝1-3r²+2r³,N_j＝(r-2r²+r³)l

N_k＝3r²-2r³,N_l＝(r³-r²)l

式中：N_i,j,k,l为平面梁单元节点自由度的形函数，r为局部坐标系下r坐标轴方向数值积分积分点坐标。

则弯曲梁单元的单元刚度矩阵[K_b]^e为：

则一般平面梁单元的弹性刚度矩阵[K_E]^e为：

则一般平面梁单元的几何刚度矩阵[K_G]^e为：

则第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵为：

另外为避免占用过多计算资源，存储该荷载增量步下计算所得的各单元弹性刚度矩阵方便后续单元内力的更新。

S53、将上述整体坐标系下的单元切线传导矩阵，保存其非0元素的位置索引及其对应位置矩阵元素的值，使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储；完成求解域在第i个荷载增量步下第j个迭代步的整体切线刚度矩阵的组装。

S54、通过置1法对第i个荷载增量步下第j个迭代步的整体切线刚度矩阵外部荷载热流率向量{F_e}进行处理。依据边界条件，在整体切线刚度矩阵及荷载向量中，将受约束自由度的行与列除主对角线置1外其余全为0，其意义为让受约束的自由度在求解时等于设定的边界值。该值可从模型离散后的文本文件将受约束的自由度和边界数值读取可得。/>

置1法具体形式如下：

S55、配置第三方库Eigen，设定线性方程组求解器为Cholesky分解求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。将第i个荷载增量步下第j个迭代步处理过后的整体切线刚度矩阵热流率向量{F}传入，调用语句compute(F)进行线性方程组的求解，求得预测阶段切线增量/>方程形式为：

S56、判断当前荷载增量步是否为首个荷载增量步，若当前增量步为首个荷载增量步时，则不需要计算弧长大小。若当前增量步不为首个荷载增量步时，需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量/>的二范数的大小，依据判定结果对算法进行加速。其中上标i为当前荷载增量步，下标1代表为预测阶段迭代步。

弧长计算公式为：ΔSⁱ＝{ΔU^i-1}^T{ΔU^i-1}式中上标i代表当前荷载增量步；

式中：{ΔU}^T上标T为转置含义，ΔU为荷载增量步位移增量；上标i代表荷载增量步步数，下标1代表预测阶段迭代步；

S57、若当前荷载增量步为首个荷载增量步，荷载因子大小为初始荷载因子大小λ_initial，记录荷载增量步预测阶段切线增量向量。若当前荷载增量步不为首个荷载增量步荷载因子/>具体的具体计算步骤如下:

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；T为转置符号。

S6、在当前构型的刚度情况下求解并更新广义刚度参数GSP的大小，并记录当前荷载增量步的GSP大小，并与上一增量步的GSP数值进行比较。并判定当前广义刚度参数数值GSP与广义刚度参数参考值GSP_ref的关系。依据判定结果选择当前构型的校正算法。并更新增量步荷载因子及位移增量的值。

具体的实现步骤为：

S61、首先依据更新的第i个荷载增量步的切线增量计算当前荷载增量步的广义刚度参数GSP，具体计算步骤如下：

式中：为首个荷载增量步预测阶段所求得的切线增量；T为转置符号，为i-1荷载增量步预测阶段所求得的切线增量；

S62、比较相邻两个荷载增量步的广义刚度参数值，若相邻两个荷载增量步的广义刚度参数正负性不同，则更新求解过程中荷载因子符号S。并且对于求解过程中，广义刚度参数GSP在不断逼近荷载极限点时，由于切线斜率不断逼近于水平，相邻增量步切线增量趋近于无穷，GSP逐渐减小。直至相邻增量步切线增量不为锐角，此时GSP正负性发生改变。当越过荷载极限点时，相邻增量步切线增量又恢复为锐角，GSP数值逐渐增加。广义刚度参数数值在荷载极限点处有如下规律：

S63、依据求解过程中越过荷载极限点的个数及相邻增量步广义刚度参数的数值大小与既定的广义刚度参数参考数值大小进行判定，依据判定结果选择合适的迭代算法对其校正阶段求解进行加速。结构初始时刻默认为加载状态，则荷载因子符号S为负时构型为卸载状态，反之结构则为加载状态。则组合法校正阶段算法的加速策略如下所示：

S64、预测阶段荷载因子计算完成，保存当前荷载增量步预测阶段切线增量向量，更新增量步位移向量,以及荷载因子λ大小。结束荷载增量步预测阶段求解。

计算公式为：

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；上标T为转置符号

S7、进入校正求解阶段，首先更新当前构型下的外荷载向量{F_e}、内力向量{F_i}、不平衡力向量其次判断当前荷载增量步是否达到收敛标准。然后依据判断结果进入迭代校正求解流程。最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小；

具体的实现步骤为：

S71、首先依据更新第i个荷载增量步下第j个迭代步的荷载因子更新当前构型下需加载的外部荷载向量/>以及总的荷载因子λ。具体计算步骤如下；

式中：i为荷载增量步；j为迭代步；为荷载因子增量步增量；/>为荷载因子迭代步增量

S72、其次遍历单元，计算第i个荷载增量步下第j个迭代步当前构型下单元的弹性矩阵乘以局部坐标系下的节点相对位移量{u}及旋转变换矩阵[T]得到整体坐标系下的各单元内力向量/>将各分量组装至整体内力向量/>完成当前构型下各单元内力的求解。

具体计算步骤如下：

式中：线性传导矩阵的计算同S52所述，此处不再赘述。下标j-1是由于预测阶段并未进行该内力的计算。

S73、更新当前构型下的残差向量判断残差向量/>是否小于既定残差。此处的残差向量/>要依据边界条件对受约束的自由度位置进行处理，自由度受位移边界的的内力为0；

第i个荷载增量步下第j个迭代步的残差向量计算步骤如下：

若则程序跳出校正阶段循环，若荷载因子λ小于λ_max时，程序执行下一荷载增量步计算。反之，完成几何非线性问题求解。若/> 则进入迭代校正阶段。/>

S74、迭代校正阶段主要是为了消除系统内外不平衡力。依据增量步判定的算法类型对位移以及残差量采用弧长法或者荷载法的校正算法进行求解。首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是位移校正向量/>以及不平衡力校正向量/>求解，然后依据弧长约束方程或荷载方程更新荷载因子大小，最后更新总的位移向量以及总荷载因子大小。

依据增量步的判定结果可将弧长法校正阶段替换为荷载控制法校正阶段进行计算，加速校正阶段。第i个荷载增量步下第j个校正阶段荷载控制法的计算过程如下：

而对于恒定弧长法其及迭代校正过程的计算步骤如下：

算法中采用的约束方程为柱面约束方程，其方程形式为：

{ΔU^i-1}^T{ΔU^i-1}＝ΔS²

式中：i代表荷载增量步；ΔU^i-1表示上一荷载增量步位移增量；

对于第i个荷载增量步下第j个迭代步有限元方程如下：

而对于第i个荷载增量步下第j个迭代步的外力载荷如下：

而在第i个荷载增量步下第j个迭代步的位移增量如下：

因此在第i个荷载增量步下第j个迭代步方程的最终形式如下：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；代表不平衡力校正向量；

因此校正阶段求解公式为：

单元及整体切线传导矩阵的计算步骤同S52，此处不再赘述。

将参考荷载外力向量{F_ref},更新后的不平衡力向量依据边界条件进行处理，对自由度被约束的位置采用置1法进行处理。处理方法同S54，此处不再赘述。/>

配置第三方库Eigen，设定线性方程组求解器为Cholesky求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。传入即可求得位移校正向量/>同理可求得不平衡力校正向量/>

弧长法迭代校正阶段荷载因子求解方法如下所示：

在荷载增量步i下第j迭代步的位移增量可表示为如下形式：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；代表不平衡力校正向量；

位移增量代入到柱面约束方程中可得以当前迭代步荷载因子/>为变量的一元二次方程。因此在校正阶段时关于荷载因子的一元二次方程为：

则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为：

式中：i代表荷载增量步步数；j代表迭代步步数；代表预测阶段切线增量；代表位移校正向量；/>代表不平衡力校正向量；

则对于弧长法校正阶段的算法可表示为：

S75、依据S63中的判定结果选择校正阶段计算方式，计算荷载因子更新增量步位移向量以及增量步荷载因子Δλ大小。直至残差满足小于预设迭代容差停止迭代。并判断总荷载因子大小λ是否大于预设的荷载因子最大值λ_max。

S8、完成几何非线性求解，收集最终位移场的位移向量结果；提取节点位移各增量步向量分量及荷载因子值，绘制荷载位移曲线；记录求解器总耗时。

具体的实现步骤为：

S81、当总荷载因子数值大小大于预设最大荷载因子大小时，跳出荷载增量步主循环，记录当前构型下总位移向量。求解器求解完成并存储当前模型各个节点的位移向量。

S82、选择梁单元模型中某一点，提取其各增量步内位移分量及荷载因子数值大小，以荷载因子数值为纵坐标，位移分量为横坐标。绘制荷载因子曲线并输出。同时依据初始构型下及最终构型，绘制单元的变形效果图。最终完成该模型的几何非线性问题求解。

本发明公开了一种高效求解几何非线性问题的组合法，该方法将更新拉格朗日法、有限元方法、非线性算法、几何非线性理论结合起来，在现有非线性求解算法的理论基础上进行创新，以计算机卓越的数值计算能力为载体来求解几何非线性问题，在保证求解结果数值精度的前提下，提升了几何非线性问题的求解速度。

表1

以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，可利用上述揭示的技术内容做出变动或修饰为等同变化的等效实施例。但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (9)

1.一种高效求解几何非线性问题的组合法，其特征在于，包括以下步骤：

S1、建立待定几何非线性分析平面梁的几何模型或者几何装配体模型；

S2、对待分析几何模型进行空间离散化，即网格划分，生成计算所需的网格数据；

S3、对划分的网格模型进行边界条件、载荷类型大小、截面材料设定；

S4、设定组合法初始荷载因子数值的大小λ_initial、迭代容差值大小t、荷载因子符号S、设定广义刚度参数参考值GSP_ref、总荷载因子数值大小λ；

S5、使用组合法针对既定荷载边界条件下的问题求解，需计算当前参考荷载外力向量{F_ref}。进入荷载增量步预测阶段求解。计算当前构型下各单元的切线刚度矩阵[K_T]、弹性刚度矩阵[K_E]、几何刚度矩阵[K_G]、各单元截面面积A、长度L、单元夹角θ、截面惯性矩I并存储、单元内力{f_i}。求解当前增量步下的切线增量。

S6、在当前构型的刚度情况下求解并更新广义刚度参数GSP的大小，并记录当前荷载增量步的GSP大小，并与上一增量步的GSP数值进行比较。并判定当前广义刚度参数数值GSP与广义刚度参数参考值GSP_ref的关系。依据判定结果选择当前构型的校正算法。并更新增量步荷载因子及位移增量的值。

S7、进入校正求解阶段，首先更新当前构型下的外荷载向量{F_e}、内力向量{F_i}、不平衡力向量其次判断当前荷载增量步是否达到收敛标准。然后依据判断结果进入迭代校正求解流程。最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小；

S8、完成几何非线性求解，收集最终位移场的位移向量结果；提取节点位移各增量步向量分量及荷载因子值，绘制荷载位移曲线；记录求解器总耗时。

2.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S1具体实现方法为：建立宏观尺度任意规模平面梁单元的几何模型或者几何装配体模型，建立一个连续的求解域。

3.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S2具体实现方法为：

S21、离散参数的设定；

结合大规模的几何连续模型的几何特征，设定几何模型的全局网格种子、部分特征几何边的局部种子，在部分几何特征较复杂的特征边处设置局部细化参数；

S22、对几何非线性问题进行数学建模，会产生偏微分方程(PDE)，数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解，其离散方程

K(u)＝F,F∈Rⁿ

式中：n是离散化系统的总自由度，K是在当前边界荷载问题下系统的线性方程组系数矩阵，F为荷载向量，可使用线性代数方法求解未知量u；

S23、建立离散化几何模型；

设定网格划分参数，使用网格划分工具开展几何模型网格划分，可实现不同形状的大规模几何模型的空间离散，将模型离散为梁单元的连续离散体，同时将模型离散后的数据信息写出至本地文件。

4.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S3具体实现方法为：

S31、设置下列材料参数，材料杨氏模量E、单元截面形状类型S_type、单元截面长度b，单元截面高度h，并将求解所需材料参数输出到本地文件；

S32、对离散后的模型施加边界条件及载荷，指定边界约束的几何位置、载荷大小以及载荷施加的几何位置。并将受约束载荷的几何位置集合信息输出到本地文件。

5.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S4具体实现方法为：

S41、改良弧长法求解时需设置初始荷载因子值的大小，荷载因子初始值会影响求解收敛性。迭代容差值t是校正阶段求解时迭代步收敛标准。荷载因子符号S初始状态默认为正号，并在求解过程中依据广义传导参数进行调整。荷载因子最大值λ_max为1，即满载情况—参考荷载向量{F_ref}等于{F_e}外力向量。预设初始荷载因子值为0.1、迭代容差值为10^-3、广义刚度参数参考值为10^-2。

6.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S5具体实现方法为：

S51、进入组合法求解器之前,需要读取模型信息，其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储。首先根据既定载荷以及边界情况，计算当前构型下的参考荷载热流率向量{F_ref}。其次从本地文件中读取既定的初始荷载因子λ_initial、迭代容差t、荷载因子最大值λ_max。读取梁单元截面材料杨氏模量E、单元截面形状类型S_type、单元截面长度高度b、h，单元夹角θ。当前构型下的初始单元内力{f_i}等于上一构型平衡时的内力。

具体计算步骤如下：

首先在第i个荷载增量步下第j个迭代步，读取既定荷载以及边界情况，将受约束自由度的边界情况及载荷大小进行存储；

其次依据有限元理论对荷载数值大小在所约束区域上进行计算。

式中：为第i个荷载增量步下第j个迭代步的整体刚度矩阵；Ω为单元域；N为平面梁单元形函数；Γ为模型离散后荷载加载区域。f为荷载数值。

S52、进入组合法求解器荷载增量步主循环，首先遍历模型单元，对离散后模型每个单元进行节点位移的提取，求解当前构型下的初始几何构型信息即求解各单元当前构型下的单元长度L，截面面积A，单元内力{f_i}。下面给出第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵的计算方法。

由虚功原理可得平面梁单元弹性刚度、几何刚度矩阵在求解域离散后的数值表达式如下：

式中：L为单元长度，A为单元截面面积，E为单元杨氏模量,I为单元截面惯性矩，P为单元轴力大小。W^L为虚功线性项，W^NL为虚功非线性项。

具体计算步骤如下：

首先计算当前平面梁单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s)下梁单元的形函数N；

其次形函数在母单元局部坐标系下对各个坐标分量求偏导，获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量

轴力杆单元的单元刚度矩阵计算方法如下：

对于轴力杆单元其形函数N为：

式中：N_i,j为平面杆单元节点自由度的形函数，r为局部坐标系下单元长度方向数值积分积分点坐标。

则杆单元的单元刚度矩阵[K_r]^e为：

对于弯曲梁单元其形函数N为：

N_i＝1-3r²+2r³,N_j＝(r-2r²+r³)l

N_k＝3r²-2r³,N_l＝(r³-r²)l

式中：N_i,j,k,l为平面梁单元节点自由度的形函数，r为局部坐标系下r坐标轴方向数值积分积分点坐标。

则弯曲梁单元的单元刚度矩阵[K_b]^e为：

则一般平面梁单元的弹性刚度矩阵[K_E]^e为：

则一般平面梁单元的几何刚度矩阵[K_G]^e为：

则第i个荷载增量步的第j个迭代步下的单元切线刚度矩阵为：

另外为避免占用过多计算资源，存储该荷载增量步下计算所得的各单元弹性刚度矩阵方便后续单元内力的更新。

S53、将上述整体坐标系下的单元切线传导矩阵，保存其非0元素的位置索引及其对应位置矩阵元素的值，使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储；完成求解域在第i个荷载增量步下第j个迭代步的整体切线刚度矩阵的组装。

S54、通过置1法对第i个荷载增量步下第j个迭代步的整体切线刚度矩阵外部荷载热流率向量{F_e}进行处理。依据边界条件，在整体切线刚度矩阵及荷载向量中，将受约束自由度的行与列除主对角线置1外其余全为0，其意义为让受约束的自由度在求解时等于设定的边界值。该值可从模型离散后的文本文件将受约束的自由度和边界数值读取可得。

置1法具体形式如下：

S55、配置第三方库Eigen，设定线性方程组求解器为Cholesky分解求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。将第i个荷载增量步下第j个迭代步处理过后的整体切线刚度矩阵热流率向量{F}传入，调用语句进行线性方程组的求解，求得预测阶段切线增量/>方程形式为：

S56、判断当前荷载增量步是否为首个荷载增量步，若当前增量步为首个荷载增量步时，则不需要计算弧长大小。若当前增量步不为首个荷载增量步时，需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量/>的二范数的大小，依据判定结果对算法进行加速。其中上标i为当前荷载增量步，下标1代表为预测阶段迭代步。

弧长计算公式为：ΔSⁱ＝{ΔU^i-1}^T{ΔU^i-1}式中上标i代表当前荷载增量步；

式中：{ΔU}^T上标T为转置含义，ΔU为荷载增量步位移增量；上标i代表荷载增量步步数，下标1代表预测阶段迭代步；

S57、若当前荷载增量步为首个荷载增量步，荷载因子大小为初始荷载因子大小λ_initial，记录荷载增量步预测阶段切线增量向量。若当前荷载增量步不为首个荷载增量步荷载因子/>具体的具体计算步骤如下:

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；T为转置符号。

7.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S6具体实现方法为：

S61、首先依据更新的第i个荷载增量步的切线增量计算当前荷载增量步的广义刚度参数GSP，具体计算步骤如下：

式中：为首个荷载增量步预测阶段所求得的切线增量；T为转置符号，为i-1荷载增量步预测阶段所求得的切线增量；

S62、比较相邻两个荷载增量步的广义刚度参数值，若相邻两个荷载增量步的广义刚度参数正负性不同，则更新求解过程中荷载因子符号S。并且对于求解过程中，广义刚度参数GSP在不断逼近荷载极限点时，由于切线斜率不断逼近于水平，相邻增量步切线增量趋近于无穷，GSP逐渐减小。直至相邻增量步切线增量不为锐角，此时GSP正负性发生改变。当越过荷载极限点时，相邻增量步切线增量又恢复为锐角，GSP数值逐渐增加。广义刚度参数数值在荷载极限点处有如下规律：

S63、依据求解过程中越过荷载极限点的个数及相邻增量步广义刚度参数的数值大小与既定的广义刚度参数参考数值大小进行判定，依据判定结果选择合适的迭代算法对其校正阶段求解进行加速。结构初始时刻默认为加载状态，则荷载因子符号S为负时构型为卸载状态，反之结构则为加载状态。则组合法校正阶段算法的加速策略如下所示：

S64、预测阶段荷载因子计算完成，保存当前荷载增量步预测阶段切线增量向量，更新增量步位移向量,以及荷载因子λ大小。结束荷载增量步预测阶段求解。

计算公式为：

式中：为预测阶段所求得的切线增量；其中i为荷载增量步；上标T为转置符号。

8.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S7具体实现方法为：

S71、首先依据更新第i个荷载增量步下第j个迭代步的荷载因子更新当前构型下需加载的外部荷载向量/>以及总的荷载因子λ。具体计算步骤如下；

式中：i为荷载增量步；j为迭代步；为荷载因子增量步增量；/>为荷载因子迭代步增量

S72、其次遍历单元，计算第i个荷载增量步下第j个迭代步当前构型下单元的弹性矩阵乘以局部坐标系下的节点相对位移量{u}及旋转变换矩阵[T]得到整体坐标系下的各单元内力向量/>将各分量组装至整体内力向量/>完成当前构型下各单元内力的求解。

具体计算步骤如下：

式中：线性传导矩阵的计算同S52所述，此处不再赘述。下标j-1是由于预测阶段并未进行该内力的计算。

S73、更新当前构型下的残差向量判断残差向量/>是否小于既定残差。此处的残差向量/>要依据边界条件对受约束的自由度位置进行处理，自由度受位移边界的的内力为0；

第i个荷载增量步下第j个迭代步的残差向量计算步骤如下：

若则程序跳出校正阶段循环，若荷载因子λ小于λ_max时，程序执行下一荷载增量步计算。反之，完成几何非线性问题求解。若/> 则进入迭代校正阶段。

S74、迭代校正阶段主要是为了消除系统内外不平衡力。依据增量步判定的算法类型对位移以及残差量采用弧长法或者荷载法的校正算法进行求解。首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是位移校正向量/>以及不平衡力校正向量/>求解，然后依据弧长约束方程或荷载方程更新荷载因子大小，最后更新总的位移向量以及总荷载因子大小。

依据增量步的判定结果可将弧长法校正阶段替换为荷载控制法校正阶段进行计算，加速校正阶段。第i个荷载增量步下第j个校正阶段荷载控制法的计算过程如下：

而对于恒定弧长法其及迭代校正过程的计算步骤如下：

算法中采用的约束方程为柱面约束方程，其方程形式为：

{ΔU^i-1}^T{ΔU^i-1}＝ΔS²

式中：i代表荷载增量步；ΔU^i-1表示上一荷载增量步位移增量；

对于第i个荷载增量步下第j个迭代步有限元方程如下：

而对于第i个荷载增量步下第j个迭代步的外力载荷如下：

而在第i个荷载增量步下第j个迭代步的位移增量如下：

因此在第i个荷载增量步下第j个迭代步方程的最终形式如下：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；/>代表不平衡力校正向量；

因此校正阶段求解公式为：

单元及整体切线传导矩阵的计算步骤同S52，此处不再赘述。

将参考荷载外力向量{F_ref},更新后的不平衡力向量依据边界条件进行处理，对自由度被约束的位置采用置1法进行处理。处理方法同S54，此处不再赘述。

配置第三方库Eigen，设定线性方程组求解器为Cholesky求解器<SimplicialLLT>，初始化乔利斯基分解求解器，设置求解器数值精度大小。传入即可求得位移校正向量/>同理可求得不平衡力校正向量/>

弧长法迭代校正阶段荷载因子求解方法如下所示：

在荷载增量步i下第j迭代步的位移增量可表示为如下形式：

式中：i代表荷载增量步步数；j代表迭代步步数；代表位移校正向量；/>代表不平衡力校正向量；

位移增量代入到柱面约束方程中可得以当前迭代步荷载因子/>为变量的一元二次方程。因此在校正阶段时关于荷载因子的一元二次方程为：

则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为：

式中：i代表荷载增量步步数；j代表迭代步步数；代表预测阶段切线增量；代表位移校正向量；/>代表不平衡力校正向量；

则对于弧长法校正阶段的算法可表示为：

S75、依据S63中的判定结果选择校正阶段计算方式，计算荷载因子更新增量步位移向量以及增量步荷载因子Δλ大小。直至残差满足小于预设迭代容差停止迭代。并判断总荷载因子大小λ是否大于预设的荷载因子最大值λ_max。

9.根据权利要求1所述的一种高效求解几何非线性问题的组合法，其特征在于，所述步骤S8具体实现方法为：

S81、当总荷载因子数值大小大于预设最大荷载因子大小时，跳出荷载增量步主循环，记录当前构型下总位移向量。求解器求解完成并存储当前模型各个节点的位移向量。

S82、选择梁单元模型中某一点，提取其各增量步内位移分量及荷载因子数值大小，以荷载因子数值为纵坐标，位移分量为横坐标。绘制荷载因子曲线并输出。同时依据初始构型下及最终构型，绘制单元的变形效果图。最终完成该模型的几何非线性问题求解。