劉永濤，朱仁慶

(江蘇科技大學船舶與海洋工程學院，江蘇鎮(zhèn)江212003)

自從文獻［1］于1996年首先提出MPS(moving particle semi-implicit，MPS)方法以來，作為一種擅長于處理自由表面大變形問題的無網(wǎng)格粒子法，該方法獲得了快速發(fā)展，并在眾多流動問題模擬中獲得了廣泛應用，例如潰壩［2-5］、波面破碎［6-7］、波浪中浮體運動［8］、甲板上浪［9-10］、液艙晃蕩［10-12］等問題.文中針對潰壩問題，采用MPS方法進行二維自由表面大變形現(xiàn)象的數(shù)值模擬與驗證.

1 液體流動的數(shù)學模型

考慮不可壓縮粘性流體假設，則液體運動的控制方程為:

連續(xù)方程

動量輸運方程

式中:μ為流體動力粘性系數(shù);Fb為質量力;u為流速;ρ為流體密度;p為壓強.

2 MPS計算模型

2.1 核函數(shù)選取

采用MPS法對計算流場進行粒子離散化后，流體粒子就擁有位置、速度、質量和壓力這些物理量，而且在模擬過程中，流體粒子將和周圍粒子產(chǎn)生相互作用并交換物理量.通常將流體粒子間的這種相互作用限定于影響半徑范圍內(nèi)，并通過定義核函數(shù)w(r)來描述該范圍內(nèi)粒子相互作用的變化規(guī)律.文中所采用的核函數(shù)為:

式中:re為影響半徑，文中取re=3.2l0，l0為初始相鄰粒子間距.

2.2 粒子密度計算

MPS法模擬過程中，隨著流體粒子的不斷移動將可能引起流體密度的變化.這種流體密度的變化通過定義某一粒子i影響域內(nèi)相鄰粒子的密集程度，即粒子數(shù)密度ni來實現(xiàn).因此，在ri處的粒子i的粒子數(shù)密度為:

式中:rj為周圍粒子j的瞬時位置.因此對不可壓流體，這意味著計算過程中流體的粒子數(shù)密度ni需要保持為初始粒子密度常數(shù)n0.

2.3 梯度算子模型

由于動量輸運方程中存在物理量的梯度項，因此定義梯度算子計算模型.考慮物理量φ，對于某相鄰的兩個粒子i和j，它們之間梯度矢量可表示為兩者物理量差值與對應矢量半徑之比:

則在粒子i的影響域內(nèi)物理量φ的梯度即為周圍物理量梯度的加權平均:式中:d為流場的空間維數(shù).

2.4 Laplace算子模型

對于動量輸運方程中物理量的Laplace算子項，采用文獻［3］所給出的公式:式中系數(shù)λ為:

2.5 控制方程的數(shù)值解法

在每一個時間步Δt，MPS法采用分步方法來求解動量方程.

1)考慮動量方程中質量力項以及粘性力項，顯式計算粒子i的速度變化量Δu*i以及臨時速度和臨時位置為:2)計算該瞬時位于臨時位置r*i的粒子數(shù)密度n*.由于在一般情況下，n*不等于初始粒子數(shù)密度n0，為滿足流體不可壓條件，因此將其修正到n0.該修正過程通過定義粒子數(shù)密度n*的修正值n′來實現(xiàn).

3)求解壓力Poisson方程得到壓力Pn+1.該方程由考慮動量方程壓力梯度項得到的速度修正量u′和粒子數(shù)密度修正值n′之間關系建立.

首先只考慮動量方程中壓力梯度項，可得速度修正量 u′:

再根據(jù)連續(xù)方程，粒子數(shù)密度修正值n′與速度修正值u′需要滿足如下關系式:

應用迭代方法求解該Laplace方程，計算得到壓力場Pn+1.

4)將新壓力場Pn+1代入式(13)計算u′，再次修正得到新時刻速度場uni+1和粒子位置rni+1.

由上述兩式可得到壓力Poisson方程:

2.6 自由表面判別

根據(jù)粒子數(shù)密度定義可知，自由表面處的流體粒子數(shù)密度會小于內(nèi)部流體粒子數(shù)密度.因此，前述MPS法分步計算過程第2)步中，若根據(jù)式(12)得到的臨時粒子，數(shù)密度n*i＜βn0，則判別其為自由表面粒子并設定壓力為零，此處根據(jù)文獻［13］取 β 為0.97.

2.7 固體邊界條件

