本發(fā)明涉及計算機(jī)圖形學(xué)流體模擬與渲染技術(shù)領(lǐng)域，尤其涉及一種流體模擬方法與系統(tǒng)。

背景技術(shù)：

流體現(xiàn)象對于真實世界的視覺感受而言非常重要。計算機(jī)圖形學(xué)領(lǐng)域為了重現(xiàn)豐富多彩的流體現(xiàn)象，開發(fā)了多種多樣的流體模擬器。

在之前的工作中，大多數(shù)的關(guān)注重點在于如何在較少障礙物的三維空間中模擬流體。而生活中最常見的流體，如空氣，水等具有較為透明的特征，其多樣化的運動現(xiàn)象在固體表面邊界附近時才在視覺上更加明顯。這即是一類表面流流體問題的研究對象：即對于任意的三維模型，研究一薄層流體在模型表面的運動。

圖形學(xué)中存在一些研究流體在二維幾何流型中運動的工作，如z.fan等人2005年發(fā)表“adaptedunstructuredlbmforflowsimulationoncurvedsurfaces”,s.auer等人2012年發(fā)表“real-timefluideffectsonsurfacesusingtheclosestpointmethod”,分別使用lbm方法和最近點方法將二維平面的零厚度流體模擬擴(kuò)展到二維流型中，獲得類似三維流體模擬的投影效果。而對于在幾何模型表面上存在的有限厚度的薄層流體的模擬，h.wang等人在2007年發(fā)表”solvinggeneralshallowwaveequationsonsurfaces”，o.azencot等人在2015年發(fā)表“functionalthinfilmsonsurfaces”,分別使用淺水波方程和基于潤滑理論的優(yōu)化方法獲得了一定的表面流效果。然而，他們的效果只能適用于以可忽略速度運動的高粘滯性流體，并且無法在任意模型表面上做到高度并行，因而速度較慢，無法實時化。

在任意幾何模型表面模擬流體還涉及到較為復(fù)雜的數(shù)值計算技巧。l.shi等人在2004年發(fā)表“inviscidandincompressiblefluidsimulationontrianglemeshes”，s.jeong等人在2013年發(fā)表“combustionwavesonthepointsetsurface”,研究幾何模型表面的移流計算的處理方式。然而前者的計算較為復(fù)雜，后者只涉及標(biāo)量場的計算，而表面流問題需要處理速度矢量場的移流。上述h.wang等人的文章提出了表面流中一種表面張力的計算方式，然而其亦存在近似精度低，假設(shè)較強(qiáng)等問題。n.thurey等人在2006年發(fā)表“animationofopenwaterphenomenawithcoupledshallowwaterandfreesurfacesimulations”,n.chentanez等人在2015年發(fā)表“coupleing3deulerian,heightfieldandparticlemethodsforinteractivesimulationoflargescaleliquidphenomena”,將淺水波模擬器與3d流體模擬器相結(jié)合，然而他們僅能處理向上表面的互相結(jié)合問題。近年來，b.ren等人在2014年發(fā)表“multiple-fluidsphsimulationusingamixturemodel”，研究三維空間中的可混多組分流體運動，但目前圖形學(xué)的相關(guān)研究中尚不涉及表面流中的多組分流體運動計算方法。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的是解決在表面流流體的模擬中，如何高效地處理方程中的復(fù)雜的數(shù)值計算問題，提供一種表面流流體運動的實時模擬方法，以便高效地恢復(fù)真實的表面流運動，并與三維模擬器結(jié)合，恢復(fù)包含多組分流體運動在內(nèi)的多樣化的表面流流體現(xiàn)象，做到任意模型上的實時表面流流體模擬。

本發(fā)明的技術(shù)方案：

一種利用淺水波方程獲得表面流流體運動的實時模擬方法，包括：

第1:利用一種基于特征的粘滯力模型，對幾何模型進(jìn)行預(yù)計算處理；

包括對幾何模型進(jìn)行平滑化并計算平滑化三角網(wǎng)格的外法向，在平滑化三角網(wǎng)格的格點上計算控制面積，粘滯力張量矩陣與局部切平面坐標(biāo)系間仿射變換矩陣；

