基于響應面法的仔豬配奶罐攪拌器數值模擬與優化

吳壯壯，牛智有，劉梅英，劉靜，李洪成

華中農業大學工學院/農業農村部智慧養殖技術重點實驗室，武漢 430070

哺乳階段是生豬養殖的關鍵環節，提高仔豬的成活率是決定生豬養殖產業的經濟效益和社會效益的關鍵指標之一［1］。滿足哺乳仔豬的營養攝取，及時對哺乳仔豬進行喂奶或補奶是提高仔豬成活率的重要措施。在仔豬自動喂奶設備中，仔豬配奶罐主要用于承載和制備奶水，攪拌器作為仔豬配奶罐的核心部件，主要用來促進奶水混合和加熱，其攪拌性能與奶水品質密切相關。在實際生產中，由于攪拌器型式、攪拌器結構參數等因素的不合理性導致攪拌性能低下的情況時有發生。因此，優化攪拌器參數對提高仔豬配奶罐的攪拌性能尤為重要，對提高仔豬存活率具有重要意義。

傳統的攪拌器研究如LDV或PIV測量技術等操作成本高、費時、費力，難以獲得攪拌器對流動介質混合的流場信息［2-3］。隨著計算機技術的快速發展和計算流體力學理論體系的逐步完善，CFD數值模擬為獲取攪拌槽內復雜的流場信息提供了新方法［4-5］。李欣欣等［6］針對單層雙槳葉采用CFD數值模擬的方法分析了單因素條件下不同攪拌器參數對混合時間、攪拌功率和單位體積混合能的影響，結果表明單位體積混合能可以作為參數優化指標。張智等［7］基于Fluent對不同攪拌器擺放角度下的流場進行了數值模擬分析，獲得了不同擺放角度對攪拌效果的影響規律。王定標等［8］利用CFD技術和PIV測量研究了雙層槳在不同位置的流場和濃度分布，得出了槳葉高度和加料位置對攪拌功率基本沒有影響，對混合時間有較大影響的結論。賈慧靈等［9］通過對圓盤渦輪式攪拌器不同葉片傾角下的濃度場數值模擬，發現不同葉片傾角對應的最佳安裝高度不同，且當安裝高度降低時，攪拌釜內的流型由徑向流轉化為軸向流。Ranade等［10-11］采用了CFD數值模擬方法對傳統的四斜葉渦輪和直葉渦輪進行了分析，結果表明該方法可以較好地反映攪拌槽內的湍流流場分布。綜上所述，雖然國內外研究人員對攪拌器已有了較深入的研究，但主要集中于流場的分析，而對溫度場的研究較少，尤其在綜合考慮多個因素及交互作用對流場和溫度場影響方面的研究缺失。因此，優化攪拌器參數，提高攪拌器的攪拌性能，建立多指標參數優化數學模型，具有重要的現實意義。

本研究在已有的研究基礎上，設計一種雙層槳葉仔豬配奶罐，采用CFD數值模擬與響應面分析相結合的研究方法，以雙層槳葉攪拌器的轉速、層間距、槳葉角度和離底距離為優化參數，以攪拌功率、混合時間和平均溫升速率為響應指標，基于四因素三水平正交仿真試驗計算結果，構建攪拌器參數與響應值之間的響應面回歸模型，分析各設計因素及其交互作用的影響，尋求最優的設計參數組合，旨在為仔豬配奶罐中攪拌器設計提供理論參考。

1 材料與方法

1.1 仔豬配奶罐的工作原理及總體結構

仔豬配奶罐在進行仔豬奶水制備時，需要用攪拌器對奶水進行混合。根據實際需求，設計了一種雙層槳葉式仔豬配奶罐，實際容積約為127 L，有效容積100 L，其主要由罐體、外夾套、攪拌裝置、加熱裝置等組成。在攪拌裝置中，選用雙層槳葉攪拌器，依靠其攪拌作用，加快奶粉的快速溶解和傳熱。在加熱裝置中，采用電加熱-夾套水浴加熱的加熱方式對奶水進行溫度控制。此外，為了保證奶水的質量安全，仔豬配奶罐整體采用食品級304不銹鋼材質。仔豬配奶罐結構示意圖和試驗臺架如圖1所示。

