強(qiáng)對(duì)流條件下使用p-VOF 方法的低溫圓管二維結(jié)霜模擬

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

(1. 清華大學(xué) 航天航空學(xué)院，熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000)

0 引 言

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

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

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

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

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

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

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

1.1 控制方程

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

1.1.1 質(zhì)量守恒方程

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

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

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

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

1.1.3 能量守恒方程

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

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

1.1.4 水蒸氣擴(kuò)散方程

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

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

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

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

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

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

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

1.2 解算方法

1.2.1 模擬流程

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

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

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

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

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

1.2.2 相界面判定方法

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

1）單元格中霜相的體積分?jǐn)?shù)介于ε和1-ε之間；

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

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

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

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

2 二維圓管動(dòng)態(tài)結(jié)霜行為預(yù)測(cè)

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

以管壁溫度250 K、來(lái)流速度10 m/s、來(lái)流溫度310 K 和相對(duì)濕度7%（對(duì)應(yīng)水蒸氣含量為0.002 7）為基準(zhǔn)狀態(tài)。二維圓管上結(jié)霜60 min 的霜層平均厚度變化情況如圖4 所示（霜層平均厚度由霜層截面積除以1/2 圓管周長(zhǎng)得到）。在結(jié)霜初期，霜層生長(zhǎng)速率較快，隨著結(jié)霜的進(jìn)行，霜層生長(zhǎng)速率逐漸變慢。結(jié)霜30 min 時(shí)的平均霜層厚度為0.77 mm，結(jié)霜60 min 時(shí)的平均厚度為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

霜層形貌隨時(shí)間變化情況如圖5 所示，來(lái)流方向?yàn)閺淖笾劣遥t色區(qū)域?yàn)樗?，藍(lán)色區(qū)域?yàn)闈窨諝庀?，黑色區(qū)域?yàn)閳A管。隨時(shí)間變化，圓管上的霜層形貌呈現(xiàn)明顯的非均勻特征，反映了霜層與圓管擾流流動(dòng)相互耦合作用的影響。在結(jié)霜初期，圓管前緣和后緣的結(jié)霜厚度較相近且較厚，圓管側(cè)后方位置的霜層相對(duì)較薄。隨著結(jié)霜的進(jìn)行，圓管前緣和后緣的霜層厚度逐漸變大，但到20 min 以后，這兩處的霜層厚度增加已不明顯。然而，圓管側(cè)后方位置的霜層一直持續(xù)生長(zhǎng)：在結(jié)霜初期，該處的霜層最薄，在10 min 時(shí)其厚度已經(jīng)和前后緣的霜層厚度基本相當(dāng)，到20 min時(shí)，已高于其他位置的霜層，在60 min 時(shí)形成明顯凸起。

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

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

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

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

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

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

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

2.2 來(lái)流速度對(duì)圓管動(dòng)態(tài)結(jié)霜行為的影響

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

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

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

不同來(lái)流速度條件下30 min 時(shí)刻的速度場(chǎng)和溫度場(chǎng)情況如圖10 所示，當(dāng)來(lái)流速度越快，霜層越薄，是因?yàn)樗俣仍娇欤瑢?duì)流換熱量越大，達(dá)到換熱平衡時(shí)，霜層熱阻更小，所以霜層越薄。來(lái)流速度越快，圓管側(cè)后方位置的溫度梯度較其他位置的溫度梯度相對(duì)越小，該處熱流密度和傳熱較其他位置也相對(duì)越小。因此，來(lái)流速度越高，圓管側(cè)后方位置的霜層凸起越相對(duì)明顯。

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

2.3 來(lái)流濕度對(duì)圓管動(dòng)態(tài)結(jié)霜行為的影響

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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