邱棟星，汪 磊，孫寶印，古 泉*

(1.廈門(mén)大學(xué)建筑與土木工程學(xué)院，福建廈門(mén)361005；2.大連理工大學(xué)建設(shè)工程學(xué)部，遼寧大連116024)

近年來(lái)，隨著我國(guó)經(jīng)濟(jì)快速增長(zhǎng)，超大型結(jié)構(gòu)的建筑越來(lái)越多.一方面，這些重大工程的投資巨大且影響深遠(yuǎn)，對(duì)結(jié)構(gòu)的安全性提出了更高的要求；另一方面，近些年來(lái)，地震災(zāi)害頻發(fā)，地震導(dǎo)致結(jié)構(gòu)破壞甚至倒塌是造成經(jīng)濟(jì)損失與人員傷亡的重要原因之一.因此，研究超大型結(jié)構(gòu)的破壞全過(guò)程和倒塌機(jī)制對(duì)于其安全性設(shè)計(jì)意義重大，針對(duì)這一問(wèn)題提出高效且準(zhǔn)確的數(shù)值計(jì)算方法非常必要.

為了深入研究結(jié)構(gòu)的地震損傷和倒塌機(jī)制，對(duì)大型結(jié)構(gòu)進(jìn)行有限元分析時(shí)，通常需要對(duì)整體結(jié)構(gòu)建立精細(xì)化數(shù)值模型[1].考慮到大型復(fù)雜結(jié)構(gòu)精細(xì)化模型的單元數(shù)目以及自由度數(shù)動(dòng)輒10萬(wàn)～100萬(wàn)，導(dǎo)致計(jì)算效率比較低.更重要的是，當(dāng)某些局部關(guān)鍵構(gòu)件進(jìn)入非線性時(shí)，其結(jié)構(gòu)整體剛度矩陣會(huì)改變.在常規(guī)的非線性有限元計(jì)算過(guò)程中，每一時(shí)步求解非線性方程組時(shí)都需要進(jìn)行整體剛度矩陣的集成和三角分解，非常耗時(shí)，嚴(yán)重降低了計(jì)算效率[2].由孫寶印等[3]提出的數(shù)值子結(jié)構(gòu)方法將進(jìn)入非線性的局部構(gòu)件從整體有限元模型中隔離出來(lái),單獨(dú)進(jìn)行精細(xì)化模擬和分析，而整體模型仍然保持線彈性，能有效地提高計(jì)算效率；林純等[4]和孫寶印等[5]基于隔離數(shù)值子結(jié)構(gòu)方法對(duì)高層建筑結(jié)構(gòu)進(jìn)行的動(dòng)力分析驗(yàn)證了這種方法的高效性和準(zhǔn)確性.

然而此方法目前仍存在不足.非線性修正力項(xiàng)通過(guò)隔離子結(jié)構(gòu)模型的靜力分析得到，即僅僅考慮了靜力修正，慣性力在主結(jié)構(gòu)中計(jì)算而不需修正.此方法不適用于動(dòng)力子結(jié)構(gòu)問(wèn)題；且使用不迭代算法的數(shù)值子結(jié)構(gòu)方法要求計(jì)算步長(zhǎng)很小，導(dǎo)致計(jì)算時(shí)間較長(zhǎng)，計(jì)算效率不高[6-9].針對(duì)這兩個(gè)問(wèn)題，本研究拓展了數(shù)值子結(jié)構(gòu)方法，將非線性修正力計(jì)算從子結(jié)構(gòu)靜力分析拓展到動(dòng)力分析問(wèn)題中；同時(shí)提出隔離子結(jié)構(gòu)的響應(yīng)預(yù)測(cè)修正方法，可顯著增加計(jì)算步長(zhǎng)，極大地提高計(jì)算效率.基于改進(jìn)的數(shù)值子結(jié)構(gòu)方法，對(duì)梁柱和平面連續(xù)構(gòu)件進(jìn)行動(dòng)力分析，通過(guò)與完整結(jié)構(gòu)的標(biāo)準(zhǔn)解進(jìn)行對(duì)比，系統(tǒng)驗(yàn)證了此改進(jìn)方法的計(jì)算精度和計(jì)算誤差.

