康 寧，黎一鍇，何 旭

Rayleigh-Taylor不穩(wěn)定性非線性特性的數(shù)值研究

康 寧，黎一鍇*，何 旭

（北京理工大學 機械與車輛學院，北京100081）

以往對于單模態(tài)Rayleigh-Taylor（RT）不穩(wěn)定性非線性特性的研究主要集中于推導和測量恒定的氣泡推迸速度上，而缺乏對液態(tài)尖釘區(qū)域非線性動力學特性的詳細分析。采用耦合的Level-Set和Volume-of-Fluid（CLSVOF）界面捕捉方法對單模態(tài)RT不穩(wěn)定性的發(fā)展過程迸行了精確的數(shù)值模擬，并利用模擬得到的壓力場和速度場信息對RT不穩(wěn)定性非線性發(fā)展階段的穩(wěn)態(tài)動力學特性迸行了分析。模擬結果表明，在液態(tài)尖釘根部由于慣性力作用而引起的水平?jīng)_擊流會在此處形成一個局部最大壓力點，由于此處慣性力與壓強梯度的平衡，位于最大壓力點附近的流動最終將達到穩(wěn)態(tài)。通過理論分析，確定了此穩(wěn)態(tài)流動中各穩(wěn)態(tài)特征參數(shù)與初始擾動波長、慣性加速度之間的關系。這些特征參數(shù)的確定有助于將經(jīng)典低速射流的相關理論擴展應用到RT不穩(wěn)定性誘導霧化的研究領域。

Rayleigh-Taylor不穩(wěn)定性；非線性；穩(wěn)態(tài)；CLSVOF方法；數(shù)值模擬

當2種不同密度的流體層組成的系統(tǒng)受到由重流體指向輕流體的加速度時，流體層間的界面是不穩(wěn)定的，將產(chǎn)生所謂的Rayleigh-Taylor（RT）不穩(wěn)定性。RT不穩(wěn)定性在天體物理［1-2］、激光核聚變［3-4］以及液體霧化［5-6］等科學和工業(yè)領域里起著非常重要的作用。在柴油機缸內(nèi)噴霧燃燒過程中，液滴表面的RT不穩(wěn)定性被認為是實現(xiàn)液滴二次霧化的主要原因之一［7-9］。因此，對RT不穩(wěn)定性的研究受到了國內(nèi)外學者的廣泛關注。

在單模態(tài)RT不穩(wěn)定性的發(fā)展初期，由于擾動振幅較擾動波長λ是一個小量，其動力學特性可以用數(shù)學上處理比較簡單的線性理論來描述。在線性理論中，擾動振幅將隨時間以指數(shù)的形式增長。在忽略表面張力和流體黏性的條件下，線性增長率可表示為為波數(shù)，A為加速度，At＝（ρ2-ρ1）／（ρ2+ρ1）為Atwood數(shù)，ρ1和ρ2分別為輕流體和重流體的密度。表面張力和流體黏性的存在都阻礙著RT不穩(wěn)定性的發(fā)展。Bellman和Pennington［11］考慮了起穩(wěn)定作用的表面張力，推導出線性增長率為為表面張力系數(shù)，可以看出存在一個臨界波數(shù)只有波數(shù)小于臨界波數(shù)的擾動才能隨時間不斷增長。Piriz等［12］考慮了起阻尼作用的流體黏性，推導出線性增長率為某二次方程的解?？紤]流體其他性質，如可壓縮性、密度梯度等對線性增長率影響規(guī)律的研究可參考文獻［13］。

當初始擾動的振幅由于RT不穩(wěn)定性指數(shù)增長到可以與波長相比較時，由于此時的非線性影響不再可以忽略，擾動的發(fā)展將偏離線性理論。實驗中可以觀察到在RT不穩(wěn)定性的非線性發(fā)展階段，界面將呈現(xiàn)不對稱的氣泡（bubble）和尖釘（spike）形狀，而且氣泡以恒定的速度向高密度的液層推進發(fā)展。Goncharov［14］基于漸進勢流模型從理論上推出了二維和三維情況下此恒定速度的表達式。

