基于CFD的植物工厂管道通风模拟及优化

吴晨溶 程瑞锋 方 慧* 杨其长 张 晨

(1.中国农业科学院 农业环境与可持续发展研究所，北京 100081;2.农业农村部设施农业节能与废弃物处理重点实验室，北京 100081)

植物工厂是指通过设施内高精度环境控制实现农作物周年连续生产的高效农业系统，利用计算机对植物发育过程的温度、湿度、光照、CO2浓度以及营养液等环境条件进行自动控制，使设施内植物生长发育不受或很少受自然条件制约的省力型生产方式[1]。由于植物工厂充分运用了现代工程、生物、环境和信息等手段，技术高度密集，近年来一直被国际上公认为是现阶段设施农业的最高发展阶段，是衡量一个国家农业技术水平的重要标志，受到世界各国的高度重视[2]。

植物工厂内部气流场的分布情况对植物生长有重要影响[3]。研究表明[4-6]植物工厂内部气流速度在0.01～1.30 m/s范围内，叶片光合作用及蒸腾作用随着气流速度的增长逐渐增强，更适宜作物生长。但在传统通风模式的植物工厂中，内部气流往往分布不均，尤其是栽培架对风速的阻挡，作物光合作用和蒸腾作用受到抑制，导致作物产生疾病，其中以烧尖现象最为普遍，严重影响作物的产量和品质，降低了植物工厂经济效益[7]。大量研究表明，增加作物冠层气流可在一定程度上改善作物冠层微环境，促进作物光合作用和蒸腾作用。Goto等[8]设计了一种在作物冠层上方的垂直通风模式，采用单根主管连接6根支管的方式，将气流从主管通向支管气孔自上而下吹入作物冠层，该技术在一定程度上改善了作物冠层气流环境。Shibata等[9]同样设计了垂直气流系统，该系统主要由1台风机、1根集气管和6根通风管组成，每根通风管分布9个气孔，风机吹出的空气经过集气管吹入通风管，最后通过气孔吹出，该方法提高了叶片的蒸腾速率并在一定程度上抑制了烧尖现象的发生。Nishikawa等[10]将作物放置于独立的旋转平台上，以2 r/min的速度匀速旋转，一定程度上打破了冠层边界层阻力，作物产量比静置时提高了20%。Hesham等[11]提出多风扇系统，将7个风扇分别安装在栽培层的前后侧，产生来自相反水平方向的气流，改善了气流环境，促进了植物工厂内作物生长，缓解作物烧尖现象。李琨等[12]设计了一套根际通风系统，空气经风机导入栽培板与营养液面中间的空气层中，经栽培板定植孔周围预留的通风孔自下而上排出，该系统在解决传统环控方式通风温控不均匀的同时，对地上部及地下部多种微环境参数调控起到了积极作用。尽管上述通风方式对气流均匀性有一定积极作用，但试验周期较长，且有的需借助复杂的机械装置来实现通风。

近年来，计算流体力学软件(Computational fluid dynamics，CFD)以周期短、成本低、准确性高的特点被广泛应用于温室和植物工厂的气流场与温度场的模拟[13]。Lim T等[14]利用CFD模拟了植物工厂内中间过道4种不同进出口位置下气流分布情况。Baek M S等[15]借助CFD在人工光植物工厂中将空调、排风扇、内部风扇的开启设置为不同案例，模拟出植物工厂内部温度场及气流场，并选出最适宜植物生长的通风案例。张芳等[16]通过CFD找到了自然通风条件下大跨度温室内温度场和气流场的分布规律，以及不同室外风速条件下通风口开度对大跨度温室温度场和气流场的影响。刘焕等[17]利用CFD找到了人工光植物工厂内最佳气流循环模式。Zhang Y等[18]设计了带有三排气孔的空气管，空气垂直向下流到作物冠层表面，并借助CFD软件模拟植物冠层的气流场，最终确定2根空气管的方案气流最为均匀。

本研究拟将植物工厂中的单层栽培架为研究对象，设计管道通风方案，并设置不同进风口风速进行正交试验，利用CFD仿真软件对不同方案下的作物冠层表面气流场进行模拟仿真，分析不同管孔直径、管孔数量及进风口风速对作物冠层气流分布的影响，确定最佳管道通风方案，以期为改善植物工厂微环境和缓解作物烧尖现象提供一定的理论参考。

1 材料与方法

1.1 试验装置

