黃志豪，李光熙，唐桂華，李小龍，范元鴻

（1 西安交通大學熱流工程與科學教育部重點實驗室，陜西西安 710049； 2 西安航天動力研究所，陜西西安 710100）

引 言

再生冷卻一般利用燃料的熱沉冷卻燃燒室壁面，然而隨著飛行速度的增加，高超聲速飛行器燃燒室壁面的散熱面臨極大挑戰(zhàn)。常規(guī)的解決方法是通過熱管理提高燃料熱沉的利用效率[1]，然而僅依靠提高燃料熱沉的利用效率無法滿足冷卻要求，需要增加冷卻劑以輔助燃料冷卻燃燒室。常用冷卻劑有正癸烷、航空煤油RP?3、水和其他高密度烴。與其他冷卻劑相比，水具有成本低、來源廣和換熱能力強的特點，作為優(yōu)質冷卻劑廣泛應用于各種換熱器。由于燃燒室壁溫高，常壓下的水容易出現(xiàn)相變，而超臨界水的密度近似液體的密度量級，同時又具有與氣體相近的擴散能力，具有非常好的流動性和傳輸性能，超臨界水十分適合作為飛行器燃燒室壁面冷卻劑。

超臨界水在核燃料堆冷卻換熱研究較為深入，Shang 等[2]數(shù)值模擬研究了直徑對超臨界水在水平圓管中流動的影響，發(fā)現(xiàn)直徑影響傳熱是由于邊界層浮力的變化，小管徑更有利于傳熱。范辰浩等[3]實驗研究了小管徑內超臨界水的傳熱惡化特性，發(fā)現(xiàn)在高熱通量時傳熱惡化的產生是由于流動層流化的發(fā)展，而管內徑向的工質溫度和物性的巨大差異促進了流動層流化的發(fā)展。Wang等[4]基于大渦模擬研究了圓形管道中超臨界水在無重力、豎直向上和豎直向下流動下的傳熱，發(fā)現(xiàn)與無重力流動相比，向下流動的傳熱略有改善，向上流動的傳熱出現(xiàn)惡化，提出由于浮力的影響而產生這些差異。

Bai 等[5]數(shù)值模擬了豎直向上圓形管道中超臨界水非均勻加熱對流動換熱的影響，發(fā)現(xiàn)與均勻加熱方式相比，非均勻加熱方式表現(xiàn)出顯著的差異，截面溫度分布沿周向呈現(xiàn)出較大的不均勻性，傳熱強化僅出現(xiàn)在局部區(qū)域。Gao 等[6]對超臨界水在非均勻加熱的水平圓管中的流動換熱進行了數(shù)值模擬，分析了超臨界水的非均勻傳熱特性及機理，并討論了側向二次流和浮升力效應對傳熱作用開始的判據，發(fā)現(xiàn)二次流與熱物性之間的相互影響對傳熱具有重要影響。許多文獻對超臨界水傳熱強化和惡化產生的機理進行了分析[7?11]，發(fā)現(xiàn)流體物性變化對超臨界水傳熱具有重要影響。文獻中對超臨界水換熱關聯(lián)式的適用性也進行了研究[12?17]，發(fā)現(xiàn)超臨界水的關聯(lián)式需要根據不同工況的換熱特點進行修正。

現(xiàn)有研究主要關注超臨界水在圓形通道內的流動傳熱，非圓結構如類矩形、類三角形和方環(huán)形等結構也存在相關研究[18?21]，而高超聲速飛行器燃燒室壁面的冷卻通道具有方形單側加熱特點，目前缺乏超臨界水在方形單側受熱通道內流動傳熱的相關研究，亟需獲得其流動傳熱規(guī)律。

本文對超臨界水在方形冷卻通道的流動換熱進行數(shù)值模擬研究，通過分析流場和超臨界水熱物性變化，闡明傳熱惡化產生和恢復的機理，并驗證超臨界流體常用傳熱關聯(lián)式在單側加熱方形通道內的適用性，最后，提出非對稱凹槽和雙通道強化結構，防止通道傳熱惡化，提升通道綜合換熱性能并進行場協(xié)同原理分析。

