夏 斌，梁新剛，徐向華，*

(1. 清華大學 航天航空學院，熱科學與動力工程教育部重點實驗室，北京 100084；2. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000)

0 引 言

英國反應發(fā)動機公司研發(fā)的協(xié)同吸氣式火箭發(fā)動機（synergistic air-breathing rocket engine, SABRE）作為一種新型預冷組合式發(fā)動機的典型代表，被認為是未來20 年內(nèi)最有可能實現(xiàn)單級入軌的動力系統(tǒng)[1]。由于SABRE 采用了先進的進氣預冷卻技術，提高了發(fā)動機進氣量、增加了壓縮機的壓比和燃油加入量，其吸氣模態(tài)的工作馬赫數(shù)上限可提高到馬赫數(shù)5[2]，能夠順利銜吸氣模態(tài)和火箭模態(tài)，從而避免了動力模態(tài)銜接的推力陷阱問題。為了保證SABRE 在高馬赫數(shù)下的進氣量和推力，需要將來流冷卻至最低-130℃[2]。然而，發(fā)動機吸氣模態(tài)工作包線內(nèi)12 km 左右的空域中含有大量水蒸氣，如不采取有效的抑霜措施，這些水蒸氣將在低溫的預冷換熱器表面結(jié)霜，并在數(shù)秒內(nèi)將換熱器堵塞[3]。因此，預冷換熱器的抑霜技術是SABRE 的核心技術之一。然而，目前這種預冷器抑霜技術的具體技術途徑未見諸報道。

雖然目前已有大量結(jié)霜抑霜的研究，但這些研究絕大多數(shù)是針對自然對流或低速流動條件且冷面溫度高于-30℃的結(jié)霜[4-5]，而SABRE 的進氣預冷器工作在來流速度很快（10 m/s 以上[1]）、冷面溫度很低的環(huán)境。由于快速來流空氣的作用，所結(jié)霜層較致密，不同于在自然對流或低速流動下所形成的晶枝結(jié)構(gòu)和稀疏霜層[6]。這種預冷器的結(jié)霜也不同于飛機結(jié)冰，飛機結(jié)冰是液態(tài)過冷水滴撞擊機翼表面凝固為冰的過程[7-8]，而預冷器結(jié)霜是空氣中的水蒸氣在低溫表面直接凝華為霜的過程。

發(fā)展類SABRE 的組合式發(fā)動機，進氣預冷器的結(jié)霜機理和抑霜技術是必須解決的難題，其中，強對流條件下低溫表面的結(jié)霜特點和規(guī)律又是其基礎。本文中將大于流速10 m/s 的流動稱為強對流條件，以區(qū)別于目前大多數(shù)低速流動條件下結(jié)霜研究中的強迫對流。目前，低溫表面上結(jié)霜的模擬方法可分為三類：第一類是使用基于實驗經(jīng)驗關系式的一維模型[9-11]，但該類方法是單向的流動影響結(jié)霜，無法求解與霜層動態(tài)變化對流場和傳熱的影響，并且無法求解復雜外形上的結(jié)霜，主要用于平板上的結(jié)霜模擬；第二類是同時分別求解濕空氣中的流場、傳熱傳質(zhì)，以及霜層中的傳熱傳質(zhì)，濕空氣與霜層之間的耦合作用通過界面上的傳熱傳質(zhì)交換實現(xiàn)[12-15]；第三類是使用多相流方法，將霜層看作流體，使用一套方程組同時求解空氣和霜層中流場和傳熱傳質(zhì)[16-20]。以上三類模型中，多相流方法不需要處理網(wǎng)格移動和重構(gòu)問題，因此適合模擬流動、傳熱、傳質(zhì)強耦合條件下的霜層生長問題。然而，已有的結(jié)霜模擬仍集中在自然對流或低速流動條件下的結(jié)霜。

