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

喬 磊，趙 璇

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

1 概 述

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

2 工程概況

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

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

3 數(shù)值模擬

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

3.1 理論背景

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

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

(1)

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

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

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

(2)

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

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

(3)

當(dāng)流體流過(guò)多孔介質(zhì)時(shí)，有3個(gè)力作用在固體基質(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中，熱量通過(guò)以下兩個(gè)過(guò)程在多孔介質(zhì)中傳遞：當(dāng)熱量由流體運(yùn)動(dòng)引起時(shí)，存在強(qiáng)制對(duì)流；當(dāng)流體運(yùn)動(dòng)由溫度變化造成的密度差引起時(shí)，存在自由對(duì)流。

Flac中，用于對(duì)流-擴(kuò)散傳熱的能量平衡方程如下：

(6)

式中：T為溫度；qT為熱通量；qw為特定流體的流量；qvT為體積熱源強(qiáng)度；ρ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(從模型邊緣到堤岸的距離)。垂直位移固定在模型底部的邊緣處，水平位移固定在模型的左右邊緣處。用于計(jì)算的材料參數(shù)見(jiàn)表1。

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

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

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

圖2 模擬中的2D模型

4 結(jié)果分析

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

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

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

圖3 氣溫和平均波浪在點(diǎn)A-D模擬的溫度

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

分析溫度可知，在點(diǎn)B和點(diǎn)C處先開(kāi)始溫度變化，因?yàn)辄c(diǎn)B和點(diǎn)C離水庫(kù)較近，傳感器在點(diǎn)A和點(diǎn)D的反應(yīng)不如點(diǎn)B和點(diǎn)C明顯。

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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