王 蘇，路德春，杜修力

（北京工業(yè)大學 城市與工程安全減災教育部重點實驗室，北京 100124）

1 引 言

土是由土顆粒、孔隙水和孔隙氣體組成的摩擦類材料，其力學特性不僅取決于土自身材料的組成成份，而且與荷載的作用密切相關，包括兩方面，一是初始應力狀態(tài)，一般為自然狀態(tài)土的自重應力狀態(tài)；二是附加荷載的應力路徑[1]。相應地，土表現(xiàn)出兩種基本力學特性，即壓硬性和剪脹性，壓硬性指土的抗剪強度和剛度隨著約束壓力的增大而增大；剪脹性指土在附加剪應力作用下產(chǎn)生體積膨脹或收縮的特性。

地下結構賦存于巖土介質(zhì)中，其地震響應與地面結構的慣性效應不同，主要受圍巖土體的變形控制。土體的變形取決于兩部分荷載，一是地震作用的附加荷載；二是土體的自重應力。數(shù)值模擬時，土體的自重應力狀態(tài)與地震荷載的靜-動力耦合作用一般通過人工邊界條件的設置實現(xiàn)，主要有3種方法：第1種方法為靜-動力統(tǒng)一人工邊界方法[2]，這種方法是在黏彈性動力人工邊界的基礎上，通過改變?nèi)斯み吔鐥l件的剛度系數(shù)來實現(xiàn)靜-動力耦合分析，但對剛度系數(shù)的改變尚缺乏明確的物理解釋。第2種方法為動力松弛法[3]，該方法認為，土體的自重應力狀態(tài)，是土體在階躍荷載作用下從原始非受載狀態(tài)發(fā)生振動達到平衡時的穩(wěn)定狀態(tài)。動力松弛法將靜力問題看作動力問題的一部分，為階躍載荷作用下動力響應的穩(wěn)定狀態(tài)，而不考慮動力過程。它的優(yōu)點是避免了直接求解聯(lián)立方程組，對任何系統(tǒng)都有確定的求解步驟，計算簡單并便于程序?qū)崿F(xiàn)。用于地下結構-圍巖系統(tǒng)地震反應分析時，動力松弛法要求土體的本構模型適用于從非受載狀態(tài)到地震荷載作用的全過程，目前還沒有統(tǒng)一描述從小應變到大應變?nèi)秶耐恋姆蔷€性本構模型[4]，若采用土的常規(guī)本構模型則不能合理地統(tǒng)一考慮土體的自重應力狀態(tài)與地震荷載作用。第3種方法是首先進行靜力分析，然后將靜力分析得到的支座反力在靜力分析步向動力分析步轉(zhuǎn)化過程中，以集中力形式施加于側向和底部邊界面支座對應的單元節(jié)點上，以此作為動力分析的初始條件，實現(xiàn)靜力邊界條件與動力邊界條件的統(tǒng)一。

本文將地下結構-圍巖系統(tǒng)的靜-動力耦合作用歸結為考慮土體自重應力狀態(tài)影響的地下結構-圍巖系統(tǒng)地震響應問題，通過靜力分析獲得初始地應力，將其作為集中力作用于單元節(jié)點上，利用杜修力等[5－6]提出的應力型黏彈性人工邊界條件，結合有限元分析軟件ABAQUS，實現(xiàn)地下結構地震響應靜-動力耦合模擬。利用本文的模擬方法，通過自定義的混凝土材料本構模型，將極限剪應變作為破壞標準，采用生死單元模擬地下結構的地震破壞過程，分析了結構埋深即土體的初始應力狀態(tài)對其抗震性能的影響。

2 地震波垂直入射黏彈性人工邊界

對近場波動問題進行數(shù)值分析，將未知的散射場或邊界外行場用總場減去自由場或邊界入射場（用上標i表示）表示，得到含外源作用的人工邊界面l節(jié)點i方向的運動方程為[6]

式中：δij為克羅內(nèi)克符號；ml為節(jié)點l的集中質(zhì)量；klikj為節(jié)點k方向j對于節(jié)點l方向i的剛度系數(shù)；clikj為節(jié)點k方向j對于節(jié)點l方向i的阻尼系數(shù)；Al為人工邊界面上l節(jié)點的影響面積；fli為入射地震波作用下l節(jié)點i方向的等效節(jié)點力；Fli為在結點l方向i處截去的無限遠場對有限近場的作用力；式（2）括號內(nèi)右邊3項分別為克服彈簧、阻尼器及介質(zhì)所需的抗力；對于三維問題，n=3，下標i和 j分別為1、2、3，分別對應于直角坐標系的x、y和z軸。Kli、Cli為l節(jié)點的人工邊界彈簧和阻尼元件參數(shù)，如圖1所示，對于三維人工邊界，彈簧-阻尼元件參數(shù)為

對于二維人工邊界，彈簧-阻尼元件參數(shù)為

式中：λ為介質(zhì)的拉梅常數(shù)；G為剪切模量；ρ為介質(zhì)密度；cp和cs分別為P波和S波的波速；r為結構幾何中心到該人工邊界點所在邊界線或面的距離；參數(shù)A和B的值分別建議為A=0.8，B=1.0。

