蘇 璞，李 鋼，余丁浩

(大連理工大學(xué)海岸和近海工程國家重點實驗室，遼寧，大連 116024)

隨著工程結(jié)構(gòu)規(guī)模和復(fù)雜程度的日益增加，非線性分析已經(jīng)成為研究結(jié)構(gòu)損傷演化規(guī)律的重要手段。材料非線性是土木工程領(lǐng)域常見的一種非線性行為，準(zhǔn)確模擬結(jié)構(gòu)材料非線性行為對結(jié)構(gòu)性能評估意義重大。

對于材料非線性問題，數(shù)值求解過程需重復(fù)合成和分解整體剛度矩陣，尤其是自由度規(guī)模比較大時，計算成本大幅增加，耗費時間通常難以接受，盡管計算機性能的不斷提高在一定程度上緩解了這一問題，但高效的非線性數(shù)值算法依然是眾多學(xué)者追求的目標(biāo)，許多分析方法不斷涌現(xiàn)，如多尺度方法[1－3]，子結(jié)構(gòu)方法[4－6]等。其中，子結(jié)構(gòu)概念最初用于求解線彈性問題[7]，Clough等[8]將其進一步擴展到非線性問題中，其核心思想是將整體結(jié)構(gòu)劃分為線彈性子結(jié)構(gòu)和非線性子結(jié)構(gòu)，通過靜力凝聚的方式消去彈性子結(jié)構(gòu)的內(nèi)部自由度，大幅降低了系統(tǒng)自由度數(shù)目，提高了計算效率，之后該方法被眾多學(xué)者進一步完善，并應(yīng)用于不同結(jié)構(gòu)的非線性分析[9－12]，但該方法需要提前預(yù)測非線性發(fā)生的區(qū)域，其發(fā)展受到一定限制，為克服其局限性，Han和Abel[13]提出了自適應(yīng)子結(jié)構(gòu)方法，該方法可根據(jù)結(jié)構(gòu)非線性狀態(tài)動態(tài)的劃分彈性和非線性子結(jié)構(gòu)，無需預(yù)判非線性發(fā)生的區(qū)域；Sheu等[14－15]基于多級剖分策略提出了可考慮非線性區(qū)域變化的子結(jié)構(gòu)分析方法。此外，有學(xué)者基于混合變分原理提出了雙重子結(jié)構(gòu)方法[16－17]，該方法以界面力為未知量，其界面方程規(guī)模一般小于或等于傳統(tǒng)的子結(jié)構(gòu)方法[18]。盡管子結(jié)構(gòu)方法降低了結(jié)構(gòu)整體剛度矩陣的階數(shù)，但在材料非線性分析過程中仍不可避免的需要對子結(jié)構(gòu)剛度矩陣進行分解運算，該過程通常耗時較多，仍會在一定程度上降低整體效率。

對于大多數(shù)工程結(jié)構(gòu)，材料非線性行為一般僅發(fā)生在局部區(qū)域，即結(jié)構(gòu)整體剛度矩陣只有少量元素發(fā)生變化，傳統(tǒng)變剛度法須對整體剛度矩陣實時更新和分解，無法充分考慮結(jié)構(gòu)的局部非線性特征。而在數(shù)學(xué)領(lǐng)域，對于僅有少量元素變化的矩陣，可通過Woodbury公式快速求逆，目前已有許多學(xué)者將該公式應(yīng)用于局部非線性問題[19－22]。Woodbury公式的使用避免了整體剛度矩陣的實時更新和分解，僅需對小規(guī)模的Schur補矩陣進行分解，進而有效提升計算效率。為了將Woodbury公式應(yīng)用于一般有限元問題，李鋼等[23]基于材料應(yīng)變分解的思想提出了隔離非線性有限元法，在該方法中，結(jié)構(gòu)的非線性特征由額外設(shè)置的非線性自由度來體現(xiàn)，而整體剛度矩陣始終保持不變，當(dāng)結(jié)構(gòu)非線性規(guī)模較小時，結(jié)構(gòu)非線性自由度數(shù)也較少，結(jié)合Woodbury公式僅需對與非線性自由度相關(guān)的小規(guī)模的Schur補矩陣進行分解，從而實現(xiàn)非線性問題的高效求解[24]。然而，Schur補矩陣通常為滿陣，當(dāng)結(jié)構(gòu)非線性規(guī)模比較大時，Schur補矩陣階數(shù)也隨之增大，對其進行分解的計算成本大幅增加，導(dǎo)致Woodbury公式的高效性受到限制。

