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

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

1 引 言

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

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

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

2 計算模型

2.1 格子Boltzmann模型

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

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

平衡態(tài)分布函數(shù)選擇如下：

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

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

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

2.2 單相自由面模型

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

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

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

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

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

在界面運動的過程中，界面格點可能變?yōu)闅庀喔顸c或者液相格點。界面的重構(gòu)則根據(jù)質(zhì)量與密度的關(guān)系來決定。若碰撞步后界面格點的流體體積分數(shù)ε>1+κ（κ為選取的小量，如κ=10-3），則界面格點轉(zhuǎn)變?yōu)橐合喔顸c，它周圍的氣相格點則相應(yīng)地變?yōu)榻缑娓顸c；若ε<-κ，則界面格點轉(zhuǎn)變?yōu)闅庀喔顸c，它周圍的液相格點也需要相應(yīng)地變?yōu)榻缑娓顸c。為了保持質(zhì)量守恒，格點類型的轉(zhuǎn)換完成后需要對多余的質(zhì)量根據(jù)界面的法向進行重新分配。

2.3 重力的處理和邊界處理

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

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

3 計算結(jié)果與分析

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

圖3－5分別給出了對應(yīng)Fr=0.85、3.4、13.6等三種工況的入水空泡的數(shù)值模擬結(jié)果。

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

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

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

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

4 結(jié)論和展望

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

[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數(shù)值模擬[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.