熱聲發動機中的交變流動換熱器特性研究

王海濤,羅靖,張麗敏,孫巖雷,陳遠航,常德鵬,陳燕燕,余國瑤,羅二倉

(1.中國科學院理化技術研究所低溫工程學重點實驗室,100190,北京;2.中國科學院大學,100049,北京)

熱聲熱機也稱回熱式熱機,因具有高效率、高可靠性、長壽命等優勢而逐漸受到重視[1-3]。換熱器是實現發動機內部工作氣體與外界熱量交換的核心部件,回熱器位于高、低溫換熱器之間,實現熱能與聲能的相互轉換。高、低溫換熱器與回熱器等部件緊密相連,工質在相鄰部件間往復穿梭。如圖1所示,在熱聲發動機的換熱器中,絕大部分工質在各交界面的突變截面處往復流動,顯著的時變慣性與“進、出口”效應使得熱聲熱機換熱器的特性必然與定常流動換熱器存在顯著不同。

圖1 熱聲發動機核心部件示意圖

由于物理過程本身的復雜性與多樣性,至今對交變流動及其換熱規律的認識與定常流相比仍然相差甚遠。理論分析工作主要基于線性或弱非線性熱聲理論,針對層流充分發展換熱器部件進行分析[4-7]。實驗與數值分析研究方面,代表性工作主要集中在20世紀80—90年代對斯特林熱機的加熱器或冷卻器的研究[8-12],以及近些年來的二維簡單板疊式換熱器的仿真與可視化[13-17]、三維瞬態仿真研究[18-20]等。前一階段的研究主要是借鑒定常流的分析方法,嘗試采用阻力系數與換熱系數來描述換熱器的特征,但所得到的表達式往往相差較大,缺乏普適性,這其中的緣由也鮮少有分析。后一階段則主要針對簡單換熱器進行CFD二維與三維模擬以及流場可視化測量,考察不同聲場強度下換熱器內的瞬態流動、換熱分布特征,以及換熱器端部的渦流問題等。這種細致分析的結果顯示出交變流動換熱器與定常流換熱器有著顯著的不同,首先是熱流密度的分布與聲場分布及系統其他部件之間的耦合特性非常強;另外,顯示出換熱器的流動與換熱特征主要受進出口效應影響。但是,這類研究未能在流場細節之外總結出規律,且受限于研究的具體熱機形式及實驗條件(針對的是簡單換熱器結構、單一聲場條件以及低功率密度下的研究),分析結果的系統性、普適性和實用性都顯得不足。

隨著大功率熱聲熱機發展的需求日益迫切,實際熱源耦合的換熱器設計越來越成為關鍵,然而對熱機內部交變流動換熱器特性的理解尚難以滿足發展需求。因此,有必要對熱聲熱機換熱器內部的換熱和阻力特征展開細致研究,明確換熱器內部換熱與流動損失特性,系統分析并總結特性規律,有效指導熱聲熱機的工程設計與應用。本文以一臺新型10 kW熱管耦合自由活塞斯特林發電機的發動機核心單元為分析對象,分別采用回熱式熱機通用設計準一維軟件Sage[21-22]和商用CFD計算軟件FLUENT[23],對高、低溫換熱器內的流動與換熱特征進行詳細分析。基于拓撲等效模型,對比分析兩種計算結果,提出適用于Sage計算的換熱器修正模型,提升Sage模型的計算準確性,指導熱聲熱機的工程設計。

1 熱管耦合發電機整體設計

1.1 熱管耦合自由活塞斯特林發電機結構

熱管耦合自由活塞斯特林發電機由實現熱聲轉換的發動機與聲電轉換的直線電機兩部分組成。發動機側需要實現熱管與發動機加熱器的結構與熱耦合,熱管有柱形直熱管(亦可大角度彎折)與異形熱管兩種,如圖2所示。柱形直熱管更為成熟,與自由活塞斯特林的耦合方式包括熱頭橫截面插入與軸向插入式兩種。常規耦合受限于發動機熱頭在橫截面與軸向長度上的過度緊湊,使得熱管冷凝段長度與發動機的尺寸適配性較差,只能用于低功率匹配,對于高功率、緊湊型發動機結構的應用存在限制。異形集成式熱管加熱器能實現發動機與熱管之間的高效耦合,發電機功率可達幾十千瓦。熱管結構與發動機的承壓壁集成,異形熱管具有顯著結構依賴特點,通用性差,結構復雜,工藝要求高,目前尚不成熟。

