粒子法中壓力振蕩問(wèn)題的研究

朱 躍， 姜?jiǎng)僖?楊星團(tuán)， 段日強(qiáng)

（清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，教育部先進(jìn)反應(yīng)堆工程與安全重點(diǎn)實(shí)驗(yàn)室，北京100084）

1 引 言

移動(dòng)粒子半隱式法 MPS（Moving Particle Semi－implicit Method）是一種基于 Lagrange粒子的無(wú)網(wǎng)格方法，最早由Koshizuka［1］提出，隨后Koshizuka等［2，3］將 MPS方法應(yīng)用于計(jì)算核能領(lǐng)域的一些熱工水力難題，如蒸汽爆炸、兩相流和液滴破裂等大變形問(wèn)題。MPS方法發(fā)展至今，一直存在壓力振蕩問(wèn)題，國(guó)內(nèi)外學(xué)者也對(duì)此問(wèn)題進(jìn)行了研究。Sueyoshi等［4］為了解決 MPS方法中的壓力振蕩問(wèn)題，對(duì)獲得的壓力進(jìn)行后處理，雖然可以得到光滑的壓力曲線，但壓力振蕩仍然存在。Hibi等［5］提出了固體壓力邊界條件和異相異性的核函數(shù)，對(duì)壓力泊松方程進(jìn)行兩次ICCG迭代求解，以修改壓力泊松方程的方式對(duì)MPS方法中的壓力振蕩進(jìn)行抑制。Khayyer等［6］采用一個(gè)新型的壓力梯度計(jì)算公式，修改壓力泊松方程源項(xiàng)，并且對(duì)壓力泊松方程的系數(shù)矩陣項(xiàng)引入弱可壓縮條件，使得MPS方法在模擬過(guò)程中的壓力振蕩得到抑制。Khayyer等［7］提出了一種拉普拉斯算子的高階精度格式，用于離散壓力泊松方程和粘性力項(xiàng)。Kondo等［8］提出了一種新型的壓力泊松方程源項(xiàng)，該源項(xiàng)包含一個(gè)主要項(xiàng)和兩個(gè)誤差補(bǔ)償項(xiàng)三個(gè)部分，并且對(duì)靜水問(wèn)題和潰壩問(wèn)題進(jìn)行了模擬，結(jié)果表明，壓力振蕩問(wèn)題得到了抑制。Lee等［9］通過(guò)改進(jìn)壓力泊松方程的源項(xiàng)、梯度模型、粒子碰撞模型以及自由表面的識(shí)別等因素，逐漸改善壓力振蕩問(wèn)題。Tamai等［10］提出了最小二乘法移動(dòng)粒子半隱式方法 LSMPS（Least squares moving particle semiimplicit method）。Shibata等［11］加入了虛擬粒子的概念，并且在此基礎(chǔ)上修改了壓力泊松方程的系數(shù)矩陣項(xiàng)和源項(xiàng)，改進(jìn)了自由表面粒子的識(shí)別，修改了壓力梯度模型，改善了MPS方法的穩(wěn)定性。

本文針對(duì)傳統(tǒng)MPS方法中存在的壓力振蕩問(wèn)題，采用一種新型的壓力泊松方程離散格式來(lái)抑制MPS方法的壓力振蕩，新型的壓力泊松方程格式包括對(duì)粒子稀疏程度不敏感的壓力系數(shù)矩陣離散格式以及能夠有效抑制壓力振蕩的混合源項(xiàng)。在核函數(shù)的選取上，選擇了二次樣條核函數(shù)，Johnson等［12］首次運(yùn)用該核函數(shù)模擬高速?zèng)_擊問(wèn)題，有效地改善了計(jì)算穩(wěn)定性。針對(duì)MPS方法在模擬過(guò)程中的粒子插值不完整問(wèn)題，提出了粒子插值不完整性的修正。最終有效地抑制了在模擬過(guò)程中的壓力振蕩。

2 數(shù)值計(jì)算方法

2．1 控制方程

對(duì)于不可壓縮流體，其連續(xù)性方程和Navier－Stokes方程可表示為

