分析土和結(jié)構(gòu)相互作用的一種實用耦合方法①

E-mail:quangu@xmu.edu.cn。

古泉， 彭伊, 曾志弘

(廈門大學(xué)建筑與土木工程學(xué)院,福建 廈門 361005)

摘要：提出一種新的數(shù)值解與解析解耦合的理論和計算方法，研究土-結(jié)構(gòu)相互作用(SSI)體系的地震動力響應(yīng)。采用大型有限元軟件OpenSees模擬復(fù)雜結(jié)構(gòu)的非線性行為，用等效線彈性頻域內(nèi)解析解模擬地基土的行為，使用時域離散遞歸方法將頻域內(nèi)的解析解轉(zhuǎn)化到時域內(nèi)，再通過子結(jié)構(gòu)邊界上力和位移的協(xié)調(diào)條件來求解。二者之間的耦合和實時數(shù)據(jù)交流通過CS集成方法來實現(xiàn)。以一個單自由度算例和一個實際工程為例，驗證此方法的精度、穩(wěn)定性和工程實用性，對比在考慮和不考慮SSI體系情況下結(jié)構(gòu)動力響應(yīng)的區(qū)別。本文所提的耦合SSI計算方法和部分研究成果可為工程設(shè)計人員提供參考。

關(guān)鍵詞：土-結(jié)構(gòu)相互作用； 子結(jié)構(gòu)法； CS方法； 數(shù)值解與解析解耦合方法； OpenSees; 非線性地震動力響應(yīng)分析

收稿日期：①2014-08-20

基金項目：國家自然科學(xué)基金重大國際(中美)合作項目(51261120376)；國家自然科學(xué)基金重大研究計劃集成項目(91315301-12)；國家自然科學(xué)基金面上項目(51578473)

作者簡介：古泉(1974-)，男，福建廈門人，博士，教授，主要從事高層建筑結(jié)構(gòu)與土結(jié)體系的非線性地震動力分析研究。

中圖分類號:TU352.12; U441.3文獻標(biāo)志碼：A

DOI:10.3969/j.issn.1000-0844.2015.03.0845

A Practical Coupling Method for Analyzing Soil-structure Interaction

GU Quan, PENG Yi, ZENG Zhi-hong

(SchoolofArchitectureandCivilEngineering,XiamenUniversity,Xiamen361005,Fujian,China)

Abstract：A novel practical method is presented for the analysis of soil-structure interaction (SSI). In this method, the structure is modeled by nonlinear FEM (OpenSees), and the soil is modeled by a time domain solution that is transformed from a frequency domain analytical solution using a discrete time domain recursive filter. The boundary conditions of force and displacement between soil and structure are satisfied using Newton’s method, and the coupling between the two substructures is based on CS integration techniques. We use structure and soil systems with a single degree of freedom, and study a real SSI example to determine the efficiency, accuracy, and applicability of the proposed method. Furthermore, we examine the differences between the conditions when considering and not considering SSI effects. This study proposes a practical method for nonlinear seismic analysis of SSI systems, and the research results provide valuable insights for engineering applications.

Key words： soil-structure interaction; substructure method; CS method; method of coupling numerical and analytical solutions; OpenSees; nonlinear seismic dynamic response analysis

0引言

在地震作用下土與結(jié)構(gòu)相互作用(SSI)能改變系統(tǒng)的固有頻率和阻尼等[1]，從而對結(jié)構(gòu)響應(yīng)產(chǎn)生重要的影響[2]?；趯嶋H工程中的觀測，即使采用相同的上部結(jié)構(gòu)，不同地區(qū)的土質(zhì)條件差異仍會使結(jié)構(gòu)的動力響應(yīng)存在很大差別[3-11]。雖然SSI普遍存在，但無法將SSI對結(jié)構(gòu)的作用得出一個一般性的結(jié)論[12-16]。在復(fù)雜的重大基礎(chǔ)工程中往往需要考慮SSI的影響，所以研究實用SSI體系分析方法，探究其對結(jié)構(gòu)安全性和運行性能的影響是很有必要的。

