杜 軍, 趙忠超, 黃 波

(1.江蘇科技大學(xué) 能源與動力工程學(xué)院, 江蘇 鎮(zhèn)江 212003)(2.哈爾濱理工大學(xué) 機械動力工程學(xué)院, 黑龍江 哈爾濱 150006)

溫室園藝在世界各地均有廣泛的發(fā)展和應(yīng)用,現(xiàn)代化農(nóng)業(yè)、林業(yè)和養(yǎng)殖業(yè)都與溫室業(yè)的發(fā)展有著密切的聯(lián)系．溫室系統(tǒng)是一個復(fù)雜的物理和生物系統(tǒng),包含許多非線性動態(tài)過程,如:動態(tài)傳熱過程、水份輸運及平衡過程、作物光合作用等過程[1]．許多研究者對這樣的過程進(jìn)行了不同類別的研究[2-4]．文獻(xiàn)[5]中的溫室氣溫數(shù)學(xué)模型,顯得有些簡單,考慮的影響因素較少,尤其沒有考慮溫室內(nèi)土壤對氣溫的相關(guān)影響．文獻(xiàn)[6]對溫室內(nèi)的不同對象建立能量平衡方程,主要給出了溫室采光面和土壤表層的能量平衡方程,但沒有求解．這種模型方法研究對象不明確,對溫室系統(tǒng)來講可以列出多個類似方程,聯(lián)立求解很困難,使問題趨向復(fù)雜化．文獻(xiàn)[7]應(yīng)用建筑學(xué)上的反應(yīng)系數(shù)法建立不同研究對象的熱平衡方程,然后聯(lián)立求解[7]．該數(shù)學(xué)模型大量涉及到溫室各內(nèi)表面溫度,而表面溫度在當(dāng)前測試條件和手段下,具有很難測量準(zhǔn)確的特點,不易檢驗．而且這樣的數(shù)學(xué)模型在計算時需要用到模擬日前3天和模擬日的室外逐時氣溫,輸入條件比較復(fù)雜．

文獻(xiàn)[8]針對不同結(jié)構(gòu)溫室對風(fēng)速引起的熱流變化進(jìn)行了研究,目的是建立有風(fēng)干燥夜間情況下的子模型．文獻(xiàn)[9]則建立了溫室內(nèi)不同對象的能量方程,包括土壤表面．除了模型研究外,很多學(xué)者研究采用多種技術(shù)給溫室供熱以抵御寒冷的室外氣候[10-12]．文獻(xiàn)[13]采用氣源熱泵在冬季的夜間給一處151 m2的實驗溫室加熱,其空間氣溫分布比較均勻,并且溫度波動較?。墨I(xiàn)[14]采用地源熱泵和潛熱儲存罐給溫室供熱,其熱泵系數(shù)達(dá)到2.3～3.8．

文中采用溫室內(nèi)氣溫和土溫的相關(guān)性數(shù)學(xué)模型,設(shè)計數(shù)值模擬流程,進(jìn)行數(shù)值求解．以北方一處冬季實際運行中的傾斜式溫室為對象進(jìn)行實驗測試及研究．結(jié)合數(shù)值模擬結(jié)果,對溫室內(nèi)溫光特性進(jìn)行研究,分析日光溫室的太陽輻射能及室內(nèi)氣溫和土溫的變化規(guī)律．

1 數(shù)學(xué)模型的選取

典型傾斜式溫室能量守恒如圖1．

圖1 溫室能量守恒Fig.1 Energy balance of a greenhouse

參照文獻(xiàn)[15],將溫室沿東西方向分割為N個小溫室,認(rèn)為小溫室內(nèi)氣溫均勻一致．以小溫室做為一個系統(tǒng),溫室內(nèi)空氣和土壤做為研究對象．從能量平衡角度出發(fā),認(rèn)為能量的凈增量(減量)等于空氣內(nèi)能的增量(減量),有陽光照射時采光面?zhèn)鳠嵯禂?shù)變化及北墻蓄熱量對氣溫的影響單獨考慮,根據(jù)能量守恒定律：