式中u為流體的速度，P為流體的壓力，t為時(shí)間，ρ為流體的密度，μ為流體的動(dòng)力粘度，G為重力加速度。

2．2 傳統(tǒng)MPS方法

2．2．1 核函數(shù)

在MPS方法中，粒子間的相互作用關(guān)系最基本的體現(xiàn)就是通過(guò)核函數(shù)來(lái)實(shí)現(xiàn)的。在MPS方法中，核函數(shù)的公式為

式中r＝｜ri－rj｜為兩個(gè)粒子之間的距離，ri為i粒子的位置，rj為i粒子支持域內(nèi)j粒子的位置，re為粒子作用域的半徑，一般取re＝2．1l0，其中l(wèi)0為初始粒子間距。

2．2．2 梯度模型

梯度模型是MPS方法中的另一個(gè)基本模型。假定∮為任意標(biāo)量，那么表示粒子之間相互作用的梯度模型可表示為

式中d為空間維度，n0為初始粒子數(shù)密度，wij＝w（｜rj－ri｜）。在 MPS法的模擬中，只有作為標(biāo)量的壓力值P參與式（4）的計(jì)算。為了計(jì)算的穩(wěn)定性，MPS方法中的壓力梯度計(jì)算一般表示為

式中 Pj為i粒子支持域內(nèi)粒子的壓力，Pimin為i粒子支持域內(nèi)所有粒子中最小的粒子壓力。

2．2．3 Laplace模型

采用Laplace微分算子度量物理量的空間擴(kuò)散。假定∮為任意標(biāo)量，那么表示粒子之間相互作用的Laplace模型可表示為

式中d為空間維度，λ定義為

引入λ是對(duì)有限范圍內(nèi)的核函數(shù)代替無(wú)限范圍的高斯函數(shù)帶來(lái)誤差的一種補(bǔ)償。Laplace模型用來(lái)離散二階導(dǎo)數(shù)項(xiàng)，在MPS法的模擬中，粘性力項(xiàng)用Laplace微分算子進(jìn)行離散。

2．2．4 壓力泊松方程

在傳統(tǒng)的MPS方法中，壓力是通過(guò)求解壓力泊松方程得到的。其中，系數(shù)矩陣項(xiàng)的離散通過(guò)Laplacian模型進(jìn)行離散，源項(xiàng)采用粒子數(shù)密度的偏差，可表示為

2．3 改進(jìn)的MPS方法

2．3．1 壓力泊松方程

在改進(jìn)的MPS方法中，本文借鑒了SPH方法的思想，SPH 方法是由 Lucy等［13，14］獨(dú)立提出來(lái)的一種粒子法，SPH方法廣泛應(yīng)用于流體動(dòng)力學(xué)問(wèn)題。SPH方法和MPS方法類(lèi)似，都屬于拉格朗日粒子法，但兩者也有著明顯的區(qū)別。在SPH方法中，流體壓力通過(guò)狀態(tài)方程計(jì)算得到，計(jì)算效率較高，在MPS方法中，流體壓力通過(guò)求解壓力泊松方程得到，計(jì)算量較大；SPH方法有著較為嚴(yán)格的數(shù)學(xué)理論基礎(chǔ)，可以變化得到不同形式的離散格式，而MPS方法則是通過(guò)粒子間的加權(quán)平均計(jì)算得到，沒(méi)有類(lèi)似于SPH方法的嚴(yán)格的數(shù)學(xué)推導(dǎo)。

在改進(jìn)的MPS方法中，對(duì)于系數(shù)矩陣項(xiàng)的離散，采用 Monaghan［15］在SPH方法中離散熱傳導(dǎo)的格式，此離散格式對(duì)于粒子分布的不均勻性不敏感。此離散格式也同時(shí)用于離散壓力泊松方程［16］，并且有效地解決了壓力解耦問(wèn)題。在源項(xiàng)的選擇方面，本文選擇能夠較好解決壓力振蕩的混合源項(xiàng)，其主要由速度散度項(xiàng)和粒子數(shù)密度偏差項(xiàng)兩部分構(gòu)成。所以得到的壓力泊松方程為

