斯特林循環(huán)中射流現(xiàn)象及其對(duì)循環(huán)特性影響

劉柏岙，王樹(shù)林，肖 剛

（能源清潔利用國(guó)家重點(diǎn)實(shí)驗(yàn)室浙江大學(xué)能源工程學(xué)院，浙江 杭州 310027）

斯特林發(fā)動(dòng)機(jī)應(yīng)用前景較為廣泛，其在余熱回收[1-2]，由斯特林發(fā)動(dòng)機(jī)構(gòu)成的微型熱電聯(lián)產(chǎn)系統(tǒng)[3-5]和碟式斯特林發(fā)電系統(tǒng)[6-8]等方面受到關(guān)注。

段晨[9]建立了考慮“不可逆因子”的斯特林循環(huán)熱力學(xué)模型，采用相似設(shè)計(jì)和多目標(biāo)優(yōu)化方法設(shè)計(jì)制造了1 臺(tái)β型斯特林發(fā)動(dòng)機(jī)。瞿凡等[10]基于1 臺(tái)β型斯特林發(fā)動(dòng)機(jī)，對(duì)加熱器進(jìn)行了輻射與燃燒2 種加熱方式的數(shù)值模擬與優(yōu)化設(shè)計(jì)。黃護(hù)林等[11]對(duì)斯特林發(fā)動(dòng)機(jī)回?zé)崞鬟M(jìn)行數(shù)值模擬，研究了長(zhǎng)徑比和孔隙率對(duì)回?zé)崞饔行约傲髯钃p失的影響。王思[12]考慮回?zé)崞鲀?nèi)部輻射換熱，對(duì)絲網(wǎng)結(jié)構(gòu)及材料物性進(jìn)行數(shù)值模擬，得到絲網(wǎng)材料比熱容的最優(yōu)值。Dellali 等人[13]對(duì)振蕩流條件下，柱陣列回?zé)崞鞯膲航颠M(jìn)行實(shí)驗(yàn)研究，推導(dǎo)了回?zé)崞髂Σ料禂?shù)的實(shí)驗(yàn)關(guān)聯(lián)式。目前，對(duì)于斯特林發(fā)動(dòng)機(jī)的研究集中在整機(jī)循環(huán)特性實(shí)驗(yàn)研究、斯特林循環(huán)分析方法、回?zé)崞髁鲃?dòng)與換熱特性研究等各個(gè)方面。

在斯特林發(fā)動(dòng)機(jī)運(yùn)行過(guò)程中，工質(zhì)在各個(gè)部件間循環(huán)流動(dòng)，流動(dòng)過(guò)程中不可避免地會(huì)發(fā)生射流現(xiàn)象，并產(chǎn)生一定的功損失和熱損失。尤其是在回?zé)崞髋c加熱器、冷卻器接口部分的間隙內(nèi)產(chǎn)生的射流現(xiàn)象，直接影響了工質(zhì)在回?zé)崞鲀?nèi)的流動(dòng)與換熱情況?；谌S數(shù)值模擬，本文研究了斯特林發(fā)動(dòng)機(jī)內(nèi)的射流現(xiàn)象對(duì)斯特林循環(huán)的影響及損失機(jī)制，并對(duì)回?zé)崞鞑糠诌M(jìn)行了結(jié)構(gòu)優(yōu)化設(shè)計(jì)，最后依托1 臺(tái)百瓦級(jí)β型斯特林發(fā)動(dòng)機(jī)，對(duì)不同工況下的模擬結(jié)果進(jìn)行了實(shí)驗(yàn)驗(yàn)證。

1 斯特林發(fā)動(dòng)機(jī)三維數(shù)值模型

基于Ansys CFD 軟件，對(duì)斯特林發(fā)動(dòng)機(jī)內(nèi)復(fù)雜的流動(dòng)換熱情況進(jìn)行三維數(shù)值模擬。使用Ansys DesignModeler 建立幾何模型，并創(chuàng)建流體計(jì)算域?；贗CEM CFD 進(jìn)行了全結(jié)構(gòu)化網(wǎng)格的生成。利用Ansys Fluent 18.0 軟件進(jìn)行數(shù)值計(jì)算。