1 數(shù)值子結(jié)構(gòu)方法及其改進(jìn)方法

數(shù)值子結(jié)構(gòu)方法將超大型結(jié)構(gòu)的非線性動(dòng)力問(wèn)題分解為兩個(gè)較簡(jiǎn)單的問(wèn)題，即適度規(guī)模的主結(jié)構(gòu)非精細(xì)化分析和少量非線性隔離子結(jié)構(gòu)的精細(xì)化分析問(wèn)題.主結(jié)構(gòu)采用線彈性非精細(xì)化模型，計(jì)算規(guī)模適度；在整個(gè)計(jì)算過(guò)程中剛度保持不變，只需要進(jìn)行一次剛度矩陣集成和三角分解，效率很高.少量非線性區(qū)域采用精細(xì)化隔離子結(jié)構(gòu)模型，規(guī)模較小，計(jì)算效率高.主結(jié)構(gòu)與子結(jié)構(gòu)之間的的邊界滿足位移邊界協(xié)調(diào)和力平衡條件.本研究中主、子結(jié)構(gòu)之間的信息傳遞使用客戶機(jī)/服務(wù)器(C/S)集成技術(shù)，可實(shí)現(xiàn)基于網(wǎng)絡(luò)通道的并行分布式計(jì)算[10].以下簡(jiǎn)要介紹數(shù)值子結(jié)構(gòu)方法的計(jì)算公式.

基于有限元的整體結(jié)構(gòu)運(yùn)動(dòng)方程為[11]：

(1)

1.1 基于子結(jié)構(gòu)靜力分析法的非線性修正力

將式(1)按Newmark-β法對(duì)時(shí)間t進(jìn)行離散，則

(2)

在式(2)的左右兩邊同時(shí)加上Kun+1項(xiàng)并整理，可得

(3)

(4)

(5)

(6)

如圖1所示，ub為子結(jié)構(gòu)邊界上的節(jié)點(diǎn)位移.假設(shè)子結(jié)構(gòu)邊界上節(jié)點(diǎn)位移是連續(xù)的，ub可根據(jù)式(7)插值得到：

ub=NuNL.

(7)

其中：N為形函數(shù)，用于插值得到子結(jié)構(gòu)連續(xù)邊界上的節(jié)點(diǎn)位移；uNL(見(jiàn)圖1)為線彈性主結(jié)構(gòu)隔離區(qū)域的節(jié)點(diǎn)位移.在每一個(gè)小的時(shí)步中，式(7)也可以被寫(xiě)成以下增量形式：

δub=NδuNL.

(8)

根據(jù)虛功原理，主結(jié)構(gòu)隔離系統(tǒng)與子結(jié)構(gòu)系統(tǒng)的虛功是相等的，即

(9)

(10)

將式(8)帶入到式(10)中可得

RNL=NTRsub,b.

(11)

由此可得隔離區(qū)域的非線性修正力在子結(jié)構(gòu)進(jìn)行精細(xì)化模擬時(shí)的計(jì)算公式：

(12)

圖1 精細(xì)化子結(jié)構(gòu)模型的非線性修正力計(jì)算方法Fig.1Calculation method of nonlinear correction force for refined substructure model

以上介紹的是子結(jié)構(gòu)模型進(jìn)行靜力分析時(shí)的非線性修正力的計(jì)算，其只考慮位移邊界u，而對(duì)子結(jié)構(gòu)的慣性力項(xiàng)并未考慮，無(wú)法準(zhǔn)確模擬子結(jié)構(gòu)的動(dòng)力問(wèn)題.

因此，以下介紹子結(jié)構(gòu)動(dòng)力分析時(shí)的非線性修正力的公式推導(dǎo).