式中γ為一個(gè)參數(shù)，取γ＝0．3。γ的選擇具有經(jīng)驗(yàn)性，并沒(méi)有固定的選擇依據(jù)。

2．3．2 核函數(shù)

在核函數(shù)的選擇上，本文采用的是二次樣條核函數(shù)，其公式和一階導(dǎo)數(shù)為

式中s＝rij／h，rij為兩個(gè)粒子之間的距離，h為光滑長(zhǎng)度。在MPS方法中，為了使二次樣條核函數(shù)適應(yīng)MPS方法，當(dāng)計(jì)算粒子數(shù)密度和梯度時(shí)，光滑長(zhǎng)度h＝1．05l0，當(dāng)計(jì)算Laplacian模型時(shí)，光滑長(zhǎng)度h＝2．0l0；αd的值在一維、二維和三維空間中分別為1／h，2／（πh2）和5／（4πh3）。

2．3．3 粒子插值不完整性的修正

在粒子法中，粒子積分值的近似等于加權(quán)累加該粒子支持域內(nèi)所有鄰居粒子的值。當(dāng)中心粒子支持域內(nèi)的粒子分布均勻時(shí)，對(duì)于該粒子積分值的插值近似是完整可靠的。粒子2支持域內(nèi)的粒子分布均勻如圖1所示，其粒子積分值的粒子插值近似是準(zhǔn)確的。但是在MPS方法模擬的過(guò)程中，粒子的分布會(huì)變得不均勻，這將使得粒子插值近似產(chǎn)生誤差。模擬過(guò)程中，粒子3支持域內(nèi)的粒子分布不均勻，所以其粒子積分值的粒子插值近似將會(huì)產(chǎn)生誤差。在自由表面邊界處的粒子，粒子支持域由自由表面截?cái)?，這將導(dǎo)致支持域內(nèi)的鄰居粒子小于流體內(nèi)部粒子，造成很大的誤差，如圖1中粒子1。

為了解決上述問(wèn)題，本文提出了一個(gè)解決粒子插值不完整性的方案?；舅枷霝?，每個(gè)粒子在初始時(shí)刻，其鄰居粒子是均勻分布的，則在其支持域內(nèi)的粒子數(shù)量是一定的。在模擬的過(guò)程中，如果支持域內(nèi)的粒子數(shù)量等于初始時(shí)刻的粒子數(shù)量，則認(rèn)為粒子分布保持均勻；否則進(jìn)行修正，使其數(shù)量等于初始時(shí)刻的粒子數(shù)量。

圖1 不同流體粒子支持域內(nèi)粒子分布示意Fig．1 Particle distribution in the support domain of different particles

由于MPS方法中，其支持域半徑的選取有兩種情況，在粒子數(shù)密度模型和梯度模型中，支持域半徑re＝2．1l0；在拉普拉斯模型中，支持域半徑re＝4．0l0。對(duì)于re＝2．1l0，粒子在初始時(shí)刻其支持域內(nèi)的鄰居粒子數(shù)量為12個(gè)，如圖2所示。并且這12個(gè)粒子均勻地分布在邊長(zhǎng)為l0的12個(gè)網(wǎng)格內(nèi)。然后進(jìn)行如下判斷。

（1）判斷粒子在支持域內(nèi)的鄰居粒子是否是12個(gè)，是則認(rèn)為該粒子插值是完整的。

（2）如果鄰居粒子小于12個(gè)，通過(guò)依次檢索每個(gè)網(wǎng)格中的粒子將其進(jìn)行補(bǔ)充。如圖3所示，在網(wǎng)格9內(nèi)沒(méi)有鄰居粒子，則通過(guò)其關(guān)于粒子i中心對(duì)稱(chēng)的鄰居粒子4來(lái)映射補(bǔ)充，其位置計(jì)算為xj9＝2xi－xj4，yj9＝2yi－yj4。壓力的計(jì)算分以下兩種情況。

（i）當(dāng)粒子9位于自由表面外時(shí)，壓力值pj9＝0。

