結(jié)構(gòu)-地基系統(tǒng)高階時域模型的動輸入方法

周文凱 陳燈紅

(1.三峽大學(xué) 土木與建筑學(xué)院,湖北 宜昌 443002;2.三峽大學(xué) 期刊社,湖北 宜昌 443002)

結(jié)構(gòu)-地基動力相互作用是結(jié)構(gòu)動力分析及結(jié)構(gòu)抗震計算中的一個重要方面.研究結(jié)構(gòu)-地基動力相互作用問題時,為考慮無限地基的輻射阻尼效應(yīng),最理想的做法是將地基范圍取得足夠大.由于動力分析中,網(wǎng)格尺寸通常受到最小波長的限制,這樣勢必會要求很大的計算機存儲量和工作量,實現(xiàn)時代價很大.因此,通常情況下在無限的地基中截斷一定范圍,這樣該系統(tǒng)分為近場的有限域和遠場的無限域部分.有限域一般采用有限元法考慮,如何準確地模擬無限域的輻射邊界條件是研究的難點.基于各種計算假設(shè)和不同的數(shù)學(xué)、力學(xué)模型,許多學(xué)者建立了多種人工邊界條件,它們可分為兩類,即局部的近似邊界和全局的精確邊界.近似邊界包括黏性邊界(viscous boundary)[1]、黏彈性邊界(viscous-spring boundary)[2-3]、多次透射公式(multi-transmitting formula)[4]、無限元(infinite element)[5]等;精確方法包括邊界元法(BEM)[6]、比例邊界有限元法(SBFEM)[7-8]等.

在動輸入方面,一般而言,有作為封閉系統(tǒng)的振動體系和作為開放系統(tǒng)的波動體系兩種方式.對于封閉系統(tǒng),將結(jié)構(gòu)的慣性力作為外力輸入,目前一般工業(yè)和民用建筑物多采用這種輸入方式.然而,對于尺寸、質(zhì)量很大,空間效應(yīng)顯著的高壩結(jié)構(gòu),傳統(tǒng)的無質(zhì)量地基假定不符合實際[9],地基的慣性和阻尼不能忽略,其動力學(xué)方程必須作為開放系統(tǒng)的波動問題來求解.

比例邊界有限元法是由 Wolf、Song于20世紀90年代提出的一種半解析數(shù)值方法.該方法在發(fā)展之初是為了解決無限域中的波動傳播問題,并具有顯著的優(yōu)勢.近年來,Song、Bazyar提出分別采用Pade近似[10]和連分式方法[11]來求解無限域的動力剛度矩陣,這兩種近似算法具有收斂范圍廣、收斂速度快等優(yōu)點,在不影響計算精度前提下其計算效率比早期的卷積積分算法有明顯提高.Birk等[12]針對多自由度系統(tǒng)高階連分式展開時可能造成剛度矩陣病態(tài)的問題,改進了文獻[11]中的連分式算法,算例表明改進后的算法明顯提高了數(shù)值穩(wěn)定性和魯棒性.陳燈紅等[13]在文獻[12]的基礎(chǔ)上發(fā)展了一種大壩-地基動力相互作用的高階時域模型.

本文在文獻[13]的研究成果基礎(chǔ)上,在時域里建立了一種結(jié)構(gòu)-地基相互作用的高階模型,借鑒黏彈性人工邊界等效力的動輸入方法[14-15],建立了與上述高階時域算法相匹配的波動輸入模型,并通過兩個數(shù)值算例驗證了文中建立方法的有效性.

1 SBFEM位移控制方程

基于彈性動力學(xué)基本方程,采用比例邊界坐標(ξ為徑向,η為環(huán)向)變換,文獻[8]詳細地推導(dǎo)出了SBFEM位移控制方程為:

其中,s為空間維數(shù),s=2或3分別表示二維或三維情況;u(ξ)為待求的位移向量.二維情況下,4個系數(shù)矩陣E0、E1、E2、M0分別表示為