(1)

式中:Gy為進(jìn)入溫室太陽凈輻射量(W);Qs為溫室維護(hù)結(jié)構(gòu)總散熱量(W);Qt為土壤表面對空氣的傳熱量(W);Qz為作物對空氣的傳熱量(W);Qf為通風(fēng)引起的空氣顯熱和潛熱量(W);Qr為供熱量(W);Ta為溫室內(nèi)空氣溫度(℃);ca為溫室內(nèi)空氣定壓比熱(J/(kg℃));va為溫室內(nèi)空氣體積(m3);ρa為溫室內(nèi)空氣密度(kg/m3);m為溫室內(nèi)空氣對太陽凈輻射吸收份額;t為時間(s).

以土壤表面為坐標(biāo)原點,垂直向下為正半軸,建立Z坐標(biāo)系．認(rèn)為土壤內(nèi)部溫度變化滿足Z方向一維熱傳導(dǎo)方程:

(2)

式中:T為溫室內(nèi)土壤溫度(℃);Z為坐標(biāo)軸;ct為溫室內(nèi)土壤比熱(J/(kg℃));ρt為溫室內(nèi)土壤密度(kg/m3);λ為溫室內(nèi)土壤導(dǎo)熱系數(shù)(W/(m℃)).

溫室內(nèi)空氣和土壤可以通過傅里葉定律建立聯(lián)系:

(3)

式中:Tt為溫室內(nèi)土壤表面溫度(℃);qa為溫室內(nèi)空氣對土壤表面熱流密度(W/m2);qy為太陽輻射對土壤表面熱流密度(W/m2).

式(1,2,3)構(gòu)成了傳熱微分方程組,參照文獻(xiàn)[15],給出該微分方程組的模擬解．

(4)

(5)

式中:η=η1+η2+η3;η1=(Aqkq+Apkp+Ackc+ρacaV+rVε)/(ρacava);η2=Atht/(ρacava);η3=Azhz/(ρacava);M=m/(ρacava);N=1/(ρacava).

2 數(shù)學(xué)模型的求解

求解思路:將式(4)中的Gy，T0，Tz，Tt和Qr離散為時間坐標(biāo)上的階躍函數(shù),即時間步長Δt內(nèi)為常數(shù)輸入量,則使問題得以簡化．如此處理后,仍不能正常求解,因求解式(4)時將產(chǎn)生1個積分常數(shù),而很難確定初值條件．因此聯(lián)立tN+1時刻和tN時刻解的表達(dá)式,采用Δt考慮時間變量得出氣溫模擬公式．對式(5)采用有限差分法在時間和一維空間坐標(biāo)上進(jìn)行離散化,推導(dǎo)出相應(yīng)的數(shù)值模擬公式．

對式(4)進(jìn)行求解，得

(6)

Ta=ce-ηt+U

(7)

式中：c為氣溫積分常數(shù).因此有

TaN-UN=ce-ηtN，TaN+1-UN+1=ce-ηtN+1

消去c得

TaN+1=e-ηΔtTaN+UN+1-e-ηΔtUN

(8)

式中：Δt為時間步長,Δt=tN+1-tN,(s).

溫室內(nèi)空氣隨著溫度的變化,含濕量也發(fā)生相應(yīng)變化,這一變化過程會吸收或放出熱量,從而導(dǎo)致已知時刻的下一時刻氣溫值產(chǎn)生一定誤差．當(dāng)正負(fù)誤差逐點積累后,有可能在某一時刻互相抵消,也有可能在某一時刻達(dá)到峰值而產(chǎn)生較大誤差．為了減小這種潛在的誤差,采用如下方法修正．

圖2中,TaN為已知時刻氣溫值,對應(yīng)含濕量為xN；TaN+1為下一時刻理論計算氣溫值,對應(yīng)含濕量為xN+1；TaZ為修正后氣溫值,對應(yīng)含濕量為xZ．