1.1 幾何模型與網(wǎng)格劃分

考慮到幾何與流場(chǎng)的對(duì)稱性，為了提高計(jì)算精度與計(jì)算效率，構(gòu)造了1/4 幾何模型如圖1所示。對(duì)照β型斯特林發(fā)動(dòng)機(jī)的設(shè)計(jì)制造參數(shù)，確定了流體域的幾何參數(shù)，并進(jìn)行了全結(jié)構(gòu)化網(wǎng)格劃分，結(jié)果如圖2所示。

圖1 斯特林發(fā)動(dòng)機(jī)幾何模型（m）Fig.1 Geometric model of Stirling engine(m)

圖2 斯特林發(fā)動(dòng)機(jī)網(wǎng)格劃分Fig.2 Mesh generation for the Stirling engine

在相同的幾何模型上，生成了不同尺寸和數(shù)量的3 種網(wǎng)格，其網(wǎng)格數(shù)分別為153.6 萬(wàn)、506.8 萬(wàn)和1 247.6 萬(wàn)。計(jì)算過(guò)程中發(fā)現(xiàn)，3 種網(wǎng)格計(jì)算所得的溫度、壓力、循環(huán)效率和循環(huán)指示功均有良好的收斂性。在其他條件相同的情況下，循環(huán)效率和循環(huán)指示功隨網(wǎng)格數(shù)目的變化如圖3所示。從圖3 可以看到，506.8 萬(wàn)網(wǎng)格與1 247.6 萬(wàn)網(wǎng)格之間循環(huán)效率相對(duì)誤差0.75%，循環(huán)輸出功相對(duì)誤差0.87%。綜合考慮計(jì)算精度與計(jì)算所需資源，本文選擇506.8 萬(wàn)網(wǎng)格進(jìn)行計(jì)算。

圖3 網(wǎng)格無(wú)關(guān)性驗(yàn)證Fig.3 Grid independence verification

1.2 三維數(shù)值模擬求解器設(shè)置

綜合考慮百瓦級(jí)β 型斯特林發(fā)動(dòng)機(jī)的幾何構(gòu)造、多孔介質(zhì)區(qū)域以及運(yùn)動(dòng)部件動(dòng)網(wǎng)格，以氮?dú)鉃楣べ|(zhì)，將求解器類型設(shè)置為瞬態(tài)、壓力基，選擇SIMPLE 算法，選擇具有相當(dāng)精度的Standardk-ε湍流模型。對(duì)應(yīng)實(shí)驗(yàn)工況參數(shù)設(shè)置三維數(shù)值模擬的邊界條件見(jiàn)表1。

表1 三維數(shù)值模擬邊界條件Tab.1 Boundary conditions for three dimensional numerical simulation

本模型是利用流體計(jì)算域邊界的運(yùn)動(dòng)，模擬實(shí)際斯特林發(fā)動(dòng)機(jī)中配氣活塞和動(dòng)力活塞的運(yùn)動(dòng)，使用動(dòng)態(tài)層方法（Laying）對(duì)邊界上的網(wǎng)格進(jìn)行更新，同時(shí)利用用戶自定義函數(shù)（UDF）指定邊界運(yùn)動(dòng)，結(jié)合斯特林發(fā)動(dòng)機(jī)菱形傳動(dòng)設(shè)計(jì)參數(shù)和DEFINE_CG_MOTION 宏，定義模型中配氣活塞和動(dòng)力活塞的運(yùn)動(dòng)。

該模型利用Fluent 軟件自帶的多孔介質(zhì)模型，模擬實(shí)際工質(zhì)在金屬絲網(wǎng)回?zé)崞鞑糠值牧鲃?dòng)狀況。多孔介質(zhì)模型簡(jiǎn)化了多孔結(jié)構(gòu)，求解時(shí)在動(dòng)量方程中增添由動(dòng)量損失項(xiàng)和慣性損失項(xiàng)組成的動(dòng)量源項(xiàng)。

