基于LBM方法的圓盤等速入水空泡的數值模擬

張 珂，顏 開，褚學森，程貫一

（中國船舶科學研究中心，江蘇 無錫214082）

1 引 言

自由面流動廣泛存在于自然界和日常生活中，其特點表現為存在可移動的自由表面。有著廣泛的工程應用背景。例如，在風浪環境中航行的船舶，其搖蕩過程中艏艉會與水面發生砰擊，嚴重時所產生的沖擊力可導致船舶局部結構的破壞；船舶液艙中的自由面隨著船身的搖擺發生晃蕩，產生對艙壁的沖擊載荷。惡劣海況中的甲板上浪也屬于自由面流動。物體入水問題是典型的自由面流動問題。空投魚雷在入水過程中，包含撞水(濺水)、流動形成、空泡打開、空泡面閉合、空泡深閉合直至空泡消失等一系列相繼發生的現象，在該過程中強烈的入水載荷有可能引起結構強度破壞問題，也可能影響內部儀器的正常使用。入水空泡將影響入水初期彈道。對入水載荷、入水空泡和入水彈道的研究一直是國內外的重要研究課題。

研究物體入水的數值模擬方法可有多種，如質點標記（MAC）法[1]、VOF方法[2]、光滑粒子流體動力學（SPH）方法[3]、邊界元法[4]等等。這些方法分別具有其自身的優勢和存在的問題。當自由面變化非常劇烈時，尤其是當氣液發生摻混時，自由面邊界變得模糊不清，這些方法往往難以奏效，必須尋求更加有效的方法。格子Boltzmann方法是一種新興的數值模擬方法，它基于統計物理學，具有獨特的粒子特性。其微觀動力學背景使得它具有許多其它基于N-S方程的數值方法所沒有的獨特優點[5]。格子Boltzmann方法的基本思路是：按照分子運動論的觀點，流場是運動粒子的宏觀效應；該方法將時間和空間完全離散，將流場劃分為格子；流體被抽象成大量的微觀粒子，并且這些微觀粒子根據某種簡單的運動規則在離散的格點上進行遷移和碰撞；粒子分布函數的演化在宏觀上反映流體的運動規律，流場的密度、速度等宏觀量可由粒子分布函數計算得到。

本文嘗試采用一種格子Boltzmann單相自由面模型[6]進行物體入水空泡流動的數值模擬研究，獲得入水空泡的變化規律。該模型的特點是忽略系統中氣相對液相的動力學影響，適用于具有大密度比的氣液兩相流動。同時這一嘗試將為在模擬更為復雜的自由面流動中采用格子Boltzmann方法打下基礎。

2 計算模型

2.1 格子Boltzmann模型

標準的格子BGK模型[7]如下式所示：

其中，i表示格子速度方向，fi（x,t）表示在t時刻x位置在i方向具有速度ci的粒子的分布函數，τ為無量綱松弛時間，為平衡態分布函數。

平衡態分布函數選擇如下：

宏觀量的速度、密度等可通過格點各方向的粒子分布函數的適當運算得到：

其中n為格點離散速度方向的個數。

本文采用D3Q19模型離散速度，其離散速度模型如圖1所示，格子聲速。 當i=0時，權重wi=1/3；當i=1～6時，權重wi=1/18；當i=7～8時，權重wi=1/36。

2.2 單相自由面模型

Thuery[6]提出的格子Boltzmann單相自由面模型適用于大密度比的氣液兩相流動，其特點是忽略了系統中氣相對液相的動力學影響。該模型中流場被劃分為三類格點：液相格點、界面格點、氣相格點。其中液相格點被液相完全充滿，氣相格點完全被氣相充滿，界面格點既包含液相也包含氣相。并且氣相格點和液相格點之間必須有界面格點存在。

每個格點都定義了流體體積分數ε＝m/ρ，其中m為格點的質量，ρ為格點密度。根據定義有：氣相格點ε=0，界面格點0<ε<1，液相格點ε＞1。界面格點質量的變化直接通過相鄰格點的粒子分布函數計算得到。對于x位置的界面格點，其相鄰格點位置為x+ciΔt，ci為離散速度。則該方向的質量流量可按下式計算：

其中下標iopp表示i的反方向。下一時間步x位置界面格點的質量為

由于氣相格點的物理量在計算中被忽略，從氣相格點遷移至液相格點的分布函數需要根據宏觀邊界條件進行重構：

其中ρA為氣體壓力對應的密度，下標iopp代表i的反向，u為x位置界面流體的速度。對于邊界處界面格點，本文采用文獻[8]中提出的使用重構的分布函數的方式處理，以獲得邊界處更為合理的自由面速度。

在界面運動的過程中，界面格點可能變為氣相格點或者液相格點。界面的重構則根據質量與密度的關系來決定。若碰撞步后界面格點的流體體積分數ε>1+κ（κ為選取的小量，如κ=10-3），則界面格點轉變為液相格點，它周圍的氣相格點則相應地變為界面格點；若ε<-κ，則界面格點轉變為氣相格點，它周圍的液相格點也需要相應地變為界面格點。為了保持質量守恒，格點類型的轉換完成后需要對多余的質量根據界面的法向進行重新分配。

2.3 重力的處理和邊界處理

由于重力對入水空泡的發展規律有重要的影響，本文引入重力的影響。重力的作用是一種額外的體積力，在格子Boltzmann方法中可以通過改變動量的方法引入重力的影響[9]。平衡態分布函數由新的速度=u+τg求得，即其中τ為松弛時間。