圖2氣溫修正示意圖
Fig.2Diagramoftemperaturecorrection

根據(jù)TaN→TaN+1理論計算吸收或放出的熱量應(yīng)該等于修正后實際吸收或放出的熱量,于是有

ρavaca(TaN+1-TaN)=ρavaca(TaZ-TaN)+

var(xZ-xN)

(9)

(10)

在式(10)中TaZ和xZ均為未知數(shù),而只有一個方程,在數(shù)學(xué)上無法求解．現(xiàn)在采用計算機編程查找的方法來解決:先比較xN和xN+1的大小,可以令x=min{xN,xN+1}．在限定范圍{xN,xN+1}內(nèi),令循環(huán)步長為0.01,即:x=x+0.01,依據(jù)x值,確定對應(yīng)的TaZ,由TaZ確定對應(yīng)的x′值,定義查找精度ε=|x-x′|≤0.01,能夠得到氣溫修正值TaZ．當(dāng)然,也可以令溫度值循環(huán),用含濕量校正．

(11)

整理式(11)得:

(12)

圖3 土壤分層Fig.3 Diagram of soil layers

如將土壤表面定義為0層,向下依次為1層、2層、3層…,則Tt=T(0),Z0=0,式(23)即為:

(13)

(14)

因此,土壤表面溫度在N時刻輸入量為

(15)

土壤某一深度處溫度基本保持不變,可假定為常數(shù)[16];文獻(xiàn)[17]也指出土壤溫度在地面以下某一深度基本恒定．事實上,這一深度及溫度隨不同地點和不同季節(jié)等條件將有一些變化,但在所研究期間和特定條件下可通過測量得到土壤恒溫Th和恒溫深度h．

土壤內(nèi)部溫度場求解過程如圖4所示,括號內(nèi)為時刻,下角標(biāo)為Z軸坐標(biāo),單位為cm．

圖4 土壤內(nèi)部溫度場求解Fig.4 Diagram of solving the temperature fields inside the soil

假定在1:00時輸入所需的初始參數(shù),令Δt=3 600 s,可以得出任意時刻氣溫和土壤溫度值．即輸入1:00時刻的氣溫Ta(1)和土溫T(1),可以得到氣溫Ta(2),Ta(3)，…，Ta(24)以及土溫T(2),T(3)，…，T(24)．

對溫室內(nèi)太陽輻射能量分配的計算方法,按照晴天狀況計算．這是由于研究對象是冬季嚴(yán)寒地區(qū)實際運行溫室,一般在陰天時溫室不會揭簾運行．而晴天時,進(jìn)入溫室的陽光多為直射輻射,方向性很強,在計算溫室內(nèi)太陽輻射能量的分配時,假設(shè)溫室各內(nèi)表面均為漫灰表面．采用角系數(shù)方法考慮溫室內(nèi)輻射太陽能的分配．

在圖5中,線段vw和zw均已知,vw為溫室總高(脊高),zw是人行過道的寬度．在三角形Δuvw中,yw=zw/cosA;Δuxy和Δuvw為相似三角形,最后得到太陽光線能夠照射到的北墻高度:

Hso=xy=vw-zw·tgL/cosA

(16)

式中：Hso為太陽照射北墻高度(m).

根據(jù)以上分析,引入不完全角系數(shù),即:采光面發(fā)出的太陽輻射量只能到達(dá)土壤和北墻沒有陰影的部分．結(jié)合圖6,晴天時角系數(shù)為

(17)

(18)

圖5 溫室陽光陰影Fig.5 Diagram of sunlight shadow inside a greenhouse

圖6 溫室內(nèi)表面橫剖Fig.6 Cross-section diagram of a greenhouse

可見光輻射能量Gyk到達(dá)土壤及北墻的輻射分量為

Etk=GykZhsXab,bc

(19)

Eqk=Gyk(H0+HzZhc)Xab,cd

(20)

