自由活塞斯特林發(fā)動機(jī)的響應(yīng)面法設(shè)計

冶文蓮 ，王小軍，劉迎文 ，閆春杰

（1.蘭州空間技術(shù)物理研究所 真空技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室，蘭州 730000；2.西安交通大學(xué)能源與動力工程學(xué)院 熱流科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室，西安 710049）

0 引言

自由活塞斯特林發(fā)動機(jī)（Free Piston Stirling Engine，F(xiàn)PSE）作為一種外燃式發(fā)動機(jī)，具有效率高、壽命長、結(jié)構(gòu)簡單、噪音低、不易磨損、自啟動等特點(diǎn)?？衫脽嵩捶秶鷱V，例如太陽能、同位素能源、生物能、余熱等[1]。其技術(shù)特征為：柔性彈簧支撐、無磨損間隙密封、動力活塞與配氣活塞通過流體傳動完成機(jī)械能的輸出。在所有發(fā)動機(jī)中，F(xiàn)PSE的循環(huán)效率最接近卡諾循環(huán)，可以達(dá)到理論卡諾循環(huán)效率的60%[2]。

近幾十年，研究人員在斯特林發(fā)動機(jī)的熱力學(xué)分析方面開展了充分研究。根據(jù)Martini[3]的命名法，斯特林發(fā)動機(jī)熱力學(xué)分析法大體可分為零級到四級五大類，包括用于斯特林裝置初步設(shè)計優(yōu)化的零級分析法，以斯密特等溫分析法[4]為代表的一級分析法，以絕熱分析為代表的二級分析法[5]，由Finkelestin提出的三級分析法，基于CFD軟件的四級分析法[6-7]。其中，三級分析法采用節(jié)點(diǎn)法求解每個節(jié)點(diǎn)處氣體工質(zhì)的質(zhì)量、動量以及能量守恒方程，各部件采用模塊化連接，可建立自由活塞斯特林發(fā)動機(jī)整機(jī)內(nèi)部氣體工質(zhì)瞬時流量、壓力等模型，能有效預(yù)測整機(jī)的輸出功率和熱效率等性能參數(shù)。Gedeon[11]編制了專門的軟件Sage，國內(nèi)外研究學(xué)者采用此軟件開展了斯特林發(fā)動機(jī)的模擬計算[8-9]，在預(yù)測斯特林發(fā)動機(jī)性能方面準(zhǔn)確性較高，已被NASA GRC廣泛采用[10-12]。四級分析法雖然在整機(jī)設(shè)計、性能優(yōu)化等方面具有很好的應(yīng)用，亦可直觀地表示整機(jī)內(nèi)部溫度場、壓力場、速度場等分布，但是在整機(jī)設(shè)計優(yōu)化方面比較耗時。因此，本文將采用Sage軟件開展FPSE的設(shè)計計算。

斯特林發(fā)動機(jī)的性能好壞主要取決于眾多幾何參數(shù)的選取，包括活塞桿直徑、活塞直徑、加熱器、回?zé)崞饕约袄鋮s器結(jié)構(gòu)參數(shù)等，評價斯特林發(fā)動機(jī)熱力學(xué)性能的指標(biāo)有輸出功率、熱效率，文中將其作為輸出參數(shù)進(jìn)行優(yōu)化設(shè)計。從國內(nèi)外相關(guān)文獻(xiàn)中可知，各個參數(shù)間具有復(fù)雜的相關(guān)性，而且前期優(yōu)化設(shè)計中更多的考慮了單因素的影響[13]。為了克服上述問題，文中引入一種新的方法-響應(yīng)面法（RSM）[14-16]，可有效預(yù)測斯特林發(fā)動機(jī)性能。由于影響發(fā)動機(jī)性能的結(jié)構(gòu)參數(shù)較多（20個及以上），為便于后期整機(jī)調(diào)試，選取活塞桿直徑、加熱器、回?zé)崞饕约袄鋮s器長度作為設(shè)計參數(shù)，輸出功率和熱效率作為響應(yīng)值，揭示發(fā)動機(jī)輸出功率和效率隨各項(xiàng)參數(shù)的變化規(guī)律。

