土體非線性對土－結構動力相互作用的影響
——以Millikan圖書館為例

梁建文，楊晶麗，韓 冰，朱 俊

（1.天津大學建筑工程學院，天津 300350；2.南京航空航天大學土木與機場工程系，江蘇南京 210016；3.北京工業(yè)大學建筑工程學院，北京 100124）

引言

地震激勵下，土體與結構之間存在土－結構動力相互作用（SSI）效應。文獻［1］采用間接邊界元法（IBEM）研究了SV波入射下，層狀半空間中剛性基礎－剪力墻結構體系的SSI，研究了場地動力特性對SSI的影響。文中在文獻［1］基礎上，對近場土體進行有限元法（FEM）模擬，遠場土體仍然采用間接邊界元法模擬，采用等效線性化方法考慮土體非線性，進一步研究土體非線性對SSI的影響。

等效線性化方法采用等效剪切模量和阻尼比隨土層動剪應變的變化近似考慮土體的非線性動力特性，具有形式簡單，計算效率高的優(yōu)點，是目前用來估計場地非線性的主要方法之一。Esteban 等［2］采用等效線性方法分析了土層的非線性動力特性對理想SDOF 結構模型的動力響應的影響，研究表明了當輸入地震波的主頻率接近考慮土體非線性后衰減的場地基頻值時，土體非線性會增加上部的結構損傷，而且相比于結構高度、基礎尺寸等參數，結構質量對SSI的影響更明顯。Massumi等［3］采用等效線性模型模擬有限元土域，分析了SSI對結構動力響應的影響，結果表明土體剛度較小時，SSI會使結構的撓度和層間位移角顯著增加；Quoc等［4］利用粘性土的等效線性骨架曲線考慮土體非線性，研究了非線性土－樁－結構動力相互作用的重要影響因素。Liu等［5］通過ABAQUS有限元軟件研究土－結構相互作用的重要影響，分別結合等效線性模型和彈塑性模型考慮土體非線性，結果表明2 類非線性數值模型的模擬結果均與振動臺試驗結果吻合。韓冰等［6］采用有限元－間接邊界元耦合法，通過等效線性化法考慮土層的非線性，研究了場地中存在透鏡體時對土－結構動力相互作用的影響。

文中采用有限元－間接邊界元（FEM－IBEM）耦合法，以著名的Millikan圖書館結構為例，計算El Centro波和Taft波從基巖面垂直入射時，不同地震波幅值，不同土體非線性參數下，結構基底剪力響應和頂部層間位移角響應，研究土體非線性對SSI的影響，得出一些有參考價值的結論。

1 模型與方法

采用FEM－IBEM 耦合法建立模型。采用有限元模擬近場土體，采用間接邊界元法（IBEM）模擬遠場土體，自動滿足無窮遠輻射條件。模擬近場土體的有限元（FEM）子域要足夠大從而保證IBEM 子域幾乎不受FEM子域散射的影響，使得遠場土體非線性可視作與自由場土體非線性等效［7］，這樣，SSI體系的等效線性迭代計算只需在有限元子域內進行，而將自由場非線性迭代后得到的土層特性參數作為IBEM 子域的計算輸入參數，避免重復建立半無限空間的整體動力剛度矩陣。在此前提下，確定較小的FEM 子域尺寸以提高計算效率。文中不考慮基礎與土體交界處可能發(fā)生的滑移等接觸問題。

1.1 模型

土－結構動力相互作用（SSI）體系的物理模型如圖1所示。

圖1 SSI模型Fig.1 SSI model

剪力墻高為H，寬為W，剛性基礎半徑為r，彈性基巖上覆蓋土層厚度為D。剪力墻與基礎之間以及基礎與土層之間都完全固接。地震波從基巖面垂直入射。

土層的特征參數包括土層質量密度ρL，阻尼比ζL和泊松比νL，以及拉梅常數λL和μL，土層剪切波速和壓縮波速分別為βL和αL。彈性基巖的特征參數包括基巖質量密度ρR，阻尼比ζR和泊松比νR，以及拉梅常數λR和μR，基巖的剪切波速和壓縮波速分別為βR和αR。

