食用菌培养室热湿耦合场的数值模拟研究

张璐杨 周 勃 孙成才 邢韶缨 王肖宇

（1.沈阳工业大学建筑与土木工程学院 沈阳 110870；2.辽宁沣知稼农业科技发展有限公司 沈阳 110870）

0 引言

21 世纪后，我国的食用菌产值比重在农业中已高居第五位，食用菌产业正在由小型企业、分散农户的生产方式逐步转向规模化的工厂化生产[1,2]。食用菌工厂化生产是在封闭的培养室内，人们自由控制建筑内部温湿度、通风和光照等条件，通过机械化、标准化的栽培达到绿色、高效的生产[3]。虽然可以通过环境控制来克服外界不利环境带来的影响，但目前为了满足市场上对食用菌的巨大需求，培养室内的培养瓶密度不断增加，导致室内温度和湿度的分布不均匀，这既增加了工厂环境调控的成本和难度，也影响了出菇的质量与产量[4,5]。因此，对食用菌培养室的热湿耦合场进行模拟分析，有助于提高食用菌成品率、降低工厂运维成本，对建立合理的空气调节系统以及配置环境调控设备具有至关重要的实际工程意义。

CFD 数值模拟不仅可以有效地弥补实验周期较长、成本较高的缺点，而且可以对农业建筑的室内环境进行可靠的分析[6]。最早将CFD 运用于农业建筑的研究可追溯至上世纪八十年代，Okushima等人模拟了无作物影响时斜顶温室的环境[7]。而随着数值仿真技术的发展，模型逐渐由二维至三维，由单一模拟至后期优化设计。有学者建立了南方连栋温室的三维CFD 模型，结合试验及仿真结果研究了不同机械通风方案下的降温效果，优化后可为温室节省约60%的能耗[8]。Boulard 等人重点考虑了作物的蒸腾、光合作用的耦合，采用CFD 方法模拟了玻璃温室内的流场分布情况[9]。这些研究说明，CFD 作为研究建筑内部环境的有效工具，在模拟温室内流场分布情况方面具有较高的可信度。

事实上，在不同的生长阶段食用菌与环境之间的显热、潜热交换并不相同[10]。关于动植物对建筑内部环境的影响，国内外学者大多选择采用CFD方法直接建模[11]、将其简化为发热板状结构[12]或者近似为多孔介质[13]三种方法。直接建模的方法虽使得结果更加准确，但对于生理结构较为复杂的食用菌而言较为繁琐。简化为发热板状结构的方法又无法准确地展现动植物对室内空气流场的影响。因此，本文选择将食用菌生长区域近似为多孔介质的方法，这既能反映食用菌对空气的阻力作用，又能增加模拟的准确度。

由于大多数学者并未将食用菌散湿量考虑进去[14,15]，但食用菌培养室内的通风过程往往涉及到温度、相对湿度、CO2浓度和光照等条件的相互作用，是一个由多因素、多物理场之间相互影响、相互耦合的复杂空气流动过程。因此，本文通过CFD三维数值模拟以及现场实测的方法，结合食用菌特定时期的发热量与散湿量，分析食用菌培养室热湿耦合场的分布规律，并提出两种培养室内气流组织的改造方案，为食用菌提供最佳的生长环境。

1 CFD 模型的建立

1.1 食用菌培养室概况

为验证模型的准确性，于2020年12月3日至12月6日在食用菌培养室进行测试，测试地点位于北京市昌平区小汤山农业科技园（40.1°N，116.2°E）。园内分布4 个大型食用菌厂房，文中以其中生育室内的金针菇培养室为原型进行研究，培养室长10m，宽10.5m，高6m，围护结构采用150mm 厚聚氨酯板。由于金针菇在生长过程中会散发大量的热量和湿量，因此需要制冷系统来保持良好的循环通风，进而维持室内的温湿度。室内每台空气调节设备的额定制冷量为14.8kW，距离屋顶0.2m，进风口直径0.42m，共4 台。图1为食用菌厂房和金针菇培养室实物图。

图1 食用菌培养室实物图Fig.1 Physical map of edible fungus culture room

在培养室内部，设置有7 排6 列6 层放置金针菇培养瓶的菇架，共有培养瓶24192 瓶。每排菇架长9m，宽1m，高3.3m。测试期间该生育室内的金针菇处于采收前期，最佳生长温度为5.5℃。培养室内需保持5.5℃左右，相对湿度90%左右。

