負壓梯度下自由面流動模擬的無網格粒子法Dirichlet邊界條件

韓沛東,劉啟新,周子棋,孫中國,席光

(西安交通大學能源與動力工程學院,710049,西安)

液滴或氣泡的破裂與融合等自由表面流動是工業(yè)過程常見的流體動力學現(xiàn)象。由于自由表面常存在劇烈運動變形,傳統(tǒng)基于網格剖分的模擬中多采用Level Set或VOF(volume of fluid)等數(shù)值方法,界面重構和尖角精確捕捉仍較復雜;基于拉格朗日框架的無網格粒子法避免了所有基于網格的限制[1]。光滑粒子流體動力學方法SPH(smoothed particle hydrodynamics)及移動粒子半隱式法MPS(moving particle semi-implicit method)等已廣泛應用于計算大變形非定常流動[2-4]。在不可壓流動計算中,與大氣連通的自由表面粒子壓力值常賦為0,可準確捕捉自由界面,但在負壓梯度(即流體壓力低于大氣壓,或界面法向壓力梯度力朝向流體內部)算例中,自由表面粒子的核函數(shù)截斷誤差會影響方法的計算精度和穩(wěn)定性[5],造成諸如粒子聚集、表面誤判、壓力求解誤差增大等問題[6],是無網格粒子法研究的難點和熱點之一。

國內外學者圍繞算子修正、粒子分布重構、表面粒子判定、自由邊界修正等方面開展了一系列改善研究。其中,Koshizuka等[7-8]采用粒子穩(wěn)定項PST提升梯度算子穩(wěn)定性;Liu等[9-10]提出用一階、二階修正矩陣提高算法在粒子分布不均勻時的精度。Khayyer等[11]基于粒子遷移技術PS(particle shifting)提出了優(yōu)化的粒子遷移技術OPS(optizmized PS),通過粒子位置遷移調整(重構),顯著提升了典型算例的收斂性;Jandaghian等[12]等提出了修正的粒子遷移技術CPS(corrected PS),與粒子碰撞模型結合,避免了粒子聚集并降低了數(shù)值噪聲。Chen等[13-16]提出概念粒子模型,開發(fā)了無需自由表面粒子判定的MPS-NSD(no surface detection)方法,獲得了進一步發(fā)展與應用;Khayyer等[11]提出了基于散度的表面粒子判定,提升了判定的準確性;朱躍等[17]提出了基于幾何法和體積法的表面判定方式。Liu等[18-19]通過加入虛粒子修正自由表面及附近粒子梯度算子;Matsunaga等[20]通過移動的曲面網格來表示可變形的自由面邊界,提升了計算精度;殷競存等[21]采用改進的ghost粒子方法處理自由表面,獲得較平滑的邊界。上述工作對負壓自由面流動的求解精度和收斂性起到一定提升作用。

然而,無網格粒子法通常對自由表面的狄利克雷邊界條件(Dirichlet boundary condition)[1]直接賦值,將自由表面粒子所處位置即視為自由面邊界,忽略了自由表面粒子不規(guī)則排布帶來的空間離散誤差。實際計算中,這一誤差并不總能通過前述技術得到完全修正[22]。

方形液滴旋轉算例[23]是檢驗負壓條件下自由面大變形計算精度與穩(wěn)定性的經典算例。相關研究大都聚焦高階壓力泊松方程和梯度、散度、拉普拉斯算子等方面,而本文利用該算例研究自由表面核截斷區(qū)域計算精度。

本文提出了一種基于虛擬邊界的狄利克雷邊界條件,并在MPS框架下通過方形液滴旋轉算例進行了驗證,該方法可顯著改善大變形復雜曲面的求解精度和穩(wěn)定性。

1 數(shù)值方法

1.1 MPS方法

MPS方法[1]用于模擬不可壓縮流動,在拉格朗日描述下,通過配點法實現(xiàn)空間離散。不可壓縮流體的控制方程包括質量守恒方程與動量守恒方程

ρ·u=0

(1)

(2)

式中:ρ為流體密度;u為流體速度;p為壓力;μ為動力黏度;f為體積力。

