自由活塞斯特林發(fā)動(dòng)機(jī)Re-1000的模擬研究

陳 曦，崔 浩

（上海理工大學(xué)能源與動(dòng)力工程學(xué)院，上海 200093）

0 引言

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

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

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

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

自由活塞斯特林發(fā)動(dòng)機(jī)（FPSE），主要包括：配氣活塞、動(dòng)力活塞、加熱器、回?zé)崞骱屠鋮s器等，物理模型如圖1所示。

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

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

表1 Re-1000自由活塞斯特林發(fā)動(dòng)機(jī)的結(jié)構(gòu)參數(shù)Table1 Structural parameters of the Re-1000 free-piston Stirling engine

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

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

2.1 氣體域求解

Sage軟件會基于所建立的發(fā)動(dòng)機(jī)模型，在空間節(jié)點(diǎn)上采用中心差分法，在時(shí)間節(jié)點(diǎn)上采用后三點(diǎn)差分法，把發(fā)動(dòng)機(jī)的控制方程離散成一個(gè)大型的非線性方程組，然后采用非線性解算器對這個(gè)方程組進(jìn)行迭代求解，得到發(fā)動(dòng)機(jī)模型的輸出參數(shù)，其求解過程非?？欤蠼馑玫臅r(shí)間一般都在幾分鐘之內(nèi)。

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

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

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

連續(xù)性方程：

動(dòng)量方程：

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

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

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

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

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

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

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

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

為了迭代求解方程，還需要得出氣體區(qū)域流動(dòng)幾何的摩擦因子f、努賽爾數(shù)Nu、軸向?qū)崧蔔k的經(jīng)驗(yàn)公式。

2.2 一維Sage模型

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

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

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

3 模擬結(jié)果及驗(yàn)證

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

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

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

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

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

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

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

4 結(jié)論

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