SPH方法對氣液兩相流自由界面運動的追蹤模擬

沈雁鳴， 陳堅強

（中國空氣動力研究與發(fā)展中心，四川 綿陽621000）

0 引 言

空氣動力學(xué)和水動力學(xué)是流體力學(xué)中兩個重要分支，二者研究對象在可壓縮性、密度、壓力等方面具有巨大差異，因此在以往的研究中聯(lián)系并不緊密。但是，隨著時代的發(fā)展和技術(shù)的進(jìn)步，二者的交叉領(lǐng)域卻越來越受到關(guān)注，其中以氣液兩相流動最為突出。氣液兩相流動是自然界中常見的流動現(xiàn)象，風(fēng)浪，風(fēng)雨、水泡、汽水等等是典型的氣液兩相流，這其中涉及到多相介質(zhì)作用、自由界面運動、粘性力和表面張力作用等問題。這些流體力學(xué)問題也是目前流體力學(xué)界關(guān)注的熱點和難點。深入研究氣液兩相流動的機理，了解多相介質(zhì)界面運動的數(shù)學(xué)和物理規(guī)律，對解決船舶工程、海洋工程、化學(xué)工程、水利工程中的實際問題，具有十分廣泛的現(xiàn)實意義。

由于氣液兩相流動中氣體和液體間相互作用使得氣液界面產(chǎn)生變形，因此進(jìn)行數(shù)值模擬時要求必須能夠準(zhǔn)確跟蹤變形界面。目前學(xué)術(shù)界已經(jīng)形成了多種基于網(wǎng)格的成熟的計算方法，比較典型的有MAC法、VOF方法、Level－set方法。這些方法均取得了較好的應(yīng)用效果，但是由于它們都是基于網(wǎng)格劃分的，在處理大變形問題時需要對網(wǎng)格劃分等過程進(jìn)行復(fù)雜操作，不方便實施，難以處理網(wǎng)格扭曲畸變現(xiàn)象，而且在處理多介質(zhì)問題時還需要進(jìn)行繁瑣的示標(biāo)記錄過程。近年來，隨著無網(wǎng)格方法的發(fā)展和應(yīng)用的日趨成熟，國際計算物理學(xué)界開始將無網(wǎng)格方法應(yīng)用到多介質(zhì)自由界面流動的數(shù)值模擬中。這類方法簡單形象，穩(wěn)定有效，在一定的改進(jìn)后能夠高精度、高效率地模擬自由界面流動，光滑粒子流體動力學(xué)方法（Smoothed Particle Hydrodynamics，SPH）就是一系列無網(wǎng)格方法中的一種。

光滑粒子流體動力學(xué)方法是一種用于求解連續(xù)介質(zhì)動力學(xué)的無網(wǎng)格方法，其核心思想來源于插值理論，即通過一種稱之為核函數(shù)的積分核進(jìn)行估值近似，從而將流體力學(xué)方程轉(zhuǎn)化為數(shù)值計算用的SPH方程組。在整個流場中，流場介質(zhì)被離散為一系列粒子，這些粒子攜帶所有的流場變量，受控制方程約束在空間中進(jìn)行一定規(guī)律的運動。因此，算法本身受應(yīng)用限制較少，理論上可以適用于任意變形問題。該方法在1977年由Gingold ＆Manahan和Lucy分別提出［1］，經(jīng)過不斷完善與改進(jìn)，已被成功運用于天體物理、流體力學(xué)、材料力學(xué)等各領(lǐng)域的計算研究中［2］。本文即利用該方法，結(jié)合微可壓縮模型（SCM），并利用了界面控制方法XSPH速度修正以及Van Der waals狀態(tài)方程修正，對典型的氣液兩相流動如二維潰壩、氣泡上浮等問題進(jìn)行了數(shù)值模擬，分析了空氣和水相互作用機理以及界面運動規(guī)律，結(jié)果表明該方法在模擬多相介質(zhì)界面運動問題準(zhǔn)確有效，可用于處理更為復(fù)雜的多相流動工程問題。

1 數(shù)學(xué)模型

1.1 SPH方法

光滑粒子動力學(xué)（SPH）方法的核心是插值［3］，宏觀變量A（r）能方便地借助于一組無序點上的值表示成積分插值計算得到，其形式為：

假設(shè)質(zhì)點所表示的物體密度為ρ（r′），則式（1）可以進(jìn)一步表示為：

積分往往可以近似為有限個數(shù)目的部分之和，因此，

這樣式（2）可以近似為：

利用核的可微分性質(zhì)很容易計算得到〈▽A〉，〈▽·A〉，〈▽×A〉的表達(dá)式，此處A（r）可以取為ρ（r），P（r）或者v（r）。流體力學(xué)控制方程 N－S方程可通過此方法作變換可得到以下方程組：