本文以子結(jié)構(gòu)和隔離非線性有限元法為基礎(chǔ)，提出基于子結(jié)構(gòu)的Woodbury非線性分析方法；該方法將整體結(jié)構(gòu)劃分為若干個子結(jié)構(gòu)的同時，也將Schur補矩陣分解為若干個子矩陣(即子結(jié)構(gòu)的Schur補矩陣)，有效降低了Schur補矩陣的規(guī)模；且整個分析過程僅需對子結(jié)構(gòu)的Schur補矩陣進行分解，改進了 Schur補矩陣階數(shù)較大時 Woodbury公式計算效率降低的不足，進一步拓展了Woodbury公式的適用范圍；依據(jù)時間復(fù)雜度理論建立效率分析模型，并通過數(shù)值算例證明了本文方法的準(zhǔn)確性和高效性。

1 基本理論

1.1 子結(jié)構(gòu)法

子結(jié)構(gòu)法是將整體結(jié)構(gòu)的求解域劃分為若干個子域，各子域?qū)?yīng)的部分被定義為子結(jié)構(gòu)，子結(jié)構(gòu)間相鄰的部分稱為界面；同時，整體結(jié)構(gòu)位移自由度隨之分為內(nèi)部自由度和界面自由度，通過靜力凝聚消去內(nèi)部自由度，即可得到一個僅與界面自由度相關(guān)的控制方程。因此，子結(jié)構(gòu)法是通過“化整為零”的方式縮減問題規(guī)模，其求解流程為：將整體結(jié)構(gòu)劃分為若干個子結(jié)構(gòu)；裝配子結(jié)構(gòu)得到界面方程；利用界面信息得到每個子結(jié)構(gòu)的響應(yīng)[8]。

假定某結(jié)構(gòu)求解域為?，邊界為Γ，如圖1(a)所示；依據(jù)子結(jié)構(gòu)法將整體結(jié)構(gòu)劃分為若干個子結(jié)構(gòu)，這里以兩個子結(jié)構(gòu)為例，如圖1(b)所示；其中：子結(jié)構(gòu)1和子結(jié)構(gòu)2對應(yīng)的求解域分別為?1和?2；Γ1和Γ2分別為?1和?2所獨有的邊界，而Γ12為?1和?2的公共界面。

圖1 子結(jié)構(gòu)方法示意Fig.1 Diagram of substructuring method

這里將由子結(jié)構(gòu)組成的系統(tǒng)定義為子結(jié)構(gòu)系統(tǒng)[25]?？紤]界面處子結(jié)構(gòu)間的相互作用力以及位移連續(xù)性條件，則子結(jié)構(gòu)系統(tǒng)的平衡方程和位移協(xié)調(diào)方程分別為[17]：

式中：s為子結(jié)構(gòu)編號；NS為子結(jié)構(gòu)個數(shù)，這里NS=2；K(s)、u(s)和F(s)分別為子結(jié)構(gòu)s的剛度矩陣、位移向量和荷載向量；λ為拉格朗日乘子；B(s)是與位移自由度相關(guān)的矩陣。由式(1)解得u(s)，并將其代入式(2)，可得到子結(jié)構(gòu)系統(tǒng)的界面方程為：

其中：為子結(jié)構(gòu)s剛度矩陣K(s)的廣義逆[16]，當(dāng)K(s)非奇異時，。式(3)通常采用共軛梯度法來求解，將λ回代到式(1)，即可得到各子結(jié)構(gòu)的位移響應(yīng)。在整個求解過程中，僅需對各子結(jié)構(gòu)剛度矩陣K(s)進行更新和分解，避免了對結(jié)構(gòu)整體剛度矩陣的相應(yīng)運算，可有效降低計算成本。

1.2 基于Woodbury公式的結(jié)構(gòu)分析方法