第2：在模擬開始時載入預(yù)計算的數(shù)據(jù)；

對于可形變模型表面上的模擬，在每個時間步開始時重新計算平滑化三角網(wǎng)格的外法向，仿射變換矩陣與粘滯力張量矩陣；

所述的可形變模型表面上的模擬及重新計算方法進(jìn)一步包括：

在可形變模型表面上的模擬中，在平滑時調(diào)整平滑網(wǎng)格頂點的位置，使得原網(wǎng)格的頂點總在平滑網(wǎng)格頂點沿其法向方向的直線上；之后，通過移動平滑網(wǎng)格頂點的位置或編輯底面高度控制原模型的形變；

第3:使用單組分淺水波方程計算單組分的表面流流體運動，或使用多組分淺水波方程計算多組分的表面流流體運動；

第4:選擇性包括與3d模擬器進(jìn)行交互；

對于包含交互的情況，在第1步中的預(yù)計算時根據(jù)模型生成粒子3d模擬器中的粒子邊界；同時將3d模擬器進(jìn)行一個時間步的計算，并在時間步的最后，計算淺水波方程中三角網(wǎng)格上的模擬器以及粒子3d模擬器間的物理量交互。

其中，

第1步所述預(yù)計算進(jìn)一步包括：

第1.1、將原幾何模型平滑化；

第1.2、對每個平滑化網(wǎng)格頂點，計算鄰域內(nèi)的一個半平整化模型(som)；

所述som的計算方法包括：

第1.2.1、對每個所考慮的平滑模型上的格點v，將其周圍鄰域內(nèi)的點沿該考慮格點的法向投影至切平面；

第1.2.2、v沿其法向?qū)?yīng)一個原模型上的點vp,對于任意原模型上的點vp',同樣可沿同一方向獲得平滑模型上的投影對應(yīng)點v’；

第1.2.3、v，v’,vp,vp'一定能夠構(gòu)成一個平面，與v的切平面交于一直線，vv’與該直線的夾角為θ，則將vp'以vp為中心在該平面內(nèi)旋轉(zhuǎn)θ角度；

第1.2.4、對所有v附近的vp'作上述三步，獲得v鄰域內(nèi)的一個“半平整化模型”(som)；

第1.3、在som上計算每個三角形的形狀運算符s(shapeoperator)；

第1.4、利用第1.3步中的計算結(jié)果獲得每個平滑網(wǎng)格頂點切平面內(nèi)的平均粗糙程度最大的方向和與其垂直的方向，以及兩個方向上的粗糙系數(shù)；

第1.5、利用第1.4步結(jié)果獲得粘滯力張量矩陣，粘滯力張量矩陣與速度矢量的矩陣乘積結(jié)果即為表面流流體運動中所受的模型產(chǎn)生的粘滯力。

第3步所述的計算單組分的表面流流體運動的方法包括：

第3.1、在每個平滑網(wǎng)格點上存儲表面流流體的速度場與深度場；

第3.2、在每個時間步開始時采用半拉格朗日法對速度場與深度場進(jìn)行移流計算；

所述的在每個時間步開始時采用半拉格朗日法對速度場與深度場進(jìn)行移流計算包括：

第3.2.1、在考慮頂點i時，對于所有周圍的頂點j，設(shè)其法向分別為ni、nj，將頂點j的3d坐標(biāo)沿方向ni×nj旋轉(zhuǎn)，使ni與旋轉(zhuǎn)后的nj方向相同；

第3.2.2、此時兩點i、j的切平面平行，計算兩點局部切平面坐標(biāo)系間的放射矩陣mij；

第3.2.3、在半拉格朗日法進(jìn)行速度插值時，即能夠直接使用矩陣mij和頂點j的速度uj插值獲得頂點i的速度ui；

第3.2.4、進(jìn)一步對深度場進(jìn)行插值，方法是直接對原深度場進(jìn)行插值；

第3.3、計算流體所受粘滯力；

第3.4、計算流體所受表面張力壓強(qiáng)；

所述計算流體所受表面張力壓強(qiáng)的方法是：

第3.4.1、判斷網(wǎng)格點處于液體內(nèi)部還是邊緣；