其中W為光滑核函數(shù)，為剪切應(yīng)力張量為剪切應(yīng)變率張量為總應(yīng)力張量，Πij為人工粘性，Hi為人工熱量，vij為粒子之間速度差。粒子的密度計算中，只考慮本相粒子相互間的插值。而動量和能量計算中可使用作用域中所有粒子進(jìn)行插值。

對于核函數(shù)W，它是兩粒子之間的距離與h的函數(shù)，其中h是核函數(shù)影響域的尺寸，本文使用Gaussian形式：

式中S＝r／h，d表示維度。

1.2 微可壓縮模型（SCM）

SPH方法是用于模擬可壓縮流動問題的，因此使用該方法模擬不可壓縮流動時，必須使用一定的方法保證可壓縮條件的實現(xiàn)。鄧小剛［4］在處理氣體低速流動問題時提出了微可壓縮模型（SCM），在此我們將該模型應(yīng)用于液體自由界面流動中，把不可壓縮流動考慮為微可壓縮流。應(yīng)用這種方法可以有效地緩解壓力計算對密度變化敏感所導(dǎo)致的計算不穩(wěn)定。

上式中ρ0為參考密度，一般取初始密度值，P0為參考壓力，對應(yīng)初始密度下的壓力值。γ在氣體時取常數(shù)值1.4，對液體水則取常數(shù)值7，這可以使微小的密度變化產(chǎn)生較大的壓力變化，從而保證密度波動保持在很小的范圍之內(nèi)。該壓力表達(dá)式使低速空氣和低速液體在自由界面流動中狀態(tài)方程實現(xiàn)了統(tǒng)一。

對于自由流動的流體，在只有重力的作用下，其表面初始密度值為ρ0，但靜水壓力值等于0。為了保證此條件的實現(xiàn)，可將（9）式變換為下面形式：

該方程將自由流動的初始條件添加到微可壓縮模型中。參考壓力值P0也只是一個常數(shù)，用于限定密度的最大改變量，為了保證密度波動不至于過大，可以取值為P0＝ρ0c2／γ。其中c為參考聲速，一般取值為c＝10vmax。

然而，重力場會對液體從上到下產(chǎn)生與高度相關(guān)的靜水壓力，在計算開始后重力的突然施加會導(dǎo)致液體的彈性振蕩，為了克服這個問題，可以給出初始密度隨高度變化的函數(shù)式［5］：

在初始狀態(tài)時通過該方式對不同高度液體密度進(jìn)行賦值，可以有效地克服因重力施加導(dǎo)致的液面高度振蕩。

1.3 固壁邊界條件與界面控制方法

對于固壁條件，在Liu［3］等的研究中，他們使用了兩種虛粒子來處理固壁邊界條件。其中第一種虛粒子設(shè)置于固定邊界，與 Monaghan［5］所使用的相似。該粒子用于對臨近邊界的粒子施加Lennard－Jones排斥力防止流體粒子穿透邊界，排斥力表達(dá)式如下：

本文對容器壁的處理使用了該方法。同時，為了更好地模擬出壁面邊界的無滑移條件，防止密度計算過程中的數(shù)值振蕩，本文還按照 Morris［6］提出的方法，在固壁邊界布置與流體粒子相同的虛粒子，參與流體粒子密度和壓力的計算，但這種虛粒子并不代表固體內(nèi)部的真實特性。

同時，為了防止計算中粒子出現(xiàn)非物理穿透，本文還使用了XSPH速度修正［7］以及Van Der waals狀態(tài)方程修正［7］。XSPH速度修正表達(dá)式為

其中＝0.5（ρi＋ρj）。該修正式可以使交界面更為清晰光滑。

Van Der waals狀態(tài)方程修正式為

上式右邊項即為Van Der waals修正項，一般只在多相流動中氣態(tài)相中添加，可以有效地控制交界面形狀，防止氣體粒子發(fā)生崩潰或者逃逸。

2 結(jié)果與分析

通過對數(shù)學(xué)模型和計算方法的理解和研究，本文使用光滑粒子流體動力學(xué)方法對常見的自由界面流動二維潰壩及氣泡上浮問題進(jìn)行了計算，比較并分析了計算結(jié)果。

2.1 潰壩問題

潰壩流動是一種典型的自由表面流動，該流動包含了水體與壁面作用產(chǎn)生的水花飛濺融合、自由表面變形等多種復(fù)雜現(xiàn)象，完整而準(zhǔn)確地模擬這些復(fù)雜現(xiàn)象在傳統(tǒng)網(wǎng)格方法中是較為困難的。