針對以上問題,本文設計出一種高效軸向直熱管自由活塞斯特林發電機,通過降低加熱器孔隙率與降低系統工作頻率,提高加熱器內工質位移幅度,從而有效增長加熱器長度,即增長熱管冷凝段長度。本文所研究的發電機結構優化設計基于Sage軟件完成,如圖3所示,為基于熱力-電磁設計完成的整機基本結構示意圖。直熱管插入如圖3(b)所示的套筒內,氣體在銅塊的矩形窄縫中進行往復流動換熱,同時為實現熱管加熱器與發電機的高效耦合,負載發電機部分如圖3(a)所示,采用對置電機的形式,實現單側發電機熱耦合的同時減小系統振動。

(a)發動機整體結構 (b)熱管加熱器局部

經設計優化后的換熱器結構參數:翅片式高溫加熱器長度為250 mm,流道寬度為1.3 mm,流道高度為10 mm,翅片平均厚度為7.986 mm,流道數為108,孔隙率為0.046;回熱器長度為52 mm,采用絲綿填充,孔隙率為0.897;室溫換熱器為管殼式,長度為70 mm,管徑為2 mm,孔隙率為0.08,流道數為768。

1.2 整機性能參數

熱管耦合發電機整機優化后的額定工況性能參數如表1所示。額定設計發電功率為10.77 kW,對應加熱器吸熱量為30 kW,冷卻器排熱量約為19 kW。

表1 熱管耦合自由活塞斯特林發電機性能參數

發動機側冷卻器、回熱器與加熱器核心段的壓力、體積流率波動幅值分布如圖4(a)所示,聲功、聲阻抗相位分布如圖4(b)所示。為適應熱管冷凝段長度,加熱器長度為250 mm,可見加熱器內明顯的壓力幅值下降與聲功降低,即加長的換熱器明顯增加了聲功損失。

(a)發動機內壓力波動與體積流率分布

基于以上熱管耦合發電機結構與發動機核心部件額定工況,本文以冷卻器、回熱器與加熱器組成的核心單元為研究對象,耦合熱聲轉換,重點分析加熱器與冷卻器內的流動與換熱特性。

2 換熱器計算模型

2.1 換熱器拓撲等效

為實現對換熱器內流動與換熱特征的系統分析,同時簡化計算過程,將三維流體計算區域通過拓撲等效簡化為二維平面結構。加熱器與冷卻器均為窄縫結構,相鄰的部件為回熱器與空腔(加熱器側連接膨脹腔,冷卻器側連接壓縮腔),回熱器為多孔介質結構,具有顯著的均勻化流動效果。因此,回熱器與冷卻器的拓撲對應性即可解耦,即將加熱器與冷卻器窄縫簡化為二維結構時,與回熱器的結構匹配只需要滿足界面上的孔隙率對應即可。由于換熱器孔隙率相對于回熱器與空腔足夠小,因此忽略換熱器窄縫三維分布對回熱器與空腔內射流與二維射流的差異。以狹長窄縫二維等效為基礎,以孔隙率最小的加熱器一條完整窄縫為最小計算單元,以窄縫寬度為基準。冷卻器原設計為管束,二維近似通過水力直徑、換熱面積及孔隙率相近原則進行等效,等效后的窄縫數量取相對加熱器鄰近整數,適當調整水力直徑以匹配相對孔隙率,回熱器寬度與空腔寬度以相對換熱器孔隙率確定,加熱器、回熱器、冷卻器的長度保持不變,空腔長度則以相對容積確定。從而將三維結構等效為二維平面結構。壓縮腔與膨脹腔容積根據CFD動活塞運動要求進行容積擴展,以壓縮腔與冷卻器界面、膨脹腔與加熱器界面上的壓力、體積流率與整機完全相同為原則,將計算邊界外延。

圖5為按照拓撲等效原則獲得的自由活塞斯特林發電機二維模型,結構參數如表2所示。基于該二維模型,可同時采用sage與CFD進行仿真計算,sage一維模型計算結果即可以與原整機設計計算結果進行對比校驗,以檢驗拓撲簡化過程的合理性,同時也可為CFD局部仿真提供邊界條件。基于CFD仿真對換熱器時變換熱特征進行詳細分析,并與sage一維模型計算結果進行對比,考察sage一維模型計算的有效性。

