大曲發(fā)酵過(guò)程中曲房環(huán)境溫度的數(shù)值分析

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

（1.四川輕化工大學(xué) 機(jī)械工程學(xué)院，四川 宜賓 644000；2.四川輕化工大學(xué) 釀酒生物技術(shù)及應(yīng)用四川省重點(diǎn)實(shí)驗(yàn)室，四川 宜賓 644000）

大曲是白酒釀造過(guò)程中的重要物質(zhì)，也是必不可少的糖化劑、發(fā)酵劑和生香劑[1-2]。大曲發(fā)酵質(zhì)量是決定白酒產(chǎn)品質(zhì)量的重要環(huán)節(jié)之一，因此釀好酒必須要有好曲[3-5]。而影響大曲發(fā)酵質(zhì)量的因素有很多，其中發(fā)酵時(shí)的大曲溫度及水分變化和曲房溫度、濕度、O2、CO2濃度是主要影響因素。目前，眾多企業(yè)與學(xué)者已開展曲房環(huán)境參數(shù)在線監(jiān)測(cè)與控制的相關(guān)研究，并取得了一定的成果。如李秀榮等[6-7]已設(shè)計(jì)出無(wú)線溫濕度檢測(cè)系統(tǒng)，實(shí)現(xiàn)了不同區(qū)域溫濕度的監(jiān)控。古貝春、三井等企業(yè)已將傳感器物聯(lián)網(wǎng)技術(shù)應(yīng)用于大曲培曲過(guò)程的溫度控制中，實(shí)現(xiàn)環(huán)境參數(shù)的實(shí)時(shí)監(jiān)測(cè)，自動(dòng)開關(guān)門窗與發(fā)送預(yù)警消息功能[8-9]。但上述成果未對(duì)環(huán)境溫度變化規(guī)律及其分布進(jìn)行深入研究，無(wú)法較好實(shí)現(xiàn)智能化控制[10]。需建立一套完整的大曲發(fā)酵曲房環(huán)境參數(shù)模型，解析大曲在上緩、中挺、后緩落不同階段的發(fā)酵變化規(guī)律。

利用傳感器物聯(lián)網(wǎng)技術(shù)獲取曲房環(huán)境溫度數(shù)據(jù)的基礎(chǔ)上，對(duì)數(shù)據(jù)進(jìn)行解析。利用計(jì)算流體力學(xué)（computational fluid dynamics，CFD）仿真軟件FLUENT對(duì)整個(gè)發(fā)酵曲房上緩期與中挺期的環(huán)境溫度場(chǎng)進(jìn)行模擬仿真，將結(jié)果與實(shí)測(cè)的曲房溫度場(chǎng)云圖進(jìn)行對(duì)比分析，優(yōu)化后的曲房環(huán)境溫度仿真模型能夠反映大曲發(fā)酵過(guò)程中的實(shí)際溫度分布及其變化規(guī)律。通過(guò)對(duì)現(xiàn)有大曲發(fā)酵過(guò)程溫度變化的研究，掌握其變化與微生物發(fā)酵的關(guān)系，為后續(xù)智能化調(diào)控提供理論依據(jù)，建立的仿真模型為后續(xù)控制系統(tǒng)提供溫度場(chǎng)變化的數(shù)學(xué)模型奠定基礎(chǔ)。

1 曲房溫度參數(shù)檢測(cè)

1.1 檢測(cè)點(diǎn)位的布置

大曲發(fā)酵過(guò)程中，曲房?jī)?nèi)空間各點(diǎn)的環(huán)境溫度不僅與大曲的發(fā)酵時(shí)間有關(guān)，而且與空間位置有關(guān)，還受大曲發(fā)酵的差異性及并房（堆曲的不同高度）的影響。

在實(shí)際測(cè)量中，充分考慮以上因素進(jìn)行測(cè)量點(diǎn)的布局，如圖1所示，在垂直方向上布置上、中、下三層，在水平面每層布置7個(gè)檢測(cè)點(diǎn)位，下層7個(gè)檢測(cè)點(diǎn)位相鄰位置分別布置有7個(gè)曲內(nèi)溫度檢測(cè)點(diǎn)。由于大曲在發(fā)酵中需要加蓋稻草用于保溫保濕，因此下層7個(gè)檢測(cè)點(diǎn)位于稻草內(nèi)。

圖1 曲房檢測(cè)點(diǎn)位分布圖Fig.1 Distribution diagram of detection point location of Daqu room

