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

（天津大學，內(nèi)燃機燃燒學國家重點實驗室，天津 300072）

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

1 前言

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

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

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

2 多相流模型

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

2.1 Eulerian模型

質(zhì)量方程為：

動量方程為：

能量方程為：

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

2.1.1 動量交換模型

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

式中，F(xiàn)eb為外部體積力；Flift為升力；Fvm為虛擬質(zhì)量力；FTD為湍流離散力；FWL為壁面潤滑力；Fdrag為曳力，是兩相間最主要的力。

升力方程為：

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

虛擬質(zhì)量力方程為：

式中，ρp為p相的密度；Cvm=0.5為虛擬質(zhì)量力系數(shù)。

壁面潤滑力方程為：

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

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

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

湍流離散力與曳力關(guān)系密切，其方程為：

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

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

表1 部分相間作用力系數(shù)模型

2.1.2 內(nèi)部換熱模型

在流場內(nèi)部，質(zhì)量交換主要取決于兩相的溫度、換熱系數(shù)和潛熱，溫度決定熱量傳遞的方向，質(zhì)量交換率為：

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

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

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

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

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

2.1.3 壁面沸騰換熱模型

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

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

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

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

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

氣泡覆蓋面積與氣泡脫離直徑、核點密度有關(guān)：

式中，dd為氣泡脫離直徑；K為經(jīng)驗常數(shù)，可由Del Valle and Kenning[14]關(guān)聯(lián)式獲得；Ja為雅克比數(shù)。

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

表2 部分壁面沸騰模型

2.2 Mixture模型

混合相質(zhì)量方程為：

混合相動量方程為：

混合相能量方程為：

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

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

滑移速度與相對速度的關(guān)系為：

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

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

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

2.3 模型篩選流程

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

圖2 多相流模型篩選流程

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

3.1 發(fā)動機模型

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

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

3.2 邊界條件

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

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

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

3.3.1 曳力系數(shù)模型比較

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

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

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

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

表4 火力面溫度測量結(jié)果 K

3.3.2 換熱模型比較

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

3.3.3 非曳力系數(shù)模型比較

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

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

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

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

圖8 Mixture模型仿真結(jié)果

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

4 冷卻通道結(jié)構(gòu)研究

4.1 交界面溫度

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

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

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

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

4.2 冷卻液流速

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

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

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

4.3 冷卻液氣相分布

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

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

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

5 結(jié)束語

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