1 數(shù)學(xué)模型

1.1 模型介紹

本文研究的自由活塞斯特林發(fā)動機(jī)結(jié)構(gòu)示意圖如圖1所示，采用β型結(jié)構(gòu)。主要由冷卻器、回?zé)崞?、加熱器、兩個活塞、直線電機(jī)等組件組成。配氣活塞桿穿過動力活塞，兩活塞通過氣體壓力的交變流動自由運(yùn)動，均由板彈簧支撐。由于配氣活塞端溫度較高，選用材料為In718，動力活塞材料為硬質(zhì)鋁合金。FPSE基本參數(shù)如表1所列。

表1 自由活塞斯特林發(fā)動機(jī)基本參數(shù)Tab.1 The technical specifications of the FPSE

圖1 自由活塞斯特林發(fā)動機(jī)結(jié)構(gòu)示意圖Fig.1 The schematic illustration of the FPSE

Sage軟件[17]是一款用于模擬、優(yōu)化斯特林循環(huán)發(fā)動機(jī)、斯特林制冷機(jī)的一維專用軟件，已在回?zé)崾桨l(fā)動機(jī)的研究中被采用。其采用可視化界面，對發(fā)動機(jī)的各個部件一一建模，各模塊之間通過質(zhì)量、能量流動進(jìn)行連接。根據(jù)圖1的斯特林發(fā)動機(jī)結(jié)構(gòu)，采用Sage對其建模，如圖2所示。模型中包括配氣活塞，動力活塞，環(huán)形翅片式加熱器和冷卻器，不銹鋼隨機(jī)纖維回?zé)崞鳌⑴蛎浨?、壓縮腔和背壓腔等模塊。此外，在建模過程中考慮了配氣活塞處穿梭損失，回?zé)崞鞅诿鎸?dǎo)熱損失、配氣活塞殼體導(dǎo)熱損失等，以及膨脹和壓縮腔部分的死容積。各個部件間根據(jù)實(shí)際狀態(tài)進(jìn)行質(zhì)量流、能量流的連接。模型中假設(shè)氦氣為理想氣體，忽略輻射漏熱以及與外界環(huán)境的對流換熱。壓縮腔和膨脹腔外壁設(shè)置為絕熱邊界條件。

圖2 自由活塞斯特林發(fā)動機(jī)Sage模型圖Fig.2 Sage model of FPSE

1.2 氣體動力方程

自由活塞斯特林發(fā)動機(jī)內(nèi)部氣體的流動和傳熱為非定常周期型可壓縮流動。為便于說明一維流動和傳熱，將空間平均和時間平均的內(nèi)部流動氣體各方程表述為[17]：

連續(xù)方程：

動量方程：

能量方程：

式中：t為時間，s；ρ為工質(zhì)氣體密度，kg/m3；A為氣體域截面積，即x方向上的截面積，m2；u為流體速度，m/s；p為壓力，Pa；e為質(zhì)量能；q為氣體軸向?qū)釤崃?，W/m2；Qw為對流換熱熱流，W；動量方程中的F是斯托克斯的應(yīng)力項(xiàng)，在一維方程中無法求解?？梢杂脫Q熱器的水力直徑dh和長度L的Darcy阻力系數(shù)f和總局部損失系數(shù)K來表示：

其中阻力系數(shù)根據(jù)不同換熱器結(jié)構(gòu)來選取，文中采用的回?zé)崞鳛椴讳P鋼隨機(jī)纖維結(jié)構(gòu)，加熱器和冷卻器為翅片結(jié)構(gòu)。其中回?zé)崞髯枇ο禂?shù)表示為：

