張凱,孫中國(guó),席光

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,710049,西安)

移動(dòng)粒子半隱式法(MPS)是一種無(wú)網(wǎng)格數(shù)值方法,相比于傳統(tǒng)網(wǎng)格剖分方法,移動(dòng)粒子半隱式法采用拉格朗日描述,在處理大變形問(wèn)題時(shí),只需追蹤粒子攜帶的位置、壓力和速度等信息,避免了傳統(tǒng)網(wǎng)格方法容易產(chǎn)生網(wǎng)格畸變以及需要頻繁網(wǎng)格重構(gòu)的情況,因而在處理具有大變形復(fù)雜邊界或移動(dòng)邊界的問(wèn)題時(shí)具有明顯優(yōu)勢(shì)[1]。近年來(lái),利用MPS法模擬楔形體入水、多孔介質(zhì)流動(dòng)等自由表面流體流動(dòng)現(xiàn)象的方法得到了廣泛的研究與應(yīng)用[2-4]。流體壓力采用隱式求解,求解精度和穩(wěn)定性將直接影響MPS算法的整體精度,壓力求解的數(shù)值振蕩容易導(dǎo)致模擬精度下降甚至不收斂,因此有效抑制MPS法在流動(dòng)模擬中的壓力振蕩對(duì)提升算法精度與穩(wěn)定性具有重要意義。

國(guó)內(nèi)外學(xué)者已從多方面開展了抑制MPS算法壓力振蕩的研究,包括壓力梯度修正、表面粒子判定、添加泊松方程源項(xiàng)等。Khayyer等通過(guò)在壓力梯度以及泊松方程中添加源項(xiàng)并允許輕微的可壓縮性,提出了改進(jìn)MPS法用來(lái)提升預(yù)測(cè)沿海結(jié)構(gòu)波浪沖擊壓力的精度[5]。Shibata等提出了在自由表面上布置虛擬粒子的方法來(lái)抑制壓力振蕩,該方法同時(shí)可模擬負(fù)壓現(xiàn)象[6]。Debadatta等通過(guò)修改壓力泊松方程源項(xiàng),并采用混合表面粒子判定技術(shù)以及改進(jìn)的固液邊界條件,獲得較穩(wěn)定的壓力求解,模擬了劇烈晃動(dòng)現(xiàn)象[7]。

上述學(xué)者對(duì)壓力振蕩的改進(jìn)主要集中在壓力梯度以及泊松方程方面,對(duì)核函數(shù)的研究相對(duì)較少,實(shí)際上,核函數(shù)數(shù)值特征對(duì)模擬的穩(wěn)定性具有重要影響。Ataie-Ashtiani等對(duì)不同核函數(shù)提高潰壩模擬穩(wěn)定性做了相關(guān)研究,發(fā)現(xiàn)緊支域半徑以及粒子大小對(duì)模擬穩(wěn)定性具有較大影響[8]。潘徐杰等發(fā)現(xiàn)在粒子相撞時(shí)具有有限值的核函數(shù)可減緩壓力振蕩,并研究了5種常見核函數(shù)對(duì)模擬的適用性[9-10]。Shobeyri發(fā)現(xiàn),與具有陡峭形狀的核函數(shù)相比,扁平形狀的核函數(shù)可提高模擬的穩(wěn)定性[11]。朱躍等采用兩種核函數(shù)模擬了潰壩問(wèn)題,與傳統(tǒng)核函數(shù)相比,采用高斯型核函數(shù)模擬得到的粒子數(shù)密度偏差波動(dòng)較小,壓力振蕩幅度較小[12]。

本文分析了MPS法核函數(shù)形狀特征對(duì)壓力求解穩(wěn)定性的影響,通過(guò)歸納總結(jié)主要特性,提出了一種可有效抑制壓力振蕩的新核函數(shù)。針對(duì)靜水壓力問(wèn)題的靜壓求解以及液體晃動(dòng)問(wèn)題的動(dòng)力學(xué)壓力計(jì)算,采用典型算例進(jìn)行驗(yàn)證。

1 數(shù)值方法

