不同初始含水率砂土水-汽-熱耦合遷移數(shù)值模擬

郭 毅，晉 華，桂金鵬

(太原理工大學(xué)水利科學(xué)與工程學(xué)院，太原 030024)

干旱半干旱地區(qū)，晝夜溫差較大，溫度波動(dòng)引起的液態(tài)水流動(dòng)和水汽輸送會(huì)對土壤水分的時(shí)空分布產(chǎn)生顯著影響[1,2]。通常情況下，干旱區(qū)近地表土壤表層的含水量很低，水汽通量成為總水分通量的重要組成部分[3-5]。而且，水汽遷移同時(shí)攜帶了大量能量，進(jìn)而影響了土壤內(nèi)的溫度分布。例如，Cahill等研究發(fā)現(xiàn)在裸土頂層2 cm處，40%到60%的熱通量傳遞是由水蒸氣遷移產(chǎn)生的[6]。土壤的水熱狀態(tài)會(huì)對發(fā)生在土壤內(nèi)物理，生物和化學(xué)過程產(chǎn)生重要影響[7]。因此，有必要對干旱區(qū)土壤水-汽-熱耦合遷移規(guī)律展開研究。

Penman[8]首次提出了基于Fick定律的土壤中水蒸氣遷移理論。隨后，Philip和de Vries[9]建立了溫度梯度和含水率梯度雙重作用下的水-汽-熱耦合遷移數(shù)學(xué)模型，即PDV模型，并引入了水汽擴(kuò)散增強(qiáng)因子解釋溫度梯度驅(qū)動(dòng)的水蒸氣通量。Milly[10]改進(jìn)了PDV模型，在模型使用基質(zhì)吸力代替了體積含水率，考慮了土壤的非均質(zhì)性。Nassar等[11]對PDV模型進(jìn)行了擴(kuò)展，提出了土壤水分、熱量和溶質(zhì)耦合作用的傳輸方程。Bittelli等[12]進(jìn)行了裸土條件下熱量，液態(tài)水和水汽傳輸?shù)膶?shí)驗(yàn)研究和數(shù)值模擬。國內(nèi)學(xué)者也在該領(lǐng)域做了大量的研究。例如，蔡樹英[13]通過室內(nèi)土柱蒸發(fā)試驗(yàn)，對PDV模型進(jìn)行了驗(yàn)證，并與等溫模型進(jìn)行了對比。Liu Wei等[14]進(jìn)行了自然環(huán)境下熱量、水分和氣體傳遞的數(shù)值模擬，分析了環(huán)境溫度、相對濕度、太陽輻射和風(fēng)速等參數(shù)對水熱遷移的影響。王華軍等[15]對砂土和壤土進(jìn)行了考慮水蒸氣遷移的數(shù)值模擬研究，發(fā)現(xiàn)兩種土質(zhì)在低飽和度時(shí)，水蒸汽遷移效應(yīng)明顯。

綜上，國內(nèi)外學(xué)者已經(jīng)進(jìn)行了很多水-汽-熱遷移的相關(guān)研究，但不同初始含水率下水蒸氣對土壤水熱遷移影響的定量研究尚未見到，因此，有必要開展該方面的研究。本文以一端有熱源影響的一維土柱為研究對象，利用Hydrus-1D軟件將非飽和土壤水分運(yùn)動(dòng)方程和熱量傳遞方程進(jìn)行了耦合迭代計(jì)算，對不同初始含水率的砂土進(jìn)行了有無水蒸氣影響對比模擬，分析水蒸氣遷移對整體水分遷移和溫度傳遞的影響，以期為進(jìn)一步研究土壤水-汽-熱耦合傳輸規(guī)律提供參考。

1 理論模型

1.1 物理模型

本次物理模型以100 cm長土柱為研究對象，對土柱左端施加恒溫?zé)嵩?，在距熱? cm處設(shè)置第1個(gè)觀測點(diǎn)，然后分別以5、10、10、15、15和20 cm間隔的布置其他6個(gè)觀測點(diǎn)，研究土柱含水率和溫度的變化規(guī)律，物理模型如圖1所示。為了簡化模型，做出如下假設(shè)：①土樣為均質(zhì)和各向同性的多孔介質(zhì)；②不考慮土壤的形變；③只考慮水平方向的水汽熱傳遞；④不考慮熱輻射，只考慮熱傳導(dǎo)和熱對流。

圖1 一維土柱水-汽-熱耦合遷移模型示意圖Fig.1 One-dimensional model of coupled water-vapor-heat migration in soil

1.2 數(shù)值模型