第3.4.2、對于液體內(nèi)部的格點，利用表面能的梯度獲得表面張力，進(jìn)而計算出表面張力壓強(qiáng)；

第3.4.3、對于液體邊緣的格點，利用表面能沿表面切平面方向的梯度，根據(jù)接觸角數(shù)據(jù)，計算出邊緣的表面張力，進(jìn)而計算出表面張力壓強(qiáng)；

第3.5、在網(wǎng)格上計算速度的散度與流體壓強(qiáng)的梯度；

第3.6、根據(jù)淺水波方程獲得速度場與深度場的變化率。

對于固體組分和液體組分的混合組分，在第3.5步和第3.6步之間還應(yīng)當(dāng)包括：

(1)：計算液體與固體間的粘滯作用力；之后：

(2)：計算固體的溶解與析出及其對流體運動造成的影響；

所述計算固體的溶解與析出及其對流體運動造成影響的方法是：

a：根據(jù)用戶指定的溶解度和溶解速率計算固體當(dāng)前時間步的溶解量，或者當(dāng)超過溶解度時的析出量；

b：根據(jù)步驟a中的數(shù)據(jù)計算液體因此損失或獲得的動量，加入到液體速度變化率的計算中。

最后根據(jù)淺水波方程獲得固體與液體兩組分的速度場與深度場的變化率。

第4步所述與3d模擬器進(jìn)行交互的方法是：

第4.1、在模擬開始前根據(jù)模型生成3d模擬器的粒子邊界；

第4.2、同時將粒子3d模擬器進(jìn)行一個時間步的計算；

第4.3、根據(jù)用戶給定的閾值計算淺水波模擬器應(yīng)向粒子3d模擬器傳送的液體質(zhì)量、位置和速度，并生成3d模擬器中的粒子；

第4.4、判斷3d模擬器中運動軌跡與淺水波模擬器模型表面相交的粒子，根據(jù)粒子質(zhì)量計算淺水波方程中的體積源項，影響下一時間步的淺水波方程模擬器計算，同時在3d模擬器中刪除該粒子。

本發(fā)明同時還提供了一種利用淺水波方程獲得表面流流體運動的實時模擬系統(tǒng)，該系統(tǒng)包括：

輸入/輸出模塊：用于模擬過程中數(shù)據(jù)的輸入輸出，包括加載模擬過程中所需的必要物理參數(shù)與用戶指定的參數(shù)，保存流體模擬器所獲得的圖片與數(shù)據(jù)結(jié)果等；

預(yù)計算模塊：用于預(yù)計算數(shù)據(jù)與需要對預(yù)計算數(shù)據(jù)重新計算情況下的重新計算；

載入與更新模塊：用于載入預(yù)計算數(shù)據(jù)，及在每時間步開始時需要對預(yù)計算數(shù)據(jù)重新計算情況下的重新計算；

移流計算模塊：用于計算3d模型表面上的單組分或多組分流體的移流；

表面張力計算模塊：用于計算單組分或多組分的表面流流體的表面張力；

流體速度與深度更新模塊：用于計算單組分或多組分流體的表面流流體速度與深度的變化率；

與3d模擬器交互模塊：用于處理與3d粒子模擬器的交互模擬。

本發(fā)明的優(yōu)點和有益效果：

采用本發(fā)明方法的主要優(yōu)勢在于提供了一個系統(tǒng)性，高效率，能夠處理非可忽略速度流體的表面流流體解決方案。其中，本方法首次提出了基于特征的粘滯模型與可處理表面溶解的表面多組分流體運動模型。本方法對于淺水波方程中移流、表面張力、與3d模擬器結(jié)合等問題的處理方式相比前人的方法也更加高效。本方法可以完全并行在gpu上，獲得實時的模擬速度。

附圖說明

圖1是本發(fā)明算法流程圖；

圖2是復(fù)雜表面上(t＝0.4、2.7、4.5、6.4、9.0s時)沿物體特征的表面流運動效果；

圖3是可形變表面的(t＝0.7、3.5、4.4、7.9、10.0s時)表面流運動效果；

