基于界面追蹤法的Kelvin-Helmholtz不穩定性直接數值模擬

過海龍，張 瑩，李文彬，李培生

(南昌大學 機電工程學院，南昌 330000)

0 概述

隨著國民經濟的快速發展，國家對生態環境的保護也越來越重視，對汽車尾氣排放的管控也越發嚴格，高效低污染的汽油機或柴油機成為內燃機目前發展的主要趨勢。其中高壓燃油噴射霧化技術對降低燃燒階段缸內顆粒物生成量是一種十分有效且重要的手段。在燃油噴射霧化[1]的過程中經常涉及到一種重要的流體動力學不穩定性：開爾文-亥姆霍茲不穩定性(Kelvin-Helmholtz instability， KHI)。文獻[2-3]中使用開爾文亥姆霍茲-瑞麗泰勒(Kelvin-Helmholtz and Rayleigh-Taylor, KH-RT)模型研究了燃油射流霧化的過程，霧化過程可大致分為液體尚未分裂和分裂后形成液滴最終霧化兩步，這兩步均會影響噴霧最終的形態，從而影響內燃機后期的燃燒過程。本文在K-H模型中研究射流尚未分裂時的運動，通過波高和線性增長率來研究液體尚未分裂的運動，為進一步研究液體的霧化過程提供一定理論基礎。

KHI的數值研究[4]一直是計算流體力學的經典問題，其中相應的數值研究方法包括線性區域和非線性區域。一般而言，用于研究線性穩定性和非線性穩定性的方法主要包括數值模擬、試驗分析及理論研究。其中，文獻[4]中通過試驗對KHI的非線性增長進行了探索，結果顯示很難通過試驗對KHI的非線性增長進行研究。文獻[5]中對相關線性穩定性進行了解析研究。文獻[6]中使用渦流模型研究KHI，研究發現KHI在后期會產生二次不穩定性。目前數值模擬已成為KHI數值研究中最常用的一類方法。

文獻[7]使用光滑粒子流體動力學方法(smoothed partical hydrodynamics， SPH)[8-11]對界面的發展進行了研究，結果顯示，KHI的增長速度與阿特伍德數之間存在正相關關系。另外，還在不同密度和不同表面張力的條件下進行了模擬并分析其對KHI的影響。

文獻[12]中使用格子玻爾茲曼方法(lattice Boltzmann method， LBM)[12-14]研究了二維KHI，其針對表面張力對二維混合層中演化的影響的研究結果表明：部分動能可以轉化為勢能。

文獻[15]中通過使用氣泡多相直接數值模擬(direct numerical simulation， DNS)對流動模型的影響因子進行了探究。研究結果顯示，該模型能夠較為準確地模擬出界面在不同初始條件下的變化。該方法具有計算量小、能快速完成數值模擬的優勢，但不能準確捕捉界面的運動情況。界面追蹤法(front tracking method， FTM)[16-22]可以準確地捕捉界面的移動，文獻[23]中于2001年通過FTM成功模擬了三維流動、二維流動和軸對稱流動。

準確跟蹤界面的變化對研究KHI的發展非常重要。本文中采用FTM從界面向上、向內及水平的發展三方面對KHI的發展進行研究。為分析和觀察界面向上的發展，本文使用卷起高度h和線性增長率γ作為參數來衡量界面的向上發展，通過界面比較圖來比較界面的橫向和向內發展。利用FTM對三組分不混溶流體的二維KHI進行了數值模擬，分別研究了中間流體層密度、黏度、重力、波數、剪切力對上、下界面KHI發展方面的影響。

1 數值方法

1.1 控制方程

對于二維不可壓縮流動，其納維-斯托克斯(Navier-Stokes, N-S)方程可表示為式(1)和式(2)。

(1)

(2)

式中，u為流場中的速度矢量；g為重力加速度；q為曲率；p為壓強；ρ和ν為非連續的密度場和黏度場；n為表面法向量；σ為界面上的表面張力系數；δ為狄拉克函數，并且其只有在界面上時才不為0；t為時間。

弗勞德準數(Fr)是一個用于確定液體流動狀態的無量綱數，Fr>1的流動稱為超臨界流動。

