單環路液氫溫區脈動熱管高充液率工況計算流體動力學（CFD）模擬

徐金柱，焦波，孫瀟，王芳，甘智華

（1 哈爾濱理工大學榮成校區機械工程系，山東榮成264316；2 浙江大學制冷與低溫研究所，浙江杭州310027）

隨著電子芯片的小型化，其局部熱通量已經高達102～103W/cm2[1]，嚴重威脅到設備的安全運行。為了解決這一問題，亟需開發一種高效的傳熱元件。脈動熱管（pulsating heat pipe，PHP）具有體積小、結構簡單、傳熱效率高等優點，一經提出就引起了研究者們的注意。它由一根毛細管彎折而成，充注完成后由于毛細力的作用工質在管內形成隨機分布的氣液塞，工作時，氣塞會由于蒸發和冷凝的作用而膨脹或縮小，推動工質在管內形成振蕩或循環流動。與傳統熱管相比，PHP啟動更快，傳熱溫差更小，在低溫超導[2]、空間應用[3]、余熱回收[4]等領域具有廣闊的發展潛力。但PHP內部為復雜的氣液兩相流，目前仍無法完全理解其運行機理[5]，并且影響PHP傳熱性能的因素較多，增加了其研究難度。

工質對PHP 的傳熱性能有很大影響，它不僅決定了PHP 的管道臨界直徑，而且不同工質PHP的換熱能力有很大不同[6]，工質的選擇取決于PHP的運行溫度和應用范圍。徐冬等[7]總結了近些年來低溫脈動熱管的實驗及理論研究進展，他們指出低溫溫區PHP 的研究起步較晚，影響傳熱性能的參數尚未被完全確定。隨著低溫超導材料（臨界溫度39K）的發現和應用，液氫溫區PHP的高效傳熱引起了研究人員的關注。2011 年，Natsume等[9]對氮、氖、氫低溫PHP 的傳熱性能進行了測試，結果表明氫工質PHP的有效熱導率高達500～3500W/(m·K)。2014 年，甘智華等[10]設計了針對液氫溫區PHP的實驗臺，以探索多種參數（彎頭數、加熱功率、充液率、傾角等）對PHP 傳熱性能的影響。他們[11]測試了充液率34.2%情況下脈動熱管的運行過程，隨后，變充液率和變加熱功率下PHP的性能得到了實驗研究[12-13]。2019年，孫瀟和甘智華等[14]在實驗中將脈動熱管的冷凝段溫度分別維持在不同數值，分析了不同的冷卻溫度對液氫溫區PHP傳熱性能的影響，發現冷凝段溫度的提高可以增加PHP 的傳熱效率，并不會對PHP 的啟動造成影響。

另一方面，研究PHP 傳熱機理的手段也得到了極大發展。隨著計算機內存的擴大，計算速度的加快及算法的優化， 計算流體動力學（computational fluid dynamics，CFD）數值模擬被大量運用于分析PHP的流動及傳熱特性。劉向東等[15]建立了3D 數值模型分析了板式脈動熱管內工質的流動方式，他們指出工質在板式PHP 存在三種流動方式：離散的氣泡流，氣塞流和長氣塞流。Pouryoussefi和Zhang[16]建立了一個2D模型對PHP進行了分析，探究了以水作為工質的脈動熱管的流動情況。他們指出，對于閉式PHP，當蒸發段的熱通量較高時，從蒸發段向冷凝段運動的液塞具有較大的速度，這個速度足夠大以至于液塞可以穿越冷凝段到達下一個蒸發段。此時，PHP內的流動即從振蕩流變為了循環流。Wang等[17]建立了2D模型對充液率分別為30%、50%和70%的PHP 進行了數值模擬，他們得出，在底部加熱模式下，隨著加熱功率的增加，不同充液率的PHP傳熱熱阻均會減小，而充液率較低的PHP具有更小的熱阻。隨后他們[18]對單環路水工質PHP 的傳熱性能進行了CFD 模擬,利用2D 模型，分析了以水為工質的PHP 在加熱功率變化時熱阻的大小，并將模擬結果與Saha等[19]所做的實驗相對比。他們指出，模擬與實驗得到的熱阻變化趨勢一致，并且它們之間的誤差為10.01%。在低溫PHP 的CFD 模擬中，陳曦等[20]在對氮工質PHP 進行模擬時發現，在PHP 穩定運行后，管內工質的循環流動并不是一直在進行，而是在運行一段時間后會經歷一段時間的停滯，然后重新運行。

