樁-土-結構相互作用分析中土動力模型的修正

鄧浩昀， 金新陽， 顧 明， 黃吉鋒

(1. 同濟大學 土木工程防災國家重點實驗室， 上海 200092； 2. 中國建筑科學研究院 建筑結構研究所， 北京 100020)

樁-土-結構相互作用已成為高層建筑、橋梁工程、核電站工程等結構抗震設計的研究重點，主要分析方法包括整體分析法和子結構法.整體分析法將上部結構、樁基礎和土作為整體進行計算，充分考慮土的非線性、樁和土的相互作用以及土和結構的慣性相互作用，比較符合實際情況.在整體分析法中，場地無限域模擬和土動力模型選擇是影響分析精度的2個主要因素.對于無限域模擬，通常在場地四周施加人工邊界來模擬遠場地基介質的輻射阻尼效應.土動力模型主要有以下2類：等效線性模型(ELM)和基于各向異性運動硬化的塑性模型[1].對于等效線性模型，假定地震過程中土的材料性質是線性的，在每一次迭代計算中土的剪切模量和阻尼比與等效剪切應變相關.與塑性模型相比，等效線性模型能夠減少計算難度，同時由于概念明確，因而得到廣泛應用.

Ishihara[2]研究發現，強震作用下基礎近場區域的土會產生較大的剪切應變，非線性強.若地基土選用等效線性模型，則結構動力分析得到的結果不夠精確.Yoshida等[3]也指出，強震作用下土的動力模型采用等效線性模型會使結構的位移峰值增大.然而，當前在進行樁-土-結構相互作用分析時，這些誤差往往被忽略.因此，在分析強震作用下樁-土-結構體系的動力響應時，需要對近場區域土的動力模型進行修正.常用的方法是將場地劃分為規則的近場區域和遠場區域，近場區域采用非線性模型，遠場區域采用線性模型.Casciati等[4]建立了樁-土-結構模型并進行動力分析，其中近場區域使用等效線性模型，遠場區域使用線性模型，驗證了該模型具有良好的精度.然而，地基土的等效剪切模量和阻尼比是通過SHAKE91自由場分析求得后導入三維有限元模型的，未考慮上部結構對地基土的影響.

本文基于已有研究提出一種土動力模型修正方法.首先，建立樁-土-結構有限元模型，其中土動力模型為等效線性模型，并通過分析得到近場區域的尺寸.其次，將近場區域土的剪切模量乘以修正系數，得到修正等效線性模型(MELM)，使修正等效線性模型的動力分析結果與摩爾-庫倫(MC)模型趨于一致.然后，對100個算例重復上述分析，得到每個模型的修正系數.最終，對修正系數進行回歸分析，獲得修正系數與結構周期和場地周期的經驗公式.

1 有限元模型

1.1 接觸

結構柱(梁單元)-筏板(實體單元)之間選用動力耦合的方法，耦合梁單元節點和對應柱底面積區域內節點的所有自由度.土(實體單元)-筏板之間選用接觸離散算法，法向為硬接觸，切向為庫倫摩擦模型.樁(梁單元)-土之間采用嵌入模型，土的動力反應將作為位移約束條件施加給樁的節點，實現樁-土動力耦合.

1.2 結構阻尼

結構阻尼采用Rayleigh阻尼，阻尼矩陣C寫為質量矩陣M和剛度矩陣K的線性比例之和，如下所示：

C=αM+βK

(1)

由瑞利阻尼理論可知，α和β可以通過下式確定：

(2)

(3)

式中：ξi、ξj為2個振型的阻尼比，取0.02；ωi取結構第1階頻率，ωj取結構第2～20階頻率的平均值[5].

1.3 重力荷載施加

已有研究[6]表明：在進行動力分析時，忽略重力將無法考慮初始應力對接觸狀態的影響，不能真實反映土與結構的受力和變形情況.本文通過2步施加重力荷載：① 固定地基土底部邊界并對整個模型施加重力加速度場(從0逐漸增加至9.8 m·s-2)，得到整個體系的應力狀態文件；② 將獲得的應力狀態文件作為初始條件施加到模型上，并同時施加重力加速度9.8 m·s-2，待整個模型的位移趨于穩定，施加地震波進行動力分析.

1.4 等效線性模型

