專利名稱：：使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體詳細(xì)運(yùn)動(dòng)的方法
技術(shù)領(lǐng)域：
：本發(fā)明涉及一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)的運(yùn)動(dòng)的方法。
背景技術(shù)：
：1介紹當(dāng)水猛烈地與固體、空氣、或水本身相互作用時(shí)，其呈現(xiàn)各種結(jié)構(gòu)，包括液滴/泡、薄水膜、以及漩渦，如圖7中所示。當(dāng)流體經(jīng)歷以高雷諾數(shù)為特征的運(yùn)動(dòng)時(shí)能夠出現(xiàn)這樣的特征，其中，雷諾數(shù)表示慣性與流體的粘度之比的相對(duì)大小。高速移動(dòng)的水是典型的高雷諾數(shù)流體。本發(fā)明探索這樣的現(xiàn)象的基于物理的模擬。假設(shè)納維爾-斯托克斯方程正確地模擬流體的物理運(yùn)動(dòng)，看似真實(shí)的模擬應(yīng)能夠再現(xiàn)水的高雷諾數(shù)行為。雖然本工作主要討論水，但也可應(yīng)用于任何流體。然而，迄今為止，尚未令人滿意地再現(xiàn)這些現(xiàn)象。本發(fā)明認(rèn)識(shí)到導(dǎo)致這種失敗的因素，并提出一種允許高雷諾數(shù)液體運(yùn)動(dòng)的真實(shí)的模擬的方法。看似不真實(shí)的高雷諾數(shù)流體的模擬與數(shù)值耗散有關(guān)。具體地說，納維爾-斯托克斯方程的離散化模擬不可避免地需要估計(jì)非網(wǎng)格點(diǎn)處的物理量。在大多數(shù)方法中，通過網(wǎng)格點(diǎn)處的物理量的值的插值來計(jì)算這樣的值。然而，由此近似引入的誤差起到類似于非物理粘度或數(shù)值耗散的作用。能夠通過使用更細(xì)的網(wǎng)格來降低這種耗散；然而，增加網(wǎng)格分辨率會(huì)將計(jì)算時(shí)間和存儲(chǔ)器需求增加至不切實(shí)際的水平。在過去的幾年中，已提出多種優(yōu)秀的技術(shù)來解決這個(gè)問題。將質(zhì)點(diǎn)引入歐拉方案可以幫助捕捉流體的慣性運(yùn)動(dòng)，并增加模擬的有效分辨率。Enright等人[2002b]介紹了質(zhì)點(diǎn)水平集方法，其通過在界面周圍散布質(zhì)點(diǎn)來增加表面追蹤的準(zhǔn)確性。Losasso等人[2004]提出了一種基于八叉樹的多水平流體解算器，其允許在諸如水表面的更感興趣的區(qū)域中的更精細(xì)分辨率的模擬。不幸的是，以上技術(shù)不能產(chǎn)生具有充分細(xì)節(jié)和真實(shí)性的高雷諾數(shù)液體。被模擬的流體表現(xiàn)出比在真實(shí)物理中更大的粘度；流體常?？雌饋砗?粘，而且通常在復(fù)雜流動(dòng)中不表現(xiàn)出小尺度特征的運(yùn)動(dòng)。申請(qǐng)人發(fā)現(xiàn)以上耗散抑制技術(shù)產(chǎn)生這種人為瑕疵的原因是它們不能有效地抑制速度耗散，即使它們顯著地降低了質(zhì)量耗散。在本發(fā)明中，申請(qǐng)人介紹了一種稱為導(dǎo)數(shù)質(zhì)點(diǎn)(derivativeparticle)的新概念，并開發(fā)了一種基于該概念的流體模擬器。申請(qǐng)人利用用于平流部分的模擬的拉爾朗日方案的無耗散性質(zhì)；對(duì)于需要模擬細(xì)節(jié)的流體區(qū)域，申請(qǐng)人使用質(zhì)點(diǎn)來求解平流步驟?；诰W(wǎng)格與基于質(zhì)點(diǎn)的模擬器之間的切換可能引入數(shù)值耗散。本發(fā)明開發(fā)了一種允許詳細(xì)流體運(yùn)動(dòng)的再現(xiàn)的新的轉(zhuǎn)換程序。此工作的一個(gè)創(chuàng)新方面在于，除了存儲(chǔ)物理量(速度和水平集值)之外，所述導(dǎo)數(shù)質(zhì)點(diǎn)還存儲(chǔ)那些量的空間導(dǎo)數(shù)，這使得能夠?qū)崿F(xiàn)非網(wǎng)格/非質(zhì)點(diǎn)位置處的物理量的更精確評(píng)估。本工作的質(zhì)點(diǎn)的使用與質(zhì)點(diǎn)水平集方法[Enright等人2002b]的不同之處在于導(dǎo)數(shù)質(zhì)點(diǎn)攜帶流體速度以及水平集值。實(shí)驗(yàn)顯示提出的模擬器準(zhǔn)確地追蹤界面。更有趣的是，提出的方法被證實(shí)顯著降低非物理阻尼，允許在真實(shí)的高雷諾數(shù)流體中發(fā)生的小尺度特征的宏觀運(yùn)動(dòng)的再現(xiàn)。本說明書的其余部分如下組織第2部分回顧先前的工作；第3部分給出所述模擬器的概述；第4部分介紹基于八叉樹的約束插值剖面(CIP)解算器；第5部分介紹導(dǎo)數(shù)質(zhì)點(diǎn)模型；第6部分報(bào)告我們的實(shí)驗(yàn)結(jié)果；第7部分總結(jié)本說明書。2先前的工作由于Foster和Metaxas[1996；1997]首先介紹了基于完全3D納維爾-斯托克斯模擬的流體動(dòng)畫技術(shù)，所以該方法已在計(jì)算機(jī)圖形學(xué)界盛行。JosStam[1999]介紹了一種稱為穩(wěn)定流體(StableFluids)的穩(wěn)定流體解算器。該解算器的平流步驟使用半拉格朗日方法[Staniforth和C^ot`e1991]來實(shí)施，其即使在使用大時(shí)步時(shí)也保持穩(wěn)定。從那時(shí)起，已有計(jì)算機(jī)圖形學(xué)中基于半拉格朗日方法的快速流體動(dòng)畫技術(shù)的積極開發(fā)。Rasmussen等人[2003]提出了一種使用2D半拉格朗日解算器來產(chǎn)生氣體的3D大尺度動(dòng)畫的技術(shù)，且Feldman等人[2003]提出了一種將基于質(zhì)點(diǎn)的燃燒模型并入穩(wěn)定流體解算器的爆炸模型。Treuille等人[2003；2004]介紹了一種生成滿足指定的關(guān)鍵幀約束的流體流動(dòng)的優(yōu)化技術(shù)。另外，F(xiàn)eldman等人提出了結(jié)合交錯(cuò)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格[Feldman等人2005a]的混合網(wǎng)格的使用，而且還提出了可變形網(wǎng)格的使用[Feldman等人2005b]，對(duì)于該方法他們必須修改半拉格朗日方法。為了模擬流體，除穩(wěn)定的納維爾-斯托克斯解算器之外，還需要用于追蹤液體表面的模型。為此，F(xiàn)oster和Fedkiw[2001]提出了一種將無質(zhì)量質(zhì)點(diǎn)引入到水平集場(chǎng)的混合表面模型。這種模型推動(dòng)了質(zhì)點(diǎn)水平集方法[Enright等人2002b]的開發(fā)，所述質(zhì)點(diǎn)水平集方法能夠以顯著的準(zhǔn)確率捕捉流體表面的動(dòng)態(tài)運(yùn)動(dòng)。Enright等人[2005]證明該質(zhì)點(diǎn)水平集方法允許半拉格朗日方案中的大時(shí)步的使用。該方法已被用于建模流體與剛性體之間的相互作用[Carlson等人2004]以及流體與可變性薄殼物體之間的相互作用[Guendelman等人2005]，而且被用于模擬粘彈性流體[Goktekin等人2004]、沙[Zhu和Bridson2005]、多相流體[Hong和Kim2005]、以及水滴[Wang等人2005]。作為基于網(wǎng)格的方法的替代，還已研究了純粹基于質(zhì)點(diǎn)的方法。Stam和Fiume[1995]使用光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)(SPH)來建模氣態(tài)現(xiàn)象。在SPH中，通過顆粒的集合來表現(xiàn)流體，且模擬器通過計(jì)算納維爾-斯托克斯方程的每一項(xiàng)來計(jì)算它們的運(yùn)動(dòng)。M¨uller等人[2003]使用SPH模型來模擬流體，并且在[M¨uller等人2005]中，他們使用該技術(shù)來模擬多相流體相互作用。Premoˇze等人[2003]將移動(dòng)質(zhì)點(diǎn)半隱式(MPS)法引入到圖形學(xué)界，這是一種在模擬流體的不可壓縮運(yùn)動(dòng)的過程中顯示出比SPH更好的性能的技術(shù)。然而，純粹基于質(zhì)點(diǎn)的方法的基本缺陷是它們?nèi)鄙俦砻孀粉櫮芰?；如果使用不適當(dāng)?shù)馁|(zhì)點(diǎn)數(shù)目，則它們可能產(chǎn)生粒狀或過度光滑表面。削弱模擬結(jié)果的視覺真實(shí)性(以及物理準(zhǔn)確性)的主要因素是速度場(chǎng)的數(shù)值耗散。為了降低模擬氣態(tài)現(xiàn)象時(shí)的耗散，F(xiàn)edkiw等人[2001]使用三次插值。它們還包括稱為渦度約束(vorticityconfinement)的附加步驟，其將速度場(chǎng)的漩度(curl)放大，產(chǎn)生煙氣運(yùn)動(dòng)中的真實(shí)的渦漩形部分。通過采用渦流(vortex)質(zhì)點(diǎn)法來執(zhí)行渦度耗散的更基于物理的預(yù)防[Selle等人2005；Park和Kim2005]。這種方法用漩度型納維爾-斯托克斯方程來工作并利用質(zhì)點(diǎn)來求解(solve)平流項(xiàng)，這導(dǎo)致渦度的有效保持。然而，在液體的情況下，粘度表現(xiàn)出更突出的作用；特別是，渦漩形運(yùn)動(dòng)存在時(shí)間短而且較少觀察到。因此，建模高雷諾數(shù)液體的行為需要在更一般的情況下工作的速度耗散抑制技術(shù)。為了降低液體模擬中的耗散，Song等人[2005]提出了一種基于CIP平流方法的技術(shù)。他們的技術(shù)用三階精度來求解速度平流。然而，通過這種方法，不能避免來自基于網(wǎng)格的平流的數(shù)值粘度。Zhu和Bridson[2005]提出了一種基于流體隱式質(zhì)點(diǎn)(FLuid-Implicit-Particle)(FLIP)方法的流體模擬技術(shù)。該FLIP方法[Brackbill和Ruppel1986]用質(zhì)點(diǎn)來求解平流步驟，而不包括網(wǎng)格的所有其它步驟。雖然這種方法中的質(zhì)點(diǎn)平流沒有數(shù)值耗散，但當(dāng)在網(wǎng)格與質(zhì)點(diǎn)之間來回傳遞速度值時(shí)，引入了插值誤差。因此，長期以來存在對(duì)使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法的需要。本發(fā)明意在解決這些問題并滿足期盼已久的需要。
發(fā)明內(nèi)容本發(fā)明設(shè)法解決現(xiàn)有技術(shù)的缺點(diǎn)。本發(fā)明的目的是提供一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法。本發(fā)明的另一目的是提供一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法，其中，用拉格朗日方案來實(shí)施平流部分。本發(fā)明的又一目的是提供一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法，其提供在平流與非平流部分的切換處轉(zhuǎn)換物理量的程序，抑制不必要的數(shù)值耗散的介入。3概述在多相流體解算器的開發(fā)中，申請(qǐng)人假設(shè)空氣和水兩者都是不可壓縮的。不可壓縮流體的納維爾-斯托克斯方程可以寫成其中u是速度，p是壓力，f是外力，μ是粘度系數(shù)，ρ是流體密度。為了準(zhǔn)確地建?？缭浇橘|(zhì)之間的界面的密度和粘度的不連續(xù)性，申請(qǐng)人采用具有可變密度的虛擬流體方法[Liu等人2000；Kang等人2000；Hong和Kim2005]。在我們的解算器中，表面追蹤是基于水平集方法。根據(jù)下式而更新水平集場(chǎng)φ符號(hào)距離條件對(duì)于納維爾-斯托克斯方程的時(shí)間積分，申請(qǐng)人采用分步法[Stam1999]，其遞增地計(jì)算(accountsfor)等式(1)中的項(xiàng)的效果和等式(2)中的質(zhì)量守恒條件。申請(qǐng)人開發(fā)的技術(shù)在這里關(guān)注于平流部分的增加。因此，如圖8所示，申請(qǐng)人將分步(fractionalsteps)分組成兩部分，非平流部分和平流部分。等式(3)的時(shí)間積分僅涉及平流部分。平流部分由網(wǎng)格平流部分和質(zhì)點(diǎn)平流部分組成。所述網(wǎng)格平流部分平流輸送根據(jù)當(dāng)前速度場(chǎng)從非平流部分計(jì)算的速度和水平集場(chǎng)。我們的方法與純粹歐拉模擬的不同之處在于我們的模擬器包括質(zhì)點(diǎn)平流部分，其用于模擬需要產(chǎn)生細(xì)節(jié)的區(qū)域。為此，申請(qǐng)人在空氣-水界面中引入了大量質(zhì)點(diǎn)。在網(wǎng)格平流部分中，平流輸送速度和平流輸送水平集是歐拉平流步驟。這些步驟的精確處理形成高雷諾數(shù)行為的成功模擬的基礎(chǔ)。對(duì)于歐拉平流，申請(qǐng)人開發(fā)了一種基于八叉樹的CIP解算器(第4部分)，其將速度和質(zhì)量耗散降低到顯著的水平。當(dāng)流體區(qū)域的平流需要使用質(zhì)點(diǎn)平流部分來模擬時(shí)，當(dāng)前在網(wǎng)格上限定的物理量(例如水平集和速度)被輸送到質(zhì)點(diǎn)。在該輸送之后，質(zhì)點(diǎn)能夠被直接地平流輸送。質(zhì)點(diǎn)平流的結(jié)果然后被傳回到網(wǎng)格。在程序中的這一點(diǎn)上，不管模擬是經(jīng)由網(wǎng)格平流部分還是質(zhì)點(diǎn)平流部分進(jìn)行的，最終速度和水平集值都存儲(chǔ)在網(wǎng)格點(diǎn)上。申請(qǐng)人在第5部分中介紹了網(wǎng)格至質(zhì)點(diǎn)和質(zhì)點(diǎn)至網(wǎng)格速度/水平集轉(zhuǎn)換程序。基于質(zhì)點(diǎn)的平流的使用和轉(zhuǎn)換程序的開發(fā)是再現(xiàn)流體運(yùn)動(dòng)的細(xì)節(jié)所必不可少的。一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法包括以下步驟a)基于八叉樹數(shù)據(jù)結(jié)構(gòu)用自適應(yīng)網(wǎng)格來建模流體；b)對(duì)于網(wǎng)格點(diǎn)處的流體速度求解不可壓縮流體的納維爾-斯托克斯方程；c)在納維爾-斯托克斯方程的時(shí)間積分中使用分步法；d)將分步分組成非平流部分和平流部分；e)根據(jù)模擬的詳細(xì)程度來選擇網(wǎng)格平流部分或質(zhì)點(diǎn)平流部分；f)對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP方法；g)對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型。當(dāng)模擬的詳細(xì)程度高時(shí)，使用對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟。對(duì)于高度感興趣的區(qū)域的細(xì)節(jié)自適應(yīng)網(wǎng)格增加了網(wǎng)格的數(shù)目。八叉樹數(shù)據(jù)結(jié)構(gòu)包括樹形數(shù)據(jù)結(jié)構(gòu)，每個(gè)內(nèi)部結(jié)點(diǎn)具有多達(dá)八個(gè)孩子。不可壓縮流體的納維爾-斯托克斯方程可以寫成其中u是速度，p是壓力，f是外力，μ是粘度系數(shù)，ρ是流體密度。該方法還可以包括使用水平集方法來追蹤流體的表面的步驟。所述水平集方法包括水平集場(chǎng)φ。根據(jù)下式來更新所述水平集場(chǎng)φ符號(hào)距離的時(shí)間積分僅涉及平流部分。網(wǎng)格平流部分平流輸送根據(jù)當(dāng)前速度場(chǎng)從非平流部分計(jì)算的速度和水平集場(chǎng)。網(wǎng)格平流部分使用歐拉平流步驟來平流輸送速度和水平集。通過基于八叉樹的CIP(約束插值剖面)法來執(zhí)行歐拉平流。對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP方法的步驟包括用基于八叉樹的CIP方法來平流輸送速度和水平集及其導(dǎo)數(shù)的步驟。直接從方程獲得平流輸送導(dǎo)數(shù)的方程，其中，φ是被平流輸送的物理量，且通過對(duì)x微分該方程來獲得空間變量x的導(dǎo)數(shù)的平流輸送方程，即求解x的導(dǎo)數(shù)的平流輸送方程的步驟包括采用使用分步法的步驟，其使用分步法的步驟包括以下步驟a)使用有限差分法來求解非平流部分b)平流輸送根據(jù)的結(jié)果；以及c)在由網(wǎng)格限定的單元(cell)的中心處進(jìn)行壓力值取樣并在結(jié)點(diǎn)處進(jìn)行包括速度、水平集及其導(dǎo)數(shù)的其它值的取樣。該方法可以進(jìn)一步包括以下步驟a)將基于八叉樹的CIP方法與自適應(yīng)網(wǎng)格相結(jié)合；以及b)通過用基于八叉樹的CIP方法模擬平流來使由八叉樹數(shù)據(jù)結(jié)構(gòu)的網(wǎng)格分辨率中的區(qū)域性變化引起的八叉樹人為瑕疵最小化。導(dǎo)數(shù)質(zhì)點(diǎn)模型使用基于質(zhì)點(diǎn)的拉格朗日方案。用和來計(jì)算網(wǎng)格點(diǎn)處的水平集及其導(dǎo)數(shù)，其中p是由導(dǎo)數(shù)質(zhì)點(diǎn)攜帶的梯度。該方法可以進(jìn)一步包括以下步驟a)在對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟之前將網(wǎng)格上限定的速度和水平集轉(zhuǎn)換成質(zhì)點(diǎn)上限定的速度和水平集(網(wǎng)格至質(zhì)點(diǎn))；以及b)在返回到對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP模型的步驟之前將質(zhì)點(diǎn)上限定的速度和水平集轉(zhuǎn)換成網(wǎng)格上限定的速度和水平集(質(zhì)點(diǎn)至網(wǎng)格)。轉(zhuǎn)換速度和水平集網(wǎng)格至質(zhì)點(diǎn)的步驟包括下述步驟，即對(duì)于速度uP的轉(zhuǎn)換使用具有用單調(diào)CIP插值來代替線性插值的FLIP方法的步驟，且其中是就在非平流部分開始之前的質(zhì)點(diǎn)速度，且是由于非平流部分而引起的Gi處的速度變化。對(duì)于速度uP、網(wǎng)格點(diǎn)G、選自G的每個(gè)象限的最近質(zhì)點(diǎn)P1、P2、P3和P4、Pi的速度以及和的x和y分量的導(dǎo)數(shù)，轉(zhuǎn)換速度和水平集質(zhì)點(diǎn)至網(wǎng)格的步驟包括以下步驟a)將的導(dǎo)數(shù)坐標(biāo)變換成從而表示沿著方向的分量且表示垂直于的分量，其中，類似地獲得的坐標(biāo)變換的導(dǎo)數(shù)b)沿著方向?qū)Σ襟E(a)中獲得的結(jié)果進(jìn)行一維單調(diào)CIP插值以便計(jì)算位置A處的uA及其方向?qū)?shù)M；c)通過應(yīng)用上述坐標(biāo)變換的反轉(zhuǎn)(inverse)來從獲得(uAx，uAy)；d)對(duì)質(zhì)點(diǎn)P3和P4執(zhí)行步驟(a)～(c)以計(jì)算B的uB和(uBx，uBy)；以及e)對(duì)步驟(c)和(d)中獲得的結(jié)果執(zhí)行y方向單調(diào)CIP插值，其給出G處的速度和導(dǎo)數(shù)。對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟包括以下步驟a)調(diào)整質(zhì)點(diǎn)的水平集值；以及b)執(zhí)行質(zhì)點(diǎn)重新播種。本發(fā)明的優(yōu)點(diǎn)是(1)所述使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法降低質(zhì)量和速度兩者的耗散，導(dǎo)致動(dòng)態(tài)流體的真實(shí)再現(xiàn)；(2)所述使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法用拉格朗日方案來實(shí)施平流部分；以及(3)所述使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法提供質(zhì)點(diǎn)至網(wǎng)格和網(wǎng)格至質(zhì)點(diǎn)程序。雖然簡(jiǎn)要地概述了本發(fā)明，但通過以下附圖、詳細(xì)說明及權(quán)利要求書，能夠獲得對(duì)本發(fā)明的更全面的理解。參照附圖將更透徹地理解本發(fā)明的這些及其它特征、方面和優(yōu)點(diǎn)。圖1是示出了根據(jù)本發(fā)明的使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法的流程圖；圖2是示出了根據(jù)本發(fā)明的方法的另一流程圖；圖3是示出了根據(jù)本發(fā)明的方法的又一流程圖；圖4是示出了網(wǎng)格至質(zhì)點(diǎn)和質(zhì)點(diǎn)至網(wǎng)格轉(zhuǎn)換的一部分流程圖；圖5是示出了網(wǎng)格至質(zhì)點(diǎn)和質(zhì)點(diǎn)至網(wǎng)格轉(zhuǎn)換的另一部分流程圖；圖6是示出了使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟的細(xì)節(jié)的一部分流程圖；圖7示出了由球的沖擊而產(chǎn)生的水的模擬結(jié)果；圖8示出了根據(jù)本發(fā)明的方法的體系結(jié)構(gòu)；圖9示出了從質(zhì)點(diǎn)速度到網(wǎng)格速度的計(jì)算；圖10示出了從2D潰壩(breaking-dam)模擬截取的截圖上部和下部序列分別示出了用質(zhì)點(diǎn)水平集方法的傳統(tǒng)線性模型和導(dǎo)數(shù)質(zhì)點(diǎn)模型產(chǎn)生的結(jié)果；圖11示出了從水滴模擬截取的截圖左圖像和右圖像分別示出了用質(zhì)點(diǎn)水平集方法的傳統(tǒng)線性模型和導(dǎo)數(shù)質(zhì)點(diǎn)模型產(chǎn)生的結(jié)果；圖12示出了從旋轉(zhuǎn)的水箱的模擬截取的截圖；以及圖13示出了從洪水模擬截取的截圖底端行的三個(gè)圖像示出了洪水的推進(jìn)(頂視圖)；頂部圖像示出了側(cè)視圖。具體實(shí)施例方式4開發(fā)基于八叉樹CIP解算器本部分描述通過將CIP方法與八叉樹數(shù)據(jù)結(jié)構(gòu)相結(jié)合而進(jìn)行的模擬器的網(wǎng)格平流部分的開發(fā)。4.1CIP方法介紹在用分步法來求解納維爾-斯托克斯方程(第3部分)的過程中，平流項(xiàng)和由于其雙曲線性質(zhì)而需要特別注意。半拉格朗日方法提供了一種穩(wěn)定方案，在該方案內(nèi)處理以上雙曲線方程[Stam1999]。然而，不幸的是，這種方法遭遇嚴(yán)重的數(shù)值耗散，該數(shù)值耗散源于用來根據(jù)相鄰網(wǎng)格點(diǎn)處的物理量而確定回推(backtracked)的非網(wǎng)格點(diǎn)處的物理量的線性插值。同時(shí)CIP方法是最初由Yabe和Aoki[1991a；1991b]提出并隨后由Yabe等人[2001]改進(jìn)的三階技術(shù)。這種方法的關(guān)鍵思想是不僅平流輸送其物理量，而且還平流輸送其導(dǎo)數(shù)。這里，問題將是如何平流輸送導(dǎo)數(shù)。Yabe和Aoki[1991a]的感興趣的觀察結(jié)果是能夠直接從初始雙曲線方程獲得用于平流輸送導(dǎo)數(shù)的方程其中，φ是正在被平流輸送的物理量。對(duì)于空間變量x微分方程(5)，申請(qǐng)人得到其可以用來預(yù)測(cè)x隨時(shí)間的演變。在求解方程(6)的過程中，申請(qǐng)人再次應(yīng)用分步法首先，申請(qǐng)人使用有限差分來求解非平流部分然后，申請(qǐng)人平流輸送根據(jù)下式的結(jié)果對(duì)CIP平流的詳細(xì)介紹可以在[Song等人2005]中找到。以上方法的二維和三維情況的擴(kuò)展存在于[Yabe和Aoki1991b]。Song等人[2005]發(fā)現(xiàn)[1991b]中給出的推廣可能引起不穩(wěn)定，并提出解決這個(gè)問題的修改。在本工作中，為了求解平流部分，申請(qǐng)人采用帶有二階Runge-Kutta時(shí)間積分的CIP方法。4.2將CIP與自適應(yīng)八叉樹網(wǎng)格相結(jié)合由Losasso等人[2004]介紹的基于八叉樹數(shù)據(jù)結(jié)構(gòu)的自適應(yīng)網(wǎng)格的使用允許在全部流體中以非齊次精度來執(zhí)行模擬。從實(shí)踐觀點(diǎn)看，這種方法是有用的，因?yàn)槠湓试S通過引入少量的附加計(jì)算和存儲(chǔ)器來以更高的精度模擬諸如表面的更感興趣的區(qū)域。這里，申請(qǐng)人采用這種八叉樹數(shù)據(jù)結(jié)構(gòu)，并提出這種自適應(yīng)方法的實(shí)際值能夠通過將其與CIP方法相結(jié)合而得到進(jìn)一步改善。Losasso等人[2004]對(duì)于半拉格朗日平流步驟使用線性插值。然而，如果申請(qǐng)人用CIP插值來代替線性插值，則平流將具有三階精度。如在Guendelman等人[2005]中一樣，申請(qǐng)人在單元中心處對(duì)壓力值進(jìn)行取樣，但是在結(jié)點(diǎn)處對(duì)所有其它量—速度、水平集及其導(dǎo)數(shù)—進(jìn)行取樣。當(dāng)單元被細(xì)化時(shí)，通過參照所有導(dǎo)數(shù)值來對(duì)新網(wǎng)格點(diǎn)處的值進(jìn)行CIP插值。我們注意到CIP方法非常好地與八叉樹數(shù)據(jù)結(jié)構(gòu)相配合，因?yàn)殡m然其具有三階空間精度，但其被限定在單一網(wǎng)格單元模板(stencil)而不是多個(gè)模板上。由于這種緊湊性，申請(qǐng)人能夠在沒有任何主要修改的情況下將CIP方法用于模擬自適應(yīng)網(wǎng)格。相反，對(duì)于被限定在多個(gè)模板上的高階方案，難以將該方法從一階擴(kuò)展到三階由于模擬自適應(yīng)地細(xì)化網(wǎng)格，所以將產(chǎn)生具有不同尺寸的單元，這使得多模板高階方案的開發(fā)成為難題。我們還將注意到一種申請(qǐng)人稱為八叉樹人為瑕疵的特征。八叉樹數(shù)據(jù)結(jié)構(gòu)的網(wǎng)格分辨率的區(qū)域性變化產(chǎn)生耗散的量的變化。質(zhì)量耗散的區(qū)域性變化不顯著。然而，對(duì)于速度耗散，區(qū)域性差異可能顯著，尤其是對(duì)于快速流體運(yùn)動(dòng)。例如，在潰壩的模擬中，水經(jīng)歷快速的橫向移動(dòng)。接近分辨率高的表面的流體經(jīng)歷少量的數(shù)值擴(kuò)散并因此而進(jìn)行快速移動(dòng)。相反，底部的流體由于低的網(wǎng)格分辨率而以更粘流體的特征的方式移動(dòng)。當(dāng)這兩類運(yùn)動(dòng)一起出現(xiàn)時(shí)，水的上部好像是在下部上爬行。包括CIP在內(nèi)的高階方案避免不了這種人為瑕疵。然后，當(dāng)使用本工作中提出的基于八叉樹的CIP解算器來模擬平流時(shí)，所述人為瑕疵較不顯著；申請(qǐng)人將這種改善歸因于數(shù)值耗散的總體降低，尤其是在低分辨率區(qū)域中。5導(dǎo)數(shù)質(zhì)點(diǎn)模型在動(dòng)量/質(zhì)量守恒方面，基于質(zhì)點(diǎn)的拉格朗日方案比基于網(wǎng)格的歐拉方案更精確。因此，申請(qǐng)人將拉格朗日方案應(yīng)用于需要更詳細(xì)地模擬的區(qū)域。然而，這種方法在表面追蹤和壓力/粘度計(jì)算方面具有其本身的局限性。因此，申請(qǐng)人僅將該方法應(yīng)用于平流部分的模擬；模擬的所有其它部分均基于歐拉方案。拉格朗日與歐拉方案之間的切換要求物理量(速度和水平集值)的網(wǎng)格至質(zhì)點(diǎn)和質(zhì)點(diǎn)至網(wǎng)格轉(zhuǎn)換。應(yīng)小心地開發(fā)用于執(zhí)行這些轉(zhuǎn)換的程序以便它們不引入不必要的數(shù)值擴(kuò)散。在本部分中，申請(qǐng)人介紹了不削弱拉格朗日平流的需要的特性的新型轉(zhuǎn)換程序，其是使得高雷諾數(shù)流體中的小尺度特征能夠再現(xiàn)的必要的組件。為簡(jiǎn)單起見，對(duì)2D進(jìn)行本部分的說明。5.1網(wǎng)格至質(zhì)點(diǎn)速度轉(zhuǎn)換在圖8中所示的平流部分的末端，(1)需要被平流輸送的速度存儲(chǔ)在網(wǎng)格點(diǎn)處，而且(2)大量質(zhì)點(diǎn)散布在網(wǎng)格上。質(zhì)點(diǎn)平流部分必須找出質(zhì)點(diǎn)的速度從而其速度和水平集值產(chǎn)生拉格朗日運(yùn)動(dòng)。假設(shè)質(zhì)點(diǎn)P處于由四個(gè)網(wǎng)格點(diǎn)G1、G2、G3與G4限定的單元中。對(duì)于i∈{1，2，3，4}，使uGi＝(uGi，vGi)為Gi處的網(wǎng)格速度。申請(qǐng)人正在尋找能夠用于P的速度uP＝(uP，vP)的公式。在單元中質(zhì)點(diǎn)方法[Harlow1963]中使用的一種可能方法是使用雙線性插值其中，wi是基于質(zhì)點(diǎn)位置P而確定的雙線性加權(quán)。不幸的是，該轉(zhuǎn)換引入嚴(yán)重的數(shù)值擴(kuò)散。在FLIP方法[Brackbill和Ruppel1986]中，僅網(wǎng)格速度的遞增部分被傳遞質(zhì)點(diǎn)速度。即，其中，是就在非平流部分開始之前的質(zhì)點(diǎn)速度，且是由于非平流部分而引起的Gi處的速度變化。這種方法已被示出顯著降低數(shù)值擴(kuò)散[Zhu和Bridson2005]。因此，在網(wǎng)格至質(zhì)點(diǎn)速度轉(zhuǎn)換程序的開發(fā)中，申請(qǐng)人使用FLIP方法，但是用單調(diào)CIP插值[Song等人2005]來代替線性插值。這不僅給出uP，而且給出5.2質(zhì)點(diǎn)至網(wǎng)格速度轉(zhuǎn)換假設(shè)每個(gè)質(zhì)點(diǎn)根據(jù)速度uP移動(dòng)，攜帶速度的空間導(dǎo)數(shù)以及速度本身。申請(qǐng)人現(xiàn)在必須開發(fā)將質(zhì)點(diǎn)平流的結(jié)果傳送到網(wǎng)格的程序。質(zhì)點(diǎn)至網(wǎng)格速度轉(zhuǎn)換比網(wǎng)格至質(zhì)點(diǎn)轉(zhuǎn)換更復(fù)雜。參照?qǐng)D9，使G為網(wǎng)格點(diǎn)。申請(qǐng)人從G的每個(gè)象限選擇最近的質(zhì)點(diǎn)，即P1、P2、P3和P4。使為Pi的速度，并使和分別為的x和y分量的導(dǎo)數(shù)。申請(qǐng)人必須找到G處的速度uG＝(uG，vG)及其空間導(dǎo)數(shù)的公式。由于所述四個(gè)質(zhì)點(diǎn)沒有成矩形地放置，所以申請(qǐng)人不能使用傳統(tǒng)的基于網(wǎng)格的CIP插值。我們的質(zhì)點(diǎn)至網(wǎng)格速度轉(zhuǎn)換由以下步驟組成。為簡(jiǎn)單起見，申請(qǐng)人僅示出x分量的計(jì)算。(1)申請(qǐng)人將的導(dǎo)數(shù)坐標(biāo)變換成，從而，表示沿方向的分量，且u1⊥|表示垂直于的分量。同樣地，申請(qǐng)人得到的坐標(biāo)變換的導(dǎo)數(shù)(2)申請(qǐng)人沿著方向?qū)Σ襟E(1)中得到的結(jié)果執(zhí)行一維單調(diào)CIP插值以便計(jì)算位置A處的uA及其方向的導(dǎo)數(shù)(3)申請(qǐng)人通過以上坐標(biāo)變換的反轉(zhuǎn)來從獲得(uAx，uAy)。(4)同樣地，申請(qǐng)人對(duì)質(zhì)點(diǎn)P3和P4執(zhí)行步驟(1)～(3)以計(jì)算B的uB和(uBx，uBy)。(5)申請(qǐng)人對(duì)步驟(3)和(4)中得到的結(jié)果執(zhí)行y方向單調(diào)CIP插值，這給出G處的速度和導(dǎo)數(shù)。這種單調(diào)插值法能夠直接擴(kuò)展到3D情況。申請(qǐng)人注意到本工作未采用一般徑向基插值方案，因?yàn)?1)他們不利用導(dǎo)數(shù)信息，而且(2)已證明難以修改方案以便它們體現(xiàn)導(dǎo)數(shù)但保持單調(diào)性。5.3水平集轉(zhuǎn)換應(yīng)不同于速度轉(zhuǎn)換來執(zhí)行水平集轉(zhuǎn)換；水平集值表示到表面的最短距離，所以，該換算將不基于插值法。這里，申請(qǐng)人使用廣泛使用的水平集轉(zhuǎn)換程序[Enright等人2002a；Enright等人2002b]的修改的且更準(zhǔn)確的形式。在初始的質(zhì)點(diǎn)水平集方法中，使用球隱函數(shù)φ(x)＝sp(|φp|-|x-xp|)來計(jì)算網(wǎng)格點(diǎn)x處的水平集值，其中，p是質(zhì)點(diǎn)的水平集，且對(duì)于正(負(fù))質(zhì)點(diǎn)，sp＝+1(-1)。在本工作中，利用存儲(chǔ)在導(dǎo)數(shù)質(zhì)點(diǎn)中的導(dǎo)數(shù)信息，申請(qǐng)人用下式來計(jì)算該網(wǎng)格點(diǎn)處的水平集及其導(dǎo)數(shù)其中p是由導(dǎo)數(shù)質(zhì)點(diǎn)攜帶的梯度。由于方程(8)是基于梯度方向的距離而不是歐氏距離，所以其產(chǎn)生比球隱函數(shù)更準(zhǔn)確的結(jié)果。除以上修改之外，該轉(zhuǎn)換程序與Enright等人[2002a]中提出的相同。6試驗(yàn)結(jié)果本發(fā)明中提出的技術(shù)在具有雙G52.5GHz處理器和5.5GB的存儲(chǔ)器的PowerMac上實(shí)施。申請(qǐng)人使用程序來模擬真實(shí)世界中產(chǎn)生高雷諾數(shù)流體行為的多種實(shí)驗(yàn)情況。在實(shí)驗(yàn)中，申請(qǐng)人用以下常數(shù)來求解由水和空氣組成的兩相流體g＝9.8m/sec2、ρwater＝1000kg/m3、μwater＝1.137×103kg/ms、ρa(bǔ)ir＝1.226kg/m3、以及μair＝1.78×105kg/ms，其中g(shù)是重力。使用移動(dòng)立方體算法來進(jìn)行水表面的提取，而且通過智能射來進(jìn)行渲染。潰壩為了比較導(dǎo)數(shù)質(zhì)點(diǎn)模型與用質(zhì)點(diǎn)水平集方法的線性模型之間的數(shù)值擴(kuò)散，申請(qǐng)人用有效分辨率1282來執(zhí)行2D潰壩測(cè)試，其中，在重力場(chǎng)釋放0.2×0.4m的水柱。圖10示出了結(jié)果，申請(qǐng)人能夠看到導(dǎo)數(shù)質(zhì)點(diǎn)產(chǎn)生擴(kuò)散較少的結(jié)果波的斷裂更劇烈，而且很好地保存渦度。落到淺水上的大塊水圖11示出了從下述模擬截取的截圖即其中大塊水落到0.05m深的蓄水池上。申請(qǐng)人對(duì)于蓄水池的底表面使用無滑動(dòng)邊界條件，這引起水在頂部快速地移動(dòng)而在底部慢速地移動(dòng)，導(dǎo)致皇冠狀飛濺。這項(xiàng)實(shí)驗(yàn)用2563的有效網(wǎng)格分辨率來執(zhí)行。還用使用質(zhì)點(diǎn)水平集方法的傳統(tǒng)線性模型來模擬相同的情景。比較證明所提出的技術(shù)產(chǎn)生更真實(shí)的結(jié)果。導(dǎo)數(shù)質(zhì)點(diǎn)模型平均每時(shí)間步長用60秒，而線性模型平均每時(shí)間步長用30秒，均包括文件輸出時(shí)間。固體球的沖擊在這項(xiàng)實(shí)驗(yàn)中，具有半徑0.15m的固體球以速度(5.0，-3.0，0.0)m/sec沖擊到水中。該沖擊產(chǎn)生圖7中所示的復(fù)雜結(jié)構(gòu)。這項(xiàng)實(shí)驗(yàn)表明所提出的技術(shù)能夠產(chǎn)生猛烈的固體-水-空氣相互作用中發(fā)生的詳細(xì)的流體運(yùn)動(dòng)和表面特征。這項(xiàng)實(shí)驗(yàn)的有效分辨率是192×128×128。旋轉(zhuǎn)水箱圖12示出了具有半滿的水的0.9×0.9×0.3m旋轉(zhuǎn)箱。在實(shí)驗(yàn)中，離心力產(chǎn)生將水推向水箱側(cè)的效果。這種模擬以96×96×32的有效網(wǎng)格分辨率來執(zhí)行。這項(xiàng)實(shí)驗(yàn)證明即使在粗分辨率的網(wǎng)格中，導(dǎo)數(shù)質(zhì)點(diǎn)模型也能夠模擬流體的復(fù)雜運(yùn)動(dòng)。洪水在這項(xiàng)實(shí)驗(yàn)中，申請(qǐng)人模擬在6m×2m×4m域上發(fā)生的中等尺度的洪水。圖13示出了從這項(xiàng)實(shí)驗(yàn)截取的截圖。當(dāng)移除堤壩時(shí)，水開始滑落臺(tái)階。由于沿路的阻礙，水進(jìn)行產(chǎn)生復(fù)雜幾何形狀的猛烈運(yùn)動(dòng)。這項(xiàng)實(shí)驗(yàn)的有限分辨率是384×128×256。圖1～6示出了本發(fā)明的流程圖。如圖1～6中所示，一種用于使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法包括以下步驟a)用基于八叉樹數(shù)據(jù)結(jié)構(gòu)的自適應(yīng)網(wǎng)格建模流體(S100)；b)對(duì)于網(wǎng)格點(diǎn)處的流體速度求解不可壓縮流體的納維爾-斯托克斯方程(S200)；c)在納維爾-斯托克斯方程的時(shí)間積分中使用分步法(S300)；d)將分步分組成非平流部分和平流部分(S400)；e)根據(jù)模擬的詳細(xì)程度來選擇網(wǎng)格平流部分或質(zhì)點(diǎn)平流部分(S500)；f)對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP方法(S600)；以及g)對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型(S700)。當(dāng)模擬的詳細(xì)程度高時(shí)使用對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟(S700)。對(duì)于高度感興趣的區(qū)域的細(xì)節(jié)自適應(yīng)網(wǎng)格增加網(wǎng)格的數(shù)目。八叉樹數(shù)據(jù)結(jié)構(gòu)包括樹形數(shù)據(jù)結(jié)構(gòu)，每個(gè)內(nèi)部結(jié)點(diǎn)具有多達(dá)八個(gè)孩子。不可壓縮流體的納維爾-斯托克斯方程可以寫成其中u是速度，p是壓力，f是外力，μ是粘度系數(shù)，ρ是流體密度。該方法可以進(jìn)一步包括使用水平集方法來追蹤流體的表面的步驟(S800)。如圖1至3中所示，所述水平集方法包括水平集場(chǎng)φ。根據(jù)下式來更新所述水平集場(chǎng)φ符號(hào)距離的時(shí)間積分僅涉及平流部分。網(wǎng)格平流部分平流輸送根據(jù)當(dāng)前速度場(chǎng)從非平流部分計(jì)算的速度和水平集場(chǎng)。網(wǎng)格平流部分使用歐拉平流步驟來平流輸送速度和水平集。通過基于八叉樹的CIP(約束插值剖面)方法來執(zhí)行所述歐拉平流。對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP方法的步驟(S600)包括用基于八叉樹CIP方法來平流輸送速度和水平集及其導(dǎo)數(shù)的步驟(S610)。直接從方程獲得用于平流輸送導(dǎo)數(shù)的方程，其中，φ是被平流輸送的物理量，且通過對(duì)于x求方程的微分來獲得空間變量x的導(dǎo)數(shù)的平流輸送方程，即如圖2和圖3中所示，求解x的導(dǎo)數(shù)的平流輸送方程的步驟(S610)包括采用使用分步法的步驟，且所述使用分步法的步驟包括以下步驟a)使用有限差分法來求解非平流部分(S620)；b)平流輸送根據(jù)的結(jié)果(S630)；以及c)在由網(wǎng)格限定的單元的中心處進(jìn)行壓力值取樣并在結(jié)點(diǎn)處進(jìn)行包括速度、水平集及其導(dǎo)數(shù)的其它值進(jìn)行取樣(S640)。如圖3中所示，該方法可以進(jìn)一步包括以下步驟a)將基于八叉樹的CIP法與自適應(yīng)網(wǎng)格相結(jié)合(S650)；以及b)通過用基于八叉樹的CIP法模擬平流來使由八叉樹數(shù)據(jù)結(jié)構(gòu)的網(wǎng)格分辨率中的區(qū)域性變化引起的八叉樹人為瑕疵最小化(S660)。導(dǎo)數(shù)質(zhì)點(diǎn)模型使用基于質(zhì)點(diǎn)的拉格朗日方案。用和來計(jì)算網(wǎng)格點(diǎn)處的水平集及其導(dǎo)數(shù)，其中p是由導(dǎo)數(shù)質(zhì)點(diǎn)攜帶的梯度。該方法可以進(jìn)一步包括以下步驟a)在對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟(S700)之前將網(wǎng)格上限定的速度和水平集轉(zhuǎn)換成質(zhì)點(diǎn)上限定的速度和水平集(S690)(網(wǎng)格至質(zhì)點(diǎn))；以及b)在返回到對(duì)于網(wǎng)格平流部分使用基于八叉樹的CIP法的步驟(S600)之前將質(zhì)點(diǎn)上限定的速度和水平集轉(zhuǎn)換成網(wǎng)格上限定的速度和水平集(S710)(質(zhì)點(diǎn)至網(wǎng)格)。轉(zhuǎn)換速度和水平集的步驟(S690)網(wǎng)格至質(zhì)點(diǎn)包括將具有用單調(diào)CIP插值代替線性插值的FLIP方法用于速度和uP的轉(zhuǎn)換的步驟(S692)，其中，是就在非平流部分開始之前的質(zhì)點(diǎn)速度，且是由于非平流部分而引起的Gi處的速度變化。如圖5中所示，對(duì)于速度uP、網(wǎng)格點(diǎn)G、選自G的每個(gè)象限的最近質(zhì)點(diǎn)P1、P2、P3和P4、Pi的速度、以及和的x和y分量的導(dǎo)數(shù)，轉(zhuǎn)換速度和水平集的步驟(S710)質(zhì)點(diǎn)至網(wǎng)格包括以下步驟a)將的導(dǎo)數(shù)坐標(biāo)變換成因此表示沿著方向的分量且u1⊥表示垂直于的分量，其中，類似地獲得的坐標(biāo)變換的導(dǎo)數(shù)(S712)；b)沿著方向?qū)Σ襟E(a)中獲得的結(jié)果進(jìn)行一維單調(diào)CIP插值以便計(jì)算位置A處的uA及其方向的導(dǎo)數(shù)(S714)；c)通過應(yīng)用以上坐標(biāo)變換的反轉(zhuǎn)來從獲得(uAx，uAy)(S716)；d)對(duì)質(zhì)點(diǎn)P3和P4執(zhí)行步驟(a)～(c)以計(jì)算B的uB和(uBx，uBy)(S718)；以及e)對(duì)步驟(c)和(d)中獲得的結(jié)果執(zhí)行y方向單調(diào)CIP插值，其給出G處的速度和導(dǎo)數(shù)(S720)。如圖6中所示，對(duì)于質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟包括以下步驟a)調(diào)整質(zhì)點(diǎn)的水平集值(S702)；以及b)執(zhí)行質(zhì)點(diǎn)重新播種(S704)。7結(jié)論在本發(fā)明中，申請(qǐng)人已介紹了一種稱為導(dǎo)數(shù)質(zhì)點(diǎn)的新概念，并提出了增加平流部分的精度的流體模擬技術(shù)。使用歐拉方案來實(shí)施模擬器的非平流部分，而使用拉格朗日方案來實(shí)施平流部分。在開發(fā)以這種方式結(jié)合的模擬器時(shí)的重要問題是如何轉(zhuǎn)換平流與非平流部分的切換處的物理量。申請(qǐng)人通過開發(fā)有效地抑制不必要的數(shù)值耗散的轉(zhuǎn)換程序而成功地克服了這個(gè)問題。在不需要執(zhí)行質(zhì)點(diǎn)平流的流體區(qū)域中，使用基于八叉樹的CIP解算器來模擬平流，這是申請(qǐng)人通過將CIP插值與八叉樹數(shù)據(jù)結(jié)構(gòu)相結(jié)合而在本工作中開發(fā)的。這種方法還用來降低數(shù)值耗散。多次實(shí)驗(yàn)的結(jié)果證明提出的技術(shù)顯著降低了速度耗散，導(dǎo)致動(dòng)態(tài)流體的真實(shí)再現(xiàn)。沒有先前的研究人員的開拓性工作，本工作將是不可能的。完全使用質(zhì)點(diǎn)來求解平流部分的思想來自PIC[Harlow1963]和FLIP[Brackbill和Ruppel1986]方法，但是對(duì)于質(zhì)點(diǎn)至網(wǎng)格和網(wǎng)格至質(zhì)點(diǎn)速度換算，申請(qǐng)人使用基于導(dǎo)數(shù)的三次插值來代替線性插值。另外，讓質(zhì)點(diǎn)攜帶水平集值的想法來自質(zhì)點(diǎn)水平集方法[Enright等人2002b]。這種模型在本工作中得到擴(kuò)展，從而，質(zhì)點(diǎn)還攜帶速度信息。此外，利用導(dǎo)數(shù)信息來自CIP方法[Yabe等人2001]。這里，申請(qǐng)人將這種方法的應(yīng)用擴(kuò)展至質(zhì)點(diǎn)，引起質(zhì)點(diǎn)至網(wǎng)格速度轉(zhuǎn)換程序的開發(fā)，該轉(zhuǎn)換程序包含CIP方法本身的擴(kuò)展非矩形構(gòu)造中的質(zhì)點(diǎn)的CIP插值。在實(shí)驗(yàn)中，強(qiáng)調(diào)再現(xiàn)細(xì)節(jié)的能力。但是，提出的技術(shù)還能夠應(yīng)用于大尺度的流體運(yùn)動(dòng)以產(chǎn)生準(zhǔn)確的結(jié)果。該技術(shù)包括增加的用于存儲(chǔ)導(dǎo)數(shù)信息的存儲(chǔ)器的量。然而，申請(qǐng)人證明能夠在當(dāng)前的PC上模擬2563網(wǎng)格，這已經(jīng)是用于描繪感興趣的流體情景的可用的分辨率。雖然已參照不同的實(shí)施例示出并描述了本發(fā)明，但本領(lǐng)域的技術(shù)人員將認(rèn)識(shí)到在不脫離權(quán)利要求書所限定的本發(fā)明的精神和范圍的情況下可以進(jìn)行形式、細(xì)節(jié)、組成以及操作上修改。參考文獻(xiàn)BRACKBILL，J.U.，ANDRUPPEL，H.M.1986.FlipAmethodforadaptivelyzoned，particle-in-cellcalculationsoffluidflowsintwodimensions.J.Comp.Phys.65，314-343.CARLSON，M.，MUCHA，R.J.，ANDTURK，G.2004.RigidfluidAnimatingtheinterplaybetweenrigidbodiesandfluid.ACMTransactionsonGraphics23，3，377-384.ENRIGHT，D.，F(xiàn)EDKIW，R.，F(xiàn)ERZIGER，J.，ANDMITCHELL，I.2002.Ahybridparticlelevelsetmethodforimprovedinterfacecapturing.J.Comp.Phys.183，83-116.ENRIGHT，D.，MARSCHNER，S.，ANDFEDKIW，R.2002.Animationandrenderingofcomplexwatersurfaces.ACMTransactionsonGraphics21，3，736-744.ENRIGHT，D.，LOSASSO，F(xiàn).，ANDFEDKIW，R.2005.Afastandaccuratesemi-lagrangianparticlelevel′setmethod.ComputersandStructures83，479-490.FEDKIW，R.，STAM，J.，ANDJENSEN，H.W.2001.Visualsimulationofsmoke.ComputerGraphics(Proc.ACMSIGGRAPH2001)35，15-22.FELDMAN，B.E.，O′BRIEN，J.F.，ANDARIKAN，0.2003.Animatingsuspendedparticleexplosions.ACMTransactionsonGraphics22，3，708-715.FELDMAN，B.E.，O′BRIEN，J.F.，ANDKLINGNER，B.M.2005.Animatinggaseswithhybridmeshes.ACMTransactionsonGraphics24，3，904-909.FELDMAN，B.E.，O′BRIEN，J.F.，KLINGNER，B.M.，ANDGOKTEKIN，T.G.2005.Fluidsindeformingmeshes.InProceedingsofEurographics/ACMSIGGRAPHSymposiumonComputerAnimation2005，255-259.FOSTER，N.，ANDFEDKIW，R.2001.Practicalanimationofliquids.ComputerGraphics(Proc.ACMSIGGRAPH2001)35，23-30.FOSTER，N.，ANDMETAXAS，D.1996.Realisticanimationofliquids.GraphicalmodelsandimageprocessingGMIP58，5，471-483.FOSTER，N.，ANDMETAXAS，D.1997.Controllingfluidanimation.InComputerGraphicsInternational97，178-188.GOKTEKIN，T.G.，BARGTEIL，A.W.，ANDO′BRIEN，J.F.2004.Amethodforanimatingviscoelasticfluids.ACMTransactionsonGraphics23，3，463-468.GUENDELMAN，E.，SELLE，A.，LOSASSO，F(xiàn).，ANDFEDKIW，R.2005.Couplingwaterandsmoketothindeformableandrigidshells.ACMTransactionsonGraphics24，3，973-981.HARLOW，F(xiàn).H.1963.Theparticle-in-cellmethodfornumericalsolutionofproblemsinfluiddynamics.InExperimentalarithmetic，high-speedcomputationsandmathematics.HONG，J.-M.，ANDKIM，C-H.2005.Discontinuousfluids.ACMTransactionsonGraphics24，3，915-920.KANG，M.，F(xiàn)EDKIW，R.，ANDLIU，X.-D.2000.Aboundary-conditioncapturingmethodformultiphaseincompressibleflow.J.Sci.Comput.15，323-360.LIU，X.-D.，F(xiàn)EDKIW，R.，ANDKANG，M.2000.Aboundaryconditioncapturingmethodforpoisson′sequationonirregulardomains.J.Comp.Phys.160，151-178.LOSASSO，F(xiàn).，GIBOU，F(xiàn).，ANDFEDKIW，R.2004.Simulatingwaterandsmokewithanoctreedatastructure.ACMTransactionsonGraphics23，3，457-462.MCNAMARA，A.，TREUILLE，A.，POPOVI′C，Z.，ANDSTAM，J.2004.Fluidcontrolusingtheadjointmethod.ACMTransactionsonGraphics23，3，449-456.M＂ULLER，M.，CHARYPAR，D.，ANDGROSS，M.2003.Particlebasedfluidsimulationforinteractiveapplications.InProceedingsofACMSIGGRAPH/EurographicsSymposiumonComputerAnimation2003，154-159.M＂ULLER，M.，SOLENTHALER，B.，KEISER，R.，ANDGROSS，M.2005.Particle-basedfluid-fluidinteraction.InProceedingsofEurographics/ACMSIGGRAPHSymposiumonComputerAnimation2005，237-244.PARK，S.I.，ANDKIM，M.J.2005.Vortexfluidforgaseousphenomena.InProceedingsofEurographics/ACMSIGGRAPHSymposiumonComputerAnimation2005，261-270.PREMCTZE，S.，TASDIZEN，T.，BIGLER，J.，LEFOHN，A.，ANDWHITAKER，R.T.2003.Particle-basedsimulationoffluids.InEurographics2003proceedings，BlackwellPublishers，401-410.RASMUSSEN，N.，NGUYEN，D.Q.，GEIGER，W.，ANDFEDKIW，R.2003.Smokesimulationforlargescalephenomena.ACMTransactionsonGraphics22，3，703-707.SELLE，A.，RASMUSSEN，N.，ANDFEDKIW，R.2005.Avortexparticlemethodforsmoke，waterandexplosions.ACMTransactionsonGraphics24，3，910-914.SONG，O.-Y.，SHIN，H.，ANDKO，H.-S.2005.Stablebutnondissipativewater.ACMTransactionsonGraphics24，1，81-97.STAM，J.，ANDFIUME，E.1995.Depictingfireandothergaseousphenomenausingdiffusionprocesses.ComputerGraphics(Proc.ACMSIGGRAPH′95)29，AnnualConferenceSeries，129-136.STAM，J.1999.Stablefluids.ComputerGraphics(Proc.ACMSIGGRAPH′99)33，AnnualConferenceSeries，121-128.STANIFORTH，A.，ANDCλOTTE，J.1991.Semi-lagrangianintegrationschemeforatmosphericmodel-areview.Mon.WeatherRev.119，12，2206-2223.TREUILLE，A.，MCNAMARA，A.，POPOVI′C，Z.，ANDSTAM，J.2003.Keyframecontrolofsmokesimulations.ACMTransactionsonGraphics22，3，716-723.WANG，H.，MUCHA，P.J.，ANDTURK，G.2005.Waterdropsonsurfaces.ACMTransactionsonGraphics24，3，921-929.YABE，T.，ANDAOKI，T.1991.Auniversalsolverforhyperbolicequationsbycubic-polynomialinterpolationi.one-dimensionalsolver.Comp.Phys.Comm.66，219-232.YABE，T.，ANDAOKI，T.1991.Auniversalsolverforhyperbolicequationsbycubic-polynomialinterpolationii.two-andthreedimensionalsolvers.Comp.Phys.Comm.66，233-242.YABE，T.，XIAO，F(xiàn).，ANDUTSUMI，T.2001.Theconstrainedinterpolationprofilemethodformultiphaseanalysis.J.Comp.Phys.169，556-593.ZHU，Y.，ANDBRIDSON，R.2005.Animatingsandasafluid.ACMTransactionsonGraphics24，3，965-972.權(quán)利要求1.一種使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體的詳細(xì)運(yùn)動(dòng)的方法，包括以下步驟a)用基于八叉樹數(shù)據(jù)結(jié)構(gòu)的自適應(yīng)網(wǎng)格來建模流體；b)對(duì)于網(wǎng)格點(diǎn)處的流體速度求解用于不可壓縮流體的納維爾-斯托克斯方程；c)在納維爾-斯托克斯方程的時(shí)間積分中使用分步法；d)將各分步分組成非平流部分和平流部分；e)根據(jù)模擬的詳細(xì)程度來選擇網(wǎng)格平流部分或質(zhì)點(diǎn)平流部分；f)對(duì)網(wǎng)格平流部分使用基于八叉樹的CIP方法；以及g)對(duì)質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型，其中，當(dāng)模擬的詳細(xì)程度高時(shí)，使用所述的對(duì)質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟。2.根據(jù)權(quán)利要求1所述的方法，其中，對(duì)于在細(xì)節(jié)上高度關(guān)注的區(qū)域，所述自適應(yīng)網(wǎng)格增加網(wǎng)格的數(shù)目。3.根據(jù)權(quán)利要求1所述的方法，其中，所述八叉樹數(shù)據(jù)結(jié)構(gòu)包括樹形數(shù)據(jù)結(jié)構(gòu)，其中，每個(gè)內(nèi)部結(jié)點(diǎn)具有多達(dá)八個(gè)的孩子。4.根據(jù)權(quán)利要求1所述的方法，其中，所述用于不可壓縮流體的納維爾-斯托克斯方程能夠?qū)懗善渲?，u是速度，p是壓力，f是外力，μ是粘度系數(shù)，ρ是流體密度。5.根據(jù)權(quán)利要求1所述的方法，還包括使用水平集方法來追蹤流體的表面的步驟，其中，所述水平集方法包括水平集場(chǎng)φ。6.根據(jù)權(quán)利要求5所述的方法，其中，根據(jù)下式來更新所述水平集場(chǎng)φ|帶有符號(hào)距離7.根據(jù)權(quán)利要求6所述的方法，其中，所述的時(shí)間積分僅涉及平流部分。8.根據(jù)權(quán)利要求6所述的方法，其中，所述網(wǎng)格平流部分根據(jù)當(dāng)前的速度場(chǎng)來平流輸送根據(jù)非平流部分計(jì)算的速度和水平集場(chǎng)。9.根據(jù)權(quán)利要求6所述的方法，其中，所述網(wǎng)格平流部分使用歐拉平流步驟來平流輸送速度和水平集。10.根據(jù)權(quán)利要求9所述的方法，其中，通過基于八叉樹的CIP(約束插值剖面)方法來執(zhí)行所述歐拉平流。11.根據(jù)權(quán)利要求1所述的方法，其中，所述的對(duì)網(wǎng)格平流部分使用基于八叉樹的CIP方法的步驟包括用基于八叉樹的CIP方法來平流輸送速度和水平集和它們的導(dǎo)數(shù)的步驟。12.根據(jù)權(quán)利要求11所述的方法，其中，直接從方程獲得用于平流輸送導(dǎo)數(shù)的方程，其中，φ是被平流輸送的物理量，且通過關(guān)于x求方程的微分來獲得空間變量x的導(dǎo)數(shù)的平流輸送方程，即13.根據(jù)權(quán)利要求12所述的方法，其中，求解x的導(dǎo)數(shù)的平流輸送方程的步驟包括使用分步法的步驟，其中，所述使用分步法的步驟包括以下步驟a)使用有限差分法來求解非平流部分b)平流輸送根據(jù)的結(jié)果；以及c)在由網(wǎng)格所限定單元的中心處對(duì)壓力值取樣，并在各結(jié)點(diǎn)處對(duì)包括速度、水平集和它們的導(dǎo)數(shù)的其它值取樣。14.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括步驟a)將基于八叉樹的CIP方法與自適應(yīng)網(wǎng)格相結(jié)合；以及b)通過用基于八叉樹的CIP方法模擬平流來最小化八叉樹數(shù)據(jù)結(jié)構(gòu)的網(wǎng)格分辨率中的區(qū)域性變化所引起的八叉樹人為瑕疵。15.根據(jù)權(quán)利要求1所述的方法，其中，所述導(dǎo)數(shù)質(zhì)點(diǎn)模型使用基于質(zhì)點(diǎn)的拉格朗日方案，其中，用下式來計(jì)算網(wǎng)格點(diǎn)處的水平集及其導(dǎo)數(shù)以及其中，p是導(dǎo)數(shù)質(zhì)點(diǎn)攜帶的梯度。16.根據(jù)權(quán)利要求1所述的方法，進(jìn)一步包括以下步驟a)在所述的對(duì)質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟之前將各網(wǎng)格上定義的速度和水平集轉(zhuǎn)換成各質(zhì)點(diǎn)上定義的速度和水平集(網(wǎng)格到質(zhì)點(diǎn))；以及b)在返回到所述的對(duì)網(wǎng)格平流部分使用基于八叉樹的CIP方法的步驟之前將各質(zhì)點(diǎn)上定義的速度和水平集轉(zhuǎn)換成各網(wǎng)格上定義的速度和水平集(質(zhì)點(diǎn)到網(wǎng)格)。17.根據(jù)權(quán)利要求16所述的方法，其中，網(wǎng)格到質(zhì)點(diǎn)的轉(zhuǎn)換速度和水平集的步驟包括對(duì)于速度uP的轉(zhuǎn)換使用具有用單調(diào)CIP插值代替線性插值的FLIP方法的步驟，其中其中是恰好在非平流部分開端之前的質(zhì)點(diǎn)速度，且是由于非平流部分引起的Gi處的速度變化。18.根據(jù)權(quán)利要求16所述的方法，其中，對(duì)于速度uP、網(wǎng)格點(diǎn)G、選自G的每個(gè)象限的最近質(zhì)點(diǎn)P1、P2、P3和P4、Pi的速度、以及和的x分量和y分量的導(dǎo)數(shù)，質(zhì)點(diǎn)到網(wǎng)格的轉(zhuǎn)換速度和水平集的步驟包括以下步驟a)將的導(dǎo)數(shù)坐標(biāo)變換成(u1‖u1⊥)，使得u1‖表示沿著方向的分量且u1⊥表示垂直于的分量，其中，類似地獲得的坐標(biāo)變換的導(dǎo)數(shù)(u2‖，u2⊥)；b)沿著方向?qū)Σ襟E(a)中獲得的結(jié)果進(jìn)行一維單調(diào)CIP插值以計(jì)算位置A處的uA及其方向的導(dǎo)數(shù)(uA‖，uA⊥)；c)通過應(yīng)用上述坐標(biāo)變換的逆變換，從(uA‖，uA⊥)獲得(uAx，uAy)；d)對(duì)質(zhì)點(diǎn)P3和P4執(zhí)行步驟(a)～(c)以計(jì)算B的uB和(uBx，uBy)，B為網(wǎng)格線和連接P3和P4的線段之間的交點(diǎn)；以及e)對(duì)步驟(c)和(d)中獲得的結(jié)果執(zhí)行y方向單調(diào)CIP插值，給出G處的速度和導(dǎo)數(shù)。19.根據(jù)權(quán)利要求15所述的方法，其中，所述的對(duì)質(zhì)點(diǎn)平流部分使用導(dǎo)數(shù)質(zhì)點(diǎn)模型的步驟包括以下步驟a)調(diào)整質(zhì)點(diǎn)的水平集值；以及b)執(zhí)行質(zhì)點(diǎn)重新播種。全文摘要使用導(dǎo)數(shù)質(zhì)點(diǎn)來模擬流體詳細(xì)運(yùn)動(dòng)的方法。本發(fā)明提供了一種方法，該方法提供了一種使用質(zhì)點(diǎn)和導(dǎo)數(shù)信息來顯著降低速度的非物理耗散的新流體模擬技術(shù)。在求解傳統(tǒng)納維爾-斯托克斯方程的過程中，所述方法用質(zhì)點(diǎn)模擬代替平流部分。當(dāng)在基于網(wǎng)格和基于質(zhì)點(diǎn)的模擬器之間切換時(shí)，必須轉(zhuǎn)換諸如水平集和速度的物理量。開發(fā)了一種利用存儲(chǔ)在質(zhì)點(diǎn)中以及網(wǎng)格點(diǎn)中的導(dǎo)數(shù)信息的新型耗散抑制轉(zhuǎn)換程序。通過多次實(shí)驗(yàn)，提出的技術(shù)可以再現(xiàn)諸如液滴/泡、薄水膜以及漩渦等高雷諾數(shù)流體的詳細(xì)的運(yùn)動(dòng)。平流中增加的精度還能夠用來在更大尺度的流體模擬中產(chǎn)生更好的結(jié)果，所述平流中增加的精度形成提出的技術(shù)的基礎(chǔ)。文檔編號(hào)G06T17/00GK101484922SQ200680054876公開日2009年7月15日申請(qǐng)日期2006年10月27日優(yōu)先權(quán)日2006年4月5日發(fā)明者宋五泳,金度燁,高亨錫申請(qǐng)人:財(cái)團(tuán)法人Seoul大學(xué)校產(chǎn)學(xué)協(xié)力財(cái)團(tuán)