基于近似Woodbury法的桿系結(jié)構(gòu)靜力穩(wěn)定性分析

靳永強，張 濤，莫振林，袁永強，李彥濱，周雄雄

(1.中國建筑第二工程局有限公司，北京 100160；2.中國建筑西南勘察設(shè)計研究院有限公司，四川 成都 610052；3.西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100)

在土木工程領(lǐng)域中，由長度比橫截面尺寸大很多的細長桿件所組成的結(jié)構(gòu)稱為桿系結(jié)構(gòu)，如剛架、桁架、網(wǎng)殼等結(jié)構(gòu)，其具有空間大、桿件布置靈活及結(jié)構(gòu)形式優(yōu)美等特點，被廣泛應(yīng)用于民用和工業(yè)等建筑中[1-3]。在外荷載作用下，桿系結(jié)構(gòu)中的部分桿件會承受較大的軸壓力，其靜力穩(wěn)定性分析是不可忽略的關(guān)鍵環(huán)節(jié)。通常，在桿系結(jié)構(gòu)的設(shè)計與分析中，強度問題(代表桿件截面所能承受的最大內(nèi)力)和穩(wěn)定性問題(表示結(jié)構(gòu)或者構(gòu)件不能以原來的平衡狀態(tài)承受附加荷載)是緊密聯(lián)系的[4-5]。因此，為了得到桿系結(jié)構(gòu)最為真實的臨界荷載值及荷載與位移全過程曲線，在靜力推覆分析時，需同時考慮材料非線性和幾何非線性。

針對桿系結(jié)構(gòu)的靜力穩(wěn)定性分析，已有研究取得了一系列成果，主要可分為兩類[6-8]：第一類是基于構(gòu)件層面的宏觀分析方法，具有建模簡單、計算量小等特點。例如，韋良[9]基于平衡微分方程并引入穩(wěn)定函數(shù)推導(dǎo)了考慮幾何非線性的梁柱單元剛度矩陣，同時利用塑性鉸模型，提出了基于QR法的鋼框架穩(wěn)定性分析方法，能夠準確、高效地獲得結(jié)構(gòu)極限承載力及變形能力；Li等[10]、靳永強等[11]基于擬力法理論建立了軸壓構(gòu)件宏觀分析模型，通過設(shè)置兩個局部塑性機制來分別模擬彎曲變形和軸向塑性伸長行為，實現(xiàn)軸壓構(gòu)件的非線性屈曲分析。第二類是基于材料層面的非線性有限元分析方法，具有計算精度良好、計算量大等特點。例如，Nguyen等[12]采用更新拉格朗日列式建立了考慮幾何非線性行為的平衡方程，并利用纖維塑性鉸法來模擬空間框架結(jié)構(gòu)的側(cè)扭穩(wěn)定性行為，可實現(xiàn)塑性變形沿桿件長度上的開展，通過數(shù)值算例驗證了該方法的準確性和高效性；葉康生等[13]、吳可偉[14]基于共旋坐標法和塑性區(qū)法提出了一種桿系結(jié)構(gòu)彈塑性大位移分析方法，將單元端部截面分割為若干小塊面積，單元內(nèi)的截面一直保持彈性狀態(tài)，快速實現(xiàn)結(jié)構(gòu)的靜力穩(wěn)定性分析。

在桿系結(jié)構(gòu)的靜力推覆分析中，其數(shù)值模擬的難點在于：失穩(wěn)瞬間，結(jié)構(gòu)剛度發(fā)生突變，承載力驟降，荷載-位移關(guān)系曲線呈現(xiàn)下降趨勢。為了越過荷載-位移曲線的極值點，通常將采用位移控制算法來施加推覆力[15]。為了克服結(jié)構(gòu)中不同部位非線性程度的差異，滿足荷載分布恒定比例關(guān)系，且提高數(shù)值求解的穩(wěn)定性，黃羽立等[16-17]提出了基于多點位移控制的推覆分析方法。該方法能夠較好地獲取荷載-位移曲線的下降段，準確追蹤失穩(wěn)后的非線性平衡路徑，具有良好的數(shù)值穩(wěn)定性優(yōu)勢。