1.1 移動(dòng)粒子半隱式法

移動(dòng)粒子半隱式法中,流體在流動(dòng)過(guò)程中遵循能量守恒方程和動(dòng)量守恒方程,矢量形式為

(1)

(2)

式中:ρ為密度;t為時(shí)間;u為速度矢量;p為壓力;μ為動(dòng)力黏度;f為重力;為梯度算子;2為拉普拉斯算子。算子式分別為

(3)

(4)

(5)

式中:d為維度;〈〉i為粒子i物理量的近似值;xij為粒子i與j之間的x向距離;wij為核函數(shù);r=|ri-rj|,ri、rj為粒子i、j的坐標(biāo);re為核函數(shù)的影響半徑,一般情況取2.1l0,l0為單個(gè)流體粒子的直徑。

粒子i的粒子數(shù)密度可表示為

(6)

對(duì)于不可壓縮流體,粒子數(shù)密度保持不變,取初始布置時(shí)的粒子數(shù)密度n0。

粒子壓力可通過(guò)求解壓力泊松方程得到

(7)

為了使計(jì)算結(jié)果更穩(wěn)定,調(diào)整泊松方程的源項(xiàng)[13],該方法可獲得更穩(wěn)定的壓力場(chǎng),即

(8)

式中γ為(0,1)之間的常值。

1.2 核函數(shù)基本要素與構(gòu)造思路

在MPS算法中,粒子間的相互作用用核函數(shù)加權(quán)表示,傳統(tǒng)核函數(shù)和加權(quán)形式參見式(5)(6)。當(dāng)粒子間距r>re時(shí),核函數(shù)值為0,忽略粒子間相互作用,核函數(shù)滿足緊支性條件并可提高計(jì)算效率。計(jì)算中物理參數(shù)值可通過(guò)核函數(shù)加權(quán)計(jì)算獲得,核函數(shù)滿足非負(fù)性條件可防止出現(xiàn)非物理現(xiàn)象。為了描述粒子間相互作用的實(shí)際情況,即粒子間距離越小,粒子間的相互作用越大,核函數(shù)須滿足單調(diào)遞減條件。

傳統(tǒng)核函數(shù)在粒子間距非常小時(shí)值為無(wú)窮大,導(dǎo)致在壓力方程建立和求解過(guò)程中容易出現(xiàn)數(shù)值振蕩現(xiàn)象。眾多學(xué)者在模擬過(guò)程中會(huì)采用不同形狀核函數(shù)以減輕模擬過(guò)程中的壓力振蕩現(xiàn)象,一種常見的核函數(shù)處理方式是在粒子相互碰撞時(shí),將核函數(shù)的值改為有限值[9-11];另外,核函數(shù)的形狀特征對(duì)壓力求解穩(wěn)定性的影響是本文關(guān)注的重點(diǎn)?；谏鲜龇治?為便于把握曲線的形狀特征,本文構(gòu)造有限值核函數(shù)及以指數(shù)函數(shù)為基礎(chǔ)的多項(xiàng)式形式核函數(shù),重點(diǎn)研究核函數(shù)的極值和形狀特征,并分別通過(guò)流體靜力學(xué)和流體動(dòng)力學(xué)中兩個(gè)典型算例,即靜水壓力求解和液體晃動(dòng)壓力求解,來(lái)分析并驗(yàn)證指數(shù)多項(xiàng)式型核函數(shù)對(duì)抑制壓力振蕩現(xiàn)象的效果。

2 靜水壓力求解

構(gòu)造兩種有限值核函數(shù),基于靜水壓力問(wèn)題模擬,對(duì)比分析不同核函數(shù)在模擬過(guò)程中對(duì)壓力求解穩(wěn)定性的影響。分別從核函數(shù)是否為有限值、核函數(shù)表征形狀特點(diǎn)以及核函數(shù)真實(shí)形狀特征3個(gè)方面,分析歸納抑制壓力振蕩的核函數(shù)特性。本文算例為二維模型,重力加速度g=-9.8 m/s2。

2.1 算例形狀及初始布置

