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

康 寧，黎一鍇，何 旭

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

康 寧，黎一鍇*，何 旭

（北京理工大學(xué) 機(jī)械與車輛學(xué)院，北京100081）

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

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

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

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

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

由于非線性在數(shù)學(xué)處理上的復(fù)雜性，雖然可以推導(dǎo)出局部個(gè)別參數(shù)（如氣泡的恒定推進(jìn)速度）的漸進(jìn)解，但是很難得到描述全局范圍動(dòng)力學(xué)特性的解析解（如流體內(nèi)的壓力場(chǎng)和速度場(chǎng)分布）。近年來，隨著計(jì)算機(jī)計(jì)算速度以及計(jì)算方法的不斷發(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ā)展特性。在國(guó)內(nèi)，葉文華等［18-20］利用高精度的通量校正傳輸（Flux Corrected Transport，F(xiàn)CT）方法對(duì)激光燒蝕RT不穩(wěn)定性進(jìn)行了模擬，研究了電子熱傳導(dǎo)燒蝕在長(zhǎng)波長(zhǎng)擾動(dòng)的非線性RT不穩(wěn)定性演變中所起的作用。程會(huì)方等［21］利用移動(dòng)粒子半隱式法（Moving-Particle Sem i-implicitmethod，MPS）模擬了二維RT不穩(wěn)定性，得到的界面形狀與實(shí)驗(yàn)結(jié)果在表觀上基本一致。

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

1 物理模型和計(jì)算方法

1．1 控制方程

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

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

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

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

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

1．2 定解條件

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

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

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

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

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

1．3 計(jì)算方法

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

1．4 數(shù)值驗(yàn)證

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

圖3 模型驗(yàn)證結(jié)果Fig.3 Results ofmodel validation

2 計(jì)算結(jié)果與討論

從RT不穩(wěn)定性線性增長(zhǎng)率的表達(dá)式可以看到，對(duì)于大Atwood數(shù)，只有在無量綱參數(shù)A*＝ρ1λ2A／σ＞4π2的條件下，表面變形的振幅才能夠增長(zhǎng)。本文對(duì)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個(gè)值進(jìn)行了模擬。所有模擬的初始速度擾動(dòng)振幅為計(jì)算在液體尖釘?shù)捻敳康竭_(dá)計(jì)算域的上邊界后停止。

2．1 非線性表面變形的整體動(dòng)力學(xué)特性

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

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

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

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

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

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

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

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

式中：c1＝1．26。

2．2 穩(wěn)態(tài)流動(dòng)的動(dòng)力學(xué)特性

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

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

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

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

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

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

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

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

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

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

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

計(jì)算結(jié)果（見圖10）也證實(shí)了這個(gè)條件。

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

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

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

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

在RT不穩(wěn)定性中，液體表面的非線性運(yùn)動(dòng)學(xué)條件為

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

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

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

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

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

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

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

其解為

圖11 表面不動(dòng)點(diǎn)所對(duì)應(yīng)的x坐標(biāo)隨時(shí)間的變化規(guī)律Fig.11 Temporal changes of x coordinate where vertical velocity vanishes on surface

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

2.2.3 最大壓力點(diǎn)處出流速度v*3和出流寬度2b*

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

（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］解茂昭.燃油噴霧場(chǎng)結(jié)構(gòu)和霧化機(jī)理［J］.力學(xué)與實(shí)踐，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)定性線性增長(zhǎng)率和非線性行為的數(shù)值研究［J］.強(qiáng)激光與粒子束，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］.強(qiáng)激光與粒子束，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)定性線性增長(zhǎng)率的預(yù)熱致穩(wěn)公式［J］.物理學(xué)報(bào)，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］程會(huì)方，段日強(qiáng)，姜?jiǎng)僖?Rayleigh-Taylor不穩(wěn)定性的MPS數(shù)值模擬［J］.核動(dòng)力工程，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)燃機(jī)流動(dòng)與燃燒。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)絡(luò)出版時(shí)間：2016-02-18 11：20

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

國(guó)家自然科學(xué)基金（51476011）；北京理工大學(xué)青年教師學(xué)術(shù)啟動(dòng)計(jì)劃（3030012261599）

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

康寧，黎一鍇，何旭.Rayleigh-Taylor不穩(wěn)定性非線性特性的數(shù)值研究［J］.北京航空航天大學(xué)學(xué)報(bào)，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