本文的計算模型如圖1所示：初始時刻水柱高度H為1m，水面寬度L為2m。整個水箱高度D為3m，長度d為5.366m，當(dāng)有空氣存在時水箱中其他位置全部充滿空氣。重力加速度g為9.8m／s2，水的動力粘性系數(shù)為0.01Pa·s，空氣為0.0001Pa·s。水柱在重力的作用下垮塌，沿著水平固壁邊界流動直到遇到前方3.366m處垂直固壁后發(fā)生撞擊。無空氣存在時水粒子數(shù)為5000個，有空氣存在時空氣粒子數(shù)為35200個。粒子間距Δx為0.02m，光滑長度h為1.33Δx，計算時間步長 Δt為0.00001s。初始密度按照公式（11）給出，代入公式（10）便可得到重力場下的初始壓力分布。狀態(tài)方程中水的聲速值為14（gH）0.5，空氣的聲速值為198（gH）0.5。

圖2是兩種模型下不同時刻水柱變形數(shù)值模擬結(jié)果，其中上面三幅圖是無空氣存在條件下0.4s、1.8s和2.2s時水面變化，下面三幅圖則為有空氣存在時對應(yīng)時刻氣水兩相自由界面變化。由圖中可以看出，計算開始后水柱在重力的作用下迅速垮塌，浪頭沿著水平固壁迅速推進(jìn)，直到遇到垂直固壁后產(chǎn)生猛烈撞擊，水頭沿著垂直壁面向上爬升并向外翻滾，翻滾產(chǎn)生的浪頭同后部水面再次發(fā)生拍擊產(chǎn)生二次浪花飛濺。通過上下對比也可以發(fā)現(xiàn)，空氣的存在一定程度上阻礙了水流的激烈變化，使得自由界面變化更為平緩一些。

圖1 潰壩流動計算模型示意圖Fig.1 Model of dam－break problem

圖2 計算開始后兩種模型不同時刻水柱變形圖Fig.2 Dam－break flow and impact against a vertical wall

為了進(jìn)一步觀察氣水兩相界面狀況，浪頭撞擊墻壁后的翻滾情況被放大后在圖3中展示。由圖中可以看出，浪頭翻滾后將一部分空氣包裹成氣泡，而空氣和水的交界面連續(xù)性和光滑度都較好。本文通過界面控制方法和狀態(tài)方程的修正，使得光滑粒子流體動力學(xué)方法可以生動地模擬出潰壩問題后浪花的翻滾、飛濺以及融合過程，對氣水兩相自由界面的追蹤模擬也較為有效。

圖3 有空氣存在時浪頭撞擊墻壁后的翻滾情況Fig.3 Detail of water impacting vertical wall

本文也給出了定量上的比較以驗證計算的準(zhǔn)確性。如圖4所示，為水頭前沿離左側(cè)固壁距離與時間關(guān)系實驗值和計算結(jié)果的比較。由圖中可以看出，計算結(jié)果中水頭的前進(jìn)距離比實驗結(jié)果略大一些，空氣的存在增加了水流的阻力，因此水頭前進(jìn)速度比無空氣時慢。SPH方法對潰壩問題的模擬在定性和定量上都是較為準(zhǔn)確有效。

2.2 氣泡上浮

氣泡上浮是另外一種典型的自由表面流動，圓形氣泡在浮力的作用下上升，上升過程中氣泡逐漸變形甚至破裂。這一過程中流體粘性力、壓阻、表面張力作用過程較為復(fù)雜。在此我們使用光滑粒子流體動力學(xué)方法對此問題進(jìn)行模擬。

圖4 潰壩問題水頭前沿離左側(cè)固壁距離與時間關(guān)系Fig.4 Time evolution of the water－front

計算模型如圖5所示：初始時刻，半徑R為0.03m的氣泡，在寬為6R，高為10R的水中靜止不動，氣泡離底面高度2R，計算開始后氣泡在浮力的作用下上浮并逐漸變形。計算中空氣與水的其他參數(shù)與潰壩模型相同。粒子間距Δx為0.0015m，光滑長度h為1.33Δx，計算時間步長Δt為0.00001s，計算中未考慮表面張力作用。

圖5 氣泡上浮計算模型示意圖Fig.5 Model of bubble rising

