張 艷，杜震宇

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

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

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

1 試驗內(nèi)容

1.1 供試日光溫室

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

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

1.2 測點布置與編號

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

圖2 南北向測點位置及編號(單位：mm)Fig.2 North-south position and number of measuring points(unit:mm)

圖3 土壤層(≤30 cm)測點位置及編號(單位：mm)Fig.3 Soil layer (≤30 cm) measuring points location and number(unit:mm)

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

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

2 理論模型

2.1 物理模型

本模型中沒有考慮室內(nèi)作物，所以認(rèn)為地面為裸地。地面的熱濕過程包括：土壤垂直方向的熱量傳遞，與室內(nèi)空氣的熱對流交換，與室內(nèi)各個面的長波輻射熱交換，吸收的太陽直射輻射、散射輻射，土壤表面水汽相變傳遞的熱濕量等熱濕交換過程。圖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為太陽輻射量，J；H為對流輻射換熱量，J；L為長波輻射換熱量，J；γ為汽化潛熱，J/kg；W為水分傳輸量，kg.

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

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

模型簡化假設(shè)：

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

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

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

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

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

連續(xù)性方程：

(2)

動量方程：

(3)

(4)

(5)

(6)

能量方程：

(7)

k方程：

(8)

ε方程：

(9)

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

(10)

(11)

參考文獻[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)格劃分比較困難，計算結(jié)果不準(zhǔn)確，而非結(jié)構(gòu)網(wǎng)格能對不規(guī)整模型進行劃分，生成高質(zhì)量的網(wǎng)格，所以選擇使用ICEM CFD非結(jié)構(gòu)網(wǎng)格劃分。計算區(qū)域網(wǎng)格劃分如圖5所示。

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

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

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

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

2.5 數(shù)值計算方法

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

2.6 邊界條件

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

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

h=18.63v0.603.

(14)

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

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

對于迎風(fēng)面：

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

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

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

土壤的原始溫度場

(15)

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

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

表3 邊界條件Table 3 Boundary conditions

3 結(jié)果與討論

3.1 模擬驗證

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

圖6 溫度場驗證Fig.6 Temperature field verification

圖7 濕度場驗證Fig.7 Moisture field verification

3.2 土壤溫度場與濕度場的分布

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

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

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

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

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

3.3 土壤不同初始體積含水率下溫度場的變化

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

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

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

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

4 結(jié)論

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

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

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

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