本文發(fā)展了一種新的多相流模擬方法—p-VOF 方法（pseudo-VOF），針對作為預冷換熱器基本單元的圓管，開展了強對流條件下的常物性霜層二維圓管動態(tài)結(jié)霜模擬研究，分析了強對流條件下低溫圓管結(jié)霜的大致規(guī)律，為預冷換熱器結(jié)霜特性研究提供了基礎方法。

1 動態(tài)結(jié)霜模擬算法

p-VOF 方法是一種多相流模擬偽VOF 方法。該方法是基于VOF 方法的整體架構(gòu)，求解簡化的質(zhì)量守恒方程，替代求解CFD 軟件中原有的體積分數(shù)方程。由于不解決原有的體積分數(shù)方程，所以稱這種方法為偽VOF 方法。

模擬進行了以下簡化和假設：1）濕空氣是具有常物性的不可壓縮牛頓流體；2）當霜面溫度低于冰點，且濕空氣中水蒸氣濃度大于霜面飽和水蒸氣濃度時，水蒸氣在霜面凝華為霜，反之亦然；3）由于強對流條件下形成的霜層致密，可忽略霜層內(nèi)部的傳質(zhì)；4）相變只在水蒸氣和霜層（冷表面）之間的接觸面上發(fā)生。

1.1 控制方程

p-VOF 方法采用兩相流模型，其中濕空氣是主相，霜是第二相，濕空氣包含干空氣和水蒸氣兩種組分。

1.1.1 質(zhì)量守恒方程

原VOF 方法中，濕空氣相和霜相兩相使用相同的速度，其質(zhì)量守恒方程（體積分數(shù)方程）如下：

式中：α為體積分數(shù)，%；ρ為密度，kg/m3；u為速度矢量，m/s；t為時間，s；下標i代表a 或f，分別表示濕空氣相或霜相；兩相的體積分數(shù)之和為1。Sm是相變引起的質(zhì)量源相或質(zhì)量匯項，kg/(m3·s)，當i取f 時，Sm即為Smf。

由于原VOF 方法將所有物相均當作流體處理，所以每種物相均有對流項。但由于實際中的霜相是不會移動的，加之濕空氣為不可壓縮的假設，所以對于結(jié)霜模擬，質(zhì)量守恒方程（1）中的對流項（方程左邊第二項）可以忽略。因此，方程（1）可以簡化為：

式中：p為壓力，Pa；μ為混合物的黏性，N·s/m2；F是體積力或動量源相，N/m3。所有的變量都定義為兩相混合，并且能夠通過各相的體積分數(shù)權重計算出相關參數(shù)。

1.1.3 能量守恒方程

能量守恒方程也由濕空氣相和霜相兩相共用：

式中：E是兩相混合物的內(nèi)能，J/kg；keff是有效熱導率，W/(m·K)；Sh是相變導致的放熱或吸熱，J/(m3·s)。

1.1.4 水蒸氣擴散方程

濕空氣相中水蒸氣組分的擴散方程如下：

式中：w為濕空氣中水蒸氣的質(zhì)量分數(shù)，即水蒸氣含量；Dw是水蒸氣擴散系數(shù)，m2/s。

1.1.5 相變質(zhì)量傳遞速率模型

按照假設，相變只發(fā)生在霜層表面或冷板表面，發(fā)生相變的地方為相界面。在相界面處，霜相相變質(zhì)量源項，也即濕空氣相中水蒸氣變?yōu)樗嗟馁|(zhì)量傳遞速率可根據(jù)動理論[21]計算得到：

式中：ρv是相界面處水蒸氣密度，kg/m3；T是相界面的溫度，K；ρs是相界面處溫度為T時對應的飽和水蒸氣密度，kg/m3；L是網(wǎng)格單元特征長度，m；Rw是水蒸氣氣體常數(shù)，J/(kg·K)；σ是校正系數(shù)，對于水的相變一般取0.03[21]。其中，除以網(wǎng)格單元特征長度L的目的是將面相變速率轉(zhuǎn)化為體積相變速率。

