辛鵬飛, 苗建印, 匡以武, 張紅星, 王 文
(1. 上海交通大學 機械與動力工程學院,上海 200240; 2. 北京空間飛行器總體設計部 空間熱控技術北京市重點實驗室,北京 100086)
隨著電子產品集成度和功率的不斷提高,電子器件的散熱壓力也逐漸增大,散熱需求達到了 1 000 W/cm2量級[1-2],泵驅兩相系統結合微通道換熱技術可以獲得較高冷卻能力,在航空航天、制冷空調、通信、燃料電池等領域為保證電子器件的正常運行都有著廣泛的應用和發展前景[3].泵驅兩相流在運行過程中存在壓力和流量失穩導致的運行可靠性問題,主要分為靜態不穩定以及動態不穩定兩種類型[4].流量漂移(Ledinegg不穩定)是靜態不穩定中經常出現的情況,它是一種系統層面的不穩定性[5].對于單個通道來說,負斜率區域導致在一定壓降范圍內,壓降與流量之間不再是一一對應的關系.對于并聯管道來說,由于進出口壓降保持一致,因而在流量漂移區間內可能存在多種不同的流量組合.流量的變化使得部分管道的出口為過冷或者過熱,從而影響換熱器的性能.
在針對多個分散的熱控單元進行液體冷卻時,多個散熱器的并聯運行不可避免,散熱器間流量的相互影響會導致其換熱能力產生一定波動特征.Akagawa等[6]觀察了并聯蒸汽發生器管道中的流量漂移現象,進行了理論分析并提出了穩定性的判定標準.Minzer等[7-8]研究了太陽能直接蒸汽發生器并聯管道中流量分配,理論分析了兩并聯管道中的流量分布及穩定區域.張炳雷等[9]在低高寬比微通道實驗中觀察到了流量漂移不穩定現象的存在,發現流量漂移易發生在大流量和低進口溫度下.Oevelen等[10-11]利用廣義特征值的方法分析大量并行微通道下的流量分配,并考慮了側壁面導熱對流量分配的影響.
在自然循環系統中,也存在多換熱器之間流量分配問題.楊瑞昌等[12]建模分析了自然循環壓水堆一回路系統內并聯換熱管內正流和倒流的流量.郝建立等[13]通過分析自然循環蒸汽發生器一回路倒流現象,發現工作于負斜率區域是流量漂移的原因.彭傳新等[14]對并聯通道兩相自然循環系統中流量漂移導致的流動失穩進行實驗研究.
文獻中針對多并聯系統中的穩定性分析采用的方法主要包括線性變換和拉普拉斯變換[6,8,15]、廣義特征值解法[10]、矩陣分析[16]等.
現有文獻中研究對象多為動力系統和供暖系統的并聯管道的流量漂移現象,而對于小型泵驅冷卻系統多個并聯換熱器件的研究相對較少.與并聯管道所不同的是,并聯熱沉間由于流動換熱結構、傳熱所導致的密度和流型變化等原因使得影響流量分配因素更多,出現的特征也有所不同.此外,絕大多數研究的工質以水為對象,在泵驅兩相液冷回路中還常采用氨和有機低沸點流體作為工質,其物性與水也有所不同.
本文針對泵驅兩相冷卻系統的并聯多蒸發器進行建模,采用氨作為冷卻工質,研究了在較高熱流密度下影響單個蒸發器流動特征曲線的因素,并在此基礎上研究了相關因素對多蒸發器系統流量分配的影響規律,分析了流量變化對于蒸發器溫度分布的影響.
多蒸發器兩相泵驅系統如圖1所示.圖中:Lin表示蒸發器進口連接管道長度;Lout表示蒸發器出口連接管道長度.換熱器尺寸參照文獻[17]中蒸發器實驗模型.為簡化計算,對多蒸發器系統進行以下假設:① 蒸發器結構相同;② 蒸發器內部各微通道間流量均勻分配;③ 外部管道內壁光滑;④ 忽略蒸發過程中飽和壓力變化對物性參數的影響;⑤ 忽略換熱器與外界環境之間的漏熱.

