焦文靜，戴傳山，王秋香

（天津大學(xué)機(jī)械學(xué)院，中低溫?zé)崮芨咝Ю媒逃恐攸c(diǎn)實(shí)驗(yàn)室，天津300072）

換熱器是將熱流體的部分熱量傳遞給冷流體的設(shè)備，廣泛地應(yīng)用于石油、化工、冶金、動(dòng)力、核能等工業(yè)部門及國防工程中。換熱器作為化工過程中必不可少的單元設(shè)備，其材料及動(dòng)力消耗占整個(gè)工藝設(shè)備消耗的30%左右［1］。換熱器的設(shè)計(jì)水平往往直接影響著工業(yè)的投資效益和能效［2］，使得換熱器的強(qiáng)化傳熱技術(shù)越來越受到重視。如何提高換熱器的性能，一直是換熱器研究和開發(fā)的目標(biāo)［3－4］。換熱器微型化，不僅僅意味著節(jié)約材料和空間，更重要的是與常規(guī)尺度的換熱器相比，實(shí)現(xiàn)相同目標(biāo)的同時(shí)，通常情況下可以大幅度的縮減所消耗的能量，具有高傳熱系數(shù)、高表面積體積比、低傳熱溫差、低流動(dòng)阻力等優(yōu)點(diǎn)，與普通換熱器相比，其單位體積傳熱系數(shù)高達(dá)幾十到幾百M(fèi)W／（m3·K），比普通換熱器要高12個(gè)數(shù)量級(jí)。器件裝置微型化的強(qiáng)大發(fā)展趨勢推動(dòng)了微型換熱器技術(shù)的迅猛發(fā)展，使其成為熱能傳遞與動(dòng)力轉(zhuǎn)換領(lǐng)域的一個(gè)重要研究方向［5］。微尺度換熱器的特征尺度在微米到亞毫米尺度范圍內(nèi)，微尺度傳熱與流動(dòng)問題機(jī)理非常復(fù)雜，需要對(duì)其進(jìn)行大量的理論模擬與實(shí)驗(yàn)研究。格子波爾茲曼方法（LBM）是20世紀(jì)80年代新興發(fā)展起來的一種數(shù)值計(jì)算方法，具有粒子特性的先天優(yōu)勢，相比其它常規(guī)的數(shù)值方法更容易處理復(fù)雜的理論模擬問題［6－7］。在LBM數(shù)值模擬中，大多采用均勻網(wǎng)格［8－10］對(duì)模型進(jìn)行剖分，這雖然有利于LBM的實(shí)現(xiàn)，但在計(jì)算精度和計(jì)算代價(jià)的問題上限制了其應(yīng)用范圍。均勻網(wǎng)格模擬的收斂速度受網(wǎng)格步長的制約：網(wǎng)格越密，精度相對(duì)較高但收斂速度越慢。為了解決這個(gè)問題，多格子波爾茲曼方法被學(xué)者提了出來［11］。對(duì)復(fù)雜流動(dòng)問題，在靠近流場邊界的區(qū)域其物理量一般變化復(fù)雜、劇烈，對(duì)這些區(qū)域應(yīng)該采用較細(xì)的網(wǎng)格來保證計(jì)算精度；而在遠(yuǎn)離邊界的區(qū)域，物理量變化相對(duì)平緩，采用較粗的網(wǎng)格可以在不影響精度的前提下大大提高收斂速度。所以，探討根據(jù)物理量變化的劇烈程度劃分成粗細(xì)不同的網(wǎng)格，對(duì)實(shí)現(xiàn)高效的LBM數(shù)值模擬計(jì)算有著重要的現(xiàn)實(shí)意義。

1 多格子LBM

這種方法最先由Yu［11］提出。多格子LBM對(duì)處理復(fù)雜邊界非常有效，能夠較好的滿足精度和高效的雙重要求。這種方法是將計(jì)算區(qū)域分為不同的塊，各個(gè)塊之間通過邊界相互交換信息，交換的信息必須滿足一定的關(guān)系，從而使質(zhì)量、動(dòng)量守恒，使作用力能連續(xù)的通過邊界。