1.2 基于子結(jié)構(gòu)動(dòng)力分析法的非線性修正力

動(dòng)力分析時(shí)把運(yùn)動(dòng)方程對(duì)時(shí)間t進(jìn)行離散，每一時(shí)步都滿足以下方程：

(13)

(14)

(15)

其中

(16)

而等效動(dòng)力剛度和改進(jìn)的非線性修正力分別為：

(17)

(18)

(19)

根據(jù)虛功原理，主結(jié)構(gòu)隔離系統(tǒng)與子結(jié)構(gòu)系統(tǒng)的虛功是相等的，即

[RNL+(MNLüNL)]TδuNL=[Rsub+

(Msub,büsub,b)]T,[Rsub,i+

(20)

其中：MNL為主結(jié)構(gòu)隔離區(qū)域的節(jié)點(diǎn)質(zhì)量，Msub為子結(jié)構(gòu)系統(tǒng)的節(jié)點(diǎn)質(zhì)量，Msub,b、Msub,i分別為子結(jié)構(gòu)系統(tǒng)的邊界節(jié)點(diǎn)和內(nèi)部節(jié)點(diǎn)質(zhì)量.由于內(nèi)部節(jié)點(diǎn)不受外力作用，其內(nèi)部節(jié)點(diǎn)的結(jié)構(gòu)抗力Rsub,i與慣性力Msub,iüsub,i的和始終為0.因此，同理可得

RNL+(MNLüNL)=NT[Rsub,b+(Msub,büsub,b)],

(21)

將式(21)代入(18)，可得子結(jié)構(gòu)動(dòng)力分析時(shí)隔離區(qū)域的非線性修正力的計(jì)算公式：

NT(Msub,büsub,b)].

(22)

1.3 算法流程簡(jiǎn)介

本研究分別采用不迭代與迭代算法對(duì)隔離數(shù)值子結(jié)構(gòu)中的動(dòng)力問(wèn)題進(jìn)行分析.其中不迭代算法的流程可分為以下幾個(gè)步驟：

1.4 預(yù)測(cè)改進(jìn)的隔離數(shù)值子結(jié)構(gòu)方法

(23)

預(yù)測(cè)改進(jìn)的隔離數(shù)值子結(jié)構(gòu)方法可分為以下幾個(gè)步驟：

2 隔離數(shù)值子結(jié)構(gòu)方法的驗(yàn)證(梁柱單元)

本節(jié)采用一根長(zhǎng)度為6 m的豎向懸臂梁對(duì)隔離數(shù)值子結(jié)構(gòu)方法進(jìn)行驗(yàn)證.梁柱模型如圖2所示.圖2(a)為主結(jié)構(gòu)模型，其包含2個(gè)梁柱單元，3個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)1采用固定約束.梁柱截面為復(fù)合截面，抗拉壓以及抗彎材料均選用應(yīng)力σ-應(yīng)變?chǔ)湃鐖D3所示的雙線性本構(gòu)模型，其中，抗拉壓剛度和抗彎剛度一樣，其彈性模量E均為5.74×106MPa，屈服點(diǎn)的應(yīng)力fy為1.47×104MPa，b為剛度比，此處取0.6.

在計(jì)算過(guò)程中，為了簡(jiǎn)化問(wèn)題，選取圖2(a)中的①單元為隔離單元(線彈性的主結(jié)構(gòu)模型，見(jiàn)圖2(a)右側(cè))，圖2(b)為精細(xì)化彈塑性子結(jié)構(gòu)模型.在主結(jié)構(gòu)的3號(hào)節(jié)點(diǎn)分別施加靜力、動(dòng)力荷載，并使用兩種隔離數(shù)值子結(jié)構(gòu)方法對(duì)該模型進(jìn)行響應(yīng)分析.

圖2 梁柱單元隔離數(shù)值子結(jié)構(gòu)方法的模型示意圖Fig.2Schematic diagram of the model for numerical substructures method of beam-column unit