能量守恒方程（4）中的能量源項Sh為相變引起的潛熱變化，在該模擬中只在相界面處為非零值：

式中：γ為相變潛熱，J/kg。

1.2 解算方法

1.2.1 模擬流程

該瞬態(tài)二維結(jié)霜模擬基于FLUENT 的VOF 方法框架通過用戶自定義函數(shù)（user defined function，UDF）實現(xiàn)。模擬的整體流程如圖1 所示，圖中n為已進行的時間步數(shù)，N為設定的總時間步數(shù)。

圖1 結(jié)霜模擬流程圖Fig. 1 Flow chart of frosting simulations

使用壓力-速度耦合解算的PISO 算法，速度使用PRESTO！方法離散。動量方程和能量方程采用二階迎風格式，水蒸氣擴散方程采用一階迎風格式。

在每個計算時間步中，可由方程（2）整理得到方程（8）。相界面上霜相的體積分數(shù)變化值，即可通過顯式求解方程（8）得到：

其中，Δt是時間步長，s；Δαf是時間步長Δt內(nèi)的霜相體積分數(shù)變化，%。由于顯式求解方程（8），穩(wěn)定計算的時間步長可提高至1 s 量級。

1.2.2 相界面判定方法

由于相變只發(fā)生在相界面，只需要對相界面上的單元格進行結(jié)霜相變計算即可。滿足如下條件的網(wǎng)格單元則判定為相界面單元（以下ε為一個極小值，例如1×10-12）：

1）單元格中霜相的體積分數(shù)介于ε和1-ε之間；

2）單元格中霜相的體積分數(shù)小于ε，但該單元格的毗鄰單元中有霜相體積分數(shù)大于1-ε的單元格；

3）單元格毗鄰低溫表面且霜相體積分數(shù)小于ε。

1.2.3 實現(xiàn)與修正方法

p-VOF 方法中霜相體積分數(shù)更新、源項計算、相界面網(wǎng)格單元搜索等功能均通過UDFs 實現(xiàn)，并關聯(lián)到FLUENT 的非穩(wěn)態(tài)計算過程中。為實現(xiàn)霜層固定不動，將霜相的黏性系數(shù)設置為一個大值（例如1×108）并將霜相的速度置零。此外，由于p-VOF 方法使用了FLUENT 的VOF 方法架構(gòu)卻未求解FLUENT原有的VOF 方程，熱導率和黏性系數(shù)不能正確傳輸用于計算，因此還需要使用UDFs 對各相的熱導率和黏性系數(shù)進行賦值，以得到正確的結(jié)霜模擬。

2 二維圓管動態(tài)結(jié)霜行為預測

對二維低溫圓管動態(tài)結(jié)霜行為進行了模擬，并比較了各種因素的影響。本研究的主要目的是確認p-VOF方法能夠用于強對流條件下低溫圓管外形的結(jié)霜模擬，且能正確反映圓管外形下流動、傳熱與霜層形貌變化的相互耦合關系。因此，卡門渦街現(xiàn)象作為次要因素而暫時忽略（捕捉卡門渦街流動的時間步長需為10-6s 量級以下，時長30 min 及以上的結(jié)霜模擬也無法承受這種巨大的時間開銷）。同時，鑒于圓管外形的對稱性并為減小計算量，使用半圓管模型網(wǎng)格進行模擬，模擬所使用的網(wǎng)格如圖2 所示，圓管半徑為5 mm。對管壁附近的網(wǎng)格進行了加密處理，如圖3所示。使用SA 湍流模型，時間步長使用1 s。模擬所使用的霜層密度為560 kg/m3、熱導率為0.2 W/(m·K)。模擬狀態(tài)如表1 所示。

圖2 結(jié)霜模擬使用的網(wǎng)格Fig. 2 Mesh for frosting simulations

圖3 管壁附近的的加密網(wǎng)格Fig. 3 Close-up view of the refined mesh around the tube