由于非線性在數(shù)學處理上的復雜性，雖然可以推導出局部個別參數(shù)（如氣泡的恒定推進速度）的漸進解，但是很難得到描述全局范圍動力學特性的解析解（如流體內(nèi)的壓力場和速度場分布）。近年來，隨著計算機計算速度以及計算方法的不斷發(fā)展，利用數(shù)值模擬研究 RT不穩(wěn)定性的研究也越來越多。Baker等［15］利用渦層方法（vortex sheetmethod）模擬了無黏、不可壓流體的二維RT不穩(wěn)定性發(fā)展過程，并研究了液層厚度的影響規(guī)律。Ramaprabhu等［16-17］通過模擬三維RT不穩(wěn)定性研究了不同Atwood數(shù)下（0.005～1）界面的發(fā)展特性。在國內(nèi)，葉文華等［18-20］利用高精度的通量校正傳輸（Flux Corrected Transport，F(xiàn)CT）方法對激光燒蝕RT不穩(wěn)定性進行了模擬，研究了電子熱傳導燒蝕在長波長擾動的非線性RT不穩(wěn)定性演變中所起的作用。程會方等［21］利用移動粒子半隱式法（Moving-Particle Sem i-implicitmethod，MPS）模擬了二維RT不穩(wěn)定性，得到的界面形狀與實驗結果在表觀上基本一致。

可以看出，以往關于RT不穩(wěn)定性非線性動力學特性的研究主要集中于氣泡，尤其是氣泡的恒定推進速度上。但是，在利用RT不穩(wěn)定性實現(xiàn)液體霧化的應用方面，尖釘?shù)膭恿W特性更加重要。恒定的氣泡推進速度意味著恒定的液體波谷表面下沉速度，進而使液體以恒定的速度流入尖釘區(qū)域，使尖釘頂部以恒定的速度破碎形成小液滴。因此，深入研究這些恒定的非線性動力學特性將有利于更好地利用RT不穩(wěn)定性實現(xiàn)霧化。本文基于數(shù)值計算得到的壓力場與速度場結果，對RT不穩(wěn)定性的非線性動力學特性，尤其是對液體霧化起重要作用的恒定特性進行了深入分析。

1 物理模型和計算方法

1．1 控制方程

本文考慮氣-液兩相流的RT不穩(wěn)定性。密度為ρ1的氣體層覆蓋在密度為ρ2（ρ2＞ρ1）的液體層上方，且整個系統(tǒng)以恒定的加速度A+g向下加速運動，g為重力加速度。本文中氣體和液體均被視為理想不可壓的牛頓流體。液體霧化過程中液-氣密度之比較大，一般在100～1 000范圍內(nèi)變化，因而典型的Atwood數(shù)At非常接近于1。為保證數(shù)值計算的穩(wěn)定性，本文選取的液-氣密度比ρ2／ρ1為100，相應的Atwood數(shù)At為0.98。將參考坐標系固定在底面固壁上，則作用于流體系統(tǒng)的有效債性加速度為A，如圖1所示。圖中：ξ為表面變形的位移。在此參考坐標系下，整個系統(tǒng)的控制方程為

式中：u為速度矢量；ρ為流體密度函數(shù)；p為流體壓強；ey為y方向上的單位向量。

圖1 物理模型示意圖Fig.1 Schematic of physicalmodel

表面張力通過連續(xù)表面力（Continuum Surface Force，CSF）模型［22］轉化為體積力Fs進入控制方程。

式中：A*表征了起不穩(wěn)定作用的債性力和起穩(wěn)定作用的表面張力之間的比值。

1．2 定解條件

本文計算在二維矩形區(qū)域內(nèi)進行，計算區(qū)域的尺寸如圖1所示。在水平方向（x方向），由于只考慮單模態(tài)擾動的RT不穩(wěn)定性，選取x方向的跨度為Lx＝λ＝2π／k。根據(jù)文獻［15］，形成尖釘?shù)牧黧w動力學行為均發(fā)生在距離表面1／k的范圍之內(nèi)，因此，本文選取液層的初始厚度h0＝λ＞1／k。雖然隨著時間的發(fā)展，液層的波谷表面會下降，但在本文所有的模擬工況中，在計算結束時液層表面的最低點到底面固壁的距離仍有0．5λ＞1／k。所以，液層厚度在本文中不會對尖釘?shù)膭恿W特性產(chǎn)生較大的影響。為給尖釘提供足夠的縱向發(fā)展空間，氣層豎直方向的跨度設為1.5λ，整個計算區(qū)域豎直方向的總跨度為Ly＝2．5λ。

在計算開始時，給予氣-液界面一個非常小的正弦豎直速度初始擾動vs（t＝0，x）＝vs0sin（kx）。計算區(qū)域的左邊界和右邊界均設為周期邊界，底部邊界設為滑移壁面，頂部邊界設為自由邊界。具體的初始和邊界條件已標示在圖1中。

