路文渊 贾 蕾

0 引言

近年来，CFD这一工具正在被用来模拟研究外部环境和自身结构对温室室内环境的影响。如何将CFD方法与具体问题很好地结合起来依然是我们当前研究的重点所在。自然通风、机械通风、温室结构等问题的CFD模拟研究均已展开，并取得了一定成就。国内学者李元哲等［1］最早于1994年运用热力学、传热学和建筑采光的基本理论，建立了日光温室微气候数学模型。2005年，佟国红等［2］对东北型日光温室温度环境进行了非稳态求解。国外自1989年，Okushima等［3］首次应用CFD研究自然通风下无植物单跨Venlo型温室的通风问题以来，也有许多学者用CFD对各种连栋温室的非稳态自然通风和热湿环境问题进行了研究。但由于日光温室是我国特有的温室结构形式，所以国外学者对其研究还甚少。

本文在对山西地区自然通风日光温室热湿环境试验研究的基础上，建立了冬季自然通风条件下该类日光温室的三维CFD模型，并对稳态条件下的温度场进行了模拟计算，用日光温室现场试验数据对该模型进行了验证。

1 模型的建立

1.1 几何模型

测试和建模用日光温室位于太原市郊区，坐北朝南，温室方位角为南偏西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轴，建立直角坐标系。该类日光温室围护结构材料的热物性参数见文献［4］，X方向截面几何尺寸及空气温度测试布点见图1。为验证模拟计算值，在X方向，相距2.3 m的三个截面上各设9个空气温度测点，图1中圆黑点是温度测点，黑点旁边数据是三个截面测点编号。

1.2 热物理模型

日光温室自然通风是包括质量、热量和动量交换，甚至还有物质相变的一个非常复杂的过程。室内外空气通过风口进出温室形成了自然通风，并作为载体参与了温室内的传热和传质过程。所以以室内外的空气作为研究对象，采用计算流体力学方法对室内外空气流动过程进行分析，就可以对日光温室自然通风过程进行全面的研究。本文主要研究冬季白天温室采用自然通风时室内温度场的稳态分布，因此，参考文献［4］作一些合理的简化和假设后，温室的热物理模型如图2所示。

图1 日光温室横断面及空气温度测试布点图（单位：m）

图2 日光温室的热物理模型

1.3 数学模型

1.3.1 控制方程

温室自然通风过程中空气的流动除室外新风与温室内的空气产生对流外，还有空气温差形成浮升力作用下的自然对流。在进行封闭腔内自然对流换热的数值计算时，为便于处理由于温差而引起的浮升力项，常常采用Boussinesq假设。按文献［5］这一假设由三部分组成:1)流体中的粘性耗散忽略不计;2)除密度外其他物性为常数;3)对密度仅考虑动量方程中与体积力有关的项，其余各项中的密度亦作为常数。

流体流动受质量、动量和能量守恒定律的支配，对于温室热环境这种带有传热的流动问题，从流动和传热问题的一般原理方程上着手进行研究，采用Boussinesq假设，建立以质量、动量与能量守恒方程为主体的方程组，联立求解。适用于本文所描述日光温室自然对流换热问题的各基本方程通用形式为［6］:

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

表1 中c1，c2，c3，cD是K—ε 湍流模型的系数;σK，σε，σH是施密特或普朗特数;μeff是有效粘性系数;Vj是速度分量;K是湍流动能;ε是湍流能量耗散率;H是比焓，H=CPT;μt是湍流粘性系数;vt是湍流动能粘性系数;gi是重力分量;CP是定压比热;φ是通用变量;β是体积膨胀系数;θ=H－H0是相对焓值，H0是参考值;c1，c2，c3，cD，σK，σε和 σH见文献［6］。

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

1.3.2 辐射传输方程

本文主要研究温室内热环境的变化规律，从能量角度研究温室内通过PVC膜接受的太阳辐射。离散坐标辐射模型可以模拟半透明介质内的辐射传递过程，所以选择离散坐标辐射模型对太阳辐射进行模拟求解。离散坐标模型求解的是从有限个立体角发出的辐射传输方程，每个立体角对应着笛卡儿坐标系下的固定方向S→。它不进行射线跟踪，相反，把辐射传输方程转化为空间坐标系下的辐射强度的输运方程，有多少个立体角方向S→，就求解多少辐射强度输运方程。该模型把沿S→方向传输的辐射方程视为某个场方程。这样，辐射传输方程化为:

入射辐射在温室内的分配和计算，远不如室外辐射那样具有成型的公式可以参考，因为形态各异的作物对日光的反射使得室内辐射的分配和计算变得相当复杂。本文对温室内太阳直射的分配及计算在吸收已有研究成果的基础上，做如下假定:1)各部分墙体的辐射量按漫灰表面辐射理论进行计算;2)忽略山墙、骨架等在室内作物及地表形成的阴影，认为各时刻通过各面入射的全部直射辐射在地表层均匀分配。采用离散坐标辐射模型对温室入射辐射进行分配计算。

1.3.3 壁面函数法

无论何种湍流模型，都是针对充分发展的湍流才有效的，也就是说，这些模型均是高Re数的湍流模型。而在壁面区，流动情况变化很大，特别是在粘性底层，湍流应力几乎不起作用，不能用以上K—ε模型来求解这个区域内的流动。目前工程计算解决这个问题的常用方法就是壁面函数法，对于湍流核心区的流动使用高Re数K—ε模型求解，而在壁面区不进行求解，直接使用半经验公式将壁面上的物理量与湍流核心区内的求解变量联系起来。这样，不需要对壁面区内的流动进行求解，就可直接得到与壁面相邻控制体积的接点变量值。在划分网格时，不需要在壁面区加密，只需要把第一个内节点布置在对数律成立的区域内，即配置到湍流充分发展的区域。壁面函数公式就好像一个桥梁，将壁面值同相邻控制体积的节点变量值联系起来。相对于采用低Re数K—ε模型，壁面函数法计算效率高，工程使用性强。所以本文也采用壁面函数法处理壁面值。

