朱曉臨，何紅虹，郭清偉，周韻若

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

流體模擬是計算機圖形學的一個重要分支；模擬固體入水是流體模擬的一個重要研究課題，有著廣泛地應用背景；例如，運動員跳水，水上飛機降落、水下航行器的投放等。其中，在固體入水的模擬過程中，流體與固體的相互作用，流體與固體交互時交界面的壓力處理、流體表面的大變形都是研究的重點和難點。

光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)方法在處理大變形和快速移動的界面具有獨特的優(yōu)勢。因此，SPH方法被廣泛應用于各種流-固耦合問題。但傳統(tǒng)SPH方法在模擬運動固體邊界時，由于固體邊界附近流體粒子支持域被截斷，流體粒子密度不連續(xù)，使得計算精度較差，會導致邊界處壓力振蕩等問題；而流體表面因為受到固體的沖擊力，流體液面較粗糙。

傳統(tǒng)SPH方法處理固體邊界的方法有邊界力法、虛粒子法和耦合力法。文獻[1]最早使用邊界力法處理固壁邊界，但該方法存在多個未知參數(shù)，需要根據(jù)具體問題進行調(diào)整；若調(diào)整不當，則可能導致流體粒子穿透固壁邊界和粒子混亂等問題。為了防止流體粒子穿透固體邊界，文獻[2]使用邊界力法對流體粒子施加較大的排斥力，并提供有限的邊界處理，給方程引入了剛度參數(shù)；排斥力只在流體粒子穿透固體邊界時起作用，粒子到邊界的距離隨時間變化，且粒子可能在邊界處出現(xiàn)異常加速，不符合事實。文獻[3]提出一種雙向的邊界力法，防止流體粒子穿透固體邊界，采用修正的SPH 壓力公式，避免了負壓；在流固發(fā)生碰撞時，對碰撞點處的相對速度施加一個邊界條件，但容易發(fā)生粒子堆積現(xiàn)象。為提高計算精度，解決邊界處流體粒子支持域被截斷的問題。文獻[4]采用虛粒子法處理固體邊界，利用固定虛粒子表征固壁邊界，除不發(fā)生位移外，虛粒子與流體粒子屬性完全相同，可出現(xiàn)少許流體粒子穿透固壁邊界的現(xiàn)象。文獻[5]在模擬水下爆炸等問題時，提出了一種鏡像生成虛粒子的方法，該方法守恒性較好、精度高，但動態(tài)生成虛粒子較復雜。文獻[6]利用虛粒子法模擬固壁邊界，直接將內(nèi)部流場壓力插值到虛粒子上，虛粒子通過壓力梯度對流場施加邊界力；但是，當流體中出現(xiàn)壓力動蕩時，插值得到的邊界虛粒子壓力值可能為負，虛粒子會對流體產(chǎn)生非正常的吸引力，導致流體粒子穿透固壁邊界。

針對邊界力法和虛粒子法的不足，文獻[7]提出耦合力法，將邊界力與虛粒子相結(jié)合施加邊界條件，防止流體粒子穿透固體邊界，提高了計算精度。該方法能夠獲得精確的光滑壓力場，適用于復雜或移動的固體邊界，但是單獨使用邊界粒子表征固壁邊界，可能會導致流體粒子穿透固體邊界。文獻[8]在模擬固體入水和浮體出水時，采用了耦合力法，對非均勻分布粒子進行密度校正和核梯度校正；但其采用高階SPH 格式，雖精度較高，但計算復雜，模擬速度較慢，實時性較差。

