瑞利-泰勒不穩(wěn)定問題的光滑粒子法模擬研究?

楊秀峰 劉謀斌

1)(DepartMent of Mechanical Engineering,Iowa State University,AMes,IA 50011,USA)

2)(北京大學(xué)工學(xué)院,北京 100187)

瑞利-泰勒不穩(wěn)定問題的光滑粒子法模擬研究?

楊秀峰1)劉謀斌2)?

1)(DepartMent of Mechanical Engineering,Iowa State University,AMes,IA 50011,USA)

2)(北京大學(xué)工學(xué)院,北京 100187)

(2017年4月1日收到;2017年6月2日收到修改稿)

提出了一種適用于模擬多相流的光滑粒子法,該方法對密度方程在交界面處的離散格式進行了修正以適應(yīng)多相流所涉及的大密度比問題,在不同相粒子之間施加了很小的排斥力以防止粒子穿透交界面,并采用了最新發(fā)展的雙曲型光滑函數(shù)以消除應(yīng)力不穩(wěn)定問題.應(yīng)用該多相流光滑粒子法模擬研究了單模態(tài)和多模態(tài)瑞利-泰勒不穩(wěn)定問題.通過與文獻中結(jié)果的對比研究表明:在模擬瑞利-泰勒不穩(wěn)定問題時,本文方法的結(jié)果明顯優(yōu)于文獻中的大部分光滑粒子法模擬結(jié)果,與Grenier等(2009 J.Comput.Phys.228 8380)的結(jié)果相當(dāng),但本文方法比Grenier等的方法簡單方便.對于單模態(tài)瑞利-泰勒不穩(wěn)定問題,研究了交界面的形態(tài),渦結(jié)構(gòu)的演化過程以及貫穿深度隨時間的變化關(guān)系.對于多模態(tài)瑞利-泰勒不穩(wěn)定問題,研究了交界面演化過程中小尺度結(jié)構(gòu)合并成大尺度結(jié)構(gòu)的過程,水平方向的平均密度隨高度的變化關(guān)系,以及貫穿深度隨時間的變化關(guān)系.

瑞利-泰勒不穩(wěn)定性,光滑粒子法,多相流

1 引 言

瑞利-泰勒不穩(wěn)定(Rayleigh-Taylor instability)問題是指當(dāng)密度大的流體處于密度小的流體之上時,在重力的作用下產(chǎn)生的一種交界面不穩(wěn)定現(xiàn)象.Rayleigh[1]最早從理論上分析了重力作用下密度不同的兩種流體的穩(wěn)定性,從理論上證明:當(dāng)上層流體的密度小于下層流體的密度時,交界面是穩(wěn)定的,擾動使交界面產(chǎn)生的振蕩會逐漸衰減;相反,如果上層流體的密度大于下層流體的密度,則交界面是不穩(wěn)定的,擾動導(dǎo)致的交界面振蕩會被迅速放大.Taylor[2]從理論上分析了重力作用下上層流體密度大于下層流體密度時交界面的穩(wěn)定性問題,指出交界面的初始擾動會隨時間以指數(shù)關(guān)系增加,并給出了相應(yīng)的變化公式.為了驗證Taylor的理論結(jié)果,LeWis[3]做了一系列的實驗,在圓柱形的裝置中觀測了不穩(wěn)定交界面的演化過程.實驗結(jié)果表明Taylor的結(jié)果適用于交界面演化的初始階段.LeWis[3]的實驗結(jié)果還表明交界面不穩(wěn)定性的發(fā)展過程可分為三個階段:首先是初始階段,擾動幅度以指數(shù)增加;隨后為過渡階段,擾動幅度的增加速度逐漸增大到最大值;最后是勻速階段,兩種流體相互貫穿的速度不再增大,貫穿深度勻速增加.

瑞利-泰勒不穩(wěn)定問題主要取決于兩種流體的密度和重力加速度.重力加速度越大,不穩(wěn)定性就越強.如果沒有重力,則不會產(chǎn)生瑞利-泰勒不穩(wěn)定現(xiàn)象.流體密度的影響一般用無量綱數(shù)——Atwood數(shù)(At)來表征:

