基于物理的不可壓縮流體模擬技術(shù)綜述①

吳德陽 唐 勇 劉浩陽 劉培藝 汪國平**

(*燕山大學(xué)信息科學(xué)與工程學(xué)院 秦皇島 066004)

(**北京大學(xué)信息科學(xué)技術(shù)學(xué)院 北京 100871)

(***河北省計算機虛擬技術(shù)與系統(tǒng)集成重點實驗室 秦皇島 066004)

(****北京市虛擬仿真與可視化工程研究中心 北京 100871)

0 引 言

近年來，隨著數(shù)字信息、圖形技術(shù)的快速發(fā)展，虛擬現(xiàn)實作為一種創(chuàng)建和體驗虛擬世界的仿真技術(shù)，能夠模擬出與自然景象極其逼真的效果，使人如同身臨其境，該技術(shù)已經(jīng)被廣泛應(yīng)用于消防、油畫生成、游戲場景、影視特效等領(lǐng)域[1-3]。而流體模擬作為計算機圖形學(xué)的一個研究分支，可以豐富虛擬場景，增強現(xiàn)實感、沉浸感，因此對于流體仿真的研究具有重要意義。隨著計算機硬件的計算能力不斷提升，流體仿真的計算速度得到了很大改善，但對于實時仿真還存在很大的提升空間。通常，許多研究者通過犧牲部分仿真效果來滿足實時性的要求，但游戲、影視特效等領(lǐng)域?qū)崟r性和真實感兩者都具有很高的要求，因此，如何快速、真實地模擬不可壓縮流體仍是計算機圖形學(xué)領(lǐng)域面臨的挑戰(zhàn)。

在單相的流體模擬中，求解不可壓縮性一直是流體模擬的研究熱點，雖然近年來隨著對流體模擬算法的不斷改進，取得了一些新穎的研究成果[4,5]，但由于受表面張力[6]、邊界和空氣壓強等因素的影響，對于求解邊界處的不可壓縮性還存在一定的困難。而在真實感上，研究者為了達到更精細的效果，常常將網(wǎng)格進行細化或增加粒子數(shù)，雖然細節(jié)效果更明顯，但與之帶來的后果是增加了算法的計算耗時。為此，有些文獻則利用GPU的并行計算能力來提高算法的實時性能，如文獻[7,8]利用GPU對流體模擬進行加速的策略。在求解不可壓縮流體的基礎(chǔ)上，流體與固體的雙向耦合成為了流固耦合的難點，相比于單種流體運動的模擬，流固耦合除了需要求解流體的不可壓之外，還需要考慮流-固界面穿透問題。兩相流或者多相流之間耦合的主要難點在于不同相之間的交界面特性，以及不同流體之間的混合、兩種化學(xué)物質(zhì)參與反應(yīng)所產(chǎn)生的效果。多相流界面主要分為液-固界面、氣-液界面、液-液界面以及氣-液-固界面4種類型，不同相之間構(gòu)成的界面也有所不同，如液相和固相構(gòu)成的界面大多為尖銳性界面，而液相與液相構(gòu)成的界面大多為彌散界面。目前處理尖銳性多相流界面的方法主要有有限體積法[9]和粒子水平集方法[10]，而描述彌散界面的方法主要包含相場法[11]和格子波爾茲曼法( lattice Boltzmann method，LBM )[12]等。目前國內(nèi)外已有部分文獻對流體仿真技術(shù)進行了總結(jié)[13-15]，為了與這些綜述文獻區(qū)分開，本文主要圍繞2003年至2019年基于物理的不可壓縮流體仿真算法進行梳理與總結(jié)，并通過比較分析，給出了當(dāng)前流體仿真技術(shù)所面臨的問題和挑戰(zhàn)。

本文第1節(jié)介紹了不可壓流體納維-斯托克斯方程(Navier-Stokes equation，N-S)的求解過程。第2節(jié)介紹了基于物理的流體仿真方法，分別為Lagrange粒子的方法、基于Eular網(wǎng)格的方法以及基于Lagrange粒子-Eular網(wǎng)格的混合方法。第3節(jié)對第2節(jié)中的物理仿真方法進行對比分析和評價。第4節(jié)提出現(xiàn)有研究存在的不足并對未來值得關(guān)注的研究方向進行初步探討。第5節(jié)總結(jié)全文。

1 不可壓流體的N-S方程求解

不可壓縮流體的運動通常由著名的N-S方程來描述，主要由質(zhì)量守恒方程、動量守恒方程、能量守恒方程3大方程組成。在不可壓流體仿真中，為了簡化計算，通常將密度設(shè)為恒定值并忽略溫度的影響，其基本形式如下[16]：

(1)

質(zhì)量守恒：▽·u=0

(2)

其中，k▽2u為粘力項，(u·▽)u為平流項，▽p為壓力項，f為外力項，式(2)為流體的不可壓條件，N-S方程求解流程如圖1所示。

圖1 N-S方程求解流程圖

具體求解步驟如下：(1)初始化速度場、壓力場等。(2)施加外力，主要為重力。(3)對速度進行平流操作，主要有前向平流和后向平流。前向平流根據(jù)當(dāng)前時刻的中心速度，向前計算下一時刻的位置，并根據(jù)下一時刻的物理量通過雙線性插值操作更新周圍4個網(wǎng)格的物理量，最后將下一時刻的物理量作為當(dāng)前時刻的物理量。后向平流根據(jù)當(dāng)前的中心速度，向后回溯Δt時刻的位置，然后根據(jù)Δt時刻位置的物理量進行雙線性插值更新周圍4個網(wǎng)格的物理量。(4) 進行壓力投影，使流體滿足不可壓條件。(5) 添加粘力。(6) 添加步驟(4)得到的壓強更新速度。

