基于速度修正的固壁邊界處理方法

朱曉臨，陳 薇

基于速度修正的固壁邊界處理方法

朱曉臨，陳 薇

(合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，安徽 合肥 230009)

固壁邊界處理方法的研究一直是流體模擬中的難點(diǎn)問題，常見的固壁邊界處理方法有邊界力法和虛粒子法。邊界力法通過對靠近邊界的流體粒子施加作用力防止其穿透邊界，但模型參數(shù)較多，力的大小難以調(diào)控，且在計(jì)算中會產(chǎn)生邊界截?cái)嗾`差問題；虛粒子法通過在邊界外生成虛粒子解決了邊界截?cái)嗾`差問題，但在處理復(fù)雜邊界問題時，由于外部的虛粒子的生成較困難，且分布不均，計(jì)算精度受到影響，出現(xiàn)粒子飛散的情況。針對2種方法存在的問題，提出一種基于速度修正的固壁邊界處理方法，無需求解邊界力或在邊界外生成虛粒子，直接利用動量方程和計(jì)算速度耗損量求出流體粒子碰到邊界后的反彈速度，大大降低了處理邊界的復(fù)雜程度，也克服了2種方法在邊界拐角處粒子不均勻采樣而導(dǎo)致的算法不穩(wěn)定的問題。模擬仿真驗(yàn)證了該方法在穩(wěn)定性、計(jì)算效率方面均較傳統(tǒng)邊界力法和虛粒子法更好；隨著粒子數(shù)的增加，該方法耗時更少、計(jì)算效率更高，對復(fù)雜場景的模擬效果更好。

流體模擬；固壁邊界；邊界力法；虛粒子法；速度修正

光滑粒子流體動力學(xué)(smooth particle hydrodynamics, SPH)方法是流體模擬，尤其是液體模擬的主要方法。流體模擬的固壁邊界處理方法主要包括邊界力法和虛粒子法2種。

邊界力法應(yīng)用一層邊界粒子來表示固壁邊界，通過該層邊界粒子對鄰近邊界的流體粒子施加作用力，防止流體粒子穿透固壁邊界。其關(guān)鍵是壁面作用力的施加模型。

