超高土石壩心墻孔壓不均勻分布特征及其機理研究

鄒德高，姜秋婷，劉京茂，金 偉，朱先文

(1.大連理工大學 海岸與近海工程國家重點實驗室，遼寧 大連 116024；2.大連理工大學 水利工程學院，遼寧 大連 116024；3.中國電建集團成都勘測設計研究院有限公司，四川 成都 610072)

1 研究背景

實測資料顯示，高土石壩施工過程中心墻內會產(chǎn)生較高的超靜孔隙水壓力(簡稱孔壓)。例如，小浪底斜心墻壩(壩高154 m)竣工時最大孔壓為1.37 MPa[1]；糯扎渡心墻壩(壩高261.5 m)滿蓄時最大孔壓為2.30 MPa[2]。這些監(jiān)測成果為發(fā)展和驗證高心墻壩安全評價方法、推動科學認識大壩工作性態(tài)提供了寶貴支撐，但已有工程普遍存在有效滲壓計數(shù)目相對較少的問題。在建的高295 m兩河口心墻壩建立了完備的監(jiān)測系統(tǒng)，獲得了豐富的孔壓監(jiān)測資料，目前最大監(jiān)測孔壓已達3.74 MPa[3]，同時發(fā)現(xiàn)孔壓呈顯著不均勻分布的特征，超出了以往的認識，引起了工程界對兩河口大壩心墻安全的擔憂。

為了分析孔壓對工程安全的影響，眾多學者提出并發(fā)展了基于有效應力原理的流固耦合分析方法，主要包括Biot流固耦合法[4-7]和不排水流固耦合法[8-9]兩種，兩者各有優(yōu)缺點。Biot流固耦合法在理論上更嚴密，但面臨著方程組系數(shù)矩陣病態(tài)問題，病態(tài)程度與時間步長、單元尺寸、滲透性等多種因素有關，并且求解效率低、易失敗[10]，相關研究成果主要集中在二維分析。近年來，許多學者[4-5]開始對土石壩進行大型三維Biot流固耦合分析，但相關研究仍較少。不排水流固耦合法雖僅適用于土體飽和度較高、孔隙流體滲透性極低的情況，但該方法可在傳統(tǒng)有限元計算軟件中植入，實現(xiàn)方便，并且求解效率高、穩(wěn)定性好。同時，大量研究[11-12]表明孔隙流體的壓縮性與飽和度密切相關。為反映飽和度對孔壓的影響，曹雪山等[13]在Biot流固耦合法中考慮了孔隙流體的壓縮性，但其使用時仍存在前述問題。此外，心墻孔壓模擬的準確性還取決于采用的本構模型。目前心墻力學特性模擬大都采用鄧肯E-B本構模型，該模型計算的水平位移較實測偏大[14-15]，這會導致心墻孔壓模擬失真。相比之下，彈塑性本構模型在模擬土體力學特性時結果更合理[16]。

綜上，為了真實反映土石壩心墻孔壓與飽和度的耦合效應及其三維分布特征，本文首先將飽和度作為初始材料參數(shù)引入到廣義塑性本構模型中，建立了不排水條件下簡化、高效的流固耦合分析方法，并驗證了該方法在模擬心墻壩填筑孔壓和變形時的適用性。基于心墻實測飽和度隨機分布規(guī)律和三維流固耦合數(shù)值模擬，系統(tǒng)研究了兩河口大壩心墻孔壓分布特征及其形成機理。本文研究可為分析超高心墻壩填筑期孔壓發(fā)展規(guī)律及其影響提供可靠的支撐手段。

2 兩河口心墻壩孔壓監(jiān)測規(guī)律

兩河口心墻壩位于四川省甘孜州雅江縣，最大壩高295 m，壩頂寬16 m，河床部位心墻底高程2580 m，壩頂高程2875 m，正常蓄水位2865 m，上下游壩坡坡度分別為1∶2.0和1∶1.9(考慮馬道的綜合壩坡)。大壩施工過程中以左岸壩頂為原點，沿壩軸線共布設5個監(jiān)測橫斷面。根據(jù)已有現(xiàn)場監(jiān)測資料[3]，選擇0+340 m橫斷面為典型斷面，心墻中滲壓計(PDB)和電磁式沉降環(huán)(DC4)埋設位置如圖1所示。截至2021年8月，大壩填筑至高程2860 m，此時心墻最大沉降為2985 mm，位于0.5H(H為壩高)，如圖2所示。

圖1 典型斷面心墻滲壓計和電磁式沉降環(huán)分布

圖2 電磁式沉降環(huán)測值沿高程分布

