大曲发酵过程中曲房环境温度的数值分析

白云松，田建平 *，黄海飞，王开铸，杨海栗，黄 丹

（1.四川轻化工大学 机械工程学院，四川 宜宾 644000；2.四川轻化工大学 酿酒生物技术及应用四川省重点实验室，四川 宜宾 644000）

大曲是白酒酿造过程中的重要物质，也是必不可少的糖化剂、发酵剂和生香剂[1-2]。大曲发酵质量是决定白酒产品质量的重要环节之一，因此酿好酒必须要有好曲[3-5]。而影响大曲发酵质量的因素有很多，其中发酵时的大曲温度及水分变化和曲房温度、湿度、O2、CO2浓度是主要影响因素。目前，众多企业与学者已开展曲房环境参数在线监测与控制的相关研究，并取得了一定的成果。如李秀荣等[6-7]已设计出无线温湿度检测系统，实现了不同区域温湿度的监控。古贝春、三井等企业已将传感器物联网技术应用于大曲培曲过程的温度控制中，实现环境参数的实时监测，自动开关门窗与发送预警消息功能[8-9]。但上述成果未对环境温度变化规律及其分布进行深入研究，无法较好实现智能化控制[10]。需建立一套完整的大曲发酵曲房环境参数模型，解析大曲在上缓、中挺、后缓落不同阶段的发酵变化规律。

利用传感器物联网技术获取曲房环境温度数据的基础上，对数据进行解析。利用计算流体力学（computational fluid dynamics，CFD）仿真软件FLUENT对整个发酵曲房上缓期与中挺期的环境温度场进行模拟仿真，将结果与实测的曲房温度场云图进行对比分析，优化后的曲房环境温度仿真模型能够反映大曲发酵过程中的实际温度分布及其变化规律。通过对现有大曲发酵过程温度变化的研究，掌握其变化与微生物发酵的关系，为后续智能化调控提供理论依据，建立的仿真模型为后续控制系统提供温度场变化的数学模型奠定基础。

1 曲房温度参数检测

1.1 检测点位的布置

大曲发酵过程中，曲房内空间各点的环境温度不仅与大曲的发酵时间有关，而且与空间位置有关，还受大曲发酵的差异性及并房（堆曲的不同高度）的影响。

在实际测量中，充分考虑以上因素进行测量点的布局，如图1所示，在垂直方向上布置上、中、下三层，在水平面每层布置7个检测点位，下层7个检测点位相邻位置分别布置有7个曲内温度检测点。由于大曲在发酵中需要加盖稻草用于保温保湿，因此下层7个检测点位于稻草内。

图1 曲房检测点位分布图Fig.1 Distribution diagram of detection point location of Daqu room

1.2 检测数据结果

冬季（外界环境温度10～15 ℃）一个发酵周期（28 d）的不同点位曲内温度随时间变化的曲线图见图2。

图2 曲内温度-时间变化曲线图Fig.2 Internal temperature-time change curve of Daqu

由图2可知，大曲发酵主要分为三个阶段：上缓期，为大曲的主要发酵期，温度在不断缓慢上升；中挺期，大曲温度基本维持在55～60 ℃，温度无较大变化；后缓落期，部分大曲还在发酵，但温度开始呈现自然缓慢下降趋势。每个阶段之间都有一个较大波动段，时间分别为第6天和第12天，此时正好为操作工人对曲块进行翻曲并房。曲内温度在不同的测量点存在一定的差异性。该差异性将影响曲房内环境温度在不同空间位置的变化。

表1～表3所示为冬季大曲发酵曲房第2天（上缓期）、第5天（上缓期）以及第8天（中挺期）三个时间段某一时刻的温度检测数据。由表中数据可以看出上缓期曲房内不同位置的曲内温度差异性较大，最大温差达到9.3 ℃（上缓期两个表格的七个点位温差最大分别为9.3 ℃与6.8 ℃）。中间区域的大曲温度最高，四周温度相对较低。曲房环境温度受稻草传热影响，上层与中层点位环境温度相对于下层点位温差较小，每层各点温差在1.5 ℃范围内。

