基于Lu 模型的土石壩溫度場對壩體滲流場的響應分析

唐友山

（江西交通職業技術學院,江西 南昌 330013）

1 研究背景

土石壩是多孔介質材料,壩體熱量通過地表氣溫以熱傳導的形式傳入壩體內部,在地下水滲流影響下,溫度場會發生改變,溫度梯度也會引起水體流動,同時,溫度的變化也會使得流體物理性質發生變化,壩體滲流場與溫度場是相互影響的耦合關系,可以通過研究壩體溫度場的變化情況來反饋壩體滲流場,從而達到監測大壩滲漏的目的[1]。關于兩場耦合的研究很多,多數研究只考慮單向耦合關系,如模型中未能考慮流體物性隨溫度的變化,非飽和區域的導熱系數的空間非均質性等[2]。土石壩溫度場受氣溫影響,在表層溫度向大壩內部傳遞時,需穿過非飽和區域才能傳遞到壩體內部,非飽和滲流會影響壩體溫度場,如土中含水率對傳熱參數的影響[3-5]。因此,不能簡單地將非飽和土體當作干燥土體或者完全飽和土體計算,需考慮土體含水率影響,本文通過引入Lu 模型來刻畫導熱系數與含水量之間的定量關系,模型中同時考慮水體的溫變效應,使得模型計算結果更加接近實際工程。

2 描述非飽和介質導熱系數的Lu 模型

Lu 模型考慮土體礦物成分及含水率影響,導熱系數在完全飽和土（λsat）和干燥土（λdry）間插值[6],其表達式為：

λdry和λsat分別為：

式中：λs為土中固體顆粒導熱系數,,q 為石英在固體顆粒中的含量,λq為石英導熱系數,λ0為其他礦物質綜合導熱系數,其中,λq=7.7 W·m-1·℃-1,λ0=2.0 W·m-1·℃-1（q＞20%）,λ0=3.0 W·m-1·℃-1（q＜20%）。

對于常規土,插值系數Ke為：

式中：χ為參數,對于粘土、壤土、砂土分別取0.58、0.9、1.05。

3 飽和非飽和滲流場與溫度場數學模型

3.1 溫度場影響下的飽和-非飽和滲流場控制方程

非飽和土滲流計算中,土體滲透系數的影響因素很多,含水率是重要因素,通常采用土水特征曲線來表征土壤含水率與滲透系數間的關系,研究中一般通過建立經驗模型來進行計算,在眾多模型中,VG 模型[7-8]具有參數物理意義明確、精度高、適應性廣的特點,因此,本文選擇VG 模型進行計算,其表達式為：

式中：h（θ）為土壤基質吸力；k（θ）為非飽和土體滲透系數；ks為飽和土滲透系數；θ為含水率；θr和θs分別為土壤殘余含水率和飽和含水率；α和n 為與土體類別有關的參數,。

VG 模型參數可通過實驗獲取,沒有實驗條件的情況下也可按照經驗取定。非飽和土水力數據庫UNSODA[9]中含有世界各地大量土壤水力參數,可以通過在該數據庫中查取水力參數,通過數據擬合得到VG 模型參數。缺少數據的情況下也可按照經驗取定。

考慮溫度對滲透系數和溫度梯度對滲流的影響,得到滲流區域Ω內任意一點p（x,y,z）的總水頭方程滿足下式：

式中：H 為總水頭；hp為土壤基質吸力；k 為土壤滲透系數；μ為給水度；Qh滲流源匯項；,為土體容水度；θ為土體含水率；n 為土體孔隙率；SS為單位貯水率；Γ為邊界；其他符號意義同上。

3.2 滲流場影響下溫度場控制方程

土石壩壩體熱量傳遞包括滲流傳熱和多孔介質熱傳導兩部分,則單位時間內流入土體單元的熱量為：

式中：λ為導熱系數；q 為熱流密度；cw為流體比熱容；ρw為流體密度；t 為時間；v 為流體流速。

流入單元體的凈熱量應等于單元體所吸收的熱量,若單元體內無熱源,則溫度區域Ω內任意一點p（x,y,z）的溫度函數T 應滿足下式：

式中：T0（x,y,z）初始時刻溫度場；S1為已知溫度邊界；T1（x,y,z）為已知溫度邊界；q（x,y,z,t）為已知熱流邊界（q=0 為絕熱邊界）；S2為已知熱流邊界；n2為S2法向；n3為S3法向；β為表面放熱系數；Ta為地面溫度；Qh問溫度場源匯項；其他符號同上。

4 考慮非飽和滲流的壩體溫度場特征

4.1 幾何模型

幾何模型為一均質土石壩斷面,壩體填筑材料為粉土；最大壩高16.2 m,壩頂寬5 m。上游坡比為1∶3,下游壩坡為1∶2.5。壩腳設堆石棱體排水,頂寬1.5 m,高3 m,排水體的材料為礫石。為了便于分析計算結果,模型中取了7 個觀察點（1#～7#）,各觀察點坐標見表2,計算模型見圖1。

圖1 均質土石壩幾何模型

表2 觀察點坐標

4.2 計算參數

壩體材料的熱力學和物理力學參數見表3。非飽和滲流參數參考文獻見表4[11]。

表3 試驗土的熱力學和物理力學參數

表4 非飽和滲流計算參數

4.3 邊界條件與初始條件

滲流場邊界設定為：上游邊界AB 與下游邊界DC 取為定水頭邊界條件,其他邊界取為零通量邊界條件。