圖1 多蒸發器熱沉并聯模型
本文計算中將單個蒸發器簡化成一個單通道換熱單元,在此基礎上考慮工質在蒸發器內的各個局部壓降,包括從管道進出蒸發器靜壓室、靜壓室的截面變化以及靜壓室進出通道等過程中引起的局部阻力.計算中邊界條件設置為出口壓力保持恒定,每個工況下進口流量恒定,同時,多蒸發器連接在同一集箱上,進出口壓差相同.
對于水平方向布置的蒸發器,其一維動量方程為
(1)
式中:G為蒸發器的質量流速;u為通道內流速;z為計算點到通道入口的距離;p為壓力;fw為壁面摩擦系數;t為時間.
對于并行蒸發器,蒸發器間的總流量保持不變,滿足:
(2)
式中:W為質量流量;N為并行蒸發器的數目.
根據假設,將動量方程沿通道長度進行積分,得到蒸發器的集總參數模型:
pout=pin-pout-Fi(Wi)
(3)
式中:L為蒸發器的等效長度;A為等效截面積;Δpa、Δpf、Δpc、Δpe分別為蒸發器內的加速壓降、摩擦壓降[18]、突縮壓降[19]、突擴壓降[20];pin、pout為通道進口、出口壓力;Fi(Wi)為蒸發器中各部分壓降之和.
對于單相區域,其摩擦壓降可表示為
(4)
Po=24(1-1.355 3β+1.946 7β2-
1.701 2β3+0.956 4β4-0.253 7β5)
(5)
式中:pf,sp為單相區摩擦壓降;Po為泊肅葉數;μ為運動黏度;ρ為流體密度;ν為比體積;Dh為水力直徑;β為通道短邊與長邊的比值.
兩相區域壓降摩擦壓降表示為
(6)
(7)
(8)
式中:pf,tp為兩相區摩擦壓降;C為Chisholm因子,在層流狀態下為5;x為干度;l表示液相;g表示氣相.
蒸發器中突縮壓降Δpc、突擴壓降Δpe分別表示為
(9)
(10)
式中:γ表示測試通道截面積與集箱截面積的比值;Cc為收縮系數;Ke為膨脹系數;xe,out為蒸發器出口干度;G1為突縮前質量流率;G2為突縮后質量流率.
同時,考慮到蒸發器通道以及壁面在熱流下溫度變化,將流體和蒸發器壁面沿流動方向分為若干個節點,如圖2所示,考慮節點間的導熱以及與流體節點間的對流換熱.圖中:Tf為流體溫度;Tw為壁面溫度;i為蒸發器中通道的編號.

圖2 蒸發器節點換熱示意圖
壁面導熱能量方程為
(11)

忽略黏性耗散以及流體間的擴散項,流體側的能量方程為
(12)
式中:Hf為流體焓值.
由于兩相區的存在,將流體的溫度與焓值的關系用分段函數來表示[11]:
(13)
式中:Tsat為飽和溫度;cl為液體比熱容;cg為氣體比熱容;Hl為飽和液體焓值;Hg為飽和氣體焓值.

(14)
式中:η0為肋系數,文中取99.12%;hch為蒸發器內的微通道高度;h為對流換熱系數,采用Bertsch等[21]提出的關聯式,其表達式為
h=hNBS+hconv,tpF
(15)
hNB為核態沸騰項;hconv, tp為對流項;S為核沸騰項的抑制因子;F為對流項因子.
對于壁面的能量方程,考慮到熱量在微通道壁面中向前后兩個方向的傳遞,對擴散項采用中心差分的格式;對于流體能量方程中的對流項,采用一階迎風格式[22].
氨是一種自然制冷劑,在微通道中的流動換熱研究較少.氨具有僅次于水的汽化潛熱,工作壓力較高,表面張力較小,全球變暖潛能值(GWP)和消耗臭氧潛能值(ODP)為0[23],在熱泵制冷[24]、航空航天[25]、環路熱管[26]等領域有著較多的應用,本文將采用液氨作為冷卻工質,選取的飽和溫度為25 ℃.
蒸發器的材料為銅,表面鍍鎳,結構參數為:wch=0.4 mm;hch=0.6 mm;L=10.8 mm;wf=0.3 mm;底層厚度wt=1 mm.采用局部加熱的形式,加熱區域長度為6 mm,位于蒸發器底部中央.
根據李雅普諾夫穩定性的理論,省略非線性方程中的高次項,利用近似一次線性系統的穩定性準則來評價其穩定性[27].并聯熱沉進出口壓力相同,可表示為
(16)
代入式(3)可得:
(17)
令Wi=Wi0+wi,其中Wi0為穩態時i蒸發器內的流量,wi為穩態時受到的擾動量.令W=[W1W2…WN]T, 則初始狀態點為W0=[W10W20…WN0]T,在W0處做泰勒展開,
(18)
其中,特征矩陣B為