在構(gòu)件層面，Wong等[18]基于位移分解原理提出了結(jié)構(gòu)非線性高效分析方法-擬力法，該方法被廣泛應(yīng)用于工程結(jié)構(gòu)的非線性分析和振動控制等領(lǐng)域[19]。為了實現(xiàn)結(jié)構(gòu)的精細化高效模擬，Li等[20]將位移分解思想拓展到材料應(yīng)變層面，提出了隔離非線性有限元法，其基本原理是將材料應(yīng)變分解為線彈性與非彈性兩部分，并與有限元相結(jié)合，形成具有剛度矩陣彈塑性分離形式的控制方程，以Woodbury公式為求解工具，實現(xiàn)了局部材料非線性問題的高效求解[21]。針對結(jié)構(gòu)出現(xiàn)大范圍材料非線性問題，Li等[22]通過引入組合近似法至Woodbury公式中，提出了近似Woodbury法來高效求解控制方程。隨后，李鋼等[23]將隔離非線性有限元法嵌入到更新拉格朗日格式中，建立了考慮雙重非線性(材料和幾何非線性)行為的梁單元控制方程，并對近似Woodbury法進行了改進，高效求解整體控制方程，實現(xiàn)了桿系結(jié)構(gòu)的動力雙重非線性分析，但在求解其靜力穩(wěn)定性問題方面還有待進一步完善。

鑒于此，本文基于多點位移控制算法和近似Woodbury法，提出了桿系結(jié)構(gòu)的靜力穩(wěn)定性高效分析方法。該方法利用Woodbury公式對多點位移控制的平衡方程進行展開，結(jié)合隔離非線性理論形成具有剛度矩陣彈塑性分離形式的控制方程，采用近似Woodbury法對其進行高效求解，在桿系結(jié)構(gòu)的靜力穩(wěn)定性分析中，避免整體剛度矩陣的實時分解。以Williams雙桿體系[22]為數(shù)值算例，驗證了本文方法在靜力穩(wěn)定性分析中的可行性和有效性。最后，通過某框架結(jié)構(gòu)的靜力穩(wěn)定性分析，論證了該方法能高效地獲得荷載-位移關(guān)系曲線的下降段，處理負剛度問題，反映結(jié)構(gòu)的承載力、穩(wěn)定性及剛度的整個變化歷程。

1 多點位移控制算法概述

多點位移控制算法是一種較為新穎的靜力推覆分析方法，通過引入位移剛性約束來保持恒定推覆側(cè)力分布，以實現(xiàn)結(jié)構(gòu)穩(wěn)定性分析[16-17]。該方法在解決結(jié)構(gòu)復(fù)雜非線性問題時，特別是處理負剛度問題具有顯著優(yōu)勢，能夠獲取結(jié)構(gòu)從彈性階段到軟化階段的全過程。對于桿系結(jié)構(gòu)，其失穩(wěn)瞬間平衡路徑可能會發(fā)生“跳躍”等情況，為了精確及完整地追蹤非線性屈曲路徑，本文采用多點位移控制進行桿系結(jié)構(gòu)的靜力推覆分析。

以圖1所示的n自由度模型比例加載為例，對多點位移控制算法進行說明。在該模型上施加荷載(F1∶F2∶…∶Fn)滿足恒定比例關(guān)系為(p1∶p2∶…∶pn=1∶2∶…∶n)，應(yīng)滿足如下約束方程[16-17]：

圖1 多點位移控制的靜力推覆分析

(1)

式中：P={p1,p2,…,pn}T，pi是第i個自由度上施加荷載的比例系數(shù)；U={u1,u2,…,un}T，ui是第i個自由度所對應(yīng)的位移；u0是新增約束方程引入自由度所對應(yīng)的位移。結(jié)構(gòu)的平衡方程可寫為：

F-Q=0

(2)