表2 等效換熱核心單元結構參數

圖5 發動機核心單元拓撲等效模型

圖5所示計算模型中,高溫加熱器的壁面溫度為823 K,低溫冷卻器的壁面溫度為333 K,其它壁面邊界設置均為絕熱。軸向計算邊界為避免開口邊界帶來的入流焓流不確定性,采用雙活塞邊界進行計算。活塞邊界即設定往復振蕩虛擬活塞邊界面,模擬活塞對熱聲核心單元的往復壓縮、膨脹效應,能量從邊界上以時均聲功形式進出系統,活塞邊界可設置絕熱,無不確定性能量,因此對于系統內部的能量分布分析精度更高。Sage一維模型可采用虛擬活塞邊界進行計算;CFD仿真中,兩端虛擬活塞邊界則必須采用動網格,為保證活塞相對平衡位置不變,還需要進行模擬充氣與活塞同步的初始過程模擬。活塞同步后,即可加載壓縮活塞與膨脹活塞運動的速度控制,活塞速度隨時間的變化函數如下

upiston(t)=ωxpiston1cos(ωt+θpiston)

(1)

式中:xpiston1、θpiston分別表示活塞面的位移幅值與相位;ω是角頻率。

2.2 Sage一維模型

Sage是由美國Gedeon等于1995年根據MS-DOS編寫開發的一款商業計算軟件,主要是針對回熱式熱機的通用設計計算軟件。其采用圖形化界面將與實際物理組件相對應的氣體流動、傳熱和其他建模細節封裝在一些特定的模型組件中,用戶通過將封裝好的組件按照特定的需求進行相互連接,即可對復雜的回熱式熱機系統建立仿真模型,利用質量、動量和能量守恒方程、本構方程將各部件的參數進行有機耦合,最終實現整機的求解與多參數協同優化。圖6給出圖5所對應的Sage一維計算模型。

ρstdy—系統載荷;Pphsr—體積位移;mGt—進出口的質量流;Qstdy—線性熱源;T0—低溫熱源;Th—高溫熱源。

回熱器模型選擇軟件內嵌的絲綿多孔介質模型,其阻力的本構方程如下

f=a1/Re+a2Rea3

(2)

式中:a1=25.7α+79.8;a2=0.146α+3.76;a3=-0.002 83α-0.074 8;α=ε/(1-ε);ε為孔隙率。

換熱特性的本構方程如下

Nu=1+b1Peb2

(3)

式中:b1=0.186α;b2=0.55;Pe為貝克來數,即雷諾數Re與普朗特數Pr的乘積。

等效模型中的換熱器均為平板模型,采用軟件內嵌的板式換熱器模型,基于線性熱聲理論平板模型獲得板式換熱器模型的阻力與換熱特征層流,湍流則借用湍流平板穩態流動特性,本構方程如下

f=0.11(ε/dh+68/Re)0.25

(4)

Nu=0.035Re0.75Pr0.33

(5)

式中:dh為絲綿等效水力直徑。

Sage一維模型中可通過設置額外的阻力系數,對模塊進行整體阻力修正,阻力修正通過定義以下壓力梯度表達式中的K實現。即計算模塊的壓降由模塊特征阻力系數f與附加可用戶自定義的阻力系數K決定。

在流動方向上的壓降表達式如下

(6)

式中:L為模型長度;ρu|u|/2為流體動能。

2.3 CFD模型

CFD仿真使用ANSYS FLUENT 2021 R1版本。采用活塞邊界動網格進行瞬態計算,發動機內氦氣的交變流動計算屬于小馬赫數可壓縮流動。將氦氣設置為理想氣體,采用基于壓力法的求解器計算,選用PISO 算法,壓力使用PRESTO!離散格式,其他均采用二階迎風離散格式,時間離散采用二階時間差分。換熱器內流速較高,以幅值雷諾數作為判斷基準,加熱器雷諾數幅值變化范圍在3 000以上,明顯高于2 300,因此可采用k-ε湍流模型。模型中計入氦氣物性隨溫度的變化,物性公式列于表3。回熱器采用熱平衡多孔介質模型,只考慮回熱器內流阻模型。阻力系數根據Sage模型中所用的阻力系數公式,通過局部擬合轉換為黏性阻力系數與慣性阻力系數表征的阻力特性,獲得FLUENT中多孔介質模型設置參數1/K(黏性阻力系數)與C2(慣性阻力系數)。