式中，Si代表i方向動(dòng)量源項(xiàng)，代表黏性阻力代表慣性阻力。引入Tanaka 等人[14]關(guān)于回?zé)崞髁髯柘禂?shù)f的經(jīng)驗(yàn)關(guān)聯(lián)式

最終計(jì)算得到主流方向上的黏性阻力系數(shù)Cf和慣性阻力系數(shù)Cj：

式中，φ為孔隙率，dh為金屬絲網(wǎng)水力直徑。

對(duì)于簡(jiǎn)化為多孔介質(zhì)模型的回?zé)崞鱾鳠徇^(guò)程，本模型基于非熱平衡假設(shè)，利用用戶自定義標(biāo)量方程（UDS）定義額外的標(biāo)量輸運(yùn)方程，并計(jì)算流固間換熱過(guò)程，其中固體能量方程如下：

式中，h為流固間對(duì)流換熱系數(shù)，ρs為固體密度，Es為固體介質(zhì)總能，ks為固體介質(zhì)導(dǎo)熱系數(shù)，Ts為固體介質(zhì)溫度，Tf為流體溫度，hfs為流固間對(duì)流換熱系數(shù)，為源項(xiàng)。

引入Gedeon 等人[15]采用瞬時(shí)局部雷諾數(shù)描述的換熱關(guān)聯(lián)式。

通過(guò)分別計(jì)算固相非穩(wěn)態(tài)項(xiàng)、固體熱擴(kuò)散相、考慮徑向熱擴(kuò)散的流體熱彌散項(xiàng)和能量方程源項(xiàng)，最終得到多孔介質(zhì)區(qū)域流固耦合換熱情況。

1.3 斯特林發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)與換熱特性

基于斯特林發(fā)動(dòng)機(jī)三維數(shù)值模擬，對(duì)特定工況（N2-1.88MPa）和曲軸某一特定角度（33.26°）下，斯特林發(fā)動(dòng)機(jī)內(nèi)部溫度場(chǎng)和流場(chǎng)進(jìn)行了研究。在曲軸轉(zhuǎn)角為33.26°下，斯特林發(fā)動(dòng)機(jī)膨脹腔體積接近最小值，正處在壓縮腔體積減小，膨脹腔體積增大，工質(zhì)由壓縮腔流向膨脹腔的所謂“冷吹”過(guò)程。

斯特林發(fā)動(dòng)機(jī)內(nèi)溫度分布如圖4所示。由圖4可見(jiàn)：膨脹腔和加熱器內(nèi)工質(zhì)溫度較高，壓縮腔和冷卻器內(nèi)溫度較低；回?zé)崞鞑糠譁囟扔蔁岫说嚼涠搜剌S向下降，同時(shí)在冷吹過(guò)程中，回?zé)崞鲀?nèi)部工質(zhì)流向加熱器被加熱。

圖4 斯特林發(fā)動(dòng)機(jī)內(nèi)溫度分布Fig.4 Temperature distributions in the Stirling engine

圖5 為斯特林發(fā)動(dòng)機(jī)內(nèi)速度流線。

圖5 斯特林發(fā)動(dòng)機(jī)內(nèi)速度流線Fig.5 Velocity streamlines in the Stirling engine

由圖5 可以看出，在“冷吹”過(guò)程中，加熱器和回?zé)崞?、冷卻器和壓縮腔之間的連接處和加熱器彎管處，工質(zhì)流速最大，膨脹腔和壓縮腔內(nèi)工質(zhì)流速較小，回?zé)崞鲀?nèi)工質(zhì)的流速最小。在“冷吹”過(guò)程中，配氣活塞和動(dòng)力活塞的運(yùn)動(dòng)迫使工質(zhì)從加熱器流向膨脹腔，由于橫截面積的變化，在加熱器和膨脹腔接口處以及冷卻器與回?zé)崞鹘涌谔幑べ|(zhì)流速迅速增大，產(chǎn)生較為明顯的射流現(xiàn)象。

2 斯特林發(fā)動(dòng)機(jī)內(nèi)射流現(xiàn)象與回?zé)崞鹘Y(jié)構(gòu)優(yōu)化

