柴油機(jī)缸內(nèi)輻射換熱三維數(shù)值模擬

付麗榮，張文平，明平劍，羅躍生

(1.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，黑龍江哈爾濱150001;2.哈爾濱工程大學(xué)理學(xué)院，黑龍江哈爾濱150001)

高負(fù)荷柴油機(jī)的輻射傳熱占總傳熱量的10%～45%［1］，直接關(guān)系到發(fā)動(dòng)機(jī)熱效率以及因傳熱引起的各種熱負(fù)荷、熱強(qiáng)度問(wèn)題;同時(shí)，輻射熱流量對(duì)燃燒性能的影響很大，因此需對(duì)缸內(nèi)輻射傳熱進(jìn)行深入研究［2-3］。柴油機(jī)缸內(nèi)輻射傳熱的研究主要基于實(shí)驗(yàn)［4-6］和數(shù)值模擬［7-9］;實(shí)驗(yàn)研究法，即用熱輻射傳感器直接測(cè)量氣缸內(nèi)輻射傳熱量或者用二色法測(cè)量氣缸內(nèi)火焰溫度，再求得輻射傳熱量。數(shù)值模擬主要利用商用軟件，如 Fluent［8］、AVL-Fire［9］。目前主要的輻射傳熱模擬方法有:區(qū)域法、蒙特卡洛法、離散傳遞法、有限體積法、無(wú)網(wǎng)格法［10］等。而有限體積法，可確保輻射能量整體守恒，對(duì)不規(guī)則邊界適應(yīng)性強(qiáng)，易處理各向異性散射，可以采用與流場(chǎng)計(jì)算相同的網(wǎng)格，在迭代計(jì)算中可以不對(duì)溫度和熱流進(jìn)行插值計(jì)算。因此本文采用有限體積法進(jìn)行編程，從而可以與通用流體力學(xué)數(shù)值模擬軟件(general transportation equation analysis，GTEA)［11］的流場(chǎng)計(jì)算采用同一套網(wǎng)格，簡(jiǎn)化了計(jì)算程序，提高了計(jì)算的通用性。

本文在前期工作基礎(chǔ)上，針對(duì)三維輻射傳熱進(jìn)行數(shù)值模擬。本文借助有限體積法計(jì)算三維封閉空腔內(nèi)的輻射傳熱，并用3個(gè)算例——六面體、圓柱體和TBD620型柴油機(jī)的ω燃燒室來(lái)驗(yàn)證此方法和程序的可靠性與準(zhǔn)確性，同時(shí)考慮了網(wǎng)格離散個(gè)數(shù)和立體角離散個(gè)數(shù)對(duì)計(jì)算結(jié)果的影響。

1 控制方程

1.1 輻射傳熱方程

對(duì)于吸收、發(fā)射、散射性灰介質(zhì)的輻射傳遞方程表達(dá)式如下［12］

式中:κ、σs分別為吸收系數(shù)和散射系數(shù)，β0=κ+σs為衰減系數(shù);I為空間位置r、方向s處的輻射強(qiáng)度;Ib為黑體輻射強(qiáng)度，Ib=6T4/π;σ為斯蒂芬-玻爾茲曼常數(shù)，為 5.67×10-8W/(M2·K4)Φ(s，si)為輻射從入射方向s到散射方向si的散射相函數(shù);Ωi是中心方向?yàn)閟i的立體角。

本文所研究的介質(zhì)是非散射的，所以方程(1)可以簡(jiǎn)化為

1.2 輻射傳熱方程的邊界條件

RTE的邊界條件定義如下

式中:εw為壁面發(fā)射率，下標(biāo)w表示壁面位置，nw為單位法向量。式(3)表示當(dāng)輻射能離開(kāi)一個(gè)光學(xué)反射面時(shí)有2部分的貢獻(xiàn)，一是高溫表面放出的輻射量，二是投射到此表面的反射量。

2 數(shù)值離散

2.1 有限體積法

為了得到離散化方程，將方程(2)對(duì)控制體ΔV和立體角ΔΩmn積分，可得ΔΩmn內(nèi)輻射能量守恒方程的有限體積表達(dá)式:

應(yīng)用Gauss定理把體積分轉(zhuǎn)換成面積分，假設(shè)輻射強(qiáng)度I和Ib在控制體和立體角內(nèi)為常數(shù)，同時(shí)，源項(xiàng)在之內(nèi)為均勻的，則方程(4)的左端和右端分別可以表示為如下形式:

將式(5)、(6)代入到方程(4)，可得

式中:

