β型斯特林發動機有限體積建模及熱力學分析研究

邱文會，閆春杰，王小軍

（蘭州空間技術物理研究所 真空技術與物理重點實驗室，蘭州 730000）

關鍵字：斯特林發動機；動網格；有限體積；渦流

0 引言

在能源短缺形勢日益嚴峻的背景下，斯特林發動機因工作效率高、運轉特性好、結構簡單、易損部件少等優點受到國內外的廣泛重視[1]。與傳統的斯特林發動機相比，β型斯特林發動機具有質輕、無潤滑油、工作可靠等優勢，因此，在航空航天、太陽能發電、工業余熱回收利用、高海拔熱電轉化中有著廣闊的應用前景[2]。

當前研究斯特林熱力學的主要方法有實驗法、理論法和數值仿真方法。實驗研究雖然是最基本也是最可靠的研究方法，但是研究周期長，費用高。理論法較為局限，只適合前期設計階段的初步分析，結果并不可靠。數值仿真方法，即高維CFD方法，可以得到完整的物理圖像，可以對關注的結構或區域進行定量分析，尤其是可以對試驗難測量區域的數據進行采集和分析，是未來研究斯特林發動機的有力工具。

張存泉等[3]對回熱器建立了相關的理論模型，在對控制方程的無量綱化簡化和處理的基礎上，推導出不可逆的回熱器熵流方程。賈明興等[4]利用流體仿真軟件建立了三維模型，分析了孔隙率的變化對回熱器流動特性的影響，模擬結果發現，在高雷諾數下，隨著孔隙率的增大壓降也會增大，雷諾數較小時，孔隙率對壓降影響較小。Mahkamov[5]使用CFD軟件建立了一個生物質能斯特林機的三維仿真模型，模擬結果表明，斯特林發動機中連接腔室之間的“死容積”對斯特林的功率有著非常重要的影響，同時，在理論分析的基礎上對斯特林機設計和結構提出了改進意見。Singh等[6]用商業流體軟件對在研的斯特林機建立了二維模型，分析了在低雷諾數條件下，努賽爾數隨雷諾數的變化情況，同時也建立了三維模型，模擬了不同絲網目數下，雷諾數對回熱器壓降的影響。Kraitong等[7]針對低溫差太陽能碟式斯特林機建立模型，使用神經網絡算法對二階分析法進行改進，利用MATLAB軟件計算求解，并使用流體軟件對斯特林整機工作過程進行模擬仿真，最后將數值計算結果與三維仿真進行對比，結果顯示二者誤差很小。

基于Fluent流體仿真軟件平臺，在理論研究的基礎上，對在研的斯特林發動機進行合理降維轉化，建立了二維仿真模型，其中動力活塞和配氣活塞的振動頻率和振幅由動網格和UDF（用戶自定函數）控制，重點關注溫度場、壓力場和流場的分布，分析熱功參數的變化對工質流動和熱功轉化的影響。

1 有限體積建模

1.1 數學模型

回熱器采用多孔介質模型，其他腔室包括加熱器、冷卻器、壓縮腔和膨脹腔采用湍流模型。仿真實驗需要對模型進行簡化，才能達到CFD軟件運行的要求，簡化和假設為：（1）熱機氣缸內部工質為理想可壓縮氣體；（2）熱機的氣缸外殼都視為絕熱壁面，不考慮與外界的對流換熱和輻射換熱；（3）只考慮工質和兩個換熱器之間的對流換熱，忽略二者之間的輻射換熱；（4）忽略配氣活塞和動力活塞與缸體之間的摩擦以及工質泄漏損失；（5）加熱器和冷卻器面積由其相對應換熱翅片與相應腔室體積比來確定。

在Fluent軟件中有不同的湍流模型可以選擇，不同的模型適用于不同的流動情況。在斯特林熱機仿真領域使用廣泛的湍流模型有k-ε模型和k-ω模型。k-ε模型是半經驗公式，假設流動為完全湍流，可以忽略分子的黏性影響。標準k-ε模型的雷諾平均Navier-Stokes方程如式（1）和式（2）。

