可壓縮流動和可壓縮多相流的數(shù)值方法研究領(lǐng)域十分廣闊，近幾十年來取得了重大成就。然而，由于它適用問題廣泛，所以目前仍舊是一個熱門的研究話題。最通用的兩相可壓縮流模型稱為7方程模型(1D)，該模型考慮了每相的密度、動量、內(nèi)能、壓力、質(zhì)量分?jǐn)?shù)(10個未知數(shù))。每個相的狀態(tài)方程與ρ, p, e和質(zhì)量分?jǐn)?shù)的約束有關(guān)，留下了7個未知數(shù)。7方程模型可被進(jìn)一步簡化為考慮單一壓力的6方程模型以及僅考慮混合速度的5方程模型?；趬毫Φ目蓧嚎s流方法與基于密度的可壓縮流方法的二分法是眾所周知的，它們各有優(yōu)缺點(diǎn)。

Saurel和Lemetayer提出了針對界面、激波和傳質(zhì)問題的多相流總體框架。本文提出一種壓力基方法對全局連續(xù)性方程進(jìn)行非守恒處理以導(dǎo)出壓力方程。

?

一、數(shù)學(xué)方程

本文提出了一種假設(shè)單一壓力、速度和溫度的四方程可壓縮多相流模型。該模型可通過求解單獨(dú)的溫度方程擴(kuò)展為5方程模型，并使用代數(shù)滑移模型擴(kuò)展為6方程版本。

?

集合平均模型

以體積分?jǐn)?shù)(αk)表示的相質(zhì)量守恒方程為:

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

這里引入混合速度平流的物質(zhì)導(dǎo)數(shù)為：

??

其中ψ是任何平流量，上標(biāo)m表示平流速度是混合速度。

質(zhì)量守恒方程可以用質(zhì)量分?jǐn)?shù)(Yk)以物質(zhì)導(dǎo)數(shù)形式重寫為:

??

?

界面跟蹤

對于界面跟蹤方法，只考慮由level-set函數(shù)表示的兩相流系統(tǒng)。沒有解出相質(zhì)量守恒方程。Level-set方程的演化控制方程如下：

??

其中φ為到界面的距離，u為界面的速度。平滑Heaviside函數(shù)Hk依賴于距離函數(shù)Hk = Hk(φ)。Hk的物質(zhì)導(dǎo)數(shù)為：

?

其中δk(φ)為選擇界面的光滑狄拉克函數(shù)。

?

Mixture動量方程

Mixture動量方程如下：

??

其中一個關(guān)鍵假設(shè)是混合粘度模型，該模型考慮了由于滑移引起的粘性動量傳遞。對于界面跟蹤方法，漂移速度項(xiàng)設(shè)為零。

?

Mixture溫度方程

通過定義顯混合焓推導(dǎo)出混合溫度方程，其中T0為共同參考溫度，Cp為混合熱容量(假設(shè)相熱容恒定)。

??

式中Φ為動能的粘性耗散，最后一項(xiàng)為P個不同分量發(fā)生相變時的潛熱貢獻(xiàn)，而混合速度對T平流的物質(zhì)導(dǎo)數(shù)可以簡化為:

??

?

系綜平均系統(tǒng)壓力方程的發(fā)展

導(dǎo)出可壓縮多相流壓力方程的起點(diǎn)是非守恒形式的Mixture連續(xù)方程：

??

取ρ = ρ(p, T, Yk)我們可以將上述方程展開為基于獨(dú)立未知量的線性化，如下：

??

為了評估約束導(dǎo)數(shù)，我們從相對于相密度和質(zhì)量分?jǐn)?shù)的密度定義開始：

??

式中ρk = ρk(p, T)。在恒定(T, Yk)下求ρ對p的導(dǎo)數(shù)：

??

這是等溫混合聲速。同樣地，常數(shù)(p, Yk)處的導(dǎo)數(shù)可以寫成：

??

恒定(p, T)時密度對Yk的導(dǎo)數(shù)為：

?

??

?

將方程(12)-(14)、方程(8)、方程(3)代入方程(10)，得到：

??

這是n相可壓縮流壓力方程的最終形式。

?

單相流的限制

對于可壓縮的單相流，方程(15)可簡化為：

??

由于理想氣體滿足狀態(tài)方程p = ρRT，所以可以導(dǎo)出：

??

利用Cp?Cv = R和Cp/Cv = γ的關(guān)系。對于理想氣體等熵流動，將方程簡化為聲波方程。

?

不可壓縮多相流的限制

在不可壓縮流動的極限下，由于ρk不是(p, T)的函數(shù)，所以可得：

????

在不存在相間滑移和相變的情況下，這給出了速度場的無散度條件。

?

界面跟蹤的壓力方程

對于level-set界面跟蹤方法，我們使用Hk作為主要變量，線性化連續(xù)性方程：

??

利用方程(5)的定義，可得：

??