Woodbury公式(詳見附錄)現(xiàn)已被廣泛應(yīng)用于高效求解結(jié)構(gòu)非線性問題[19－22]，該類方法可統(tǒng)一稱為基于 Woodbury公式的結(jié)構(gòu)分析方法(以下簡稱Woodbury法)。隔離非線性有限元法作為一種高效的Woodbury法，其以有限元基本理論為基礎(chǔ)，具有通用性好，適應(yīng)范圍廣的特點。

隔離非線性有限元法的核心思想是將材料應(yīng)變分解為線性和非線性兩部分。以圖2所示的單軸應(yīng)力-應(yīng)變關(guān)系作簡要介紹，假設(shè)某個增量步內(nèi)，材料非線性狀態(tài)由A點到B點，應(yīng)變增量Δε根據(jù)初始初始彈性模量Ee可分解為：

式中：Δε?為線彈性應(yīng)變增量；Δε?為非線性應(yīng)變增量。與應(yīng)變增量Δε對應(yīng)的應(yīng)力增量Δσ可表示為：

當(dāng)考慮多維材料本構(gòu)關(guān)系時，上述應(yīng)變分解過程同樣適用，此時式(4)和式(5)可相應(yīng)寫為：

式中：Δε、Δε?和 Δε?分別為應(yīng)變向量、線彈性應(yīng)變向量和非線性應(yīng)變向量；Δσ為應(yīng)力向量；De為材料初始彈性本構(gòu)矩陣。

圖2 應(yīng)變非線性分解Fig.2 Nonlinear decomposition of material strain

基于上述應(yīng)變分解思想，結(jié)合虛功原理可得到隔離非線性有限元法的控制方程為[23]：

式中：Ke為結(jié)構(gòu)的初始彈性剛度矩陣；Kr為系數(shù)矩陣，代表結(jié)構(gòu)非線性自由度與節(jié)點力之間的關(guān)系；Krr為分塊對角矩陣，體現(xiàn)結(jié)構(gòu)材料非線性變形信息；Δu和分別為結(jié)構(gòu)的結(jié)點位移向量和非線性應(yīng)變向量；ΔF為結(jié)構(gòu)結(jié)點荷載向量。采用Woodbury公式求解上述控制方程，可得：

2 基于子結(jié)構(gòu)的Woodbury法

2.1 控制方程

依據(jù)有限元法對圖1(a)所示結(jié)構(gòu)的求解域進行離散化，如圖3(a)所示；采用子結(jié)構(gòu)法將整體結(jié)構(gòu)劃分為若干個子結(jié)構(gòu)，這里以NS=2為例，相應(yīng)子結(jié)構(gòu)系統(tǒng)如圖3(b)所示。圖3中深色部分表示結(jié)構(gòu)在外荷載作用下可能發(fā)生非線性變形的區(qū)域。各子結(jié)構(gòu)結(jié)點分為界面結(jié)點和內(nèi)部結(jié)點，子結(jié)構(gòu)s的結(jié)點位移自由度N(s)相應(yīng)表示為：

圖3 子結(jié)構(gòu)剖分示意Fig.3 Diagram of substructuring partition

將1.2節(jié)中應(yīng)變分解思想應(yīng)用于各子結(jié)構(gòu)，并考慮子結(jié)構(gòu)之間的變形協(xié)調(diào)關(guān)系，可得子結(jié)構(gòu)系統(tǒng)的平衡方程和位移協(xié)調(diào)方程分別為：

式中：

由式(3)、式(15)、式(16)、式(17)和式(18)可知，子結(jié)構(gòu)法在求解界面方程時，需對子結(jié)構(gòu)剛度矩陣K(s)進行更新的分解，而基于子結(jié)構(gòu)的Woodbury法僅需對小規(guī)模的矩陣進行分解。

2.2 方程求解

本文控制方程的迭代求解方案基于 Newton-Raphson(簡稱N-R)格式，每一個N-R迭代步需依次求解式(16)和式(14)，其中式(16)采用共軛梯度法求解，具體迭代格式如圖4所示。

由圖4可知，在共軛梯度法求解Δλ的過程中，由于合成ΔFB和計算乘積KBpk涉及到子結(jié)構(gòu)Schur補矩陣的分解或回代運算，計算量較大；其余基本是矩陣和向量或向量和向量的乘積運算，計算量相對較小，可忽略不計。

