基于COMSOL Multiphysics的某水庫均質(zhì)壩滲流與穩(wěn)定性模擬分析研究

李宏偉

(烏魯木齊市水利勘測設(shè)計院(有限責(zé)任公司)，烏魯木齊 830001)

1 概 述

我國水資源分布面積較廣[1-3]，如何有效利用水資源是水利工作者夜以繼日思考的問題。為此，水利工程技術(shù)人員考慮修建水庫大壩等水工建筑物，用以水資源調(diào)度[4-5]。這些水工建筑物給人類的生存發(fā)展提供重要推動力，但其安全性一直是業(yè)內(nèi)研究考慮的重要方面。已有設(shè)計人員與學(xué)者基于室內(nèi)水工模型試驗(yàn)，設(shè)計安全穩(wěn)定性較高的水工建筑，提出水庫大壩類水工結(jié)構(gòu)安全設(shè)計施工的要點(diǎn)[6-9]。還有一些學(xué)者借助數(shù)值有限元軟件，模擬計算水工建筑物在長期運(yùn)營荷載下安全穩(wěn)定性，準(zhǔn)確預(yù)判壩體失穩(wěn)[10-12]。不論采用何種手段，其目的均是為了保證水工建筑物的安全狀態(tài)。針對具體的工程案例，采用數(shù)值模擬軟件，綜合分析獲得滲流與結(jié)構(gòu)穩(wěn)定性，為工程施工提供重要的理論參考。

2 滲流與穩(wěn)定性分析理論

2.1 滲流理論

滲流基本理論方程為達(dá)西公式，其表達(dá)式如下：

(1)

式中流速v與水力坡降J成正比關(guān)系。

基于達(dá)西基本理論，聯(lián)系巖土材料各向異性特點(diǎn)，獲得各方向上達(dá)西滲流微分表達(dá)式：

(2)

針對不可壓縮流體，獲得容納黏滯參數(shù)的應(yīng)力方程為：

(3)

而在多孔介質(zhì)中滲流運(yùn)動與孔隙度有關(guān)，運(yùn)動方程可表述為：

(4)

圖1 單元體中滲流示意圖

基于質(zhì)量守恒定律，單元體水質(zhì)量變化曲線為：

(5)

式中：n為介質(zhì)孔隙率；ρ為流體密度；V為單元體體積。

聯(lián)系滲流量公式，并考慮流體不可壓縮性，可得到：

(6)

式(6)即為剛度無限大的介質(zhì)中流體運(yùn)動微分方程，任意時刻任意點(diǎn)的流量值為零。基于二維流體方程，推廣至三維維度內(nèi)，穩(wěn)定滲流微分方程為：

(7)

當(dāng)介質(zhì)為各向同性時，微分方程基于Laplace變換，可得：

(8)

若流體運(yùn)動處于非穩(wěn)定狀態(tài)，則滲流運(yùn)動勢必需考慮介質(zhì)狀態(tài)與含水層，運(yùn)動微分方程可表述為：

(9)

當(dāng)介質(zhì)材料為各項同性時，式(9)可簡化得到：

(10)

不論是穩(wěn)定滲流或是非穩(wěn)定，運(yùn)動微分方程的求解主要基于初始條件與邊界條件，獲得符合條件的收斂值，有限元計算收斂值即是依靠計算機(jī)強(qiáng)大計算能力，獲得精度解。在前述分析基礎(chǔ)上，滲流運(yùn)動方程基本表達(dá)式為：

(11)

式中：h=h(x，y，z)為水頭函數(shù)關(guān)系式；Ss為蓄水量。

基于變分與等參量變換，獲得單元滲流矩陣表達(dá)式：

(12)

聯(lián)系滲流基本微分方程與單元體水質(zhì)量本構(gòu)，等參變換有：

(13)

利用高斯積分與簡化的雅可比矩陣，積分解為：

(14)

式中：φ(ξ，η，ζ)為積分函數(shù)在高斯積分點(diǎn)處的最高精度函數(shù)值。

通過微分方程解與滲流邊界條件聯(lián)系，即可獲得滲流場中滲流速度、水頭壓力等特征參數(shù)值。

2.2 穩(wěn)定性理論

COMSOL Multiphysics有限元軟件計算平臺計算水庫壩坡穩(wěn)定性可采用強(qiáng)度折減法，以邊坡失穩(wěn)和有限元計算不收斂同時出現(xiàn)，作為解精度要求的依據(jù)。強(qiáng)度折減法本質(zhì)上是以M-C強(qiáng)度準(zhǔn)則為內(nèi)在依據(jù)，公式如下：

(15)

其中：

(16)

基于有限元軟件計算壩坡直至不收斂時，巖土體抗剪參數(shù)臨界值及安全穩(wěn)定系數(shù)作為分析重要基礎(chǔ)，為評價壩坡安全穩(wěn)定性提供參量。

3 工程概況