1.2 几何模型与网格划分

根据金针菇培养室的室内布置与几何尺寸，本文未考虑测试期间处于关闭状态的加湿设备，仅考虑每层菇架中的金针菇培养瓶区域，且培养瓶之间的细小缝隙忽略不计。培养室密封性良好，仅通过风机进行通风，因此将整个金针菇培养室内的空气作为计算域，在ANSYS ICEM 中建立三维模型。图2为培养室几何结构图，图3为室内空气流动过程的示意图。

图2 培养室几何结构图Fig.2 Geometry of the incubation chamber

图3 室内空气流动过程Fig.3 Indoor air flow process

在利用数值计算解决室内流动与传热问题时，网格的划分将直接影响数值模拟计算的速度与准确度。本研究选用非结构化四面体网格对所选的食用菌培养室进行划分，根据食用菌培养室对环境精度的要求，对风口处以及近壁面区域等流场变化梯度较大的地方进行网格加密处理。

为得到网格独立性分析结果，对该培养室几何模型划分了6 种网格模型并进行网格独立性检验。本文的主要目的是为了解决食用菌培养室的热湿耦合场分布的问题，因此选择将整个培养室以及食用菌生长区域的平均温度、平均相对湿度的计算结果作为网格独立性的检验标准。通过计算达到稳态后，得到的最终结果如图4所示。

图4 网格独立性检验结果Fig.4 Grid independence test results

从图4中可以看出，随着网格数量以不同的程度逐渐增加，各组结果的变化趋势基本趋于稳定。其中，当网格数量小于150 万时，计算结果仍存在微小变化，而当网格数量增加到150 万以上时，各组的计算结果已基本保持不变。因此可以认为当网格数量在150 万以上时，模型的计算结果与网格数量无关。为提高计算效率，本文选择网格数相对较少的网格模型进行模拟计算，共生成266688 个节点、1667615 个网格，网格质量符合要求。网格划分结果如图5所示。

图5 网格的划分Fig.5 Grid partition and encryption

1.3 数学模型

1.3.1 基本控制方程

食用菌培养室内的空气流动是以三大方程（质量守恒方程、能量守恒方程以及动量守恒方程）来描述的，若流动中包含其余组分，则还要满足组分守恒定律。本文将室内气体看作水蒸气与空气的混合，并视为定常不可压缩流体，三大方程与组分方程可表述为式（1）形式[16,17]：

式中：ρ为密度，kg/m ；Iφ为广义扩散系数；u为速度矢量，m/s；φ为通用变量，即速度、温度等；Sφ为广义源项。

1.3.2 湍流模型

在自然环境以及许多工程设备中流体的流动常为湍流流动，食用菌培养室内气流将不可避免的发生湍流现象，因此将室内空气流动看作湍流流动。本文选择标准k-ε模型来求解培养室内的气流流动，与之对应的输运方程为[18-19]：

式中：Gb为由浮力引起的湍动能k的产生项；YM为可压湍流中脉动扩张的贡献；Sk、Sε为源项；Gk为由平均速度引起的湍动能k的产生项；σε、σk、G1ε、G2ε、G3ε为经验常数，通常σε= 1.3、σk=1.0、G1ε=1.44、G2ε=1.92、G3ε=0.09。

1.3.3 多孔介质模型

为了描述食用菌对培养室内热湿耦合场的影响，应尽可能使模型与实际情况相符合，但若根据食用菌本身的实际大小来建立模型，不仅会增大网格数量、影响收敛速度，而且会影响模拟的准确度。因此目前在模拟障碍物对室内环境的影响时，大多采用将其简化为多孔介质的方法，通过设定阻力系数来表示对空气流动的阻力作用[20]。多孔介质模型是在动量方程中添加了一个与速度有关的源项[21]：

式中：Si为动量方程源项；v为速度值；μ为动力黏度；D、C 为指定的矩阵；等号右边第一项为黏性损失，第二项为惯性损失。

对于简单均匀的多孔介质，上式可改写为：

式中：1/α为黏性阻力系数；C2为惯性阻力系数。

1.4 培养室内环境参数的测定与方法

