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

?

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

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

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

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

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

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

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

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

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

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

1 相關(guān)工作

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

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

1.2 邊界處理

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

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

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

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

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

2 本文算法

2.1 基本方法

2.1.1 Navier-Stokes方程

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

2.1.2 SPH方法

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

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

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

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

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

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

由此得到平滑后的張力為

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

2.2 邊界處理

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

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

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

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

2.2.1 高斯積分法則

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

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

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

2.2.2 密度修正

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

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

因此，固體對第個粒子貢獻的密度為