表1 第二天温度数据（上缓期）Table 1 Temperature data of the second day (slow rising period)

表2 第五天温度数据（上缓期）Table 2 Temperature data of the fifth day (slow rising period)

表3 第八天温度数据（中挺期）Table 3 Temperature data of the eighth day (intermediate maintaining period)

2 理论模型

以某酒厂曲房房间尺寸为计算、仿真模型。大曲发酵的传热理论模型包括大曲、稻草、墙壁与空气的自然对流换热控制方程，以及稻草与墙壁内部的传热控制方程。

2.1 对流传热控制方程

将空气在曲房内的流动视为三维稳态不可压缩气体的层流流动，基于Boussinesq假设的重力影响，曲房内空气在曲块与门窗的温差下做低速流动，忽略黏性耗散，密度ρ为常数，则曲房内的自然对流传热的控制方程[11-13]见式（1）～（5）：

连续方程：

动量方程：

能量方程：

其中：u、v、w分别为x、y、z三个方向的分速度，m/s；ρ是气体的密度，kg/m3；p为流体的压强，Pa；μ为流体的动力黏度，kg/（s·m）；g为重力加速度，m/s2；β=是流体的体膨胀系数，1/K；t为流体温度，K；Δt为温度差，K；是流体的热扩散率，m2/s；Cp为流体的定压比热容，kJ/（kg·K）。

2.2 固体传热控制方程

FLUENT在计算固体域的热传递时使用的能量方程形式[14-15]：

式中：ρ为密度，kg/m3；h为显焓，kJ；k为热导率，W（m·k）；T为温度，K；Sh为体积热源。

3 曲房仿真模型

3.1 曲房仿真模型的建立

由于实际检测到的曲内温度存在较大的差异性，曲房为对称结构，且检测点位均位于曲房一侧或对称线上，可将点位1、点位3、点位5和点位6对称到曲房另一侧，仿真分析中热源个数由实测点位数的7个变为11个。

采用SolidWorks三维软件建立发酵曲房三维模型。由于房间结构相对较为复杂，在保证相关物理量准确的条件下，对模型的门、窗、壁面、曲块等结构进行简化处理。所简化上缓期与中挺期的三维模型如图3、图4所示。根据检测点位的分布将曲堆简化分为11个区域。每个区域对应一个检测点位曲内温度数据。

利用ANSYS meshing进行网格划分，流体区域的网格划分方式采用四面体结构，固体区域的网格划分方式采用扫掠法。Mesh后的上缓期有限元模型有292 541个节点，149805个单元，中挺期有限元模型有402740个节点，214 296个单元，中挺期有限元模型图如图5所示。

图3 上缓期（单层曲）曲房简化模型Fig.3 Simplified model of Daqu room in the period of slow rising(single-layer Daqu)

图4 中挺期（三层曲）曲房简化模型Fig.4 Simplified model of Daqu room in the period of intermediate maintaining (three layers of Daqu)

图5 中挺期曲房有限元模型Fig.5 Finite element model of Daqu room in the period of intermediate maintaining

3.2 仿真条件的设置

利用ANSYS fluent15.0对有限元模型进行求解计算。在fluent软件中，仿真条件的设置直接影响仿真的结果。一般情况下，受到环境的影响和摆放位置的不同，不同大曲温度参数是不同的。本文选取曲块作为热源，不同区域的大曲内温度作为初始边界条件，模型属于大曲、稻草与空气的对流换热模型以及固体内部传热模型。且大曲发酵温度上升在0.1～0.3 ℃/h之间。在短时间范围内，可以将曲房内温度视为稳态，所以该仿真采用稳态仿真。

另外，本研究中影响曲房温度分布的最重要因素是曲房内稻草的导热系数。常态下，稻草的物理性质如下：密度100 kg/m3，比热容2 010 J/（kg·K），导热系数0.04 J/（m·s·K）。而在大曲发酵曲房这种特殊情况下，稻草受到湿度、温度以及自身状态的改变，其物理性质均发生了变化。在稳态情况下，稻草密度与比热容的大小对仿真结果的影响较小而导热系数的大小决定最终的仿真结果。经过多次的数值模拟，并将与实际曲房的环境温度数据比较后，不断对模型参数进行优化、修正。确定上缓期的导热系数为0.001 2～0.001 6 J/（m·s·K），中挺期的导热系数为0.001 8 J/（m·s·K）。