1 物理模型及數(shù)值模擬方法

1.1 計算模型

航空發(fā)動機冷卻通道緊貼高溫壁面，冷卻介質在多個微小通道內流動，考慮模型的周期性，研究的模型可以簡化為單側加熱的方通道，如圖1所示。方通道當量直徑為D，壁厚a=2 mm，方形通道加熱長度為Lt，為得到較均勻的速度分布并防止回流，模擬過程中增加長度Ld=Lu=0.1 m 的入口段和出口段。

圖1 方形通道示意圖Fig.1 Schematic of square channel

為提升通道換熱性能，提出了多種強化通道。圖2（a）～（c）分別表示對稱凹槽、倒角凹槽和非對稱倒角凹槽的結構。其中加熱長度Lt=0.5 m，凹槽半徑R=2 mm，槽深d=1 mm，凹槽間距為25 mm，凹槽數(shù)目為16?？紤]入口處換熱較強，無須強化，第一個凹槽布置在距離入口50 mm 處，其余凹槽均勻分布于內壁面靠近加熱面處。凹槽加倒角結構如圖2（b）、（c）所示，其中圓角半徑r=1 mm，r1=2 mm，r2=0.6 mm。非對稱倒角結構由曲線與6 個點（P0，P1，P2，P3，P4，P5）確定。此外，也研究了非對稱凹槽最低點偏下游以及偏上游對通道流動換熱的影響。

圖2 凹槽結構示意圖Fig.2 Schematic of grooves

雙通道結構的布置是以1 mm 厚度的壁面平均分割流體域，考慮兩種布置形式，即上下布置（下層靠近加熱端）與并排布置。二者均為對稱結構，且結構參數(shù)一致。

1.2 邊界條件與求解方法

采用質量流量入口、壓力出口與均勻熱流邊界條件，計算工況如表1 所示。基于商業(yè)軟件Fluent求解，速度與壓力的耦合采用SIMPLE 算法求解，控制方程中對流項的離散采用二階迎風格式，擴散項采用QUICK 格式。當連續(xù)性方程、動量方程和湍流方程殘差值小于1.0×10?5，能量方程的殘差值小于1.0×10?8，且不隨計算步數(shù)變化時認為結果收斂。超臨界水熱物性通過調用NIST數(shù)據庫獲得。

表1 計算工況參數(shù)Table 1 Conditions of numerical simulation

1.3 模型驗證

在靠近固體和流體區(qū)域的內壁附近采用精細網格，確保無量綱距離y+小于1，以滿足SSTk?ω模型的要求。經網格的獨立性考核，最終采用627 萬的網格數(shù)量進行計算。

由于目前缺乏對方形通道中超臨界水傳熱的實驗研究，首先對超臨界水在圓形管道中豎直向上流動的數(shù)值模擬結果與Yamagata等[22]的實驗結果進行對比，實驗段由AISI 316型不銹鋼管光管制成，內徑為7.5 mm，壁厚為2 mm。實驗工況為：壓力P=24.5 MPa，質量流 速G=1260 kg·m?2·s?1，熱 通量q=233 kW·m?2。如圖3 所示，計算結果與文獻中的數(shù)據吻合很好，最大誤差為0.59%，表明本計算模型能夠較好地處理超臨界水的流動傳熱現(xiàn)象。

圖3 實驗結果與數(shù)值模擬結果的比較Fig.3 Comparison between available experimental and present numerical results

2 結果分析與討論

2.1 傳熱惡化機理分析

本節(jié)對超臨界水在豎直向上流動時出現(xiàn)傳熱惡化的工況（case 1）進行了分析。由于研究對象為超臨界水的豎直向上流動，而此時流動方向與重力方向相反，二次流僅僅依靠流體密度差產生，影響十分微小，因此不考慮二次流動的影響。

Koshizuka 等[23]對傳熱惡化時的表面?zhèn)鳠嵯禂?shù)做出如下定義:

其中hDB為Dittus?Boelter 關聯(lián)式[24]計算的表面?zhèn)鳠嵯禂?shù)。

