邱文會(huì)，閆春杰，王小軍

（蘭州空間技術(shù)物理研究所 真空技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室，蘭州 730000）

關(guān)鍵字：斯特林發(fā)動(dòng)機(jī)；動(dòng)網(wǎng)格；有限體積；渦流

0 引言

在能源短缺形勢(shì)日益嚴(yán)峻的背景下，斯特林發(fā)動(dòng)機(jī)因工作效率高、運(yùn)轉(zhuǎn)特性好、結(jié)構(gòu)簡(jiǎn)單、易損部件少等優(yōu)點(diǎn)受到國(guó)內(nèi)外的廣泛重視[1]。與傳統(tǒng)的斯特林發(fā)動(dòng)機(jī)相比，β型斯特林發(fā)動(dòng)機(jī)具有質(zhì)輕、無(wú)潤(rùn)滑油、工作可靠等優(yōu)勢(shì)，因此，在航空航天、太陽(yáng)能發(fā)電、工業(yè)余熱回收利用、高海拔熱電轉(zhuǎn)化中有著廣闊的應(yīng)用前景[2]。

當(dāng)前研究斯特林熱力學(xué)的主要方法有實(shí)驗(yàn)法、理論法和數(shù)值仿真方法。實(shí)驗(yàn)研究雖然是最基本也是最可靠的研究方法，但是研究周期長(zhǎng)，費(fèi)用高。理論法較為局限，只適合前期設(shè)計(jì)階段的初步分析，結(jié)果并不可靠。數(shù)值仿真方法，即高維CFD方法，可以得到完整的物理圖像，可以對(duì)關(guān)注的結(jié)構(gòu)或區(qū)域進(jìn)行定量分析，尤其是可以對(duì)試驗(yàn)難測(cè)量區(qū)域的數(shù)據(jù)進(jìn)行采集和分析，是未來(lái)研究斯特林發(fā)動(dòng)機(jī)的有力工具。

張存泉等[3]對(duì)回?zé)崞鹘⒘讼嚓P(guān)的理論模型，在對(duì)控制方程的無(wú)量綱化簡(jiǎn)化和處理的基礎(chǔ)上，推導(dǎo)出不可逆的回?zé)崞黛亓鞣匠?。賈明興等[4]利用流體仿真軟件建立了三維模型，分析了孔隙率的變化對(duì)回?zé)崞髁鲃?dòng)特性的影響，模擬結(jié)果發(fā)現(xiàn)，在高雷諾數(shù)下，隨著孔隙率的增大壓降也會(huì)增大，雷諾數(shù)較小時(shí)，孔隙率對(duì)壓降影響較小。Mahkamov[5]使用CFD軟件建立了一個(gè)生物質(zhì)能斯特林機(jī)的三維仿真模型，模擬結(jié)果表明，斯特林發(fā)動(dòng)機(jī)中連接腔室之間的“死容積”對(duì)斯特林的功率有著非常重要的影響，同時(shí)，在理論分析的基礎(chǔ)上對(duì)斯特林機(jī)設(shè)計(jì)和結(jié)構(gòu)提出了改進(jìn)意見。Singh等[6]用商業(yè)流體軟件對(duì)在研的斯特林機(jī)建立了二維模型，分析了在低雷諾數(shù)條件下，努賽爾數(shù)隨雷諾數(shù)的變化情況，同時(shí)也建立了三維模型，模擬了不同絲網(wǎng)目數(shù)下，雷諾數(shù)對(duì)回?zé)崞鲏航档挠绊?。Kraitong等[7]針對(duì)低溫差太陽(yáng)能碟式斯特林機(jī)建立模型，使用神經(jīng)網(wǎng)絡(luò)算法對(duì)二階分析法進(jìn)行改進(jìn)，利用MATLAB軟件計(jì)算求解，并使用流體軟件對(duì)斯特林整機(jī)工作過(guò)程進(jìn)行模擬仿真，最后將數(shù)值計(jì)算結(jié)果與三維仿真進(jìn)行對(duì)比，結(jié)果顯示二者誤差很小。

基于Fluent流體仿真軟件平臺(tái)，在理論研究的基礎(chǔ)上，對(duì)在研的斯特林發(fā)動(dòng)機(jī)進(jìn)行合理降維轉(zhuǎn)化，建立了二維仿真模型，其中動(dòng)力活塞和配氣活塞的振動(dòng)頻率和振幅由動(dòng)網(wǎng)格和UDF（用戶自定函數(shù)）控制，重點(diǎn)關(guān)注溫度場(chǎng)、壓力場(chǎng)和流場(chǎng)的分布，分析熱功參數(shù)的變化對(duì)工質(zhì)流動(dòng)和熱功轉(zhuǎn)化的影響。