(3)

式中，U3為下層流體水平速度；L為計算域的長。

另外，對于多相流界面曲率較大的情況，本文使用韋伯數(We)來對流體運動進行分析。韋伯數與慣性力之間呈正相關關系。

(4)

式中，ρ3為下層流體密度。

本文中使用雷諾數(Re)這一無量綱數來表示流體流動狀態。

(5)

式中，v3為下層流體豎直方向速度。

通過在方程(1)中對壓力項進行分離，并引入瞬時速度得到方程(6)。

(6)

1.2 數值模擬方法

本文采用界面追蹤法(FTM)進行數值模擬，并采用二階精度投影方法對N-S方程進行求解。界面由相互連接的標記點指示，該界面用于更新每個網格點的密度和黏度。

1.2.1 界面追蹤和標記

本文將界面的演化與固定網格相結合，然后使用雙線性插值法在網格和界面之間進行信息交換。式中的權重可以解釋為面積分數(圖1)，且其數學表達式見式(7)。

(7)

圖1 權重系數值表示圖

(8)

本文中使用簡單的顯式一階時間積分，并通過式(9)獲得下一個時刻的界面位置。

(9)

1.2.2 表面張力

表面張力是模擬不混溶多相流的重要參數。其單位面積在界面上的力可以通過式(10)得到。

Fσ=σqn

(10)

本文對單位界面網格的表面張力進行了分析。對于二維不混溶流動，可得：

(11)

式中，s為界面元素；τ為切向量。由此，界面上的表面張力方程可以寫成式(12)形式。

Fσ=σqnδ(x-xf)

(12)

式中，δ(x-xf) 僅在處于界面上時為“1”，否則均為“0”。然后，可以得到單元界面網格上的表面張力表達式如式(13)所示。

(13)

(14)

式中，μ為流體的黏度。

2 問題描述及分析

本文采用長為L、高為H的計算域，如圖2所示。為簡單起見，該計算域的高度等于其長度。首先將周期邊界條件應用于計算域的左右兩側，并將固壁邊界條件施加于該計算域的上下兩側。在該模型中，使用下層流體的密度、速度及計算域的長度作為無量綱尺度來模擬研究界面演化。

圖2 KHI物理模型

該流場中的密度分布為ρ1<ρ2<ρ3(下標1、2、3分別表示上、中、下層流體)。相鄰流層流體的水平運動方向相反，初始時水平速度Un(n=1,2,3)的絕對值均為1。在上述所有模擬中，U1=U3=1，U2=-1。然后將數值模擬的垂直速度Vn(n=1,2,3)設為0。

在該數值模擬中，在初始時刻對界面施加正弦擾動。定義函數f(x,y;y0,ε)如式(15)所示。

(15)

式中，ε為速度梯度的厚度；y0為界面的初始位置。

將兩個形式為A=A0eikx+ω t的正弦擾動作用于界面。式中，A為擾動幅度；ω為波的角頻率；k為波數；A0為初始擾動幅度。當既有表面張力又存在重力時，該特征值條件可由公式(16)給出。

(16)

數值增長速率可以寫成式(17)形式。

(17)

由于理查森數(Ri)的定義是浮力與剪切力之比，因此Ri可以表示為式(18)和式(19)形式。

(18)

(19)

式中，Rit為上界面的Ri;Rib為下界面的Ri;σt為上界面表面張力系數;σb為下界面表面張力系數。

2.1 物理模型驗證

為驗證FTM的有效性，本文將不同初始擾動下得到的界面波高結果與文獻[24]得到的波高進行了比較。本文所用的模擬條件均與參考文獻[24]條件相同,即：Rit=0.047，Rib=0.134，A0分別為0.01、0.02及0.05。圖3為本文和參考文獻[24]在不同界面初始擾動下界面卷起高度隨無量綱時間的變化，可以看出本文的模擬結果與參考文獻中的分析結果高度一致。

圖3 參考文獻[24]模擬結果與本文模擬結果對比圖

2.2 網格收斂研究