由于本文采用的LBM模型為弱可壓縮模型，因此在物體入水過程中會產生一系列的壓力波。這些壓力波在計算邊界處會發生反射，而產生的反射波會對流場的壓力分布產生影響，甚至影響自由面的形狀。為了減弱該影響，本文的入口邊界采用了基于插值的疊加格式[10]。四周的流場邊界則采用了充分發展邊界處理格式。圓盤處采用半步長反彈格式。

3 計算結果與分析

本文進行了不同Froude數下圓盤垂直等速入水的全三維數值模擬，并將計算結果與Bergmann等人[4]的實驗結果進行了比較。實驗中，入水圓盤的半徑為30mm，在三種工況中，平板分別以0.5m/s、1m/s和2m/s速度入水，則相應的Froude（定義為Fr=U2/gR）數分別為0.85，3.4和13.6。數值模擬計算示意圖如圖2所示，計算域的網格數為150×150×180。采用流體運動而圓盤靜止的方式進行圓盤垂直等速入水的模擬，計算中圓盤半徑R=15，初始水深H=40，重力加速度g=0.000 005，邊界入口速度U等于圓盤入水速度。

圖3－5分別給出了對應Fr=0.85、3.4、13.6等三種工況的入水空泡的數值模擬結果。

圖6為Bergmann等人給出的對應Fr=0.85、3.4、13.6等三種工況（分別對應圖6（a）、（b）、（c））的實驗照片[4]）。該實驗與現有的其它圓盤入水實驗的區別是其底部通過連桿嚴格控制了圓盤入水的速度。實驗照片中白色的線條系文獻[4]作者采用邊界元方法得到的計算結果，時間τ表示該時刻到空泡閉合所需的時間。

將圖3－5與圖6的相應試驗結果相比較，可見：在Fr=0.85工況下由于Froude數相對較小，入水空泡在離水面不遠處立即閉合，這與實驗過程相一致。隨著Froude數的增大，入水空泡的閉合深度逐漸增加。數值模擬結果也復現了這一現象。總體上看，數值模擬結果與試驗結果較為吻合。圖5的數值模擬結果中，液面附近未形成如實驗照片圖6（b）所示的冠狀結構，也未見液面的破碎、液滴的飛濺現象。原因可能是在本文所使用的單機計算性能的限制下，數值模擬中圓盤的半徑的格子密度不夠大，還不足以捕捉到液面附近精細的拓撲結果變化。

圖7給出了本文計算的入水空泡閉合時空泡相對閉合深度Zcoll/R與Froude數平方根的關系，并將計算結果與Bergmann等人[4]的實驗結果進行了比較，其中Zcoll為入水空泡閉合深度。

從圖7可見，本文數值模擬得到的空泡相對閉合深度Zcoll/R與Froude數平方根基本呈線性關系，這與實驗結果相一致；同時，其斜率也與實驗結果擬合得到的斜率相一致。計算得到的空泡相對閉合深度總體上比實驗值略為偏大，可能的原因是在模擬過程中，由于單機計算條件的限制，計算域的尺寸不夠大而產生的影響。

4 結論和展望

本文嘗試采用格子Boltzmann單相自由面模型數值模擬了圓盤垂直等速入水空泡的發展過程；計算結果給出了不同Froude數下圓盤入水空泡隨時間的變化過程；研究了入水空泡相對閉合深度與Froude數之間的關系；數值計算結果與相關實驗結果吻合較好。表明格子Boltzmann單相自由面模型可以用來研究出水空泡的發展過程問題，為其模擬更為復雜的自由面流動問題打下了基礎。

[1]陳九錫,顏 開.用MAC方法計算平頭物體垂直等速入水空泡[J].空氣動力學學報,1986,4(1):47-55.

[2]Kleefsman K M T,Fekken G,Veldman A E P,et al.A Volume-of-Fluid based simulation method for wave impact problems[J].Journal of Computational Physics,2005,206:363-393.

[3]龔 凱,劉 樺.圓盤垂直入水的SPH數值模擬[C].成都:第二十屆全國水動力學研討會文集,2007.

[4]Raymond Bergmann,Devaray van der Meer,Gekle,et al.Controlled impact of a disk on a water surface:Cavity dynamics[J].Journal of Fluid Mechanics,2009,633:381-409.

[5]Luo L S.The lattice-gas and lattice Boltzmann methods:past,present and future[C]//Proceeding of International Conference on Applied Computational Fluid Dynamics.Beijing,China,2000.

[6]Nils Thürey.Physically based animation of free surface flows with the Lattice Boltzmann Method[D].PHD Thesis,University of Erlangen-Nuremberg,2007.

[7]Qian Y,d’Humières D,Lallemand P.Lattice BGK models for Navier-Stokes equation[J].Europhysics Letters,1992,17(6),479-484.

[8]湯 波,李俊峰等.帶自由面流體運動的單相格子Boltzmann方法模擬[J].清華大學學報（自然科學版）,2008,48(11):34-360.

[9]Buick J M.Lattice Boltzmann methods in interfacial wave modeling[D].PHD Thesis,University of Edinburgh,U.K.,1997.

[10]Yu D Z,Luo L S,Wer S Y.Viscous flow computation with the method of lattice Boltzmann equation[J].Progress in Aerospace Sciences,2003,39:329-367.