針對流體表面較粗糙的問題，文獻[9]采用各向異性核函數(shù)重構(gòu)更為光滑、平整的流體表面，真實地展現(xiàn)出流體表面的細節(jié)；但由于確定粒子各向異性核函數(shù)非常復雜，該方法重構(gòu)流體表面的速度較慢。針對傳統(tǒng)各向異性核函數(shù)構(gòu)造流體表面速度慢的缺陷，文獻[10]提出流體內(nèi)部粒子采用各向同性核函數(shù)，邊界粒子采用各向異性核函數(shù)，提高了流體表面繪制的速度，但模擬的逼真度不夠。文獻[11]根據(jù)流體粒子的特征值比值將粒子分為近表面粒子和內(nèi)部粒子，表面重構(gòu)時近表面粒子參與計算,內(nèi)部粒子根據(jù)鄰居粒子的數(shù)量直接對色函數(shù)賦值；該方法計算效率與流體粒子數(shù)有關，粒子數(shù)增多時，計算速度較慢。文獻[12]使用δ+-SPH方法，通過對自由液面粒子位移的修正，對液面的壓力和光滑度進行處理；采用自適應粒子細化算法，解決粒子無序的問題；但是靠近邊界的粒子容易產(chǎn)生負壓。

本文針對固體邊界附近流體粒子支持域被截斷，流體粒子密度不連續(xù)，導致流體粒子在固體邊界處出現(xiàn)穿透、堆積以及流體表面變形難以處理等問題；通過改進壓力計算方法，使流體粒子在固體邊界附近正常流動，采用耦合力法，改進軟斥力，對流體粒子進行排斥，可以很好的穩(wěn)定交互界面壓力場；再對流體表面粒子位置進行校正，使流體表面粒子分布均勻，液面流動自然，提升了模擬效果。

1 SPH方法簡介

1.1 SPH方法

SPH方法是一種自適應、無網(wǎng)格的拉格朗日型粒子法。由LUCY[13]與GINGOLD 和MONAGHAN[14]在解決天體物理學時提出，后被廣泛應用于流體力學和固體力學中。通過離散帶有物理屬性值的粒子，由核近似和粒子近似計算得到標量A在位置r處的對應值為

其中，mj為粒子j的質(zhì)量；rj為位置；ρj為密度；Aj為在rj處的場量；函數(shù)W(r,h)是半徑為h的光滑核函數(shù)。由式(1)可知，其梯度和拉普拉斯算子只影響核函數(shù)，即

其中，符號?為梯度算子；?2為拉普拉斯算子。

核函數(shù)是SPH方法的核心，本文關于不同作用力的計算，采用文獻[15]中提出的核函數(shù)。

1.2 控制方程

將固體看作受剛體運動約束的流體，統(tǒng)一納入Navier-Stokes方程中求解?；谌蹩蓧嚎sSPH方法，動量方程和連續(xù)性方程為

其中，ρ為密度；v為速度；p為壓強；μ為黏度系數(shù)；F為額外力。

密度計算為

壓強由下列狀態(tài)方程計算得到

針對固體邊界附近流體粒子支持域被截斷的問題，為提高模擬方法的穩(wěn)定性和解決固體邊界處流體粒子密度不連續(xù)，流體粒子i的密度由周圍所有鄰居粒子密度之和計算得到。

本文采用文獻[16]的方法進行密度校正，將邊界粒子密度引入流固相互作用力的計算中。邊界粒子的質(zhì)量統(tǒng)一為mb，從而使邊界粒子的體積不依賴于質(zhì)量，保證了流體密度的穩(wěn)定性，邊界粒子的體積為

2 本文工作

在固體入水實驗中，需要對3 種邊界進行處理，即運動的固體邊界，模擬實驗需設置水槽兩側(cè)的靜止固壁邊界、以及流體表面的自由流動液面邊界。針對傳統(tǒng)方法的不足，需對上述3 種不同的邊界分別進行處理。

2.1 運動固體邊界的處理

通過改進耦合力法，結(jié)合邊界力易于計算以及虛粒子可以消除固體邊界附近流體粒子密度不連續(xù)的優(yōu)勢，處理運動的固體邊界。如圖1 所示，固體邊界粒子(紅色)，有自己的質(zhì)量和密度，在固體邊界粒子上布置排斥力，用于給流體粒子施加排斥力防止?jié)B透，運動固體邊界內(nèi)布置2 層虛粒子 (黑色)，虛粒子之間的距離設置為0.5h，用于彌補流體到達固體附近時，流體粒子的支持域h被截斷的問題，虛粒子基本屬性與流體粒子相同，速度與固體相同。

