移動(dòng)粒子半隱式法的大渦模擬研究與應(yīng)用

潘徐杰，張懷新

（上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200030）

1 引 言

移動(dòng)粒子半隱式法（MPS）由Koshizuka[1-2]最早于1995年提出的一種用于計(jì)算不可壓縮流體運(yùn)動(dòng)的無網(wǎng)格方法，是一種完全Lagrange的無網(wǎng)格數(shù)值方法。MPS通過各代表流體或邊界條件粒子的性質(zhì)來表現(xiàn)出流場(chǎng)的時(shí)－空物理特征。在MPS法中，控制方程采用Lagrange形式的Navier-Stokes方程，方程中不出現(xiàn)對(duì)流項(xiàng)，每個(gè)計(jì)算部分為顯式和隱式兩步，第一步考慮重力與邊界條件的影響，對(duì)粒子的速度和坐標(biāo)進(jìn)行修正，第二步由第一步引起的粒子數(shù)密度變化為條件，建立壓力Poisson方程，求得壓力后對(duì)粒子進(jìn)行第二次修正，使得粒子的數(shù)密度回到初始狀態(tài)，恰好滿足了不可壓縮條件。在MPS中，控制方程各項(xiàng)由粒子與周圍粒子的關(guān)系得到，方程通過這種關(guān)系進(jìn)行離散。同傳統(tǒng)的網(wǎng)格法比較，MPS由于控制方程中不存在對(duì)流項(xiàng)，避免了可能存在的數(shù)值耗散；因使用可移動(dòng)的粒子來代表流體，使得其在大變形自由表面的模擬中能夠較為滿意地表現(xiàn)出自由表面的形狀；并且粒子布置容易性使得該方法在各種復(fù)雜邊界條件中均能得心應(yīng)手。

由于MPS廣泛的適用性和良好的前景，該方法在被提出后[1-2]已被運(yùn)用到許多領(lǐng)域中。1999年，Gotoh[3]用MPS法來模擬波浪的傳播問題，討論了擋板在帶有固定不可穿透斜面的水槽中突然運(yùn)動(dòng)問題以及可穿透斜面的流體粒子運(yùn)動(dòng)過程問題、模擬了帶有垂直墻的破碎波問題；從2001年至2003年，Sueyoshi[4-6]在海洋工程領(lǐng)域使用MPS法做了一系列的研究。包括潰壩模型、晃蕩模型和具有斜坡長(zhǎng)底的潰壩模型的模擬[4]，不具有和具有消波板的潰壩模型的模擬[5]，以及具有液艙和不具有液艙的二維浮體在由數(shù)值水池模擬出的惡劣海況下的運(yùn)動(dòng)[6]，其中還包括了破損情況下的二維浮體運(yùn)動(dòng)等等；2004年，Hibi[7]則將重點(diǎn)放在了由MPS法自身引起的壓力振蕩現(xiàn)象，專門研究解決壓力振蕩問題，提出了四種不同的方法來緩解壓力問題。其中包括新的底部流體粒子壓力計(jì)算方式、異相異性的核函數(shù)、兩次的ICCG迭代、以及修改的壓力Poisson方程；2005年，Shao[8]使用帶有大渦湍流模式的MPS法和SPH法模擬同一個(gè)潰壩模型；另外，同年Gotoh[9]用MPS模擬了水躍問題；2007年，Alam[10]使用MPS模擬了楔形體入水問題，對(duì)于水面的破碎又分考慮表面張力和不考慮表面張力兩種情況；另外，2008年，Sueyoshi[14]使用MPS模擬了甲板上浪問題。

本文使用MPS法，結(jié)合Smagorinsky渦粘模型，模擬了半滿液箱的橫搖晃蕩運(yùn)動(dòng)。針對(duì)由MPS自身引起的壓力震蕩現(xiàn)象，選擇了帶權(quán)重的面積－時(shí)間平均法，將壓力震蕩現(xiàn)象緩解至可實(shí)際接受的范圍內(nèi)；在自由表面判別方法上，針對(duì)原方法容易造成非自由表面粒子被誤判這一缺陷，提出了一種基于鄰居搜索的混合判別方法，該方法在MPS法容易出錯(cuò)的階段使用鄰居搜索，能夠準(zhǔn)別地判別自由表面并減少誤判情況的產(chǎn)生。

2 控制方程

MPS使用Lagrangian形式的Navier-Stokes方程作為控制方程，過濾后的控制方程有如下形式：