至此，由Boussinesq假设、基本控制方程组、湍流模型、壁面函数法、离散坐标辐射模型等构成了自然通风日光温室完整的热环境数学模型。利用FLUENT对方程组的离散、求解及相关参数的设定见文献［4］。

2 方程组的求解

2.1 网格的划分

本文选取温室中端部长度为9.2 m的实测封闭空间作为计算区域。流动与传热问题数值计算结果最终的精度及计算过程的效率，主要取决于所生成的网格和所采用的算法。由于该计算区域体积较大，为了保证计算的精度，计算域采用对适应不规则区域离散划分的非结构化非均匀网格，在入口和出口附近流场变化梯度较大的区域进行网格加密处理，并进行网格独立性考核，选择结果均匀性较好的网格为163 993个体网格，生成34 468个节点。

2.2 边界条件的确定

CFD数值模拟以室内外空气作为研究对象，其余如外界气象条件、温室围护结构和室内土壤都作为数值计算的边界条件来进行处理。

2.2.1 外界气象条件

室外风速和风向以速度进口边界条件输入，根据室外风速和风向的实测值以速度矢量的形式输入。入口速度为0.3 m/s，入口气流温度292.3 K;出口为出流边界;流场初始温度305 K。

2.2.2 围护结构

在稳态数值计算中，对后墙和覆盖材料用实测值给定，分别为301.5 K和310.2 K。由于计算和测试区域是被聚氯乙烯膜与邻室隔断，所以“山墙”处理为绝热，后屋面也处理为绝热。土壤层表面温度采用实测303 K值给定。

3 三维稳态CFD模拟与验证

选取某测试日13:00时的测试值进行模拟和对比，温度场计算结果见图3。由图3a)可知:温室内温度场分布的梯度非常明显，尤其是前屋面近壁区的等温面非常狭小，温差较大。从X方向截面上可以观察到明显的温度分层，除近壁区和入口处，温度分布随着空间高度的升高而增大，热浮力作用明显。通过Y方向的断面可以更好地观察温度场在竖直方向的分布情况(见图3b))，且明显看出沿Y轴方向入口处温度最低，从入口到出口的截面方向温度逐渐升高。根据Bot等［7］的研究，在风速低于0.5 m/s时热浮力是温室通风的主要驱动力，本次模拟计算时考虑了热浮力的影响，室内风速小于入口0.3 m/s，这与上述结论相一致。图3c)反映了不同高度水平面上温度场分布情况，基本上呈中心对称，与对称的计算域及边界条件设置吻合。

图4中line-4，line-5和 line-6分别表示 Z=240 cm，160 cm，80 cm高度不同Y的温度值。由图4可知:在靠近后墙处室内温度从高到低依次升高，但差距不大，这是由于后屋面和卷起的保温被遮挡太阳辐射，靠近后墙处不同高度接受的太阳辐射强度不等所致;从距离后墙200 cm处开始，室内不同高度温度变化趋势又变为从高到低依次降低，且差距较大，这是热空气上浮和进风口温度较低共同影响的结果。这些均与实测温度场值变化趋势相同。不同的变化趋势在中间必有一个温度沿高度几乎不变的交汇点，该点为距离后墙1.20 m左右，与后屋水平投影长度1.38 m相当。

图5是从模拟结果中截取的测试点的温度值与实测值的比较。由图5可知，模拟结果较为准确地反映了温室内不同点温度变化趋势，整体比试验数据偏小。最大误差为5.2℃，且都出现在最上部靠近膜边界处的测点，最小误差为1℃，出现在底部紧邻后墙的测点。

图3 计算的温度分布云图

由于温室自然通风是在室外自然风引起的风压和室内外温度差引起的热压共同作用下产生的室内外空气通过风口流入和流出温室的运动，因自然通风口附近的气流运动十分复杂，风口处的空气流动方向和大小也随时在变化，数值模拟很难完全满足这种边界条件，且模拟过程中未考虑植物对室内温度场和速度场的影响，这些都可能是造成模拟值比实测值偏小的原因。

图4 X=460 cm截面不同高度的温度值比较

图5 模拟值与测试值的比较

4 结论及建议

虽然数值模拟的温室内温度分布与现场实测结果在数值上存在差异，但除边界点外，最大误差在4.5℃以内，而且从空间分布的总体趋势来看是一致的，说明对温室热环境的数值模拟是成功的，所建立的CFD模型及其边界条件是有效的。但日光温室热惰性大，影响因素多，不可能处于热稳定状态，CFD作为一种模拟仿真工具，随着计算机技术的飞速发展，这一方法也得到了长足进步，今后应该考虑植物对流动和换热的影响，并考虑非稳态模拟。

［1］ 李元哲，吴德让，于 竹.日光温室微气候的模拟与实验研究［J］.农业工程学报，1994，10(1):130-136.

［2］ 佟国红，李保明.日光温室温度环境非稳态模拟求解方法初探［A］.北京:中国农业工程学会学术年会论文集第五分册［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］ 贾 蕾.自然通风日光温室热环境的试验研究与数值分析［D］.太原:太原理工大学，2009.

［5］ 陶文铨.数值传热学［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.