圖1 三維黏彈性人工邊界Fig.1 Three-dimensional viscoelastic artificial boundary

3 靜-動力耦合邊界條件模擬方法

通過靜力分析求得土體與結構各單元的內(nèi)力及支座反力，將所得內(nèi)力作為初始內(nèi)力施加到各單元上，邊界處采用黏彈性人工邊界條件，并將靜力分析得到的支座反力施加到黏彈性邊界的節(jié)點上，再進行一次靜力分析，此次所得結果便是地下結構-圍巖系統(tǒng)的初始應力狀態(tài)，并且無初始變形，通常稱為地應力平衡。將初始地應力作為動力分析的初始狀態(tài)，之后再進行地震響應分析，即可得到考慮土體初始自重應力狀態(tài)影響的地下結構地震響應。

基于有限元分析軟件 ABAQUS，編制了靜-動力耦合人工邊界模擬子程序，程序流程如圖2所示。為驗證本文提出的靜-動力耦合邊界模擬方法以及相應的子程序，采用數(shù)值算例進行驗證。

圖2 靜-動力耦合人工邊界施加流程圖Fig.2 Flow chart of artificial boundary setting instatic-dynamic coupling

算例有限元模型如圖3所示，土體區(qū)域為250 m×250 m×250 m，單元為10 m×10 m×10 m的六面體實體單元。土體的材料參數(shù)分別為，密度ρ=2 g/cm3，彈性模量E=117 MPa，泊松比υ=0.3，剪切波速cS=150 m/s。輸入的脈沖為垂直入射S波，時間步長為0.005 s。位移與速度時程曲線如圖4所示。

圖3 有限元模型及黏彈性人工邊界Fig.3 Finite element model and viscoelastic artificial boundary

圖4 位移與速度時程曲線Fig.4 Displacement and velocity time histories

靜力分析后自重應力達到了K0固結狀態(tài)，豎向應力符合 σv=ρg z 關系，z為距地表的埋深；水平向應力符合 σh= K0σv關系， K0= υ (1 - υ)。圖5為靜力分析后土體的豎向位移云圖，各位置處的豎向位移均小于在4×10-5m，該結果取決于土體的彈性模量E與柏松比υ，忽略圖5所示的豎向固結變形對地下結構的地震反應對分析幾乎無影響。圖 6是采用黏彈性邊界后土體底面和地表位移數(shù)值解與解析解的比較，表明本文設置的黏彈性邊界條件可較好地模擬邊界處的能量輻射，并實現(xiàn)地震動的輸入。

圖5 初始時刻土體豎向位移Fig.5 Vertical displacement of soil at initial time

圖6 底邊界與地表位移時程解析解與數(shù)值計算結果的比較Fig.6 Bottom boundary and surface displacement time analytical solutions and numerical calculation results

4 地下結構-圍巖系統(tǒng)靜-動力耦合模擬

地下結構的埋深不同，主要反映了結構埋深處土體的自重應力狀態(tài)不同，埋深越深，土體的自重應力越大。本文分別模擬了結構頂板距地表的埋深為10、15、20 m共3個二維算例。地下結構-圍巖系統(tǒng)有限元模型中，結構模型水平向長為20 m，高為6 m，內(nèi)有2根截面寬為0.4 m的柱子，頂板及左右墻厚為0.4 m，底板厚為0.6 m；土體范圍為結構兩側各為20 m，結構底板距計算區(qū)域下邊界20 m，3個算例結構頂板距地表分別為10、15、20 m。土體單元網(wǎng)格為 1 m×1 m；結構大部分單元網(wǎng)格為0.3 m×0.2 m，局部細劃，結構埋深10 m的有限元模型如圖7所示。土體與結構均采用四面體實體單元，土-結構相互作用界面采用接觸面對罰函數(shù)法，摩擦系數(shù)取0.6。

圖7 地下結構-圍巖系統(tǒng)數(shù)值計算模型Fig.7 Numerical calculation model of underground structure-rock system

地震動輸入為修正后的 Kobe波，位移峰值與速度峰值均放大2倍，時間步長為0.005 s，輸入的地震波只考慮水平x向SV波，位移時程和速度時程曲線如圖8所示。

土體采用 Drucker-Prager本構模型，模型參數(shù)如表1所示。鋼筋混凝土結構等效為均一材料，并基于ABAQUS自定義彈塑性損傷模型[7]，損傷因子D與等效塑性應變?yōu)橹笖?shù)關系，表達式為

圖8 地震動的位移時程與速度時程曲線Fig.8 Displacement and velocity time histories of seismic waves

表1 土的物理力學參數(shù)Table 1 Physical and mechanical parameters of soil

采用生死單元法模擬地下結構-圍巖系統(tǒng)的地震響應。利用ωD作為單元刪除的判斷準則，ωD為隨著損傷因子D單調(diào)增加的狀態(tài)變量，表示單元的損傷程度，始終 ΔωD≥0，當ωD=1時，單元被“刪除”?！皠h除”后的單元，其材料的變形模量趨于0，與未刪除的相鄰單元保持變形協(xié)調(diào)，位移連續(xù)。在有限元模型中對“刪除”的單元進行隱藏，不進行網(wǎng)格重新劃分，是一種假想的刪除。數(shù)值模擬結果如圖9～11所示，DUCTCRT即為ωD的值。