剪力墻的特征參數包括沿縱向單位長度的質量Μb，泊松比νb和阻尼比ζb，以及剪力墻基頻fb和等效剪切波速βb= 4fbH。基礎沿縱向單位長度的質量為M0。

地震作用下基礎的位移響應包括水平平動位移dx，豎向平動位移dz和關于圓心O的轉角θ。剪力墻沿水平方向和豎向的絕對位移分別為dbx和dbz，結構頂部相對于底部（或基礎）的相對位移，即剪力墻自身發(fā)生的彈性位移分別為drelx和drelz。

FEM－IBEM 耦合法模型如圖2 所示。其中，F(xiàn)EM 子域由近場土體組成，子域寬度2*m1，深度m2；IBEM子域由遠場土體和彈性基巖組成。將土層進一步劃分為n個子層，子層厚度dj（j=1，2，…，n）。IBEM 子域與FEM子域的界面稱為邊界S1，土－結構界面稱為邊界S2。FEM子域需要全部劃分有限單元，并將位于邊界S1和S2上的單元稱為耦合單元，其余為內部單元，而IBEM 子域和結構只需在各自邊界S1和S2處離散邊界單元即可。沿共同邊界S1，需要使IBEM 子域的邊界單元數與FEM 子域的耦合單元數相同，并且單元節(jié)點重合。沿共同邊界S2，F(xiàn)EM子域的耦合單元數與結構在S2上的邊界單元數相同，并且單元節(jié)點重疊。同時，邊界S1，S2上各對應單元響應之間應滿足以下條件：

圖2 FEM－IBEM耦合法計算模型示意圖Fig.2 FEM－IBEM model

式中，上標“－”和“+”分別代表邊界單元和耦合單元。單元位移u和應力t分別在各自子域獨立求解得到。另外，在邊界S1和S2上引入虛擬均布荷載P1，P2和P3，可看作是分別施加在IBEM 子域，F(xiàn)EM 子域和結構上的外荷載。而通過式（1）～式（3）建立的聯(lián)系最終可求解出荷載P1，P2和P3，這也是采用耦合法求解整個體系地震響應的核心步驟。

1.2 IBEM子域的求解方法

通常，在間接邊界元法中，遠場任意位置的響應可認為是該點自由場響應與散射場響應的疊加。首先，自由場響應u′（x，z）和t′（x，z）可通過直接剛度法計算得到：

式中，S為場地整體剛度矩陣，由每一子層的動力剛度矩陣Ej（j=1，2，…，n）和基巖剛度矩陣ER按一定順序集整而成；外荷載僅在基巖面由入射波引起，因此，荷載向量那么，通過方程（4）可解出各土層界面上的位移當第j層土的上下界面位移已知后，即確定了（z=0，z=dl）處的邊界條件，就可以在該土層的局部坐標系下求解土層中任意位置上的自由場位移響應u′(x,z)和應力響應t′(x,z)了。

求解散射場響應時，需要在邊界S1上離散N個邊界單元并引入格林影響函數為位移格林函數，其物理意義為在第j個邊界單元上依次施加單位水平均布荷載，單位豎向均布荷載時第i層土中任意點的位移響應。其求解過程為：先固定荷載作用層上下界面，根據該土層波數域中的運動方程得到上下界面節(jié)點位移，應力的特解和齊解；再取消固定，將節(jié)點反力作為外荷載施加于荷載作用層的上下界面，由式（4）求得場地中各土層的界面位移，那么任意i土層內任意位置的位移響應ui(k′,z)也相應可解（k′為散射場波數），注意荷載作用層的響應還應增加上一步的方程齊解和特解部分；最后對ui(k′,z)進行傅里葉變換得到第i層土在空間域中的位移響應即為其位移格林函數。同樣，場地的應力格林函數可由本構關系得到。

將格林函數與邊界上虛擬均布荷載P1相乘即為散射場響應，其中單元荷載向量可認為是FEM子域施加在IBEM子域上的荷載(j= 1,2,…,N)，上標“H”、“V”分別代表水平方向和豎向的荷載。