1.2 檢測(cè)數(shù)據(jù)結(jié)果

冬季（外界環(huán)境溫度10～15 ℃）一個(gè)發(fā)酵周期（28 d）的不同點(diǎn)位曲內(nèi)溫度隨時(shí)間變化的曲線圖見圖2。

圖2 曲內(nèi)溫度-時(shí)間變化曲線圖Fig.2 Internal temperature-time change curve of Daqu

由圖2可知，大曲發(fā)酵主要分為三個(gè)階段：上緩期，為大曲的主要發(fā)酵期，溫度在不斷緩慢上升；中挺期，大曲溫度基本維持在55～60 ℃，溫度無(wú)較大變化；后緩落期，部分大曲還在發(fā)酵，但溫度開始呈現(xiàn)自然緩慢下降趨勢(shì)。每個(gè)階段之間都有一個(gè)較大波動(dòng)段，時(shí)間分別為第6天和第12天，此時(shí)正好為操作工人對(duì)曲塊進(jìn)行翻曲并房。曲內(nèi)溫度在不同的測(cè)量點(diǎn)存在一定的差異性。該差異性將影響曲房?jī)?nèi)環(huán)境溫度在不同空間位置的變化。

表1～表3所示為冬季大曲發(fā)酵曲房第2天（上緩期）、第5天（上緩期）以及第8天（中挺期）三個(gè)時(shí)間段某一時(shí)刻的溫度檢測(cè)數(shù)據(jù)。由表中數(shù)據(jù)可以看出上緩期曲房?jī)?nèi)不同位置的曲內(nèi)溫度差異性較大，最大溫差達(dá)到9.3 ℃（上緩期兩個(gè)表格的七個(gè)點(diǎn)位溫差最大分別為9.3 ℃與6.8 ℃）。中間區(qū)域的大曲溫度最高，四周溫度相對(duì)較低。曲房環(huán)境溫度受稻草傳熱影響，上層與中層點(diǎn)位環(huán)境溫度相對(duì)于下層點(diǎn)位溫差較小，每層各點(diǎn)溫差在1.5 ℃范圍內(nèi)。

表1 第二天溫度數(shù)據(jù)（上緩期）Table 1 Temperature data of the second day (slow rising period)

表2 第五天溫度數(shù)據(jù)（上緩期）Table 2 Temperature data of the fifth day (slow rising period)

表3 第八天溫度數(shù)據(jù)（中挺期）Table 3 Temperature data of the eighth day (intermediate maintaining period)

2 理論模型

以某酒廠曲房房間尺寸為計(jì)算、仿真模型。大曲發(fā)酵的傳熱理論模型包括大曲、稻草、墻壁與空氣的自然對(duì)流換熱控制方程，以及稻草與墻壁內(nèi)部的傳熱控制方程。

2.1 對(duì)流傳熱控制方程

將空氣在曲房?jī)?nèi)的流動(dòng)視為三維穩(wěn)態(tài)不可壓縮氣體的層流流動(dòng)，基于Boussinesq假設(shè)的重力影響，曲房?jī)?nèi)空氣在曲塊與門窗的溫差下做低速流動(dòng)，忽略黏性耗散，密度ρ為常數(shù)，則曲房?jī)?nèi)的自然對(duì)流傳熱的控制方程[11-13]見式（1）～（5）：

連續(xù)方程：

動(dòng)量方程：

能量方程：

其中：u、v、w分別為x、y、z三個(gè)方向的分速度，m/s；ρ是氣體的密度，kg/m3；p為流體的壓強(qiáng)，Pa；μ為流體的動(dòng)力黏度，kg/（s·m）；g為重力加速度，m/s2；β=是流體的體膨脹系數(shù)，1/K；t為流體溫度，K；Δt為溫度差，K；是流體的熱擴(kuò)散率，m2/s；Cp為流體的定壓比熱容，kJ/（kg·K）。

2.2 固體傳熱控制方程

FLUENT在計(jì)算固體域的熱傳遞時(shí)使用的能量方程形式[14-15]：

式中：ρ為密度，kg/m3；h為顯焓，kJ；k為熱導(dǎo)率，W（m·k）；T為溫度，K；Sh為體積熱源。

3 曲房仿真模型

3.1 曲房仿真模型的建立