圖6是氣泡上浮過程不同時刻界面變化的計算結(jié)果，其中左邊是0.7s時氣泡在水中的位置以及形狀，右邊則為氣泡最初狀態(tài)、0.2s、0.5s、0.7s直到0.9s時氣泡形狀計算結(jié)果。由圖中可以看出，氣泡上浮過程中，由于頭部阻力影響，圓形氣泡開始扁平，尾部出現(xiàn)凹陷，隨后凹陷不斷加劇，并朝兩端擴展，發(fā)展出蘑菇狀凹陷，整個氣泡則變?yōu)轳R鞍形，下端兩個鞍腳向內(nèi)翻轉(zhuǎn)，不斷拉長變細(xì)最后斷裂，分離成為兩個小氣泡。頭部則收縮成為一個前拱后平的大氣泡，兩個小氣泡跟隨在其尾部朝液面繼續(xù)上升。

圖7則是使用SPH方法的計算結(jié)果和參考文獻(xiàn)使用Level－set計算結(jié)果的比較。由圖中可以看出，本文計算結(jié)果與level－set計算結(jié)果［6］在氣泡形狀和氣泡上升位置都較為吻合。SPH方法在兩相自由運動界面的追蹤上較為清晰光滑，而Level－Set方法在氣泡分裂過程中出現(xiàn)了幾個小氣泡。

圖6 氣泡上浮過程不同時刻界面變化Fig.6 Time evolution of bubble deformation

圖7 SPH計算結(jié)果與參考文獻(xiàn)Level－Set計算結(jié)果比較Fig.7 Rising－bubble deformation，the present SPH method results compared with the Level－Set solution

圖8是氣泡頭部位置及速度隨時間變化曲線。由氣泡頭部位移變化曲線看來，氣泡頭部在0.3s之前做加速上升，0.3s以后基本做線性運動。由圖上速度曲線可以看出，在0.2s以前氣泡由漸變?yōu)闄E圓的過程中，氣泡上升速度慢慢提高，但增幅減小。0.2s以后隨著氣泡尾部凹陷，氣泡阻力增大，逐漸與浮力和重力達(dá)到平衡，上升速度穩(wěn)定在1m／s左右。

圖8 氣泡頭部位置及速度隨時間變化曲線Fig.8 Time evolutions of position and velocity of bubble head

3 結(jié) 論

本文利用光滑粒子流體動力學(xué)方法，結(jié)合微可壓縮模型（SCM），對典型的氣液兩相流動如二維潰壩、氣泡上浮等問題進(jìn)行了數(shù)值模擬。為了克服計算振蕩，防止界面崩潰，本文還利用了界面控制方法XSPH速度修正以及Van Der waals狀態(tài)方程修正。模擬結(jié)果和實驗結(jié)果以及Level－Set方法計算結(jié)果進(jìn)行了比較，分析后得到了以下結(jié)論：

潰壩問題中，空氣的存在一定程度上阻礙了水流的激烈變化，使得自由界面變化更為平緩一些。因此在當(dāng)水氣界面相對速度較大時，不能只考慮水介質(zhì)，空氣的存在產(chǎn)生的氣墊效應(yīng)也必須考慮。

氣泡上浮問題中，氣泡上升過程不斷變形，界面發(fā)生不斷彎曲變形直至斷裂，斷裂后又可以重新組成新的氣泡單元，這一界面的變化過程是較為符合物理實際的。

因此，本文所用光滑粒子流體動力學(xué)方法和微可壓縮模型（SCM）在模擬氣液兩相介質(zhì)界面運動問題中十分有效，所引入XSPH速度修正以及Van Der waals狀態(tài)方程修正能夠有效克服計算振蕩，控制界面變化。通過進(jìn)一步的改進(jìn)，可利用該方法處理更為復(fù)雜的多相流動工程問題。

［1］GINGOLD R A ，MONAGHAN J J.Smoothed particle hydrodynamics：theory and application to non－spherical stars.monthly noticesr［J］.AstronomySoc.，1977.181：375－389.

［2］MONAGHAN J J.An introduction to SPH［J］.ComputerPhysicsCommunica－tions，1988，48：89－96.

［3］LIU G R，LIU M B.Smoothed Particle hydrodynamics：a mesh－free particle method［M］. World Scientific，2003.

［4］DENG X G，ZHUANG F G，MAO M L.On low mach number perfect gas flow calculations［R］.AIAA Paper 99－3317，1999.

［5］MONAGHAN J J.Simulating free surface flows with SPH［J］.JournalofComputationalPhysics，1994，110：399－406.

［6］MORRIS J P，MONAGHAN J J.A switch to reduce SPH viscosity［J］.JournalofComputationalPhysics，1997，136：41－50.

［7］COLAGROSSI A，LANDRINI M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics［J］.JournalofComputationalPhysics，2003，191：448－475.

［8］OGER G，DORING M，ALESSANDRINI B，et al.Twodimensional SPH simulations of wedge water entries［J］.JournalofComputationalPhysics，213（2006）：803－822.