2 基于物理的不可壓流體模擬方法

基于物理的不可壓流體仿真方法主要由3大類構(gòu)成，分別為Lagrange粒子法、Eular網(wǎng)格法和Lagrange-Eular混合法。Lagrange粒子法主要以光滑粒子動力學(xué)(smoothed particle hydrodynamics,SPH)法、基于位置的流體法(position based fluid,PBF)和LBM法為代表；Eular網(wǎng)格主要涉及高度場網(wǎng)格、水平集方法和相場法；Lagrange-Eular混合法主要以質(zhì)點網(wǎng)格法(particle in cell，PIC)和流體隱式粒子法(fluid implicit particle，F(xiàn)LIP)為代表，不可壓流體仿真方法框架如圖2所示。

圖2 流體仿真方法框架

2.1 Lagrange粒子法

Lagrange粒子法的基本思想是將連續(xù)的流體離散成一堆攜帶物理量(密度、速度、質(zhì)量等)的粒子，然后通過計算所有粒子之間的相互作用力來描述流體質(zhì)點隨時間變化的規(guī)律。本節(jié)主要介紹光滑流體動力學(xué)和基于位置的流體模擬方法以及這2種方法的應(yīng)用。

2.1.1 SPH法

在Lagrange粒子法中，SPH是流體模擬方法最為典型的方法之一。由于該方法實現(xiàn)簡單、且能處理拓撲變化較大的邊界效果，如噴濺現(xiàn)象等，因此被廣泛應(yīng)用于單種流體模擬[17]、流固耦合[18]、多相流[19]等，其基本模型如圖3所示。

圖3 SPH粒子模型

SPH的基本思想是通過在位置x處對光滑核半徑r范圍內(nèi)的n個粒子累加求和，得到該處的物理屬性的累積量[20]：

(3)

模擬的大致步驟如下：(1) 初始化粒子信息，如位置等。(2) 利用式(3)求出粒子的密度ρ。(3) 根據(jù)上一步驟得到的粒子密度ρ計算粒子的壓強進而推導(dǎo)出由半徑r范圍內(nèi)鄰居粒子對該粒子施加的壓力F。(4) 根據(jù)步驟(3)得到的壓力F利用牛頓第二定律求出粒子的加速度a。(5) 利用粒子的加速度a求出粒子的速度。(6) 根據(jù)速度更新粒子的位置。

(1) 單種流體的仿真

最初的SPH由Stam[20]引入計算機圖形學(xué)，用于模擬火焰等現(xiàn)象。隨后Müller等人[21]將SPH用于流體模擬，模擬的流體效果相比于Eular網(wǎng)格法的真實感更強，其能模擬水珠飛濺的效果，如圖4所示，但需要離線渲染，且該方法采用理想氣體關(guān)聯(lián)壓力和密度，由于與實際的流體壓力和密度不同，理想氣體關(guān)聯(lián)壓強會導(dǎo)致模擬的流體存在較高的壓縮性。

圖4 SPH粒子法仿真效果[21]

微可壓SPH法(weakly compressible SPH，WCSPH)[22]，針對文獻[21]的高壓縮性問題，Becker等人[22]引入一種新的表面張力模型，可以避免控制方程中導(dǎo)數(shù)項的計算，同時為了防止數(shù)值耗散，在模型中添加了人工粘度來提高數(shù)值的穩(wěn)定性。WCSPH很好地改善了傳統(tǒng)SPH的不可壓問題，但需要用較小的時間步長來維持模擬的穩(wěn)定性。

預(yù)測校正SPH法(predictive corrective incompressible SPH，PCISPH)[23]，為了克服WCSPH對時間步長的敏感性問題，Solenthaler等人[23]通過預(yù)測-矯正粒子的位置和密度，減少上一幀與下一幀之間的密度差，并關(guān)聯(lián)密度與壓強的關(guān)系，實現(xiàn)不可壓縮性，相比WCSPH方案，其時間步長更長，數(shù)值求解更穩(wěn)定。

局部泊松SPH(local Poisson SPH， LPSPH)[24]，He等人[24]為了彌補WCSPH大密度誤差和PCISPH全局求解泊松方程耗時的不足，利用局部作用機理，將泊松方程轉(zhuǎn)換成積分形式，然后通過離散化將積分方程變換為局部壓力積分域，減少求解泊松方程的時間，同時為了減少密度誤差，將局部壓力集成到預(yù)測矯正器中，可以避免由于密度波動引起的不穩(wěn)定問題。

隱式不可壓SPH法(implicit incompressible SPH， IISPH)[25]，與WCSPH和PCISPH相比，LPSPH在時間步長和收斂速度上具有很大的改進，但LPSPH受模擬場景規(guī)模的限制，只能用于小規(guī)模流體的模擬。為此，Ihmsen等人[25]通過離散化壓力泊松方程，并在最終的速度更新中使用壓力項來表示，可以有效地改進大場景下的時間步長，且能使密度差降低至0.01%。He等人[26]提出一種基于SPH的自由表面流中小尺度薄片特征的魯棒模擬方法，利用擴散模型將自由表面建模成非常薄的區(qū)域，即在界面區(qū)域內(nèi)的粒子表示成1，界面外則表示為0，界面處則是從0到1光滑過渡，這種建模方法可以在界面粒子較少的情況下，仍然可以魯棒地表示流體拉伸過程中的薄片特征。Bender等人[27]提出一種無散度的SPH方法(divergence-free SPH，DFSPH)，利用壓強對密度進行 2 次調(diào)整， 使流體實現(xiàn)不可壓縮性。Band等人[28]將偏微分方程離散化引入IISPH的壓力求解方程中，用于求解流體的壓力加速度，可以有效地解決IISPH的壓力震蕩問題。Li等人[29]利用Bender等人[27]提出的DFSPH對流體的不可壓縮性進行求解N-S方程，通過求解密度約束來修正粒子的位置，以保證粒子密度在相對恒定狀態(tài)。