方程（3）是Gradient模型，方程（4）是Laplacian模型，在上述兩方程中，d表示空間的維數(shù)，n0是初始粒子數(shù)密度，矢量xi和矢量xj是粒子i和粒子j的坐標(biāo)，方程（4）中的λ可由方程（5）求得，粒子數(shù)密度n0可由方程（6）求得：

在方程（3）、（4）、（5）和（6）中的w（r）為核函數(shù)，可由方程（7）求得：

在核函數(shù)中，re是核函數(shù)的作用距離，r是兩粒子間的距離。MPS中的壓力，可由壓力Poisson方程（8）求得：

將方程（2）中的亞格子應(yīng)力采用分子粘性形式，可寫為：

上式中Cs為Smagorinsky系數(shù)，Δ為MPS初始粒子間距。有了上述的方程，MPS法的計(jì)算過程可由圖1所示。

同一般MPS算法不同的是，考慮大渦模擬的MPS多了亞格子應(yīng)力需要考慮，亞格子應(yīng)力同體積力與邊界條件一起考慮，如下式所示：

上式中Δui*是中間速度值，Δt為時(shí)間增量。而（9）式中的項(xiàng)則并入壓力項(xiàng)，同壓力泊松方程一同求解，其他過程均同一般的MPS算法一致[13]。

3 基于鄰居搜索的混合自由表面判別

原MPS法中使用粒子數(shù)密度的變化來識(shí)別自由表面，如方程（11）所示：

式中＜n*＞i在某個(gè)時(shí)間步下的粒子i的數(shù)密度，MPS理論推薦β可以選取介于0.8到0.99之間的數(shù)。

在MPS模擬中，當(dāng)粒子處于自由表面處時(shí)，粒子周圍的鄰居較少，所以其數(shù)密度也較少，可以方便地識(shí)別為自由表面。在MPS的模擬過程中，自由表面的識(shí)別是在考慮了體積力和邊界條件的第一次粒子速度位置修正后進(jìn)行的，然后自由表面粒子壓力賦0后作為初始條件參與壓力泊松方程的求解。MPS理論要求在每一個(gè)模擬時(shí)間步之后，各粒子的數(shù)密度回歸為初始粒子數(shù)密度，相當(dāng)于（11）式中數(shù)密度參數(shù)β為1。但是自由表面的識(shí)別在第一次修正之后，粒子的位置發(fā)生了變化，也就是參數(shù)β值發(fā)生了變化，而MPS僅能保證在第二次修正后β值回復(fù)為1，所以在自由表面識(shí)別的過程中，將有非自由表面的粒子因β值略小而被誤判為自由表面粒子，從而壓力賦0參與泊松方程的求解，給模擬的結(jié)果帶來誤差，這種情況在流體運(yùn)動(dòng)比較劇烈時(shí)尤為明顯。

本文使用一種基于鄰居搜索的自由表面判定方法，鄰居搜索區(qū)域同Gradient模型方程（3）中的鄰居搜索距離一致，那么模擬時(shí)就不需要進(jìn)行額外的鄰居搜索工作?？紤]到并非只是最表層粒子屬于自由表面，實(shí)際模擬中靠近表層的粒子都可以作為自由表面處理，所以在鄰居搜索時(shí)，只考慮粒子周圍1.0倍粒子間距到2.1倍粒子間距之間的區(qū)域，將該區(qū)域分成八分，如有兩相鄰區(qū)域中沒有鄰居，則認(rèn)為該粒子屬于自由表面，反之則是一般粒子。

由于對(duì)每個(gè)粒子進(jìn)行鄰居搜索以判定自由表面較耗費(fèi)資源，而實(shí)際所需要的工作只是要消除因粒子數(shù)密度參數(shù)略小而被誤判的情況發(fā)生，所以模擬中進(jìn)行鄰居搜索的對(duì)象僅為按方程（11）所選出的粒子。另外，鑒于發(fā)生拍擊破碎的情況，可能會(huì)出現(xiàn)粒子不存在兩相鄰區(qū)域中沒有鄰居但確實(shí)屬于自由表面的情形，可將粒子的數(shù)密度要明顯小于一般粒子的直接歸并為自由表面。所以在使用鄰居搜索法時(shí)，只需要對(duì)部分粒子進(jìn)行識(shí)別，這些粒子數(shù)密度條件可滿足以下條件：

式中A為粒子數(shù)密度判定下限，由于這種自由表面判別方式結(jié)合傳統(tǒng)的MPS自由表面判別方式和鄰居搜索法，所以稱為混合自由表面判別法。

4 壓力震蕩現(xiàn)象與緩解

