日光溫室穩(wěn)態(tài)溫度場(chǎng)的數(shù)值模擬

路文淵 賈 蕾

0 引言

近年來，CFD這一工具正在被用來模擬研究外部環(huán)境和自身結(jié)構(gòu)對(duì)溫室室內(nèi)環(huán)境的影響。如何將CFD方法與具體問題很好地結(jié)合起來依然是我們當(dāng)前研究的重點(diǎn)所在。自然通風(fēng)、機(jī)械通風(fēng)、溫室結(jié)構(gòu)等問題的CFD模擬研究均已展開，并取得了一定成就。國(guó)內(nèi)學(xué)者李元哲等［1］最早于1994年運(yùn)用熱力學(xué)、傳熱學(xué)和建筑采光的基本理論，建立了日光溫室微氣候數(shù)學(xué)模型。2005年，佟國(guó)紅等［2］對(duì)東北型日光溫室溫度環(huán)境進(jìn)行了非穩(wěn)態(tài)求解。國(guó)外自1989年，Okushima等［3］首次應(yīng)用CFD研究自然通風(fēng)下無植物單跨Venlo型溫室的通風(fēng)問題以來，也有許多學(xué)者用CFD對(duì)各種連棟溫室的非穩(wěn)態(tài)自然通風(fēng)和熱濕環(huán)境問題進(jìn)行了研究。但由于日光溫室是我國(guó)特有的溫室結(jié)構(gòu)形式，所以國(guó)外學(xué)者對(duì)其研究還甚少。

本文在對(duì)山西地區(qū)自然通風(fēng)日光溫室熱濕環(huán)境試驗(yàn)研究的基礎(chǔ)上，建立了冬季自然通風(fēng)條件下該類日光溫室的三維CFD模型，并對(duì)穩(wěn)態(tài)條件下的溫度場(chǎng)進(jìn)行了模擬計(jì)算，用日光溫室現(xiàn)場(chǎng)試驗(yàn)數(shù)據(jù)對(duì)該模型進(jìn)行了驗(yàn)證。

1 模型的建立

1.1 幾何模型

測(cè)試和建模用日光溫室位于太原市郊區(qū)，坐北朝南，溫室方位角為南偏西5°，脊高3.3 m，北墻高2.3 m;后墻和山墻均為石灰漿砌實(shí)心粘土紅磚，厚度分別為0.80 m和0.37 m;后屋面厚度0.27 m(石棉瓦加紅磚)，水平投影長(zhǎng)度1.38 m，仰角40°;前屋面塑料薄膜使用三層共擠聚乙烯流滴性的PVC膜，厚度為0.002 m，夜間使用保溫棉被保溫;溫室跨度為8 m，長(zhǎng)度為56 m。分別以溫室長(zhǎng)、寬和高三個(gè)方向?yàn)閄軸，Y軸和Z軸，建立直角坐標(biāo)系。該類日光溫室圍護(hù)結(jié)構(gòu)材料的熱物性參數(shù)見文獻(xiàn)［4］，X方向截面幾何尺寸及空氣溫度測(cè)試布點(diǎn)見圖1。為驗(yàn)證模擬計(jì)算值，在X方向，相距2.3 m的三個(gè)截面上各設(shè)9個(gè)空氣溫度測(cè)點(diǎn)，圖1中圓黑點(diǎn)是溫度測(cè)點(diǎn)，黑點(diǎn)旁邊數(shù)據(jù)是三個(gè)截面測(cè)點(diǎn)編號(hào)。

1.2 熱物理模型

日光溫室自然通風(fēng)是包括質(zhì)量、熱量和動(dòng)量交換，甚至還有物質(zhì)相變的一個(gè)非常復(fù)雜的過程。室內(nèi)外空氣通過風(fēng)口進(jìn)出溫室形成了自然通風(fēng)，并作為載體參與了溫室內(nèi)的傳熱和傳質(zhì)過程。所以以室內(nèi)外的空氣作為研究對(duì)象，采用計(jì)算流體力學(xué)方法對(duì)室內(nèi)外空氣流動(dòng)過程進(jìn)行分析，就可以對(duì)日光溫室自然通風(fēng)過程進(jìn)行全面的研究。本文主要研究冬季白天溫室采用自然通風(fēng)時(shí)室內(nèi)溫度場(chǎng)的穩(wěn)態(tài)分布，因此，參考文獻(xiàn)［4］作一些合理的簡(jiǎn)化和假設(shè)后，溫室的熱物理模型如圖2所示。

圖1 日光溫室橫斷面及空氣溫度測(cè)試布點(diǎn)圖（單位：m）

圖2 日光溫室的熱物理模型

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

1.3.1 控制方程

