固體入水的邊界壓力處理及表面粒子位置的修正

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

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

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

光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)方法在處理大變形和快速移動(dòng)的界面具有獨(dú)特的優(yōu)勢(shì)。因此，SPH方法被廣泛應(yīng)用于各種流-固耦合問題。但傳統(tǒng)SPH方法在模擬運(yùn)動(dòng)固體邊界時(shí)，由于固體邊界附近流體粒子支持域被截?cái)?，流體粒子密度不連續(xù)，使得計(jì)算精度較差，會(huì)導(dǎo)致邊界處壓力振蕩等問題；而流體表面因?yàn)槭艿焦腆w的沖擊力，流體液面較粗糙。

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

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

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

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

1 SPH方法簡(jiǎn)介

1.1 SPH方法

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

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

其中，符號(hào)?為梯度算子；?2為拉普拉斯算子。

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

1.2 控制方程

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

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

密度計(jì)算為

壓強(qiáng)由下列狀態(tài)方程計(jì)算得到

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

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

2 本文工作

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

2.1 運(yùn)動(dòng)固體邊界的處理

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

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

流體粒子i的壓力由周圍鄰居粒子j(j為流體粒子或固體粒子)的壓力通過插值得到，如圖2(a)所示。運(yùn)動(dòng)固體邊界粒子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)

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

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

相應(yīng)地，邊界附近流體粒子黏性力的計(jì)算同樣考慮固體粒子的貢獻(xiàn)，改進(jìn)的流體粒子黏性力為

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

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

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

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

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

2.2 靜止固壁邊界的處理

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

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

圖3 固壁邊界虛粒子分布示意圖 (藍(lán)色為流體粒子，黑色為固壁邊界粒子，白色為虛粒子) 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 流體表面粒子的位置修正

在固體入水問題中，流體表面的流體粒子與運(yùn)動(dòng)固體的相互作用通常與流體表面的變化和破裂有關(guān)。因?yàn)楣腆w粒子對(duì)流體粒子施加了排斥力，流體表面會(huì)比未施加排斥力的表面更加粗糙，所以需要對(duì)流體表面粒子位置進(jìn)行修正。

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

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

本文采用文獻(xiàn)[18]中提出的距離場(chǎng)法搜索表面粒子。水槽的固壁邊界外設(shè)置的虛粒子剛好補(bǔ)充了左右邊界的搜索范圍，而運(yùn)動(dòng)固體邊界的虛粒子設(shè)置不滿足2h距離，不需要單獨(dú)進(jìn)行交界面粒子的搜索，故搜索出來的流體粒子均為表面流體粒子，減少了計(jì)算量。

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

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

圖5 流體粒子i方向速度示意圖(藍(lán)色為流體粒子，黑色為固體粒子) 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時(shí)刻，表面流體粒子的速度為

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

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

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

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

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

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

二維圓柱置于空中1 m 處，分別以6.5 m/s 和20 m/s 的速度下落。實(shí)驗(yàn)初始設(shè)置如圖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 速度入水的實(shí)驗(yàn)結(jié)果

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

表1 二維圓柱以6.5 m/s 速度入水各方法在 不同幀數(shù)下的用時(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 幀的模擬效果對(duì)比 ((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 幀的模擬效果對(duì)比 ((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 幀的模擬效果對(duì)比 ((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 幀的模擬效果對(duì)比 ((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 幀的模擬效果對(duì)比 ((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 幀的模擬效果對(duì)比 ((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ù)下的用時(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 小球入水渲染實(shí)驗(yàn)

本文分別渲染了以6.5 m/s 和20 m/s 速度下落的小球入水實(shí)驗(yàn)，小球半徑為0.04 m，密度為2 600 kg/m3。并給出了不同速度下小球處于同一高度的實(shí)驗(yàn)效果圖(圖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)

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

4 結(jié) 論

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