祝 笛

（江蘇科技大學(xué)計算機(jī)學(xué)院，江蘇鎮(zhèn)江 212003）

0 引言

流體場景的模擬一直是計算機(jī)圖形學(xué)領(lǐng)域的研究熱點。氣泡作為日常生活中一種常見的流體形態(tài)，其模擬被廣泛地應(yīng)用于影視特效、游戲動畫、軍事仿真等領(lǐng)域。由于氣泡運動狀態(tài)多變且存在氣液交互作用，如何真實、高效地模擬氣泡具有重要研究意義。

計算流體動力學(xué)（Computational Fluid Dynamics，CFD）［1］是以物理模型為基礎(chǔ)，基于物理規(guī)律描述流體的運動過程，可逼真地模擬出各種自然流體，因此成為國內(nèi)外計算機(jī)圖形學(xué)領(lǐng)域的研究熱點。隨著計算機(jī)技術(shù)的發(fā)展，基于數(shù)值求解的CFD 成為流體仿真研究的主流方法，目前主要分為歐拉網(wǎng)格法和拉格朗日粒子法，其中拉格朗日粒子法是將流體離散化為大量粒子，通過求解單個粒子的物理量進(jìn)行流體仿真，可有效避免流體變形時網(wǎng)格扭曲造成的精度破壞，成為流體高精度仿真的主要方法。

1 相關(guān)研究

光滑粒子流體動力學(xué)（Smoothed Particle Hydrodynamics，SPH）是近年來發(fā)展起來的一種拉格朗日粒子法，作為一種無網(wǎng)格的流體模擬方法，其運用大量攜帶物理屬性的粒子進(jìn)行流體建模［2-3］，能夠精確模擬出流體運動的細(xì)節(jié)，有效增強了流體模擬的真實感，成為目前流體模擬的主要方法。例如，Monaghan［4］采用SPH 方法進(jìn)行自由表面流動的仿真；Müller 等［5］采用SPH 方法實現(xiàn)了在交互式幀率的基礎(chǔ)上對可壓縮流體的模擬。

氣泡作為液體的一種常見細(xì)節(jié)特征，其模擬效果對于流體仿真的現(xiàn)實感和直觀性具有重要影響，許多學(xué)者為改善氣泡模擬效果進(jìn)行了探索。例如，Kim 等［6］將粒子與網(wǎng)格相結(jié)合，對氣體和液體之間的數(shù)值求解劃分階段，采用空間平均方法分散氣泡流實現(xiàn)了微型氣泡復(fù)雜場景的模擬；武小龍等［7］結(jié)合流體體積法（Volume of Fluid）和Skin-Onion 法模擬了氣泡的生成和運動，但僅模擬了單個氣泡的獨立運動，并未考慮鄰近氣泡之間的交互影響；Müller等［8］采用SPH 方法建立了密度比為10∶1的多相流模型，將浮力引入粒子受力模型，對氣泡和液體的交互進(jìn)行模擬，但氣泡和液體交界處存在錯誤的密度估值，影響了流體壓強的計算；Cleary 等［9］提出一種基于耦合離散氣泡的SPH流體模型，可獨立模擬氣泡和液體。該模型將每個氣泡描述為一個離散的實體，采用SPH 方法模擬液體，然后在氣泡與其周圍液體之間引入一個拉力模型，實現(xiàn)了氣泡與液體之間交互作用的模擬。此外，氣泡在液體中的運動表現(xiàn)出各種復(fù)雜特征，如變形、上浮、聚集等現(xiàn)象，這些運動細(xì)節(jié)的模擬可以在一定程度上增強流體交互模擬的真實感，許多學(xué)者進(jìn)行了相關(guān)研究。例如，Ming 等［10］采用SPH 方法建立了多相流模型，模擬了復(fù)雜氣泡上升運動的現(xiàn)象，通過優(yōu)化壓力和表面張力解決了小尺度、高密度比的氣泡上升問題；Zheng 等［11］建立了一個基于連續(xù)流體模擬的框架，引入多相流界面跟蹤的區(qū)域水平集、半隱式表面張力模型等方法，模擬出了逼真的氣泡動態(tài)效果；Gao 等［12］將溫度因素引入流體模型，建立了一種流體受熱的簡化模型，將潛熱模型與相關(guān)物理規(guī)律相結(jié)合，實現(xiàn)了熱條件下氣泡的產(chǎn)生和交互模擬；Arai 等［13］采用SPH 和流體體積的方法模擬了氣泡上升的形變，成功捕獲了氣泡形變細(xì)節(jié)的瞬間形態(tài)。