現(xiàn)階段，處理土-結(jié)構(gòu)相互作用問題的常用方法主要有有限元整體分析方法和子結(jié)構(gòu)法[12-16]。有限元法對于大型模型中土體的模擬往往會涉及到較多自由度和大型剛度矩陣的處理，其復(fù)雜的運算可能會導(dǎo)致冗長的計算時間并產(chǎn)生大量的計算成本。子結(jié)構(gòu)法是縮減自由度數(shù)目的一種有效方法，將整體SSI體系分為兩個或者多個子結(jié)構(gòu)來處理，對縮減后的運動方程進行求解，但各個子結(jié)構(gòu)之間的循環(huán)迭代和數(shù)據(jù)交流使得此方法應(yīng)用于非線性結(jié)構(gòu)體系時會受到一定局限。

本文在子結(jié)構(gòu)法的基礎(chǔ)上，提出數(shù)值解與解析解耦合的分析方法。該方法結(jié)合有限元分析和基于時域離散遞歸方法的頻域內(nèi)解析解，既能滿足復(fù)雜的非線性結(jié)構(gòu)體系的精確分析，又能實現(xiàn)土體的簡單快速模擬，對于大規(guī)模的線性或者非線性土-結(jié)構(gòu)相互作用體系的數(shù)值分析和模擬具有較大的應(yīng)用潛力。

1基于子結(jié)構(gòu)法的耦合方法

假定土-結(jié)構(gòu)相互作用體系由線彈性半無限大土、剛性基礎(chǔ)、上部結(jié)構(gòu)組成，結(jié)構(gòu)可以進入非線性。上部結(jié)構(gòu)用有限元軟件OpenSees求解[17]，土的頻域內(nèi)解析解通過時域離散遞歸方法可轉(zhuǎn)化到時域[16]，通過整個體系不同子結(jié)構(gòu)邊界上力和位移的協(xié)調(diào)條件來求解，并用一種高效率和實用的集成技術(shù)(即CS技術(shù))來實現(xiàn)有限元軟件和土-結(jié)構(gòu)相互作用框架的耦合。

1.1土-結(jié)構(gòu)相互作用體系頻域解析解的時域表示

剛性基礎(chǔ)的總位移為：

其中ug為自由場的輸入位移；us為剛性基礎(chǔ)施加在土上的廣義力和力矩所產(chǎn)生的額外位移。位移us可由下式計算得到：

其中C(ω)為土中剛性基礎(chǔ)的柔度矩陣，該矩陣為頻率的函數(shù)，與基礎(chǔ)的幾何特性和土體的性質(zhì)有關(guān)[4]，柔度矩陣求逆即為頻域內(nèi)剛度矩陣(即C(ω)=K(ω)-1)。對于不同的土壤類型和剛性基礎(chǔ)，可在文獻中找到該剛度矩陣的計算方法[4,13]。Fs為剛性基礎(chǔ)作用在土體上的廣義力，可由下式計算得到：

M0為剛性基礎(chǔ)的質(zhì)量矩陣;Fb為上部結(jié)構(gòu)作用在剛性基礎(chǔ)的廣義力，基于本文所提出的耦合方法，時域內(nèi)的廣義力矩陣Fb可由有限元軟件計算得到，即通過OpenSees分析得到。

為了計算土-結(jié)構(gòu)相互作用體系的動力時程反應(yīng)，需要對系統(tǒng)在時域內(nèi)進行分析，因此需將C(ω)在頻域內(nèi)的解析解轉(zhuǎn)化到時域內(nèi)。本文將文獻中已有的單變量的離散遞歸方法[16]擴展到多變量的時域離散遞歸方法，通過以下例子來說明該方法。

在一個二維的均勻半無限地基中，正方形剛性地基的柔度矩陣在時域內(nèi)表示為：

其中a為基礎(chǔ)寬度的一半；uB、θB分別是基礎(chǔ)的水平位移和轉(zhuǎn)角；kij、Cij(i, j=H,M分別代表水平和轉(zhuǎn)動位移)分別為頻域內(nèi)的剛度和柔度函數(shù)；Hs、Ms分別為作用在土上的水平力和彎矩。通過時域離散遞歸方法，用加權(quán)最小二乘法對用一系列無量綱頻率表示的柔度矩陣C(ω)中每一項進行最優(yōu)化模擬。因此C(ω)在時域內(nèi)可通過位移序列和力序列之間的特定關(guān)系式來表示：

