堤壩滲漏及結構破壞條件下滲流傳熱特性的數值分析

蔣中明,楊江寅,黃湘宜,廖峻慧,肖喆臻

(1.長沙理工大學水利與環境工程學院,湖南 長沙 410114; 2.長沙理工大學洞庭湖水環境治理與生態修復湖南省重點實驗室,湖南 長沙 410114)

滲透破壞既是威脅堤壩安全的主要問題,也是堤壩破壞的主要形式[1]。準確探測滲透破壞的產生及擴展,可為防汛搶險的決策提供強有力的科學依據。近年來,以溫度為對象的堤壩滲漏隱患探測技術成為研究熱點[2-4]。美國、德國等國家最早將溫度作為天然示蹤劑,通過捕捉溫度變化來研究堤壩滲漏情況[5]。該方法于20世紀80年代引入國內,后被稱為溫度示蹤滲漏監測技術[6]。早期溫度示蹤滲漏監測技術需在水工建筑物或其基礎內埋設大量熱敏溫度計,由于埋設點有限,其對溫度場的監測不夠準確[7]。隨后興起的分布式光纖溫度測量技術克服了前者的不足,但需在堤壩內埋設大量光纖,工作量依然巨大,且在已建工程中難以實現。目前,紅外熱成像法憑借其非接觸、無污染、速度快、可連續掃描等優點,成為堤壩滲漏探測技術研究的重點發展方向[8]。

深入分析滲漏條件下堤壩內部溫度場的變化規律,探明滲漏出口處溫度的演變特性是紅外熱成像堤壩滲漏探測技術的基礎。國內外研究人員針對滲漏條件下堤壩溫度場變化規律等方面開展了大量的研究,并取得了豐碩的研究成果。理論分析方面,羅日洪等[9]建立了穩定滲流影響下堤壩滲漏熱流耦合模型,分析了壩體中滲流場對溫度場的影響;蘇懷智等[8]認為當滲流存在時,水體遷移將迫使土體溫度與水溫相適應,從而引起溫度場的局部不規則變化。室內試驗方面,倪楓[10]通過物理模型試驗研究了不同入滲水頭、滲漏通道及入滲水溫等條件下堤壩溫度場與滲流場的變化規律;馬佳佳等[11]開展了多種工況下土壩滲漏的紅外熱像探測試驗;Song等[12]采用水熱耦合分析對堤壩滲流狀況進行監測,提出溫度可成為判斷堤壩前期滲流狀況的有效工具。數值模擬方面,陳建生等[13]提出了利用虛擬熱源法研究壩基裂隙巖體中存在的集中滲漏通道,通過溫度場來判定鉆孔周圍是否存在滲漏,并定量計算出滲漏量等參數;Yousefi等[14]以Shamil大壩為對象,研究了壩體內滲流對壓力及溫度場的影響,表明溫度比壓力更能反映滲流狀態;張啟義等[15]通過對雙層堤基異常滲漏險情下的溫度場進行計算分析,在滲漏通道出口處觀測到明顯的低溫區;Cheng等[16]基于滲流傳熱理論,采用數值模擬法對土石壩滲漏過程中的傳熱特性進行分析,提出孔隙內滲流穩定后,飽和區與非飽和區之間會形成溫差,進而可通過溫度判斷壩體內浸潤線位置。

堤壩滲漏是一個十分復雜的非穩態過程,盡管堤壩工程滲漏問題已從單一場的研究向多場耦合方向發展,但在數值模擬方面依舊存在考慮的環境條件和邊界荷載不夠全面等問題。為此,本文在現有研究基礎上,采用數值模擬方法,考慮堤壩工程實際邊界情況,研究堤壩滲漏乃至結構破壞過程中溫度場的演變特性,以期探明堤壩滲漏及結構破壞與溫度場演變的內在聯系,為基于紅外熱成像的堤壩滲漏隱患探測技術提供一條有益的實現途徑。

1 數值模擬理論

1.1 多孔介質滲流傳熱理論

多孔介質滲流傳熱過程中涉及復雜的耦合作用。滲流方面,堤內及地下水的流動符合達西定律。假設流體和固體不可壓縮,流體運動方程可參考文獻[17]。傳熱方面,堤壩滲流是一個十分緩慢的過程,可忽略土體與滲水之間的能量交換,基于局部熱平衡理論的對流-擴散熱傳導能量平衡方程可參考文獻[18]。

1.2 多孔介質熱流固耦合理論

多孔介質中土體、水及溫度之間相互影響?；诙嗫捉橘|線彈性假設,當產生非等溫滲流時,土體本構方程表述為[19]

