鄒德高，劉京茂，汪玉冰，李俊超，李 多，陳 濤，王 鋒

（1.大連理工大學(xué)水利工程學(xué)院，遼寧省大連市 116024；2.大連理工大學(xué)海岸和近海工程國家重點實驗室，遼寧省大連市 116024；3.浙江大學(xué)軟弱土與環(huán)境土工教育部重點實驗室，浙江杭州市 310058；4.四川華能瀘定水電有限公司，四川省成都市 610072；5.中國電建集團成都勘測設(shè)計研究院有限公司，四川省成都市 610072）

0 引言

土石壩因其良好的地基適應(yīng)性，是深厚覆蓋層上建壩的首選壩型[1]。根據(jù)我國水電能源戰(zhàn)略布局的需求，我國建設(shè)和規(guī)劃了一批坐落在覆蓋層上的高土石壩工程，例如：瀑布溝（覆蓋層深度h=70m）、長河壩（h=76.5m）、阿爾塔什（h=90m）、冶勒（h=420m）、旁多（h=424m）等[1，2]，這些高壩很多位于地震頻發(fā)、地質(zhì)復(fù)雜的強震區(qū)且覆蓋層存在易液化土層，強震作用下覆蓋層地基變形和穩(wěn)定是評價深厚覆蓋層上高土石壩抗震安全的重要內(nèi)容。

有限元數(shù)值分析方法是目前研究深厚覆蓋層上高土石壩抗震安全的重要手段。有限元數(shù)值分析主要是通過計算獲得土工構(gòu)筑物的地震加速度、動應(yīng)力或振動孔壓等變量的空間分布，從而為土工構(gòu)筑物的抗震穩(wěn)定性評價或地基液化分析提供依據(jù)。有限元分析土體或壩體的地震動力響應(yīng)時，土的動力非線性性質(zhì)可以采用等價黏彈性模型及彈塑性模型，相應(yīng)地將這些本構(gòu)模型與有限元結(jié)合，發(fā)展了等價線性分析方法以及彈塑性動力分析方法[3]。其中基于等價黏彈性模型的等價線性分析方法理論簡單、應(yīng)用方便，而且在參數(shù)確定和應(yīng)用方面積累了比較豐富的試驗資料和工程經(jīng)驗，已為工程界普遍接受；另一種基于彈塑性模型的非線性分析方法能夠較好地模擬土體的實際反應(yīng)，并能夠直接計算土工構(gòu)筑物的地震變形和孔壓發(fā)展過程，在理論上更為嚴密和合理[4，5]。但已有的商業(yè)軟件集成的大都是理想彈塑性等經(jīng)典彈塑性本構(gòu)模型[6]，這類本構(gòu)模型在屈服面內(nèi)一般認為是彈性的，即不發(fā)生塑性累積變形的，并且往復(fù)循環(huán)加載時也無法計算塑性變形，因此這類型模型無法反映土體卸載體縮引起的孔壓增加以及循環(huán)往復(fù)加載引起的孔隙水壓力增加（剪縮引起的正孔壓）或減?。裘浺鸬呢摽讐海6]，導(dǎo)致這類模型難以合理的描述地基變形和孔壓發(fā)展規(guī)律，不利于準確評價深厚覆蓋層上高土石壩的抗震安全，因此必須開發(fā)適用于深厚覆蓋層上高土石壩抗震安全評價的分析軟件和本構(gòu)模型。

本文基于大連理工大學(xué)自主開發(fā)的巖土工程高性能非線性有限元分析軟件GEODYNA，考慮飽和多孔介質(zhì)Biot動力固結(jié)理論，采用大連理工大學(xué)靜動統(tǒng)一實用化彈塑性本構(gòu)模型[9]，對地震作用下西部某水電站覆蓋層地基加速度和孔隙水壓力進行了數(shù)值模擬，驗證了開發(fā)的本構(gòu)模型和軟件的合理性和可靠性。研究成果可為評價深厚覆蓋層上高土石壩抗震安全提供可靠的評價技術(shù)。

1 離心機振動臺試驗

1.1 模型簡介

試驗以西部某水電站工程為研究對象，試驗研究采用浙江大學(xué)研制的土工離心機[10]。試驗?zāi)Ｐ拖洳捎茂B環(huán)式，其內(nèi)部尺寸為：740 mm×340 mm×425 mm（長度 × 寬度 × 高度），試驗時離心加速度為50g。由于覆蓋層地基斷面材料比較復(fù)雜，為了便于離心機試驗研究，對地基材料分區(qū)進行了簡化，簡化后的地基共2層（見圖1）。試驗上層礫石層選用配置的福建標準砂，厚度為100 mm，試驗下層砂土選用現(xiàn)場砂，厚度為300mm。模型內(nèi)部均設(shè)置了9個水平向加速度傳感器（編號A1-A9）和9個孔壓計（編號為P1-P9），加速度傳感器和孔壓計的位置相同（見圖1）。