圖1 仔豬配奶罐結構示意圖（A）和試驗臺架（B）Fig.1 Structural diagram of piglet milk preparation tank （A） and test bench （B）

1.2 仔豬配奶罐數值模擬

1）物理模型構建。雙層槳葉仔豬配奶罐簡化模型如圖2所示，配奶罐直徑為500 mm，高度為600 mm，初選攪拌器直徑為350 mm，攪拌器軸徑30 mm，槳葉寬度35 mm，槳葉角度為90°，攪拌器離底距離為150 mm，攪拌器層間距為200 mm，轉速為60 r/min。

圖2 雙層槳葉攪拌器簡化模型Fig.2 Simplified model of double-layer paddle agitator

2）網格無關性分析。使用Mesh網格劃分軟件對導入的模型進行網格劃分。在網格劃分時采用結構化網格與非結構化網格相結合的方式，對形狀規則的靜止區域采用六面體結構化網格劃分，對形狀復雜的槳葉區采用空間適應能力較強的非結構化四面體網格，并對槳葉及槳葉區流體進行局部加密。為了更好地反映流場狀態以及提高數值模擬結果精度，進行網格無關性分析。對配奶罐內5個不同速度監測點進行速度分析，網格數量從5萬增加到36萬，分別進行CFD數值模擬，根據分析結果選擇28萬作為數值模擬的網格數量。

3）邊界條件設置。采用多重坐標系法（multiple reference frame， MFR）解決配奶罐內運動區域和靜止區域的交互問題，將計算域劃分為包含槳葉運動的槳葉區和除了槳葉區以外的靜止區，2個區域通過Interface面進行數據交換。配奶罐的內壁面、下表面設定為無滑移壁面，上表面設定為對稱邊界。設定計算域的流體介質為奶水，奶水密度為1 062 kg/m3，黏度為8 MPa·s，比熱容為3.5 kJ/（kg·℃），導熱系數為0.45 W/（m?K）。

1.3 模擬方法與模型驗證

1）模擬方法。通過設置不同參數模擬不同的環境進行流場和溫度場的數值分析。流場模擬采取穩態計算方法，采用標準k-ε雙方程湍流模型，以SIM?PLE算法作為壓力速度耦合方式，收斂精度設定為10-4。待穩態流場收斂后，在配奶罐內的某一點加入示蹤劑進行混合時間的模擬，將此穩態流場作為初始條件進行示蹤劑濃度場的瞬態模擬，開啟組分傳輸模型，僅激活組分傳輸項，打開示蹤劑方程，關閉其他方程，設置殘差收斂標準為10-7。溫度場模擬采用瞬態計算方法，為了簡化計算，將配奶罐內壁面和下表面設定為恒溫熱源360 K（86.85 ℃），槳葉區和靜止區初始溫度設定為300 K（26.85 ℃），能量方程殘差設置為10-8，仿真時間300 s，示蹤劑加料點及監測點和溫度監測點如圖3所示。

圖3 示蹤劑加料點及濃度監測點（A）、溫度監測點（B）Fig.3 Tracer feeding point and concentration monitoring point（A），temperature monitoring points（B）

2）仿真模型驗證。參考文獻［12］驗證方法，將功率準數模擬值與功率準數理論值進行比較。功率準數模擬值、攪拌功率P、功率準數理論值及雷諾數（Re）可由公式（1）～（7）計算得到［13-14］，力矩可由數值模擬得到，在進行理論值計算時，雙層槳葉攪拌器的功率準數可以近似看作為單層槳的功率準數乘以攪拌槳的層數［15］。經計算分析得知，功率準數理論值普遍比模擬值大，兩者相差較小，但整體變化趨勢相同，這是因為CFD數值模擬是基于流體各向同性的特性計算得到的，而實際中流體特性是各向異性的，因此，可以使用此模型進行數值模擬分析。

式（1）～（7）中：ω為角速度，rad/s；N為攪拌轉速，r/s；ρ為攪拌介質密度，kg/m3；b為槳葉寬度，mm；d為攪拌器直徑，mm；θ為槳葉角度，（°）；H為液面高度，mm；D為罐體內徑，mm；μ為攪拌介質黏度，Pa?s。