式中,B1、B2為應(yīng)變-位移矩陣;D為彈性矩陣;J為ξ=1邊界上的Jacobian矩陣;N(η)為形函數(shù);ρ為質(zhì)量密度.

2 高階時域算法

2.1 近場有限域

結(jié)構(gòu)-地基系統(tǒng)的動力相互作用示意如圖1所示,近場有限域動力學(xué)方程可表示為

式中,M、C、K分別表示有限域的質(zhì)量矩陣、阻尼矩陣和剛度矩陣;下標s、b分別為結(jié)構(gòu)除交界面以外的自由度和結(jié)構(gòu)-地基交界面上的自由度;Ps、Pb、Rb分別為結(jié)構(gòu)所受外力、地基所受外力、結(jié)構(gòu)-地基系統(tǒng)的相互作用力.

圖1 土-結(jié)構(gòu)相互作用系統(tǒng)

2.2 結(jié)構(gòu)-地基耦合系統(tǒng)的運動方程

將以u(ξ)為未知量的二階偏微分方程式(1)轉(zhuǎn)化為以無限域動力剛度矩陣S∞(ω)為未知量的比例邊界有限元控制方程,然后將S∞(ω)采用連分式展開求解[13,16],根據(jù)邊界ξ=1上的力-位移平衡條件,建立了無限域時域分析的邊界條件,即

將近場有限域和遠場無限域通過結(jié)構(gòu)-地基交界面上的相互作用力向量Rb(如圖1)進行耦合[13,17].對于線性系統(tǒng),聯(lián)合式(3)、(4),并按照自由度順序?qū)?yīng)相加并展開,得到結(jié)構(gòu)-地基耦合系統(tǒng)的動力學(xué)方程為:

其中,耦合系統(tǒng)的剛度、阻尼、質(zhì)量矩陣Kc、Cc、Mc,以及待求向量dc、荷載向量fc分別表示為

從上式可以看出,耦合系統(tǒng)的剛度矩陣由K∞及高階項(階數(shù)i=1,2,…,M H)組成,阻尼矩陣由C∞及高階項組成,該式即為高階邊界條件.若上式中i=0,則邊界條件退化為剛度為K∞、阻尼為C∞的彈簧-阻尼器單元,即黏彈性邊界,其只有一階精度.

3 動輸入方法

由式(4)可看出,該邊界實質(zhì)上是由若干彈簧-阻尼器并聯(lián)組成的高階邊界條件,其性質(zhì)與黏彈性人工邊界相似,能夠吸收由近場有限域向遠場無限域散射的外行波.針對外源波動輸入問題時,需要在人工邊界處施加自由場荷載,則人工邊界上節(jié)點的等效荷載為

式中,kb、cb分別為式(4)中表示無限域剛度矩陣Ku、阻尼矩陣Cu的子矩陣;分別為自由場位移和速度矢量;Pf為自由場應(yīng)力矢量,可根據(jù)下式計算得到

其中,l、m為邊界外法線的方向余弦;應(yīng)力σ可由彈性理論求得,即

式中,λ、G為 Lame常數(shù),并有ρ為介質(zhì)的質(zhì)量密度,cs、cp分別為剪切波、膨脹波的波速.

當(dāng)從底部邊界垂直入射P波時,底部的邊界條件為u=0,ν=ν0(t),則在高度h處,自由場位移、速度可表示為入射波和反射波的疊加,即

式中,H為底部邊界到自由表面的距離.

對于底部邊界,其邊界條件為h=0,{l,m}T={0,-1}T,將式(8)～(10)代入式(7)中,等效荷載Pb可表示為

當(dāng)從底邊界垂直入射SV波時,底部的邊界條件為u=u0(t),ν=0,則等效荷載Pb計算為

對于左、右兩側(cè)邊界的等效荷載Pb也可按類似的方法求出.

4 數(shù)值算例

4.1 算例1