(19)
令ΔW=W-W0=Keλt,代入式(11)可得:
(20)
式中:λ為特征根;m為等效慣性系數(L/A);K為系數矩陣.
當系統特征矩陣B的特征根λ全為負值時,系統是穩定的;當至少有一個特征根λ>0時,系統是不穩定;當僅存在零實部的特征根,則不能判斷原系統的穩定性.
圖3表示了計算的流程圖.圖中:Tw,new為迭代后壁面溫度;Tf,new為迭代后流體溫度;Hf,new為迭代后流體焓值;RH為流體焓值迭代設定殘差;RT為溫度迭代設定殘差.對于單個蒸發器,沿流動方向將蒸發器分為若干個節點,先通過Gauss-Seidel方法迭代計算蒸發器的壁面和流體的溫度分布,再根據每個控制體內流體的溫度求解對應壓降,結合蒸發器的局部壓降得到系統的流動特征曲線.并根據單個蒸發器的流動特征曲線,得到系統的特征曲線以及蒸發器間的流量分配情況.

圖3 計算邏輯圖
在進口溫度為5 ℃時,對蒸發器的壁面溫度進行網格無關性檢驗,如圖4所示.在收斂標準一致的情況下,隨著劃分網格數目Nc的增加,其溫度分布逐漸趨于一致.為保證計算精度以及減小計算量,沿流動方向將劃分216個一維網格.

圖4 網格無關性檢驗
將模型計算結果與Miglani等[28]以水為工質、通道尺寸為1 mm×1 mm×55 mm的微通道流量分配實驗數據對比,總體積流量Q=9.3 mL/min,入口溫度Tin=88.5 ℃,飽和溫度Tsat=100.8 ℃,如圖5所示.圖中:Wtot為總流量;Wi/Wtot為流量比;PT為總加熱功率.可以發現,模型計算結果與實驗模擬結果比較吻合,最大偏差為3.2%.綜上所述,此模型滿足計算要求.

圖5 模型驗證
與常規通道相似,進口過冷度ΔTin、熱流密度Q″、飽和溫度Tsat以及通道水力直徑Dh都會影響蒸發器的流動特征曲線,如圖6所示.圖中:Δp為蒸發器壓降;灰色虛線表示負斜率區域的消失的邊界.圖6(a)中,隨著進口過冷度的減少,負斜率區域傾斜程度逐漸降低,負斜率區域對應的長度先增加后減小,進口過冷度增加到一定數值后,負斜率區域消失.加熱功率的影響如圖6(b)所示,隨著加熱功率的增大,負斜率區域從無到有,傾斜程度逐漸增大,且對應的負斜率區間長度逐漸增大.在進口過冷度一定時,飽和溫度的影響如圖6(c)所示,飽和溫度的增加會降低工質間的密度差,從而降低通道內的孔隙率,使得兩相摩擦壓降降低[29],飽和溫度的提高會降低負斜率區域的長度和傾斜程度,提高穩定性.通道尺寸的影響如圖6(d)所示,微通道水力直徑Dh的減小會增加整體的壓降,卻減小了負斜率區間的傾斜程度,原因在于增加了通道進口的局部壓降,類似于節流的效果,可以提高系統的穩定性.可以發現,虛線所代表的負斜率消失對應的質量流速隨著進口過冷度的減小、加熱功率的增大、飽和溫度的增大以及微通道水力直徑減小而增大,而飽和溫度的增大對應著汽化潛熱的降低,這與Kuang 等[30]得到的理論表達式變化趨勢相同.

圖6 熱力參數和通道尺寸對蒸發器流動特征曲線的影響
在相同的過冷度15 ℃和加熱功率下,水和氨的流動特征曲線如圖7所示,可以看出,以水為工質的整體壓降大于以氨為工質,且水的負斜率區域的范圍和傾斜程度遠大于氨.由于水的液氣密度比為 1 600,遠大于氨的77,導致在兩相出現區域水的壓降比較大.在相同的過冷度下,由于氨的比熱容較大,使得氨的沸騰起始點小于水.但在純液體工況下,由于相同質量流速下,氨的流速較大,使得其壓降大于水.

