母德強　賈洪振　姜振海

（長春工業(yè)大學(xué))

?

板式換熱器換熱區(qū)的數(shù)值計算

母德強*賈洪振姜振海

（長春工業(yè)大學(xué))

摘要根據(jù)BR0.6板式換熱器板片的實體結(jié)構(gòu)，應(yīng)用CAD軟件建立了板式換熱器換熱區(qū)三維模型；采用ANSYS軟件中的ICEM模塊將流動區(qū)域劃分為兩種類型網(wǎng)格，主流區(qū)為四面體非結(jié)構(gòu)網(wǎng)格，靠近壁面處為邊界層網(wǎng)格；選取低雷諾數(shù)κ-ε模型進行數(shù)值計算。運用Fluent軟件分析了板片熱流道內(nèi)壓力、速度和溫度的分布以及湍動能變化狀況。結(jié)果表明：出口溫度為300 K，流體流動狀態(tài)為十字交叉流，兩相鄰?fù)膭幽軜O值點間距與人字形波紋板片接觸點兩相鄰點間距在數(shù)值上相一致；低雷諾數(shù)κ-ε模型克服了壁面函數(shù)法在靠近壁面區(qū)半經(jīng)驗公式的不足。

關(guān)鍵詞板式換熱器換熱區(qū)邊界層網(wǎng)格低雷諾數(shù)κ-ε模型數(shù)值計算

*母德強，男，1961年生，教授，博士生導(dǎo)師。長春市，130012。

0　引言

板式換熱器是一種由一系列具有一定波紋形狀的金屬片疊裝而成的熱量交換裝置，具有結(jié)構(gòu)緊湊、換熱效率高等特點，因而廣泛應(yīng)用于諸多工業(yè)部門，并已成為一種重要的換熱設(shè)備。我國板式換熱器的制造企業(yè)僅有四五十家，年產(chǎn)各種規(guī)格的板式換熱器2000臺以上[1]。

但我國許多企業(yè)具有自主知識產(chǎn)權(quán)的產(chǎn)品較少，開發(fā)新產(chǎn)品時主要通過實物測試來研究板式換熱器的流體流動特性，這不僅浪費大量的物力、財力，而且也無法詳細(xì)獲得其流體內(nèi)部的特性。由傳熱基本方程

可知，影響換熱量的主要因素是與板式換熱器結(jié)構(gòu)相關(guān)的因素，包括傳熱系數(shù)K、傳熱面積A以及與傳熱單元和流動狀況相關(guān)的對數(shù)平均溫差Δtm[2]。目前，主要通過提高傳熱系數(shù)K，即采用強化傳熱的方法來增大換熱量。因此，應(yīng)用商用流體軟件進行數(shù)值分析成了換熱器設(shè)計的常規(guī)方法。查閱資料發(fā)現(xiàn)，采用RNG κ-ε模型加壁面函數(shù)法[2]對人字形板式換熱器流體特性進行數(shù)值分析者居多，但此方法在近壁面處的處理不夠精確。

隨著計算機技術(shù)的發(fā)展，網(wǎng)格劃分也越來越細(xì)化，特別是計算機計算能力的飛速提升為低雷諾數(shù)κ-ε模型的使用提供了空間。鑒于此，本文使用低雷諾數(shù)κ-ε模型，采用商用流體軟件Fluent對換熱區(qū)流場進行仿真計算，并根據(jù)計算結(jié)果確定影響換熱器性能的參數(shù)。

1　理論基礎(chǔ)

1.1數(shù)值計算所依據(jù)的數(shù)學(xué)模型

數(shù)值計算所依據(jù)的數(shù)學(xué)模型假設(shè)如下[3- 4]：

（1）工作介質(zhì)為不可壓縮的牛頓流體；

（2）重力和由于密度差異引起的浮升力忽略不計；

（3）由于換熱器內(nèi)流體流速較低，忽略流體流動時的黏性耗散所產(chǎn)生的熱效應(yīng)。

1.2基本方程

基本方程主要有連續(xù)性方程、動量方程和能量守恒方程。

（1）連續(xù)性方程

式中u、v、w——分別為x、y、z方向上的速度分量。

（2）動量方程

式中i——表示i方向；

Ui——i方向上的速度分量；

ρ——流體密度；

