朱曉臨，汪歡歡，殷競(jìng)存

?

SPH流固耦合模擬自適應(yīng)采樣固壁虛粒子邊界處理方法

朱曉臨，汪歡歡，殷競(jìng)存

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

固壁虛粒子邊界處理方法是流體模擬中一種主要邊界處理方法，但其不能確保流體粒子不穿透固體邊界，并且計(jì)算量較大。為防止流體粒子穿透邊界，在邊界附近設(shè)置一個(gè)阻尼區(qū)，阻尼區(qū)內(nèi)的流體粒子被邊界施加一個(gè)彈性力和一個(gè)和流體粒子運(yùn)動(dòng)速度方向相反的阻尼力，使得邊界附近流體粒子更加穩(wěn)定。為減少計(jì)算量，提出兩種邊界粒子自適應(yīng)采樣法：一種是依據(jù)邊界周圍粒子數(shù)目的不同，邊界粒子自適應(yīng)地采樣質(zhì)量不同的大小粒子；另一種是依據(jù)邊界周圍粒子數(shù)目的不同，邊界粒子自適應(yīng)的采樣不同層數(shù)的相同質(zhì)量粒子。與傳統(tǒng)的固體邊界粒子采樣方法相比，該方法減少了邊界粒子數(shù)目，加快了模擬速度，節(jié)省了計(jì)算機(jī)內(nèi)存，基于GPU加速技術(shù)實(shí)現(xiàn)的三維流體模擬，能夠進(jìn)行實(shí)時(shí)交互。

固壁虛粒子；阻尼區(qū)；彈性力；阻尼力；邊界粒子自適應(yīng)采樣；GPU加速技術(shù)

流體和固體的交互在日常生活中是普遍存在的，例如石頭掉進(jìn)水里、衣服浮在水面上、冰塊的融化、海綿吸水等。流固耦合動(dòng)力學(xué)涉及到流體動(dòng)力學(xué)和固體動(dòng)力學(xué)，對(duì)于流固耦合問(wèn)題中的固體研究總體上可以分為剛體和可變形的物體；流體主要分為氣體與液體。從總體來(lái)看，流固耦合問(wèn)題按照耦合機(jī)理可以分為兩類：第一類是指流體與固體的相互作用只發(fā)生在流體與固體的交界面上，即固體不具有滲水性，如石頭和水的作用；第二類是指流體區(qū)域與固體區(qū)域部分或全部重疊在一起，如滲流問(wèn)題。在計(jì)算機(jī)圖形學(xué)領(lǐng)域，流固耦合問(wèn)題的模擬主要是基于物理的偏微分方程求解方法的數(shù)值模擬。主要的流固耦合數(shù)值模擬方法有有限元法、邊界元法、光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)方法等。然而，有限元法需要全局離散，導(dǎo)致問(wèn)題的自由度和原始信息量大；邊界元法只在邊界上離散問(wèn)題域，降低了問(wèn)題的復(fù)雜度；有限元法和邊界元法都是歐拉網(wǎng)格方法，難以模擬具有較大變形的流固耦合自然現(xiàn)象。SPH方法是一種拉格朗日無(wú)網(wǎng)格插值方法，能方便地處理具有大變形流固耦合問(wèn)題。對(duì)于SPH方法進(jìn)行流固耦合模擬，邊界處理是研究的重點(diǎn)和難點(diǎn)。對(duì)于SPH流固耦合的邊界處理，目前主要有4類方法：基于距離的排斥力方法、固壁虛粒子法、虛粒子法、預(yù)測(cè)矯正法前3類方法是基于固體分別給予流體不同的力而實(shí)現(xiàn)耦合效應(yīng)；第4類方法是計(jì)算機(jī)圖形學(xué)領(lǐng)域中常用的正法?；诰嚯x的排斥力方法計(jì)算的固體流體之間的力不是物理力，導(dǎo)致模擬時(shí)粒子堆積問(wèn)題；固壁虛粒子法由于將固體看成由流體粒子組成，能避免了粒子堆積問(wèn)題，但計(jì)算量較大，并且不能確保流體粒子不穿透固體邊界；虛粒子法由于需要計(jì)算固體邊界處的法向量，并且動(dòng)態(tài)的生成邊界粒子，計(jì)算量很大，對(duì)于沒有確定形狀的邊界很難計(jì)算邊界法向量，也不能保證粒子不穿透邊界；預(yù)測(cè)矯正法是對(duì)每個(gè)流體粒子進(jìn)行預(yù)測(cè)位置，對(duì)于即將發(fā)生穿透固體的流體粒子進(jìn)行位置以及速度矯正，流體與固體之間沒有相互作用力，沒有物理規(guī)律支持。目前，對(duì)于SPH流固耦合的邊界處理主要使用固壁虛粒子法。