流體被表征為具有一定質量的離散粒子,粒子攜帶位置、速度、壓力等物理量信息,通過粒子間相互作用,離散控制方程中的各項微分算子,導出梯度算子和拉普拉斯算子的光滑近似式。

MPS方法采用核函數(shù)模型描述粒子間作用,并憑此離散控制方程,本文采用如下核函數(shù)[6]

(3)

式中:w(r)為核函數(shù)值;r=|ri-rj|為兩粒子i、j之間距離;re為光滑半徑,其值為2.1倍粒子直徑l0。

原始MPS方法中,負壓會帶來數(shù)值不穩(wěn)定[23],通常將負壓調整為零或更換更穩(wěn)定的梯度公式[24]來實現(xiàn)僅有斥力的壓力梯度計算。當計算存在負壓與吸引力時,上述解決方式不再適用,因此本文采用帶修正矩陣的梯度公式[23]

〈φ〉i=

(4)

式中:Cij為在引入的一個修正項,形式為一個矩陣的逆矩陣。在二維問題中,Cij可表示為

(5)

拉普拉斯算子計算公式為

(6)

式中:φ為任意標量函數(shù);d為維度(二維d=2);n0為粒子數(shù)密度常數(shù)。在忽略連續(xù)域內核函數(shù)截斷誤差的情況下,拉普拉斯算子可簡化為[9]

(7)

(8)

1.2 優(yōu)化的粒子遷移技術

在使用上述模型時,粒子傾向于沿流線分散,會導致計算不穩(wěn)定,常需采用OPS技術[11]來減小粒子分布的不均勻性。參考菲克第一定律,將高濃度粒子向低濃度粒子區(qū)域進行微小移動,即

(9)

(10)

式中:Cshift≤0.5為遷移系數(shù);d0=l0;Ci為粒子分布濃度的梯度。

自由表面粒子以及鄰近表面的內部粒子僅保留沿自由面切向的遷移分量為

(11)

(12)

(13)

粒子遷移后,憑借原流場參數(shù)更新速度場,即

ui′=ui+ui·δrii′

(14)

OPS技術使計算穩(wěn)定性明顯提升,但同時也導致表面粒子斥力失效[12],在計算含負壓的自由面流動時,仍然會影響收斂性。

1.3 表面粒子判定

若要保證壓力泊松方程計算結果準確,精確地判定自由表面粒子并施加Dirichlet邊界條件至關重要。為保證表面粒子應判盡判,同時又避免內部粒子被誤判為表面粒子,本文采用了多種判定方法組合的形式。

傳統(tǒng)MPS方法[1]將粒子數(shù)密度小于αn0的粒子視為表面粒子,但常出現(xiàn)內部粒子誤判為表面粒子的情況。將粒子位置矢量rij的散度·r作為評判標準[11],對于粒子均勻分布時的內部粒子該值為2。文獻[16]提出虛擬光源和光幕的概念,將粒子i視為光源,鄰近粒子視為遮擋物,將未遮擋面積占比Rill作為判定標準。

本文將上述方法相結合,對于粒子數(shù)密度小于0.99n0或·r<1.9的粒子進行光源遮擋判定,若Rill<1/6則視為表面粒子。若粒子數(shù)密度小于0.85n0或·r<1.5,則將該粒子強制判定為表面粒子。對于一個表面粒子,若其影響域內沒有內部粒子,則將其視為游離粒子,不參與自由表面的計算。判定結果如圖1所示。

圖1 自由表面附近粒子分布及虛擬邊界示意圖

2 虛擬邊界

2.1 機理分析

傳統(tǒng)MPS方法將自由表面粒子視為實際邊界,并在施加狄利克雷邊界條件時將其壓力賦值為0,如圖1所示。針對不可壓縮流體,被自由表面粒子包裹的內部流體粒子的粒子數(shù)密度n保持恒定值n0,且在每個時間步粒子流動發(fā)生位移后,都會修正為n0。粒子間的相對位置決定了粒子間的相互作用,進而影響流動。

若自由表面粒子位置分布不夠平滑,其用于描述實際自由面的形狀便會出現(xiàn)誤差,進而影響計算的精度和穩(wěn)定性。單個自由表面粒子沿表面垂直方向發(fā)生偏移時,一般會產生以下兩種情況。