求解式(16)得到Δλ后，將其回代到式(14)即可得到各子結(jié)構(gòu)的位移 Δu(s)。對比式(9)和式(14)可知，傳統(tǒng)Woodbury法需分解整體結(jié)構(gòu)的Schur補矩陣Ksch，而基于子結(jié)構(gòu)的 Woodbury法僅需分解各子結(jié)構(gòu)的 Schur補矩陣(s=1,2,…,NS)。圖5給出了Ksch和的具體形式，由圖5可知，的階數(shù)與整體結(jié)構(gòu)非線性自由度數(shù)相關(guān)，而的階數(shù)僅與子結(jié)構(gòu)s的非線性自由度數(shù)相關(guān)，因此，的階數(shù)一般低于Ksch的階數(shù)，且的階數(shù)隨子結(jié)構(gòu)數(shù)的增多而降低。

圖4 共軛梯度法迭代格式Fig.4 Iterative format of conjugate gradient method

圖5 Schur補矩陣的階數(shù)示意Fig.5 Dimension of Schur complement matrix

此外，通過對比式(14)、式(15)、式(17)和式(18)可知，回代求解位移 Δu(s)、計算乘積KBpk和合成ΔFB的計算過程相差不大。例如對于乘積KBpk，其計算過程從右到左依次執(zhí)行，首先，計算向量pk與矩陣的乘積，即：

圖6 qk 具體計算步驟Fig.6 Detailed calculation steps of qk

由此，便完成了乘積KBpk的計算(上述變量均為中間變量)。回代求解位移Δu(s)和合成ΔFB均可分解成上述類似的步驟進行運算。值得注意的是，由于計算ΔFB的過程中已經(jīng)對矩陣進行了分解，因此后續(xù)在計算KBpk和回代求解位移 Δu(s)時，只需要對矩陣進行回代即可。

3 時間復(fù)雜度效率模型

時間復(fù)雜度作為一種僅與算法本身相關(guān)的效率評價方法，具體是指算法執(zhí)行過程中所需的計算工作量，能相對客觀地評價一個算法的效率[26－27]。

根據(jù)2.2節(jié)論述可知，基于子結(jié)構(gòu)的Woodbury法在求解過程中，合成ΔFB、計算KBpk和回代求解位移 Δu(s)三項由于涉及到子結(jié)構(gòu) Schur補矩陣的分解或回代運算，計算量較大，其余運算計算成本相對較小，可忽略不計。此外，上述三項的計算過程也相似，以乘積KBpk為例，分別統(tǒng)計其計算過程每一步驟的時間復(fù)雜度并匯總，即可得到計算KBpk的時間復(fù)雜度，見表1。合成ΔFB和回代求解位移 Δu(s)的時間復(fù)雜度也可通過類似的方式得到，本文不再贅述。將上述三項的時間復(fù)雜度匯總并忽略低階項，即可得到基于子結(jié)構(gòu)的 Woodbury法的時間復(fù)雜度，見表2。

表1 計算KBpk的時間復(fù)雜度Table 1 Time complexity of calculating KBpk

結(jié)合式(22)，則基于子結(jié)構(gòu)的Woodbury法的時間復(fù)雜度TPro可表示為：

由式(23)可知，界面方程共軛梯度法求解迭代次數(shù)NCG和非線性自由度數(shù)是影響基于子結(jié)構(gòu)的Woodbury法時間復(fù)雜度的重要參數(shù)，NCG和越大，該方法的時間復(fù)雜度越高，即效率越低。

4 算例

4.1 工程概況

某7層鋼框架結(jié)構(gòu)，詳細(xì)尺寸如圖7(a)和圖7(b)所示；框架梁柱截面均為工字型，具體信息如圖7(c)所示。每根梁和柱沿長度方向被劃分為5個單元，共計2940個單元，其類型為3結(jié)點Timoshenko梁，結(jié)點自由度總數(shù)為 15456。所有單元截面的纖維數(shù)量均為31根，上、下翼緣各劃分11根纖維，中間腹板劃分9根纖維。在梁、柱節(jié)點處布置集中質(zhì)量，分布情況見表3。鋼材初始彈性模量為Ee=200 GPa，屈服應(yīng)力為σy=300 MPa，應(yīng)力-應(yīng)變關(guān)系為線性硬化模型，硬化模量為Eh=0.02Ee。地震波選取臺灣集集地震波，輸入方向為平面雙向，兩個方向的峰值加速度均調(diào)幅到1.2g。時間間隔為0.01 s，總分析時長為30 s。