1.2.1 水分運(yùn)動(dòng)方程

非飽和土壤水分傳輸有兩種方式，一種是液態(tài)水的傳輸，一種是氣態(tài)水即水蒸氣的傳輸，驅(qū)動(dòng)力有基質(zhì)勢梯度和溫度梯度，水分遷移方程可用如下方程表述[16]:

(1)

式中：θ是總體積含水率，θ=θl+θv；θl是液態(tài)水體積含水率；θv是水蒸氣體積含水率；KLh是基質(zhì)勢梯度驅(qū)動(dòng)的液態(tài)水等溫水力傳導(dǎo)率；KLT是溫度梯度驅(qū)動(dòng)的液態(tài)水非等溫水力傳導(dǎo)率；Kvh是基質(zhì)勢梯度驅(qū)動(dòng)的水蒸氣等溫水力傳導(dǎo)率；KvT是溫度梯度驅(qū)動(dòng)的水蒸氣非等溫水力傳導(dǎo)率；T是溫度；t是時(shí)間；h是負(fù)壓水頭。x表示距熱源的距離，原點(diǎn)位于土柱最左端，向右為正。

土壤水分特征曲線是描述負(fù)壓水頭與含水率的關(guān)系曲線，對分析和評價(jià)土壤水分的保持和運(yùn)動(dòng)有重要的意義。該特征曲線和液態(tài)水等溫水力傳導(dǎo)率KLh可由VG模型表述[17]。

(2)

m=1-1/n,n>1

(3)

(4)

KLT，Kvh和KvT由下列經(jīng)驗(yàn)公式求得[19-21]：

(6)

(7)

(8)

式中：GWT是表征土壤水分特征曲線的溫度依賴性的參數(shù)，取值7[22]，γ是土壤水的表面張力，γ=75.6-0.1425T-2.38×10-4T2；γ0是25 ℃時(shí)候土壤水表面張力；D是水蒸氣擴(kuò)散度；ρw是液態(tài)水密度；ρvs是飽和水汽密度；Hr是相對濕度；M是水分子的摩爾質(zhì)量；g是重力加速度；R是氣體常數(shù)；η是水蒸氣擴(kuò)散增強(qiáng)因子[23]。

水蒸氣擴(kuò)散度可以定義為:

D=τθaDa

(9)

式中：τ是土體孔隙曲率系數(shù)；Da是自由大氣擴(kuò)散系數(shù)；θa是被空氣所填充的孔隙率。τ和Da可用下列給出的公式求得[24]。

(10)

(11)

ρvs、Hr、η分別由下列經(jīng)驗(yàn)公式求得[24,10,23]：

ρvs=

(12)

(13)

(14)

式中：fc是土壤中黏土的質(zhì)量分?jǐn)?shù)。

1.2.2 熱量傳遞方程

非飽和土壤中熱量傳輸分為固體骨架熱傳導(dǎo)、液態(tài)水對流換熱、水蒸氣對流換熱和液態(tài)水的蒸發(fā)潛熱四部分，基于能量守恒定律，熱量傳輸可用如下公式表述[16]：

(15)

土壤表觀比熱容和表觀熱導(dǎo)率可以下列經(jīng)驗(yàn)公式求得[25,26]：

CP(θ)=Cnθn+Cwθw+Cvθv

(16)

λ(θ)=λ0(θ)+βtcw|q|

(17)

式中：θ是體積分?jǐn)?shù)；下標(biāo)n,w,v分別代表固液氣三相；βt為熱彌散系數(shù)；λ0(θ)為多孔介質(zhì)熱導(dǎo)率，可由公式λ0(θ)=b1+b2θ+b3θ0.5[27]求得。

1.2.3 初始條件與邊界條件

初始條件：在初始狀態(tài)下，土柱內(nèi)各處溫度、含水率相等，即：

T(x,t)=T0t=0

(18)

θ(x,t)=θ0t=0

(19)

邊界條件：在x=0處

T(x,t)=TH

(20)

(21)

式中：TH為熱源溫度。

在x=L處

T(x,t)=Tc

(22)

(23)

式中：Tc為環(huán)境溫度。

式(1)、式(15)與上述初始條件和邊界條件構(gòu)成非飽和土壤一維水平水-汽-熱耦合遷移數(shù)學(xué)模型。

1.3 數(shù)值計(jì)算

非飽和土壤水-汽-熱遷移方程組具有非線性和耦合性，無法求得解析值，故利用Hydrus-1D軟件對水熱傳遞方程進(jìn)行離散化，進(jìn)而求得數(shù)值解。模型求解區(qū)域?yàn)?00 cm長土柱，空間計(jì)算步長為1 cm，模擬時(shí)間為1 440 min，初始計(jì)算時(shí)間步長為1 min，最小計(jì)算時(shí)間步長為0.01 min，最大計(jì)算時(shí)間步長為5 min。

