日光溫室內(nèi)淺層土壤溫濕度場(chǎng)的數(shù)值模擬

張 艷，杜震宇

(太原理工大學(xué) 環(huán)境科學(xué)與工程學(xué)院，太原 030024)

近年來(lái)，日光溫室在我國(guó)農(nóng)業(yè)發(fā)展中取得了顯著成就，已成為我國(guó)農(nóng)業(yè)的主要生產(chǎn)形式，但日光溫室內(nèi)熱濕環(huán)境比較復(fù)雜，室外氣象條件、圍護(hù)結(jié)構(gòu)、土壤物性參數(shù)、作物的生長(zhǎng)等都會(huì)影響室內(nèi)溫濕度的調(diào)控，從而影響室內(nèi)作物的產(chǎn)量[1-3]。由于土壤的蓄放熱對(duì)溫室的熱濕環(huán)境有著決定性的影響[4]，因此有必要對(duì)土壤溫濕度場(chǎng)進(jìn)行研究，為研發(fā)低成本、運(yùn)行費(fèi)用低、節(jié)能的溫濕度調(diào)控裝置奠定基礎(chǔ)，以提高日光溫室的熱濕環(huán)境調(diào)控技術(shù)。國(guó)內(nèi)外學(xué)者對(duì)于日光溫室熱濕環(huán)境已做了部分研究。SHARMA、TIWARI et al[5-6]等利用Runge-Kutta方法分析溫室室內(nèi)空氣熱濕交換規(guī)律，建立了考慮導(dǎo)熱、對(duì)流等作用的能量平衡方程，但未給出土壤區(qū)域的相關(guān)方程。FAYER[7]提出了土壤熱濕傳遞的理論模型，但忽略了土壤溫度梯度對(duì)水分傳遞的影響。楊艷超[8]對(duì)山東省日光溫室微氣候條件的模擬研究表明日光溫室內(nèi)潛熱、顯熱的變化趨勢(shì)與太陽(yáng)總輻射相一致。劉宏[9]對(duì)日光溫室土壤-空氣換熱器的換熱特性進(jìn)行了試驗(yàn)研究，測(cè)試了不同空氣流速換熱管內(nèi)空氣溫濕度和換熱管周圍的土壤溫度，得出了最佳入口空氣流速。李鈺楠[10]使用Visual Basic6.0編寫了土壤熱濕耦合迭代求解軟件，探究了土壤物性等初始計(jì)算參數(shù)對(duì)土壤-空氣換熱器埋地?fù)Q熱管周圍土壤溫濕度場(chǎng)的影響。范毅等[11]對(duì)日光溫室環(huán)境下土壤空氣換熱器的換熱特性進(jìn)行了研究，分析了試驗(yàn)工況下土壤空氣換熱器的動(dòng)態(tài)換熱過(guò)程及系統(tǒng)性能變化規(guī)律。YENER et al[12]對(duì)自然條件下土耳其不同區(qū)域不同深度土壤溫度進(jìn)行理論與試驗(yàn)研究，分析比較了土壤溫度與空氣溫度的變化規(guī)律。MAHDAVI et al[13]分析比較了自然條件下有秸稈覆蓋與裸土淺層土壤蒸發(fā)量的變化。

由此可見(jiàn)，國(guó)內(nèi)外學(xué)者做了大量有關(guān)溫室室內(nèi)溫濕度、熱環(huán)境以及土壤-空氣換熱器埋管周圍土壤溫濕度場(chǎng)的研究，對(duì)于淺層土壤的研究多為裸露在室外的土壤，且忽略了濕度場(chǎng)對(duì)溫度場(chǎng)的影響，對(duì)日光溫室內(nèi)淺層土壤溫濕度場(chǎng)的理論與試驗(yàn)研究都很少。本文利用Fluent15.0對(duì)日光溫室淺層土壤(注：≤100 cm深度)溫濕度場(chǎng)進(jìn)行三維非穩(wěn)態(tài)數(shù)值模擬，在此基礎(chǔ)上利用該模型研究了大尺度內(nèi)淺層土壤溫濕度場(chǎng)的變化規(guī)律，旨在為研發(fā)日光溫室熱濕環(huán)境調(diào)控裝置提供科學(xué)依據(jù)。

1 試驗(yàn)內(nèi)容

1.1 供試日光溫室

該試驗(yàn)所用日光溫室位于太原市小店區(qū)，坐北朝南，方位角為南偏西5 °，南北跨度9.7 m，東西凈長(zhǎng)50 m，后墻高度3.1 m.南面采光面為塑料薄膜，在膜的最底部和頂部有可啟閉的通風(fēng)口，可以將熱濕空氣排出[14]。溫室結(jié)構(gòu)尺寸如圖1所示。