1 有限體積建模

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

回?zé)崞鞑捎枚嗫捉橘|(zhì)模型，其他腔室包括加熱器、冷卻器、壓縮腔和膨脹腔采用湍流模型。仿真實(shí)驗(yàn)需要對(duì)模型進(jìn)行簡(jiǎn)化，才能達(dá)到CFD軟件運(yùn)行的要求，簡(jiǎn)化和假設(shè)為：（1）熱機(jī)氣缸內(nèi)部工質(zhì)為理想可壓縮氣體；（2）熱機(jī)的氣缸外殼都視為絕熱壁面，不考慮與外界的對(duì)流換熱和輻射換熱；（3）只考慮工質(zhì)和兩個(gè)換熱器之間的對(duì)流換熱，忽略二者之間的輻射換熱；（4）忽略配氣活塞和動(dòng)力活塞與缸體之間的摩擦以及工質(zhì)泄漏損失；（5）加熱器和冷卻器面積由其相對(duì)應(yīng)換熱翅片與相應(yīng)腔室體積比來(lái)確定。

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

標(biāo)準(zhǔn)k-ε模型的湍動(dòng)能k和湍流耗散ε的輸運(yùn)方程如式（3）和式（4）：

式中：Gk為層流速度梯度產(chǎn)生的湍動(dòng)能項(xiàng)；Gb為浮力產(chǎn)生的湍動(dòng)能項(xiàng)；YM為湍流膨脹到全局流程中對(duì)耗散項(xiàng)的貢獻(xiàn)值，C為常量；σk和σε為k方程和ε方程的湍流Prandtl數(shù)；Sk和Sε為用戶定義的湍動(dòng)能項(xiàng)和湍流耗散源項(xiàng)。

多孔介質(zhì)模型的控制方程式是在上述公式基礎(chǔ)上進(jìn)行改進(jìn)得到的，在質(zhì)量方程的流體密度項(xiàng)上都乘上一個(gè)系數(shù)φ，即孔隙率。將動(dòng)量損失項(xiàng)St添加到動(dòng)量方程中，St運(yùn)算如式（5）：

回?zé)崞髦械哪芰糠匠虅t應(yīng)用局部熱平衡方程[8]：

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

1.2 幾何模型和參數(shù)設(shè)置

1.2.1 幾何模型和邊界處理

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

仿真模型邊界參數(shù)需要根據(jù)實(shí)際情況來(lái)設(shè)置。如圖1所示，膨脹腔和壓縮腔外邊界選擇絕熱壁面，加熱器和冷卻器選擇固定溫度外壁和熱耦合屬性，即換熱器壁面參與對(duì)流換熱，回?zé)崞鞑捎枚嗫捉橘|(zhì)模型，其相對(duì)獨(dú)立的模塊不需要單獨(dú)進(jìn)行設(shè)置，實(shí)際上膨脹腔溫度較高時(shí)會(huì)通過(guò)對(duì)流和輻射散熱，如果考慮這部分建模，會(huì)大幅增加模型的復(fù)雜程度和計(jì)算時(shí)間，而且求解方法也有適用的場(chǎng)合，如果采用輻射求解模型就會(huì)嚴(yán)重影響仿真結(jié)果的準(zhǔn)確性和解的收斂性，實(shí)際上輻射換熱量相比工質(zhì)對(duì)流換熱的熱量也是非常小的，所以在整個(gè)模型中不考慮發(fā)動(dòng)機(jī)外壁面，內(nèi)壁面間的輻射換熱。

圖1 斯特林機(jī)計(jì)算區(qū)域圖Fig.1 Computational domain of Stirling engine

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

1.2.2 網(wǎng)格劃分和求解器設(shè)置

網(wǎng)格由Design Model模塊進(jìn)行劃分，其中回?zé)崞鞑糠謩澐譃榻Y(jié)構(gòu)化四邊形網(wǎng)格，因?yàn)榕錃饣钊蛣?dòng)力活塞動(dòng)網(wǎng)格的設(shè)置，所以其他腔室劃分為三角形網(wǎng)格，有利于動(dòng)網(wǎng)格的更新和重構(gòu)。動(dòng)網(wǎng)格更新的方法主要有彈性光滑（Spring Smoothing）、鋪層（Layering）和局部重構(gòu)（Local Remeshing）。在文中使用鋪層和局部重構(gòu)的組合方式。

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