通過兩塊系統(tǒng)來說明多格子LBM方法［12］。兩種格子的空間步長比率m定義如式 （1）。

Δxc、Δxf分別是粗網(wǎng)格和細(xì)網(wǎng)格的空間步長。流體的黏度由式 （2）求出。

為了滿足黏度為常數(shù)，粗網(wǎng)格的松弛時(shí)間τc和細(xì)網(wǎng)格的松弛時(shí)間τf必須滿足式 （3）。

同時(shí)為了保證各變量在通過粗、細(xì)網(wǎng)格邊界時(shí)滿足連續(xù)性，粗、細(xì)網(wǎng)格的密度分布函數(shù)在碰撞后必須滿足式 （4）、式 （5）。

圖1 粗、細(xì)網(wǎng)格邊界結(jié)構(gòu)圖

從圖中可以看出，粗網(wǎng)格與細(xì)網(wǎng)格在邊界上互相交錯(cuò)。值得注意的是，在細(xì)網(wǎng)格邊界MN上，黑點(diǎn)標(biāo)記的格子點(diǎn)是沒有原始信息的，由白點(diǎn)標(biāo)記的格子通過采用3次樣條插值的空間差分得到式（6）。

想要進(jìn)行合理準(zhǔn)確的模擬計(jì)算，粗、細(xì)網(wǎng)格的梳理速度必須相同，也就是說細(xì)網(wǎng)格進(jìn)行m步梳理，粗網(wǎng)格只進(jìn)行一步梳理，因此在細(xì)網(wǎng)格邊界MN上需要用時(shí)間插值方法從而得出細(xì)網(wǎng)格在f（t（n＋1／m），MN）時(shí)的值。同樣我們采用三點(diǎn)拉格朗日插值法得到式 （7）。

對(duì)于速度場與溫度場的模擬，這些處理方法同樣適用。

2 驗(yàn)證程序正確性

有自然對(duì)流換熱的模擬較為復(fù)雜，對(duì)于自編程序是否能準(zhǔn)確地研究含有自然對(duì)流的換熱情況，需要驗(yàn)證程序的正確性。所以，可以通過模擬正方形空間內(nèi)單熱圓管自然對(duì)流這一經(jīng)典算例并與前人結(jié)果進(jìn)行對(duì)比以驗(yàn)證程序的正確性。作者模擬了瑞利數(shù)Ra＝1×104，Ra＝5×104和Ra＝1×1053種情況下熱管表面平均努塞爾數(shù)、流線、等溫線并與已有文獻(xiàn)［13－15］進(jìn)行了對(duì)比。

圖2是熱管在方形空間內(nèi)的換熱的多格子LBM模擬理論模型。正方形計(jì)算區(qū)域采用粗細(xì)兩種網(wǎng)格，藍(lán)色線代表粗網(wǎng)格，A－B－C－D是粗網(wǎng)格的內(nèi)邊界，紅色線代表細(xì)網(wǎng)格，a－b－c－d是細(xì)網(wǎng)格的外邊界。計(jì)算區(qū)域?yàn)?.5×2.5，劃分成50×50的粗網(wǎng)格，步長Δxc＝Δyc為0.05以及內(nèi)部50×50細(xì)網(wǎng)格步長Δxf＝Δyf為0.025，圓管直徑為1。對(duì)于速度場，粗細(xì)網(wǎng)格計(jì)算區(qū)域中所有節(jié)點(diǎn)初始速度為0，四周均采用無滑移邊界條件。對(duì)于溫度場，圓管為恒壁溫邊界條件，量綱為1溫度為1，四周采用恒壁溫邊界條件，量綱為1溫度為0，流場中各節(jié)點(diǎn)的初始無量綱溫度均為0。模型中熱管的直徑設(shè)為d，正方形區(qū)域邊長定義為H，同時(shí)定義d／H為管徑空間比，在這里取d／H＝0.4，為了方便與已有文獻(xiàn)做比較，瑞利數(shù)Ra定義如式 （8），努塞爾數(shù)Nu定義如式（9）、式（10）。