1.4 數(shù)值模擬方案

基于非飽和土壤水-汽-熱耦合遷移數(shù)值模型，設(shè)置以下6種模擬情景，模擬研究不同初始含水率下水蒸氣對土壤整體水分遷移和溫度變化的影響，模擬方案如表1所示。為了保證模擬結(jié)果適用的廣泛性，所選土質(zhì)為Hydrus軟件內(nèi)置的砂土，砂土VG模型參數(shù)如表2所示。

表1 模擬方案表Tab.1 Simulation scheme table

表2 砂土VG模型參數(shù)Tab.2 Parameters of VG model for sand

2 模擬結(jié)果與分析

考慮到砂土含水率和溫度在距熱源較近測點(diǎn)變化劇烈，在離熱源較遠(yuǎn)的測點(diǎn)變化較小，因此本文選擇在3和8 cm測點(diǎn)的溫濕度動(dòng)態(tài)變化進(jìn)行研究。同時(shí)，由于含水率的絕對變化量相對較小，為方便對比水蒸汽在非飽和土壤水汽熱遷移中的貢獻(xiàn)，對含水率和溫度變化均進(jìn)行了無量綱處理，具體公式如下[16]：

無量綱土壤含水率：

(26)

無量綱土壤溫度：

(27)

2.1 含水率變化模擬結(jié)果與分析

圖2是初始含水率分別為0.08、0.11、0.13、0.15、0.20和0.30 cm3/cm3時(shí)砂土含水率隨時(shí)間變化曲線?？梢钥吹?，無論是否考慮水蒸氣遷移，在距熱源3和8 cm處均會(huì)形成一個(gè)短暫的含水率峰值，該峰值現(xiàn)象的出現(xiàn)與陳紅兵、王華軍等[28,29]人的研究結(jié)果一致，這是因?yàn)樵跍囟忍荻群秃侍荻入p重作用下，水分產(chǎn)生了遷移，在該點(diǎn)處形成了短暫的水分滯留；在考慮水蒸氣遷移后，含水率變化更加劇烈，表現(xiàn)為含水率峰值與模擬結(jié)束時(shí)刻土壤最終含水率之間的差值變大。下面以模擬結(jié)束時(shí)刻距熱源3 cm處的最終含水率為例，分析不同初始含水率下水蒸氣遷移對整體水分遷移的影響。當(dāng)砂土初始含水率為0.08 cm3/cm3時(shí)，僅考慮液態(tài)水遷移，含水率減少0.039，考慮水蒸氣遷移，含水率減少0.154，含水率的減少量增幅為294.87%；當(dāng)砂土初始含水率為0.11 cm3/cm3，僅考慮液態(tài)水遷移，含水率減少0.091，考慮水蒸氣遷移后，含水率減少0.130，含水率減少量增幅為42.86%；當(dāng)砂土初始含水率為0.13 cm3/cm3時(shí)，僅考慮液態(tài)水遷移，含水率減少0.133，考慮水蒸氣遷移，含水率減少0.148，含水率減少量增幅為11.28%；當(dāng)砂土初始含水率為0.15 cm3/cm3，僅考慮液態(tài)水遷移，測點(diǎn)含水率減少0.159，考慮水蒸氣遷移，含水率減少0.169，含水率減少量增幅為6.29%；當(dāng)砂土初始含水率為0.20 cm3/cm3，僅考慮液態(tài)水遷移，含水率減少0.208，考慮水蒸氣遷移后，含水率減少0.209，含水率減少量增幅為0.48%；當(dāng)砂土初始含水率為0.30 cm3/cm3，僅考慮液態(tài)水遷移，含水率減少0.222，考慮水蒸氣遷移，含水率減少0.222，含水率減少量增幅為0%。進(jìn)一步地，圖3是不同初始含水率下，考慮水蒸氣遷移后3 cm處含水率減少量的增幅圖?？梢钥吹剑撼跏己试降停魵鈱ν寥雷罱K含水率的影響越大，初始含水率為0.08 cm3/cm3時(shí)，考慮水蒸氣后，砂土含水率減少量增幅極大，這表明水蒸氣遷移在整體水分遷移中占據(jù)了主導(dǎo)作用，隨初始含水率的增加到0.11 cm3/cm3，砂土含水率減少量增幅迅速減小，當(dāng)初始含水率增加到0.20 cm3/cm3時(shí)，水蒸氣遷移對整體水分遷移的影響可以忽略。