培养室内部测点的位置如图6所示。为反映食用菌生长区域的环境状态，按菇架的摆放位置共布置18 个测点，沿着培养室的走廊方向，距地面1.1m布置9 个测点，距地面2.2m 布置9 个测点，每个测点均选择在培养室走廊的中间位置。室内温湿度的测量选用泰仕TES1361 记忆式温湿度表，温度测量精度为±0.8℃、测量范围为-20～60℃，相对湿度测量精度为±3%、测量范围为10%～95%。室内的风速由天跃TY9900 数字微风仪测得，测量范围为0～10m/s，精度为±0.2 m/s。壁面表面温度采用福禄克FLUKE 63 红外温度计，测量范围为-32～535℃、精度为±1℃。测量期间每隔1 小时记录一次数据，各参数每次均测量三次，取平均值作为最终值进行计算。各测点的空间位置如表1所示。

表1 各测点的空间坐标（单位：m）Table 1 Spatial coordinates of each measuring points

图6 测点的布置Fig.6 Layout of measuring points

1.5 计算方法与边界条件

采用基于有限体积法的CFD 软件ANSYS FLUENT 进行计算。为简化域内流场结构，将食用菌的生长区域设为多孔介质区域，并通过速度与压降的拟合关系确定瓶体阻力系数[22]。在模拟过程中由于金针菇呼吸、蒸腾作用所释放的热量与湿量以多孔介质区域的能量源项的形式输入，将培养瓶区域设为热源，根据工厂每瓶金针菇散热量0.15W 的条件，得到多孔介质区域内热源为32W/m 。金针菇的散湿量通过下式计算[23]。

式中：m为排风质量流量，kg/s；Δd为培养室进、出口空气含湿量之差，g/kg。

将围护结构设为“wall”壁面边界条件，多孔介质区域与室内其余空气流域间的边界设为“interior”内部边界条件，空调设备的进风口设为“velocity-inlet”速度进口边界条件，回风口选择设置为“pressure-outlet”压力出口边界条件，其余壁面均认为绝热。其中的参数条件均按现场实测数据进行设定，边界条件的设置如表2所示。

表2 边界条件与参数设置Table 2 Boundary conditions and parameter settings

计算时采用三维稳态方法求解，同时打开能量方程和组分运输方程来展现室内热湿的相互耦合。湍流模型选用标准k-ε模型，选择“Full Buoyancy Effects”选项，近壁面区域采用标准壁面函数进行计算。离散格式选择二阶迎风格式，压力-速度耦合项选用SIMPLEC 半隐式算法求解。收敛残差设置为0.001，其中的能量项为10-6，最终经过约1800次迭代达到收敛。

2 结果分析与模型验证

2.1 数值模拟结果分析

根据以上设置进行计算，为反映室内空气流场的分布情况，选择各测点所在的垂直切面Y=2.6m、Y=5.6m、Y=8.6m 来对结果进行分析，各切面的温度场、湿度场及速度场云图如图7所示。由图7（a）中可以看出，菇架上方温度分布较为均匀，约5.8℃，食用菌生长区域温度较高，最低处约5.51℃，最高处约6.17℃，温差约0.66℃，培养室内最高温度约为6.2℃。三个切面的平均温度约为5.86℃，最高约6.17℃。沿气流流动方向，培养室中心区域的温度高于四周。

室内相对湿度场云图如图7（b）所示，室内相对湿度场的分布与温度场分布基本一致，温度高处区域相对湿度较高，培养室内部食用菌生长区域最低相对湿度约为90.42%。并且结合速度场云图7（c）来看，由于金针菇培养瓶的摆放比较密集，当低温气体经进风口流入房间后，无法均匀地穿过整个金针菇生长区域，每层培养瓶间所产生的热量和湿量难以带出，导致在培养室中部的培养瓶附近温度与相对湿度较高，生长区域温差较大。因此，有必要对食用菌培养室的通风设计进行改造优化。

图7 原始模型模拟结果Fig.7 Simulation results of the original model

2.2 模型的验证