某水庫為當(dāng)?shù)匦钏{(diào)度重要樞紐，滿足當(dāng)?shù)毓まr(nóng)業(yè)基本發(fā)展需求，并在旱季提供水資源輸送渠道，灌溉農(nóng)田面積超過2.667×104hm2，水庫庫容為500×104m3，有一條苗河提供水源，該河流流域面積超過10 km2，河流淤泥沉積懸浮沙過多時，影響到下游水力發(fā)電后，該水庫可作為反調(diào)節(jié)水利工程，引水調(diào)節(jié)河流水質(zhì)。蓄水工程包括泄洪道、涵閘室輸水孔與消力池等水工建筑，水庫正常水位設(shè)計為920.3 m，洪水位為921.3 m，泄洪道等水利設(shè)施均為Ⅳ級建筑，調(diào)度水平均流量為1.2 m3/s。壩體主軸線設(shè)計修建為與河道正交，設(shè)計高程為922 m，壩頂寬度為2.5 m，上下游壩體坡度為0.5和0.4，壩身屬均質(zhì)土。

現(xiàn)場地質(zhì)踏勘表明，壩址場地內(nèi)地質(zhì)構(gòu)造屬平靜狀態(tài)，無較活躍斷層帶與斷裂帶，但可見部分出露地表巖石節(jié)理發(fā)育較豐富，夾有風(fēng)化碎屑巖層。依據(jù)鉆孔資料得知，場地內(nèi)覆蓋層包括人工填土與粉質(zhì)黏土等，人工填土以種植土為主，韌性較差，夾有第四系河流搬運(yùn)作用形成的礫石，粒徑為2～10 mm，磨圓度較高，由于經(jīng)過壓實(shí)作用，含水量較低，僅為15%。區(qū)域內(nèi)地形地貌以河漫灘沖積平原為特點(diǎn)，兩側(cè)岸坡依附在山體中，山體坡度為25°～41°，巖層以半風(fēng)化花崗巖為主，中粗粒結(jié)構(gòu)。基巖為砂礫巖，下伏另有二長花崗巖與之不整合接觸，標(biāo)準(zhǔn)承載力超過300 kPa。為保證壩體安全穩(wěn)定，水庫沿線岸坡均已完成清淤清障，提升壩基承載力。

4 均質(zhì)壩滲流與穩(wěn)定模擬計算

4.1 模型建立與研究工況

均質(zhì)壩體斷面圖見圖2(a)。本文將基于斷面形態(tài)與巖土材料建立幾何模型，后以sat格式形式導(dǎo)入至COMSOL計算平臺中，幾何模型見圖2(b)。

圖2 斷面圖與幾何模型

以單元網(wǎng)格劃分幾何模型，獲得三角形單元網(wǎng)格為主的數(shù)值模型，共劃分出單元網(wǎng)格2 512個，節(jié)點(diǎn)數(shù)2 165個，單元質(zhì)量最小為0.45，見圖3(a)。圖3(b)為模型中單元質(zhì)量直方圖，直方值較高區(qū)域集中于后半部分，即壩體滲流與穩(wěn)定性可在求解之時獲得較高精度解。

圖3 數(shù)值模型與直方圖

所有邊界荷載參數(shù)均參照實(shí)際工程勘察資料報告，邊界約束荷載與工況有關(guān)，滲流場以正常蓄水位、水位回落期開展；穩(wěn)定性研究工況與壩體修建以及水庫運(yùn)營有關(guān)，故以穩(wěn)定滲流期、水位回落期及壩體施工期開展計算分析。

4.2 滲流分析

水位回落曲線與時間關(guān)系曲線見圖4。在COMSOL中將定義荷載與時間相關(guān)函數(shù)，分析水位回落曲線特征可知，水位陡降位于第24 d，速率為8 m/d。筆者認(rèn)為水位回落后半階段屬驟降工況，應(yīng)考慮對上游壩坡滲流與穩(wěn)定性影響。

圖4 水位回落曲線與時間關(guān)系曲線

基于上述分析分別獲得正常蓄水位與水位回落期兩工況下滲流場特征參數(shù)解，見圖5、圖6。從圖5、圖6中可看出，正常蓄水位穩(wěn)定滲流狀態(tài)下浸潤面并未蔓延至壩頂區(qū)域，迎水側(cè)壩身大部分區(qū)域滲透壓力均為正值；從滲透水壓力分布云圖亦可看出，下游側(cè)背水側(cè)及壩踵等區(qū)域滲透壓力為負(fù)值以下，最小為-0.2 MPa，最大滲透水壓力為1.4 MPa，位于上游迎水側(cè)壩基上覆蓋層等區(qū)域，壩坡與下游側(cè)滲透水壓力差符合滲流特征，初步判斷上游壩坡安全。結(jié)合壩體滲透水壓力矢量分布圖可知，滲透方向在壩體內(nèi)部由上游側(cè)指向下游側(cè)壩基，另基于COMSOL后處理分析獲得壩體滲流水頭等勢線，可看出上游側(cè)逐漸遞減至下游側(cè)，壩頂中心線等勢值為50%，水頭遞減規(guī)律顯著。正常蓄水位壩體最大水力坡降值為2.89，位于底部溝槽處，浸潤線上部起始點(diǎn)處水力坡降值為0.08。綜合正常蓄水位各特征參數(shù)值可認(rèn)為該均值壩體在正常蓄水位下滲流處于安全狀態(tài)，壩體不會出現(xiàn)滲透破壞等現(xiàn)象。