在固體邊界上設置一層邊界粒子來防止流體粒子穿越固體邊界.由于靠近邊界處粒子數(shù)密度將會小于內(nèi)部流體粒子數(shù)密度，并被誤判為自由表面粒子，為此，在固體邊界外部另外平行布置若干層位置固定的虛擬粒子［2］，文中取為3層，具體布置如圖1，該類粒子僅參與局部粒子數(shù)密度計算，而不參與壓強和流速計算.

圖1 邊界處粒子布置示意圖Fig.1 Schematic view of the layout of boundary particles

3 算例

3.1 潰壩模擬1

流動域尺寸:長 L=0.8 m，高度 H=0.6 m，初始靜止水柱長度 a=0.175 m，高度h=0.35 m，被擋板限制于流域左側.當擋板打開后，水柱將在重力作用下向右側流動，并在不同時刻流經(jīng)底板的不同水平位置處，到達右壁面后在慣性作用下繼續(xù)沿壁面上沖.為研究粒子直徑對計算結果的影響，取粒子直徑 l0分別為0.004，0.006，0.008 m.

圖2 自由表面的試驗結果［14］與數(shù)值模擬結果對比Fig.2 Comparison of free surface configurations between experimental results and numerical results

圖3 波面流經(jīng)底板水平位置的模擬結果對比Fig.3 Comparison of numerical results for the water front positions along the tank bottom

圖4 波面流經(jīng)底板水平位置的試驗以及模擬結果(l0=0.004 m)對比Fig.4 Experimental and numerical results(l0=0.004 m)for the water front positions along the tank bottom

對上述潰壩過程的考察，主要關注各瞬時的自由表面大變形特征，相關對應結果可用圖2～4表示.圖2為瞬時自由表面結果，顯示了在不同時刻自由表面的數(shù)值模擬結果(右圖，對比3種粒子直徑流動模擬結果可知較小粒子直徑對復雜局部流動能有較精細地刻畫，本處僅提供粒子直徑等于0.004 m 的計算結果)與試驗結果(左圖，Hu［14］的試驗)的對比情況.從4個典型時刻的對比圖可以看出，自由表面的數(shù)值模擬結果與試驗結果較為一致.圖3為3種不同粒子直徑條件下不同時刻底板處波面水平位置的數(shù)值模擬對比，可見隨著粒子直徑減小，數(shù)值模擬結果趨于穩(wěn)定，文獻［4］也有類似結論.因此，選取較小的流體粒子直徑會有較好的計算結果，但壓力Poisson方程的求解也更耗時.綜合權衡計算耗時以及數(shù)值模擬精度，文中選取相對較小的粒子直徑.圖4為波面時歷結果，顯示了在不同時刻運動波面流經(jīng)底板水平位置時的對比，可見 MPS 模擬結果(l0=0.004 m)與 Hu［14］的試驗結果在無量綱時刻小于2.0時吻合較好，在這之后兩者略有分離，并且MPS模擬結果略大于試驗結果.

3.2 潰壩模擬2

流動域尺寸:長 L=0.584 m，高度 H=0.438 m;初始靜止水柱長度a=0.146 m，高度h=0.292 m，位于流動域左側.在距離左壁面0.292 m的底板處固定一方形木塊，其長度為0.024 m，高度為0.048 m.當擋板打開后，水柱向右側流動并經(jīng)過木塊位置處，波面破碎后繼續(xù)沖向右壁面.為較精確細致地反映波面流經(jīng)木塊后的破碎變形狀況，本算例選取粒子直徑為較小值，即l0=0.002 m，流動域共計有流體粒子10 658個，邊界粒子1 022個，邊界外三層虛擬粒子3 102個.

圖5 自由表面的試驗結果［15-16］與數(shù)值模擬結果對比Fig.5 Comparison of free surface configurations betweenexperimental results［15 -16］ and numerical results

在本算例中，由于水柱流經(jīng)方木塊受阻并產(chǎn)生劇烈沖擊后，沿木塊頂部快速躍起形成舌狀，此后部分流體在慣性作用下繼續(xù)向右沖擊右壁面，其余部分受重力作用下落并沖擊底板.該過程中自由表面產(chǎn)生顯著的破碎、翻卷、合并、沖頂?shù)鹊湫偷拇笞冃芜\動(圖5).對上述過程的真實有效模擬是重點也是難點.在圖5中，從0.1～0.5s4個不同時刻，對比自由表面的數(shù)值模擬結果(右圖)與試驗結果［15-16］(左圖)可以看出，0.1～0.3s各瞬時大變形自由表面的MPS數(shù)值模擬效果較為真實，僅在0.5s時刻兩者有明顯差別，表現(xiàn)在試驗結果中流動域右下側出現(xiàn)大體積氣體卷入現(xiàn)象，而模擬結果流體下落并堆積于右下側，同時由于粒子法本身原因出現(xiàn)大量散落的流體粒子.

