可液化地基地下結構地震反應特征簡化有效應力分析

趙 凱, 王秋哲, 王彥臻, 莊海洋, 陳國興

(南京工業大學 巖土工程研究所， 南京 210009)

大量震害現象表明地下結構在強地震作用下可能會出現嚴重的震害以及次生災害[1-8]。由于地下工程結構受周圍土體約束，土-結構動力相互作用關系十分復雜。采用動力有效應力法分析可液化地基地下結構的地震反應特性，可以更科學地解釋強地震動作用下土-地下結構體系的運動模式及破壞機理。

目前國內外學者對可液化地基地下結構動力有效應力分析方法進行了系統研究，歸納起來主要有兩大類：① 不排水解耦有效應力法：假定在動力作用期間土體中的孔隙水封閉在土骨架中，不向外排出。在總應力分析法的基礎上，增加殘余孔隙水壓力或殘余變形的計算。陳國興等[9-13]基于ABAQUS軟件平臺，實現了修正Davidenkov本構模型與修正Byrne孔壓增量模型的有效應力算法，通過地下結構非線性精細化建模和高效的有限元數值模擬平臺，系統地研究了輸入地震動特性、場地條對復雜地鐵地下車站結構體系非線性地震反應特征、結構損傷與空間效應特性的影響。② 基于Biot動力固結理論的雙相介質耦合場統一算法：即采用一定的數值積分格式直接對Biot動力方程進行逐步積分計算的方法。Prevost等[14-15]提出了一種將Biot方程應用于飽和多孔介質動力分析的有限元計算方法，并應用于可液化地基中地下結構的地震反應分析。Khoshnoudian等[16]采用考慮運動硬化的動彈塑性模型描述土的特性，基于Biot部分動力固結模型，數值模擬了可液化地基中內襯隧道的地震反應特性。黃雨等[17]采用了Armstrong-Frederick 型非線性運動硬化準則和非關聯流動準則來描述動荷載作用下砂土的循環活動性以及液化強度，基于Biot動力固結理論對液化場地樁基礎的地震反應進行了三維有效應力分析。王剛等[18]基于能夠模擬飽和砂土液化后大應變響應的彈塑性循環本構模型，采用DIANA-SWANDYNE II對阪神地震中破壞的大開車站進行了有效應力分析，說明了考慮液化大變形的土與地下結構動力相互作用分析方法及其有效性。鄒炎等[19]采用Yang等[20]提出的砂土液化大變形本構模型，基于OpenSees對飽和砂土場地的地震反應進行了非線性動力有限元分析，闡述了地基的液化機理。上述研究中均采用復雜的黏彈塑性本構模型(如邊界面模型和嵌套面模型)來描述土的動應力-應變關系。這將會導致Biot動力固結方程求解難度大，數值計算收斂困難等難題，且對于使用者的要求較高[21-22]。

大量試驗研究表明，Masing型非線性動本構模型能較為理想地描述各類土體的非線性動力特性，模型參數少且物理意義明確，從而使得該類模型成為應用最為廣泛的黏彈性動本構模型。然而，該類動本構模型僅提供了一維非線性剪應力-剪應變關系，且不能夠描述循環剪切和壓縮作用引起的塑性體應變[23]。為此，本文給出了適宜于將一維應力-應變關系向三維空間擴展的等效剪應變算法及加卸載判據，并基于Masing法則構造了三維應力空間中的Davidenkov本構模型；同時補充了筆者提出的剪應力和正應力差耦合作用下剪切-體積應變耦合模型作為Biot動力固結方程中計算超靜孔隙水壓力的源項；隨后基于FLAC3D軟件平臺，建立了適宜于分析可液化地基中土-地下結構相互作用問題的有效應力算法。最終利用該算法初步研究了可液化地基中地下結構周圍場地的地震反應規律。

1 可液化場地動力有效應力分析方法

1.1 控制方程

結合Biot動力固結控制方程，考慮地震荷載作用下飽和土體積變化的相容條件，建立土骨架-孔隙流體兩相介質動力有效應力分析方法。假設均質各向同性多孔介質，孔隙水的運動符合Darcy定律，土骨架和孔隙水是可壓縮的，控制方程如下

(1)

(2)

式中：?為拉普拉斯算子；上標“·”為對時間求導；p為超孔壓；ns為孔隙率；k為滲透系數；u為土單元位移矢量；ε為土單元應變張量；G為動剪切模量；v為泊松比；γw為孔隙流體重度；Kf為流體壓縮系數。

(3)

(4)

雙向耦合剪切條件下剪切-體積應變耦合模型表達式如下

(5)

λ=RES-RESt

(6)

1.2 循環本構方程

采用Masing法則構造三維應力空間中的Davidenkov本構模型描述土骨架的非線性和滯回性[28-29]，如圖1所示。初始骨架曲線的表達式為