式中：ui、Fi分別為位移和廣義力矩陣中的一個元素；aj、bj、cj、dj、ej、hj為待定的離散過濾系數(shù)；m、k、l和r、s、n分別代表力和位移序列中各項的標(biāo)號。當(dāng)前時步輸出位移是當(dāng)前時步和過去時步的力與過去時步的位移的線性組合。F(t-r)、u(t-r)分別代表第t-r步的廣義力和位移。利用傅里葉變換，定義z=eiωΔ，其中Δ為兩個連續(xù)時間步之間的時間間隔，則可得到：

其中：

系數(shù)aj、bj、cj、dj、ej和hj可利用加權(quán)最小二乘法最小化以下殘差得到：

其中:W(ω)為權(quán)重，在本文中取統(tǒng)一值。系數(shù)求出之后即可通過式(5)解出時域內(nèi)當(dāng)前時步剛性基礎(chǔ)的位移。

1.2耦合子結(jié)構(gòu)的隱式和顯式算法

先介紹一個隱式算法，用于計算土-結(jié)構(gòu)相互作用體系的動力響應(yīng)。如圖1所示，在計算每一時步n時執(zhí)行以下操作：

(1) 假設(shè)當(dāng)前步的基礎(chǔ)位移為un，即為迭代過程中的變量。一旦un已知，則可求出當(dāng)前時步的其他變量(如：力、位移)。上一收斂時步(即第n-1時步)的基礎(chǔ)位移un-1可作為試算的初值un。

(2) 對于上部結(jié)構(gòu)，通過指定基礎(chǔ)位移un通過多點激勵方法計算其動力響應(yīng)，由有限元軟件OpenSees計算可得出基礎(chǔ)反力。

(3) 利用上部結(jié)構(gòu)作用于基礎(chǔ)上的力Fb,n、剛性基礎(chǔ)的質(zhì)量和位移，通過式(3)可以求出剛性基礎(chǔ)作用在土體上的廣義力Fs,n。其中當(dāng)前步的速度和加速度可以用當(dāng)前步和前一時步的位移、前一步的速度和加速度計算得到，比如用Newmark-Beta法插值得到。

(4) 通過式(2)，利用離散遞歸方法可求出額外位移us,n，即除了自由場輸入位移ug,n之外土的額外位移。從圖1中可看出，從步驟(2)～(4)可歸納為如下方程式：

其中：un即為剛性基礎(chǔ)的位移;us,n為土的額外位移。

圖1　土-結(jié)構(gòu)相互作用體系的隱式和顯式算法流程 Fig.1　Flowchart of the implicit algorithm and explicit 　　　 algorithm for SSI system

(5) 利用剛性基礎(chǔ)和土之間的位移協(xié)調(diào)條件，計算當(dāng)前時步基礎(chǔ)和土位移的殘差。由于已從步驟(2)～(4)中得出土的實際位移un=us,n+ug,n,=ψ(un)+ug,n，則該殘差可表示為：

如果殘差值小于用戶規(guī)定的誤差界限則該時步收斂，終止循環(huán)，進而轉(zhuǎn)入計算下一時步，即第n+1步，重復(fù)步驟(1)～(5)。如果達不到收斂誤差的要求，則通過牛頓更新法則給出新的試算值un，重復(fù)步驟(2)～(5)以繼續(xù)迭代直至收斂。

顯式算法與隱式算法有兩個不同：(1)步驟(1)中當(dāng)前步的基礎(chǔ)位移取為un=un-1,s+un,g(即用上一時步中土的額外位移代替當(dāng)前步的值)；(2)省去了步驟(5)，即不用檢查收斂性。當(dāng)精度要求相同時，顯式算法比隱式算法要省去很多個時間步的計算，因此對于大規(guī)模土-結(jié)構(gòu)相互作用體系來說，顯式算法是一個很好的選擇。

1.3CS方法的基本原理

本文中有限元軟件OpenSees和土-結(jié)構(gòu)相互作用體系之間簡單、高效的數(shù)據(jù)交流非常重要。如圖1所示在整個土-結(jié)構(gòu)相互作用體系的流程中，每個迭代步都需要多次讓OpenSees計算上部結(jié)構(gòu)的動力響應(yīng)，而土-結(jié)構(gòu)相互作用體系的非線性響應(yīng)計算是依賴于時間歷史的，這需要讓OpenSees常駐內(nèi)存。