(1)當壓力梯度力朝向自由面外部時(如圖1中F所示),若該粒子產生向外的偏移,壓力梯度力減小,粒子之間斥力減小,該粒子向外偏移的趨勢也會隨之減小。反之,若該粒子產生向內的偏移,梯度力增大,粒子之間斥力增加,將遏制該粒子繼續(xù)向內偏移。此時誤差被自然修正,不會積累。

(2)當壓力梯度力朝向自由面內部時(如圖1中F′所示),若該粒子產生向外的偏移,壓力梯度力減小,誤差將逐漸擴大,該粒子會飄離自由表面,如圖2(a)所示。反之,若粒子產生向內的偏移,粒子靠近,梯度力增大,可能導致局部高壓無法抵消梯度力增大帶來的吸引力,引起粒子繼續(xù)向自由面內部擠壓,產生非物理的局部高壓,并最終導致計算發(fā)散,如圖2(b)所示。

(a)粒子飄散 (b)局部高壓失穩(wěn)

2.2 方法原理

本文提出一種基于虛擬邊界的狄利克雷邊界條件賦值方法,其基本原理為在原計算的邊界分布基礎上構建表征真實表面邊界的虛擬表面邊界,通過計算表面粒子與虛擬邊界的偏差量來對狄利克雷邊界條件賦值進行修正。其中,實際的自由表面邊界由表面粒子及其鄰域搜索范圍內的其他表面粒子連接而成,局部常為一條短曲線段,整個自由表面由許多這樣的小曲線段聯(lián)接組成。

每個短曲線段用二次曲線擬合,表面粒子位于擬合曲線上。自由表面的法向量與擬合曲線的對稱軸重合,則擬合二次曲線y軸的朝向為該表面粒子的法向量方向。

2.3 算法流程

狄利克雷邊界條件改進算法流程如圖3所示。

圖3 狄利克雷邊界條件改進算法流程

對于某自由表面粒子i,搜索其附近的其他自由表面粒子位置,同時計算i粒子的法向量,本文采用式(12)計算修正的法向量,采用坐標旋轉公式

(15)

坐標旋轉變換關系如圖4所示,以i粒子為原點,對附近表面粒子進行坐標旋轉,使y軸與法向量重合,θ為法向量與y軸的夾角。

圖4 坐標旋轉變換關系

粒子坐標旋轉后,采用最小二乘法求解待擬合二次曲線的各項系數(shù),并聯(lián)立求解各項系數(shù),即可得擬合二階曲線方程

(16)

Q函數(shù)極值點對各項系數(shù)的導數(shù)為0,即

(17)

(18)

通過以上二階多項式擬合,計算得到偏差最小的擬合曲線,即虛擬邊界。憑借虛擬邊界方程計算得出邊界到該粒子的偏移矢量(0,b0)。并將偏移矢量沿-θ方向旋轉回原始笛卡爾坐標,得到i粒子到虛擬邊界的距離矢量δi。

最終,通過該表面粒子與此曲線的偏差和已求得的流場參數(shù)梯度進行插值計算,得出該表面粒子應當賦予的修正參數(shù)值

φi=φs+φi·δi

(19)

式中:φi為i粒子賦予的參數(shù)值;φs為狄利克雷邊界條件預設的參數(shù)值;φi為i粒子處參數(shù)φ的梯度。

為保障二階多項式的擬合精度,一般情況下,需要5～7個數(shù)據點。本文采用re≥2.5l0作為該模型搜索半徑,鄰近表面粒子與加上i粒子自身,滿足上述擬合條件。

若出現(xiàn)i粒子附近的表面粒子分布與二次曲線分布差異較大的情況,可先對結果的方差進行校核,并在對狄利克雷邊界條件賦值時,對b0的大小進行限制,令

b0=min(b0,δmax)

(20)

δmax取值影響自由表面粒子在虛擬邊界附近松散分布程度,本研究中取值為0.6l0。實際模擬過程中可根據具體算例需求進行調整,推薦的取值范圍為0.1l0～0.8l0。

此外,上述狄利克雷邊界條件的賦值方法也可用于壓力出口邊界,采用改進的狄利克雷邊界條件對表面粒子進行壓力進行賦值

pi=ps+pi·δi

(21)

