基于力學-熱學-水文地質耦合的堤防穩定性數值模擬

喬 磊，趙 璇

(德州市水利局 丁東水庫運行維護中心，山東 德州 253000)

1 概 述

山東省德州市的丁東水庫沿線布設了防洪堤，防洪堤是保護地區免受洪水侵襲的一種常用方法[1]。這種堤防是比較通用的一種巖土結構，但有的堤防較為薄弱，對人民的生活和財產而言存在一定危險，所以監測堤防穩定性也成為當下重點關注的項目[2-3]。如今，研究人員正在尋找一種通用有效的方法來監測堤防穩定性，并預測其破壞的時間和地點[4-5]。許多項目嘗試使用熱測量，即用光纖特性來測量堤防溫度或其他儀器來測量、估計堤防中流體的流量。ISMOP項目是基于溫度和孔隙壓力傳感器建立的一個復雜又危險的預測系統。在實際測量之前，先對水庫內水位變化過程中發生的反應進行數值模擬，并得出相應結果，同時為了檢驗洪水波對堤防穩定性的影響程度需要耦合數值模型。

2 工程概況

丁東水庫位于山東省德州市陵城區丁莊鄉境內，西北距德州市25.0 km，是一座中型圍壩引水式平原水庫，屬于黃河中下游沖擊平原河間洼地。巖性中上部為沙壤土夾裂隙黏土，下部為極細砂、細砂層。壩基地層巖性為第四系全新統沖積層，壩基地層巖性為第四系全新統沖擊堆積層(a1Q4)，夾湖積堆薄層(1Q4)，自上而下分為五大層，10個亞層。堤頂公路為瀝青混凝土路，圍壩部分全長11 636 m，凈寬6 m，隔壩部分全長1 400 m，凈寬8.0 m，面層鋪設厚度為3 cm。現根據大壩由西至東方向的橫截面構建模型，模型見圖1。

圖1 堤壩地質模型

3 數值模擬

本文采用Itasca Flac 2D 7.0軟件進行數值模擬。該軟件使用有限差分法求解巖土工程問題。

3.1 理論背景

用達西公式表示地下水-力學耦合流體的運輸：

qw=-K?(P-ρwg·x)

(1)

式中：K為流體流動系數(或滲透率)；ρw為流體密度；g為重力。

該方程中的流體密度與溫度變化相關，可表示為如下形式:

ρw=ρ0[1-βf(T-T0)]

(2)

式中：T0為參考溫度；βf為流體體積的熱膨脹。

滲透率K與水力傳導率kH的關系為：

(3)

當流體流過多孔介質時，有3個力作用在固體基質上，即固體重量、浮力和阻力或滲透力。模擬中都考慮了這些力。

用傅里葉公式描述Flac中的熱傳遞：

qT=-kT?T

(4)

其中，有效導熱系數kT是由兩部分構成，分別為流體kwT導熱系數和固體ksT導熱系數，可由下式表示：

kT=ksT+nSkwT

(5)

在Flac中，熱量通過以下兩個過程在多孔介質中傳遞：當熱量由流體運動引起時，存在強制對流；當流體運動由溫度變化造成的密度差引起時，存在自由對流。

Flac中，用于對流-擴散傳熱的能量平衡方程如下：

(6)

式中：T為溫度；qT為熱通量；qw為特定流體的流量；qvT為體積熱源強度；ρ0為流體的參考密度；cw為流體的比熱容；cT為有效比熱容。

cT可定義為如下形式：

cT=ρdCv+nSρ0cw

(7)

式中：ρd和Cv分別為固體物質的體積密度和體積比熱容；n為孔隙度；S為飽和度。

3.2 數值模擬

數值模擬建立924×70的方形網格，尺寸為0.1 m。建立的模型深度比地形低2.5 m，為了避免邊界效應，偏移量為20 m(從模型邊緣到堤岸的距離)。垂直位移固定在模型底部的邊緣處，水平位移固定在模型的左右邊緣處。用于計算的材料參數見表1。

表1 建模中使用的力學、熱學和水文地質特性

再根據以下4個步驟進行建模：①僅計算地質(無堤防)的力；②計算地質和堤防的力；③1 h計算步驟：首先進行流體傳遞，然后傳送熱量，最后使水力機械平衡。每一步完成之后，水位和氣溫的值將隨之變化；④安全系數的計算。

