矩方法及其在水凝結研究中的應用與發展

羅喜勝, 曹 赟, 秦豐華

(中國科學技術大學 近代力學系, 合肥 230026)

0 引 言

在流體力學領域，經常會遇到大量微小粒子在流體中運動的問題。這類問題中，往往還伴隨著粒子數量和粒子尺寸分布的改變。比較典型的例子：水的同質凝結問題。對于這類問題，需要求解粒子的總量平衡方程(Population Balance Equation，PBE)。如果用拉格朗日法求解，需要追蹤大量的粒子軌跡，計算量龐大，很難運用到實際問題當中。

矩方法，又可以稱作內部坐標矩方法，通過求解內部坐標的各階矩的輸運方程，反應粒子數量及尺寸分布的變化。這里的內部坐標，可以是粒子的半徑、表面積、體積等等。以球型粒子，內部坐標取粒子半徑為例：零階矩表示粒子數量，一階矩表示粒子的半徑和，二階矩、三階矩分別和粒子的表面積和、體積和為同一量綱(之間相差一個常數)。根據內部坐標的各階矩，并對粒子分布函數進行適當假設，可以反推出粒子數量和粒子尺寸分布的信息。相比于運用拉格朗日描述的求解方法，矩方法大大簡化了計算量，使得工程應用成為可能。

1 矩方法

1.1 傳統矩方法

矩方法最早是Hulburt[1]于1964年從經典統計力學角度引入的，并且應用此方法，討論了水滴的凝結增長問題。下面對傳統矩方法做一個簡單介紹。首先，定義水蒸氣凝結生成的液滴在相空間的數密度f(ξ,t)，即分布函數。其中t為時間，ξ為可以描述諸如液滴尺寸、空間分布等狀態特性的相空間坐標：

ξ=[ξ1,ξ2,ξ3,…](1)

相空間坐標速度可以寫成：

進而得到了關于分布函數f(ξ,t)的控制方程：

其中h(ξ,t)為相空間中液滴數密度的增量。對于實際的三維空間中的相變問題，方程(3)退化成：

其中:xi為笛卡爾坐標;vi為笛卡爾坐標中液滴的速度;ri為內部坐標;Gj(c,T,r)為液滴的增長率;B(c,T,r)為成核率;c、T是水蒸氣濃度和溫度，描述了影響凝結成核與增長的環境因素。假設液滴為球型，內部坐標取液滴半徑，引入內部坐標的k階矩：

聯立方程(4)和(5)，得到了內部坐標矩的輸運方程：

觀察方程(6)，一般來說，成核率由成核理論給出，在經典成核理論中[2-5]，成核率與液滴的尺寸分布無關，因此等式右邊并不涉及具體的積分運算。等式左邊第三項為液滴增長項，由方程(6)可知，液滴增長率必須滿足與半徑成線性關系，否則k階矩方程中還將包含非k階矩，使得方程組無法封閉求解。然而，液滴的增長率并非嚴格滿足液滴半徑的線性關系，甚至相差很多，這一限制也成為傳統矩方法的重要缺陷。在實際應用中，Hill[6]提出了與半徑無關的平均增長率的概念，并且把它應用到細長管中凝結的研究中。其后的學者基本沿用這一處理方法，此方法雖然在數學上并不嚴格，但是得到的數值結果與實驗吻合較好[7-8]。

另一方面，可以通過給出粒子分布函數f(ξ,t)的具體形式來封閉方程，并且f(ξ,t)中的未知量要關于內部變量獨立。目前比較常用的分布函數為對數正態分布，假設粒子是球型，內部坐標取粒子半徑，其形式如下：

其中：rg是幾何平均半徑，σ(t)幾何標準差，N(t)是粒子總數。f(r,t)中存在3個未知量，因此至少需要3個矩變量來建立矩方程。f(r,t)中的未知量和矩變量的關系如下：

需要指出的是，f(ξ,t)的形式并不唯一，其他滿足條件的分布也是可行的。

1.2 求積矩方法和直接求積矩方法

前文提到，傳統矩方法要求液滴增長率只能是液滴半徑的線性關系，因此對許多應用都帶來限制。為了解決這一問題，McGraw[9]發展出了求積矩方法(Quadrature Method of Moments, QMOM)。這一方法的核心思想是通過高斯求積法處理各階矩，使積分形式轉變成求和形式。根據這一思想，各階矩和與矩相關的增長公式變形為如下形式：

其中:ri是高斯點，wi是權重。高斯點取得越多，計算精度越高，同時計算量也越大。綜合考慮計算量與計算精度，實際應用中一般取3個高斯點。具體計算時，首先根據P-B算法(product-difference algorithm)[10]構建一個三對角矩陣，求解該矩陣的特征值和特征向量，進一步求解得到高斯點與該點對應的權重。對于此方法的效率和穩定性，學者也進行了驗證[11]。

雖然求積矩方法最初的目的是為了解決液滴增長率的限制問題，但它的意義卻不止于此。首先，應用該方法避免了對分布函數f(ξ,t)的任何假設，f(ξ,t)的形式不需要顯式知道。其次，高斯求積思想的引入，間接體現出矩與內部坐標之間的關系，為矩方法的進一步發展開闊了思路。在實際應用中，求積矩方法不僅可以用來研究液滴的增長問題，也被用來研究凝結，液滴聚合、破碎等問題[12-13]。

求積矩方法也存在著一定的局限性。求積矩方法追蹤的仍然是各階矩，這一點與傳統矩方法并無本質區別。因此，求積矩方法不能處理相間、內部坐標間有相互作用的問題。另一方面，由于需要計算矩，求積矩方法在推廣到多組分問題時，計算變得十分復雜。因而，在求積矩方法的基礎上，Fox等[14]發展了直接求積矩方法(Direct Quadrature Method of Moments, DQMOM)，分布函數f(ξ;x,t)通過Dirac函數表示出來：

其中:N表示求積點數，NS表示組分數。將f(ξ;x,t)代入PBE中，得到直接求積矩方法的控制方程。由于直接求積矩方法追蹤的是高斯點和對應的權重，各高斯點擁有單獨的相速度。因此，此方法可以體現各組分、各相以及各高斯點所代表的相坐標之間的差別。

直接求積矩方法的應用，可以參考Fox等的文章[15-16]。