式中:Etk為可見光到達(dá)土壤的輻射分量(W);Eqk為可見光到達(dá)北墻的輻射分量(W);Hz為作物平均高度,測量值(m);H0為無遮擋高度,H0=Hso-Hz,如果Hso

作物吸收的一次可見光輻射分量為

Ezk=Gyk-Etk-Eqk

(21)

式中：Ezk為作物吸收的可見光輻射分量(W).

作物水平方向和垂直方向投影孔隙率可以根據(jù)作物在水平方向和垂直方向的陰影率求得．作物水平方向和垂直方向的陰影率定義為作物在兩個方向上一定面積內(nèi),投影面積所占的百分比如圖7所示．所以得到:

Zys=Szs/Ss

(22)

Zyc=Szc/Sc

(23)

Zhs=1-Zys

(24)

Zhc=1-Zyc

(25)

式中:Zys為作物水平方向陰影率;Zyc為作物垂直方向陰影率;Szs為作物水平方向投影面積(m2);Szc為作物垂直方向投影面積(m2);Ss為水平方向一定面積(m2);Sc為垂直方向一定面積(m2).

圖7 作物投影Fig.7 Projection diagram of a plant

類似可求紅外光輻射能量Gyh到達(dá)土壤及北墻的輻射分量為:

Eth=KvGyhZhsXab,bc

(26)

Eqh=KvGyh(H0+HzZhc)Xab,cd

(27)

式中:Eth為紅外光到達(dá)土壤的輻射分量(W);Eqh為紅外光到達(dá)北墻的輻射分量(W);KV為溫室內(nèi)氣體輻射減弱率,KV=1-αa;αa為溫室內(nèi)氣體吸收率.

溫室內(nèi)氣體吸收的一次紅外光輻射分量為

Eah=Gyh-Eth-Eqh

(28)

式中:Eah為溫室內(nèi)氣體吸收的紅外光輻射分量(W).

因此,只計一次反射輻射,可得最終土壤及北墻對太陽輻射的凈吸收量、作物對可見光輻射的凈吸收量、溫室內(nèi)氣體對紅外光輻射的凈吸收量:

Et=αt[Etk+Eth+ρq(Eqk+EqhKv)ZhsXcd,bc]

(29)

Eq=αq[Eqk+Eqh+ρt(Etk+EthKv)ZhcXbc,cd]

(30)

Ez=Ezk+ρqEqk(1-ZhsXcd,bc)+

ρtEtk(1-ZhcXbc,cd)g

(31)

Ea=Eah+αa(ρqEqh+ρtEth)g

(32)

式中:Et為土壤對太陽輻射的凈吸收量(W);Eq為北墻對太陽輻射的凈吸收量(W);Ez為作物對可見光輻射的凈吸收量(W);Ea為溫室內(nèi)氣體對紅外光輻射的凈吸收量(W).

近似取可見光輻射量等于紅外光輻射量,代入以上各式整理得:

(33)

(34)

(1-ZhsXcd,bc)Xab,cd+

ρtZhs(1-ZhcXbc,cd)Xab,bc]

(35)

(36)

3 實驗測試及分析研究

實驗溫室位于北緯45°41′的哈爾濱市遠(yuǎn)郊,溫室東西延長106 m,溫室內(nèi)占地面積800 m2,種植550株油桃樹,每株間距1 m．溫室保溫草簾較薄,平均為4.6 cm厚,而且北坡由木板加稻草組成,透風(fēng)量很大,整體保溫防寒性能極差．該溫室為典型的土木結(jié)構(gòu)溫室,其結(jié)構(gòu)及供熱方式如圖8．

圖8 實驗溫室結(jié)構(gòu)Fig.8 Structure diagram of a testing greenhouse

3.1 實驗測點布置

進(jìn)行實驗時,沿溫室長度方向取6個截面,對每個截面上3點氣溫和8點土溫進(jìn)行逐時測試,共計66個測點,如圖9,10所示．所有測點(包括室外氣溫)的溫度傳感器為Cu100銅電阻,溫度采集用WSC-411p型數(shù)字顯示儀,其精度為±0.5%Fs,分辨率為0.1℃．