其他设置如下：

仿真模型启动能量设置（Energy），粘性方程（Viscous）设置采用默认层流（Laminar）。

曲房工作环境是封闭腔内自然对流换热，房间内流体介质为空气，空气的密度项设置为boussinesq，密度大小为1.2 kg/m3曲房内固体材料参数汇总见表4。

表4 曲房固体材料参数Table 4 Parameters of solid materials in Daqu room

（3）对房间内传热的热边界的设置采用流-固耦合传热边界面。在导入装配体的时候，玻璃、木门、稻草、大曲与房间内流体相接触的表面为流体与固体耦合的传热壁面，需要将其接触面的边界条件（boundary conditions）设置为耦合面（couple）。玻璃、木门与外界环境的接触面设置为对流（convection），温度为10 ℃，传热系数如表4所示。曲房的工作温度设定为10 ℃。

（4）不同区域热源温度对应表1～表3中曲内温度数据。

3.3 求解设置

求解方法选用在压力与速度的耦合（pressure-velocity coupling）算法中的Simple算法，梯度采用最小二乘单元，其他物理量采用二阶迎风格式的收敛标准，默认松弛因子不变，残差设置将能量项更改1e-5，其他保持默认不变。

4 结果对比分析

根据上述设置过程，最终曲房三个阶段仿真结果如图6～图8所示。分析图示温度分布可知，曲房环境温度分布规律基本相同，温度梯度分布均匀，温度最高区域出现在稻草下，为热源。温度最低在两侧门、窗附近，并向曲房内呈梯度升高。受到固体（稻草）传热的影响，在热源附近有一小部分区域温度梯度下降较快。受到自然对流换热影响，曲房中、上层环境温度在垂直方向由下往上呈缓慢下降趋势。

图6 第二天曲房仿真温度云图Fig.6 Simulation temperature cloud chart of Daqu room on the second day

由图6可知，曲房内稻草外环境温度小部分区域温度达到26.2 ℃，大部分区域温度在20.3～24.7 ℃。与实测最高温度22.1 ℃相差2.6 ℃，误差为11.7%，与实测最低温度20.7 ℃相差0.4 ℃，误差为1.9%。

图7 第五天曲房仿真温度云图Fig.7 Simulation temperature cloud chart of Daqu room on the fifth day

由图7可知，曲房内稻草外环境温度小部分区域温度达到35.1 ℃，大部分区域温度在26.1～32.8 ℃。与实测最高温度30.6 ℃相差2.2 ℃，误差为7.1%，与实测最低温度29.2 ℃相差3.1 ℃，误差为11.8%。

图8 第八天曲房仿真温度云图Fig.8 Simulation temperature cloud chart of Daqu room on the eighth day

由图8可知，曲房内稻草外环境温度小部分区域温度达到40.5 ℃，大部分区域温度在28.9～38.2 ℃。与实测最高温度34.2 ℃相差4 ℃，误差为11.6%，与实测最低温度31.8 ℃相差2.9 ℃，误差为9.1%。结果说明曲房的对流换热、传热计算具有足够的精度。

5 结论

对已采集到的大曲发酵曲房实时温度数据进行解析，确立了大曲发酵的三个阶段。利用已采集到的实时温度数据对上缓期和中挺期的大曲发酵过程中曲房环境温度进行数值模拟，得到曲房环境温度仿真模型。解析了曲房温度的分布规律并确定上升期与稳定期的稻草导热系数分别为0.001 2～0.001 6 J/（m·s·K）与0.001 8 J/（m·s·K）。仿真结果与实测结果的最大误差为11.8%，说明曲房的对流换热、传热计算具有足够的精度。为进一步研究智能化控制曲房整体环境温度奠定基础。