基于光滑粒子的氣泡實時仿真

祝 笛

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

0 引言

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

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

1 相關研究

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

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

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

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

2 基于SPH的流體建模

2.1 SPH原理

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

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

光滑核函數定義了粒子相鄰的支持域（見圖1），決定了計算結果的準確度。該函數通過將鄰域粒子的相應物理量進行加權平均，對場量作平滑處理。

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

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

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

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

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

2.2 控制方程

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

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

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

2.3 SPH流體方程

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

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

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

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

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

3 氣泡仿真

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

3.1 氣泡自由上浮

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

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

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

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

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

3.2 鄰域搜索

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

3.2.1 哈希表搜索法

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

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

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

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

3.2.2 區(qū)間劃分

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

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

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

3.3 氣泡懸浮

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

Table 1 Particle properties表1 粒子屬性

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

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

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

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

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

Step 3：對縱坐標y進行判斷：若y＞ymax，則令y=liquid.y，vy=liquid.vy，Fbuoyancy=-mg。

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

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

4 仿真結果及分析

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

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

Fig.5 Comparison of memory usage圖5 內存占用情況比較

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

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

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

5 結語

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