圖1 模型試驗情況Figure 1 Illustration of model test information

離心機振動臺試驗輸入地震是硬梁包水電站場地譜地震波，輸入地震波僅為水平向。如圖2所示為A0測點實測的輸入加速度時程（模型試驗中實測值為其50倍）。在后文數(shù)值模擬中，以實測地震波作為計算輸入地震動。

圖2 模型試驗中實測水平地震波Figure 2 Time history curve of seismic acceleration

1.2 計算模型簡介

如圖3所示是本次計算采用的計算網(wǎng)格圖（基于模型實際尺寸建模），單元采用四邊形等參單元。在后文進行動力試驗時，地震動采用底部輸入，為了模擬疊環(huán)式模型箱，對地基兩側(cè)節(jié)點的位移自由度綁定進行模擬分析。試驗過程中地基土均為飽和土，為了模擬孔壓的時程變化過程，因此在后文進行動力分析時采用基于動力固結(jié)方程和彈塑性模型的有效應(yīng)力方法進行數(shù)值模擬。

圖3 模型試驗有限元網(wǎng)格（單位：mm）Figure 3 Finite element model （Unit：mm）

2 數(shù)值模擬方法與理論

2.1 動力流固耦合分析方法簡介

基于Biot動力固結(jié)方程的有效應(yīng)力分析方法是分析飽和土液化變形問題的有效途徑。將Biot動力固結(jié)方程進行空間域離散，并引入系統(tǒng)阻尼并寫成矩陣形式[11，12]，可得：

式中：M——土體的質(zhì)量矩陣；

C——阻尼矩陣；

Q——耦合矩陣；

f(1)——土體的荷載向量；

Mf——流體的質(zhì)量矩陣；

S——流體的壓縮矩陣；

H——流體的滲透矩陣；

f(2)——流體的荷載向量。

本文在計算阻尼矩陣時阻尼比取0.05。大連理工大學(xué)在飽和多孔介質(zhì)彈塑性動力框架上，將巖土工程有限元分析過程進行了類的抽象和封裝，集成了土工構(gòu)筑物和地基的動力有效應(yīng)力分析方法，自主開發(fā)了巖土工程高性能非線性有限元分析軟件GEODYNA。

2.2 地基土材料本構(gòu)簡介

土體本構(gòu)模型的合理性直接影響數(shù)值模擬結(jié)果的可靠性。大連理工大學(xué)在壓力相關(guān)和狀態(tài)相關(guān)靜動統(tǒng)一本構(gòu)模型的基礎(chǔ)上[13，14]，針對土體關(guān)鍵力學(xué)特性—應(yīng)力歷史相關(guān)性模擬存在的問題對塑性模量以及本構(gòu)模型實用性等進行了進一步的發(fā)展[13，14]。改進模型基于比例記憶理論實現(xiàn)了不規(guī)則加載過程中重要應(yīng)力歷史的記憶，更合理的描述了復(fù)雜循環(huán)加載條件下土體應(yīng)力應(yīng)變關(guān)系。該模型參數(shù)的物理意義大都明確，且可根據(jù)常規(guī)試驗直接確定。作者等（2018）[9]采用已有復(fù)雜應(yīng)力路徑條件下（包括，等小主應(yīng)力、等平均主應(yīng)力、等大主應(yīng)力、等應(yīng)力比條件下的單向循環(huán)和雙向循環(huán)）的試驗結(jié)果對模型進行了驗證，結(jié)果表明，改進模型可以較好地反映不規(guī)則加載條件下土體的循環(huán)滯回、循環(huán)硬化等變形特性，尤其是可以考慮復(fù)雜循環(huán)不規(guī)則加載歷史對后期加載變形的影響。同時基于第十五屆國際大壩數(shù)值分析標準研討會發(fā)布的Menta面板堆石壩的動力反應(yīng)分析案例，驗證了本構(gòu)模型在模擬筑壩堆石料力學(xué)特性方面的合理性[15]。提出模型同樣適用于砂土，表1給出了提出模型的本構(gòu)模型參數(shù)，其是根據(jù)離心機模型試驗中現(xiàn)場砂和福建標準砂土的三軸試驗成果標定獲得。

表1 本構(gòu)模型參數(shù)Table 1 Model parameters of soil