溫室自然通風(fēng)過程中空氣的流動(dòng)除室外新風(fēng)與溫室內(nèi)的空氣產(chǎn)生對(duì)流外，還有空氣溫差形成浮升力作用下的自然對(duì)流。在進(jìn)行封閉腔內(nèi)自然對(duì)流換熱的數(shù)值計(jì)算時(shí)，為便于處理由于溫差而引起的浮升力項(xiàng)，常常采用Boussinesq假設(shè)。按文獻(xiàn)［5］這一假設(shè)由三部分組成:1)流體中的粘性耗散忽略不計(jì);2)除密度外其他物性為常數(shù);3)對(duì)密度僅考慮動(dòng)量方程中與體積力有關(guān)的項(xiàng)，其余各項(xiàng)中的密度亦作為常數(shù)。

流體流動(dòng)受質(zhì)量、動(dòng)量和能量守恒定律的支配，對(duì)于溫室熱環(huán)境這種帶有傳熱的流動(dòng)問題，從流動(dòng)和傳熱問題的一般原理方程上著手進(jìn)行研究，采用Boussinesq假設(shè)，建立以質(zhì)量、動(dòng)量與能量守恒方程為主體的方程組，聯(lián)立求解。適用于本文所描述日光溫室自然對(duì)流換熱問題的各基本方程通用形式為［6］:

其中，φ，Γφ，eff，SN和 SB給在表1 中。

表1 中c1，c2，c3，cD是K—ε 湍流模型的系數(shù);σK，σε，σH是施密特或普朗特?cái)?shù);μeff是有效粘性系數(shù);Vj是速度分量;K是湍流動(dòng)能;ε是湍流能量耗散率;H是比焓，H=CPT;μt是湍流粘性系數(shù);vt是湍流動(dòng)能粘性系數(shù);gi是重力分量;CP是定壓比熱;φ是通用變量;β是體積膨脹系數(shù);θ=H－H0是相對(duì)焓值，H0是參考值;c1，c2，c3，cD，σK，σε和 σH見文獻(xiàn)［6］。

表1 φ，Γφ，eff，SN 和 SB

1.3.2 輻射傳輸方程

本文主要研究溫室內(nèi)熱環(huán)境的變化規(guī)律，從能量角度研究溫室內(nèi)通過PVC膜接受的太陽輻射。離散坐標(biāo)輻射模型可以模擬半透明介質(zhì)內(nèi)的輻射傳遞過程，所以選擇離散坐標(biāo)輻射模型對(duì)太陽輻射進(jìn)行模擬求解。離散坐標(biāo)模型求解的是從有限個(gè)立體角發(fā)出的輻射傳輸方程，每個(gè)立體角對(duì)應(yīng)著笛卡兒坐標(biāo)系下的固定方向S→。它不進(jìn)行射線跟蹤，相反，把輻射傳輸方程轉(zhuǎn)化為空間坐標(biāo)系下的輻射強(qiáng)度的輸運(yùn)方程，有多少個(gè)立體角方向S→，就求解多少輻射強(qiáng)度輸運(yùn)方程。該模型把沿S→方向傳輸?shù)妮椛浞匠桃暈槟硞€(gè)場(chǎng)方程。這樣，輻射傳輸方程化為:

入射輻射在溫室內(nèi)的分配和計(jì)算，遠(yuǎn)不如室外輻射那樣具有成型的公式可以參考，因?yàn)樾螒B(tài)各異的作物對(duì)日光的反射使得室內(nèi)輻射的分配和計(jì)算變得相當(dāng)復(fù)雜。本文對(duì)溫室內(nèi)太陽直射的分配及計(jì)算在吸收已有研究成果的基礎(chǔ)上，做如下假定:1)各部分墻體的輻射量按漫灰表面輻射理論進(jìn)行計(jì)算;2)忽略山墻、骨架等在室內(nèi)作物及地表形成的陰影，認(rèn)為各時(shí)刻通過各面入射的全部直射輻射在地表層均勻分配。采用離散坐標(biāo)輻射模型對(duì)溫室入射輻射進(jìn)行分配計(jì)算。

1.3.3 壁面函數(shù)法

