结构非线性分析中，通常仅有少数构件进入非线性，其余大部分构件仍保持弹性。针对局部非线性问题，许多学者提出了高效分析方法，如子结构方法[1]、多尺度方法[2]、拟力法[3－4]以及基于Woodbury的非线性分析方法[5－6]等。Woodbury公式是对仅有少数元素变化的矩阵快速求逆的数学方法，最初被应用在结构重分析领域中。重分析[7－8]是指当结构参数(构件形状、尺寸、材料等)改变时，通过利用结构修改前的信息对结构进行快速分析的方法，该方法避免了对结构重复的完全分析，可高效精确地求解修改后结构的响应。目前，重分析方法已被应用到不同力学问题中，如结构设计[9]、动静态分析[10－11]、灵敏度分析[12－13]以及非线性分析[14－17]等。对于结构局部非线性问题，每个增量步中切线刚度矩阵仅少数与非线性相关的元素发生改变，可表示为弹性刚度矩阵与其低秩修正矩阵的和，符合Woodbury公式的应用条件，因此利用弹性刚度矩阵的逆矩阵可高效求解切线刚度矩阵的逆，避免了传统非线性分析中对切线刚度矩阵的重新合成和分解，提高了局部非线性分析的效率。例如，隔离非线性有限元方法是基于拟力法[3,18－19]的思想，并与有限元方法结合实现了结构精细化分析[6,20－22]，该方法通过对材料应变进行弹塑性分解，在每个增量步中将结构切线刚度矩阵等效成弹性刚度矩阵与塑性相关的修正矩阵的和，因此可利用Woodbury公式求解结构响应，实现了局部非线性的高效分析。目前对基于Woodbury非线性方法的效率评价[5,23]仅从增量步中线性方程组的求解进行分析，而通常在增量步中还需要迭代计算消除线性化引入的误差，总体的分析效率和非线性迭代算法的性能密切相关。若迭代算法收敛性差或收敛速度慢，分析计算代价高的问题仍无法避免。因此，研究基于Woodbury的非线性方法的迭代求解性能是必要的。
在传统结构非线性分析中，Newton-Raphson(N-R)法是求解非线性方程的经典迭代算法，该方法具有二阶收敛性，由于每个迭代步中都要对结构切线刚度矩阵更新和分解，当非线性程度较高时，迭代次数增多，尤其对于大规模结构会显著降低分析效率。而修正牛顿法[24]仅在增量步的第一个迭代步中更新和分解刚度矩阵，有效降低了每个迭代步的计算量，但该方法仅为线性收敛，且仅对弱非线性问题适用，对于刚度矩阵实时变化大的强非线性问题，该方法的收敛速率慢甚至不收敛。拟牛顿法[22]的收敛性介于N-R法和修正牛顿法之间，具有超线性收敛性，该方法在每个迭代步中对原有的柔度矩阵修正而不必重新合成和分解切线刚度矩阵，也不必沿用已有的刚度矩阵，但缺点是求解柔度矩阵的过程可能会破坏刚度矩阵的带状特性，使迭代计算偏离平衡路径[25]。上述3种迭代方法是结构非线性分析中最常用的计算方法，Bathe等[26]等在早期分析了修正牛顿法和拟牛顿法中Broyden-Fletcher-Goldfarb-Shanno(BFGS)方法的求解性能，结合算例和相关经验给出如下建议：在结构刚度下降(材料软化)的情况下，适合选用修正牛顿法；在复杂材料非线性(混凝土材料)或结构刚度缓慢增加的情况下，适合选用BFGS方法。但该文献中仅通过定性分析给出上述结论，并没有对算法效率进行理论定量分析。由于非线性分析过程的复杂性和耗时性，已有的迭代算法都有各自的适用性和局限性，因此许多学者为改进结构非线性分析的精度和效率，针对不同非线性问题提出了更有效的迭代算法。Scott和Fenves[25]基于Krylov子空间法提出了加速的牛顿算法求解非线性平衡方程，该方法通过用低秩最小二乘分析法改善了结构非线性平衡路径的搜索方向，对框架结构动力连续倒塌等强非线性问题有良好的求解性能，算例结果表明该方法比N-R法的收敛半径大，收敛速度更快。Kim等[27]结合神经网络法提出了一种预测校正方法进行结构非线性分析，其主要思想是先用神经网络算法对迭代初始点进行预测，再用N-R法从该预测点迭代求解，提高了收敛速度，且该方法适用于求解几何非线性等复杂变形的非线性问题。Saffari等[28－29]提出了用于结构非线性分析的同伦摄动法，算例结果表明该方法可以有效地减少钢框架结构非线性分析的迭代步数和计算时间。张令心等[30]提出了改进的直角N-R法，采用切线刚度矩阵和法向刚度矩阵混合的法来增加位移的增量步长，在不增加额外计算量的情况下，提高了N-R法的迭代效率，但是该方法的稳定性和通用性还需要进一步研究和验证。
本文对基于Woodbury非线性方法的迭代求解进行研究，根据该方法的迭代特点，采用N-R法、修正牛顿法、3阶两点法、4阶两点法、三点法五种算法分别求解平衡方程，并通过理论和算例对比分析各迭代算法的计算性能。以隔离非线性有限元法为例，介绍了5种迭代方法求解其平衡方程的迭代过程，并利用时间复杂度理论建立了5种迭代算法求解该方法平衡方程的时间复杂度分析模型，对迭代方法的效率进行量化对比。通过2个采用力控制加载的结构非线性分析算例，对比了5种迭代算法求解基于Woodbury非线性方法平衡方程的效率和精度，分析了各算法适用的非线性问题。最后，对5种算法的综合性能进行评价，为基于Woodbury非线性方法的迭代算法选用提供参考建议。
1 基于Woodbury非线性方法及迭代求解
1.1 Woodbury公式在结构非线性分析中的应用
非线性分析中，刚度矩阵随位移变化而改变，此时结构的力与变形关系是非线性的。对于一个自由度数为n的结构，整体的非线性平衡方程可表达为[31]：
式中：X(n×1阶)为位移响应；F(n×1阶)为外荷载向量；
为结构节点恢复力向量。
通常，可用增量法求解式(1)，将作用于结构的荷载分为若干增量，如图1所示。每次施加一个荷载增量并求解出其对应的位移响应，再逐步施加荷载直至达到相应的荷载水平。在某个荷载增量步中，刚度矩阵取为上一级荷载下的结构刚度矩阵，且在该增量步中保持不变，即将非线性方程式(1)进行线性化，则第k+1步的增量方程为：
式中：
和ΔRk分别为荷载增量和恢复力增量。通过线性化处理，恢复力增量可近似表达为：
图1 非线性方程增量求解示意图
Fig.1 Incremental solution of the nonlinear equation
将式(3)代入式(2)，整理可得：
式中，Kt(n×n阶)为第k个增量步的切线刚度矩阵。
在传统非线性分析中，通常采用直接法求解式(4)，得到每个增量步的位移向量。为了提高求解效率，此时可将切线刚度矩阵拆分为弹性刚度矩阵Ke(n×n阶)与产生非线性相关的修正矩阵ΔK(n×n阶)的和，即：
局部非线性问题中，结构进入非线性的范围很小，因此ΔK是一个稀疏矩阵，通常可写成2个低秩矩阵的乘积，即
均为n×p阶矩阵。则用Woodbury公式[5]求解切线刚度矩阵的逆为：
式中：
称为舒尔补矩阵。可以看出，Woodbury公式将n×n阶矩阵的求逆运算转换成p×p阶舒尔补矩阵的求逆以及矩阵和向量之间的运算，当p≪n时，该公式可显著提高矩阵求逆的效率。
基于Woodbury的非线性分析方法的核心思想是用式(6)求解式(4)中每个增量步的结构位移响应，但每种方法的区别在于矩阵ΔK的构造方式不同。本文以隔离非线性有限元方法[20]为例介绍该类方法的求解过程。隔离非线性有限元方法的基本思想是将材料应变分解为弹性应变与塑性应变，图2为建立其控制方程的流程图。通过材料应变分解引入的塑性变形成为附加未知量，借鉴传统有限元中建立单元位移插值函数(矩阵N)的思想，在单元内设置塑性插值点得到塑性应变插值函数(矩阵C)来反映塑性变形在单元中的分布形式，此时结构的非线性行为可用塑性变形及与其相关的局部塑性矩阵表示。根据虚功原理和塑性插值点处的内力平衡条件可建立单元的控制方程，再对单元集成得到每个增量步中结构整体的控制方程[6]：
式中：
为与塑性变形相关的矩阵，两者均是稀疏矩阵，d为塑性插值点处产生塑性的自由度数；ΔX(n×1阶)为结构节点位移增量向量；ΔF(n×1阶)为结构外荷载增量向量；
(d×1阶)为塑性变形增量向量。
图2 建立隔离非线性有限元法控制方程的流程图
Fig.2 Flow chart of establishing the governing equation of the Inelasticity-Separated Finite Element Method
对式(7)整理，消去
可得：
式(8)的系数矩阵与结构切线刚度矩阵Kt是等价的，即
矩阵
为弹性刚度矩阵Ke的修正矩阵ΔK，因此用Woodbury公式[6]求解位移响应为：
通常采用从右至左的顺序[23]求解式(9)，主要计算为矩阵
的分解及矩阵和向量的相乘运算。式(9)的效率分析表明[23]，当结构产生局部非线性时，其渐进时间复杂度为O(mn)(m为刚度矩阵的带宽)，而传统LDLT分解法的渐进时间复杂度为O(m2n)，因此用Woodbury公式求解结构响应具有高效性。
1.2 非线性方程的迭代求解
通过1.1节的理论可知，隔离非线性有限元方法应用Woodbury公式求解增量步的控制方程，避免了传统方法中对切线刚度矩阵的重复合成和分解，提高了每个增量步的计算效率。但是，在结构非线性分析中，由于引入了式(3)的线性近似关系，得到的结果可能有较大的误差，通常在每个增量步中需要迭代计算，使位移结果满足精度要求。因此，迭代算法的计算性能也会决定非线性分析的效率。本节选用不同的迭代算法求解隔离非线性有限元法的平衡方程，研究其迭代性能。选用的原则是算法要符合该方法的求解特点，即迭代算法应基于线性化的思想并且对切线刚度矩阵进行分解，否则不适用于求解其非线性平衡方程。例如：非线性共轭梯度法[32]不是线性化方法以及拟牛顿法[33]直接得到柔度矩阵而无需对切线刚度矩阵分解，这2种方法无法求解隔离非线性有限元法的平衡方程。该文选用N-R法、修正牛顿法、3阶两点法、4阶两点法和三点法5种非线性迭代算法，介绍了各算法求解隔离非线性有限元法平衡方程的迭代过程，其中前两种方法为目前非线性分析中常用的迭代算法，其余三种方法是数学中的迭代算法，因此将上述5种方法分为两类介绍。
1)N-R法和修正牛顿法
用迭代方法求解结构非线性平衡方程式(1)的本质是在每个增量步中消除线性化近似模型和非线性模型之间不平衡力的过程，区别在于不同方法的迭代策略不同。其中，N-R法是非线性分析中常用的迭代算法，具有2阶收敛性。该方法迭代方程为[34]：
式中：下标i表示迭代步，为了表达方便，省略了增量步上标k；不平衡力向量为
在隔离非线性有限元法中，切线刚度矩阵
可用式(9)中的Woodbury公式求解位移增量。由1.1节的分析可知，对于局部非线性问题，式(9)的计算量随结构自由度数n增多而增加[23]，因此若每个迭代步都重新求解式(9)，对于大规模结构会显著降低非线性分析的效率。
修正牛顿法可改善N-R法的上述问题，避免了每个迭代步均更新和分解切线刚度矩阵，而更新频率可根据结构不同的非线性程度确定，即非线性较强更新频率高，非线性弱更新频率低，通常在每个增量步的第一个迭代步更新和分解切线刚度矩阵即可。该方法的迭代方程为[34]：
式中，
为一个增量步的初始刚度矩阵。
此时，用Woodbury公式求解式(11)的线性方程组时仅需要矩阵和向量的相乘运算及回代计算，降低了每个迭代步的计算量，但该方法的收敛性为线性收敛。
2)两点法和三点法
两点法和三点法属于数学中的多点迭代方法，求解过程中无需计算高阶导数。与N-R法相同，2种方法均为求解非线性方程的线性化方法，但在每个迭代步中需多次求解线性方程组。其中，第一次求解是对迭代初值的预测，其余是对预测值的校正，目的是通过预测校正的过程更好逼近真实解。两点法包括3阶两点法和4阶两点法2种形式，求解式(1)的迭代格式可统一写为[35]：
式中：f(X)=F-R(X)，Xj 为一个迭代步的中间变量。若式(12)中2个线性方程组的雅可比矩阵相等，即f′(Xj)=f′(Xi)，则该迭代格式具有3阶收敛性，称为3阶两点法；若2个线性方程组的雅可比矩阵不相等，则该迭代格式具有4阶收敛性，称为4阶两点法，但计算2次雅可比矩阵会显著增加计算量。为此，Chun[36]引入一个关于ti的多项式实函数g(ti)来用第一个方程组的f′(Xi)近似求解f′(Xj)，此时f′(Xj)≈f′(Xi)/g(ti)，该迭代格式仍具有4阶收敛性。
将式(1)代入式(12)式中：整理可得两点法的迭代方程分别为：
3阶两点法：
4阶两点法：
则该迭代步的总位移为Xi+1=Xi+ΔXj+ΔXj0。其中，ΔXj=Xj-Xi ，ΔXj0=Xi+1-Xj 均表示位移增量向量；Δfi=F-Ri、Δfj=F-Rj均表示每次求解的不平衡力增量向量；由式(3)可知，
即结构切线刚度矩阵为节点力对位移的导数。多项式函数g(ti)取为：
图3 两点法迭代计算示意图
Fig.3 Iterative calculation of the two-point method
式(13)和式(14)的计算流程如图3所示。隔离非线性有限元方法中，切线刚度矩阵为[Kt]i=Ke-
可用Woodbury求解式(13)和式(14)的线性方程组，由于迭代格式中仅需更新一次切线刚度矩阵，因此求解中只分解第一个方程组的舒尔补矩阵，第二个方程组进行相应的回代计算即可。
Džunić等[37]提出的三点法与4阶两点法具有相似的迭代格式，但三点法中每个迭代步需求解三次线性方程组，具有8阶收敛性。其迭代格式为：
式中：Xj和Xj0是该迭代步的中间变量；q(ti,si)为关于si、ti的多项式实函数[38]。
同样地，将式(1)代入式(16)中，整理可得到三点法的迭代方程为：
则该迭代步的总位移Xi+1=Xi+ΔXj+ΔXj0+ΔXj1。式(17)中，ΔXj0=Xj0-Xj、ΔXj1=Xi+1-Xj0表示位移增量向量；Δfj=F-Rj 、Δfj0=F-Rj0表示不平衡力向量。多项式函数g(ti)取值同式(15)，q(ti,si)取为：
式(17)的计算流程如图4所示。同样地，用Woodbury公式计算式(17)时，只需一次舒尔补矩阵分解运算，其余2个方程组仅需相应的回代计算。
图4 三点法迭代计算示意图
Fig.4 Iterative calculation of the three-point method
2 迭代算法的求解性能分析
2.1 迭代效率的影响因素
迭代算法效率的影响因素主要有以下3个方面[39]：① 算法的收敛速度；② 算法一个迭代步的计算量；③ 算法的总计算量。显然，仅从某个方面衡量迭代算法效率是不全面的。本文从上述3个方面对比分析N-R法、修正牛顿法、3阶两点法、4阶两点法以及三点法的计算效率。其中，5种方法的收敛速度和一个迭代步的运算对比结果如表1所示。
表1 5种方法的收敛速度和一个迭代步的运算对比
Table 1 Convergence rate and operations in each iteration of the five algorithms
对比项 方法N-R法修正牛顿法 3阶两点法 4阶两点法三点法收敛阶数 2 1 3 4 8矩阵分解次数1 0 1 1 1矩阵回代次数1 1 2 2 3
由表1可知，算法的收敛速度由快到慢依次为：三点法、4阶两点法、3阶两点法、N-R法、修正牛顿法。通过对比各迭代算法一个迭代步的运算可定性判断其计算量高低。其中，三点法最高；4阶两点法和3阶两点法次之，虽然2种方法的矩阵分解及回代次数相同，但4阶两点法需额外计算一个多项式函数，因此其计算量比3阶两点法略高；其次是N-R法；修正牛顿法最低。
迭代算法的总计算量等于总迭代步数乘以一个迭代步的计算量。由上述分析可知，收敛速度快的算法通常一个迭代步的计算量高，收敛速度慢的算法一个迭代步的计算量低。因此，迭代算法的总计算量无法定性判断，需要根据迭代格式和每个迭代步中线性方程组的规模及求解方式进行量化分析。本文利用时间复杂度理论建立各迭代算法时间复杂度模型来量化分析5种迭代算法求解隔离非线性有限元方法平衡方程的计算效率。
2.2 迭代算法时间复杂度模型及对比分析
算法时间复杂度理论是常用的算法效率定量评价方法，该方法的分析结果不依赖于计算机性能等因素的影响，可客观判断算法效率。其基本思想是将算法执行步骤拆分成加法、减法、乘法及除法四则基本运算，分别统计每种运算的执行次数，再乘以各自的执行时间，即可得到算法的时间复杂度函数。假设计算机执行四则基本运算(加减乘除)的时间相同，则时间复杂度函数的计算公式为[40]：
式中：T(a)为算法的时间复杂度函数；a表示问题的规模；t为四则运算的执行时间(一般取单位时间)；ei(a)(i=1,2,3,4)分别为加减乘除的操作次数。
根据上述的时间复杂度理论，分别建立用N-R法(式(10))、修正牛顿法(式(11))、3阶两点法(式(13))、4阶两点法(式(14))和三点法(式(17))求解隔离非线性有限元平衡方程的时间复杂度分析模型，结果如表2所示。其中，iN、iMN、iL3、iL4及iS分别为N-R法、修正牛顿法、3阶两点法、4阶两点法及三点法的迭代步数；m为弹性刚度矩阵Ke的带宽；α和β分别为稀疏矩阵K′和K′每列最大非零元素个数，和单元类型相关；γ为修正牛顿法中增量步数与总迭代步数的比值，其倒数表示每个增量步中平均迭代步数，反映了刚度矩阵更新的频率。
比较其余4种迭代算法和N-R法的时间复杂度模型可知，若某个迭代算法的时间复杂度小于N-R法，则有：
式中，Tj (j=1,2,3,4)分别代表其余4种迭代算法的时间复杂度TS、TL4、TL3、TMN。
表2 5种迭代算法的时间复杂度分析模型
Table 2 Time complexity analysis models of the five iterative algorithms
迭代算法 时间复杂度模型Ti+d+αd+βd+mn+dm■--■■2 3 d N-R法2=344113 3 2 3 NN ■■■Ti+d+d+βd+d+22 MNMN=γd 3γ■■■3 4γ 2γ24α修正牛顿法4mndm+---■■331 d2■Ti+dαd+βdmn+ndm+■3■+--2■d 3阶两点法2=5+882222 3 2 3 L3L3■■■2 Ti+dαd+βdmn+ndm+■3■=5+8862221 3 2 3■+--d 4阶两点法2 L4L4■■■3■■三点法Ti+dαd+βdmn+ndm+■--=7+1212123131 3 2 3 d 2 2 SS+■■■
虽然3阶两点法、4阶两点法以及三点法的理论收敛速度高于N-R法，但是一个迭代步的计算量也高于N-R法。因此与N-R法相比，3种迭代算法的迭代步数降低比例必须大于一定值(该值为至少降低的比例值)才能使总的时间复杂度低于N-R法(即满足式(20))，因此本文将这个比例值的下界pmin作为量化指标，令Tj=TN并代入表2中三种算法的总时间复杂度，可得：
式中：pmin-L3、pmin-L4、pmin-S分别为3阶两点法、4阶两点法和三点法迭代步数降低比例最小值；
分别为N-R法、3阶两点法、4阶两点法以及三点法每个迭代步的时间复杂度，即表2中括号中的部分。
显然，迭代步数降低比例最小值pmin越小，3种算法的总时间复杂度越易低于N-R法，算法越易保证效率优势。由式(21)可以看出，pmin与结构自由度数n和产生塑性自由度数d相关。为了直观地表示该值的变化趋势，图5给出了在不同塑性自由度比l=d/n的情况下，迭代步数降低比例的最小值pmin-L3、pmin-L4、pmin-S随自由度数n的变化曲线。其中，取结构刚度矩阵的半带宽为
从图5可以看出，结构自由度数n取值为102～106，除了n较小时，3种方法(三点法、3阶两点法、4阶两点法)的迭代步数降低比例最小值pmin的量级很低，在10-2～10-4。由于三点法每个迭代步的时间复杂度最高，该方法迭代步数降低比例最小值pmin-S最大。而3阶两点法和4阶两点法每个迭代步的时间复杂度相差很小，因此pmin-L3与pmin-L4的曲线无显著差别。3种方法的pmin均随自由度数n的增大而降低，与N-R法相比，3种方法效率的变化趋势相同，即n越大，3种算法效率优势越大，计算性能越好；若n较小，则3种算法不易发挥效率优势。当结构自由度数n相同时，随着塑性自由度比l增大，3种方法的pmin均逐渐减小，3种算法效率优势越显著。
表3给出了几组结构参数下3种方法的pmin取值。可以看出，当自由度数n＜103时，3种方法的pmin值较大，最大值达到0.68，即该方法的迭代步数至少比N-R法降低68%才能保证其效率高于N-R法；当n＞103时，三种方法的pmin显著降低，特别是塑性自由度比例l较大时，最低量级达到10-4～10-5，即3种方法的迭代步数只需比N-R法降低10-2%～10-3%，则效率高于N-R法，效率优势显著。因此，一般情况下，当n＜103且l较小，pmin较大，此时N-R法的效率更高；当n＞103或l较大，pmin显著降低，此时3种算法的效率更高。
图5 3种迭代算法pmin随n的变化图
Fig.5 Relationship between pmin and n
表3 3种迭代算法pmin值
Table 3 pmin of the three iterative algorithms
参数pmin迭代方法l=0.05 n=102 n=103 n=104 n=106 3阶两点法 0.50 0.46 0.14 3.11×10-4 4阶两点法 0.52 0.47 0.15 3.12×10-4三点法 0.68 0.64 0.25 6.24×10-4参数pmin迭代方法l=0.10 n=102 n=103 n=104 n=106 3阶两点法 0.49 0.31 0.028 8.40×10-5 4阶两点法 0.50 0.32 0.028 8.40×10-5三点法 0.67 0.48 0.055 1.68×10-4
与两点法和三点法不同，修正牛顿法的理论收敛速度低于N-R法，而每个迭代步的时间复杂度也低于N-R法，即该方法降低了N-R法每个迭代步的时间复杂度，但增加了迭代步数。若要使修正牛顿法的时间复杂度满足式(20)，则该方法迭代步数增加的比例值必须要小于一定值(该值为至多增加的比例值)，同样令TMN=TN 可得到该比例值的上界pmax-MN：
式中：
表示修正牛顿法每个迭代步的时间复杂度。pmax越大，与N-R法相比，修正牛顿法越易保证效率优势。图6给出了不同γ情况下，迭代步数增加比例的最大值pmax随自由度数n的变化曲线。
图6中，自由度数n的取值为102~106，γ分别取0.4和0.6，即修正牛顿法中每个增量步的平均迭代步数分别为2.5和1.67。从图6(a)~图6(b)可以看出，当结构产生的塑性自由度比例l相同时，pmax随n的增大而增大，即结构自由度数n越大，修正牛顿法的效率优势越大；当n相同时，pmax随l的增大而增大，即塑性自由度比例l越大，修正牛顿法的效率优势越大。此外，对比图6(a)、图6(b)，γ越小，即矩阵更新频率越低，pmax越大，修正牛顿法越易有效率优势。
表4列出了在不同结构参数和γ值时的pmax。可以看出，当n≤103且l较小时，pmax值很小，最小量级仅有10-3，说明此时只要修正牛顿法的迭代步数超过N-R法，其效率将低于N-R法；当n>103或l较大时，pmax逐渐增大，且当n=106，此时在不同l情况下，pmax均趋于相同值，分别趋于1.5和0.67，说明只要修正牛顿法增加的迭代步数不超过N-R法的1.5倍和0.67倍，其效率高于N-R法。因此，通常情况下，当n＜103且l较小，pmax较小，此时N-R法的效率更高；当n＞103或l较大，pmax显著提高，此时修正牛顿法的效率更高，且当n达到一定值时，修正牛顿法的效率优势会趋于稳定。
图6 修正牛顿法的pmax随n变化图
Fig.6 Relationship between pmax and n
表4 修正牛顿法的pmax值
Table 4 pmax of the modified Newton method
pmax参数l=0.05 n=102 n=103 n=104 n=106 γ= 5.4×10-3 9.56×10-2 0.9945 1.4988 0.4 γ= 3.6×10-3 6.17×10-2 0.4979 0.6663 0.6 pmax参数l=0.10 n=102 n=103 n=104 n=106 γ= 0.0307 0.4969 1.3952 1.4997 0.4 γ= 0.0203 0.2842 0.6349 0.6666 0.6
2.3 迭代算法性能的综合评价指标
当迭代算法可以得到收敛解时，用2.2节中建立的时间复杂度模型对比各迭代算法的效率才是有意义的。结构非线性分析过程中，对于不同的非线性问题或不同的计算条件，迭代算法有可能得不到收敛解。因此，对迭代算法求解性能评价，还应根据算法收敛性、计算误差及效率等方面综合判断。本文的评价方法为：首先，对算例中迭代算法求解性能进行排序，排序的原则是当各算法的计算误差在规定范围(通常为5%)内时，依据总时间复杂度由低到高排序，若算法没有得到收敛解，则不参与排序。其次，根据每个算例中算法排序结果，计算算法性能的综合评价指标[41]：
式中：h为评价的算法个数，本文中h=5；p为计算工况数；Qij表示各迭代算法的排序次数，i=1～5分别代表N-R法、修正牛顿法、3阶两点法、4阶两点法以及三点法；j=1～5代表算法排序。如Q11表示N-R法排序第1的次数。
从式(23)可以看出，Gi越大，算法综合性能越好。若Q11=p，表明N-R法在所有计算工况下性能排序均是第一，此时G1=100，综合性能最好。然而，式(23)并没有计算算法得不到收敛解的情况，此时可理解为若算法在某工况下不收敛，则此工况对算法综合指标的贡献为零。
3 算例分析
本节通过2个算例，共4种荷载工况，对5种迭代算法进行精度效率对比和性能综合评价。求解精度评价指标为误差，将N-R法的解看作精确解。效率评价指标根据迭代效率影响因素得到，其中，收敛速度用总迭代步数或平均迭代步数衡量；一个迭代步的计算量用平均时间复杂度(等于所有计算步的总时间复杂度与计算步数的比值)衡量；总的计算量用总时间复杂度(等于平均时间复杂度乘以总迭代步数)衡量。最后，根据每种工况下的算法排序，用式(23)计算各算法的综合性能指标。
3.1 空间钢框架非线性分析
5层空间钢框架，层高3.6 m，跨度6 m，构件截面尺寸如图7所示。有限元模型中每个构件划分3个单元，共1815个单元，1480个节点，8610个自由度。首层质量为81×104 kg，2层～4层质量为72×104 kg，5层质量为63×104 kg。钢材本构关系采用Clough随动硬化模型，弹性模量为Ee=2.0×105 MPa，屈服强度为fy=400 MPa，屈服后刚度系数为0.05。对该结构分别进行非线性动力时程分析及静力非线性分析，考虑结构的P-Δ效应。两种荷载工况分别为：① 施加x向、y向地震动，x方向的峰值加速度为3.5 g，y方向为3.15 g。步长为0.01 s，计算2000步。② 施加倒三角单调荷载(图7所示)，采用力加载方案，力增量为ΔP=400 N，计算2000步。
结构产生的塑性自由度数目曲线和顶点位移曲线分别如图8、图9所示，迭代算法的各项指标对比结果如表5所示。工况1是对该框架结构进行动力时程分析，对比结果为：
1)精度方面：在此工况下，N-R法、4阶两点法以及三点法得到收敛解，从图9(a)可以看出，3种方法顶点位移时程曲线吻合很好，与N-R法相比，其余2种方法的误差很小，分别为0.089%和0.11%。由于输入的地震动强度大，结构产生的塑性自由度达到总自由度数的9%左右，且地震荷载属于随机动力荷载，其作用形式复杂，无规律，易引起结构产生较复杂的非线性行为，因此切线刚度矩阵更新频率低的修正牛顿法和3阶两点法计算发散。
2)效率方面：从总迭代步数或每个增量步的平均迭代步数可以看出，三点法收敛速度最快、其次是4阶两点法，N-R法最慢；N-R法的平均时间复杂度最少，三点法最高，但两者相差少于5%；由于4阶两点法和三点法的迭代步数比N-R法显著降低，实际降低比例分别为45.22%和59.65%，远超过保证其效率优势的迭代步数降低比例最小值pmin(4.74%)和(9.06%)，因此4阶两点法和三点法的总时间复杂度比N-R法显著降低，分别降低约43.97%和55.78%。
综上，动力荷载工况中，三点法的迭代步数显著少于其他方法，弥补了每个迭代步时间复杂度高的缺点，在保证精度的条件下效率最高，因此三点法的计算性能最好。
图7 5层空间钢框架
Fig.7 5-story space steel frame structure
图8 塑性自由度数随加载步变化曲线
Fig.8 The curve of the inelastic degrees of freedom
图9 顶点位移随加载步变化曲线
Fig.9 The curve of the top displacement
表5 迭代算法对比
Table 5 Comparison of the iterative algorithms
注：表中的降低百分比均是与N-R法相比的计算结果，计算公式为｜λi-λN｜/λN,λi和λN分别为某种算法和N-R法的所对比的指标值。
方法对比指标迭代步数 平均迭代步数 实际降低(↓)/提高(↑)比例/(%)理论pmin或pmax/(%)平均时间复杂度总时间复杂度 顶点位移误差/(%)排序N-R法 6612 3.306 0 0 8 3.462×10 2.290×1012 0 3工况14阶两点法 3622 1.811 45.22↓ pmin=4.74 8 3.543×10 1.283×1012 0.0891 2三点法 2668 1.334 59.65↓ pmin=9.06 8 3.624×10 9.668×1011 0.11 1 N-R法 4255 2.123 0 0 8 9.734×10 4.142×1012 0 5修正牛顿法 4244 2.122 0.26↓ pmax=103 8 4.634×10 12 1.967×10 5 7.53×10- 1工况23阶两点法 2756 1.378 35.23↓ pmin=3.85 8 9.8275×1012 2.708×10 4 1.648×10- 4 4阶两点法 2581 1.29 39.3↓ pmin=3.87 8 9.8279×1012 2.537×10 5 7.96×10- 3三点法 2020 1.01 52.53↓ pmin=7.46 8 9.922×10 12 2.00×10 4 1.668×10- 2
工况2是对该结构施加静力单调荷载，产生的塑性自由度数平均值为823，约占总位移自由度数的9.56%。对比结果为：
1)精度方面：5种迭代算法计算结果均收敛，且位移曲线与N-R方法吻合很好(图9(b))，误差量级非常小，均在10-4和10-5之间，证明4种迭代算法的精度与N-R法相当。
2)效率方面：迭代步数最少的是三点法，其次是4阶两点法和3阶两点法，3种方法均少于N-R法，这与理论收敛速度相吻合。而修正牛顿法的迭代步数也比N-R法略小，此时收敛速度和N-R法相当。平均时间复杂度最低的是修正牛顿法，最高的是三点法；由于3阶两点法、4阶两点法和三点法的迭代步数比N-R法显著降低，均远远超过相应的降低比例最小值pmin，因此总时间复杂度均低于N-R法，分别降低约34.6%、38.7%以及51.7%；修正牛顿法的总时间复杂度也比N-R法降低约52.5%。
综上，在静力单调荷载工况下，修正牛顿法的迭代步数没有显著增加，在保证精度的条件下其效率最高，因此该方法的计算性能最好。
3.2 平面薄钢板静力非线性分析
承受均布荷载的薄钢板，尺寸和形状如图10所示，板厚0.1 m。有限元模型采用三节点平面单元，共离散成4538个单元，2390个结点，4698个位移自由度。钢板的本构模型采用Clough随动强化弹塑性模型，弹性模量为Ee=2.0×105 MPa，屈服强度为fy=400 MPa，屈服后刚度系数为0.02，泊松比为0.2。对该薄板进行静力非线性分析，共2种荷载工况：① 施加静力单调荷载，加载模式为变步长加载，计算2000步。在钢板的弹性阶段，前5个计算步，力增量为ΔP=90 kN/m，其余计算步的力增量为ΔP=0.78 kN/m。② 正弦加载，表达式为P(t)=Pmax sin(t)，Pmax=2000 kN/m，计算2000个荷载步，
图10 平面薄钢板
Fig.10 Plain steel plate
图11和图12分别为产生塑性自由度数目变化曲线和薄板边缘点(图10)的位移随加载步的变化曲线(位移方向为负方向)，各迭代算法的对比指标如表6所示。工况1是对平面薄钢板结构进行静力单调加载，产生塑性自由度平均值约为515，约占总位移自由度数的11%。对比结果为：
1)精度方面：除修正牛顿法外，其余迭代方法均得到收敛解，且各算法位移曲线与N-R法吻合很好，误差量级仅为10-4。产生这种现象的原因为：修正牛顿法的切线刚度矩阵更新频率最低，该方法的收敛性易受结构非线性程度和加载步长等因素影响，收敛性较差。
2)效率方面：三点法的收敛速度最快，N-R法最慢。3阶两点法、4阶两点法及三点法的迭代步数无明显差别，均比N-R法降低约50%，远超过相应的降低比例最小值pmin为5.97%、6.01%和11.35%，因此3种算法总时间复杂度均比N-R法低一个量级，计算效率显著高于N-R法。
图11 塑性自由度数曲线
Fig.11 The curve of the inelastic degrees of freedom
图12 位移随加载步变化曲线
Fig.12 The curve of the displacement
表6 迭代算法对比
Table 6 Comparison of the iterative algorithms
方法对比指标迭代步数 平均迭代步数 实际降低(↓)/提高(↑)比例/(%)理论min/(%)p平均时间复杂度总时间复杂度 顶点位移误差/(%)排序N-R法 4032 2.016 0 0 8 4.5736×10 1.844×1012 0 4 3阶两点法 2008 1 50.2↓ 5.97 8 4.6147×10 11 9.266×10 4工况12.83×10- 2 4阶两点法 2007 1 50.2↓ 6.01 8 4.6149×10 11 9.262×10 4 4.90×10- 1三点法 2004 1 50.3↓ 11.35 8 4.6563×10 11 9.331×10 4 1.13×10- 3 N-R法 4221 2.11 0 0 8 7.683×10 3.243×1012 0 3工况24阶两点法 2186 1.09 48.2↓ 4.2 8 7.731×10 1.690×1012 0.0174 2三点法 2105 1.05 50.1↓ 8.06 8 7.780×10 12 1.638×10 3 9.17×10- 1
综上，三点法的迭代步数与4阶两点法和3阶两点法相差甚微，无明显优势，此时三点法的总时间复杂度高于2种方法。因此，在此工况下，4阶两点法在保证精度的条件下效率最高，该方法的计算性能最好。
工况2是对平面薄钢板结构进行静力正弦往复加载，产生塑性自由度平均值约为601，约占总位移自由度数的12.8%。对比结果为：
1)精度方面：N-R法、4阶两点法和三点法得到收敛解，且三种方法位移曲线吻合较好。与N-R法相比，4阶两点法和三点法的误差量级仅有10-3～10-2。由于N-R法、4阶两点法以及三点法对切线刚度矩阵实时更新，与结构实际的刚度矩阵更接近，因此在往复加载的复杂条件下收敛性较好。而3阶两点法和修正牛顿法中切线刚度矩阵更新频率较其他3种方法低，在此工况下收敛性差，无法得到收敛解。
2)效率方面：三点法的迭代步数最少，4阶两点法次之，均少于N-R法，即三点法收敛速度最快，N-R法最慢；4阶两点法和三点法的迭代步数比N-R法分别降低约48.2%和50.1%，均远超过其理论降低比例最小值pmin，因此总时间复杂度远低于N-R法，分别降低约47.9%和49.5%。
综上，正弦往复加载条件下，三点法在保证精度的条件下效率最高，因此该方法的计算性能最好。
3.3 迭代算法性能综合评价
通过对第3.1节～第3.2节2个算例中5种迭代算法性能的对比分析，可以看出：在结构非线性分析中，该文选用的5种迭代算法并没有绝对的最优算法，都有各自的适用性和优势。为了综合评价各算法的计算性能，考虑迭代算法是否收敛、收敛速度、总时间复杂度以及误差等因素，根据式(23)可得到5种迭代算法的综合评价指标Gi，如表7所示。
表7 5种迭代算法综合评价指标Gi
Table 7 Comprehensive performance indexes Gi of the five iterative algorithms
注：j=0代表没有得到收敛解的次数。
Qij方法j=0 j=1 j=2 j=3 j=4 j=5Gi N-R法(i=1)0 0 0 2 1 1 45修正牛顿法(i=2)3 1 0 0 0 0 25 3阶两点法(i=3)2 0 1 0 1 0 30 4阶两点法(i=4)0 1 2 1 0 0 80三点法(i=5)0 2 1 1 0 0 85
从表7可以看出，N-R法、4阶两点法及三点法在所有荷载工况下均得到了收敛解，收敛性较好3阶两点法有2次无收敛解，即动力加载及静力往复加载条件2种工况，收敛性较差；而修正牛顿法有3次无收敛解，分别为动力、静力往复和单调加载条件3种工况，该方法只有在简单加载条件且步长较小时才能得到收敛解，收敛性最差。由式(23)可知，无收敛解的情况对综合指标Gi的贡献为0，因此无收敛解的工况越多，Gi越小。在得到收敛解的情况下，算法排序高的情况越多，Gi越大。综合指标Gi的结果表明：三点法的综合性能最好，4阶两点法次之，N-R法第3，3阶两点法和修正牛顿法的性能较差，分别位于第4和第5。
4 结论
本文用N-R法、修正牛顿法、3阶两点法、4阶两点法及三点法5种迭代算法求解基于Woodbury非线性方法的平衡方程，利用算法时间复杂度理论建立了5种方法的时间复杂度分析模型。通过对比其余4种方法和N-R法的时间复杂度，得到3阶两点法、4阶两点法和三点法的迭代步数降低比例的最小值pmin以及修正牛顿法迭代步数增加比例最大值pmax，并将其作为量化指标对比了4种迭代算法与N-R法的效率。最后，通过数值算例对各迭代算法的收敛性、收敛速度、时间复杂度以及误差等进行分析，并综合对比了各迭代算法的计算性能。得到结论如下：
(1)选用的迭代算法可有效地求解基于Woodbury非线性方法的平衡方程，在得到收敛解的情况下，其余4种方法的计算结果与N-R法之间的误差量级很小，满足精度要求。其中，N-R法、4阶两点法及三点法的收敛性较好；3阶两点法的收敛性较差，修正牛顿法的收敛性最差。
(2)通过理论和算例分析，对基于Woodbury非线性方法平衡方程的迭代算法选用给出如下建议：一般情况下，对于小规模结构(自由度数n＜1000)且产生塑性自由度数较少(塑性自由度比例l≤5%)，选用N-R法求解的效率最高；对于中大型规模结构(n＞1000)或产生塑性自由度数较多(l＞5%)，在动力或静力往复等复杂的加载条件下，选用三点法求解的效率最高，而在静力单调的简单加载条件且步长较小时，选用修正牛顿法求解的效率最高。
(3)通过计算5种迭代算法求解基于Woodbury非线性方法平衡方程的综合性能指标，可得算法综合性能由高到低排序为：三点法、4阶两点法、N-R法、3阶两点法和修正牛顿法。因此，在通常情况下，三点法和4阶两点法的计算性能优于N-R法，而3阶两点法和修正牛顿法的计算性能不如N-R法。
[1]孙宝印, 古泉, 张沛洲, 等.钢筋混凝土框架结构弹塑性数值子结构分析方法 [J].工程力学, 2016, 33(5)：44－49.Sun Baoyin, Gu Quan, Zhang Peizhou, et al.Elastoplastic numerical substructures method of reinforced concrete frame structures [J].Engineering Mechanics, 2016, 33(5)： 44－49.(in Chinese)
[2]Kim D J, Duarte C A, Proenca S P.A generalized finite element method with global-local enrichment functions for confined plasticity problems [J].Computational Mechanics, 2012, 50(5)： 563－578.
[3]Lin T H.Theory of inelastic structures [M].New York：John Wiley & Sons, 1968： 43－55.
[4]Li G, Zhang Y, Li H N.Nonlinear seismic analysis ofreinforced concrete frames using the force analogy method [J].Earthquake Engineering & Structural Dynamics, 2014, 43(14)： 2115－2134.
[5]Akgün M A, Garcelon J H, Haftka R T.Fast exact linear and non-linear structural reanalysis and the Sherman-Morrison-Woodbury formulas [J].International Journal for Numerical Methods in Engineering, 2001, 50(7)：1587－1606.
[6]Li G, Yu D.Efficient inelasticity-separated finite element method for material nonlinearity analysis [J].Journal of Engineering Mechanics, 2018, 144(4)： 04018008-1－04018008-11.
[7]高国强.车身结构设计中快速重分析方法研究[D].长沙： 湖南大学, 2015： 16－30.Gao Guoqiang.The research of fast reanalysis algorithm for structural design of automobile body [D].Changsha：Hunan University, 2015： 16－30.(in Chinese)
[8]王琥, 种浩, 高国强, 等.重分析方法研究进展及展望[J].工程力学, 2017, 34(5)： 1－16.Wang Hu, Chong Hao, Gao Guoqiang, et al.Review of advances and outlook in reanalysis methods [J].Engineering Mechanics, 2017, 34(5)： 1－16.(in Chinese)
[9]宋琦, 杨任, 陈璞.一种新的结构修改算法及其在工程设计中的应用[J].工程力学, 2016, 33(7)：1－6.Song Qi, Yang Ren, Chen Pu.A new algorithm for structural modifications and its applications [J].Engineering Mechanics, 2016, 33(7)： 1－6.(in Chinese)
[10]Huang C, Verchery G.An exact structural static reanalysis method [J].Communications in numerical methods in engineering, 1997, 13(2)： 103－112.
[11]Deng L, Ghosn M.Pseudoforce method for nonlinear analysis and reanalysis of structural systems [J].Journal of Structural Engineering, 2001, 127(5)： 570－578.
[12]Kirsch U.Reanalysis and sensitivity reanalysis by combined approximations [J].Structural and Multidisciplinary Optimization, 2009, 40(1)： 1－15.
[13]李书, 冯太华, 范绪箕.结构静力问题的重特征值灵敏度分析 [J].工程力学, 1996, 13(4)： 97－104.Li Shu, Feng Taihua, Fan Xuji.Sensitivity analysis with repeated eigenvalues in structural static problem [J].Engineering Mechanics, 1996, 13(4)： 97－104.(in Chinese)
[14]Kirsch U, Bogomolni M.Nonlinear and dynamic structural analysis using combined approximations [J].Computers & structures, 2007, 85(10)： 566－578.
[15]Kirsch U, Bogomolni M, Sheinman I.Nonlinear dynamic reanalysis of structures by combined approximations [J].Computer Methods in Applied Mechanics &Engineering, 2006, 195： 4420－4432.
[16]Hurtado J E.Reanalysis of linear and nonlinear structures using iterated shanks transformation [J].Computer Methods in Applied Mechanics & Engineering, 2002,191： 4215－4229.
[17]Huang G, Chen X, Yang Z, et al.Exact analysis and reanalysis methods for structures with nonlinear boundary conditions [J].Computers & Structures, 2018,198： 12－22.
[18]Wong K K F, Yang R.Inelastic dynamic response of structures using force analogy method [J].Journal of Engineering Mechanics, 1999, 125(10)： 1190－1199.
[19]靳永强, 李钢, 李宏男.基于拟力法的钢支撑非线性滞回行为模拟[J].工程力学, 2017, 34(10)： 139－148.Jin Yongqiang, Li Gang, Li Hongnan.Numerical simulation of steel brace hysteretic behavior based on the force analogy method [J].Engineering Mechanics, 2017,34(10)： 139－148.(in Chinese)
[20]李钢, 余丁浩, 李宏男.基于拟力法的纤维梁有限元非线性分析方法 [J].建筑结构学报, 2016, 37(9)： 61－68.Li Gang, Yu Dinghao, Li Hongnan.Nonlinear fiber beam element analysis based on force analogy method [J].Journal of Building Structures, 2016, 37(9)： 61－68.(in Chinese)
[21]Li G, Yu D, Li H N.Seismic response analysis of reinforced concrete frames using inelasticity-separated fibre beam-column model [J].Earthquake Engineering &Structural Dynamics, 2018, 47(5)：1291－1308.
[22]Yu D H, Li G, Li H N.Improved woodbury solution method for nonlinear analysis with high-rank modifications based on a sparse approximation approach[J].Journal of Engineering Mechanics, 2018, 144(11)：04018103-1－04018103-13.
[23]李钢，贾硕，余丁浩.基于算法复杂度理论的拟力法计算效率评价 [J].计算力学学报, 2018, 35(2)： 129－137.Li Gang, Jia Shuo, Yu Dinghao.The efficiency evaluation of force analogy method based on the algorithm complexity theory [J].Chinese Journal of Computational Mechanics, 2018, 35(2)： 129－137.(in Chinese)
[24]何政, 欧进萍.钢筋混凝土结构非线性分析 [M].哈尔滨： 哈尔滨工业大学出版社, 2007： 168－177.He Zheng, Ou Jinping.Nonlinear analysis of reinforced concrete structures [M].Harbin： Harbin Institute of Technology Press, 2007： 168－177.(in Chinese)
[25]Scott M H, Fenves G L.Krylov subspace accelerated newton algorithm： Application to dynamic progressive collapse simulation of frames [J].Journal of Structural Engineering, 2010, 136(5)： 473－480.
[26]Bathe K J, Cimento A P.Some practical procedures for the solution of nonlinear finite element equations [J].Computer Methods in Applied Mechanics &Engineering, 1980, 22(1)： 59－85.
[27]Kim J H, Kim Y H.A predictor-corrector method for structural nonlinear analysis [J].Computer Methods in Applied Mechanics & Engineering, 2001, 191(8-10)：959－974.
[28]Saffari H, Mansouri I, Bagheripour M H, et al.Elastoplastic analysis of steel plane frames using homotopy perturbation method [J].Journal of Constructional Steel Research, 2012, 70(2)： 350－357.
[29]Mansouri I, Saffari H.A fast hybrid algorithm for nonlinear analysis of structures [J].Asian Journal of Civil Engineering, 2014, 15(2)： 213－230.
[30]Zhang L X, Hu M Y, Wu P C.Application analysis of an improved N-R iterative method about finite element numerical calculation [J].Key Engineering Materials Trans Tech Publications, 2011, 450： 518－521.
[31]王勖成.有限单元法[M].北京： 清华大学出版社,2003： 547－555.Wang Xucheng.Finite element method [M].Beijing：Tsinghua University Press, 2003： 547－555.(in Chinese)
[32]Suarjana M.Conjugate gradient method for linear and nonlinear structural analysis on sequential and parallel computers [D].USA： Stanford University, 1994： 1－30.
[33]De Borst R, Crisfield M A, Remmers J J C, et al.Nonlinear finite element analysis of solids and structures[M].New York： John Wiley & Sons, 2012： 30－46.
[34]Bathe K J.Finite element method [M].New Jersey：Prentice Hall, 1996： 755－759.
[35]Ruiz Á A M, Argyros I K.Two-step newton methods [J].Journal of Complexity, 2014, 30(4)： 533－553.
[36]Chun C.Some fourth-order iterative methods for solving nonlinear equations [J].Applied Mathematics &Computation, 2008, 195(2)： 454－459.
[37]Džunić J, Petković M S, Petković L D.A family of optimal three-point methods for solving nonlinear equations using two parametric functions [J].Applied Mathematics & Computation, 2011, 217(19)： 7612－7619.
[38]Petkovi M S, Petkovi L D.Families of optimal multipoint methods for solving nonlinear equations： a survey [J].Applicable Analysis & Discrete Mathematics,2010, 4(1)： 1－22.
[39]Ralston A, Rabinowitz P.A first course in numerical analysis [M].Courier Corporation, 2001： 25－29.
[40]王晓东.计算机算法设计与分析[M].第3版.北京：电子工业出版社, 2007： 1－3.Wang Xiaodong.The design and analysis of computer algorithms [M].3rd ed.Beijing： Publishing House of Electronics Industry, 2007： 1－3.(in Chinese)
[41]Rezaiee-Pajand M, Estiri H.Mixing dynamic relaxation method with load factor and displacement increments [J].Computers & Structures, 2016, 168： 78－91.