σij=2Gεij+λεkk+αPδij-3βTbKb(T-T0)δij

(1)

式中:σ為應力;ε為應變;i、j、k代表坐標軸方向;G為剪切模量;λ為拉梅常數;α為比奧系數;P為孔隙水壓力;δ為Kronecher符號;βTb為多孔介質熱膨脹系數;Kb為多孔介質體積模量;T為溫度;T0為參考溫度。

1.3 土體結構損傷分析理論

沈珠江[20]基于連續介質力學基本概念,提出了考慮黏土結構破壞損傷過程的損傷模型,后被廣泛應用于工程界。該土體損傷模型通過定義損傷系數,對各項力學及流體計算參數進行折減,以反映土體的損傷程度。損傷系數和損傷過程中任意時刻土體參數的損傷值為[21]

ω=1-e-(aεv+bεs)

(2)

E=ωE0+(1-ω)Ed

(3)

K=ωK0+(1-ω)Kd

(4)

式中:ω為損傷系數;a、b為擬合參數;εv為單元體積應變;εs為單元剪切應變;E、K分別為損傷后的彈性模量和滲透系數;E0為單元初始時刻彈性模量;Ed為單元定義損傷彈性模量;K0為單元初始時刻滲透系數;Kd為單元定義損傷滲透系數。

利用內嵌FISH語言實現基于FLAC3D的土體損傷計算,具體步驟如下:

a.讀取單元初始參數,并保存在單元額外變量中,作為后續損傷計算的基準。

b.讀取各單元的體積應變及剪切應變,采用式(2)計算各個單元的損傷系數值,并保存于第二個單元額外變量中。

c.根據各個單元的損傷系數值,采用式(3)(4)分別計算損傷過程中各個單元的參數值,并保存于第三個單元額外變量中。

d.運用FISH的賦值功能,將第三個單元額外變量賦值給模型的對應參數,由此實現對單元參數的損傷。

e.通過FISH的循環調用功能,將上述步驟反復執行,以模擬計算過程中土體單元實時損傷的效果。

2 模型有效性驗證

2.1 飽和多孔介質流體熱對流計算

飽和多孔介質中,流體流動時會攜帶熱量而發生強制對流傳熱,流體靜止時會由于內部溫度差異導致密度不均勻,進而發生自由對流傳熱。Zhao等[22]提出了一種求解穩態飽和多孔介質流體自由對流問題的漸進逼近方法,解決了采用傳統有限元方法難以得到準確對流解的弊端。本文通過引入文獻[22]中飽和多孔介質中流體的熱對流計算,驗證FLAC3D軟件求解多孔介質中滲流傳熱問題的有效性及準確性。文獻[22]計算模型及監測點布置如圖1所示,模型高度為10km,寬度為10km,厚度為1m,共計網格節點169個,單元總數144個,監測點5個。

圖1 熱對流計算模型(單位:km)

滲流計算邊界:左右兩側及上下兩側均為不透水邊界。傳熱計算邊界:頂面為固定溫度邊界,大小為20℃;底面也是固定溫度邊界,其值為220℃;左右兩側均為絕熱邊界。具體計算參數如下:滲透率為1×10-8m/s,孔隙率為0.1,導熱率為0.335W/(m·K),比熱容為803.2J/(kg·K),熱膨脹系數為7×10-5K-1。圖2為初始溫度分布和初始孔隙壓力分布圖。

圖2 初始溫度和初始孔隙壓力分布

給模型施加順時針旋轉的流速,驅動流體發生強制對流。隨著強制對流現象的產生,多孔介質內的溫度發生改變。當瑞利數大于臨界值時,自由對流現象產生。最后,將初始設置的傾斜加速度去除,由流體內部溫度差驅動流體發生自由對流。飽和多孔介質中流體自由對流的溫度分布穩態解如圖3所示(圖3(b)為采用文獻[22]的漸進逼近方法計算得到的解析解)。

圖3 穩態溫度分布(單位:℃)

通過對比圖2(a)及圖3(a)可知,在自由對流過程中,流體攜帶熱量運動,極大地改變了整體溫度分布;對比圖3(a)及圖3(b)可知,基于FLAC3D軟件計算得到的溫度數值解與采用文獻[22]方法計算得到的解析解幾乎吻合。為進一步判斷兩者擬合度,提取監測點溫度分布并進行比對,結果見表1。由表1可知,數值解與解析解最大溫度差值出現在1號監測點,為3.4℃,誤差僅為3.89%。因此,基于FLAC3D軟件計算得到的數值解與解析解高度吻合,滲流傳熱計算的有效性及準確性得以驗證。

