趙 靜,許立群,安泳鋼,王洪玉,唐 勇,*

(1.燕山大學(xué) 河北省計算機(jī)虛擬技術(shù)與系統(tǒng)集成重點(diǎn)實(shí)驗(yàn)室,河北 秦皇島 066004;2.燕山大學(xué) 信息科學(xué)與工程學(xué)院,河北 秦皇島 066004;3.東北大學(xué)秦皇島分校,河北 秦皇島 066004)

0 引言

在虛擬世界中模擬自然現(xiàn)象是計算機(jī)圖形學(xué)中一個重要研究領(lǐng)域,而流固交互現(xiàn)象因其流體拓?fù)浣Y(jié)構(gòu)變化劇烈、非線性固體模型種類多樣、流固接觸面邊界條件復(fù)雜以及整體方案需隱式求解等特點(diǎn),使之成為圖形學(xué)中最充滿挑戰(zhàn)的一個分支,這也吸引了國內(nèi)外很多研究者對流固交互技術(shù)展開深入研究。目前流固交互技術(shù)已經(jīng)廣泛應(yīng)用于力學(xué)分析、虛擬現(xiàn)實(shí)、游戲影視等領(lǐng)域。由于虛擬空間中的流體運(yùn)動通過求解偏微分方程獲得,固體的形變效果由其應(yīng)力應(yīng)變關(guān)系決定,因此目前研究者主要從流體和固體的表示方法和離散化方案、拓?fù)浣Y(jié)構(gòu)、邊界條件以及加速策略等幾個角度對流固交互過程中涉及的技術(shù)進(jìn)行探索和研究。

如今,在基于物理的真實(shí)場景模擬中對流固交互技術(shù)的需求正越來越強(qiáng)烈,如模擬填充空氣的橡膠輪胎、液壓過程、空中飛行的飛機(jī)或者水中漂浮的船等。除了自然環(huán)境的真實(shí)感重現(xiàn),在虛擬手術(shù)[1-2]、數(shù)字裝配[3]以及軟體機(jī)器人[4]等領(lǐng)域也都應(yīng)用了大量的流固交互相關(guān)內(nèi)容。

由于流固交互在實(shí)現(xiàn)過程中會使用到不同的表示方法,這使得同步處理交互過程中的邊界條件充滿挑戰(zhàn)。最初,通常使用分區(qū)方案處理流固交互過程,即在流體和固體間單獨(dú)地迭代計算,并使用其中一個材料的屬性定義另一個材料的邊界條件[5-6]。然而這種簡單直接的計算方法很容易出現(xiàn)數(shù)值穩(wěn)定性問題,需要使用較小的時間步長來解決。因此,很多研究者開始探索一種整體方案,Zarifi等[7]以及Akbay等[8]先后證明通過整體方案可以使用更大的計算時間步長并獲得更好的系統(tǒng)穩(wěn)定性。不同方案的選擇受流體和固體所使用的表示方法影響。由于固體的靜止?fàn)顟B(tài)形狀很容易定義,因此通常使用拉格朗日網(wǎng)格表示法對其狀態(tài)進(jìn)行更新[9];相反,流體由于沒有靜止?fàn)顟B(tài)形狀,因此通常使用歐拉網(wǎng)格表示法進(jìn)行模擬[10]。這使得在流體和固體交界面以相同的法向速度移動需要額外的處理方案,Batty等[11]通過將力離散到歐拉網(wǎng)格表面來近似計算不規(guī)則表面,而Zarifi等[7]則通過更加復(fù)雜的cut-cell方法,以使邊界結(jié)果更加準(zhǔn)確。為了避免不同表示方法帶來的耦合問題,很多研究者開始探索使用相同的方法來表示流體和固體[12-13]。另一方面,非線性固體模型對流固交互的真實(shí)性也有很大的影響。通常很多研究者會使用線性模型近似表示非線性固體模型,并通過迭代收斂得到非線性模型的解[14],但這也會占用大量的計算性能。而混合粒子-網(wǎng)格方法的提出[15-17]很好地緩解了這一問題。

可以看出,目前在基于物理的流固交互技術(shù)方面的研究主要從流體和固體所使用的表示方法出發(fā),依據(jù)分區(qū)方案和整體方案的設(shè)計特點(diǎn),結(jié)合真實(shí)性和效率問題優(yōu)化流體運(yùn)動方程的解算過程并選擇合適的固體模型;其次,結(jié)合不同表示方法的特點(diǎn)選用特定的離散化方案,針對性地改進(jìn)當(dāng)前表示方法下實(shí)現(xiàn)效果;最后,從流固交互數(shù)值計算的加速策略中進(jìn)一步提升系統(tǒng)的魯棒性和效率。

1 流體運(yùn)動方程

在基于物理的流固交互中,流體的運(yùn)動方程通過一個偏微分方程表示,即Navier-Stokes方程(NS方程)[10]

其中,ρ為流體密度,u(x,t)是空間中流體速度,ut是的縮寫,p表示壓力,μ為運(yùn)動粘度,f表示外力,如重力,是流體密度的物質(zhì)導(dǎo)數(shù)。通常數(shù)學(xué)界傾Δ向于使Δ用ΔΔ表示拉普拉斯算子,而工程界使用2或·表示。NS方程從數(shù)學(xué)含義上表示牛頓流體的動量守恒(式(1))和質(zhì)量守恒(式(2))。在特定情況下,也會增加關(guān)于氣壓、溫度以及密度相關(guān)的狀態(tài)方程[18]實(shí)現(xiàn)更豐富的流體效果。NS方程中主要考慮兩點(diǎn):1)應(yīng)用到流體運(yùn)動的牛頓第二定理;2)假設(shè)流體中的應(yīng)力為擴(kuò)散黏性項(xiàng)和壓力項(xiàng)兩部分。通常,對于非黏性流體,會將NS方程中的黏性項(xiàng)μΔu移除,并稱之為歐拉方程。而對于不可壓縮流體,則,即式(2)可以簡化為。

2 固體模型

