楊秀峰 劉謀斌

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

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

瑞利-泰勒不穩(wěn)定問(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日收到修改稿)

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

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

1 引 言

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

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

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

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

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

2 數(shù)值方法

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

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

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

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

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

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

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

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

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

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

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

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

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

控制方程(4)和(5)不封閉,可通過(guò)如下?tīng)顟B(tài)方程來(lái)計(jì)算壓強(qiáng):

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

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

對(duì)于粒子的位置移動(dòng),采用XSPH運(yùn)動(dòng)方式[27]:

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

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

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

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

為了驗(yàn)證本文提出的多相流光滑粒子法方法,圖2對(duì)比了本文的模擬結(jié)果和文獻(xiàn)中的幾種結(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]的交界面不清晰,交界面附近不同種類(lèi)的流體粒子相互穿透的現(xiàn)象非常明顯,而且局部的復(fù)雜變化沒(méi)能模擬出來(lái).Monaghan和Rafiee[34]的交界面比較清晰,但局部的復(fù)雜變化也沒(méi)能模擬出來(lái).相對(duì)而言,本文的結(jié)果與Level-set和Grenier等[29]的結(jié)果不但交界面清晰,而且局部細(xì)節(jié)豐富.

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

圖3 (網(wǎng)刊彩色)三種空間分辨率的模擬結(jié)果對(duì)比: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對(duì)比了本文方法采用三種分辨率模擬的結(jié)果.從圖3可以看出,隨著分辨率的增加,交界面變得更加光滑,而且細(xì)節(jié)增多,但是這種變化并不顯著.因此,分辨率會(huì)影響模擬結(jié)果,但并不是導(dǎo)致圖2中不同模擬結(jié)果的主要原因.圖3中的顏色表示流體壓強(qiáng),可以看出,本文方法得到的壓強(qiáng)分布不但整體上光滑,而且交界面和壁面附近的壓強(qiáng)也是連續(xù)的,沒(méi)有出現(xiàn)異?，F(xiàn)象.本文模擬結(jié)果的交界面提取方法見(jiàn)文獻(xiàn)[36],該方法先以粒子為頂點(diǎn)構(gòu)建三角形網(wǎng)格,然后在網(wǎng)格上提取交界面.

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

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

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

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

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

圖9 多模態(tài)瑞利-泰勒不穩(wěn)定交界面的貫穿深度隨時(shí)間的變化(HH表示上層重的流體向下進(jìn)入輕的流體的貫穿深度;HL表示下層輕的流體向上進(jìn)入重的流體的貫穿深度)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時(shí),流體的貫穿深度相對(duì)較小,交界面形狀也比較規(guī)則,所以平均密度的變化比較規(guī)則.到t=5時(shí),密度出現(xiàn)了不同尺度的變化,這是因?yàn)榻唤缑娴男〗Y(jié)構(gòu)已逐漸合并成不同尺度的大結(jié)構(gòu),而這種結(jié)構(gòu)尺度的變化也體現(xiàn)在了平均密度的變化上.無(wú)規(guī)則的尺度變化逐漸加劇,到t=7時(shí),已經(jīng)完全看不出初始階段的規(guī)則變化.

圖9給出了兩種流體相互貫穿的深度隨時(shí)間的變化情況.從圖9可以看出,兩條線的形狀相似.初始階段,貫穿深度的增加稍快,而后緩慢增加.在t=4附近,穿透距離加速增加,而后幾乎線性增加.這是因?yàn)?原先的小“泡”開(kāi)始合并成大“泡”,而大“泡”的運(yùn)動(dòng)速度大于小“泡”.當(dāng)貫穿深度到達(dá)0.9附近時(shí),因?yàn)樯舷聝啥吮诿娴淖饔?貫穿深度趨于穩(wěn)定,貫穿速度迅速降低.

5 結(jié) 論

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

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

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):11302237,U 1530110)資助的課題.

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

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

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