圖4 表示了數(shù)值模擬與Dittus?Boelter 關聯(lián)式計算的內壁面局部表面?zhèn)鳠嵯禂?shù)的比值隨流體平均焓的變化，可以發(fā)現(xiàn)傳熱惡化的起始點A、傳熱惡化點B 和傳熱惡化恢復點C，這是后續(xù)分析的三個基準點。

圖4 豎直向上與無重力流動內壁面h/hDB沿程變化Fig.4 h/hDB on the inner wall with and without gravity

需要注意的是，當流體的焓值過高時，會再次出現(xiàn)h/hDB＜0.3 的現(xiàn)象，這是由于流體溫度靠近擬臨界溫度，比定壓熱容增加劇烈，Prandtl 數(shù)過大，導致hDB偏高，因此前面ABC 三點才是傳熱惡化產生的區(qū)域。

圖5、圖6展示了傳熱惡化區(qū)域A、B 和C 點的流場及物性分布，分別從邊界層厚度與近壁區(qū)湍動能兩個方面來解釋傳熱惡化產生和恢復的原因。為防止局部高溫出現(xiàn)，主要關注單側加熱下的受熱面，故圖5、圖6 主要表示了加熱面單側的一半流場與物性分布。其中，A 到B 的過程被定義為傳熱惡化加劇的過程，由圖5（a）、（b）可以看出，從A到B點黏性底層的急劇升溫導致流體動力黏度增加，動力黏度增加導致邊界層中黏性底層增厚，如圖5（c）所示。由于黏性底層中傳熱方式以導熱為主，而在傳熱惡化區(qū)，與緩沖層與對數(shù)層相比，黏性底層內熱導率整體處于較低的水平[25]，綜合影響下增加了黏性底層整體熱阻，降低了邊界層流體的流動換熱能力。而對于B 到C 過程，黏性底層內流體動力黏度降低，邊界層厚度恢復到與A點相似的水平，從而使C 點換熱能力提升，傳熱惡化恢復到正常傳熱狀態(tài)。

圖5 邊界層厚度的影響Fig.5 Effect of boundary layer thickness

圖6（a）可以看到，相比于A、C兩點，B點處黏性底層與緩沖層內密度均處于較低水平，緩沖層內尤為明顯，因此從A到B點的過程中，低密度的流體在浮力作用下加速上升，緩沖層內出現(xiàn)明顯的熱加速現(xiàn)象，該現(xiàn)象使得B 點處緩沖層與對數(shù)層內的速度分布比較平緩，從而使得該處湍動能明顯降低，導致傳熱惡化發(fā)生，如圖6（b）、（c）所示。而從B 到C點過程中，對數(shù)層流體密度進一步減小，流體速度增加，形成M 形速度分布，增強了湍流的強度，導致湍動能增加，從而使傳熱惡化恢復到正常傳熱水平。

圖6 近壁區(qū)湍動能的影響Fig.6 Effect of near?wall turbulent kinetic energy

綜上所述，導致傳熱惡化的原因為以下兩點：（1）黏性底層內黏度增大使得邊界層局部增厚；（2）黏性底層與緩沖層內密度減小導致湍動能降低。上述耦合影響導致了傳熱惡化的出現(xiàn)。

2.2 常用關聯(lián)式的驗證

文獻中存在較多超臨界流體換熱的關聯(lián)式，主要是在傳統(tǒng)Dittus?Boelter 關聯(lián)式的形式上引入壁面物性參數(shù)修正，以提高關聯(lián)式精度。選擇了6 個常用于超臨界流體換熱的關聯(lián)式以驗證其對方通道單側加熱超臨界水的適用性。式（2）～式（7）分別來 自 于Gupta 等[26]、Mokry 等[27]、Jackson[28]、Bishop等[29]、Kim等[30]和Fan等[31]，關聯(lián)式中計算的物性是根據模擬得到的流體平均溫度與相應固體壁面溫度，由NIST數(shù)據庫查出對應的物性參數(shù)。