本文引入绝对误差（absolute error,AE）、平均相对误差（average relative error,ARE）、最大相对误差（maximum relative error,MRE）、均方根误差（root mean square error,RMSE）[24]，来分别对图6中各测点温度、相对湿度、速度的模拟值与实测值进行比较，来验证模型的准确性。图8显示该工况下各测点温度的模拟值与实测值的吻合良好，绝对误差在0.04～0.38℃之间，RMSE 约为0.25℃，MRE为6.1%，ARE 为3.6%，整个培养室内的平均温度为5.81℃，各测点的平均温度约为5.85℃，符合该工况下培养室内部环境的控制要求。并且可以看出各测点相对湿度的模拟值与实测值的绝对误差在0.25%～3.36%之间，RMSE 约为1.65%，MRE 为3.8%，ARE 为1.5%，整个培养室内的平均相对湿度约为90.83%，各测点的平均相对湿度约为91.04%。

图8 温度场、相对湿度场模拟结果对比Fig.8 Comparison of simulation results of temperature field and relative humidity field

各测点的速度模拟值与实测值的对比如图9所示，模拟结果显示各测点的平均风速约为0.32m/s，绝对误差在之间0～0.19m/s。结果表明，温度场、相对湿度场与速度场的分布较为一致，结果虽存在一定的误差，但整体的精确度较高，可以证明模型有效，能够比较准确地模拟食用菌培养室的热湿耦合场的分布情况。

图9 速度场模拟结果对比Fig.9 Comparison of speed field simulation results

2.3 食用菌散湿量分析

为准确地反映散湿量对培养室环境的影响，分别计算考虑食用菌散湿和非散湿工况下室内相对湿度场的变化，其余条件均与原始模型相同，将两工况下室内相对湿度场的模拟结果进行对比，如表3所示。试验结果表明，考虑食用菌散湿量的工况下室内的相对湿度接近90%，相对湿度差仅2.11%，与实测结果较为吻合，证明该工况能够达到食用菌的生长条件。若不考虑由于食用菌的蒸腾作用所造成的散湿现象，将促使室内相对湿度梯度增大，室内相对湿度差高达16.22%，计算结果与现场实测结果出现了较大偏差，表明该工况下的培养室相对湿度环境不利于食用菌的发育。

表3 食用菌散湿对室内相对湿度的影响Table 3 Effect of edible fungus on indoor relative humidity

由于食用菌培养室为密闭的高湿空间，相对湿度过低或变化梯度过大均不利于食用菌的生长，因此在数值模拟时理应将食用菌的散湿量考虑进去，避免在对其培养室热湿耦合场进行分析时产生较大偏差。

3 改造方案数值模拟

3.1 改造方案的提出

由以上模拟结果可知，食用菌培养室内的热湿耦合场的分布并不均匀，尤其在室内的中下部分，空气流动较差，在部分区域热量和湿量出现大量的堆积，无法达到食用菌最佳的生长环境要求，而温差过大易产生培养室出菇不同步等问题。

因此，需要对该培养室的气流组织进行优化设计，本文提出如图10 所示的两种改造方案，来寻求解决上述问题的方法。

图10 各方案示意图Fig.10 Schematic diagram of each program

方案一：在送风口处安装水平向下30°的百叶，通过改变送风角度，使气流向培养室下方传送，其余条件不变。

方案二：改变送风口和回风口位置，并设置导流板，尽可能使气流均匀地到达培养室底部，其余条件不变。

3.2 方案模拟结果

3.2.1 温度场、相对湿度场的改造对比

采用和原始模型相同的计算方法，对上述两种改进方案进行数值模拟，结果如图11、图12 所示。原始模型中空气经进风口喷出，绝大多数热量被墙体吸收，与外界环境进行换热。因此方案一选择利用角度向下的百叶来改变送风角度，目的是为了减少热量损耗，使大量的冷空气直接进入菇架，来达到为食用菌生长区域降温的目的。在培养室内部，平均温度约为5.87℃，食用菌生长区域平均温度约为5.97℃，最低约5.64℃，最高约6.25℃，温差约为0.61℃。其中Y=2.6m、Y=5.6m、Y=8.6m 切面的平均温度约为5.91℃，最高约6.25℃。室内平均相对湿度为91%，生长区域最低相对湿度约为90.17%，相对湿度云图与温度云图分布一致，但由于改变了送风角度，导致房间两侧菇架的上部区域表面风速升高，厂家需调整菇架的摆放位置来避免顶部食用菌过于干燥的现象。此外，从Y=5.6m 处的计算结果可以发现，方案一仍难以解决培养室中部热量与湿量大量堆积的现象。

