改進的降階牛頓迭代數(shù)值子結(jié)構(gòu)方法

孫寶印，申 偉，孫天舒，歐進萍,4

(1. 河海大學土木與交通學院，江蘇，南京 210098；2. 廣西大學土木建筑工程學院，廣西，南寧 530004；3. 大連理工大學土木工程學院，遼寧，大連 116024；4. 哈爾濱工業(yè)大學土木工程學院，黑龍江，哈爾濱 150090)

計算效率和模擬精度是高層建筑等重大工程地震彈塑性全過程分析的關(guān)鍵技術(shù)瓶頸。一方面，此類結(jié)構(gòu)的規(guī)模龐大，為了精確刻畫其非線性行為，有必要針對整體結(jié)構(gòu)建立精細化的有限元模型，自由度動輒數(shù)十上百萬，結(jié)構(gòu)剛度矩陣的集成和分解耗時嚴重。更為關(guān)鍵的是，結(jié)構(gòu)一旦進入非線性狀態(tài)，傳統(tǒng)的變剛度法需要反復(fù)迭代更新結(jié)構(gòu)剛度，計算效率更是成為棘手難題。另一方面，高層結(jié)構(gòu)地震倒塌破壞主要是由局部關(guān)鍵構(gòu)件嚴重損傷導(dǎo)致，而這部分關(guān)鍵構(gòu)件的非線性變形特征具有顯著的多樣性和復(fù)雜性，對模型精細程度提出了較高的要求。實際工程中為了滿足計算效率要求，無法針對整體結(jié)構(gòu)建立精細模型，進而難以保證局部關(guān)鍵構(gòu)件的模擬精度。由此可見，計算效率和模擬精度之間相互矛盾、制約?，F(xiàn)有的計算方法較難有效兼顧這類大規(guī)模結(jié)構(gòu)非線性分析的效率和精度，因此，發(fā)展更加高效、精確的非線性分析理論和計算方法十分必要。

高層結(jié)構(gòu)在地震作用下呈現(xiàn)局部非線性特征，如何充分利用這種局部化特征來提高大規(guī)模結(jié)構(gòu)非線性分析的計算效率，是近幾十年來國內(nèi)外學者關(guān)注的重點。Clough 等[1]利用子結(jié)構(gòu)方法對局部非線性的大規(guī)模結(jié)構(gòu)系統(tǒng)進行動力分析。將整體結(jié)構(gòu)劃分為線彈性和塑性區(qū)域：一方面，通過靜力凝聚消去線彈性區(qū)域內(nèi)部自由度，降低求解矩陣規(guī)模來提高計算效率；另一方面，提高塑性區(qū)模型精細程度來有效兼顧模擬精度。Wong 等[2]提出擬力法，定義能夠描述結(jié)構(gòu)構(gòu)件非線性行為的塑性自由度，建立相應(yīng)的塑性機制，來模擬結(jié)構(gòu)的非線性行為。并采用初始剛度對狀態(tài)方程進行迭代求解，利用顯式求解算法，可實現(xiàn)高效、穩(wěn)定的結(jié)構(gòu)非線性分析。隨后，Li 等[3-5]在擬力法的基礎(chǔ)上作了進一步改進，將材料應(yīng)變分解為線彈性應(yīng)變和非線性應(yīng)變兩部分，并在單元內(nèi)構(gòu)造非線性應(yīng)變的分布場模型，提出了單元非線性變形機制的模擬方法。再結(jié)合虛功原理，建立具有隔離非線性特征的結(jié)構(gòu)控制方程，在該方程迭代分析時僅需對一個小規(guī)模的Schur 補矩陣進行迭代更新，從而能夠顯著提高結(jié)構(gòu)非線性分析的計算效率。此外，陸新征等[6]、Sun 等[7-8]以及陳宇等[9]從多尺度建模角度，對整體結(jié)構(gòu)進行差異化建模分析，在保證模擬精度的前提下，可大幅降低結(jié)構(gòu)模型自由度。