隨著基于物理流固交互技術(shù)的不斷發(fā)展以及固體本構(gòu)模型的引入,在交互過程中使用的固體模型也越來越真實(shí)和豐富。根據(jù)固體的應(yīng)力應(yīng)變關(guān)系,大體上可以將其分為兩大類:理想的不能形變的剛體以及可以進(jìn)行彈性和塑性形變的軟體或形變體。

2.1 剛體

剛體是流固交互中最常見的固體模型。通常自然界中的任何固體都會發(fā)生一定的形變,因此剛體是流體交互中的一種理想模型。在這個理想的限定下,剛體所有點(diǎn)的位置可以通過6個自由度來表示。首先,定義剛體的質(zhì)心為空間中的一點(diǎn),獲得表示位置的3個自由度信息。其次,以剛體的質(zhì)心為定位點(diǎn),就可以得到剛體內(nèi)任意一點(diǎn)基于質(zhì)心的3個自由度信息。利用剛體中所有點(diǎn)彼此之間保持固定距離的約束關(guān)系,就可以得出剛體運(yùn)動的規(guī)律。

Baraff等[19]介紹了完整的剛體模擬方法。通過將剛體的質(zhì)心視為坐標(biāo)原點(diǎn),以此方便地將坐標(biāo)系固定到剛體上進(jìn)行運(yùn)動分析。通常,剛體上的點(diǎn)相對于質(zhì)點(diǎn)的坐標(biāo)系稱為物體空間,而將定義所有對象的公共坐標(biāo)系稱為世界空間。剛體的物體空間和世界空間如圖1所示。

圖1 物體空間(左)和世界空間(右)Fig.1 Object space(left)and world space(right)

在時刻t,當(dāng)剛體中的某一點(diǎn)在物體空間r0處,剛體質(zhì)心在世界空間x(t),而剛體的方向可以通過一個關(guān)于質(zhì)心的旋轉(zhuǎn)矩陣R(t)描述。那么,該點(diǎn)的世界空間坐標(biāo)可以表示為

2.2 形變體

2.2.1 線性彈性模型

在基于物理的流固交互場景中,線性模型是常見的非線性彈性理論的簡化,描述了物體如何形變以及如何將外部力轉(zhuǎn)化為內(nèi)部應(yīng)力。線性彈性模型的基礎(chǔ)線性化假設(shè)是:1)無限小的應(yīng)變或小的形變,2)應(yīng)變和應(yīng)力之間是一種線性關(guān)系,3)僅適用于不會產(chǎn)生屈服的應(yīng)力狀態(tài)。

線性彈性模型的邊界條件通過3個關(guān)于線動量平衡的張量偏微分方程和6個最小應(yīng)變-位移關(guān)系控制。整個偏微分方程系統(tǒng)由一組線性算術(shù)本構(gòu)關(guān)系組成。線性彈性模型的控制方程[20]為

式(3)為由牛頓第二定律推導(dǎo)出的形變體運(yùn)動方程,式(4)為應(yīng)變-位移方程,式(5)為由胡克定理推導(dǎo)出的形變體本構(gòu)方程,對于彈性材料,胡克定律表示材料的運(yùn)動行為并關(guān)聯(lián)應(yīng)力和應(yīng)變的關(guān)系。其中,σ是柯西應(yīng)力張量,ε是無限小應(yīng)變張量,d為位移向量,C表示剛度Δ張量,Fbody是單位體積所受的力,ρ為質(zhì)量密度,表示Nabla算子,d¨表示位移關(guān)于時間的二階導(dǎo)數(shù)。

2.2.2 超彈性模型

隨著流固交互研究的不斷發(fā)展,研究者們發(fā)現(xiàn)線性彈性模型并不能準(zhǔn)確地表示彈性材料的真實(shí)行為。對于可以產(chǎn)生較大形變的材料,其應(yīng)力-應(yīng)變關(guān)系可以定義為非線性的、各向同性的、或不可壓縮的且通常依賴于材料特有的應(yīng)變率。而超彈性模型提供一種表示材料應(yīng)力-應(yīng)變關(guān)系的方法[21]。

最簡單的超彈性模型是Saint Venant-Kirchhoff模型,它僅僅將線性模型幾何擴(kuò)展到非線性領(lǐng)域。該模型存在一般形式和各向同性形式,分別如下:

其中,S為 Piola-Kirchhoff第二應(yīng)力張量,E為Lagrangian-Green應(yīng)力張量,C表示剛度張量,λ和μ為 Lamé常數(shù),I為單位矩陣。Saint Venant-Kirchhoff模型的應(yīng)變-能量密度函數(shù)為

其他常見的非線性模型還有Neo-Hookean模型[22]和Fixed Corotated模型[23]。

3 表示方法和空間離散化

3.1 歐拉和拉格朗日表示法

在流固交互中,不同的表示方法對流固交互的實(shí)現(xiàn)效果有直接的影響。在經(jīng)典場論中[24],拉格朗日表示法是指觀察者跟隨在空間中隨時間移動的單獨(dú)的流體包來觀察流體運(yùn)動的一種方式。繪制流體包隨時間運(yùn)動的位置可以獲得該包的運(yùn)動蹤跡。拉格朗日規(guī)范可以視為坐在船上并隨河水漂流。歐拉表示法是指聚焦在流體隨時間流動的空間中某一指定位置上來觀察流體運(yùn)動的方法。歐拉規(guī)范可以看作在河岸邊看河水流過某一固定的位置。

拉格朗日和歐拉表示法有時可以不嚴(yán)謹(jǐn)?shù)乇硎鰹槔窭嗜蘸蜌W拉參考系。然而,通常拉格朗日或歐拉表示法可以應(yīng)用于任何觀察者參考系中,并且可以用于基于任何坐標(biāo)系統(tǒng)的參考系中。

在計算流體動力學(xué)或計算機(jī)圖形學(xué)中,歐拉模擬使用固定的網(wǎng)格,而拉格朗日方式以隨速度場移動的模擬節(jié)點(diǎn)為特征。

