路文淵 賈 蕾

0 引言

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

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

1 模型的建立

1.1 幾何模型

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

1.2 熱物理模型

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

圖1 日光溫室橫斷面及空氣溫度測試布點圖（單位：m）

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

1.3 數(shù)學模型

1.3.1 控制方程

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

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

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

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

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

1.3.2 輻射傳輸方程

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

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

1.3.3 壁面函數(shù)法

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

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

2 方程組的求解

2.1 網(wǎng)格的劃分

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

2.2 邊界條件的確定

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

2.2.1 外界氣象條件

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

2.2.2 圍護結構

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

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

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

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

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

圖3 計算的溫度分布云圖

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

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

圖5 模擬值與測試值的比較

4 結論及建議

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

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

［2］ 佟國紅，李保明.日光溫室溫度環(huán)境非穩(wěn)態(tài)模擬求解方法初探［A］.北京:中國農(nóng)業(yè)工程學會學術年會論文集第五分冊［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］ 賈 蕾.自然通風日光溫室熱環(huán)境的試驗研究與數(shù)值分析［D］.太原:太原理工大學，2009.

［5］ 陶文銓.數(shù)值傳熱學［M］.西安:西安交通大學出版社，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.