傳統(tǒng)的子結(jié)構(gòu)方法、擬力法以及多尺度方法均需要預(yù)判結(jié)構(gòu)的塑性區(qū)域，并且模型網(wǎng)格一旦確認，分析過程中較難更改。為了規(guī)避這些缺陷，孫寶印等[10-13]提出數(shù)值子結(jié)構(gòu)方法(如圖1所示)，將原始高層結(jié)構(gòu)的非線性分析轉(zhuǎn)化為主結(jié)構(gòu)的等效線彈性分析和屈服構(gòu)件的隔離子結(jié)構(gòu)非線性分析。在整個地震彈塑性時程分析過程中，由線彈性主結(jié)構(gòu)模型計算得到整體結(jié)構(gòu)的動力響應(yīng)，根據(jù)結(jié)構(gòu)全局位移實時判斷構(gòu)件的屈服狀態(tài)。將構(gòu)件屈服隔離并建立精細化的子結(jié)構(gòu)模型，由子結(jié)構(gòu)非線性分析計算得到屈服構(gòu)件的真實內(nèi)力，并將其傳回主結(jié)構(gòu)。主結(jié)構(gòu)將轉(zhuǎn)化得到的屈服構(gòu)件等效外荷載作為外力項施加在屈服構(gòu)件端部節(jié)點上，最后，再由主結(jié)構(gòu)進行等效線彈性分析得到整體結(jié)構(gòu)的真實響應(yīng)。經(jīng)研究表明：該方法可以有效兼顧大規(guī)模結(jié)構(gòu)非線性分析的計算效率和模擬精度。

由于上述數(shù)值子結(jié)構(gòu)方法采用初始彈性剛度進行非線性迭代分析，仍存在一定不足。于是，本文擬對該方法作進一步改進，提出一種降階的牛頓迭代數(shù)值子結(jié)構(gòu)方法。首先，簡介數(shù)值子結(jié)構(gòu)方法基本原理及其不足；其次，提出降階的牛頓迭代算法，并介紹屈服構(gòu)件的隔離子結(jié)構(gòu)非線性分析；最后，通過一平面15 層3 跨鋼結(jié)構(gòu)的地震彈塑性時程分析來驗證所提方法的精度和效率。

1 數(shù)值子結(jié)構(gòu)方法基本原理及不足

1.1 數(shù)值子結(jié)構(gòu)方法

根據(jù)數(shù)值子結(jié)構(gòu)方法基本原理[10-13]，結(jié)構(gòu)運動方程可以表示為：

式中：U、和分別為結(jié)構(gòu)的位移、速度和加速度；M、C0、K0、F和Fc分別為結(jié)構(gòu)的質(zhì)量矩陣、初始阻尼矩陣、初始剛度矩陣、外力和非線性修正力，分別表示為：

式中：Ne和Ns分別為總單元數(shù)和屈服單元數(shù)；B(e)為用于檢索結(jié)構(gòu)中單元所在位置的布爾矩陣；為屈服單元的非線性修正力，可表示為：

式中：ke和ue分別為屈服單元的初始彈性剛度和節(jié)點位移；re為由隔離子結(jié)構(gòu)計算得到的屈服單元內(nèi)力。

式(1)等式左邊是一個線彈性系統(tǒng)，結(jié)構(gòu)非線性是由等式右邊的非線性修正力Fc來考慮。對于局部非線性結(jié)構(gòu)系統(tǒng)，屈服單元數(shù)Ns＜＜Ne遠小于總單元數(shù)Ne，則非線性修正力Fc的計算量較小。

1.2 數(shù)值子結(jié)構(gòu)方法的不足