其中ρH和ρL分別表示重的流體和輕的流體的密度.Atwood數(shù)的取值范圍是0≤At≤1.當(dāng)兩種流體的密度相等時(At=0),不存在瑞利-泰勒不穩(wěn)定問題.當(dāng)兩種流體密度差增大時,Atwood數(shù)也增大,交界面的不穩(wěn)定性增強.當(dāng)輕的流體不存在時(At=1),即上層流體置于真空之上,瑞利-泰勒不穩(wěn)定問題就成了流體的自由下落問題.對于二維問題,當(dāng)Atwood數(shù)較小時,相互貫穿的流體會形成兩個反向旋轉(zhuǎn)的渦.當(dāng)Atwood數(shù)增大時,則渦的尺寸會減小.當(dāng)Atwood數(shù)大到一定值時,上層密度大的流體穿過下層密度小的流體迅速下落,而密度小的流體則會迅速上浮[4].另外,流體的表面張力和黏性等因素也會影響瑞利-泰勒不穩(wěn)定性,因為它們會對流體的相互貫穿起到一定程度的阻礙作用[5].

瑞利-泰勒不穩(wěn)定問題不僅具有科學(xué)上的研究意義,而且具有應(yīng)用上的研究價值[6?8].因此,在Rayleigh和Taylor之后,這一問題成了流體交界面穩(wěn)定性研究的一個重要的基本問題.近幾十年,研究者們發(fā)展了許多模型來預(yù)測瑞利-泰勒不穩(wěn)定問題中兩種流體相互貫穿的深度[9?13].然而,到目前為止,還沒有一個模型能夠與實驗結(jié)果完全相符[10,14].對于交界面形狀的演化過程,更是沒有模型能夠預(yù)測.隨著計算流體力學(xué)的快速發(fā)展,數(shù)值模擬成了研究瑞利-泰勒不穩(wěn)定問題的一個重要手段[15?19].

本文的目的是給出一種適用于模擬多相流的光滑粒子法(smoothed particle hydrodynaMics,SPH method),并應(yīng)用該方法模擬研究單模態(tài)和多模態(tài)瑞利-泰勒不穩(wěn)定流動的演化過程.

2 數(shù)值方法

光滑粒子法是一種無網(wǎng)格粒子類方法[20?24].在模擬流體問題時,光滑粒子法將連續(xù)的流體離散成可運動的粒子,粒子之間通過光滑函數(shù)進行關(guān)聯(lián),然后將流體的控制方程離散成相應(yīng)的粒子方程進行數(shù)值求解.

在光滑粒子法中,任一函數(shù)f在粒子a上的值可以離散成如下求和形式:

其中腳標(biāo)a和b表示粒子,r表示位置矢量,W表示光滑函數(shù)(也稱為核函數(shù)),V表示粒子的體積.參數(shù)h表示光滑長度,用于控制光滑函數(shù)W的寬度,把上式中的求和范圍限制在有限的區(qū)域內(nèi).本文使用的光滑函數(shù)為雙曲型核函數(shù)[25,26]:

其中s=r/h,r表示兩個粒子之間的距離.該光滑函數(shù)能夠有效地消除光滑粒子法在模擬流體問題時潛在的應(yīng)力不穩(wěn)定問題[25,26].傳統(tǒng)的光滑粒子法一般采用鐘形光滑函數(shù)(如三次樣條函數(shù)和高斯函數(shù)等).在粒子附近很小的一個區(qū)域內(nèi),當(dāng)壓強不變時,鐘形光滑函數(shù)會導(dǎo)致粒子之間的排斥力隨著粒子間距的減小而減小,從而導(dǎo)致粒子在這個區(qū)域內(nèi)聚集,即應(yīng)力不穩(wěn)定現(xiàn)象.本文使用的雙曲型光滑函數(shù)不存在這種問題,因為當(dāng)壓強不變時,排斥力會隨著粒子間距的減小而增大,從而防止粒子聚集,消除應(yīng)力不穩(wěn)定問題.

本文使用納維-斯托克斯方程作為流體的控制方程,其拉格朗日形式為:

其中ρ,u,p,μ分別表示流體的密度、速度、壓強和黏性系數(shù);g表示重力加速度.

在光滑粒子法中,密度方程(4)可以離散成如下形式:

其中?aWab表示以粒子a為中心的光滑函數(shù)在粒子b處的梯度.對于密度差較大的多相流,方程(6)會高估密度大的流體粒子對密度小的流體粒子的貢獻.因此,在計算交界面處的大密度粒子對小密度粒子的作用時,采用如下形式的密度方程:

其中腳標(biāo)H和L分別表示密度大的流體和密度小的流體.與方程(6)相比,上式可以有效地處理大密度比多相流的交界面,降低交界面處的計算誤差,并有助于得到清晰光滑的交界面.需要注意的是,對于同相粒子之間的相互作用,以及小密度粒子對大密度粒子的作用,仍使用方程(6)計算.

動量方程(5)中的壓力梯度項可離散成(9)式得到的黏性一般會比真實的物理黏性稍小,可以通過人工黏性進行補償.本文采用如下形式的人工黏性[27]:

其中參數(shù)α用于控制人工黏性的大小,本文取值為α=0.02.人工黏性(11)僅在相鄰粒子相互靠近時起作用.另外,人工黏性還能在一定程度上降低粒子的非物理振動[27].

為了防止交界面兩側(cè)的粒子相互穿透,并得到光滑清晰的交界面,當(dāng)不同相的兩個粒子相互作用時,在它們之間施加如下形式的排斥力[28?30]:

其中參數(shù)εR的取值范圍為0—0.1,本文取值為0.02.該排斥力的形式實際上是由壓強梯度項的離散形式(8)式演變而來,不同的是,該斥力的值遠小于壓強梯度,而且只在不同相的粒子之間起作用,即在交界面兩側(cè)的粒子之間施加排斥力,從而抑制粒子穿透交界面.另外,該排斥力在一定程度上也能起到表面張力的作用.

控制方程(4)和(5)不封閉,可通過如下狀態(tài)方程來計算壓強:

其中c表示數(shù)值聲速,ρr表示參考密度.

在光滑粒子法中,壓強通常會產(chǎn)生較大的非物理波動.為了降低壓強波動,一般可以通過密度重置來實現(xiàn).本文采用如下形式的密度重置方程[31],每20個時間步修正一次密度:

對于粒子的位置移動,采用XSPH運動方式[27]:

其中參數(shù)ε的取值范圍為0—1,本文取值為0.5.(16)式的求和范圍是同相的相鄰粒子.另外,在求解密度方程時,也改用速度u′[32],其他方程中的速度仍采用原速度u.

3 單模態(tài)瑞利-泰勒不穩(wěn)定問題

圖1是單模態(tài)瑞利-泰勒不穩(wěn)定問題的示意圖,其中計算區(qū)域的高度為H=2,寬度為L=1,交界面初始形狀為y=1?0.15 sin(2πx),上層流體密度為ρH=1.8,下層流體密度為ρL=1.0,重力加速度為g=1.兩種流體的動力黏性系數(shù)相等:νH= νL=0.0025.Atwood數(shù)為At=2/7.所有邊界都是剛性無滑移固體壁面條件.初始時刻,整個系統(tǒng)處于靜止?fàn)顟B(tài).然后,在重力的作用下,上層重的流體開始向下運動,而下層輕的流體則向上運動.

圖1 單模態(tài)瑞利-泰勒不穩(wěn)定問題的初始設(shè)置Fig.1.Initial set-up for single-Mode Ray leigh-Tay lor instability p rob lem.

為了驗證本文提出的多相流光滑粒子法方法,圖2對比了本文的模擬結(jié)果和文獻中的幾種結(jié)果.圖2中本文模擬結(jié)果采用的空間分辨率為200×400,其他結(jié)果采用的空間分辨率分別是Levelset 312×624[29];Grenier等[29]300×600;Chen等[33]200×400,Monaghan和Rafiee[34]150×300;Hu和AdaMs[35]60×120.圖2表明各種模擬結(jié)果中交界面的整體結(jié)構(gòu)大致相似,但局部結(jié)構(gòu)差別較大,例如Chen等[33]與Hu和Adams[35]的交界面不清晰,交界面附近不同種類的流體粒子相互穿透的現(xiàn)象非常明顯,而且局部的復(fù)雜變化沒能模擬出來.Monaghan和Rafiee[34]的交界面比較清晰,但局部的復(fù)雜變化也沒能模擬出來.相對而言,本文的結(jié)果與Level-set和Grenier等[29]的結(jié)果不但交界面清晰,而且局部細節(jié)豐富.

圖2 (網(wǎng)刊彩色)本文結(jié)果與文獻中結(jié)果的對比(t=5)Fig.2.(color online)CoMparison of the resu lts of this paper and literature(t=5).