圖3為心墻底部滲壓計實測孔壓時程以及填筑和水位的變化過程?？梢钥吹剑?1)不同位置的滲壓計孔壓發(fā)展是同步的，并且與填筑高程變化呈正相關關系；(2)2021年5月至8月水位增長108 m，實測孔壓僅增長0.1 MPa左右；(3)實測最大孔壓為3.74 MPa，明顯高于同類工程(見研究背景部分)。說明填筑期上游庫水還未來得及滲入心墻，高孔壓主要由填筑重力引起的，因此本文分析中并未考慮上游水頭入滲的影響。

圖3 心墻滲壓計實測孔壓時程

圖4為心墻滲壓計測值沿高程的分布。如圖4(a)所示，2019年8月不同橫斷面同一高程孔壓的相對差值達52%；同一橫斷面上下游對稱位置的滲壓計測值相對差值達91%。如圖4(b)(c)，孔壓隨深度呈非線性遞增規(guī)律；不同高程孔壓最大值所在位置沿高程是往復變化的，部分高程位于心墻邊緣處。說明兩河口大壩心墻孔壓呈顯著不均勻分布特征。

圖4 心墻滲壓計孔壓測值

3 基于廣義塑性本構模型的不排水流固耦合分析法及驗證

3.1 方法介紹為了反映兩河口大壩心墻孔壓與飽和度的耦合效應及其三維分布特征，首先將飽和度作為材料參數(shù)引入到廣義塑性本構模型中，建立了不排水條件下簡化、高效的流固耦合分析方法。具體實現(xiàn)過程如下。

(1)心墻有效應力-應變關系模擬。土石壩中筑壩料表現(xiàn)出明顯的彈塑性變形特性，為了更合理地描述筑壩料的力學特性，鄒德高等[17]考慮筑壩料的壓力相關性和循環(huán)滯回特性，對原始廣義塑性本構模型進行了改進，并成功應用于多座土石壩的靜、動力反應分析[18-19]。因此，本文采用改進后廣義塑性本構模型模擬心墻有效應力與應變的關系，兩者關系表達為：

dσ′=Dep：dε

(1)

dε=dεe+dεp

(2)

(3)

式中：dεe為彈性應變增量；dεp為塑性應變增量；Dep為彈塑性矩陣；De為彈性矩陣；ng為塑性流動方向；n為加載方向矢量；H為塑性模量；下標“L”和“U”分別表示加載和卸載。該模型的詳細介紹本文不再贅述，模型參數(shù)及其物理意義詳細簡介見文獻[20]。

(2)心墻孔壓-應變關系模擬。土石壩心墻土填筑飽和度通常在80%～93%[13]，氣體以封閉氣泡的形式存在于水中，孔隙流體的壓縮性由氣體和水共同決定。因此，將孔隙水氣作為一種可壓縮混合流體進行模擬時，孔壓和應變關系可描述為[8]：

dpw=Dw：dε

(4)

(5)

(6)

式中：Dw為孔隙流體的剛度矩陣；dpw為孔壓增量；dε為應變增量；Kaw為孔隙流體的體積壓縮模量；I為元素都為1的3×3矩陣；O為3×3的零矩陣。

已有研究[21-22]表明，對于土體中液相相連通、氣相不連通的氣封閉狀態(tài)(飽和度大于80%左右)，可忽略基質吸力的影響，則pw=pair，pair為孔隙氣體壓力。在該狀態(tài)下，土體中孔隙水氣的體積壓縮模量可以表示為[23]：

(7)

式中：n為孔隙率；Kf為純水的壓縮模量，Pa；pw為孔壓，Pa；pa為1個大氣壓，約等于1×105Pa；S為當前飽和度；h為空氣在水中的體積溶解系數(shù)，一般取值0.02。

隨著pair(等于pw)的增加，水中的氣泡會進一步溶解到水中，土體飽和度會逐級增大。依據(jù)Boyle定律和Henry定律，孔隙氣體壓力變化條件下當前飽和度S與初始飽和度S0關系可以表示為[24]：

(8)

(3)流固耦合方法的實現(xiàn)。土石壩心墻的滲透系數(shù)極小[13，25]，在填筑施工過程中流體滲流速度十分緩慢，此時土與孔隙流體呈同步變形的狀態(tài)[8]。根據(jù)有效應力原理，將式(1)和(4)相加，可得心墻總應力增量與應變增量關系：

dσ=(Dep+Dw)：dε

(9)

從上式可以看到，只需在已有的本構模型中增加Dw并記錄高斯點孔壓即可實現(xiàn)心墻孔壓的分析，并且當孔隙流體模量取0時，可自動退化為式(1)。由于有限元分析時外力是與總應力平衡的，因此上式很容易在單相介質中實現(xiàn)，并且無需修改應力積分和迭代求解流程。更為重要的是，因本文流固耦合方法本質是單相問題求解，因此避免了Biot流固耦合分析矩陣病態(tài)、求解速度慢的問題。

