短時(shí)凍區(qū)氣候下土坡淺層水熱狀況的數(shù)值分析

闕 云，陳曉鵬，榮 耀

(1.福州大學(xué) 土木工程學(xué)院，福建 福州350108;2.江西農(nóng)業(yè)大學(xué) 工學(xué)院，江西 南昌330045)

0 引 言

福建北部既有殘積土邊坡廣泛在每年降雨量很少的冬季出現(xiàn)淺層垮塌現(xiàn)象，嚴(yán)重影響公路正常運(yùn)營(yíng)，但破壞原因始終存在困惑?，F(xiàn)場(chǎng)調(diào)查表明，該區(qū)域內(nèi)土質(zhì)邊坡冬季失穩(wěn)原因與地處短時(shí)凍區(qū)所引發(fā)的邊坡淺層短時(shí)反復(fù)凍融息息相關(guān)。因此，開展短時(shí)凍區(qū)氣候下的土坡淺層水熱狀況研究具有重要意義。事實(shí)上，國(guó)內(nèi)外眾多學(xué)者針對(duì)凍土工程相關(guān)水熱力變化展開了大量研究，取得了豐碩成果［1-9］，如:Harlan［1］(1973)提出了土凍結(jié)過程中水、熱遷移相互作用的耦合模型并得到了廣泛應(yīng)用。熱力學(xué)模型是Duquennoi 在1985 年首次提出［2］。何敏等［4］在飽和凍土多孔多相介質(zhì)理論的基礎(chǔ)上，建立了飽和凍土溫度場(chǎng)、水分場(chǎng)、變形場(chǎng)耦合分析數(shù)值模型，并進(jìn)行了驗(yàn)證。宋存牛［5］對(duì)國(guó)內(nèi)外在水熱力耦合作用機(jī)理、耦合作用方程方面的研究進(jìn)行了綜評(píng)，指出了各種研究結(jié)果之間差異、不同凍脹模式建立的依據(jù)、研究方法以及適用條件。然而以上研究多數(shù)集中于季節(jié)及多年凍土，針對(duì)短時(shí)凍土區(qū)邊坡淺層水熱變化研究還較少。相比與多年凍土與季節(jié)凍土，南方短時(shí)凍土無論是其深度、凍結(jié)及熱融過程都不盡相同?；诖耍疚囊罁?jù)南方短時(shí)凍區(qū)氣候特點(diǎn)，針對(duì)短時(shí)凍區(qū)內(nèi)殘積土邊坡淺層水熱耦合進(jìn)行數(shù)值分析，探究土坡淺層水熱遷移狀況。

1 水熱耦合模型及程序編制

1.1 水熱耦合模型

水熱耦合模型采用Harlan 模型［6］:

式中，θ、θi分別為含水率和含冰率;D( θ) 、K(θ)分別為水分?jǐn)U散率和導(dǎo)水系數(shù);ρi、ρl分別為冰、水的密度;C、λ、Lf分別為比熱容、導(dǎo)熱系數(shù)及凍融潛熱;T、t 分別為溫度及時(shí)間。其中未凍水含量計(jì)算式通過測(cè)溫法室內(nèi)試驗(yàn)測(cè)得:θu=0.519T-1.18。

1.2 有限差分模型

基于交替隱式差分(ADI)的隱式差分方法對(duì)式(1)進(jìn)行求解，采用MATLAB 語言編制程序，差分模型如圖1 所示，其中坡高為7.5 m，坡底寬5 m。

1.3 邊界條件

基于第二類邊界條件確定坡面邊界，主要考慮蒸發(fā)、太陽輻射、凈輻射、對(duì)流換熱、蒸發(fā)潛熱。坡底坡高邊界采用隔水隔溫邊界。其中蒸發(fā)強(qiáng)度采用Campbell 模型［10］，如式(2)。

式中，E 為蒸發(fā)強(qiáng)度;ρvs、ρva分別為地表、空氣水汽密度;rv為阻力系數(shù)。

熱邊界主要由蒸發(fā)潛熱、太陽總輻射、凈輻射及對(duì)流換熱累加而得。其中蒸發(fā)潛熱采用式(3)［8］計(jì)算。

式中，L=2.49463×109-2.247×106Ts;Ts為土壤表層溫度。

太陽輻射采用HIPERPAVⅡ模型［11］計(jì)算，如式(4)。