表3 氦氣物性參數

2.4 CFD二維模型網格與時間尺度無關性檢驗

基于額定工況,進行加密網格與時間尺度的計算,同時對壓力、體積流率、溫度、聲功、換熱量的波動與時均值進行監測,波動量幅值、相位及時均量均通過FFT分析獲得。圖7與圖8給出了典型參數的影響。其中網格尺度影響相對較小,網格量大于19 800后,對計算結果影響可忽略,因此后續計算網格均采用19 800這套網格模型。

圖7 壓力波動幅度延程分布隨網格數的變化

圖8 不同時間步長下時均能量守恒性隨計算時間的變化

時間尺度對波動參數幅度的影響較小,但是對時均能量的影響較大。由于計算區域為封閉系統,活塞邊界的時均聲功差異應該等于高、低溫換熱器的時均換熱量差異,但CFD二維模型的時均能量計算發現存在明顯的能量不平衡,且該能量不平衡性直接與時間尺度相關。如圖8中黑色線代表的每周期計算500個點時,系統有600 W左右的不平衡;每周期計算達到5 000個點時,系統能量不平衡降低到170 W左右。繼續減小時間步長,會導致計算量過大,相較于10 kW量級的系統典型時均能流,計算誤差可接受,本文選擇一個周期計算5 000步。對小時間步長的要求意味著熱聲系統中時均能量的準確求解較為困難,這主要是因為能量的主要特征為波動量,而時均值為多參數耦合值,相較于波動量為高階小量,且與參數間的相位關系直接相關,因而計算精度對時間尺度要求非常高。

2.5 額定工況Sage模型與CFD模型計算對比

表4給出了基于圖5的Sage模型與CFD模型計算結果對比。其中:Pcomp表示壓縮腔入口壓力一階幅值;Pexp1表示膨脹腔1入口壓力一階幅值;PVchxout表示冷卻器出口截面時均聲功;PVhhxin表示加熱器入口截面時均聲功;Qchx表示冷卻器壁面時均換熱量;Qhhx表示加熱器壁面時均換熱量。兩種模型計算均采用相同活塞邊界,因此壓縮腔與膨脹腔活塞界面處的體積流率完全相同。二者計算的時均能量相當,差異在10%附近,但壓力差異相對較大,即內部阻力特性差異較大。

1、20世紀50年代初的探索期是箏樂演奏技術發展的一個分水嶺，打破了“右手演奏，左手和聲”的觀念。這個時期的代表作《慶豐年》。

表4 兩種模型相同截面計算參數的對比

3 換熱器特性分析

為詳細分析發動機換熱器的換熱特性,基于固定結構參數與額定工況,保持其他參數不變,進行不同位移幅度與不同聲場相位下的換熱器內流動與換熱特性分析。其中位移幅值的變化是通過等比例改變活塞速度來實現的,從額定工況的0.1倍變化到2倍,改變位移幅度時,計算區域內的聲場相位分布變化非常小。聲場相位的變化通過保持額定工況位移幅度不變、改變兩個活塞之間的相位差(在0～2π范圍內變化)實現。

3.1 換熱器內換熱分布特性

如圖9所示的額定工況下換熱分布對比可知,CFD二維模型與Sage一維模型的計算對時均溫差與時均熱流密度qw分布的捕捉基本相似,數值相近,說明Sage一維模型可以合理地獲得換熱器內時均換熱分布與集中特征。但是,與CFD二維模型計算結果仍然存在明顯差異,該差異與溫度場的局部二維特性有關,也與流阻損失產生的場能量分布差異有關。計算聲場下,換熱器表現為靠近回熱器一側換熱更為顯著,由于為等壁溫條件,因此對應的時均溫差也更大。

(a)額定工況下換熱器內時均溫差分布