圖2 單熱管在方形空間內(nèi)換熱的多格子LBM模擬理論模型

瑞利數(shù)Ra

局部Nu數(shù)

表面平均Nu數(shù)

圖3給出了Ra＝104情況下溫度場和流場的模擬結(jié)果對(duì)比，圖3（a）、（b）是Peng［12］等模擬的溫度場和流場，圖3（c）、（d）是王秋香等［14－15］采用均勻網(wǎng)格的模擬結(jié)果，圖3（e）是采用多格子LBM的模擬結(jié)果，值得一提的是，計(jì)算時(shí)間相比均勻網(wǎng)格模擬節(jié)約了30%40%。

表1給出了平均努塞爾數(shù)Nu模擬結(jié)果的對(duì)比，本文模擬結(jié)果與Peng等［13］模擬結(jié)果的相對(duì)誤差在0.2%2.6%，與王秋香等均勻網(wǎng)格的模擬結(jié)果誤差在0.5%1.1%，相對(duì)誤差在可接受范圍內(nèi)，充分說明了多格子LBM在大大節(jié)約計(jì)算時(shí)間的同時(shí)能夠保證計(jì)算的準(zhǔn)確性。

圖3 溫度場和流場的模擬結(jié)果對(duì)比

表1 平均努塞爾數(shù)Nu模擬結(jié)果的對(duì)比

3 冷熱微圓管自然對(duì)流理論模型

圖4 冷熱微圓管自然對(duì)流換熱的多格子LBM模擬理論模型

圖4是冷熱微圓管自然對(duì)流換熱的多格子LBM模擬理論模型。在模擬過程中設(shè)定雷諾數(shù)Re為0，研究純自然對(duì)流情況下（具體是Ra＝103時(shí)的換熱情況），冷熱管中心連線與浮升力呈不同夾角對(duì)整體傳熱性能的影響。

在整個(gè)模擬計(jì)算中，冷熱微圓管直徑相同，量綱為1長度為1，計(jì)算區(qū)域20×20，劃分成外圍200×200的粗網(wǎng)格，步長Δxc＝Δyc為0.1，以及內(nèi)部200×200細(xì)網(wǎng)格步長Δxf＝Δyf為0.05。對(duì)應(yīng)于真實(shí)模型中：冷熱微圓管外管徑為1.0 mm，管間距為1.5d，也就是1.5 mm。速度場中，粗細(xì)網(wǎng)格計(jì)算區(qū)域中所有節(jié)點(diǎn)初始速度為0，四周采用無滑移邊界條件。溫度場中，冷熱圓管為恒壁溫邊界條件，量綱為1，溫度為1和 2，四周為恒壁溫邊界條件，量綱為1，溫度為1.5，流場中各節(jié)點(diǎn)的初始溫度均為量綱為1，溫度1.5。值得注意的是，與驗(yàn)證程序不同，這里以圓管的直徑d為特征尺寸來定義瑞利數(shù)Ra以及微管表面平均Nu。

4 模擬結(jié)果與分析

圖5（a）給出了冷熱管中心連線與浮升力垂直時(shí)的溫度場流場的分布圖。從圖5（a）中可以看出，這種布管方式下，浮升力使被加熱流體迅速上升，被冷卻流體快速下沉，形成穩(wěn)定的以冷熱管中心連線為軸的上下對(duì)稱的溫度場分布。冷熱管之間互有影響，熱管的溫度場向右傾斜 （冷管方向），冷管的溫度場向左傾斜（熱管方向）。

從流線分布可以看出，貼著熱管壁面被加熱的流體迅速上升，至上壁面分成兩部分，一部分順時(shí)針沿右壁面下行，至冷熱管中心線水平位置后，流向冷管；而另一部分逆時(shí)針沿左壁下降，至冷熱管中心水平位置后，流向熱管。圍繞冷管的流線與熱管類似，但方向相反，基本上為中心對(duì)稱。