(2) 流固耦合仿真

由于SPH在處理邊界時具有較大的靈活性，該方法除了用于模擬單種流體的運動效果，也被廣泛用于流固耦合的模擬。如Shao等人[30]提出一種流固耦合算法，在固體邊界處給流體粒子施加一個相反的作用力，防止流體粒子穿透固體邊界，解決流固界面的穿透問題，如圖5所示，該方法能實現(xiàn)流-固雙向耦合、湍流效果，但由于在邊界施加這種排斥力并不屬于物理力，因此在模擬時常常發(fā)生流體粒子堆積現(xiàn)象。Gao等人[31]使用基于位置的動力學(xué)求解的位置約束和考慮固體邊界粒子對流體粒子的相對貢獻，有效解決了粒子堆積的問題。Gissler等人[32]將SPH用于模擬流體與剛體之間的耦合，在該界面中，流體壓力求解器每次迭代都會更新剛體顆粒的速度，與SPH相比可以支持更大的時間步長，減少了模擬流體所需的計算時間。

圖5 SPH粒子法流固耦合效果[30]

通過對上述的SPH及其演變方法分析發(fā)現(xiàn)，SPH在流體模擬中具有如下優(yōu)勢：具有較大的靈活性，能夠模擬不同類型流體，如粘性流體、不可壓流體與流體之間的交互效果； SPH只需要計算鄰域內(nèi)局部粒子的信息，不需要通過求解全局的N-S方程即可得到流體的運動，因此在小規(guī)模流體中具有較高的計算效率。

2.1.2 PBF法

雖然SPH方法具有簡單和易于實現(xiàn)的優(yōu)勢，但由于SPH對鄰居粒子的密度較為敏感，通過直接求解壓強來強制不可壓縮性是比較耗時的，因此在粒子較多時無法達到實時的效果。然而在游戲和電影等場景中，往往對實時性要求較高，而且大多數(shù)算法通常需要設(shè)置小的時間步長來維持系統(tǒng)的穩(wěn)定性，對于大規(guī)模的場景無法達到大時間步長的模擬，為此需要一種既能保證真實感又能達到實時的建模方法。

Müller等人[33]提出了一種基于牛頓第二運動定律的基于位置的動力學(xué)方法(position based dynamics，PBD)，基本思想是通過攜帶質(zhì)量、位置、速度等物理量的N個粒子和M個約束來表示一個物體，并通過歐拉積分預(yù)測粒子的位置，同時將內(nèi)力建模成一種約束，然后通過約束投影到最終的位置。由于采用顯式積分求解約束方程，因此可以有效地消除模擬過程中的不穩(wěn)定性，但該方法只用于布料的運動建模。隨后，Macklin等人[34]基于PBD的思想，提出一種基于位置的流體方法PBF，為了強制不可壓縮性，PBF通過使用狀態(tài)方程建立當(dāng)前粒子位置與鄰居粒子位置之間的密度約束函數(shù)：

(4)

其中，ρ0為初始密度，ρi則通過SPH的密度估算器[22]求得，因此避免了SPH先計算力再由力計算粒子位置所帶來的拉伸不穩(wěn)定問題，如圖6所示。Zhang等人[35]通過對輸入的模型進行采樣并自動生成控制粒子，用于控制流體生成的形狀，該方法在傳統(tǒng)的PBF上增加了流體形狀的可控性，能夠使工程師使用少量的控制參數(shù)去生成與目標形狀的流體效果，而且能夠滿足實時性的要求。

圖6 PBF粒子法仿真效果[34]

常規(guī)的PBF在模擬包含大量粒子的流體時，收斂速度慢的主要原因在于PBF為了求解流體的不可壓縮性時，需要通過迭代不斷地調(diào)整流體粒子的位置。為此，K?ster等人[36]提出一種自適應(yīng)PBF方法，根據(jù)每個粒子的水平細節(jié)信息自適應(yīng)地調(diào)整求解器的迭代次數(shù)，以保證穩(wěn)定性與運算效率的平衡。

PBF的優(yōu)勢在于其能夠?qū)崟r模擬流體運動，同時還能允許較大的時間步長，解決了傳統(tǒng)SPH方法不能實時模擬和時間步長的限制。

2.1.3 LBM法

LBM法是一種用于求解N-S方程的新型數(shù)值方法，其基本思想是使用隨機運動的微觀粒子表示流體中的分子或宏觀流體的一部分，并通過粒子流動和粒子碰撞完成動量的交換，LBM的基本模型如下[37,38]：

(5)

