牛金帝，張西文，邱 宇，扈 萍，陳子俊

（濟(jì)南大學(xué) 土木建筑學(xué)院, 濟(jì)南 250022）

0 引言

自17 世紀(jì)至今，山東境內(nèi)及周邊地區(qū)已發(fā)生15 余次強(qiáng)地震。多數(shù)土石壩建于20 世紀(jì)80年代，由于當(dāng)時(shí)技術(shù)水平落后，修建土石壩時(shí)并未考慮壩基砂土液化的影響，導(dǎo)致土石壩在地震發(fā)生時(shí)發(fā)生較大危害，出現(xiàn)滑坡、震陷等問題，因此進(jìn)行土石壩動(dòng)力響應(yīng)研究是非常有必要的。

土體的動(dòng)力分析方法根據(jù)本構(gòu)模型可分為基于等價(jià)粘彈性的等效線性分析方法和基于(黏)彈塑性模型的非線性分析方法。等效線性分析方法是用一個(gè)等效的剪切模量和阻尼比代替所有不同應(yīng)變幅值下的剪切模量和阻尼比，將非線性問題轉(zhuǎn)化為線性問題，利用頻域線性波動(dòng)方法求解[1]。但是等效線性分析方法存在不能反映土層中真實(shí)的運(yùn)動(dòng)過程，在地震動(dòng)輸入較大時(shí)計(jì)算誤差較大，可能出現(xiàn)死循環(huán)現(xiàn)象等問題。若想得到真實(shí)的土體地震反應(yīng)，宜采用基于(黏)彈塑性模型的非線性分析方法。非線性分析方法將土體視作為（黏）彈塑性變形材料，模型由初始加荷曲線、移動(dòng)的骨干曲線和開放的滯回圈組成。該模型能更好地模擬殘余應(yīng)變，進(jìn)行動(dòng)力分析時(shí)可以直接計(jì)算殘余變形，在動(dòng)力分析中可以隨時(shí)計(jì)算切線模量并進(jìn)行非線性計(jì)算，這樣得到的動(dòng)力響應(yīng)過程能夠更好地接近實(shí)際情況[2]。