無論何種湍流模型，都是針對(duì)充分發(fā)展的湍流才有效的，也就是說，這些模型均是高Re數(shù)的湍流模型。而在壁面區(qū)，流動(dòng)情況變化很大，特別是在粘性底層，湍流應(yīng)力幾乎不起作用，不能用以上K—ε模型來求解這個(gè)區(qū)域內(nèi)的流動(dòng)。目前工程計(jì)算解決這個(gè)問題的常用方法就是壁面函數(shù)法，對(duì)于湍流核心區(qū)的流動(dòng)使用高Re數(shù)K—ε模型求解，而在壁面區(qū)不進(jìn)行求解，直接使用半經(jīng)驗(yàn)公式將壁面上的物理量與湍流核心區(qū)內(nèi)的求解變量聯(lián)系起來。這樣，不需要對(duì)壁面區(qū)內(nèi)的流動(dòng)進(jìn)行求解，就可直接得到與壁面相鄰控制體積的接點(diǎn)變量值。在劃分網(wǎng)格時(shí)，不需要在壁面區(qū)加密，只需要把第一個(gè)內(nèi)節(jié)點(diǎn)布置在對(duì)數(shù)律成立的區(qū)域內(nèi)，即配置到湍流充分發(fā)展的區(qū)域。壁面函數(shù)公式就好像一個(gè)橋梁，將壁面值同相鄰控制體積的節(jié)點(diǎn)變量值聯(lián)系起來。相對(duì)于采用低Re數(shù)K—ε模型，壁面函數(shù)法計(jì)算效率高，工程使用性強(qiáng)。所以本文也采用壁面函數(shù)法處理壁面值。

至此，由Boussinesq假設(shè)、基本控制方程組、湍流模型、壁面函數(shù)法、離散坐標(biāo)輻射模型等構(gòu)成了自然通風(fēng)日光溫室完整的熱環(huán)境數(shù)學(xué)模型。利用FLUENT對(duì)方程組的離散、求解及相關(guān)參數(shù)的設(shè)定見文獻(xiàn)［4］。

2 方程組的求解

2.1 網(wǎng)格的劃分

本文選取溫室中端部長(zhǎng)度為9.2 m的實(shí)測(cè)封閉空間作為計(jì)算區(qū)域。流動(dòng)與傳熱問題數(shù)值計(jì)算結(jié)果最終的精度及計(jì)算過程的效率，主要取決于所生成的網(wǎng)格和所采用的算法。由于該計(jì)算區(qū)域體積較大，為了保證計(jì)算的精度，計(jì)算域采用對(duì)適應(yīng)不規(guī)則區(qū)域離散劃分的非結(jié)構(gòu)化非均勻網(wǎng)格，在入口和出口附近流場(chǎng)變化梯度較大的區(qū)域進(jìn)行網(wǎng)格加密處理，并進(jìn)行網(wǎng)格獨(dú)立性考核，選擇結(jié)果均勻性較好的網(wǎng)格為163 993個(gè)體網(wǎng)格，生成34 468個(gè)節(jié)點(diǎn)。

2.2 邊界條件的確定

CFD數(shù)值模擬以室內(nèi)外空氣作為研究對(duì)象，其余如外界氣象條件、溫室圍護(hù)結(jié)構(gòu)和室內(nèi)土壤都作為數(shù)值計(jì)算的邊界條件來進(jìn)行處理。

2.2.1 外界氣象條件

室外風(fēng)速和風(fēng)向以速度進(jìn)口邊界條件輸入，根據(jù)室外風(fēng)速和風(fēng)向的實(shí)測(cè)值以速度矢量的形式輸入。入口速度為0.3 m/s，入口氣流溫度292.3 K;出口為出流邊界;流場(chǎng)初始溫度305 K。

2.2.2 圍護(hù)結(jié)構(gòu)

在穩(wěn)態(tài)數(shù)值計(jì)算中，對(duì)后墻和覆蓋材料用實(shí)測(cè)值給定，分別為301.5 K和310.2 K。由于計(jì)算和測(cè)試區(qū)域是被聚氯乙烯膜與鄰室隔斷，所以“山墻”處理為絕熱，后屋面也處理為絕熱。土壤層表面溫度采用實(shí)測(cè)303 K值給定。

3 三維穩(wěn)態(tài)CFD模擬與驗(yàn)證

選取某測(cè)試日13:00時(shí)的測(cè)試值進(jìn)行模擬和對(duì)比，溫度場(chǎng)計(jì)算結(jié)果見圖3。由圖3a)可知:溫室內(nèi)溫度場(chǎng)分布的梯度非常明顯，尤其是前屋面近壁區(qū)的等溫面非常狹小，溫差較大。從X方向截面上可以觀察到明顯的溫度分層，除近壁區(qū)和入口處，溫度分布隨著空間高度的升高而增大，熱浮力作用明顯。通過Y方向的斷面可以更好地觀察溫度場(chǎng)在豎直方向的分布情況(見圖3b))，且明顯看出沿Y軸方向入口處溫度最低，從入口到出口的截面方向溫度逐漸升高。根據(jù)Bot等［7］的研究，在風(fēng)速低于0.5 m/s時(shí)熱浮力是溫室通風(fēng)的主要驅(qū)動(dòng)力，本次模擬計(jì)算時(shí)考慮了熱浮力的影響，室內(nèi)風(fēng)速小于入口0.3 m/s，這與上述結(jié)論相一致。圖3c)反映了不同高度水平面上溫度場(chǎng)分布情況，基本上呈中心對(duì)稱，與對(duì)稱的計(jì)算域及邊界條件設(shè)置吻合。