圖9 測試截面分布Fig.9 Distribution of testing sections

圖10 測點分布Fig.10 Distribution of testing points

3.2 實驗結(jié)果及模擬研究

運行工況如圖11，12所示,圖11中含義數(shù)值1代表停爐,2代表揭簾,黑色線段為各自區(qū)間．一般來講,在嚴(yán)重陰天或下雪等不利的天氣情況以及積雪未清理等特殊情況發(fā)生時,溫室運行中不適宜揭簾．圖12表明,由于受停電、跑水、壓爐以及受熱面積灰等多種因素影響,鍋爐供熱系統(tǒng)的運行功率與設(shè)計的平均功率相差較大．

圖11 運行工況Fig.11 Operating conditions

圖12 供熱功率Fig.12 Heating power

氣溫采用平均氣溫值,土溫采用5cm土溫值,實驗結(jié)果與模擬結(jié)果如圖13，14，15所示．

圖13 進(jìn)入溫室太陽凈輻射量Fig.13 Net solar radiation inside a greenhouse

圖14 模擬氣溫Fig.14 Simulated air temperatures

圖15 模擬土溫Fig.15 Simulated soil temperatures

由圖13可以看出,在連續(xù)的3天內(nèi),進(jìn)入溫室的太陽凈輻射量在同一時間段內(nèi)相差不大,一般在12∶00～13∶00太陽輻射最強,這段時間前后太陽輻射逐漸減弱．由圖14和15可見,模擬結(jié)果與測試數(shù)據(jù)擬合得很好,氣溫和土溫模擬誤差分別在±1.54℃,±1.42℃以內(nèi)．氣溫隨太陽輻射量增加升溫幅度較大,在中午可以達(dá)到20 ℃左右,在無太陽輻射且不供熱時,氣溫在3～4 h內(nèi)可降到0℃．5 cm土層溫度受太陽輻射量影響也很大,最大升溫幅度達(dá)到5 ℃,在蓋簾2～3 h內(nèi)5 cm土層降溫幅度稍大,連續(xù)幾日有陽光照射時,對整體土層溫度的升高很有利．因此,合理掌握好揭簾和蓋簾時間,最大限度利用太陽能,是溫室經(jīng)濟(jì)運行的關(guān)鍵措施．

3.3 實驗及模擬結(jié)果分析

按照圖8的整體程序流程圖,參照文獻(xiàn)[15，18]進(jìn)行了程序計算,得到了進(jìn)入溫室的凈輻射量以及氣溫和土溫的模擬值．對于氣溫和土溫的相對分布均勻主要受到兩個條件的影響:一是供熱時采用了大熱管供熱系統(tǒng),溫度沿溫室長度方向的均勻性較好;二是計算時室外氣溫和太陽輻射采用了逐時測量值,以便較好驗證模型的準(zhǔn)確性．對于計算的氣溫和土溫絕對誤差似乎不大,但是由于個別時點的溫度絕對值較小,所以相對誤差超過了100 %,但不會影響整體評價分析的效果．

4 結(jié)論

1) 氣溫隨時間沿著東西向變化,在同一天內(nèi)氣溫隨太陽輻射量增加升溫幅度較大,在中午可以達(dá)到20 ℃左右;

2) 土溫隨著土壤深度變化,土壤溫度分層趨勢明顯,在80 cm土層基本達(dá)到恒定溫度6.3 ℃,5 cm土層溫度日變化較大;

3) 土壤溫度與室內(nèi)氣溫相互影響,但各層土壤溫度受室內(nèi)氣溫影響遠(yuǎn)不如受陽光照射影響大;

4) 在連續(xù)的3天內(nèi),進(jìn)入溫室的太陽凈輻射量在同一時間段內(nèi)相差不大,一般在12 ∶00～13 ∶00太陽輻射最強．

參考文獻(xiàn)(References)