1994年，MONAGHAN[1]最先提出通過施加邊界力的方式進(jìn)行固壁邊界處理，并提出蘭納-瓊斯(Lennard-Jones，L-J)勢函數(shù)作為作用力的施加模型，但該模型中存在許多未知參數(shù)，需要根據(jù)具體問題進(jìn)行調(diào)整，調(diào)整不當(dāng)可能會導(dǎo)致流體粒子穿透邊界或粒子飛亂等問題。2004年，MüLLER等[2]基于L-J勢函數(shù)實(shí)現(xiàn)了在微可壓縮情況下，將SPH方法用于模擬流體與可變形固體的相互作用，并且在作用力較大時，流體粒子不會穿透固體邊界，但需要更小的時間步長來確保模擬的穩(wěn)定性。2007年，HARADA等[3]針對邊界粒子不足的問題，提出一種邊界粒子密度權(quán)重函數(shù)，將邊界粒子的密度加入到流體密度的計(jì)算中，從而減小了因邊界粒子的不足造成的計(jì)算誤差。2009年，BECKER等[4]以直接力的方法處理固壁邊界，同時采用預(yù)測-校正技術(shù)計(jì)算粒子的速度和位置。上述方法僅適合規(guī)則固體邊界的情況，且在計(jì)算鄰近邊界粒子的屬性時由于邊界外測粒子的缺失導(dǎo)致計(jì)算會產(chǎn)生截?cái)嗾`差。

2009年，MONAGHAN和KAJTAR[5]采用一種簡化的邊界力計(jì)算公式克服了L-J勢函數(shù)方法計(jì)算邊界力參數(shù)較多的不足，但力的大小和作用范圍仍難以控制。2012年，LIU等[6]提出了一種新型邊界力模型，并將邊界粒子與虛粒子結(jié)合施加邊界條件，但若單獨(dú)使用一層邊界粒子表示固壁邊界，仍會出現(xiàn)流體粒子穿透邊界的問題。邊界力法可以方便地表征復(fù)雜的邊界，且實(shí)現(xiàn)方法相對簡單，易于編程。然而，主要缺點(diǎn)為邊界力的范圍和強(qiáng)度難以控制，若邊界力太大，則會導(dǎo)致邊界附近粒子屬性的數(shù)值計(jì)算波動較大造成粒子飛亂；反之則將導(dǎo)致流體粒子穿透邊界。

虛粒子法應(yīng)用多層虛粒子表征固壁邊界，通過插值或鏡像方式將流體粒子的部分屬性(質(zhì)量、密度、壓強(qiáng)等)賦予虛粒子，虛粒子與流體粒子通過流體力學(xué)控制方程相互作用。虛粒子法的關(guān)鍵是其屬性值的獲取及與流體粒子的相互作用。

1996年，RANDLES和LIBERSKY[7]提出通過鏡像方式生成虛粒子進(jìn)行邊界處理，如果流體粒子到邊界的距離小于其支持域半徑，那么在固壁外就會自動生成鏡像虛粒子，流體粒子和鏡像虛粒子的速度大小相等但方向相反，其他物理量相同。1997年，MORRIS和MONAGHAN[8]使用插值虛粒子法模擬SPH流體與固壁邊界的相互作用。2002年，LIU等[9]提出在固壁邊界處布置2種虛粒子，第一種是固定在邊界上，且對靠近邊界的流體粒子施加一個排斥力，從而防止流體粒子穿透邊界；第二種是由鄰近邊界的流體粒子對稱產(chǎn)生的，并賦予其相應(yīng)物理量，從而使得邊界附近的流體粒子在進(jìn)行積分計(jì)算時不會產(chǎn)生截?cái)嗾`差。2006年，HU和ADAMS[10]使用鏡像虛粒子法模擬流體在規(guī)則固體邊界條件下的運(yùn)動，虛粒子設(shè)置在流體粒子作用域與邊界重疊的部分，并且虛粒子與流體粒子具有相同的密度、質(zhì)量、壓力和黏度，將虛粒子看作是流體粒子進(jìn)行處理。上述方法只適用于規(guī)則固壁直邊界，對于復(fù)雜邊界情況，虛粒子的生成過于復(fù)雜，計(jì)算量較大。

2012年，SCHECHTER和BRIDSON[11]提出了一種新鏡像粒子法來模擬流體的自由表面和進(jìn)行邊界處理。該方法在液體和固體之間以及液體和空氣之間產(chǎn)生鏡像粒子，并使用液體粒子通過插值運(yùn)算推導(dǎo)出鏡像粒子的屬性值。該方法解決了存在于流體自由表面附近的非自然張力視覺形變問題和在固壁邊界處無法真實(shí)模擬的流體的黏附現(xiàn)象。2015年，劉虎等[12]提出了應(yīng)用一系列的虛粒子來表征固壁邊界的處理方法，當(dāng)求解方程式時，讓虛粒子具有流體粒子的部分屬性(密度、質(zhì)量、壓力等)，求解其連續(xù)方程和狀態(tài)方程，并將邊界問題耦合到方程的求解過程中，無需提出顯式的邊界條件。該方法的優(yōu)點(diǎn)是無需對邊界進(jìn)行特殊的處理，缺點(diǎn)是通過求解連續(xù)方程和狀態(tài)方程得到的虛擬粒子的壓力，與真實(shí)壓力之間存在誤差，導(dǎo)致模擬的不穩(wěn)定，出現(xiàn)粒子飛散的情況。雖然虛粒子法的守恒性較好，精度也較高，但在生成自由表面的虛粒子時過于復(fù)雜，計(jì)算時間大大增加，固壁邊界形狀復(fù)雜時更加如此。

總體而言，邊界力法適合處理復(fù)雜的固體邊界，然而模型的參數(shù)較多，力的大小難以調(diào)控，并且在計(jì)算過程中會產(chǎn)生邊界截?cái)嗾`差問題；虛粒子法能夠很好地解決邊界截?cái)嗾`差，但在處理復(fù)雜邊界問題時，邊界外部的虛粒子的生成較困難，計(jì)算精度和穩(wěn)定性難以保證，會出現(xiàn)粒子飛亂甚至程序運(yùn)行崩潰的情況。