本文采用均勻交錯網(wǎng)格對計算區(qū)域進行離散化。為確定網(wǎng)格的尺寸，筆者研究了典型模擬工況下（A*＝8.8π3）不同網(wǎng)格（Δ＝Δx＝Δy＝λ／16，λ／32，λ／40，λ／48，λ／56，λ／64）對主要非線性動力學參數(shù)（恒定波谷沉降速度 V和波谷振幅 δ0）的影響，如圖2所示。結果顯示，當小于λ／48之后，網(wǎng)格尺寸對主要的非線性動力學特性影響非常小。本文數(shù)值模擬采用的網(wǎng)格尺寸為Δ＝Δx＝Δy＝λ／64。時間步長Δt受2個條件的限制，即自由表面和表面張力波在一個時間步長內(nèi)不能移動超過一個網(wǎng)格尺寸的距離。

式中：Cr為Courant數(shù)，本文中取為0.25；umax和vmax分別為每個計算時刻水平速度和豎直速度的最大值。

圖2 網(wǎng)格尺寸對波谷沉降速度V和振幅δ0的影響Fig.2 Effects of grid sizes on trough descending velocity V and amplitudeδ0

1．3 計算方法

氣-液RT不穩(wěn)定性模擬的關鍵是氣-液界面的捕捉方法。本文采用的計算方法建立在耦合Level-Set和Volume-of-Fluid（CLSVOF）界面捕捉方法基礎之上［23］。氣-液界面通過一條等LS函數(shù)線φ＝0隱式表征，而VOF方法作為LS方法的補充以保證液體質量守恒［24］。數(shù)值計算方法的細節(jié)可參考文獻［25］，作者利用同樣的數(shù)值計算方法成功研究了Faraday不穩(wěn)定性中具體的動力學特性。

1．4 數(shù)值驗證

為驗證本文所用數(shù)值方法以及程序代碼的正確性，將利用本文代碼的計算結果與Baker等［26］利用渦層方法進行模擬的結果進行比較。圖3給出了初始擾動振幅為0.1λ條件下，尖釘頂點和氣泡頂點位移隨時間的變化規(guī)律。圖中還給出了線性理論解作為參考?？梢钥闯觯瑹o論在線性發(fā)展階段還是非線性發(fā)展階段，利用本文程序代碼計算的結果與Baker等［26］的計算結果符合得很好。

圖3 模型驗證結果Fig.3 Results ofmodel validation

2 計算結果與討論

從RT不穩(wěn)定性線性增長率的表達式可以看到，對于大Atwood數(shù)，只有在無量綱參數(shù)A*＝ρ1λ2A／σ＞4π2的條件下，表面變形的振幅才能夠增長。本文對A*＝2.4π3，3.2π3，4.0π3，4.8π3，5.6π3，6.4π3，7.2π3，8.0π3和8.8π3等9個值進行了模擬。所有模擬的初始速度擾動振幅為計算在液體尖釘?shù)捻敳康竭_計算域的上邊界后停止。

2．1 非線性表面變形的整體動力學特性

圖4給出了RT不穩(wěn)定性非線性階段從液層表面形成尖釘?shù)陌l(fā)展過程。圖中所有不同時刻界面形狀曲線都是在一個工況下計算得到的?？梢钥闯?，在表面波峰的根部會形成一個局部最大壓強點pmax，而在這個最大壓強點以上則會形成一個動力學上獨立于底部液層的尖釘區(qū)域。這種動力學特性與Faraday不穩(wěn)定性非常類似［25］。由于本文關注的是Atwood數(shù)很大的情況，因此不會在尖釘頂部出現(xiàn)Kelvin Helmhotz（KH）不穩(wěn)定性所引起的卷曲（roll-up）結構［16，26］。

當把形成局部最大壓力點之后不同時刻表面波谷的最低點重疊在一起時，如圖4（b）所示，可以看到如下現(xiàn)象：①不同時刻波谷區(qū)域的液面保持著振幅為δ0的正弦曲線形狀（如圖4（b）中虛正弦曲線所示）；②在波峰區(qū)域可以看到不同時刻的最大壓力點（如圖4（b）中的小圓圈所示）重合在一起，而且最大壓力點到虛正弦波峰頂點的距離η也保持恒定的值。

