陳 曦，崔 浩

（上海理工大學能源與動力工程學院，上海 200093）

0 引言

1816年，蘇格蘭牧師Robert Stirling發(fā)明了斯特林發(fā)動機（又稱熱氣機），是一種外燃、閉式循環(huán)、往復活塞式的能量轉換裝置[1]。20世紀60年代，俄亥俄大學機械工程學院教授William Beale改進斯特林發(fā)動機結構，發(fā)明了自由活塞斯特林發(fā)動機（FPSE），具有效率高、壽命長、結構簡單、噪聲低、不易磨損和可自啟動等優(yōu)點[2]。

按Martini的命名規(guī)則，斯特林發(fā)動機熱力學分析方法主要分為五類：零級分析法、一級分析法、二級分析法、三級分析法和四級分析法。零級分析法，即實驗模型，斯特林發(fā)動機輸出功率的估算公式來源于大量實驗數(shù)據(jù)。目前為止主要有三種：Malmo公式法、指示功率法和Beale數(shù)法[3]。一級分析法又稱為等溫分析法，最主要的假設是熱腔和冷腔內(nèi)工質(zhì)的循環(huán)溫度恒定，通過理論推導一般可得到解析解，進而可考慮各種熱損失和流動損失。一級分析法首先由Schmidt[4]完成。Schmidt并沒考慮熱損失和流動損失，因此，在實際應用時需要對Schmidt分析結果進行修正，修正系數(shù)一般為0.45～0.55[5]。一級分析法存在比較大的理論誤差，一般只適用于定性分析。二級分析法是對每一個腔體內(nèi)質(zhì)量方程、能量方程和氣體狀態(tài)方程進行求解，進而可計算各種功率損失和熱損失，最終得到斯特林發(fā)動機的輸出功率和熱效率。二級分析法有三類：等溫模型、絕熱模型和半絕熱模型。等溫模型、絕熱模型更接近實際情況，而且方程簡單易懂，計算簡便，是二級分析方法中最常用的模型[6]，但是二級分析法中是假設功率損失和熱損失之間沒有關聯(lián)，故而該方法并沒有對損失機理進行深入研究。Finkelestin[7]首先提出三級分析法，該方法考慮了氣體在各腔體內(nèi)的一維流動，采用節(jié)點法列出每個節(jié)點處工質(zhì)的質(zhì)量守恒、動量守恒和能量守恒偏微分方程，采用數(shù)值求解，建立斯特林發(fā)動機瞬時能量和工質(zhì)流動的模型，精確地模擬斯特林發(fā)動機輸出功率和熱效率?？紤]各損失之間的相互影響，有助于對損失機理的研究。Gedeon[8]將三級分析法編成軟件GLIMPS，然后又改進GLIMPS的代碼編制了Sage軟件[9]，三級分析法解決了一級分析法和二級分析法的空間誤差問題，得到了最廣泛的發(fā)展和應用，是目前發(fā)展最成熟的斯特林循環(huán)分析法。四級分析法即CFD分析，考慮了在一維模型中忽略的橫向梯度與其他多維現(xiàn)象，并可以得到斯特林發(fā)動機工作過程中各種損失之間的相互影響，有助于更深入地了解斯特林發(fā)動機的傳熱特性和工質(zhì)流動特性。

采用三級分析軟件Sage，結合了Re-1000的結構參數(shù)和運行參數(shù)，建立了Re-1000的一維Sage模型，并將Sage模擬結果與實驗數(shù)據(jù)進行對比。

1 自由活塞斯特林發(fā)動機及參數(shù)

自由活塞斯特林發(fā)動機（FPSE），主要包括：配氣活塞、動力活塞、加熱器、回熱器和冷卻器等，物理模型如圖1所示。

圖1 自由活塞斯特林發(fā)動機的工作原理圖Fig.1 The principle diagram of a free-piston Stirling engine

1979年，美國Sunpower公司和俄亥俄大學等設計并制造了Re-1000斯特林發(fā)動機，以供NASA Lewis研究中心研究使用[10]，斯特林發(fā)動機的主要結構參數(shù)如表1所列。Re-1000機器為自由活塞斯特林發(fā)動機，經(jīng)過多年的機器實驗，科研人員記錄了大量數(shù)據(jù)用于評估自由活塞斯特林發(fā)動機的熱力學性能，如工作腔平均壓力，加熱器溫度，冷卻器溫度，回熱器空隙率等參數(shù)對發(fā)動機性能的影響[11-12]。這些實驗數(shù)據(jù)也可用于驗證編寫計算機程序的準確性。NASA報告TM-82999、TM-83407和TM-87126給出了該機器的原始實驗數(shù)據(jù)。為驗證模擬的準確性，將模擬結果與Re-1000機器實驗數(shù)據(jù)進行對比。