在氣泡仿真中存在氣泡數(shù)量較多且計算消耗較大的問題，如何提升仿真計算效率成為流體研究的重要內(nèi)容。例如，段興鋒等［14］利用三維空間網(wǎng)格對整個模擬區(qū)域進(jìn)行均勻劃分，采用并行前綴求和與并行計數(shù)排序進(jìn)行鄰域粒子的查找，同時利用CUDA 并行加速，提升了仿真計算效率；龍廳［15］在鏈表搜索的基礎(chǔ)上優(yōu)化了存放粒子單元的尺寸大小，極大地提高了三維情況下粒子的搜索效率；葉靜然等［16］將粒子的樹形結(jié)構(gòu)存儲于GPU 上，并直接利用樹形搜索相鄰粒子，提高了模擬速度并仿真了管道水流。

然而，以上研究并沒有涵蓋氣泡仿真的全部問題。例如，由于內(nèi)部空氣密度小于液體密度，氣泡在液體中生成后呈自由上浮運動狀態(tài)。在運動過程中，氣泡表現(xiàn)出形態(tài)變化、融合、懸浮及破裂等多種復(fù)雜現(xiàn)象，且氣泡與液體之間存在交互作用。因此，本文采用基于SPH 的多相模型［17］，通過求解納維—斯托克斯［18］（Navier-Stokes Equation，N-S）方程分別計算氣泡和液體的物理量以進(jìn)行流體模擬，然后加入作用力實現(xiàn)兩相流體的交互。在氣泡上升至液體表面后，氣泡會懸浮于液面，并且與液面的起伏相耦合。為此，本文采用基于液體粒子搜索的氣液耦合方法［19］，借助鄰域液體的運動狀態(tài)更新氣泡的相關(guān)屬性，實現(xiàn)氣液運動起伏的耦合。此外，在液體模型中存在粒子數(shù)量較大且運動狀態(tài)多變的情況，搜索鄰域液體粒子時存在搜索空間復(fù)雜度較高且計算消耗較大的問題。為降低搜索空間的復(fù)雜度，提升仿真實時性，本文提出一種基于區(qū)間劃分結(jié)合哈希表的鄰域搜索方法，從3 個維度對液體空間進(jìn)行區(qū)間劃分，再將劃分出的區(qū)間進(jìn)行網(wǎng)格編號并映射至哈希表，計算并更新粒子索引，縮減鄰域搜索范圍。該法在保證仿真視覺效果的同時，有效降低了時間和空間復(fù)雜度，減少了仿真過程的計算消耗，提升了氣泡仿真的實時性。

2 基于SPH的流體建模

2.1 SPH原理

SPH 是基于插值理論的方法，將流體離散為帶有密度、速度、質(zhì)量等一系列物理量的大量粒子，流體空間內(nèi)任意位置的粒子屬性可以通過其鄰域采樣點的屬性插值計算得到，計算公式為：

式中，i表示所求粒子，j表示鄰域內(nèi)的所有粒子，Ai、Aj分別表示粒子i、j的物理屬性，mj表示粒子j的質(zhì)量，ρj表示粒子j的密度，W(xij，h)表示光滑核函數(shù)，xij表示xi-xj，xi、xj分別為粒子i、j的位置，h為該函數(shù)的支持半徑。

光滑核函數(shù)定義了粒子相鄰的支持域（見圖1），決定了計算結(jié)果的準(zhǔn)確度。該函數(shù)通過將鄰域粒子的相應(yīng)物理量進(jìn)行加權(quán)平均，對場量作平滑處理。

Fig.1 Particle adjacent support domain圖1 粒子相鄰支持域

本文在流體建模時參照文獻(xiàn)［8］，采用WPloy6核函數(shù)計算密度，表示為：

采用Wspiky(r，h)核函數(shù)計算壓力，表示為：

采用Wviscosity(r，h)核函數(shù)計算粘滯力，表示為：

上述核函數(shù)在粒子距離較近時能夠產(chǎn)生較大排斥力，避免出現(xiàn)粒子成簇的情況。此外，在復(fù)雜交互的仿真中，其可使粘滯力計算更加穩(wěn)定。

2.2 控制方程

N-S 方程一般用于不可壓縮流體的描述，遵循流體運動的質(zhì)量守恒和動量守恒定律。該方程的一般形式為：

