基于FTM方法的二維Kelvin-Helmholtz不穩(wěn)定性數(shù)值模擬

， ， ， ，

(南昌大學(xué) 機電工程學(xué)院，南昌 330031)

1 引 言

K-H(Kelvin-Helmholtz instability)[1]不穩(wěn)定性是自由剪切流的無粘不穩(wěn)定性。一般認(rèn)為，不互溶的混合流體中，上層流體與下層流體間的速度差異足夠大時，界面上存在的微小擾動會使界面變形，出現(xiàn)Kelvin-Helmholtz不穩(wěn)定性。

K-H不穩(wěn)定性的數(shù)值研究一直是計算流體力學(xué)的經(jīng)典問題。晶格-玻爾茲曼方法LBM[2-5]LBM(Lattice -Boltzmann Method)能通過耦合分子間的相互作用正確地模擬界面的運動，其與界面運動和瑞利-泰勒不穩(wěn)定性[6]的研究類似，對于不混溶的兩相混合層的研究，LBM多相流模擬有良好的前景。

Lee等[7]用具有高分辨率和高精度的相場法(Phase-field Method)模擬多組分流體的K-H不穩(wěn)定性，并研究了韋伯?dāng)?shù)、密度比和速度差對K-H不穩(wěn)定性的影響。相場法的優(yōu)點在于無需進行顯性界面追蹤，并且能自動處理多相物理界面拓?fù)浣Y(jié)構(gòu)的改變，計算的時間步長也可以更大，在處理多組分流體時有獨特優(yōu)勢，然而傳統(tǒng)相場法計算的精度卻有待提高。

Shadloo等[8]用平滑顆粒流體動力學(xué)(SPH)方法模擬K-H不穩(wěn)定性，并研究了理查德森數(shù)和密度比的影響，發(fā)現(xiàn)K-H不穩(wěn)定性的增長主要受理查德森數(shù)的影響，而不僅是穩(wěn)定力的作用。Tryggvason等[9]利用渦片單元法[10-12]研究了相對薄的渦量層K-H不穩(wěn)定性的大振幅階段。在對非常薄的渦量層K-H不穩(wěn)定性的大振幅演化求解中，將無粘性、正規(guī)化的渦片單元模型和全粘性納維斯托克斯求解方案進行比較。結(jié)果表明，非粘性模型的零正規(guī)化標(biāo)度的極限、極限的高雷諾數(shù)和較小的初始厚度的粘性計算基本一致。

FTM[13-14]通過顯示跟蹤相分界面，具有以高精度高分辨率捕捉界面結(jié)構(gòu)拓?fù)渥兓膬?yōu)點。本文利用該優(yōu)點模擬了二維下流體的K-H不穩(wěn)定性，通過選取不同的速度梯度層厚度、速度差、邊界條件和理查德森數(shù)研究了K-H不穩(wěn)定性的界面形態(tài)演化的過程。此外，將Neumann邊界條件和Dirichlet邊界條件下的渦量分布進行對比，從而判斷兩種邊界條件對渦量分布的影響程度。

2 數(shù)值方法和數(shù)學(xué)模型

2.1 控制方程

使用顯式的有限差分方法求解控制方程，二維不可壓縮流動的動量方程為

σκδ(x-xf)

(1)

式中xf為界面位置,ρ和μ分別為非連續(xù)的密度場和黏度場,κ為平均曲率,δ只有在界面上時才不為0,其他參數(shù)遵循習(xí)慣約定。

2.2 密度場和黏度場重構(gòu)

不互溶、不可壓縮的流體在界面兩側(cè)保持各自的特性，因此在界面上存在物性階躍。FTM算法中，通過賦予界面一定厚度的過渡區(qū)，讓流體的特性在這個區(qū)域內(nèi)實現(xiàn)平滑過渡。當(dāng)界面移動時，界面上的密度和黏度分布隨之發(fā)生變化，故引入Heaviside函數(shù)[15]來表征這種變化，

(2)

(3)

(4)

2.3 界面追蹤

界面的運動通常是與固定網(wǎng)格結(jié)合在一起，采用雙線性插值法實現(xiàn)界面與網(wǎng)格間的信息交流，其表達(dá)式為

(5)

2.4 表面張力

界面上表面張力的計算表達(dá)式為