綜上所述，液氫溫區脈動熱管的實驗研究已經積累了一定的研究成果，CFD數值模擬也被廣泛地應用到PHP 理論分析中，同時有效性得到驗證，但尚未發現針對液氫溫區PHP 的CFD 模擬研究，啟動階段及運行階段脈動熱管管內的壓力變化，以及壓力變化與脈動熱管流動方式的關系，在以往的CFD 數值模擬中均未發現相關報道。由于氫的密度、汽化潛熱遠小于常溫工質，飽和壓力對溫度的導數dpsat/dT高于常溫工質一個數量級[14]，導致管內壓力隨溫度的變化更加劇烈，在傳熱性能及運行機理也與常溫脈動熱管有所差異，因此開展低溫液氫溫區脈動熱管的研究對于全面理解脈動熱管傳熱機理與流動規律具有重大意義。由于低溫工況下實驗難度大、耗時長，可視化實驗尤其困難，利用CFD模擬可以獲得PHP 的傳熱性能和內部流動特性。基于此，在本文作者課題組已有的研究基礎之上[5,10-13,21]，本文針對單環路液氫溫區PHP建立了二維數值模型，對80%充液率工況進行了模擬，并與實驗對比驗證了模型的有效性，分析了在啟動過程中PHP 內工質的流動與管內的壓力、溫度的關系，探索了穩定運行階段PHP 中溫度及壓力的振蕩規律，并討論了蒸發段的加熱量對PHP 傳熱性能的影響，旨在發展液氫溫區PHP理論分析模型，為指導工程應用提供依據。

1 控制方程

在本研究中，采用VOF 方法追蹤PHP 中氣液界面。VOF方法[22]適用于兩種或兩種以上不相溶的流體。因此，根據計算單元中體積分數的大小可以將控制體分為兩種情況：一種是純流體，αk=1；另一種是混合流體，αk≠1。在PHP 中只存在氣液兩相，當控制體中只存在氣相或液相時，控制方程與單相流的控制方程相同。當控制體中存在相界面時，控制體中氣相與液相的體積分數之和為1，如式(1)所示。

流動及傳熱的控制方程見式(2)～式(5)[23]。

連續性方程

動量守恒方程

能量守恒方程

動量方程和能量方程所求解出的速度場和溫度場被氣相和液相流體共享[24]，因此混合相的物性被用于動量方程與能量方程中，見式(6)。

在PHP 運行過程中，發生在氣液界面上的傳熱和傳質是引起界面變化根本原因[25]。為了計算相變過程中傳質量和傳熱量，Lee 等[26]提出了一個傳質模型，開發了一套完整的CFD 代碼，并將其連接到控制方程中。確定了傳質量后，相變過程中的傳熱量由傳質量與汽化潛熱的乘積而計算[式(7)～式(10)]。

2 物理模型和邊界條件

PHP管徑必須足夠小，才能使工質在管道中呈氣液塞分布，此時，脈動熱管的臨界直徑Dcrit需滿足以下條件[式(11)][27]。

密度和表面張力是關于溫度的函數，那么PHP的臨界直徑同樣依賴于工質的運行溫度，圖1表示的是氫工質脈動熱管的臨界直徑與運行溫度的關系，在實驗中取脈動熱管的內徑為2.3mm，此時對應的最高工作溫度為28.6K。當運行溫度低于此溫度時，該值均小于臨界直徑，滿足PHP運行條件。

圖1 臨界直徑與運行溫度的關系

由式(11)可知，重力和表面張力在內徑的選擇上起到了重要作用，由于本文所建立的是二維模型，相當于將物理模型沿紙面拉伸1m[24]。為了比較兩者的差別，如表1所示，對二維模型與三維模型的彎液面進行了分析，假設微元體高度為ΔL。