圖3 雙線性本構(gòu)模型Fig.3Bilinear constitutive model

2.1 梁柱單元靜力荷載下彈塑性驗(yàn)證

圖4 梁柱單元標(biāo)準(zhǔn)模型Fig.4Standard model of beam-column unit

圖6 梁柱單元?jiǎng)恿奢d下節(jié)點(diǎn)4的位移響應(yīng)Fig.6Displacement response of node 4 under dynamic load of beam-column unit

圖5為梁柱單元的4號(hào)節(jié)點(diǎn)在靜力荷載作用下水平位移.由圖可知，采用迭代與不迭代算法的隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果與標(biāo)準(zhǔn)模型的解保持一致；當(dāng)荷載加載完成時(shí)(t=1.0 s)，單元①已進(jìn)入塑性階段，標(biāo)準(zhǔn)模型的4號(hào)節(jié)點(diǎn)的水平位移為3.879 35 cm，而采用迭代隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果也為3.879 35 cm，與標(biāo)準(zhǔn)解相同；采用不迭代隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果為3.866 88 cm，與標(biāo)準(zhǔn)解的相對(duì)誤差為0.3%.

圖5 梁柱單元靜力荷載下節(jié)點(diǎn)4的位移響應(yīng)Fig.5Displacement response of node 4 under static load of beam-column element

2.2 梁柱單元?jiǎng)恿奢d下彈塑性驗(yàn)證

當(dāng)對(duì)模型進(jìn)行動(dòng)力響應(yīng)分析時(shí)，由于需要考慮慣性項(xiàng)，所以對(duì)模型的精細(xì)化模擬需要對(duì)質(zhì)量進(jìn)行重新分配.本節(jié)中所采取的質(zhì)量矩陣為集中質(zhì)量矩陣，質(zhì)量分配的規(guī)則如圖2和圖4所示.節(jié)點(diǎn)4受到單點(diǎn)激勵(lì)，輸入的波形選用的是1989年美國(guó)洛馬-普雷塔的Loma Prieta地震波，放大倍數(shù)為300倍.

如圖6所示，對(duì)該模型進(jìn)行動(dòng)力分析時(shí)，靜力迭代算法的非線性修正力計(jì)算是基于子結(jié)構(gòu)靜力分析得到的，其計(jì)算結(jié)果與標(biāo)準(zhǔn)解的誤差為5%.將子結(jié)構(gòu)從靜力分析拓展至動(dòng)力分析后，迭代的隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果與標(biāo)準(zhǔn)解保持一致；不迭代算法的計(jì)算步長(zhǎng)為10-3s時(shí)，計(jì)算過(guò)程中的位移發(fā)散，計(jì)算無(wú)法完成；縮小計(jì)算時(shí)間步長(zhǎng)后，不迭代算法與標(biāo)準(zhǔn)解在彈性階段吻合較好，誤差保持在1%以內(nèi)，但是當(dāng)單元1進(jìn)入塑性階段，計(jì)算結(jié)果并不理想，誤差達(dá)到35%，此時(shí)，不迭代算法的計(jì)算時(shí)間步長(zhǎng)為10-4s，計(jì)算成本較大.為了加快不迭代算法的計(jì)算效率及計(jì)算精度，本文提出的預(yù)測(cè)改進(jìn)的不迭代隔離數(shù)值子結(jié)構(gòu)方法將時(shí)間步長(zhǎng)增大為10-3s，相比于普通的不迭代算法的計(jì)算效率提高了10倍，且計(jì)算結(jié)果與標(biāo)準(zhǔn)解的誤差僅為0.81%，極大地提高了不迭代算法在模擬結(jié)構(gòu)動(dòng)力響應(yīng)分析時(shí)的計(jì)算精度.

3 隔離數(shù)值子結(jié)構(gòu)方法的驗(yàn)證(平面四節(jié)點(diǎn)單元)