p——壓強分布；

μ——流體動力黏度。

（3）能量守恒方程

式中T——溫度;

α——流體擴散率。

1.3模型選取

計算流體在近壁面處的流動主要有兩種方法，即壁面函數(shù)法和低雷諾數(shù)模型[2]。用壁面函數(shù)法時，第一個節(jié)點布置在湍流區(qū)內(nèi)，通過壁面函數(shù)直接求解湍流區(qū)內(nèi)的數(shù)值。這個方法忽略了黏性底層，得不到黏性底層的精確數(shù)值結(jié)果。而低雷諾數(shù)模型則是一種直接覆蓋黏性底層的計算，在近壁面區(qū)設(shè)置了加密的邊界層網(wǎng)格，因此其計算量較大。需要強調(diào)的是，在流道內(nèi)流體處于完全湍流的狀態(tài)下，運用低雷諾數(shù)模型和標(biāo)準(zhǔn)κ-ε模型所得計算結(jié)果基本上是一致的。因此本文選用低雷諾數(shù)κ-ε模型。其中，κ方程和ε方程分別為：

式中κ——流體湍動能；

ε——湍流耗散率；

Gκ——由平均速度梯度引起的湍動能κ的產(chǎn)生項；

ρ——流體密度；

μt——在溫度t時對應(yīng)的流體動力黏度；

n——壁面法向量坐標(biāo)；實際計算時，法向坐標(biāo)n可近似取為x、y、z中任意一個；

u——與壁面平行的流速；

ui——流速在i方向的分量；

C1ε、C2ε、Cμ——經(jīng)驗常數(shù)；

σκ、σε——分別為湍動能κ和耗散率ε對應(yīng)的普朗特數(shù)（Pr）。

本文中，常數(shù)C1ε=1.44、C2ε=1.90、Cμ=0.99、σκ=1.0、σε=1.2。式（6）中，等式右邊最后一項是低Re數(shù)模型區(qū)別于高Re數(shù)模型的項。阻尼函數(shù)f1、f2和fμ的引入，實際上是對標(biāo)準(zhǔn)κ-ε模型中的系數(shù)C1ε、C2ε和Cμ進行修正。各系數(shù)計算式為[2]

f1≈1.0

Re為湍流雷諾數(shù)。

2　建立模型

換熱區(qū)是換熱器中用于熱量交換的主要部分，對換熱器的換熱效率起決定作用。本文根據(jù)BR0.6人字形板式換熱器換熱區(qū)的實體結(jié)構(gòu)（見圖1），利用PROE 5.0軟件繪制了換熱區(qū)幾何模型。其參數(shù)為：波紋高度h=3 mm，波紋間距λ=9 mm，波紋傾角β=60°。

圖1　板式換熱器換熱板實體結(jié)構(gòu)

（1）網(wǎng)格的劃分

網(wǎng)格劃分如圖2所示。

圖2　帶有邊界層網(wǎng)格的四面體網(wǎng)格

采用ANSYS軟件中的ICEM劃分網(wǎng)格。根據(jù)板式換熱器板間距較小、流道復(fù)雜多變的特點，不同流動區(qū)域使用不同類型網(wǎng)格。主流區(qū)網(wǎng)格劃分采用非結(jié)構(gòu)化四面體網(wǎng)格，靠近壁面處采用三棱柱類型的邊界層網(wǎng)格，充分保證了靠近壁面處雷諾數(shù)較低的情況或者是邊界層里面存在部分層流及轉(zhuǎn)捩狀況[4]。劃分的網(wǎng)格最大尺寸小于0.3 mm，網(wǎng)格數(shù)量為30萬個。

（2）計算方式

采用的主要分析軟件為ANSYS軟件下的Fluent模塊，計算采用分離變量隱式法求解，速度和壓力耦合采用SIMPLE方法，二階精度的迎風(fēng)格式離散[5]。

（3）邊界條件

進口采用速度入口，入口溫度為350 K。出口邊界條件為壓力出口邊界條件，設(shè)定為0.101 3 MPa[6]。外部邊界條件設(shè)為無滑移速度邊界條件u=0，溫度分布服從絕熱邊界條件，t/n｜ω=0。

3　數(shù)值模擬結(jié)果及分析

