自由活塞斯特林發動機Re-1000的模擬研究

陳 曦，崔 浩

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

0 引言

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

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

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

1 自由活塞斯特林發動機及參數

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

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

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

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

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

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

2.1 氣體域求解

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

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

圖2 氣體區域模型示意圖Fig.2 Gas area model

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

連續性方程：

動量方程：

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

（1）絲網氣體區域（Matrix gas）

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

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

（2）氣缸氣體區域（Cylinder-space gas）

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

（3）管道氣體區域（Duct gas）

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

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

2.2 一維Sage模型

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

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

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

3 模擬結果及驗證

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

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

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

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

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

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

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

4 結論

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