式中：a1=25.7α+79.8；a2=0.146α+3.76；a3=-0.002 83α-0.074 08；α=

加熱器和冷卻器阻力系數(shù)：

式中：φ為回?zé)崞骺紫堵?；Re為雷諾數(shù)；ε為表面不規(guī)則物體的平均高度。

在絲網(wǎng)氣體區(qū)域中，能量方程中的Qw由通過截面單位長度的熱流表示：

式中：Nu為努賽爾數(shù)；kg為氣體的導(dǎo)熱系數(shù)，W/（m·K）；Ax為單位長度的流體和表面接觸面積，m2；Tw-T為計算流體通過表面和平均截面間的溫差，K。不同類型的氣體區(qū)域由于流動換熱特性不同，Qw有不同的表達(dá)，應(yīng)根據(jù)實(shí)際情況進(jìn)行計算。

對上述氣體動力方程式（1）～（3）在均勻分布的網(wǎng)格點(diǎn)上進(jìn)行離散求解。在這些網(wǎng)格點(diǎn)上，變量ρ、ρuA及ρe通過隱式函數(shù)進(jìn)行插值，然后在一系列網(wǎng)格上進(jìn)行求解，再求解密度ρ以滿足整個系統(tǒng)的動量和能量守恒。

2 響應(yīng)面優(yōu)化方法

響應(yīng)面優(yōu)化法（Response Surface Methodology，RSM），也叫響應(yīng)面曲面法，主要通過對指定設(shè)計空間內(nèi)樣本點(diǎn)的集合進(jìn)行有限的試驗(yàn)設(shè)計，采用多元線性回歸的方法，將試驗(yàn)中的各個因素及其交叉因素進(jìn)行多項(xiàng)式擬合，再根據(jù)函數(shù)的響應(yīng)面等值線和回歸方程等確定試驗(yàn)的較優(yōu)組合。在工程優(yōu)化設(shè)計中，應(yīng)用響應(yīng)面法可以得到響應(yīng)值與輸入變量之間的變化關(guān)系。RSM主要包括Box-Behnken（BBD）設(shè)計、均勻外殼設(shè)計（Uniform Shell Design，USD）和中心組合設(shè)計（Central Composite Design，CCD）三種，BBD是RSM方法中常用的試驗(yàn)設(shè)計方法，適用于2～5個因素的試驗(yàn)，文中采用BBD優(yōu)化法設(shè)計試驗(yàn)點(diǎn)。采用的RSM方法主要用于評價部分關(guān)鍵結(jié)構(gòu)參數(shù)對斯特林發(fā)動機(jī)輸出功率和熱效率的影響。RSM的實(shí)施過程如下：

首先，根據(jù)優(yōu)化需求定義因素的水平、范圍和響應(yīng)值。其次，基于BBD優(yōu)化方法進(jìn)行相應(yīng)的RSM設(shè)計得出試驗(yàn)組合，將這些組合值代入Sage模型中進(jìn)行計算得到響應(yīng)值。RSM法根據(jù)試驗(yàn)安排及計算結(jié)果得出二次回歸擬合方程，同時對擬合結(jié)果進(jìn)行方差分析（ANOVA）。根據(jù)方差分析結(jié)果判斷出回歸方程擬合度是否良好，若擬合度良好則認(rèn)為模型可靠性高，從而進(jìn)行結(jié)果分析。若回歸方程擬合度較差，則需要進(jìn)一步添加試驗(yàn)點(diǎn)，對其重新設(shè)計。采用RSM建立結(jié)構(gòu)參數(shù)與自由活塞斯特林發(fā)動機(jī)輸出功率和熱效率之間的關(guān)聯(lián)式。其中，二次回歸模型描述如式（6）：

式中：y為FPSE系統(tǒng)的響應(yīng)值；xi和xj為隨機(jī)變量；β0、βi、βii和（i=0，1，2，…，k；j=0，1，2，…，k）為待定系數(shù)，由樣本點(diǎn)迭代確定；ε為統(tǒng)計誤差。