1.4 攪拌評價指標

攪拌評價指標主要包括攪拌功率、混合時間、平均溫升速率。理想狀態下攪拌功率正好為攪拌作業功率，在攪拌系統中應合理控制攪拌功率，避免過大或過小?；旌蠒r間是指從攪拌開始到罐內液體理化特性參數不存在明顯差異時的時間，國際上通常用95%的規則確定混合時間，即當一個或多個監測點達到最終穩定濃度的±5%所用的時間為混合時間。平均溫升速率是平均溫度與加熱時間的比值，能夠反映出在相同加熱時間下不同攪拌器參數對溫度的影響程度，平均溫升速率越大，加熱效率越高。

1.5 響應面法設計

本研究采用BBD法進行響應面試驗設計。在響應面分析中采用二階多項式模型，分別構建設計因素與響應指標之間的函數關系。選用的二階多項式模型基函數為

式（8）中，y為響應指標；xi、xj為設計變量；k為設計變量個數；β0為回歸方程常數；βi、βii、βij分別為回歸方程的線性偏移系數、二階偏移系數和交互系數。本研究基于雙層槳葉仔豬配奶罐，使用Design-Ex?pert8.0.6軟件以攪拌器轉速（X1）、層間距（X2）、槳葉角度（X3）和離底距離（X4）為設計優化參數，以攪拌功率（Y1）、混合時間（Y2）和平均溫升速率（Y3）為響應指標，設計四因素三水平正交仿真試驗，共設計出29組試驗。按照試驗組的攪拌器參數分別構建不同的仿真模型并展開CFD數值模擬分析。試驗因素水平如表1所示。

表1 試驗因素及水平Table 1 Test factors and levels

2 結果與分析

2.1 仿真模型構建與網格劃分

采用Design Modeler構建仿真模型，使用Mesh進行網格劃分，仿真模型及網格劃分示意圖如圖4所示。

圖4 仿真模型（A）及網格劃分（B）Fig.4 Simulation model （A） and meshing （B）

2.2 響應面模型建立與分析

1）響應面模型建立。BBD響應面試驗設計方案及數值模擬結果如表2所示，二階多項式響應面回歸模型分別為：

表2 攪拌器參數響應曲面試驗設計與結果Table 2 Experimental design and calculation results of agitator parameter response surface

2）響應面模型方差分析。為了檢驗響應面回歸模型擬合的準確性，采用方差分析法對設計因素與響應指標間的響應面模型分別進行顯著性分析，以復相關系數R2和校正決定系數R2adj評價響應面回歸模型的擬合效果，其值越接近于1，表示擬合效果越好［16］。響應面回歸模型誤差分析如表3所示，響應面方差分析如表4、表5和表6所示，可以看出，攪拌功率、混合時間和平均溫升速率的響應面模型P值均小于0.000 1，表明響應面回歸模型極顯著，具有統計學意義。攪拌功率、混合時間和平均溫升速率響應面模型的復相關系數R2分別為0.991 7、0.988 1、0.947 2，校 正 決 定 系 數R2adj分 別 為0.983 4、0.976 1、0.894 4，表明擬合得到的響應面回歸方程擬合程度高，準確性強，預測值與實際值間具有高度相關性，可以用該模型對攪拌器相關指標進行分析及預測。

表3 回歸模型誤差分析Table 3 Analysis of error in regression model

由表4可知，4個設計因素中X1、X3對攪拌功率Y1響應面模型的影響均為極顯著，交互項中X1X3影響極顯著，二次項中X12影響極顯著，X32影響顯著，其余項不顯著，設計參數影響大小順序為X1>X3>X4>X2。由表5可知，在4個設計因素中X1、X2、X3和X4對混合時間Y2響應面模型的影響均為極顯著；交互項中X2X3、X1X2和X1X4影響極顯著，X1X3影響顯著；二次項中X12、X42影響極顯著，其余項不顯著，設計因素影響大小順序為X1>X4>X2>X3。由表6可知，4個設計因素中X1、X2和X3對平均溫升速率Y3響應面模型影響顯著，交互項中均不顯著，二次項中X32、X42影響極顯著，其余項不顯著，設計因素影響大小順序為X1>X3>X2>X4。