以文獻[14]中15 層3 跨的平面鋼框架結(jié)構(gòu)為研究對象。采用不動點迭代數(shù)值子結(jié)構(gòu)方法[10-13]對其進行地震彈塑性時程分析時，主結(jié)構(gòu)中每根梁柱構(gòu)件由30 個基于歐拉-伯努利梁理論的彈性單元模擬，分析過程中根據(jù)截面最外側(cè)應(yīng)變是否超過材料屈服應(yīng)變來判斷構(gòu)件的彈塑性狀態(tài)[15]，屈服單元由隔離子結(jié)構(gòu)進行非線性分析。在隔離子結(jié)構(gòu)中，屈服單元采用基于歐拉-伯努利梁理論的分布式塑性鉸纖維單元模擬，鋼材采用OpenSees中Steel02 本構(gòu)模型模擬，材料彈性模量E=2.0×106MPa，硬化比b=0.1。在不同工況和收斂誤差下，不動點迭代數(shù)值子結(jié)構(gòu)方法的迭代次數(shù)時程曲線如圖2 所示，迭代信息如表1 所示?？梢钥闯觯斒諗空`差為1.0×10-5mm 時，8 度和9 度罕遇地震作用下，數(shù)值子結(jié)構(gòu)方法需要最多24 步和26 步收斂，平均迭代次數(shù)分別為4.6 和6.9；而當收斂誤差減小至1.0×10-8mm 時，最多需要41 步和44 步，平均迭代次數(shù)分別為7.1 和11.1。當荷載水平增大或收斂誤差減小，不動點迭代數(shù)值子結(jié)構(gòu)方法所需的迭代收斂次數(shù)均有所增加，尤其是當收斂條件變得苛刻，收斂速率顯著降低。于是，下文將對該方法作進一步改進，以提高其收斂速率。

圖2 不動點迭代數(shù)值子結(jié)構(gòu)方法迭代次數(shù)時程曲線Fig. 2 Iteration number time-history curves using the numerical substructure method based on fixed-point iteration

表1 不動點迭代數(shù)值子結(jié)構(gòu)方法迭代信息統(tǒng)計表Table 1 Table of iteration information using the numerical substructure method based on fixed-point iteration

2 降階的牛頓迭代算法

經(jīng)過時間離散后，第n時步結(jié)構(gòu)的等效運動方程滿足：

不失一般性，以Newmark-β 逐步積分法為例，第n時步速度和位移分別表示為：

式中， γ 和 β為控制參數(shù)。本文采用平均加速度法，則 γ 和 β取值分別為0.5 和0.25。

由式(5)可以得到第n時步速度和加速度：

式 中：c1=γ/(βΔt)；c2=1/β(Δt)2；a1=1-γ/β；a2=(1-γ/2β)Δt；a3=-1/βΔt；a4=1-1/2β。

將式(6)代入式(4)，得到結(jié)構(gòu)的等效平衡方程：

式中，K和Pn分別為結(jié)構(gòu)等效初始剛度和外力：

由式(7)可以得到結(jié)構(gòu)的位移場：

定義與屈服單元有關(guān)的自由度為塑性自由度，則塑性自由度位移us,n可以表示為：

式中，Bs,n為檢索塑性自由度位置的布爾矩陣。

將式(10)代入式(11)，塑性自由度位移殘差φ(us,n)可以表示為：

采用牛頓迭代算法，已知當前塑性自由度位移us,n和殘差φ(us,n)，得到更新后的塑性自由度位移：

式中，Hs,n為與塑性自由度有關(guān)的迭代矩陣：

3 隔離子結(jié)構(gòu)非線性分析

考慮到高層結(jié)構(gòu)在地震作用下發(fā)生倒塌破壞，局部關(guān)鍵構(gòu)件非線性程度較高，而大部分屈服構(gòu)件的非線性程度較低(稱之為一般屈服構(gòu)件)。為了滿足模擬精度需求，一般屈服構(gòu)件和少數(shù)關(guān)鍵構(gòu)件所采用的隔離子結(jié)構(gòu)模型的精細程度不同。為了簡化起見，本節(jié)以梁柱構(gòu)件為例說明。

3.1 一般屈服梁柱構(gòu)件

針對一般屈服梁柱構(gòu)件，纖維模型可以較為準確模擬其非線性行為。由于主、子結(jié)構(gòu)中梁柱構(gòu)件均為桿系模型，它們的邊界自由度一致，因此，由隔離子結(jié)構(gòu)可直接計算得到屈服構(gòu)件內(nèi)力和切線剛度，這里不作詳細贅述。

3.2 關(guān)鍵梁柱構(gòu)件

以如圖3(a)所示的平面框架結(jié)構(gòu)為例，將主結(jié)構(gòu)中關(guān)鍵柱構(gòu)件隔離并建立精細化的子結(jié)構(gòu)模型(例如：實體單元)，如圖3(b)所示。為了消除剛體位移的影響，將主結(jié)構(gòu)中屈服構(gòu)件對應(yīng)的彈性單元節(jié)點位移ue轉(zhuǎn)換得到局部坐標位移ve：

圖3 關(guān)鍵構(gòu)件隔離子結(jié)構(gòu)精細化模型邊界處理Fig. 3 Boundary treatment of isolated refined substructures for key components

式中：上標i為精細化子結(jié)構(gòu)第i個邊界；轉(zhuǎn)換矩陣T可以表示為：

式中， α為桿系單元軸向與全局坐標系X軸之間的夾角。

如圖3(c)所示，在局部坐標系下，根據(jù)平截面假定，精細化子結(jié)構(gòu)中第i個邊界中第j個節(jié)點的位移滿足以下轉(zhuǎn)換關(guān)系：

將式(16)代入式(18)，并集成得到全局坐標系下精細化子結(jié)構(gòu)的邊界位移轉(zhuǎn)換關(guān)系：

式中，C為全局坐標系下子結(jié)構(gòu)的邊界轉(zhuǎn)換矩陣：

式中，m1和m2分別為第1 和第2 個邊界的節(jié)點數(shù)。

根據(jù)主、子結(jié)構(gòu)邊界反力所作虛功相等，圖3(d)中邊界反力之間滿足：

將式(20)代入式(22)，得到屈服構(gòu)件內(nèi)力re：

再將式(23)代入式(3)，得到屈服構(gòu)件的非線性修正力=keue-CTqb。

隔離子結(jié)構(gòu)的增量平衡方程可以表示為：

式中：下標i和b分別為隔離子結(jié)構(gòu)內(nèi)部和邊界自由度；ks和q分別為隔離子結(jié)構(gòu)的剛度和外力；Δ為增量。

將式(20)代入式(24)，并將式(24)第二個等式左乘CT，再通過靜力凝聚，則有：

4 分析流程

采用數(shù)值子結(jié)構(gòu)方法進行結(jié)構(gòu)地震彈塑性時程分析時，具體分析流程如下：

1)初始化。令n=0、Fn=0 、=0 、=0、Bs,n=0 以及=0。根據(jù)式(8)集成等效剛度矩陣K，計算并保存其逆矩陣K-1。根據(jù)總時步NS P，得到荷載增量 dF=F/NS P。

2)令n=n+1、i=0、=Bs,n-1以 及=，計算當前時步外力Fn=Fn-1+dF，并由式(9)計算得到等效外力Pn。根據(jù)式(10)計算得到主結(jié)構(gòu)位移，更新整體結(jié)構(gòu)屈服構(gòu)件單元狀態(tài)。若屈服單元數(shù)為零，則系統(tǒng)收斂，跳到第7)步；否則，更新布爾矩陣，由式(11)得到塑性自由度位移，進入下一步。

5)根據(jù)下式判斷塑性自由度系統(tǒng)的收斂性。

5 算例分析

為了驗證本文提出的改進的降階牛頓迭代數(shù)值子結(jié)構(gòu)方法(下文簡稱本文方法)的計算效率和模擬精度，將其與傳統(tǒng)的牛頓算法(下文簡稱傳統(tǒng)方法)進行對比分析，同樣以上述15 層3 跨的平面鋼框架結(jié)構(gòu)為例。為避免平臺的差異影響對比分析的客觀性，本文算例均在OpenSees 平臺中進行，采用BandGeneral 求解器，計算硬件條件是Interl Xeon W-2155CPU，主頻是3.3 GHz，內(nèi)存為64 GB。

5.1 屈服單元隔離非線性分析

進行動力時程分析之前，將豎向重力荷載分10 步逐漸加到結(jié)構(gòu)上。圖4 給出8 度和9 度罕遇地震作用下整體結(jié)構(gòu)中隔離單元占比時程曲線。從圖4 中可以看出，結(jié)構(gòu)最終屈服單元數(shù)占比分別是4.4%和9.6%；8 度罕遇地震作用下，1.62 s時刻結(jié)構(gòu)進入彈塑性狀態(tài)，直到12 s 左右屈服單元數(shù)達到最大；9 度罕遇地震作用下，1.46 s 時刻進入屈服狀態(tài)，直到12 s 左右塑性區(qū)不再擴展。

圖4 隔離單元所占比例時程曲線Fig. 4 Isolated element percentage time-history curves

圖5 給出該平面鋼結(jié)構(gòu)在8 度和9 度罕遇地震作用下的塑性區(qū)發(fā)展歷程及大小，圖中數(shù)值1～6 表示有1 個～6 個屈服單元。由圖5 可知，8 度罕遇地震下，1 層～4 層梁端先進入屈服狀態(tài)，隨后，5 層～9 層梁端出現(xiàn)塑性區(qū)，其中，1 層～4 層梁端塑性區(qū)面積較大；9 度罕遇地震下，1 層～6 層梁端先進入屈服狀態(tài)，然后，7 層～11 層梁端和底層柱端出現(xiàn)塑性區(qū)，其中，1 層～5 層梁端塑性區(qū)擴展面積較大，底層柱端塑性區(qū)面積相對較小。

圖5 塑性區(qū)發(fā)展歷程及大小Fig. 5 Development history and size of plastic zone

5.2 準確性和可靠性分析

作為對比，采用傳統(tǒng)方法分析時，該結(jié)構(gòu)梁柱構(gòu)件同樣劃分30 個非線性纖維梁單元。傳統(tǒng)方法和本文方法得到的結(jié)構(gòu)頂點位移和基底剪力時程曲線如圖6～圖9 所示，可以看出，兩種方法的計算結(jié)果幾乎一致。

圖6 8 度罕遇地震鋼結(jié)構(gòu)頂層位移時程曲線Fig. 6 Time-history curves of top displacement for the steel structure under 8 degree rare earthquakes

圖7 9 度罕遇地震鋼結(jié)構(gòu)頂層位移時程曲線Fig. 7 Time-history curves of top displacement for the steel structure under 9 degree rare earthquakes

圖8 8 度罕遇地震鋼結(jié)構(gòu)基底剪力時程曲線Fig. 8 Time-history curves of base reaction for the steelstructure under 8 degree rare earthquakes

圖9 9 度罕遇地震鋼結(jié)構(gòu)基底剪力時程曲線Fig. 9 Time-history curves of base reaction for the steelstructure under 9 degree rare earthquakes

進一步，表2 對比兩種方法的最大頂點位移和最大基底剪力以及它們的相對誤差。由表2 可知，8 度和9 度罕遇地震作用下，結(jié)構(gòu)最大頂層位移和基底剪力的相對誤差均不超過0.2%。由式(28)計算得到頂層位移時程曲線的均方差分別為0.14%和0.04%，基底剪力時程曲線的均方差分別為0.02%和0.01%。進一步驗證了本文方法的準確性、可靠性。

表2 傳統(tǒng)方法與本文方法的誤差Table 2 Errors using the traditional and proposed methods

式中：N為總數(shù)據(jù)量；Xi和Yi為第i時刻響應(yīng)值。

5.3 計算效率分析

8 度和9 度罕遇地震作用下，傳統(tǒng)方法和本文方法的迭代次數(shù)時程曲線分別如圖10 和圖11 所示。此外，傳統(tǒng)方法、本文方法以及不動點迭代數(shù)值子結(jié)構(gòu)方法(簡稱原方法)的計算信息如表3所示。由圖10 和圖11 可知，本文方法的迭代收斂次數(shù)略高于傳統(tǒng)方法，當?shù)卣鸷奢d水平增大或迭代收斂誤差減小，本文方法的迭代收斂次數(shù)變化較小。由表3 可以看出，當收斂誤差為1.0×10-5mm時，8 度和9 度罕遇地震作用下，傳統(tǒng)方法最多分別需要4 步和5 步收斂，平均迭代次數(shù)分別為2.7和3.1，本文方法最多分別需要6 步和7 步收斂，平均迭代次數(shù)分別為3.2 和3.7，而原方法最多分別需要24 步和26 步收斂，平均迭代次數(shù)分別為4.6和6.9；當收斂誤差為1.0×10-8mm 時，8 度和9 度罕遇地震作用下，傳統(tǒng)方法最多分別需要5 步和6 步收斂，平均迭代次數(shù)分別為3.2 和3.5，本文方法最多分別需要6 步和7 步收斂，平均迭代次數(shù)分別為3.9 和4.1，而原方法最多分別需要41 步和44 步收斂，平均迭代次數(shù)分別為7.1 和11.1。由此可見，本文方法所需的荷載步最大收斂次數(shù)和平均迭代次數(shù)略大于傳統(tǒng)方法，但明顯小于文獻[10 - 13]中的不動點迭代數(shù)值子結(jié)構(gòu)方法。

表3 不同方法迭代信息統(tǒng)計表Table 3 Table of iteration information using various methods

圖10 8 度罕遇地震迭代步數(shù)時程曲線Fig. 10 Iteration number time-history curves under an 8-degree rare earthquake

圖11 9 度罕遇地震迭代步數(shù)時程曲線Fig. 11 Iteration number time-history curves under a 9-degree rare earthquake

當收斂誤差為1.0×10-5mm 時，8 度和9 度罕遇地震下傳統(tǒng)方法和本文方法的計算信息統(tǒng)計表如表4 所示。可以看出：傳統(tǒng)方法分別需要對整體結(jié)構(gòu)剛度K進行7282 次和8267 次集成和分解，計算耗時分別為313 min 和398 min，而本文方法只需在初始化過程中進行一次集成和分解，非線性分析過程中僅需對Hs矩陣進行8409 次和9743 次集成和分解，計算耗時分別為157 min和332 min。

表4 傳統(tǒng)方法與本文方法的計算信息統(tǒng)計表Table 4 Computation information using the traditional and proposed methods

圖12 給出了傳統(tǒng)方法和本文方法求解矩陣維數(shù)的時程曲線。從圖12 可以看出：隨著地震荷載水平的增大，塑性自由度增多，迭代更新和分解的Hs矩陣規(guī)模變大；8 度和9 度罕遇地震下，本文方法求解的最大Hs矩陣的維數(shù)分別為666 和1713，而傳統(tǒng)方法求解的K矩陣維數(shù)高達9315。由此可見，結(jié)構(gòu)非線性程度越低，本文方法求解的Hs矩陣規(guī)模越小，計算效率越高。

圖12 矩陣運算規(guī)模時程曲線Fig. 12 Matrix operation size time-history curves

基于變剛度的傳統(tǒng)方法需要重復(fù)集成和分解整體結(jié)構(gòu)的剛度矩陣K，以及更新所有構(gòu)件單元的狀態(tài)，而采用本文方法進行非線性分析時，僅對剛度矩陣K進行一次運算，迭代過程中只更新規(guī)模較小的Hs矩陣，且只對數(shù)量較少的屈服構(gòu)件單元進行狀態(tài)更新?？梢酝茰y，針對局部非線性的大規(guī)模結(jié)構(gòu)系統(tǒng)的彈塑性分析，本文方法較傳統(tǒng)方法具有顯著的計算優(yōu)勢。

6 結(jié)論

本文對文獻[10 - 13]中不動點迭代數(shù)值子結(jié)構(gòu)算法作進一步改進，提出了降階的牛頓迭代數(shù)值子結(jié)構(gòu)方法，并利用該方法對一平面15 層3 跨鋼框架結(jié)構(gòu)進行地震彈塑性時程分析，得出以下結(jié)論：

(1)改進的降階牛頓迭代數(shù)值子結(jié)構(gòu)方法與傳統(tǒng)牛頓算法得到的結(jié)構(gòu)頂層位移和基底剪力時程曲線一致，從而驗證了本文方法的準確性和可靠性。

(2)本文方法收斂速率遠高于不動點迭代數(shù)值子結(jié)構(gòu)方法，接近傳統(tǒng)牛頓算法的二次收斂。當?shù)卣鸷奢d水平增加或迭代收斂條件變得苛刻，不動點迭代數(shù)值子結(jié)構(gòu)方法的收斂次數(shù)明顯增多，而本文方法卻沒有顯著變化。

(3)本文方法的計算效率高于傳統(tǒng)牛頓算法。對于局部非線性結(jié)構(gòu)系統(tǒng)，結(jié)構(gòu)模型塑性自由度規(guī)模較小，本文方法求解的Hs矩陣維度遠小于傳統(tǒng)方法求解的K矩陣，計算效率較高，且隨著結(jié)構(gòu)規(guī)模的增加，計算效率優(yōu)勢將更加明顯。

(4)數(shù)值子結(jié)構(gòu)方法可以有效提高局部非線性結(jié)構(gòu)系統(tǒng)的模擬精度。隔離屈服構(gòu)件單元并建立精細化子結(jié)構(gòu)有限元模型，既可以精確捕捉局部屈服構(gòu)件的非線性行為，又可以有效揭示整體結(jié)構(gòu)的失效路徑和塑性分布。