基于Fluent的棉花烘房流場均勻性優化*

錢睿，吳達科，楊明金

(西南大學工程技術學院，重慶市，400715)

0 引言

干燥是使農產品可長期保存的重要手段，在傳統的氣流上升式[1]熱泵熱風干燥工藝中，熱泵將熱空氣通過管道送入烘房對農產品直接進行干燥，或進入烘房后穿過普通勻風板進行烘干，如此干燥過程結束后出現不同位置農產品干燥過度或干燥不足的情況，整體干燥品質較差，直接影響到農產品的后續儲存。因此，對此類烘房進行結構優化勢在必行。但在實際改進中，依靠經驗在現場改變烘房結構的方法并不可靠，且會造成大量的資源浪費。利用計算流體力學軟件模擬烘房內部氣流流動情況可以準確暴露內部結構存在的“缺陷”，從而有針對性地進行優化改進，將大大減少人力物力消耗。

近些年來，CFD技術在噴霧干燥領域[2-5]、熱風干燥領域[5-8]以及干燥參數優化[9-11]等方面應用發展迅速。國內外學者借助CFD技術對于農產品熱泵烘房的結構優化進行了廣泛深入的研究。Babu等[12]利用Fluent軟件比較了四種不同托盤結構葉片干燥室的干燥性能，確定了最佳優化方案，大大提高了葉片干燥時的傳熱傳質效率；王嘉麟等[13]針對固定床式花生莢果烘房設計了一種擋板組合式勻風機構，很好地改善了干燥均勻性問題；王振文等[14]利用Fluent提出在菊花烘房設置導流板可以使氣流分布更加均勻；劉瑞等[15]提出在菊花熱泵烘房進風口設置引風罩，借助于Fluent研究了引風罩最佳角度與長度，改善了整體空間的氣流均勻性。上述研究針對不同農產品或不同類型烘房的優化均有一定效果，但目前對于氣流上升式熱泵烘房內部進行結構優化，以改善內部流場均勻性的研究少之又少。

本文針對已研發的通用型氣流上升式熱泵烘房存在的農產品干燥不均勻問題，利用CFD技術研究烘房空載與滿載下的內部氣流分布情況，探尋影響熱空氣分布的結構因素，設計9組不同角度的導流葉片與不同變方孔徑勻風板的組合優化方案，探討最佳組合方式，為后續實際干燥工藝改進提供技術參考。

1 烘房建模及模擬方法

1.1 三維模型建立及網格劃分

本文所用熱泵是以熱泵烘干雙機V102-X為原型，主要由熱泵供熱系統、溫控系統、烘房熱風循環系統等組成。烘房主要由底部勻風區域、物料干燥區域、回風區域組成。熱空氣所需熱量由空氣源熱泵及輔助電加熱產生。具體烘房尺寸參數如表1所示。

表1 烘房主要參數Tab. 1 Drying room size parameters

此空氣源熱泵干燥的完整過程是：空氣在熱泵系統中經過冷凝器被加熱，通過管道進入烘房對物料進行干燥，熱空氣帶走物料中的水分，再通過回風口返回熱泵系統進行除濕加熱，如此循環。利用Solidworks所建立烘房模型以及在ICEM中網格劃分如圖1所示。

(a) 空載烘房三維圖 (b) 空載烘房網格

本文以棉花作為待烘干物料，由于實際烘房采用的是整體烘干，因此滿載時將棉花層設置為多孔介質區域，具體尺寸為1 620 mm×1 080 mm×1 267 mm。

1.2 基本假設

1) 烘房材料具有良好的保溫效果以及絕熱性能，可看作對外熱損較少，視為絕熱。

2) 實際干燥過程，流體包含空氣與水蒸氣兩種物態，但由于本文重點研究烘房內流場的分布，所以將流體看作干熱空氣。

3) 熱空氣在烘房內流動時，壓強與溫度變化小，密度變化很小，因此可視熱空氣為不可壓縮流體[16]。

1.3 數學模型

1.3.1 湍流模型