其中，式(5)的左邊表示LBM的streaming操作，右邊表示collision操作。雖然Eular網(wǎng)格法和Lagrange粒子在模擬流體的宏觀運動上具有很好的靈活性和準確性，但對于流體的微觀運動卻很難捕獲。相反，LBM方法則通過微觀分子來描述流體的宏觀運動狀態(tài)，避免了傳統(tǒng)數(shù)值方法的復(fù)雜性和精度問題，因此，近年來LBM方法在游戲、工程和物理仿真領(lǐng)域都有廣泛的應(yīng)用。在模擬大場景時，通常使用基于淺水方程的LBM，如Judice等人[39]利用LBM方法用于模擬游戲的流體動畫，并將LBM應(yīng)用于GPU，因此，該方法能夠?qū)崟r模擬大規(guī)模的流體動畫場景，而LBM的簡單高效特性在設(shè)計游戲場景中具有很好的優(yōu)勢。Bauza等人[40]提出一種基于LBM的實時流體交互方法，該方法結(jié)合淺水方程實現(xiàn)動畫與人的實時交互。然而在工程應(yīng)用中，傳統(tǒng)的LBM在模擬多相流效果時，存在密度比低、雷諾數(shù)小的問題，為此許多研究者改進了粒子概率分布函數(shù)f，如Li等人[41]與Wu等人[42]利用LBM模擬多相流效果，在蒸氣點附近利用平均插補的方法校正概率分布函數(shù)f，能夠模擬密度比大于500的液滴運動和液滴震蕩效果。

基于上述分析，LBM方法具有如下特點：

? 計算局部性好，易于擴展到GPU并行框架。

? LBM具有程序易實施、計算速度快的特點。

? 由于LBM基于微觀理論，因此能夠很好地描述多相流及復(fù)雜的邊界效果。

2.2 Eular網(wǎng)格法

基于Lagrange粒子法的流體模擬方法雖然具有靈活、方便的優(yōu)勢，并且SPH經(jīng)過近幾年的不斷改進，已經(jīng)在效果和效率上有很大的提升。與Lagrange粒子法不同，基于Eular網(wǎng)格的流體模擬方法的核心思想是將模擬的流體固定在網(wǎng)格上，然后通過計算流過網(wǎng)格節(jié)點的物理量來描述流體的運動，在整個仿真過程中，網(wǎng)格不隨流體的運動而運動。

2.2.1 高度場網(wǎng)格法

高度場網(wǎng)格法是基于淺水方程的一種流體模擬方法，其基本思想是忽略垂直方向的動力學(xué)問題，將流體產(chǎn)生的水波高度值存儲于二維的網(wǎng)格中，進而降低求解N-S方程和數(shù)值解的復(fù)雜度，淺水方程的基本形式如下：

(6)

(7)

其中，u為速度，h為水體的深度，s為地形高度。從式(6)和(7)可以看出，水波的流動只與2個參數(shù)有關(guān)，一個是高度信息，另一個則是流體的速度，并且這2個公式只在x方向上進行求解，去除了水深方向的信息，簡化了計算的復(fù)雜度，這也是基于高度場網(wǎng)格的主要優(yōu)勢。

Layton等人[43]提出一種用于產(chǎn)生水波的高效淺水方程方法，該方法通過非線性淺水方程來預(yù)測隨時間變化的水波高度，實時地模擬了水波效果。由于該方法使用Semi-Lagrange方法[44]求解淺水方程，因此可以達到無條件穩(wěn)定的大時間步長。隨后，Wang等人[45]將該方法用于求解固體與水波之間的雙向耦合，能夠模擬表面波紋、水滴、流體與固體的耦合效果，但無法模擬波浪運動過程中的破碎效果。為此，Thurey等人[46]提出一種基于淺水方程的實時破碎波浪仿真方法，首先用淺水方程模擬高度場網(wǎng)格，其次，為了產(chǎn)生波浪的破碎效果，通過檢測高度場中的陡峭波區(qū)域并用線段進行標記，然后在波峰處產(chǎn)生粒子，用于模擬波浪破碎的效果，如圖7所示。

圖7 高度場網(wǎng)格法仿真效果[46]

基于Thurey等人[46]的方法存在粒子偽影問題，Chentanez等人[47]直接在高度場網(wǎng)格產(chǎn)生噴霧粒子、飛濺粒子和泡沫粒子，并與網(wǎng)格進行信息傳遞，因此產(chǎn)生的破碎波現(xiàn)象更符合物理運動，但所描述的波長不得小于網(wǎng)格的大小，否則會導(dǎo)致生成的破碎波丟失部分細節(jié)。Nielsen等人[48]使用基于FFT的波浪模擬來生成高分辨率的紋理，用于替換渲染網(wǎng)格，產(chǎn)生的波浪細節(jié)效果比文獻[47]的更逼真，且實時性更好。

為了進一步增強大規(guī)模場景中的流體細節(jié)效果，Chentanez等人[49]通過耦合3D Eular網(wǎng)格、粒子法和高度場，將高度場網(wǎng)格耦合在3D Eular網(wǎng)格中用于產(chǎn)生流體運動過程中的水波，而在3D網(wǎng)格表面區(qū)域添加流體在飛濺過程中所需要的細節(jié)粒子，該方法可以快速地模擬大規(guī)模場景下的波浪破碎效果。邵緒強等人[50]基于淺水方程提出一種大規(guī)模實時的、穩(wěn)定的Eular網(wǎng)格法，對產(chǎn)生毛刺現(xiàn)象的高度場值進行平滑處理，而針對細節(jié)部分，其采用的是將高度場與SPH粒子進行結(jié)合，并通過隔點采樣的方法降低計算復(fù)雜度，相比文獻[49]的細節(jié)效果更加真實。

綜上所述，基于高度場網(wǎng)格法的流體模擬方法相比粒子法在實時性能上具有較大的提升優(yōu)勢，由于高度場網(wǎng)格是基于淺水方程[51]，同時去除了垂直方向的信息，計算效率較高，因此特別適合模擬大規(guī)模場景。