在受到地震荷載作用時，土的應力-應變關系表現為復雜的滯回曲線.等效線性模型采用等效剪切模量G和等效阻尼比λ來反映土的滯回曲線，G和λ為應變γ的函數.Hardin等[7]根據試驗資料提供了如下經驗公式：

(4)

式中：Gmax和γd分別為最大剪切模量和參考應變，可根據試驗數據進行確定.λ在工程中通常采用如下經驗公式[8]：

(5)

最大阻尼比λmax可通過試驗或經驗公式得到.本文中G/Gmax-γ和λ-γ曲線采用孫靜[9]對國內不同地區的土樣進行共振柱試驗的結果.

Abaqus軟件中，等效線性模型可以通過Kelvin模型來反映土體在周期載荷下的滯回性.該模型由線彈性彈簧和阻尼裝置并聯組成，應力-應變關系為

(6)

式中：E和η分別為模量系數和黏滯系數；σ為應力；ε為應變.

將一維模型推廣至三維模型的情況，如下所示：

(7)

(8)

式中：ν和μ為拉梅系數；σii和σij分別為正應力和剪應力；εii和γij分別為正應變和剪應變張量；εV=εii+εjj+εkk，ην=2νξ/ω，ημ=2μξ/ω，其中ξ為阻尼比，ω為結構基頻.

等效線性模型子程序在Abaqus軟件中的開發及使用流程如圖1所示.文獻[1]中指出，等效線性模型中，等效剪切應變可取動力作用過程中最大應變與折減系數0.65的乘積.

圖1 等效線性模型子程序流程圖Fig.1 Subroutine flow chart of equivalent linear model

1.5 人工邊界及場地寬度

黏彈性人工邊界是當前使用較多的一種局部人工邊界.劉晶波等[10]基于三維球面波動方程推導出黏彈性邊界的法向和切向人工邊界條件，并將其等效為人工截斷邊界上連續分布的并聯彈簧-阻尼系統.在Abaqus軟件中可以通過在邊界節點施加彈簧單元和阻尼器單元來實現，并設置相應的剛度系數Ki和阻尼系數Ci，地震動的輸入方法可以通過施加等效荷載來實現[11].

場地范圍的選擇是樁-土-結構相互作用分析中的重要問題.若截取場地范圍過大，則計算時間太長；反之，散射波在人工邊界上產生的反射會對計算結果產生較大影響.已有研究對場地的截取范圍結論不一[12-13]，需要進一步分析場地范圍對結構動力響應的影響.

本文結構底部尺寸為20 m×20 m，高度為20 m，對4種場地范圍的樁-土-結構模型進行動力分析，得到結構加速度峰值及相對位移峰值(見表1).結果表明：當場地寬度(100 m×100 m)為結構寬度的5倍時，結構的相對位移及加速度與400 m×400 m場地結果的誤差分別為1.42%和1.46%，基本可以消除側向邊界對結構動力響應的影響.

表1 不同場地范圍結構的動力響應峰值Tab.1 Maximum value of structural dynamic response in different site dimensions

本文研究的是一維地震波輸入(x-z截面)，故y-z截面方向的寬度可以適當減小，有利于提高計算效率.在表1選定的場地基礎上，增加了2組場地進行了對比分析(見表2).結果表明：60 m×100 m、80 m×100 m的結果與100 m×100 m的結果十分接近，故最終選用的場地范圍為60 m×100 m.

2 模型參數

2.1 結構參數

本文研究了5層、10層、15層和20層的四跨鋼框架.層高3 m，跨度5 m，構件密度7 850 kg·m-3.樓板厚度在0.1～0.2 m之間.5層和10層結構使用筏板基礎， 15層和20層結構使用樁筏基礎.筏板尺寸為22 m×22 m×1 m，群樁按5×5的方式排列，間距為5 m，樁直徑為0.5 m.

表2 yz截面不同場地范圍結構的動力響應峰值Tab.2 Maximum value of structural dynamic response in different y-z site dimensions

2.2 場地類型及參數

根據《建筑抗震設計規范》，建筑場地類別分為Ⅰ、Ⅱ、Ⅲ、Ⅳ 4類，樁-土-結構相互作用對軟土地基上的結構影響較大，故本文只分析Ⅱ、Ⅲ、Ⅳ 3類場地，其中Ⅱ類和Ⅲ類場地各2種.各場地的分層情況和力學參數如表3所示，泊松比均為0.35.