由于自身算法的原因，MPS中的壓力存在震蕩現(xiàn)象[11-13]。在這方面，Hibi[7]和Sueyoshi[5]都進(jìn)行過研究。Sueyoshi對(duì)浮體的非線性運(yùn)動(dòng)模擬中的壓力震蕩現(xiàn)象，使用了固定的面積—時(shí)間平均方法，緩解壓力震蕩并取得了很顯著的效果。而Hibi的多種方法雖然有一定的效果，但是卻在不同程度上復(fù)雜化了MPS算法，并且其效果同面積—時(shí)間平均方法相比，沒有絕對(duì)的提高，所以本文采用面積—時(shí)間平均方法。而在實(shí)際模擬過程中發(fā)現(xiàn)，數(shù)十倍于壓力正常值的壓力畸點(diǎn)較為頻繁地出現(xiàn)在靠近邊界處。這是因?yàn)樵诨问幍哪M中，邊界具有速度，在考慮了邊界條件的修正后，邊界粒子可能與流體粒子距離過近而導(dǎo)致其數(shù)密度過大，從而由方程（8）可看出該粒子的壓力值將過大。從而本文使用帶權(quán)重的壓力平均方式，當(dāng)輸出粒子點(diǎn)周圍有壓力畸點(diǎn)出現(xiàn)時(shí)，該平均方式可以一定程度上緩解受壓力畸點(diǎn)的影響。

5 MPS的大渦數(shù)值模擬

首先選取潰壩模型進(jìn)行模擬，模型的尺寸如圖3所示。在模擬中可以發(fā)現(xiàn)，基于鄰居搜索的自由表面判定方法能夠較為準(zhǔn)確地捕捉到自由表面粒子（如圖5所示），同原MPS的自由表面判定方法相比（如圖6所示），明顯地減少了自由表面誤判的產(chǎn)生，使用了混合自由表面判定條件的大渦MPS法，能夠準(zhǔn)確地模擬出潰壩的過程。

接著選取晃蕩模型進(jìn)行模擬，模型尺寸如圖4所示，液箱做強(qiáng)迫的橫搖正弦運(yùn)動(dòng)，其幅度為4°，橫搖的頻率為0.89Hz，橫搖的中心位于艙底中心。而模擬中的三個(gè)壓力輸出值分別放置于左壁的1/2高處、3/4高處與左壁頂部。從模擬結(jié)果可以看出，使用了混合自由表面判定條件后，模擬運(yùn)行順利，混合自由表面判定方法能夠準(zhǔn)確地捕捉自由表面粒子。由于在模擬中，除了使用混合自由表面判定條件外，還使用了區(qū)域限定。根據(jù)文獻(xiàn)[12]的研究結(jié)果，可以大約估算出自由表面將會(huì)出現(xiàn)的位置，對(duì)自由表面不會(huì)出現(xiàn)的位置劃為限定區(qū)域，對(duì)這部分區(qū)域里面的粒子不進(jìn)行自由表面判定而直接認(rèn)定為屬于非自由表面粒子，這樣做不但再次減小了非自由表面粒子被誤判的可能性，而且只對(duì)部分粒子進(jìn)行自由表面判定，還增加了模擬的效率，在本次模擬中，被限定的區(qū)域設(shè)置在液箱的半高以下。從模擬的結(jié)果上看，不多的自由表面誤判情況都集中在較靠近自由表面附近的區(qū)域，而在液箱中部和底部始終沒有誤判產(chǎn)生，如圖7所示，由于晃蕩模型中的液位較潰壩模型相比要深，可能會(huì)有偶爾的非自由表面粒子被誤判為自由表面的情況出現(xiàn)在液箱底部的情形，但使用了區(qū)域限制之后，則不會(huì)出現(xiàn)這種情況。在模擬的（1/4）T和（3/4）T時(shí)，液艙中的液體發(fā)生拍擊現(xiàn)象，與實(shí)驗(yàn)結(jié)論一致[13]。當(dāng)左壁頂部發(fā)生拍擊時(shí)，左壁各壓力點(diǎn)處的壓力均處于峰值，與拍擊基本處于同一時(shí)刻，這一點(diǎn)也與實(shí)驗(yàn)結(jié)論一致[13]。從圖7中可以看出，MPS在模擬自由表面的運(yùn)動(dòng)方面，有著出色的表現(xiàn)。模擬的壓力結(jié)果采用面積－時(shí)間平均法進(jìn)行處理，面積為壓力輸出點(diǎn)同臨近點(diǎn)按方程（13）進(jìn)行，其中n在本文取值為5，在時(shí)間平均中，時(shí)間段越長(zhǎng)，效果則越好[12]，但是過長(zhǎng)的時(shí)間可能將別的時(shí)段的壓力平均進(jìn)來并影響結(jié)果，所以本文選取時(shí)間段為0.012秒，約等于本文模擬周期1.1236秒的1%。從圖8可以看出，平均后的壓力無論在表現(xiàn)規(guī)律和極值上，都同實(shí)驗(yàn)值保持基本一致[13]。