2.2.2 界面模型法

界面捕獲方法最常用的主要有流體體積法(volume of fluid，VOF)、水平集方法(level set，LS)以及相場法(phase field,PF)，VOF法通常使用流體體積函數(shù)F計算網(wǎng)格單元中流體的體積分數(shù)，然而在一個兩相流系統(tǒng)中，每個單元格都有如下3種可能：

當(dāng)F=0時，表示該網(wǎng)格單元被流體1所占據(jù)。

當(dāng)F=1時，表示該網(wǎng)格單元被流體2所占據(jù)。

當(dāng)0

并通過控制方程來跟蹤兩相界面變化：

(8)

其中u為速度，由于采用較為精確的數(shù)值方法計算網(wǎng)格單元的流體體積分數(shù)，因此可以保證質(zhì)量守恒的特性，但由于VOF定義的界面是不連續(xù)的，因此很難精確地計算界面的曲率和法線等物理信息。

(1) LS法

另一種界面捕獲方法LS法，將界面表示成連續(xù)光滑的零厚度的水平集函數(shù)φ，然后定義界面內(nèi)側(cè)的流體φ<0，在界面外側(cè)的流體φ>0，在界面層上的φ=0，LS模型如下:

(9)

從LS模型可以看出，LS方法定義的界面是光滑的零水平集函數(shù)，這種光滑特性可以準確地計算界面的曲率和法線等物理信息。然而在水平集方法中最大的挑戰(zhàn)是質(zhì)量守恒問題，由于隨著界面的演化，界面層的質(zhì)量不斷發(fā)生變化，為了保證界面在演化過程中質(zhì)量守恒，研究者們提出了許多守恒的水平集算法用于減少界面在演化過程中的質(zhì)量丟失。如Olsson等人[52]提出一種守恒的LS方法，該方法使用守恒的中間步對水平集函數(shù)進行平流，可以保證界面的輪廓和厚度是恒定的，不足的是該方法只用于2維的數(shù)值驗證。Losasso等人[53]提出一種SPH和LS雙向耦合的流體模擬方法，通過耦合SPH和LS法用于建模流體的流動，模擬的流體更符合實際的物理運動，但該方法在液體表面生成虛擬粒子，增加了系統(tǒng)開銷?；谖墨I[52]的工作，文獻[54,55]通過添加額外項修改重新初始化水平集函數(shù)來保持尖銳界面的質(zhì)量守恒。文獻[56]提出新的Heaviside函數(shù)求解水平集來保證質(zhì)量守恒，雖然通過重新初始化水平集函數(shù)可以有效解決質(zhì)量不守恒問題，但會降低算法的收斂速度。

在上述的改進算法中，大都關(guān)注解決LS方法的質(zhì)量守恒問題，而在常見的流體現(xiàn)象中，由于多相流在運動過程中，界面的厚度和形狀不斷發(fā)生變化，因此如何控制這些界面特性也是LS方法所需要解決的問題?；诖耍琒antos等人[57]利用區(qū)域水平集法跟蹤多相流的界面，并通過局部矯正來最小化總體積誤差。相比于原有的LS方法能夠處理更多不同流體間的界面，并且在界面層使用空氣或者水進行填充用于控制界面的厚度，這種方式會造成一定的人工偽影現(xiàn)象。Balczar等人[58]針對流體變形問題提出一種新的多標記水平集方法，利用多個水平集函數(shù)φi表示混合流體中每個子區(qū)域Ωi:Ωd={Ω1, Ω2,…, Ωn}，其中Ωn表示被劃分成區(qū)域個數(shù)，因此可以在相同的體積中控制不同的界面運動，同時將表面張力轉(zhuǎn)換為一種體積力，避免了傳統(tǒng)水平集方法的數(shù)值凝結(jié)問題。

(2) VOF-LS耦合方法

由于LS方法能夠很好地捕獲界面的光滑特性，VOF方法可以保證質(zhì)量守恒，將兩者耦合一起可以同時獲得LS方法和VOF方法的優(yōu)勢，為此，近些年來，許多研究者將兩者進行耦合用于處理更加復(fù)雜的界面演化現(xiàn)象。如文獻[59,60]提出一種耦合水平集合流體體積法(coupled level set and volume of fluid，CLSVOF)方案，該類方法主要通過水平集函數(shù)計算界面的表面張力和曲率等，CLSVOF方法能夠保證質(zhì)量守恒的情況下準確地計算界面的曲率和法向等信息。文獻[61,62]提出一種流體體積與水平集耦合法(volume-of-fluid-level-set，VOSET)，該方法直接利用幾何方法計算水平集界面，可以有效解決界面因震蕩引起的不連續(xù)問題。由于具有質(zhì)量守恒和計算精確的優(yōu)點，因此被廣泛用于捕獲不同多相流界面。

水平集方法及其改進的方案雖然能保證質(zhì)量守恒和準確計算曲率等信息，但是該方法仍存在以下2個局限性：

? 由于LS方法定義多相流界面是光滑的零水平集函數(shù)，界面的厚度為零，通常只被用于兩相流的模擬，且該方法難以表現(xiàn)出多相流混合過程中的擴散特性。

? 為了更加精確地表示界面的零厚度特性，LS方法通常通過網(wǎng)格細化來表示界面的細節(jié)特征度，因此增加了計算時間和復(fù)雜度。

(3) 相場法

