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

郭 毅，晉 華，桂金鵬

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

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

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

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

1 理論模型

1.1 物理模型

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

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

1.2 數值模型

1.2.1 水分運動方程

非飽和土壤水分傳輸有兩種方式，一種是液態水的傳輸，一種是氣態水即水蒸氣的傳輸，驅動力有基質勢梯度和溫度梯度，水分遷移方程可用如下方程表述[16]:

(1)

式中：θ是總體積含水率，θ=θl+θv；θl是液態水體積含水率；θv是水蒸氣體積含水率；KLh是基質勢梯度驅動的液態水等溫水力傳導率；KLT是溫度梯度驅動的液態水非等溫水力傳導率；Kvh是基質勢梯度驅動的水蒸氣等溫水力傳導率；KvT是溫度梯度驅動的水蒸氣非等溫水力傳導率；T是溫度；t是時間；h是負壓水頭。x表示距熱源的距離，原點位于土柱最左端，向右為正。

土壤水分特征曲線是描述負壓水頭與含水率的關系曲線，對分析和評價土壤水分的保持和運動有重要的意義。該特征曲線和液態水等溫水力傳導率KLh可由VG模型表述[17]。

(2)

m=1-1/n,n>1

(3)

(4)

KLT，Kvh和KvT由下列經驗公式求得[19-21]：

(6)

(7)

(8)

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

水蒸氣擴散度可以定義為:

D=τθaDa

(9)

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

(10)

(11)

ρvs、Hr、η分別由下列經驗公式求得[24,10,23]：

ρvs=

(12)

(13)

(14)

式中：fc是土壤中黏土的質量分數。

1.2.2 熱量傳遞方程

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

(15)

土壤表觀比熱容和表觀熱導率可以下列經驗公式求得[25,26]：

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

(16)

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

(17)

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

1.2.3 初始條件與邊界條件

初始條件：在初始狀態下，土柱內各處溫度、含水率相等，即：

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為環境溫度。

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

1.3 數值計算

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

1.4 數值模擬方案

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

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

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

2 模擬結果與分析

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

無量綱土壤含水率：

(26)

無量綱土壤溫度：

(27)

2.1 含水率變化模擬結果與分析

圖2是初始含水率分別為0.08、0.11、0.13、0.15、0.20和0.30 cm3/cm3時砂土含水率隨時間變化曲線。可以看到，無論是否考慮水蒸氣遷移，在距熱源3和8 cm處均會形成一個短暫的含水率峰值，該峰值現象的出現與陳紅兵、王華軍等[28,29]人的研究結果一致，這是因為在溫度梯度和含水率梯度雙重作用下，水分產生了遷移，在該點處形成了短暫的水分滯留；在考慮水蒸氣遷移后，含水率變化更加劇烈，表現為含水率峰值與模擬結束時刻土壤最終含水率之間的差值變大。下面以模擬結束時刻距熱源3 cm處的最終含水率為例，分析不同初始含水率下水蒸氣遷移對整體水分遷移的影響。當砂土初始含水率為0.08 cm3/cm3時，僅考慮液態水遷移，含水率減少0.039，考慮水蒸氣遷移，含水率減少0.154，含水率的減少量增幅為294.87%；當砂土初始含水率為0.11 cm3/cm3，僅考慮液態水遷移，含水率減少0.091，考慮水蒸氣遷移后，含水率減少0.130，含水率減少量增幅為42.86%；當砂土初始含水率為0.13 cm3/cm3時，僅考慮液態水遷移，含水率減少0.133，考慮水蒸氣遷移，含水率減少0.148，含水率減少量增幅為11.28%；當砂土初始含水率為0.15 cm3/cm3，僅考慮液態水遷移，測點含水率減少0.159，考慮水蒸氣遷移，含水率減少0.169，含水率減少量增幅為6.29%；當砂土初始含水率為0.20 cm3/cm3，僅考慮液態水遷移，含水率減少0.208，考慮水蒸氣遷移后，含水率減少0.209，含水率減少量增幅為0.48%；當砂土初始含水率為0.30 cm3/cm3，僅考慮液態水遷移，含水率減少0.222，考慮水蒸氣遷移，含水率減少0.222，含水率減少量增幅為0%。進一步地，圖3是不同初始含水率下，考慮水蒸氣遷移后3 cm處含水率減少量的增幅圖。可以看到：初始含水率越低，水蒸氣對土壤最終含水率的影響越大，初始含水率為0.08 cm3/cm3時，考慮水蒸氣后，砂土含水率減少量增幅極大，這表明水蒸氣遷移在整體水分遷移中占據了主導作用，隨初始含水率的增加到0.11 cm3/cm3，砂土含水率減少量增幅迅速減小，當初始含水率增加到0.20 cm3/cm3時，水蒸氣遷移對整體水分遷移的影響可以忽略。

圖2 不同初始含水率條件下砂土含水率隨時間變化曲線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 溫度變化模擬結果與分析

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

圖4 不同初始含水率條件下砂土溫度隨時間變化曲線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 結 論

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

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

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

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