圖7展示了各個關聯(lián)式的驗證結果，綜合來看，Gupta、Kim 和Jackson 關聯(lián)式與計算數(shù)據誤差最大，大部分結果的驗證誤差均大于30%。Bishop 關聯(lián)式有部分點的驗證誤差大于30%，Mokry 關聯(lián)式在小熱通量和大質量流速時與計算數(shù)據吻合較好，驗證誤差在30%以內，但其他工況與計算數(shù)據的符合較差，驗證誤差均大于30%。對于Fan 關聯(lián)式，除高質量流速（case 9）和大管徑（case 10）的工況的部分數(shù)據外，絕大多數(shù)數(shù)據驗證結果誤差均在30%以內，預測精度較高。根據2.1 節(jié)的分析，傳熱惡化的發(fā)生主要由近壁區(qū)黏度和密度的變化導致的，而Fan關聯(lián)式引入了壁面密度和黏度的修正，因此該關聯(lián)式精度較高，推薦使用Fan 關聯(lián)式來預測單側加熱方形通道內超臨界水的換熱性能。

圖7 關聯(lián)式計算與數(shù)值模擬結果的比較Fig.7 Comparison between the correlations and present numerical results

2.3 不同結構強化換熱效果比較

為防止豎直向上流動工況出現(xiàn)傳熱惡化，本節(jié)對比了2.1 節(jié)所述多種強化換熱結構的強化換熱效果。圖8 為G=500 kg·m?2·s?1時不同強化換熱結構的加熱底面平均溫度隨加熱長度的變化?？梢钥闯?，對傳熱惡化工況，5 種強化換熱結構均可降低壁面溫度，其中，四種凹槽結構和并排的雙通道均可以避免傳熱惡化，顯著降低壁溫。并排雙通道、非對稱倒角下游凹槽、倒角凹槽、非對稱倒角上游凹槽和對稱凹槽，在光滑通道傳熱惡化點處的壁溫下降幅度分別為18.55%、16.12%、14.76%、14.65%和13.29%，強化換熱效果為：并排雙通道＞非對稱倒角下游凹槽＞倒角凹槽＞非對稱倒角上游凹槽＞對稱凹槽，上下雙通道強化換熱效果最差，不能避免傳熱惡化的產生。

圖8 加熱底面平均溫度隨加熱長度的變化Fig.8 Average temperature on the heated surface against heating length

為綜合評價強化換熱結構的換熱能力和阻力損失，引入等泵功條件下綜合換熱評價因子PEC。

圖9、圖10 分別表示了豎直向上流動不同強化換熱結構Nu/Nu0與PEC 隨質量流速的變化?？梢钥闯龇菍ΨQ倒角下游凹槽綜合換熱性能最優(yōu)，非對稱倒角下游凹槽在質量流速G=500～1000 kg·m?2·s?1時，PEC 有微小的下降，但維持在一個較高的水平，且Nu/Nu0比值維持較高水平，說明在該流量區(qū)間，隨著流量增加非對稱倒角下游凹槽對流體的擾動增加，相較于光滑管，強化傳熱效果能維持較高水平。對稱凹槽與倒角凹槽結構也具有一定的綜合強化換熱性能，并排雙通道僅僅在質量流速G=500 kg·m?2·s?1時有一定的強化效果，其他情況下PEC＜1，綜合換熱性能較差。除了上下通道，其他通道在G=500 kg·m?2·s?1時PEC 能達到最大數(shù)值，是因為在該流量下，強化結構對流體的擾動影響最大，Nu/Nu0達到最大值，此時強化傳熱效果最佳。上下雙通道的Nu/Nu0較小，強化換熱效果有限，其PEC 在所有質量流速下均小于1，綜合換熱性能不如光管。

圖9 不同強化換熱結構的Nu/Nu0Fig.9 Nu/Nu0 of different enhanced structures

圖10 不同強化換熱結構的PECFig.10 PEC of different enhanced structures

綜上所述，在q=1000 kW·m?2、G=200～1500 kg·m?2·s?1的工況下，雙通道結構換熱效果有限，不適用于超臨界水的強化換熱；凹槽結構均具有一定的強化換熱效果，其中非對稱倒角下游凹槽結構強化換熱效果最佳。

2.4 強化換熱機理分析

