基于Ren模型的土石壩集中滲漏流熱耦合研究

倪 楓 任 杰,2 張文兵 王 帆

（1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710048；2.水利部堤防安全與病害防治工程技術(shù)研究中心,鄭州 450003；3.河海大學(xué) 水利水電學(xué)院,南京 210098）

1 研究背景

土石壩滲漏一直以來(lái)都是土石壩隱患診斷的重點(diǎn)問(wèn)題,對(duì)于土石壩的滲漏監(jiān)測(cè)主要方法有同位素示蹤法、電阻率成像法、滲流熱監(jiān)測(cè)技術(shù)以及彈性波探測(cè)等,這幾種方法各自有優(yōu)缺點(diǎn)[1].同位素示蹤雖然可以準(zhǔn)確地示蹤出土石壩滲漏軌跡,但是布置造價(jià)過(guò)高,難以實(shí)現(xiàn);電阻率成像法等方法缺乏試驗(yàn)支撐,并且圖像難以解譯,如果應(yīng)用到實(shí)際工程中比較困難.相比之下溫度作為天然的示蹤劑可以更方便地應(yīng)用到實(shí)際工程中,因此滲流熱監(jiān)測(cè)技術(shù)是目前研究土石壩滲漏的一個(gè)重要手段.滲流熱監(jiān)測(cè)技術(shù)的理論基礎(chǔ)是巖土體中滲流場(chǎng)與溫度場(chǎng)的耦合作用[2],該技術(shù)的原理是介質(zhì)中的水分運(yùn)移控制著土石壩溫度場(chǎng)的分布,壩體溫度場(chǎng)的時(shí)空分布反過(guò)來(lái)能夠反映壩體滲流場(chǎng)的特征.

土壤導(dǎo)熱系數(shù)是土壤熱力學(xué)特性的一個(gè)重要指標(biāo),是影響土壤傳熱過(guò)程中溫度分布和能量變化的主要因素,也是研究其物理過(guò)程,如流-熱耦合、物質(zhì)傳遞、氣體擴(kuò)散等的基礎(chǔ)[3-4].Lu等[5]根據(jù)前人的研究,在2007年將土壤大致劃分為粗土和細(xì)土,并引入了和砂粒土含量相關(guān)的經(jīng)驗(yàn)常數(shù).2015年,Lu和Dong等[6]提出了閉合土體的熱導(dǎo)率經(jīng)驗(yàn)方程,該方程將土壤水分特征曲線與導(dǎo)熱系數(shù)結(jié)合起來(lái),可以較準(zhǔn)確地預(yù)測(cè)不同類型土體的導(dǎo)熱系數(shù).考慮到有機(jī)質(zhì)對(duì)熱導(dǎo)率的影響和土壤顆粒組分對(duì)形狀因子的影響,Ren等[7]于2019年提出一個(gè)新的經(jīng)驗(yàn)?zāi)Ｐ?與之前的熱傳導(dǎo)系數(shù)模型相比,Ren（2019）模型充分反映了土體有機(jī)質(zhì)成分和顆粒成分對(duì)形狀因子的影響.

本文基于COMSOL Multiphysics軟件建立Ren（2019）模型的土石壩流-熱耦合數(shù)學(xué)模型,以典型二維心墻土石壩的數(shù)值模擬作為算例,將模擬結(jié)果與吳志偉等[8]基于Lu（2007）模型的土石壩流熱耦合模型進(jìn)行對(duì)比,在此基礎(chǔ)上模擬了土石壩心墻在不同集中滲漏情況下溫度場(chǎng)和滲流場(chǎng)的變化,研究結(jié)果可為新的導(dǎo)熱系數(shù)在滲流熱監(jiān)測(cè)技術(shù)的應(yīng)用提供參考.

2 數(shù)學(xué)模型

2.1 滲流場(chǎng)模型

滲流場(chǎng)采用理查茲方程表示[9-10]:

式中:ρw為水的密度;g為重力加速度;Se為土體相對(duì)飽和度;Ss為彈性貯水率;p為壓強(qiáng);θ為含水量;ks為介質(zhì)飽和滲透系數(shù);kr（θ）為非飽和帶相對(duì)滲透率[0≤kr（θ）≤1];μ（T）為水的動(dòng)力黏度[10];z代表需要計(jì)算位置的高程;Cm為土體容水度,代表水分特征曲線斜率的負(fù)倒數(shù);Qm代表水流源匯項(xiàng);?代表拉普拉斯算子.

用van Genuchten模型描述非飽和帶土壤水分特征曲線[10]:

式中:θr為殘余含水率;θs為飽和含水率.其中Se為:

式中:Se為土壤相對(duì)飽和度;α和nv為VG 模型參數(shù);hp為壓力水頭;α是水分特征曲線進(jìn)氣值的倒數(shù);nv是分特征曲線坡度的指示參數(shù),m＝1－1/nv,在非飽和帶壓力水頭hp等于負(fù)壓水頭hc[11].

土體容水度Cm和相對(duì)滲透率kr分別采用下述經(jīng)驗(yàn)公式表示:

2.2 溫度場(chǎng)模型

溫度場(chǎng)采用熱對(duì)流傳熱方程描述[12]:

式中:T為溫度;Cw、Cs分別為流體和土體的比熱容;λ為有效導(dǎo)熱系數(shù);DH為水動(dòng)力彌散系數(shù);u為平均速度;Qs為熱量源匯項(xiàng).

COMSOL中自帶計(jì)算傳熱模塊是用熱量運(yùn)移方程進(jìn)行計(jì)算的,用PDE 模塊來(lái)代替COMSOL 自帶的傳熱模型進(jìn)行計(jì)算,可靈活修改導(dǎo)熱系數(shù)經(jīng)驗(yàn)?zāi)Ｐ?

COMSOL中提供的系數(shù)型PDE方程如下[14]:

式中:ea、da、c、α、γ、β、a、q和h均為自定義系數(shù);Ω為求解域,Ω的偏微分是其外邊界;n是外法線方向;f、g、r分別為求解域和邊界上的源項(xiàng).

2.3 Ren（2019）土體導(dǎo)熱系數(shù)模型

Ren等[7]考慮有機(jī)質(zhì)含量對(duì)導(dǎo)熱系數(shù)的影響及土壤粒組成分對(duì)形狀因子的影響,提出一個(gè)新的導(dǎo)熱系數(shù)經(jīng)驗(yàn)?zāi)Ｐ?

式中:λsat和λdry與砂土含量和干土容重（γd）有關(guān),當(dāng)砂含量0＜Csand＜1 時(shí),干土容重11＜γd＜20.γd（k N/m3）和ρb（g/cm3）之間的關(guān)系式為:

式中:g＝9.8 m2·s－1,是重力加速度.

Ke與θ的新關(guān)系如下:

式中:α和β是λ（θ）的形狀因子,它們的取值會(huì)影響導(dǎo)熱系數(shù)曲線的斜率和增長(zhǎng)速度,表達(dá)式如下:

式中:Com為有機(jī)質(zhì)的質(zhì)量比（g·kg－1）;Cclay、Csilt、Csand分別表示黏粒、粉粒、砂粒的質(zhì)量分?jǐn)?shù)（%）.

3 實(shí)例應(yīng)用

3.1 研究區(qū)域概況

選取湖北省黃岡市大山背水庫(kù)進(jìn)行土石壩的集中滲漏研究[15].大山背土石壩工程的基本概況如下:壩體為黏土心墻代料組合壩,壩頂長(zhǎng)134 m、寬7 m;壩頂高程108.12 m,壩高23.2 m;心墻頂高程106.32 m,頂寬2.0 m,邊坡為1∶0.15.壩體填筑材料假設(shè)為粉質(zhì)黏土,心墻材料為黏土,壩體排水體采用褥墊式水平排水,排水褥墊層高2 m,寬22 m,褥墊層材料為礫石.為了說(shuō)明不同滲漏方式對(duì)壩體溫度場(chǎng)和滲流場(chǎng)的影響,本文假設(shè)在某一特定時(shí)刻,心墻形成0.02 m 寬的裂縫,并假設(shè)心墻裂縫分別存在3個(gè)不同高度（6.10、10.10、14.10 m）,如圖1所示.