式中:ps為壓力邊界給出的壓力值;pi可取i粒子的壓力梯度或其上游2l0距離處的粒子的壓力梯度。僅對距離壓力出口一定范圍內的表面粒子進行上述修正,可獲得更平滑的壓力分布。

將上述方法推廣至三維時,需采用最小二乘法對粒子i附近的表面粒子進行曲面擬合,擬合二階曲面時須采用的待定系數(shù)將由10個縮減至6個。

3 數(shù)值模擬驗證

3.1 典型算例及設置

方形液滴旋轉算例[12]是二維方形液滴在無粘和無表面張力的情況下以幾何中心為轉軸發(fā)生轉動,在離心力作用下方塊四角向外拉伸形成彎曲尖突,中心部分形成明顯負壓,四邊發(fā)生顯著內縮的劇烈形變過程。該算例中自由面平滑、尖突清晰、壓力由外向內形成負壓且成圓形等特征是檢驗負壓計算和算法穩(wěn)定性的經典算例。

圖5 方形液滴旋轉算例

3.2 模擬結果

改進后方形液滴旋轉算例壓力云圖及與文獻結果比較如圖6所示,其中方框為初始液滴位置,斜線為各角點的理論運行軌跡。圖6(a)～(d)展示了0.5～2.0 s方形液滴的變形過程,與文獻[18]中計算結果相比,形態(tài)和壓力均吻合較好。但相比文獻,使用偽粒子重整表面技術和更高精度的復雜模型,本文模型及實現(xiàn)過程更加簡單快捷。

(a)t=0.5 s (b)t=1.0 s (c)t=1.5 s (d)t=2.0 s (e)文獻結果

圖7展示了添加修正矩陣和粒子遷移技術的方形液滴旋轉計算結果,以及在此基礎上添加了本文改進方法后的計算結果。結果表明,修正矩陣和改進的粒子遷移技術可在一定程度上避免拉伸失穩(wěn)誘發(fā)的表面粒子聚集,但無法完全避免梯度力不足導致的粒子飄散問題。

采用本文提出的改進方法后(如圖7(b)所示),表面粒子飄散現(xiàn)象得到明顯抑制,液滴外形演變過程與理想結果相符,液滴內部的壓力分布也變得更合理。t=2.0 s時,液滴中心壓力為-54.53 Pa,換算為無量綱壓力約為-0.055,與文獻[11,18,22]計算結果相符。但四邊內凹處邊界曲線仍存在明顯毛刺,考慮與粒子搜索半徑選取相關。圖7(c)展示了調整擬合曲線搜索半徑由re=2.1l0增大為4l0后,曲線擬合數(shù)據由5～7個點進一步擴充至7～9個點,曲線擬合精度提升,邊界毛刺現(xiàn)象基本消失。

(a)改進前 (b)改進后re=2.1l0 (c)改進后re=4l0

4 結 論

本文針對負壓梯度條件下自由面大變形流動問題,通過構建虛擬邊界改進了拉格朗日框架下的狄利克雷邊界條件,模擬了方形液滴旋轉典型算例,并得到以下結論。

(1)通過采用最小二乘法對自由表面位置進行二階多項式擬合,構建虛擬邊界,憑借虛擬邊界與表面粒子的偏移矢量修正粒子賦值,改進了狄利克雷邊界條件,顯著提高了自由表面核截斷區(qū)域的計算精度與穩(wěn)定性。

(2)采用上述方法模擬方形液滴旋轉算例,無需更高階精度的復雜模型便可獲得較為理想的結果。大變形曲面外緣粒子因梯度力偏差而導致的粒子飄離表面情況得到顯著抑制,液滴內部壓力分布更接近預期。驗證了方法的有效性與準確性。

(3)通過對上述方法中擬合偏移量b0的誤差進行約束和修正,以及增加虛擬邊界擬合曲線的搜索范圍至re=4l0,可顯著減小虛擬邊界誤差。

本文提出的狄利克雷邊界條件改進方法為準確捕捉負壓梯度自由面大變形流動的界面、提升復雜界面流動計算的穩(wěn)定性提供了新的思路。未來將繼續(xù)擴展至三維情況并與相界面、表面張力等數(shù)值模型耦合。