表1 結(jié)霜模擬狀態(tài)表Table 1 Parameters for frosting simulations

2.1 低溫圓管動態(tài)結(jié)霜行為

以管壁溫度250 K、來流速度10 m/s、來流溫度310 K 和相對濕度7%（對應水蒸氣含量為0.002 7）為基準狀態(tài)。二維圓管上結(jié)霜60 min 的霜層平均厚度變化情況如圖4 所示（霜層平均厚度由霜層截面積除以1/2 圓管周長得到）。在結(jié)霜初期，霜層生長速率較快，隨著結(jié)霜的進行，霜層生長速率逐漸變慢。結(jié)霜30 min 時的平均霜層厚度為0.77 mm，結(jié)霜60 min 時的平均厚度為0.93 mm，即后30 min 的霜層平均厚度僅增加了0.16 mm。

圖4 圓管表面霜層平均厚度變化情況Fig. 4 Variation of the average thickness of the frost layer with time on the tube surface

霜層形貌隨時間變化情況如圖5 所示，來流方向為從左至右，紅色區(qū)域為霜相，藍色區(qū)域為濕空氣相，黑色區(qū)域為圓管。隨時間變化，圓管上的霜層形貌呈現(xiàn)明顯的非均勻特征，反映了霜層與圓管擾流流動相互耦合作用的影響。在結(jié)霜初期，圓管前緣和后緣的結(jié)霜厚度較相近且較厚，圓管側(cè)后方位置的霜層相對較薄。隨著結(jié)霜的進行，圓管前緣和后緣的霜層厚度逐漸變大，但到20 min 以后，這兩處的霜層厚度增加已不明顯。然而，圓管側(cè)后方位置的霜層一直持續(xù)生長：在結(jié)霜初期，該處的霜層最薄，在10 min 時其厚度已經(jīng)和前后緣的霜層厚度基本相當，到20 min時，已高于其他位置的霜層，在60 min 時形成明顯凸起。

圖5 圓管表面霜層形貌隨時間變化情況Fig. 5 Evolution of frost layer morphology on tube surface

圖6 為不同時刻的圓管附近的速度場和溫度場變化圖，其中，上半部分為帶流線的速度云圖，下半部分為帶流線的溫度云圖。從速度云圖中可以看到，在圓管前緣位置處為流動滯止區(qū)，在圓管的側(cè)面靠前的位置有一個高速區(qū)，氣流在此位置有一定的加速。在圓管側(cè)后方位置，開始發(fā)生流動分離，從此處位置開始在圓管后部產(chǎn)生明顯的回流區(qū)，在回流區(qū)內(nèi)的空氣流速較低。隨著結(jié)霜進行，霜層的增厚，圓管后部的回流區(qū)逐漸變大。從溫度云圖中可知，在圓管前緣的溫度梯度大，后緣的溫度梯度稍小，在側(cè)后方的溫度邊界層呈凸出的形式，該部分的溫度梯度最小，這與圓管上流動的滯止、分離和回流直接相關。

圖6 圓管附近速度場（上半部分）和溫度場（下半部分）隨時間變化情況Fig. 6 Variations of velocity (upper) and temperature (lower)contours with time near the tube

圖7 是圓管附近濕空氣中水蒸氣含量變化情況。該圖中，圓管上的藍色部分為霜層，霜層內(nèi)的濕空氣相體積分數(shù)為0，因此水蒸氣含量也相應地為0。在5 min 時，圓管側(cè)后方位置附近的水蒸氣含量最低，隨著結(jié)霜的進行，霜層附近的水蒸氣含量逐漸升高，在60 min 時，除了霜面附近極少區(qū)域有低水蒸氣含量的地方，霜層附近其他位置的水蒸氣含量基本都已達到了來流水蒸氣含量的0.002 7。

圖7 圓管附近水蒸氣含量隨時間變化情況Fig. 7 Variation of the water vapor content with time near the tube