本文涉及到高階時域算法和波動輸入方法兩個關(guān)鍵問題,首先為了驗證高階時域算法的正確性,考慮如圖2所示的相互作用問題,取e=b=1.近場有限域、遠場無限域地基的材料參數(shù)分別取為G1=4G=4,ρ1=ρ=1,ν1=0.25;G2=G=1,ρ2=ρ=1,ν2=0.25.施加在頂部自由面上的均布荷載P(t)如圖3所示.

圖2 耦合系統(tǒng)

圖3 施加的均布荷載

本文算法中,耦合系統(tǒng)采用四結(jié)點四邊形單元離散,劃分72個四結(jié)點單元,共91個結(jié)點,其中截斷邊界采用24個兩結(jié)點線單元離散,相似中心選在O點,如圖4所示.計算總時間為20b/cs,積分步長為t=0.02b/cs.

圖4 耦合系統(tǒng)的網(wǎng)格圖

為了說明本文算法的精確性,采用有限元擴展網(wǎng)格進行了對比分析.擴展網(wǎng)格的計算區(qū)域大小取為對稱的20b×20b,采用八結(jié)點四邊形單元離散,劃分了6 400個單元共19 521個結(jié)點,該模型中截斷邊界采用固定邊界條件.

為了對比說明本文算法的優(yōu)勢,同時進行了有限元法結(jié)合黏彈性邊界對該問題的分析[18].計算區(qū)域取為4b×2b,有限元網(wǎng)格采用八結(jié)點四邊形單元離散,劃分了288個單元共937個結(jié)點.

圖5為觀測點O、A點的無量綱豎向位移時程.由圖5可知,本方法(M H=6)的結(jié)果與有限元擴展網(wǎng)格的結(jié)果在無量綱時間5～10 s之間略有偏移;當(dāng)M H增大到8時,兩者吻合得很好;而黏彈性邊界在無量綱時間5～10 s之間偏移較大,需要進一步擴大計算范圍才能提高計算精度.

圖5 O、A點豎向位移時程

4.2 算例2

其次,為驗證上述波動輸入法的正確性,采用底部垂直入射SV波的均勻彈性半空間進行驗證.介質(zhì)的力學(xué)參數(shù)為:彈模E=13 230 MPa,泊松比ν=0.25,密度ρ=2 700 kg/m3,波速cs=1 400 m/s,cp=2 425 m/s,入射波的位移方程為:

計算區(qū)域的范圍為762 m×381 m,有限元單元尺寸為15.4 m×15.4 m,單元數(shù)為450,其中截斷邊界采用60個兩結(jié)點線單元離散,結(jié)點數(shù)為496.計算總時間為2 s,時間步長取為0.01 s.經(jīng)過計算得到系統(tǒng)頂部和底部各點的水平向位移時程,如圖6所示.

圖6 水平向位移時程

可以看出:頂部自由表面的最大值(t=0.544 s)為入射SV波幅值的2倍;底部位移時程的前半段(t≤0.544 s)是入射波引起的位移,后半段(0.544 s≤t≤1.089 s)是反射波引起的位移.上述數(shù)值解與理論解吻合良好,從而說明本文建立的高階算法及其動輸入方法是正確的、精確的.

5 結(jié) 語

采用有限元法模擬有限域,采用基于SBFEM的高階透射邊界模擬遠場無限域,構(gòu)建了一種結(jié)構(gòu)-地基動力相互作用的高階時域模型及相適配的動輸入機制.借鑒黏彈性人工邊界中等效力的動輸入機制,將截斷邊界上的總波場分解為散射場和自由場,其中局部場地效應(yīng)引起的散射場由高階人工邊界條件自動滿足,無局部場地效應(yīng)影響的自由場由彈性理論求出.通過兩個數(shù)值算例驗證了建立的高階時域模型及相應(yīng)的動輸入機制的正確性、精確性,為復(fù)雜結(jié)構(gòu)的地震動輸入奠定基礎(chǔ).