圖5給出了模擬工況A*＝4.8π3下，波谷最低點的垂直坐標（常被稱為氣泡頂點）、表面變形中性點的垂直坐標和局部最大壓力點的垂直坐標隨時間的變化規(guī)律曲線。作為參考，圖5也給出了線性理論中、表面變形中性點與波谷最低點的距離以及最大壓力點與虛正弦波峰頂點的距離隨時間的變化曲線。其他模擬工況下的各曲線變化規(guī)律與圖5相似。

從圖5中可以看出，在發(fā)展的初始階段，表面變形嚴格地遵循線性理論解：

圖4 液層表面的動力學特征參數(shù)Fig.4 Features of jet formation from a liquid layer

圖5 和隨時間的變化曲線Fig.5 Temporal changing curves oand

隨著表面振幅的不斷增大，非線性因素變得逐漸明顯，從而降低了表面振幅的增長速度，最終，波谷以恒定的速度V*向下沉降。在線性理論中，中性表面的高度會保持在平衡位置，即然而由于非線性的影響在某個時刻之后開始以與波谷沉降相同的速度下降，因而使波谷正弦表面的振幅保持在一個固定值上。可以為界將整個發(fā)展過程分為2個階段：線性發(fā)展階段（t*＜t*1）和非線性發(fā)展階段從圖5中還可以看出和的下降速度一致，這表明最大壓力點與波谷表面作為一個整體以恒定的速度向下沉降。值得注意的是，直到計算結束，液體表面最低點到底面固壁的距離仍然保持在以上，因此，圖5中所示的所有穩(wěn)態(tài)動力學特性與底面的固壁邊界條件無關。

圖6 與無量綱參數(shù)A*的關系Fig.6 Dependence ofandη*on dimensionless parameter A*

式中：c1＝1．26。

2．2 穩(wěn)態(tài)流動的動力學特性

2.2.1 波谷表面的恒定沉降速度V*

從圖4（b）和圖5可以看出，在非線性發(fā)展階段，整個波谷表面都以相同的速度 V*向下沉降。這個速度就是RT不穩(wěn)定性研究中提及的氣泡推進速度。圖7給出了模擬工況A*＝3.2π3，4.8π3，6.4π3和8.0π3下波谷最低點的無量綱豎直方向速度隨時間的變化規(guī)律?？梢钥吹剿心M工況下的值最終都會穩(wěn)定在一個常數(shù)，即

圖7 不同工況下隨時間的變化規(guī)律Fig.7 Temporal changes ofunder different cases

Baker等［26］利用渦層方法得到了與式（6）相同的數(shù)值。漸進勢流模型理論［14］給出At＝1的條件下二維氣泡頂點推進速度為。這個理論結果與式（6）也吻合得很好。從第2.2節(jié)的討論將會看到，V表征了液體通過最大壓力點進入自由尖釘區(qū)域的速度，因此其是研究RT不穩(wěn)定性誘導霧化過程中的一個重要參數(shù)。本節(jié)將基于數(shù)值計算得到的詳細壓力場和速度場信息，對RT不穩(wěn)定性非線性發(fā)展的物理過程進行分析，并建立新的模型計算波谷表面沉降速度V。

沿著波谷中心線（x*＝0．75），由于對稱性，水平方向的速度u*＝0，因此y方向上液體區(qū)域的無量綱動量方程可表示為

圖8（a）給出了模擬工況A*＝4.8π3下非線性發(fā)展階段不同時刻t*＝3.41，3.66，3.90，4.14沿波谷中心線的無量綱壓強分布規(guī)律?？梢钥吹?，隨著波谷表面的沉降，無量綱壓強相應地增大。當把圖8（a）中的縱坐標y*替換為之后（見圖8（b）），可以看到所有的壓力曲線幾乎重合在一起，表明隨時間變化的壓力函數(shù)p*（x*，y*，t*）在非線性發(fā)展階段可以寫為2個變量的函數(shù)p*（x*，ζ*）。另外，從圖4和圖5已經(jīng)知道，在非線性發(fā)展階段，整個正弦表面（包括正弦波谷和虛正弦波峰）以及最大壓力點都以相同的穩(wěn)定速度沉降。當站在波谷表面最低點來觀察時，最大壓力點以下的液體流動為穩(wěn)態(tài)流動，因此豎直方向速度可以表示為v*（x*，y*，t*）＝v*（x*，ζ*），則式（7）可以化簡為