不可壓縮流體的流動滿足質量守恒方程以及動量守恒方程。湍流模型選擇標準k-ε模型，此模型即在湍動能k方程的基礎上，加入一個湍動耗散率ε的方程，形成k-ε雙方程模型[17]。

Gb-ρε-YM+Sk

(1)

(2)

(3)

(4)

(5)

式中：μt——普朗特假定湍動粘度；

μ——流體動力粘性系數，取值為1.789 4×10-5；

ui、uj——湍流速度在第i、j個方向上的分量；

Gk——平均速度梯度引起的湍動能產生；

Gb——浮力影響引起的湍動能產生；

YM——過渡的擴散所產生的波動，對于不可壓縮流體YM=0；

Sk、Sε——用戶自定義源項；

C1ε、C2ε、C3ε——經驗常數，Fluent中默認取1.44、1.922、0.09；

σk、σε——湍動能和湍動耗散率對應的普朗特數，Fluent默認取1.0、1.3；

gi——重力加速度在第i個方向上的分量；

Prt——湍動普朗特數，默認取0.85；

Cμ——經驗常數；

T——流體溫度；

β——熱膨脹系數；

ρ——流體密度，取1.225 kg/m2；

t——時間。

1.3.2 多孔介質模型

滿載時空氣在物料層之間流動會受到粘性阻力以及慣性阻力的影響，因此本文將棉花層設置為六層多孔介質區域[18-19]。在Fluent中多孔介質區域需要確定孔隙率、粘性阻力系數以及慣性阻力系數。多孔介質的模擬即在標準流體流動方程中加入了動量源項[20]。源項由粘性損失項和慣性損失項兩部分組成。

(6)

式中：Si——i方向上的(x，y，z)動量源項；

D、C——指定的系數矩陣；

|v|——速度大小。

對于簡單均質多孔介質情況，可將公式簡化為式(7)。

(7)

式中：α——滲透性系數；

C2——慣性阻力系數。

α與C2通過式(8)和式(9)求出。

(8)

(9)

Dp為多孔介質的當量直徑，本文取為棉花的平均直徑50 mm，棉花孔隙率φ通過孔隙率實驗[21]測得棉花的自身密度ρ1以及堆積密度ρ2，然后根據式(10)計算，即可得到孔隙率約為70%。

(10)

1.4 邊界條件

利用風速計采集入口9個點的風速及空氣溫度，采集數據如表2所示。

表2 入口速度溫度采集數據Tab. 2 Inlet velocity temperature collection data

求得平均速度和平均溫度分別約為10.1 m/s和42 ℃，作為進風口的邊界條件；出風口設定為壓力出口，表壓設置為101 325 Pa；進出口的湍流強度以及水力直徑通過式(11)和式(12)計算。

進出口的湍流強度

(11)

(12)

式中：Re——雷諾數；

d——圓管直徑，在本模型中進風口為矩形，d可認為等于水力直徑L。

L=4A/P

(13)

式中：A——進出風口的截面面積；

P——進出風口的截面周長。

具體邊界條件設置參數如表3所示。

表3 邊界條件參數Tab. 3 Boundary condition parameters

1.5 評價指標

1) 速度云圖以及矢量圖。在CFD仿真后處理中，通過速度云圖[22-23]可以直觀地反映流場均勻性以及分布特性，而風速過低時可引入速度矢量圖來分析流場。

2) 相對標準偏差CV。為比較流場均勻性的改善程度，引入了相對標準偏差CV[24-25]。

(14)

式中：σv——平面內所有監測點風速的標準差；

vi——各監測點風速，i=1，2，3…；

n——平面內均勻分布的監測點數目。

相對標準偏差CV越大，表示烘房內部的流場越不均勻；CV越小，則表示流場越均勻。

3) 均勻性指數γv。均勻性指數[26]是基于統計偏差定義，可全面反映截面的流體速度分布特性。γv越大，表示流體流動均勻性越好。

(15)

1.6 Fluent仿真可靠性驗證

為驗證所建立CFD烘房模型以及仿真結果的正確性與可靠性，在實際烘房中三個高度平面上按圖2所示每個平面取35個風速監測點，共105個點，利用風速傳感器測量各點風速，并記錄相對應點的仿真值，將實際測量值與仿真值進行對比，得到圖3。