圖1 二維心墻土石壩結(jié)構(gòu)示意圖

3.2 模型邊界與具體參數(shù)

3.2.1 邊界條件

滲流場(chǎng)方面,壩體閃游定水頭為104.92 m,邊界為A-E-D,壩體下游定水頭86.92 m,邊界為B-C,其余邊界均為零通量邊界.

溫度場(chǎng)方面,由于湖北省屬于華中地區(qū)且本研究對(duì)象壩高在23.2 m,根據(jù)文獻(xiàn)[8,16]的數(shù)據(jù)分析和論證,庫(kù)水溫度隨高程變化基本可忽略,因此簡(jiǎn)化邊界A-E溫度按當(dāng)?shù)厮畮?kù)的平均水溫取10℃;A-B作為絕熱邊界;B-C-D-E作為大氣邊界,溫度取日平均氣溫.當(dāng)?shù)販囟入S時(shí)間周期性變化的傅里葉級(jí)數(shù)函數(shù)[14]如下:

式中:T0為平均溫度;k為傅里葉級(jí)數(shù)的階數(shù);Ak和Bk為傅里葉級(jí)數(shù)的系數(shù);L為周期.

這里采用湖北省黃岡市的氣溫?cái)?shù)據(jù)（數(shù)據(jù)來(lái)源于中國(guó)氣象數(shù)據(jù)網(wǎng)站）并根據(jù)傅里葉函數(shù)來(lái)取5階傅里葉級(jí)數(shù)就可以滿足精度要求.擬合得氣溫回歸模型系數(shù)為:T0＝16.37、A1＝－10.28、A2＝－0.541 7、A3＝0.266 7、A4＝－0.291 7、A5＝－0.088 22、B1＝－7.133、B2＝－0.216 5、B3＝0.3、B4＝－0.043 3、B5＝0.083 44.將系數(shù)代入式（15）即可求得溫度隨時(shí)間的周期性規(guī)律.

3.2.2 模型參數(shù)

結(jié)合當(dāng)?shù)貕误w材料并參考文獻(xiàn)[15],給出模型滲流場(chǎng)及溫度場(chǎng)的計(jì)算參數(shù),見表1.區(qū)域1、2、3分別代表圖1中的壩體、心墻、排水褥墊層,3個(gè)區(qū)域的材料主體分別是粉質(zhì)黏土、黏土、礫石,心墻裂縫區(qū)域的飽和水力傳導(dǎo)率參照排水褥墊層取10－2m/s.

表1 心墻土石壩流-熱耦合模型計(jì)算參數(shù)

4 結(jié)果與討論

4.1 模型對(duì)比分析

為了開展基于Ren（2019）模型的土石壩流-熱耦合研究,有必要驗(yàn)證該模型的可靠性和準(zhǔn)確性,因此參照吳志偉等[8]進(jìn)行的基于Lu模型土石壩集中滲漏研究,在與之模型邊界條件和參數(shù)都相同的情況下模擬基于Ren（2019）模型的土石壩集中滲漏研究.基本工況為水庫(kù)從0 d時(shí)刻開始蓄水,在0～2 190 d之間處于蓄水之后的正常工作期,隨后在第2 190 d時(shí)心墻裂縫在距壩底10.08 m 處開始發(fā)生集中滲漏,并持續(xù)200 d.因此設(shè)置第一個(gè)計(jì)算區(qū)段為2 190 d,之后發(fā)生集中滲漏,計(jì)算時(shí)間為200 d,取時(shí)間步長(zhǎng)1 d.通過(guò)模擬計(jì)算出心墻裂縫出口處的溫度和滲透流速變化,將結(jié)果與吳志偉等[8]計(jì)算出的結(jié)果進(jìn)行對(duì)比如圖2所示.