間接邊界元法的詳細計算方法見文獻［1］和［8］。

因此，根據間接邊界元法，IBEM子域在共同邊界S1上的邊界單元的平均位移響應和平均應力響應如下所示：

1.3 FEM子域的求解方法

近場有限單元法子域內的土體需要按四節(jié)點等參單元劃分網格。FEM 子域中沿S1邊界劃分的耦合單元數為N，沿S2邊界劃分的耦合單元數為M。根據有限單元法理論中平面應變問題的求解思路，集整所有土單元的單元剛度矩陣，建立FEM子域的整體動力剛度矩陣KF。KF可分區(qū)表示為：

因此地震作用下，由PF引起的FEM子域的動力響應為：

式中：A為將單元均布荷載轉化為單元等效節(jié)點荷載的轉換矩陣；uu為內部節(jié)點的位移向量；us為耦合邊界節(jié)點的位移向量。

1.4 結構響應的求解方法

結構頂部的位移響應為［9］：

結構的基底剪力Fbx為：

1.5 耦合方程

由式（5），式（8）和式（9）給出的沿S1、S2的單元位移和應力響應，式（1）～式（3）可具體化為：

由式（16），可將P3寫作用P2替代的關系式如下：

然后將式（17）代入式（14）和式（15），有

聯(lián)立式（18）和式（19）并寫成矩陣方程形式，即為FEM－IBEM耦合法的總耦合方程，如下：

由耦合方程式（20）最終可解得間接邊界元法子域的虛擬荷載向量P1和有限單元法子域的虛擬荷載向量P2，再由式（17）可得到虛擬荷載向量P3。最后，通過式（11）～式（13）可求解結構的地震響應。

1.6 等效線性化迭代計算

等效線性化的計算流程如下。根據土層的初始剪切模量G0和阻尼比ζ進行首次運算，將土單元經過前一次完整動力計算后的整個地震持時內的最大剪應變γmax，乘以等效系數0.65得到單元等效剪應變γeff；再分別根據等效剪切模量比、阻尼比與等效應變的關系曲線G/G0－γ，ζ－γ得到土單元的新剪切模量值Gi和新阻尼比值ζi，并將Gi和ζi作為開始本次動力迭代計算的輸入參數。若本次迭代計算后得到的單元等效剪應變與前一次的相差|γi-γi-1|/γi-1在某允許值σ范圍內，即認為該單元收斂［10］，否則，繼續(xù)求取新的Gi+1和ζi+1進入下一次動力計算。文中取σ=0.05，并以某次迭代有95%以上的單元收斂作為收斂條件。將最后一次動力迭代計算后得到的系統(tǒng)響應作為完成非線性計算的最終輸出結果。

2 方法驗證

圖3為文中耦合方法和文獻［1］的間接邊界元法的計算結果對比。基礎寬度r=10 m，土層厚度D=30 m，40 m，土層和基巖參數為ρR=ρL=1 900 kg/m3，ζL=0.05，ζR=0.02，νL=νR= 1 3，βL=200 m/s，βR=400 m/s，1 000 m/s。上部結構H=20 m，W=20 m，βb=200 m/s，ε=βLH βbr=2，ζb=0，νb= 1 3，MbM0=2，M0MS=1。以基巖面垂直入射SV 波時的基礎位移和結構頂部的相對位移為例進行了對比驗證。可以看出，文中計算結果與文獻［1］的結果相吻合，證明了文中方法的有效性。

圖3 垂直入射SV波時基礎位移響應和頂部相對位移響應與文獻［1］的結果對比Fig.3 Comparison of the foundation displacement and the relative displacement at the top of the structure for vertical incidence of SV-wave with the results in the paper［1］

圖3 （續(xù)）Fig.3 （Continued）

3 算例分析

本節(jié)以著名的Millikan圖書館為例進行算例分析，具體參數見文獻［11－12］。