表1 監測點穩態溫度值

2.2 飽和土熱固結計算

飽和土在熱固結過程中受外部荷載、孔隙壓力、溫度荷載的同時作用,伴隨著孔隙水的排出、溫度的擴散,土體內部的有效應力變化較為復雜。白冰[23]基于飽和多孔介質熱-水-力完全耦合的控制方程,通過有限傅立葉變換及其逆變換,推導了土柱內部溫度、孔隙壓力和位移演化過程的解析表達式,并對一維熱彈性固結模型進行解析求解;孫致學等[24]利用該解析解進行了熱流固耦合模型的數值驗證。為驗證基于FLAC3D軟件求解熱流固耦合問題的有效性及準確性,本文引入文獻[24]中經典一維飽和土的熱彈性固結算例進行分析。計算模型及監測點布置情況如圖4所示,模型高度為7m,寬度為3m,厚度為1m,共計網格節點1240個,單元總數570個,監測點4個,分別設置在z坐標為0、3.5、5.0、7.0m處。

圖4 熱彈性固結計算模型

模型初始孔隙壓力為10kPa,初始溫度為10℃。力學計算邊界:模型頂面受向下的均布荷載作用,大小為10kPa;模型左右兩側為法向位移約束;底面為固定約束。滲流計算邊界:左右兩側及底面均為不透水邊界;頂面為自由滲漏邊界。傳熱計算邊界:頂面受固定溫度荷載作用,大小為60℃;左右兩側及底面均為絕熱邊界。具體計算參數如下:彈性模量為60MPa,泊松比為0.4,滲透率為4×10-6m/s,孔隙率為0.2,導熱率為0.836W/(m·K),比熱容為167.2J/(kg·K),熱膨脹系數為3×10-7K-1,密度為2000kg/m3,流體密度為1000kg/m3。

圖5為土體不同位置處豎向位移、孔隙壓力及溫度的解析解和數值解對比圖,可見基于FLAC3D軟件計算得到的數值解與文獻[24]中的解析解高度吻合,其有效性及準確性得以驗證。

3 模型應用與結果分析

3.1 工程背景

根據湖南省水利水電勘測設計研究總院2021年4月的《護城垸工程地質勘察報告》,護城垸位于洞庭湖區華容縣,西南臨藕池河東支,東靠華容河,垸內總面積364.6km2,總人口37.8萬人,是洞庭湖區重點堤垸。該垸一線大堤分華容河段和藕池河段,總長82.816km,其堤防級別為Ⅱ級,建筑物級別為Ⅱ級。其中,華容河堤段樁號為70+030～76+790的斷面堤身土質量較差,密實程度不均。1998年、2016年汛期該堤段堤身、堤腳出現散浸及滲水險情。

3.2 計算模型

以華容河堤段樁號為71+590的斷面為研究對象,建立堤防計算網格模型。堤頂高程為35.97m,寬6.12m,堤身高6.47m,內坡比為1︰4,外坡比為1∶2.58。堤身為素填土,堤基由粉質黏土、淤泥粉質黏土及粉細砂層組成。計算模型堤身及堤基采用八節點六面體單元劃分,共計網格節點數27882個,單元總數13640個。計算網格及分組情況如圖6所示,計算參數見表2。

圖6 計算網格及模型分組

為深入分析下游堤基及地表處滲透破壞及其傳熱特征,設置7個監測點用于監測各狀態變量隨時間變化的規律。測點1位于下游堤腳處,測點2、3、4位于下游堤腳附近地表處,測點5、6、7位于下游地表以下0.75m處的堤基表層,具體編號及位置如圖7所示。

圖7 模型監測點位置

3.3 計算方案及邊界條件

洞庭湖水位受季節及氣候的影響,水位高程變化復雜。為模擬堤防滲流實際情況,擬定計算工況為天然工況、初始工況及高水位工況。其中天然工況及初始工況為穩態計算工況,高水位工況為瞬態計算工況。

a.天然工況。為模擬堤防枯水期滲流及溫度分布,擬定堤防天然工況為上游無水的情況。滲流方面,根據實測資料擬定地下水位。傳熱方面,模型底部溫度為20℃,頂部由于受到氣溫的作用,擬定溫度為25℃。

b.初始工況。滲流方面,擬定上游水位為2020年汛期華容河水文站平均水位,高程為30.8m。傳熱方面,擬定堤頂、堤防背水側及下游地表溫度為2020年夏季岳陽市平均氣溫,大小為28℃;擬定上游湖水表面溫度為2020年汛期洞庭湖表層平均水溫[25],大小為25℃;上游堤坡及堤基表面溫度與水溫一致,從上至下逐級遞減,堤基表面溫度為22.5℃;模型底部溫度保持不變,大小為20℃。