表2 不同網(wǎng)格下監(jiān)控點(diǎn)物理量幅值Tab.2 Physical quantities range of monitoring points in different grids

壓力耦合采用經(jīng)典的Simple算法，相比Simpler算法，雖然占用內(nèi)存較大，但計(jì)算速度快，魯棒性好。離散格式選擇二階迎風(fēng)格式，時(shí)間的離散采用全隱式格式。收斂的判斷標(biāo)準(zhǔn)為能量方程的殘差小于10-6，其他物理量的殘差均小于10-4，判斷計(jì)算達(dá)到穩(wěn)定的標(biāo)準(zhǔn)是圖1中各監(jiān)控點(diǎn)監(jiān)測(cè)的溫度、壓力、流速等數(shù)據(jù)在最后3個(gè)循環(huán)周期內(nèi)變化幅度不超過(guò)1%，即可認(rèn)定仿真達(dá)到穩(wěn)定。

2 模擬結(jié)果與分析

通過(guò)數(shù)值模擬得到了斯特林發(fā)動(dòng)機(jī)結(jié)構(gòu)內(nèi)溫度場(chǎng)、壓力損失和速度分布等情況，并分析了充氣壓力和加熱器溫度的變化對(duì)工質(zhì)流動(dòng)和做功的影響。

2.1 溫度分布

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

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

回?zé)崞鬏S向溫度在一個(gè)周期內(nèi)近似為線性分布，并且從低溫端到高溫端溫度逐漸升高，但任何一個(gè)徑向截面溫度都是隨著工質(zhì)流動(dòng)而變化的。由于回?zé)崞鲀?nèi)部密度分布不均，壓力非單調(diào)分布，導(dǎo)致熱端溫升快、溫降慢，冷端溫降快、溫升慢。

2.2 壓力分布

圖3是一個(gè)周期內(nèi)斯特林發(fā)動(dòng)機(jī)的壓力分布云圖，在3/4T時(shí)，即配氣活塞位于平衡位置向上止點(diǎn)運(yùn)動(dòng)，動(dòng)力活塞位于下止點(diǎn)時(shí)，膨脹腔和壓縮腔壓差可以達(dá)到12 kPa，而在1T時(shí)即循環(huán)開始，壓差只有3 kPa左右。從圖中可以明顯看出，壓縮腔和膨脹腔之前的壓降主要來(lái)源于回?zé)崞?，在大部分時(shí)間內(nèi)，軸向上回?zé)崞鞫即嬖诿黠@的壓力梯度，由于加熱器和冷卻器流道的存在，所以換熱器兩端也存在一定的壓差。應(yīng)用后處理工具可以得到相關(guān)壓力變化均值，圖4是3/4T回?zé)崞鬏S向壓力分布圖，3/4T時(shí)換熱器、回?zé)崞骱屠鋮s器兩端的壓降值分別為12 463 Pa、432 Pa和876 Pa。可以看出，與回?zé)崞鲀啥藟航迪啾龋訜崞骱屠鋮s器的兩端壓降很小。

圖3 一個(gè)周期內(nèi)斯特林發(fā)動(dòng)機(jī)的壓力分布云圖Fig.3 Pressure distribution of a cycle stirling engine

圖4 3/4T回?zé)崞鬏S向壓力分布圖Fig.4 Axial pressure distribution of 3/4T regenerator

2.3 速度分布

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

2.4 充氣壓力對(duì)質(zhì)量流率和做功的影響

在不改變其他參數(shù)的條件下，研究了充氣壓力（8×105Pa、12×105Pa、16×105Pa、20×105Pa）對(duì)回?zé)崞魃舷逻吔缳|(zhì)量流率和膨脹腔與壓縮腔PV功及輸出功的影響。如圖6所示，隨著充氣壓力的增大，充入的工質(zhì)質(zhì)量也增加，回?zé)崞鲀啥速|(zhì)量流率幅值有明顯的增大，且流率增幅與壓力增幅呈近似線性關(guān)系，從圖中也可以看出，回?zé)崞鳠岫说馁|(zhì)量流率與回?zé)崞骼涠说馁|(zhì)量流率相比較小，主要是因?yàn)闊岫藴囟群芨撸べ|(zhì)密度小，所以質(zhì)量流量相對(duì)較小。

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

