方腔內冷熱微管自然對流換熱的多格子LBM研究

焦文靜，戴傳山，王秋香

（天津大學機械學院，中低溫熱能高效利用教育部重點實驗室，天津300072）

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

1 多格子LBM

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

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

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

為了滿足黏度為常數，粗網格的松弛時間τc和細網格的松弛時間τf必須滿足式 （3）。

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

圖1 粗、細網格邊界結構圖

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

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

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

2 驗證程序正確性

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

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

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

瑞利數Ra

局部Nu數

表面平均Nu數

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

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

圖3 溫度場和流場的模擬結果對比

表1 平均努塞爾數Nu模擬結果的對比

3 冷熱微圓管自然對流理論模型

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

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

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

4 模擬結果與分析

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

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

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

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

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

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

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

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

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

表2 不同布管方式下平均Nu對比

5 結 論

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

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

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

符 號 說 明

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

Δx——網格空間步長，量綱為1

τ——網格的松弛時間，量綱為1

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

Ra——瑞利數，量綱為1

f——波爾茲曼密度分布函數

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

β——體積膨脹系數，K－1

α——熱擴散系數，m2／s

Nu——努塞爾數，量綱為1下角標

f——細網格

c——粗網格

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

［2］ 華賁，仵浩，劉二恒.基于火用經濟評價的換熱器最優傳熱溫差化［J］.化工進展，2009，28（7）：1142－1146.

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

［4］ 陸應生，陳墓玲，潘寧忠，等.強化傳熱元件與高效換熱器研究進展［J］.化工進展，1998，1：46－48.

［5］ 戴傳山，王秋香，孫平樂.污垢對微管換熱器用于發電系統的影響［J］.化工進展，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］.天津：天津大學機械學院，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］ 王秋香，戴傳山，李琪.流向對管束外混合對流影響的數值模擬［J］.工程熱物理學報，2011，32（7）：1221－1224.

［15］ 王秋香.微管換熱器傳熱與流動性能的實驗研究 ［D］.天津：天津大學機械學院，2010.