式中，ρ表示流體密度，v表示流體速度，p表示流體壓強，μ表示粘度系數(shù)；-(v· ?)v表示流體的對流加速力，該力場使得流體運動時的空間變化產(chǎn)生加速度；-?p表示流體壓力，該力場使得流體向壓強較小的區(qū)域流動；μ?2v表示粘滯力，該力場是由于流體自身粘性而產(chǎn)生的阻力，其大小受到流體粘度系數(shù)與速度差的影響；f表示單位體積流體的外力場，如重力場。

采用恒定數(shù)量的粒子描述流體運動時，單個粒子的質(zhì)量同樣恒定，因此粒子模型滿足式（5）中的質(zhì)量守恒方程。此外，粒子的移動與流體一致，不存在對流時的數(shù)值耗散，因此式（5）中的對流項可以忽略不計。流體受力可分為壓力、粘滯力和外力3 種，這3 種力場決定了粒子速度的瞬時變化量。式（5）的控制方程可簡化為：

2.3 SPH流體方程

求解N-S 方程時，需將其轉(zhuǎn)換為SPH 方法的插值形式。根據(jù)式（1）可知，在位置x處的粒子密度為：

在流體仿真中，通過粒子密度可以計算出理想氣體狀態(tài)下的壓力場為：

式中，k表示氣體常量，在溫度恒定時k值不變；ρ0表示流體的靜止密度?；赟PH 原理將粒子所受的各個力場表示為物理屬性的插值模式，根據(jù)式（1）和式（6）可以得到：

根據(jù)粒子的加速度ai可以計算出時間步長內(nèi)粒子i的速度與位置分別為：

式中，vi(t+Δt)和xi(t+Δt)分別為更新后粒子i的速度與位置，Δt為時間步長。

3 氣泡仿真

氣泡在整個生命周期中可以分為3 種狀態(tài)，分別為氣泡產(chǎn)生、氣泡上浮運動以及氣泡到達(dá)液面懸浮后破裂。在自然條件下，可將液體中產(chǎn)生的氣泡當(dāng)作一種特殊的流體，氣泡與液體之間可以視為彼此獨立且相互作用的關(guān)系。本文基于SPH 方法分別對液體與氣泡建立仿真模型，在該流體模型的基礎(chǔ)上對氣泡在液體中的運動狀態(tài)進(jìn)行模擬。

3.1 氣泡自由上浮

對液體粒子的區(qū)域進(jìn)行網(wǎng)格劃分，在液體底部設(shè)置隨機(jī)生成氣泡的區(qū)域Ω(X+?，Y，Z+φ)，其中?和φ分別為浮動范圍。每個時間步長內(nèi)，在該區(qū)域內(nèi)隨機(jī)位置上生成一定數(shù)量的氣泡，在氣泡初生時設(shè)置氣泡粒子的大小、生命周期和初始速度等屬性值。

氣泡粒子生成后作上升運動，由于受到各種作用力的影響，粒子的運動速度不斷發(fā)生改變。根據(jù)文獻(xiàn)［5］研究結(jié)果可知，除SPH 流體模型中的壓力和重力外，氣泡運動還受到凝聚力、拉力和浮力的影響。在液體空間中，位置相近的氣泡會相互聚集，凝聚力控制著鄰近氣泡粒子的凝聚效果，該力表示為：

拉力使得氣泡粒子與液體粒子之間相互作用，實現(xiàn)兩相流體交互的效果，該力表示為：

在氣泡未到達(dá)液體表面時，浮力推動了氣泡粒子整體呈現(xiàn)上升運動，該力表示為：

由于氣泡運動時還受到壓力等作用力，結(jié)合式（10）和式（15）可以得到氣泡粒子i在所受合力作用下產(chǎn)生的加速度，再根據(jù)式（11）更新氣泡粒子的狀態(tài)信息。每個時間步長內(nèi)大量氣泡粒子的狀態(tài)持續(xù)更新，從而模擬出氣泡上浮運動。

3.2 鄰域搜索

在懸浮狀態(tài)下，氣泡與液面運動的耦合需要基于對最鄰近液體粒子的搜索，其實質(zhì)是計算并比較各液體粒子與該氣泡粒子之間的距離差值。然而液體模型中存在粒子數(shù)量較大且運動狀態(tài)多變的情況，搜索鄰域液體粒子時存在一定的空間復(fù)雜度。為降低搜索空間的復(fù)雜度，提升仿真實時性，本文提出一種基于區(qū)間劃分結(jié)合哈希表的鄰域搜索方法。

3.2.1 哈希表搜索法