對比霜層形貌、溫度、速度場以及水蒸氣含量變化情況，圓管側(cè)后方位置在結(jié)霜初期霜層較薄是因為此處速度梯度小、溫度梯度小、水蒸氣含量低，導致此處的傳熱較弱、水蒸氣傳質(zhì)較慢，由此表現(xiàn)為側(cè)后方位置霜層生長較慢、霜層較薄。而圓管后緣處，雖然濕空氣流速較低，但是由于該處的流動方向是朝向圓管后緣的，使得該處的溫度梯度和水蒸氣含量相對于圓管側(cè)后方位置更大，霜層生長速率更高，使得圓管前緣和后緣處的霜層厚度無未明顯差異。

隨著結(jié)霜的進行，圓管側(cè)后方位置的霜層逐漸凸起高于其他位置的霜層。其原因是圓管側(cè)后方位置的空氣流速低且溫度梯度小，雖然結(jié)霜速率小，霜層增長慢，但此處的熱流密度更小、達到傳熱平衡時的霜層更厚。

2.2 來流速度對圓管動態(tài)結(jié)霜行為的影響

不同來流速度條件下霜層平均厚度變化如圖8所示，在結(jié)霜初期，來流速度越快，霜層生長速率越大，相同時間內(nèi)霜層的平均厚度越大。但是，來流速度越大，結(jié)霜速率明顯減慢的時間越早，霜層的最終厚度也越薄。不同來流速度條件下霜層形貌變化情況如圖9 所示，當來流速度越高，圓管上的霜層越薄，側(cè)后方的霜層凸出越明顯。

圖8 不同來流速度下霜層厚度變化情況Fig. 8 Variations of the average frost layer thickness with time on the tube surface under different air velocities

圖9 不同來流速度下圓管霜層形貌變化情況Fig. 9 Temporal evolutions of the frost layer morphology on the tube surface under different air velocities

不同來流速度條件下30 min 時刻的速度場和溫度場情況如圖10 所示，當來流速度越快，霜層越薄，是因為速度越快，對流換熱量越大，達到換熱平衡時，霜層熱阻更小，所以霜層越薄。來流速度越快，圓管側(cè)后方位置的溫度梯度較其他位置的溫度梯度相對越小，該處熱流密度和傳熱較其他位置也相對越小。因此，來流速度越高，圓管側(cè)后方位置的霜層凸起越相對明顯。

圖10 不同來流速度下速度場（上半部分）與溫度場（下半部分）分布情況（30 min）Fig. 10 Velocity (upper) and temperature (lower) contours under different air velocities at 30 min

2.3 來流濕度對圓管動態(tài)結(jié)霜行為的影響

不同濕度條件下的平均霜層厚度變化情況如圖11 所示，相對濕度（RH）分別為5%、7%和9%。隨著來流濕度的增大，霜層生長速率加快，相同時間內(nèi)霜層的平均厚度越大。這是因為來流濕度越大，作為結(jié)霜驅(qū)動力的水蒸氣濃度差也越大，這使得水蒸氣相變速率越大，結(jié)霜速率越大，從而霜層更厚。

圖11 不同來流濕度下平均霜層厚度變化情況Fig. 11 Variations of the average frost layer thickness with time under different air humidity

不同來流濕度條件下，霜層形貌隨時間的變化情況如圖12 所示。圓管上霜層的形貌，隨著來流濕度的增大，呈現(xiàn)出圓管各處霜層均增厚的情形。這是因為來流濕度的增大，只是通過式（6）來增大結(jié)霜速率，對流動和傳熱沒有直接影響，而這種濕度增大對結(jié)霜速率的增益，對圓管各處結(jié)霜的影響較一致。

圖12 不同來流濕度下霜層形貌變化情況Fig. 12 Temporal evolution of the frost layer morphology under different air humidity