圖4是由藝術(shù)家根據(jù)表面顏色梯度人為指定t獲得的(t＝2.0、4.0s時)藝術(shù)性效果；

圖5是在上述基礎(chǔ)上(t＝2.0、4.0s時)多組分表面流模擬效果示意圖；

圖6是(t＝1.8、6.5、8.3、12.3s時)與3d模擬器結(jié)合的模擬效果；

圖7是(t＝4.6、6.1s時)不同風(fēng)速環(huán)境下的模擬效果；

圖8是獲得頂點v鄰域內(nèi)一個“半平整化模型”(som)示意圖；

圖9是控制面積示意圖，即連接各三角形質(zhì)心的陰影區(qū)域。

圖10是利用淺水波方程獲得表面流流體運動的實時模擬系統(tǒng)框圖。

具體實施方式

一、如圖10所示利用淺水波方程獲得表面流流體運動的實時模擬系統(tǒng)，該系統(tǒng)包括：

輸入/輸出模塊：用于模擬過程中數(shù)據(jù)的輸入輸出，包括加載模擬過程中所需的必要物理參數(shù)與用戶指定的參數(shù)，保存流體模擬器所獲得的圖片與數(shù)據(jù)結(jié)果等；

預(yù)計算模塊：用于預(yù)計算數(shù)據(jù)與需要對預(yù)計算數(shù)據(jù)重新計算情況下的重新計算；

移流計算模塊：用于計算3d模型表面上的單組分或多組分流體的移流；

表面張力計算模塊：用于計算單組分或多組分的表面流流體的表面張力；

流體速度與深度更新模塊：用于計算單組分或多組分流體的表面流流體速度與深度的變化率；

與3d模擬器交互模塊：用于處理與3d粒子模擬器的交互模擬。

二、利用淺水波方程獲得表面流流體運動的實時模擬方法

本發(fā)明實時模擬方法的算法流程(參見圖1)步驟具體如下：

a：利用一種基于特征的粘滯力模型，預(yù)計算所有可以預(yù)計算的數(shù)據(jù)并存儲。

其中，粘滯力張量可以由藝術(shù)家人為指定獲得想要的藝術(shù)效果。對于存在與3d模擬器交互的情況，根據(jù)模型生成3d模擬器中的粒子邊界。

在本發(fā)明方法中，為了成功模擬出真實世界中流體會沿著表面溝壑等細(xì)節(jié)特征流動的現(xiàn)象，提出了一種基于特征的粘滯力模型如下：

首先，將原模型平滑化，本方法的后續(xù)模擬將進(jìn)行在該平滑后的模型上。為了恢復(fù)平滑過程中損失的原模型細(xì)節(jié)信息，在每個平滑模型的格點上將計算一個粘滯力張量t。方法如下：

1.對每個所考慮的平滑模型上的格點v，將其周圍鄰域內(nèi)的點沿該考慮格點的法向投影至切平面；

2.v沿其法向?qū)?yīng)一個原模型上的點vp,對于任意原模型上的點vp',同樣可沿同一方向獲得平滑模型上的投影對應(yīng)點v’；

3.v，v’,vp,vp'一定可以構(gòu)成一個平面，與v的切平面交于一直線，vv’與該直線的夾角為θ，則將vp′以vp為中心在該平面內(nèi)旋轉(zhuǎn)θ角度；

4.對所有v附近的vp′作上述三步，獲得v鄰域內(nèi)的一個“半平整化模型”(som)。如圖8所示。

5.在該som上計算每個三角面片的形狀運算符(shapeoperator)s；

6.對每個v，計算如下等式：

其中，投影相交集iv由所有沿v法向方向投影后，與v在平滑表面上控制面積相交的som上的三角面片組成。a與s為三角面片的面積與形狀運算符。x為v切平面上的方向矢量。該等式的意義在于找到v切平面上“平均粗糙程度最大”的方向與該方向上的粗糙系數(shù)τmax。控制面積如圖9所示，為連接各三角形質(zhì)心的陰影區(qū)域。

7.計算與垂直的方向并計算下式：

8.設(shè)其中e為局部切平面坐標(biāo)系中的坐標(biāo)軸矢量，m為仿射矩陣。則t＝μmdiag{τmax,τ⊥}m^-1

