地下結(jié)構(gòu)地震破壞靜-動(dòng)力耦合模擬研究

王 蘇，路德春，杜修力

（北京工業(yè)大學(xué) 城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100124）

1 引 言

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

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

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

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

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

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

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

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

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

3 靜-動(dòng)力耦合邊界條件模擬方法

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

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

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

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

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

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

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

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

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

4 地下結(jié)構(gòu)-圍巖系統(tǒng)靜-動(dòng)力耦合模擬

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

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

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

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

圖8 地震動(dòng)的位移時(shí)程與速度時(shí)程曲線Fig.8 Displacement and velocity time histories of seismic waves

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

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

圖9 埋深10 m時(shí)地下結(jié)構(gòu)的地震破壞過程Fig.9 Earthquake failure process of underground structure with the depth of 10 m

圖11 埋深20 m時(shí)地下結(jié)構(gòu)的地震破壞過程Fig.11 Earthquake failure process of underground structure with the depth of 20 m

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

圖12 不同埋深時(shí)結(jié)構(gòu)的地震破壞過程Fig.12 Seismic damage process of structure at different buried depths

5 結(jié) 論

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

（1）本文提出的土體自重應(yīng)力的分析思路以及相應(yīng)的人工邊界條件的設(shè)置方法，為地下結(jié)構(gòu)的抗震性能分析與評(píng)價(jià)提供了一種有利工具。

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

[1]路德春, 姚仰平. 砂土的應(yīng)力路徑本構(gòu)模型[J]. 力學(xué)學(xué)報(bào), 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]劉晶波, 李彬. 三維黏彈性靜-動(dòng)力統(tǒng)一人工邊界[J].中國科學(xué): 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]杜修力. 工程波動(dòng)理論和方法[M]. 北京: 科學(xué)出版社,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ìn)廷. 近場(chǎng)波動(dòng)模擬的一種應(yīng)力邊界條件[J]. 力學(xué)學(xué)報(bào), 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]杜修力. 局部解耦的時(shí)域波分析方法[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]王蘇. 地下結(jié)構(gòu)地震破壞過程彈塑性模擬研究[D]. 北京: 北京工業(yè)大學(xué), 2012.