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

(1.大連理工大學(xué) 海岸與近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.大連理工大學(xué) 水利工程學(xué)院，遼寧 大連 116024；3.中國(guó)電建集團(tuán)成都勘測(cè)設(shè)計(jì)研究院有限公司，四川 成都 610072)

1 研究背景

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

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

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

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

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

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

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

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

圖3 心墻滲壓計(jì)實(shí)測(cè)孔壓時(shí)程

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

圖4 心墻滲壓計(jì)孔壓測(cè)值

3 基于廣義塑性本構(gòu)模型的不排水流固耦合分析法及驗(yàn)證

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

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

dσ′=Dep：dε

(1)

dε=dεe+dεp

(2)

(3)

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

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

dpw=Dw：dε

(4)

(5)

(6)

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

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

(7)

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

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

(8)

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

dσ=(Dep+Dw)：dε

(9)

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

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

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

圖5 心墻沉降和孔壓時(shí)程

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

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

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

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

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

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

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

圖7 心墻實(shí)測(cè)初始飽和度分布

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

表2 計(jì)算工況

4.3 結(jié)果分析

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

圖8 心墻孔壓時(shí)程(滲壓計(jì)PDB-60)

圖9 心墻孔壓分布

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

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

4.3.2 工況3 (1)孔壓時(shí)程和分布。圖11給出了滲壓計(jì)PDB-60位置實(shí)測(cè)和工況3中7組計(jì)算孔壓時(shí)程的對(duì)比?？梢缘玫剑河?jì)算孔壓時(shí)程與初始飽和度密切相關(guān)，總體上初始飽和度越大計(jì)算孔壓越大，但兩者關(guān)系呈非線性變化關(guān)系；當(dāng)S0在92%～93%時(shí)，計(jì)算孔壓時(shí)程與實(shí)測(cè)基本吻合，因此PDB-60位置的初始飽和度約為92%～93%。

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

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

圖12 心墻孔壓分布散點(diǎn)圖

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

(2)沉降變形分布。圖13給出了工況3心墻沉降沿高程的分布。7組數(shù)值模擬沉降完全重合，并與工況2的結(jié)果基本一致，最大差值僅為3%左右。說(shuō)明保持平均飽和度不變，考慮飽和度隨機(jī)性分布對(duì)心墻沉降變形影響不大，這是因?yàn)榇髩纬两凳峭馏w變形的累積效果，與土體平均飽和度相關(guān)，受材料的局部特性影響較小。

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

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

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

(4)應(yīng)力水平分布。土體的受力破壞程度主要與其應(yīng)力水平相關(guān)(或與破壞面的距離有關(guān))。圖15為工況3心墻中部一列單元的應(yīng)力水平(應(yīng)力水平=當(dāng)前應(yīng)力比/破壞應(yīng)力比，等于1時(shí)表示破壞，等于0時(shí)表示等壓狀態(tài))。可以看到：考慮飽和度的隨機(jī)性對(duì)心墻應(yīng)力水平分布也具有較大影響，不同數(shù)值模擬中同一高程心墻應(yīng)力水平具有明顯的差別，但心墻應(yīng)力水平均低于0.85；此外，雖然心墻底部中間孔壓最大，但該位置應(yīng)力水平低于0.45。說(shuō)明兩河口大壩心墻高孔壓狀態(tài)下的力學(xué)性態(tài)是良好的。

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

5 結(jié)論

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

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

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