其中μ為粘滯系數(shù)。為計算方便，在實現(xiàn)中用戶可直接指定摩擦系數(shù)γ＝μ/ρ用于代入(2)式后的計算。

9.最后，γb＝t·u

以上基于特征的粘滯力模型將可以獲得自然的流體沿表面特征運動的模擬效果。

在上面所述的模擬方法中，對于非形變模型上的模擬，原始模型的平滑化，各類仿射矩陣計算，各頂點的控制面積與粘滯力張量t均使用一個預(yù)計算步驟計算并存儲，在之后復(fù)數(shù)的模擬過程中可以直接調(diào)用，無需重復(fù)計算。

b：在模擬開始時載入上述預(yù)計算的數(shù)據(jù)，對可形變模型，在每個時間步開始時重新計算必須更新的數(shù)據(jù)。具體而言，我們在預(yù)計算的平滑化時調(diào)整平滑網(wǎng)格頂點的位置，使得原網(wǎng)格的頂點總在平滑網(wǎng)格頂點沿其法向方向的直線上。之后，我們通過移動平滑網(wǎng)格頂點的位置或編輯底面高度控制原模型的形變。這樣，僅需重新計算平滑網(wǎng)格的法向，仿射矩陣與粘滯力張量。我們的方法采用后面這種較為高效的手段模擬可形變模型上的表面流運動。

c：采用淺水波公式計算單組分或多組分的表面流流體運動，具體項(移流，表面張力，粘滯力等)的技術(shù)實現(xiàn)方法具體說明如下。

在本發(fā)明中，采用淺水波方程作為基本模擬方法，其基本方程為：

在上述兩式中，第(1)式為表面流流體沿表面法向方向的深度的變化方程(質(zhì)量方程)，第(2)式為流體的速度變化方程(運動方程)。其中d為流體沿表面法向方向的深度，u為速度，p為壓強(qiáng)，由靜水壓ρgh與表面張力壓強(qiáng)ps相加而得，ρ為密度，γ＝γw-γb為風(fēng)帶來的粘滯力γw與底面帶來的粘滯力γb之差，aext為外力等導(dǎo)致的加速度，q為體積源。

為風(fēng)所帶來的粘滯力計算公式，γw為用戶設(shè)定的系數(shù)，uw為風(fēng)速。

對于(1)(2)兩式的移流項，本方法中采用半拉格朗日法處理。對于矢量場的移流快速插值，本方法提出使用如下方式：

(1)在考慮頂點i時，對于所有周圍的頂點j，設(shè)其法向分別為ni，nj，將頂點j的3d坐標(biāo)沿方向ni×nj旋轉(zhuǎn)，使ni與旋轉(zhuǎn)后的nj方向相同。

(2)此時兩點i，j的切平面平行，可以計算兩點局部切平面坐標(biāo)系間的放射矩陣mij。

(3)在半拉格朗日法進(jìn)行速度插值時，即可直接使用矩陣mij和頂點j的速度uj插值獲得頂點i的速度ui。

對于頂點v處表面張力的計算，我們基于表面能觀點使用如下等式：

求和對所有包括所考慮的頂點v的流體表面的三角面片進(jìn)行，式中的梯度表示僅對v的坐標(biāo)進(jìn)行的偏導(dǎo)算子。σ為表面張力系數(shù)。

對液體與固體之間的接觸角，設(shè)其為θc，計算

σlg＝σ,

則

其中三個求和分別對所有包含v的液體表面，液體-固體交面，固體-氣體交面進(jìn)行。其中的梯度算子代表沿v的切平面的二維偏導(dǎo)。

最后，對液體內(nèi)部的頂點，使用作為表面張力壓強(qiáng)，對液體邊界上的點，使用作為表面張力壓強(qiáng)。其中面積a的垂直角標(biāo)分別代表垂直于液體表面的投影面積與垂直與v的切平面上液體邊界線法向的投影面積。該表面張力壓強(qiáng)的正負(fù)號可根據(jù)(5)(6)兩式左邊的力矢量與模型外法向是否同向確定：同向為負(fù)，異向為正。