3.2 方法驗證基于兩河口心墻壩二維有限元模型，采用Biot流固耦合法對上述不排水流固耦合法的適用性進行了驗證。兩種方法中材料參數(shù)均相同，壩體填筑過程與實際工程保持一致(共歷時6 a)，材料分區(qū)、填筑模擬等詳細資料可見第4部分三維分析介紹。不排水流固耦合法中僅心墻采用不排水方法模擬，Biot流固耦合法中心墻兩側設定為排水邊界。已有研究[1，26-28]表明，試驗得到滲透系數(shù)會明顯高估現(xiàn)場孔壓的消散。為了與實際孔壓累積和消散趨勢保持一致，開展了不同心墻滲透系數(shù)的計算分析，結果表明，當滲透系數(shù)取1×10-10m/s時，計算的孔壓消散發(fā)生的時間與實測時間是基本一致的(見圖5滲透系數(shù)1×10-9m/s和1×10-10m/s的Biot流固耦合計算結果對比)，因此可認為心墻的滲透系數(shù)是1×10-10m/s。其值明顯小于試驗值，這與現(xiàn)場心墻非飽和、封閉氣泡對孔隙水以及反濾料對細顆粒的阻滯作用有助于提高心墻的抗?jié)B能力等因素有關。

圖5也給出了不排水流固耦合法與Biot流固耦合法(心墻滲透系數(shù)為1×10-10m/s)計算結果的對比。結果表明，Biot流固耦合法心墻計算沉降略大于不排水流固耦合法計算沉降，最大差值為9%；孔壓則略小于不排水流固耦合法計算結果，最大差值為7.5%左右。綜上，兩種分析方法的計算結果差別較小，不排水流固耦合法可用于兩河口心墻壩填筑期孔壓和變形模擬。

圖5 心墻沉降和孔壓時程

4 大壩三維靜力有限元分析和討論

4.1 有限元模型與材料參數(shù)

4.1.1 有限元模型 本次研究采用如圖6所示的有限元模型。為模擬兩河口心墻壩施工過程及水位變化，大壩分為49層填筑，其中圍堰9層，壩體40層。采用實際建設過程的蓄水方案[3]，邊填筑邊蓄水，庫水位分47步蓄至壩頂以下10 m。有限元模型采用了六面體等參單元和比例邊界多面體單元模擬，單元總數(shù)為57 754個，結點總數(shù)為61 230個，自由度總數(shù)為183 690個。計算采用大連理工大學抗震研究所開發(fā)的巖土工程非線性有限元分析程序GEODYNA[29]。

圖6 大壩有限元網(wǎng)格

4.1.2 材料參數(shù) 大壩計算模型包括圍堰、堆石料、過渡料和心墻料等分區(qū)(見圖6)。其中，心墻料參數(shù)根據(jù)室內試驗[30]確定，如表1所示。

表1 廣義塑性本構模型參數(shù)

4.2 計算工況兩河口心墻壩施工過程中對心墻料含水量進行了監(jiān)測，本文依據(jù)現(xiàn)場檢測提供的160個數(shù)據(jù)，對心墻初始飽和度進行了分析。如圖7所示，心墻實測初始飽和度S0基本分布在80%～100%，滿足均值88%，方差0.001的正態(tài)分布N(88%，0.001)。

圖7 心墻實測初始飽和度分布

為全面分析心墻飽和度對孔壓和變形的影響，開展了如表2所示工況的比較研究。其中工況1：S0為100%；工況2：S0為88%；工況3：S0按N(88%，0.001)隨機取值?？紤]到監(jiān)測條件有限，無法監(jiān)測每個位置的飽和度，為了充分考慮初始飽和度的隨機性，工況3進行了7組數(shù)值模擬。

表2 計算工況

4.3 結果分析

4.3.1 工況1和工況2 圖8對比了工況1和工況2中滲壓計PDB-60位置計算與實測孔壓時程，兩種工況中計算孔壓與實測值均差別較大，最大差值分別為43%和-46%。圖9對比了兩個工況下心墻計算和實測孔壓分布，可以得到：工況1較實測偏大，工況2較實測偏小，并且兩種工況計算孔壓離散性都較小，而實測值離散性較大。說明假定初始飽和度為常數(shù)難以反映實際孔壓量值大小和孔壓分布的顯著不均勻性。

圖8 心墻孔壓時程(滲壓計PDB-60)

圖9 心墻孔壓分布

圖10給出了心墻沉降沿高程的分布圖。工況1心墻最大沉降為3235 mm，位于0.63H，與實測最大沉降量值(2985 mm)和位置(0.5H)均存在較大的差別；并且0.4H以下計算沉降較實測值偏小，0.4H以上較實測值偏大。工況2心墻最大沉降為2729 mm，位于0.53H，與實測吻合較好；并且計算沉降與實測結果沿壩高分布規(guī)律是基本一致的，說明反映現(xiàn)場平均飽和度影響的變形計算結果與實際更吻合。同時對比3.2節(jié)二維計算沉降可知，不考慮大壩的三維空間效應會顯著高估兩河口大壩的計算變形。