由于實(shí)際檢測(cè)到的曲內(nèi)溫度存在較大的差異性，曲房為對(duì)稱結(jié)構(gòu)，且檢測(cè)點(diǎn)位均位于曲房一側(cè)或?qū)ΨQ線上，可將點(diǎn)位1、點(diǎn)位3、點(diǎn)位5和點(diǎn)位6對(duì)稱到曲房另一側(cè)，仿真分析中熱源個(gè)數(shù)由實(shí)測(cè)點(diǎn)位數(shù)的7個(gè)變?yōu)?1個(gè)。

采用SolidWorks三維軟件建立發(fā)酵曲房三維模型。由于房間結(jié)構(gòu)相對(duì)較為復(fù)雜，在保證相關(guān)物理量準(zhǔn)確的條件下，對(duì)模型的門、窗、壁面、曲塊等結(jié)構(gòu)進(jìn)行簡(jiǎn)化處理。所簡(jiǎn)化上緩期與中挺期的三維模型如圖3、圖4所示。根據(jù)檢測(cè)點(diǎn)位的分布將曲堆簡(jiǎn)化分為11個(gè)區(qū)域。每個(gè)區(qū)域?qū)?yīng)一個(gè)檢測(cè)點(diǎn)位曲內(nèi)溫度數(shù)據(jù)。

利用ANSYS meshing進(jìn)行網(wǎng)格劃分，流體區(qū)域的網(wǎng)格劃分方式采用四面體結(jié)構(gòu)，固體區(qū)域的網(wǎng)格劃分方式采用掃掠法。Mesh后的上緩期有限元模型有292 541個(gè)節(jié)點(diǎn)，149805個(gè)單元，中挺期有限元模型有402740個(gè)節(jié)點(diǎn)，214 296個(gè)單元，中挺期有限元模型圖如圖5所示。

圖3 上緩期（單層曲）曲房簡(jiǎn)化模型Fig.3 Simplified model of Daqu room in the period of slow rising(single-layer Daqu)

圖4 中挺期（三層曲）曲房簡(jiǎn)化模型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 仿真條件的設(shè)置

利用ANSYS fluent15.0對(duì)有限元模型進(jìn)行求解計(jì)算。在fluent軟件中，仿真條件的設(shè)置直接影響仿真的結(jié)果。一般情況下，受到環(huán)境的影響和擺放位置的不同，不同大曲溫度參數(shù)是不同的。本文選取曲塊作為熱源，不同區(qū)域的大曲內(nèi)溫度作為初始邊界條件，模型屬于大曲、稻草與空氣的對(duì)流換熱模型以及固體內(nèi)部傳熱模型。且大曲發(fā)酵溫度上升在0.1～0.3 ℃/h之間。在短時(shí)間范圍內(nèi)，可以將曲房?jī)?nèi)溫度視為穩(wěn)態(tài)，所以該仿真采用穩(wěn)態(tài)仿真。

另外，本研究中影響曲房溫度分布的最重要因素是曲房?jī)?nèi)稻草的導(dǎo)熱系數(shù)。常態(tài)下，稻草的物理性質(zhì)如下：密度100 kg/m3，比熱容2 010 J/（kg·K），導(dǎo)熱系數(shù)0.04 J/（m·s·K）。而在大曲發(fā)酵曲房這種特殊情況下，稻草受到濕度、溫度以及自身狀態(tài)的改變，其物理性質(zhì)均發(fā)生了變化。在穩(wěn)態(tài)情況下，稻草密度與比熱容的大小對(duì)仿真結(jié)果的影響較小而導(dǎo)熱系數(shù)的大小決定最終的仿真結(jié)果。經(jīng)過(guò)多次的數(shù)值模擬，并將與實(shí)際曲房的環(huán)境溫度數(shù)據(jù)比較后，不斷對(duì)模型參數(shù)進(jìn)行優(yōu)化、修正。確定上緩期的導(dǎo)熱系數(shù)為0.001 2～0.001 6 J/（m·s·K），中挺期的導(dǎo)熱系數(shù)為0.001 8 J/（m·s·K）。

其他設(shè)置如下：

仿真模型啟動(dòng)能量設(shè)置（Energy），粘性方程（Viscous）設(shè)置采用默認(rèn)層流（Laminar）。

曲房工作環(huán)境是封閉腔內(nèi)自然對(duì)流換熱，房間內(nèi)流體介質(zhì)為空氣，空氣的密度項(xiàng)設(shè)置為boussinesq，密度大小為1.2 kg/m3曲房?jī)?nèi)固體材料參數(shù)匯總見表4。