為了實現(xiàn)OpenSees和SSI程序之間的協(xié)同工作，采用CS算法實現(xiàn)他們之間的集成、數(shù)據(jù)交流和實時通訊[18]。在土-結(jié)構(gòu)相互作用體系中插入一個客戶端，該客戶端是一個C++對象，此時Open-Sees將有限元模型儲存在其內(nèi)存中并作為一個服務(wù)器接收和執(zhí)行客戶端的指令，通過該客戶端即可實現(xiàn)OpenSees和土-結(jié)構(gòu)相互作用體系之間的實時數(shù)據(jù)交流。文獻[18]中有CS方法建立鏈接的詳細(xì)介紹。

2數(shù)值算例

2.1單自由度土-結(jié)構(gòu)相互作用體系

以一個均勻的彈性半空間的方形剛性表面地基和1層框架結(jié)構(gòu)組成的SSI系統(tǒng)為例。上部結(jié)構(gòu)高度為36.3m，基礎(chǔ)長度為25.2m，等效柱截面寬度為6.24m，樓面頂層的集中質(zhì)量為7.550N，梁、柱的彈性模量分別為2×1014N/m2和2.7×1010N/m2，阻尼比為1%。線性組合的梁、柱截面大小為b×h=6.24m×6.24m，梁上均布荷載為122.5kN/m，在結(jié)構(gòu)動力計算中低頻振型對動力響應(yīng)的影響要比高頻振型大得多，往往只需要求結(jié)構(gòu)體系的第一頻率和第二頻率，求得的一階頻率、二階頻率分別為2.33Hz和13.95Hz。土的材料參數(shù)：密度為1 874kg/m3;剪切波速為400 m/s;遲滯阻尼比為ξβ=0.001，ξα=0.000 5；泊松比為1/3。

將一組振幅為0.1m，頻率從1～4Hz的正弦波作為自由場地震動輸入，得到歸一化位移頻率響應(yīng)曲線(即輸出值除以輸入正弦波的幅值)，在本算例中最佳時間步長為0.003 9s(圖2)，數(shù)值解和解析解獲得的共振頻率完全一致(即1.76Hz)，初步驗證了該算法的正確性。

圖2　時間步長為0.003 9 s時解析和耦合法所求 　　　出歸一化的位移頻率響應(yīng)函數(shù)的比較 Fig.2　Comparison of normalized displacement frequency 　　　 responses using analytical and coupling methods 　　　 (time step is 0.003 9 second)

圖3　正弦波的幅值為0.001 m時隱式算法和顯式算 　　　法所求出非線性土-結(jié)構(gòu)相互作用體系的歸一 　　　化位移頻率響應(yīng)函數(shù)的比較 Fig.3　Comparison of normalized displacement frequency 　　　 response functions of nonlinear SSI system using implicit 　　　 and explicit algorithms (sine wave’s amplitude is 0.001 m)

將算例中的梁、柱單元改為非線性。結(jié)構(gòu)抗彎屈服彎矩設(shè)定為109 N·m，時間步長取0.003 9 s，選取振幅為0.001 m的正弦波。如圖3所示為隱式算法和顯示算法結(jié)果的比較，二者具有相同的共振頻率(即1.76 Hz)，雖然二者的峰值位移響應(yīng)有差別，但地震動力響應(yīng)時程的誤差很小，該誤差在實際工程中屬于可接受的范圍內(nèi)。

2.2美國加州理工大學(xué)Millikan圖書館的非線性地震動力響應(yīng)

基于本文所提出的耦合算法計算美國加州理工大學(xué)Millikan圖書館的非線性地震動力響應(yīng)(圖4)，以驗證該方法對土-結(jié)構(gòu)相互作用問題的適用性。

圖4　美國加州理工大學(xué)Millikan 圖書館 Fig.4　The Millikan library in California Institute of Technology