溫度場邊界條件設定為：上游邊界AB 與下游邊界AD和DC 取為已知溫度邊界條件,壩底取為絕熱邊界。

4.4 計算模型

計算模型見圖1,模型中采用Lu 模型來描述非飽和滲流區域土體的導熱系數與飽和度的定量關系,由于comsol 軟件自帶的多孔介質熱傳導模塊默認的導熱系數為體積平均模型,需修改模塊,引用Lu 模型。首先需定義λdry、λsat和λs,修改軟件中熱導系數方程。同時,模型中引入水的動力粘度和水密度與溫度間的函數關系,待后續的模型參數定義時引用。

模型中其他參數的取值見表3 和表4。滲流場邊界條件為：上游邊界AB 水頭為13.2 m,下游邊界DC 水頭為1 m,其余邊界取為零通量邊界；溫度場邊界條件為：上游壩面A 處溫度設為18℃,上游壩面溫度沿水深呈線性變化,庫底溫度穩定為5℃。下游邊界AC 溫度為20℃,BC 邊界取為絕熱邊界,初始溫度取為10℃。

4.5 計算結果分析

（1）液體體積分數對巖土體熱參數的影響

圖2 為大壩壓力水頭分布圖,圖中,壓力水頭為0 m 的等水頭線即為壩體浸潤線,非飽和區域土體飽和度會影響壩體材料的導熱性。圖3 為壩體非飽和區域的飽和度分布圖。由圖3 可知,作用在土體上的基質吸力使得地下水向上運移,浸潤線附近及以下土體飽和度為 1.0,越靠近壩頂土體飽和度越低,壩體頂部附近的相對飽和度達到 0.5,棱體上邊界的相對飽和度最低,為0.12,這是因為棱體材料為礫石,土體孔隙通道大,基質吸力低,水分向上遷移困難。飽和度在非飽和區域內非線性變化。

圖2 壓力水頭分布圖

圖3 相對飽和度分布圖

非飽和帶土體的導熱系數與其含水量的關系是由 Lu 模型給出的,見圖4,導熱系數隨液體體積分數的增大而增大。圖5 為壩體非飽和土導熱系數分布圖。由圖5 可知,浸潤線以下的飽和土的導熱系數為1.81 W·m-1·K-1,壩體導熱系數最大值出現在下游水位以下的排水棱體處,為2.92 W·m-1·K-1,這是因為棱體材料為礫石,介質石英含量高,且水位以下巖體飽和,故導熱系數最大。最小值出現在排水棱體上邊界,為1.48 W·m-1·K-1,這是因為棱體上邊界含水量最低,介質孔隙空氣含量大,由于空氣的導熱系數遠遠低于水體導熱系數,阻礙了熱量在巖體中的傳遞,故導熱系數最小。最小值為最大值的50.7%。

圖4 非飽和土導熱系數與液體體積分數關系圖

圖5 非飽和土導熱系數分布圖

（2）不同熱參數計算方法的對比分析

在進行多孔介質滲流傳熱的數值模擬研究中,非飽和巖土體的導熱系數取值,常用以下兩種方法：（1）常量法,導熱系數取值按飽和巖土體或某一含水量狀態給出；（2）變量法,導熱系數按固、液、氣三相的體積分數確定。本文所采用的Lu 模型采用變量法,是一個半經驗半理論模型,通過在comsol 中修改等效熱傳導系數方程來進行計算。在前面計算的基礎上,通過引用不同的熱參數計算方法來說明Lu 模型的有效性,工況假設見表5。表中為水體導熱系數,取為0.58 W·m-1·K-1,為空氣導熱系數,取為0.024 W·m-1·K-1,按照Lu 模型給出的公式計算,其余參數同上。

表5 不同熱參數計算方法的對比分析工況

運用上述5 種方法分別計算得到的導熱系數垂向分布見圖6。由圖6 可知,Lu 模型、指數模型和體積平均模型在壩體垂向上均受含水率的影響,導熱系數隨高程的增大而減小,越靠近壩頂,導熱系數越低,當高程低于11 m 時,壩體導熱系數基本穩定不變,這是因為越接近壩頂,壩體填土含水率越低,故導熱系數也越低,當高程低于11 m 時,土體飽和導熱系數穩定不變,由圖4 也可以看出,土體導熱系數隨含水率的增大而增大。

圖6 不同熱參數計算方法的非飽和土導熱系數垂向分布

5 種方法分別計算得到的導熱系數中,采用體積平均模型計算得到的導熱系數最大,為2.35 W·m-1·K-1～2.46 W·m-1·K-1,超過了飽和土體的導熱系數,同時變幅也最小,說明采用體積平均模型計算得到的土體導熱系數與實際偏差較大；采用Lu 模型和指數模型計算所得到的導熱系數分布形態相近,計算值很接近,均介于干燥土體導熱系數值與飽和土體導熱系數值之間,說明兩種方法的計算值是合理的,且兩種方法計算所得的導熱系數受含水率影響很大,能夠反應含水率與熱傳導系數的非線性關系。因此,采用Lu 模型可以很好地刻畫導熱系數與飽和度之間的定量關系,評價不同飽和度巖土體的導熱系數,可為非飽和巖土體傳熱模型提供理論基礎。

5 結論

（1）通過建立土石壩飽和-非飽和滲流場與溫度場全耦合模型,可以為開展兩場耦合數值模擬提供理論基礎,從而建立更加精細的數值模型。

（2）Lu 模型能夠精確的描述土體導熱系數與含水率間的定量關系,相比于常用的體積平均模型,Lu 導熱系數模型計算結果更加合理。