建立一個(gè)三維直角坐標(biāo)系，見(jiàn)圖1.以水平向西為Z軸正方向，以水平向南為X軸正方向，以垂直于地表指向天空為Y軸正方向，以溫室地表面東西向與南北向中部位置的交點(diǎn)為坐標(biāo)原點(diǎn)。

1.2 測(cè)點(diǎn)布置與編號(hào)

供試溫室沒(méi)有種植作物，地面為裸地。試驗(yàn)測(cè)點(diǎn)布置在溫室中部的位置，如圖2所示。為方便討論，將試驗(yàn)所布置的三個(gè)測(cè)點(diǎn)從北到南依次命名為測(cè)點(diǎn)1、測(cè)點(diǎn)2、測(cè)點(diǎn)3.在圖2所示測(cè)點(diǎn)位置距地表30 cm范圍內(nèi)的土壤層內(nèi)埋設(shè)土壤溫濕度傳感器，測(cè)點(diǎn)的具體位置及編號(hào)如圖3所示[15]。

圖2 南北向測(cè)點(diǎn)位置及編號(hào)(單位：mm)Fig.2 North-south position and number of measuring points(unit:mm)

圖3 土壤層(≤30 cm)測(cè)點(diǎn)位置及編號(hào)(單位：mm)Fig.3 Soil layer (≤30 cm) measuring points location and number(unit:mm)

1.3 數(shù)據(jù)采集與整理

試驗(yàn)數(shù)據(jù)的采集時(shí)間為2017年3月11日—3月25日。通過(guò)傳感器采集土壤溫濕度，自動(dòng)儲(chǔ)存至數(shù)據(jù)采集儀。存儲(chǔ)的數(shù)據(jù)可通過(guò)USB端口導(dǎo)入計(jì)算機(jī)內(nèi)，進(jìn)行數(shù)據(jù)整理與分析。

2 理論模型

2.1 物理模型

本模型中沒(méi)有考慮室內(nèi)作物，所以認(rèn)為地面為裸地。地面的熱濕過(guò)程包括：土壤垂直方向的熱量傳遞，與室內(nèi)空氣的熱對(duì)流交換，與室內(nèi)各個(gè)面的長(zhǎng)波輻射熱交換，吸收的太陽(yáng)直射輻射、散射輻射，土壤表面水汽相變傳遞的熱濕量等熱濕交換過(guò)程。圖4為日光溫室地表能量收支示意圖。

圖4 日光溫室地表能量收支示意圖Fig.4 Sketch of the surface heat balance of sunlight greenhouse

溫室內(nèi)地表面能量平衡模型方程如下：

Q+I+H+L-γW=0 .

(1)

式中：Q為熱傳導(dǎo)熱量，J；I為太陽(yáng)輻射量，J；H為對(duì)流輻射換熱量，J；L為長(zhǎng)波輻射換熱量，J；γ為汽化潛熱，J/kg；W為水分傳輸量，kg.

本文對(duì)實(shí)際模型進(jìn)行簡(jiǎn)化和假設(shè)，建立了淺層土壤溫濕度場(chǎng)的三維非穩(wěn)態(tài)數(shù)值模型。采用CFD數(shù)值計(jì)算方法，考慮了多變的太陽(yáng)輻射和室外環(huán)境條件，用UDS定義土壤含水率，編寫濕方程的UDF程序，導(dǎo)入Fluent中對(duì)日光溫室淺層土壤溫濕度場(chǎng)進(jìn)行數(shù)值求解。

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

模型簡(jiǎn)化假設(shè)：

1) 采用Boussinesq假設(shè)，以便于處理由于溫差而引起的浮升力項(xiàng)；

2) 將溫室圍護(hù)結(jié)構(gòu)材料視為常物性的固體材料，土壤視為多孔介質(zhì)；

3) 不考慮土壤中水分的宏觀流動(dòng)，忽略重力對(duì)土壤內(nèi)水分的影響；

4) 土壤中的水蒸氣與溫室內(nèi)的濕空氣遵循理想氣體狀態(tài)方程。

空氣區(qū)域控制方程包括如下方程。

連續(xù)性方程：

(2)

動(dòng)量方程：

(3)

(4)

(5)

(6)

能量方程：

(7)

k方程：

(8)

ε方程：

(9)

土壤區(qū)域的濕遷移方程：

(10)

(11)

參考文獻(xiàn)[16]給出的DT、Dθ的擬合公式：