式中:F={F1,F2,…,Fn}T，是關(guān)于位移U={u1,u2,…,un}T的非線性函數(shù)，表示結(jié)構(gòu)的內(nèi)力；Q={Q1,Q2,…,Qn}T，表示結(jié)構(gòu)所受的外荷載。式(1)是方程組(2)的約束方程，可基于拉格朗日乘子法構(gòu)建新的方程組，其形式如下[16-17]：

(3)

式中:λ為拉格朗日乘子，表示力加載系數(shù)。為了求解位移U和λ，可采用牛頓迭代法求解非線性方程組(3)：

(4)

其中:i表示第i個迭代步；Kt表示第i迭代步的結(jié)構(gòu)切線剛度矩陣；C0的表達式為：

(5)

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

2.1 基于多點位移控制的隔離非線性方程

為了便于推導(dǎo)，將式(4)寫為增量形式，經(jīng)整理可得：

(6)

式中:ΔU和Δλ分別表示第i迭代步至第i+1迭代步的位移增量和加載系數(shù)增量；ΔR表示迭代過程中產(chǎn)生的不平衡力(等于外荷載向量和當前迭代步的結(jié)構(gòu)內(nèi)力向量之差，即ΔR=Q-F)。對上述增量方程直接求解時，需要實時分解包含切線剛度Kt和荷載比例系數(shù)向量P組成的矩陣，這將花費大量的時間成本。為此，可基于Woodbury公式，對方程(6)中矩陣求逆部分進行展開：

(7)

將式(7)代入式(6)，經(jīng)整理可得：

ΔU=ΔUR+ΔλΔUP

(8)

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

(9)

式中：ΔUR和ΔUP分別表示不平衡力ΔR和荷載比例系數(shù)向量P所產(chǎn)生的位移增量，具體表達式分別如下：

(10)

(11)

對比式(6)和式(8)可知：在求解增量位移ΔU時，式(6)需要對包含切線剛度Kt和荷載比例系數(shù)向量P組成的矩陣進行分解，而式(8)僅需分解切線剛度矩陣Kt。根據(jù)文獻[23]，可將式(10)和式(11)表達成具有隔離非線性特征的控制方程，形式分別如下：

(12)

(13)

本文采用更新拉格朗日格式描述幾何非線性行為，因此，上述兩個方程中的所有矩陣和向量均是以t時刻的已知構(gòu)型作為參考系。由式(12)和式(13)可知，多點位移控制的平衡方程(式(6))可基于Woodbury公式進行轉(zhuǎn)換，將切線剛度矩陣Kt和荷載比例系數(shù)加載向量P拆分開來，分別形成由不平衡力ΔR和荷載比例系數(shù)加載向量P作為外荷載的平衡方程(式(10)和式(11))，再結(jié)合隔離非線性的基本理論，將式(10)和式(11)表達成具有隔離非線性特征的控制方程(式(12)和式(13))。

2.2 高效求解方法

根據(jù)文獻[23]可知，具備隔離非線性特征的控制方程(式(12)和式(13))可采用近似Woodbury法進行高效求解，其核心思想是：采用Woodbury公式對控制方程進行精確展開，假定實時變化的整體剛度矩陣在某一時間段內(nèi)保持不變，與舒爾補矩陣相關(guān)的線性方程采用組合近似法高效求解，在桿系結(jié)構(gòu)的靜力穩(wěn)定性分析中，能避免矩陣的實時更新和分解運算，可顯著提高計算效率。在某段時間區(qū)域t∈[Tj,Tj+1]，基于近似Woodbury法求解式(12)和式(13)的具體表達式分別如下：

(14)

(15)

表1 混合近似法的時間復(fù)雜度

其中:計算次序(3)中的項目可采用組合近似法進行求解，具體求解流程可參見文獻[23]。

2.3 近似矩陣更新準則

為了能夠消除矩陣近似所引起的誤差，可建立近似矩陣更新準則對近似矩陣進行“適時”更新，即[23]：

φ=

(16)