影響斯特林發(fā)動機(jī)熱力學(xué)性能的參數(shù)比較多，包括三個換熱器結(jié)構(gòu)參數(shù)、冷熱端溫度和充氣壓力等。國內(nèi)外研究者對這方面也做過相關(guān)研究[18，19]。在本次RSM試驗(yàn)中，將配氣活塞桿直徑、加熱器、回?zé)崞饕约袄鋮s器的長度作為因變量，主要是便于后期整機(jī)熱力學(xué)性能的調(diào)試，若要改動其他參數(shù)，勢必給整機(jī)結(jié)構(gòu)的調(diào)整帶來諸多不便?；?zé)崞髟趯?shí)際設(shè)計計算過程中，配氣活塞長度也隨換熱器長度發(fā)生變化。各參數(shù)的范圍和水平主要基于實(shí)際工程設(shè)計，其分別由A、B、C和D表示，輸出功率和熱效率為響應(yīng)值，分別用Wnet和Eff表示。其中，活塞桿直徑范圍為0.004～0.008 m，加熱器長度范圍0.010～0.025 m，回?zé)崞鏖L度范圍0.030～0.045 m，冷卻器長度范圍0.015～0.030 m，如表2所列。

表2 各因素水平和范圍Tab.2 Ranges and levels in BBD

3 結(jié)果分析

根據(jù)BBD試驗(yàn)設(shè)計安排試驗(yàn)點(diǎn)、各因素和水平，隨機(jī)得出試驗(yàn)總次數(shù)為29次。將不同結(jié)構(gòu)參數(shù)組合的因變量輸入至Sage模型中，計算得出不同配氣活塞桿直徑、加熱器、冷卻器以及回?zé)崞鏖L度下斯特林發(fā)動機(jī)的輸出功率和熱效率值，結(jié)果如表3所列。表中第1列表示試驗(yàn)次數(shù)，后4列表示試驗(yàn)條件，最后2列為試驗(yàn)結(jié)果。Dr表示活塞桿直徑，Lh、Lr和Lc分別表示加熱器、回?zé)崞骱屠鋮s器長度。對表3中數(shù)據(jù)進(jìn)行回歸分析，得到四因素對輸出功率和熱效率的二次多項(xiàng)回歸擬合方程：

表3 試驗(yàn)安排矩陣Tab.3 Design of experimental matrix

3.1ANOVA分析

基于BBD算法，式（7）和式（8）給出了輸出功率和熱效率的二次模型回歸方程。在響應(yīng)面優(yōu)化分析中二次回歸方程的準(zhǔn)確性和顯著性由方差分析檢驗(yàn)。表4、表5分別給出輸出功率和熱效率響應(yīng)面優(yōu)化模型的方差分析結(jié)果。平方和反映所優(yōu)化試驗(yàn)數(shù)據(jù)的總體波動性，方差由每部分的平方和除以自由度得到，F(xiàn)值反映模型對響應(yīng)值的影響，P值用以判斷因素影響相應(yīng)的顯著性，P值小于0.05表示因素對響應(yīng)影響很大。R2是判斷回歸方程有效性的一個重要參數(shù)，其值越接近于1，表明模型準(zhǔn)確度較高。修正擬合系數(shù)R2adj反映模型能夠解釋的響應(yīng)值變化百分比，預(yù)測擬合系數(shù)R2Pred與修正擬合系數(shù)R2adj越接近于1，表明模型越顯著。信噪比用來衡量有效數(shù)據(jù)和干擾之間的比例關(guān)系，通常認(rèn)為其值大于4，說明有效數(shù)據(jù)充分，模型可靠，能夠準(zhǔn)確預(yù)測響應(yīng)值。由表4可知，R2為0.999 9，R2adj為0.999 8，說明模型的預(yù)測值和實(shí)際值吻合度較好。還可得知，活塞桿直徑、加熱器、回?zé)崞骱屠鋮s器長度及AC、AD、BC、BD、CD、A2、B2、C2和D2的P值都小于0.05，說明參數(shù)對輸出功率有顯著影響，而其他參數(shù)如AB的P值大于0.05，說明其影響較小。由表5可得，R2為0.998 5，R2adj為0.997 1，R2Pred為0.999 4說明吻合度較好。除了AB和AC的P值大于0.05外，其他參數(shù)對熱效率都有顯著的作用。由表4和表5得知，BC、BD及CD兩兩參數(shù)對輸出功率和熱效率都有顯著影響，因此作為重點(diǎn)因素進(jìn)行分析。