表4 攪拌功率響應曲面二次全模型方差分析Table 4 Quadratic full model variance analysis of mixing power response surface

表5 混合時間響應面二次全模型方差分析Table 5 Mixed time response surface quadratic full mod?el analysis of variance

表6 平均溫升速率響應面二次全模型方差分析Table 6 Quadratic full model variance analysis of mean temperature rise rate response surface

3）響應面交互作用影響。由響應面模型方差分析可知，攪拌功率響應面模型中X1X3交互作用影響顯著，混合時間響應面模型中X2X3、X1X3、X1X4和X1X2交互作用影響顯著，平均溫升速率響應面模型交互作用影響均不顯著，因此在分析因素交互作用影響時著重考慮設計參數顯著交互項對攪拌功率和混合時間的影響。X1X3交互作用對攪拌功率的影響如圖5所示，由圖5可知，攪拌功率的大小主要受槳葉角度和轉速的影響，與槳葉層間距和離底距離基本沒有影響。在設計因素水平取值范圍內攪拌功率隨著槳葉角度和轉速的增大而增大，轉速對攪拌功率的影響最大。原因是隨著槳葉角度的增加，攪拌器槳葉與攪拌介質的有效接觸面積增加，攪拌介質對攪拌器槳葉的反作用力增大，引起攪拌力矩的增加，因此導致攪拌功率增大。而轉速的增加會直接引起轉矩的增加，導致攪拌功率增大。

圖5 X1X3交互作用對攪拌功率的影響Fig.5 Effect of X1X3 interaction on stirring power

由圖6A可知，當層間距一定時，混合時間隨槳葉角度的增大而增大，層間距處于較低水平時（L=100 mm），槳葉角度對混合時間的影響相對較小，最優槳葉角度范圍為30°～48°。當槳葉角度一定時，混合時間隨層間距的增大而增大，槳葉角度為90°時，層間距對混合時間的影響較大，最優的層間距范圍為100～175 mm，混合時間在槳葉角度為30°，層間距在100 mm時最小。由圖6B可知，當槳葉角度一定時，混合時間隨轉速的增大而減小，轉速處于30～70 r/min時，隨著轉速的增大，混合時間整體降低最多。轉速處于70～100 r/min時，混合時間變化趨勢趨于平緩且相對較小，最優的轉速范圍為70～100 r/min。當轉速一定時，混合時間隨槳葉角度的減小而減小，混合時間在轉速為100 r/min，槳葉角度為30°時達到最小。由圖6C可知，當離底距離一定時，混合時間隨轉速增大而減小。轉速為85～100 r/min時，轉速的增大對混合時間變化趨勢影響較小，且在此轉速范圍內混合時間隨離底距離的增大先增大后減小。由圖6D可知，當層間距一定時，混合時間隨轉速的增大而減小。轉速處于30～70 r/min時，隨著轉速的增大，混合時間整體降低最多。轉速處于72～100 r/min時，混合時間變化平穩。當轉速一定時，混合時間隨層間距的減小而減小，當轉速為100 r/min時，層間距的變化對混合時間影響最小，混合時間在轉速為100 r/min，層間距為100 mm時達到最小。

圖6 交互作用對混合時間的影響Fig.6 Effect of interaction on mixing time

4）響應面模型驗證。為了驗證轉速、層間距、槳葉角度和離底距離4個設計因素在取值范圍內變化時，反映響應指標回歸模型預測的準確性，在設計因素取值范圍內進行隨機取值，開展10組附加試驗進行預測模型的精度驗證，附加試驗組如表7所示。

分別按照表7進行仿真模型的構建，采取與表2響應面試驗組相同的操作屬性開展數值模擬分析并得到數值模擬值。將表7設計因素的取值分別代入式（9）～（11）計算得到模型預測值，模型預測值與數值模擬值對比分析結果如圖7所示。計算得知攪拌功率、混合時間和平均溫升速率的模型預測值與數值模擬值的平均相對誤差分別為9.13%、8.41%、4.08%，表明所建立的模型可以對攪拌功率、混合時間與平均溫升速率進行預測。

圖7 預測值與模擬值對比Fig.7 Comparison between predicted and simulated values