1.3 近期矩方法的發展

1.3.1 有限區域試探函數矩方法

有限區域試探函數矩方法(Finite size domain Complete set of trial functions Method Of Moments, FCMOM)[17-18]，通過一系列正交的基函數來近似粒子分布函數f(ξ,t)：

其中，Φn(ξ)是正交基函數，Cn(t)是膨脹因子?？紤]到實際問題中，粒子的尺寸是有限的，可以將分布函數轉換到有限區域[ξmin,ξmax]求解。并且可以通過無量綱化，進一步將求解域化簡到標準無量綱區域[-1,1]上。此時，原問題變成動邊界問題。[-1,1]上無量綱分布函數表示為：

進而可以得到系數的cn(t)表達式：

方程(14)中，忽略階數為負值的矩。

與其他方法相比，有限區域試探函數矩方法可以得到分布函數f(ξ,t)的近似表達式。

1.3.2 調節因子矩方法

在求積矩方法中，由于矩的數值跨度太大，會造成求解矩陣的困難。另一方面，在直接求積矩方法中，往往由于源項過大，會造成負的權值點。針對這一問題，有學者[19]提出調節因子矩方法(Adjustable MOM)。

方程(15)中，P為調節因子，通過改變P的大小，可以對矩的數值范圍進行調節。

1.3.3 擴展矩方法

Falola等[20]發展了一套擴展矩方法(EMOM)，將矩用一系列加權求和表示：

其中ni為ξi處的粒子數密度。更高階的表示，可以寫成：

這種方法可以看成是直接求積矩方法的特例，但是計算上更便捷。

1.3.4 擴展直接求積矩方法

擴展直接求積矩方法(EQMOM)的思想與FCMOM方法較為相近，可以得到粒子分布函數的具體表達形式。對比公式(13)，粒子分布函數改寫為[21]：

其中δσ(ξ;ξα)不再表示δ函數而是一個分布，如Gamma分布：

δσ(ξ;ξα)也可以用對數正態分布表示[22]。

1.3.5 處理非圓粒子

一般的矩方法中，大多假設粒子為圓形。Alexiadis等[23]引入了粒子形狀因子，研究了處理非圓粒子的矩方法。并且將此方法應用于粒子聚合問題的研究，取得了較好的結果。

2 經典矩方法在水凝結研究中的應用

如前所述，在液滴增長率與半徑成非線性關系時，矩輸運方程(6)在數學上不封閉。經典矩方法主要指根據Hill[6]提出的液滴平均增長率概念，解決封閉求解問題發展而來的矩方法。通過此方法得到的矩方程具有結構簡單，計算量小的優點。此外，經典矩方法還可以方便的與常用的商業軟件如FLUENT等結合，處理復雜流動中的水蒸氣相變問題[24]。國內外學者應用矩方法研究了水蒸氣凝結在燃燒加熱高超聲速風洞[25-26]、機翼氣動特性[27-30]、汽輪機能量轉換[31-32]、大氣中氣溶膠顆粒形成與傳輸[33-36]等眾多領域開展了大量的研究工作。下面將對經典矩方法做一簡要介紹，并且應用經典矩方法，對Prandtl-Meyer膨脹流動中的凝結問題進行了研究。

2.1 經典矩方法

k=0,1,2,...(20)

其中:

若所考慮的多相體系中，凝結液滴很小，氣液兩相的差異可以忽略，將多相流動介質作為混合物處理，ρ為混合物密度，方程(6)中液滴運動速度也替換為流動速度。與μn定義于單位體積略有不同，矩Qn定義為單位混合物質量的相關量，如零階矩Q0為單位質量混合物中的液滴數密度，而Q3為單位質量混合物中液滴所占的體積，與熱力學參數直接相關。

一般來說，所求液滴尺寸分布的矩的階數越高，矩方程個數越多，反映液滴尺寸分布的信息也越多。但考慮到具體計算的效率，一般只取低階矩進行計算。經典矩方法采取四個低階矩方程作為相變控制方程：

2.2 凝結矩方法的建立

在建立矩方程時，沒有考慮液相與氣相間的速度滑移，將液滴與氣體當作一個整體進行考慮。因此完整的控制方程包括了混合相的控制方程(不考慮黏性時為歐拉方程，考慮黏性時為納維-斯托克斯方程)和矩方程。將相變控制方程與二維氣體動力學Euler方程結合，就可以得到整個含有相變過程的高速流動控制方程組:

其中:U為變量，F和G為通量，S為源項。分別有：

其中E表示單位質量混合物的總能量。由分布函數的定義，可知單位體積混合物中包含的液滴總體積：

由此可以計算液態水的質量分數：

ρl為液態水的密度?？刂品匠讨械娜A矩方程也可改寫：

上式實際上是液態水質量的對流方程。上述方程組不封閉，還需要添加氣體狀態方程(由于液相質量占比很小，因此忽略液相體積對壓強的影響)和壓力、溫度修正方程。氣體狀態方程在很多文獻中都有涉及，此處不做贅述。

通過求解矩的演化過程，得到各階矩隨時間的變化情況，再對液滴分布函數作適當假設，就可以推導出液滴尺寸分布的信息。

2.3 伴隨同質凝結的Prandtl-Meyer膨脹流動

超聲速氣流經過一個轉角，氣體膨脹，馬赫數上升，壓力、密度、溫度降低，形成一列中心稀疏波。在二維情況下，Prandtl和Meyer首次給出了流動的解析解，被收入大多數空氣動力學教材中，因此這一列稀疏波又被稱為Prandtl-Meyer(以下簡稱PM)膨脹波。當氣體中包含水蒸氣時，沿流線氣體膨脹，壓強降低，水蒸氣分壓也隨之降低；同時，溫度降低，水蒸氣飽和蒸氣壓降低。由于水蒸氣飽和蒸氣壓與溫度呈冪指數關系，從而導致水蒸氣飽和蒸氣壓更快的降低，使得水蒸氣飽和度S(水蒸氣分壓與飽和蒸氣壓之比)沿流線逐步升高。當飽和度S大于1時，就可能發生凝結。