圖7 水和氨的流動特征曲線
蒸發器進出口連接管長度對系統流動特征也會產生影響,如圖8所示.當進口管道長度增加時,負斜率對應的范圍以及傾斜程度明顯減少,主要因為進口單相段壓降是流量的單調函數,能夠提高系統的穩定性;當出口管道長度增加時,負斜率對應的區間以及傾斜程度會明顯增加,原因是出口管道多為兩相混合物,其壓降與流量之間的函數關系是非線性的,從而使得系統的穩定性降低.圖9表示在連接管道總長為220 mm時,進口管道內的壓降Δpin、出口管道內的壓降Δpout與總壓降Δptot的比值.當出口長度較長時,管道內的壓降會大于蒸發器內部的壓降,而出口管道流動特性曲線中存在的負斜率區域會增加系統的不穩定性;進口長度較長時,進口管道壓降較大的占比會提高系統的穩定性.當進口管道長度增加到一定數值后,可以抵消蒸發器產生的負斜率區間,但會增加系統整體的壓降和占用的空間.因此,對于特定的多熱源冷卻系統,合理的安排蒸發器的位置有助于提高系統的穩定性.

圖8 進出口管道長度對流動特征曲線的影響

圖9 進出口連接管道內壓降與蒸發器內壓降比值
由于蒸發器流動特性曲線中負斜率區域的存在,在滿足質量守恒與動量守恒的前提下,蒸發器內部的流量分配有著多種的可能性.當進口過冷度ΔTin=20 ℃,總質量流速Gtot=1 500 kg/(m2·s)時,蒸發器中的流量分配存在3種情況,圖10中兩個換熱單元的流量壓降曲線相交于各自的負斜率區域,其對應的質量流速和始終為 1 500 kg/(m2·s),兩條流動特征曲線的交點即為可能的流量分布情況.圖中:Δp1為蒸發器1壓降;Δp2為蒸發器2壓降;x1為蒸發器1出口干度.當兩個蒸發器完全相同時,其流量組合為(1 101/399) kg/(m2·s)或(750/750) kg/(m2·s).將工況點代入式(20)計算特征矩陣的特征根,工況點1、3特征根實部為負值,相對穩態流量變化量ΔW在擾動后逐漸降低,是穩定工作點;工況點2為正值,ΔW在擾動后逐漸增大,為不穩定工作點.圖中:Wa表示沸騰起始的質量流速;Wb表示負斜率區域消失的質量流速.從而將流動特征曲線分為3個區域,III 區位于單相液區,II 區以及 I 區大部分位于兩相區,單相氣區在整體流量變化中占據的范圍較小.

圖10 兩并聯熱沉間的流量漂移
流量漂移的存在使得部分蒸發器得到更多的流量,從而出口為單相狀態;部分蒸發器得到較少的流量,出口流體為高干度甚至過熱狀態.如上圖初始狀態兩個蒸發器流體的出口干度均為0.02,當發生流量漂移后,蒸發器流體的出口干度變為0.12或0.兩相區流體的干度對流體對流換熱系數有著比較大的影響,因此流量的變化會對蒸發器整體溫度分布有著比較大的影響.


圖11 流量漂移前后蒸發器壁面溫度、流體溫度以及單位長度對流換熱量的變化
通過對上述進行分析可以發現,在部分區域加熱時,發生流量漂移之后對于蒸發器整體溫度分布來說,不一定是不利的.如上述情況下,流量漂移增加了一個蒸發器整體換熱系數,使其壁面溫度降低,而另一個蒸發器壁面溫度不會發生明顯的變化.
圖12為大范圍質量流速下兩個蒸發器由于流量漂移而出現的系統壓降和流量分配特征,進口過冷度為20 ℃、加熱功率為500 W/cm2,其中紅色虛線表示不穩定狀態,藍色實線表示穩定狀態,并表示出上文討論的工況點.可以發現,在此參數下不穩定的工況點主要出現在兩個蒸發器都運行在負斜率區域時,這是由于蒸發器在此參數下的流動特征曲線負斜率區域較為平緩.將圖10中的沸騰起始流量Wa和負斜率區域消失流量Wb表示在圖12中,假設蒸發器2中的流量始終高于蒸發器1,其流量分布為帶有灰色陰影的曲線.在沸騰起始流量Wa之前,蒸發器內部呈現均勻分配狀態;流量繼續降低導致流量分配不均的情況出現,蒸發器2在 III 區隨總流量減小而增大到極值W2,max,蒸發器1在 II 區流量降低直到Wb.總流量繼續降低,蒸發器2在 II 區流量從極值W2, max降低到Wa,蒸發器1在 I 區流量降低到極小值W1,min.隨后,總流量的降低使得蒸發器1的流量增加、蒸發器2的流量減低,最后都達到負斜率消失點Wb,此后兩蒸發器間保持均勻分配的狀態.