式中，qs為太陽輻射吸收熱流，W/m2;αs為太陽輻射的吸收率;If為考慮太陽在24 h 中的輻射角度的輻射強(qiáng)度系數(shù)，一般假設(shè)其服從正弦分布，晚上為零;qsolar為太陽輻射峰值。

凈輻射計(jì)算依據(jù)Stefan-Boltzman 定律確定［6］，其中當(dāng)表層結(jié)冰時(shí)運(yùn)用式(5)計(jì)算，表層不結(jié)冰時(shí)由式(6)計(jì)算。

式中，Lsp、Ls分別為雪層和土壤的表層凈長(zhǎng)波輻射量;TK為絕對(duì)溫度;εac、σ 分別為熱輻射系數(shù)(土壤取值為0.65，積雪取值為0.8)、常數(shù)(5.6697×10-8W/m2·K4)。

對(duì)流換熱［9］采用式(7)計(jì)算，對(duì)流換熱系數(shù)由式(8)確定。

式中，hc為對(duì)流換熱系數(shù);C 為常數(shù)，取決于形狀和熱流條件，當(dāng)?shù)乇頊囟雀哂诳諝鈺r(shí)取1.79，當(dāng)?shù)乇頊囟鹊陀诳諝鈺r(shí)取0.89;w 為風(fēng)速，m/s;Ta為空氣溫度;Ts為坡表溫度。

選取2013 年12 月27 日福建省武夷山市省道303 線K364+900 處現(xiàn)場(chǎng)實(shí)測(cè)氣象數(shù)據(jù)作為本文模擬基準(zhǔn)氣象邊界條件，具體如表1 所示。

表1 氣象邊界條件參數(shù)Tab.1 Meteorological parameters of the boundary conditions

1.4 計(jì)算參數(shù)

土水特征曲線采用濾紙法確定，并運(yùn)用(Van Genuchten 模型，以下簡(jiǎn)稱VG 模型)擬合，VG 模型方程［12］如式(9)所示。

式中，θ 為土壤含水率;θr、θs分別為殘余含水率以及飽和含水率(由試驗(yàn)測(cè)得分別為0.10、0.38);H 為土壤負(fù)壓;α、n、m(m=1-1/n)為擬合參數(shù)，由試驗(yàn)結(jié)果擬合可得α=0.0163、n=1.7445。

未凍土導(dǎo)水系數(shù)K(θ)與凍土導(dǎo)水系數(shù)Ki(θ)分別由式(10)和式(11)［3］計(jì)算。

式中，Se為有效飽和度;Ks為飽和導(dǎo)水率，由室內(nèi)試驗(yàn)可得土壤飽和導(dǎo)水率為0.002 m/h;I 為經(jīng)驗(yàn)常數(shù)，可近似表示為I=1010θi。

未凍土水分?jǐn)U散系數(shù)由式(12)計(jì)算，凍土水分?jǐn)U散系數(shù)Di(θ)由式(13)［3］計(jì)算。

土壤熱參數(shù)取值:未凍土骨料顆粒比熱為0.84 kJ/(kg·K);凍土骨料顆粒比熱為0.77 kJ/(kg·K);土壤中水和冰的比熱分別取4.182 kJ/(kg·K)、2.09 kJ/(kg·K)［11］，凍土導(dǎo)熱系數(shù)λf與未凍土導(dǎo)熱系數(shù)λu分布按式(14)、(15)［13］計(jì)算得到。

式中，rd為土質(zhì)干密度

1.5 初始條件

模擬時(shí)設(shè)定初始含水率為0.25，由室內(nèi)凍結(jié)溫度試驗(yàn)得，該含水率下土壤凍結(jié)溫度約為-1.9 ℃。由現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果［14］確定初始溫度分布為:

式中，x 為坡面水平深度(m)，y 為坡面垂直高度(m)。

2 計(jì)算結(jié)果分析

2.1 基準(zhǔn)條件下邊坡淺層水熱狀況

依據(jù)表1 邊界條件，通過數(shù)值模擬可得邊坡坡中(坡高4 m)淺層溫度及含冰率變化如圖2 所示。由圖2 可見:當(dāng)氣溫由零上逐漸降低至零下時(shí)，邊坡淺層由表至內(nèi)逐層凍結(jié)，當(dāng)日最低氣溫達(dá)至-4.3℃時(shí)，土坡淺層最大凍深約為5 cm。整個(gè)凍期下，土坡淺層結(jié)冰時(shí)間段為4:00 ～8:20，原因是由于7:00以后由于輻射影響淺層冰開始消融。邊坡淺層含冰率峰值約為0.14，出現(xiàn)在7:00 時(shí)刻，位置并非出現(xiàn)在邊坡表層，而是在淺層2 cm 左右位置。