本節(jié)對凹槽結構強化換熱的機理進行詳細分析。圖11 為G=500 kg·m?2·s?1時四種凹槽結構附近切面的流線圖，可以看到凹槽區(qū)域產生明顯的旋渦，導致該區(qū)域的流動分離，從而破壞邊界層的發(fā)展，增加流體的擾動。流體在凹槽下游位置與加熱壁面再接觸，使得該處邊界層變薄，從而強化換熱。由于旋渦區(qū)域僅出現(xiàn)在凹槽結構內部，對主流流動影響較小，因此不會帶來較大的壓降。值得注意的是，凹槽內部的旋渦雖然增強了流體的擾動，但由于其速度較小，因此換熱能力較弱，導致凹槽內部區(qū)域換熱效果較差，強化的部位只發(fā)生在凹槽下游的再附著點附近。因此，在保證流動分離的前提下，應盡可能減小旋渦強度以降低阻力。對稱凹槽結構由于壁面結構的突變，產生了較強的流動分離，導致在凹槽區(qū)域出現(xiàn)了一大一小兩個旋渦，如圖11（a）所示。較大的旋渦可以使主流產生流動分離從而強化換熱，但較小的旋渦對換熱并沒有積極作用，反而由于增強了流體的擾動而增大該處的流動阻力。與對稱凹槽相比，倒角凹槽壁面附近流線更為順滑，消除了凹槽區(qū)域較小的旋渦，使得流阻降低。同時，在凹槽下游部位，再附著區(qū)域由尖角變?yōu)槠交^渡的圓角，相當于增加了強化區(qū)域，因此倒角凹槽比對稱凹槽具有更高的傳熱性能。非對稱倒角上游凹槽結構上游陡峭而下游平緩，產生了一大一小兩個旋渦，由于有利于強化傳熱的旋渦產生于下游，而下游的區(qū)域較為狹窄，所以非對稱倒角上游強化傳熱效果一般。非對稱倒角下游凹槽結構壁面下游陡峭而上游平緩，增強了流體的擾動，產生了一大三小的四個旋渦，大的旋渦對倒角凹槽結構影響范圍更大，強化換熱效果更為明顯，小的旋渦則增加了流動阻力，綜合來看大的旋渦影響范圍更廣，因此非對稱倒角凹槽結構強化換熱效果更好，推薦使用非對稱倒角凹槽結構以改善傳熱惡化。

圖11 凹槽附近剖面流線圖Fig.11 Streamlines around various grooves

為進一步揭示傳熱強化機理，引入Guo 等[32]提出的場協(xié)同原理。一方面是溫度梯度場與流場的協(xié)同關系，即溫度梯度與速度的協(xié)同角β，協(xié)同角β越小，對流傳熱能力越好；另一方面是速度矢量與速度梯度矢量的協(xié)同關系，即速度與速度梯度的協(xié)同角α，協(xié)同角α越大，流動阻力越小。

圖12、圖13 均為沿著流動的中心面剖面圖，分別表示G=500 kg·m?2·s?1時對稱凹槽、倒角凹槽和非對稱倒角凹槽在Z=0.25 m 附近的協(xié)同角β與α的分布?？梢钥闯觯瑓f(xié)同角β值較小的范圍：非對稱倒角下游凹槽＞倒角凹槽＞對稱凹槽＞非對稱倒角上游凹槽，說明非對稱倒角下游凹槽結構的傳熱能力最強。此外，對稱凹槽結構中，β值較小的區(qū)域主要集中在凹槽內部，但由于該區(qū)域流速較小，換熱能力較差，因此對主流換熱能力的強化效果有限。而倒角凹槽與非對稱凹槽結構中β值較小的區(qū)域主要集中在凹槽與主流區(qū)域交界處，該處的換熱強化有利于凹槽區(qū)域和主流之間的熱量交換，且非對稱凹槽β值較小的區(qū)域向主流區(qū)沿伸幅度大于倒角凹槽，因此，非對稱倒角凹槽結構更有利于傳熱的強化。

圖12 凹槽附近的協(xié)同角β云圖Fig.12 Distribution of synergy angle β around various grooves

圖13 凹槽附近的協(xié)同角α云圖Fig.13 Distribution of synergy angle α around various grooves