圖12 兩并聯熱沉的流量分配
為了評價蒸發器間的流量分配,引入無量綱流量βi[32]和流量分配不均勻系數φ:
(21)
(22)
式中:n為工況點數目;βavg為n個工況點的平均無量綱流量.


圖13 不同加熱功率下的兩蒸發器流量漂移
當進口過冷度逐漸降低時,所有工況的不均勻系數φ隨之降低,如圖14所示.因為隨著進口過冷度的降低,蒸發器流動特性曲線中的負斜率區域斜率減小,從而降低了工況的不均勻程度.當進口過冷度達到15 ℃后,不穩定區間只會出現在均勻分配的區域,不穩定工況點的不均勻系數變為0.不均勻分配程度隨過冷度減小逐漸降低,直至完全消失.

圖14 不同過冷度下兩蒸發器的不均勻系數和不穩定區間占比
由于進出口的管道長度會影響蒸發器的流動特征曲線,所以蒸發器的布置對系統的流量分布有著比較大的影響.圖15表示在進口過冷度為20 ℃、加熱功率為500 W/cm2時,當總長度保持不變,采用流量分配不均勻系數φ表示蒸發器布置位置對系統流量分布的影響.當蒸發器之間的距離發生變化時,其流量分配情況也存在一定的差異.定義無量綱長度di=d/Ltot,表示蒸發器間距在總長度所占的比例.可以發現,在總長度為220 mm時,隨著di增大,所有工況點的不均勻系數呈現逐漸增大的趨勢,這是由于隨著蒸發器之間距離的增大,其流動特征曲線的差異逐漸增大,導致蒸發器的運行狀態更加偏離均分的情況,進口管道長的蒸發器內,會得到更多的流量.穩定工況點除di為0.636外,其余均呈現上升的趨勢,而不均勻系數的下降是因為不均勻分配工況中穩定點所占比例減小導致的.此外,隨著蒸發器之間距離的增大,其最大流量比也會增加,在這些工況點發生流量漂移后的危害也會增大.

圖15 不均勻系數和最大流量比與兩蒸發器間距離的關系
本文針對并聯微通道熱沉系統進行建模,研究在部分區域加熱情況下蒸發器熱沉系統中的流量分配問題,并分析了不同參數對單個蒸發器以及系統的影響情況.為了更好地表示流量分配對蒸發器系統的影響,對蒸發器的溫度分布進行了計算.主要結論如下:
蒸發器系統中的靜態不穩定性與其流動特征曲線的負斜率區域相關.進口過冷度的增加、加熱功率增加以及飽和溫度的降低都會增加流動特征曲線的負斜率傾斜程度,引起系統的不穩定;通道尺寸的降低可以增加進口壓降,從而降低負斜率區域傾斜程度,但會增加系統整體壓降.
蒸發器的結構參數對不穩定性也會有較大的影響.在總長度保持一定時,增加蒸發器進口連接管道所占比例能夠提高系統的穩定性,甚至可以消除負斜率區域;增加出口連接管道所占比例會加劇系統的不穩定性.合理安排蒸發器的位置有利于提高系統的穩定性.
當蒸發器系統運行在不穩定區間時,受到擾動便會產生流量漂移,從而改變蒸發器間的流量分配情況.通過表示出不同流量下蒸發器溫度的分布情況,可以發現,在部分加熱的條件下,產生流量漂移不一定會惡化系統的對流換熱.流量降低使得部分兩相區出現在加熱區域內,換熱系數的增大降低了蒸發器壁面的最高溫度.而得到更多流量的蒸發器由于全部是單相換熱,其平均溫度會大于初始狀態.
兩個蒸發器加熱功率存在差異時,低加熱功率的蒸發器會得到更多的流量,隨著加熱功率差異的增大,其不均勻程度也會增加;進口過冷度的減小會降低系統整體的不均勻分配程度,由于兩相區長度的增加,其不穩定區間占比先會有所上升,最后負斜率區域消失,不穩定區間也消失;不同的進出口管道長度配置會影響系統流量的分配傾向,進口長度占比大的通道會得到更多的流量,隨著蒸發器間的無量綱間距增大,其不均勻程度也會增加,發生流量漂移后的危害也會增加.