τ=G·γ=Gmax·γ·[1-H(γ)]

(7)

其中

(8)

圖1 修正Davidenkov應力-應變關系示意圖Fig.1 Schematic illustration of hypsteretic stress-straincurve of modified Denvidenkov model

由Masing法則構造的滯回曲線表達式為

(9)

式中：τ、γ分別為剪應力和剪應變；Gmax為最大動剪切模量；A、B和γc為土的試驗參數；γc為加卸載轉折點處的應變。

Masing型非線性動本構模型的核心是時變剪切模量的計算，其準確性與等效剪應變的選取密切相關，即等效剪應變算法的選取是本構關系由一維空間向三維空間擴展的關鍵問題。以應變偏量的第二不變量描述的等效剪應變能綜合反映三維空間中土體的應力-應變關系，而直接采用應變偏量的第二不變量作為等效剪應變時，其數值只能描述一維應力-應變關系曲線中應變大于零的部分，需將等效剪應變改由增量形式通過加卸載判斷后再次計算[30]

(10)

(11)

(12)

(13)

1.3 土骨架模型與孔壓增量模型的耦合方式

動力有效應力分析方法可以合理地考慮土骨架與孔隙水兩相介質相互作用以及場地液化過程中土體力學性質的變化。如果在動荷載施加過程中孔隙流體的排水通道受阻，對于同一土骨架模型，可以采用上一狀態的骨架曲線來預測當前狀態的骨架曲線，即初始剪切模量隨結構孔壓(亦即隨有效應力變化)的累積而減小。

(14)

2 模型驗證

Martin等基于干砂的常幅循環直剪試驗，建立了剪應變幅值相關的體應變增量模型。Byrne[31]將Martin等采用多項式給出的體應變增量表達式改為指數函數形式，該模型僅有兩個參數且物理意義也更為明確，廣泛用于地下結構地震反應分析[32-33]。

(15)

式中，參數C1和C2可由修正標準貫入錘擊數Nl或相對密度Dr(以百分數表示)確定

C1=7 600(Dr)-2.5或C1=8.7(N1)-1.25

(16)

C2=0.4/C1

(17)

圖2 不同分析模型預測值與實測值的對比Fig.2 Comparison between the measured and predicted excess pore pressure time-histories based on different models

表1 模型計算參數Tab.1 Mechanical parameters of the soil element

圖2為本文模型、修正Byrne模型以及室內試驗獲得的超孔壓增長時程曲線。由圖可知，采用修正Byrne模型能夠較好地預測一維剪應力路徑下超孔壓的發展，而在正應力差和剪應力耦合作用下(即橢圓形和圓形應力路徑)的預測結果遠小于試驗結果。筆者建立的模型預測結果與試驗數據均吻合較好，能夠反映應力路徑變化對超孔壓發展的影響規律。這是由于修正Byrne模型由直剪試驗的數據建立，難以反映正應力差和剪應力的耦合變化引起土單元的超孔壓發展規律。飽和砂土的孔壓增長規律與應力路徑密切相關，對于RCS相同而a/b不同的雙向耦合剪切試驗，對試樣所施加應力路徑不一致，因此受到的剪切效果也不相同。值得一提的是，應力孔壓(超孔壓波動值)隨著循環周數增加逐步放大，這是由于結構孔壓的增長引起土骨架模量衰減導致的。

圖3以Chen等給出的應變控制不排水循環三軸加載試驗驗證本文模型的適用性。可以看出，本文模型預測的結果對于不同相對密度，不同加載類型的試驗結果均較為理想。

圖3 應變控制不排水循環三軸試驗預測值與實測值對比

3 算 例

3.1 計算模型建立

選取某砂質場地矩形截面隧道為算例，土層總厚度為60 m，底面為不透水基巖，隧道為三跨結構，總寬L= 40 m，邊跨寬16 m，中跨寬8 m，高為10 m。頂板、底板及側墻厚度均取1 m，隧道埋深為1 m，結構采用C45混凝土。模型如圖4所示，采用四節點等應變有限差分單元，網格大小約為0.25 m × 0.25 m，土-結相互作用接觸面為法向硬接觸，接觸面切向采用摩擦接觸。土-結接觸面法向剛度kn= 350 GPa，剪切剛度ks= 350 GPa，摩擦系數μ=0.7。場地的剪切波速Vs= 203 m/s，滲透系數K= 1×10-5cm/s，飽和重度γ= 19.8 kN/m3，臨界阻尼比D取值為5%。模型計算參數見表1。

圖4 模型網格示意圖Fig.4 Grid layout of computation domain in detail