本方法還提供了表面多組分流體運動的模擬方式，即液體中溶解固體的模擬，公式如下：

角標(biāo)l,s分別代表液體、固體的物理量。

其中γl＝γbl+γdrag，由以下兩式計算：

γbl＝λγlt·ul-κ(λ-1)γst

γdrag＝κ(λ-1)cd(ul-us)|ul-us|

其中

而

定義當(dāng)前溶解量qs＝ep(dl+ds)|ul|，則液體因固體溶解的動量損失

其中ep為用戶設(shè)置的溶解速率。

在模擬中，本方法設(shè)置一個溶解率上限，在模擬時檢查ds與dl的比值不超過該上限。對于ds超過的部分，將其從中減去并黏附回模型表面上使得底面高度b相應(yīng)增加同樣的量。

d：對于包含3d模擬器交互的情況，同時將3d模擬器進(jìn)行一個時間步，在時間步的最后，處理兩個模擬器間的物理量交互。

本方法中，利用體積源q處理與3d模擬器的結(jié)合問題。當(dāng)某點i的液體深度d大于用戶指定的∈1時，我們將d設(shè)為用戶指定的∈2，并將總計體積為(d-∈2)ai的液體粒子加入sph三維模擬器，其中ai為i的控制面積。該sph粒子的初始位置為其中vi是i的世界位置坐標(biāo)，b是淺水波模型中的底面高度，n為該處的外法向。該粒子的初速為其中分子上為當(dāng)前時間步計算出的液體深度變化。于此同時，我們還搜索i點直接相連的頂點，若其深度大于用戶指定的(稍小的)值η1，用與上面同樣的后續(xù)步驟將值設(shè)置為用戶指定的η2。

反過來，當(dāng)一個sph粒子運動至模型表面上時，我們認(rèn)為其被表面流體吸收。此時在i點(1)式中，q在一個時間步內(nèi)臨時增加同時刪除該sph粒子。

圖2至圖7分別給出了采用本發(fā)明方法進(jìn)行實時模擬的效果，其中各參數(shù)均被設(shè)置如下：

γ＝1.3,σ＝1.0,κ＝2.0,cd＝10.0,ep＝0.005

1＝3.0,2＝1.6,η1＝2.4,η2＝2.1。

圖2是復(fù)雜表面上(t＝0.4、2.7、4.5、6.4、9.0s時)沿物體特征的表面流運動效果；所使用的模型頂點數(shù)約44000，模擬速度為每幀7.44ms。

圖3是可形變表面(t＝0.7、3.5、4.4、7.9、10.0s時)的表面流運動效果；模型定點數(shù)23200，模擬速度41.36ms/幀，其中用于形變表面信息重新計算的時間為28.92ms/幀。

圖4是由藝術(shù)家根據(jù)表面顏色梯度人為指定t獲得的(t＝2.0、4.0s時)藝術(shù)性效果；其中底面顏色的梯度方向被指定為最大粗糙程度的方向，相應(yīng)粗糙系數(shù)數(shù)值為顏色梯度值，且該粗糙系數(shù)被指定為5倍于垂直方向的粗糙系數(shù)；摩擦系數(shù)被指定為80網(wǎng)格頂點數(shù)52200，速度9.94ms/幀。

圖5是在上述基礎(chǔ)上(t＝2.0、4.0s時)多組分表面流模擬效果示意圖；其中，固液密度比為2：1，速度15.24ms/幀。

圖6是(t＝1.8、6.5、8.3、12.3s時)與3d模擬器結(jié)合的模擬效果；模型包含23200個頂點，粒子模擬器中共有129000個粒子。模擬時間為使用淺水波的三角網(wǎng)格上模擬器20.52ms/幀，粒子模擬器26.73ms/幀，合計47.25ms/幀。

圖7是(t＝4.6、6.1s時)不同風(fēng)速環(huán)境下的模擬效果；圖中由上到下為：無風(fēng)；風(fēng)迎面吹來；風(fēng)由左手方向吹來的結(jié)果；使用的風(fēng)速場由另外的現(xiàn)有三維網(wǎng)格氣體模擬器事先生成。