在歐拉規(guī)范中,流場用位置x和時間t的函數(shù)表示。例如,流場速度可以通過函數(shù)u(x,t)表示。同樣,流場的壓力和密度可以分別表示為p(x,t)和ρ(x,t)。通過這些函數(shù)的行為就可以確定整個流場。在拉格朗日規(guī)范中,將整個流體分割成一個個單獨(dú)的流體包,這些流體包隨時間移動,且由某個與時間無關(guān)的向量場x0標(biāo)簽化。通常,x0表示初始時刻t0流體包的質(zhì)心位置,并以這種方式來分析物體隨時間推進(jìn)可能的變化。在拉格朗日規(guī)范描述中,流場通過在給定時刻t流體包x0的位置函數(shù)X(x0,t)表示。通過追蹤這些隨時間移動的流體包的特性和運(yùn)動可以確定整個流場的特性。

以圖2為例,在歐拉規(guī)范中,在煙囪上部位置O放置一個溫度計,記錄與時間t有關(guān)的函數(shù)T。當(dāng)不同的煙霧粒子流經(jīng)點(diǎn)O的時候,溫度計的值T會發(fā)生變化。給定O點(diǎn)溫度為T(x0,y0,z0,t),那么整個場中的溫度可以表示為T(x,y,z,t),其中x,y,z與時間t無關(guān)。然而,在拉格朗日規(guī)范中,溫度計則被依附在煙霧粒子A上,有TA=TA(t)。如果知道每個粒子隨時間變化的位置函數(shù),就可以將拉格朗日規(guī)范轉(zhuǎn)化為歐拉規(guī)范。

圖2 流場表示方法示意圖Fig.2 Illustration of flow field specification

根據(jù)上面的陳述,這兩種規(guī)范有如下關(guān)系:

等式兩邊描述了在時刻t,標(biāo)簽為x0的流體包或粒子的速度。

3.2 空間離散化

NS方程是一種偏微分方程,除了一些較為簡單的線性情況外,大多數(shù)情況下很難得到解析解。由于計算機(jī)的內(nèi)存空間有限,所以對這些問題的離散化是必要的,進(jìn)而可以通過數(shù)值方法獲得偏微分方程的近似解。常見的空間離散化方法有:有限差分法、有限元法、有限體積法、譜元法、格子玻爾茲曼方法、邊界元方法等。本文僅對常用的有限差分方法和有限元法進(jìn)行簡單的介紹。

由于有限差分方法和有限元方法相對容易實(shí)現(xiàn),目前有限差分方法[25]和有限元方法[26]已經(jīng)是求解偏微分方程數(shù)值解的最常見的方法。

4 時間離散化

在流固交互計算中除了進(jìn)行空間離散化,還必須對時間進(jìn)行離散化。為了形成動畫,通常需要以固定的頻率按時間序列計算下一時刻的輸出結(jié)果。盡管可以選擇不同尺度的時間步長,但是通常采用一個固定的時間步長來簡化問題并避免偽影。雖然,目前在實(shí)際的模擬中每一幀仍需要計算多個時間步,但是大多數(shù)最新的動畫系統(tǒng)中已經(jīng)盡可能地減少每幀的時間步數(shù)。

4.1 顯式時間積分

顯式時間積分是指在更新時刻t+Δt的狀態(tài)的時候,僅使用時刻t的量計算。常見的歐拉積分器

就是一種顯式積分器,式(8)中,通過時刻t的位置x(t)和速度u(x,t)顯式地計算時刻t+Δt的位置x(t+Δt)。還有很多類似的積分方法,但具有不同的特性。在研究此類積分器的時候,重要的是,對于特定應(yīng)用而言重要的屬性,可能對于計算機(jī)動畫而言并不重要,并且看起來非常相似的積分技巧在實(shí)際中可能會有非常不同的表現(xiàn)。目前常見的顯式積分器包括:梯形法積分、中點(diǎn)法積分以及辛歐拉積分器。

4.2 隱式時間積分

顯式時間積分通過系統(tǒng)中當(dāng)前時刻的狀態(tài)計算后一時刻的系統(tǒng)狀態(tài),而隱式時間積分通過解算包含當(dāng)前時刻系統(tǒng)狀態(tài)和后一時刻系統(tǒng)狀態(tài)的方程,以獲得最終的解。隱式方法比顯式方法實(shí)現(xiàn)難度更大,但是在解決某些剛性問題的時候,顯式時間積分需要很小的時間步長來維持?jǐn)?shù)值穩(wěn)定性。為了解決這類問題并獲得更高的精度,通過隱式方法就可以在較大時間步長的條件下以較少的計算量實(shí)現(xiàn)。

由于隱式方法不能直接處理所有類型的差分算子,因此在使用隱式方法的時候需要利用算子分割方法,這是隱式方法的主要局限性。

5 流固交互方法

根據(jù)交互場景中流體和固體所使用表示方法的差別,流固交互方法大致可以分為四類:歐拉流體和拉格朗日固體的交互、拉格朗日流體和拉格朗日固體的交互、歐拉流體和歐拉固體的交互以及混合的粒子網(wǎng)格方法,每種交互類型中根據(jù)具體實(shí)現(xiàn)方法的不同,又可以進(jìn)一步劃分。

5.1 歐拉流體和拉格朗日固體交互

5.1.1 MAC法

歐拉流體方法最早應(yīng)用于計算流體力學(xué)領(lǐng)域,1965年,Harlow等[27]使用 MAC(Marker and Cell)方法和有限差分形式的NS方程,進(jìn)行了不可壓縮流體和剛體的交互模擬,并實(shí)現(xiàn)了滑動和非滑動邊界條件。隨后,Foster等[28]首次將MAC交錯網(wǎng)格引入計算機(jī)圖形學(xué),通過解算二維或三維空間的NS方程,實(shí)現(xiàn)一種可控的復(fù)雜行為流固交互技術(shù),其中固體通過拉格朗日網(wǎng)格法表示。但是,該方法中流體只能實(shí)現(xiàn)與網(wǎng)格對齊的不可移動的固體交互,處理曲線邊界會出現(xiàn)體素化偽影。為了解決參數(shù)化曲線邊界偽影問題和移動邊界條件問題,他們[29]又使用半拉格朗日方法和水平集方法,分別計算法向和切向速度分量,通過速度場獲取隱式表面。但是這個方法只能處理在一個網(wǎng)格節(jié)點(diǎn)中存在一個多邊形面的情況,當(dāng)網(wǎng)格分辨率較低的時候只能按照多個邊的平均法線處理。Rasmussen等[30]提出使用粒子水平集的方法,結(jié)合移動網(wǎng)格窗體技術(shù),實(shí)現(xiàn)了更加豐富的流體飛濺效果。但是流固交互邊界的體素化偽影問題依舊存在。