圖1 粒子分布初始狀態(tài) Fig.1 Initial state of the particle distribution

流體粒子i的壓力由周圍鄰居粒子j(j為流體粒子或固體粒子)的壓力通過插值得到，如圖2(a)所示。運動固體邊界粒子b的壓力只受到流體粒子的影響，通過鄰居流體粒子的壓力插值得到，如圖2(b)所示。

圖2 黑色是固體粒子，白色是流體粒子((a)流體粒子受力方式；(b)固體粒子受力方式) Fig.2 Black is solid particles and white is fluid particles ((a) Forced method of fluid particle;(b) Forced method of solid particle)

采用文獻[17]中壓力和黏性力的計算公式，則流體粒子i的壓力和黏性力分別為

流體壓力很大程度上受固體粒子的質(zhì)量影響，因此本文在流體粒子壓力計算中加入固體粒子質(zhì)量及密度。因為固體質(zhì)量與流體質(zhì)量差異較大，直接加入固體質(zhì)量的計算，容易產(chǎn)生壓力振蕩，所以本文在計算壓力時，對于固體邊界粒子質(zhì)量的貢獻選擇使用體積ibV代替，使邊界壓力計算不依賴于固體質(zhì)量，但又受到固體邊界粒子的影響，可以有效提高數(shù)值穩(wěn)定性，避免邊界處壓力振蕩引起的流體粒子在固體邊界附近分離或波動的現(xiàn)象。對于流體粒子i，改進的壓力計算式為

相應地，邊界附近流體粒子黏性力的計算同樣考慮固體粒子的貢獻，改進的流體粒子黏性力為

其中，mb為固體的質(zhì)量；g為重力加速度。

本文對文獻[9]中的耦合力法進行改進，原固體邊界粒子對流體粒子施加的排斥力模型為

由于原始耦合力法僅通過邊界粒子表征固體邊界，邊界粒子不具有任何物理屬性，會使流體粒子穿透邊界。而在運動固體與流體的交互過程中，流場的運動受固體的影響較大。本文在計算排斥力時，添加固體粒子屬性的計算，可以更好地平衡排斥力的大小，減輕壓力擾動。對接近固體邊界的流體粒子進行有限距離的排斥，有效防止流體粒子穿透固體邊界。針對運動固體入水問題，改進軟斥力模型為

其中，cs為色函數(shù)；χ,f (η)均為可變參數(shù)系數(shù)，通過粒子之間的位置和核半徑的關系，判斷粒子之間距離的函數(shù)[9]，即

其中，xij為流體粒子i到固體粒子j之間的向量；rij為2 個粒子之間的距離，Δd為2 個粒子的初始距離。

2.2 靜止固壁邊界的處理

對于固體落入靜水，水槽兩邊固壁邊界防穿透的設置，采用虛粒子法處理。在此，由于固體落入水槽底部時，固體粒子推開流體粒子后，流體粒子支持域被截斷，水槽的底部邊界處會出現(xiàn)流體粒子與固體邊界分離現(xiàn)象，故在水槽底部設置2 層虛粒子，用于固體到達水槽底部時，補充底部邊界處流體粒子的支持域。

在固壁邊界外2h(h為核半徑)距離內(nèi)鋪設虛粒子，方便下一步對流體表面粒子位置進行修正時，搜索表面修正區(qū)域內(nèi)的流體粒子，補充兩側(cè)流體粒子的支持域，避免搜索出兩側(cè)的流體粒子，節(jié)省計算。固壁邊界外虛粒子設置如圖3 所示，圓內(nèi)是支持域半徑為2h的流體粒子i的支持域。

圖3 固壁邊界虛粒子分布示意圖 (藍色為流體粒子，黑色為固壁邊界粒子，白色為虛粒子) Fig.3 Schematic diagram of the distribution of virtual particles at solid wall boundary (Blue is fluid particles,black is solid boundary particles and white is virtual particles)