本節(jié)中的模型如圖7所示，圖7(a)為完全彈性的主結(jié)構(gòu)模型，取單元①作為隔離數(shù)值子結(jié)構(gòu)(如圖7(b))，并將其精細(xì)化，如圖7(c).主結(jié)構(gòu)共有8個(gè)節(jié)點(diǎn)，3個(gè)單元，每個(gè)節(jié)點(diǎn)2個(gè)自由度.結(jié)構(gòu)采用四節(jié)點(diǎn)單元，每個(gè)單元大小為3 m×3 m. 1和2節(jié)點(diǎn)為固定端.本算例為平面應(yīng)變問(wèn)題，單元厚度取0.1 m，分別使用兩種隔離數(shù)值子結(jié)構(gòu)方法對(duì)模型在靜力以及動(dòng)力荷載作用下的結(jié)構(gòu)響應(yīng)進(jìn)行分析.

圖7 四節(jié)點(diǎn)單元隔離數(shù)值子結(jié)構(gòu)方法的模型示意圖Fig.7Schematic diagram of the model for numerical substructures method of four-node unit

3.1 平面四節(jié)點(diǎn)單元靜力荷載下的彈塑性驗(yàn)證

圖8 四節(jié)點(diǎn)單元隔離數(shù)值子結(jié)構(gòu)方法的標(biāo)準(zhǔn)模型示意圖Fig.8Schematic diagram of the standard model for numerical substructure method of four-node unit

為了驗(yàn)證兩種隔離數(shù)值子結(jié)構(gòu)方法在平面四節(jié)點(diǎn)單元上的正確性，本節(jié)中參照隔離數(shù)值子結(jié)構(gòu)方法的模型網(wǎng)格尺寸進(jìn)行完整結(jié)構(gòu)的建模(見(jiàn)圖8)，并將其作為標(biāo)準(zhǔn)模型.由于在標(biāo)準(zhǔn)模型的求解過(guò)程中需要使節(jié)點(diǎn)9、10、12、13滿足一定的位移約束條件，即在求解等效動(dòng)力平衡方程時(shí)需增加一個(gè)約束矩陣T，則動(dòng)力方程為：

約束矩陣T可由節(jié)點(diǎn)9、10、12、13與相鄰節(jié)點(diǎn)的位置關(guān)系確定.例如在本節(jié)中，3、4、13號(hào)節(jié)點(diǎn)的位移滿足以下關(guān)系：

圖9 平面四節(jié)點(diǎn)單元靜力荷載下節(jié)點(diǎn)8的位移響應(yīng)Fig.9Displacement response of node 8 under static load of planar four-node element

由圖9可知，當(dāng)靜力荷載加載完成時(shí)(即t=1.0 s)，標(biāo)準(zhǔn)模型的8號(hào)節(jié)點(diǎn)水平位移的標(biāo)準(zhǔn)解為55.064 4 cm，而采用迭代的隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果為55.064 4 cm，與標(biāo)準(zhǔn)解相同；而不迭代算法的計(jì)算結(jié)果為53.649 7 cm，與標(biāo)準(zhǔn)解的相對(duì)誤差為2.6%.

3.2 平面四節(jié)點(diǎn)單元?jiǎng)恿奢d下的彈塑性驗(yàn)證

在動(dòng)力荷載下精細(xì)化的子結(jié)構(gòu)模型需要考慮慣性項(xiàng)，需要對(duì)節(jié)點(diǎn)質(zhì)量進(jìn)行重新分配，其分配規(guī)則如圖7和8所示.在標(biāo)準(zhǔn)模型中，5、6、7、8號(hào)節(jié)點(diǎn)的質(zhì)量均為50 kg，3和4號(hào)節(jié)點(diǎn)的質(zhì)量均為100/3 kg，10、11、12、13號(hào)節(jié)點(diǎn)的質(zhì)量均為25/3 kg.7號(hào)節(jié)點(diǎn)受到單點(diǎn)激勵(lì)，輸入的波形選用的是1989年美國(guó)洛馬-普雷塔的Loma Prieta地震波，放大倍數(shù)為1 600倍.