5.1.2 ALE法

1974年,Hirt等[31]首次提出了ALE(Arbitrary Lagrangian-Eulerian)方法。Feldman等[32]使用基于速度的非結(jié)構(gòu)化四面體網(wǎng)格表示歐拉流體,利用基于半拉格朗日的ALE方法,進(jìn)行表面追溯。但是流體形變過程中的網(wǎng)格重建操作會消耗大量的計算機(jī)性能,同時數(shù)值耗散現(xiàn)象明顯。Klingner等[33]利用非結(jié)構(gòu)化的自適應(yīng)網(wǎng)格改進(jìn)上述數(shù)值耗散問題,但是網(wǎng)格重建問題依然存在。

5.1.3 Rigid Fluid法

Carlson等[34]通過拉格朗日乘子定義固體內(nèi)部速度場并約束固體運(yùn)動,從而將固體視為歐拉流體,實(shí)現(xiàn)流體和固體間的雙向耦合交互。但是這種方法在進(jìn)行流體和很薄的固體交互時,由于流體會透過固體約束,發(fā)生漏液現(xiàn)象。Guendelman等[6]改進(jìn)了上述方法,通過使用低維度的三角形表面,處理無限薄的固體模型,消除漏液現(xiàn)象。

5.1.4 cut-cell法

Roble等[35]使用cut-cell方法,將歐拉網(wǎng)格裁切成和不規(guī)則物體邊界一致的形狀,對規(guī)則的有限差分方法進(jìn)行修改,結(jié)合相應(yīng)的流體散度計算方法,使體素化網(wǎng)格可以更好地匹配物體細(xì)節(jié)。Batty等[11]將壓力投影方法改寫為動能最小化問題,通過計算使離散動能最小的離散壓力值,保證生成對稱正半定的線性系統(tǒng),使解算速度提升兩個數(shù)量級。但是都不能處理形變體問題。Robinson-Mosher等[36]使用固體速度作為邊界條件作用于流體,并用流體的壓力作用的固體,實(shí)現(xiàn)形變體的流體耦合處理,之后他們[37]通過分析形變結(jié)構(gòu)的阻尼矩陣,獲得對稱正定線性系統(tǒng),改進(jìn)了之前的數(shù)值不穩(wěn)定問題,但都不能實(shí)現(xiàn)流固間自由滑動。Zarifi等[7](圖3)結(jié)合Chentanez等[38]提出的強(qiáng)耦合流固交互整體系統(tǒng)方案,使用cut-cell離散方案,實(shí)現(xiàn)自由滑動邊界條件,對流體使用有限體積離散化,同時對固體使用有限元方法,有效避免了基于網(wǎng)格方法出現(xiàn)的階梯狀體素化偽影問題。

圖3 使用cut-cell方法的流固交互Fig.3 Fluid-solid interaction by cut-cell method

5.1.5 流函數(shù)法

Ando等[39]使用流函數(shù)(Stream Function)代替壓力投影方法,實(shí)現(xiàn)了一種無條件無散度的窄帶流固交互方法。利用Helmholtz-Hodge分解可以直接計算出最終的速度,但是這種方法并不能很好地處理可壓縮流體。

5.1.6 小結(jié)

歐拉表示法可以容易地實(shí)現(xiàn)流體的不可壓縮性,MAC方法通過將速度和壓力投影到網(wǎng)格上,使得空間離散化過程很方便,而其體素化偽影問題也可以通過cut-cell法得到緩解,但是歐拉流體數(shù)值耗散問題嚴(yán)重。同時,使用歐拉流體和拉格朗日固體交互很難實(shí)現(xiàn)整體系統(tǒng)方案,這使得時間步長和系統(tǒng)穩(wěn)定性都受到一定的限制。

5.2 拉格朗日流體和拉格朗日固體交互

5.2.1 DSC法

Misztal等[40]基于拉格朗日形變單純復(fù)形(Deformable Simplicial Complex,DSC)方法,使用非結(jié)構(gòu)化四面體網(wǎng)格和有限元方法,將NS方程的解特征化為二階優(yōu)化問題。通過偽可壓縮方程求解穩(wěn)定的壓力值,降低數(shù)值耗散現(xiàn)象,但網(wǎng)格重建和倒置問題突出。Clausen等[12]設(shè)計了一個局部的三維網(wǎng)格修復(fù)算法緩解大形變情況下網(wǎng)格倒置問題,可以魯棒地分割或合并四面體網(wǎng)格。同時在四面體元網(wǎng)格上存儲信息,并使用線性基函數(shù)進(jìn)行插值。但是這種局部方法對時間步長要求嚴(yán)格,網(wǎng)格重建問題仍需占據(jù)大量的計算時間。

5.2.2 SPH法

由于基于拉格朗日網(wǎng)格的方法需要頻繁地處理網(wǎng)格的拓?fù)渥兓?因此Solenthaler等[41]引入了局部參考形的定義,利用光滑粒子流體動力學(xué)(Smoothed Particle Hydrodynamics,SPH)和表面重建技術(shù),實(shí)現(xiàn)了融化和凝固效果。Akinci等[42]實(shí)現(xiàn)了SPH流體和任意形狀的剛體的交互。通過邊界粒子采樣固體表面,解決不規(guī)則邊界采樣問題。而且粒子方法可以自然地模擬小尺度細(xì)節(jié)。但是上述兩個方法在處理可壓縮流體邊界的時候會出現(xiàn)穿透問題。Peer等[43]直接從形變梯度上獲取旋轉(zhuǎn)信息,使用一階連續(xù)的內(nèi)核梯度修正算法進(jìn)行旋轉(zhuǎn)估計,結(jié)合不可壓縮SPH方法,簡化了邊界條件計算過程。但是由于加入了更大的顯式摩擦力項(xiàng),導(dǎo)致時間步長被限制。Gissler等[44](圖4)對固體表面進(jìn)行粒子采樣,將流體迭代解算器和固體SPH解算器相連接,通過使用大的時間步長顯著地減少計算時間,并能和不同的SPH解算器相結(jié)合。