圖書館結構參數H=44 m，W=25 m，r=12.5 m，fb=2.33 Hz，νb= 1 3，ζb=0.01，MbM0=7.5，M0MS=0.2。圖書館場地參數如表1 所示。沿彈性基巖面垂直入射的地震波分別選用幅值為0.1 g，0.2 g的El Centro 波和Taft 波，波型為SV 波。地震波的加速度時程曲線如圖4（a）和（b）所示。對上述4 個工況分別進行線性和非線性分析。其中，非線性分析時各工況下又分開采用了2類非線性參數，一次計算過程中所有土層均使用同一組土動力參數。2組土層動剪切模量、阻尼比與剪應變的關系曲線如圖5所示，取自文獻［13］。

圖4 基巖輸入地震波的加速度時程曲線Fig.4 Acceleration time-history curves of the seismic waves

圖5 兩組等效線性G/G0－γ和ζ－γ關系曲線Fig.5 Two types of the relation curves of G/G0－γ and ζ－γ

表1 Millikan圖書館場地參數［11］Table 1 Parameters of the Millikan library building site［11］

為了方便描述考慮土體非線性后非線性SSI響應與線性SSI響應之間的差異，定義2個參數如下：

式中：χF表示考慮土體非線性后的結構基底剪力響應較線性響應的變化幅度；χD表示考慮土體非線性后的結構頂部層間位移角響應較線性響應的變化幅度（Db=drelx/H）。

圖6 和圖7 分別給出了Taft 波和El Centro 波入射下，結構基底剪力Fbx非線性響應與線性響應的比較。由圖可見，在考慮土體非線性后，F(xiàn)bx時程響應峰值均有不同程度的減少。當基巖輸入幅值為0.1 g 的Taft 波時，采用兩類非線性參數（Ⅰ和Ⅱ）進行等效線性計算后的基底剪力分別比線性情況時降低了|χF|Ⅰ=8.4%，|χF|Ⅱ=20.54%；當基巖輸入幅值為0.2 g 的Taft 波 時，|χF|Ⅰ=24.18%，|χF|Ⅱ=43.08%。當基巖輸入0.1 g 幅值的El Centro 波時，F(xiàn)bx非線性響應的降低幅度分別為|χF|Ⅰ=22.13%，|χF|Ⅱ=31.75%；當基巖輸入0.2 g 幅值的El Centro 波時，|χF|Ⅰ=44.53%，|χF|Ⅱ=54.28%。顯然，采用非線性更強的Ⅱ類參數計算后的基底剪力下降幅度|χF|Ⅱ均大于Ⅰ類情況時的下降幅度|χF|Ⅰ；而且地震波幅值為0.2 g時的基底剪力下降幅度均比0.1 g時更明顯。

圖6 基巖輸入Taft波情況結構基底剪力Fbx時程曲線Fig.6 Time-history curves of the base shear for Taft wave input from bedrock

圖7 基巖輸入El Centro波情況結構基底剪力Fbx時程曲線Fig.7 Time-history curves of the base shear for El Centro wave input from bedrock

同樣，圖8 和圖9 分別給出了Taft波和El Centro 波入射下，結構頂部層間位移角Db的非線性響應與線性響應的比較。由圖可見，在考慮土體非線性后，層間位移角Db的時程響應峰值均有不同程度的減少。當基巖輸入幅值為0.1 g 的Taft 波時，層間位移角Db的非線性響應降低幅度分別為|χD|Ⅰ=6.65%，|χD|Ⅱ=15.71%；當基巖輸入幅值為0.2 g 的Taft 波時，|χD|Ⅰ=17.22%，|χD|Ⅱ=34.89%。當基巖輸入0.1 g 幅值的El Centro 波時，Db非線性響應的降低幅度分別為|χD|Ⅰ=14.32%，|χD|Ⅱ=24.95%；當基巖輸入0.2 g 幅值的El Centro 波時，|χD|Ⅰ=40.93%，|χD|Ⅱ=53.31%。顯然，對于結構的層間位移角響應，非線性響應規(guī)律與基底剪力類似，土層的非線性越強，地震波幅值越大，非線性響應越明顯。

圖8 基巖輸入Taft波情況層間位移角時程曲線Fig.8 Time-history curves of the structure drift for Taft wave input from bedrock