1 相關(guān)工作

早期的流體模擬主要通過(guò)參數(shù)建模來(lái)模擬，文獻(xiàn)[1]采用對(duì)波浪進(jìn)行函數(shù)建模方法模擬波浪效果，但并非基于物理規(guī)律計(jì)算，效果不真實(shí)。FOSTER和METAXAS[2]首次通過(guò)對(duì)流體運(yùn)動(dòng)控制方程數(shù)值求解進(jìn)行流體模擬。STAM和FIUME[3]首次將SPH方法引入到計(jì)算機(jī)圖形學(xué)領(lǐng)域，模擬氣體與火焰效果。隨后DESBRUN和GASCUEL[4]應(yīng)用SPH方法模擬可變形物體；MULLER等[5]運(yùn)用SPH方法模擬流體。SPH方法在流體模擬領(lǐng)域得到快速發(fā)展，并進(jìn)一步發(fā)展到流固耦合問(wèn)題中。在SPH方法進(jìn)行流固耦合模擬時(shí)，研究的難點(diǎn)和重點(diǎn)在于固體邊界的處理。MONAGHAN等[6]最先提出通過(guò)固體邊界上的一系列位置點(diǎn)對(duì)靠近固體邊界的流體粒子施加一個(gè)排斥力以防止流體粒子穿透發(fā)生。后來(lái)很多學(xué)者對(duì)此進(jìn)行研究，并提出了一系列的排斥力公式。HU和ADAMS[7]基于虛粒子方法實(shí)現(xiàn)了固體邊界是平面條件下SPH流固耦合的穩(wěn)定模擬；SOLENTHALER等[8]對(duì)固體邊界進(jìn)行粒子采樣，通過(guò)提出一種統(tǒng)一化的粒子模型，即固壁虛粒子法，在流體粒子的物理屬性值計(jì)算過(guò)程中，將固體邊界粒子看作流體粒子，參與流體計(jì)算，以保證流固交界面附近流體粒子的密度和壓強(qiáng)分布均勻，模擬出液體固化和固體融化的流固耦合效果。在固壁虛粒子方法的基礎(chǔ)上，SUN和HAN[9]對(duì)每個(gè)流體粒子實(shí)行偵察機(jī)制，將即將發(fā)生穿透剛體的流體粒子速度改為為剛體速度，流體粒子位置改為邊界點(diǎn)，但是局部效果不真實(shí)。AKINCI等[10]通過(guò)考慮固體邊界粒子對(duì)流體粒子屬性計(jì)算的相對(duì)貢獻(xiàn)，有效解決了邊界粒子分布不均勻?qū)е碌倪吔缌Ｗ佣逊e，并模擬了流體與剛性體狀、面狀及線狀物體的耦合。以上固壁虛粒子法都不能保證流體粒子不穿透固壁邊界，且計(jì)算量較大。

本文針對(duì)SPH流固耦合邊界處理最常用的固壁虛粒子方法進(jìn)行改進(jìn)，提出SPH流固耦合模擬自適應(yīng)采樣固壁虛粒子邊界處理方法，即為阻尼區(qū)結(jié)合自適應(yīng)邊界粒子采樣邊界處理方法。