表7 附加試驗組Table 7 Additional test groups

2.3 優化結果分析

以攪拌功率最低、混合時間最低和平均溫升速率最高為優化目標，運用 Design Expert 軟件對建立的響應指標的二次全因素響應面回歸模型進行最優參數求解，約束條件為：（1）目標函數：minY1；minY2；maxY3；（2）變量區間：30≤X1≤100，100≤X2≤250，30≤X3≤90，50≤X4≤200。優 化 后 得到的各因素最優參數組合為：攪拌器轉速80 r/min，層間距170 mm，槳葉角度30°，離底距離100 mm。根據最優參數結果構建仿真模型，優化前后響應指標對比分析如表8所示，可以看出，與優化前方案相比，攪拌功率減小27.08%，混合時間減小70.15%，平均溫升速率提升9.57%。

表8 優化前后響應指標對比Table 8 Comparison of response indexes before and after optimization

為了進一步驗證優化方案與初選方案的優劣性，分析優化前后流場湍流動能云圖和溫度分布云圖。由圖8A、B可知，優化后的流場湍流強度明顯優于優化前，且湍流分布更加均勻，尤其在攪拌器的底部和近軸區湍流動能分布死區明顯減小，因此優化后的模型流場分布要優于優化前的模型。由圖8C、D可知，優化前的溫度分布在罐體上端形成了較大的低溫區，高溫區主要分布在罐體底部，而優化后的整體溫度分布均勻性明顯優于優化前，一方面是因為優化后的模型轉速的增加引起整個流體區域流速的增加，促進了流體間的熱量交換，另一方面是因為優化后的攪拌器模型槳葉形狀由平直葉型變為折葉型，平直葉型以徑向流和切向流為主，而折葉型是以徑向流、切向流和軸向流組合的混合流動，促進了近壁面溫度和罐體底部的溫度向罐體上端傳熱，使整體溫度分布更加均勻?？傮w來看，響應面優化后的攪拌器，功率損耗減少，加熱效率提高，混合強度增強，溫度分布更加均勻，更適用于仔豬配奶罐。

圖8 優化前后湍流動能分布云圖和溫度分布云圖對比Fig.8 Comparison of turbulence kinetic energy distribution and temperature distribution before and after optimization

3 討論

本研究以仔豬喂奶裝置中雙層槳葉配奶罐為研究對象，以攪拌器槳葉轉速、層間距、槳葉角度和離底距離為設計變量，攪拌功率、混合時間和平均溫升速率為響應指標，基于CFD數值模擬結果建立了設計變量與響應指標間的二次多項式預測模型，分析了各因素及其交互作用對各響應指標的影響。研究表明，轉速和槳葉角度對攪拌功率的影響極顯著，轉速、層間距、槳葉角度和離底距離對混合時間影響極顯著，轉速對平均溫升速率影響極顯著，層間距和槳葉角度對平均溫升速率影響顯著，其中轉速是最大的影響因子。攪拌功率響應面模型中轉速與槳葉角度交互作用影響顯著，混合時間響應面模型中轉速與層間距交互作用、轉速與槳葉角度交互作用、轉速與離底距離交互作用和層間距與槳葉角度交互作用影響顯著，平均溫升速率響應面模型交互作用影響均不顯著。由響應面回歸模型的誤差分析及附加試驗分析可知，基于二次多項式構建的攪拌功率、混合時間和平均溫升速率的預測模型的復相關系數R2分別為0.991 7、0.988 1、0.947 2，校正決定系數R2adj分別為0.983 4、0.976 1、0.894 4，模型預測值與數值模擬值的平均相對誤差分別為9.13%、8.41%、4.08%，表明預測模型具有較好的可靠性?；陧憫娣▋灮玫降膬灮桨概c初選方案相比，攪拌功率減小27.08%，混合時間減小70.15%，平均溫升速率提升9.57%，湍流動能分布和溫度分布相較優化前的模型有明顯改善。本研究基于響應面分析和數值模擬所得到的攪拌器參數，相較優化前攪拌性能得到明顯提升。但是目前僅對攪拌器進行了仿真優化，尚未進行試驗研究，后續研究可以對此方面進行設計補充，綜合仿真結果和試驗結果進一步優化攪拌器。