表1 Re-1000自由活塞斯特林發(fā)動機的結構參數(shù)Table1 Structural parameters of the Re-1000 free-piston Stirling engine

2 模擬軟件及一維模擬模型

Sage軟件中斯特林模塊是斯特林整機模擬軟件[13-14]，可以模擬不同類型的斯特林循環(huán)的發(fā)動機與制冷機。不論是那種斯特林機械，都是由換熱固體壁面、氣體區(qū)域、圓筒、換熱器、活塞等部件組成，各部件通過合適的熱流邊界、力的作用面相互耦合。主要利用斯特林模塊來模擬分析斯特林發(fā)動機。

2.1 氣體域求解

Sage軟件會基于所建立的發(fā)動機模型，在空間節(jié)點上采用中心差分法，在時間節(jié)點上采用后三點差分法，把發(fā)動機的控制方程離散成一個大型的非線性方程組，然后采用非線性解算器對這個方程組進行迭代求解，得到發(fā)動機模型的輸出參數(shù)，其求解過程非常快，求解所用的時間一般都在幾分鐘之內(nèi)。

Sage模擬的核心是氣體域的求解。氣體區(qū)域可以簡化為具有長度L、流動截面積A、濕周Sx的一維矩形域，如圖2所示。質(zhì)量流率m（t）僅通過x正負方向兩端邊界與其他氣體區(qū)域相連，熱流q（x，t）僅在z方向負邊界與熱固體壁面相連。在軟件運行求解中，氣體區(qū)域僅在x方向進行離散，每個固體區(qū)域內(nèi)都采用均勻網(wǎng)格劃分，其中網(wǎng)格數(shù)的劃分可根據(jù)各部件內(nèi)動態(tài)參數(shù)變化的劇烈程度自行設定。且在x方向上為變截面，即A（x）。當z方向上表面隨位移變化時，該截面積為時間的函數(shù)，即A（x，t）。

圖2 氣體區(qū)域模型示意圖Fig.2 Gas area model

假定控制容積進出口邊界固定，允許側邊界為時均流截面，且側邊界為非滲透，則流量僅通過進出口邊界而不通過側邊界。忽略該氣體區(qū)域的體積力，則該控制容積的氣體控制方程可以簡化為：

連續(xù)性方程：

動量方程：

能量方程：式中：t為時間；ρ為氣體的密度；u為速度；e為質(zhì)量能（包括內(nèi)能與動能）；A為氣體區(qū)域的截面積；p為壓力；F為黏性壓力梯度；q為氣體軸向導熱熱流；Qw為對流換熱熱流。由于發(fā)動機不同部件中的氣體的流動換熱特性不同，將氣體區(qū)域分為三類：多孔絲網(wǎng)換熱流動區(qū)域（Matrix gas）、氣缸內(nèi)活塞運動導致的變?nèi)莘e腔體內(nèi)的氣體區(qū)域（Cylinder-space gas）和管道流動氣體區(qū)域（Duct gas）。Qw表示單位長度氣體區(qū)域通過z軸負邊界表面與固體壁面間的換熱量，不同類型的氣體區(qū)域由于流動換熱的特性不同，Qw的計算公式也不同。

（1）絲網(wǎng)氣體區(qū)域（Matrix gas）

因為絲網(wǎng)區(qū)域中氣體的熱滲透深度大于絲網(wǎng)絲徑，因而氣體和絲網(wǎng)間的換熱與壁面溫度變化同相。

式中：Nu為努賽爾數(shù)；k為氣體導熱系數(shù)；dh為水力直徑；Sx為濕周（單位長度的濕面積）；（Tw-T）為z負向表面與平均截面間的溫差。

（2）氣缸氣體區(qū)域（Cylinder-space gas）

在氣缸氣體區(qū)域中，氣體和絲網(wǎng)間的換熱與壁面溫度變化存在相位差，引入努賽爾數(shù)的復數(shù)形式，其計算可以通過線性時均層流理論得到。

（3）管道氣體區(qū)域（Duct gas）

管道區(qū)域內(nèi)的壓縮和對流都會引起氣體的溫度波動，前者是由于壓力變化，后者是由于流動過程中的溫度梯度。壓縮引起的溫度變化大約是對流引起的溫度變化的兩倍。兩種不同機理引起的溫度波動造成了努賽爾數(shù)計算的差異。

為了迭代求解方程，還需要得出氣體區(qū)域流動幾何的摩擦因子f、努賽爾數(shù)Nu、軸向導熱率Nk的經(jīng)驗公式。

2.2 一維Sage模型

Sage軟件是通過非常有邏輯的樹狀結構來對發(fā)動機進行建模的。其建模采用圖形化界面，可選取合適的模塊來模擬斯特林發(fā)動機的不同部件。雖然各部件的結構和工作機理不同，但在模擬中都可以看成工質(zhì)在各種形狀結構的流道內(nèi)的換熱和流動。將不同固體模塊和氣體模塊進行組合，即可實現(xiàn)斯特林發(fā)動機不同部件的模擬。不同模塊間根據(jù)實際情況再進行質(zhì)量流、能量流連接，最終完成整機建模。