圖6 不同充氣壓力下回?zé)崞鲀啥速|(zhì)量流率曲線Fig.6 Mass flow curve of regenerator under different charged pressure

圖7是壓縮腔、膨脹腔和輸出凈功隨充氣壓力的變化曲線，可以發(fā)現(xiàn)，隨著充氣壓力的增加，膨脹腔和壓縮腔PV功都在增加，因?yàn)榕蛎浨粌?nèi)工質(zhì)平均壓力高于壓縮腔工質(zhì)平均壓力，在活塞行程和頻率不變的情況下，一個(gè)循環(huán)內(nèi)PV功更多，所以凈功也

隨之增加。

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

2.5 加熱器溫度對(duì)質(zhì)量流率和做功的影響

對(duì)整機(jī)性能影響最大的是加熱器溫度。當(dāng)其他參數(shù)保持不變時(shí)，改變加熱器的溫度，研究了回?zé)崞魃舷逻吔缳|(zhì)量流率，膨脹腔與壓縮腔PV功和凈輸出功的變化情況。

由圖8可以看出，加熱器溫度越高，通過(guò)加熱器上邊界的質(zhì)量流率越低，主要是因?yàn)榧訜崞鳒囟茸兏撸涓浇患訜岬墓べ|(zhì)的比容越大，流過(guò)加熱器上邊界的質(zhì)量流率就越低。在整個(gè)腔室內(nèi)所充工質(zhì)質(zhì)量不變的情況下，冷卻器附近工質(zhì)的密度相對(duì)較大，流過(guò)冷卻器下邊界的質(zhì)量流量越大。如圖9所示，隨著加熱器熱端溫度的升高，膨脹腔內(nèi)工質(zhì)壓力也隨之增大，PV功增加，壓縮腔的PV功則會(huì)減少，輸出功增加，由此可以看出，加熱器溫度越高，斯特林發(fā)動(dòng)機(jī)膨脹腔PV功和凈輸出功都會(huì)明顯增大。

圖8 不同加熱器溫度下回?zé)崞鲀啥速|(zhì)量流率曲線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 結(jié)論

利用FLUENT軟件對(duì)β型自由活塞式斯特林發(fā)動(dòng)機(jī)進(jìn)行二維仿真，通過(guò)數(shù)學(xué)模型的建立和求解，在仿真實(shí)驗(yàn)穩(wěn)定后，得到了一個(gè)周期內(nèi)不同時(shí)間節(jié)點(diǎn)，熱機(jī)內(nèi)部工質(zhì)溫度場(chǎng)，壓力場(chǎng)和速度場(chǎng)的變化情況，通過(guò)監(jiān)控點(diǎn)的設(shè)置，重點(diǎn)關(guān)注了回?zé)崞魃舷逻吔缳|(zhì)量流率，各個(gè)腔室PV功和凈輸出功等變化，可以得到結(jié)論：

（1）加熱器和換熱器溫度、壓力等分布相對(duì)均勻，回?zé)崞鞫瞬苛鲌?chǎng)變化劇烈，腔室底部流體熱質(zhì)交換不足。回?zé)崞髦泄べ|(zhì)存在明顯的溫度、壓力梯度，且沿軸向近似為線性分布，與單向流相比，回?zé)崞鲀?nèi)工質(zhì)的流動(dòng)存在復(fù)雜的非單調(diào)非對(duì)稱的振蕩流。

（2）壓縮腔、膨脹腔中溫度和壓力分布較為均勻，但是壓力在回?zé)崞鲀啥思扒皇揖植繀^(qū)域分布不均，對(duì)照流場(chǎng)分布云圖可知，在這些區(qū)域存在明顯的渦流，這些渦流會(huì)增大流阻耗散能量，影響熱機(jī)整體熱效率。

（3）充氣壓力的增大導(dǎo)致腔室內(nèi)質(zhì)量流量增大，膨脹腔和壓縮腔工質(zhì)壓力也隨之增大，膨脹腔PV功增加的更多，凈功變大；隨著加熱器溫度的提高，換熱器與工質(zhì)換熱更加劇烈，回?zé)崞鲏航到档?，膨脹腔工質(zhì)PV功增大，壓縮腔PV功降低，所以增加加熱器溫度可以明顯增加斯特林發(fā)動(dòng)機(jī)的熱功轉(zhuǎn)換效率，但同時(shí)也要考慮加熱器材料的耐熱性能和缸體材料的耐腐蝕性能。