基于近似Woodbury法的杆系结构静力稳定性分析

靳永强，张 涛，莫振林，袁永强，李彦滨，周雄雄

(1.中国建筑第二工程局有限公司，北京 100160；2.中国建筑西南勘察设计研究院有限公司，四川 成都 610052；3.西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

在土木工程领域中，由长度比横截面尺寸大很多的细长杆件所组成的结构称为杆系结构，如刚架、桁架、网壳等结构，其具有空间大、杆件布置灵活及结构形式优美等特点，被广泛应用于民用和工业等建筑中[1-3]。在外荷载作用下，杆系结构中的部分杆件会承受较大的轴压力，其静力稳定性分析是不可忽略的关键环节。通常，在杆系结构的设计与分析中，强度问题(代表杆件截面所能承受的最大内力)和稳定性问题(表示结构或者构件不能以原来的平衡状态承受附加荷载)是紧密联系的[4-5]。因此，为了得到杆系结构最为真实的临界荷载值及荷载与位移全过程曲线，在静力推覆分析时，需同时考虑材料非线性和几何非线性。

针对杆系结构的静力稳定性分析，已有研究取得了一系列成果，主要可分为两类[6-8]：第一类是基于构件层面的宏观分析方法，具有建模简单、计算量小等特点。例如，韦良[9]基于平衡微分方程并引入稳定函数推导了考虑几何非线性的梁柱单元刚度矩阵，同时利用塑性铰模型，提出了基于QR法的钢框架稳定性分析方法，能够准确、高效地获得结构极限承载力及变形能力；Li等[10]、靳永强等[11]基于拟力法理论建立了轴压构件宏观分析模型，通过设置两个局部塑性机制来分别模拟弯曲变形和轴向塑性伸长行为，实现轴压构件的非线性屈曲分析。第二类是基于材料层面的非线性有限元分析方法，具有计算精度良好、计算量大等特点。例如，Nguyen等[12]采用更新拉格朗日列式建立了考虑几何非线性行为的平衡方程，并利用纤维塑性铰法来模拟空间框架结构的侧扭稳定性行为，可实现塑性变形沿杆件长度上的开展，通过数值算例验证了该方法的准确性和高效性；叶康生等[13]、吴可伟[14]基于共旋坐标法和塑性区法提出了一种杆系结构弹塑性大位移分析方法，将单元端部截面分割为若干小块面积，单元内的截面一直保持弹性状态，快速实现结构的静力稳定性分析。

在杆系结构的静力推覆分析中，其数值模拟的难点在于：失稳瞬间，结构刚度发生突变，承载力骤降，荷载-位移关系曲线呈现下降趋势。为了越过荷载-位移曲线的极值点，通常将采用位移控制算法来施加推覆力[15]。为了克服结构中不同部位非线性程度的差异，满足荷载分布恒定比例关系，且提高数值求解的稳定性，黄羽立等[16-17]提出了基于多点位移控制的推覆分析方法。该方法能够较好地获取荷载-位移曲线的下降段，准确追踪失稳后的非线性平衡路径，具有良好的数值稳定性优势。

在构件层面，Wong等[18]基于位移分解原理提出了结构非线性高效分析方法-拟力法，该方法被广泛应用于工程结构的非线性分析和振动控制等领域[19]。为了实现结构的精细化高效模拟，Li等[20]将位移分解思想拓展到材料应变层面，提出了隔离非线性有限元法，其基本原理是将材料应变分解为线弹性与非弹性两部分，并与有限元相结合，形成具有刚度矩阵弹塑性分离形式的控制方程，以Woodbury公式为求解工具，实现了局部材料非线性问题的高效求解[21]。针对结构出现大范围材料非线性问题，Li等[22]通过引入组合近似法至Woodbury公式中，提出了近似Woodbury法来高效求解控制方程。随后，李钢等[23]将隔离非线性有限元法嵌入到更新拉格朗日格式中，建立了考虑双重非线性(材料和几何非线性)行为的梁单元控制方程，并对近似Woodbury法进行了改进，高效求解整体控制方程，实现了杆系结构的动力双重非线性分析，但在求解其静力稳定性问题方面还有待进一步完善。

鉴于此，本文基于多点位移控制算法和近似Woodbury法，提出了杆系结构的静力稳定性高效分析方法。该方法利用Woodbury公式对多点位移控制的平衡方程进行展开，结合隔离非线性理论形成具有刚度矩阵弹塑性分离形式的控制方程，采用近似Woodbury法对其进行高效求解，在杆系结构的静力稳定性分析中，避免整体刚度矩阵的实时分解。以Williams双杆体系[22]为数值算例，验证了本文方法在静力稳定性分析中的可行性和有效性。最后，通过某框架结构的静力稳定性分析，论证了该方法能高效地获得荷载-位移关系曲线的下降段，处理负刚度问题，反映结构的承载力、稳定性及刚度的整个变化历程。

1 多点位移控制算法概述

多点位移控制算法是一种较为新颖的静力推覆分析方法，通过引入位移刚性约束来保持恒定推覆侧力分布，以实现结构稳定性分析[16-17]。该方法在解决结构复杂非线性问题时，特别是处理负刚度问题具有显著优势，能够获取结构从弹性阶段到软化阶段的全过程。对于杆系结构，其失稳瞬间平衡路径可能会发生“跳跃”等情况，为了精确及完整地追踪非线性屈曲路径，本文采用多点位移控制进行杆系结构的静力推覆分析。

以图1所示的n自由度模型比例加载为例，对多点位移控制算法进行说明。在该模型上施加荷载(F1∶F2∶…∶Fn)满足恒定比例关系为(p1∶p2∶…∶pn=1∶2∶…∶n)，应满足如下约束方程[16-17]：

图1 多点位移控制的静力推覆分析

(1)

式中：P={p1,p2,…,pn}T，pi是第i个自由度上施加荷载的比例系数；U={u1,u2,…,un}T，ui是第i个自由度所对应的位移；u0是新增约束方程引入自由度所对应的位移。结构的平衡方程可写为：

F-Q=0

(2)

式中:F={F1,F2,…,Fn}T，是关于位移U={u1,u2,…,un}T的非线性函数，表示结构的内力；Q={Q1,Q2,…,Qn}T，表示结构所受的外荷载。式(1)是方程组(2)的约束方程，可基于拉格朗日乘子法构建新的方程组，其形式如下[16-17]：

(3)

式中:λ为拉格朗日乘子，表示力加载系数。为了求解位移U和λ，可采用牛顿迭代法求解非线性方程组(3)：

(4)

其中:i表示第i个迭代步；Kt表示第i迭代步的结构切线刚度矩阵；C0的表达式为：

(5)

2 基于近似Woodbury法的多点位移控制推覆分析

2.1 基于多点位移控制的隔离非线性方程

为了便于推导，将式(4)写为增量形式，经整理可得：

(6)

式中:ΔU和Δλ分别表示第i迭代步至第i+1迭代步的位移增量和加载系数增量；ΔR表示迭代过程中产生的不平衡力(等于外荷载向量和当前迭代步的结构内力向量之差，即ΔR=Q-F)。对上述增量方程直接求解时，需要实时分解包含切线刚度Kt和荷载比例系数向量P组成的矩阵，这将花费大量的时间成本。为此，可基于Woodbury公式，对方程(6)中矩阵求逆部分进行展开：

(7)

将式(7)代入式(6)，经整理可得：

ΔU=ΔUR+ΔλΔUP

(8)

Δλ=(C0-PTΔUR)/(PTΔUP)

(9)

式中：ΔUR和ΔUP分别表示不平衡力ΔR和荷载比例系数向量P所产生的位移增量，具体表达式分别如下：

(10)

(11)

对比式(6)和式(8)可知：在求解增量位移ΔU时，式(6)需要对包含切线刚度Kt和荷载比例系数向量P组成的矩阵进行分解，而式(8)仅需分解切线刚度矩阵Kt。根据文献[23]，可将式(10)和式(11)表达成具有隔离非线性特征的控制方程，形式分别如下：

(12)

(13)

本文采用更新拉格朗日格式描述几何非线性行为，因此，上述两个方程中的所有矩阵和向量均是以t时刻的已知构型作为参考系。由式(12)和式(13)可知，多点位移控制的平衡方程(式(6))可基于Woodbury公式进行转换，将切线刚度矩阵Kt和荷载比例系数加载向量P拆分开来，分别形成由不平衡力ΔR和荷载比例系数加载向量P作为外荷载的平衡方程(式(10)和式(11))，再结合隔离非线性的基本理论，将式(10)和式(11)表达成具有隔离非线性特征的控制方程(式(12)和式(13))。

2.2 高效求解方法

根据文献[23]可知，具备隔离非线性特征的控制方程(式(12)和式(13))可采用近似Woodbury法进行高效求解，其核心思想是：采用Woodbury公式对控制方程进行精确展开，假定实时变化的整体刚度矩阵在某一时间段内保持不变，与舒尔补矩阵相关的线性方程采用组合近似法高效求解，在杆系结构的静力稳定性分析中，能避免矩阵的实时更新和分解运算，可显著提高计算效率。在某段时间区域t∈[Tj,Tj+1]，基于近似Woodbury法求解式(12)和式(13)的具体表达式分别如下：

(14)

(15)

表1 混合近似法的时间复杂度

其中:计算次序(3)中的项目可采用组合近似法进行求解，具体求解流程可参见文献[23]。

2.3 近似矩阵更新准则

为了能够消除矩阵近似所引起的误差，可建立近似矩阵更新准则对近似矩阵进行“适时”更新，即[23]：

φ=

(16)

其中:|·|表示绝对值；φ表示由矩阵近似所产生的相对误差；Δλ(t)表示t时刻的加载系数增量；ε0为用户自定义参数，用于控制矩阵更新频次[23]。矩阵Kt(t)表示结构的切线刚度，对方程(12)进行消元，可得：

(17)

针对某时间段(t∈[Tj,Tj+1])，当式(16)成立时，则矩阵近似所引起的误差(见式(14)和式(15))可在迭代过程中进行消除。当t=Tj+1时，式(16)不成立，表明：矩阵近似所引起的误差较大，需将式(14)和式(15)中的矩阵Keg(Tj)和Kf(Tj)进行更新，即分别用矩阵Keg(Tj+1)和Kf(Tj+1)进行替换。因此，通过建立近似矩阵更新准则，可“适时”地对近似矩阵进行更新，使本文方法在满足计算精度的前提下提高计算效率。

2.4 求解流程

针对杆系结构的静力稳定性分析，本文方法的求解流程如图2所示。

图2 本文方法求解流程图

2.5 效率分析

时间复杂度是一种客观反映算法效率的工具，能够定量地评估算法本身在执行过程中所需要的理论时间，不依赖于计算机性能、编程技巧等外在因素[24-25]。就式(10)和式(12)而言，当式(10)采用传统有限元法(矩阵分解采用LDLT法)进行求解，式(12)采用近似的Woodbury法进行求解时，可依据文献[23]分别将上述两种方法的时间复杂度列于表2。

表2 三种方法的时间复杂度函数(t∈[Tj, Tj+1])

3 数值分析

3.1 模型验证

以经典的Williams双杆体系[26]为数值算例(如图3所示)，验证本文方法在解决负刚度问题时的有效性。结构杆件均为φ114 mm×4 mm的钢管，材料本构关系为理想弹塑性，屈服应力为240 MPa，弹性模量为210 GPa。每根杆件划分成10个单元(3结点Timoshenko纤维梁单元)；杆件截面被分割成24个小面积(24根纤维)，环向分成12等份，径向分成2等份。

图3 Williams 双杆体系

图4给出了本文方法与沈世钊等[26]模拟的竖向荷载-位移曲线，两者吻合较好。其中，本文方法和沈世钊等[26]所获得的临界荷载值分别为15.55 kN和 15.56 kN，表明：本文方法在进行稳定性分析时可精确获得结构的荷载-位移全过程曲线，能追踪发生“跳跃”失稳的平衡路径。

图4 Williams 双杆体系的荷载-位移曲线

3.2 数值算例

以某3层三跨钢框架结构为算例，进行静力稳定性分析，其立面和平面图如图5所示。该结构层高均为3.9 m，X方向为3跨(边跨跨度为6 m，中跨跨度为9 m)，Z方向为4跨(跨度均为6 m)，梁柱截面尺寸如图5所示。

图5 3层三跨钢框架结构

各杆件均采用3结点Timoshenko纤维梁单元，柱划分6个单元，边跨梁划分8个单元，中跨梁划分12个单元，有限元模型总计1 068个单元，总自由度为5 850，构件的截面纤维数均为19根。钢材本构关系采用双线性随动强化模型，其中弹性模量Ee=206 GPa，屈服强度σy=345 MPa，屈服后刚度Et=0.02Ee。恒荷载集中施加于梁柱节点，大小均为750 kN。采用倒三角侧向荷载模式对结构进行水平推覆(见图5)；CA法的最大基向量个数取5。

图6给出了静力推覆荷载作用下结构第①榀的基底剪力和顶点位移全过程曲线，可见采用本文方法能够较好地追踪结构失稳后的下降段曲线，即软化过程，可以获得稳定的临界荷载值及对应的侧向位移，分别为1 411 kN和0.441 m。图6也绘制了仅考虑材料非线性行为时第①榀的基底剪力和顶点位移曲线，对比结果表明：几何非线性行为导致结构失稳，在进行稳定性分析时，需同时考虑结构的材料和几何非线性(双重非线性)行为，以获得实际的极限承载力。图7给出了静力推覆分析过程中各层推覆力和顶点位移关系曲线，可以看出各层施加的推覆力严格符合等比例荷载模式，满足多点位移控制算法的荷载分布恒定比例关系。

图6 第①榀基底剪力-顶点位移曲线对比图

图7 各层推覆力荷载值-顶点位移曲线图

在静力推覆分析过程中，结构的塑性自由度数(p)时程曲线如图8所示。其中，被激活的最大塑性自由度数为1 725，占总自由度数的29.5%，每个增量步的平均塑性自由度为1 360.5(总增量步为1 788)，占总自由度数的23.3%。

图8 塑性自由度时程曲线

为了验证本文方法的高效性，图9给出了本文方法和传统有限元法(采用LDLT分解法)的时间复杂度时程曲线，其中，荷载增量步的时间复杂度等于每个迭代步所需时间复杂度的总和。由图9可知，传统有限元法在每个增量步的时间复杂度仅随迭代步数而发生变化，本文方法的时间复杂度与迭代步数、基向量个数、所激活的塑性自由度数及近似矩阵更新次数均相关，在执行自适应策略时，近似矩阵的更新会导致该增量步的时间复杂度变大。在整个静力推覆分析过程中，矩阵更新的次数为58，与总增量步数(1 788)的比例为1∶31，表明：近似矩阵更新准则(见式(16))能较好地控制矩阵更新次数，在保证迭代收敛的前提下使计算效率最优。此外，迭代步数也会影响算法的计算成本，表3统计了本文方法和传统有限元法的总迭代步数，可知：本文方法所需的迭代步数要多于传统有限元法，即本文方法的收敛性要低于传统有限元法，原因在于：近似Woodbury法中近似矩阵所引起的误差需要更多迭代步进行消除。同时，表3统计了上述两种方法的总时间复杂度和每个荷载增量步的平均时间复杂度，本文方法和传统有限元法的总时间复杂度分别为0.98×1011和5.10×1011，平均到每个增量步的时间复杂度分别为0.55×108和2.85×108，结果表明：在计算效率方面，本文方法比传统有限元法快5.18倍，验证了该方法能高效分析杆系结构的静力稳定性问题。

图9 时间复杂度时程曲线

表3 总迭代步数和时间复杂度对比

4 结 论

针对杆系结构的静力稳定性问题，本文基于多点位移控制算法和近似Woodbury公式，提出了高效的静力推覆分析方法。利用Woodbury公式对多点位移控制的平衡方程进行展开，结合隔离非线性理论形成了具有刚度矩阵弹塑性分离形式的控制方程，并采用近似Woodbury法对控制方程进行高效求解。最后，将该方法应用于某框架结构的静力稳定性分析中，并得到如下结论：

(1)本文方法是以Woodbury公式为理论基础，将多点位移控制的平衡方程变为具有隔离非线性特性的控制方程，实现刚度矩阵弹塑性分离，并采用近似Woodbury法高效求解控制方程。

(2)本文方法能够高效分析杆系结构的静力稳定性问题，可追踪失稳后的非线性平衡路径，获得下降段曲线，处理负刚度问题，能够准确获得荷载-位移全过程曲线，完整反映结构的承载力、稳定性及刚度的整个变化历程。