分析土和結構相互作用的一種實用耦合方法①

分析土和結構相互作用的一種實用耦合方法①

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

古泉， 彭伊, 曾志弘

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

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

關鍵詞：土-結構相互作用； 子結構法； CS方法； 數值解與解析解耦合方法； OpenSees; 非線性地震動力響應分析

收稿日期：①2014-08-20

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

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

中圖分類號:TU352.12; U441.3文獻標志碼：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引言

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

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

本文在子結構法的基礎上，提出數值解與解析解耦合的分析方法。該方法結合有限元分析和基于時域離散遞歸方法的頻域內解析解，既能滿足復雜的非線性結構體系的精確分析，又能實現土體的簡單快速模擬，對于大規模的線性或者非線性土-結構相互作用體系的數值分析和模擬具有較大的應用潛力。

1基于子結構法的耦合方法

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

1.1土-結構相互作用體系頻域解析解的時域表示

剛性基礎的總位移為：

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

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

M0為剛性基礎的質量矩陣;Fb為上部結構作用在剛性基礎的廣義力，基于本文所提出的耦合方法，時域內的廣義力矩陣Fb可由有限元軟件計算得到，即通過OpenSees分析得到。

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

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

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

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

其中：

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

其中:W(ω)為權重，在本文中取統一值。系數求出之后即可通過式(5)解出時域內當前時步剛性基礎的位移。

1.2耦合子結構的隱式和顯式算法

先介紹一個隱式算法，用于計算土-結構相互作用體系的動力響應。如圖1所示，在計算每一時步n時執行以下操作：

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

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

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

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

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

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

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

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

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

1.3CS方法的基本原理

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

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

2數值算例

2.1單自由度土-結構相互作用體系

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

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

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

圖3　正弦波的幅值為0.001 m時隱式算法和顯式算 　　　法所求出非線性土-結構相互作用體系的歸一 　　　化位移頻率響應函數的比較 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)

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

2.2美國加州理工大學Millikan圖書館的非線性地震動力響應

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

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

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

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

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

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

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

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

圖8　結構不同位置的彎矩曲率響應 Fig.8　Moment-curvature response at different position 　　　 of the structure

3結論

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

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

參考文獻(References)

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

[2]StewartJP，FenvesGL，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，FenvesGL.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.