圖2 心墻裂縫出口處滲透流速及溫度曲線對(duì)比圖

從圖2可以看出,基于兩個(gè)不同導(dǎo)熱系數(shù)模型計(jì)算出的心墻裂縫出口處溫度和滲透流速變化趨勢(shì)接近一致,兩個(gè)模型在第2 190 d之前的心墻裂縫出口處溫度變化趨勢(shì)基本一致;發(fā)生集中滲漏以后,兩個(gè)模型心墻裂縫出口處溫度變化基本一致,在經(jīng)過(guò)40 d左右溫度同時(shí)達(dá)到最低,最低溫度與庫(kù)水溫度基本一致,之后受大氣溫度影響逐漸上升;裂縫出口處滲透流速,兩個(gè)模型在第2 190 d之前基本都接近于0 m/s,在第2 190 d時(shí)突變后的滲透流速同屬于一個(gè)數(shù)量級(jí)10－3m/s.吳志偉等[8]計(jì)算出的心墻裂縫出口處滲透流速在第2 190 d突增后一直維持在1.33×10－3m/s,基于Ren模型計(jì)算出的心墻裂縫出口處滲透流速在第2190 d突增到1.48×10－3m/s,之后200 d逐漸降低穩(wěn)定到1.31×10－3m/s,滲透流速在突增之后逐漸降低的原因是發(fā)生集中滲漏之后心墻上游面水頭逐漸降低,后滲透流速達(dá)到穩(wěn)定.

通過(guò)對(duì)上述兩個(gè)模型計(jì)算結(jié)果的綜合分析,得出基于Ren（2019）模型的土石壩流-熱耦合模型可以準(zhǔn)確地反映土石壩溫度場(chǎng)和滲流場(chǎng)的變化過(guò)程.

4.2 不同集中滲漏條件下土石壩溫度場(chǎng)及滲流場(chǎng)的變化特征

為研究不同集中滲漏條件下土石壩溫度場(chǎng)及滲流場(chǎng)的變化規(guī)律及特征,建立基于Ren（2019）模型的土石壩飽和-非飽和流-熱耦合模型.設(shè)置土石壩在正常運(yùn)行8a（2920 d）之后發(fā)生集中滲漏,假設(shè)土石壩心墻分別在距壩底6.08、10.08、14.08 m 形成水平貫穿裂縫,對(duì)應(yīng)工況1、2、3,集中滲漏計(jì)算時(shí)間為300 d,時(shí)間步長(zhǎng)1 d.

各工況以及無(wú)集中滲漏的壩體滲流場(chǎng)如圖3所示,圖中紅線代表壩體的浸潤(rùn)線.由圖3（a）可知,浸潤(rùn)線在穿過(guò)心墻后急劇降低,之后基本與排水褥墊層持平,說(shuō)明此時(shí)壩體的滲流庫(kù)水基本都能從排水褥墊層排到下游,有效地保證了土石壩的安全運(yùn)行.由圖3（b）、（c）、（d）與圖3（a）對(duì)比可知,當(dāng)心墻發(fā)生集中滲漏時(shí),浸潤(rùn)線在穿過(guò)心墻下游表面時(shí)會(huì)隨著心墻裂縫的高度變化而發(fā)生變化,心墻裂縫距壩底越高,壩體浸潤(rùn)線在穿過(guò)心墻下游表面時(shí)的位置越高.各工況以及無(wú)集中滲漏的壩體溫度場(chǎng)如圖4所示.在沒有集中滲漏的情況下,壩體溫度場(chǎng)是從心墻前庫(kù)水的低溫區(qū)逐步向下游的高溫區(qū)過(guò)渡,等溫線分布比較均勻,而在集中滲漏之后,心墻前庫(kù)水通過(guò)圖4（b）、（c）和（d）所示的不同高度的滲漏通道分別進(jìn)入心墻下游非飽和帶,等溫線在心墻裂縫進(jìn)口處較為密集,且低溫庫(kù)水滲過(guò)集中滲漏通道后使心墻后壩體溫度降低,相應(yīng)的大氣溫度對(duì)壩體溫度影響的區(qū)域減少.