圖10中不同位移幅度下的時均熱流密度分布變化表明,隨著換熱器內位移幅度低于額定工況,換熱器內出現無效換熱區間(即時間換熱為0),超過額定工況后,換熱器換熱強度分布出現近等比例抬升。Sage一維模型計算與CFD二維模型計算結果結論一致,說明換熱器有效長度必須與其內部工質位移幅度相匹配,位移過小,換熱器有效換熱區域變小,印證了熱聲系統中的換熱器進、出口效應本質與本征緊湊性特征。

圖10 不同位移幅度下時均熱流密度分布的CFD二維模型計算結果

圖11將時均熱流密度與時均氣固溫差代入努塞特數計算式,獲得時均等效Nu(Nu0),冷卻器中的Nu0顯著高于加熱器,且延程分布特性更為顯著,對速度也更為敏感。這進一步說明了熱聲發動機中換熱器換熱特征具有顯著的局部特征。

圖11 不同位移幅度下Nu0的CFD二維模型計算結果

鑒于Sage一維模型能較好地捕捉換熱器內的換熱特征,且計算速度顯著快于CFD二維模型,因此以下采用Sage一維模型計算不同聲場相位時的換熱器內換熱特性。如圖12所示,調節計算邊界上的活塞相位差(Δθ),得到的各工況下聲場相位分布圖。圖13給出相應的時均聲功分布情況,可見核心單元內的能流分布發生了顯著變化。

圖12 不同邊界活塞相位差下的聲阻抗角分布

圖13 不同邊界活塞相位差下的時均聲功分布

由圖14可知,雖然加熱器與冷卻器外壁面溫度保持不變,即加熱器外壁面為823 K,冷卻器外壁面為333 K,但聲場的變化不僅改變了熱流密度分布特性,還決定了換熱器是吸熱還是放熱。因此,熱聲系統中的換熱器換熱特性是聲場能流驅動的,由時均焓流決定時均換熱,而不是溫差驅動換熱。時均溫差是時均能流變化產生時均換熱的結果而不是換熱的驅動力,這與穩態流動完全不同。

圖14 不同邊界活塞相位差下的時均熱流密度分布

3.2 阻力特性分析

圖15給出了額定工況下冷卻器、回熱器、加熱器中壓力、體積流率、聲阻抗相位角以及時均聲功的分布對比。其中最大的差異是壓力波動在回熱器與換熱器的交界面上的壓力差異。這是因為Sage作為一維模型無法直接獲得復雜流動變化帶來的局部損失,且已有換熱器模塊不能針對部件間的相對結構特征進行自動局部阻力調整,而是需要人為設定相關損失。CFD二維模型仿真可以獲得突變截面流場特性的變化,包括射流的影響。從計算結果可以看出,由于加熱器孔隙率非常小,在回熱器與加熱器界面上產生了顯著的一階波動壓力損失,時均能量也可看到明顯的聲功損失。換熱器與兩側腔體的連接處損失則相對小得多。這是因為,加熱器的高速射流在回熱器這一高阻力部件中產生的損失要比射流在空腔中顯著得多。這是換熱器計算模型中Sage一維模型與CFD二維模型的最大差異來源。

(a)壓力

(a)壓力

3.3 Sage模型阻力修正

基于CFD模型計算結果,可對換熱器與回熱器的突變界面處的流阻損失進行特征擬合,反代入Sage模型計算模型進行修正,將修正結果再次與CFD模型計算結果進行對比,即可判斷修正模型的適用性。另一方面,基于穩態突變界面局部損失系數與Swift計算方法亦可對突變界面特性進行修正。以下將針對這兩類方法進行對比驗證。

由Swift方法對熱聲系統中噴射泵相關理論分析[24]可知,氣體往復穿梭變截面時,半個周期經歷流道擴張,半個周期經歷流道收縮,變截面處的局部阻力系數K應具有時間依賴性,即為時間的函數,因此

(7)

使用正弦速度

U(t)=|U1|sinωt

(8)

將式(8)代入式(7),可求得平均壓降為

(9)

則其微分形式為

(10)

根據穩態流動中面積比與管道截面突變局部阻力系數表[25],計算模型面積比,即可獲得交變流動界面等效阻力系數值。