合適的網格對于提高模擬準確性及研究效率均具有重要意義。為此，本文中針對三組分流體的網格收斂問題進行了分析，并對其進行了數值模擬試驗，以期找出最合適的網格數量。所有模擬計算均是基于64×64、128×128、256×256、512×512的網格劃分情況下進行，且模擬初始參數均為Rit=0.057，Rib=0.142，A0=0.02。通過不同網格數下界面的比較圖來選出最適合的網格密度。圖4為不同網格密度在無量綱時間為0.5時的界面比較圖。從圖4中可以明顯地看出，一旦網格數增加到128×128 或以上，界面明顯趨向于合并成如圖4所示的單一曲線，因此選擇128×128的網格進行后續研究。

圖4 網格收斂研究

3 數值分析及其結果

3.1 中間流體層的影響

處于中間位置的流體同時與上層流體和下層流體相接觸，其在三組分不混溶流體的運動中起著重要作用，因此本文中針對中間層流體密度和黏度對三組分不混溶流體KHI的影響分別進行了模擬。

3.1.1 中間層流體密度的影響

眾所周知，KHI通?？梢栽?個方向上發展，即向上發展、向內發展和橫向發展。其中，向內和橫向發展的KHI現象較為明顯。為了能更好地研究界面向上的發展，將對上下界面的卷起高度(h)和增長速率(γ)進行分析研究。圖5為卷起高度示意圖。

圖5 卷起高度示意圖

不同無量綱中間流體密度時上下界面卷起高度和增長速率隨無量綱時間的變化如圖6所示。隨著時間的推移，上下界面的卷起高度會逐漸增大，且卷起高度和增長速率均表現為先增大后減小，這表明界面的卷起和向上發展均有一個上限。除此之外，本研究中還發現，下界面卷起的高度會隨中間層密度的增加而增加。這是因為中間層的慣性力增大，使得下界面變得越來越不穩定，而這一作用與動能的增加具有相同的效果。如圖6(a)所示，在界面發展的后期，在Ri較大的情況下的卷起高度要比Ri較小時大。這是因為與Ri較小的情況相比，Ri較大的情況下卷起高度變化較慢。這表明Ri的增大不僅能夠抑制界面的演化，而且具有相對較大的勢能還能夠起到抑制界面退化的作用。

圖6 不同中間流體密度時上下界面卷起高度和增長速率 隨無量綱時間的變化

圖7為不同密度的中間流體層在不同無量綱時間下的界面比較圖。由圖7(a)、圖7(b)可知，上界面KHI水平方向的發展隨著中間層密度的增加而減弱，而下界面KHI水平方向的發展隨中間層密度的增加而加強。這是因為中間層流體密度的增加對上界面是起到穩定作用，對下界面卻會加劇界面不穩定性的發展。由圖7(c)可知，上下界面KHI向內發展的程度均隨著中間流體密度的增加而增加。

圖7 不同密度的中間流體層在不同無量綱時間下的 界面比較圖

3.1.2 中間層流體黏度的影響

在三組分不混溶流體體系中，中間流體的黏度會對整個體系的界面演化產生重要作用，因此本節中單獨考慮了不同黏度的中間流體的黏性流動。

圖8為中間層流體具有不同無量綱黏度時上下界面的卷起高度和增長速率隨無量綱時間的變化。從圖中可以看出，所有情況下的卷起高度和增長速率都呈現出先增大后減小的趨勢，并且隨著黏滯力的增大，卷起高度和增長速率都會減少，這表明黏滯力主要抑制界面的變化。可以看出，增加中間層動態黏度的效果相當于由于動能降低而增加Ri，即中間流體層的黏度對上下界面的卷起高度具有抑制作用。

圖8 中間層流體具有不同黏度時上下界面的卷起高度和 增長速率隨無量綱時間的變化

圖9為不同中間流體層黏度在不同無量綱時間下的界面比較圖。如圖9所示，隨著中間流體層黏度的增加，上下界面的卷起運動被抑制。因此，中間層黏度的增加能夠抑制上、下界面KHI的向上和向內發展。

圖9 不同中間流體層黏度在不同無量綱時間下的 界面比較圖

3.2 Fr的影響

重力是維持界面穩定的重要因素之一，本文中對不同Fr下黏性流體的界面演化情況進行了模擬。