2.1 膨脹腔內(nèi)射流現(xiàn)象及其影響

在“冷吹”過(guò)程中，加熱器內(nèi)高溫工質(zhì)進(jìn)入膨脹腔內(nèi)并產(chǎn)生較為明顯的射流現(xiàn)象，膨脹腔內(nèi)密度分布如圖6所示。由圖6 可見(jiàn)，膨脹腔內(nèi)密度分布高度不均勻。由射流現(xiàn)象導(dǎo)致的密度分布，與斯特林發(fā)動(dòng)機(jī)內(nèi)部平均壓力密切相關(guān)。模擬結(jié)果顯示，在整個(gè)斯特林循環(huán)過(guò)程中的平均壓力越大（1.41、1.88、2.36 MPa）參與做功的工質(zhì)質(zhì)量越大，同時(shí)平均壓力的增大使密度分布的不均勻性更加明顯，膨脹腔內(nèi)的工質(zhì)密度差也越大（0.372、0.507、0.642 kg/m3）。平均壓力的增大顯著增強(qiáng)了加熱器高溫工質(zhì)射流進(jìn)入膨脹腔內(nèi)的傳熱過(guò)程，增強(qiáng)了冷熱工質(zhì)混合作用。

圖6 膨脹腔內(nèi)密度分布Fig.6 Density distributions in the expansion chamber

對(duì)于膨脹腔內(nèi)的射流現(xiàn)象，另一個(gè)重要的影響因素則是轉(zhuǎn)速。轉(zhuǎn)速的改變直接影響了斯特林發(fā)動(dòng)機(jī)工作腔中工質(zhì)的流速，膨脹腔內(nèi)工質(zhì)速度矢量分布如圖7所示。

圖7 膨脹腔內(nèi)工質(zhì)速度矢量分布Fig.7 Velocity vectors in the expansion chamber

由圖7 可見(jiàn)，在加熱器與膨脹腔接口部分，工質(zhì)流速最大，同時(shí)在這部分高速工質(zhì)和膨脹腔壁面之間，存在速度較小的漩渦。由模擬結(jié)果得到，隨著轉(zhuǎn)速的增大（600、800、1 000 r/min），由加熱器流向膨脹腔的工質(zhì)的流速不斷增大（最大流速分別為7.64、10.04、12.44 m/s），工質(zhì)的射流沖擊作用更加強(qiáng)烈，產(chǎn)生的漩渦及能量耗散也更多。

2.2 回?zé)崞黜敹碎g隙內(nèi)射流現(xiàn)象及回?zé)崞鹘Y(jié)構(gòu)優(yōu)化

在回?zé)崞髋c加熱器、冷卻器的連接部分存在一定大小的間隙。以“冷吹”過(guò)程為例，工質(zhì)由冷卻器流向回?zé)崞?，在絲網(wǎng)頂面間隙內(nèi)發(fā)生射流現(xiàn)象（圖8）。

圖8 回?zé)崞鲀?nèi)射流現(xiàn)象示意Fig.8 Jet impingement in the regenerator

定義回?zé)崞鏖g隙長(zhǎng)度比λ為

式中，d為絲網(wǎng)頂面間隙長(zhǎng)度，L為回?zé)崞鏖L(zhǎng)度。當(dāng)前斯特林發(fā)動(dòng)機(jī)回?zé)崞鏖g隙長(zhǎng)度比λ設(shè)計(jì)值為0.074。為進(jìn)一步研究回?zé)崞黜敹碎g隙內(nèi)射流沖擊對(duì)斯特林循環(huán)特性的影響，分別設(shè)置6 組不同λ值0、0.037、0.074、0.130、0.167 和0.259，進(jìn)行三維數(shù)值模擬。不同回?zé)崞鏖g隙長(zhǎng)度比下循環(huán)效率變化如圖9所示。

圖9 不同回?zé)崞鏖g隙長(zhǎng)度比下循環(huán)效率變化Fig.9 Changes of the cycle efficiency at different ratios of regenerator gap length

