劉永濤，朱仁慶

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

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

1 液體流動(dòng)的數(shù)學(xué)模型

考慮不可壓縮粘性流體假設(shè)，則液體運(yùn)動(dòng)的控制方程為:

連續(xù)方程

動(dòng)量輸運(yùn)方程

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

2 MPS計(jì)算模型

2.1 核函數(shù)選取

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

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

2.2 粒子密度計(jì)算

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

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

2.3 梯度算子模型

由于動(dòng)量輸運(yùn)方程中存在物理量的梯度項(xiàng)，因此定義梯度算子計(jì)算模型.考慮物理量φ，對(duì)于某相鄰的兩個(gè)粒子i和j，它們之間梯度矢量可表示為兩者物理量差值與對(duì)應(yīng)矢量半徑之比:

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

2.4 Laplace算子模型

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

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

在每一個(gè)時(shí)間步Δt，MPS法采用分步方法來求解動(dòng)量方程.

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

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

首先只考慮動(dòng)量方程中壓力梯度項(xiàng)，可得速度修正量 u′:

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

應(yīng)用迭代方法求解該Laplace方程，計(jì)算得到壓力場Pn+1.

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

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

2.6 自由表面判別

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

2.7 固體邊界條件

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

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

3 算例

3.1 潰壩模擬1

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

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

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

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

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

3.2 潰壩模擬2

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

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

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

4 結(jié)論

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

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］徐剛，段文洋.自由面流動(dòng)模擬的MPS算法研究［J］.哈爾濱工程大學(xué)學(xué)報(bào)，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］.水動(dòng)力學(xué)研究與進(jìn)展，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］.江蘇科技大學(xué)學(xué)報(bào):自然科學(xué)版，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.