根據(jù)Re-1000的結構參數(shù)和運行參數(shù)，建立了Re-1000的一維Sage模型，模型根目錄如圖3所示。Sage建模中各部件的模型是以分層樹狀的結構有邏輯地連接在一起的。該模型有幾點基本假設[15-16]：（1）工質(zhì)氣體為氦氣，并采用理想氣體模型；（2）熱物性參數(shù)為溫度的函數(shù)；（3）氣體流動和傳熱為一維；（4）回熱器內(nèi)填料不可壓縮，各截面上填料空隙率不變；（5）忽略輻射漏熱和壁面與環(huán)境的對流換熱。

圖3 Re-1000機器的Sage模型圖Fig.3 Sage model of Re-1000 machine

3 模擬結果及驗證

根據(jù)自由活塞斯特林發(fā)動機樣機Re-1000的結構參數(shù)和運行參數(shù)，利用Sage軟件模擬計算斯特林機器的輸出功率和熱效率，并與Re-1000實驗數(shù)據(jù)進行對比。模擬時工質(zhì)為He，冷卻溫度和加熱溫度保持不變，分別為298 K和873 K，工作頻率為30.2 Hz。輸出功率和熱效率隨平均壓力的變化如圖4所示?？梢钥闯觯敵龉β?、熱效率的模擬值和實驗值呈相同的變化趨勢。隨著壓力的增大，指示功在增大。而熱效率隨著壓力的增大先增大后減小，即存在最佳充氣壓力，使得自由活塞斯特林發(fā)動機的熱效率最大。輸出功率模擬值與實驗值的誤差在15%左右，熱效率模擬值比實驗值大4%。

產(chǎn)生誤差的主要原因包括兩個方面，一方面加熱器和冷卻器的內(nèi)外壁溫差沒有考慮。Re-1000樣機實驗數(shù)據(jù)中給出的是加熱器和冷卻器外壁溫度，并未給出內(nèi)壁溫度。因換熱器內(nèi)外壁存在溫差，這將導致加熱器內(nèi)壁溫度低于外壁溫度，冷卻器內(nèi)壁溫度高于外壁溫度，即模擬中機器冷熱源溫差要比實際溫差要大，這將導致指示功和熱效率模擬結果偏大。另一方面，樣機參數(shù)并不全面，間隙密封尺寸并沒有給出具體數(shù)值等，這些因素都會導致模擬值與實驗值存在偏差。

圖4 輸出功率和熱效率隨充氣壓力的變化曲線Fig.4 Output power and thermal efficiency vs.pressure

輸出功率和熱效率隨加熱溫度的變化如圖5所示。此時冷卻溫度為298 K，充氣壓力為7 MPa，工質(zhì)為He，工作頻率為30.2 Hz。從圖5中可以看出，輸出功率、熱效率的模擬值和實驗值呈相同的變化趨勢。隨著加熱溫度的升高，輸出功率和熱效率均增大；輸出功率的模擬值和實驗值誤差15%左右，熱效率的模擬值比實驗值大了4%。

圖5 輸出功率和熱效率隨加熱溫度的變化曲線Fig.5 Output power and thermal efficiency vs.heating temperature

輸出功率和熱效率隨冷卻溫度的變化如圖6所示。此時加熱溫度為873 K，充氣壓力為7 MPa，工質(zhì)為He，工作頻率為30.2 Hz。從圖6中可以看出，輸出功率、熱效率的模擬值和實驗值也呈相同的變化趨勢。隨著冷卻溫度的升高，輸出功率和熱效率均減?。惠敵龉β实哪M值和實驗值誤差15%左右，熱效率的模擬值比實驗值大了4%。主要原因是冷熱源溫差對自由活塞斯特林發(fā)動機的輸出功率和熱效率有重要影響，對于Re-1000機器而言，冷熱源溫差越大，機器的輸出功率和熱效率越大。

圖6 輸出功率和熱效率隨冷卻溫度的變化曲線Fig.6 Output power and thermal efficiency vs.heat dissipation temperature

4 結論

根據(jù)Re-1000結構參數(shù)和運行參數(shù)，采用三級分析軟件Sage，建立了Re-1000自由活塞斯特林發(fā)動機的一維模型，并將Sage模擬結果與實驗數(shù)據(jù)進行對比，模擬結果顯示：輸出功率、熱效率的模擬值和實驗值呈相同的變化趨勢。隨著加熱溫度的升高，輸出功率和熱效率均增大；隨著冷卻溫度的升高，輸出功率和熱效率均減小。輸出功率模擬值與實驗值誤差在15%左右，熱效率模擬值比實驗值大4%。