（ii）當(dāng)粒子9在流體內(nèi)部時(shí)，壓力值

式中 下標(biāo)j9k為新生成的粒子9支持域內(nèi)的鄰居粒子。

圖2 初始時(shí)刻re＝2．1l0時(shí)支持域內(nèi)粒子分布情況Fig．2 Particle distribution in the support domain ofre＝2．1l0condition at initial instant

圖3 粒子的補(bǔ)充Fig．3 Supplement of particles

（3）最后檢查鄰居粒子數(shù)是否超過(guò)12個(gè)，如果超過(guò)12個(gè)，則將多余的粒子進(jìn)行刪除。刪除原則如下，如圖4所示，在網(wǎng)格2中出現(xiàn)了兩個(gè)粒子j1和j2，對(duì)粒子進(jìn)行粒子數(shù)密度j1和j2計(jì)算，并且選擇刪除偏離初始粒子數(shù)密度值n0較大的粒子，偏離值Δn計(jì)算如下，re＝4．0l0時(shí)，粒子在其支持域內(nèi)的鄰居粒子數(shù)量為48個(gè)，其計(jì)算過(guò)程與上述情況類(lèi)似。需要說(shuō)明的是，上述生成的補(bǔ)充粒子并不是在實(shí)際的模擬過(guò)程中生成新粒子，而是在梯度模型、粒子數(shù)密度模型和Laplacian模型的插值計(jì)算過(guò)程中加入補(bǔ)充粒子的信息，以達(dá)到在MPS方法模擬過(guò)程中粒子均勻分布的效果。

3 數(shù)值模擬

3．1 潰壩模型的計(jì)算

潰壩模型是一個(gè)經(jīng)典的研究自由表面大變形問(wèn)題的模型，本文選其作為研究MPS方法中壓力振蕩的第一個(gè)算例。Hu等［17］對(duì)潰壩問(wèn)題進(jìn)行了實(shí)驗(yàn)研究，并且測(cè)得潰壩過(guò)程中的壓力－時(shí)間曲線。為了與文獻(xiàn)［17］的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較，本文采用的潰壩模型尺寸與文獻(xiàn)［17］相同，其潰壩模型如圖5所示，A點(diǎn)為壓力監(jiān)測(cè)點(diǎn)。潰壩內(nèi)流體介質(zhì)為水，計(jì)算的粒子數(shù)量為7592個(gè)。設(shè)計(jì)了三種計(jì)算方案列入表1。

3．2 潰壩模型的計(jì)算結(jié)果及分析

圖4 粒子的刪除Fig．4 Deletion of particles

圖5 潰壩計(jì)算模型Fig．5 Sketch of the dam－break problem

表1 在模擬過(guò)程中的不同計(jì)算方案設(shè)置Tab．1 Condition of equations and schemes used in the simulation

圖6為A點(diǎn)在不同計(jì)算方案下的壓力－時(shí)間對(duì)比。從圖6（a）可以看出，傳統(tǒng)MPS方法模擬潰壩問(wèn)題時(shí)，在監(jiān)測(cè)點(diǎn)A處出現(xiàn)了劇烈振蕩。而采用改進(jìn)方案1時(shí)，相對(duì)于傳統(tǒng)的MPS方法，其壓力振蕩幅度減小了很多，并且其壓力振蕩的頻率也較低。在t＝0．35s～0．42s和t＝0．7s～0．95s時(shí)間內(nèi)，采用傳統(tǒng)MPS方法出現(xiàn)了明顯大于改進(jìn)方案1的壓力振蕩幅度，這兩段時(shí)間內(nèi)非物理壓力波動(dòng)的原因不同。在t＝0．35s～0．42s時(shí)間內(nèi)，潰壩流體沖擊右側(cè)豎直壁面，所以A點(diǎn)的壓力振蕩幅度非常大；而在t＝0．7s～0．95s這段時(shí)間內(nèi)，是波浪翻卷后又一次沖擊水面造成A點(diǎn)出現(xiàn)劇烈的壓力振蕩。因?yàn)閭鹘y(tǒng)的MPS方法中源項(xiàng)采用了粒子數(shù)密度變化率，在劇烈沖擊過(guò)程中，粒子的疏密程度將發(fā)生劇烈的變化，所以粒子數(shù)密度變化也將非常敏感，直接導(dǎo)致了傳統(tǒng)MPS方法源項(xiàng)值變化十分劇烈，從而其壓力振蕩的幅度會(huì)比較劇烈。而改進(jìn)方案1的源項(xiàng)采用了混合速度散度和粒子數(shù)密度變化率的混合源項(xiàng)?；旌显错?xiàng)中粒子的速度散度并不會(huì)因?yàn)榱Ｗ拥氖杳艹潭茸兓鴦×易兓?，所以其壓力振蕩幅度比較小。所以相對(duì)于傳統(tǒng)MPS方法，采用新的壓力泊松方程形式，能夠很好地抑制壓力振動(dòng)的幅度和頻率。