表4 輸出功率的方差分析Tab.4ANOVAfor output power

表5 熱效率的方差分析Tab.5 ANOVAfor thermal efficiency

圖3為響應(yīng)面（RSM）預(yù)測的輸出功率、熱效率和Sage模型計算值之間的關(guān)系圖。由圖可知，大部分設(shè)計點(diǎn)與圖中斜線十分接近，說明模型的模擬結(jié)果能夠有效預(yù)測斯特林發(fā)動機(jī)的熱力學(xué)性能。

圖3 響應(yīng)面預(yù)測結(jié)果與Sage計算結(jié)果圖Fig.3 RSM predictions versus Sage model's results

3.2 各因素之間耦合影響分析

為研究四個參數(shù)對自由活塞斯特林發(fā)動機(jī)性能的相互影響，RSM分析法提供多因素的二維等高線圖和三維表面分析，顯示了兩個因素之間的相互作用對輸出功率和熱效率的影響。圖4～9分別給出了BC、BD和CD（對應(yīng)于表4和表5的重要因素）對輸出功率和熱效率的相互影響。

圖4和圖5為加熱器和回?zé)崞鏖L度的綜合效應(yīng)對FPSE輸出功率和熱效率的影響，此時活塞桿直徑為0.006 m，冷卻器長度為0.022 5 m。由圖4（a）可知，當(dāng)加熱器長度增加至0.025 m，回?zé)崞鏖L度最短時的輸出功率最大。由圖4（b）得知，加熱器長度和回?zé)崞鏖L度間的相互作用對輸出功率影響較大，即兩者長度增加輸出功率均線性降低。從圖5（a）和（b）所示的等高線和三維表面圖中可看出，回?zé)崞鏖L度增加至最大時，熱效率隨加熱器長度增加而非線性增加。這是因?yàn)榧訜崞骱突責(zé)崞鏖L度的增加會導(dǎo)致流阻損失的增加，使得輸出功率隨之降低，但在模型中引起的回?zé)崞魍鈿け诿鎸?dǎo)熱損失降低。因此，針對實(shí)際機(jī)型，需選擇合適的回?zé)崞骱图訜崞鏖L度。

圖4 加熱器與回?zé)崞鏖L度對輸出功率的耦合影響圖Fig.4 The combined effect of length of heater and regenerator on output power

圖6和圖7顯示加熱器和冷卻器長度的綜合效應(yīng)對輸出功率和熱效率的影響，此時活塞桿直徑為0.006 m，回?zé)崞鏖L度為0.037 5 m。由圖6（a）中可清晰地看出，左下角區(qū)域?yàn)檩敵龉β实淖畲笾?，該值大?00 W，此時的加熱器長度為0.015 m，冷卻器長度為0.01 m。由圖6（b）可知，加熱器長度由0.025 m降低至0.01 m時，輸出功率隨冷卻器長度的增加而減小，并呈現(xiàn)線性趨勢。從圖7（a）和（b）所示的等高線和三維表面圖中可以看出，加熱器長度增加至最大時，熱效率隨冷卻器長度先增加，達(dá)到最大值后降低。這是因?yàn)楫?dāng)冷卻器長度比較短時，壓降損失影響不明顯，輸出功率降低亦不明顯，熱效率增加。當(dāng)換熱器長度增加到一定值時，由流動阻力引起的壓降損失對發(fā)動機(jī)性能的影響變大。因此，發(fā)動機(jī)性能隨冷卻器長度的增加有最佳值。