圖9 基巖輸入El Centro波情況層間位移角時程曲線Fig.9 Time-history curves of the structure drift for El Centro wave input from bedrock

圖10和圖11分別給出了Taft波和El Centro波入射下，土體非線性（I和II）和線性情況下結構頂部絕對加速度Ab的傅里葉譜比較。由圖可見，考慮土體非線性后，系統(tǒng)頻率明顯降低。當基巖輸入幅值為0.1 g的Taft波時，SSI系統(tǒng)頻率由2.03 Hz降到了1.54 Hz；基巖輸入幅值為0.2 g的Taft波時，SSI系統(tǒng)頻率由2.03 Hz降到了1.17 Hz。當基巖輸入0.1 g幅值的El Centro波時，考慮土體非線性后，SSI系統(tǒng)頻率由1.95 Hz降到了1.46 Hz；輸入0.2 g幅值的El Centro波時，SSI系統(tǒng)頻率由1.95 Hz降到了1.17 Hz。顯然，地震波幅值越大，系統(tǒng)頻率降低幅度越大。

圖10 基巖輸入Taft波時結構頂部絕對加速度時程曲線及其傅里葉譜Fig.10 Time-history curves of the absolute acceleration at the top of the structure and their Fourier spectra for Taft wave

圖10 （續(xù)）Fig.10 （Continued）

圖11 基巖輸入El Centro波時結構頂部絕對加速度時程曲線及其傅里葉譜Fig.11 Time-history curves of the absolute acceleration at the top of the structure and their Fourier spectra for El Centro wave

4 結論

文中建立了層狀半空間中SSI體系的二維FEM－IBEM 耦合法模型，采用等效線性化方法同時考慮近場和遠場土體的非線性，以著名的Millikan圖書館為例，計算了Taft波和El Centro波垂直入射情況下，分別采用2 類非線性參數時的結構基底剪力和層間位移角。最后通過與線性SSI 響應比較，研究了土體非線性對SSI的影響。主要結論如下：

（1）土體非線性對SSI 有明顯影響：考慮土體非線性后，結構基底剪力響應和層間位移角均比線性情況時有明顯降低；結構系統(tǒng)頻率均明顯低于線性情況。文中算例中，基巖輸入Taft波時基底剪力最小和最大降幅分別為8.4%和43%；層間位移角最小和最大降幅分別為6.7%和35%。基巖輸入El Centro 波時基底剪力的最小和最大降幅分別為22%和54%；層間位移角最小和最大降幅分別為14.3%和53.3%。

（2）土體的非線性越強，對非線性SSI 的影響越明顯。文中算例中，采用非線性更強的Ⅱ類土體參數時的非線性SSI 響應衰減幅度均比采用Ⅰ類參數時更明顯。基巖輸入幅值為0.1 g 的Taft 波時，采用兩類非線性參數（Ⅰ和Ⅱ）的基底剪力的降幅分別為8.4%和20.5%，層間位移角的降幅分別為6.7%和15.7%；基巖輸入幅值為0.2 g的Taft波時，基底剪力的降幅分別為24%和43%，層間位移角的降幅分別為17%和35%；基巖輸入幅值為0.1g 的El Centro 波時，采用兩類非線性參數（Ⅰ和Ⅱ）的基底剪力的降幅分別為22%和32%，層間位移角的降幅分別為14.3%和25%；基巖輸入幅值為0.2 g的El Centro波時，基底剪力的降幅分別為44.5%和54%，層間位移角的降幅分別為41%和53.3%。

（3）地震波入射幅值對非線性SSI 有明顯影響，且地震波幅值越大，非線性SSI 越明顯。文中算例中，基巖輸入幅值為0.2 g 的Taft 波或El Centro 波時，結構基底剪力，層間位移角以及系統(tǒng)頻率的非線性衰減幅度均比幅值為0.1 g時更明顯。基巖輸入Taft波（0.1 g和0.2 g）時系統(tǒng)頻率的降幅分別為24%和42%；基巖輸入El Centro波（0.1 g和0.2 g）時系統(tǒng)頻率的降幅分別為25%和40%。