3 計算與試驗結(jié)果對比分析

3.1 地震加速度

以下數(shù)值模擬和試驗結(jié)果，均是根據(jù)縮尺規(guī)則將模型結(jié)果等效到原型。圖4給出了A2（距離模型箱底部的高度為0.15m）、A5（距離模型箱底部的高度為0.25m）、A8（距離模型箱底部的高度為0.35m）三個測點處的加速度時程與實測的對比。從圖4中可以看出：埋深最深測點吻合的最好，埋深最淺的測點吻合較差；位于下側(cè)現(xiàn)場砂土層中的A2、A5測點，計算和實測加速度幅值大小變化規(guī)律基本吻合；位于上側(cè)福建標準砂土層的A8測點，在約6s以前計算和實測加速度幅值大小基本一致，在6s以后計算加速度幅值明顯小于實測值。

圖4 計算和實測加速度時程對比Figure 4 Comparison between measured and simulated acceleration time histories

圖5給出了A0、A8點計算和實測傅里葉譜的結(jié)果，可以看到：實際輸入地震動的主頻約為1.6Hz，其與模型基頻（震前約2.0Hz）存在一定的差異，地震過程中因孔壓上升導(dǎo)致地基模量降低導(dǎo)致模型基頻略有降低，震后的基頻約為1.7Hz。在頻率1.7Hz位置（與地震動主頻接近），計算和實測的傅里葉譜幅值均出現(xiàn)了峰值，但兩者傅里葉譜值相差不大。在1.0Hz附近區(qū)域，計算的傅里葉譜極值大于實測值，這與計算和實測的加速度時程峰值規(guī)律并不相同。在大于3.0Hz的范圍內(nèi)，計算的傅里葉幅值要明顯的小于實測值。綜上可知，計算加速度的峰值較小是因為高頻含量較小導(dǎo)致的。造成這種差異的原因可能一方面與數(shù)值分析中未考慮模型箱底部和側(cè)邊界地震波反射有關(guān)[16]；另一方面與數(shù)值分析模擬接近液化低有效應(yīng)力時土體阻尼比大小選擇有關(guān)；此外有學(xué)者認為這種現(xiàn)象與土體的剪脹有關(guān)，但從后文孔壓時程可以看到，計算孔壓和實測值吻合較好[17]。限于篇幅，本文并未對上述因素進行數(shù)值分析探討。

圖5 實測和計算傅里葉譜（單位：g）Figure 5 Comparison between measured and simulated acceleration Fourier spectrum （Unit：g ）

3.2 孔隙水壓力

圖6 給出了P2、P5和P8測點的計算和試驗孔壓發(fā)展過程對比。從圖6中可以看出：位于下側(cè)現(xiàn)場砂土層的P2測點計算和實測孔壓值發(fā)展過程和量值均比較接近；位于下側(cè)現(xiàn)場砂土層的P5測點實測孔壓略小于計算值，但兩者的變化規(guī)律是一致的；位于上側(cè)福建標準砂土層的P8測點實測孔壓略大于計算值，但孔壓累積和消散的過程是基本一致的?？傮w來講，實測與計算的孔壓無論在量值和發(fā)展規(guī)律上都吻合較好。

圖6 孔壓時程對比Figure 6 Comparison between measured and simulated pore water pressure time histories

4 結(jié)論

深厚覆蓋層上建設(shè)高土石壩工程是壩工界正在面臨的難題之一，其中強震作用下覆蓋層地基變形和穩(wěn)定是深厚覆蓋層上高土石壩的抗震安全評價重要內(nèi)容。本文應(yīng)用大連理工大學(xué)自主研發(fā)的巖土工程高性能非線性有限元分析軟件GEODYNA和靜動統(tǒng)一實用化彈塑性本構(gòu)模型，對西部某水電站覆蓋層地基離心機試驗進行了動力彈塑性有效應(yīng)力數(shù)值分析。結(jié)果表明：計算和實測加速度幅值大小變化規(guī)律吻合較好；計算的孔隙水壓力隨地震時間的上升斜率、幅度大小變化及消散趨勢與試驗數(shù)據(jù)基本一致，說明本文采用的本構(gòu)模型能合理描述地震作用下地基土的孔壓發(fā)展規(guī)律。

本文研究成果驗證了開發(fā)的本構(gòu)模型和軟件的合理性和可靠性，在此基礎(chǔ)上聯(lián)合已開發(fā)的土—界面—結(jié)構(gòu)體系跨尺度計算理論和方法[18-20]，可為實際深厚覆蓋層高土石壩地震安全評價和抗震措施設(shè)計提供精細化數(shù)值分析工具。