圖5 加熱器與回?zé)崞鏖L度對熱效率的耦合影響圖Fig.5 The combined effect of length of heater and regenerator on thermal efficiency

圖6 加熱器與冷卻器長度對輸出功率的耦合影響圖Fig.6 The combined effect of length of heater and cooler on output power

圖7 加熱器與冷卻器長度對熱效率的耦合影響圖Fig.7 The combined effect of length of heater and cooler on thermal efficiency

圖8和圖9顯示了回?zé)崞骱屠鋮s器長度的綜合效應(yīng)對輸出功率和熱效率的影響，此時活塞桿直徑為0.006 m，加熱器長度為0.015 m。由圖8（a）中可知，等高線圖中左下角區(qū)域的輸出功率值最大，該值大于110 W，此時冷卻器的長度為0.015 m，回?zé)崞鏖L度為0.03 m。圖8（b）三維曲面圖中，回?zé)崞鏖L度增加時，輸出功率隨冷卻器長度的增加而降低，而且回?zé)崞鏖L度的增加使得輸出功率的降低幅度更大。由圖9（a）和（b）所示的等高線和三維表面圖得知，回?zé)崞鏖L度增加至0.045 m時，熱效率隨冷卻器長度呈二次曲線形式增加。因此，綜合考慮發(fā)動機(jī)性能，應(yīng)該根據(jù)實(shí)際結(jié)構(gòu)設(shè)計以及后續(xù)加工工藝等因素選擇合適長度的回?zé)崞骱屠鋮s器。

圖8 回?zé)崞髋c冷卻器長度對輸出功率的耦合影響圖Fig.8 The combined effect of length of regenerator and cooler on output power

圖9 回?zé)崞髋c冷卻器長度對熱功率的耦合影響圖Fig.9 The combined effect of length of regenerator and cooler on thermal efficiency

4 結(jié)論

響應(yīng)面是一種綜合試驗(yàn)設(shè)計和數(shù)學(xué)建模的優(yōu)化方法，可有效減少試驗(yàn)次數(shù)，給出比較直觀、簡單易懂的等高線圖和三維立體圖，并可有效考察影響因素之間的交互作用。通過采用RSM揭示了自由活塞斯特林發(fā)動機(jī)換熱器以及配氣活塞桿結(jié)構(gòu)參數(shù)對整機(jī)熱力學(xué)性能的耦合影響規(guī)律。得出結(jié)論如下：

（1）基于BBD方法安排了29次試驗(yàn)矩陣，建立輸出功率和熱效率的二次回歸模型。通過對回歸方程的方差分析和模型診斷，發(fā)現(xiàn)所選的四個結(jié)構(gòu)參數(shù)對斯特林發(fā)動機(jī)性能具有顯著影響。響應(yīng)面預(yù)測模型與Sage計算結(jié)果具有高度一致性，有效證明了模型的有效性和可靠性。

（2）由RSM的二維等高線圖和三維表面圖，探明了兩個因素之間的相互作用對輸出功率和熱效率的影響規(guī)律。結(jié)果顯示：輸出功率隨著回?zé)崞?、加熱器以及冷卻器長度的增加而降低，熱效率隨之增加而增加，冷卻器長度對熱效率有最佳值，較長的換熱器長度會對整機(jī)熱效率產(chǎn)生積極影響。在以后對FPSE的改進(jìn)設(shè)計中，需合理選擇其長度，以優(yōu)化發(fā)動機(jī)性能并使結(jié)構(gòu)更緊湊。本文對于有效設(shè)計斯特林發(fā)動機(jī)結(jié)構(gòu)參數(shù)以及預(yù)測其性能具有良好的指導(dǎo)意義。