伴隨凝結的PM流動，在科學實驗和工業生產中也是常見的現象。例如，在超聲速飛機機翼背風面，超聲速氣流膨脹，速度加快，溫度降低。如果來流中水蒸氣含量較高，水蒸氣有可能發生凝結。凝結放出潛熱，使局部壓力提高，形成各種波系結構(后文稱凝結波或凝結激波)，改變了機翼上的壓力分布，對機翼的受力產生影響。而在工業生產中，蒸汽輪機中的凝結問題也可以簡化成伴隨凝結的PM流動問題。一方面，凝結放熱，影響流場結構，進而影響蒸汽輪機的效率。另一方面，在這種有限的區域內，由于凝結的自我抑制性，凝結波可能會產生自激振蕩，給葉片施加一個周期性的載荷，影響葉片的使用壽命，同時產生噪聲。

對于伴隨凝結的PM流動問題，已有不少學者開展了相關研究。Smith[37]對典型的二維超聲速流動繞過尖角膨脹形成的PM波中的凝結問題進行了實驗研究，發現凝結區域向來流方向彎曲。這表明之前Steffen[38]研究此類問題時提出的假設“從拐點發出的射線上放熱量相同”并不恰當。Kurshakov等[39]同樣進行了實驗研究，發現了膨脹扇中凝結導致的激波。凝結激波并不與拐點相連，而是在距拐點一定距離的區域形成。Frank[40]則首次觀測到了凝結激波的非定常振蕩現象，并對伴隨凝結的定常和非定常流動進行了比較。Delale等采用漸近方法分析了PM流動中的凝結現象，并分為沒有凝結激波的亞臨界情況[[41]和存在凝結激波的超臨界情況[42]。結果表明：亞臨界時凝結區域向來流方向彎曲，這與Smith的實驗結果吻合；超臨界時給出了凝結激波的起始位置，并不與拐點相接，也符合Kurshakov的實驗觀測結果。在具體的應用方面，White等[43]對蒸汽輪機葉片處的凝結問題進行了實驗研究并進行了數值模擬。實驗結果與數值結果吻合得較好。Kim等[44]通過實驗和數值研究，提出了一種控制凝結的多孔結構，可減弱膨脹扇中凝結激波強度，且能完全抑制凝結激波的振蕩。

PM流動膨脹會導致凝結現象發生，凝結釋放熱量將改變流場結構，流場變化進一步影響凝結的進行，兩者耦合在一起，使得流動更為復雜。在一定參數范圍內還可能產生振蕩現象，呈現出強烈的非定常效應。這使得對伴隨凝結的PM流動的分析變得相當困難。有界區域內的繞角PM流動與噴管中單純的膨脹流動都可以看成是氣流經過一列簡單波進行膨脹，流動特征存在相似之處。對細長噴管中的凝結問題以及凝結波的振蕩模式都有相應的研究[6-7,45]，而對于伴隨凝結的繞角膨脹的PM 流動，構型不同導致空間演化更為復雜，正如前人所說[46-48]，仍需要進一步研究，深入認識其規律。數值模擬是研究凝結與流動耦合作用的一種有效手段，本節采用經典矩方法對PM流動中水蒸氣同質凝結對流動的影響以及凝結放熱引起的凝結激波進行了研究。

流動模型依據超聲速G2噴管進行設計，如圖1所示。噴管喉道上游的收縮段與原始的超聲速G2噴管收縮段相同。而在喉道下游的擴張段，上壁面是一個與噴管軸線平行的水平面，下壁面是與噴管軸線向外偏轉30°角的斜面，拐點在喉道下游0.072 cm的位置，整個擴張段延伸到喉道下游6.5 cm處。超聲速來流由二維噴管產生，在噴管喉道下游繞角膨脹，形成PM流動。氣流經過短暫的膨脹，在拐點處略微超過聲速，當地馬赫數約為1.035。通過逐漸提高來流中水蒸氣的飽和度，研究不同條件下凝結波的產生和運動形態。流動介質為氮氣和水蒸氣的混合氣體，僅考慮水蒸氣的同質凝結。計算采用自適應網格，根據前人應用經典矩方法研究凝結問題的經驗[8]，自適應加密后網格最小尺度達到0.1 mm量級，即可很好捕捉到各種凝結流場的結構。

圖1 伴隨凝結的PM流動計算示意圖Fig.1 Sketch of Prandtl-Meyer flow with condensation

初始時刻，在喉道下游設置一個間斷面，間斷面左側上游設置高壓，右側下游設置低壓。計算開始后，間斷面處產生向上游運動的膨脹波和向下游運動的激波，并逐漸在噴管和膨脹扇中匹配出穩定的流場。如果直接啟動凝結計算，凝結放熱產生額外的凝結波，干擾了流動的匹配過程，使得流場很難穩定下來。因此，對于每個算例，首先模擬了未發生凝結情況下的流場，當流動達到穩定時再開啟凝結計算。由于來流條件是一定的，通過不同方式建立起來的流場，最終都會趨向同一個結果。

首先，對來流飽和度較小的情況進行了模擬。噴管入口處的溫度Tini=330 K，壓強pini=1 atm，水蒸氣飽和度Sini=0.45。此來流條件下，噴管中水蒸氣飽和度仍然小于1，經過拐點膨脹后，才會發生凝結。作為比較，同時計算了相同來流條件下，不啟動凝結的情況。結果如圖2所示，圖中箭頭指向流動方向?？梢姰斄鲃又邪殡S有水蒸氣凝結時，溫度場不再相對于拐點呈射線狀分布，而是由于凝結潛熱的釋放，在膨脹扇中出現一個高溫區域；在距離拐點較遠處，由于凝結激波的影響，溫度場還會發生一定程度的扭曲。凝結對壓力場的影響相對較小，壓力等值線仍大致呈射線狀分布，說明PM 流動的膨脹作用相較于水蒸氣凝結的減質添熱效應對壓力場的影響更大。同時還可以看出壁面對壓力場的影響：壓力等值線趨向垂直于上壁面，當啟動凝結計算后，這種趨勢更加明顯。

圖2 啟動凝結(實線)與不啟動凝結(虛線)時溫度、壓力分布對比Fig.2 The distribution of temperature and pressure in PM flow with (solid line) and without (dashed line) condensation

從有凝結時數值紋影圖(圖3(a))中可以清晰地看到由水蒸氣凝結放熱產生的激波，激波從膨脹扇中發出，向來流方向彎曲，同時又距離拐點一定距離。這些特征都與前人的理論分析和實驗結果吻合。流動膨脹過程會帶來密度變小，越靠近拐點，膨脹越劇烈，因此紋影圖中拐點附近顯示為高密度梯度區域；激波同樣會導致密度變化，波后密度升高，這與膨脹過程正相反。凝結放熱隨著距拐點距離縮小而逐漸集中，最終形成激波且強度逐漸加強；進一步靠近拐點，流動膨脹作用增強并逐步超過凝結激波作用，激波逐漸減弱并消失，拐點附近主要體現為膨脹帶來的密度快速降低。這也是Kurshakov等觀測到凝結激波不從拐點發出的物理原因。

(a) 拐角為尖點

(b) 拐角為圓角

從圖3(a)中還可以看到，在拐點下游附近、靠近斜面的區域也存在一道激波。這是由于拐點附近膨脹過快，凝結還來不及充分發展，氣流已流過膨脹扇。然而在膨脹扇之后，水蒸氣飽和度依然很高，水蒸氣會繼續凝結放熱，形成了該激波。實際上，拐點在數學上是一個奇點，該點附近膨脹率和冷卻率都趨向于無窮大。目前水蒸氣同質凝結過程中使用的穩態成核理論和基于準靜態的液滴增長理論在此流動條件下都存在問題。基于此，設計了另一個轉角模型，用半徑r=2mm的圓角替換了原拐角尖點，數值紋影如圖3(b)所示?？梢钥吹皆诰嚯x拐點較遠的地方，流動形態和凝結激波位置并無顯著差別；但在轉角附近，凝結激波在距離壁面很近的地方還可以清晰的看到。這說明圓拐角確實弱化了此處的極端冷卻率和膨脹率。同時，由于冷卻率的降低，轉角下游壁面處的激波形態也發生了改變。與尖點拐角相比，差別僅出現在拐角附近區域，說明拐角構型對凝結的影響具有局部性，在大尺度情形下可以忽略其影響。另外，真實流動在壁面附近存在邊界層，邊界層中流體速度變慢，膨脹率和冷卻率必然不會很大，邊界層對伴隨凝結的PM流動的影響還需要進一步討論。

當來流水蒸氣含量較高時，會出現凝結激波振蕩現象。氣流繞過拐角偏轉膨脹后發生凝結，釋放大量潛熱，促使膨脹扇中出現凝結激波；波后溫度升高，凝結減弱，并促使凝結激波向上游運動。而逐漸減弱的凝結過程釋放熱量降低，導致激波強度逐漸減弱，最終在喉道附近消失。此時，流場恢復到了無凝結狀態，重復之前凝結、激波相互影響的過程，呈現出凝結激波的振蕩現象。如果提高來流中水蒸氣含量，也即提高入口處的水蒸氣飽和度Sini，凝結激波的振蕩頻率也會相應提高。為了統計激波振蕩頻率規律，保持入口來流壓力不變pini=1 atm，計算了兩種來流溫度Tini=310 K、330 K下，振蕩頻率隨來流水蒸氣飽和度的變化關系，如圖4所示。來流溫度相同時，凝結激波的振蕩頻率與水蒸氣飽和度基本成線性關系。來流溫度越高，相同飽和度下水蒸氣含量越高，凝結量越大，凝結激波越強，振蕩頻率越高、隨水蒸氣飽和度上升得越快。

圖4 凝結激波的振蕩頻率隨來流溫度和來流水蒸氣飽和度的變化關系Fig.4 Oscillation frequency of condensation wave in dependency on the temperature and saturate ratio of incoming flow

3 矩方法在水凝結研究中的新發展

3.1 異質凝結矩方法

作為另一種主要的凝結現象，異質凝結也在生產和生活中占有重要地位。水蒸氣異質凝結與同質凝結的主要區別在于形成凝結核的過程不同，同質成核是指水蒸氣分子在一定條件下自發團聚成凝結核，而異質成核是指水蒸氣在其他物質表面的凝結成核。生活中常見的雨、霧等自然現象，以及運用異質凝結清除顆粒等，都是水蒸氣異質凝結的典型例子。目前關于異質凝結的研究主要集中在凝結的機理上，有關凝結與流動相互作用的研究還比較少。如前所述，矩方法將流動與凝結過程作為一個完整體系建立控制方程，正適合研究凝結與流動的相互作用。本節將經典矩方法進行拓展，應用到異質凝結領域，進而對激波管中的異質凝結問題進行參數研究。

3.1.1 異質凝結矩方法的建立

與同質凝結一樣，異質凝結也包含成核與增長兩個子過程。異質成核中，研究較多的是水蒸氣在雜質粒子上的成核過程。當雜質粒子完全被液態水包裹，則進入增長過程，雜質粒子的影響減弱以致于可以忽略。換言之，異質成核生成的液滴與同質成核生成的液滴在增長過程上沒有差別。因此，異質凝結矩方法的建立主要是將異質成核過程整合到矩方程中。

關于異質成核過程的研究，有很長的歷史。最早Volmer [49]通過動理論給出了不可溶物質平板表面上的異質成核率公式，Fletcher [50]將其推廣到不可溶球形粒子表面，Kotake等 [51]引入了表面張力的影響。Girshick等 [5]提出的應用在同質成核理論中的自洽平衡態液滴數密度分布也被應用于異質成核理論。最近，范煜等 [52]綜合了上述理論并添加了線張力的影響，改進了異質成核公式。當前也有不少學者在分子尺度研究異質成核過程的特征與機理 [53-56]。雜質粒子上異質成核過程十分復雜，大致包含3個階段：首先，受固體粒子的影響，水蒸氣分子優先在顆粒表面凝結生成大量離散的液簇；這些液簇逐漸增長、碰并，數量減少、體積變大，形成附著于顆粒表面的球冠狀液滴；最終所有液滴合并為完整的液膜，完全覆蓋顆粒表面，異質成核結束。此時的雜質粒子被稱為活化粒子，活化粒子可以被看作液滴，隨后的增長過程也與液滴的增長相同。

然而在宏觀流動計算方法中完整模擬成核過程是非常困難的，本文采用“瞬間活化”假設對整個異質成核過程進行近似[57]。由同質成核理論可知，只有當液滴半徑大于臨界半徑r*時液滴才會凝結增長。既然異質、同質凝結在增長階段沒有差別，這一結論也適用于異質成核生成的液滴。如果忽略液膜厚度，小于臨界半徑的粒子即使表面上被覆蓋一層液膜，此液膜也會因為蒸發強于凝結而趨向消失。因此假設只有半徑大于臨界半徑的顆粒才能活化，且忽略具體成核過程，直接形成液膜并增長(圖5)。此處的臨界半徑指的是根據當前的熱力學狀態，由同質凝結理論計算得到的臨界半徑。對于半徑在百納米以下的粒子，根據范煜等[52]提出的成核公式可知，粒子活化的時間很短，表面覆蓋的液膜很薄，“瞬間活化”假設是合適的。

雜質粒子一般不是均勻的顆粒，假設粒子為圓球形，半徑分布滿足函數g(r)。隨著流動狀態的改變，部分粒子被活化并發生凝結，活化粒子半徑分布h(r,r*)不僅與粒子半徑有關，也與臨界半徑r*有關。由“瞬間活化”假設可得

圖5 “瞬間活化”過程示意圖Fig.5 Schematic of instantaneous-wetting model

流場發展過程中，臨界半徑隨時間變化r*=r*(t)，活化粒子分布也隨時間變化，則成核率Jhe=?h/?t，從而矩輸運方程(6)可寫為

k=0,1,2,…(29)

其中Qk仍然是關于液滴尺寸分布的k階矩。異質凝結中液滴尺寸分布與同質凝結有無差別，以及雜質顆粒尺寸分布與液滴分布的關系，還是需要探討的問題。上式與同質凝結矩方程(20)類似，主要差別僅在于成核過程項:

該項主要由雜質粒子尺寸分布及臨界半徑決定。能夠將雜質粒子尺寸分布的影響也包含到計算模型中，也是異質矩方法的特點之一。

異質凝結矩方法的控制方程也由歐拉方程和矩方程(29)中前四階構成，與同質凝結矩方法的控制方程(23)、(24)形式相同，僅需將源項S替換為:

由于異質凝結的三階矩也包含了雜質顆粒的體積，三階矩方程不能替換為液滴質量對流方程。由于活化顆粒的尺寸分布函數已知，可以很容易的計算出已凝結雜質顆粒的總體積

因此，凝結液態水質量分數

與同質凝結矩方法類似，異質凝結矩方法也將液相與氣相當作均勻混合物處理，不考慮相間作用。因此，凝結放出的潛熱的影響，需要通過壓力修正與溫度修正來體現。其修正過程與同質凝結矩方法類似，此處不再贅述。

3.1.2 激波管中的異質凝結

我們用拓展的矩方法對激波管中的異質凝結問題進行了研究。如圖6所示，激波管總長L=1 m，其中充滿氮氣和水蒸氣的混合氣體，氣體中均勻分散有不可溶的惰性顆粒，充當異質凝結的凝結核。激波管中間D處設有一道隔膜，將整個激波管分為高壓區(high-pressure section，HPS)和低壓區(low-pressure section，LPS)，壓強分別為1 atm和0.5 atm；溫度相同，均為295 K，水蒸氣飽和度分別為0.9和0.45。 計算域邊界均采用固壁條件。在P1處設有探針，記錄各參量隨時間的變化情況。

圖6 封閉激波管算例示意圖Fig.6 Sketch ofheterogeneous condensation in a closed shock tube

為了對比分析凝結對流動的影響，分別計算了考慮與不考慮異質凝結兩種情況，不考慮凝結僅關閉凝結源項計算即可。圖7給出了兩種情況下激波管軸線上的密度在x-t平面內的等值線圖。

在零時刻，激波管中間的隔膜被打開，高壓段和低壓段聯通，產生了一道右行的激波(SW)和接觸間斷(CS)。激波過后，低壓段密度升高。同時，一道膨脹波開始向左側的高壓段傳播，并導致高壓段密度降低。之后，膨脹波在高壓段左側壁面處反射，形成反射膨脹波(RE)，進一步降低了高壓段密度。另一方面，右行激波也在低壓段的右側壁面反射，形成反射激波(RS)，進一步壓縮了低壓段的混合氣體。反射激波與接觸間斷相交，產生一道左行的透射反射激波(TRS)和一道右行弱激波(S1)。透射反射激波又會和反射膨脹波相交，并在高壓段左側壁面再次反射。由此可見，激波管問題包含了多種非線性波系間的相互作用。從圖7(a)中可以看出，各種波系以及不同波系間的相互作用都被很好的模擬出來，這也從另一方面說明了此程序在模擬可壓縮流動中，有很高的可靠性。當可壓縮流動中還伴隨有異質凝結時，凝結與不同波系之間又會產生相互影響。從圖7(b)中可以明顯看出，激波管中的波系結構變得更加復雜。這主要由凝結時潛熱的放出和液滴蒸發時從周圍氣體吸熱量所導致的。

(a) 不考慮異質凝結

(b) 考慮異質凝結

SW：激波，CS：接觸間斷，EF：膨脹波，RS：反射激波，TRS：透射反射激波，S1：接觸間斷上反射的弱激波，RE：反射膨脹波，HC1、HC2：異質凝結區域，EV1、EV2、EV3：蒸發區域。

圖7 激波管軸線上的密度在x-t平面內的云圖
Fig.7 Density contours for the closed shock tube problemin a space-time diagram without condensation (a)and with heterogeneous condensation (b)

對比有無異質凝結的密度x-t圖可以看出，激波管破膜后，膨脹波經過高壓段，使高壓段中混合氣體溫度降低，水蒸氣飽和度升高，異質凝結發生(HC1)。當膨脹波在高壓段左側壁面反射后，高壓段溫度繼續降低，異質凝結再次發生，在圖7(b)中用HC2標記出來。凝結放出潛熱，被氣體吸收，導致局部壓力升高，產生壓力波，并影響了膨脹波EF和反射膨脹波RE，使其形狀產生畸變。而右行激波經過反射和透射，透射反射激波(TRS)進入凝結區域，激波后溫度升高，水蒸氣飽和度降低，液滴開始蒸發(EV1)。而后，透射反射激波穿過反射膨脹波的膨脹扇，抑制了膨脹扇中的凝結，導致膨脹扇發生扭曲，并使反射膨脹波后的液滴開始蒸發，形成第二個蒸發區EV2。最后透射反射激波經過高壓段左側壁面的再次反射，使高壓段溫度進一步提高，形成了第三個蒸發區EV3。 蒸發過程中，液滴吸收氣體的熱量，在每個蒸發區都會形成一系列膨脹波，這一點與同質相變的結果相似[7, 58]。而與同質相變的不同之處在于，同質凝結往往發生在水蒸氣飽和度較高的條件下(理論上水蒸氣飽和度超過1就會發生同質凝結，但顯著的同質凝結往往發生在水蒸氣飽和度3以上的環境中)。這就導致同質凝結的強度較大，放熱較為集中，形成顯著的凝結波甚至是凝結激波，反觀異質凝結就較為平緩。但對于同樣的初始條件，不同凝結過程的放熱總量是相同的，異質凝結對波系的影響也不容忽視。

3.2 相間滑移矩方法

經典矩方法已被證明在處理凝結問題時高效、實用。然而，如前所述，不論處理的問題是同質凝結還是異質凝結，都有一個共同的基本假設，即忽略液相與氣相之間的相互作用，將凝結生成的液滴和氣體當作均勻混合物處理。一般來說，在液滴較小、流動較為平穩的情況下，液滴與氣體熱力學狀態、流動狀態均較為一致，即相間作用不顯著，該假設是成立的。當流動變化劇烈(如激波、旋渦等)、液滴較大(如燃燒加熱風洞中的同質凝結液滴尺寸可增長到微米級)時，氣相與液相之間除了水蒸氣凝結、液滴蒸發帶來的相變潛熱釋放、吸收之外，還存在著由相間速度差異所帶來的動量、能量交換以及由溫度差異所帶來的熱量傳導。在某些條件下，兩相間存在著強烈的相互作用，或相間作用不強但持續時間長，相間差異被積累放大。如水蒸氣在旋渦中的凝結問題[58]，圖8給出了伴隨水蒸氣同質凝結的旋渦流動的消光結果。從圖中可以看出，旋渦中心出現一個“白點”，是因為液滴跟隨氣流持續旋轉，在離心力的作用下離開旋渦中心，并在邊緣處聚集。這一現象通過經典凝結矩方法無法重現，并且由于旋渦中心溫度最低，凝結最劇烈，經典矩方法模擬給出此處的液態水量是最多的，得到完全相反的結論。由此可見，在此條件下，氣液兩相間相互作用時不可忽略的。

圖8 伴隨凝結啟動渦的實驗消光結果[58]Fig.8 Light-extinction image of vortex shedding in Ludwieg tube[58]

為了考慮相間相互作用，一些學者將氣液兩相區分開，通過拉格朗日法追蹤液相顆粒運動，氣相仍用歐拉方法進行描述[59]。這種方法直觀、可靠，但凝結生成的液滴數量極大，拉格朗日法計算量巨大，即使依靠經驗對液滴進行分組，仍然需要追蹤大量的液滴簇路徑，大計算量是該方法最主要的瓶頸。也有學者針對矩方法進行改進，Marchisio等[14]在求積矩方法的基礎上，提出了直接求積矩方法。由于每個積分點都有獨立的相速度，因此通過添加必要的代數關系，就可以得到相間的相互作用。Brown等[60]在研究氣溶膠問題時，多分散膠體系統用矩方程進行描述，其平均速度包括了慣性速度、擴散速度和膠體熱泳速度的修正，考慮了相間的質量和熱量交換，并發展出一套計算程序。但是，沒有考慮相間的動量交換，不能描述液相對氣相的影響。

本節在經典矩方法的基礎上，將氣液兩相分開處理，建立了一套考慮兩相間質量、動量和能量交換的凝結矩方法。此方法與一般的雙流體模型方法[61-62]類似，氣液兩相都采用歐拉方法進行描述，相間具有速度滑移。液相平均速度與液滴尺寸分布無關。由于各階矩描述的是液滴尺寸分布的信息，很自然的矩方程的對流速度采用液相的平均速度。與直接求積矩方法相比，此方法具有經典矩方法物理意義明確的特點，并且更容易通過計算流體力學方法實現。而與一般的雙流體模型法相比，此方法并非簡單的對液相參數進行統計平均，可以提供更多液相的信息。

3.2.1 相間滑移矩方法的建立

與經典矩方法不同，由于考慮了相間的速度滑移，控制方程包含了三個部分：氣相控制方程、液相控制方程以及矩方程。對于無粘問題，氣相控制方程為歐拉方程；液相為彌散于氣相中的液滴，忽略其對壓力的影響，僅考慮液相的對流輸運過程。凝結過程則由矩方程描述，需要注意的是經典矩方法中液滴與氣相具有相同的速度，考慮速度滑移效應時，矩方程中對流速度應是液相速度。水蒸氣凝結帶來的質量交換和潛熱釋放、相間速度滑移引起的動量和能量交換都通過源項實現。因此，二維無粘情況下相間滑移矩方法控制方程的各項如下：

其中下標“g”和“d”分別代表氣相(gas)和液相(droplet)。α是氣相中水蒸氣的質量分數，三階矩方程替換為水蒸氣質量守恒方程。其它符號與前文一致。

源項Si表示由于相變帶來的各類相間交換，是相間滑移矩方法的核心。S1、S5、S9分別是水蒸氣相變導致的氣相、液相及水蒸氣質量變化源項。水蒸氣在成核過程生成臨界半徑r*液滴，消耗水蒸氣帶來等量的液相質量增長率:

而液滴增長過程同樣帶來液相質量增長率:

因此相變帶來的總質量變化率Sm=Sm1+Sm2，從而

S1=S9=-S5=-Sm(37)

動量源項S2、S3、S6、S7包含伴隨相變質量轉換的動量交換和相間速度差異產生黏性力帶來的動量交換，故:

S2=-S6=-fx-Sm1ug-Sm2ud

S3=-S7=-fy-Sm1vg-Sm2vd(38)

對于剛成核的液滴，認為其速度為當地氣相的速度。而液滴增長帶來的液相增加質量，認為其速度為當地的液相平均速度。因此，成核帶來的動量變化根據當地氣相速度計算，液滴增長帶來的動量變化根據當地的液相平均速度計算。f=(fx,fy)表示相間速度滑移帶來的黏性作用力，其表達式如下：

能量源項S4、S8包含相間速度滑移帶來的能量交換、相間溫度差帶來的熱傳導和凝結帶來的潛熱釋放：

S4=-fxud-fyvd-q-Sm1Ec1-Sm2Ec2

其中，氣液兩相間熱交換率[64]:

其中，Tw1、Tw2分別是成核過程和增長過程液滴表面的濕球溫度[65]。成核時，凝結核較小，認為整個核的溫度都是Tw1；增長過程中的液滴內部平均溫度為液相溫度Td，而表面濕球溫度Tw2是新增長液膜的溫度。一般來說，Tw1、Tw2和Td都不相同。凝結釋放的潛熱L是溫度的函數[8]，這里假設潛熱全部被氣相吸收，因此潛熱項只出現在氣相能量方程的源項中。另外，由于成核和增長過程中生成的液態水具有不同的速度，攜帶的動能也不同，分別用液相平均速度和氣相速度計算。

與經典矩方法相比，通過上述源項使得相間滑移矩方法對同質成核和液滴增長的模擬更接近真實物理過程。

矩方程的源項S10、S11、S12與經典矩方法并無差別，為了完整重寫如下：

S12=J(43)

同樣的，上述方程組還需要用氣體狀態方程進行封閉。如果忽略液相對氣相壓力的影響，則氣體狀態方程與之前并無不同。由于各相擁有單獨的能量方程，并且水蒸氣的相變潛熱已經在能量方程的源項中體現出來，因此與經典矩方法不同，這里不再需要進行壓力和溫度的修正。

3.2.2 伴隨凝結的旋渦運動

為了考察相間滑移的影響，采用前述相間滑移矩方法計算了伴隨凝結的旋渦運動。

首先討論較為基礎的蘭金(Rankine)復合渦中的凝結問題。蘭金渦由速度線性增長的渦核和二維位勢渦組合而成，速度分布如下：

其中，R為渦核半徑，Ω為旋轉角速度。這一速度分布與氣象學中的熱帶氣旋類似，常被作為研究熱帶氣旋的初步理論模型。需要指出的是，蘭金渦一般被作為不可壓旋渦來研究的，但對于可壓縮流動，仍然可以借用其速度分布。其壓強分布，可以通過可壓縮流動的伯努利方程求解。

計算中流動介質為氮氣和水蒸氣的混合氣體，通過調整二者質量比實現不同的初始飽和度。渦核半徑R=2 cm，轉動角速度Ω=14 870 s-1，這里角速度設置較大，以產生足夠強的旋渦。選取半徑100 cm的圓形計算域，如圖9所示，渦核位于計算域中心。由于計算域半徑遠大于渦核半徑，邊界條件的影響可以忽略。由于流動中沒有強間斷，最小網格尺度可略大，設為0.5 mm。遠場壓強p0=100 kPa、溫度T0=298 K。

計算初始時全場賦值為未發生凝結的穩定蘭金渦流場，水蒸氣處于過飽和狀態。而后啟動計算，凝結隨之發生。作為示例，圖10給出了遠場水蒸氣飽和度S0=1.2時旋渦中心的液相密度云圖，顏色越深密度越大。旋渦外圍由于飽和度不是太高，并沒有同質凝結發生。隨著向旋渦中心靠近，壓力、溫度均降低，凝結也逐步變得劇烈，液相密度增加。同時，凝結放出潛熱，產生凝結波向外側傳播。早期凝結波較強，可能會抑制外側凝結的產生，這從t=10 μs的液相密度分布可以清晰地看到。凝結生成的液滴，由于相間滑移并不嚴格跟隨氣體運動，而是在離心力的作用下沿徑向向外遷移，導致旋渦中心液相密度降低。在t=410 μs之后，液相密度分布圖顯示，旋渦中心出現“白點”，并且隨時間增大。

另外，還發現一個有趣的現象，隨著時間發展，液相密度分布逐漸出現分層結構，形成一系列同心環。在渦核內部，氣流旋轉速度隨半徑近似線性變快，液滴在向外遷移的過程中被氣流帶動加速，受到更大的離心力，進一步促進向外遷移，這也是中心出現“白點”的原因；而在渦核外部，氣流速度隨著半徑降低，液滴則被氣流帶動減速，離心力減小，向外遷移擴散速度變慢。渦核內產生的液滴將在渦核邊緣聚集，形成內部第一道暗紋。在距離旋渦中心較遠的區域，水蒸氣飽和度雖然高于1，但不足以支持顯著的同質凝結成核過程，當內部液滴遷移到此區域時，沒有消耗的水蒸氣會促使液滴繼續凝結增長。這使得旋渦外側的液滴尺寸要比內側大的多，而大液滴跟隨性更差，相間滑移作用更強，徑向遷移更快。同時，凝結釋放潛熱使得溫度升高、飽和度降低，抑制液滴增長，降低徑向遷移。在相間滑移與凝結自我抑制的共同作用下，液滴在局部位置聚集，液相密度分布呈現環狀結構。

隨后，應用相間滑移矩方法對啟動渦中的凝結問題進行了數值研究。參照路德維希(Ludwieg)管中斜劈誘導旋渦中的凝結實驗[8]設置計算條件，如圖11所示。計算采用氮氣和水蒸氣的混合氣體。計算域長110 cm，高10 cm。斜劈尺寸如圖所示，左邊界位于x=17 cm處。隔膜位于x=90 cm處，左側為高壓區，壓力100 kPa，水蒸氣飽和度S=0.88；右側為低壓區，壓強為22 kPa，水蒸氣組分含量與高壓區相同。初始溫度296 K。計算網格自適應加密，最大加密4層，加密后最小網格尺寸約為0.1 mm。

圖9 蘭金渦算例網格示意圖Fig.9 Computation domain and mesh for Rankine vortex flow with condensation

圖10 蘭金渦中心區域液相密度演化, S0=1.2Fig.10 Variation of liquid phase density distribution in Rankine vortex, S0=1.2

計算啟動后破膜，間斷面處產生向高壓段運動的膨脹波。當膨脹波到達斜劈后，在尖角處誘導產生旋渦。隨著時間推進，旋渦不斷增強，旋渦中心的溫度不斷降低，飽和度升高，隨后發生凝結。圖12給出了不同時刻，旋渦中的氣相溫度和液相密度分布云圖。從圖中可以看出，當旋渦中的溫度下降到250 K左右時(t=2.6 ms)，水蒸氣開始發生同質凝結。此時，旋渦外的溫度要高于旋渦內的溫度。隨著時間的推進，旋渦緩慢向右運動，強度逐漸上升，渦心凝結加劇(t=2.8 ms)，旋渦內的液相密度開始升高。伴隨著凝結潛熱的放出，渦心溫度迅速升高，但還是略低于外部溫度。此時溫度最低處位于剪切層內。在t=3.0 ms時刻，旋渦中心過飽和的水蒸氣已經凝結完全，渦心的溫度隨著旋渦強度的增加開始緩慢下降。凝結主要發生在旋渦邊緣，可以觀察到，在障礙物尖角右側溫度最高。在t=3.2 ms時刻，膨脹波中的凝結波與旋渦接觸，兩者作用的結果使得旋渦右側出現凝結，此處溫度升高。而旋渦中的溫度繼續緩慢降低。從液相密度圖上還可以發現，剪切層上出現微弱的不穩定現象。在t=3.4 ms時刻，凝結波與旋渦作用，加劇了剪切層上的擾動。渦心的溫度由于凝結波的影響開始回升。在t=3.6 ms時刻，凝結波基本掃過旋渦，渦心溫度進一步上升。剪切層上的界面不穩定性發展出小的渦結構，并且伴隨著大渦的運動被卷入旋渦內部。對比不同時刻的液相密度分布可以發現，旋渦中水蒸氣凝結完全后，渦心的液相密度就開始逐漸降低。液滴隨著氣流持續旋轉，在離心力的作用下向旋渦外側運動。這說明相間滑移的影響在伴隨凝結的旋渦運動中不可忽略。

圖11 斜劈誘導啟動渦算例示意圖Fig.11 Sketch of vortex shedding in Ludwieg tube with water vapor condensation

圖12 啟動渦中氣相溫度場(左)及液相密度場(右)的演化過程Fig.12 Distribution of gas phase temperature (left) and liquid phase density (right) during the vortex shedding process

圖13給出了凝結產生前期(t=2.8 ms)和后期(t=3.75 ms)啟動渦中液相流線結果?？梢钥闯?，液相流線從旋渦中心發出，盤旋向外。這再次表明液滴在離心力的作用下，逐漸離開旋渦中心。前期凝結量不大，液滴僅在旋渦范圍內運動。液相流線終止于斜劈壁面上。隨著流動發展，液滴被甩出旋渦，液相流線向下游延伸。同時還可以發現，斜劈尖端的脫渦會對流線產生擾動。觀察流線的疏密程度可以發現，靠近旋渦中心，流線較密，而在旋渦邊緣，流線間距較大。這說明液滴在離開旋渦的過程中，向外速度逐漸加快，與前述蘭金渦中液滴的遷移特征也是一致的。

圖14給出了t=3.75 ms時刻，實驗和不同計算方法得到的液相密度分布結果。經典矩方法的計算模擬結果顯示，旋渦中心的液相密度最高。這是因為旋渦中心凝結最為劇烈，液滴生成后完全跟隨氣流運動，無法離開旋渦中心，與實驗結果是不符的。而相間滑移矩方法模擬能準確重現消光實驗中觀測到的旋渦中心的“白點”，同時與紋影實驗結果一致，可以捕捉從斜劈尖端分離的渦結構。

圖13 啟動渦中凝結產生前期(t=2.8 ms)和后期(t=3.75 ms)液相流線示意圖Fig.13 Streamlines of the liquid phase in the starting vortex at different time

(b) 紋影實驗

(c) 經典矩方法

(d) 相間滑移矩方法

也應該注意到，與實驗結果相比，相間滑移矩方法仍還存在一些不足。實驗觀測表明液態水會向旋渦外圍聚集，形成類似“旋臂”結構，兩種數值方法都沒有模擬出這一結構。另一方面，實驗結果顯示，不僅在渦心，旋渦中整體的液相密度都較低，而數值結果中旋渦中的液相密度偏高。分析原因，可能是由于沒有考慮液滴碰并效應造成的。如果考慮液滴碰并，理論上液滴數密度將減少，平均半徑將增大，液滴相對氣流的跟隨性變弱，相間差距更加明顯。液滴將更快的離開旋渦，并且更容易在旋渦邊緣聚集。

4 總結與展望

本文首先介紹了矩方法的發展過程，包括傳統矩方法、求積矩方法和直接求積矩方法。傳統矩方法需要對分布函數的具體形式作出假設。求積矩方法克服了傳統矩方法液滴增長率必須是半徑線性關系的限制(這一限制實在內部坐標取液滴半徑條件下做出的)，并且不需要知道分布函數的具體形式。而直接求積矩方法中，雖然分布函數顯式給出，但是不需要對其具體形式作出額外假設。傳統矩方法，計算過程較為簡單，各項物理含義清楚明了。求積矩方法需要進行矩陣計算，計算過程較為復雜，但仍然追蹤各階矩。直接求積矩方法，直接追蹤高斯點和對應權值，可以反映相間、高斯點對應的內部坐標間的相互作用，這是前兩者所不具備的。

隨后，介紹了矩方法的擴展，包括推廣到異質凝結以及考慮相間相互作用。通過引入“瞬間活化”假設，簡化異質成核模型，僅需修改矩方程的源項，即可得到伴隨異質凝結的流動問題的控制方程。應用該方法對激波管中異質凝結問題進行了參數研究，加深了對異質凝結和流動相互作用的認識。為了描述氣、液相之間的相互作用，將氣液兩相分開處理，建立了一套考慮兩相間質量、動量和能量交換的凝結矩方法。應用此方法研究了伴隨同質凝結的旋渦運動問題，結果表明相間滑移矩方法能正確的捕捉到實驗觀測發現的“白點”現象，即凝結生成的液滴在離心力作用下離開旋渦中心，使得該區域液相密度極低。

近期，矩方法的發展主要在于具體問題的應用；彌補之前方法的不足；反推出分布函數的具體形式；處理非圓粒子等。今后，矩方法的發展還將主要集中在處理具體問題上，如：粒子的破碎、聚合，各相間的相互作用等方面。