圖2 測量點位置

(a) Z=0.1 m處風速仿真值與測量值 (b) Z=0.8 m處風速仿真值與測量值 (c) Z=1.5 m處風速仿真值與測量值

由圖3可知，三個平面上各35個監測點實測值與仿真值存在一定誤差，但總體變化趨勢基本一致，經計算得到實測值與仿真值的平均相對誤差約為11.7%。在實際工程應用中，要求相對誤差不應高于20%，因此本文的Fluent仿真中所建立CFD模型網格滿足模擬要求，所選取求解器模型適用此烘房模擬，模擬結果符合實際情況，可知此CFD模型是具有可靠性的，可用于后續優化仿真方案的驗證分析。

2 初始烘房模擬結果

2.1 初始空載烘房流場分析

對初始空載烘房仿真后，截取Z=0.1 m處平面的速度云圖觀察熱空氣進入烘房后在勻風區的分布情況，如圖4(a)所示，烘房下部的勻風區熱空氣主要集中在中間部分，無法向兩側完全散開，呈現出勻風區兩側的風速明顯低于中間風速，并在兩側形成對稱的幾處較大渦旋，初步分析是由于進風口寬度有限所致，這些渦旋會導致對應位置的烘干區部分區域空氣流速很低，但本文中熱泵的送風口寬度不可改變，所以需采取方案優化烘房進風口。

(a) Z=0.1 m處速度云圖

由于內部空間較大，具有一定初始速度的熱空氣進入烘房后速度逐漸減弱并雜亂分散，受兩側渦旋影響，在中間位置有一小段熱空氣流速突然增大。為觀察延X軸方向實際干燥區前后風場分布情況，截取了Y=0.6 m處速度云圖，如圖4(b)所示，具有一定壓強的熱空氣進入烘房后直接沖擊至烘房后壁面，導致很大一部分熱空氣從鄰近后壁面的圓孔進入烘干區，在靠近進風口處，通過勻風板進入干燥區的熱空氣較少，后壁面附近風量明顯多于進風口側，并且由于進風口與出風口在同一水平面，導致不可避免地出現很大的渦旋。

另外，截取烘干區Z方向上的五個截面，并在各個截面選取20個風速監測點，處理得到初始結構下的烘干區五個截面的平均速度、相對標準偏差CV以及均勻性系數 的具體數值如表4所示，烘干區從下至上平均速度呈上升趨勢，越靠近出風口風速越高，且各截面的相對標準偏差CV都比較高，均勻性系數均處于較低水平，流場整體均勻性較差。

表4 初始結構烘干區XY面流場評價指標Tab. 4 Evaluation index of XY surface flow field in drying zone of initial structure

2.2 初始滿載烘房流場分析

為觀察滿載時烘干區速度場情況，截取了XZ平面上Y=0.6 m處截面，得到圖5所示速度云圖。

(a) 滿載Y=0.6 m速度云圖

由于多孔介質層之間的粘性阻力及慣性阻力的存在，導致風速過低，無法直觀地觀察到物料層之間的風場情況，因此引入速度矢量圖來代替速度云圖。從速度矢量圖可以看出，具有一定初速度的熱空氣進入烘房后，從勻風板前后部進入烘干區熱空氣較多，由烘房中間部分進入物料層的熱空氣較少，這是由于物料層與烘房壁面間存在空隙，熱空氣沿空隙流向出風口。物料層由下至上，熱空氣始終受粘性阻力及慣性阻力影響，流動逐漸減少，有較多速度“缺口”，在物料層右上角出現風速“空白區”，是由于烘房內外壓強差的存在導致附近熱空氣直接流向出風口。

3 烘房結構優化設計

為改善烘房由于進風口寬度有限導致的勻風區兩側風場不均的問題，在進風口處設置圖6所示葉片用以導流，將進來的熱空氣充分擴散，減少勻風區渦旋情況。為使仿真快速有效地進行，在確定仿真方案之前，進行了葉片角度從15°～60°以15°為梯度差的預仿真，觀察四種角度葉片對流場所造成的影響，結果顯示：超過30°的葉片對熱空氣分散效果較差，因此最終選取了角度的三個梯度進行仿真分析，分別為10°、15°、20°。