Fσ=σκnδ(x-xf)

(6)

式中δ(x-xf)為狄拉克函數(shù)，在界面上為1，其余都為0。

取單位面元上的表面張力為研究對象，對于二維流動有

κn=?t/?s

(7)

因此得到單位界面元上的表面張力為

(8)

式中t為界面上的切向量，s為界面元，l為界面上的節(jié)點。

2.5 數(shù)值模型

h(x,y;a,y0,ampl,ε)= tanh{[y-y0-

ampl*sin(aπx)]/ε}

(9)

擾動形式為由函數(shù)A=A0e- i k x +n t賦予界面一個初始擾動，其中k為波數(shù)，t為時間，當(dāng)考慮重力和表面張力時，根據(jù)文獻[1]的結(jié)果，可表示為

(10)

式中 若等號右側(cè)第二項的結(jié)果為非負(fù)數(shù)，擾動穩(wěn)定，否則，擾動將會隨時間而增長。當(dāng)擾動不穩(wěn)定時，等號右側(cè)第二項可以寫為

(11)

經(jīng)過簡單的數(shù)學(xué)變換，可以得到式(12)，Ri為理查德森數(shù)，在物理學(xué)中，理查德森數(shù)用于表征流體勢能與其動能的比值，勢能包含了重力項和表面張力項。

(12)

圖1 流體布局示意圖

Fig.1 Schematic on fluid layout

(13)

(14)

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

3.1速度梯度層厚度的影響

在不考慮速度梯度層的厚度時，上下兩種流體的速度相等，即|u1|=|u2|，實際過程中界面交界處存在速度的過渡區(qū)，從而會產(chǎn)生速度梯度層，速度梯度層的厚度會影響K-H不穩(wěn)定性的發(fā)展，如圖2所示。ε=0.10時，受速度梯度層的影響，界面向內(nèi)發(fā)展的趨勢明顯受到抑制，更多的是向右發(fā)展，渦量分布也比較散。隨著ε減小，速度梯度層變薄，對K-H不穩(wěn)定性的影響減弱，K-H不穩(wěn)定性發(fā)生的時間也提前。當(dāng)ε≤0.04時，界面形態(tài)和渦量分布已經(jīng)非常相似，由速度梯度層厚度造成的影響基本可以忽略。因此本文在余下的數(shù)值模擬中，均選取ε=0.02時的速度梯度層厚度為基準(zhǔn)。

3.2 速度差的影響

在不考慮切向剪切力的情況下，界面平衡時，剪切力和表面張力滿足

μ(?u/?y)=σκ

(15)

式中μ為流體粘度,κ為表面界面曲率，在本節(jié)中初始條件為

u(x,y)=h(x,y;2,0.5,0.01,0.02)

v(x,y)=0

(16)

不考慮梯度層的厚度，上下兩種流體的速度大小相同，使用2×|u|來衡量速度梯度的大小。從式(14)可知，當(dāng)界面處的速度梯度不足以破壞表面張力維持的穩(wěn)定時，界面的穩(wěn)定將不會打破。圖3分別為t=2.0 s時刻速度梯度為2×|0.4u|， 2×|0.8u|， 2×|1.2u|和2×|1.6u|四種情況下界面的形態(tài)(a,b,c和d)和流場渦量(a1，b1，c1和d1)分布?？梢钥闯觯俣炔钤酱?，界面上的剪切力越大，當(dāng)剪切力的作用大于表面張力和重力的作用時，界面上的擾動就會發(fā)展為K-H不穩(wěn)定性。對比圖3的四種情況可以發(fā)現(xiàn)，界面的不穩(wěn)定發(fā)生的時間會因為速度梯度的增加而提前，并且同一時刻的界面變形及渦量的擴散強度均會隨著速度梯度的增大而增大。

3.3 邊界條件的影響

對比兩種邊界條件的影響，得到圖4的界面形態(tài)的演化和圖5渦量隨時間的變化，兩種邊界條件的初始條件均如式(16)所示。