從以上量綱分析可知，在相同表面張力的作用下，二維模型的內徑為三維模型內徑的一半。因此在模擬中內徑為實驗內徑的二分之一至1.15mm，汪健生等[28]在對水工質脈動熱管模擬時采用了相同的內徑處理方式。

表1 二維模型與三維模型的區別

由以上分析，本文建立了單環路PHP的二維數值模型，如圖2(a)所示，垂直放置并采用底部加熱模式。蒸發段、絕熱段和冷凝段的長度分別為51.5 mm、100mm和51.5mm。工質為氫，充液率為80%。

圖2 物理模型

在本研究中，作出了以下假設：①忽略壁面厚度；②氣體密度遵循理想氣體狀態方程，其余物性為溫度的函數；③飽和溫度為壓力的函數；④流動為層流。為了保證網格的質量和改善計算精度，在模擬中采用ICEM 生成結構化網格。圖2(b)展示了部分蒸發段的網格示意圖，近壁面處的網格被加密。

根據實際的物理過程，整個模擬過程分為3個階段。

（1）初始階段 實際過程中工質充注到PHP后，由于毛細作用，工質在PHP 內會隨機形成交替分布的氣塞和液塞。在模擬時，為了得到這種氣液塞隨機分布狀態，初始化后，將所有壁面邊界條件設置為冷凝段溫度（19K）。初始時，管內工質的溫度也與此相同，因此在初始化時，將管內壓力設置為初始溫度對應狀態下的飽和壓力。

（2）啟動階段 管內形成相對穩定的氣液塞分布后，蒸發段、絕熱段和冷凝段的邊界條件分別為：恒熱通量、絕熱和恒溫。維持冷凝段的溫度為19K，蒸發段的熱通量為加熱量與蒸發段的表面積之比。

（3）穩定階段 當蒸發段溫度在小范圍周期振蕩即認為PHP 已經達到穩定階段。模擬中在啟動成功并且溫度達到穩定后，逐漸提高蒸發段的加熱量，以分析其對PHP傳熱性能的影響。

本模擬中，所有控制方程由基于有限體積法的軟件Fluent 求解。PISO（pressure implicit splitting of operator）算法用于迭代離散方程，壓力的差值格式為PRESTO!格式，動量和能量方程的迭代采用二階迎風差分格式。模擬中測試結果顯示采用默認的松弛因子就可以獲得收斂的結果。用于計算的時間步長為10-4s，收斂標準分別是：連續性方程5×10-4，速度10-4，能量方程10-7。

3 結果與討論

3.1 模型的有效性驗證

為了檢驗網格獨立性，驗證數值結果的可靠性，采用了4種不同的網格，并計算得出了相應的網格數對應的熱阻，如圖3所示。考慮數值模擬的計算精度，在本研究中網格數量取為58772。

圖3 網格獨立性驗證（Q=0.5W）

為了驗證數值模型的有效性，將模擬結果與文獻[14]中的實驗結果進行了對比，實驗中采用了彎頭數N=2的脈動熱管，在本文中對單環路脈動熱管進行了數值模擬，為了比較兩者的傳熱性能，將加熱功率減小到文獻中的1/2，以此保證每個彎頭的熱通量相同。PHP的傳熱性能用熱阻來說明，總熱阻的表達如式(12)所示。

圖4給出了模擬與實驗熱阻的對比，可以看出兩者具有良好的一致性，最大的熱阻誤差不超過15%，這是因為在模擬過程中會不可避免地由于幾何建模、假設條件等引入數值誤差，并且模擬與實驗的邊界條件不可能完全相同，導致模擬與計算存在一定誤差。從圖4中可以看出，在充液率為80%的實驗工況下，隨著加熱功率的增加，脈動熱管的熱阻逐漸減小，但熱阻變化率隨著加熱功率的增加而降低，CFD數值模擬得到了與之相接近的變化趨勢。這是因為當加熱功率增加后，蒸發段的工質更容易發生核態沸騰，液膜蒸發速度加快，相變速率提高，導致脈動熱管內工質的振蕩頻率提高，從而提高了傳熱性能降低了熱阻。隨著加熱功率的增加，熱阻逐漸趨于穩定。這是由于脈動熱管流動需要克服的流動阻力，由于充液率較高，管內工質的流動更接近單相流，熱阻隨著加熱功率的變化更容易達到穩定。