方程(20)中的最后一項(xiàng)為由于相變導(dǎo)致的源項(xiàng)。需要注意的是，流體體積模量的體積加權(quán)求和并不能得到由系綜平均方程得到的Wallis聲速。

利用這一公式，提出一種處理不可壓縮和可壓縮多相流的綜合平均和界面跟蹤的單一方法。非守恒方程的優(yōu)點(diǎn)在于它是不可壓縮level-set方法的自然推廣，而不可壓縮level-set方法本身就是非守恒的。這種方法不受接觸不連續(xù)點(diǎn)附近偽壓力振蕩的影響。方程(15)和(20)的右側(cè)表示不可壓縮項(xiàng)對方程的貢獻(xiàn)，而可壓縮項(xiàng)對方程貢獻(xiàn)置于左側(cè)。此外，該方法還提出了流體的正確體積模量(恒壓和恒溫的組合)，可用于使用等溫或等熵導(dǎo)數(shù)的壓力基方法。其他現(xiàn)象如體加熱和相變的影響也被考慮在內(nèi)。

?

二、數(shù)值離散化

該模型被實(shí)現(xiàn)到積鼎VirtualFlow軟件中。VirtualFlow是一個專門用于多相流建模的有限體積CFD求解器。除壓力方程采用一階迎風(fēng)外，所有方程的對流項(xiàng)均采用HLPA格式離散化，對非定常問題采用一階歐拉時間步進(jìn)。

?

三、驗(yàn) 證

單相流Joukowsky問題

這里給出閥門突然關(guān)閉的Joukowsky問題。問題由(圖1)一個充滿長度為L = 1m的水域組成，該區(qū)域具有固定的壓力邊界，壓力為P0 = 1bar。初始化流速為V0 = 0.0666m/s。根據(jù)Modified Tait狀態(tài)方程，水的密度為ρ = 1000 kg/(m^3)，聲速為c = 1481.47m/s。根據(jù)理論解，期望壓力比為(P0 +ρCV0)/P0 = 2。

??

圖1:Joukowsky問題的草圖

?

區(qū)域封閉端(閥)處的壓力演化如圖2所示，并與解析解進(jìn)行了對比。在這種情況下，計算域被離散為100個控制體。壓力峰值預(yù)測準(zhǔn)確，但隨著時間的推移，由于數(shù)值擴(kuò)散，波形變得模糊。精細(xì)化網(wǎng)格(200和400單元)的結(jié)果如圖3所示，表明它們逐步收斂。

??

圖2:閥門處壓力演化對比圖

??

圖3:隨著網(wǎng)格細(xì)化，最后周期內(nèi)閥門壓力變化情況

?

多相流

這里介紹了該方法在界面流和混合流中的應(yīng)用。在這兩種情況下，結(jié)果已與完善的解析解進(jìn)行了對比。

?

兩相混合物中的聲速

為了驗(yàn)證該模型對不同氣液混合物的聲速預(yù)測是否正確，我們初始化了一個長度為0.1 m的一維計算域，其中氣體的體積分?jǐn)?shù)分別為1、0.97、0.5、0.1、0，范圍從純氣體到純液體。在計算域的一端設(shè)置一個振幅為1000 Pa的初始壓力脈沖，測量信號到達(dá)另一端的時間。使用64個控制體對計算域進(jìn)行離散化。氣相被認(rèn)為是蒸汽，并被建模為理想氣體，其分子量為18 g/mol，熱容比為1.2675。用Modified Tait狀態(tài)方程對液體進(jìn)行建模。系統(tǒng)壓力設(shè)置為1bar，溫度初始化為400k。

多相混合物中的聲速由Wallis聲速給出，為:

??

在所有情況下，遠(yuǎn)端的壓力信號如圖4所示。到達(dá)速度最快的是純水，其次是純空氣，根據(jù)方程(21)的預(yù)測，之后依次到達(dá)的是體積分?jǐn)?shù)分別為0.97、0.1和0.5的工況。將數(shù)值結(jié)果計算的聲速與表1中的解析值進(jìn)行比較。結(jié)果與分析值吻合良好(在3%以內(nèi))。

??

圖4:到達(dá)遠(yuǎn)端的壓力信號

?

表1:混合聲速分析結(jié)果對比

??

?

?

本文模擬了入射壓力波在空氣-水界面上的傳播和反射，即:(a)壓力波從空氣側(cè)入射，(b)壓力波從水側(cè)到達(dá)。透射和反射系數(shù)由聲阻抗給出:

??

其中z = ρc為聲阻抗，下標(biāo)l和r表示左(入射)和右流體。圖5給出了問題描述，其中空氣被建模為分子量為28 g/mol，熱容量比為1.4的理想氣體，水的模型使用Modified Tait狀態(tài)方程。

??

圖5:聲波透射與反射問題示意圖

