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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

(4)

(5)

(6)

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

ub=NuNL.

(7)

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

δ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ì)化模擬時的計算公式：

(12)

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

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

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

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

動力分析時把運動方程對時間t進(jìn)行離散，每一時步都滿足以下方程：

(13)

(14)

(15)

其中

(16)

而等效動力剛度和改進(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é)點質(zhì)量，Msub為子結(jié)構(gòu)系統(tǒng)的節(jié)點質(zhì)量，Msub,b、Msub,i分別為子結(jié)構(gòu)系統(tǒng)的邊界節(jié)點和內(nèi)部節(jié)點質(zhì)量.由于內(nèi)部節(jié)點不受外力作用，其內(nèi)部節(jié)點的結(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)動力分析時隔離區(qū)域的非線性修正力的計算公式：

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

(22)

1.3 算法流程簡介

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

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

(23)

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

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

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

在計算過程中，為了簡化問題，選取圖2(a)中的①單元為隔離單元(線彈性的主結(jié)構(gòu)模型，見圖2(a)右側(cè))，圖2(b)為精細(xì)化彈塑性子結(jié)構(gòu)模型.在主結(jié)構(gòu)的3號節(jié)點分別施加靜力、動力荷載，并使用兩種隔離數(shù)值子結(jié)構(gòu)方法對該模型進(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 梁柱單元靜力荷載下彈塑性驗證

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

圖6 梁柱單元動力荷載下節(jié)點4的位移響應(yīng)Fig.6Displacement response of node 4 under dynamic load of beam-column unit

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

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

2.2 梁柱單元動力荷載下彈塑性驗證

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

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

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

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

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

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

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

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

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

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

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

3.2 平面四節(jié)點單元動力荷載下的彈塑性驗證

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

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

4 結(jié) 論

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

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

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

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