該圖書館高43.9m，寬21m，利用OpenSees模擬上部結(jié)構(gòu)，采用基于位移的非線性梁/柱單元，該單元具有雙線性的非線性彎矩-曲率關(guān)系。柱單元的參數(shù)如下：屈服彎矩為1.15×108 N·m，EI=1.40×1011N·m2，屈服后剛度與初始剛度的比值為0.01。梁用彈性梁柱單元模擬,EA=1.29×1011N。基礎(chǔ)假設(shè)為剛性，土的材料參數(shù)與單自由度算例一致。該結(jié)構(gòu)的一階固有頻率為2.33Hz。

由耦合法求得歸一化的位移頻率響應(yīng)函數(shù)如圖5所示。通過將簡諧地震激勵的幅度設(shè)置為非常小的值，該土-結(jié)構(gòu)相互作用體系仍處于線彈性狀態(tài)。用該方法求出結(jié)構(gòu)的共振頻率為1.82Hz，非常接近原位強制振動測試結(jié)果1.90Hz[20]。結(jié)構(gòu)頂部的絕對位移值為38.8，轉(zhuǎn)動和平動位移幅值分別為10.71和1.95。

圖5　耦合法求得Millikan 圖書館考慮土-結(jié)構(gòu) 　　　 相互作用的歸一化的位移頻率響應(yīng)函數(shù) Fig.5　 The normalized displacement frequency response 　　　 functions of the Millikan library considering SSI 　　　 using the coupling method

地震激勵改為SanFernando地震，結(jié)構(gòu)明顯屈服并產(chǎn)生強非線性行為。為了研究土-結(jié)構(gòu)相互作用對結(jié)構(gòu)的影響，將考慮和不考慮土-結(jié)構(gòu)相互作用下結(jié)構(gòu)的響應(yīng)時程進行對比。各層的位移和加速度包絡(luò)線如圖6所示，從圖中得出當(dāng)考慮土-結(jié)構(gòu)相互作用時各層的最大位移均比僅考慮結(jié)構(gòu)系統(tǒng)的最大位移值大，并且隨著高度的增加該差別會越來越大。土-結(jié)構(gòu)相互作用體系頂層的最大位移為0.262 1m，只考慮結(jié)構(gòu)體系頂層的最大位移為0.207 5m。相反,土-結(jié)構(gòu)相互作用體系最大加速度要小于僅考慮結(jié)構(gòu)體系的加速度[在圖7(c)中也可以觀測到該現(xiàn)象]。

圖6　考慮和不考慮土-結(jié)構(gòu)相互作用下各層的 　　　位移和加速度包絡(luò)線 Fig.6　Enveloping lines of displacement and acceleration 　　　 at each storey with and without SSI effect

結(jié)構(gòu)的響應(yīng)時程如圖7所示，考慮和不考慮土-結(jié)構(gòu)相互作用下結(jié)構(gòu)不同位置的彎矩曲率響應(yīng)如圖8所示。從結(jié)果中看出，土-結(jié)構(gòu)相互作用對于結(jié)構(gòu)系統(tǒng)有不利的影響。隱式和顯式算法的結(jié)果非常相近，該誤差在工程中屬于可接受范圍，說明該顯式算法可以用于現(xiàn)實的大尺度工程中。

圖8　結(jié)構(gòu)不同位置的彎矩曲率響應(yīng) Fig.8　Moment-curvature response at different position 　　　 of the structure

3結(jié)論

提出一個實用、高效的耦合解析解和數(shù)值解的算法，用于土-結(jié)構(gòu)相互作用體系的地震分析。該方法即利用有限元軟件在大規(guī)模線性/非線性結(jié)構(gòu)系統(tǒng)模擬中的強大分析功能，并利用時域離散遞歸方法對半空間土體的高效分析，因此可能應(yīng)用到解決大型民用基礎(chǔ)設(shè)施的實際土結(jié)相互作用問題中。

將該方法應(yīng)用到一個單自由度土-結(jié)構(gòu)相互作用體系和美國加州理工大學(xué)Millikan圖書館實際工程算例中，驗證其正確性。此方法可應(yīng)用到大規(guī)模土-結(jié)構(gòu)相互作用體系的非線性地震分析中，并可用于分析是否考慮土的相互作用的區(qū)別。本文提出的方法實用、高效，為分析大規(guī)模非線性土-結(jié)構(gòu)相互作用體系的動力時程響應(yīng)提供了重要工具，該研究的部分結(jié)果對工程實踐有一定指導(dǎo)意義。