圖6 A點(diǎn)處不同MPS計(jì)算方案以及文獻(xiàn)［17］實(shí)驗(yàn)數(shù)據(jù)的壓力－時(shí)間對(duì)比Fig．6 Comparison of time histories of pressure at Apoint by different schemes and experimental data by Ref．［17］

從圖6（b）可以看出，采用粒子插值不完整性修正后，A點(diǎn)處的壓力振蕩幅度進(jìn)一步降低，并且壓力值與文獻(xiàn)［17］的實(shí)驗(yàn)值非常接近。實(shí)際上文獻(xiàn)［17］的實(shí)驗(yàn)值是經(jīng)過(guò)平均化處理后的結(jié)果，改進(jìn)方案2與文獻(xiàn)［17］的實(shí)驗(yàn)結(jié)果非常吻合。在t＝0．35s時(shí)，改進(jìn)方案2計(jì)算得到的壓力峰值與文獻(xiàn)［17］的實(shí)驗(yàn)值接近，在t＝0．75s時(shí)，出現(xiàn)了與文獻(xiàn)［17］實(shí)驗(yàn)值類(lèi)似的壓力峰值，這是由于波浪翻卷后又一次沖擊水面導(dǎo)致A點(diǎn)處產(chǎn)生的壓力峰值。圖7所示為A點(diǎn)處兩種計(jì)算方案下源項(xiàng)值和粒子數(shù)密度變化率隨時(shí)間的變化。從圖7（b）可以看出，相對(duì)于沒(méi)有采用粒子插值不完整修正的改進(jìn)方案1，采用粒子插值不完整性修正后的改進(jìn)方案2在A點(diǎn)處的粒子數(shù)密度變化率明顯降低，并且沒(méi)有類(lèi)似于改進(jìn)方案1中t＝0．41s和t＝0．84s時(shí)偏離度很大的峰值，從而在圖7（a）中其源項(xiàng)值的波動(dòng)幅度也遠(yuǎn)小于改進(jìn)方案1，進(jìn)而其壓力振蕩的幅度減小。結(jié)果表明，采用粒子插值不完整修正，能夠進(jìn)一步抑制MPS方法在模擬過(guò)程中的壓力振蕩現(xiàn)象。需要說(shuō)明的是，雖然采用改進(jìn)方案2時(shí)壓力振蕩仍然存在，但其振蕩的幅度已經(jīng)在可接受范圍內(nèi)，并且因?yàn)榱Ｗ拥目臻g分布不連續(xù)，MPS方法不可避免會(huì)出現(xiàn)一定程度的壓力振蕩。