圖5 浸潤線(正常蓄水位)

圖6 滲流特征參數(shù)分布圖(正常蓄水位)

為研究水位回落期滲流場特征，給出壩體各個時間段內(nèi)水力坡降在壩坡高程上分布，見圖7，圖7中每一條曲線即為一個時間段。分析水力坡降在水位回落過程中表現(xiàn)可知，非穩(wěn)定滲流場中水力坡降值隨水位回落，逐漸增大，在第12 d～18 d每下降1 m高程，水力坡降增大0.009；當(dāng)?shù)?8 d～25 d，每下降1 m高程，水力坡降增大0.034，速率增大近1個量級，即水力坡降增大速率與降落速度有關(guān)。

圖7 各時間段水力坡降在壩體上分布曲線

在t=27 d后，水力坡降值超過0.53，即越過安全臨界水力坡降值，壩坡滲流安全性受到挑戰(zhàn)，為此利用COMSOL提取出t=27 d后滲流特征參數(shù)分布，見圖8、圖9。

圖8 浸潤線(水位回落期)

圖9 滲流特征參數(shù)分布圖(水位回落期)

從圖8、圖9中可分析看出，浸潤線最高點(diǎn)高程為917.5 m，逼近正常蓄水位高程；底部溝槽處滲透水壓力接近于零，最大滲透水壓力位于上游側(cè)壩基處，但壩體內(nèi)部滲透水壓力分布較大，下游側(cè)壩基處滲透水壓力由增大趨勢。由滲流矢量方向可看出，滲流方向指向分為兩類：上游側(cè)壩坡與壩體內(nèi)部，即上游側(cè)壩坡滲透壓力會在一定程度促進(jìn)壩坡滲流運(yùn)動，導(dǎo)致壩坡滲漏量增加或滲透破壞。最大水力坡降值為5.61，處于上游側(cè)壩坡內(nèi)，相比正常蓄水位下增大94.1%，水位回落過程中非穩(wěn)定滲流場各特征參數(shù)種種表現(xiàn)表明，該工況下壩坡內(nèi)滲流活動處于危險區(qū)間，應(yīng)考慮針對上游壩坡增加防滲措施結(jié)構(gòu)，在壩身表面鋪設(shè)土工防滲膜，合理調(diào)控水庫水位回落，保證壩體滲流安全。

4.3 穩(wěn)定性分析

穩(wěn)定性分析按照前述理論基礎(chǔ)，借助COMSOL有限元軟件開展正常蓄水位、水位回落(正常蓄水位回落至死水期)、施工期3個工況下壩坡穩(wěn)定性分析，獲得3種工況下壩坡失穩(wěn)滑弧云圖，見圖10。其中，水位回落工況共采用有效應(yīng)力強(qiáng)度折減法與最小應(yīng)力組合兩種計算方案。正常蓄水期與施工期壩坡模擬計算出的滑弧均出現(xiàn)在下游側(cè)，而水位回落期滑弧出現(xiàn)在上游側(cè)，表明水位回落期壩坡失穩(wěn)首先出現(xiàn)在迎水側(cè)。經(jīng)計算獲得3種工況下安全穩(wěn)定系數(shù)，見表1，3種工況下安全系數(shù)分別為1.54、1.33和1.73。對比規(guī)范要求臨界最小安全系數(shù)，3種工況下均處于安全狀態(tài)。

圖10 壩坡失穩(wěn)滑弧云圖

表1 3種工況下安全穩(wěn)定系數(shù)

5 結(jié) 論

引入三維滲流分析理論與強(qiáng)度折減方法計算穩(wěn)定性，基于某水庫均質(zhì)壩體工程資料分析，借助COMSOL Multiphysics數(shù)值軟件建立分析模型，研究不同工況下壩體滲流與穩(wěn)定性，結(jié)論如下：

1) 獲得了正常蓄水位下壩坡的上下游側(cè)滲透壓力差分布、浸潤線、水頭等勢線均由上游側(cè)遞減至下游側(cè)特征，最大水力坡降值為2.89，位于底部溝槽處，滲流處于安全，壩體不會出現(xiàn)滲透破壞。

2) 水位回落期中水力坡降隨水位回落逐漸增大，且增大速率與降落水位有關(guān)；浸潤線最高點(diǎn)高程為917.5 m，滲流方向指向上游側(cè)壩坡與壩體內(nèi)部；最大水力坡降值為5.61，處于上游側(cè)壩坡內(nèi)，壩體滲流處于危險狀態(tài)，應(yīng)考慮布設(shè)防滲結(jié)構(gòu)。

3) 分析了正常蓄水位、水位回落期、施工期3個工況下壩坡危險滑弧，得到3種工況下安全系數(shù)分別為1.54、1.33和1.73，均高于臨界最小安全系數(shù)，壩體穩(wěn)定性處于安全狀態(tài)。