圖6 進風口優化設計

針對初始烘房延X軸方向烘干區前后風場不均的問題，并為提高進入烘干區的熱空氣流速，將原有勻風板上的圓孔改為變徑方孔，如圖7所示。

(a) A類勻風板 (b) B類勻風板

由于變徑方孔的存在，接近后壁面處勻風板上的方孔徑較小，上升至烘干區的熱空氣隨之變少，下方壓強增大，促使熱空氣向勻風板中前部移動，使得勻風區的熱空氣在三個區域內能較均勻地上升至烘干區。從前到后，A類勻風板孔徑為45 mm，35 mm，25 mm；B類勻風板孔徑為40 mm，30 mm，20 mm；C類勻風板孔徑為35 mm，25 mm，15 mm。

將上述的兩種優化方法中進行組合，設計出9種方案，如表5所示。

表5 仿真方案Tab. 5 Simulation schemes

4 結構優化模擬結果

4.1 優化后空載烘房流場分析

為比較優化前后烘房內部橫向風場變化，分別在原始結構以及9種優化方案得到的仿真結果中截取了XZ平面上Y=0.6 m處的速度云圖，并在此截面烘干區范圍內選取24個風速監測點，計算出各優化方案在此截面上的平均速度、相對標準偏差以及均勻性系數，如圖8所示。

9種優化方案在XZ平面上較原始結構均大幅度提高了平均速度；除20°葉片下的3種優化方案外，另外6種方案也均大幅度降低了相對標準偏差并提高了均勻性系數。

為反映各方案烘干區不同高度平面的氣流狀態，截取了Z=0.6、0.85、1.1、1.35、1.6 m這5個平面的速度云圖。在這5個平面上各選取了20個風速監測點，經處理后得到平均速度以及相對標準偏差的點線圖如圖9所示。

從圖9中不同優化方案與原始狀態下的截面平均速度和相對標準偏差對比可以看出，所有優化方案在XY平面上均不同程度提高了平均速度并降低了相對標準偏差。圖9(a)中，在葉片10°的情況下，3類變孔徑勻風板在提高整體平均速度方面差異不是很大，但圖9(d)中B類勻風板的相對標準偏差整體處于下降趨勢，且大幅度降低，而勻風板C與A的相對標準偏差在不同截面的波動較大，因此在10°葉片組中最優選為B勻風板。同理，在葉片15°以及20°的情況下，分別是勻風板B和勻風板C優化效果較同組其他方案效果更優。

(a) 9種方案Y=0.6 m處平均風速 (b) 9種方案Y=0.6 m處相對標準偏差 (c) 9種方案Y=0.6 m處均勻性系數

(a) 10°葉片下各方案平均速度 (b) 15°葉片下各方案平均速度 (c) 20°葉片下各方案平均速度

(d) 10°葉片下各方案相對標準偏差 (e) 15°葉片下各方案相對標準偏差 (f) 20°葉片下各方案相對標準偏差

4.2 確定最優方案

將上述各葉片下最優方案進一步對比，如圖10所示。10°葉片與勻風板B在XY截面上的平均速度比另外兩種組合較原始結構高出更多，但綜合相對標準偏差，15°葉片和勻風板B一直處于下降趨勢，比較穩定。

(a) 各組最優方案平均速度對比 (b) 各組最優方案相對標準偏差對比 (c) 各組最優方案均勻性系數對比

將所監測風速點進一步處理得到了各組最優方案的均勻性系數對比，可知15°葉片和勻風板B的組合下，各平面的均勻性系數均處于較好狀態，且隨著高度的上升均勻性逐漸增加。綜合上述各評價指標，最終選定15°葉片和勻風板B為最優組合。