參考文獻(References)

[1]SafakE.DetectionandIdentificationofSoil-structureInteractioninBuildingsfromVibrationRecordings[J].JournalofStructuralEngineering, 1995, 121(5): 899-906.

[2]StewartJP，F(xiàn)envesGL，SeedRB.SeismicSoil-structureInteractioninBuildings.I:AnalyticalMethods[J].JournalofGeotechnicalandGeoenvironmentalEngineering,1999,125(1):26-37.

[3]JenningsPC，BielakJ.DynamicsofBuilding-soilInteraction[J].BulletinoftheSeismologicalSocietyofAmerica,1973,63：9-48.

[4]LucoJE.LinearSoil-StructureInteraction[R].ReportUCRL-15 721.Livermore，California:LawrenceLivermoreNationalLaboratory,1980.

[5]LucoJE,TrifunacMD,WongHL.IsolationofSoil-structureInteractionEffectsbyFfull-scaleForcedVibrationTests[J].EarthquakeEngineeringandStructuralDynamic,1988,16:1-21.

[6]MerrittGH,HousnerGW.EffectsofFoundationComplianceonEarthquakeStressesinMultistoryBuildings[J].BullSeismSocAmer,1954,44(4):551-569.

[7]GuQ.PerformanceandRiskAssessmentofSoil-StructureInteractionSystemsBasedonFiniteElementReliabilityMethods[J].MathematicalProblemsinEngineering,2014(2):1-16.

[8]FatahiB,TabatabaiefarH,SamaliB.PerformanceBasedAssessmentofDynamicSoil-structureInteractionEffectsonSeismicResponseofBuildingFrames[J].Geo-Risk,2014:344-351.

[9]TrifunacMD.DiscussionofSeismicSoil-structureInteractioninBuildings.I:AnalyticalMethods;II:Empiricalfindings[J].JournalofGeotechnicalandGeoenvironmentalEngineering,2000,126(7):668-670.

[10]VeletsosAS,MeekJW.DynamicBehaviorofBuilding-foundationSystems[J].JournalofEarthquakeEngineeringandStructuralDynamics,1974,3(2):121-138.

[11]WongHL,TrifunacMD,LucoJE.AComparisonofSoilStructureInteractionCalculationswithResultsofaFull-scaleForcedVibrationTest[J].SoilDynamicsandEarthquakeEngineering,1988,7(1):22-31.

[12]MylonakisG,GazetasG.SeismicSoil-structureInteraction:beneficialorDetrimental?[J].JournalofEarthquakeEngineering,2000,4(3):277-301.

[13]JeremicB,KunnathS,XiongF.InfluenceofSoil-foundation-structureInteractiononSeismicResponseofTheI-880Viaduct[J].EngStructures, 2004,26:391-402.

[14]TangY,ZhangJ.ProbabilisticSeismicDemandAnalysisofaSlenderRCShearWallConsideringSoil-structureInteractionEffects[J].EngineeringStructures, 2011, 33(1):218-229.

[15]StewartJP,SeedRB,FenvesGL.SeismicSoil-structureInteractioninBuildings.II:EmpiricalFindings[J].JournalofGeotechnical&GeoenvironmentalEngineering,ASCE,1999,128:38-48.

[16]SafakE.Time-domainRepresentationofFrequency-dependentFoundationImpedanceFunctions[J].SoilDynamicsandEarthquakeEngineering,2006,26(1): 65-70.

[17]MazzoniS，McKennaF，F(xiàn)envesGL.OpenSeesCommandLanguageManual[R].PacificEarthquakeEngineeringCenter,UniversityofCalifornia,Berkeley,Calif,USA,2005,http://opensees.berkeley.edu/.

[18]GuQ,OzcelikO.IntegratingOpenSeesWithOtherSoftware——WithApplicationtoCouplingProblemsinCivilEngineering[J].StructuralEngineeringMechanics，2011,40 (1):85-103.

[19]LennartL.SystemIdentification:TheoryfortheUser.1998.

[20]JohnFC,ModernDigitalSeismology-Instrumentation,andSmallAmplitudeStudiesintheEngineeringWorld[D].Pasadena,California:CaliforniaInstituteofTechnology,2004.