哈希表搜索法是一種高效的鏈表法［20］，該法將液體粒子空間劃分為若干網(wǎng)格單元，然后將網(wǎng)格位置與哈希表相映射，并將粒子信息存儲于哈希表的索引中。首先，需要將液體空間均勻劃分為若干個邊長為l的網(wǎng)格，設(shè)置邊長為光滑核半徑（l=h），并為每個網(wǎng)格賦予編號；然后，利用哈希函數(shù)將粒子的坐標(biāo)(i，j，k)存儲到哈希表相應(yīng)的索引位置，根據(jù)網(wǎng)格號從哈希表中獲取粒子索引。由于粒子的位置是根據(jù)固定的網(wǎng)格邊長進(jìn)行離散的，粒子的位置(x，y，z)在存入索引時需要除以網(wǎng)格邊長并取整，即i=x/l，j=y/l，k=z/l。

當(dāng)搜索粒子i的鄰域關(guān)系時，哈希表搜索法只需要根據(jù)i的位置找到距離最近的網(wǎng)格點，即在二維空間中僅需查找3 × 3 的網(wǎng)格，三維空間中查找3 × 3 × 3 的網(wǎng)格便可以搜索到鄰域粒子，而無需遍歷計算所有粒子。鄰域搜索范圍如圖2所示。

Fig.2 Neighborhood search range圖2 鄰域搜索范圍

當(dāng)空間區(qū)域內(nèi)有N個粒子時，如果進(jìn)行全局鄰域搜索，那么需要計算的次數(shù)為N(N-1)，時間T(n)=n(n-1)=n2-n，其時間復(fù)雜度為O(n2)。將區(qū)域劃分為n個網(wǎng)格后，單個網(wǎng)格中的粒子為N/n個，每個粒子鄰域搜索的運算次數(shù)為27*N/n，那么所有粒子的運算次數(shù)為N*(27*N/n)，時 間T(n)=(27*N/n)n，其時間復(fù)雜度為O(n)。

3.2.2 區(qū)間劃分

流體的狀態(tài)信息在每一幀都會更新，僅采用哈希表搜索懸浮氣泡的鄰域液體粒子時需要時刻更新該氣泡粒子和液體粒子的哈希索引。然而氣液起伏耦合僅需要查找最鄰近的液體粒子，可以將氣泡所在液面位置作為切入點，通過劃分液體區(qū)間縮減搜索空間。

假設(shè)計算域Ω 有N個液體粒子，以懸浮氣泡的位置(x，y，z)為中心，以光滑核半徑h為劃分距離。首先，在Χ方向上作劃分，將|xi-h|和|xi+h|區(qū)間內(nèi)的液體粒子搜索至集合A[C(xk)，...，C(xs)]；然后，在Y 方向上對集合A[C(xk)，...，C(xs)]中的粒子進(jìn)行判斷，將在|yi-h|和|yi+h|區(qū)間內(nèi)的粒子更新至集合B[C(xk)，...，C(xs)]。同理，以|zi-h| 和|zi+h| 為區(qū)間劃分集合B[C(xk)，...，C(xs)]，可得到結(jié)合D[C(xk)，...，C(xs)]。液體粒子所在空間經(jīng)劃分后，更新范圍可以從Ω 縮減至集合D；最后，將集合D作為網(wǎng)格區(qū)域映射到哈希表，計算并更新粒子索引，區(qū)間劃分結(jié)構(gòu)如圖3 所示。由此，在仿真中，每幀更新時僅需要計算D區(qū)間內(nèi)的液體粒子，可有效降低計算消耗和空間復(fù)雜度。

Fig.3 Neighborhood search range after interval division圖3 區(qū)間劃分后鄰域搜索范圍

3.3 氣泡懸浮

氣泡懸浮于液面一小段時間后會破裂消失。氣泡薄膜破裂后可以分散成大量小氣泡，但由于其物理模型的復(fù)雜性，本文采用文獻(xiàn)［6］中的方法，通過一個簡易模型模擬氣泡在液面的浮動和破裂，在達(dá)到視覺效果的同時簡化了大量數(shù)值計算。在定義粒子屬性時采用相同的結(jié)構(gòu)體，該結(jié)構(gòu)體存儲了密度、速度、生命周期及位置等一系列屬性值，具體如表1所示。

Table 1 Particle properties表1 粒子屬性