圖9 埋深10 m時地下結構的地震破壞過程Fig.9 Earthquake failure process of underground structure with the depth of 10 m

圖11 埋深20 m時地下結構的地震破壞過程Fig.11 Earthquake failure process of underground structure with the depth of 20 m

對于埋深10 m的結構，其地震破壞首先從中柱與底板的連接處開始，隨著塑性應變的增大，損傷進一步積累，另一中柱與底板的連接處也開始破壞，最后是側墻與底板連接處破壞。對于埋深15 m的結構，其地震破壞的規(guī)律與埋深10 m時相似，區(qū)別在于埋深15 m時結構破壞的時間整體早于埋深10 m時結構破壞的時間。對于埋深20 m的結構，其地震破壞同樣首先從中柱與底板的連接處開始，然后是中柱與頂板的連接處破壞，最后是邊墻與底板連接處、邊墻與頂板連接處破壞，每一破壞處的破壞區(qū)域以及結構整體的破壞程度均顯著大于埋深10 m和15 m時，結構地震破壞的時間與埋深10 m時接近。將3個埋深的地下結構破壞過程整理為圖12，埋深15 m時的地下結構單元初始破壞的時刻在2.252 s，比埋深10 m和埋深20 m的地下結構單元初始破壞的時刻早，而且結構破壞的時刻也比另外兩種埋深早，說明此地下結構形式，在地震作用下存在一個最不利埋深，在該埋深時結構最容易發(fā)生地震破壞。埋深20 m時的地下結構單元初始破壞的時刻最晚，表明地下結構埋深超過最不利埋深后，隨著埋深的加深結構越不容易發(fā)生地震破壞；另一方面，埋深20 m時的地下結構中柱與頂板的連接處以及側墻與頂板的連接處均發(fā)生了地震破壞，表明深埋地下結構一旦發(fā)生地震破壞，其破壞程度隨著埋深的加深而加重。

圖12 不同埋深時結構的地震破壞過程Fig.12 Seismic damage process of structure at different buried depths

5 結 論

本文從土的力學特性分析入手，將荷載作用分為自重荷載與附加荷載兩部分，并將自重應力作為附加荷載作用的初始狀態(tài)，該思路與當前研究土的力學特性與本構模型的建模方法一致。進而將地下結構-圍巖系統(tǒng)的靜-動力耦合作用歸結為考慮土體自重應力狀態(tài)影響的地下結構-圍巖系統(tǒng)地震響應問題，分別進行靜力與動力分析，通過人工邊界條件的設置實現(xiàn)了地下結構地震響應靜-動力耦合模擬。

（1）本文提出的土體自重應力的分析思路以及相應的人工邊界條件的設置方法，為地下結構的抗震性能分析與評價提供了一種有利工具。

（2）初始地應力顯著地影響地下結構-圍巖系統(tǒng)的動力響應，地下結構對于抗震性能而言存在一個最不利埋深，處于最不利埋深時最容易發(fā)生地震破壞；埋深超過最不利埋深時，埋深越深地下結果越不容易發(fā)生地震破壞，一旦發(fā)生地震破壞，結構的埋深越深，其地震破壞的程度越重。

[1]路德春, 姚仰平. 砂土的應力路徑本構模型[J]. 力學學報, 2005, 37(4): 451－459.LU De-chun, YAO Yang-ping. Constitutive model of sand considering complex stress paths[J]. Acta Mechanica Sinica, 2005, 37(4): 451－459.

[2]劉晶波, 李彬. 三維黏彈性靜-動力統(tǒng)一人工邊界[J].中國科學: E 輯, 2005, 35(9): 966－980.LIU Jing-bo, LI Bin. 3D viscoelastic static-dynamic artificial boundary[J]. Science in China: Serial E, 2005,35(9): 966－980.

[3]杜修力. 工程波動理論和方法[M]. 北京: 科學出版社,2009.

[4]PAPADIMITRIOU A G, BOUCKOVALAS G D.Plasticity model for sand under small and large cyclic strains: A multiaxial formulation[J]. Soil Dynamics and Earthquake Engineering, 2002, 22(3): 191－204.

[5]杜修力, 趙密, 王進廷. 近場波動模擬的一種應力邊界條件[J]. 力學學報, 2006, 38(1): 49－56.DU Xiu-li, ZHAO Mi, WANG Jin-ting. A stress artificial boundary in FEA for near-field wave problem[J]. Acta Mechanica Sinica, 2006, 38(1): 49－56.

[6]杜修力. 局部解耦的時域波分析方法[J]. 世界地震工程, 2000, 16(3): 22－26.DU Xiu-li. A partially decoupling analytical method for wave propagation problems in time domain[J]. World Information on Earthquake Engineering, 2000, 16(3):22－26.

[7]王蘇. 地下結構地震破壞過程彈塑性模擬研究[D]. 北京: 北京工業(yè)大學, 2012.