圖4 使用SPH方法的流固交互Fig.4 Fluid-solid interaction by SPH

SPH方法在處理鄰近搜索的時候需要消耗大量的計算性能,而且不能很好地構(gòu)造隱式系統(tǒng),并行計算很難實(shí)現(xiàn)。

5.2.3 PBD法

Macklin等[45]使用基于位置的動力學(xué)(Position Based Dynamics,PBD),提出一個易于并行的約束解算器,可以實(shí)時地呈現(xiàn)動畫效果,兼顧簡單性和魯棒性。同時該方法改進(jìn)了摩擦模型,增加了動摩擦和靜摩擦效果。但是PBD沒有使用基于物理的本構(gòu)模型,真實(shí)性不足。

5.2.4 小結(jié)

拉格朗日方法常用的空間離散化拓?fù)浣Y(jié)構(gòu)主要為拉格朗日網(wǎng)格和拉格朗日粒子。拉格朗日網(wǎng)格(如DSC)方法,很好地解決了歐拉網(wǎng)格體素化偽影問題,但是其在處理大形變的時候網(wǎng)格重建和網(wǎng)格倒置問題十分明顯,需要使用小的時間步長和特定的優(yōu)化方法處理,對計算機(jī)性能要求較高;拉格朗日粒子(如SPH、PBD)方法,解決了網(wǎng)格拓?fù)渥兓瘑栴},并能實(shí)現(xiàn)次網(wǎng)格細(xì)節(jié),但是單純的粒子法需要昂貴的鄰近粒子搜索操作,這進(jìn)一步提高了對計算機(jī)性能的要求。不過,拉格朗日流體和固體間的交互便于構(gòu)建整體方案,使整個系統(tǒng)的真實(shí)性更明顯。

5.3 歐拉流體和歐拉固體交互

5.3.1 IBM法

Peskin等[46]在生物流體力學(xué)領(lǐng)域,使用浸入邊界法(Immersed Boundary Method,IBM)實(shí)現(xiàn)了流體和固體的交互。在歐拉視角的笛卡爾網(wǎng)格上對浸入邊界方程進(jìn)行空間離散化,并通過二階龍格庫塔方法實(shí)現(xiàn)時間離散化。隨后,Levin等[47]將這個方法擴(kuò)展到彈性固體,并且可以直接模擬體積化的數(shù)據(jù)模型。解決了對非線性形變的碰撞檢測和處理問題,并消除了黏性偽影問題。由于使用歐拉網(wǎng)格離散化,易于GPU并行計算。但是,該方式對時間步長要求較高,同時不適用于無限薄物體的模擬。Teng等[13]結(jié)合上述歐拉固體表示法和相應(yīng)解算器,實(shí)現(xiàn)了完全基于歐拉方法的不可壓縮流體和形變體的雙向耦合。通過使用隱式積分器,降低對時間步長的要求。但是,這種方法只能處理單向流現(xiàn)象。

5.3.2 EOL法

Sueda等[48]結(jié)合拉格朗日方法和歐拉方法的特點(diǎn),提出了一種包含簡化節(jié)點(diǎn)和歐拉節(jié)點(diǎn)的Eulerian-on-Lagrangian(EOL)方法。該方法可以魯棒、高效、準(zhǔn)確地模擬包含大量約束的系統(tǒng)。隨后,Fan等[49]利用EOL方法不需要顯式離散化的特點(diǎn),將高頻形變的情況用歐拉方法處理,并提出一個基于最小二乘法的接觸感知解算器。在保證系統(tǒng)穩(wěn)定的同時使時間步長最大。但是由于顯式地處理彈性力,導(dǎo)致系統(tǒng)可能出現(xiàn)不穩(wěn)定現(xiàn)象。接著,他們[50]對人體肌肉使用EOL方法處理體積守恒和大形變問題,實(shí)現(xiàn)了肌肉骨骼等關(guān)節(jié)組織的近距離接觸模擬。該方法可以直接通過MRI數(shù)據(jù)建立肌肉模型,但只能實(shí)現(xiàn)較簡單的肌肉運(yùn)動效果。Sachdeva等[51]利用EOL離散方法,消除了準(zhǔn)靜態(tài)方程中不必要的縱向自由度,同時維持橫向運(yùn)動,實(shí)現(xiàn)了對薄/細(xì)東西的模擬。該方法簡化了肌腱和骨頭間的碰撞處理問題,在較大時間步長的前提下實(shí)現(xiàn)了高度剛性的肌腱模擬。但是該方法對參數(shù)過于敏感,且實(shí)現(xiàn)的效果較為簡單。Weidner等[52]將EOL擴(kuò)展到處理三維空間中的二維物體,模擬布料和尖銳集合體的交互,通過在布料中插入EOL頂點(diǎn),實(shí)現(xiàn)尖銳交互點(diǎn)布料的精準(zhǔn)彎曲效果,同時避免了過約束導(dǎo)致的鎖定問題和抖動問題,但是這種方法網(wǎng)格重建效率并不高。

Li等[53]使用歐拉固體的模擬方法,提出了一個現(xiàn)存動畫管線的后處理技術(shù)。通過建模二維空間的任意拓?fù)浣Y(jié)構(gòu)的超彈性膜,可以方便地實(shí)現(xiàn)覆蓋到任意物體上的皮膚或薄膜效果。利用類似于紋理貼圖的方法,可以使用低分辨率網(wǎng)格實(shí)現(xiàn)高分辨率的表面細(xì)節(jié)。

5.3.3 小結(jié)

歐拉固體繼承了歐拉表示法的優(yōu)勢,IBM方法可以直接模擬體積化的數(shù)據(jù),并且便于GPU并行計算;EOL方法雖然是在拉格朗日空間使用歐拉節(jié)點(diǎn),但是其本質(zhì)仍是歐拉方法特性。歐拉方法可以隱式地進(jìn)行空間離散化,不過,在處理彈性固體的時候,歐拉方法對時間步長要求較高。