圖4 實驗熱阻與模擬熱阻對比

3.2 初始階段PHP內壓力分布特點

當工質充注到PHP 后，由于毛細作用管內的工質會形成隨機分布的氣塞與液塞，圖5為初始階段管內氣體體積分數分布及壓力分布云圖。從整體圖上看，PHP內部的壓力存在一個梯度，高壓區位于蒸發段，最高壓力出現在蒸發段的底部，兩個通道沿著冷凝段方向壓力逐漸下降，這是由于重力在初始階段的氣液分布與壓力分布中起到了重要作用。在局部區域，壓力波動下降，高壓與低壓相互交錯，圖6展示了軸線方向上[P-Q是圖2(a)中蒸發段右側的直管段部分]壓力與對應位置的體積分數的變化。可以發現壓力的變化和氣液分布存在一定的相關性，高壓區域對應著氣塞位置，低壓區域對應液塞位置。氣塞和液塞間存在清晰的相界面，界面之間存在一個明顯的壓力跳躍[29]，這是由于表面張力和液膜彎液面引起的。初始階段的氣液分布與壓力分布為PHP的啟動提供了基礎。

圖5 初始階段氣體體積分數分布云圖及壓力云圖

圖6 初始階段部分蒸發段的壓力分布

3.3 啟動階段PHP的流動分析

脈動熱管的啟動與管內壓力的變化息息相關，在啟動階段，蒸發段與冷凝段依舊由于重力的影響而存在壓差，但由于液體的蒸發和冷凝使PHP 內壓力變化的量級遠大于此，因此可以忽略蒸發段和冷凝段的壓差，通過監測脈動熱管蒸發段內點6#流體的壓力[6#的位置如圖2(a)所示]，從而得到PHP管內壓力隨時間的變化曲線，如圖7 所示。同時，在模擬中得到了啟動階段管內氣體的體積分數變化，可以表明工質在不同時刻的流動狀態，由此分析壓力變化對工質流動的影響。在本研究中工質經歷了兩次循環（圖7中ABC和CDE）完成啟動。在每次循環中，工質都經歷了以下3個過程：通道內同時向上流動、單方向循環流動和逆向流動。

圖7 PHP內流體壓力波動曲線（Q=0.27W）

經過初始階段的計算，脈動熱管內形成了隨機分布的氣塞和液塞。從圖7中A點開始對PHP的蒸發段施加恒熱通量邊界條件，脈動熱管開始進入啟動階段。蒸發段的壁面被加熱，內部的流體發生核態沸騰，一些汽化核心隨機形成，產生小氣泡，原來處于蒸發段的氣塞也會隨著液膜蒸發而膨脹，由于PHP 蒸發段的加熱面是對稱布置的，蒸發段兩根通道內的氣塞會由于膨脹作用同時向冷凝段運動，圖8(a)展示了蒸發段右側通道P-Q 的氣塞變化。與此同時，在冷凝段中由于PHP 壓力升高導致飽和溫度升高，氣塞溫度低于飽和溫度，氣塞開始冷凝，逐漸縮小為小氣泡，原來的小氣泡冷凝為液相，如圖8(b)所示。此階段工質在通道內同時向上流動。