為了解決水平集法界面擴散特性的局限性，Shen等人[63]提出一種用于模擬多相流擴散運動的相場模型，由于界面的厚度是一層薄而非零的過渡區(qū)域，因此在模擬過程中，這種擴散界面可以捕獲更多的流體運動細節(jié)。描述彌散界面的方法主要有2種：一種是基于守恒場的相場方法(cahn-hilliadrd,CH)，另一種是基于非守恒場的相場法(allen-cahn,AC)。基于CH相場模型一般用擴散動力學(xué)表示相場濃度隨空間和時間的變化，通常用于描述相與相之間演化過程中的守恒場變量，如兩相流中的混合流體中的濃度場(也稱為質(zhì)量分數(shù))等，起初用于材料領(lǐng)域的枝晶、共晶和包晶的生長[64]，后被廣泛用于多相流領(lǐng)域。CH相場模型的控制方程為

(10)

(11)

式(10)左邊第1項表示相場濃度隨時間的變化，第2項表示多相流中平流對相場φ的影響，第3項表示擴散項。其中u為速度場；M是有關(guān)相場的遷移率，通常只與相場濃度相關(guān)，如文獻[65]的M=1-φ2，或者為常數(shù)，即M=1，關(guān)于遷移的設(shè)置可以參考文獻[66]；μ為化學(xué)勢， 當(dāng)達到平衡時，化學(xué)勢必須滿足μ=0。相場模型用于控制多相流中交界面的創(chuàng)建和演化，而不可壓流體的流動則由N-S方程控制，將相場方程和N-S方程耦合一起得到著名的“模型H”[67]。

He等人[68]基于相場理論提出一種用于界面控制的流體混合方法，該方法主要包含2部分：第1部分是利用網(wǎng)格法求解N-S方程，獲得流體平流的速度；第2部分在擴散方程中引入擴散因子和銳化因子，j=αjs+βjd，其中α為銳化因子，β為擴散因子，用于控制界面的彌漫效果。如圖8所示，生成的流體界面為彌漫界面，與LS方法相比，這種方法的優(yōu)點是多相流混合過程中界面的厚度是可控的，擴散界面的細節(jié)比較突出。Li等人[69]在原有的CH方程中添加矯正項使得兩相之間的界面滿足光滑、均勻分布的特性，該方法簡化了復(fù)雜的網(wǎng)格拓撲，并且捕獲的界面效果要優(yōu)于原始的CH相場模型。文獻[70,71]提出一種用于相場模型的全自適應(yīng)能量穩(wěn)定的數(shù)值模擬方法，文獻[70]通過引用自適應(yīng)網(wǎng)格對擴散界面區(qū)域進行細化，可以使得擴散界面變得非常薄，同時引入改進的歐拉方法用于提高二階空間有限差分方法的精度， 相比于基于常規(guī)網(wǎng)格的相場法，精度更高，但由于界面處網(wǎng)格是自適應(yīng)細化的，因此當(dāng)界面厚度較小時，需要較大的計算開銷。

圖8 相場法仿真效果[68]

CH相場模型雖然在描述流體擴散的彌散界面具有較好的表達優(yōu)勢，但由于CH模型涉及4階空間導(dǎo)數(shù)，因此需要昂貴的計算資源，相比CH相場模型，AC相場模型的計算只涉及2階的空間導(dǎo)數(shù)，因此實現(xiàn)較為容易，AC相場模型通常用于描述非保守場變量的物理量， 則2維情況下非守恒的AC相場模型可表示為

(12)

由于傳統(tǒng)的AC相場模型不能保持質(zhì)量守恒，因此許多研究者通過引入矯正因子[72-74]來解決AC的不守恒性，守恒的AC方程可寫成如下公式：

(13)

其中λ(t)為引入的矯正因子，雖然文獻[72-74]引入的矯正因子可以保證相場在演化過程中保持質(zhì)量守恒，但是引入的矯正因子會使得模型過于依賴參數(shù)的設(shè)置，同時由時間依賴的對流-擴散方程存在高梯度或者不連續(xù)的缺點，會造成界面的局部震蕩問題。針對這一問題，Bazilevs等人[75]利用相場方程的自身殘差對原有的模型進行修改，同時將不連續(xù)捕獲YZβ-DC算子添加到有限元公式中，使標量保持在預(yù)期的物理范圍內(nèi)，解決了傳統(tǒng)相場模型對參數(shù)的依賴性。

求解多相流動方程除了在不守恒的AC相場方程中引入矯正因子外，還需通過細化網(wǎng)格來提高彌散界面的精度，如文獻[76,77]通過對均勻網(wǎng)格進一步細化用于求解多相流的彌散界面，相比非細化網(wǎng)格的精度要高。

通過對比LS模型和CH相場模型可以發(fā)現(xiàn)，2個模型的定義很相似，其最主要的區(qū)別在于LS法在界面處的φ=0，而CH模型-1<φ<1，同時CH模型的控制方程比LS法多了一項擴散因子。在描述銳化界面特性上，LS法相比于CH模型和AC模型使用更廣泛，LS法在計算曲率和法向上具有較好的優(yōu)勢，同時LS法描述的界面是零厚度光滑函數(shù)，界面尖銳性可以直接通過求解LS模型即可。

2.3 基于Lagrange粒子-Eular網(wǎng)格的混合方法

