大曲發酵過程中曲房環境溫度的數值分析

白云松，田建平 *，黃海飛，王開鑄，楊海栗，黃 丹

（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%，說明曲房的對流換熱、傳熱計算具有足夠的精度。為進一步研究智能化控制曲房整體環境溫度奠定基礎。