標準k-ε模型的湍動能k和湍流耗散ε的輸運方程如式（3）和式（4）：

式中：Gk為層流速度梯度產生的湍動能項；Gb為浮力產生的湍動能項；YM為湍流膨脹到全局流程中對耗散項的貢獻值，C為常量；σk和σε為k方程和ε方程的湍流Prandtl數；Sk和Sε為用戶定義的湍動能項和湍流耗散源項。

多孔介質模型的控制方程式是在上述公式基礎上進行改進得到的，在質量方程的流體密度項上都乘上一個系數φ，即孔隙率。將動量損失項St添加到動量方程中，St運算如式（5）：

回熱器中的能量方程則應用局部熱平衡方程[8]：

式中：K為Darcy滲透率；Cf是Forchheimer慣性阻力系數；C和β-1分別是多孔模型中設置的慣性阻力系數和黏性阻力系數。

1.2 幾何模型和參數設置

1.2.1 幾何模型和邊界處理

對于有動網格且運行為非穩態的計算模型，如果直接建立三維模型，不僅困難，而且計算周期長，所以取其半剖面應用Design Model模塊建立二維模型。首先根據長度位置參數進行草圖的繪制，接著生成草圖中的各個面域，其中加熱器和冷卻器模型采用面域的布爾運算得到，最后是邊界的命名。進入Mesh模塊，對不同位置的面域之間使用界面（Interface）來連接，尤其是換熱器與流域之間的界面需要設置熱耦合（Coupled）屬性，才能實現換熱器與工質之間的換熱，主要參數如圖1所示。在加熱器和冷卻器中填充矩形長條代表換熱器，且矩形面積與空白面積的比值與實際結構中換熱器與空體積之比相等，從而實現了三維實體向二維模型的轉化。具體計算條件及邊界運動的設置如表1所列。

仿真模型邊界參數需要根據實際情況來設置。如圖1所示，膨脹腔和壓縮腔外邊界選擇絕熱壁面，加熱器和冷卻器選擇固定溫度外壁和熱耦合屬性，即換熱器壁面參與對流換熱，回熱器采用多孔介質模型，其相對獨立的模塊不需要單獨進行設置，實際上膨脹腔溫度較高時會通過對流和輻射散熱，如果考慮這部分建模，會大幅增加模型的復雜程度和計算時間，而且求解方法也有適用的場合，如果采用輻射求解模型就會嚴重影響仿真結果的準確性和解的收斂性，實際上輻射換熱量相比工質對流換熱的熱量也是非常小的，所以在整個模型中不考慮發動機外壁面，內壁面間的輻射換熱。

圖1 斯特林機計算區域圖Fig.1 Computational domain of Stirling engine

表1 斯特林仿真初始值及邊界條件Tab.1 Initial value and boundary conditions of stirling simulation

1.2.2 網格劃分和求解器設置

網格由Design Model模塊進行劃分，其中回熱器部分劃分為結構化四邊形網格，因為配氣活塞和動力活塞動網格的設置，所以其他腔室劃分為三角形網格，有利于動網格的更新和重構。動網格更新的方法主要有彈性光滑（Spring Smoothing）、鋪層（Layering）和局部重構（Local Remeshing）。在文中使用鋪層和局部重構的組合方式。

對于物理模型，網格劃分越密，精度也就越高，但是也會顯著增加內存需求和計算量，所以在仿真實驗前進行網格無關性驗證非常重要。本文劃分了四組網格，監控點設置如圖1所示，選取膨脹監控點的溫度振幅和壓力振幅作為監控量，各組網格對應的物理量幅值列于表2中，可以看出，第二組相比第一組網格監控點物理量的相對誤差較大，第三組和第四組相比第二組相對誤差都在2%以內，所以在滿足計算精度的要求下，考慮計算成本，選擇第四組網格，網格數量為43 522個。經過前期仿真實驗，坍塌因子（Collapse Factor）與分裂因子（Split Factor）分別設置為0.3和0.5時網格重構效果很好，且有利于計算結果的收斂。為了確保模型更新動網格時不發生負體積，要求在一個時間步長（Time Step）緊貼壁面的網格寬度要大于網格移動的距離。因為緊貼壁面的網格高度為0.1 mm，所以在模擬計算中取時間步長為0.000 01 s。