圖10 壩軸線沉降沿高程的分布

4.3.2 工況3 (1)孔壓時程和分布。圖11給出了滲壓計PDB-60位置實測和工況3中7組計算孔壓時程的對比。可以得到：計算孔壓時程與初始飽和度密切相關，總體上初始飽和度越大計算孔壓越大，但兩者關系呈非線性變化關系；當S0在92%～93%時，計算孔壓時程與實測基本吻合，因此PDB-60位置的初始飽和度約為92%～93%。

圖11 工況3中7組數(shù)值模擬的孔壓時程(滲壓計PDB-60)

圖12對比了心墻實測孔壓與工況3計算結果沿壩高的分布，限于篇幅僅給出了2組計算結果。如圖12所示，3和7號數(shù)值模擬計算孔壓分布與監(jiān)測值基本一致(圖12(a)(b))，7組數(shù)值模擬孔壓基本能夠覆蓋實測孔壓分布區(qū)域(圖12(c))，說明考慮心墻飽和度隨機性后，可數(shù)值再現(xiàn)心墻孔壓不均勻分布規(guī)律。此外，不同數(shù)值模擬中孔壓最大值在3.5～4.7 MPa，均位于心墻底部，但計算最大值大于監(jiān)測最大值(3.74 MPa)，說明增加滲壓計數(shù)量可能會監(jiān)測到更大的孔壓值。但高孔壓也間接說明了心墻防滲性能良好，如果心墻滲透系數(shù)很大(見3.2節(jié))，孔壓會迅速滲透轉移，不會出現(xiàn)局部高孔壓的現(xiàn)象。

圖12 心墻孔壓分布散點圖

綜上，結合工況1和工況2計算結果可知，兩河口大壩的心墻孔壓的不均勻分布特征是心墻飽和度隨機性導致的。

(2)沉降變形分布。圖13給出了工況3心墻沉降沿高程的分布。7組數(shù)值模擬沉降完全重合，并與工況2的結果基本一致，最大差值僅為3%左右。說明保持平均飽和度不變，考慮飽和度隨機性分布對心墻沉降變形影響不大，這是因為大壩沉降是土體變形的累積效果，與土體平均飽和度相關，受材料的局部特性影響較小。

圖13 工況3中壩軸線沉降沿高程分布(未竣工)

(3)主應力分布。圖14為工況3心墻上游側單元的有效小主應力沿高程的分布，為簡潔，圖中僅給出了1號、3號和7號數(shù)值模擬的結果?？梢钥吹剑翰煌瑪?shù)值模擬中同一高程心墻有效小主應力具有明顯的差別，但有效小主應力均大于0，未出現(xiàn)拉應力。

圖14 工況3心墻上游側有效小主應力(未竣工)

(4)應力水平分布。土體的受力破壞程度主要與其應力水平相關(或與破壞面的距離有關)。圖15為工況3心墻中部一列單元的應力水平(應力水平=當前應力比/破壞應力比，等于1時表示破壞，等于0時表示等壓狀態(tài))?？梢钥吹剑嚎紤]飽和度的隨機性對心墻應力水平分布也具有較大影響，不同數(shù)值模擬中同一高程心墻應力水平具有明顯的差別，但心墻應力水平均低于0.85；此外，雖然心墻底部中間孔壓最大，但該位置應力水平低于0.45。說明兩河口大壩心墻高孔壓狀態(tài)下的力學性態(tài)是良好的。

圖15 工況3心墻應力水平(未竣工)

5 結論

(1)將飽和度作為初始材料參數(shù)引入廣義塑性本構模型，并考慮了飽和度隨著圍壓變化的過程，建立了不排水條件下簡化、高效的三維流固耦合分析方法，為分析高土石壩施工期心墻孔壓分布規(guī)律及其影響提供了有效的工具。

(2)心墻飽和度初始值顯著影響其孔壓，飽和度越大孔壓越大。當心墻初始飽和度以實測正態(tài)分布隨機取值時，數(shù)值分析再現(xiàn)了兩河口大壩心墻實測的高孔壓以及孔壓明顯不均勻分布的特征，闡明了孔壓分布規(guī)律與飽和度隨機性是直接相關的。

(3)飽和度隨機性對心墻沉降變形的影響低于3%，對心墻的局部應力和應力水平影響較為明顯，但心墻有效小主應力均大于0，應力水平均低于0.85，且底部中間高孔壓處的應力水平低于0.45。因此，飽和度隨機性導致的高孔壓以及孔壓不均勻分布不會影響兩河口大壩心墻的安全。