隨著蒸發段氣塞液膜的不斷蒸發，氣泡合并，蒸發段中氣塞逐漸變長，脈動熱管內的壓力也不斷增加。2.4s時脈動熱管內的壓力與氣體體積分數分布如圖9(a)所示，可以發現蒸發段左右兩側通道內的壓力明顯不同，左側通道內的壓力低于右側通道的壓力。這是因為汽化核心的產生位置與蒸發作用都具有一定的隨機性，蒸發段的氣塞數量、大小因此會產生差異，氣液分布狀態發生了變化，PHP蒸發段左右兩側通道出現壓差。當壓差足夠克服重力和剪切力時，管內工質開始向同一方向的流動，如圖9(b)所示，在此工況中，PHP內的工質最開始沿順時針方向流動。而此時PHP 整體的壓力仍在增加，直到3.3s附近，氣塞即將運動到冷凝段，壓力達到最高值（圖7中B點）。如圖9(b)所示，此階段工質的流動方式為單方向循環流動。

PHP的壓力最高值對應飽和溫度的最高值，蒸發段的工質需要達到更高的溫度才能發生核態沸騰與蒸發，而此時工質剛剛從冷凝段到達蒸發段，需要較長的加熱時間，因此在3.5～4.0s 期間蒸發段幾乎全部為液相，無法為工質運動提供必要的壓力。由于慣性作用，工質會繼續向前流動一段時間后轉變流動方向，如圖10 所示。隨著氣塞在運動中冷凝使管內壓力降低，飽和溫度降低，蒸發段的工質重新開始發生核態沸騰，產生小氣泡，壓力開始回升（圖7 中C 點）。由于工質的逆向流動，使蒸發段右側通道的工質在蒸發段的加熱時間比左側通道的工質加熱時間長，因此它會比左側通道更早發生核態沸騰，產生更多的氣泡，如圖10 中4.2s所示。此時，工質的流動方式為逆向流動。

隨著加熱的進行，蒸發段左右兩側出現核態沸騰現象，同時冷凝段的氣塞也相繼液化，如圖11所示。PHP兩根通道中的工質重新開始同時向上流動。由于右側通道比左側通道具有更多氣塞，工質比上一次更早形成單方向循環流，氣塞更快地到達冷凝段，使得管內壓力的峰值（圖7中D點）低于第一次的峰值（圖7 中B 點），并在隨后工質改變流動方向，完成第二次循環。隨著循環的進行，壓力的振蕩趨于穩定，PHP進入穩定階段。

圖8 氣塞的變化

圖9 工質沿順時針流動開始時壓力與氣體體積分數云圖及體積分數隨時間的變化

圖10 PHP氣體體積分數變化及流動方向轉變

3.4 穩定階段PHP的流動及傳熱特性

在穩定運行過程中，脈動熱管中的流動主要為順時針循環流，這種流動方式使通道出現上升管和下降管，工質順時針流動時，左側通道為上升管而右側通道為下降管。工質的溫度與氣體體積分數分布如圖12 所示，可以發現上升管中流體的溫度始終高于下降管中流體的溫度。這是因為下降管中的流體剛經過冷凝段的冷卻，流體溫度較低；而上升管中的流體在蒸發段經過了充分加熱，在上升管發生核態沸騰，產生氣泡，氣泡相互合并形成氣塞，將工質分隔為氣液塞分布狀態，到達冷凝段后冷凝。

圖11 PHP內氣體體積分數變化（4.3～5.9s）

工質的溫度變化是判斷PHP 是否成功啟動的重要指標。PHP 內的溫度變化如圖13 所示，溫度振蕩的峰值均已在圖中標出。