表6是優化后烘干區XY不同高度平面的平均速度、相對標準偏差CV以及均勻性系數的具體數據，與表4對比，各平面平均速度提高了0.231 m/s，平均相對標準偏差由60.51%降至38.87%，平均均勻性系數由0.76提升至0.85，烘干區由下至上，風速逐漸增大，風場愈發均勻。

表6 優化后烘干區XY平面流場評價指標Tab. 6 Evaluation index of XY plane flow field in drying zone after optimization

4.3 最優方案下滿載烘房流場分析

優化前后對比如圖11、圖12所示。

(a) 優化前Y=0.35 m監測點風速 (b) 優化前Y=0.6 m監測點風速 (c) 優化前Y=0.85 m監測點風速

(d) 優化后Y=0.35 m監測點風速 (e) 優化后Y=0.6 m監測點風速 (f) 優化后Y=0.85 m監測點風速

(a) 優化前Z=0.5 m監測點風速 (b) 優化前Z=1.0 m監測點風速 (c) 優化前Z=1.5 m監測點風速

(d) 優化后Z=0.5 m監測點風速 (e) 優化后Z=1.0 m監測點風速 (f) 優化后Z=1.5 m監測點風速

確定優化方案后，加入多孔介質并進行模擬仿真。為比較優化前后物料層內XZ不同位置平面的風場變化，使用ISO CLIP工具選定了Y=0.35、0.6、0.85 m三平面上處于物料層內部的部分截面，并各沿進風方向以及垂直進風方向選取了300個風速記錄點，分別記為x1,x2。由圖11可知，優化前后三個物料層內部截面的風場變化趨勢相似；另外，原始烘房下在靠近進風口以及出風口側風速過高，而其他位置的速度偏低，平均速度僅為0.087 m/s，優化后提升至0.2 m/s；原始烘房下的三維曲面圖波動較大，三截面的平均風速相對標準偏差為44.8%，優化后降至25.1%，平均均勻性系數由0.84提升至0.91。

為了從烘干區XY不同高度平面上分析優化前后流場變化，同樣使用ISO CLIP工具選定了Z=0.5、1.0、1.5 m三平面上處于物料層內部的部分截面，各取了300個風速記錄點后繪制三維曲面圖如圖12所示。原始烘房下內部中間位置熱空氣流速很低，平均速度僅0.096 m/s，優化后被提高至0.23 m/s；另外，原始烘房在Z=0.5 m和1.5 m處速度波動很大，是因為這兩處平面靠近進出風口側，而Z=1.0 m由于處于物料層中部，其上各點風速極低，變化范圍小，而優化后的烘房風場變化顯然更加平穩，經結構優化后平均相對標準偏差由30.1%降至24.32%，且優化后三平面的速度均勻性系數均在0.9以上，平均均勻性系數由0.88提升至0.94。

5 結論

本文基于Fluent軟件對氣流上升式熱泵烘房流場進行了模擬仿真，在驗證了仿真模型的可靠性與正確性的基礎下完成了結構優化。

1) 通過在進風口設置不同角度導流葉片，并且將勻風板上原有的統一直徑圓孔改為變徑方孔，共設計了9種優化方案組合。通過比較速度云圖、速度矢量圖、平均速度、相對標準偏差以及均勻性系數，最終確定15°導流葉片與B類勻風板的組合為最優組合。

2) 優化后的烘房相較于初始烘房，烘干區平均速度提高了0.231 m/s，烘干區XY不同高度平面的平均相對標準偏差CV由60.51%降至38.87%，平均均勻性系數由0.76提升至0.85；滿載時烘干區物料層范圍內XZ不同位置平面的平均速度0.087 m/s增強至0.2 m/s，平均相對標準偏差CV由44.8%降至25.1%，平均均勻性系數由0.84提升至0.91；烘干區物料層范圍內XY不同高度平面的平均速度由0.096 m/s增強至0.23 m/s，平均相對標準偏差CV由30.1%降至24.32%，平均均勻性系數由0.88提升至0.94，即優化后的烘房流場整體低速區域大大減少，流場均勻性有效得到改善，將更有利于實際烘干，從而提高產品的最終干燥效果。本研究可為氣流上升式熱泵烘房的優化設計提供理論方法及技術方案。