張俊紅 徐喆軒 王靜超 徐天舒 林澤峰 馬梁

（天津大學，內燃機燃燒學國家重點實驗室，天津 300072）

主題詞：發(fā)動機 過冷沸騰 多相流 相間模型 冷卻結構

1 前言

發(fā)動機冷卻系統(tǒng)直接影響發(fā)動機的熱負荷水平，其性能的優(yōu)劣對整機的動力性、經濟性和排放性能具有重要影響。發(fā)動機冷卻水腔結構復雜，內部結構不可見，導致設計難度較大，常用的試驗方法難以應用，因此采用計算流體力學（Computational Fluid Dynamics，CFD）對其進行優(yōu)化設計受到研究人員的廣泛關注。傳統(tǒng)無相變CFD的方法由于采用顯熱代替潛熱的方式，準確性不高，而考慮相變的多相流CFD模擬可以表征潛熱吸熱的過程，準確地模擬流體溫度場。

傅松[1]分別采用對流傳熱修正傳熱系數的方法和均勻多相流沸騰傳熱模型建立了用于缸蓋冷卻通道內沸騰傳熱計算的修正算法；Arash[2]采用兩相的Mixture多相流模型對柴油機冷卻通道進行模擬，得到冷卻通道內壓力、溫度、體積分數的分布；谷芳[3]利用流體體積（Volume of Fluid，VOF）兩相流方法建立過冷沸騰模型，采用T型管和鑄鋁缸蓋進行試驗驗證，模型適應性很高；花仕洋[4]采用歐拉（Eulerian）多相流模型對單缸發(fā)動機缸蓋進行模擬，與無相變對流傳熱結果對比，Eulerian模型計算結果更準確；張俊紅[5]采用混合（Mixture）多相流模型簡化復雜的過冷沸騰問題，提高了計算準確性。上述發(fā)動機冷卻通道的傳熱仿真研究都將多相流與對流傳熱進行比較，以證明仿真的準確性，但由于研究中采用的多相流模型多種多樣，因而無法說明哪一種模型更適合發(fā)動機冷卻通道的過冷沸騰模擬。

為準確選擇多相流模型，本文針對某發(fā)動機冷卻通道中的過冷沸騰現象，建立多種過冷沸騰模型，通過模擬驗證和比較建立準確高效的多相流模型，并分析不同結構發(fā)動機冷卻通道的流動換熱現象。

2 多相流模型

換熱設備中宏觀多相流模擬常用Eulerian模型和Mixture模型，這兩個模型同樣應用于發(fā)動機冷卻通道的多相流模擬。

2.1 Eulerian模型

質量方程為：

動量方程為：

能量方程為：

式中，αq、ρq分別為q相的體積分數和密度；m pq、mqp分別為p相到q相和q相到p相的質量交換；up、uq分別為p相和q相的速度；uq、uqp分別為p相相對于q相、q相相對于p相的速度；τq為剪切張力；g為重力；∑F為相間作用力之和；hq為q相的焓；?qq為熱通量；Sq為能量的源項；Qpq為換熱量；hpq、hqp為潛熱；n為相的總數目。

2.1.1 動量交換模型

兩相間的作用力十分復雜，一般將其分解：

式中，Feb為外部體積力；Flift為升力；Fvm為虛擬質量力；FTD為湍流離散力；FWL為壁面潤滑力；Fdrag為曳力，是兩相間最主要的力。

升力方程為：

式中，Cl為升力系數；αg為氣相體積分數。

虛擬質量力方程為：

式中，ρp為p相的密度；Cvm=0.5為虛擬質量力系數。

壁面潤滑力方程為：

式中，CWL為壁面潤滑力系數；nw為指向遠離墻壁方向的法向單位距離。

流體對流體的曳力模型方程為：

式中，Δu為兩相間的速度差；CD為曳力系數；db為氣泡直徑。

湍流離散力與曳力關系密切，其方程為：

式中，μtq為動力粘度；δpq=0.75；αp為p相的體積分數。

相間作用力除本身參數外，主要由各力系數決定，部分相間作用力系數模型如表1所示。

表1 部分相間作用力系數模型

2.1.2 內部換熱模型

在流場內部，質量交換主要取決于兩相的溫度、換熱系數和潛熱，溫度決定熱量傳遞的方向，質量交換率為：

式中，αi為界面濃度；hinter為界面處的換熱系數；Q為總熱量；QT為表熱；Tsat為飽和溫度；Tl為液相溫度；hlg為潛熱。

常用的界面換熱系數模型有Ranz-Marshall[12]、Tomiyama[10]、Hughmark[13]模型等。Ranz-Marshall引用雷諾數和普朗特（Prandtl）數來表征相間努塞爾數Nu：