圖10為不同Fr時上下界面卷起高度和增長速率隨無量綱時間的變化。從圖10中可以看出，上下界面的卷起高度隨Fr的減小呈減小趨勢，不同模擬條件下的增長率具有相同的變化趨勢。具體而言，當t在0.2～0.9的范圍內時，上下界面增長速率會先增加，并在t=0.5時達到峰值然后減小，并且上界面Ri越大，上界面的卷起高度就越小，這表明表面張力和重力在界面演化過程中起著穩定的作用。

圖10 不同Fr下卷起高度和增長速率隨無量綱時間的變化

圖11為不同Fr時在不同無量綱時間下的界面比較圖。如圖11所示隨著Fr的減小，上界面的向內演化受到抑制。然而，由于底層流體的流體黏度較大，因此不同Fr時下界面向內卷曲程度的差別不大。

圖11 不同Fr下在不同無量綱時間下的界面比較圖

3.3 表面張力和波數的影響

重力和表面張力都是維持界面穩定的主要因素。表面張力的方向與界面相切，能夠抑制上下界面的向上和向內演化。不同無量綱表面張力系數下的上下界面卷起高度和增長速率隨無量綱時間的變化圖12所示。由圖12可見，上下界面各自的卷起高度曲線和上界面增長速率曲線隨著表面張力的改變卻幾乎完全吻合，這表明當波數足夠小時，重力在界面向上演化中能夠起到更為重要的作用。

圖12 不同表面張力系數下的卷起高度和增長速率隨無量綱 時間的變化

圖13為不同表面張力系數下不同無量綱時間的界面比較圖。從圖13中可以看出，隨著表面張力的增加，界面向內演化會受到抑制。根據理論KHI公式可以看出，重力對界面演化的影響會隨著波數的增加而減小。這意味著波數越多，表面張力對KHI的影響就越明顯。

圖14為在不同波數下無量綱時間為0.8的云圖。通過對比云圖可以看出，界面隨著波數的增加而拉長。圖15為不同波數下的上下界面卷起高度和增長速率隨無量綱時間的變化。如圖15所示，當波數為4時，界面卷起高度較小，卷起高度和增長速率隨波數的增加而提前達到峰值。

圖14 無量綱時間t=0.8時不同波數下的云圖

圖15 不同波數下的上下界面卷起高度和增長速率隨無量綱時間t的變化

3.4 剪切力的影響

當界面上的剪切作用能夠抵消來自重力和表面張力的穩定作用時，流體界面便會發生卷曲，即剪切力是引發KHI的主要因素。本文中分別就一些動能相同但剪切力不同的情況進行了模擬，且通過改變界面厚度來改變剪切力。

圖16為不同無量綱剪切力作用下上下界面卷起高度和增長速率隨無量綱時間的變化。從圖中可以看出，隨著厚度的增加，界面的卷起高度和增長速率均呈減小趨勢。圖17為不同速度梯度厚度在不同無量綱時間的界面比較圖。從圖中可以看出，隨著剪切力的減小，界面向內演化也會受到一定的抑制。綜上，KHI的發展受到三個方向的抑制。當界面厚度大于擾動波長的0.03倍時，界面后期不會演化成貓眼狀。

圖16 不同剪切力作用下上下界面卷起高度和增長速率隨無量綱時間的變化

圖17 不同速度梯度厚度在不同無量綱時間的界面比較圖

4 結論

(1) 三組分不混溶流體中間流體層密度越大，下界面的卷起高度和增長速率越大，而上界面水平方向的KHI發展越快。Ri的增大不僅能夠抑制界面的演化，而且還能夠限制界面的退化。

(2) 增加三組分不混溶流體相互接觸兩層流體中的一層流體的黏度能夠抑制界面的向內演化，并且三組分不混溶流體中間流體層的黏度越小，界面的卷起高度和增大速率就越大。

(3)Fr與卷起高度、數值增長速率之間呈正相關關系。隨著表面張力系數的增大，界面向內演化將受到抑制。

(4) 當波數大于4時，不會觸發典型形式的KHI的發生。

(5) 剪切力越大，界面的卷起高度和增長速率就越大，且剪切力對KHI有著顯著的影響。