圖8是傳統(tǒng)MPS方法、改進(jìn)方案1和改進(jìn)方案2模擬潰壩問(wèn)題在不同時(shí)刻的壓力場(chǎng)分布對(duì)比?？梢钥闯觯瑐鹘y(tǒng)MPS方法的壓力場(chǎng)空間分布極其不穩(wěn)定，出現(xiàn)了很多局部壓力驟增和驟減的區(qū)域，這是由于傳統(tǒng)MPS方法會(huì)產(chǎn)生壓力振蕩。而改進(jìn)方案1的空間壓力場(chǎng)相對(duì)于傳統(tǒng)MPS方法光滑許多，其原因是多方面的。如其選擇了對(duì)粒子稀疏程度不敏感的系數(shù)矩陣離散格式，能夠較好抑制壓力振蕩的二次樣條核函數(shù)以及混合的源項(xiàng)。但改進(jìn)方案1的空間壓力分布仍有一些變化不太光滑的區(qū)域。改進(jìn)方案2的空間壓力分布最為光滑，這是因?yàn)楦倪M(jìn)方案2在改進(jìn)方案1的基礎(chǔ)上，又增加了粒子插值不完整性的修正，從而有效地抑制了模擬過(guò)程中壓力變化驟增的區(qū)域，其空間壓力分布最為光滑。

圖7 A點(diǎn)處兩種方案下源項(xiàng)值和粒子數(shù)密度變化率隨時(shí)間的變化Fig．7 Comparison of time histories of source term value and change rate of particle number density at Apoint by different schemes

圖8 三種不同方案在不同時(shí)刻潰壩壓力分布對(duì)比Fig．8 Comparison of pressure field of dam－break with different MPS methods

圖9對(duì)比了傳統(tǒng)MPS方法和改進(jìn)方案2在潰壩模擬中的速度場(chǎng)?？梢钥闯?，兩種方法模擬潰壩時(shí)，獲得的速度場(chǎng)有著很大的不同。不同時(shí)刻，傳統(tǒng)MPS方法計(jì)算獲得的速度場(chǎng)在0m～0．6m位置處存在一個(gè)非常紊亂的區(qū)域，并且自由表面也不光滑。而采用改進(jìn)方案2所獲得的速度場(chǎng)非常光滑。在流體撞擊壁面前，流體從0m～1．2m位置處速度均勻地增加，流體前端的速度最大，而在尾部的流體保持靜止，并且流體的自由表面也非常光滑。通過(guò)對(duì)比可以看出，改進(jìn)方案2相比于傳統(tǒng)MPS方法可以獲得更為準(zhǔn)確和光滑的速度場(chǎng)。

3．3 計(jì)算時(shí)間對(duì)比

表2為三種方案在粒子數(shù)量分別為4592，7592和12328時(shí)的CPU計(jì)算時(shí)間，其中CPU計(jì)算時(shí)間指的是每循環(huán)10個(gè)時(shí)間步所需的時(shí)間?？梢钥闯?，隨著粒子數(shù)量的增加，三種方案的CPU計(jì)算時(shí)間都會(huì)增加。并且傳統(tǒng)MPS方法所需的時(shí)間最短；改進(jìn)方案1由于采用了新型的壓力泊松方程離散格式，其計(jì)算所需的時(shí)間有所增加；改進(jìn)方案2由于在改進(jìn)方案1的基礎(chǔ)上，增加了粒子插值完整性的修正，所需的計(jì)算時(shí)間最多。總的來(lái)說(shuō)，改進(jìn)的MPS方法并沒(méi)有大幅度地增加計(jì)算時(shí)間，為了得到較光滑的壓力場(chǎng)，這部分計(jì)算時(shí)間的增加是值得的。

表2 不同方案的CPU計(jì)算時(shí)間對(duì)比Tab．2 Comparison of computational time between three different schemes

3．4 靜水模型的計(jì)算

采用靜水模型檢驗(yàn)改進(jìn)方案2抑制壓力振蕩的效果。圖10為模擬區(qū)域，液面高度為0．48m，液體寬度為0．36m，容器總高度為0．6m。在流體底部中心位置處設(shè)置壓力監(jiān)測(cè)點(diǎn)P。本文靜水壓力的計(jì)算公式可表示為

式中P為該點(diǎn)處流體的壓力，ρ為流體密度，水的密度ρ＝1000kg／m3，g為重力加速度，取g＝9．8m／s2，ΔH為該點(diǎn)距離流體液面的垂直高度差。靜水模型中，粒子數(shù)量為12392。

3．5 靜水模型計(jì)算結(jié)果及分析

