劉柏岙，王樹林，肖 剛

（能源清潔利用國家重點實驗室浙江大學能源工程學院，浙江 杭州 310027）

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

段晨[9]建立了考慮“不可逆因子”的斯特林循環(huán)熱力學模型，采用相似設計和多目標優(yōu)化方法設計制造了1 臺β型斯特林發(fā)動機。瞿凡等[10]基于1 臺β型斯特林發(fā)動機，對加熱器進行了輻射與燃燒2 種加熱方式的數(shù)值模擬與優(yōu)化設計。黃護林等[11]對斯特林發(fā)動機回熱器進行數(shù)值模擬，研究了長徑比和孔隙率對回熱器有效性及流阻損失的影響。王思[12]考慮回熱器內(nèi)部輻射換熱，對絲網(wǎng)結(jié)構(gòu)及材料物性進行數(shù)值模擬，得到絲網(wǎng)材料比熱容的最優(yōu)值。Dellali 等人[13]對振蕩流條件下，柱陣列回熱器的壓降進行實驗研究，推導了回熱器摩擦系數(shù)的實驗關(guān)聯(lián)式。目前，對于斯特林發(fā)動機的研究集中在整機循環(huán)特性實驗研究、斯特林循環(huán)分析方法、回熱器流動與換熱特性研究等各個方面。

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

1 斯特林發(fā)動機三維數(shù)值模型

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

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

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

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

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

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

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

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

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

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

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

該模型利用Fluent 軟件自帶的多孔介質(zhì)模型，模擬實際工質(zhì)在金屬絲網(wǎng)回熱器部分的流動狀況。多孔介質(zhì)模型簡化了多孔結(jié)構(gòu)，求解時在動量方程中增添由動量損失項和慣性損失項組成的動量源項。

式中，Si代表i方向動量源項，代表黏性阻力代表慣性阻力。引入Tanaka 等人[14]關(guān)于回熱器流阻系數(shù)f的經(jīng)驗關(guān)聯(lián)式

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

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

對于簡化為多孔介質(zhì)模型的回熱器傳熱過程，本模型基于非熱平衡假設，利用用戶自定義標量方程（UDS）定義額外的標量輸運方程，并計算流固間換熱過程，其中固體能量方程如下：

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

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

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

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

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

斯特林發(fā)動機內(nèi)溫度分布如圖4所示。由圖4可見：膨脹腔和加熱器內(nèi)工質(zhì)溫度較高，壓縮腔和冷卻器內(nèi)溫度較低；回熱器部分溫度由熱端到冷端沿軸向下降，同時在冷吹過程中，回熱器內(nèi)部工質(zhì)流向加熱器被加熱。

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

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

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

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

2 斯特林發(fā)動機內(nèi)射流現(xiàn)象與回熱器結(jié)構(gòu)優(yōu)化

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

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

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

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

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

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

2.2 回熱器頂端間隙內(nèi)射流現(xiàn)象及回熱器結(jié)構(gòu)優(yōu)化

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

圖8 回熱器內(nèi)射流現(xiàn)象示意Fig.8 Jet impingement in the regenerator

定義回熱器間隙長度比λ為

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

圖9 不同回熱器間隙長度比下循環(huán)效率變化Fig.9 Changes of the cycle efficiency at different ratios of regenerator gap length

由圖9 可見，當λ值范圍為0.037～0.167 時循環(huán)效率較大，與回熱器頂端無間隙的情況（λ=0）相比，當λ=0.074 時，循環(huán)效率提高了8.32%。

不同回熱器間隙長度比下工質(zhì)流速分布如圖10所示。

圖10 不同回熱器間隙長度比下工質(zhì)流速分布Fig.10 Velocity distributions at different ratios of regenerator gap length

由圖10 可見：在冷吹過程中，當λ=0，即回熱器頂端無間隙、無射流現(xiàn)象時，工質(zhì)直接從冷卻器流入回熱器，并只流過回熱器的一小部分基體，使回熱器的回熱效果受到一定影響（循環(huán)效率η=9.50%）；當λ=0.074 時，工質(zhì)在回熱器頂端間隙內(nèi)發(fā)生射流現(xiàn)象，在較廣的區(qū)域內(nèi)流入回熱器（循環(huán)效率η=10.29%）；當回熱器頂端間隙進一步增大，λ=0.259 時，在回熱器頂端間隙內(nèi)出現(xiàn)了較為明顯的漩渦，工質(zhì)做功能力下降（循環(huán)效率η為8.90%）。

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

3 百瓦級β型斯特林發(fā)動機實驗驗證

3.1 實驗系統(tǒng)組成及測點布置

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

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

以氮氣為工質(zhì)，采取鹽浴方式加熱并使用冷卻水進行冷卻，在一定轉(zhuǎn)速和充氣壓力范圍內(nèi)進行實驗。具體實驗設計工況見表2。

表2 實驗臺實驗工況Tab.2 Experimental working conditions for test bench

3.2 模擬結(jié)果實驗驗證

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

3）回熱器頂端間隙內(nèi)射流現(xiàn)象影響了工質(zhì)在回熱器內(nèi)的流動與換熱，回熱器間隙長度比λ存在1 個最佳范圍，使斯特林發(fā)動機循環(huán)效率最大。在本研究中λ由0 增至0.037～0.130 時，對應的最高循環(huán)效率比提高了8.32%。