圖4中l(wèi)ine-4，line-5和 line-6分別表示 Z=240 cm，160 cm，80 cm高度不同Y的溫度值。由圖4可知:在靠近后墻處室內(nèi)溫度從高到低依次升高，但差距不大，這是由于后屋面和卷起的保溫被遮擋太陽輻射，靠近后墻處不同高度接受的太陽輻射強(qiáng)度不等所致;從距離后墻200 cm處開始，室內(nèi)不同高度溫度變化趨勢(shì)又變?yōu)閺母叩降鸵来谓档?，且差距較大，這是熱空氣上浮和進(jìn)風(fēng)口溫度較低共同影響的結(jié)果。這些均與實(shí)測(cè)溫度場(chǎng)值變化趨勢(shì)相同。不同的變化趨勢(shì)在中間必有一個(gè)溫度沿高度幾乎不變的交匯點(diǎn)，該點(diǎn)為距離后墻1.20 m左右，與后屋水平投影長(zhǎng)度1.38 m相當(dāng)。

圖5是從模擬結(jié)果中截取的測(cè)試點(diǎn)的溫度值與實(shí)測(cè)值的比較。由圖5可知，模擬結(jié)果較為準(zhǔn)確地反映了溫室內(nèi)不同點(diǎn)溫度變化趨勢(shì)，整體比試驗(yàn)數(shù)據(jù)偏小。最大誤差為5.2℃，且都出現(xiàn)在最上部靠近膜邊界處的測(cè)點(diǎn)，最小誤差為1℃，出現(xiàn)在底部緊鄰后墻的測(cè)點(diǎn)。

圖3 計(jì)算的溫度分布云圖

由于溫室自然通風(fēng)是在室外自然風(fēng)引起的風(fēng)壓和室內(nèi)外溫度差引起的熱壓共同作用下產(chǎn)生的室內(nèi)外空氣通過風(fēng)口流入和流出溫室的運(yùn)動(dòng)，因自然通風(fēng)口附近的氣流運(yùn)動(dòng)十分復(fù)雜，風(fēng)口處的空氣流動(dòng)方向和大小也隨時(shí)在變化，數(shù)值模擬很難完全滿足這種邊界條件，且模擬過程中未考慮植物對(duì)室內(nèi)溫度場(chǎng)和速度場(chǎng)的影響，這些都可能是造成模擬值比實(shí)測(cè)值偏小的原因。

圖4 X=460 cm截面不同高度的溫度值比較

圖5 模擬值與測(cè)試值的比較

4 結(jié)論及建議

雖然數(shù)值模擬的溫室內(nèi)溫度分布與現(xiàn)場(chǎng)實(shí)測(cè)結(jié)果在數(shù)值上存在差異，但除邊界點(diǎn)外，最大誤差在4.5℃以內(nèi)，而且從空間分布的總體趨勢(shì)來看是一致的，說明對(duì)溫室熱環(huán)境的數(shù)值模擬是成功的，所建立的CFD模型及其邊界條件是有效的。但日光溫室熱惰性大，影響因素多，不可能處于熱穩(wěn)定狀態(tài)，CFD作為一種模擬仿真工具，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，這一方法也得到了長(zhǎng)足進(jìn)步，今后應(yīng)該考慮植物對(duì)流動(dòng)和換熱的影響，并考慮非穩(wěn)態(tài)模擬。

［1］ 李元哲，吳德讓，于 竹.日光溫室微氣候的模擬與實(shí)驗(yàn)研究［J］.農(nóng)業(yè)工程學(xué)報(bào)，1994，10(1):130-136.

［2］ 佟國(guó)紅，李保明.日光溫室溫度環(huán)境非穩(wěn)態(tài)模擬求解方法初探［A］.北京:中國(guó)農(nóng)業(yè)工程學(xué)會(huì)學(xué)術(shù)年會(huì)論文集第五分冊(cè)［C］.2005:95-99.

［3］ Okushima L，Sase S，Nara M.A support system for natural ventilation design of greenhouses based on computational aerodynamics［J］.Acta Horticulturae，1989(284):129-136.

［4］ 賈 蕾.自然通風(fēng)日光溫室熱環(huán)境的試驗(yàn)研究與數(shù)值分析［D］.太原:太原理工大學(xué)，2009.

［5］ 陶文銓.數(shù)值傳熱學(xué)［M］.西安:西安交通大學(xué)出版社，2001.

［6］ Qingyan Chen.The methematical foundation of the CHAMPION SGE computer code(revision)，1987(3):119-121.

［7］ Bot G.P.A.Greenhouse climate:from physical proecesses to a dynamic model［D］.Doetoral Thesis:Agricultural University of Wageninge，1983.