圖3 不同時(shí)刻壩體壓力水頭分布圖（單位:m）

圖4 不同集中滲漏條件下壩體溫度場(chǎng)的變化（單位:℃）

比較上述4種情況的土石壩溫度場(chǎng)和滲流場(chǎng),得出壩體溫度場(chǎng)分布在一定程度上受集中滲漏的影響,而且不同高度的集中滲漏通道對(duì)壩體溫度場(chǎng)的影響不同,可以通過(guò)測(cè)溫來(lái)反演滲流場(chǎng)去確定集中滲漏通道的具體位置.

4.3 心墻裂縫出口處滲透流速及溫度變化規(guī)律分析

圖5（a）為不同高度的心墻裂縫出口處滲透流速變化對(duì)比圖.由圖5可以看出,不同高度的心墻裂縫出口處滲透流速在2 920 d之前一直保持接近0 m/s,在2 920 d時(shí)3種工況出口處滲透流速發(fā)生突變,且都達(dá)到10－3m/s以上,之后隨著時(shí)間推移逐漸降低并在300d后基本全部達(dá)到穩(wěn)定.從圖5可以看出對(duì)不同高度心墻裂縫突變后的滲透流速大小不同,且心墻裂縫越靠近壩底,滲透流速突變?cè)酱笄疫_(dá)到穩(wěn)定之后的滲透流速也越大.

圖5（b）為不同高度的心墻裂縫出口處溫度變化對(duì)比圖.由圖5（b）可以看出,不同高度的心墻裂縫出口處溫度在發(fā)生集中滲漏的瞬間都發(fā)生了溫度突降.通過(guò)分析圖4中無(wú)滲漏情況的壩體溫度場(chǎng),可以得出工況1、2、3在集中滲漏發(fā)生前的溫度分別為14.4℃、11.5℃、11.3℃.當(dāng)集中滲漏在3種工況分別發(fā)生時(shí),對(duì)于工況1來(lái)說(shuō),裂縫出口處溫度突降以后隨著時(shí)間推移一直降低到10℃左右,而后溫度隨大氣溫度變化而變化;對(duì)于工況2和工況3來(lái)說(shuō),裂縫出口處溫度突降以后隨著時(shí)間推移先升高后降低到10℃左右,而后溫度也隨大氣溫度變化成正弦變化.工況2和工況3在裂縫出口處溫度在突降以后先升高后降低是由于外界大氣溫度在壩體內(nèi)傳遞具有滯后性造成的.

圖5 不同高度心墻裂縫出口處溫度與滲透流速變化對(duì)比

綜上,當(dāng)心墻出現(xiàn)裂縫時(shí),心墻下游表面的溫度和滲透流速都會(huì)發(fā)生突變.裂縫距壩底高度越高,發(fā)生集中滲漏時(shí)溫差變化越大,滲透流速變化越小,裂縫距壩底高度越低,發(fā)生集中滲漏時(shí)滲透流速變化越大,溫差變化越小.

5 結(jié)論

1）基于COMSOL Multiphysics軟件的Ren（2019）模型適用于土石壩飽和-非飽和流-熱耦合研究,可以準(zhǔn)確地描述土石壩溫度場(chǎng)和滲流場(chǎng)的動(dòng)態(tài)變化過(guò)程.

2）集中滲漏通過(guò)改變滲流場(chǎng)來(lái)改變土石壩溫度場(chǎng)的分布,可以通過(guò)溫度場(chǎng)去反演滲流場(chǎng)去確定集中滲漏通道的具體位置.

3）心墻裂縫距壩底高度越高,發(fā)生集中滲漏時(shí)溫差變化越大,滲透流速變化越小.在實(shí)際工程中可以在壩體內(nèi)部鋪設(shè)光纖來(lái)監(jiān)測(cè)壩體的溫度場(chǎng)變化,從而精確地確定發(fā)生集中滲漏的位置,為實(shí)際工程安全提供保障.