2.3 流體表面粒子的位置修正

在固體入水問題中，流體表面的流體粒子與運動固體的相互作用通常與流體表面的變化和破裂有關。因為固體粒子對流體粒子施加了排斥力，流體表面會比未施加排斥力的表面更加粗糙，所以需要對流體表面粒子位置進行修正。

當固體與流體交互時，搜索表面的流體粒子，同時也搜索出流-固交互界面的粒子。最終流體表面和流-固交互界面2h(h是核半徑)距離內(nèi)的所有流體粒子，形成一個修正區(qū)域[12](圖4)。并通過修正流體表面粒子的法向速度，將凸出表面的流體粒子當前時刻的法向速度設為零，保留粒子切線方向的速度，再利用當前時刻的法向速度，更新流體粒子下一時刻的速度，修正表面流體粒子的位置，提升流體表面流動效果。

圖4 修正區(qū)域示意圖 Fig.4 Schematic diagram of correction area

本文采用文獻[18]中提出的距離場法搜索表面粒子。水槽的固壁邊界外設置的虛粒子剛好補充了左右邊界的搜索范圍，而運動固體邊界的虛粒子設置不滿足2h距離，不需要單獨進行交界面粒子的搜索，故搜索出來的流體粒子均為表面流體粒子，減少了計算量。

計算修正區(qū)域內(nèi)流體粒子相對表面的法線向量為

其中，ni是提取的自由表面流體粒子i的法向。若|ni(r)|＞2h，則粒子i為表面流體粒子；否則為內(nèi)部流體粒子。搜索出自由表面流體粒子后，令粒子i當前時刻法向速度為零，保持切線方向速度，流體粒子i方向速度如圖5 所示。

圖5 流體粒子i方向速度示意圖(藍色為流體粒子，黑色為固體粒子) Fig.5 Schematic diagram of the velocity in the direction of fluid particle i(Blue is a fluid particle and black is a solid particle)

于是，在t時刻，表面流體粒子的速度為

修正位置只修正表面飛散的流體粒子當前時刻的位置，而t+1 時刻的速度計算中，利用t時刻的法向速度繼續(xù)更新，即

更新下一步流體表面粒子的位置為

其中，vi是流體粒子i在t時刻的速度；vi+1是流體粒子i在t+1 時刻的速度。

3 實驗結(jié)果與分析

實驗在Windows 10 平臺進行，開發(fā)環(huán)境為visual studio 2013，渲染實驗在Unity2018 上完成。實驗配置為Intel(R) Core(TM) i5-5200U，2.20 GHz CPU，4 G 內(nèi)存，NVIDIA Geforce GT 620 顯卡。

通過經(jīng)典的二維圓柱入水實驗的模擬，實現(xiàn)了邊界力法、虛粒子法、耦合力法及本文方法以不同速度下落的入水場景，驗證了本文方法的有效性。實驗時間步長取0.05 s，粒子支持域半徑h取0.04 m。流體粒子個數(shù)為5 200 個，二維圓柱粒子個數(shù)為180 個。流體粒子初始密度取1 000 kg/m3，固體粒子密度取2 600 kg/m3。根據(jù)不同幀數(shù)下的對比，表明本文方法優(yōu)于其他3 種方法。

二維圓柱置于空中1 m 處，分別以6.5 m/s 和20 m/s 的速度下落。實驗初始設置如圖6 所示。

圖6 二維圓柱入水初始位置示意圖 Fig.6 Schematic diagram of the initial position of the two-dimensional cylinder entering the water

3.1 二維圓柱以6.5 m/s 和20 m/s 速度入水的實驗結(jié)果