算例模擬區(qū)域和粒子初始布置如圖1所示,其中水池底面中點(diǎn)p1為壓力監(jiān)測(cè)點(diǎn),l0=0.008 m,液體部分離散為2 700個(gè)粒子,壁面粒子布置為3層。p1點(diǎn)壓力的解析解為

p=ρgh

(9)

式中h為水池深度。

(a)幾何模型 (b)粒子初始離散模型圖1 靜水問(wèn)題模型

2.2 討論及結(jié)果分析

2.2.1 有限值核函數(shù) 分別采用3種不同核函數(shù)模擬靜水壓力問(wèn)題,其中w1(r)為傳統(tǒng)核函數(shù),在粒子間距為0時(shí)核函數(shù)值為無(wú)窮大,w2(r)為有限值核函數(shù)[14],w3(r)為本文構(gòu)造的指數(shù)型多項(xiàng)式核函數(shù),為有限值且在整個(gè)支持域內(nèi)無(wú)奇點(diǎn),即

(10)

(11)

3種核函數(shù)的曲線如圖2所示。3種不同核函數(shù)模擬靜水問(wèn)題時(shí),監(jiān)測(cè)點(diǎn)p1點(diǎn)的壓力隨時(shí)間變化的曲線如圖3所示。由圖3可知:傳統(tǒng)核函數(shù)w1(r)壓力振蕩比較明顯,最終計(jì)算不收斂;核函數(shù)w2(r)雖在粒子相互碰撞時(shí)采用了有限值,但未能明顯抑制壓力振蕩;新構(gòu)造核函數(shù)w3(r)壓力振蕩得到明顯抑制,且模擬結(jié)果與解析解吻合較好。

圖2 核函數(shù)w1(r)、w2(r)、w3(r)曲線

圖3 3種核函數(shù)模擬靜水問(wèn)題時(shí)p1點(diǎn)壓力隨時(shí)間的變化

采用w1(r)時(shí),當(dāng)相鄰粒子間距較小時(shí),其核函數(shù)值將迅速增大為較大數(shù)值(趨近無(wú)窮大),造成相鄰粒子壓力巨大,發(fā)生粒子飛濺彈開的數(shù)值現(xiàn)象,這與實(shí)際物理過(guò)程不符。當(dāng)采用有限值核函數(shù)w2(r)時(shí),壓力振蕩仍然較嚴(yán)重,這是由于當(dāng)粒子間距很小時(shí),核函數(shù)值隨粒子間距變化依然比較敏感(大斜率),這使得粒子間的作用關(guān)系在間距較小時(shí)呈現(xiàn)出不穩(wěn)定的趨勢(shì),成為誘發(fā)壓力振蕩的重要因素。

本文提出變化率相對(duì)較小的w3(r),在粒子間距較小時(shí),其核函數(shù)值增大至合理范圍,與周邊粒子及前一時(shí)刻自身的核函數(shù)值相差較小,不會(huì)產(chǎn)生過(guò)大壓力,粒子運(yùn)動(dòng)相對(duì)穩(wěn)定,避免了劇烈飛濺。結(jié)果分析表明,壓力求解穩(wěn)定的核函數(shù)須在滿足粒子間距與核函數(shù)值成反比例關(guān)系的基礎(chǔ)上,同時(shí)應(yīng)具備在粒子間距較近,即r/re為0～0.1時(shí)為有限值且變化率較小的形狀特征。

2.2.2 核函數(shù)表征形狀特點(diǎn) 為便于探究核函數(shù)曲線的形狀特點(diǎn),尤其是[0,1]區(qū)域兩端曲線斜率對(duì)壓力振蕩的影響,本文構(gòu)造核函數(shù)采用w(r)=exp(-a(r/re)2)-exp(-a)的形式,其中a為常數(shù),當(dāng)r/re=1時(shí)核函數(shù)值為0,保證了核函數(shù)的連續(xù)性。a取不同的值時(shí),核函數(shù)曲線具有不同的形狀特點(diǎn),a越小核函數(shù)形狀越呈現(xiàn)出扁平狀,核函數(shù)值在0.9