2 SPH流固耦合建模

2.1 SPH方法

SPH方法是一種積分插值的方法，將模擬區(qū)域離散成一系列的帶有物理屬性值的粒子。這些物理屬性值包括粒子的速度、位置、所受的力、壓強(qiáng)和溫度等。模擬區(qū)域中位置處的物理屬性值()可以通過(guò)包含位置的支持域中的粒子屬性值進(jìn)行插值計(jì)算，即

其中，m為第粒子的質(zhì)量；A為第個(gè)粒子的屬性值；為第個(gè)粒子的密度；函數(shù)(-x,)為光滑核函數(shù)；x為第個(gè)粒子的位置；為光滑核函數(shù)的支持域半徑。

2.2 SPH流體固體建模

因?yàn)榱黧w的運(yùn)動(dòng)規(guī)律遵守Navier-Stokes方程，所以流體的運(yùn)動(dòng)模擬就是求解Navier-Stokes方程，即

因?yàn)楸疚牡牧黧w模擬中含有自由面，而自由面附近處流體粒子周圍的流體粒子數(shù)目偏少，若采用式(5)進(jìn)行密度計(jì)算，將導(dǎo)致自由面附近處流體粒子密度偏小，所以本文采用修正的密度計(jì)算公式，即

壓強(qiáng)計(jì)算公式采用不可壓縮性條件，即

因?yàn)榱黧w粒子對(duì)之間的壓力和粘性力滿足牛頓第二定律，所以本文采用對(duì)稱化的方式進(jìn)行壓力和粘性力計(jì)算，即

(1) 密度核函數(shù)

(2) 壓強(qiáng)核函數(shù)

(3) 粘性核函數(shù)

2.3 邊界粒子自適應(yīng)采樣

傳統(tǒng)固壁虛粒子方法的固體邊界粒子采樣通常是在模擬的初始時(shí)刻就給出所有邊界粒子的分布，粒子的層數(shù)通常是兩層粒子，并且邊界粒子數(shù)目以及分布不會(huì)隨著模擬的進(jìn)行而改變。在每一個(gè)模擬時(shí)間步，所有邊界粒子都要進(jìn)行全部物理屬性值的計(jì)算，計(jì)算量較大。隨著模擬的進(jìn)行，流體并不總是和固體有接觸，因此，沒有和流體區(qū)域接觸的固體邊界的計(jì)算是多余的。本文提出兩種邊界粒子自適應(yīng)采樣的方法：一種是依據(jù)邊界周圍粒子數(shù)目的不同，邊界粒子自適應(yīng)地采樣質(zhì)量不同的大小粒子，稱之為方法一；另一種是依據(jù)邊界周圍粒子數(shù)目的不同，邊界粒子自適應(yīng)的采樣不同層數(shù)的相同質(zhì)量粒子，稱之為方法二。

2.3.1 傳統(tǒng)的固壁虛粒子法

圖1 固壁虛粒子法邊界粒子采樣

2.3.2 邊界粒子自適應(yīng)采樣方法一

圖2 邊界粒子采樣一

2.3.3 邊界粒子自適應(yīng)采樣方法二

圖3 邊界粒子采樣二

2.4 邊界彈簧阻尼力模型