5.4 混合粒子-網(wǎng)格法

5.4.1 SPH粒子-網(wǎng)格法

Hong等[54]使用SPH結(jié)合歐拉網(wǎng)格,基于SPH渦度約束方法實(shí)現(xiàn)了多相流氣泡和水的交互,并利用新的氣泡模型呈現(xiàn)出次網(wǎng)格氣泡細(xì)節(jié)。但對時間步長要求較高。Patkar等[55]改進(jìn)上述方法,基于Rayleigh-Plesset方程,實(shí)現(xiàn)了氣泡振蕩和合并現(xiàn)象的模擬。利用新的基于渦度的氣泡播種機(jī)制,模擬交互過程產(chǎn)生的氣泡,并實(shí)現(xiàn)了氣泡的可壓縮性和不同尺度氣泡的無縫過渡。Losasso等[56]通過區(qū)分聚集和擴(kuò)散水域,并分別使用水平集方法和SPH粒子模擬。但是,這種方法在計算粒子密度時存在一定誤差,而且可能產(chǎn)生數(shù)值噪聲偽影。Lee等[57]使用 SPH逃逸粒子并結(jié)合BFECC(Back and Forth Error Compensation and Correction)和CIP(Constrained Interpolation Profile)方法確保了二階精度的位置更新。Gao等[58]使用歐拉網(wǎng)格和SPH粒子實(shí)現(xiàn)了高速和低速氣體流。但是這種方法在物理上并不嚴(yán)格的能量守恒。Raveendran等[59]通過對局部SPH粒子進(jìn)行密度修正,并利用粗網(wǎng)格上的泊松解算器計算無散度的速度場。這種方法消除了純網(wǎng)格方法的體素化偽影,并加快了系統(tǒng)收斂速度。Zhu等[60]使用SPH-FLIP方法進(jìn)行流體模擬,可以生成并保留湍流細(xì)節(jié),但是這種方法的效率并不是太高。

5.4.2 Voronoi/Power diagrams法

Sin等[61]基于拉格朗日粒子方法,通過在采樣點(diǎn)上計算Voronoi來處理壓力投影。使用散度定理確保Voronoi單元無散度,并處理了復(fù)雜的邊界條件。這種方法靈活地結(jié)合基于點(diǎn)和基于網(wǎng)格的方法,保證了質(zhì)量守恒。Brochu等[62]提出了一個簡化的Voronoi/Delaunay網(wǎng)格速度插值方案。使用壓力采樣方法,結(jié)合基于網(wǎng)格的曲率估計和流體邊界條件,準(zhǔn)確地獲得表面張力。這種方法雖然消除了顯示表面追溯方法產(chǎn)生的噪聲,但是其他偽影問題依然存在。De Goes等[63]使用拉格朗日冪粒子法,將流體粒子視為體積包,并通過Power diagrams時間序列來描述流體運(yùn)動。通過基于網(wǎng)格的壓力投影,準(zhǔn)確地計算空間中流體體積,有效減少數(shù)值阻尼和粒子漂移問題。

5.4.3 FLIP/PIC法

1995年,Sulsky等[64]將FLIP/PIC擴(kuò)展到固體力學(xué),使用基于物質(zhì)點(diǎn)和計算網(wǎng)格兩種表示法,消除拉格朗日網(wǎng)格的大形變糾纏問題,并利用網(wǎng)格計算空間梯度變化,充分展現(xiàn)了物質(zhì)點(diǎn)法的優(yōu)勢。Zhu等[65]通過改進(jìn)水模擬器參數(shù),實(shí)現(xiàn)對沙子的模擬。利用粒子準(zhǔn)確地進(jìn)行表面追蹤,并使用網(wǎng)格執(zhí)行不可壓縮性和邊界條件。Hong等[66]使用FLIP方法,基于深度和雷諾數(shù),實(shí)現(xiàn)了不可壓縮流體的模擬,并量化了流體不同區(qū)域的形變能力。McAdams等[67]基于FLIP流體解算器,實(shí)現(xiàn)了頭發(fā)的批量交互和體積守恒計算。雖然,基于拉格朗日方法的自碰撞檢測可以準(zhǔn)確地處理頭發(fā)間的精細(xì)交互,但是需要較高的網(wǎng)格分辨率,且存在角動量耗散問題。Han等[68]基于物質(zhì)點(diǎn)法實(shí)現(xiàn)頭發(fā)、線股等共維模型的摩擦接觸和自碰撞處理,利用拉格朗日自由度減少偽影問題,但出現(xiàn)線股壓縮現(xiàn)象時會導(dǎo)致能量不守恒,且不能很好地模擬剛體運(yùn)動。Stomakhin等[69]基于連續(xù)介質(zhì)力學(xué),提出彈塑性本構(gòu)模型,可以自動地處理雪的自碰撞和斷裂效果,并使用半隱式方案可以處理大范圍的剛度問題。但是,雪的參數(shù)需要手動調(diào)整,且碰撞處理為非滑動邊界條件。隨后,他們[70]提出了流體本構(gòu)模型的分割方法,使用物質(zhì)點(diǎn)法模擬熱傳遞、熔化和凝固現(xiàn)象,并提出增強(qiáng)型MPM方法,實(shí)現(xiàn)流體不可壓縮性,但仍不能實(shí)現(xiàn)自由滑動邊界條件。Yan等[71]基于Noyes-Whitney方程,實(shí)現(xiàn)了流體間混溶和不混溶的多流體方案,但僅實(shí)現(xiàn)了弱可壓縮流體。Hu等[17]使用移動最小二乘法和CPIC方法,利用顏色距離場實(shí)現(xiàn)不連續(xù)速度場切割,但無法處理次網(wǎng)格細(xì)節(jié)。Fang等[72]使用接觸面離散積分點(diǎn)方法,改進(jìn)了傳統(tǒng)物質(zhì)點(diǎn)法非滑動邊界條件,實(shí)現(xiàn)非線性彈性固體和不可壓縮流體間的自由滑動。但在處理固體和固體交互時會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。Daviet等[73]利用材料流變學(xué),并結(jié)合體積摩擦和Drucker-Prager屈服準(zhǔn)則,模擬真實(shí)自由流動的顆粒狀物質(zhì)。但只能實(shí)現(xiàn)均質(zhì)材料的模擬。Tampubolon等[74]使用連續(xù)混合理論模擬滑坡或泥石流等濕沙現(xiàn)象,利用動量交互項(xiàng)完成多個類型物質(zhì)間的動量轉(zhuǎn)移并保證動量守恒,但不能實(shí)現(xiàn)毛細(xì)現(xiàn)象。Gao等[15]改進(jìn)了上述方法,使用Drucker-Prager彈塑性流體準(zhǔn)則,實(shí)現(xiàn)多網(wǎng)格流體預(yù)處理解算器。利用新的粒子密度估計方法,實(shí)現(xiàn)多孔材料和水的交互。但由于使用顯式積分方案,在相對速度過大時會出現(xiàn)數(shù)值偽影。Fei等[75](圖5)依靠混合理論更新速度場和飽和度變量,并充分考慮布料中多孔材料的微觀結(jié)構(gòu),實(shí)現(xiàn)了布料對水的吸收和滴落效果。但是這種方法需要較高的網(wǎng)格分辨率和較小的時間步長。