c.高水位工況。據華容河水文站資料記載,華容河2020年汛期高水位持續時間長達40d以上。因此擬定計算時間為50d,以模擬長時高水位作用下堤防的異常滲漏情況。滲流方面,擬定堤防上游水位為華容河大堤保證水位,高程為35.5m。傳熱方面,擬定大氣溫度為30.6℃,堤頂及下游堤坡、地表處均為30.6℃的對流換熱邊界;上游堤基表面溫度不變,為22.5℃;堤坡溫度從下向上逐級遞增,水面處溫度為28℃。

3.4 計算結果分析

3.4.1滲流場分布特性

圖8為不同工況孔隙壓力分布圖。天然工況下,在流體自重作用下,孔隙壓力呈分層分布,從上至下孔隙壓力逐級遞增,水位線以下堤基土部位為飽和狀態,以上為非飽和狀態?？紫秹毫ψ畲笾禐?.12MPa,出現在堤基底部。初始工況下,在穩定水位的作用下,滲水進入了堤身及堤基內部,并與地下水連通?？紫秹毫Τ史謱臃植?且上游側整體較下游側大。在穩定高水位的作用下,堤內浸潤線抬升,堤身非飽和區大幅減小,堤防整體孔隙壓力左側大于右側,呈分層分布,孔隙壓力最大值增加至0.14MPa,分布于上游側堤基底部。

圖8 不同工況孔隙壓力分布(單位:MPa)

圖9為高水位工況水力梯度分布圖。堤防整體水力梯度左側大于右側,水力梯度最大值為2.4,分布于上游堤坡及堤腳處。此時下游堤腳處水力梯度為0.25。

3.4.2溫度場分布特性

圖10為不同工況溫度分布圖。天然工況下,在熱傳導作用下溫度場呈分層分布,且由表及里逐級遞減。溫度最大值分布在堤防表面,大小為25℃;初始工況下,在夏季高溫作用下,堤防表面持續升溫,溫度最大值達到28℃。在低溫水體的強制對流作用下,溫度最小值為22.5℃,分布于上游側地基內部。由于水溫低于堤基溫度,導致上游處溫度由表及里逐漸升高,堤內溫度低于大氣溫度,導致下游處溫度由表及里逐漸降低。高水位工況下,在強制對流作用下,下游堤腳處溫度逐漸下降。在較高氣溫作用下,堤頂及下游堤坡處溫度逐漸升高。隨著滲流量的不斷增加,下游堤腳處溫度梯度逐漸增大,與周圍地表及堤坡處溫差較大。此時下游堤坡處溫度為29.5℃,堤腳處為28.1℃。

3.4.3損傷場分布特性

圖11為根據式(2)得到的高水位工況損傷系數及塑性區分布圖。圖11(a)表明,上游堤身大面積部位損傷系數值達到1.0,代表該處損傷情況較為嚴重。下游堤腳處損傷系數大小為0.14。圖11(b)表明,塑性區與損傷系數分布位置與規律較為一致。上游側堤身內部出現連貫的塑性區分布,此時堤身存在較大的安全隱患。

圖11 高水位工況損傷系數和塑性區分布

3.4.4溫度場與滲流場關聯關系分析

由于不同位置處滲流場及溫度場的分布形態均具有一定的差異性,難以把握滲流場與溫度場的內在聯系,需要針對性選取下游堤基及地表處測點,進行各狀態變量的歷時分析。

圖12為堤防不同測點的孔隙壓力、水力梯度及溫度歷時曲線。圖12(a)表明,由于下游堤腳及附近地表直接暴露在大氣中,滲水的及時排出使得孔隙壓力得到釋放,導致下游堤腳及附近地表處測點的孔隙壓力曲線均低于堤基表層測點。隨著時間的推移,下游堤腳及附近地表測點的孔隙壓力隨時間先急劇減小后緩慢減小,堤基表層測點的孔隙壓力逐漸減小。

圖12 各狀態變量歷時曲線

圖12(b)表明,下游堤腳及附近地表處與大氣接觸,壓力梯度較大,導致該處測點的水力梯度曲線均高于堤基表層測點。從整體上看,堤基表層測點的水力梯度隨時間緩慢增大,下游堤腳及附近地表的水力梯度呈先急劇增大后緩慢減小的趨勢。出現這種現象的原因是下游堤腳及附近地表處孔隙壓力消散較快,較大的壓力梯度使得水力梯度急劇增大;隨著孔隙壓力趨于穩定,下游堤腳及附近地表測點的水力梯度逐步減小。據《護城垸工程地質勘察報告》,堤基表層土體的臨界水力梯度為0.32。隨著水力梯度超過臨界值,下游堤腳及附近地表處不同位置先后出現局部滲透破壞。首先是測點2,隨之是測點3、4,最后是測點1。