DT=-2.349 5×10-5·θ6+2.037×10-5·θ5- 6.694 6×10-6·θ4+1.043 8×10-6·θ3- 7.869 9×10-8·θ2+2.735+10-9·θ+ 5.403 9×10-12.

(12)

Dθ=2.217 6×10-6(θ6·e-7) exp(-2.09×10-3·T) .

(13)

其中，e為土壤孔隙率，本文取砂土孔隙率0.37.

2.3 幾何模型

利用ICEM CFD 15.0軟件建立日光溫室和土壤的幾何模型，土壤厚度100 cm.鑒于溫室結(jié)構(gòu)的不規(guī)整性，結(jié)構(gòu)化網(wǎng)格劃分比較困難，計(jì)算結(jié)果不準(zhǔn)確，而非結(jié)構(gòu)網(wǎng)格能對(duì)不規(guī)整模型進(jìn)行劃分，生成高質(zhì)量的網(wǎng)格，所以選擇使用ICEM CFD非結(jié)構(gòu)網(wǎng)格劃分。計(jì)算區(qū)域網(wǎng)格劃分如圖5所示。

圖5 模型網(wǎng)格劃分Fig.5 Model mesh generation

2.4 網(wǎng)格獨(dú)立性考核

數(shù)值模擬結(jié)果的準(zhǔn)確性與網(wǎng)格數(shù)目有很大關(guān)系，網(wǎng)格數(shù)目過(guò)多會(huì)導(dǎo)致計(jì)算機(jī)計(jì)算時(shí)間過(guò)長(zhǎng)，數(shù)目太少會(huì)影響計(jì)算結(jié)果的準(zhǔn)確性，因此采用數(shù)值模擬計(jì)算進(jìn)行網(wǎng)格獨(dú)立性考核，本次模擬共選取了5套網(wǎng)格，網(wǎng)格數(shù)量分別為：515 720、1 032 885、1 456 720、1 852 662和2 216 874.如表1所示，選取同一時(shí)刻測(cè)點(diǎn)1的溫度值進(jìn)行獨(dú)立性考核，網(wǎng)格數(shù)1 852 662測(cè)點(diǎn)1溫度值為17.9 ℃，網(wǎng)格數(shù)2 216 874測(cè)點(diǎn)1溫度值為17.8 ℃，相對(duì)誤差小于1%，考慮到計(jì)算時(shí)長(zhǎng)以及模擬結(jié)果的準(zhǔn)確性，本文選取網(wǎng)格數(shù)量1 852 662來(lái)進(jìn)行數(shù)值模擬。

表1 網(wǎng)格獨(dú)立性考核Table 1 Mesh independence test

2.5 數(shù)值計(jì)算方法

控制方程的求解采用有限容積法，采用SIMPLE算法處理壓力與速度的耦合關(guān)系，對(duì)流項(xiàng)的離散格式采用三階精度的QUICK格式。數(shù)值模擬采用RNGk-ε湍流模型，非平衡壁面函數(shù)，DO輻射模型。用UDS定義土壤含水率，編寫濕方程的UDF程序，導(dǎo)入Fluent中進(jìn)行數(shù)值求解。模擬計(jì)算從3月11日18時(shí)00分起，時(shí)間步長(zhǎng)300 s，連續(xù)計(jì)算2 d.

2.6 邊界條件

模型的計(jì)算域?yàn)槿展鉁厥覂?nèi)的空氣區(qū)域與深度為100 cm的土壤區(qū)域，地表面既與溫室內(nèi)空氣進(jìn)行傳熱傳質(zhì)，也與土壤層存在溫度與水分的傳遞，因此將地表面設(shè)為流固耦合面，溫室圍護(hù)結(jié)構(gòu)和塑料薄膜設(shè)為第三類邊界條件，材料的物性參數(shù)如表2所示，將來(lái)風(fēng)面設(shè)為velocity-inlet邊界條件，風(fēng)速為2.5 m/s，將實(shí)測(cè)的室外空氣溫度回歸擬合出空氣溫度隨時(shí)間變化的函數(shù)，將公式編入U(xiǎn)DF，設(shè)為入口溫度邊界條件，出口設(shè)為outflow.將土壤的原始溫度場(chǎng)編入U(xiǎn)DF，設(shè)為土壤原始溫度，作為土壤的邊界條件。將土壤深100 cm處的含濕量測(cè)量值設(shè)為土壤的初始體積含水率，即0.2 m3/m3.

圍護(hù)結(jié)構(gòu)與室外空氣的對(duì)流換熱系數(shù)與室外風(fēng)速有關(guān)，具體計(jì)算式為[17]：