试验单层栽培架位于中国农业科学院农业环境与可持续发展研究所植物工厂内，单层栽培架规格(长×宽×高)为1.4 m×0.7 m×0.4 m，采用穿孔聚乙烯管道进行通风，管道直径为0.05 m，试验作物为生菜(Lactucasativacv. Tiberius)。本试验只考虑作物冠层和冠层上部空气层，灯板到栽培板的距离一般为0.3 m，成熟期的生菜平均株高为0.1 m[19]，作物冠层和冠层上部空气层的高度分别设为0.1 m和0.2 m，管道架在栽培架外侧，管道气孔水平方向的中心线与空气层和作物冠层相交切线重合。在Workbench中的Geometry模块中创建单层植物栽培架模型(图1)，以模型东北角为原点O(0,0,0)正南方向为X轴正向，垂直向上为Y轴正向，正西方向为Z轴正向。

1.空气层；2.作物冠层；3.管道1.Air interior; 2.Crop interior; 3.Tube图1 试验单层栽培架结构图Fig.1 Schematic diagram of single shelf cultivation frame

1.2 试验设置

在温室内使用穿孔聚乙烯管道来分配空气已变得越来越普遍，但在穿孔管道设计中，气孔的直径、数量、形状、间距等因素都会影响空气流通率和气流速度。在Well C的试验中得出，气孔截面积之和与管道口面积之比小于1.5时风速排放更为均匀[20]。为满足合适的管孔径比，本研究通过计算设置了3组方案，所有方案的管道长度都为1.4 m，管道其余参数见表1，并在Workbench中的Geometry模块中创建了3组方案的管道模型(图2)。通风管道的气流由进风口进入管道后经气孔流入作物冠层内部，进风口风速分别设为2、3、4和5 m/s进行模拟。Kitaya等[21]研究表明，作物冠层表面的适宜风速为0.3～1.0 m/s，所以在本研究中，按风速将植物冠层表面划分为3个区域：停滞区(≤0.3 m/s)、适宜区(0.3～1.0 m/s)和高速区(≥1.0 m/s)。试验所用泡沫栽培板尺寸(长×宽×高)为0.7 m×0.65 m×0.014 m，栽培密度17株/m2，每个栽培板16个定植孔(4行4列)，定植孔直径为0.03 m，选用8个定植孔定植生菜，在栽培板上呈三角形分布。

表1 不同管道通风方案的管道基本参数Table 1 Basic parameters of different tube ventilation cases

图2 不同管道通风方案的管道模型Fig.2 Tube models of different tube ventilation cases

1.3 测点布置

试验中选取X=0.35 m截面，在Y=0.5、1.0、1.5、2.0、2.5 m高度处各均匀布置4个测点，共20个测点(P1，P2，…，P20)；选取Z=0.4 m截面，在Y=0.5、1.0、1.5、2.0、2.5 m高度处各均匀布置2个测点，共10个测点(P21，P22，…，P30)。2个截面共30个测点测量风速值进行模拟验证(图3)。气流速度测量使用热线风速仪(Climomaster 6501-BG，日本加野麦克斯公司)，风速测量范围为0.01～5.00 m/s，风速测量公差为标准值的±2%，工作温度范围为-20～70 ℃，相对湿度范围为2%～98%，气压范围为-5～5 kPa。

P1，P2，…，P30为风速测点P1，P2，…，P30 are the measuring points of air velocity图3 X=0.35 m截面(a)和Z=0.4 m截面(b)风速测点布置Fig.3 Air velocity measuring points of sections of X=0.35 m (a) and Z=0.4 m (b)

2 数学模型

2.1 CFD建模与网格划分

利用Workbench中的Geometry模块和ANSYS Mesh模块对栽培装置进行几何建模和网格划分，采用FLUENT 19.2对栽培区域进行模拟计算。在Geometry模块中创建与单层栽培架原型相同尺寸的栽培装置模型，包括作物冠层、空气层和管道3个部分。将构建的模型输入到ANSYS Mesh模块中，利用Proximity and Curvature方法进行网格划分(图4)，管道部分网格进行局部加密，3组方案模型的网格平均偏斜度和最大偏斜度分别为0.732 5和0.74、0.747 5和0.76、0.722 5和0.74，网格数为465 453～502 345(表2)，网格质量均为优[22]。

图4 栽培装置网格划分图Fig.4 Mesh in single shelf cultivation frame

表2 不同管道通风方案的网格参数Table 2 Mesh quality of different tube ventilation cases

2.2 控制方程

将栽培层内的通风过程视为三维、稳态、不可压缩的湍流流动过程；湍流模型选用k-ε模型(realizable)；使用有限体积法对控制方程进行离散；用SIMPLE算法(压力耦合方程组的半隐式方法)计算求解压力耦合方程；气流数值计算按照流体力学守恒数学公式计算，即连续性方程、动量方程、能量方程，符合质量、动量、能量守恒定律[23]，公式如下：

