P2PSand 和Finn 本構模型在水庫土石壩液化模擬中對比

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

（濟南大學 土木建筑學院, 濟南 250022）

0 引言

自17 世紀至今，山東境內及周邊地區已發生15 余次強地震。多數土石壩建于20 世紀80年代，由于當時技術水平落后，修建土石壩時并未考慮壩基砂土液化的影響，導致土石壩在地震發生時發生較大危害，出現滑坡、震陷等問題，因此進行土石壩動力響應研究是非常有必要的。

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

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

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

1 工程概況

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

表1 材料參數Table 1 Material parameters

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

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

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

2 土的本構模型

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

P2PSand 模型由19 個參數控制，其中輸入參數3 個，其余參數均為默認值。

Finn 模型是Martin 等根據試驗提出的，可用于解決土在循環荷載作用下體積應變以及孔隙水壓力的變化規律問題[20]。該模型的實質為在摩爾庫倫模型的基礎上增加了動孔壓的上升模式，并假定動孔壓的上升與塑性體積應變增量相關。

Byrne 在實驗數據基礎上對Finn 模型[21]進行了簡化，簡化后的計算公式為:

式中：Δεvd為砂土體積剪應變增量；εvd為砂土體積剪應變；γ為剪應變；C1、C2為模型參數。經驗計算公式為[21]

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

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

3 動力響應對比計算分析

3.1 壩基液化分析

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

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

圖3 0.1 g 地震波下相對密實度為0.45 時不同本構模型的超孔壓比云圖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 地震波下相對密實度為0.45 時不同本構模型的超孔壓比云圖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 地震波下相對密實度為0.45 時不同本構模型的超孔壓比云圖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 時不同密實度下監測點A 超孔壓比時程曲線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 不同本構模型下水庫土石壩變形云圖可以看出，兩種本構模型的水平位移均大于豎向位移，且背風坡的水平位移大于迎風坡的水平位移。由圖8 中加速度峰值為0.2 g 時不同密實度下監測點B的豎向位移、C的水平位移可知，P2PSand 本構模型下的水平位移遠大于Finn 本構下的水平位移而P2PSand 本構模型下的豎向位移略大于Finn 本構模型下的豎向位移。相對密實度為0.3 和0.45時，P2PSand 本構模型下的水平位移約為Finn 本構模型下水平位移的3 倍，豎向位移約為2 倍；相對密實度為0.7 和0.9 時，P2PSand 本構模型下的水平位移約為Finn 本構模型下水平位移的4 倍，豎向位移約為2.5 倍。由圖6 可知，P2PSand 本構下的液化速度與程度均大于Finn 模型，因此P2PSand 本構下的位移也大于Finn 模型。

圖7 不同本構模型下水庫土石壩位移變形云圖（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 時不同密實度下監測點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 穩定性分析

由圖9 可知，水庫土石壩在兩種本構模型下壩坡均出現兩條明顯的滑動面，穩定性均降低，且處于P2PSand 本構模型下的水庫土石壩最大剪切應變增量約為Finn 本構模型下水庫土石壩最大剪切應變增量的2 倍。P2PSand 本構模型下的水庫土石壩最大剪切應變增量最大值位于壩基左側，Finn 本構模型下水庫土石壩最大剪切應變增量值位于壩頂。

4 結語

通過對兩種不同本構的數值模擬可得：

1）地震作用下壩基可液化粉砂層會發生局部液化，兩種本構下的液化規律相同，P2PSand 本構模型下局部液化區域為重度液化，而Finn 本構模型下液化區域為輕微液化。

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

3）不同本構模型下剪切應變最大值出現的部位不同，均出現兩條明顯的滑動面，穩定性降低。