式中，Re為相對雷諾數；Pr為普朗特數。

Tomiyama提出了適用于較低雷諾數流動情況下的換熱模型：

Hughmark拓展了Ranz-Marshall變化范圍：

2.1.3 壁面沸騰換熱模型

由于兩相的運動，流體與壁面形成3種換熱形式：氣相或液相對壁面的對流換熱，其中氣相對流換熱程度根據氣相份額決定，在核態(tài)沸騰階段基本可以忽略；由液相轉化成氣相的潛熱；由于氣泡離開壁面形成的氣液交替與壁面換熱的現象，簡稱“淬火”換熱。3種換熱形式如圖1所示，根據這3種換熱形式得到的壁面換熱量計算式為：

式中，Qw為壁面處總熱量；Qcon為對流換熱熱量；Qe為氣相潛熱熱量；Qq為“淬火”換熱熱量。

圖1 壁面換熱形式和氣泡運動行為

大量研究以壁面處的氣泡脫離直徑、氣泡脫離頻率和氣泡核點密度來描述壁面處的沸騰換熱程度：

式中，hcon為對流換熱系數；Tw為壁面溫度；Ab、(1-Ab)分別為氣相和液相覆蓋的面積；mb為氣泡在脫離時刻的質量；Nw為核點密度；kl為液體的導熱；f為氣泡脫離頻率；λl為擴散率。

氣泡覆蓋面積與氣泡脫離直徑、核點密度有關：

式中，dd為氣泡脫離直徑；K為經驗常數，可由Del Valle and Kenning[14]關聯式獲得；Ja為雅克比數。

氣泡脫離直徑、氣泡脫離頻率和氣泡核點密度的部分模型如表2所示。

表2 部分壁面沸騰模型

2.2 Mixture模型

混合相質量方程為：

混合相動量方程為：

混合相能量方程為：

式中，vm為質量加權平均速度；ρm為混合密度；μm為混合相的動力粘度；F為相間力；αk為k相體積分數；ρk為k相密度；vk為k相速度；hk為k相焓；vdr,k為二次相的漂移速度；keff為有效導熱系數；T為溫度；SE為其他的體積熱源。

Mixture模型是對混合相進行求解，再通過滑移速度求解各相。

滑移速度與相對速度的關系為：

式中，vpq、vqk分別為p相和k相對于主相q的相對速度；ck為k相的質量分數；a為加速度；η為湍流擴散率。

Mixture模型采用平均界面換熱系數簡化描述多相流的換熱過程（包括過冷沸騰），其中Lee[20]模型以熱為機理描述蒸發(fā)和冷凝過程：

式中，mlv、mvl分別為液相轉為氣相、氣相轉為液相的質量；Ccoe為蒸發(fā)或冷凝頻率；αl、αv分別為液相、氣相的體積分數；ρl、ρv分別為液相、氣相的密度；Tl為液相溫度；M為流體摩爾質量；R為氣體常數；hfg為潛熱。

2.3 模型篩選流程

對于多相流模型中的子模型，忽略適用范圍，約有幾百種配合方式。多相流模型在流場模擬應用上的首要問題是如何準確地選取模型。在本文涉及的發(fā)動機冷卻通道內流場數值模擬中，按模型的重要程度逐步進行篩選，得到最終的多相流模型，多相流模型篩選流程如圖2所示。

圖2 多相流模型篩選流程

3 發(fā)動機過冷沸騰仿真

3.1 發(fā)動機模型

針對某直列6缸4沖程柴油機，建立其機體和冷卻通道的三維幾何模型，并簡化部分結構生成網格，如圖3所示。方案1為原始冷卻通道網格模型，用來比較多相流的模擬結果；方案2為在原始冷卻通道的基礎上，將位于發(fā)動機中間的冷卻通道入口改到單獨一側，保持出口位置不變，考察流動方向引起的冷卻變化；方案3為在原始冷卻通道的基礎上，切斷缸體冷卻通道與缸蓋冷卻通道的聯系，將冷卻通道改為分離式，缸蓋冷卻通道入口設在兩缸冷卻通道之間，在原出口基礎上兩側各增加1個小出口，缸體冷卻通道入口位置不變，出口設置在缸體兩側，考察缸蓋缸體分離冷卻的效果。方案1～3作為對比進行發(fā)動機冷卻通道結構的研究。

圖3 發(fā)動機網格模型

3.2 邊界條件

根據柴油機廠商提供的熱邊界條件，考慮發(fā)動機缸套換熱的不均勻性，將缸套劃分為上、中、下3個面，并將熱邊界條件按面積權重進行修正，計算得到熱邊界條件如表3所示。冷卻液入口溫度為363 K，入口流量為15.1 m3/h。模擬中將發(fā)動機缸體缸蓋和冷卻液作為整體，采用流固耦合方法計算固液能量交換，采用k-ε湍流模型進行各相內的湍流計算。