由圖9 可見(jiàn)，當(dāng)λ值范圍為0.037～0.167 時(shí)循環(huán)效率較大，與回?zé)崞黜敹藷o(wú)間隙的情況（λ=0）相比，當(dāng)λ=0.074 時(shí)，循環(huán)效率提高了8.32%。

不同回?zé)崞鏖g隙長(zhǎng)度比下工質(zhì)流速分布如圖10所示。

圖10 不同回?zé)崞鏖g隙長(zhǎng)度比下工質(zhì)流速分布Fig.10 Velocity distributions at different ratios of regenerator gap length

由圖10 可見(jiàn)：在冷吹過(guò)程中，當(dāng)λ=0，即回?zé)崞黜敹藷o(wú)間隙、無(wú)射流現(xiàn)象時(shí)，工質(zhì)直接從冷卻器流入回?zé)崞?，并只流過(guò)回?zé)崞鞯囊恍〔糠只w，使回?zé)崞鞯幕責(zé)嵝Ч艿揭欢ㄓ绊懀ㄑh(huán)效率η=9.50%）；當(dāng)λ=0.074 時(shí)，工質(zhì)在回?zé)崞黜敹碎g隙內(nèi)發(fā)生射流現(xiàn)象，在較廣的區(qū)域內(nèi)流入回?zé)崞鳎ㄑh(huán)效率η=10.29%）；當(dāng)回?zé)崞黜敹碎g隙進(jìn)一步增大，λ=0.259 時(shí)，在回?zé)崞黜敹碎g隙內(nèi)出現(xiàn)了較為明顯的漩渦，工質(zhì)做功能力下降（循環(huán)效率η為8.90%）。

回?zé)崞黜敹碎g隙內(nèi)射流現(xiàn)象的存在，一方面使工質(zhì)更加充分地流入回?zé)崞?，提高了回?zé)嵝?；另一方面使間隙內(nèi)產(chǎn)生漩渦，造成工質(zhì)做功能力的下降。因此，以提高循環(huán)效率作為回?zé)崞鹘Y(jié)構(gòu)優(yōu)化目標(biāo)時(shí)，該斯特林發(fā)動(dòng)機(jī)存在最佳回?zé)崞鏖g隙長(zhǎng)度比λ，其值范圍為0.037～0.130。

3 百瓦級(jí)β型斯特林發(fā)動(dòng)機(jī)實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)系統(tǒng)組成及測(cè)點(diǎn)布置

實(shí)驗(yàn)系統(tǒng)由加熱裝置、百瓦級(jí)β型斯特林發(fā)動(dòng)機(jī)主體以及數(shù)據(jù)采集裝置3 部分組成。百瓦級(jí)β型斯特林發(fā)動(dòng)機(jī)主要由飛輪、菱形傳動(dòng)機(jī)構(gòu)、動(dòng)力活塞、配氣活塞、冷卻器、回?zé)崞骱图訜崞鹘M成。實(shí)驗(yàn)測(cè)點(diǎn)布置及三維結(jié)構(gòu)如圖11所示。

圖11 β型斯特林發(fā)動(dòng)機(jī)結(jié)構(gòu)示意Fig.11 Schematic diagram of β type Stirling engine

以氮?dú)鉃楣べ|(zhì)，采取鹽浴方式加熱并使用冷卻水進(jìn)行冷卻，在一定轉(zhuǎn)速和充氣壓力范圍內(nèi)進(jìn)行實(shí)驗(yàn)。具體實(shí)驗(yàn)設(shè)計(jì)工況見(jiàn)表2。

表2 實(shí)驗(yàn)臺(tái)實(shí)驗(yàn)工況Tab.2 Experimental working conditions for test bench

3.2 模擬結(jié)果實(shí)驗(yàn)驗(yàn)證