圖3為板式換熱器換熱區(qū)流道靜壓力分布，壓力降分布從右下端進口往左上端出口呈現(xiàn)階梯狀分布，壓力梯度分布比較均勻。由圖3可見，進口端（在彩圖中呈橙紅色）壓力明顯較大，出口端（藍(lán)色）壓力梯度較小，這表明人字形換熱器流體壓力損失較大。而入口段左右端壓力存在不相等情況，這是因為其入口面積不同。

圖3　換熱區(qū)流道內(nèi)壓力分布

圖4和圖5為板式換熱器換熱區(qū)流道內(nèi)速度分布和速度矢量圖。由圖4可見，流體流動狀態(tài)主要呈現(xiàn)十字交叉流，流體在板片上沿著溝槽流動，到達邊緣處被反射到另一流道內(nèi)。這種流動有助于增強湍流流動[4]，加強換熱效果，同時也使得流道內(nèi)壓力降加大。所以換熱器在考慮換熱效率的同時，也要兼顧壓力損失。

圖6為換熱器流道內(nèi)溫度分布圖，進口溫度為350 K，出口溫度為300 K。圖6中右上側(cè)為熱流體進口，左下側(cè)為出口，可以看出溫度在進口區(qū)域沿流體流動方向變化較劇烈。在板片接觸點后面出現(xiàn)高溫區(qū)域，其原因是流體在板片接觸點處發(fā)生激烈的湍流運動。由于板片接觸點的阻擋，在接觸點后面湍流顯著減小。除此以外，在進口端還存在著溫差，這是由于流體分布不均勻而引起了溫度分布不均勻。

圖4　換熱區(qū)流道內(nèi)速度分布

圖5　換熱區(qū)流道內(nèi)速度矢量分布

圖6　換熱區(qū)流道內(nèi)溫度分布

圖7清楚地顯示了湍動能的變化情況，數(shù)值大小表示湍流的激烈程度。在左側(cè)區(qū)域入口（黑色），湍動能在圖中基本不存在，但出口（綠色）的湍動能清晰地反映在圖上右側(cè)部分，并且數(shù)值都相對較小。這是因為入口設(shè)置為等流速入口，出口設(shè)置為壓力出口，而且出口不受流道形狀的限制，所以出口處湍流變化不再那么激烈。人們最關(guān)心的流體變化情況，在圖7中的分布呈現(xiàn)脈沖狀，脈沖間距比較均勻，其間距值為0.006 m左右，與兩板片接觸點的行間距相同。湍動能值反映了湍流的激烈程度，圖7中的脈沖分布處為流體在板片接觸點區(qū)域的湍流情況，此處流體發(fā)生突然轉(zhuǎn)向，湍動能值發(fā)生突變。

圖7　換熱區(qū)流道內(nèi)湍動能分布

4　模擬數(shù)據(jù)及后處理

（1）兩反向疊加的波紋板片，其兩相鄰接觸點在板片長度方向上的投影間距為：

式中L——流體進出口間的波紋板長度；

n——沿板片長度L方向任意相鄰兩行板片 接觸點個數(shù)的和（圖6中白色斑點為 板片接觸點），n=13；

h——兩相鄰接觸點沿L方向的最短矩離。計算得出h=6.15 mm，與圖7湍動能的極值點間距值基本保持一致。

（2）流體的對流換熱系數(shù)可以用相似準(zhǔn)則方程表示[7]：

式中u——流道間的流速；

de——當(dāng)量直徑；

式中b——板片有效寬度；

s——板間距，s=3 mm；

Pr——普朗特準(zhǔn)數(shù)；

——流體的運動黏度，本文中取值為

ν＝1.01×10- 6m2/s

α——流體的熱擴散率。

式中，Re的適用范圍是1280≤Re≤3830。計算結(jié)果相關(guān)性為99.5%[8]。

設(shè)置入口速度取值為0.2、0.3、0.4、0.5、0.6 m/s，得出對應(yīng)速度下的壓力降為5.12、6.96、8.36、11.52、14.95 kPa，Re值為1288、1782、2376、2970、3564。

由圖8可知，壓力降隨雷諾數(shù)的增大而增大，且變化率有增大趨勢。由此可見，在增大流速的同時，帶來的是Re和壓力降的增大。