圖5（b）給出了冷熱管中心連線與浮升力成45°的溫度場流場分布圖。從圖5（b）中可以看出，這種布管方式下，浮升力使被加熱流體迅速上升，被冷卻流體快速下沉，形成穩(wěn)定的以冷熱管連線中心為原點(diǎn)的對(duì)稱的溫度場分布。相比冷熱管中心連線與浮升力垂直布管方式，冷熱流體分別在靠近底角點(diǎn)和頂角點(diǎn)聚集，換熱效果略差于前者。

圖5 冷熱管不同布管方式的溫度場流場

圖6 不同布管方式下熱管、冷管的局部Nu圖

從流線分布可以看出，貼著熱管壁面被加熱的流體迅速上升，至上角點(diǎn)附近分成兩股，一路繼續(xù)順時(shí)針沿右上壁面下行，在水平方向繞回構(gòu)成回路。另一股逆時(shí)針方向流動(dòng)沿左上壁面下行，水平方向繞回構(gòu)成回路。而冷管周圍被冷卻的流體略偏左向下流動(dòng)，在底角點(diǎn)分成兩股，一路逆時(shí)針沿有下壁面上行，行至在水平方向受來流熱流體阻礙，沿水平方向流回同時(shí)與熱流體進(jìn)行換熱。另一路沿左下壁面順時(shí)針上行，同樣行至在水平方向受來流熱流體阻礙，沿水平方向流回同時(shí)與熱流體進(jìn)行換熱。這種布管方式被加熱的流體與被冷卻的流體之間沒有進(jìn)行充分的摻混，換熱不充分。

圖5（c）、圖5（d）分別是冷熱管中心連線與浮升力平行熱上冷下以及熱下冷上的溫度場流場對(duì)比圖。熱管在上冷管在下時(shí)，形成以冷熱管連線中心上下左右對(duì)稱的溫度場流場分布。浮升力使被加熱的向上流動(dòng)，被冷卻的向下匯集。

熱管在下冷管在上時(shí)，熱管周圍被加熱的流體在浮升力作用下向上運(yùn)動(dòng)，但是在這過程中，受到正上方冷管的阻擋，必須要繞過冷管邊緣才能繼續(xù)上升。同樣，冷管周圍被冷卻的流體在向下運(yùn)動(dòng)的過程中受到熱管的阻礙，也必須繞過熱管壁面才能繼續(xù)下行。這就增大了被加熱被冷卻流體之間換熱的阻力，使得整體換熱效果較差。圖6給出了不同布管方式下熱管、冷管表面的局部Nu圖。

表2給出了不同布管方式下冷熱管表面平均Nu以及冷熱管整體平均Nu的對(duì)比。結(jié)果表明，冷熱管連線與浮升力垂直時(shí)冷熱管的平均Nu最大，表示換熱最好，熱下冷上布管換熱最差。

表2 不同布管方式下平均Nu對(duì)比

5 結(jié) 論

（1）冷熱管中心連線與浮升力垂直的布管方式，浮升力使被加熱流體迅速上升，被冷卻流體快速下沉，形成穩(wěn)定的以冷熱管中心連線為軸的上下對(duì)稱的溫度場分布。冷熱流體之間換熱充分，整體換熱效果最好。冷熱管中心連線與浮升力成45°布管方式的溫度場流場相比冷熱管中心連線與浮升力垂直布管方式的溫度場流場，冷熱流體主要在靠近底角點(diǎn)和頂角點(diǎn)聚集，換熱略差于前者。

（2）熱管在上冷管在下的布管方式，形成以冷熱管連線中心上下左右對(duì)稱的溫度場流場分布。浮升力使被加熱流體向上流動(dòng)，被冷卻的流體向下匯集。冷熱流體之間換熱較好。