基于網(wǎng)格法的流體模擬由于不需要考慮粒子與鄰居粒子之間的信息傳遞，因此計算效率相對較高。然而由于網(wǎng)格法是將問題域固定在網(wǎng)格上，因此網(wǎng)格的大小決定著流體細節(jié)質(zhì)量。對于粗網(wǎng)格模擬的效果往往較為粗糙，細節(jié)表現(xiàn)不夠，但可以達到實時模擬，而細網(wǎng)格模擬的細節(jié)效果更為逼真，但需要較大的計算開銷。相反，SPH則通過將流體域建模成離散的粒子，更擅長流體的細節(jié)模擬，尤其在邊界處的細節(jié)，但當(dāng)粒子數(shù)增多時，計算較為耗時。由于網(wǎng)格法和粒子法都各有不同的優(yōu)勢和不足，因此，將這2種方法結(jié)合起來，可以利用兩者的優(yōu)勢，使模擬的流體效果更加符合物理的運動。Lagrange粒子-Eular網(wǎng)格混合法的基本思想是在網(wǎng)格中分布一定的攜帶物理量的流體粒子，首先將粒子的速度插值到網(wǎng)格邊中心進行計算，在進行平流時再將網(wǎng)格速度插值到粒子上，完成速度的更新操作。

常見的Lagrange粒子/Eular網(wǎng)格混合法主要包含PIC法和FLIP法。基于標記網(wǎng)格(marker-and-cell， MAC)[81]的插值思想，兩者都是將粒子信息插值到網(wǎng)格上，這種插值方法的過程如圖9所示。具體步驟：(1) 初始化粒子信息；(2) 將速度插值到網(wǎng)格上，并在網(wǎng)格上對壓力進行求解；(3) 將速度插回粒子，并更新粒子速度； (4) 更新粒子的位置。

圖9 粒子/網(wǎng)格插值過程

2.3.1 PIC法

PIC方法是將粒子與網(wǎng)格結(jié)合的一種數(shù)值方法，原始的PIC方法[82,83]將流體表示成攜帶質(zhì)量和位置信息的粒子，然后直接在網(wǎng)格上計算速度，將下一時間步的網(wǎng)格速度插值回當(dāng)前粒子的速度，在穩(wěn)定性上具有較好的優(yōu)勢，但是由于信息在粒子與網(wǎng)格之間傳輸會存在丟失信息問題，因此產(chǎn)生高耗散的現(xiàn)象。

Edwards等人[84]為了提高PIC的數(shù)值精度，提出一種高階精確的PIC方法，先是利用高階移動最小二乘法近似網(wǎng)格上的粒子，然后通過有限差分法[85]來估算網(wǎng)格上的物理量，最后將這些物理量從網(wǎng)格插回粒子，并通過在粒子上添加正則化項來消除PIC的噪聲，該方法具有較低的數(shù)值耗散。雖然文獻[84]能達到較高數(shù)值精度，但還存在一定的數(shù)值耗散問題，為此，Jiang等人[86]提出一種仿射式的PIC法(affine particle-in-cell method,APIC)，將速度和位移的增量從網(wǎng)格插值到粒子，并過濾從網(wǎng)格到粒子的傳輸信息，該方法不但能消除原始PIC的人工耗散問題，而且還能克服FLIP的不穩(wěn)定性問題。Hammerquist等人[87]通過拓展原有的PIC方法，提出一種擴展的PIC方法(extended PIC，XPIC)，該方法在速度中添加平滑項來降低FLIP的噪聲，并通過最小化新速度與平滑速度之間的誤差來解決PIC的數(shù)值耗散問題。

2.3.2 FLIP法

利用PIC方法能夠獲得較好的穩(wěn)定性能，但存在較高的數(shù)值耗散問題，為此，Brackbill等人[88]提出一種用于低耗散的FLIP法，用于解決PIC數(shù)值耗散問題。Zhu等人[89]將FLIP方法引入計算機圖形學(xué)，與PIC方法不同，F(xiàn)LIP是將上一時刻和當(dāng)前時刻的網(wǎng)格速度的差值疊加到當(dāng)前粒子的速度上，因此FLIP相比于PIC模擬的浪花效果更加明顯，同時能避免PIC方法的數(shù)值耗散問題。

Boyd等人[90]拓展了文獻[89]的FLIP方法，提出一種用于模擬兩相流的MultiFLIP方法，首先在常規(guī)網(wǎng)格上求解N-S方程獲得平流速度，然后利用粒子水平集方法建模不同相之間的界面，最后關(guān)聯(lián)相位與MultiFLIP粒子，由于FLIP需要在流體域上保持粒子，因此用粒子重建的液體表面更為自然。Yang等人[91]通過自適應(yīng)FLIP方法計算每個流體粒子的密度、壓強更新粒子的速度和位置，并采用線性混合PIC和FLIP的速度，最后通過深度圖直接在網(wǎng)格周圍尋找彌散粒子，增強了復(fù)雜流體間相互作用的細節(jié)。Ferstl等人[92]基于文獻[91]的重采樣思想提出一種窄邊帶的FLIP方法(narrow band-FLIP，NB-FLIP)，通過對表面固定窄帶區(qū)域不斷重采樣粒子，以控制表面窄帶的寬度，該方法不但減少了模擬域的FLIP粒子，還解決了PIC法的耗散與FLIP法的不穩(wěn)定性問題。鄒長軍等人[93]利用NB-FLIP方法模擬多相流現(xiàn)象，模擬單相流時與文獻[86]得到的仿真效果相似，但是在處理兩種流體界面時，在不同相交界處的界面明顯分離現(xiàn)象。