區(qū)域尺寸設(shè)為1 m，空氣-水界面位于中心。系統(tǒng)壓力設(shè)置為1bar，溫度初始化為300k。振幅為1000 Pa的壓力脈沖初始化在距離空氣或水相界面0.45 m處。本算例采用level-set界面跟蹤方法。

解析透射和反射系數(shù)如表2所示，并給出了相應(yīng)的數(shù)值結(jié)果。這個問題很有趣，因?yàn)樵诳諝?水的情況下，反射波的振幅與入射波的振幅相同，而透射波的振幅被放大了2倍。在波來自高阻抗流體(水-空氣)的情況下，反射波改變了相位，幾乎沒有傳輸任何東西?？梢姡蓧嚎s界面流動模型能夠較準(zhǔn)確地再現(xiàn)這種效應(yīng)。不同波的壓力信號如圖6所示。

圖6:入射波、反射波和透射波的壓力信號，上圖為空氣-水，下圖為水-空氣。

?

表2:解析結(jié)果與數(shù)值結(jié)果的比較

??

?

四、應(yīng) 用

飽和系統(tǒng)的壓力傳播

上一節(jié)提出的問題被擴(kuò)展到包括1bar時蒸汽和水之間的傳質(zhì)。溫度初始化為飽和溫度372.756 K。潛熱規(guī)定為2.2 MJ。選擇蒸汽初始體積分?jǐn)?shù)= 0.5的情況。采用基于Ranz-Marshall相關(guān)的分散相傳質(zhì)模型。傳質(zhì)速率(kg/(m^3))為:

??

本文中的數(shù)值模擬不考慮相間滑移。氣泡/液滴直徑通過直接輸入來確定。通過指定更低的dp值，系統(tǒng)可以被迫進(jìn)入相平衡狀態(tài)。這里模擬了三個不同的工況，即20毫米、2毫米和20微米。

?

圖7:蒸汽-水飽和系統(tǒng)聲速變化

?

飽和系統(tǒng)的聲速如圖7所示?？梢钥闯觯瑢τ? bar的系統(tǒng)，聲速可以低至1 m/s。這里模擬的情況用紅點(diǎn)表示。

三種情況下遠(yuǎn)端壓力信號如圖8所示。對于較大的氣泡/液滴直徑，聲速幾乎不受影響(≈21 m/s)，接近圖7所示的數(shù)據(jù)。當(dāng)dp = 20微米時，壓力信號衰減明顯，到達(dá)遠(yuǎn)端需要更長的時間。聲速估計為1.96 m/s，這與圖7中水-蒸汽熱力學(xué)平衡的數(shù)據(jù)非常接近。

??

圖8：飽和蒸汽-水系統(tǒng)到達(dá)計算域遠(yuǎn)端的壓力信號

?

噴嘴的臨界流動

作為可壓縮多相流模型的另一個應(yīng)用，使用上一節(jié)描述的傳質(zhì)模型(方程(24))來解決閥門和噴嘴中的臨界流動問題。二維平面噴管示意圖如圖9所示。

??

圖9:噴嘴幾何形狀

?

入口條件：p0 = 3mpa, T0 = 507 K，蒸汽體積分?jǐn)?shù)為0.1的飽和蒸汽-水。通過降低出口壓力，找到臨界流動條件。采用無滑移的可壓縮均勻混合模型，模擬了進(jìn)出口壓力比分別為0.7、0.5、0.4和0.3的情況。傳質(zhì)模型的氣泡直徑設(shè)置為100微米。根據(jù)水的Antoine方程，計算飽和溫度作為壓力的函數(shù)。利用穩(wěn)態(tài)解算器進(jìn)行仿真，使其收斂。

圖10顯示了噴嘴內(nèi)蒸汽體積分?jǐn)?shù)的變化情況。Flashing導(dǎo)致體積分?jǐn)?shù)從入口的0.1增加到噴嘴阻塞前的0.7左右。在多相流情況下，聲速非線性特性極強(qiáng)的，所謂的聲線(黑色輪廓線)并不一定出現(xiàn)在喉部，如圖10中較低的兩個壓強(qiáng)所示。

圖10:不同出口壓力比下的空隙率等高線

?

確認(rèn)臨界流量的唯一方法是觀察質(zhì)量流量，并驗(yàn)證它們不會隨著出口壓力的降低而進(jìn)一步變化。由表3所示的質(zhì)量流量可以看出，在壓力比為0.3時，流量是臨界的。橫截面積的寬度為1cm。

表3:通過噴嘴的質(zhì)量流量。

??

?

?

結(jié) 論

本文是由積鼎科技提出的一種非保守可壓縮多相流模型，該模型是不可壓縮多相流模型的自然推廣，適用于界面跟蹤和系綜平均兩種方法。用這種單壓力、單溫度模型精確模擬了幾個問題。盡管它比標(biāo)準(zhǔn)的五方程模型更有局限性，但它仍然可以應(yīng)用于包括沖擊和相變在內(nèi)的廣泛問題。