圖8 不同時刻波谷中心線無量綱壓強的分布規(guī)律Fig.8 Dimensionless pressure distribution of trough center line at different time

圖9給出了模擬工況A*＝4.8π3下非線性發(fā)展階段不同時刻t*＝3.41，3.66，3.90，4.14沿波峰中心線（x*＝0．25）的無量綱壓強分布規(guī)律?？梢钥吹?，除了尖釘頂端區(qū)域由于頂部收縮引起的較大表面壓強外，在最大壓力點（實心黑點）以上液體尖釘區(qū)域內(nèi)的壓力會緩和到與周圍氣體壓力相一致。而由于本文中的Atwood數(shù)接近于1，液-氣密度比較大，因此在這個液體尖釘區(qū)域內(nèi)的壓強p*及其豎直方向梯度?p*／?y*（＝?p*／?ζ*）都接近于0。比較圖9（a）與圖8（a）可以看到，在接近底部壁面區(qū)域，沿波峰中心線的壓強分布與沿波谷中心線的壓強分布一致，即

式（7）和式（8）在波峰中心線上最大壓力點以下區(qū)域仍然成立。將式（8）沿波峰中心線對變量ζ*從（對應y*＝0）到（對應）進行積分，得

圖9 不同時刻波峰中心線處無量綱壓強的分布規(guī)律Fig.9 Dimensionless pressure distribution of crest center line at different time

另外，當考慮到非線性發(fā)展階段虛正弦波峰表面與正弦波谷表面以一個整體向下運動，由于質量守恒，則在虛正弦波峰頂點處的豎直方向速度應滿足條件：

計算結果（見圖10）也證實了這個條件。

圖10 不同工況下隨時間的變化規(guī)律Fig.10 Temporal changes ofunder different cases

將式（10）、式（12）～式（14）代入式（11），得

數(shù)值結果（見圖6）顯示δ*0和無量綱參數(shù)A*無關，而且保持在恒定值0.125。筆者將在第2.2.2節(jié)中對如何推導得到δ*0進行討論。從式（15）得到的表面恒定沉降速度V*＝0.236，此數(shù)值與直接模擬得到的結果（見式（6））以及基于漸進勢流模型的理論研究結果［14］吻合得很好。

除了波谷表面的恒定沉降速度V*之外，另外一個與穩(wěn)定流動相關的重要特征量就是數(shù)值與無量綱參數(shù)A*無關的正弦表面振幅，其表征了液體尖釘根部、波谷表面附近所形成穩(wěn)定流的厚度。本節(jié)將根據(jù)壓力場和速度場的計算結果建立模型，推導出的表達式。

在RT不穩(wěn)定性中，液體表面的非線性運動學條件為

式中：h*（x*，t*）為表面高度函數(shù)；和分別為表面速度的水平和豎直方向的分量。

在線性理論中，式（16）等號右側第2項非線性項會被忽略。在線性發(fā)展階段，由于表面各點速度的水平分量和表面梯度?h*／?x*都很小，因此這種忽略是合理的。但在非線性發(fā)展階段，由于與?h*／?x*的乘積不再是一個小值，因而在分析中必需考慮此項的影響。

觀察組男性33例，女性27例，平均年齡為（65.95±5.89）歲；疾病類型：7例為膽囊炎患者，19例為闌尾炎患者，14例為腹膜炎患者，12例為下肢靜脈曲張患者，8例為乳腺纖維瘤患者。對照組男性32例，女性28例，平均年齡為（65.16±5.54）歲；疾病類型：6例為膽囊炎患者，20例為闌尾炎患者，15例為腹膜炎患者，11例為下肢靜脈曲張患者，8例為乳腺纖維瘤患者。兩組患者一般資料差異不具有統(tǒng)計學意義（P＞0.05）。

圖11給出了模擬工況A*＝4.8π3下表面上對應點的x坐標隨時間的變化規(guī)律?？梢钥吹剑谡麄€發(fā)展過程中始終是在中性點上。如前所述，在非線性發(fā)展階段，波谷表面可以近似成振幅為的正弦曲線，因此，波谷區(qū)域（x*≥0．5）的表面高度可表示為在中性點上，式（16）可化簡為

在模擬工況A*＝4.8π3下，如圖12所示和位于同一水平面上波峰中心線位置的豎直方向速度大小相同。此關系對其他模擬工況同樣成立。