在本文中，利用循環(huán)指示功Win，循環(huán)輸出功率P和循環(huán)效率η3 個(gè)指標(biāo)來(lái)衡量斯特林循環(huán)特性。Win代表1 個(gè)循環(huán)周期內(nèi)，即曲軸旋轉(zhuǎn)1 周，工質(zhì)完成1 次膨脹壓縮的過(guò)程中，斯特林發(fā)動(dòng)機(jī)對(duì)外輸出功。pi代表1 個(gè)循環(huán)周期內(nèi)采集到的壓力信號(hào)，Vi代表對(duì)應(yīng)的工作腔體積，n代表1 個(gè)循環(huán)周期內(nèi)采集到的數(shù)據(jù)組數(shù)。在1 個(gè)循環(huán)周期內(nèi)，對(duì)n組壓力平均值與對(duì)應(yīng)的體積差值的乘積進(jìn)行累加求和，得到Win。

式中P代表在單位時(shí)間內(nèi)，斯特林發(fā)動(dòng)機(jī)的對(duì)外輸出功，可根據(jù)轉(zhuǎn)速υ和Win得到。

η綜合衡量了斯特林循環(huán)過(guò)程中流動(dòng)與換熱損失，利用冷卻功率Qco和循環(huán)輸出功W可以得到[16]。

式中，c為冷卻水的比熱容，q為冷卻水流量，ΔTco為冷卻器進(jìn)出口水的溫差。

在實(shí)驗(yàn)條件下數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比如圖12、圖13所示。由圖12、圖13 可見(jiàn)：在轉(zhuǎn)速相同且壓力不同的情況下，Win的相對(duì)誤差范圍為3.00%～5.78%，η的相對(duì)誤差范圍為5.89%～9.91%；在壓力相同且轉(zhuǎn)速不同的情況下，Win的三維數(shù)值模擬結(jié)果與實(shí)驗(yàn)值的相對(duì)誤差范圍為 3.58%～7.72%，η的相對(duì)誤差范圍為5.48%～9.89%?？梢?jiàn)，模擬結(jié)果與實(shí)驗(yàn)結(jié)果一致性較好。

圖12 不同壓力下數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比Fig.12 Comparison between the numerical model results and experimental results at different pressures

圖13 不同轉(zhuǎn)速下數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比Fig.13 Comparison between the numerical model results and experimental results at different rotating speeds

由圖12、圖13 還可見(jiàn)：在實(shí)驗(yàn)條件下，壓力增大，循環(huán)指示功和循環(huán)效率越大；轉(zhuǎn)速增大，循環(huán)指示功減小，循環(huán)效率增大但增長(zhǎng)幅度減小。

4 結(jié)論

1）三維模型與實(shí)驗(yàn)結(jié)果具有較好的一致性，其結(jié)果均顯示，對(duì)于百瓦級(jí)β型斯特林發(fā)動(dòng)機(jī)，在一定范圍內(nèi)增大平均壓力，可增大斯特林發(fā)動(dòng)機(jī)中工質(zhì)質(zhì)量，同時(shí)增強(qiáng)射流現(xiàn)象所導(dǎo)致的工質(zhì)分布不均勻，以及冷熱流體混合和傳熱過(guò)程，從而增大循環(huán)指示功，提升循環(huán)效率。增大斯特林發(fā)動(dòng)機(jī)轉(zhuǎn)速，會(huì)增大工質(zhì)流速，使得工質(zhì)射流過(guò)程更加明顯，漩渦產(chǎn)生的能量耗散更大，從而減小循環(huán)指示功，增大循環(huán)輸出功，循環(huán)效率增大但增幅較小。

2）在工質(zhì)由壓縮腔流向膨脹腔的“冷吹”過(guò)程中，加熱器內(nèi)工質(zhì)進(jìn)入膨脹腔發(fā)生射流現(xiàn)象，冷熱工質(zhì)混合導(dǎo)致膨脹腔溫度分布不均勻，并產(chǎn)生沖擊作用和漩渦，使工質(zhì)做功能力下降。

3）回?zé)崞黜敹碎g隙內(nèi)射流現(xiàn)象影響了工質(zhì)在回?zé)崞鲀?nèi)的流動(dòng)與換熱，回?zé)崞鏖g隙長(zhǎng)度比λ存在1 個(gè)最佳范圍，使斯特林發(fā)動(dòng)機(jī)循環(huán)效率最大。在本研究中λ由0 增至0.037～0.130 時(shí)，對(duì)應(yīng)的最高循環(huán)效率比提高了8.32%。