在有關(guān)巖土地震的工程實(shí)踐中，基于(黏)彈塑性模型的非線性分析方法被廣泛運(yùn)用于模擬砂類材料的液化現(xiàn)象，例如MD97 模型[3]、DM04 模型[4]，以及在此基礎(chǔ)上拓展的Sanisand 模型[5]、UBCsand模型[6]、PM4sand 模型[7]等。DM04 模型需要高準(zhǔn)確性的實(shí)驗(yàn)室測試數(shù)據(jù)且不具備液化后剛度退化的機(jī)制，因此在其基礎(chǔ)上進(jìn)行改造并建立了雙塑性面砂土模型P2PSand 模型。該模型具有DM04 模型理論的穩(wěn)健性和簡明性，以及UBCsand 模型和PM4sand 模型面向?qū)嵺`的特性。P2Psand 適用于一般的三維條件，P2PSand 模型參數(shù)適用于標(biāo)準(zhǔn)循環(huán)阻力場砂土，與半經(jīng)驗(yàn)過程中的循環(huán)阻力圖兼容，且該模型允許用戶指定自定義參數(shù)，而非默認(rèn)材料參數(shù)[8-10]。

土體是非線性材料，地震動(dòng)力越大，則土體的非線性特征越強(qiáng)，因此研究土石壩動(dòng)力響應(yīng)應(yīng)采取非線性分析方法，雙塑性面砂土模型Practical Twosurface Plastic Sand (簡稱P2PSand）模型[11-12]與Finn本構(gòu)模型均采用非線性分析方法?；赑2PSand 本構(gòu)模型，Zhao C 等[13-14]在DM04 模型上修改了一些公式，并保留了一般的三維公式，模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)以及不同初始和加載條件下的經(jīng)驗(yàn)關(guān)系吻合較好，從而提高了模型的性能；吳宏[15]等采用三維數(shù)值方法研究穿越不同密實(shí)度狀態(tài)飽和砂土地層的盾構(gòu)隧道的地震響應(yīng)規(guī)律?；贔inn 本構(gòu)模型，胡南雄等[16]對某中型水庫大壩壩基可液化土層進(jìn)行數(shù)值模擬，分析了不同庫水位作用下壩基可液化土層的孔隙水壓力變化過程及壩體典型節(jié)點(diǎn)的動(dòng)力響應(yīng)；嚴(yán)祖文等[17]分析壩基土及大壩整體抗震穩(wěn)定，并根據(jù)研究成果提供合理可行的防震措施、供水庫加固設(shè)計(jì)與施工應(yīng)用；李伯根等[18]對大后溝尾礦堆積壩進(jìn)行地震動(dòng)力響應(yīng)分析并進(jìn)行液化評價(jià)，分析其動(dòng)力穩(wěn)定性，以指導(dǎo)尾礦庫的安全運(yùn)行。P2PSand 本構(gòu)模型是于2019年新開發(fā)的本構(gòu)，國內(nèi)對該本構(gòu)的應(yīng)用研究較少，因此本文將采用FLAC3D 7.0 內(nèi)置P2PSand 本構(gòu)模型與Finn 本構(gòu)模型進(jìn)行數(shù)值模擬，從而得出不同本構(gòu)模型下的土石壩動(dòng)力響應(yīng)特征，找出土石壩最薄弱區(qū)域，為土石壩抗震除險(xiǎn)加固提供理論依據(jù)。

1 工程概況

本文計(jì)算模型以山東某水庫土石壩為背景，該水庫土石的模型及監(jiān)測點(diǎn)分布如圖1 所示。該大壩底端固定，兩側(cè)水平方向固定，豎直方向自由，地震波由底部輸入。計(jì)算監(jiān)測變量為監(jiān)測點(diǎn)A的超孔壓比，監(jiān)測點(diǎn)B的豎向位移，監(jiān)測點(diǎn)C的水平位移等。該水庫土石壩材料參數(shù)如表1 所示。

表1 材料參數(shù)Table 1 Material parameters

圖1 水庫土石壩模型及監(jiān)測點(diǎn)分布示意圖Fig. 1 Reservoir earth-rock dam model and monitoring site distribution diagram

模型邊界設(shè)置為自由場邊界，以此來減少邊界對地震波的反射；材料阻尼設(shè)置為0.1571；選用FLAC3D 建模，設(shè)置土層的材料參數(shù)和邊界條件，于土石壩底部輸入不同峰值強(qiáng)度的地震波，為了保證數(shù)值模擬的計(jì)算速度和精度，通過濾波把大于10 Hz的高頻分量過濾掉，將濾波后的地震動(dòng)作為基底輸入地震動(dòng)。以峰值為0.2 g 的地震波為例，利用SeismoSignal 軟件對原始地震波進(jìn)行濾波和基線校正，并在模型底部水平施加處理后的地震波（圖2）。其中，在FLAC3D 分別設(shè)置P2PSand 和Finn 本構(gòu)模型模擬地震砂土液化現(xiàn)象，在土石壩底部施加峰值強(qiáng)度為0.1 g、0.2 g、0.3 g 強(qiáng)度的地震波，g 為重力加速度，可液化粉砂層相對密實(shí)度依次為0.30、0.45、0.70、0.90。

圖2 利用SeismoSignal 軟件對原始地震波進(jìn)行濾波和基線校正處理后的地震波Fig. 2 Seismic wave after filtering and baseline correction of original seismic wave using SeismoSignal software

2 土的本構(gòu)模型

在實(shí)際工程應(yīng)用中，所采用的本構(gòu)模型不僅要能夠準(zhǔn)確表達(dá)土的應(yīng)力-應(yīng)變關(guān)系，而且還必須具有形式簡單、參數(shù)易于確定等特點(diǎn)。P2PSand 本構(gòu)模型和Finn 本構(gòu)模型同時(shí)具有模型形式簡單與參數(shù)易于確定的優(yōu)勢，其參數(shù)由相對密實(shí)度決定；P2PSand 本構(gòu)模型和Finn 本構(gòu)模型的其他參數(shù)由與相對密實(shí)度有關(guān)的經(jīng)驗(yàn)公式取得，相對密實(shí)度在工程中容易獲取。FLAC3D[19]在進(jìn)行動(dòng)力分析時(shí)所采用的方法為基于(黏)彈塑性模型的非線性分析方法，遵循Biot 流體-機(jī)械耦合孔隙力學(xué)理論的應(yīng)用。

P2PSand 模型由19 個(gè)參數(shù)控制，其中輸入?yún)?shù)3 個(gè)，其余參數(shù)均為默認(rèn)值。

Finn 模型是Martin 等根據(jù)試驗(yàn)提出的，可用于解決土在循環(huán)荷載作用下體積應(yīng)變以及孔隙水壓力的變化規(guī)律問題[20]。該模型的實(shí)質(zhì)為在摩爾庫倫模型的基礎(chǔ)上增加了動(dòng)孔壓的上升模式，并假定動(dòng)孔壓的上升與塑性體積應(yīng)變增量相關(guān)。

Byrne 在實(shí)驗(yàn)數(shù)據(jù)基礎(chǔ)上對Finn 模型[21]進(jìn)行了簡化，簡化后的計(jì)算公式為:

式中：Δεvd為砂土體積剪應(yīng)變增量；εvd為砂土體積剪應(yīng)變；γ為剪應(yīng)變；C1、C2為模型參數(shù)。經(jīng)驗(yàn)計(jì)算公式為[21]

表2 給出可液化粉砂層的P2PSand 本構(gòu)模型與Finn 本構(gòu)模型參數(shù)。由表2 可知，2 個(gè)模型參數(shù)循環(huán)因子Kc、彈性剪切模量Gr、模型參數(shù)C1、C2均由相對密實(shí)度dr確定。

表2 可液化粉砂層的P2PSand 本構(gòu)模型與Finn 本構(gòu)模型參數(shù)Table 2 P2PSand constitutive model and Finn constitutive model parameters of liquefiable silty sand layer

3 動(dòng)力響應(yīng)對比計(jì)算分析

3.1 壩基液化分析

為了更直觀地反映水庫土石壩砂土液化的程度，定義超孔壓比excess pore water pressure ratio（epwpr）描述壩基液化程度，對監(jiān)測點(diǎn)A的超孔壓比進(jìn)行監(jiān)測，計(jì)算公式如下：

當(dāng)epwpr=1時(shí)，說明土體完全液化。工程上當(dāng)epwpr＞0.6時(shí)，飽和砂土地基將產(chǎn)生失效。因此，本文將epwpr=0.6的單元定義為準(zhǔn)液化單元，這也是考慮到土石壩的穩(wěn)定性而做的一種安全儲(chǔ)備。圖3～5 分別為加速度峰值為0.1 g、0.2 g、0.3 g 地震波輸入、相對密實(shí)度為0.45 時(shí)不同本構(gòu)模型的超孔壓比云圖，圖6 為加速度峰值為0.2 g 時(shí)不同密實(shí)度下監(jiān)測點(diǎn)B超孔壓比時(shí)程曲線圖。由圖3～6 可知，當(dāng)加速度峰值為0.1 g 時(shí)，水庫土石壩壩基兩側(cè)發(fā)生液化，壩體下部壩基并未發(fā)生液化，隨著加速的峰值增大，水庫土石壩壩體下部壩基出現(xiàn)砂土液化現(xiàn)象?；赑2PSand 本構(gòu)下的可液化粉砂層的地震砂土液化程度與范圍略大于Finn 本構(gòu)模型下的程度與范圍。兩種本構(gòu)模型下，可液化粉砂層液化的規(guī)律相同：t＜4 s 時(shí)，超孔壓比增長較??；4 s＜t＜16 s時(shí)，超孔壓比增長逐漸增大；當(dāng)t＞16 s 時(shí)，超孔壓比增長趨于穩(wěn)定。由圖3～5 可以看出，兩種本構(gòu)模型下，可液化粉砂層液化的面積與范圍相差無幾，但液化程度相差較大；P2PSand 本構(gòu)模型下壩基兩側(cè)超孔壓比較大數(shù)值已接近1.0，處于完全液化，而Finn 本構(gòu)模型下壩基兩側(cè)超孔壓比最大數(shù)值約為0.7，處于中等液化。由于P2PSand 模型采用了DM04模型理論的穩(wěn)健性和簡明性，以及UBCsand 模型和PM4sand 模型面向?qū)嵺`的特性，因此計(jì)算結(jié)果更貼近于實(shí)際。