不同來流濕度條件下30 min 時刻的速度場和溫度場如圖13 所示，當來流濕度越大，圓管上霜層越厚，由此導致的阻擋來流的實際圓管半徑越大，使得霜層后部的回流區(qū)域更大。當來流濕度越大，圓管上霜層越厚，霜層內(nèi)部的溫度梯度越小，但霜層外部的流場的溫度分布無明顯差異。

圖13 不同來流濕度下速度場（上半部分）和溫度場（下半部分）分布情況（30 min）Fig. 13 Velocity (upper) and temperature (lower) contours under different air humidity at 30 min

2.4 管壁溫度對圓管動態(tài)結(jié)霜行為的影響

不同管壁溫度的霜層平均厚度變化曲線如圖14所示，圓管表面溫度越低，霜層生長速率越快，相同時間內(nèi)霜層的平均厚度越大。

不同管壁溫度下霜層形貌變化情況如圖15 所示，管壁溫度越低，圓管上的霜層生長越快。管壁溫度越低，在結(jié)霜初期圓管側(cè)后方的霜層凹陷越明顯，但隨著結(jié)霜進行，側(cè)后方位置的霜層均會逐漸凸起并高于其他位置的霜層厚度。在側(cè)后方的位置會出現(xiàn)霜層有空隙的形貌，且管壁溫度越低，這種現(xiàn)象越明顯。發(fā)生這種現(xiàn)象的原因可能是，管壁溫度低、結(jié)霜速率大，一旦有擾動形成小的孔隙，孔隙周圍的霜層仍會迅速生長，并導致孔隙周圍的水蒸氣含量維持較低水平，這樣孔隙處的結(jié)霜速率將更小，這就使得孔隙特征一直保持下來而不會被填滿。這種現(xiàn)象是由于水蒸氣供給能力相對較小和水蒸氣凝華能力相對較強（溫度相對較低）的結(jié)霜環(huán)境所導致的。這種水蒸氣供給慢和相變轉(zhuǎn)化（凝華）能力強的情況，與雪花或自然對流下結(jié)霜的晶枝生長的原理類似。

圖15 不同壁面溫度下圓管表面平均霜層形貌變化情況Fig. 15 Temporal evolutions of the frost layer morphology under different tube surface temperatures

不同管壁溫度條件下30 min 時刻的速度場和溫度場如圖16 所示。從圖中可看出，當管壁溫度越低，霜層表面溫度也越低，式（6）中霜面溫度所對應的飽和水蒸氣密度越低，從而水蒸氣相變速率越大、結(jié)霜越快。此外，當管壁溫度越低，圓管上霜層越厚，由此導致的阻擋來流的實際半徑越大，使得霜層后部的回流區(qū)域越大。

圖16 不同管壁溫度時的速度場(上半部分)和溫度場（下半部分）分布情況（30 min）Fig. 16 Velocity (upper) and temperature (lower) contours under different tube surface temperatures at 30 min

3 結(jié) 論

本文p-VOF 方法可用于強對流條件下低溫圓管外形的二維動態(tài)結(jié)霜模擬，模擬結(jié)果正確反映了圓管上傳熱、傳質(zhì)、流動與霜層生長之間的耦合關系，主要結(jié)論如下：

1）氣流流經(jīng)結(jié)霜的圓管時會發(fā)生流動分離，流動分離起始點位于圓管側(cè)后方位置的霜面，圓管側(cè)后方位置的霜層在初始結(jié)霜階段較薄，隨著結(jié)霜進行，圓管側(cè)后方位置的霜層將會高于其他位置的霜層。圓管前緣和后緣位置的霜層在初始結(jié)霜階段生長較快，隨著結(jié)霜進行，其霜層生長速率減慢。

2）在常物性霜層條件下，來流濕度越大或管壁溫度越低，霜層生長速率越快，霜層越厚；來流速度越大，在結(jié)霜初期的霜層生長速率越大，但霜層生長速率明顯減慢的時間越早，最終的霜層也越薄。