圖2 基準(zhǔn)模擬計(jì)算結(jié)果Fig.2 Benchmark simulation results

2.2 日最低氣溫對(duì)淺層水熱狀況的影響

依據(jù)閩北地區(qū)歷史氣溫資料，選取日最低氣溫分別為-7 ℃、-5 ℃、-4.3 ℃，其他參數(shù)與氣象邊界參見基準(zhǔn)條件(見表1)，計(jì)算結(jié)果如圖3 至圖4 所示(其中日最低氣溫為-4.3 ℃的計(jì)算結(jié)果見圖2)。

圖3 日最低氣溫為－7 ℃時(shí)計(jì)算結(jié)果Fig.3 Calculation results when the daily minimum temperature is －7 ℃

圖4 日最低氣溫為－5 ℃時(shí)計(jì)算結(jié)果Fig.4 Calculation results when the daily minimum temperature is －5 ℃

由圖2、圖3 和圖4 可見:三種日最低氣溫條件下，邊坡淺層最大凍深分別約為11 cm、7 cm 及5 cm。當(dāng)日最低氣溫為-7 ℃時(shí)，其凍結(jié)時(shí)段為2:00 ～11:00，其余兩種工況凍結(jié)時(shí)間較為接近，均在4:00 ～8:20。三種工況下，邊坡淺層含冰率峰值分別為0.17、0.15、0.14，峰值出現(xiàn)時(shí)刻均在7:00 ～7:20，位置依次在淺層7 cm、4 cm、2 cm，由此可見最大含冰率并非出現(xiàn)在表層，而是隨著凍深的增大而逐漸加深。

2.3 風(fēng)速對(duì)邊坡淺層水熱狀況的影響

依據(jù)閩北地區(qū)歷史氣溫資料，選取平均風(fēng)速分別為10 m/s、5.6 m/s、2 m/s，其他參數(shù)與氣象邊界參照基準(zhǔn)條件(見表1)。計(jì)算結(jié)果如圖5、圖6 所示(其中平均風(fēng)速為5.6 m/s 的計(jì)算結(jié)果見圖2)。

圖5 平均風(fēng)速為10 m·s－1時(shí)計(jì)算結(jié)果Fig.5 Calculation results when the average wind velocity is 10 m·s－1

圖6 平均風(fēng)速為2 m·s－1時(shí)計(jì)算結(jié)果Fig.6 Calculation results when the average wind velocity is 2 m·s－1

由圖2、圖5 和圖6 可見:三種風(fēng)速條件下，邊坡淺層最大凍深分別約為7 cm、5 cm 及5 cm。從凍結(jié)時(shí)間段中可以看出，平均風(fēng)速為10 m/s 時(shí)，凍結(jié)時(shí)段為4:00 ～9:20，平均風(fēng)速為2 m/s 時(shí)，凍結(jié)時(shí)段為5:00 ～7:40，說明較強(qiáng)的風(fēng)速有利于延長(zhǎng)土坡淺層結(jié)冰時(shí)間。當(dāng)平均風(fēng)速為10 m/s 時(shí)，含冰率峰值為0.14，當(dāng)平均風(fēng)速為2 m/s 時(shí)，含冰率峰值為0.12，兩種工況下含冰率峰值出現(xiàn)位置均出現(xiàn)在2 cm 左右，說明風(fēng)速對(duì)含冰率峰值及其出現(xiàn)位置影響不大。

2.4 輻射對(duì)邊坡淺層水熱狀況的影響

依據(jù)閩北地區(qū)現(xiàn)場(chǎng)輻射監(jiān)測(cè)結(jié)果，分別選取輻射峰值為200 W/m2、539 W/m2、800 W/m2，其他參數(shù)與氣象邊界參照基準(zhǔn)條件(表1)。計(jì)算結(jié)果如圖7、圖8 所示(其中輻射峰值為539 W/m2的計(jì)算結(jié)果見圖2)。

圖7 輻射峰值為200 W·m－2計(jì)算結(jié)果Fig.7 Calculation results when the radiation peak is 200 W·m－2

圖8 輻射峰值為800 W·m－2計(jì)算結(jié)果Fig.9 Calculation results when the radiation peak is 800 W·m－2