表3 場地土力學參數Tab.3 Mechanic parameters of soil at different sites

已有文獻[14]中給出了場地周期的計算公式，如下所示：

(9)

式中：N為場地土分層數；hi和vi分別為第i層厚度和剪切波速；Hi為第i層中點的深度.

2.3 地震波的篩選及反演

本文通過初選(基于臺站與地震信息[15])和復選(基于設計反應譜的雙頻段控制[16])2個步驟篩選符合條件的地震波.

根據各場地的地震反應譜及初選參數從美國太平洋工程地震研究中心數據庫(PEER)選出30條地震波，通過文獻[15]將我國的地震波參數轉換成PEER地震波初選參數，如表4所示.

表4 地震波初選參數Tab.4 Primary selection parameters of the earthquake

復選是將初選的30條地震波調幅至6.2 m·s-2(9度罕遇)，分別記錄30條地震波的反應譜在[0.1，Tg]和[T1-0.2，T1+0.5]2個區間的值；隨后，從30條地震波中選出與目標場地反應譜在2個區間內平均誤差最小的5條，用于后續的時程分析.其中，Tg為場地特征周期，T1為結構第1階周期.

本文采用等效線性水平成層土動力響應分析程序SHAKE91對地表地震波進行反演，該程序基于一維波動理論[12]可以求解出地表層(第N層)和基巖層(第1層)的位移傳遞函數QN,1，隨后可以確定地基土的底邊界(基巖面)的地震波時程(位移、速度、加速度)，如下所示：

(10)

3 近場土修正方法

3.1 地基土修正區域

已有文獻[2]表明：強震作用下，地基土的剪切應變超過0.9%時，不宜使用等效線性模型.選取100個算例(4種結構、5類場地、5條地震波)計算各樁-土-結構模型超過0.9%剪切應變的地基土區域，并對其中一個算例進行分析(見圖2).水平地震作用下，在基礎附近地基土的剪切應變明顯大于其他區域，并且數值大于0.9%，需要對該區域的等效線性模型進行修正.其他算例的結果與此算例的分布基本一致.

圖2 近場區域土的剪切應變(20層結構，Ⅲ類(1)場地)Fig.2 Shear strain in the near-field soil (20-storey building, soil type Ⅲ(1))

基于上述工作進一步確定地基土修正區域的尺寸，可將場地土分成若干區域(SS1～SS5).x-z截面和y-z截面的分區如圖3所示.

a xz截面b yz截面

圖3場地土分區

Fig.3Soilpartition

通過動力分析，求出圖3中各區域最大剪切應變的平均值，選取場地土半側區域的結果進行對比.由圖4可以發現：SS1-1、SS2-1以及SS4-1區域最大剪切應變的平均值大于0.9%，需要對該區域等效剪切模量進行修正.將SS1-1、SS2-1以及SS4-1組合起來，可以得到以下修正區域的范圍：寬度為結構寬度的1.5倍，深度為結構寬度的0.5倍.通過對其他算例的分析，發現修正區域的范圍和此算例基本一致.

a xz截面b yz截面

圖4各區域最大剪切應變

Fig.4Maximumshearstrainofeachregion

3.2 修正區域剪切模量

在樁-土-結構模型修正區域中分別采用以下3種土動力模型：① 等效線性模型；② 修正等效線性模型，修正系數k=Gnew/Gold，其中Gnew和Gold分別為修正后和修正前近場區域土的剪切模量；③ 摩爾-庫倫模型.對每個模型進行動力分析，調節修正等效線性模型的修正系數，直至修正等效線性模型與摩爾-庫倫模型中上部結構的動力響應趨于一致.

采用上述方法對100個算例進行分析，并對其中4個算例進行對比，分析不同地基土模型中，上部結構水平位移的比值(見圖5).結果表明：

(1) 由圖5a和圖5c對比得出，隨著場地周期增大，場地變軟，等效線性模型與摩爾-庫倫模型中結構水平位移的誤差越大.

(2) 由圖5a和圖5d對比得出，結構高度越大， 等效線性模型與摩爾-庫倫模型中結構水平位移的誤差越大.

(3) 由圖5a和圖5b對比得出，對相同的樁-土-結構模型，不同地震波作用下位移比ΔELM/ΔMC的值非常接近，表明同一模型下不同地震波對動力分析結果影響很小.