將空間區(qū)域離散為互不重疊的四面體或三棱柱單元，所有變量信息存儲(chǔ)在單元中心，如圖1(a)的點(diǎn)1、2、3、4所示。將4π角域離散為互不重疊的控制立體角 ( Nθ×Nφ)，把方向矢量smn定義為立體角的中心方向，θ是0～π變化的極角，φ是0～2π變化的方位角，如圖1(b)所示?？刂平潜黄椒殖搔う萴=θm+-θm-= π/Nθ，Δφn= φn+-φn-=2π/Nφ。當(dāng)立體角是確定時(shí)，方向權(quán)的符號(hào)決定了控制體表面的輻射能是流進(jìn)還是流出，它由下式表示:

本文采用一階迎風(fēng)格式將單元界面上的輻射強(qiáng)度與單元中心的輻射強(qiáng)度關(guān)聯(lián)起來(lái)，則它們的關(guān)系表示如下

式中:

將式(13)代入式(7)，則式(7)可寫(xiě)成如下線(xiàn)性方程組:

式中:

圖1 混合網(wǎng)格和立體角示意圖Fig.1 Diagram of hybrid grids and the solid angle

2.2 邊界條件離散

邊界條件可以離散成如下形式

式中:

2.3 壁面輻射熱流密度

壁面輻射熱流密度表示為

3 數(shù)值計(jì)算結(jié)果及分析

本文的輻射換熱數(shù)值模擬是應(yīng)用GTEA軟件，在文獻(xiàn)［13］基礎(chǔ)上，為了驗(yàn)證本文的UFVM方法，并且在三維輻射問(wèn)題上檢驗(yàn)所開(kāi)發(fā)的計(jì)算程序的可靠性，本文選擇2個(gè)基本的幾何體，即六面體、圓柱體。

3.1 不規(guī)則六面體輻射模擬

驗(yàn)證算例1是對(duì)不規(guī)則六面體區(qū)域內(nèi)的輻射進(jìn)行仿真并與文獻(xiàn)中結(jié)果進(jìn)行對(duì)比。幾何體結(jié)構(gòu)的建立和參數(shù)的設(shè)定與文獻(xiàn)［12］中的數(shù)據(jù)一致。如圖2(a)所示，其中z0=1 m，所有邊界為黑體面，溫度為0，壁面發(fā)射率為εw=1，內(nèi)部介質(zhì)溫度為100 K。空間離散為六面體單元，見(jiàn)圖 2(b)，(Nx×Ny×Nz)= (27×27×27)，立體角劃分為( Nθ×Nφ)= ( 8×16)，吸收系數(shù)分別采用10、1和0.1 m-1進(jìn)行計(jì)算，不考慮散射，計(jì)算底面的無(wú)量綱熱流量分布。本文將計(jì)算結(jié)果(如圖3所示)與精確解和數(shù)值解［12］進(jìn)行對(duì)比，可以看到本文的計(jì)算結(jié)果與精確解吻合很好。

圖2 六面體模型與網(wǎng)格Fig.2 A hexahedral enclosure and mesh

由圖3可見(jiàn)，當(dāng)吸收系數(shù)是10 m-1時(shí)，無(wú)量綱熱流在底面中心處的值接近于1.0，主要是由于底面此時(shí)只受到周?chē)橘|(zhì)的影響;在兩邊處熱流迅速減小，因?yàn)榇藭r(shí)兩邊受到邊界冷壁面的影響。當(dāng)吸收系數(shù)為0.1 m-1時(shí)，熱流急劇下降，這是因?yàn)榻橘|(zhì)被冷邊界冷卻的結(jié)果。

圖3 底面上AB線(xiàn)段無(wú)量綱熱流比較Fig.3 Comparison of radiative wall heat flux along line AB on the bottom wall of a quadrilateral enclosure

圖4 網(wǎng)格的影響Fig.4 The effect of different discrete grids

圖5 網(wǎng)格和立體角劃分均不同時(shí)的誤差Fig.5 The error of different both grids and solid angles