如圖10所示，在動(dòng)力荷載作用下，非線性修正力的計(jì)算是基于子結(jié)構(gòu)靜力分析得到的靜力迭代算法，其計(jì)算結(jié)果與標(biāo)準(zhǔn)解的相對(duì)誤差為7.61%.將子結(jié)構(gòu)從靜力分析拓展至動(dòng)力分析后，采用迭代的隔離數(shù)值子結(jié)構(gòu)方法的計(jì)算結(jié)果與標(biāo)準(zhǔn)解的計(jì)算結(jié)果基本保持一致.當(dāng)t=9.88 s時(shí)，標(biāo)準(zhǔn)模型中8號(hào)節(jié)點(diǎn)的水平位移為-40.593 4 cm，而采用迭代算法時(shí)，其水平位移為-40.625 3 cm，相對(duì)誤差為0.08%.采用普通的不迭代隔離數(shù)值子結(jié)構(gòu)方法，當(dāng)計(jì)算步長(zhǎng)為10-3s時(shí)，計(jì)算過(guò)程中的位移發(fā)散，計(jì)算無(wú)法完成；當(dāng)縮小計(jì)算步長(zhǎng)至10-5s，其計(jì)算結(jié)果為-40.731 8 cm，與標(biāo)準(zhǔn)解的相對(duì)誤差為0.34%；使用預(yù)測(cè)改進(jìn)的不迭代隔離數(shù)值子結(jié)構(gòu)方法后，其計(jì)算步長(zhǎng)可增大為10-3s，計(jì)算結(jié)果為-40.712 3 cm，計(jì)算效率提升了100倍，而相對(duì)誤差僅為0.29%.

4 結(jié) 論

本研究通過(guò)數(shù)值算例將改進(jìn)的數(shù)值子結(jié)構(gòu)方法與完整結(jié)構(gòu)的標(biāo)準(zhǔn)解進(jìn)行對(duì)比，系統(tǒng)驗(yàn)證了數(shù)值子結(jié)構(gòu)方法的計(jì)算精度和計(jì)算誤差，得到了如下結(jié)論：

1) 在結(jié)構(gòu)受到靜力荷載時(shí)，迭代隔離數(shù)值子結(jié)構(gòu)方法的位移響應(yīng)與標(biāo)準(zhǔn)模型的解始終保持一致，計(jì)算精度可達(dá)1 μm；而不迭代隔離數(shù)值子結(jié)構(gòu)方法在結(jié)構(gòu)處于彈性和弱非線性時(shí)與標(biāo)準(zhǔn)模型的解吻合得較好，但是當(dāng)結(jié)構(gòu)進(jìn)入強(qiáng)非線性階段時(shí)，誤差增大.

2) 在結(jié)構(gòu)受到動(dòng)力荷載時(shí)，迭代隔離數(shù)值子結(jié)構(gòu)方法的位移響應(yīng)與標(biāo)準(zhǔn)模型的解基本保持一致；不迭代隔離數(shù)值子結(jié)構(gòu)方法在結(jié)構(gòu)全部處于彈性和弱非線性階段時(shí)，吻合較好；當(dāng)子結(jié)構(gòu)進(jìn)入強(qiáng)非線性后，誤差逐漸增大，甚至不收斂.

3) 由于迭代算法計(jì)算成本較高，不迭代算法的計(jì)算結(jié)果在子結(jié)構(gòu)進(jìn)入強(qiáng)非線性后誤差較大，且步長(zhǎng)較小，效率較低.所以，本研究提出的隔離子結(jié)構(gòu)的響應(yīng)預(yù)測(cè)修正方法，能夠顯著增大不迭代隔離數(shù)值子結(jié)構(gòu)方法的步長(zhǎng)，在保證計(jì)算精度的同時(shí)，極大地提高了計(jì)算效率.