圖4中，空心圓表示的是流體粒子。固體邊界通過(guò)自適應(yīng)固體邊界粒子表示之后，固體粒子和流體粒子進(jìn)行統(tǒng)一的連續(xù)性方程和動(dòng)量方程計(jì)算，雖然解決了邊界截?cái)嗾`差問(wèn)題，但不能保證固體的非穿透性。本文在固體粒子和流體粒子進(jìn)行統(tǒng)一的連續(xù)性方程和動(dòng)量方程計(jì)算的基礎(chǔ)上，在固壁邊界附近設(shè)置一個(gè)阻尼區(qū)；當(dāng)流體粒子進(jìn)入事先設(shè)置的阻尼區(qū)之后，邊界會(huì)給流體粒子一個(gè)彈簧阻尼力，防止粒子穿透邊界。傳統(tǒng)的邊界力是關(guān)于距離的非線性力，此種類型的邊界力會(huì)產(chǎn)生邊界擾動(dòng)，邊界附近粒子不穩(wěn)定，且一旦流體粒子穿透邊界，就不再會(huì)進(jìn)入流體區(qū)域。本文通過(guò)對(duì)進(jìn)入阻尼區(qū)的流體粒子施加一個(gè)類似于彈簧的彈性力；彈性力是一種關(guān)于距離的線性力，所以理論上比傳統(tǒng)的非線性力穩(wěn)定?？拷吔绲牧黧w粒子速度可能很快，為了使得邊界處的流體粒子更加穩(wěn)定，通過(guò)施加一個(gè)和流體粒子速度相反的力，使得靠近邊界具有較快速度的流體粒子速度降下來(lái)，使得邊界附近的粒子更加的穩(wěn)定，即

圖4 彈簧阻尼模型

2.5 本文自適應(yīng)流固耦合模擬算法步驟

步驟1.流體粒子初始布局。

步驟2.自適應(yīng)邊界粒子采樣。

3 實(shí)驗(yàn)結(jié)果

本節(jié)通過(guò)三維潰壩模擬來(lái)驗(yàn)證本文提出的方法。為了更清晰、直觀地演示流體粒子穿透固體邊界的情形，采用三維潰壩模擬的二維剖面進(jìn)行演示，并與傳統(tǒng)的固壁虛粒子法進(jìn)行比較，包括模擬時(shí)間的比較。最后給出三維潰壩模擬效果的渲染圖，并與固壁虛粒子法進(jìn)行比較。

3.1 二維潰壩模擬

二維潰壩實(shí)驗(yàn)是在windows 7 平臺(tái)上進(jìn)行模擬，實(shí)驗(yàn)中的流體粒子數(shù)是1 200，固壁虛粒子的粒子數(shù)為200，實(shí)驗(yàn)開發(fā)環(huán)境為visual studio 2013，實(shí)驗(yàn)配置為Intel (R) Core (TM) i5-2400，3.10 GHz CPU，4 G內(nèi)存，NVIDIA Geforce GT 620顯卡。

圖5中紅色粒子表示的是流體粒子，藍(lán)色粒子表示的是固體邊界粒子。從圖5可以看出，固壁虛粒子邊界方法中的流體粒子穿透了固體邊界，而本文的兩種自適應(yīng)邊界粒子采樣方法沒有粒子穿透邊界。

從圖6中可以看出，本文的兩種方法所得到的壓力更加穩(wěn)定，過(guò)渡更加平滑(表1)，本文方法模擬的精度要高，模擬速度更快。

圖5 二維潰壩不同方法模擬圖

3.2 三維潰壩模擬

三維潰壩實(shí)驗(yàn)是在windows 7平臺(tái)上進(jìn)行模擬，實(shí)驗(yàn)中的流體粒子數(shù)是16 384，固壁虛粒子的粒子數(shù)為2 048，實(shí)驗(yàn)開發(fā)環(huán)境為visual studio 2013，以及cuda8.0進(jìn)行加速計(jì)算，最后的渲染部分本實(shí)驗(yàn)采用了povray進(jìn)行渲染，實(shí)驗(yàn)配置為Intel (R) Core (TM) i5-2400，3.10 GHz CPU、4 G內(nèi)存，NVIDIA Geforce GT 620顯卡。

從圖7可以看出，本文的兩種自適應(yīng)邊界粒子采樣方法的模擬結(jié)果和固壁虛粒子邊界法的模擬效果相當(dāng)，其比固壁虛粒子方法模擬速度更快(表2)。

圖6 二維潰壩不同方法模擬的壓力云圖

表1 二維潰壩不同方法模擬時(shí)間比較

圖7 三維潰壩不同方法模擬渲染圖

表2 三維潰壩不同方法模擬時(shí)間比較

4 結(jié) 論

本文提出一種改進(jìn)的固壁虛粒子邊界處理方法，通過(guò)在邊界附近設(shè)置一個(gè)阻尼區(qū)，分別采用兩種邊界粒子自適應(yīng)采樣方法，不但有效地防止了粒子穿透固體邊界，而且減少了計(jì)算量和內(nèi)存。利用GPU加速技術(shù)，通過(guò)統(tǒng)一化的粒子模型并運(yùn)用SPH方法實(shí)現(xiàn)了固體和流體的耦合模擬。

[1] PEACHEY D R. Modeling waves and surf [J]. ACM Computer Graphics, 1986, 20(4): 65-74.

[2] FOSTER N, METAXAS D. Relastic animation of liquids [J]. Graphical Models and Image Processing, 1996, 58(5): 471-483.

[3] STAM J, FIUME E. Depicting fire and other gaseous phenomena using diffusion processes [C]//Proceedings of 22ndAnnual Conference on Computer Graphics and Interactive Techiniques. New York: ACM Press, 1995: 129-136.

[4] DESBRUN M, GASCUEL M P. Smoothed particles: a new paradigm for animating highly deformable bodies [C]//Proceedings of Eurographics Workshop on Computer Animation and Simulation’ 96, Berlin: Springer, 1996: 61-76.

[5] MULLER M, CHARYPAR D, GROSS M. Particle-based fluid simulation for interactive applications [C]//Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Goslar: Eurographics Association, 2003: 154-159.

[6] MONAGHAN J J, THOMPSON M C, HOURIGAN K. Simulating free surface flow with SPH [J]. Journal of Computational Physics, 1994, 110(2): 399-406.

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

[8] SOLENTHALER B, SCHLAFLI J, PAJAROLA R. A unified particle model for fluid-solid interactions [J]. Computer Animation and Virtual Worlds, 2007, 18(1): 69-82.

[9] SUN H Q, HAN J Q. Particle-based realistic simulation of fluid-solid interaction [J]. Computer Animation and Virture World, 2010: 21(6): 589-595.

[10] AKINCI N, IHMSEN M, AKINCI G, et al. Versatile rigid-fluid coupling for incompressible SPH [J]. ACM Transactions on Graphics, 2012, 31(4): 62: 1-62: 8.

Adaptive Sampling Method for Solid Wall Virtual Particle Boundary Handling in SPH Fluid Solid Coupling Simulation

ZHU Xiaolin, WANG Huanhuan, YIN Jingcun

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

The method for solid wall virtual particle boundary handling is a main boundary handling method in fluid simulation, while it cannot ensure that no fluid penetrates the solid boundary, and requires a huge amount of computation. To prevent fluid particles from penetrating the boundary, in this paper a damp zone is set up near the boundary, and each fluid particle in this zone will be applied by the boundary with an elastic force and a damping force that is opposite to the moving velocity of the fluid particle. Therefore, the particles near the boundary become more stabilized. To reduce the amount of computation, this paper proposes two adaptive sampling methods for boundary particles based on different numbers of particles around the boundary: one is to adaptively sample particles of different masses; the other is to adaptively sample particles of the same mass in different layers. Compared with the traditional method for solid boundary particle sampling, these two adaptive sampling methods reduce the number of boundary particles, accelerate the simulation speed and save the computer memory. In addition, all of the 3D fluid simulations in this paper are based on the GPU acceleration technique, and can interact in real time.

solid wall virtual particle; damping zone; elastic force; damping force; adaptive sampling of boundary particles; GPU acceleration technique

TP 391

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

A

2095-302X(2018)03-0381-08

2017-01-14；

2017-05-11

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

朱曉臨(1964-)，男，安徽池州人，教授，博士。主要研究方向?yàn)閿?shù)值逼近、CAGD、圖形圖像處理。E-mail：zxl_hfut@126.com