場地兩側邊界采用豎向約束、水平向加彈簧阻尼器的黏彈性邊界，場地底部基巖面僅施加豎向約束。假定地震動為基巖面豎直向上傳播的SH波，由模型底部邊界輸入[34]。采用地震工程中常用的Kobe波作為輸入地震動，將輸入地震波的持續時間截取為30 s。截取后的地震波加速度時程及傅氏譜如圖5所示，在輸入地震動時將峰值調幅為0.1g，在計算過程中，僅考慮地震波單向垂直入射。

圖5 輸入地震動原始加速度時程及傅里葉譜Fig.5 Original acceleration time-history and Fourier spectraof the input Kobe wave

3.2 結果分析

圖6為場地不同位置處的土單元應力路徑((σz-σx)/2 ～τzx)。由圖可見，遠場P1(距離結構側墻100 m)在地震作用過程中處于正應力不變條件下規則的往復水平剪切狀態，應力路徑為直線型。而結構周圍場地土單元體經歷的應力路徑則明顯不同，表現為正應力差與剪應力連續變化，P2(距離結構側墻30 m)的應力路徑為斜橢圓型，而P3(距離結構側墻3 m)為斜直線型。這主要是由于地震波垂直向上傳播引起遠場土體純剪應力變化，當地震波傳播至在土-結構接觸界面時發生反射與透射現象，使地下結構周圍場地土的受力狀態更為復雜，從而在土單元體中形成比遠場土單元體更為復雜的動應力路徑，同時產生動剪應力和動正應力，兩者不同的耦合作用形成了斜橢圓型和斜直線型應力路徑。

圖7給出了隧道結構周圍場地不同位置處的超孔壓增長時程曲線，并將本文模型和修正Byrne模型的結果進行了對比。從圖中可以看出，本文模型和修正Byrne模型計算的遠場P1的超孔壓增長時程曲線基本一致，而在結構周圍P2，P3位置處，本文模型較修正Byrne模型計算得出的超孔壓發展更快，最終的超孔壓值也更大，且超孔壓波動性(即應力孔壓)也更顯著。

圖6 隧道周圍場地動應力路徑Fig.6 Predicted stress paths of the soil elements around the tunnel

圖7 隧道周圍場地超孔壓發展時程Fig.7 Time histories of excess pore pressure around the tunnel

圖8描述了隧道結構周圍場地的液化進程。由圖可見，與修正Byrne模型相比，本文模型計算的土層液化范圍明顯更大。結構周圍出現液化，尤其結構底板標高處的土層液化范圍較大，在底板角度位置出現了明顯的液化集中區，且頂板和底板最先發生液化，隨后液化區域沿著結構周圍逐漸擴展形成整片區域。隧道結構的出現改變了自由場的液化特性。距離結構較遠區域，液化最先發生在場地表面，并逐漸向下擴展。本文模型和Byrne模型計算得出的最大液化深度達到了隧道結構底板下4 m和2 m。在距離結構100 m處，在地震作用結束后，液化深度相同，均達到1.5 m。

圖8 隧道周圍場地漸進液化云圖Fig.8 Development of liquefaction zone around the tunnel

結合前文室內試驗結果可知，剪應力和正應力差的耦合剪切作用顯著加大了結構孔壓的累積速度，加劇了土骨架模量的衰退，進而導致液化區域的擴大。綜合圖6，7和8可以看出，基于直剪試驗數據建立的修正Byrne模型僅能夠較好地反映遠場P1的動力響應，而低估了隧道結構周圍P2和P3的超孔壓發展和液化進程。這與3.1節對于室內試驗模擬的規律吻合較好。

4 結 論

本文采用適宜于將一維應力-應變關系向三維空間擴展的等效剪應變算法及加卸載判據，基于Masing法則構造了三維應力空間中的Davidenkov本構模型，補充了剪應力和正應力差耦合作用下剪切-體積應變耦合模型作為Biot動力固結方程中計算超靜孔隙水壓力的源項，基于FLAC3D軟件平臺，建立了適宜于分析可液化地基中土-地下結構相互作用問題的有效應力算法，研究了地下結構所處可液化地基的地震反應規律，可以得到以下主要結論：

(1) 飽和砂土液化特性與循環加載模式和應力路徑密切相關。對飽和砂土不排水循環扭剪試驗模擬結果表明：相較于單一循環直剪試驗結果建立的修正Byrne模型，本文方法可以更為合理表征復雜應力路徑下飽和砂土的孔壓發展規律及液化過程。

(2) 土-地下結構動力相互作用顯著改變了地下結構周圍場地地震應力以及所形成的應力路徑。地震波垂直向上傳播引起遠場土體保持正應力不變條件下規則的往復水平剪應力，而地下結構周圍土體處于往復剪應力和正應力差耦合剪切狀態。

(3) 結合對于室內試驗模擬的規律可知，采用修正Byrne模型可以較好預測地下結構遠場的動力響應，卻低估了結構周圍場地的超孔壓累積速度和液化區域，本文方法能夠較好反映土-地下結構動力相互作用對結構周圍場地動力響應的影響。