表2 不同網格下監控點物理量幅值Tab.2 Physical quantities range of monitoring points in different grids

壓力耦合采用經典的Simple算法，相比Simpler算法，雖然占用內存較大，但計算速度快，魯棒性好。離散格式選擇二階迎風格式，時間的離散采用全隱式格式。收斂的判斷標準為能量方程的殘差小于10-6，其他物理量的殘差均小于10-4，判斷計算達到穩定的標準是圖1中各監控點監測的溫度、壓力、流速等數據在最后3個循環周期內變化幅度不超過1%，即可認定仿真達到穩定。

2 模擬結果與分析

通過數值模擬得到了斯特林發動機結構內溫度場、壓力損失和速度分布等情況，并分析了充氣壓力和加熱器溫度的變化對工質流動和做功的影響。

2.1 溫度分布

溫度分布云圖如圖2所示，可以看出，在一個周期內，膨脹腔與加熱器部分溫度始終較高，而壓縮腔和冷卻器的溫度都接近于仿真冷卻器的設置溫度。膨脹腔內部工質溫度沒有完全達到加熱器的溫度，即在一個周期內工質在膨脹腔中并不能完全吸熱就參與下一個周期的做功循環，尤其是靠近對稱軸附近的工質，由于渦旋的存在，無法靠近加熱器流體完成熱交換，導致溫升不夠，從而降低整體熱效率。

圖2 溫度分布云圖Fig.2 Temperature distribution

回熱器軸向溫度在一個周期內近似為線性分布，并且從低溫端到高溫端溫度逐漸升高，但任何一個徑向截面溫度都是隨著工質流動而變化的。由于回熱器內部密度分布不均，壓力非單調分布，導致熱端溫升快、溫降慢，冷端溫降快、溫升慢。

2.2 壓力分布

圖3是一個周期內斯特林發動機的壓力分布云圖，在3/4T時，即配氣活塞位于平衡位置向上止點運動，動力活塞位于下止點時，膨脹腔和壓縮腔壓差可以達到12 kPa，而在1T時即循環開始，壓差只有3 kPa左右。從圖中可以明顯看出，壓縮腔和膨脹腔之前的壓降主要來源于回熱器，在大部分時間內，軸向上回熱器都存在明顯的壓力梯度，由于加熱器和冷卻器流道的存在，所以換熱器兩端也存在一定的壓差。應用后處理工具可以得到相關壓力變化均值，圖4是3/4T回熱器軸向壓力分布圖，3/4T時換熱器、回熱器和冷卻器兩端的壓降值分別為12 463 Pa、432 Pa和876 Pa?？梢钥闯?，與回熱器兩端壓降相比，加熱器和冷卻器的兩端壓降很小。

圖3 一個周期內斯特林發動機的壓力分布云圖Fig.3 Pressure distribution of a cycle stirling engine

圖4 3/4T回熱器軸向壓力分布圖Fig.4 Axial pressure distribution of 3/4T regenerator

2.3 速度分布

速度分布云圖如圖5所示，從圖（a）和（c）中可以看出，流速較快的區域是加熱器和冷卻器與腔室連接處，在相應位置出現了渦流，導致工質之間熱質交換效率降低，對照溫度云圖可以發現，膨脹腔底部溫度較低的區域同時也是工質流速最慢的區域。在圖（b）和（d）中，膨脹腔和壓縮腔內部出現了2個及以上的渦流，這種流阻消耗能量的同時降低了工質做功能力，所以在進行斯特林發動機設計時，就應充分考慮渦流可能出現的情況，盡量使用圓角代替直角設計，換熱器和腔室連接區域增加導流葉片設計，調節活塞振幅和頻率避免狹長和太寬流域的存在。

2.4 充氣壓力對質量流率和做功的影響