圖3 (網(wǎng)刊彩色)三種空間分辨率的模擬結(jié)果對比:200×400(左),300×600(中),400×800(右)Fig.3.(color on line)CoMparison of siMu lating resu lts using three diff erent resolutions:200×400(left),300×600(Midd le),400×800(right).

需要指出的是,圖2中各種模擬結(jié)果采用的空間分辨率雖然并不完全相同,但差別不大,只有Hu和AdaMs[35]的分辨率明顯偏低.圖3對比了本文方法采用三種分辨率模擬的結(jié)果.從圖3可以看出,隨著分辨率的增加,交界面變得更加光滑,而且細節(jié)增多,但是這種變化并不顯著.因此,分辨率會影響模擬結(jié)果,但并不是導(dǎo)致圖2中不同模擬結(jié)果的主要原因.圖3中的顏色表示流體壓強,可以看出,本文方法得到的壓強分布不但整體上光滑,而且交界面和壁面附近的壓強也是連續(xù)的,沒有出現(xiàn)異常現(xiàn)象.本文模擬結(jié)果的交界面提取方法見文獻[36],該方法先以粒子為頂點構(gòu)建三角形網(wǎng)格,然后在網(wǎng)格上提取交界面.

圖4給出了時間t=1,3,5和7時的交界面形態(tài).圖5給出了與圖4相對應(yīng)的幾個時刻的整個流場中渦的結(jié)構(gòu).從圖4可以看出,交界面從初始時刻的簡單形狀逐漸演化成非常復(fù)雜的形狀.從圖5可以看出,流場中的渦結(jié)構(gòu)也隨著交界面形態(tài)的變化而變得逐漸復(fù)雜,且渦的數(shù)量也隨之增加.流體的運動產(chǎn)生渦,同時,渦也反過來影響流體的運動.初始階段,上層流體從左側(cè)向下運動,下層流體則從右側(cè)向上運動,形成一個逆時針的大渦,并在左右壁面附近形成兩個順時針的小渦(見圖5,t=1).這三個渦的中心幾乎在同一條直線上.隨著兩種流體的相互貫穿,在三個渦的作用下,交界面形成了方向相反的兩個“蘑菇”狀結(jié)構(gòu)(見圖4和圖5,t=3).三個渦都位于“蘑菇”結(jié)構(gòu)的邊沿處,并逐漸增大.同時,在右上角新形成一個順時針的小渦.由于這個小渦與右中部的渦相互靠近,并且方向相同,最終被吞并,形成一個大渦.在右側(cè)渦的擠壓作用下,中間的大渦在左上方分裂出一個小渦.與此同時,左側(cè)原有的渦向下運動,并逐漸增大,占據(jù)了下方的位置,而且也分裂出一個小渦(見圖5,t=5).值得注意的是,盡管原先的三個渦都改變了位置,且經(jīng)歷了合并或分裂過程,但它們的中心位置依然保持在幾乎同一條直線上.最后,到t=7時,因為兩種流體在底部的剪切作用,導(dǎo)致了Kelvin-Helmholtz不穩(wěn)定現(xiàn)象,在底部形成了三個逆時針旋轉(zhuǎn)的小渦.

圖4 (網(wǎng)刊彩色)單模態(tài)瑞利-泰勒不穩(wěn)定交界面的演化過程(t=1,3,5,7)Fig.4.(color on line)Evolu tion of interface for single-Mode Ray leigh-Tay lor instability p rob lem(t=1,3,5,7).

圖5 (網(wǎng)刊彩色)與圖4對應(yīng)時刻的渦結(jié)構(gòu)的演化過程(t=1,3,5,7)Fig.5.(color on line)Evolution of vortices at the saMe tiMe of Fig.4(t=1,3,5,7).

圖6 單模態(tài)瑞利-泰勒不穩(wěn)定交界面貫穿深度隨時間的變化(HH表示上層重的流體向下進入輕的流體的貫穿深度;HL表示下層輕的流體向上進入重的流體的貫穿深度)Fig.6.Penetration distance as a function of tiMe for single-Mode Ray leigh-Tay lor instability interface(HHis the distance of the heavy fl uid penetrating into the light fl uid;HLis the distance of the light fl uid penetrating into the heavy fl uid).