圖8　雷諾數(shù)與壓力降的關(guān)系曲線

圖9　雷諾數(shù)與努塞爾特數(shù)的關(guān)系曲線

圖9是Re數(shù)和Nu數(shù)的變化關(guān)系曲線，Nu隨 Re的變化經(jīng)歷了從剛開始“緩慢”到后來的“較快”再到“緩慢”的過程。Re在1800~3000范圍內(nèi)，Nu的變化較敏感。而雷諾數(shù)Re的變化與流體速度又有一定的關(guān)系，通過改變?nèi)肟谒俣瓤梢愿淖兝字Z數(shù)值。因此對于特定參數(shù)的板式換熱器，設(shè)定合適的流體速度對換熱是有幫助的。

5　結(jié)論

（1）應(yīng)用ANSYS中的ICEM繪制帶有邊界層的網(wǎng)格，采用低雷諾數(shù)κ-ε模型對板式換熱器進行數(shù)值計算。該方法能對壁面區(qū)內(nèi)部進行細(xì)致的分析，尤其是在黏性底層內(nèi)，可對流體分子的黏性作用進行充分考慮。

（2）通過軟件分析和計算，得出了換熱器流體壓力場、速度場、溫度場的分布情況。發(fā)現(xiàn)流體壓力、速度分布比較均勻，其數(shù)值隨流動依次遞減；流體的雷諾數(shù)與努塞爾特數(shù)、壓力降呈線性分布，增大雷諾數(shù)，努塞爾特數(shù)和壓力降也增大；湍動能極值點間距與波紋板接觸點間距在數(shù)值上保持一致。

（3）存在的不足之處：①由于時間和實驗設(shè)備的限制，未做實驗平臺驗證。②未進行定量分析波紋高度、波紋間距和波紋傾角在不同系數(shù)下的響應(yīng)值。

參考文獻

[1]徐志明，王月明，張仲彬.板式換熱器性能的數(shù)值模擬[J] .動力工程學(xué)報，2011，31 (3) :198-202.

[2]張師帥.計算流體動力學(xué)及其應(yīng)用[M] .武漢：華中科技大學(xué)出版社，2011.

[3]曲寧.板式換熱器流動與傳熱分析[D].濟南：山東大學(xué)，2005.

[4]張冠敏.復(fù)合波紋板式換熱器強化傳熱機理及傳熱特性研究[D] .濟南：山東大學(xué)，2006.

[5]丁源，王清.Ansys Icem CFD從入門到精通[M] .北京：清華大學(xué)出版社，2013.

[6]崔立祺.人字形板式換熱器強化傳熱研究及場協(xié)同分析[D] .杭州：浙江大學(xué)，2008.

[7]程寶華，李先瑞.板式換熱器及換熱裝置技術(shù)應(yīng)用手冊[M] .北京：中國建筑工業(yè)出版社，2009.

[8]楊世銘.傳熱學(xué)[M] .北京：人民教育出版社，1980：140-148.

設(shè)計與計算

Numerical Calculation of Heat Transfer Region in Plate Heat Exchanger

Mu Deqiang Jia Hongzhen Jiang Zhenhai

Abstract：According to the entity structure of the plate in BR0.6 plate heat exchanger, the three-dimensional model of the heat transfer region is established through CAD software.The flow region is classified into two types of grids through the ICEM module of ANSYS software, among which, the main flow region is the tetrahedral unstructured grid while the region near the wall is the boundary layer grid.Meanwhile, the low Reynolds κ-ε Model is selected for the numerical calculation.The distribution of the internal pressure, velocity and temperature as well as the variability of the turbulent kinetic energy in the hot runner of the plate are analyzed via Fluent software.It's showed that the spacing between two adjacent turbulent kinetic energy extremes is numerically consistent with the value that between the two adjacent contact points of the V-shape corrugated plate when the outlet temperature is 300 K.Moreover, the defect of the semi-rational formula from the wall-function method near the wall is revised by applying the low Reynolds κ-ε Model.

Key words：Plate heat exchanger; Heat transfer region; Boundary layer grid; Low Reynolds;κ-ε model; Numerical calculation

收稿日期：（2015-06-16）

中圖分類號TQ 051.5

DOI：10.16759/j.cnki.issn.1007- 7251.2016.02.001