圖4給出了吸收系數(shù)為1 m-1時(shí)，立體角離散個(gè)數(shù)均為 ( Nθ×Nφ)= ( 4×8)，空間區(qū)域分別離散為19 683個(gè)單元和12 167個(gè)單元時(shí)計(jì)算結(jié)果與精確解的相對(duì)誤差分布。由圖可知網(wǎng)格加密后的相對(duì)誤差不但沒(méi)有明顯的變小，反而在某些位置還有所變大。當(dāng)空間離散為12 167個(gè)單元，立體角離散個(gè)數(shù)為 ( Nθ×Nφ)= ( 4×8)時(shí)，與單元個(gè)數(shù)為 19 683，立體角離散個(gè)數(shù)為 ( Nθ×Nφ)= ( 8×16)時(shí)，計(jì)算結(jié)果與精確解的相對(duì)誤差分析見(jiàn)圖5。由圖可看出當(dāng)網(wǎng)格加密且立體角離散個(gè)數(shù)增多時(shí)，計(jì)算誤差明顯減小，計(jì)算結(jié)果更接近于精確解。因此，在數(shù)值計(jì)算中，空間離散的網(wǎng)格并不是越密越好，而是要和立體角離散的個(gè)數(shù)相對(duì)應(yīng)，即如果空間網(wǎng)格離散的個(gè)數(shù)增加，那么相應(yīng)的立體角離散的個(gè)數(shù)也要增加，才可以得到更高精度的數(shù)值解。

3.2 圓柱體輻射模擬

驗(yàn)證算例2是半徑為1m，z軸方向長(zhǎng)度為2m的圓柱體。邊界為黑體，溫度為0，壁面發(fā)射率為εw=1，內(nèi)部氣體溫度為100 K，圓柱體和空間網(wǎng)格見(jiàn)圖6。在吸收系數(shù)分別為0.1，1.0和5.0m-1情況下，計(jì)算圓柱體右側(cè)邊界的無(wú)量綱輻射熱流量，空間域離散為74 203個(gè)四面體單元，立體角劃分為 ( Nθ×Nφ)= ( 8×16)。本文的計(jì)算結(jié)果與文獻(xiàn)［14］提供的精確解和文獻(xiàn)［15-16］的數(shù)值解進(jìn)行對(duì)比，結(jié)果如圖7所示，可以看出它與精確解和數(shù)值解吻合很好。

圖6 圓柱體和網(wǎng)格Fig.6 The cylinder enclosure and mesh

圖7 無(wú)量綱輻射熱流量的對(duì)比Fig.7 Comparison of dimensionless radiative heat flux along the lateral side of the cylinder

圖7顯示了在3種吸收系數(shù)0.1、1.0和5.0 m-1下右側(cè)邊界的無(wú)量綱輻射熱流量，當(dāng)吸收系數(shù)是5.0 m-1時(shí)，到達(dá)右側(cè)壁面中間處的輻射熱流密度接近于介質(zhì)的黑體輻射強(qiáng)度，這是由于壁面的輻射強(qiáng)度只受封閉壁面周?chē)臒釟怏w的發(fā)射本領(lǐng)影響。然而，在兩邊角處，輻射熱流迅速減小，因?yàn)槭苤車(chē)敳亢偷撞坷浔诿嬗绊?。但?dāng)吸收系數(shù)是0.1 m-1時(shí)，介質(zhì)的輻射本領(lǐng)是弱的，并且無(wú)量綱輻射熱流量大幅度下降，這是因?yàn)槭芷渌浔诿婧凸鈱W(xué)薄氣體可忽略的自身衰減的深遠(yuǎn)影響。

4 ω型燃燒室輻射模擬

驗(yàn)證算例3是TBD620型柴油機(jī)ω型燃燒室。發(fā)動(dòng)機(jī)的基本參數(shù):氣缸數(shù)為16，沖程為195 mm，缸徑為170 mm，連桿長(zhǎng)度為350 mm，壓縮比為13.5。在燃燒結(jié)束至排氣開(kāi)始前這一階段，此時(shí)活塞在下止點(diǎn)位置，缸蓋、缸壁和活塞溫度均為350 K，假設(shè)氣體平均溫度為800 K;利用GAMBIT軟件劃分網(wǎng)格，通過(guò)軟件將網(wǎng)格文件轉(zhuǎn)換為CGNS格式，并讀入本文開(kāi)發(fā)的求解器;空間域內(nèi)采用混合網(wǎng)格，氣缸內(nèi)為三棱柱非結(jié)構(gòu)化網(wǎng)格，ω型燃燒室內(nèi)采用四面體網(wǎng)格，網(wǎng)格總數(shù)為63 227，圖8(a)、(b)分別顯示了計(jì)算網(wǎng)格的主視和俯視圖;立體角劃分為( Nθ×Nφ)= ( 8×16)。計(jì)算燃燒室縱向中心截面邊界線(xiàn)左側(cè)、下側(cè)和右側(cè)的輻射熱流密度，圖9中s為0～200、200～432.7、432.7～632.7mm 區(qū)間分別代表邊界線(xiàn)左側(cè)、下側(cè)和右側(cè)。