圖3 0.1 g 地震波下相對密實(shí)度為0.45 時(shí)不同本構(gòu)模型的超孔壓比云圖Fig. 3 The cloud image of over pore pressure ratio under 0.1 g seismic wave action for different constitutive models with the relative density of 0.45

圖4 0.2 g 地震波下相對密實(shí)度為0.45 時(shí)不同本構(gòu)模型的超孔壓比云圖Fig. 4 The cloud image of over pore pressure ratio under 0.2 g seismic wave action for different constitutive models with the relative density of 0.45

圖5 0.3 g 地震波下相對密實(shí)度為0.45 時(shí)不同本構(gòu)模型的超孔壓比云圖Fig. 5 The cloud image of over pore pressure ratio under 0.3 g seismic wave action for different constitutive models with the relative density of 0.45

圖6 加速度峰值為0.2 g 時(shí)不同密實(shí)度下監(jiān)測點(diǎn)A 超孔壓比時(shí)程曲線Fig. 6 The time history curve of excess pore pressure ratio for monitoring point A under different density is obtained when the peak acceleration is 0.2 g

3.2 壩體變形分析

由圖7 不同本構(gòu)模型下水庫土石壩變形云圖可以看出，兩種本構(gòu)模型的水平位移均大于豎向位移，且背風(fēng)坡的水平位移大于迎風(fēng)坡的水平位移。由圖8 中加速度峰值為0.2 g 時(shí)不同密實(shí)度下監(jiān)測點(diǎn)B的豎向位移、C的水平位移可知，P2PSand 本構(gòu)模型下的水平位移遠(yuǎn)大于Finn 本構(gòu)下的水平位移而P2PSand 本構(gòu)模型下的豎向位移略大于Finn 本構(gòu)模型下的豎向位移。相對密實(shí)度為0.3 和0.45時(shí)，P2PSand 本構(gòu)模型下的水平位移約為Finn 本構(gòu)模型下水平位移的3 倍，豎向位移約為2 倍；相對密實(shí)度為0.7 和0.9 時(shí)，P2PSand 本構(gòu)模型下的水平位移約為Finn 本構(gòu)模型下水平位移的4 倍，豎向位移約為2.5 倍。由圖6 可知，P2PSand 本構(gòu)下的液化速度與程度均大于Finn 模型，因此P2PSand 本構(gòu)下的位移也大于Finn 模型。