如圖9所示，虛正弦波峰頂部附近（實心三角形）豎直方向的壓力梯度，因此，沿波峰中心線（x*＝0．25）在和的范圍內(nèi)，式（11）中等號右側的壓力項-?p*／?ζ*可以忽略，即

將式（14）和式（15）代入式（19），得

其解為

圖11 表面不動點所對應的x坐標隨時間的變化規(guī)律Fig.11 Temporal changes of x coordinate where vertical velocity vanishes on surface

圖12 與隨時間的變化規(guī)律Fig.12 Temporal changes ofand

2.2.3 最大壓力點處出流速度v*3和出流寬度2b*

從第2.2.1節(jié)可以看出，在RT不穩(wěn)定性中，尖釘?shù)撞啃纬傻淖畲髩毫c將液體分為2個區(qū)域：位于最大壓力點以下的穩(wěn)態(tài)流動區(qū)和位于最大壓力點以上的自由流動區(qū)。其中，自由流動區(qū)可以看成是在重力作用下液體從噴嘴向下噴出的低速射流問題，而這個最大壓力點的位置可以看成是射流問題中噴嘴出口的位置。因此，如果可以確定最大壓力點處的出流速度和出流寬度2b*，則射流問題中的研究成果［27-30］就能應用到RT不穩(wěn)定性誘導霧化領域來。這是以往研究中所忽略的問題。因此，本節(jié)將分析如何確定最大壓力點處的出流速度和出流寬度2b*。

第2.2.2節(jié)中討論過在虛正弦表面波峰頂點附近的豎直方向壓力梯度可以忽略，因此，從式（18）可以得到在虛正弦表面波峰頂點處（x*＝）的速度梯度?v*／?y*（＝?v*／?ζ*）：

將式（5）代入式（24），得

從式（5）和式（25）可以看出，當A*→∞，且。因此，當無量綱參數(shù)A*足夠大時，虛正弦表面波峰頂點和最大壓力點重合，出流速度為2V*。

圖13給出了最大壓力點處出流半寬b*以及虛正弦表面波峰頂點處出流半寬a*與無量綱參數(shù)A*的函數(shù)關系?？梢钥吹剑琣*與A*的取值無關，一直保持為常數(shù)c3＝0．173，而b*則隨著A*的增大而增大。同樣，如果A*足夠大，則

圖13 a*、b*與無量綱參數(shù)A*的關系Fig.13 Dependence of a*and b*on dimensionless parameter A*

3 結 論

本文利用CLSVOF界面捕捉方法對大Atwood數(shù)下二維單模態(tài)RT不穩(wěn)定性進行了數(shù)值模擬，重點研究了RT不穩(wěn)定性尖釘區(qū)域非線性發(fā)展階段的穩(wěn)態(tài)動力學特性。通過對計算結果的分析可以得到如下結論：

1）在液體尖釘根部由于債性力作用而引起的水平?jīng)_擊流會在此處形成一個局部最大壓力點，這個最大壓力點將使位于其上部的尖釘區(qū)域的發(fā)展獨立于位于其下部的液層。由于此處債性力與壓強梯度的平衡，位于最大壓力點附近的流動最終將達到穩(wěn)態(tài)。

2）在穩(wěn)態(tài)流動中，液體波谷表面可用一條正弦曲線表征，其振幅只與初始擾動的波長有關。波谷表面連同最大壓力點將作為一個整體以恒定的速度向下運動，其數(shù)值與債性加速度和波長乘積的平方根成正比。

3）得到了最大壓力點處的出流速度和出流寬度與無量綱參數(shù)的關系。通過確定這些穩(wěn)態(tài)特征量有助于將經(jīng)典低速射流的相關理論擴展應用到RT不穩(wěn)定性誘導霧化的研究領域。

（References）

［1］ARNETTW D，BAHCALL JN，KIRSHNER R P，et al.Supernova 1987A［J］.Annual Review of Astronomy and Astrophysics，1989，27（2）：629-700.

［2］NORMAN M L，SMARR L，SMITH M D，et al.Hydrodynamic formation of twin-exhaust jets［J］.Astrophysical Journal，1981，247（1）：52-58.

［3］EVANSR G，BENNETT A J，PERT G J.Rayleigh-Taylor instabilities in laser accelerated targets［J］.Physical Review Letters，1982，49（22）：1639-1642.

［4］LINDL JD，MCCRORY R L，CAMPBELL EM.Progress toward ignition and burn propagation in inertial confinement fusion［J］.Physics Today，1992，45（9）：32-40.