表3 發(fā)動機熱邊界條件

3.3 發(fā)動機過冷沸騰仿真結果

3.3.1 曳力系數模型比較

Eulerian多相流模型中選擇Tolubinsky的氣泡脫離直徑模型、Cole的氣泡脫離頻率模型、Lemmer的氣泡核點密度模型組成壁面沸騰換熱模型，再選擇Ranz-Mar?shall的兩相間換熱系數模型，并將這些模型作為基準模型，以考察不同曳力系數模型在發(fā)動機冷卻通道中模擬的準確性。不同曳力系數模型的火力面溫度仿真結果如圖4所示，從溫度分布可看出，除Grace模型外，其余模型冷卻液入口附近的火力面溫度最高，冷卻液出口附近的火力面溫度最低，這與冷卻通道內流量分布一致。

圖4 采用不同曳力系數模型的發(fā)動機火力面溫度仿真結果

為選取最合適的曳力系數模型，利用熱電偶對該柴油發(fā)動機火力面溫度進行測量，測點布置在柴油機火力面和冷卻通道出、入口處，如圖5所示，火力面溫度測量結果如表4所示。將火力面溫度分布（圖4）與試驗測量值（表4）進行對比可發(fā)現，采用Schiller and Naumann曳力系數模型的仿真結果與試驗測量值最吻合，采用Grace曳力模型計算的溫度值過低，采用Tomiyama曳力模型和Ishii曳力模型的計算溫度值均過高。因此，本文選取Schiller and Naumann模型作為發(fā)動機冷卻通道過冷沸騰模擬的曳力系數模型。

圖5 發(fā)動機溫度測點布置

表4 火力面溫度測量結果 K

3.3.2 換熱模型比較

在采用Schiller and Naumann曳力系數模型的基礎上，選取不同界面換熱系數模型、氣泡脫離直徑和氣泡核點密度計算發(fā)動機火力面溫度，部分結果如圖6所示。整體而言，界面換熱系數模型對發(fā)動機和冷卻液耦合溫度場影響更大，表現為溫度結果浮動更明顯；采用Kocamustafaogullari氣泡核點密度模型時3個火力面溫度值相差不大，第2缸的預測結果十分接近試驗值，但整體誤差比基準模型大；采用Tomiyama的界面換熱系數模型的仿真結果大于試驗值。

3.3.3 非曳力系數模型比較

考慮Moraga的升力系數模型、Antal的壁面潤滑力模型等非曳力的相間作用力，在基準模型基礎上進行發(fā)動機過冷沸騰模擬，結果如圖7所示。對比圖7a與圖4a、圖7b與圖6c可發(fā)現，模擬結果的溫度差變化很小，這說明是否考慮除曳力之外的相間作用力，對高熱流的發(fā)動機及冷卻通道的模擬準確性影響不明顯。

圖6 采用不同換熱模型的發(fā)動機火力面溫度仿真結果

圖7 考慮非曳力的相間作用力的仿真結果

Mixture模型的計算結果如圖8所示，使用該模型的仿真結果與試驗值相比偏高，相對誤差明顯大于使用Eulerian模型得到的仿真結果，說明該模型的仿真精確度較低。

圖8 Mixture模型仿真結果

綜上所述，在進行發(fā)動機冷卻液過冷沸騰仿真時，采用Tolubinsky的氣泡脫離直徑模型、Cole的氣泡脫離頻率模型、Lemmer的氣泡核點密度模型、Ranz-Marshall的兩相間換熱系數模型以及Schiller and Naumann曳力系數模型得到的仿真結果與試驗測試結果最接近，使用上述子模型對發(fā)動機冷卻液過冷沸騰進行仿真既可以節(jié)約計算成本，又能夠保證仿真的準確性。

4 冷卻通道結構研究

4.1 交界面溫度

方案2中冷卻液入口條件與原始條件相同，即入口溫度為363 K，入口流量為15.1 m3/h。方案3中缸體的入口條件與原始條件相同，但缸蓋冷卻液入口條件有所變化，每個入口的冷卻液流量為原流量的1/5，保持總流量不變。3種結構冷卻通道的固液交界面溫度如圖9所示。

圖9 不同冷卻通道結構的發(fā)動機固液交界面溫度

由圖9可以看出，在方案1和方案3中，固液交界面或冷卻液表面溫度主要受冷卻液溫度影響，最小值為363 K，極少部分區(qū)域溫度較低。方案2中出現大面積低溫區(qū)域，該區(qū)域交界面溫度與缸體外表面溫度接近，為300 K左右；其余部分溫度與冷卻液溫度（363 K）相同。