[1] Young P C, Lees M J.Simplicity out of complexity in glasshouse climate modelling[J].ActaHortic, 1996,406:15-28.

[2] Bouchouicha M, Lafont F, Balmat J F.Neural networks, fuzzy logic and genetic algorithms for greenhouse identification[C]∥In: 2ndinternationalconferenceonelectro-technicalandautomaticcontrol(JTEA).Tunisia:[s.n.], 2002:356-2.

[3] Pessel N, Duplaix J, Balmat J F, et al.Data analysis for Neuro-Fuzzy model approach[C]∥SOFA′05IEEEinternationalworkshoponsoftcomputingapplications.Szeged-Arad:[s.n.], August, 2005: 44-50.

[4] Balas M, Duplaix J, Balas S.Modelling the heat flow of the greenhouses[C]∥SOFA′05IEEEinternationalworkshoponsoftcomputingapplications.Szeged-Arad:[s.n.], August, 2005: 37-43.

[5] Schmidt U.Greenhouse climate control with a combine model of greenhouse and plant by using online measurement of leaf temperature and transpiration[J].ActaHorticulturae, 1996,406:89-98.

[6] Gabriel P A, Manuel P, Alberto V.Mediterranean greenhouse energy balance[J].ActaHortic, 1998,456:375-382

[7] 郭慧卿,李振海,張振武.日光溫室溫度環(huán)境動態(tài)模擬[J].沈陽農(nóng)業(yè)大學(xué)學(xué)報,1994(4):438-443.

Guo Huiqin, Li Zhenhai, Zhang Zhenwu.Dynamic simulation of temperatures in a greenhouse[J].JornalofShenyangAgricultureUniversity,1994(4):438-443.(in Chinese)

[8] Balas Marius M, Duplaix Jean, Bouchouicha Moez, et al.Structural modeling of the wind′s influence over the heat flow of the greenhouses[J].JIntellFuzzySyst, 2008,19:29-40.

[9] Pita G P, Manuel P, Alberto V.Mediterranean greenhouse energy balance[J].ActaHortic, 1998,456:375-382.

[10] Ruijs M A, Vermeulen P M.Simulation model for determining the effects of farm specific aspects on energy consumption in glasshouses[J].ActaHortic, 1998,456:515-521.

[11] Nielsen B, Madsen H.Identification of a continuous time model of the heat dynamics of a greenhouse[J].ActaHortic, 1996,406:39-48.

[12] Miguel A F, Silva A M.Modelling of transport phenomena near greenhouse screens[J].ActaHortic, 1998,456:371-374.

[13] Tong Y, Kozai T, Nishioka N, et al.Greenhouse heating using heat pumps with a high coefficient of performance (COP) [J].BiosystEng, 2010,106:405-411.

[14] Benli Huseyin.Energetic performance analysis of a ground-source heat pump system with latent heat storage for a greenhouse heating[J].EnergConversManage, 2011,52:581-589.

[15] Du Jun, Bansal P, Huang B.Simulation model of a greenhouse with a heat-pipe heating system[J].AppliedEnergy, 2012,93:268-276.

[16] Tap R F,Van L G, Willigenburg.Experimental results of receding horizon optimal control of greenhouse climate[J].ActaHortic, 1996,406:229-238.

[17] Levit H J,Piacentini R D.Greenhouses heating by simultaneous use of thermal screens and underground pipe[J].ActaHortic, 1994,357:335-342.

[18] 趙忠超,杜軍,史自強,等.乘員負(fù)荷分布方式對客室內(nèi)氣流組織的影響[J].江蘇科技大學(xué)學(xué)報：自然科學(xué)版,2009,23(6):520-526.

Zhao Zhongchao,Du Jun,Shi Ziqiang， et al.Effects of the airflow with different distributions of heat load generated by passengers in passenger compartment[J].JournalofJiangsuUniversityofScienceandTechnology:NaturalScienceEdition,2009,23(6):520-526.(in Chinese)