连续性方程

(1)

动量方程

(2)

(3)

(4)

能量方程

(5)

式中：ρ为密度，kg/m3；t为时间，s；v为速度矢量；u，υ，ω分别为流体质点速度在x，y，z方向上的分量；μ为黏性系数；S为源项；T为温度，K；k为导热系数，W/(m·K)；Cp为比热容，J/(kg·K)；p为压力，Pa。

2.3 参数设置及计算方法

将植物栽培装置模型进行简化，作物冠层上部空气层设为空气介质，作物冠层设为多孔介质，黏性阻力系数和惯性阻力系数分别为25和1.3[19]；气流由进风口进入管道内部通过气孔吹入作物冠层。管道进风口设置为速度入口(Velocity-inlet)，进风速度分别为2、3、4、5 m/s；栽培装置四周气流出口设置为压力出口(Pressure-outlet)，压力为0 Pa，所有材料的热物性参数见表3。选用SIMPLE算法进行压力速度耦合运算。空间离散设置中，梯度项选用最小二乘法；动量、能量和粘性项都选用一阶迎风格式，使其更快收敛，将能量项的松弛因子设置为10-6，动量项和粘性项皆设置为10-3，迭代步数设为300步。

表3 材料的热物性参数Table 3 Thermal physical parameters of materials

3 CFD仿真结果

3.1 CFD模型验证

将30个测点的风速模拟值与实测值进行对比(图5)，平均绝对误差为0.02 m/s，均方根误差为0.024，其中误差最大的测点是第9个测点，实测风速值为0.15 m/s，CFD模拟风速值为0.085 m/s，绝对误差为0.065 m/s。总体上看，CFD风速模拟值与风速实测值吻合度良好。因此本试验所构建的单层栽培架实际内部气流场分布情况与模拟情况吻合良好，能够进行数据分析与和后续的模拟试验。

图5 单层栽培架内风速实测值与模拟值对比图Fig.5 Comparison between measured and simulated air velocity in the single shelf cultivation frame

3.2 6孔通风管道作物冠层气流分布

对6孔通风管道方案进行模拟，作物冠层气流分布云图见图6。气孔附近风速最大，作物冠层平面气孔位置呈6个气流束，风速沿气流束逐渐减弱。除6个气流束所在位置外，作物冠层平面其他区域风速较小，气流在气孔与气孔之间存在较大面积的停滞区，且在冠层平面边缘管道两端位置存在部分零风速区域，整体风速并不均匀。随着进风口风速增大，作物冠层高速区面积也逐渐增大，停滞区面积逐渐减小。在进气速度为2、3、4和5 m/s时，作物冠层平面风速平均值分别为0.28、0.43、0.56和0.66 m/s，适宜区的面积比分别为32.7%、40.1%、37.9%和34.2%，其中当进风口风速为3 m/s时，适宜区面积比最大为40.1%(表4)。

图6 6孔通风管道不同进风口风速(v)下作物冠层表面风速分布云图Fig.6 The contours of canopy surface air velocity under different inlet air velocity of 6 holes tube

表4 6孔通风管道不同进风口风速下作物冠层表面各区域的面积比Table 4 The area percentage of each area at the canopy surface under different inlet air velocity of 6 holes tube

3.3 12孔通风管道作物冠层气流分布

对12孔通风管道方案进行模拟，结果见图7：最大风速仍分布在气孔附近，当进风口风速为 2 m/s 时，作物冠层平面气孔位置呈12个气流束且风速较小，与之相比，当进风口风速为3、4和5 m/s时，气流分布较均匀，随着进风口风速的逐渐增大，气孔间的停滞区面积减小，作物冠层适宜区的面积增大，但冠层平面边缘管道两端位置仍存在零风速区域。作物冠层平面风速平均值分别为0.34、0.52、0.71和0.89 m/s，适宜区的面积比分别为47.0%、68.4%、62.2%和52.9%，其中当进风口风速为3 m/s时，面积比最大为68.4%(表5)。

图7 12孔通风管道不同进风口风速(v)下作物冠层表面风速分布云图Fig.7 Contours of canopy surface air velocity under different inlet air velocity of 12 holes tube

表5 12孔通风管道不同进风口风速下作物冠层表面各区域的面积比Table 5 The area percentage of each area at the canopy surface under different inlet air velocity of 12 holes tube

3.4 25孔通风管道作物冠层气流分布