圖2為根據上述內容建立的帶有測量點的2D模型。這些點的位置與壩內溫度和孔隙壓力傳感器的實際位置較接近。點A-D在高于地形0.4 m的地方，兩邊距離堤防交界點1 m左右。

圖2 模擬中的2D模型

4 結果分析

對不同的波浪參數進行17次數值模擬：一次模擬的是平均波浪參數，其余16次模擬只有一個參數發生變化，將各工況列于表2：水位增加對應的測試時間是從0到3.5 m；高水位3.5 m時對應的測試時間；水位下降對應的測試時間是從3.5到0.1 m；低水位0.1 m時對應的測試時間。

表2 數值模擬中波浪階段的時間 /h

安全系數(FoS)是在每個階段的波浪結束時進行計算，此參數用于評估洪水過程中堤防的穩定性。當初始土壤溫度為8℃、蓄水池內的水溫固定在11.04℃時，空氣溫度見圖3，水位變化見圖4。

圖3 氣溫和平均波浪在點A-D模擬的溫度

圖4 平均波浪參數模型的水位和位移長度

分析溫度可知，在點B和點C處先開始溫度變化，因為點B和點C離水庫較近，傳感器在點A和點D的反應不如點B和點C明顯。

由圖5可知，任意時間段，熱反應都滯后于孔隙壓力傳感器的反應，這是因為防洪堤在溫度較低的情況下，在入滲過程中發生了冷卻作用。在波浪的第一階段，1 h內垂直和水平方向的位移增加。達到高水位(3.5 m)后，點B和點C的位移開始減少，點A和點D位移繼續增加，因為點A和點D距水庫較點B和點C更遠。當水位下降時，所有點的增加速度都下降。但在低水位(0.1 m)時仍有一些位移，主要是由于堤內水流運動形成的。

圖5 基于平均波浪和水位的點A-D的孔隙壓力

圖6至圖9顯示了1 h內每階段的位移、孔隙壓力和溫度都具有平均波浪。由圖6-圖9中可以清楚地看到孔隙壓力異常，即孔隙含水飽和度不為零，與溫度變化或位移變化的區域不對應。

圖6所示為第一階段，水位增加至3.5 m后。不對稱(右)堤防和水庫底部位移值最大，水庫底部附近的孔隙壓力在減少，與空氣或水接觸的區域溫度在升高。

圖6 水位上升結束48 h后的位移、孔隙壓力和溫度

圖7所示為水庫高水位階段之后。孔隙壓力等值線幾乎為橢圓形；最大水流在不透水層上；最大位移在水庫底部、也在右邊堤防的頂部；溫度變化小于水的滲透范圍。

圖7 高水位結束108 h后的位移、孔隙壓力和溫度

圖8所示的水位下降之后，降低了位移和孔隙壓力的值，只有溫度超過8℃以上的區域才會增大，這些變化是由空氣-土壤對流和水滲透而造成的。這個結論對于最后一個階段，即圖9所示的低水位而言也是正確的，位移接近于零，因此水越過防洪堤頂部的堤防層。

圖8 水位停止下降192 h后的位移、孔隙壓力和溫度

圖9 數值模擬結束270 h后的位移、孔隙壓力和溫度

在所有數值模擬中，安全系數都是在每個階段結束時計算。不同波形的安全系數見表3，不同波浪階段的時間見表2，FoS的值按每個階段最短到最長的時間順序算出。對于所有測試過的模型，堤壩是穩定的(FoS>1)；當水庫出水速度最快時，FoS值最低，這個階段使堤壩穩定性從3.5降低到3.0。

表3 不同波形的安全系數

第一階段，增加水位時的持續時間對試驗結束后的FoS值(變化范圍為3.181到3.193)幾乎沒有影響。高水位(從3.14到3.25)和低水位(從3.09到3.26)的最終FoS值差異最大，長時間的低水位使堤防穩定狀態更接近試驗前的狀態。

5 結 論

根據本文研究內容，結論可分為以下兩個方面：

數值模擬結果表明，可以同時計算熱力場和水力場的數值，可以用熱傳感器代替昂貴的孔隙壓力傳感器，但傳感器需放置在不受溫度變化影響的位置。此外，熱傳感器對入滲水的反應要慢于孔隙壓力傳感器，是因為入滲過程中土壤和水之間發生了能量轉換。

穩定性分析結果表明，水位下降，堤防的穩定性隨之降低，安全系數從3.5降低至3.0。對所有測試過的洪水波參數來說，堤防是非常穩定的，任一階段的FoS差異都不大，均小于5%。