過(guò)冷沸騰多相流模擬及發(fā)動(dòng)機(jī)冷卻結(jié)構(gòu)研究＊

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

（天津大學(xué)，內(nèi)燃機(jī)燃燒學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072）

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

1 前言

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

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

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

2 多相流模型

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

2.1 Eulerian模型

質(zhì)量方程為：

動(dòng)量方程為：

能量方程為：

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

2.1.1 動(dòng)量交換模型

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

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

升力方程為：

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

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

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

壁面潤(rùn)滑力方程為：

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

流體對(duì)流體的曳力模型方程為：

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

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

式中，μtq為動(dòng)力粘度；δpq=0.75；αp為p相的體積分?jǐn)?shù)。

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

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

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

在流場(chǎng)內(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ù)來(lái)表征相間努塞爾數(shù)Nu：

式中，Re為相對(duì)雷諾數(shù)；Pr為普朗特?cái)?shù)。

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

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

2.1.3 壁面沸騰換熱模型

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

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

圖1 壁面換熱形式和氣泡運(yùn)動(dòng)行為

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

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

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

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

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

表2 部分壁面沸騰模型

2.2 Mixture模型

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

混合相動(dòng)量方程為：

混合相能量方程為：

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

Mixture模型是對(duì)混合相進(jìn)行求解，再通過(guò)滑移速度求解各相。

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

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

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

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

2.3 模型篩選流程

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

圖2 多相流模型篩選流程

3 發(fā)動(dòng)機(jī)過(guò)冷沸騰仿真

3.1 發(fā)動(dòng)機(jī)模型

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

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

3.2 邊界條件

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

表3 發(fā)動(dòng)機(jī)熱邊界條件

3.3 發(fā)動(dòng)機(jī)過(guò)冷沸騰仿真結(jié)果

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

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

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

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

圖5 發(fā)動(dòng)機(jī)溫度測(cè)點(diǎn)布置

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

3.3.2 換熱模型比較

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

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

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

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

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

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

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

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

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

4.1 交界面溫度

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

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

由圖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，部分位置達(dá)到460 K。上述區(qū)域溫度都達(dá)到了冷卻液的沸點(diǎn)，而冷卻通道外側(cè)表面溫度又低于冷卻液沸點(diǎn)，說(shuō)明冷卻液在發(fā)動(dòng)機(jī)冷卻通道內(nèi)發(fā)生過(guò)冷沸騰現(xiàn)象。方案2采用一側(cè)入口流動(dòng)的方式，冷卻液主要由一側(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)冷卻液流動(dòng)減弱區(qū)域（第5缸）與圖9中高溫交界面所在區(qū)域一致，均在同一缸；由圖10b可知，冷卻液流動(dòng)減弱區(qū)域與圖9中缸體交界面高溫區(qū)域相同，但在同一缸的缸蓋區(qū)域沒(méi)有受到明顯影響而產(chǎn)生高溫；由圖10c可知，缸蓋內(nèi)冷卻液流動(dòng)更加均勻，這與交界面溫度分布一致，但分離后的缸體中間位置的冷卻通道內(nèi)冷卻液流動(dòng)效果較差，其中高溫區(qū)域即由此引起。

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

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

4.3 冷卻液氣相分布

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

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

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

5 結(jié)束語(yǔ)

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