圖6給出了兩種流體的貫穿深度隨時間的變化過程.圖6中縱坐標(biāo)0表示高度方向的中點,而1則表示上下兩端的壁面.從圖6可以看出,初始階段t＜1時,兩種流體離開中點的距離幾乎同步增加.在t=1之后,重的流體向下運動的速度明顯快于輕的流體向上運動的速度.在t=4之后,重的流體已經(jīng)比較接近底部壁面,因此運動速度迅速降低;而輕的流體因為運動速度較慢,一直到t=6之后才開始降低速度.在t=7時,兩種流體已非常接近兩端壁面.根據(jù)LeWis[3]的實驗和Layzer[37]的理論結(jié)果,在不考慮壁面影響的情況下,瑞利-泰勒不穩(wěn)定導(dǎo)致的貫穿深度在初始階段隨時間成指數(shù)增加,穩(wěn)定后隨時間成勻速增加關(guān)系.在圖6中,貫穿深度小于0.3時,約為指數(shù)增加階段,貫穿深度為0.3至0.8時,約為勻速增加階段.因此,本文的數(shù)值結(jié)果與LeWis[3]的實驗結(jié)論及Layzer[37]的理論結(jié)果相符.對于貫穿深度大于0.8時,因兩端壁面的影響,變?yōu)闇p速階段.

4 多模態(tài)瑞利-泰勒不穩(wěn)定問題

多模態(tài)瑞利-泰勒不穩(wěn)定問題的初始設(shè)置與單模態(tài)的設(shè)置圖1相似,高度為H=2,寬度為L=1.不同的是,多模態(tài)問題的交界面形狀為y=1?0.05sin(πx/5),上層流體密度為ρH=3,下層流體密度為ρL=1,重力加速度為g=1.兩種流體的動力黏性系數(shù)相等:νH=νL=0.001.Atwood數(shù)為At=0.5.上下兩端為剛性無滑移固體壁面條件,左右邊界為周期邊界條件.

圖7 (網(wǎng)刊彩色)多模態(tài)瑞利-泰勒不穩(wěn)定交界面的演化過程Fig.7.(color on line)Evolution of interface forMulti-Mode Ray leigh-Tay lor instability p rob lem.

圖8 平均密度隨高度的變化(t=3,5,7)Fig.8.Average density as a function of height(t=3,5,7).

圖7給出了多模態(tài)瑞利-泰勒不穩(wěn)定交界面的演化過程.初始階段,初始擾動的波形被逐漸拉長,相鄰結(jié)構(gòu)的影響很小 (見圖7,t=1).隨著貫穿深度的增加,相鄰結(jié)構(gòu)相互影響,形成蘑菇狀結(jié)構(gòu)(見圖7,t=2和3).隨著蘑菇狀結(jié)構(gòu)的增大,相鄰結(jié)構(gòu)的相互作用也加劇,交界面整體結(jié)構(gòu)的對稱性被打破,有的蘑菇狀結(jié)構(gòu)被拉長,有的蘑菇狀結(jié)構(gòu)被壓扁(見圖7,t=4).然后小尺度結(jié)構(gòu)相互合并,形成不同尺度的大結(jié)構(gòu) (圖7,t=5—9).交界面的形態(tài)也隨之變得異常復(fù)雜.

圖9 多模態(tài)瑞利-泰勒不穩(wěn)定交界面的貫穿深度隨時間的變化(HH表示上層重的流體向下進入輕的流體的貫穿深度;HL表示下層輕的流體向上進入重的流體的貫穿深度)Fig.9.Penetration d istance as a function of tiMe for Mu lti-Mode Ray leigh-Tay lor instability interface(HHis the d istance of the heavy fl uid penetrating into the light fl uid;HLis the distance of the light fl uid penetrating into the heavy fl uid).

圖8是水平方向的平均密度隨高度的變化情況.在t=3時,流體的貫穿深度相對較小,交界面形狀也比較規(guī)則,所以平均密度的變化比較規(guī)則.到t=5時,密度出現(xiàn)了不同尺度的變化,這是因為交界面的小結(jié)構(gòu)已逐漸合并成不同尺度的大結(jié)構(gòu),而這種結(jié)構(gòu)尺度的變化也體現(xiàn)在了平均密度的變化上.無規(guī)則的尺度變化逐漸加劇,到t=7時,已經(jīng)完全看不出初始階段的規(guī)則變化.