图11 方案一模拟结果Fig.11 Simulation results of scheme 1

图12 方案二模拟结果Fig.12 Simulation results of scheme 2

方案二中改变风机的位置并在屋顶放置导流板，其目的是为了提高冷空气到达培养室地面时的速度，便于促进室内底部空气的流动，尽可能地解决培养室中下部区域热量、湿量堆积的现象。培养室中平均温度约为5.81℃，食用菌生长区域平均温度约为5.89℃，最低约5.59℃，最高6.1℃，温差为0.51℃。其中Y=2.6m、Y=5.6m、Y=8.6m 切面的平均温度约为5.83℃，最高约6.12℃。室内平均相对湿度约为90.85%，生长区域最低相对湿度约为89.92%。

模拟结果表明，两种方案都可在一定程度上缓解原始模型中培养室底部流速过低、热量与湿量大量堆积的问题。其中方案一中食用菌生长区域温差为0.61℃，最低相对湿度为90.17%，与原始模型相比变化较小。方案二中食用菌生长区域的温差约为0.51℃，最低相对湿度约为90.3%，温度场的优化效果更佳。方案一虽无法有效保障培养室中部的温湿度分布满足金针菇的最佳生长需求，但其优点在于布置更加简便，相比方案二而言大大节省了工厂的改造成本。两种方案的设计改造均有效降低了食用菌生长区域的温差可为食用菌提供更加适宜的生长环境。

3.2.2 气流组织均匀性分析

为准确地反映出食用菌培养室内气流组织的均匀性，引入不均匀系数对气流组织的均匀性进行量化分析[25,26]。选取18 个测点来对培养室的温度、速度的均匀性进行衡量，如式（7）、（8）所示：

式中：n为测点的个数；kt、ku分别为温度不均匀系数和速度不均匀系数；、分别为n个测点的平均温度和平均速度。

显然，不均匀系数越小，说明气流组织分布的均匀性越好[27]，各方案的不均匀系数如表4所示。可以看出，改变送风角度（方案一）使温度不均匀系数降低了22.7%，培养室内的温度场分布更加均匀，但由于气流到达食用菌生长区域时速度较大，相比原始模型而言，方案一中速度场的均匀性较差。而改变送风角度并设置导流板（方案二）使室内速度不均匀系数降低了11.7%，虽可以使气流更加均匀地到达培养室底部，解决了原始模型中室内底部流速过低、温度过高的现象，但培养室内部墙体侧的温度有所升高，与原始模型相比，温度场分布较不均匀。总体而言，两方案均在不同程度上改善培养室内气流组织分布不均的现象。

表4 各方案气流不均匀性对比Table 4 Comparison of air flow unevenness of each schemes

4 结论

（1）选择标准k-ε模型，并结合特定时间食用菌的发热量与散湿量，将食用菌生长区域简化为多孔介质。利用现场实测数据以及经验数据作为边界条件，对食用菌培养室内部热湿耦合场进行数值模拟计算，温度场的平均相对误差为3.6%，相对湿度场的平均相对误差为1.5%，验证了该模型的可行性，证明食用菌培养室热湿耦合场的模拟结果具有较高的准确性。

（2）分别对考虑食用菌散湿和未考虑散湿工况下培养室内部的相对湿度进行计算，未考虑散湿时室内相对湿度差为16.22%，相对湿度变化梯度较大，模拟结果与实测结果相比出现较大偏差。考虑散湿时相对湿度差为2.11%，与实测结果较为吻合。为保证结果的准确度，模拟分析食用菌培养室热湿耦合场时理应将食用菌散湿量考虑在内。

（3）对食用菌培养室提出两种改造方案，通过CFD 数值模拟得出改造后的温度场、相对湿度场以及速度场的分布规律。根据模拟结果，方案一选择布置百叶来改变送风角度，降低食用菌生长区域的温差的效果虽有不足，但室内温度不均匀系数降低了22.7%，改善培养室内的温度场分布均匀性的效果更佳，但培养室中部的热量与湿量仍难以散去。方案二选择改变送风角度并设置导流板，速度不均匀系数降低了11.7%，食用菌生长区域的温差降至0.51℃，降低食用菌生长区域温差的效果更加明显，更有利于解决出菇不同步的问题，但相比方案一而言，其维护与改造的成本较大。