K-H不穩(wěn)定性是由界面上的擾動引起的，一般有三個發(fā)展方向，分別為向上、向內(nèi)和向右(由U1的速度方向決定)，而最明顯的是向內(nèi)發(fā)展，圖4 的結(jié)果也印證了這點。兩種邊界條件下，盡管初始條件相同，同一時刻Neumann邊界下界面的變形和卷曲程度明顯比Dirichlet邊界下的劇烈，隨著界面上擾動的增長，這個現(xiàn)象愈加明顯。Neumann邊界條件下，K-H不穩(wěn)定性非常明顯，表現(xiàn)為界面向內(nèi)卷的圈數(shù)更多，卷起的液膜更薄，而在Dirichlet邊內(nèi)卷起的圈數(shù)明顯減少，卷起的液膜也變厚。受邊界條件的影響，抑制最明顯的是界面向內(nèi)發(fā)展的過程，然而界面形態(tài)的演化過程卻非常相似，這表明邊界條件并不能完全抑制K-H不穩(wěn)定性的發(fā)展。

圖2t=2.0 s時刻不同速度梯度層厚度的界面形態(tài)和渦量分布

Fig.2 Interface pattern and vorticity distribution with different thickness of velocity gradient layers whent=2.0 s

不同邊界條件的渦量隨時間的變化如圖5所示，初始時刻的渦量均勻分布在兩種流體的界面上，在方向相反的初始速度的推動下，上層流體主動向下運動，下層流體主動向上運動，渦量逐漸卷入中心，兩種流體互相滲透形成漩渦，隨著渦量在中心的堆積，界面開始變得陡峭，擾動的振幅也逐漸增長，界面形態(tài)的轉(zhuǎn)變與渦量變化一致。Neumann邊界條件下，渦量向兩邊擴散得比較少，在方向相反的初始速度作用下迅速向中心聚集形成漩渦，當(dāng)漩渦達(dá)到一定強度時，渦量才開始往兩邊擴散；而在Dirichlet邊界條件下，渦量在K-H不穩(wěn)定性發(fā)展的初期已經(jīng)開始往兩邊擴散，因此界面中心形成的漩渦強度更弱，使得K-H不穩(wěn)定性向內(nèi)的發(fā)展受到很大影響。

圖3 速度差分別為0.4u,0.8u,1.2u和1.6u時的界面形態(tài)和渦量分布

Fig.3 Interface pattern and vorticity distribution for velocity difference is 0.4u,0.8u,1.2uand 1.6urespectively

圖4 界面在不同邊界條件下隨時間的形態(tài)演變

Fig.4 Pattern evolution of interface at different boundary conditions with time

圖5 界面在不同邊界條件下的渦量隨時間變化情況

Fig.5 Vorticity evolution of interface at different boundary conditions with time

圖6t=2.0 s時不同密度條件下表面張力系數(shù)分別為0.000，0.005，0.010和0.050時界面形態(tài)的演化

Fig.6 Pattern evolution of interface att=2.0 s for surface tension is 0.000,0.005,0.010 and 0.050 respectively with different density ratio

3.4 理查德森數(shù)的影響

4 結(jié) 論

本文應(yīng)用界面跟蹤法對兩相流的K-H不穩(wěn)定性進行數(shù)值模擬，模擬得到不同的速度梯度層的厚度、不同的速度差、不同的邊界條件和不同理查德森數(shù)對兩相流的K-H不穩(wěn)定性的影響，進而得出以下結(jié)論。

(1) 速度梯度層的厚度與界面在水平分量中的移動速度呈正相關(guān)，此外卷起的程度與界面在水平分量中的移動速度亦成正相關(guān)。

(2) 初始水平速度差與界面處的卷起程度有關(guān)，初始水平速度差越大，界面卷起越多，同時其擾動增長速度也越快，K-H不穩(wěn)定性的特征形式變得更加明顯。

(3) 在Neumann邊界條件和Dirichlet邊界條件下，界面的擾動發(fā)展情況不同，且在Neumann邊界條件下的界面擾動發(fā)展更快。

(4) 理查德森數(shù)對K-H不穩(wěn)定性發(fā)生的時間和界面的累積具有顯著影響。當(dāng)理查德森數(shù)較小時，K-H不穩(wěn)定的現(xiàn)象非常明顯。隨著理查德森數(shù)的增加，界面的卷起受到抑制。

(5) 相比于相場法[1]，本文使用的界面追蹤法具有較高的計算精度，數(shù)值模擬的結(jié)果更加精確。