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

(南昌大學(xué) 機(jī)電工程學(xué)院，南昌 330000)

0 概述

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

KHI的數(shù)值研究[4]一直是計(jì)算流體力學(xué)的經(jīng)典問題，其中相應(yīng)的數(shù)值研究方法包括線性區(qū)域和非線性區(qū)域。一般而言，用于研究線性穩(wěn)定性和非線性穩(wěn)定性的方法主要包括數(shù)值模擬、試驗(yàn)分析及理論研究。其中，文獻(xiàn)[4]中通過試驗(yàn)對(duì)KHI的非線性增長(zhǎng)進(jìn)行了探索，結(jié)果顯示很難通過試驗(yàn)對(duì)KHI的非線性增長(zhǎng)進(jìn)行研究。文獻(xiàn)[5]中對(duì)相關(guān)線性穩(wěn)定性進(jìn)行了解析研究。文獻(xiàn)[6]中使用渦流模型研究KHI，研究發(fā)現(xiàn)KHI在后期會(huì)產(chǎn)生二次不穩(wěn)定性。目前數(shù)值模擬已成為KHI數(shù)值研究中最常用的一類方法。

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

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

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

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

1 數(shù)值方法

1.1 控制方程

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

(1)

(2)

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

弗勞德準(zhǔn)數(shù)(Fr)是一個(gè)用于確定液體流動(dòng)狀態(tài)的無量綱數(shù)，F(xiàn)r>1的流動(dòng)稱為超臨界流動(dòng)。

(3)

式中，U3為下層流體水平速度；L為計(jì)算域的長(zhǎng)。

另外，對(duì)于多相流界面曲率較大的情況，本文使用韋伯?dāng)?shù)(We)來對(duì)流體運(yùn)動(dòng)進(jìn)行分析。韋伯?dāng)?shù)與慣性力之間呈正相關(guān)關(guān)系。

(4)

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

本文中使用雷諾數(shù)(Re)這一無量綱數(shù)來表示流體流動(dòng)狀態(tài)。

(5)

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

通過在方程(1)中對(duì)壓力項(xiàng)進(jìn)行分離，并引入瞬時(shí)速度得到方程(6)。

(6)

1.2 數(shù)值模擬方法

本文采用界面追蹤法(FTM)進(jìn)行數(shù)值模擬，并采用二階精度投影方法對(duì)N-S方程進(jìn)行求解。界面由相互連接的標(biāo)記點(diǎn)指示，該界面用于更新每個(gè)網(wǎng)格點(diǎn)的密度和黏度。

1.2.1 界面追蹤和標(biāo)記

本文將界面的演化與固定網(wǎng)格相結(jié)合，然后使用雙線性插值法在網(wǎng)格和界面之間進(jìn)行信息交換。式中的權(quán)重可以解釋為面積分?jǐn)?shù)(圖1)，且其數(shù)學(xué)表達(dá)式見式(7)。

(7)

圖1 權(quán)重系數(shù)值表示圖

(8)

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

(9)

1.2.2 表面張力

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

Fσ=σqn

(10)

本文對(duì)單位界面網(wǎng)格的表面張力進(jìn)行了分析。對(duì)于二維不混溶流動(dòng)，可得：

(11)

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

Fσ=σqnδ(x-xf)

(12)

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

(13)

(14)

式中，μ為流體的黏度。

2 問題描述及分析

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

圖2 KHI物理模型

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

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

(15)

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

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

(16)

數(shù)值增長(zhǎng)速率可以寫成式(17)形式。

(17)

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

(18)

(19)

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

2.1 物理模型驗(yàn)證

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

圖3 參考文獻(xiàn)[24]模擬結(jié)果與本文模擬結(jié)果對(duì)比圖

2.2 網(wǎng)格收斂研究

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

圖4 網(wǎng)格收斂研究

3 數(shù)值分析及其結(jié)果

3.1 中間流體層的影響

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

3.1.1 中間層流體密度的影響

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

圖5 卷起高度示意圖

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

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

圖7為不同密度的中間流體層在不同無量綱時(shí)間下的界面比較圖。由圖7(a)、圖7(b)可知，上界面KHI水平方向的發(fā)展隨著中間層密度的增加而減弱，而下界面KHI水平方向的發(fā)展隨中間層密度的增加而加強(qiáng)。這是因?yàn)橹虚g層流體密度的增加對(duì)上界面是起到穩(wěn)定作用，對(duì)下界面卻會(huì)加劇界面不穩(wěn)定性的發(fā)展。由圖7(c)可知，上下界面KHI向內(nèi)發(fā)展的程度均隨著中間流體密度的增加而增加。

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

3.1.2 中間層流體黏度的影響

在三組分不混溶流體體系中，中間流體的黏度會(huì)對(duì)整個(gè)體系的界面演化產(chǎn)生重要作用，因此本節(jié)中單獨(dú)考慮了不同黏度的中間流體的黏性流動(dòng)。

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

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

圖9為不同中間流體層黏度在不同無量綱時(shí)間下的界面比較圖。如圖9所示，隨著中間流體層黏度的增加，上下界面的卷起運(yùn)動(dòng)被抑制。因此，中間層黏度的增加能夠抑制上、下界面KHI的向上和向內(nèi)發(fā)展。

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

3.2 Fr的影響

重力是維持界面穩(wěn)定的重要因素之一，本文中對(duì)不同F(xiàn)r下黏性流體的界面演化情況進(jìn)行了模擬。

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

圖10 不同F(xiàn)r下卷起高度和增長(zhǎng)速率隨無量綱時(shí)間的變化

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

圖11 不同F(xiàn)r下在不同無量綱時(shí)間下的界面比較圖

3.3 表面張力和波數(shù)的影響

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

圖12 不同表面張力系數(shù)下的卷起高度和增長(zhǎng)速率隨無量綱 時(shí)間的變化

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

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

圖14 無量綱時(shí)間t=0.8時(shí)不同波數(shù)下的云圖

圖15 不同波數(shù)下的上下界面卷起高度和增長(zhǎng)速率隨無量綱時(shí)間t的變化

3.4 剪切力的影響

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

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

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

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

4 結(jié)論

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

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

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

(4) 當(dāng)波數(shù)大于4時(shí)，不會(huì)觸發(fā)典型形式的KHI的發(fā)生。

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