基于SPH的雨滴打擊不規(guī)則邊界的模擬方法

劉嘉瑞，楊 猛，吳佳澤，楊 剛

?

基于SPH的雨滴打擊不規(guī)則邊界的模擬方法

劉嘉瑞1，楊 猛1，吳佳澤2，楊 剛1

(1. 北京林業(yè)大學(xué)信息學(xué)院，北京 100083；2. 深圳市六聯(lián)科技有限公司，廣東 深圳 518109)

為實現(xiàn)對雨滴打擊樹枝等不規(guī)則邊界過程的模擬，研究了流體粒子在網(wǎng)格表示的固體邊界處的受力情況，提出了一種不需要粒子采樣的邊界受力模擬方法。采用高斯積分法則對網(wǎng)格模型的三角面片進(jìn)行積分，并就此對固液邊界的粒子的密度進(jìn)行修正，以積分的方法對固體邊界處的壓力、粘性力等參數(shù)進(jìn)行計算，從而保證邊界粒子受力的連續(xù)性。同時，還提出了一種吸引力模型，用來控制粒子在沿著物體表面滑落時的運動。實驗結(jié)果表明，該方法在模擬水滴鋪展、收縮、沿著邊界流動等現(xiàn)象時達(dá)到了較為真實的效果。

SPH；高斯積分法則；邊界條件；流體

流體的模擬一般分為網(wǎng)格法和粒子法兩大類，其中粒子法由于其計算簡單，易于并行化的特點，越來越受到大家的關(guān)注。在近年的流體領(lǐng)域，各種流體模擬方法層出不窮，并都取得了不錯的成果。然而，在目前的研究中，流體的模擬大多側(cè)重于對流體本身算法的改進(jìn)。在模擬大規(guī)模的場景中，邊界力的模擬方法對效果影響不大，而研究固液耦合的情況時，由于剛體一般體積較小，學(xué)者們側(cè)重剛體和流體的相互作用，因此邊界條件的領(lǐng)域一直沒有得到足夠的發(fā)展。

在微觀流體，如雨滴等的模擬中，由于液滴本身較小，極易受外力影響，因此其形態(tài)變化主要受邊界影響。復(fù)雜邊界的表示形式主要有網(wǎng)格和粒子兩種，粒子一般為算法采樣得到，并且有一次性采樣和動態(tài)生成采樣粒子兩種方法，其優(yōu)點是便于修正粒子，保證了邊界粒子力的連續(xù)性，但是采樣算法較為復(fù)雜，并且難以在可變形的邊界方向進(jìn)一步發(fā)展。

本文受MüLLER等[1]提出的三角形積分方法的啟發(fā)，使用高斯光滑函數(shù)實現(xiàn)三角形面積上力的積分。但不同于其直接施加作用力的方法，本文參考了一些基于邊界采樣粒子的受力模型，并在此基礎(chǔ)上提出了一種新的對粒子各參數(shù)進(jìn)行修正的方法。

同時，參考AKINCI等[2]的方法，在修正邊界的同時提出了一個吸引力模型，使得流體沿著模型流動，達(dá)到了更真實的效果，如圖1所示。

圖1 雨滴打擊不規(guī)則邊界效果圖

1 相關(guān)工作

1.1 網(wǎng)格法和粒子法

關(guān)于流體的模擬，主流的兩種方法為網(wǎng)格法和粒子法；數(shù)值計算網(wǎng)格大體可以分為不隨流體運動而改變形狀的歐拉網(wǎng)格與隨著流體的運動而改變形狀的拉格朗日網(wǎng)格兩種[3]。文獻(xiàn)[4-5]在1977年分別提出光滑粒子法(smoothed particle hydrodynamics，SPH)，但最初多用于求解天體物理，文獻(xiàn)[6]提出用SPH模擬流體的方法，避免了網(wǎng)格法在計算過程中會出現(xiàn)網(wǎng)格纏結(jié)和扭曲等問題，后人的研究大多在此方法基礎(chǔ)上進(jìn)行改進(jìn)。

1.2 邊界處理

目前有很多學(xué)者做過關(guān)于流體和固體邊界的研究：MONAGHAN[7]第一次用一些線性的粒子表示運動邊界，從而計算流體在邊界的受力。BECKER等[8]提出了用粒子表示整個固體模型，然后接觸粒子直接給流體施加力的方法，較為簡單容易實施；LIU等[9]使用虛粒子(virtual particles)提供斥力，但這事實上仍然是一種直接施加力的方法。

AKINCI等[10]于2012年提出了用粒子表示的模型上修正密度、壓強(qiáng)、粘性力等的方法，獲得了較好的效果。2013年，又提出了一個較好的計算張力的模型，并提出了施加額外粘附力的方法以保證流體在固體邊界的表面運動[2]。

SCHECHTER和BRIDSON[11]提出了魂粒子(ghost particle)的方法，在粒子周圍的空氣和剛體表面鋪上虛擬的粒子，用來修正邊界粒子的密度，并模擬剛體和流體力的作用，獲得了較為真實的效果。但這些方法都需要將三角面片轉(zhuǎn)化為粒子的表示形式，不適合大型規(guī)模的剛體交互。