針對上述問題，本文提出一種基于速度修正的固壁邊界處理方法，無需在邊界外生成虛粒子或求解邊界力，直接利用動量方程和計(jì)算速度耗損量求出流體粒子遇到邊界后的反彈速度，大大降低了處理邊界的復(fù)雜程度，也克服了虛粒子法和邊界力法在拐角處粒子不均勻采樣而導(dǎo)致的算法不穩(wěn)定的問題。模擬仿真驗(yàn)證了本文方法在穩(wěn)定性、計(jì)算效率方面均比傳統(tǒng)方法更好；隨著粒子數(shù)的增加，本文方法的耗時增長慢于傳統(tǒng)方法，計(jì)算效率高的優(yōu)勢更明顯；對復(fù)雜場景的模擬效果更好。

1 SPH方法簡介

1.1 SPH基本公式

在SPH方法中，函數(shù)()定義為

其中，為光滑函數(shù)的影響和支持范圍即光滑長度；則稱為光滑核函數(shù)(smoothing kernel function)或光滑函數(shù)(smoothing function)，簡稱為核函數(shù)(kernel)。

本文采用三次樣條的核函數(shù)

其中，=/；α為函數(shù)在維空間中的歸一化系數(shù)，其在一維、二維和三維空間中的取值分別為1/，15/(7π2)，3/(2π3)。

通過將式(3)離散化，用光滑函數(shù)的支持域內(nèi)所有粒子疊加求和的形式表示任一點(diǎn)處的粒子近似式為

其中，為粒子的支持域內(nèi)的粒子總數(shù)；m，，分別為支持域范圍內(nèi)粒子的質(zhì)量、密度和位置；(–,)為粒子對粒子產(chǎn)生影響的光滑函數(shù)，且該函數(shù)與光滑長度有關(guān)。

1.2 控制方程

本文用到的是拉格朗日形式的Navier-Stokes方程(簡稱N-S方程)

本文采用SPH方法的計(jì)算過程流程如圖1所示。

圖1 SPH方法的流程圖

使用SPH方法求解N-S方程，對于任一粒子，其密度為

求得粒子密度后，還需要計(jì)算粒子壓強(qiáng)，進(jìn)而求得粒子的壓力。

利用泰特方程求壓強(qiáng)，即

其中，為常數(shù)，一般取=7；為參數(shù)，用于限制密度的最大改變量；0為參照密度。

SPH方法中求解粒子所受的壓力和黏性力的表達(dá)式分別為

流體粒子的速度、位置等信息是根據(jù)其所受合力計(jì)算所得，利用式(10)和式(11)，將N-S方程中的動量守恒方程式(7)轉(zhuǎn)化為

2 基于速度修正的固壁邊界處理方法

文獻(xiàn)[1]最早提出應(yīng)用邊界力的方式施加流固邊界條件，并提出一種基于距離的L-J勢函數(shù)的方法施加作用力

其中，系數(shù)1一般取12，2一般取4；的取值一般為流體最大速度的平方量級；0為光滑核半徑；r為流體粒子到固體邊界上位置點(diǎn)的距離；為流體粒子位置到固體邊界上位置點(diǎn)之間的向量。