［5］BEALE J C.Modeling spray atomization with the Kelvin-Helmholtz／Rayleigh Taylor hybrid model［J］.Atomization and Sprays，1999，9（6）：623-650.

［6］KONG SC，SENECAL P K，REITZ R D.Developments in spray modeling in diesel and direct-injection gasoline engines［J］. Oil＆Gas Science＆Technology，1999，54（2）：197-204.

［7］HSIANG L P，F(xiàn)AETH G M.Near-limit drop deformation and secondary breakup［J］.International Journal of Multiphase Flow，1992，18（5）：635-652.

［8］LEE C H，REITZ R D.An experimental study of the effect of gas density on the distortion and breakupmechanism of drops in high speed gas stream［J］.International Journal of Multiphase Flow，2000，26（2）：229-244.

［9］解茂昭.燃油噴霧場結構和霧化機理［J］.力學與實踐，1990，12（4）：9-15.

XIE M Z.The structure of fuel spray field and themechanism of atomization［J］.Mechanics in Engineering，1990，12（4）：9-15（in Chinese）.

［10］TAYLOR G.The instability of liquid surfaces when accelerated in a direction perpendicular to their planes.I［J］.Proceedings of the Royal Society of London Series A，Mathematical and Physical Sciences，1950，201（1065）：192-196.

［11］BELLMAN R，PENNINGTON R H.Effects of surface tension and viscosity on Taylor instability［J］.Quarterly of App liedMathematics，1953，12（2）：151-162.

［12］PIRIZ A R，CORTáZAR O D，CELA JJL，et al.The Rayleigh-Taylor instability［J］.American Journal of Physics，2006，74（12）：1095-1098.

［13］SHARP D H.An overview of Rayleigh-Taylor instability［J］. Physica D：Nonlinear Phenomena，1984，12（1-3）：3-10.

［14］GONCHAROV V N.Analytical model of nonlinear，singlemode，classical Rayleigh-Taylor instability at arbitrary Atwood numbers［J］.Physical Review Letters，2002，88（13）：134502.

［15］BAKER G R，MCCRORY R L，VERDON C P，et al.Rayleigh-Taylor instability of fluid layers［J］.Journal of Fluid Mechanics，1987，178：161-175.

［16］RAMAPRABHU P，DIMONTE G，WOODWARD P，et al.The late-time dynam ics of the single-mode Rayleigh-Taylor instability［J］.Physics of Fluids，2012，24（7）：074107.

［17］RAMAPRABHU P，DIMONTE G，YOUNG Y N，et al.Limits of the potential flow approach to the single-mode Rayleigh-Taylor problem［J］.Physical Review E，2006，74（6）：202-212.

［18］葉文華.激光燒蝕RT不穩(wěn)定性線性增長率和非線性行為的數(shù)值研究［J］.強激光與粒子束，1998，10（4）：567-572.

YE W H.Numerical studies of linear growth rates and nonlinear evolution of laser ablative Rayleigh-Taylor instability［J］.High Power Laser and Particle Beams，1998，10（4）：567-572（in Chinese）.

［19］葉文華，張維巖，陳光南，等.激光燒蝕瑞利-泰勒不穩(wěn)定性數(shù)值研究［J］.強激光與粒子束，1999，11（5）：613-618.

YE W H，ZHANG W Y，CHEN G N，et al.Numerical study of laser ablative Rayleigh-Taylor instability［J］.High Power Laser and Particle Beams，1999，11（5）：613-618（in Chinese）.

［20］葉文華，張維巖，賀賢土.燒蝕瑞利-泰勒不穩(wěn)定性線性增長率的預熱致穩(wěn)公式［J］.物理學報，2000，49（4）：762-767.

YEW H，ZHANG W Y，HE X T.Preheating stabilization formula of linear growth rate for ablative Rayleigh-Taylor instability［J］.Acta Physica Sinica，2000，49（4）：762-767（in Chinese）.

［21］程會方，段日強，姜勝耀.Rayleigh-Taylor不穩(wěn)定性的MPS數(shù)值模擬［J］.核動力工程，2010，31（s1）：123-126.

CHENG H F，DUAN R Q，JIANG SY.Numerical simulation of Rayleigh-Taylor instability with MPSmethod［J］.Nuclear Power Engineering，2010，31（s1）：123-126（in Chinese）.

［22］BRACKBILL J U，KOTHE D B，ZEMACH C.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100（2）：335-354.