圖2 不同初始含水率條件下砂土含水率隨時(shí)間變化曲線Fig.2 The moisture variation curve of sand with time under different initial moisture content

圖3 砂土不同初始含水率下考慮水蒸氣遷移后含水率減少量增幅圖Fig.3 The increase of water content reduction of sand after considering water vapor migration under different initial water content

2.2 溫度變化模擬結(jié)果與分析

圖4是熱源溫度為50 ℃，初始含水率分別為0.08、0.11、0.13、0.15、0.20和0.30 cm3/cm3時(shí)溫度隨時(shí)間變化曲線。對比有無水蒸氣對溫度模擬結(jié)果的影響，這里以模擬120 min后距熱源3 cm處的溫度變化為例進(jìn)行分析：當(dāng)初始含水率為0.08 cm3/cm3時(shí)，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，模擬點(diǎn)處溫度從0.800增加到0.828，增幅為3.5%；當(dāng)初始含水率為0.11 cm3/cm3，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，模擬點(diǎn)處溫度從0.801增加到0.823，增幅為2.75%；初始含水率為0.13 cm3/cm3，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，模擬點(diǎn)處溫度從0.799增加到0.819，增幅為2.50%；初始含水率為0.15 cm3/cm3，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，溫度從0.798上升到0.816，增幅為2.26%；含水率為0.20 cm3/cm3，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，溫度從0.796上升到0.808，增幅為1.51%；含水率為0.30 cm3/cm3，與僅有液態(tài)水遷移相比，考慮水蒸氣遷移后，模擬點(diǎn)處溫度從0.792上升到0.796，增幅為0.51%。進(jìn)一步地，圖5是不同初始含水率下，考慮水蒸氣遷移后3 cm處砂土在120 min時(shí)刻溫度增幅圖?？梢钥吹剑核魵獾倪w移加快了砂土溫度上升的速度，且隨著初始含水率增加，水蒸氣對溫度提升幅度不斷減小。砂土中水分遷移和熱量傳遞是個(gè)耦合的過程，砂土溫度變化的趨勢與含水率變化是密切相關(guān)，當(dāng)初始含水率較少時(shí)，土體孔隙的氣相成分所占比例較大，水蒸氣的遷移較其在高含水率中的遷移阻力較小，水蒸氣遷移對土壤含水率影響很大，故砂土溫度上升過程差異較大。隨著含水率增加，水蒸氣遷移對整體水分遷移貢獻(xiàn)越來越小，考慮和不考慮水蒸水的模擬差別越來越小。當(dāng)含水率增加到0.20 cm3/cm3時(shí)，模擬點(diǎn)處砂土溫度的變化幾乎一致，可以忽略水蒸氣遷移對溫度變化的影響。

圖4 不同初始含水率條件下砂土溫度隨時(shí)間變化曲線Fig.4 The temperature variation curve of sand with time under different initial moisture content

圖5 砂土不同初始含水率考慮水蒸氣遷移后溫度增幅圖Fig.5 The temperature increase diagram of sand with different initial moisture content after considering water vapor migration

3 結(jié) 論

本文基于Hydrus-1D軟件構(gòu)建了熱源作用下一維非飽和土壤水-汽-熱耦合遷移模型，進(jìn)行了多種初始含水率下有無水蒸氣的對比模擬，可以得到以下結(jié)論。

(1)與僅考慮液態(tài)水遷移相比，考慮水蒸氣遷移后會(huì)導(dǎo)致含水率遷移加劇，當(dāng)初始含水率為0.08 cm3/cm3時(shí)，水蒸氣遷移在整體水分遷移中占主導(dǎo)作用，隨著初始含水率增加0.11 cm3/cm3，水蒸氣遷移對整體水分遷移影響大幅減小，當(dāng)初始含水率繼續(xù)增加到0.20 cm3/cm3，水蒸氣對整體水分遷移的效果可忽略。

(2)與僅考慮液態(tài)水遷移相比，考慮水蒸氣遷移后會(huì)導(dǎo)致土壤溫度升高速度提升，且初始含水率越低，水蒸氣對溫度的提升幅度越大，隨著初始含水率不斷增加，水蒸氣對土壤溫度的影響不斷減小直至消失。

(3)在對低含水率砂土進(jìn)行水熱遷移模擬時(shí)，有必要將水蒸氣因素考慮進(jìn)去提高模擬精度，而對于含水率較高的砂土，可僅考慮液態(tài)水遷移以簡化計(jì)算。