圖12(c)表明,在大氣高溫條件下,下游堤腳及附近地表處測點的溫度曲線均高于堤基表層測點。堤基表層土體未直接與大氣接觸,受氣溫影響較小,因此溫度的變化主要來源于內部水體的強制對流傳熱,且隨著時間的推移,溫度呈下降趨勢。下游堤腳及附近地表土體受到滲水強制對流傳熱及大氣高溫的雙重影響,該處溫度變化趨勢較為復雜。首先,由于短時間內下游堤腳及附近地表處還未產生較大流量的滲水,在大氣高溫條件下,該處溫度升高;隨著滲流現象的產生,滲水的強制對流傳熱效應顯著,下游堤腳及附近地表處溫度急劇降低;隨著該處孔隙壓力的消散,流速緩慢減小,強制對流傳熱效應逐漸衰弱。最終在滲水強制對流傳熱及大氣高溫的共同作用下,下游堤腳及附近地表處的溫度趨于平衡。

對比圖12(b)和圖12(c)可知,隨著堤基表層不同位置先后發生局部滲透破壞,各測點的溫度值出現差異,差值約為0.18℃。因此可通過捕捉不同位置處的溫度差來預測局部滲透破壞的發生。值得一提的是,此時溫度差較小只能代表下游的局部滲透破壞,無法代表堤防整體滲透破壞或結構破壞的程度。

3.4.5溫度場與堤防結構破壞表象關聯關系分析

為進一步探明溫度與堤防結構破壞表象的關聯關系,采用FLAC3D軟件內置的強度折減法進行堤防抗滑穩定性分析。圖13為高水位工況下,通過強度折減計算得到的不同時刻的極限狀態最大剪應變分布云圖。由圖13可知,隨著時間的推移,最大剪應變數值不斷增加。第40天前,堤防潛在滑裂面位于下游側。隨著上游堤身內部損傷面積不斷擴大,第40天后堤防潛在滑裂面位置發生轉移,上游堤坡發生滑坡風險較大。

圖13 極限狀態最大剪應變

圖14為堤防抗滑穩定安全系數歷時曲線。在長時高水位作用下,堤身及堤基內部逐漸發生滲透破壞,堤防安全系數緩慢減小。從第40天開始,堤防安全系數急劇減小。根據GB 50286—2013《堤防工程設計規范》,Ⅱ級堤防正常運行條件下的安全系數允許值為1.35。由圖14可知,前期堤防安全系數為1.53,處于穩定狀態。第47天降低至允許值,此時發生大面積滑坡乃至潰堤的風險較大。通過對比可知,不同測點間溫度差于第45天后突增,最大可達0.45℃,此時堤防整體安全系數急劇減小。

圖14 安全系數歷時曲線

綜上所述,無論是下游堤腳處發生局部滲透破壞,或是上游堤身內部發生剪切破壞乃至滑坡,都會提前在下游堤腳及附近地表不同位置處出現明顯的溫度差值。隨著滲透破壞的發展、上游堤身內部塑性區的延伸,下游堤腳及附近地表處的溫度差值會逐漸增大。

4 結 論

a.堤壩滲漏過程中,下游堤基及地表處的溫度演變有階段性特征。首先,短時間內下游堤腳及附近地表處還未產生較大流量的滲水,在較高氣溫條件下,該處溫度升高;隨著滲漏現象的產生,滲水的強制對流傳熱效應顯著,下游堤腳及附近地表處溫度急劇降低;隨著該處孔隙壓力的消散,流速緩慢減小,強制對流傳熱效應逐漸衰弱;最終下游堤腳及附近地表處的溫度趨于平衡。

b.隨著下游堤腳及地表處水力梯度的不斷增大,不同位置處先后發生局部滲透破壞,出現較小溫度差,約0.18℃;隨著時間的推移,堤身內部塑性區不斷延伸,堤防安全穩定性逐漸下降。堤防安全穩定系數低于臨界值前,下游地表不同位置的溫度差進一步增大,達到0.45℃。

c.無論是下游堤腳處發生局部滲透破壞,或是上游堤身內部發生剪切破壞乃至滑坡,在下游堤腳及附近地表不同位置處均會提前出溫度差值。隨著滲透破壞的發展,該溫度差值會逐漸增大。