h=18.63v0.603.

(14)

式中：v為圍護(hù)結(jié)構(gòu)外表面的實(shí)際風(fēng)速，其值大小與大氣風(fēng)速v有關(guān)。

對(duì)于背風(fēng)面：v=0.3+0.05v.

對(duì)于迎風(fēng)面：

當(dāng)v≥2 m/s時(shí)，v=0.25v.

當(dāng)v<2 m/s時(shí)，v=0.5.

太原地區(qū)的主導(dǎo)風(fēng)向?yàn)槲鞅憋L(fēng)，平均風(fēng)速2.5 m/s，圍護(hù)結(jié)構(gòu)的對(duì)流換熱系數(shù)以及邊界條件的設(shè)置見(jiàn)表3.

土壤的原始溫度場(chǎng)

(15)

土壤的原始溫度場(chǎng)由崔良衛(wèi)等[18]計(jì)算得出，Ω為溫度波波動(dòng)頻率，7.17×10-4rad/h；a為導(dǎo)溫系數(shù)，m2/h；td為地表年平均溫度，試驗(yàn)地取11.3 ℃；Ad為地表溫度波振幅，試驗(yàn)地取17.75 ℃；y為地層深度，m；t0為土壤原始溫度，℃.

表2 材料的物性參數(shù)Table 2 Material properties

表3 邊界條件Table 3 Boundary conditions

3 結(jié)果與討論

3.1 模擬驗(yàn)證

測(cè)點(diǎn)1、測(cè)點(diǎn)2、測(cè)點(diǎn)3土壤溫度的模擬值與試驗(yàn)值的對(duì)比結(jié)果如圖6所示，測(cè)點(diǎn)1、測(cè)點(diǎn)2、測(cè)點(diǎn)3土壤含水率的模擬值與試驗(yàn)值的對(duì)比結(jié)果如圖7所示，測(cè)點(diǎn)1(Y=0)、測(cè)點(diǎn)2(Y=0)、測(cè)點(diǎn)3(Y=0)模擬值與試驗(yàn)值吻合較好，溫度平均相對(duì)誤差小于3.8%，含水率平均相對(duì)誤差小于2.6%，在工程測(cè)量的誤差允許范圍內(nèi)。故模型具有一定的合理性，數(shù)值模擬方法具有可靠性。

圖6 溫度場(chǎng)驗(yàn)證Fig.6 Temperature field verification

圖7 濕度場(chǎng)驗(yàn)證Fig.7 Moisture field verification

3.2 土壤溫度場(chǎng)與濕度場(chǎng)的分布

圖8為數(shù)值模擬出不同深度土壤層溫度隨時(shí)間的動(dòng)態(tài)變化圖，從圖中可以看出，不同深度的淺層土壤(注：≤100 cm深度)溫度變化趨勢(shì)一致，都隨時(shí)間呈現(xiàn)周期性的變化。隨著土壤深度的增加，溫度波的振幅逐漸減小，延遲時(shí)間逐漸增大。當(dāng)?shù)乇頊囟雀哂谕寥缹訙囟葧r(shí)，熱量從地表向地下傳遞，土壤層溫度高于地表溫度時(shí)，熱量從地下向地表傳遞。不同深度的土壤地溫變化幅度均表現(xiàn)出5 cm>10 cm>15 cm>20 cm>30 cm>40 cm>50 cm>100 cm的變化規(guī)律，即越靠近地表土壤層溫度日較差越大，0～30 cm土壤層溫度在一天內(nèi)的波動(dòng)振幅為19.7～7.2 ℃，土壤層溫度受氣象條件和環(huán)境變化影響較大，土壤溫度波動(dòng)顯著；30～50 cm土壤層溫度在一天內(nèi)的波動(dòng)振幅為3.6～0.4 ℃，土壤層溫度受氣象條件和環(huán)境變化影響較小，土壤溫度波動(dòng)不太明顯；50～100 cm土壤層溫度在一天內(nèi)的波動(dòng)振幅小于0.5 ℃，土壤層溫度處于比較穩(wěn)定的狀態(tài)，基本不受氣象條件和環(huán)境變化的影響。根據(jù)不同深度淺層土壤溫度梯度的變化規(guī)律，可以把0～100 cm的土層分為三個(gè)特征層：多變層(0～30 cm)，緩變層(30～50 cm)，均穩(wěn)層(50～100 cm).隨著土壤深度增加，最高氣溫與最低氣溫出現(xiàn)的時(shí)間總體上也相對(duì)滯后。

圖8 土壤溫度場(chǎng)(≤100 cm)Fig.8 Soil temperature field (≤100 cm)