由圖13（a）可以看出，對稱凹槽上游和下游部位存在較大的低α區(qū)域，這是由于對稱凹槽結構尖角的存在，流體在流動分離點和再附著點附近產生比較強烈的流動方向轉變[圖11（a）]，使得該處流動阻力嚴重增加。對凹槽結構采用倒角后，可以使流動更平滑，避免較強的流動分離及流動方向轉變，從而顯著改善這兩個區(qū)域的協(xié)同性，使局部阻力顯著降低。采用倒角的結構中低α區(qū)域影響區(qū)域：非對稱倒角下游凹槽＞倒角凹槽＞非對稱倒角上游凹槽，說明非對稱倒角下游凹槽結構流動阻力最大，非對稱倒角下游凹槽結構由于下游壁面變陡，增強了流體的擾動，強化換熱的同時也增加了局部阻力，非對稱倒角上游凹槽雖然降低了局部阻力但強化換熱效果一般，綜合PEC 考慮，非對稱倒角下游凹槽結構強化換熱效果最佳。

3 結 論

建立了超臨界水在單側加熱方通道內流動傳熱計算模型，基于數(shù)值模擬研究，分析了超臨界流體常用換熱關聯(lián)式對超臨界水在單側加熱方管中的適用性，計算結果可進一步為發(fā)動機再生冷卻通道的設計提供參考。同時詳細分析了傳熱惡化產生的機理，比較不同強化換熱結構（雙通道和凹槽）的優(yōu)劣及機理分析。主要結論如下。

（1）引起單側加熱方管內超臨界水流動傳熱惡化主要有以下兩點原因：①黏性底層內黏度增大使得邊界層局部增厚；②黏性底層與緩沖層內密度減小導致湍動能降低。上述耦合影響導致傳熱惡化的出現(xiàn)。

（2）綜合比較了六種關聯(lián)式，發(fā)現(xiàn)Fan 關聯(lián)式與計算結果符合程度最高。傳熱惡化的發(fā)生主要由近壁區(qū)黏度和密度的變化導致，而Fan 關聯(lián)式引入了壁面密度和黏度的修正，因此推薦采用Fan 關聯(lián)式預測單側加熱方形通道內超臨界水的換熱性能。

（3）在傳熱惡化產生的區(qū)域，綜合對比了雙通道和凹槽結構的強化傳熱效果，發(fā)現(xiàn)雙通道結構強化傳熱效果有限，不適用于超臨界水的強化傳熱；凹槽結構具有較好的強化傳熱效果，且引起的壓力損失小，其中非對稱倒角下游凹槽結構強化傳熱效果最佳。

（4）場協(xié)同原理分析發(fā)現(xiàn)，對于非對稱倒角下游的凹槽結構，其β值較小的區(qū)域占比大，且低α區(qū)域增加不明顯，說明其強化傳熱的同時，流動阻力增加不顯著，因此非對稱倒角下游凹槽結構綜合傳熱性能最佳。

符 號 說 明

cp——流體比定壓熱容，J·kg?1·K?1

f,f0——分別為強化結構沿程阻力系數(shù)及光滑通道沿程阻力系數(shù)

G——超臨界水質量流速，kg·m?2·s?1

H——流體焓，kJ·kg?1

h——方通道內壁面表面?zhèn)鳠嵯禂?shù)，W·m?2·K?1

k——湍動能，m2·s?2

Nub,Nu,Nu0——分別為Nusselt 數(shù)、強化結構Nusselt 數(shù)、光滑通道Nusselt數(shù)

P——超臨界水壓力，Pa

Pr——Prandtl數(shù)

qw——方通道加熱側外壁熱通量，kW·m?2

Re——Reynolds數(shù)

Tin——超臨界水入口溫度，K

v——軸向速度，m·s?1

Z——流體流過加熱段的長度，m

α——速度與速度梯度的場協(xié)同角，(°)

β——溫度梯度與速度的場協(xié)同角，(°)

λ——流體熱導率，W·m?1·K?1

μ——流體動力黏度，Pa·s

ρ——流體密度，kg·m?3

下角標

b——流體

w——壁面