其中，Twe為蒸發段壁面的平均溫度；Tev和Tcv分別為蒸發段和冷凝段流體的平均溫度，由式(13)和式(14)計算[1#為蒸發段軸線上的點，位置如圖2(a)所示]。

圖12 穩定階段PHP內工質的溫度分布及體積分數分布

圖13 PHP溫度隨時間的變化（Q=0.27W）

與圖7壓力變化曲線對比可以發現，蒸發段壁面溫度的變化與壓力的變化趨勢相同，但在啟動階段溫度的峰值所對應的時間（2.8s 和5.3s）明顯早于壓力峰值所對應的時間（3.3s 和5.8s）。開始時，PHP 的溫度振蕩較大，主要是因為剛開始PHP 中工質運動速度較小，蒸發段的氣塞持續蒸發，部分液膜被蒸干，氣塞直接接觸壁面，導致氣塞及蒸發段的壁面溫度急劇升高。而當工質開始循環流動后（2.4s后），絕熱段的流體進入蒸發段，使蒸發段的壁面溫度開始急劇下降，隨著工質的不斷流動，壁面溫度也開始趨于穩定。值得一提的是，蒸發段流體的平均溫度則可能會低于冷凝段流體的平均溫度。這是由于流體速度較高，從冷凝段流向蒸發段的冷流體在蒸發段尚未被充分加熱的原因。

3.5 加熱功率對PHP蒸發段溫度波動的影響

在模擬中改變蒸發段的熱通量，使加熱量逐漸從0.27W 增加到1W，得到了PHP 溫度隨加熱量的變化曲線，如圖14 所示。在0.27W 加熱量下，蒸發段壁面溫度有規律地振蕩，當加熱量升高到0.5W 時，壁面溫度小范圍升高，隨后穩定，而在0.5W 和0.75W 加熱功率下，PHP 的壁面溫度振幅很小，在功率增加到1W時，振幅突然增加，而當加熱量增加到1.25W 后，壁面溫度振幅逐漸發散，溫差增加，傳熱逐漸惡化。

圖14 脈動熱管溫度隨時間的變化曲線

從圖14 中可以發現，在穩定運行階段，溫度的振蕩有周期性，對不同功率下穩定運行階段的溫度進行快速傅里葉變換（FFT），得到的信號強度與頻率的關系如圖15 所示。從圖中可以發現各個功率下FFT分析均存在一個顯著的峰值，對應的頻率如圖中所示。這說明PHP 在穩定運行階段存在顯著的周期性，運行頻率隨著加熱功率的升高而升高，而當加熱功率達到一定值時，運行頻率趨于穩定。通過與熱阻的變化（圖4）進行對比，可以發現，當加熱功率較低時，溫度振蕩的周期長，頻率小，熱阻大，隨著加熱功率的升高，溫度振蕩頻率加快并趨于穩定，PHP的傳熱熱阻隨著振蕩頻率的加快而減小，并最終趨于穩定。

4 結論

本文對單環路液氫溫區PHP 建立了2D 模型，利用CFD 數值模擬研究方法，分析了其在80%充液率下的流動及傳熱特性，研究了不同階段PHP內壓力分布與溫度變化。結果表明，CFD數值模擬結果與實驗結果誤差小于15%，脈動熱管的熱阻隨著加熱功率的增加先減小而后不變。初始階段PHP內氣塞與液塞交替形成，由于重力影響，PHP內壓力存在明顯的梯度，局部位置由于氣液黏性的不同，導致壓力呈高低交錯分布。當有熱量輸入后，工質經歷了兩次通道內同時向上流動、單向循環流和逆向流動的過程完成了啟動進入穩定階段，此時PHP的主要流動形式是順時針循環流，蒸發段壁面溫度的變化與壓力的變化趨勢相同，在穩定階段呈有規律的周期波動。不同加熱功率下，溫度振蕩的頻率不同，它隨著加熱功率的增加先增大而后趨于穩定。

圖15 不同加熱功率下穩定運行階段蒸發段壁面溫度的FFT分析（FR=80%）

符號說明

d—— 直徑，m

E—— 比焓，kJ/kg

F—— 表面張力引起的體積力源項，N/m2

Fg—— 重力，N

Fs—— 表面張力，N

g—— 當地的重力加速度，m/s2

k—— 熱導率，W/(m·K)

p—— 壓力，Pa

Q—— 輸入功率，W

R—— 熱阻，K/W

S—— 源項，kg/(m3·s)或kJ/(m3·s)

T—— 溫度，K

α—— 體積分數

β—— 時間松弛因子，s-1

ρ—— 密度，kg/m3

μ—— 黏度，Pa·s

σ—— 表面張力系數，N/m

下角標

cv—— 冷凝段流體

ev—— 蒸發段流體

E—— 能量

k—— 第k相

l—— 液相

M—— 質量

sat—— 飽和狀態

v—— 氣相

wc—— 冷凝段壁面

we—— 蒸發段壁面