表4 曲房固體材料參數(shù)Table 4 Parameters of solid materials in Daqu room

（3）對(duì)房間內(nèi)傳熱的熱邊界的設(shè)置采用流-固耦合傳熱邊界面。在導(dǎo)入裝配體的時(shí)候，玻璃、木門、稻草、大曲與房間內(nèi)流體相接觸的表面為流體與固體耦合的傳熱壁面，需要將其接觸面的邊界條件（boundary conditions）設(shè)置為耦合面（couple）。玻璃、木門與外界環(huán)境的接觸面設(shè)置為對(duì)流（convection），溫度為10 ℃，傳熱系數(shù)如表4所示。曲房的工作溫度設(shè)定為10 ℃。

（4）不同區(qū)域熱源溫度對(duì)應(yīng)表1～表3中曲內(nèi)溫度數(shù)據(jù)。

3.3 求解設(shè)置

求解方法選用在壓力與速度的耦合（pressure-velocity coupling）算法中的Simple算法，梯度采用最小二乘單元，其他物理量采用二階迎風(fēng)格式的收斂標(biāo)準(zhǔn)，默認(rèn)松弛因子不變，殘差設(shè)置將能量項(xiàng)更改1e-5，其他保持默認(rèn)不變。

4 結(jié)果對(duì)比分析

根據(jù)上述設(shè)置過(guò)程，最終曲房三個(gè)階段仿真結(jié)果如圖6～圖8所示。分析圖示溫度分布可知，曲房環(huán)境溫度分布規(guī)律基本相同，溫度梯度分布均勻，溫度最高區(qū)域出現(xiàn)在稻草下，為熱源。溫度最低在兩側(cè)門、窗附近，并向曲房?jī)?nèi)呈梯度升高。受到固體（稻草）傳熱的影響，在熱源附近有一小部分區(qū)域溫度梯度下降較快。受到自然對(duì)流換熱影響，曲房中、上層環(huán)境溫度在垂直方向由下往上呈緩慢下降趨勢(shì)。

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

由圖6可知，曲房?jī)?nèi)稻草外環(huán)境溫度小部分區(qū)域溫度達(dá)到26.2 ℃，大部分區(qū)域溫度在20.3～24.7 ℃。與實(shí)測(cè)最高溫度22.1 ℃相差2.6 ℃，誤差為11.7%，與實(shí)測(cè)最低溫度20.7 ℃相差0.4 ℃，誤差為1.9%。

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

由圖7可知，曲房?jī)?nèi)稻草外環(huán)境溫度小部分區(qū)域溫度達(dá)到35.1 ℃，大部分區(qū)域溫度在26.1～32.8 ℃。與實(shí)測(cè)最高溫度30.6 ℃相差2.2 ℃，誤差為7.1%，與實(shí)測(cè)最低溫度29.2 ℃相差3.1 ℃，誤差為11.8%。

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

由圖8可知，曲房?jī)?nèi)稻草外環(huán)境溫度小部分區(qū)域溫度達(dá)到40.5 ℃，大部分區(qū)域溫度在28.9～38.2 ℃。與實(shí)測(cè)最高溫度34.2 ℃相差4 ℃，誤差為11.6%，與實(shí)測(cè)最低溫度31.8 ℃相差2.9 ℃，誤差為9.1%。結(jié)果說(shuō)明曲房的對(duì)流換熱、傳熱計(jì)算具有足夠的精度。

5 結(jié)論

對(duì)已采集到的大曲發(fā)酵曲房實(shí)時(shí)溫度數(shù)據(jù)進(jìn)行解析，確立了大曲發(fā)酵的三個(gè)階段。利用已采集到的實(shí)時(shí)溫度數(shù)據(jù)對(duì)上緩期和中挺期的大曲發(fā)酵過(guò)程中曲房環(huán)境溫度進(jìn)行數(shù)值模擬，得到曲房環(huán)境溫度仿真模型。解析了曲房溫度的分布規(guī)律并確定上升期與穩(wěn)定期的稻草導(dǎo)熱系數(shù)分別為0.001 2～0.001 6 J/（m·s·K）與0.001 8 J/（m·s·K）。仿真結(jié)果與實(shí)測(cè)結(jié)果的最大誤差為11.8%，說(shuō)明曲房的對(duì)流換熱、傳熱計(jì)算具有足夠的精度。為進(jìn)一步研究智能化控制曲房整體環(huán)境溫度奠定基礎(chǔ)。