表2 基于子結(jié)構(gòu)的Woodbury法的時間復(fù)雜度Table 2 Time complexity of Woodbury approach based on substructuring method

圖7 7層7跨鋼框架Fig.7 Seven-bay seven-story steel frame

表3 梁柱節(jié)點集中質(zhì)量Table 3 Joint masses of structure

4.2 結(jié)果分析

為驗證子結(jié)構(gòu)數(shù)對本文方法計算結(jié)果的影響，根據(jù)圖8所示的剖分方案將整體結(jié)構(gòu)劃分為2個子結(jié)構(gòu)，4個子結(jié)構(gòu)和8個子結(jié)構(gòu)三種工況，并分別進行非線性動力時程分析，最終得到三種工況下的頂層位移時程曲線如圖9(a)所示。結(jié)果表明：基于子結(jié)構(gòu)的Woodbury法三種工況的計算結(jié)果基本一致，其相對誤差在萬分之一以下，即子結(jié)構(gòu)數(shù)對該方法計算結(jié)果基本沒有影響。

為進一步驗證本文方法的準(zhǔn)確性，分別采用Woodbury法和基于子結(jié)構(gòu)的Woodbury法對該框架結(jié)構(gòu)進行非線性動力時程分析，兩種方法計算得到的頂層加速度和位移時程曲線分別如圖9(b)和圖9(c)所示。結(jié)果表明：兩種方法的頂層時程曲線基本重合，說明基于子結(jié)構(gòu)的Woodbury法的計算精度與Woodbury法基本一致。此外，在單元的非線性狀態(tài)判定過程中，可確定某一單元所有非線性插值點處不為0的非線性應(yīng)變個數(shù)，此即為該單元的非線性自由度數(shù)；通過統(tǒng)計某一分析步所有單元的非線性自由度數(shù)之和，即可得到該分析步結(jié)構(gòu)非線性自由度數(shù)，圖10給出了所有分析步結(jié)構(gòu)非線性自由度數(shù)曲線；為直觀反映結(jié)構(gòu)非線性規(guī)模，定義非線性占比γ，其為結(jié)構(gòu)非線性自由度數(shù)和結(jié)點位移自由度數(shù)的比值，即γ=Nr/N。由圖10可知，結(jié)構(gòu)最大非線性自由度數(shù)為1995，此時非線性占比為γ≈ 1 2.9%，而在整個分析過程中，結(jié)構(gòu)平均非線性自由度數(shù)為 371，相應(yīng)的平均非線性占比為γ≈ 2.4%。圖11給出了基于子結(jié)構(gòu)的 Woodbury法三種工況下的共軛梯度法迭代次數(shù)曲線，由圖11可知，隨著子結(jié)構(gòu)數(shù)的增加，共軛梯度法迭代次數(shù)也相應(yīng)增加，其中，相比2個子結(jié)構(gòu)而言，4個子結(jié)構(gòu)的迭代次數(shù)僅在非線性占比較高的分析步增加明顯，而8個子結(jié)構(gòu)的迭代次數(shù)則幾乎在每一分析步均有明顯增長。

圖8 子結(jié)構(gòu)剖分方案Fig.8 Partitioning schemes of substructures

圖9 頂層時程曲線Fig.9 Time-history curves of top floor

圖10 非線性自由度曲線Fig.10 Curves of nonlinear degrees of freedom

圖11 共軛梯度法迭代次數(shù)曲線Fig.11 Curves of iterations of conjugate gradient method

4.3 效率評價

考慮剛度矩陣帶寬與結(jié)構(gòu)結(jié)點自由度數(shù)的關(guān)系，Woodbury法一次N-R迭代求解的時間復(fù)雜度TIS為[27]：