（3）熱管在下冷管在上的布管方式，熱管周圍被加熱的流體在浮升力作用下向上運(yùn)動(dòng)，在這過程中，受到正上方冷管的阻擋，必須要繞過冷管邊緣才能繼續(xù)上升。同樣，冷管周圍被冷卻的流體在向下運(yùn)動(dòng)的過程中受到熱管的阻礙，也必須繞過熱管壁面才能繼續(xù)下行。這就增大了被加熱被冷卻流體之間的換熱阻力，使得整體換熱效果最差。

符 號(hào) 說 明

m——空間步長比率，量綱為1

Δx——網(wǎng)格空間步長，量綱為1

τ——網(wǎng)格的松弛時(shí)間，量綱為1

v——流體的黏度，量綱為1

Ra——瑞利數(shù)，量綱為1

f——波爾茲曼密度分布函數(shù)

d／H——為管徑空間比，量綱為1

β——體積膨脹系數(shù)，K－1

α——熱擴(kuò)散系數(shù)，m2／s

Nu——努塞爾數(shù)，量綱為1下角標(biāo)

f——細(xì)網(wǎng)格

c——粗網(wǎng)格

［1］ 蘇尚美，張亞男，成方園，等.微通道換熱器的特性分析及其應(yīng)用前景［J］.區(qū)域供熱，2007，5：34－38.

［2］ 華賁，仵浩，劉二恒.基于火用經(jīng)濟(jì)評(píng)價(jià)的換熱器最優(yōu)傳熱溫差化［J］.化工進(jìn)展，2009，28（7）：1142－1146.

［3］ 支浩，湯慧萍，朱紀(jì)磊.換熱器的研究發(fā)展現(xiàn)狀 ［J］.化工進(jìn)展，2009，28（s）：338－342.

［4］ 陸應(yīng)生，陳墓玲，潘寧忠，等.強(qiáng)化傳熱元件與高效換熱器研究進(jìn)展［J］.化工進(jìn)展，1998，1：46－48.

［5］ 戴傳山，王秋香，孫平樂.污垢對(duì)微管換熱器用于發(fā)電系統(tǒng)的影響［J］.化工進(jìn)展，2009，28：371－375.

［6］Yan Y Y，Zu Y Q.Numerical simulation of heat transfer and fluid flow past a rotating isothermal cylinder－A LBM approach［J］.International Journal of Heat and MassTransfer，2008，51：2519－2536.

［7］Shu C，Liu N，Chew Y T.A novel immersed boundary velocity correction－lattice Boltzmann method and its application to simulate flow past a circular cylinder［J］.Journal of Computational Physics，2007，226：1607－1622.

［8］Lima E S，Silveira N A，Damasceno J.Numerical simulation of two－dimensional flows over a circular cylinder using the immersed boundary method［J］.Comput.Phys.，2003，189：351.

［9］Feng Z G，Michaelides E.The immersed boundary－lattice Bol－tzmann method for solving fluid－particles interaction problems［J］.Comput.Phys.，2004，195：602.

［10］Feng Z G，Michaelides E.Proteus：a direct forcing method in the simulations of particulate flows［J］.Comput.Phys.，2005，202：20.

［11］Sui Y，Chew Y T，Roy P，et al.A hybrid immersedboundary and multi－block lattice Boltzmann method for simulating fluid and moving－boundaries interactions［J］.Int.Numer.Meth.Fluids，2007，53：1727－1754.

［12］ 葛鑫.微管外流動(dòng)換熱的玻爾茲曼數(shù)值模擬研究 ［D］.天津：天津大學(xué)機(jī)械學(xué)院，2009.

［13］Peng Y，Shu C，Chew Y T.Lattice kinetic scheme for the incompressible viscous thermal flows on arbitrary meshes［J］.Physical Review，2004，E 69，016703.

［14］ 王秋香，戴傳山，李琪.流向?qū)苁饣旌蠈?duì)流影響的數(shù)值模擬［J］.工程熱物理學(xué)報(bào)，2011，32（7）：1221－1224.

［15］ 王秋香.微管換熱器傳熱與流動(dòng)性能的實(shí)驗(yàn)研究 ［D］.天津：天津大學(xué)機(jī)械學(xué)院，2010.