其中:|·|表示絕對值；φ表示由矩陣近似所產(chǎn)生的相對誤差；Δλ(t)表示t時刻的加載系數(shù)增量；ε0為用戶自定義參數(shù)，用于控制矩陣更新頻次[23]。矩陣Kt(t)表示結(jié)構(gòu)的切線剛度，對方程(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 求解流程

針對桿系結(jié)構(gòu)的靜力穩(wěn)定性分析，本文方法的求解流程如圖2所示。

圖2 本文方法求解流程圖

2.5 效率分析

時間復(fù)雜度是一種客觀反映算法效率的工具，能夠定量地評估算法本身在執(zhí)行過程中所需要的理論時間，不依賴于計算機性能、編程技巧等外在因素[24-25]。就式(10)和式(12)而言，當式(10)采用傳統(tǒng)有限元法(矩陣分解采用LDLT法)進行求解，式(12)采用近似的Woodbury法進行求解時，可依據(jù)文獻[23]分別將上述兩種方法的時間復(fù)雜度列于表2。

表2 三種方法的時間復(fù)雜度函數(shù)(t∈[Tj, Tj+1])

3 數(shù)值分析

3.1 模型驗證

以經(jīng)典的Williams雙桿體系[26]為數(shù)值算例(如圖3所示)，驗證本文方法在解決負剛度問題時的有效性。結(jié)構(gòu)桿件均為φ114 mm×4 mm的鋼管，材料本構(gòu)關(guān)系為理想彈塑性，屈服應(yīng)力為240 MPa，彈性模量為210 GPa。每根桿件劃分成10個單元(3結(jié)點Timoshenko纖維梁單元)；桿件截面被分割成24個小面積(24根纖維)，環(huán)向分成12等份，徑向分成2等份。

圖3 Williams 雙桿體系

圖4給出了本文方法與沈世釗等[26]模擬的豎向荷載-位移曲線，兩者吻合較好。其中，本文方法和沈世釗等[26]所獲得的臨界荷載值分別為15.55 kN和 15.56 kN，表明：本文方法在進行穩(wěn)定性分析時可精確獲得結(jié)構(gòu)的荷載-位移全過程曲線，能追蹤發(fā)生“跳躍”失穩(wěn)的平衡路徑。

圖4 Williams 雙桿體系的荷載-位移曲線

3.2 數(shù)值算例

以某3層三跨鋼框架結(jié)構(gòu)為算例，進行靜力穩(wěn)定性分析，其立面和平面圖如圖5所示。該結(jié)構(gòu)層高均為3.9 m，X方向為3跨(邊跨跨度為6 m，中跨跨度為9 m)，Z方向為4跨(跨度均為6 m)，梁柱截面尺寸如圖5所示。

圖5 3層三跨鋼框架結(jié)構(gòu)

各桿件均采用3結(jié)點Timoshenko纖維梁單元，柱劃分6個單元，邊跨梁劃分8個單元，中跨梁劃分12個單元，有限元模型總計1 068個單元，總自由度為5 850，構(gòu)件的截面纖維數(shù)均為19根。鋼材本構(gòu)關(guān)系采用雙線性隨動強化模型，其中彈性模量Ee=206 GPa，屈服強度σy=345 MPa，屈服后剛度Et=0.02Ee。恒荷載集中施加于梁柱節(jié)點，大小均為750 kN。采用倒三角側(cè)向荷載模式對結(jié)構(gòu)進行水平推覆(見圖5)；CA法的最大基向量個數(shù)取5。

圖6給出了靜力推覆荷載作用下結(jié)構(gòu)第①榀的基底剪力和頂點位移全過程曲線，可見采用本文方法能夠較好地追蹤結(jié)構(gòu)失穩(wěn)后的下降段曲線，即軟化過程，可以獲得穩(wěn)定的臨界荷載值及對應(yīng)的側(cè)向位移，分別為1 411 kN和0.441 m。圖6也繪制了僅考慮材料非線性行為時第①榀的基底剪力和頂點位移曲線，對比結(jié)果表明：幾何非線性行為導(dǎo)致結(jié)構(gòu)失穩(wěn)，在進行穩(wěn)定性分析時，需同時考慮結(jié)構(gòu)的材料和幾何非線性(雙重非線性)行為，以獲得實際的極限承載力。圖7給出了靜力推覆分析過程中各層推覆力和頂點位移關(guān)系曲線，可以看出各層施加的推覆力嚴格符合等比例荷載模式，滿足多點位移控制算法的荷載分布恒定比例關(guān)系。

圖6 第①榀基底剪力-頂點位移曲線對比圖

圖7 各層推覆力荷載值-頂點位移曲線圖

在靜力推覆分析過程中，結(jié)構(gòu)的塑性自由度數(shù)(p)時程曲線如圖8所示。其中，被激活的最大塑性自由度數(shù)為1 725，占總自由度數(shù)的29.5%，每個增量步的平均塑性自由度為1 360.5(總增量步為1 788)，占總自由度數(shù)的23.3%。

圖8 塑性自由度時程曲線

為了驗證本文方法的高效性，圖9給出了本文方法和傳統(tǒng)有限元法(采用LDLT分解法)的時間復(fù)雜度時程曲線，其中，荷載增量步的時間復(fù)雜度等于每個迭代步所需時間復(fù)雜度的總和。由圖9可知，傳統(tǒng)有限元法在每個增量步的時間復(fù)雜度僅隨迭代步數(shù)而發(fā)生變化，本文方法的時間復(fù)雜度與迭代步數(shù)、基向量個數(shù)、所激活的塑性自由度數(shù)及近似矩陣更新次數(shù)均相關(guān)，在執(zhí)行自適應(yīng)策略時，近似矩陣的更新會導(dǎo)致該增量步的時間復(fù)雜度變大。在整個靜力推覆分析過程中，矩陣更新的次數(shù)為58，與總增量步數(shù)(1 788)的比例為1∶31，表明：近似矩陣更新準則(見式(16))能較好地控制矩陣更新次數(shù)，在保證迭代收斂的前提下使計算效率最優(yōu)。此外，迭代步數(shù)也會影響算法的計算成本，表3統(tǒng)計了本文方法和傳統(tǒng)有限元法的總迭代步數(shù)，可知：本文方法所需的迭代步數(shù)要多于傳統(tǒng)有限元法，即本文方法的收斂性要低于傳統(tǒng)有限元法，原因在于：近似Woodbury法中近似矩陣所引起的誤差需要更多迭代步進行消除。同時，表3統(tǒng)計了上述兩種方法的總時間復(fù)雜度和每個荷載增量步的平均時間復(fù)雜度，本文方法和傳統(tǒng)有限元法的總時間復(fù)雜度分別為0.98×1011和5.10×1011，平均到每個增量步的時間復(fù)雜度分別為0.55×108和2.85×108，結(jié)果表明：在計算效率方面，本文方法比傳統(tǒng)有限元法快5.18倍，驗證了該方法能高效分析桿系結(jié)構(gòu)的靜力穩(wěn)定性問題。

圖9 時間復(fù)雜度時程曲線

表3 總迭代步數(shù)和時間復(fù)雜度對比

4 結(jié) 論

針對桿系結(jié)構(gòu)的靜力穩(wěn)定性問題，本文基于多點位移控制算法和近似Woodbury公式，提出了高效的靜力推覆分析方法。利用Woodbury公式對多點位移控制的平衡方程進行展開，結(jié)合隔離非線性理論形成了具有剛度矩陣彈塑性分離形式的控制方程，并采用近似Woodbury法對控制方程進行高效求解。最后，將該方法應(yīng)用于某框架結(jié)構(gòu)的靜力穩(wěn)定性分析中，并得到如下結(jié)論：

(1)本文方法是以Woodbury公式為理論基礎(chǔ)，將多點位移控制的平衡方程變?yōu)榫哂懈綦x非線性特性的控制方程，實現(xiàn)剛度矩陣彈塑性分離，并采用近似Woodbury法高效求解控制方程。

(2)本文方法能夠高效分析桿系結(jié)構(gòu)的靜力穩(wěn)定性問題，可追蹤失穩(wěn)后的非線性平衡路徑，獲得下降段曲線，處理負剛度問題，能夠準確獲得荷載-位移全過程曲線，完整反映結(jié)構(gòu)的承載力、穩(wěn)定性及剛度的整個變化歷程。