25孔通风管道方案模拟结果(图8)表明，当进风口风速为2 m/s时，作物冠层平面气孔位置呈较明显的气流束，与之相比，当进风口风速为3、4、5 m/s 时，风速分布没有气流束，而是呈片状梯度分布，气流分布较均匀，随着进风口风速的逐渐增大，气孔间停滞区面积降低，冠层平面边缘管道两端零风速面积也相应减小。作物冠层平面风速平均值分别为0.36、0.55、0.74和0.94 m/s，适宜区的面积比分别为54.4%、86.1%、85.8%和75.0%，其中当进风口风速为3 m/s时，面积比最大为86.1%(表6)。

图8 25孔通风管道不同进风口风速(v)下作物冠层表面风速分布云图Fig.8 Contours of canopy surface air velocity under different inlet air velocity of 25 holes tube

表6 25孔通风管道不同进风口风速下作物冠层表面各区域的面积比Table 6 The area percentage of each area at the canopy surface under different inlet air velocity of 25 holes tube

3.5 气流场对比分析

由冠层平面风速分布云图(图6～8)对比可知，气孔孔径大数量少则风速范围差异大，气孔孔径小数量多则风速范围差异小。当风速为5 m/s时，25孔通风管道作物冠层平面风速平均值最高，为0.94 m/s，适宜区的面积比为75.0%，而6孔和12孔通风管道分别为34.2%和52.9%；当风速为4 m/s 时，25孔通风管道作物冠层平面风速平均值最高，为0.74 m/s，适宜区的面积比为85.8%，而6孔和12孔通风管道分别为37.9%和62.2%；当风速为3 m/s时，25孔通风管道作物冠层平面风速平均值最高，为0.55 m/s，适宜区的面积比为86.1%，而6孔和12孔通风管道分别为40.1%和68.4%；当风速为2 m/s时，25孔通风管道作物冠层平面风速平均值最高，为0.36 m/s，适宜区的面积比为54.4%，而6孔和12孔通风管道分别为32.7%和47.0%。3种管道在进风口风速逐渐降低的情况下，作物冠层平面平均风速也随之逐渐降低，但适宜区的面积比逐渐上升，在进风口风速为 3 m/s 时均达到最高值，在进风口风速为2 m/s时，适宜区的面积比均下降。综合考虑，进风口风速为3 m/s 的25孔管道通风方案最佳(图9)。

图9 不同管道通风方案在不同进风口风速下作物冠层表面适宜区的面积比Fig.9 Percentage of suitable area at canopy surface of different tube ventilation cases under different air velocity

4 结论与讨论

本研究利用CFD模拟了植物工厂内单层栽培装置作物冠层内部的气流分布情况，并将模拟值与实测值进行对比，风速平均绝对误差为0.02 m/s，均方根误差为0.024，CFD风速模拟值与实测值吻合，模拟效果良好。

基于植物工厂内作物冠层管道通风气孔的孔径与数量设计了3种(6孔、12孔和25孔)管道通风模式，4种不同的进风口风速(2、3、4和5 m/s)进行正交试验。在Hesham的试验中，传统通风模式冠层表面的平均风速为0.23 m/s，停滞区面积比为70%[11]。本试验所有管道通风方案的平均风速均高于传统通风模式，停滞区面积占比均低于传统通风模式，即所有管道通风方案的平均风速及气流均匀性均优于传统通风模式。当进风口风速为3 m/s时，25孔通风管道冠层平面适宜区的面积比为86.1%，为最佳管道通风方案。

但当进风口风速为4 m/s时，25孔通风管道冠层平面适宜区的面积比为85.8%，停滞区为8.7%，而高速区面积比为5.5%，与进风口为3 m/s时的面积比进行对比，差值分别为0.3%、2.5%和-2.8%；由此可看出，当进风口风速为4 m/s时，冠层平面适宜区和停滞区的百分比都有所减小，但差值较小，虽然高速区的面积比有所增加，但都聚集在气孔附近，对作物冠层产生影响较小，所以下一步试验将进风口风速分别设为3和4 m/s进行模拟。

本试验只研究了在固定管道高度和角度条件下，不同进风口风速、孔径和孔数对作物冠层平面气流分布的影响，但在实际生产中LED灯的散热会对栽培装置内部的气流产生一定影响。因此，下一步试验将嵌入LED灯热源模型，研究管道与冠层栽培板的高度以及气孔的倾斜度对作物冠层平面和冠层内部气流分布的影响，从而更加精确地模拟栽培装置内部环境，最终确定最优的管道通风方案。