［23］SUSSMAN M，PUCKETT E G.A coupled level set and volumeof-fluid method for computing 3D and axisymmetric incompressible two-phase flows［J］.Journal of Computational Physics，2000，162（2）：301-337.

［24］VAN DER PIJL S P，SEGAL A，VUIK C，et al.A mass-conserving level-setmethod formodeling ofmulti-phase flows［J］.International Journal for Numerical Methods in Fluids，2005，47（4）：339-361.

［25］LI Y，UMEMURA A.Two-dimensional numerical investigation on the dynamics of ligament formation by Faraday instability［J］.International Journal of Multiphase Flow，2014，60（2）：64-75.

［26］BAKER G R，MEIRON D I，ORSZAG S A.Vortex simulations of the Rayleigh-Taylor instability［J］.Physics of Fluids，1980，23（8）：1485-1490.

［27］CLANET C，LASHERAS JC.Transition from dripping to jetting［J］.Journal of Fluid Mechanics，1999，383：307-326.

［28］EGGERS J，VILLERMAUX E.Physics of liquid jets［J］.Reports on Progress in Physics，2008，71（3）：509-514.

［29］SCHULKES R M SM.The evolution and bifurcation of a pendant drop［J］.Journal of Fluid Mechanics，1994，278：83-100.

［30］UMEMURA A.Self-destabilizingmechanism of a laminar inviscid liquid jet issuing from a circular nozzle［J］.Physical Review E，2011，83（4）：046307.

Tel.：13811862628

E-mail：liyikai＠bit.edu.cn

Num erical study on non linear characteristics of Rayleigh-Taylor instability

KANG Ning，LIYikai*，HE Xu

（School of Mechanical Engineering，Beijing Institute of Technology，Beijing 100081，China）

The research on the nonlinear dynam ics of Rayleigh-Taylor（RT）beforemainly focused on deducing and measuring the constant penetration velocity of the bubble and had little detailed analysis of the nonlinear dynam ic characteristics in the liquid spike region.An accurate numerical simulation of the single-mode RT instability was carried out based on the coupled Level-Set and Volume-of-Fluid（CLSVOF）interface capturing method.The detailed information on the pressure fields and velocity fields was obtained.In addition，the steady-state dynamic characteristics in the nonlinear development stage were analyzed.Simulation results show that a localmaximum pressure point which is caused by the horizontal impinging flow with the action of inertial force appears at the root of the spike.The dependence of the different characteristic parameters of the steady flow on the initial perturbation wavelength and the inertial acceleration is determ ined.This work may extend the relevant classical theories of the low speed jet to the RT instability inducing atomization field.

Rayleigh-Taylor instability；nonlinearity；steady state；CLSVOFmethod；numerical simulation

2015-10-15；Accep ted：2016-01-22；Pub lished online：2016-02-18 11：20

s：National Natural Science Foundation of China（51476011）；Beijing Institute of Technology Research Fund Program for Young Scholars（3030012261599）

O363.2

A

1001-5965（2016）10-2059-10

康寧 男，博士研究生。主要研究方向：內(nèi)燃機流動與燃燒。E-mail：46054832＠qq.com

黎一鍇 男，博士，講師。主要研究方向：氣-液兩相流的數(shù)值模擬。

http：∥bhxb.buaa.edu.cn jbuaa＠buaa.edu.cn

DO I：10.13700／j.bh.1001-5965.2015.0667

2015-10-15；錄用日期：2016-01-22；網(wǎng)絡出版時間：2016-02-18 11：20

www.cnki.net／kcms／detail／11.2625.V.20160218.1120.003.htm l

國家自然科學基金（51476011）；北京理工大學青年教師學術啟動計劃（3030012261599）

*通訊作者：Tel.：13811862628 E-mail：liyikai＠bit.edu.cn

康寧，黎一鍇，何旭.Rayleigh-Taylor不穩(wěn)定性非線性特性的數(shù)值研究［J］.北京航空航天大學學報，2016，42（10）：2059-2068.KANG N，LIY K，HE X.Numerical study on nonlinear characteristics ofRayleigh-Taylor instability［J］.Journal of Beijing University of Aeronautics and Astronautics，2016，42（10）：2059-2068（in Chinese）.

URL：www.cnki.net／kcms／detail／11.2625.V.20160218.1120.003.htm l

*Correspond ing au thor.Tel.：13811862628 E-mail：liyikai＠bit.edu.cn