從圖7～12 可以看出，本文方法與其他3 種方法相比，流體粒子與固體交互界面較光滑。邊界力 法流體表面較粗糙且流體粒子與固體粒子之間存在較大空隙，即流體因受力較大，導致流體粒子與固體分離程度較大；虛粒子法在整個模擬過程中有少許流體粒子穿透固體邊界的現(xiàn)象；耦合力法中，流體與固體的交互界面較粗糙，且邊界附近出現(xiàn)壓力振蕩。本文方法在整個模擬過程中較穩(wěn)定，固體邊界處流體粒子流動較自然，界面分離清晰，無穿透或分離程度較大的現(xiàn)象。從表1和2 可以看出，本文方法的模擬速度較其他3 種方法快。

表1 二維圓柱以6.5 m/s 速度入水各方法在 不同幀數(shù)下的用時比較(s) Table 1 Comparison of the time consumption of the two-dimensional cylinder entering the water at a speed of 6.5 m/s under different frames (s)

圖7 6.5 m/s 速度下各方法在第65 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.7 Comparison of simulation effects of each method at frame 65 at a speed of 6.5 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

圖8 6.5 m/s 速度下各方法在第87 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.8 Comparison of simulation effects of each method at frame 87 at a speed of 6.5 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

圖9 6.5 m/s 速度下各方法在第115 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.9 Comparison of simulation effects of each method at frame 115 at a speed of 6.5 m/s ((a) Boundary force method;(b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

圖10 20 m/s 速度下各方法在第25 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.10 Comparison of simulation effects of each method at frame 25 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Ethod of this article)

圖11 20 m/s 速度下各方法在第28 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.11 Comparison of simulation effects of each method at frame 28 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

圖12 20 m/s 速度下各方法在第31 幀的模擬效果對比 ((a)邊界力法；(b)虛粒子法；(c)耦合力法；(d)本文方法) Fig.12 Comparison of simulation effects of each method at frame 31 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

表2 二維圓柱以20 m/s 速度入水各方法在 不同幀數(shù)下的用時比較 (s) Table 2 Comparison of the time consumption of the two-dimensional cylinder entering the water at a speed of 20 m/s under different frames (s)

3.2 小球入水渲染實驗

本文分別渲染了以6.5 m/s 和20 m/s 速度下落的小球入水實驗，小球半徑為0.04 m，密度為2 600 kg/m3。并給出了不同速度下小球處于同一高度的實驗效果圖(圖13～16)。

圖13 小球位于0.45 m 高處((a) 6.5 m/s 速度第36 幀；(b) 20 m/s 速度第25 幀) Fig.13 The ball is at a height of 0.45 m ((a) Frame 36 at 6.5 m/s;(b) Frame 25 at 20 m/s)

圖14 小球位于0.20 m 高處 ((a) 6.5 m/s 速度第44 幀；(b) 20 m/s 速度第28 幀) Fig.14 The ball is at a height of 0.20 m ((a) Frame 44 at 6.5 m/s;(b) Frame 28 at 20 m/s)

圖15 小球落至底部((a) 6.5 m/s 速度第75 幀；(b) 20 m/s 速度第31 幀) Fig.15 The ball falls to the bottom ((a) Frame 75 at 6.5 m/s;(b) Frame 31 at 20 m/s)

圖16 小球落至底部((a) 6.5 m/s 速度第80 幀；(b) 20 m/s 速度第65 幀) Fig.16 The ball falls to the bottom ((a) Frame 80 at 6.5 m/s;(b) Frame 65 at 20 m/s)

從渲染實驗效果圖可以看出，本文方法在模擬固體入水時，可以呈現(xiàn)出固體落水的完整過程，且模擬效果真實自然，細節(jié)展現(xiàn)較好。

4 結(jié) 論

本文通過對運動固體邊界壓力計算方法和耦合力法的改進，有效補充了邊界處流體粒子的支持域、穩(wěn)定壓力場，防止流體粒子穿透固體邊界等問題。進一步通過對流體表面粒子位置的修正，對流體液面進行優(yōu)化，提升了模擬效果。通過模擬固體入水的二維實驗和渲染實驗，對提出的方法進行了驗證，結(jié)果表明了本文方法的可行性與有效性。