圖8 計(jì)算網(wǎng)格Fig.8 Computational grid

圖9 s的定義Fig.9 Definition of s

圖10顯示的是工質(zhì)吸收系數(shù)對(duì)缸內(nèi)縱向中心截面的左側(cè)、下側(cè)和右側(cè)邊界無(wú)量綱輻射熱流密度的影響。從432.7～632.7 mm區(qū)間可以看出3種吸收系數(shù)的熱流量圖都是對(duì)稱(chēng)的，與算例2的圓柱體模型結(jié)果圖7對(duì)比可知，算例3的底面的ω型燃燒室與算例2底面是圓形的效果一樣，所得的側(cè)面熱流量圖都是對(duì)稱(chēng)的，即底面的形狀對(duì)側(cè)邊的輻射熱流量沒(méi)有影響。在0和200 mm周?chē)?，由于兩邊角處受頂部和底部冷壁面影響，所以熱流量值變小。?dāng)吸收系數(shù)為0.1 m-1時(shí)，工質(zhì)相當(dāng)于是透明的，所以各位置的熱流量幾乎一樣大。

圖10 不同吸收系數(shù)下的無(wú)量綱輻射熱流量的對(duì)比Fig.10 Comparison of dimension less radiative heat flux of different absorption coefficients

圖11顯示的是吸收系數(shù)為10.0 m-1時(shí)，工質(zhì)溫度分別為600、800、1 000 K時(shí)缸壁無(wú)量綱輻射熱流量。介質(zhì)溫度對(duì)中心部位的熱流量影響較大，對(duì)邊角處的影響較小;邊角處由于受冷壁面的影響，熱流量值變小。

圖11 介質(zhì)不同溫度下的無(wú)量綱輻射熱流量的對(duì)比Fig.11 Comparison of dimensionless radiative heat flux of medium of different temperatures

5 結(jié)論

本文通過(guò)六面體、圓柱體及ω型燃燒室的數(shù)值模擬，得出如下結(jié)論:

1)本文方法的數(shù)值解與精確解和文獻(xiàn)中的數(shù)值解吻合較好，可見(jiàn)，本文的非結(jié)構(gòu)有限體積法和求解器GTEA適用于三維輻射換熱的數(shù)值模擬。

2)不同的網(wǎng)格劃分和立體角劃分對(duì)計(jì)算結(jié)果有相應(yīng)的影響，空間離散的網(wǎng)格并不是越密越好，而是要和立體角離散的個(gè)數(shù)相對(duì)應(yīng)，即空間網(wǎng)格離散的個(gè)數(shù)要與立體角離散的個(gè)數(shù)同時(shí)增加，空間離散產(chǎn)生的假散射與角度離散產(chǎn)生的射線(xiàn)效應(yīng)在一定程度上才會(huì)互相抵消，所以可以得到更高精度的數(shù)值解。

3)可將程序應(yīng)用于柴油機(jī)燃燒室的輻射換熱模擬，燃燒室缸壁的輻射熱流量會(huì)受到缸內(nèi)工質(zhì)的吸收系數(shù)和溫度的影響;無(wú)量綱輻射熱流密度隨著工質(zhì)吸收系數(shù)的增大而增大，工質(zhì)溫度對(duì)中心部位的熱流量影響較大，對(duì)邊角處的影響較小。為后續(xù)考慮碳黑顆粒的輻射傳熱對(duì)柴油機(jī)工作性能的影響打下了基礎(chǔ)。

［1］沈季勝，沈瑜銘，陳紅巖，等.柴油機(jī)氣缸內(nèi)傳熱模型研究歷程與展望［J］.內(nèi)燃機(jī)工程，1999，20(1):72-77.SHEN Jisheng，SHEN Yuming，CHEN Hongyan，et al.Prospect and research process on radiated heat transfer models in cylinder of diesel engine［J］.Chinese Internal Combustion Engine Engineering，1999，20(1):72-77.

［2］嚴(yán)兆大，沈季勝，劉震濤，等.用蒙特卡洛法計(jì)算柴油機(jī)缸內(nèi)多元熱輻射［J］.內(nèi)燃機(jī)學(xué)報(bào)，2000，18(2):133-136.YAN Zhaoda，SHEN Jisheng，LIU Zhentao，et al.Calculation of the radiation heat transfer in cylinder of diesel engine by the Monte-Carlo method［J］.Transactions of CSICE，2000，18(2):133-136.