由圖2、圖7 和圖8 可見:三種輻射條件下，土坡淺層凍深依次為6 cm、5 cm、5 cm，說明輻射對(duì)土坡淺層凍深影響不大。凍結(jié)時(shí)間上，三種工況下土坡淺層凍結(jié)時(shí)間段依次為:4:00 ～11:00、4:00 ～8:20、4:00 ～7:40，說明雖然輻射對(duì)起凍時(shí)間幾乎無影響，但對(duì)土坡淺層冰保存時(shí)間有較大影響。三種輻射條件下，含冰率峰值均為0.14，且均出現(xiàn)在淺層2 cm 位置，說明輻射對(duì)于土坡淺層凍結(jié)主要影響其凍結(jié)時(shí)間，對(duì)含冰率峰值并沒有太大影響。

3 結(jié) 論

①福建典型短時(shí)凍區(qū)氣候下，土坡淺層最大的凍深約為10 cm 以內(nèi)，淺層一般約在4:00 開始凍結(jié)，含冰率峰值一般出現(xiàn)在表層以下。短時(shí)凍區(qū)土質(zhì)邊坡凍融失穩(wěn)分析時(shí)重點(diǎn)關(guān)注淺層10 cm 左右穩(wěn)定性。

②日最低氣溫對(duì)土坡淺層凍深有較明顯影響，風(fēng)速及太陽輻射對(duì)土坡淺層凍深影響不大。

③日最低氣溫及風(fēng)速對(duì)于土坡淺層凍結(jié)時(shí)間影響較大，二者的增大均會(huì)導(dǎo)致土坡淺層凍結(jié)時(shí)間的延長(zhǎng)，此外，低輻射也能夠促進(jìn)土坡淺層凍結(jié)時(shí)間的延長(zhǎng)。

④日最低氣溫對(duì)于含冰率峰值有較顯著影響;風(fēng)速及太陽輻射對(duì)于土坡淺層含冰率峰值影響不大。

［1］ HARLAN R L.Analysis of coupled heat-fluid transport in partially frozen soil［J］.Water Resoure Reserch，1973，9(5):1314-1323.

［2］ TAYLOR G S，LUTHIN J N.A model for coupled heat and moisture transfer during soil freezing［J］.Canadian Geotechnical Journal，1978，15(4):548-555.

［3］ 毛雪松，李寧，王秉綱，等.多年凍土路基水-熱-力耦合理論模型及數(shù)值模擬［J］.長(zhǎng)安大學(xué)學(xué)報(bào):自然科學(xué)版，2006，26(4):16-20.

［4］ 何敏，李寧，劉乃飛.飽和凍土水熱力耦合模型解析及驗(yàn)證［J］.巖土工程學(xué)報(bào)，2012，34(10):1058-1065.

［5］ 宋存牛.凍融過程中土體水熱力耦合作用理論和模型研究進(jìn)展［J］.冰川凍土，2010，32(5):982-987.

［6］ 鄭秀清，樊貴盛，邢述彥.水分在季節(jié)性非飽和凍融土壤中的運(yùn)動(dòng)［M］.北京:地質(zhì)出版社，2002.

［7］ 李雙洋，張明義，黃志軍，等.地鐵雙線隧道施工人工凍結(jié)水熱力數(shù)值分析［J］.交通運(yùn)輸工程學(xué)報(bào)，2009，9(1):67-73.

［8］ 沈科偉，榮傳新，余文文.查干淖爾礦副井凍結(jié)溫度場(chǎng)有限元分析［J］.廣西大學(xué)學(xué)報(bào):自然科學(xué)版，2011，36(4):647-652.

［9］ 程學(xué)磊，李順群，孫世娟，等.圍壓和負(fù)溫對(duì)凍土強(qiáng)度和微結(jié)構(gòu)的影響研究［J］.廣西大學(xué)學(xué)報(bào):自然科學(xué)版，2014，39(1):95-104.

［10］BRISTOW K L，HORTON R.Modeling the impact of partial surface mulch on soil heat and water flow［J］.Theoretical and Applied Climatology，1996，54(1-2):85-98.

［11］陳祖鑫.福建省花崗巖殘積土路堤溫濕耦合效應(yīng)研究［D］.福州:福州大學(xué)土木工程學(xué)院，2012.

［12］VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5):892-898.

［13］徐斅祖，王家澄，張立新.凍土物理學(xué)［M］.北京:科學(xué)出版社，2010.

［14］陳曉鵬.短時(shí)凍土區(qū)殘積土邊坡水熱耦合機(jī)理的研究［D］.福州:福州大學(xué)土木工程學(xué)院，2014.