圖9給出了兩種流體相互貫穿的深度隨時間的變化情況.從圖9可以看出,兩條線的形狀相似.初始階段,貫穿深度的增加稍快,而后緩慢增加.在t=4附近,穿透距離加速增加,而后幾乎線性增加.這是因為,原先的小“泡”開始合并成大“泡”,而大“泡”的運動速度大于小“泡”.當(dāng)貫穿深度到達0.9附近時,因為上下兩端壁面的作用,貫穿深度趨于穩(wěn)定,貫穿速度迅速降低.

5 結(jié) 論

本文在傳統(tǒng)光滑粒子法的基礎(chǔ)上發(fā)展了一種適用于模擬多相流的光滑粒子法,并應(yīng)用該方法模擬研究了單模態(tài)和多模態(tài)瑞利-泰勒不穩(wěn)定問題.該方法對密度方程在交界面處的離散格式進行了修正以適應(yīng)多相流所涉及的大密度比情況,在不同相粒子之間施加了很小的排斥力以防止粒子穿透交界面,并采用了最新發(fā)展的雙曲型光滑函數(shù)以消除應(yīng)力不穩(wěn)定問題.為了驗證本文發(fā)展的多相流光滑粒子法,將本文與文獻中的幾種模擬結(jié)果進行了對比.結(jié)果表明,在模擬瑞利-泰勒不穩(wěn)定問題時,本文方法的結(jié)果明顯優(yōu)于文獻中的大部分光滑粒子法模擬結(jié)果,與文獻[29]的結(jié)果相當(dāng),但本文方法比文獻[29]的方法更簡單方便.

對于單模態(tài)瑞利-泰勒不穩(wěn)定問題,給出了交界面的形態(tài)及渦結(jié)構(gòu)的演化過程,并研究了貫穿深度與時間的關(guān)系.貫穿深度隨時間的變化過程與實驗及理論結(jié)果相符.對于多模態(tài)瑞利-泰勒不穩(wěn)定問題,研究了交界面演化過程中小尺度結(jié)構(gòu)合并成大尺度結(jié)構(gòu)的過程.在此過程中,水平方向的平均密度隨高度的變化由初始時的規(guī)則分布逐漸演化成無規(guī)則分布.小尺度結(jié)構(gòu)向大尺度結(jié)構(gòu)的演化也導(dǎo)致了貫穿深度隨時間的增加速度由線性改變成非線性.合并過程結(jié)束后,貫穿深度再次與時間成線性關(guān)系,但是增加速度明顯大于合并前的增加速度.

[1]Ray leigh L 1883 Proc.Lond.Math.Soc.14 170

[2]Tay lor G 1950 Proc.Roy.Soc.A:Math.Phys.201 192

[3]LeWis D J 1950 Proc.Roy.Soc.A:Math.Phys.202 81

[4]Tryggvason G 1988 J.CoMput.Phys.75 253

[5]Banerjee R,K an jilal S 2015 J.Pure App l.Ind.Phys.5 73

[6]Sharp D H 1984 Physica D 12 3

[7]K ilkenny J,G lendinning S,Haan S,HamMel B,Lind l J,Mun ro D,ReMington B,Weber S,Knauer J,Verdon C 1994 Phys.P lasMas 1 1379

[8]Huang C S,Kelley M,Hysell D 1993 J.Geophys.Res.-Space Physics 98 15631

[9]A lon U,Hecht J,O fer D,Shvarts D 1995 Phys.Rev.Lett.74 534

[10]DiMonte G 2000 Phys.P lasMas 7 2255

[11]RaMshaWJ D 1998 Phys.Rev.E 58 5834

[12]G limMJ,Saltz D,Sharp D H 1998 Phys.Rev.Lett.80 712

[13]Cheng B,G limMJ,Sharp D 2002 Phys.Rev.E 66 036312

[14]Zhang Y S,He Z W,Gao F J,Li X L,T ian B L 2016 Phys.Rev.E 93 063102

[15]He X,Chen S,Zhang R 1999 J.CoMput.Phys.152 642

[16]Kadau K,Barber J L,GerMann T C,Holian B L,A lder B J 2010 Phil.Trans.R.Soc.A 368 1547

[17]RaMap rabhu P,K arkhanis V,Banerjee R,Varshochi H,K han M,LaWrie A 2016 Phys.Rev.E 93 013118

[18]Sagert I,Howell J,Staber A,Strother T,Colbry D,Bauer W2015 Phys.Rev.E 92 013009