綜上分析，由于PIC方法在穩(wěn)定性能上優(yōu)于FLIP方法，因此較適合用于處理大形變問題，但會存在較高的數(shù)值耗散問題，其主要通過2種方式來降低數(shù)值耗散的影響:(1)引入正則化項；(2)過濾粒子與網(wǎng)格之間傳遞的信息。相比之下，由于FLIP在粒子與網(wǎng)格之間傳遞的是2個時間步速度的變化值，采用疊加的方式更新速度，而不是用下一時間步的速度覆蓋當(dāng)前粒子的速度，因此不需要考慮PIC的數(shù)值問題。

3 對比分析

基于物理的流體模擬方法的分析與比較結(jié)果見表1，對比指標主要包含算法的分類、算法名稱、發(fā)表年份、技術(shù)特點、效率以及最大的仿真場景等。

從表1可以看出，根據(jù)算法的仿真場景和技術(shù)特點，基于粒子的方法從傳統(tǒng)的SPH方法演化至IISPH方法，不可壓縮性都比較好，然而在效率上，由于受到粒子數(shù)的限制，該類型算法在模擬小規(guī)模流體流動場景上具有較好的仿真效果，且大都以離線渲染為主，無法達到實時模擬。而與Lagrange粒子法相比，Eular網(wǎng)格法由于不需要對粒子進行計算和存儲，可用于解決流體仿真的實時性問題，并且大都用于大規(guī)模海浪場景的構(gòu)建。Level-Set、CLSVOF、VOSET這3種類型方法主要用于建模兩相或多相流尖銳性界面的數(shù)值建模，該類算法模擬的多相流界面具有足夠的光滑特性，且能夠準確計算表面張力；而CH相場法和AC相場法主要用于建模多相流的彌散界面，由于該類方法是基于菲克擴散定理的，因此在流體流動過程中能夠準確捕獲到精細的界面細節(jié)。

表1 流體仿真算法對比

4 存在的問題及未來展望

綜上分析，雖然近年來基于物理的不可壓縮流體模擬仿真領(lǐng)域取得了突出的成果，但在實時性、復(fù)雜流體耦合等方面仍存在一些問題，存在的不足及未來可能的研究方向如下。

4.1 大規(guī)模場景的實時繪制

在目前的流體模擬算法中，大多算法都能實時模擬小規(guī)模流體，然而對大規(guī)模場景的流體模擬卻很難達到實時，雖然使用高度場網(wǎng)格法模擬在實時性上具有一定的優(yōu)勢，但其模擬的大都為光滑流體表面，在細節(jié)效果上仍很難達到逼真的程度，同時對于航海類的大規(guī)模場景，往往需要處理的是更多的流體粒子或者更大的網(wǎng)格。因此，如何將應(yīng)用于小規(guī)模流體的仿真算法擴展于大規(guī)模場景，并能實時處理流體與固體的耦合關(guān)系，對于計算機仿真具有重要的研究意義。

4.2 增強邊界細節(jié)模擬

本文介紹的流體模擬方法中，基于粒子法的流體模擬為了表現(xiàn)更加精細的細節(jié)效果，往往需要增加更多的流體粒子，但這會造成較大時間和計算的開銷，同時，當(dāng)邊界處的粒子較為稀疏時，會導(dǎo)致拉伸不穩(wěn)定。而對于網(wǎng)格法，很難處理拓撲變形嚴重的邊界，同樣，如果通過細化網(wǎng)格來增加邊界的細節(jié)效果，會損失一定的實時性能，因此如何保證流體模擬過程中的實時性與細節(jié)效果之間的平衡是未來的面臨的一項挑戰(zhàn)。

4.3 多相流界面控制

多相流模擬可以為虛擬場景增加更豐富的細節(jié)特征，然而現(xiàn)有的多相流仿真算法為了簡化模擬過程，大多忽略溫度、氣流等外界因素對多相流界面演化的影響，而這些外界因素往往對多相流的效果具有很大影響，同時常規(guī)的多相流算法難以控制界面的運動以及厚度等，造成界面細節(jié)丟失的現(xiàn)象，因此在模擬多相流時包含一定的外界因素和控制多相流界面的生成會給多相流模擬仿真領(lǐng)域帶來新的研究方向。

4.4 多種類型方法的耦合

目前針對多種方法融合的模型較少，大多為粒子與網(wǎng)格之間的融合，如文獻[46]在高度場網(wǎng)格的基礎(chǔ)上添加粒子，文獻[51]提出的仿射FLIP方法，文獻[52]的擴展FLIP方法，在效果上具有很大的提升。但這些方法并沒有解決不同方法耦合銜接問題，尤其在2維網(wǎng)格與3維網(wǎng)格的耦合，如果處理不好可能會導(dǎo)致如文獻[46]的人工偽影問題，因此未來仍需要提出一些新的耦合模型來解決現(xiàn)有耦合模型銜接問題。

5 結(jié) 論

如今，基于物理的流體模擬已經(jīng)成為增強現(xiàn)實領(lǐng)域的主流方法，尤其是基于粒子法和網(wǎng)格法的仿真方法，為實時、逼真的流體動畫場景提供了技術(shù)支持。近年來這些方法在特定領(lǐng)域問題上已經(jīng)取得了很好的研究成果，如SPH方法在小規(guī)模流體仿真領(lǐng)域上的應(yīng)用，高度場網(wǎng)格法解決大規(guī)模實時性問題等。本文在分析基于物理的不可壓縮流體模擬方法存在的問題的基礎(chǔ)上，對現(xiàn)有的流體模擬方法進行分類和綜述，將有助于理解和研究流體仿真方法中的關(guān)鍵技術(shù)。