在不改變其他參數的條件下，研究了充氣壓力（8×105Pa、12×105Pa、16×105Pa、20×105Pa）對回熱器上下邊界質量流率和膨脹腔與壓縮腔PV功及輸出功的影響。如圖6所示，隨著充氣壓力的增大，充入的工質質量也增加，回熱器兩端質量流率幅值有明顯的增大，且流率增幅與壓力增幅呈近似線性關系，從圖中也可以看出，回熱器熱端的質量流率與回熱器冷端的質量流率相比較小，主要是因為熱端溫度很高，工質密度小，所以質量流量相對較小。

圖5 速度分布云圖Fig.5 Velocity distribution

圖6 不同充氣壓力下回熱器兩端質量流率曲線Fig.6 Mass flow curve of regenerator under different charged pressure

圖7是壓縮腔、膨脹腔和輸出凈功隨充氣壓力的變化曲線，可以發現，隨著充氣壓力的增加，膨脹腔和壓縮腔PV功都在增加，因為膨脹腔內工質平均壓力高于壓縮腔工質平均壓力，在活塞行程和頻率不變的情況下，一個循環內PV功更多，所以凈功也

隨之增加。

圖7 不同充氣壓力下各腔室PV功曲線Fig.7 PV power curves of each chamber under different charged pressure

2.5 加熱器溫度對質量流率和做功的影響

對整機性能影響最大的是加熱器溫度。當其他參數保持不變時，改變加熱器的溫度，研究了回熱器上下邊界質量流率，膨脹腔與壓縮腔PV功和凈輸出功的變化情況。

由圖8可以看出，加熱器溫度越高，通過加熱器上邊界的質量流率越低，主要是因為加熱器溫度變高，其附近被加熱的工質的比容越大，流過加熱器上邊界的質量流率就越低。在整個腔室內所充工質質量不變的情況下，冷卻器附近工質的密度相對較大，流過冷卻器下邊界的質量流量越大。如圖9所示，隨著加熱器熱端溫度的升高，膨脹腔內工質壓力也隨之增大，PV功增加，壓縮腔的PV功則會減少，輸出功增加，由此可以看出，加熱器溫度越高，斯特林發動機膨脹腔PV功和凈輸出功都會明顯增大。

圖8 不同加熱器溫度下回熱器兩端質量流率曲線Fig.8 Mass flow curve of regenerator under different heater temperature

圖9 不同加熱器溫度下各腔室PV功曲線Fig.9 PV power curve of each chamber under different heater temperature

3 結論

利用FLUENT軟件對β型自由活塞式斯特林發動機進行二維仿真，通過數學模型的建立和求解，在仿真實驗穩定后，得到了一個周期內不同時間節點，熱機內部工質溫度場，壓力場和速度場的變化情況，通過監控點的設置，重點關注了回熱器上下邊界質量流率，各個腔室PV功和凈輸出功等變化，可以得到結論：

（1）加熱器和換熱器溫度、壓力等分布相對均勻，回熱器端部流場變化劇烈，腔室底部流體熱質交換不足。回熱器中工質存在明顯的溫度、壓力梯度，且沿軸向近似為線性分布，與單向流相比，回熱器內工質的流動存在復雜的非單調非對稱的振蕩流。

（2）壓縮腔、膨脹腔中溫度和壓力分布較為均勻，但是壓力在回熱器兩端及腔室局部區域分布不均，對照流場分布云圖可知，在這些區域存在明顯的渦流，這些渦流會增大流阻耗散能量，影響熱機整體熱效率。

（3）充氣壓力的增大導致腔室內質量流量增大，膨脹腔和壓縮腔工質壓力也隨之增大，膨脹腔PV功增加的更多，凈功變大；隨著加熱器溫度的提高，換熱器與工質換熱更加劇烈，回熱器壓降降低，膨脹腔工質PV功增大，壓縮腔PV功降低，所以增加加熱器溫度可以明顯增加斯特林發動機的熱功轉換效率，但同時也要考慮加熱器材料的耐熱性能和缸體材料的耐腐蝕性能。