圖9 傳統(tǒng)MPS方法和改進(jìn)方案2模擬潰壩問(wèn)題時(shí)在不同時(shí)刻速度場(chǎng)的對(duì)比Fig．9 Comparison of velocity field of dam－break between traditional MPS method and modified MPS method 2

圖11 為不同MPS計(jì)算方案P點(diǎn)處的壓力－時(shí)間對(duì)比。從圖11（a）可以看出，原始的MPS方法模擬靜水問(wèn)題出現(xiàn)了劇烈振蕩，振蕩的幅度最大可以達(dá)到20000Pa。采用改進(jìn)方案2時(shí)，P點(diǎn)的壓力和解析解接近。但是P點(diǎn)的壓力有略微的振蕩現(xiàn)象，這是因?yàn)榱黧w粒子主要受到重力和壓力的作用，但是在計(jì)算的過(guò)程中，很難保證所有粒子受到的這兩種力時(shí)刻都保持相等，所以P點(diǎn)監(jiān)測(cè)得到的壓力有輕微的振蕩。從圖11（b）可以看出，采用改進(jìn)方案2模擬靜水問(wèn)題時(shí)，其壓力的振幅控制在400Pa以?xún)?nèi)，壓力振蕩相對(duì)于理論解的誤差保持在8%，在一個(gè)可接受范圍。圖12為傳統(tǒng)MPS方法和改進(jìn)方案2模擬靜水問(wèn)題時(shí)的壓力分布對(duì)比，可以看出，t＝10s時(shí)，傳統(tǒng)MPS方法模擬靜水問(wèn)題時(shí)，空間壓力分布非常不均勻，出現(xiàn)了非物理的壓力驟增和壓力驟減的區(qū)域，在自由表面處的粒子也處于劇烈運(yùn)動(dòng)狀態(tài)。而采用改進(jìn)方案2模擬靜水問(wèn)題時(shí)，流體的空間壓力分布非常光滑，底部流體壓力最大，并且隨著位置的升高，壓力逐漸減小，表面處的壓力為0。說(shuō)明采用改進(jìn)方案2模擬靜水問(wèn)題時(shí)，能夠得到光滑的空間壓力分布，并且其壓力振蕩得到了很好的抑制。

圖10 靜水模型示意圖Fig．10 Sketch of the hydrostatic problem

圖11 不同MPS計(jì)算方案P點(diǎn)處的壓力－時(shí)間對(duì)比Fig．11 Comparison of time histories of hydrostatic pressure at Ppoint by different schemes

圖12 不同方案模擬靜水問(wèn)題的壓力分布對(duì)比Fig．12 Comparison of pressure field of hydrostatic problem with different MPS methods

4 結(jié) 論

本文針對(duì)傳統(tǒng)MPS方法中壓力振蕩問(wèn)題，采用了一種新型的壓力泊松方程離散格式來(lái)緩解MPS方法在模擬過(guò)程中的壓力振蕩，新型的壓力泊松方程格式包括對(duì)粒子稀疏程度不敏感的壓力系數(shù)矩陣離散格式，以及能夠有效抑制壓力振蕩的混合源項(xiàng)。選擇了二次樣條核函數(shù)，該核函數(shù)能夠有效地改善計(jì)算的穩(wěn)定性。并且針對(duì)MPS方法在模擬過(guò)程中粒子插值不完整問(wèn)題，提出了粒子插值不完整性的修正。通過(guò)采用上述方法，對(duì)經(jīng)典的潰壩問(wèn)題以及靜水問(wèn)題進(jìn)行了模擬計(jì)算。結(jié)果表明，采用新型的壓力泊松方程格式、二次樣條核函數(shù)和粒子插值不完整性修正后，計(jì)算得到的空間壓力場(chǎng)變得更為光滑，并且與實(shí)際的物理現(xiàn)象吻合。通過(guò)對(duì)模擬過(guò)程中檢測(cè)點(diǎn)壓力的采集，所得的壓力－時(shí)間曲線與實(shí)驗(yàn)值和理論值十分接近，本文的改進(jìn)方法有效地抑制了MPS方法在模擬過(guò)程中的壓力振蕩。