文獻(xiàn)[1]提出了用插值方法表示三角面片的方法，獲得了較好的交互效果，在此基礎(chǔ)上，YANG等[12]在模擬血液和血管的交互時使用了代理粒子(proxy particle)，通過動態(tài)的生成代理粒子，實現(xiàn)了血管和血液的交互。但是由于粒子邊界連續(xù)性問題并沒有被解決，這兩種方法均在真實性方面有所欠缺。

本文解決了SPH形式的流體在三角網(wǎng)格表示邊界的情況下的參數(shù)修正和受力模擬問題，保證了粒子在邊界的受力的連續(xù)性。

2 本文算法

2.1 基本方法

2.1.1 Navier-Stokes方程

在流體力學(xué)中，常用式(1)的經(jīng)典運動模型——納維斯托克斯方程(Navier-Stokes equations)去描述粒子的運動，即

2.1.2 SPH方法

SPH的本質(zhì)是一種插值[14]。SPH采用離散的粒子描述連續(xù)分布的流體，每個粒子攜帶了流體各種性質(zhì)[5-6]。任一宏觀變量(如密度、壓力溫度等)能方便地借助于一組無序點上的值表示成積分插值，其形式為

由于粒子的位置、溫度和壓力等特性是一系列離散化的值，所以常被寫成離散化的形式

通過周圍粒子的質(zhì)量求和的方式得到粒子的密度為

壓力在流體的計算中使得流體的密度趨向于均勻，其大小一般由粒子和周圍粒子的密度決定。壓力場可看做壓強(qiáng)的梯度。由式(5)帶入可得壓強(qiáng)的計算公式，為了實現(xiàn)力的對稱性，一般取壓強(qiáng)和周圍粒子的平均值，求粒子所受壓力為

粘性力部分負(fù)責(zé)阻滯粒子的運動，使得粒子的運動受到周圍粒子的約束，避免粒子完全散開。邊界處的粘性力使得流體沿著表面運動，與壓力推導(dǎo)相似的推導(dǎo)方法可得

對于大型流體的模擬，如海浪、瀑布等，往往忽略表面張力對粒子的作用；然而在考慮小的液滴，如雨滴時，張力承擔(dān)了不可或缺的作用。表面張力的微觀機(jī)理是分子間引力，力的方向的計算為

由此得到平滑后的張力為

其中，為張力的系數(shù)，在實驗中得到。

2.2 邊界處理

在現(xiàn)有的算法中，網(wǎng)格表示的邊界傾向于直接用數(shù)學(xué)的方法建立出相對作用力的模型，如圖2(a)所示。其計算簡單，方便計算軟體邊界的變形力，但是由于本身方法的局限，對流體本身的邊界模擬效果并不突出。固體采樣的方法如圖2(b)所示，使得流體邊界的受力得到了較好的模擬，但是采樣的過程增加了計算的成本，并且不適合模型較大的情況。高斯三角形積分方法如圖2(c)所示，其結(jié)合了前兩者的優(yōu)點，既避免了復(fù)雜的采樣過程，又保證了流體邊界受力的連續(xù)性與合理性。

當(dāng)粒子和剛體的距離小于一定閾值時，則認(rèn)為粒子和剛體發(fā)生了交互，此時需要對邊界的密度壓力粘性力等進(jìn)行修正。

本文通過高斯三角形七點積分法則(seven point rule)計算出一定范圍內(nèi)的三角形對粒子的力貢獻(xiàn)，使得在不需要用特殊方法對邊界進(jìn)行粒子采樣的情況下對邊界的密度和受力進(jìn)行修正。假設(shè)邊界密度和流體初始密度相同，從而通過積分點計算得到固體邊界對粒子密度的貢獻(xiàn)。同理，通過計算三角形上積分點對粒子的力，近似積分得出固體邊界對粒子的作用力。

圖2 邊界受力方法比較示意圖

2.2.1 高斯積分法則

高斯積分是一種在[–1, 1]區(qū)域上的近似積分方法，點高斯積分就是在定義域上取個采樣點，計算其在每一點的函數(shù)值，并且以一定權(quán)重相加，作為最后的積分結(jié)果，不在[–1, 1]范圍內(nèi)的函數(shù)往往要映射到這個范圍上計算，即

三角形的高斯積分方法是對高斯積分的一個擴(kuò)展，其將三角形面積上的積分用幾個積分點的加權(quán)和表示，從而避免了復(fù)雜的計算(圖3)。假設(shè)所有的三角形的坐標(biāo)都是(0, 0)，(1, 0)，(0, 1)這3個點，普通的三角形都映射到這個單位三角形上，并在這個三角形上取個采樣點進(jìn)行計算，最后加權(quán)求和得到近似的積分。其形式為

其中，為普通三角形的面積。本文延續(xù)了Muller采用的七點積分方法[1]，并采用了ghost SPH中修正邊界粒子參數(shù)的方法，將固體邊界貢獻(xiàn)的密度和壓強(qiáng)等通過這種積分方式計算得到，然后加到流體粒子上進(jìn)行修正。

2.2.2 密度修正

由密度計算公式得粒子的密度等于周圍粒子的質(zhì)量的加權(quán)平均。對于周圍存在三角形的粒子，計算每個三角形的積分點，并得出這個三角形的質(zhì)量貢獻(xiàn)。

圖3 高斯七點積分方法示意圖

因此，固體對第個粒子貢獻(xiàn)的密度為