大型結構局部非線性問題的數值子結構方法的驗證與改進

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

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

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

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

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

1 數值子結構方法及其改進方法

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

基于有限元的整體結構運動方程為[11]：

(1)

1.1 基于子結構靜力分析法的非線性修正力

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

(2)

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

(3)

(4)

(5)

(6)

如圖1所示，ub為子結構邊界上的節點位移.假設子結構邊界上節點位移是連續的，ub可根據式(7)插值得到：

ub=NuNL.

(7)

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

δub=NδuNL.

(8)

根據虛功原理，主結構隔離系統與子結構系統的虛功是相等的，即

(9)

(10)

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

RNL=NTRsub,b.

(11)

由此可得隔離區域的非線性修正力在子結構進行精細化模擬時的計算公式：

(12)

圖1 精細化子結構模型的非線性修正力計算方法Fig.1Calculation method of nonlinear correction force for refined substructure model

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

因此，以下介紹子結構動力分析時的非線性修正力的公式推導.

1.2 基于子結構動力分析法的非線性修正力

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

(13)

(14)

(15)

其中

(16)

而等效動力剛度和改進的非線性修正力分別為：

(17)

(18)

(19)

根據虛功原理，主結構隔離系統與子結構系統的虛功是相等的，即

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

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

(20)

其中：MNL為主結構隔離區域的節點質量，Msub為子結構系統的節點質量，Msub,b、Msub,i分別為子結構系統的邊界節點和內部節點質量.由于內部節點不受外力作用，其內部節點的結構抗力Rsub,i與慣性力Msub,iüsub,i的和始終為0.因此，同理可得

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

(21)

將式(21)代入(18)，可得子結構動力分析時隔離區域的非線性修正力的計算公式：

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

(22)

1.3 算法流程簡介

本研究分別采用不迭代與迭代算法對隔離數值子結構中的動力問題進行分析.其中不迭代算法的流程可分為以下幾個步驟：

1.4 預測改進的隔離數值子結構方法

(23)

預測改進的隔離數值子結構方法可分為以下幾個步驟：

2 隔離數值子結構方法的驗證(梁柱單元)

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

在計算過程中，為了簡化問題，選取圖2(a)中的①單元為隔離單元(線彈性的主結構模型，見圖2(a)右側)，圖2(b)為精細化彈塑性子結構模型.在主結構的3號節點分別施加靜力、動力荷載，并使用兩種隔離數值子結構方法對該模型進行響應分析.

圖2 梁柱單元隔離數值子結構方法的模型示意圖Fig.2Schematic diagram of the model for numerical substructures method of beam-column unit

圖3 雙線性本構模型Fig.3Bilinear constitutive model

2.1 梁柱單元靜力荷載下彈塑性驗證

圖4 梁柱單元標準模型Fig.4Standard model of beam-column unit

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

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

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

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

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

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

3 隔離數值子結構方法的驗證(平面四節點單元)

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

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

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

圖8 四節點單元隔離數值子結構方法的標準模型示意圖Fig.8Schematic diagram of the standard model for numerical substructure method of four-node unit

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

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

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

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

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

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

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

4 結 論

本研究通過數值算例將改進的數值子結構方法與完整結構的標準解進行對比，系統驗證了數值子結構方法的計算精度和計算誤差，得到了如下結論：

1) 在結構受到靜力荷載時，迭代隔離數值子結構方法的位移響應與標準模型的解始終保持一致，計算精度可達1 μm；而不迭代隔離數值子結構方法在結構處于彈性和弱非線性時與標準模型的解吻合得較好，但是當結構進入強非線性階段時，誤差增大.

2) 在結構受到動力荷載時，迭代隔離數值子結構方法的位移響應與標準模型的解基本保持一致；不迭代隔離數值子結構方法在結構全部處于彈性和弱非線性階段時，吻合較好；當子結構進入強非線性后，誤差逐漸增大，甚至不收斂.

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