文獻(xiàn)[7]在處理直壁邊界問題時，能夠很好地解決邊界截?cái)嗾`差，但在處理圓弧表面、傾斜壁面等邊界問題時，邊界外部的鏡像粒子的生成較困難，并且鏡像粒子分布不均，計(jì)算精度也會受影響，出現(xiàn)粒子躍出甚至程序崩潰的情況。

目前，已有的邊界處理方法都很難兼顧計(jì)算精度、邊界復(fù)雜度以及計(jì)算效率。本文提出一種新的基于速度修正的固壁邊界處理方法：先在邊界上設(shè)置一個阻尼區(qū)，當(dāng)流體粒子運(yùn)動到阻尼區(qū)內(nèi)時，且流體粒子速度與邊界法向夾角大于π/2時，將該區(qū)域分為4個小區(qū)域，并分別對流體速度作不同處理：將流體粒子速度從，坐標(biāo)系轉(zhuǎn)化為，坐標(biāo)系表示(其中，分別為邊界的法向與切向)；然后用流體粒子與相關(guān)的邊界粒子的相對位置計(jì)算速度耗損量；之后將不同區(qū)域的流體粒子反彈到相應(yīng)區(qū)域；最后將處理后的流體粒子速度轉(zhuǎn)回，坐標(biāo)系。具體步驟如下：

(1) 添加阻尼區(qū)。以邊界斜坡為例(圖2)：在斜坡邊界(黑色線條)上方添加阻尼區(qū)邊界(藍(lán)色線條)，阻尼區(qū)的寬度為，其數(shù)值與粒子的支持域長度(即核函數(shù)的光滑長度)相同。

圖2 邊界阻尼區(qū)示意圖

(2) 流體粒子到邊界的距離。在邊界上等距分布粒子，邊界粒子間距與初始流體粒子間距有關(guān)。流體粒子到邊界的距離定義為：流體粒子到邊界粒子的2個最短距離的平均值作為該流體粒子到邊界的距離，如圖3所示。

圖3 流體粒子到邊界粒子的距離

(3) 劃分區(qū)域。當(dāng)流體粒子在阻尼區(qū)內(nèi)，且其速度方向與邊界法向夾角大于π/2時，將該區(qū)域分為①、②、③、④4個小區(qū)域，如圖4所示，其中，①~④為第1~4個小區(qū)域。

依本文方法，區(qū)域①中的粒子將反彈到區(qū)域⑧中；區(qū)域②中的粒子反彈到區(qū)域⑦中；區(qū)域③的粒子反彈到區(qū)域⑥中；區(qū)域④的粒子反彈到區(qū)域⑤中。在流體粒子碰到邊界反彈的過程中，其反彈速度會受到周圍流體粒子的影響，故需要計(jì)算速度的耗損量。

(4) 計(jì)算速度耗損量。在流體粒子周邊的邊界粒子數(shù)為，記為1,2,···,，如圖5所示。

圖5 流體粒子與邊界粒子的夾角

由圖5得

其中，為流體粒子到邊界的距離(式(14))，r為粒子到粒子α的距離(=1,2,···)。

計(jì)算粒子速度耗損量為

當(dāng)粒子到達(dá)邊界時會發(fā)生反彈，各種復(fù)雜邊界均可用類似的方法進(jìn)行處理。

將本文方法與傳統(tǒng)邊界力法和虛粒子法進(jìn)行對比，可以發(fā)現(xiàn)：①在穩(wěn)定性方面，邊界力法是通過基于距離的L-J函數(shù)來計(jì)算鄰近邊界的流體粒子受到的排斥力，當(dāng)流體粒子運(yùn)動到拐角處時，因受到來自相同距離邊界粒子施加的力的影響導(dǎo)致粒子的反彈速度的計(jì)算不穩(wěn)定；虛粒子法在處理復(fù)雜邊界問題時，邊界外部的鏡像虛粒子的生成較困難，且分布不均勻，計(jì)算精度受影響。本文方法可對邊界的形狀做檢測，通過計(jì)算速度耗損量對速度進(jìn)行修正，可以更好地控制粒子在拐角處運(yùn)動的穩(wěn)定性，克服了傳統(tǒng)方法在邊界拐角處粒子不均勻采樣而導(dǎo)致的算法不穩(wěn)定的問題。②在計(jì)算效率方面，與邊界力法相比，本文方法無需使用較多參數(shù)且通過較復(fù)雜的邊界力模型來計(jì)算力，而是直接對速度進(jìn)行修正，通過計(jì)算速度耗損量求出流體粒子遇到邊界后的反彈速度，減小了計(jì)算量；與虛粒子法相比，本文方法無需在邊界外部生成虛粒子，降低了處理邊界的復(fù)雜程度，尤其在處理復(fù)雜邊界問題時，虛粒子的生成較困難，且在計(jì)算流體粒子的速度及位置信息時，虛粒子也參與其中，使得計(jì)算量大大增加。因此，本文方法在穩(wěn)定性和計(jì)算效率方面均要優(yōu)于傳統(tǒng)的邊界力法和虛粒子法。

3 實(shí) 驗(yàn)

3.1 二維潰壩模擬

以二維潰壩實(shí)驗(yàn)?zāi)M斜坡邊界和弧形邊界場景，驗(yàn)證本文方法在穩(wěn)定性和計(jì)算效率方面的有效性。相關(guān)參數(shù)：流體粒子的初始密度為1 000 kg/m3，粒子初始間距為0.002 92 m，時間步長為10×10–3s，粒子的支持域長度=0.006 m。

(a) 邊界力法(b) 虛粒子法(c) 本文方法

由圖6可知，圖6(a)為邊界力法，穩(wěn)定性略差，粒子撞到邊界后，上爬升過程中出現(xiàn)輕微的粒子飛散現(xiàn)象；圖6(b)為虛粒子法，穩(wěn)定性差，粒子撞到邊界后，上爬升過程中粒子飛散情況明顯；圖6(c)為本文方法，在粒子撞到右側(cè)邊界后，爬升時的穩(wěn)定性好，沒有出現(xiàn)明顯的粒子飛散情況，與實(shí)際中水撞到固壁邊界后爬升的效果更為接近。

上述3種方法模擬實(shí)驗(yàn)進(jìn)行流體粒子總數(shù)分別為2 500個、5 000個和7 500個粒子，每組實(shí)驗(yàn)進(jìn)行5次取運(yùn)行到2 000幀時的平均值作為運(yùn)行時間，模擬二維潰壩斜坡邊界運(yùn)行時間的結(jié)果見表1。

表1 3種方法模擬二維潰壩斜坡邊界的運(yùn)行時間

由表1和圖7可以看出，本文方法明顯比傳統(tǒng)的邊界力法和虛粒子法的計(jì)算效率高，并且隨著粒子數(shù)的增加，耗時增長也比其他2種方法慢，計(jì)算效率高的優(yōu)勢更明顯。

圖7 隨著粒子數(shù)的增加，邊界力法、虛粒子法和本文方法耗時增長對比

同表1模擬實(shí)驗(yàn)條件，3種方法模擬二維潰壩弧形邊界的運(yùn)行時間見表2。

由圖8和圖9及表2可知，本文方法在穩(wěn)定性和計(jì)算效率方面要優(yōu)于傳統(tǒng)邊界力法及虛粒子法，并且隨著粒子數(shù)的增加，本文方法的耗時增長慢于其他2種方法，計(jì)算效率更高。

(a) 邊界力法(b) 虛粒子法(c) 本文方法

表2 3種方法模擬二維潰壩弧形邊界的運(yùn)行時間

3.2 復(fù)雜場景模擬

圖10和圖11為3種方法模擬液體在漏斗和弓形水管的運(yùn)行場景以及邊界處理效果。

由圖10和圖11可知，圖10(c)為本文方法模擬漏斗的場景，上壁未出現(xiàn)明顯的粒子飛散情況；圖11(c)為模擬弓形水箱場景，流體流動的整體穩(wěn)定性較好，且粒子分布更均勻。說明本文方法在模擬復(fù)雜場景時的穩(wěn)定性更好，適用范圍更廣。

圖9 隨著粒子數(shù)的增加，邊界力法、虛粒子法和本文方法耗時增長對比

(a) 邊界力法(b) 虛粒子法(c) 本文方法

(a) 邊界力法(b) 虛粒子法(c) 本文方法

4 結(jié) 論

針對傳統(tǒng)的邊界力法和虛粒子法的不足，本文提出一種基于速度修正的固壁邊界處理方法，既不需要求解邊界力也無需在邊界外生成虛粒子，直接利用動量方程和計(jì)算速度耗損量求出流體粒子遇到邊界后的反彈速度，大大降低了處理邊界的復(fù)雜程度，也克服了邊界力法和虛粒子法在邊界拐角處粒子不均勻采樣而導(dǎo)致的算法不穩(wěn)定的問題。模擬仿真結(jié)果表明本文方法的穩(wěn)定性、計(jì)算效率均比傳統(tǒng)的2種方法更好，且隨著粒子數(shù)的增加，本文方法的耗時增長也比上述2種方法慢，計(jì)算效率高的優(yōu)勢更明顯；通過模擬一些復(fù)雜場景，本文方法未出現(xiàn)粒子飛散情況，穩(wěn)定性好，適用范圍更廣。

[1] MONAGHAN J J. Simulating free surface flows with SPH [J]. Journal of Computational Physics, 1994, 110(2): 399-406.

[2] MüLLER M, SCHIRM S, TESCHNER M, et al. Interaction of fluids with deformable solids [J]. Computer Animation and Virtual Worlds, 2004, 15(34): 159-171.

[3] HARADA T, KOSHIZUKA S, KAWAGUCHI Y. Smoothed particle hydrodynamics in complex shapes [J/OL].[2018-10-11].https://dl.acm.org/citation.cfm?id=2614375.

[4] BECKER M, TESSENDORF H, TESCHNER M. Direct forcing for Lagrangian rigid-fluid coupling [J]. IEEE Transactions on Visualization and Computer Graphics, 2009, 15(3): 493-503.

[5] MONAGHAN J J, KAJTAR J B. SPH particle boundary forces for arbitrary boundaries [J]. Computer Physics Communications, 2009, 180(10): 1811-1820.

[6] LIU M B, SHAO J R, CHANG J Z. On the treatment of solid boundary in smoothed particle hydrodynamics [J]. Science China (Technological Sciences), 2012, 55(1): 244-254.

[7] RANDLES P W, LIBERSKY L D. Smoothed particle hydrodynamics: Some recent improvements and applications [J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 375-408.

[8] MORRIS J P, MONAGHAN J J. A switch to reduce SPH viscosity [J]. Journal of Computational Physics, 1997, 136(1): 41-50.

[9] LIU M B, LIU G R, LAM K Y. Investigations into water mitigation using a meshless particle method [J]. Shock Waves, 2002, 12(3): 181-195.

[10] HU X Y, ADAMS N A. A multi-phase SPH method for macroscopic and mesoscopic flows [J]. Journal of Computational Physics, 2006, 213(2): 844-861.

[11] SCHECHTER H, BRIDSON R. Ghost SPH for animating water [J]. ACM Transactions on Graphics, 2012, 31(4): 1-8.

[12] 劉虎, 強(qiáng)洪夫, 陳福振, 等. 一種新型光滑粒子動力學(xué)固壁邊界施加模型[J]. 物理學(xué)報(bào), 2015, 64(9): 384-397.

Treatment of Solid Boundary Based on Velocity Correction

ZHU Xiao-lin, CHEN Wei

(School of Mathematics, Hefei University of Technology, Hefei Anhui 230009, China)

The study of solid boundary treatment method has been a difficult problem in fluid simulation. The common methods of solid boundary treatment are boundary force method and virtual particle method. The boundary force method prevents fluid particles from penetrating the boundary by exerting force on the fluid particles near the boundary, but this method has too many parameters and the force is difficult to adjust, and the boundary truncation error will occur in the calculation. The virtual particle method solves the problem of boundary truncation error by generating virtual particles outside the boundary. However, when dealing with complex boundary problems, the generation of virtual particles is difficult, and the calculation accuracy will be affected due to the uneven distribution of virtual particles, which leads to particles to disperse. To solve these problems, the paper presents a new method for the treatment of solid boundary based on velocity correction, which does not need to solve the boundary force or generate virtual particles outside the boundary. The momentum equation and the velocity consumption are used directly to calculate the rebound velocity of the fluid particles when they hit the boundary, which greatly reduces the complexity of the boundary treatment. It also overcomes the problem of the instability of the boundary force method and the virtual particle method caused by the uneven sampling of particles at the corner of the boundary. Simulation results show that the proposed method is more stable and more efficient than traditional methods above, and with the increase of the number of particles, the time consuming of this method is also slower than that of the two methods, and the advantage of high computational efficiency is more obvious and the simulation effect of complex scene is also better.

fluid simulation; solid boundary; boundary force method; virtual particle method; velocity correction

TP 391.9

10.11996/JG.j.2095-302X.2019040637

A

2095-302X(2019)04-0637-07

2019-01-30；

定稿日期：2019-03-20

國家自然科學(xué)基金項(xiàng)目(61272024)

朱曉臨(1964-)，男，安徽池州人，教授，博士，碩士生導(dǎo)師。主要研究方向?yàn)閿?shù)值逼近、計(jì)算機(jī)圖形學(xué)、CAGD、圖形圖像處理等。 E-mail：zxl_hfut@126.com