(4) 在強震作用下，采用修正等效線性模型得到的結構水平位移，比等效線性模型更接近摩爾-庫倫模型的結果(ΔMELM/ΔMC的值比ΔELM/ΔMC更接近1)，故本文針對等效線性模型提出的修正方法是合適的.

a 20層-Ⅲ類(1)場地-地震波1

b 20層-Ⅲ類(1)場地-地震波2

c 20層-Ⅱ類(1)場地

d 5層-Ⅲ類(1)場地

圖5不同模型結構水平位移對比

Fig.5Comparisonofstructurallateraldisplacementamongdifferentmodels

3.3 近場區域修正系數及回歸公式

采用上述方法對100個算例進行動力分析，并分別得到相應的修正系數，以20層結構的分析結果(見表5)為例，可以發現相同樁-土-結構模型在不同地震波作用下的修正系數非常接近，這是因為在第2.3節中，每一組地震波都是基于相同的目標地震反應譜選出來的.

表5 20層結構不同場地修正系數Tab.5 Modified coefficients of 20-storey building in different sites

圖6的點為每一個樁-土-結構模型修正系數的平均值，采用回歸分析對每一條曲線進行擬合，得到便于工程使用的經驗公式.

各曲線的經驗公式可統一表示為

k=(a+bT)T2+(c+dT)T1+e

(11)

式中：a、b、c、d、e為參數，根據不同的場地條件取相應值(見表6).修正系數的計算值和真實值的平均誤差在2%以內.

圖6 剪切模量修正系數Fig.6 Modified coefficients of shear modulus表6 經驗公式參數Tab.6 Parameters of the empirical equation

場地條件abcdeⅡ類-0.045 30.090 10.216 1-0.151 41.075 1Ⅲ類0.042 70.101 00.118 3-0.243 01.189 2Ⅳ類00.321 1-0.402 7-0.101 11.525 3

4 試驗算例

為驗證本文提出的修正方法和修正系數經驗公式，采用已有試驗算例[17]進行驗證.選用文獻中的模型P4-A2，計算模型如圖7所示，基巖輸入地震波為El-Centro波，調幅峰值為0.2g.

分別采用摩爾-庫倫模型、修正等效線性模型和等效線性模型土進行分析，得到結構頂部的加速度峰值如表7所示.3種有限元模型計算結果與試驗結果非常接近，誤差分別為2.16%、3.24%、3.96%.這是因為入射地震波峰值僅為0.2g，未能使地基土產生較強的非線性，等效線性模型和非線性模型結果趨于一致，但是分析結果驗證了本文有限元模型的正確性.

圖7 算例的三維有限元模型Fig.7 Three-dimensional finite element model of the example

表7 結構頂部加速度峰值Tab.7 Peak values of acceleration at the top of structure

隨后將基巖地震波調幅至0.62g，從基巖重新輸入，分別采用摩爾-庫倫模型、修正等效線性模型和等效線性模型進行分析，得到結構的水平位移和層間剪力，如圖8所示.結果表明：對于試驗算例，采用修正等效線性模型得到的水平位移和層間剪力與摩爾-庫倫模型的平均誤差分別為2.39%和3.57%，小于與等效線性模型的平均誤差(6.45%和9.83%)，說明了本文提出修正方法的正確性.

a 水平位移b 層間剪力

圖8不同模型結構水平位移和層間剪力對比

Fig.8Comparisonofstructurallateraldisplacementandstoreyshearamongdifferentmodels

5 結論

(1) 強震作用下，樁-土-結構模型的基礎附近會出現強非線性區域，其寬度為超出結構寬度的50%，深度為結構寬度的50%.

(2) 本文分析了100個算例，得到樁-土-結構模型近場區域的修正系數，范圍在1.03至1.81之間，并且隨著場地變軟，結構變高，修正系數增大.

(3) 本文對所有算例的修正系數進行回歸分析，得到修正系數的經驗公式，公式計算值和真實值的誤差不超過2%.

(4) 在進行樁-土-結構相互作用分析時，相比等效線性模型而言，采用修正等效線性模型得到結構的動力響應與摩爾-庫侖模型的結果更為接近，故本文提出的修正方法和經驗公式具有較好的適用性和正確性.