淺層土壤中不僅存在水蒸汽的冷凝蒸發(fā)，還存在多孔介質(zhì)傳熱傳質(zhì)的Soret效應(yīng)和Dufour效應(yīng)[19]，重力、達(dá)西阻力、土壤溫度以及濕分梯度引起的毛細(xì)力的共同作用，圖9為數(shù)值模擬出不同深度土壤層的含水率隨時(shí)間的動(dòng)態(tài)變化圖，分析圖可知，不同深度的土壤濕度變化趨勢(shì)一致，都隨時(shí)間呈現(xiàn)周期性的變化。隨著土壤深度的增加，含水率的變化幅度逐漸減小，延遲時(shí)間逐漸增大。到50 cm深度的土壤層含水率基本不隨環(huán)境溫濕度的變化而變

圖9 土壤濕度場(chǎng)(≤100 cm)Fig.9 Soil moisture field (≤100 cm)

化，趨于恒定。當(dāng)溫室內(nèi)的空氣溫度高于地表面土壤溫度時(shí)，地表土壤中的水分吸收空氣熱量蒸發(fā)，地表面溫度低于空氣的露點(diǎn)溫度時(shí)，水蒸汽遇冷凝結(jié)成水滴，水蒸汽冷凝放熱，地表面吸熱。

3.3 土壤不同初始體積含水率下溫度場(chǎng)的變化

圖10為數(shù)值模擬出的土壤不同初始含水率條件下淺層土壤(注：≤100 cm深度)溫度場(chǎng)的變化。取0:00測(cè)點(diǎn)1的結(jié)果圖，可以看出，在其他初始條件相同的情況下，不同土壤初始體積含水率下土壤溫度隨深度的變化規(guī)律一致。土壤初始體積含水率越低，淺層土壤層的溫度梯度越大。當(dāng)土壤初始體積含水率為0 m3/m3時(shí)，土壤層溫度梯度最大；當(dāng)土壤初始體積含水率為0.3 m3/m3時(shí)，土壤層溫度梯度最小。這是因?yàn)椋涸谕寥莱跏俭w積含水率較高的情況下，土壤熱容量較大，溫度傳播得較慢，而在土壤初始體積含水率較低的情況下，土壤熱容量較小，溫度傳播得較快。

圖10 不同土壤初始體積含水率下土壤溫度變化Fig.10 Variations of soil temperature under different soil volume moisture

圖11為取12:00測(cè)點(diǎn)1的結(jié)果圖，可以看出，在其他初始條件相同的情況下，土壤初始體積含水率越低，沿Y方向的土壤層溫降梯度越大。當(dāng)土壤體積含水率為0 m3/m3時(shí)，沿Y方向的土壤層溫降梯度最大；當(dāng)土壤初始體積含水率為0.3 m3/m3時(shí)，沿Y方向的土壤層溫降梯度最小。同樣驗(yàn)證了：土壤初始體積含水率越低，土壤熱容量越小，溫度傳播得越快，而在土壤初始體積含水率較高的情況下，土壤熱容量較大，溫度傳播得較慢。

圖11 不同土壤初始體積含水率下土壤溫度變化Fig.11 Variations of soil temperature under different soil volume moisture

4 結(jié)論

本文建立了淺層土壤溫濕度場(chǎng)的三維非穩(wěn)態(tài)數(shù)值模型，對(duì)日光溫室淺層土壤(注：≤100 cm深度)溫濕度場(chǎng)以及土壤不同初始體積含水率下溫度場(chǎng)的變化進(jìn)行了模擬，得出以下結(jié)論：

1) 日光溫室內(nèi)淺層土壤不同深度的溫濕度變化趨勢(shì)一致，都隨時(shí)間呈現(xiàn)周期性的變化，隨著土壤深度的增加，溫度波的振幅逐漸減小，含水率的變化幅度也逐漸減小，延遲時(shí)間逐漸增大， 最高氣溫與最低氣溫出現(xiàn)的時(shí)間總體上也相對(duì)滯后，50 cm深度土壤層溫濕度基本保持恒定。

2) 土壤初始體積含水率越低，沿Y方向的土壤層溫降梯度越大，溫度傳播得越快；土壤初始體積含水率越高，沿Y方向的土壤層溫降梯度越小，溫度傳播得越慢。

3) 根據(jù)不同深度淺層土壤溫度梯度的變化規(guī)律，可以把0～100 cm的土層分為三個(gè)特征層：多變層(0～30 cm)，緩變層(30～50 cm)，均穩(wěn)層(50～100 cm).