基于力學(xué)-熱學(xué)-水文地質(zhì)耦合的堤防穩(wěn)定性數(shù)值模擬

喬 磊，趙 璇

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

1 概 述

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

2 工程概況

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

圖1 堤壩地質(zhì)模型

3 數(shù)值模擬

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

3.1 理論背景

用達(dá)西公式表示地下水-力學(xué)耦合流體的運輸：

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

(1)

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

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

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

(2)

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

滲透率K與水力傳導(dǎo)率kH的關(guān)系為：

(3)

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

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

qT=-kT?T

(4)

其中，有效導(dǎo)熱系數(shù)kT是由兩部分構(gòu)成，分別為流體kwT導(dǎo)熱系數(shù)和固體ksT導(dǎo)熱系數(shù)，可由下式表示：

kT=ksT+nSkwT

(5)

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

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

(6)

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

cT可定義為如下形式：

cT=ρdCv+nSρ0cw

(7)

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

3.2 數(shù)值模擬

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

表1 建模中使用的力學(xué)、熱學(xué)和水文地質(zhì)特性

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

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

圖2 模擬中的2D模型

4 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 數(shù)值模擬結(jié)束270 h后的位移、孔隙壓力和溫度

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

表3 不同波形的安全系數(shù)

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

5 結(jié) 論

根據(jù)本文研究內(nèi)容，結(jié)論可分為以下兩個方面：

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

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