4 結論

針對兩種潰壩問題，文中應用MPS法對自由表面大變形現(xiàn)象進行了數(shù)值模擬，對比分析了不同粒子直徑大小對數(shù)值模擬結果的影響，并與試驗結果進行了對比驗證.結果表明:選取較小粒子直徑有助于較精細地模擬大變形局部流動現(xiàn)象，且提高模擬計算結果的穩(wěn)定性;MPS方法對翻卷、破碎、沖頂?shù)茸杂杀砻娲笞冃螁栴}具有優(yōu)勢，模擬效果較好.

References)

［1］Koshizuka S，Oka Y.Moving-particle semi-implicit method for fragmentation of incompressible fluid［J］.Nuclear Science and Engineering，1996，123(3):421 -434.

［2］Koshizuka S，Tamako H，Oka Y.A particle method for incompressible viscous flow with fluid fragmentation［J］.Computational Fluid Dynamics Journal，1995，4(1):29 -46.

［3］Shakibaeinia A，Jin Y C.A weakly compressible MPS method for modeling of open-boundary free-surface flow［J］.International Journal for Numerical Methods in Fluids，2010，63(10):1208 -1232.

［4］徐剛，段文洋.自由面流動模擬的MPS算法研究［J］.哈爾濱工程大學學報，2008，29(6):539-543.Xu Gang，Duan Wenyang.An investigation of MPS algorithm for free surface flow simulation［J］.Journal of Harbin Engineering University，2008，29(6):539 -543.(in Chinese)

［5］張雨新，萬德成.用MPS方法數(shù)值模擬低充水液艙的晃蕩［J］.水動力學研究與進展，2012，27(1):100-107.Zhang Yuxin，Wan Decheng.Numerical simulation of liq-uid sloshing in low-filling tank by MPS［J］.Chinese Journal of Hydrodynamics，2012，27(1):100 -107.

［6］ Koshizuka S，Nobe A，Oka Y.Numerical analysis of breaking waves using the moving particle semi-implicit method［J］.International Journal for Numerical Methods in Fluids，1998，26(7):751 -769.

［7］Gotoh H，Sakai T.Key issues in the particle method for computation of wave breaking［J］.Coastal Engineering，2006，53(2):171 -179.

［8］Shibata K，Koshizuka S，Sakai M，et al.Lagrangian simulations of ship-wave interactions in rough seas［J］.O-cean Engineering，2012，42:13 -25.

［9］Shibata K，Koshizuka S.Numerical analysis of shipping water impact on a deck using a particle method［J］.O-cean Engineering，2007，34(9):585 -593.

［10］Shibata K，Koshizuka S，Tanizawa K.Three-dimensional numerical analysis of shipping water onto a moving ship using a particle method［J］.Journal of Marine Science and Technology，2009，14(2):214 -227.

［11］Pan Xujie，Zhang Huaixin，Lu Yuntao.Numerical simulation of viscous liquid sloshing by moving-particle semiimplicit method［J］.Journal of Marine Science and Application，2008，7(3):184 -189.

［12］龔少軍，姚震球.基于粒子法的液艙共振晃蕩現(xiàn)象研究［J］.江蘇科技大學學報:自然科學版，2010，24(6):534-538.Gong Shaojun，Yao Zhenqiu.Study on resonance sloshing by particle method［J］.Journal of Jiangsu University of Science and Technology:Natural Science Edition，2010，24(6):534 -538.(in Chinese)

［13］Khayyer A，Gotoh H.Enhancement of stability and accuracy of the moving particle semi-implicit method［J］.Journal of Computational Physics，2011，230(8):3093 -3118.

［14］Hu Changhong，Sueyoshi M.Numerical simulation and experiment on dam break problem［J］.Journal of Marine Science and Application，2010，9(2):109 -114.

［15］Martin J C，Moyce W J.An experimental study of the collapse of liquid columns on a rigid horizontal plane［J］.R Soc，1952，244(882):312 -324.

［16］Koshizuka S，Oka Y，Tamako H.A particle method for calculating splashing of incompressible viscous fluid［C］∥In:Proceedings of the International Conference，Mathematics and Computations，Reactor Physics，and Environmental Analyses.Washington，USA:American Nuclear Society，1995:1514-1521.