采用核函數(shù)w4(r)(a=3)、w5(r)(a=6)和w6(r)分別進(jìn)行數(shù)值模擬,如圖4所示,監(jiān)測(cè)點(diǎn)p1歷時(shí)10 s,壓力隨時(shí)間變化的曲線如圖5所示。核函數(shù)w4(r)中a=3時(shí),壓力振蕩較明顯;核函數(shù)w5(r)中a=6時(shí),壓力與解析解一致,隨時(shí)間分布均勻且較穩(wěn)定。

圖4 核函數(shù)w4(r)、w5(r)、w6(r)曲線

圖5 p1點(diǎn)壓力隨時(shí)間的變化

由圖4、圖5對(duì)比可知,核函數(shù)值在r/re=1附近變化率較大時(shí),容易出現(xiàn)壓力求解不穩(wěn)定。若粒子間距較大時(shí)核函數(shù)值較大,粒子相互接近的趨勢(shì)將會(huì)受到抑制,因此在r/re=1附近,核函數(shù)值應(yīng)始終維持在0附近,核函數(shù)變化率不應(yīng)過(guò)大。壓力求解穩(wěn)定的核函數(shù)具有在r/re=0處為有限值,且在r/re=0,1兩端處應(yīng)盡量保持平緩變化的形狀特點(diǎn)。

2.2.3 核函數(shù)真實(shí)形狀特征 構(gòu)造核函數(shù)采用wm(r)=bexp(-6(r/re)2)-bexp(-6)的形式,其中b取不同值時(shí),核函數(shù)wm(r)在粒子間距為0時(shí)有限值的大小不同。當(dāng)b分別取1、5、9時(shí),核函數(shù)w7(r)、w8(r)、w9(r)曲線如圖6所示,核函數(shù)與其對(duì)應(yīng)粒子數(shù)密度常數(shù)的比值wn(r)=wm(r)/n0?？芍猵1點(diǎn)壓力與b無(wú)關(guān),b取不同數(shù)值時(shí),p1點(diǎn)壓力完全一致。

圖6 核函數(shù)w7(r)～w12(r)曲線

SPH算法中光滑函數(shù)的構(gòu)造需要滿足歸一性條件[15],而MPS算法并不需要強(qiáng)制滿足,物理量參數(shù)的近似值可表示為

(12)

MPS算法中粒子數(shù)密度作為規(guī)范化因子,對(duì)粒子各作用模型進(jìn)行規(guī)范化處理。

當(dāng)b取3種不同數(shù)值時(shí),wn(r)所表達(dá)的3條曲線完全重合,模擬得到包含p1點(diǎn)在內(nèi)的壓力值分布完全一致。核函數(shù)與粒子數(shù)密度常數(shù)的比值在壓力求解時(shí)起重要作用,盡管核函數(shù)曲線形狀和極值不同,但wn(r)的曲線能真實(shí)反映粒子間的相互作用。

3 液體晃動(dòng)求解

3.1 初始布置

液體晃動(dòng)動(dòng)壓?jiǎn)栴}模擬區(qū)域和初始布置如圖7所示,其中液箱在水平方向上做正弦運(yùn)動(dòng),液箱運(yùn)動(dòng)規(guī)律為vx=0.02×5.39sin(5.39t),p2點(diǎn)為液箱右側(cè)壁一點(diǎn),p2點(diǎn)與液箱底部距離為55 mm,粒子直徑l0=0.005 m,布置8 000個(gè)液體粒子,壁面粒子布置為3層。

3.2 討論及結(jié)果分析

通過(guò)靜水壓力模擬發(fā)現(xiàn),核函數(shù)在[0,1]區(qū)間曲線兩端附近數(shù)值變化較平緩時(shí),可有效改善壓力求解振蕩現(xiàn)象;當(dāng)核函數(shù)在r/re=0.8附近數(shù)值過(guò)小時(shí),系統(tǒng)動(dòng)力學(xué)計(jì)算會(huì)出現(xiàn)偏差。為了兼顧壓力求解穩(wěn)定性與保證系統(tǒng)的動(dòng)力學(xué)性能,對(duì)核函數(shù)做進(jìn)一步改進(jìn),可提高r/re=0.8附近核函數(shù)值,新的分段指數(shù)多項(xiàng)式型核函數(shù)為