3種方案溫度較高區(qū)域均出現在缸蓋冷卻通道近火力面處，以及缸體冷卻通道近燃燒室區(qū)域上部邊緣。方案1的最高溫度出現在第5缸的缸體冷卻通道近燃燒室區(qū)域上部邊緣，為525 K，缸體冷卻通道內側表面平均溫度約為430 K。缸蓋冷卻通道的最高溫度接近500 K，位于第4缸、第5缸缸蓋，內側平均溫度約為420 K，部分位置達到460 K。上述區(qū)域溫度都達到了冷卻液的沸點，而冷卻通道外側表面溫度又低于冷卻液沸點，說明冷卻液在發(fā)動機冷卻通道內發(fā)生過冷沸騰現象。方案2采用一側入口流動的方式，冷卻液主要由一側流向另一側，在流量相同的情況下，方案2大部分區(qū)域流速大于方案1，因此冷卻液與缸體缸蓋的交界面溫度低于方案1。方案3與方案1缸蓋冷卻通道流量相同，但溫度分布更均勻，在第3和第4缸體冷卻通道內出現高溫區(qū)，溫度向兩側逐缸遞減。

4.2 冷卻液流速

3種方案的流速分布如圖10所示，方案1和方案2呈現與之前分析一致的冷卻液流速差別。由圖10a可知，在缸蓋冷卻通道內冷卻液流動減弱區(qū)域（第5缸）與圖9中高溫交界面所在區(qū)域一致，均在同一缸；由圖10b可知，冷卻液流動減弱區(qū)域與圖9中缸體交界面高溫區(qū)域相同，但在同一缸的缸蓋區(qū)域沒有受到明顯影響而產生高溫；由圖10c可知，缸蓋內冷卻液流動更加均勻，這與交界面溫度分布一致，但分離后的缸體中間位置的冷卻通道內冷卻液流動效果較差，其中高溫區(qū)域即由此引起。

圖10 不同冷卻通道結構的發(fā)動機冷卻液流速

3種方案中，冷卻液最大流速均在出口位置，其中方案3由于增加了2個小出口，冷卻液的流速減小。分離式的冷卻通道減少了上、下冷卻通道的連通，需要通過增加冷卻液出、入口使冷卻液均勻分布，如方案3缸蓋冷卻通道，否則會使流動出現“死區(qū)”，如方案3缸體冷卻通道。串聯模式的發(fā)動機冷卻通道雖然有較多的連接口，但也產生冷卻液流動減弱區(qū)，如圖10中方案1的缸蓋冷卻通道和方案2的缸體冷卻通道。方案1采用對稱的3缸串聯模式，方案2采用6缸串聯模式。

4.3 冷卻液氣相分布

過冷沸騰產生的氣體分布如圖11所示，方案1和方案3中都有較大面積的氣相覆蓋區(qū)產生，而方案2氣相覆蓋區(qū)面積很小。氣相覆蓋區(qū)均產生在高熱流量（燃燒室附近）區(qū)域，外側氣相分布為0，這與過冷沸騰現象吻合。方案1和方案3分別在第5缸缸蓋和第3缸、第4缸缸體產生比其他缸更多氣相，與流速減弱區(qū)域出現的位置相同。方案2中由于冷卻液流速較高，氣相分布區(qū)產生較少。

圖11 不同冷卻通道結構的發(fā)動機冷卻液氣相分布

相對于方案3，方案1的缸蓋冷卻通道氣相分布更廣，這也是由冷卻液流動不均勻引起的。方案1的缸體冷卻通道由中間向兩側氣相分布逐漸增加，而方案3則是在中間2缸集中更多的氣相。方案1由于冷卻液溫度逐漸上升，使冷卻液到達飽和溫度所需要的熱量降低，更多富余的熱量轉化為液體汽化所需的潛熱；方案3由于流動緩慢進而導致冷卻液加熱時間過長。

5 結束語

本文針對發(fā)動機冷卻通道過冷沸騰多相流模擬中采用模型的復雜性和多樣性問題，通過一系列仿真模擬和試驗驗證表明：發(fā)動機冷卻液過冷沸騰多相流模擬中，采用Eulerian模型比Mixture模型模擬結果更加準確；在高入口溫度和高熱流量的條件下，在Eulerian模型中采用Tolubinsky的氣泡脫離直徑模型、Cole的氣泡脫離頻率模型、Lemmer的氣泡核點密度模型、Ranz-Marshall的兩相間換熱系數模型、Schiller and Naumann曳力系數模型，可以得到準確的模擬結果。通過模擬分析可知，分離式冷卻通道結構設計可以使溫度分布更均勻，但需要更多的出入口設計，連體式的冷卻通道結構需要注重上、下冷卻通道的串聯和出、入口位置的選擇。