根據CFD模型計算所得壓力波動分布進行等效阻力系數修正嘗試了兩種方法:一種是針對波動壓力幅值與速度幅值進行修正;另一種是針對波動壓力與速度的有效值進行修正。其中波動壓力幅值修正通過下式計算出CFD二維模型的波動幅值壓降與速度幅值之間的關系,再計算出總阻力系數Ktotal

(11)

(12)

對于有效值計算法則,根據下式計算CFD二維模型的換熱器壓力與速度瞬時值,再計算出總阻力系數Ktotal

(13)

通過以上方法,針對冷卻器與加熱器在不同速度幅值下的總阻力系數進行計算。Swift計算方法獲得的阻力系數與速度無關,CFD模型數據處理獲得兩種阻力系數是速度的函數。將以上3種修正方法進行對比發現,采用有效壓降法獲得的修正最為有效。基于有效壓力修正方法使得Sage模型波動壓降與聲功損耗均高于CFD模型計算結果,但差異存在比例關系。因此,基于有效壓力的修正方法需再引入一個阻力系數修正因子

φ≈0.7

(14)

冷卻器修正后的有效總阻力系數如下

(15)

加熱器修正后的有效總阻力系數如下

(16)

以CFD模型計算為基準,各種修正方法所得加熱器與冷卻器的波動壓降以及時均聲功的誤差見圖17。

(a)冷卻器波動壓降誤差

由圖17可知,經修正因子修正后的有效阻力系數可使Sage一維模型與CFD二維模型計算的換熱器兩端壓差計算差異縮小到20%以內,即Sage一維模型的阻力特性的精度有所提高。換熱器聲功損失修正精度長加熱器比短冷卻器明顯差,這是因為長加熱器的分布特性更難用端部集總壓差簡單修正得到。如表5所示,整體性能對比表明,Sage一維模型修正后的結果明顯向CFD二維模型計算結果靠近。證明通過有效壓降與有效速度獲得的總阻力系數可用于修正Sage一維模型的阻力特性,可有效提高工程設計精度。

表5 修正前后Sage模型與CFD模型參數對比

基于以上分析,可以建立一套提升實際復雜結構發動機的分析效率與精度的方法。首先利用CFD仿真對真實換熱器與回熱器組成的核心單元結構進行給定聲場下的流動阻力特性分析,由于流阻分析的計算時間穩定較快,可以不需要等待系統熱容與時均換熱穩定即可進行分析,因而即使對于復雜三維模型的分析計算同樣具有工程可行性。所得到的換熱器阻力特征可用于修正Sage一維模型中的對應阻力系數。Sage一維模型中換熱特性的計算已經可以較為準確地捕捉換熱分布與集中特性,因而無需依賴三維CFD仿真模型的熱穩定計算。時均換熱特性的CFD計算由于對時間尺度的高要求,以及較大系統熱容穩定時間尺度之間存在的顯著矛盾,使得完全依賴三維CFD仿真進行換熱器完整特性分析的計算時間過長而難具工程指導意義。

4 結 論

本文基于10 kW熱管耦合斯特林發電機整機設計,進行發動機核心部件的拓撲等效,分別采用Sage一維模型與CFD二維模型兩種計算方法,對核心單元的換熱與阻力特性進行模擬研究,通過比較兩種模型的計算結果,得出以下結論。

(1)對換熱器內的局部換熱特征以及集總換熱特征計算結果對比,Sage一維模型與CFD二維模型具有相近結果,Sage一維模型的分析結果具有工程指導意義。換熱器換熱特性主要由聲場驅動,這可能是換熱器換熱特征分析中,對本構方程的準確性要求相對較低的原因。

(2)對換熱器內阻力特性的計算結果對比表明,Sage一維模型與CFD二維模型計算結果相差較大,換熱器射流導致的回熱器界面顯著損失是產生這一差異的主要原因。因此,為提高Sage一維計算結果的精度,需對利用CFD對實際換熱器與回熱器界面流阻進行仿真計算,提取阻力特性后在Sage一維模型中進行修正,即可較好地提高Sage一維模型對實際系統的預測準確度。

(3)進一步明確了熱聲熱機換熱器特性的顯著空間分布特征以及工況依賴特征,這些特征在熱聲換熱器中具有相似性和普遍性。但是,對于特征的量化又因其空間與工況敏感性而難以像穩態流動那樣獲得通用的本構方程,因而熱聲熱機換熱器需要針對具體設計進行特性分析。