根據(jù)非線性分析過程中每一次 N-R迭代步結(jié)構(gòu)非線性自由度數(shù)(圖(10))以及求解界面方程所需的共軛梯度法迭代次數(shù)(圖(11))，結(jié)合式(24)和式(23)可統(tǒng)計出 Woodbury法和基于子結(jié)構(gòu)的Woodbury法每一N-R迭代步的時間復(fù)雜度，表4給出了兩種方法的最大時間復(fù)雜度和平均時間復(fù)雜度，圖12給出了基于子結(jié)構(gòu)的Woodbury法三種工況下的時間復(fù)雜度曲線。

表4 兩種方法的時間復(fù)雜度對比Table 4 Time complexity comparison of two methods

圖12 時間復(fù)雜度曲線Fig.12 Curves of time complexity

由表4可知，基于子結(jié)構(gòu)的Woodbury法的最大時間復(fù)雜度和平均時間復(fù)雜度均明顯低于Woodbury法。值得注意的是，隨著劃分子結(jié)構(gòu)數(shù)的增多，基于子結(jié)構(gòu)的Woodbury法的最大時間復(fù)雜度逐漸降低；而平均時間復(fù)雜度雖然也有降低的趨勢，但8個子結(jié)構(gòu)時的平均時間復(fù)雜度要高于4個子結(jié)構(gòu)，其主要原因可結(jié)合圖10~圖12來說明。通過對比圖10~圖12可知，基于子結(jié)構(gòu)的Woodbury法的時間復(fù)雜度隨分析步的變化規(guī)律與非線性自由度數(shù)和共軛梯度法迭代次數(shù)(以下簡稱迭代次數(shù))類似，這再一次說明非線性自由度數(shù)和迭代次數(shù)是影響本文算法時間復(fù)雜度的兩個關(guān)鍵因素；通過進一步分析兩個因素對每一分析步時間復(fù)雜度的影響比重可知，在非線性占比較高的分析步，非線性自由度數(shù)起控制作用，此即基于子結(jié)構(gòu)的Woodbury法的最大時間復(fù)雜度隨著子結(jié)構(gòu)數(shù)的增多而遞減的主要原因；而在非線性占比較低的分析步，迭代次數(shù)起控制作用，此即大多數(shù)分析步中8個子結(jié)構(gòu)的時間復(fù)雜度高于4個子結(jié)構(gòu)的主要原因。

5 結(jié)論

本文基于子結(jié)構(gòu)法有效降低了Schur補矩陣的規(guī)模，同時將Woodbury公式應(yīng)用于各子結(jié)構(gòu)，僅需對子結(jié)構(gòu)Schur補矩陣進行分解，避免了整體結(jié)構(gòu)Schur補矩陣的相應(yīng)運算；依據(jù)時間復(fù)雜度理論建立效率模型，并通過數(shù)值算例驗證了基于子結(jié)構(gòu)的Woodbury法的準(zhǔn)確性和高效性。得到的主要結(jié)論如下：

(1)與 Woodbury法相比，基于子結(jié)構(gòu)的Woodbury法在保證計算精度的前提下進一步提高了Woodbury公式的計算效率，拓寬了其應(yīng)用范圍；

(2)子結(jié)構(gòu)數(shù)對基于子結(jié)構(gòu)的 Woodbury法的計算精度影響極小，但其對該方法計算效率影響較大。為保證本文方法具有較好的計算性能，可根據(jù)“界面結(jié)點自由度數(shù)與結(jié)構(gòu)結(jié)點自由度數(shù)的比值不超過10%”的經(jīng)驗準(zhǔn)則來選取合適的子結(jié)構(gòu)數(shù)。

附錄：Woodbury公式

Woodbury公式是用來對僅有少量元素變化(低秩修正)的矩陣快速求逆的數(shù)學(xué)工具。假設(shè)某個n×n階的矩陣A可逆，且其逆已知，A經(jīng)低秩修正后得到矩陣B，即B=A+WVT，其中W和V均為n×p階的矩陣，且p遠(yuǎn)小于n。若矩陣A和(I+VTA-1W)可逆，則矩陣B的逆可通過Woodbury公式進行高效求解[19]，即：

由于A-1已知，式(A1)的主要計算量集中于一個p×p階的矩陣(I+VTA-1W)的合成和求逆，且低秩修正時，該矩陣的階數(shù)p遠(yuǎn)小于矩陣B的階數(shù)n(即p<