圖5 物質(zhì)點(diǎn)法實(shí)現(xiàn)水透過毛巾Fig.5 Water through towel by material point method

5.4.4 小結(jié)

混合粒子網(wǎng)格方法因其結(jié)合了歐拉網(wǎng)格和拉格朗日粒子的優(yōu)勢,在近些年吸引了越來越多研究者的關(guān)注。SPH方法通過背景網(wǎng)格可以提高鄰近搜索速度,但也引入了額外的計算復(fù)雜度;維諾圖的引入雖然使速度計算更加方便,但是通過估計的方法導(dǎo)致整個系統(tǒng)的物理性和真實(shí)性缺失;物質(zhì)點(diǎn)法通過插值函數(shù)將粒子和網(wǎng)格聯(lián)系起來,使用網(wǎng)格計算壓力,并利用粒子進(jìn)行攜帶物質(zhì)信息,可以實(shí)現(xiàn)多種物質(zhì)的模型,但是插值函數(shù)的引入也引起了額外的數(shù)值耗散問題。

5.5 對比說明

表1總結(jié)了各類典型的流固交互方法的具體實(shí)現(xiàn)和規(guī)模效率信息。其中,詳細(xì)介紹了文獻(xiàn)中的具體實(shí)現(xiàn)方法和表示方法信息,并對空間離散方案和時間積分器進(jìn)行對比。由于各個文獻(xiàn)中使用的硬件信息差別很大,因此通過規(guī)模和效率兩個維度并不能完全論證不同文獻(xiàn)中算法的優(yōu)劣。不過,從表中仍可以看出,基于物理的流固交互方法很難達(dá)到實(shí)時性,尤其在使用粒子法的文獻(xiàn)中,為了追求更加真實(shí)的交互效果,通常都需要幾分鐘才能計算一幀圖像,這也是未來需要優(yōu)化處理的研究方向。

表1 各種流固交互技術(shù)對比Tab.1 Comparison of fluid-solid interaction methods

6 碰撞檢測方法

碰撞檢測是流固交互中必要的處理步驟,根據(jù)處理的先后順序,通常分為兩步:碰撞檢測和碰撞響應(yīng)。碰撞檢測過程一般基于距離進(jìn)行碰撞判斷,當(dāng)兩個物體出現(xiàn)交點(diǎn)或物質(zhì)點(diǎn)發(fā)生穿透,則判定碰撞發(fā)生;而碰撞響應(yīng)通常基于反射定律或特定的邊界條件進(jìn)行處理。Kelager等[76]詳細(xì)介紹了SPH方法中的碰撞檢測邊界條件,由于SPH粒子在邊界處鄰近粒子密度降低,導(dǎo)致邊界處粒子出現(xiàn)過高的黏性。Fujisawa等[77]通過修改SPH鄰近粒子密度計算方法,改善邊界黏性聚集問題。Sifakis等[78]將碰撞響應(yīng)設(shè)置為非線性約束,解決復(fù)雜的、多層布料的碰撞問題。該方案可以產(chǎn)生一個對稱正定的線性系統(tǒng),比傳統(tǒng)貪婪算法速度提高兩個數(shù)量級,但是該方法并沒有優(yōu)化碰撞檢測方法,僅提高了計算速度。McAdams等[79]提出一種可以在不規(guī)則體素化域上處理任意混合的狄利克雷邊界條件和諾依曼邊界條件的純幾何方案。Müller等[80]通過在物體間插入空氣的方法,利用單邊約束條件自動、魯棒地處理糾纏問題和布料模擬中的穿透問題。

7 加速策略

7.1 自適應(yīng)方法

在流固交互領(lǐng)域,自適應(yīng)方法因其局部優(yōu)化的特性,一直都是研究者們的重點(diǎn)研究領(lǐng)域。自適應(yīng)方法將主要的計算工作集中在視覺重點(diǎn)區(qū)域,通過對計算資源的合理分配,達(dá)到加速的目的。通常自適應(yīng)方法的類型和不同的表示方法相關(guān)聯(lián),包括自適應(yīng)網(wǎng)格方法、自適應(yīng)粒子方法、自適應(yīng)插值函數(shù)方法、自適應(yīng)核函數(shù)方法等等。Losasso等[81]提出一種不受限的八叉樹數(shù)據(jù)結(jié)構(gòu)和相應(yīng)的離散方法,并獲得對稱正定的線性系統(tǒng)。Ando等[82]使用空間自適應(yīng)的FLIP方法,減少壓力解算的復(fù)雜度,增加系統(tǒng)魯棒性。但這種方法會產(chǎn)生動量偽影。Ferstl等[83]使用自適應(yīng)八叉樹網(wǎng)格和六面體有限元離散化方法,利用幾何多網(wǎng)格解算器,獲得高效的收斂率。Setaluri等[84]基于硬件加速機(jī)制和虛擬存儲管理系統(tǒng),提出金字塔狀的稀疏均勻網(wǎng)格數(shù)據(jù)結(jié)構(gòu),通過避免間接的內(nèi)存訪問并利用自適應(yīng)多網(wǎng)格預(yù)處理共軛梯度解算器,加速壓力解算過程。Aanjaneya等[85]通過使用冪圖和金字塔式的稀疏均勻網(wǎng)格,增強(qiáng)泊松算子的規(guī)范性。Gao等[86]提出一種自適應(yīng)廣義插值方法,引入GIMP范式和離散化框架,實(shí)現(xiàn)高效的內(nèi)存讀取和并行計算。Chentanez等[87]將空間分為tall cells和cubic cells,簡化數(shù)據(jù)結(jié)構(gòu)和算法,便于系統(tǒng)并行計算。Goswami等[88]基于并行分塊技術(shù),將WCSPH流體分割成不同的計算域,獲得了兩倍的速度提升。