(13)

核函數(shù)改進(jìn)前后典型時(shí)刻液箱內(nèi)自由液面模擬結(jié)果與實(shí)驗(yàn)的對(duì)比如圖8所示。改進(jìn)前流體壓力值計(jì)算存在偏差且有瞬態(tài)波動(dòng),如圖中A、B、C區(qū)域所示,同時(shí)流體內(nèi)部出現(xiàn)較多壓力零散點(diǎn);改進(jìn)后壓力振蕩明顯減弱,壓力分層分布更加合理,壓力零散點(diǎn)基本消失,如圖中A′、B′、C′區(qū)域所示。與w1(r)相比,w6(r)在壓力求解過(guò)程中具有更好的穩(wěn)定性,液體在晃動(dòng)過(guò)程中自由液面的變化規(guī)律與實(shí)驗(yàn)結(jié)果更加吻合。

(a)幾何模型

(b)粒子初始離散模型圖7 液體晃動(dòng)模型

改進(jìn)前后液體晃動(dòng)算例監(jiān)測(cè)點(diǎn)p2的壓力計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果[16]對(duì)比如圖9所示。采用改進(jìn)前傳統(tǒng)核函數(shù),模擬結(jié)果存在較大振蕩,造成較大誤差,雙峰現(xiàn)象較模糊;改進(jìn)后壓力振蕩被抑制且成功捕獲雙峰現(xiàn)象,與實(shí)驗(yàn)結(jié)果吻合較好。壓力峰值計(jì)算較為準(zhǔn)確,谷值較實(shí)驗(yàn)值略大,需進(jìn)一步分析液體晃動(dòng)中數(shù)值耗散產(chǎn)生的影響。指數(shù)多項(xiàng)式型核函數(shù)在整個(gè)區(qū)間內(nèi)連續(xù)且最大值為有限值,函數(shù)構(gòu)造簡(jiǎn)單易于實(shí)現(xiàn),可有效抑制壓力振蕩現(xiàn)象,與傳統(tǒng)核函數(shù)相比,壓力求解具有更好的穩(wěn)定性,可廣泛適用于普通流體的模擬。

t=nT+0.2T t=nT+0.5T t=nT+0.7T(a)實(shí)驗(yàn)結(jié)果[16]

t=6.06 s t=6.41 s t=6.64 s(b)w1(r)模擬得到的壓力分布云圖

t=6.06 s t=6.41 s t=6.64 s(c)w6(r)模擬得到的壓力分布云圖圖8 不同時(shí)刻液箱內(nèi)自由液面實(shí)驗(yàn)結(jié)果與模擬結(jié)果對(duì)比

圖9 改進(jìn)前后p2點(diǎn)壓力變化的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果[16]比較

4 結(jié) 論

本文分析并研究了核函數(shù)曲線形狀特點(diǎn)對(duì)MPS法流體壓力求解穩(wěn)定性的影響,構(gòu)造了指數(shù)多項(xiàng)式型核函數(shù),通過(guò)模擬靜水壓力及液體晃動(dòng)典型問(wèn)題,驗(yàn)證了所提核函數(shù)可準(zhǔn)確模擬液體的壓力分布及監(jiān)測(cè)點(diǎn)的壓力變化,并可有效抑制求解過(guò)程中的壓力振蕩現(xiàn)象。研究發(fā)現(xiàn),核函數(shù)的最大值為有限值且r/re在[0,1]區(qū)間曲線兩端附近數(shù)值變化平緩時(shí),可有效改善壓力求解振蕩現(xiàn)象,在r/re=0.8附近核函數(shù)值過(guò)小會(huì)影響系統(tǒng)的動(dòng)力學(xué)性能。此外,核函數(shù)的極值及曲線形狀僅是表征特點(diǎn),核函數(shù)與粒子數(shù)密度的比值可真實(shí)表征粒子間的相互作用,在壓力分析穩(wěn)定性中起重要作用。