圖7 不同本構(gòu)模型下水庫土石壩位移變形云圖（0.2 g，dr=0.45，t=20s）Fig. 7 Displacement and deformation cloud map of reservoir earth-rock dam under different constitutive models（0.2g,dr=0.45，t=20s）

圖8 加速度峰值為0.2 g 時(shí)不同密實(shí)度下監(jiān)測點(diǎn)B 的豎向位移C 的水平位移Fig. 8 When the peak acceleration is 0.2 g the vertical displacement of the monitoring point A and the horizontal displacement of the monitoring point B under different compactness are obtained

3.3 穩(wěn)定性分析

由圖9 可知，水庫土石壩在兩種本構(gòu)模型下壩坡均出現(xiàn)兩條明顯的滑動(dòng)面，穩(wěn)定性均降低，且處于P2PSand 本構(gòu)模型下的水庫土石壩最大剪切應(yīng)變增量約為Finn 本構(gòu)模型下水庫土石壩最大剪切應(yīng)變增量的2 倍。P2PSand 本構(gòu)模型下的水庫土石壩最大剪切應(yīng)變增量最大值位于壩基左側(cè)，F(xiàn)inn 本構(gòu)模型下水庫土石壩最大剪切應(yīng)變增量值位于壩頂。

4 結(jié)語

通過對兩種不同本構(gòu)的數(shù)值模擬可得：

1）地震作用下壩基可液化粉砂層會(huì)發(fā)生局部液化，兩種本構(gòu)下的液化規(guī)律相同，P2PSand 本構(gòu)模型下局部液化區(qū)域?yàn)橹囟纫夯?，而Finn 本構(gòu)模型下液化區(qū)域?yàn)檩p微液化。

2）Finn 本構(gòu)模型下的水平位移與豎向位移均小于P2PSand 本構(gòu)模型下的水平與豎向位移，豎向位移在兩種本構(gòu)下相差較小。

3）不同本構(gòu)模型下剪切應(yīng)變最大值出現(xiàn)的部位不同，均出現(xiàn)兩條明顯的滑動(dòng)面，穩(wěn)定性降低。