［3］熊仕濤，陳國(guó)華.柴油機(jī)缸內(nèi)傳熱計(jì)算［J］.內(nèi)燃機(jī)學(xué)報(bào)，2001，19(3):215-218.XIONG Shitao，CHEN Guohua.Calculation of in-cylinder heat transfer of diesel engine［J］.Transactions of CSICE，2001，19(3):215-218.

［4］潘克煜，周龍保，楊中樂(lè).直噴式柴油機(jī)氣缸內(nèi)輻射傳熱的試驗(yàn)研究［J］.內(nèi)燃機(jī)工程，1996，17(4):59-67.PAN Keyu，ZHOU Longbao，YANG Zhonglu.Experimental investigation of in-cylinder heat radiation in diesel engine［J］.Chinese Internal Combustion Engine Engineering，1996，17(4):59-67.

［5］SAID A M，BUTTSWORTH D R，YUSAF T F.A review of radiation heat transfer measurement for diesel engines using the two-colour method［C］//Proceeding of 3rd International Conference on Energy and Environment.Malaysia，2009:202-207.

［6］高天，潘克煜，陳飛.直噴式柴油機(jī)缸內(nèi)氣體輻射傳熱的研究［J］.西安交通大學(xué)學(xué)報(bào)，2006，40(3):302-306.GAO Tian，PAN Keyu，Chen Fei.Investigation on in-cylinder radiation from gases in DI diesel engine［J］.Journal of Xi＇an Jiaotong University，2006，40(3):302-306.

［7］宋有生，聶宇宏.燃燒室型線(xiàn)對(duì)缸內(nèi)燃燒和輻射傳熱影響的計(jì)算分析［J］.科學(xué)技術(shù)與工程，2008，8(3):760-763.SONG Yousheng，NIE Yuhong.Calculation investigation on influence of combustion chamber pr of ile to combustion process and radiative heat transfer［J］.Science Technology and Engineering，2008，8(3):760-763.

［8］姜曉光，聶宇宏，謝凱弘，等.基于Fluent的柴油機(jī)缸內(nèi)輻射換熱分析［J］.科學(xué)技術(shù)與工程，2008，8(15):4099-4103.JIANG Xiaoguang，NIE Yuhong，XIE Kaihong，et al.Diesel engine cylinder internal radiation heat transfer analysis based on Fluent［J］.Science Technology and Engineering，2008，8(15):4099-4103.

［9］劉健，呂繼組，計(jì)時(shí)鳴.柴油機(jī)缸內(nèi)輻射換熱的多維數(shù)值模擬研究［J］.機(jī)械工程學(xué)報(bào)，2009，45(12):311-317.LIU Jian，LYU Jizu，JIShiming.Study onmulti-dimensional numerical simulation of in-cylinder radiation heat transfer of diesel engine［J］.Journal of Mechanical Engineering，2009，45(12):311-317.

［10］SADAT H，WANG C A，DEZ V L.Meshless method for solving coupled radiative and conductive heat transfer in complexmulti-dimensional geometries［J］.Applied Mathematics and Computation，2012，218:10211-10225.

［11］明平劍.基于非結(jié)構(gòu)化網(wǎng)格氣液兩相流數(shù)值方法及并行計(jì)算研究與軟件開(kāi)發(fā)［D］.哈爾濱:哈爾濱工程大學(xué)，2008:16-72.MING Pingjian.Development of numerical modeling for gas-liquid two-phase flows based on unstructured grids and parallell computing［D］.Harbin:Harbin Engineering University，2008:16-72.

［12］BEAk SW，KIM M Y，KIM JS.Nonorthogonal finite-volume solutions of radiative heat transfer in a three-dimensional enclosure［J］.Numer Heat Transfer B，1998，34:419-437.

［13］MING Pingjian，ZHANG Wenping.Numerical simulation of natural convection and radiation heat transfer in 2D enclosure on hybrid grids［J］.Numerical Heat Transfer:Part B，2012，59:116-137.

［14］KIM C，KIM M Y，YU M J，et al.Unstructured polygonal finite volume solution of radiative heat transfer in a complex axisymmetric enclosure［J］.Numerical Heat Transfer:Part B.2010，57:227-239.

［15］MURTHY JY ，MATHUR SR.Radiative heat transfer in axisymmetric geometries using an unstructured finite-volume method［J］.Numer Heat Transfer B，1998，33:397-416.

［16］KIM M Y.Assessment of the axisymmetric radiative heat transfer in a cylindrical enclosure with the finite volume method［J］.Int JHeat Mass Transfer，2008，51:5144-5153.