6 結(jié) 論

本文使用結(jié)合Smagorinsky渦粘模型的大渦移動(dòng)粒子半隱式法模擬了水柱潰壩過程和液箱橫搖晃蕩過程。針對(duì)原MPS法在模擬中普遍存在的非自由表面粒子被誤判為自由表面情況，提出了一種混合自由表面判別法；針對(duì)MPS法的壓力震蕩現(xiàn)象，使用了面積－時(shí)間平均法進(jìn)行處理。從模擬的結(jié)果中可以看出，混合自由表面判定條件能夠很好地工作，在模擬中對(duì)被誤判的非自由表面粒子修正的效果明顯，基于面積－時(shí)間平均的壓力震蕩緩解方法效果明顯，可以將壓力震蕩現(xiàn)象緩解到工程實(shí)際的運(yùn)用范疇，另外上述方法同基于大渦模擬的MPS法配合良好，模擬順利，并且其在模擬自由表面的運(yùn)動(dòng)方面，與傳統(tǒng)的網(wǎng)格比更加方便，并且表現(xiàn)也更加出色。

[1]Koshizuka S,Tamako H,Oka Y.A particle method for incompressible viscous flow with fluid fragmentation[J].Journal of Computational Fluid Dynamics,1995,4:29-46.

[2]Koshizuka S,Oka Y.Moving-Particle Semi-implicit Method for fragmentation of incompressible fluid[J].Nuclear Science and Engineering,1996,123:421-434.

[3]Gotoh H,Sakai T.Largrangian simulation of breaking waves using particle method[J].Coast Engineering Journal,1999,41:303-326.

[4]Sueyoshi M,Naito S.A study of nonlinear fluid phenomena with particle method(Part1)[J].Journal of Kansai Society Naval Architects,2001,236:191-198.(in Japanese)

[5]Sueyoshi M,Naito S.A study of nonlinear fluid phenomena with particle method(Part2)[J].Journal of Kansai Society Naval Architects,2002,237:181-186.(in Japanese)

[6]Sueyoshi M,Naito S.A study of nonlinear fluid phenomena with particle method(Part3)[J].Journal of Kansai Society Naval Architects,2003,239:81-86.(in Japanese)

[7]Hibi S,Yabushita K.A study on reduction of unusual pressure fluctuation of MPS method[J].Journal of Kansai Society Naval Architects,2004,241:125-131.(in Japanese)

[8]Shao S D,Gotoh H.Turbulence particle models for tracking free surfaces[J].Journal of Hydraulic Research,2005,43:276-289.

[9]Gotoh H,et al.Lagrangian particle method for simulation of wave overtopping on a vertical seawall[J].Coastal Engineering Journal,2005,47:157-181.

[10]Alam A,Kai H,Suzuki K.Two-dimensional numerical simulation of water splash phenomena with and without surface tension[J].Journal of Marine Science and Technology,2007,12:59-71.

[11]Pan X J,Zhang H X,Lu Y T.Moving-Particle Semi-implicit Method for vortex patterns and roll damping of 2D ship sections[J].China Ocean Engineering,2008,22:399-407.

[12]潘徐杰,張懷新.移動(dòng)粒子半隱式法晃蕩模擬中的壓力震蕩現(xiàn)象研究[J].水動(dòng)力學(xué)研究與進(jìn)展,2008,23(4):453-463.Pan X J,Zhang H X.A study on pressure fluctuation in sloshing simulation of Moving-Particle Semi-implicit Method[J].Journal of Hydrodynamic,2008,23(4):453-463.(in Chinese)

[13]潘徐杰,張懷新.用移動(dòng)粒子半隱式法模擬液艙橫搖晃蕩現(xiàn)象[J].上海交通大學(xué)學(xué)報(bào),2008,42(11):1904-1907.Pan X J,Zhang H X.Moving-Partical Semi-implicit Method for simulation of liquid sloshing on roll motion[J].Journal of Shanghai Jiaotong university,2008,42(11):1904-1907.(in Chinese)

[14]Sueyoshi M,Kashiwagi M,Naito S.Numerical simulation of wave-induced nonlinear motions of a two-dimensional floating body by the Moving Particle Semi-implicit Method[J].Journal of Marine Science and Technology,2008,13:85-94.