為了模擬氣泡在液面的持續(xù)懸浮狀態(tài)，以及氣泡與液面起伏的耦合，需要基于鄰域液體粒子更新氣泡的運動狀態(tài)，實現(xiàn)氣泡懸浮的視覺效果。在每次計算中，將位置最高的水粒子縱坐標(biāo)記為ymax，氣泡粒子的縱坐標(biāo)大于ymax時視為到達(dá)液面，記為懸浮狀態(tài)。將氣泡粒子的縱坐標(biāo)y和垂直速度vy修改為最近液體粒子liquid相應(yīng)的屬性值，并將其浮力修改為與重力值相等但方向相反的力，使得氣泡懸浮于水面。此外，設(shè)置氣泡粒子的生命值life，在每次迭代中遞減粒子生命值。當(dāng)粒子生命值小于零時刪除該粒子，并在生成區(qū)域新增粒子，實現(xiàn)氣泡的消亡和新生。其仿真步驟如下：

輸入：粒子初始化物理屬性

輸出：粒子位置及狀態(tài)變化的動畫

Step 1：搜索并存儲鄰域液體粒子。

Step 2：計算粒子加速度ai，更新粒子速度、位置和生命值。

Step 3：對縱坐標(biāo)y進(jìn)行判斷：若y＞ymax，則令y=liquid.y，vy=liquid.vy，F(xiàn)buoyancy=-mg。

Step 4：對生命值lifei進(jìn)行判斷：若lifei＜0，清除該氣泡粒子。

Step 5：循環(huán)上述步驟直至結(jié)束。

4 仿真結(jié)果及分析

實驗環(huán)境為Intel（R）Core（TM）i5-8265U CPU，8G RAM，顯卡為NVIDIA GeForce MX250。采用Visual Studio2019 和Unity2019 進(jìn)行建模與渲染。實驗中，設(shè)置粒子質(zhì)量為0.02g，液體粒子間距為0.8m，系數(shù)取值為kc=12，kb=3.5，kd=80，kmax=1.5。在同一場景中取不同粒子數(shù)量，將本文方法與全局鄰域搜索、哈希表搜索法的平均幀頻、內(nèi)存占用及粒子運動狀態(tài)進(jìn)行比較分析，結(jié)果分別如圖4、圖5、圖6所示（彩圖掃OSID 碼可見，下同）。

由圖4、圖5 可以看出，在粒子數(shù)量相同的情況下，3 種方法的內(nèi)存占用差別很小，而本文方法在仿真中的平均幀頻最高。隨著粒子數(shù)量的增多，本文方法的平均幀頻降速較為緩慢，受粒子數(shù)量的影響較小。由圖6 可以看出，3 種方法均可模擬氣泡的運動狀態(tài)，本文方法和哈希搜索法在60s 后氣泡到達(dá)液面并呈懸浮狀態(tài)，而全局搜索方法在100s后氣泡開始懸浮。

Fig.5 Comparison of memory usage圖5 內(nèi)存占用情況比較

Fig.6 Particle state comparison圖6 粒子狀態(tài)比較

圖7 展示了采用本文方法實現(xiàn)的氣泡在液體中運動的仿真場景，從左至右依次為第30 幀、第70 幀和第110幀，氣泡從液體底部生成后逐漸自由上浮運動，在其到達(dá)液面后呈懸浮狀態(tài)并與液面起伏耦合?？梢钥闯?，本文方法實現(xiàn)了氣泡在液體中運動狀態(tài)的仿真，有效提升了計算效率，保證了仿真的實時性。

Fig.7 Bubble motion simulation in the liquid圖7 氣泡在液體中的運動仿真

5 結(jié)語

本文基于SPH 方法建立了氣泡和液體的粒子模型，分析粒子受力模型并更新其狀態(tài)信息，對流體運動進(jìn)行模擬。在模擬氣泡懸浮時，借助最鄰近液體粒子的運動狀態(tài)更新氣泡相關(guān)屬性，實現(xiàn)了氣液運動起伏的耦合。同時采用液體計算域進(jìn)行區(qū)間劃分，對縮減的計算域進(jìn)行網(wǎng)格劃分，利用哈希表建立了粒子位置與網(wǎng)格的關(guān)系，降低了鄰域粒子搜索的空間維度，縮減了搜索和計算時間。通過與已有研究結(jié)果進(jìn)行比較，發(fā)現(xiàn)本文方法在實現(xiàn)氣泡運動仿真的同時有效提高了鄰域搜索效率，保證了氣泡仿真的實時性。本文方法對于煙、云等其他流體的仿真研究具有很好的借鑒價值。然而，本文目前僅對模擬的計算效率進(jìn)行了優(yōu)化，后續(xù)將會對氣泡融合及形變模擬進(jìn)行深入研究。