[19]Liang H,Li Q,Shi B,Chai Z 2016 Phys.Rev.E 93 033113

[20]Lucy L B 1977 Astron.J.82 1013

[21]G ingold R A,Monaghan J J 1977 Mon.Not.R.Astron.Soc.181 375

[22]Yang X,Peng S,Liu M,Shao J 2012 In t.J.CoMp.Meth.-Sing 9 1240002

[23]Yang X F,Peng S L,Liu MB 2014 Appl.Math.Model 38 3822

[24]Yang X,Dai L,Kong SC 2017 Proc.CoMbust.Inst.36 2393

[25]Yang X F,Liu MB 2012 Acta Phys.Sin.61 224701(in Chinese)[楊秀峰,劉謀斌 2012物理學(xué)報 61 224701]

[26]Yang X F,Liu MB,Peng S 2014 CoMput.F luids 92 199

[27]Monaghan J J 1992 Ann.Rev.Astron.Astrophys.30 543

[28]Monaghan J J 2000 J.CoMput.Phys.159 290

[29]G renier N,Antuono M,Colagrossi A,Le TouzéD,A lessand rini B 2009 J.CoMput.Phys.228 8380

[30]Yang X,Liu M2013 Sci.China:Phys.Mech.Astron.56 315

[31]Bonet J,Lok T S 1999 CoMput.Methods App l.Mech.Engrg.180 97

[32]Colagrossi A,Land rini M2003 J.CoMput.Phys.191 448

[33]Chen Z,Zong Z,Liu M,Zou L,Li H,Shu C 2015 J.CoMpu t.Phys.283 169

[34]Monaghan J,Rafiee A 2013 In t.J.NuMerical Mech.F luids 71 537

[35]Hu X Y,AdaMs N A 2009 J.CoMput.Phys.228 2082

[36]Yang X F,Liu MB 2016 Chin.J.CoMput.Mech.33 594(in Chinese)[楊秀峰,劉謀斌 2016計算力學(xué)學(xué)報 33 594]

[37]Layzer D 1955 Astrophys.J.122 1

PACS:47.20.–k,47.11.–j,68.05.–nDOI:10.7498/aps.66.164701

*Pro ject supported by the National Natural Science Foundation of China(G rant Nos.11302237,U 1530110).

?Corresponding author.E-Mail:mb liu@pku.edu.cn

N uMerical study of Ray leigh-Tay lor instab ility by using sMoothed particle hyd rodynaMics?

Yang Xiu-Feng1)Liu Mou-Bin2)?

1)(DepartMent ofMechanical Engineering,Iowa State University,AMes,IA 50011,USA)
2)(BIC-ESAT,College of Engineering,Peking University,Beijing 100187,China)

1 Ap ril 2017;revised Manuscrip t

2 June 2017)

In this paper,we present a sMoothed particle hydrodynaMics(SPH)Method for Modeling multiphase floWs.The multiphase SPH method includes a corrective discretization scheme for density app roximation around the fl uid interface to treat large density ratio,a sMall repulsive force between particles froMdiff erent phases to p revent particles froMunphysically penetrating fluid interface,and a neWly-developed hyperbolic-shaped kernel function to remove possible stress instability.This mu ltiphase SPH Method is then used to study the single-and multi-Mode Rayleigh-Taylor instability p robleMs.A coMparison between our results With the resu lts froMexisting literature shoWs that our results are obviously better than most available results froMother SPH simulations.The present results are close to those by G renier et al.while the present multiphase SPH Method is siMp ler and easier to iMp leMent than that in the work by G renier et al.(G renier,et al.2009 J.CoMput.Phys.228 8380).For the single-Mode Rayleigh-Tay lor instability,the evolutions of the interface pattern and vortex structures,and the penetration depth each as a function of time are investigated.For themu lti-Mode Rayleigh-Tay lor instability,theMerging of sMall structures into a large structure during the evolution of the interface is studied.The horizontal average density and the penetration each as a function of height are also studied.

Rayleigh-Taylor instability,smoothed particle hydrodynaMics,multiphase flow

10.7498/aps.66.164701

?國家自然科學(xué)基金(批準號:11302237,U 1530110)資助的課題.

?通信作者.E-Mail:Mb liu@pku.edu.cn

?2017中國物理學(xué)會C h inese P hysica l Society

http://Wu lixb.iphy.ac.cn