7.2 高速解算器

在解算流固交互過程中,很多系統(tǒng)求解精確解的過程十分困難,甚至不存在實(shí)際的解。為了計算方程的近似根,通常使用牛頓迭代法求解。Mazhar等[89]提出一種加速的投影梯度下降方法,對比傳統(tǒng)的高斯賽德爾和雅可比方法,收斂速度更快,且易于并行。Aanjaneya等[90]基于域分解理論和舒爾補(bǔ)碼方法,使用克利羅夫計算器代替?zhèn)鹘y(tǒng)的泊松解算器,利用multigrid v-cycle預(yù)處理器加速共軛梯度求解。

7.3 硬件并行加速

研究者為了追求更加逼真的流固交互效果,在交互系統(tǒng)中使用的粒子和網(wǎng)格數(shù)量也越來越多,因此許多研究者都考慮使用并行技術(shù)加速整個系統(tǒng)的計算效率。早期的解決方法通常是使用基于 CPU的 SIMD指令[91]或 OpenMP[92]進(jìn)行加速。但受限于CPU的設(shè)計架構(gòu),其在處理大規(guī)模并行加速場景中優(yōu)勢并不明顯;相反近年來,由于GPU的帶寬優(yōu)勢,越來越多的研究者基于GPU架構(gòu)設(shè)計并行加速策略。Wu等[93]使用FLIP方法結(jié)合GPU并行加速策略,實(shí)現(xiàn)了千萬級粒子的流體模擬場景。Hu等[94]提出基于GPU加速的Taichi圖形學(xué)編程語言,可以直接選擇不同的并行加速架構(gòu)。

8 研究趨勢預(yù)測

近年來,隨著計算機(jī)硬件的不斷發(fā)展,研究人員也可以從更多的角度分析并優(yōu)化流固交互中的數(shù)值方法。結(jié)合近些年的研究現(xiàn)狀,未來可能從以下幾個方面推進(jìn)。

1)大規(guī)模場景的實(shí)時交互。

目前,在小規(guī)模場景中,基于物理的流固交互展現(xiàn)出十分真實(shí)的模擬效果,然而在較大規(guī)模場景中,單純的基于物理的流固交互方法并不能很好地得到應(yīng)用。由于神經(jīng)網(wǎng)格等機(jī)器學(xué)習(xí)算法的流行,將基于物理的流固交互技術(shù)和機(jī)器學(xué)習(xí)相結(jié)合將是未來的熱點(diǎn)話題。此外,基于高度場和緩坡方程研究大規(guī)模流固交互場景的實(shí)時模擬技術(shù)也具有重要的意義。

2)基于成熟模型,開發(fā)市場化應(yīng)用。

現(xiàn)有的流固交互技術(shù),已經(jīng)可以實(shí)現(xiàn)非常逼真的交互場景。目前,加拿大的sideFX公司、美國的NVIDIA,以及中國的澤森科工都開始在圖形學(xué)可視化編程軟件領(lǐng)域進(jìn)行探索以降低影視特效、游戲制作、醫(yī)療仿真等非專業(yè)圖形學(xué)人員的使用門檻。未來的研究重點(diǎn)將進(jìn)一步提高軟件的擴(kuò)展性、交互性以及便利性。

3)結(jié)合計算機(jī)硬件,改進(jìn)數(shù)據(jù)結(jié)構(gòu)。

通過結(jié)合計算機(jī)硬件讀寫特點(diǎn),改進(jìn)流固交互數(shù)據(jù)存儲形式,在計算機(jī)帶寬不變的情況下,速度將得到有效的提高。傳統(tǒng)的計算機(jī)軟件,在存儲優(yōu)化這方面主要依靠編程人員的個人技巧,在實(shí)現(xiàn)相同的場景中,不同人開發(fā)的代碼效率也各不相同。特別是在流固交互場景中,百萬級甚至千萬級的粒子數(shù)量,使得單次讀寫微秒級的差別也會放大百萬倍。因此如何使數(shù)據(jù)存儲結(jié)構(gòu)更方便地讀寫和并行計算,將是未來流固交互領(lǐng)域的一個重要研究課題。

9 總結(jié)

本文系統(tǒng)地介紹了近年來流固交互技術(shù)相關(guān)研究工作。為了讓人們在計算機(jī)虛擬世界中感受到真實(shí)的交互場景,流固交互技術(shù)已經(jīng)發(fā)展了幾十年,但仍有很多方面需要進(jìn)一步探索。

基于物理的流固交互技術(shù)因選擇的表示方法不同,處理方案也各不相同,但各種方法都存在優(yōu)勢和劣勢。基于歐拉網(wǎng)格的方法在離散化、壓力解算和并行計算方面更有優(yōu)勢,基于粒子的方法在表現(xiàn)流體細(xì)節(jié)、大形變處理方面更方便,因此基于粒子-網(wǎng)格的混合方法得到了大量研究者的重視。然而,基于物理的流固交互技術(shù)在大規(guī)模場景搭建、實(shí)時交互、市場化應(yīng)用等方面仍面臨著很多問題。

本文從計算機(jī)圖形學(xué)領(lǐng)域中流固交互文獻(xiàn)所使用的表示方法角度,總結(jié)了目前相關(guān)研究中通常使用的具體實(shí)現(xiàn)方法和技術(shù),有助于理解基于物理的流固交互中的數(shù)值方法和關(guān)鍵技術(shù)。