偏微分方程數(shù)值計(jì)算在虛擬現(xiàn)實(shí)中應(yīng)用與研究

鄒長軍，尹 勇，李海江，唐 皇

(大連海事大學(xué) 航海動(dòng)態(tài)仿真和控制實(shí)驗(yàn)室，遼寧 大連 116021)

偏微分方程數(shù)值計(jì)算在虛擬現(xiàn)實(shí)中應(yīng)用與研究

鄒長軍，尹 勇，李海江，唐 皇

(大連海事大學(xué) 航海動(dòng)態(tài)仿真和控制實(shí)驗(yàn)室，遼寧 大連 116021)

隨著虛擬現(xiàn)實(shí)（Virtual Reality，VR）相關(guān)硬件的發(fā)展以及軟件技術(shù)的日趨成熟，VR技術(shù)正在進(jìn)入人們的日常生活。然而不可否認(rèn)，VR的核心問題之一——物理真實(shí)感依然沒有很好的解決。其原因大多都是由于解算模型控制方程耗費(fèi)了大量的時(shí)間，從而迫使人們采用犧牲精度的方法以滿足實(shí)時(shí)性的要求。而許多的物理模型的解算都?xì)w結(jié)于偏微分方程求解，因此如何準(zhǔn)確高效快速的求解偏微分方程（PDE）和PDE方程組對(duì)于提升虛擬現(xiàn)實(shí)系統(tǒng)的物理真實(shí)感有著至關(guān)重要的作用。本文從求解典型PDE出發(fā)，建立相應(yīng)的PDE求解器，并對(duì)計(jì)算結(jié)果進(jìn)行驗(yàn)證。最后將該思想應(yīng)用于計(jì)算流體力學(xué)領(lǐng)域，通過對(duì)流體力學(xué)中的NS方程和水力學(xué)中的淺水波方程的求解，得到計(jì)算區(qū)域的速度場(chǎng)或高度場(chǎng)，并且對(duì)結(jié)果進(jìn)行驗(yàn)證。結(jié)果表明，該方法具有較高的可信度。

虛擬現(xiàn)實(shí)；PDE；數(shù)值計(jì)算；物理真實(shí)感

0 引 言

虛擬現(xiàn)實(shí)（Virtual Reality，VR）技術(shù)是一種可以創(chuàng)建和體驗(yàn)虛擬世界的計(jì)算機(jī)仿真系統(tǒng)，它利用計(jì)算機(jī)生成一種模擬環(huán)境，是一種多源信息融合的交互式的三維動(dòng)態(tài)視景和實(shí)體行為的系統(tǒng)仿真，并使用戶沉浸到該環(huán)境中。近年來VR軟件硬件技術(shù)得到了長足的發(fā)展，消費(fèi)級(jí)的VR可穿戴設(shè)備也已經(jīng)問世。如Oculus Rift DK2、索尼Project Morpheus、三星Gear VR等都是目前較好的虛擬現(xiàn)實(shí)裝置。2016年2月教育部正式批準(zhǔn)了100個(gè)國家級(jí)虛擬仿真實(shí)驗(yàn)教學(xué)中心的建設(shè)。作者認(rèn)為此舉意義非凡，不但是從國家層面對(duì)VR技術(shù)發(fā)展趨勢(shì)的肯定，同時(shí)必將大大推進(jìn)VR技術(shù)在教育及其他相關(guān)領(lǐng)域的大發(fā)展。

與此同時(shí)，2016年03月05日德意志銀行發(fā)布了VR發(fā)展報(bào)告稱2020年VR市場(chǎng)規(guī)模將達(dá)70億美元，認(rèn)為虛擬現(xiàn)實(shí)如2007年的智能手機(jī)，一切都在有條不紊地進(jìn)行著。報(bào)告還指出，隨著大眾的需求遠(yuǎn)超預(yù)期，虛擬現(xiàn)實(shí)顯得越發(fā)真實(shí)，盡管這項(xiàng)技術(shù)會(huì)進(jìn)一步占據(jù)市場(chǎng)，但挑戰(zhàn)依然存在，而包括Facebook、谷歌在內(nèi)的科技巨頭都明確表示將繼續(xù)探索。德意志銀行最后再次強(qiáng)調(diào)，進(jìn)一步市場(chǎng)化的關(guān)鍵并不是硬件，而是內(nèi)容。其中最重要的是，創(chuàng)建沉浸式的內(nèi)容需要合適的故事情節(jié)，而這正是目前的一大軟肋。

該報(bào)告從一方面肯定了VR大發(fā)展的時(shí)代即將到來，強(qiáng)調(diào)了VR的核心在于內(nèi)容，而不是硬件設(shè)備，同時(shí)也不否認(rèn)VR技術(shù)在進(jìn)入人們的日常生活還有很長的一段路要走。這其中，VR技術(shù)中物理真實(shí)感的嚴(yán)重不足就是制約VR發(fā)展的關(guān)鍵因素之一。在創(chuàng)建高質(zhì)量VR內(nèi)容的時(shí)候，除了高沉浸感之外，人們還需要考慮場(chǎng)景中的物理真實(shí)感。目前國內(nèi)外眾多的學(xué)者就物理真實(shí)感及相關(guān)課題展開了廣泛而深入的研究。如William Mell[1]建立了基于N-S方程和燃燒控制方程的草場(chǎng)火傳播的數(shù)學(xué)模型，該模型考慮了風(fēng)、火源寬度等因素對(duì)火場(chǎng)擴(kuò)散的影響，并通過實(shí)驗(yàn)對(duì)進(jìn)行了驗(yàn)證模型。柳有權(quán)[2]采用GPU加速策略，對(duì)帶有復(fù)雜邊界的流體進(jìn)行實(shí)時(shí)三維模擬，該算法對(duì)于實(shí)時(shí)流體模擬具有重要的意義。王長波[3]針對(duì)自由表面流體的模擬，提出一種基于Lattice Boltzmann（LBM）的高效建模和繪制的方法。其基于淺水方程的LBM模型進(jìn)行流體建模及表面高度場(chǎng)計(jì)算，并提出一種基于Marching Cubes和自由表面算法結(jié)合的方法來抽取流體的表面，隨后采用考慮移動(dòng)障礙物的外力疊加機(jī)制和自適應(yīng)加密算法來進(jìn)行流體交互及表面的網(wǎng)格重構(gòu)。最后采用硬件加速技術(shù)實(shí)現(xiàn)了不同自由表面流體的繪制，如溪流、水池淺水流、洪水水淹等真實(shí)感效果。王順利[4]針對(duì)以往艦首浪三維可視化中僅從實(shí)際觀察出發(fā)建立動(dòng)態(tài)模型，不能反映艦首浪物理運(yùn)動(dòng)規(guī)律的問題，提出一種基于物理模型的艦首浪三維可視化方法。該方法將艦首浪物理模型應(yīng)用到三維可視化中，采用邊界元方法計(jì)算得到艦首浪外形數(shù)據(jù)，利用粒子系統(tǒng)技術(shù)建立了三維動(dòng)態(tài)模型。JIANG[5]采用質(zhì)點(diǎn)方法對(duì)流體和固體的模擬。采用MPM（Material Point Method）方法 建立固體的變形模型，并針對(duì)熱傳導(dǎo)，固體融化和凝固過程進(jìn)行模擬，取得了較好的可視化效果。

從已有的研究中可看出，虛擬現(xiàn)實(shí)技術(shù)的物理真實(shí)感一直都是一項(xiàng)具有挑戰(zhàn)性的課題。

通過分析不難發(fā)現(xiàn)，影響物理真實(shí)感的因素主要包括：物理模型近似程度，數(shù)值計(jì)算方法的精度，數(shù)值計(jì)算的效率，硬件的性能等，渲染引擎的效率等。為了獲得更加真實(shí)的物理真實(shí)感，作者從數(shù)值計(jì)算方法的精度和通用性方面進(jìn)行了研究。針對(duì)PDE模型方程的時(shí)間導(dǎo)數(shù)項(xiàng)、梯度項(xiàng)、散度項(xiàng)及旋度項(xiàng)等，分別采用相應(yīng)的數(shù)值離散算子。這樣，在解算不同PDE時(shí)，可以采用統(tǒng)一的離散算子接口，從而實(shí)現(xiàn)PDE解算方法的快速構(gòu)建，并且其中的離散方法可以根據(jù)用戶的需要進(jìn)行擴(kuò)充和完善。這樣可以實(shí)現(xiàn)離散算子和求解器的重復(fù)利用，可以不必反復(fù)測(cè)試，大大節(jié)省人力物力成本。

本文首先對(duì)有限體積方法及其在PDE中的應(yīng)用進(jìn)行介紹，然后從典型的偏微分模型方程出發(fā)，采用開源計(jì)算流體力學(xué)庫OpenFoam[6]，建立針對(duì)不同類型PDE的求解器，并進(jìn)行計(jì)算結(jié)果的驗(yàn)證。然后采用PISO算法，將該FVM離散方法應(yīng)用于計(jì)算流體力學(xué)NS方程組，淺水方程組的計(jì)算中，設(shè)計(jì)相應(yīng)的求解器，并對(duì)計(jì)算結(jié)果進(jìn)行驗(yàn)證。

1 有限體積方法（FVM）和PDE的數(shù)值計(jì)算中的應(yīng)用

1.1 PDE方程的分類

偏微分方程的求解方法取決于偏微分方程的類型。因此，不同類型的PDE對(duì)應(yīng)著不同的數(shù)值解法，求解時(shí)所需要的初值條件和邊值條件是不盡相同的。

以一個(gè)2階PDE為例：

引入記號(hào)

1.2 PDE邊界條件

很多物理問題通常都給出以下3種形式的邊界條件的一種：

在進(jìn)行數(shù)值計(jì)算時(shí)，必須根據(jù)特定問題進(jìn)行邊界條件或初始條件的指定，然后才能開始數(shù)值計(jì)算。否則可能無法得到準(zhǔn)確的結(jié)果。

1.3 有限體積法

目前數(shù)值離散方法主要有：有限差分方法，有限元方法和有限體積法。每種方法都會(huì)把一個(gè)簡單的PDE方程轉(zhuǎn)換成相應(yīng)的數(shù)值近似方程。其中有限差分是基于泰勒級(jí)數(shù)來建立方程的數(shù)據(jù)庫，把變量的導(dǎo)數(shù)變?yōu)榭臻g和時(shí)間上不同點(diǎn)之間的差分方法。

在有限元方法中，PDE所使用的計(jì)算域被分割成有限個(gè)子域元。假定在每個(gè)元上因變量按一種規(guī)律變化，而且這種分片拼接的圖像來描述和建立變量在整個(gè)域上的變化情況。

有限體積方法是當(dāng)前最流行的在各種CFD軟件中使用的數(shù)值離散方法。這種方法在某些方面與有限差分相似，但也體現(xiàn)出了從有限元法中吸取的一些特性。

1.4 有限體積法在PDE求解中的應(yīng)用

FVM法使用的規(guī)則網(wǎng)格提高了計(jì)算效率。近年來為了計(jì)算復(fù)雜的幾何形狀，發(fā)展了使用不規(guī)則網(wǎng)格的FVM。有限體積法結(jié)合了FDM和FEM方法的優(yōu)點(diǎn)，因此在PDE求解中有著廣泛的應(yīng)用。如：Takeshi TSUKAMOTO[7]采用FVM離散方法，基于非規(guī)則網(wǎng)格進(jìn)行VOF（volume fraction of fluid）模擬，提出了一套適合自由液面模擬的算法。F. Hermeline等[8]采用有限體積法對(duì)二維麥克斯韋方程進(jìn)行離散求解，并就算法的穩(wěn)定性、復(fù)雜度與其他離散方法進(jìn)行了比較。Zhang Ming Liang[9]提出了一種基于非結(jié)構(gòu)四叉樹網(wǎng)格模型的有限體積法算法，該算法具有健壯的、準(zhǔn)確的并且高效的特性，并成功用該算法對(duì)潰壩問題進(jìn)行模擬。BORIS ANDREIANOV[10]提出了一種收斂的反應(yīng)-擴(kuò)散模型的FVM算法。并對(duì)該FVM算法解存在性和收斂性進(jìn)行了理論分析，并對(duì)解的存在性進(jìn)行了證明。Hrvoje Jasak[6]通過Volume Of Fluid（VOF）方法獲得自由液面進(jìn)行計(jì)算，進(jìn)而采用有限體積法對(duì)船體阻力進(jìn)行了計(jì)算

2 典型PDE的數(shù)值計(jì)算及驗(yàn)證

選取具有代表性的拋物線型PDE進(jìn)行解算和驗(yàn)證，其余類型方程求解器的設(shè)計(jì)方法類似。

算例驗(yàn)證與分析：

對(duì)于熱擴(kuò)散方程

3 FVM方法在計(jì)算流體力學(xué)中的應(yīng)用

為了對(duì)不可壓縮流體求解器進(jìn)行驗(yàn)證，文章選用流體力學(xué)的典型問題方腔驅(qū)動(dòng)問題和潰壩問題進(jìn)行驗(yàn)證。通過采用PISO算法實(shí)現(xiàn)流體力學(xué)運(yùn)動(dòng)控制方程的求解。

3.1 PISO算法介紹

PISO表示有分裂算子的壓力隱式算法，由Issa于1986年提出，它最初作為對(duì)非定?？蓧嚎s流提出的一種非迭代的壓力-速度計(jì)算方法，已經(jīng)成功應(yīng)用于求解定常狀態(tài)的迭代問題。PISO包含1個(gè)預(yù)測(cè)過程和2個(gè)校正過程，可以把它看作是Simple的擴(kuò)展，即對(duì)Simple算法再多一個(gè)修正過程以達(dá)到改善他的目的。PISO算法對(duì)壓力修正方程進(jìn)行2次求解，為了計(jì)算壓力修正方程的源項(xiàng)，需增加存儲(chǔ)量。

3.2 NS方程數(shù)值計(jì)算驗(yàn)證

方腔驅(qū)動(dòng)流的數(shù)值模擬在實(shí)際工程中有著廣泛的應(yīng)用，很多重要的理論和工程研究都是用此模型來分析的，尤其在氣象、航運(yùn)、機(jī)械、采礦等領(lǐng)域應(yīng)用更為廣泛。圖2所示，腔體的長度和寬度都為單位長度，網(wǎng)格劃分為100×100，上邊界為固定速度邊界，其余3面都為固壁邊界。

對(duì)于不可壓縮流體有控制方程：

針對(duì)該控制方程設(shè)計(jì)了專門的求解器，其中非穩(wěn)態(tài)項(xiàng)采用隱式Crank-Nicolson格式，梯度項(xiàng)、散度項(xiàng)采用高斯離散格式。

求解器形式如下式所示：

式中：div為散度，grad為梯度。分別模擬量雷諾數(shù)為400，800，1 000時(shí)的流動(dòng)。如圖3所示為雷諾數(shù)為1 000時(shí)中軸線上速度分布，可以看出，該計(jì)算結(jié)果與K Poochinapan[11]的計(jì)算結(jié)果高度吻合。

3.3 淺水波方程驗(yàn)證

淺水方程是水力學(xué)中的典型問題，有淺水波運(yùn)動(dòng)控制方程：

相應(yīng)的求解器如下式：

式中：hU為速度；div為散度；grad為梯度；h0為水深；h為水位；src為源項(xiàng)。

邊界條件：初始條件和邊界條件是求解水體水流和水質(zhì)方程所必須的條件。一般對(duì)于恒定流的計(jì)算，初始時(shí)刻的變量值可以定義為0，通常稱為冷啟動(dòng)，它不會(huì)影響最終的計(jì)算結(jié)果。對(duì)于非恒定流的計(jì)算，也可以給定冷啟動(dòng)條件，經(jīng)過一段時(shí)間的計(jì)算和調(diào)節(jié)，冷啟動(dòng)帶來的影響可以被消除，計(jì)算也能滿足精度要求[9]。

1）理想潰壩模擬

一般來說，潰壩洪水波的研究可分為數(shù)學(xué)模型、物理模型和兩者結(jié)合等3種類型[9]。由于數(shù)值技術(shù)的優(yōu)越性，操作簡單，僅需給定地形條件、邊界條件就能復(fù)演水波運(yùn)動(dòng)的過程，所以受到越來越多的關(guān)注。同時(shí)也由于強(qiáng)降水、地震等自然災(zāi)害，戰(zhàn)爭(zhēng)等因素導(dǎo)致水位的急劇升高，如果水流越過河堤和壩體，會(huì)導(dǎo)致堤身破壞形成潰口，引起下游河道的沖刷。潰壩所導(dǎo)致的后果極其嚴(yán)重，如1959年法國的Malpasset拱壩潰決，1975年河南板橋水庫土石潰決等，都造成了巨大的生命和財(cái)產(chǎn)損失，因此有必要對(duì)潰壩問題進(jìn)行深入細(xì)致的研究。

理想潰壩問題是典型的非線性流的例子，經(jīng)常用來檢驗(yàn)數(shù)值格式對(duì)強(qiáng)間斷的捕捉能力，并且具有解析解，便于數(shù)值結(jié)果的對(duì)比。

2）二維部分潰壩數(shù)值模擬

目前對(duì)天然河流大多采用平面二維數(shù)學(xué)模型進(jìn)行模擬，并成功地解決了許多工程問題。平面二維數(shù)序模型沿垂向均勻分布考慮，忽略了垂向流速及加速度，水壓力按靜壓分布。

自Fennema和Chaudhry提出二維矩形部分潰壩模型驗(yàn)證數(shù)值算法處理激波捕捉的能力以來，許多學(xué)者研究了此模型，用來驗(yàn)證各種算法。數(shù)值模型的設(shè)置如圖6所示，計(jì)算區(qū)域?yàn)?00×200 m的矩形區(qū)域，中間有一擋水壩將區(qū)域一分為二，壩寬10 m，初始時(shí)刻，上下游水面靜止，上游水深為10 m，下游水深為5 m。某一時(shí)刻，擋水壩突然開有75 m寬的非對(duì)稱缺口，壩體瞬間部分潰決，數(shù)值計(jì)算中，空間步長為5 m，時(shí)間步長0.2 s。圖7所示為7.2 s時(shí)的水面輪廓線。該結(jié)果與其他學(xué)者（圖8）都得到了相近的數(shù)值解從中可以看出，該結(jié)果具有較高的準(zhǔn)確性。

4 結(jié) 語

文章從提升虛擬現(xiàn)實(shí)真實(shí)感的角度出發(fā)，介紹了針對(duì)復(fù)雜PDE數(shù)學(xué)模型的通用離散數(shù)值計(jì)算框架，該框架能夠被重復(fù)利用，并且能夠根據(jù)用戶的需求進(jìn)行定制。而要提升虛擬場(chǎng)景中的物理真實(shí)感，不可避免的要從提升數(shù)值解算效率的角度進(jìn)行研究。文章通過對(duì)典型PDE方程的數(shù)值計(jì)算，模擬了典型的物理模型問題，如拋物線型的熱擴(kuò)散PDE方程，并通過與已有文獻(xiàn)的模擬結(jié)果和實(shí)驗(yàn)結(jié)果進(jìn)行了驗(yàn)證。其次，針對(duì)計(jì)算流體力學(xué)的典型問題，方腔驅(qū)動(dòng)和一維理想潰壩，二維部分理想潰壩問題進(jìn)行的模擬，并對(duì)計(jì)算結(jié)果進(jìn)行了驗(yàn)證，結(jié)果顯示該方法具有較高的可信度。

同時(shí)在后續(xù)的研究中希望能夠?qū)?shù)值計(jì)算的結(jié)果引入到虛擬現(xiàn)實(shí)系統(tǒng)中，將該成果應(yīng)用于虛擬現(xiàn)實(shí)領(lǐng)域系統(tǒng)中，充分發(fā)揮出計(jì)算機(jī)在基于物理數(shù)值計(jì)算方面的價(jià)值。既滿足虛擬現(xiàn)實(shí)系統(tǒng)的實(shí)時(shí)性的要求，又具有較高的準(zhǔn)確性，更加逼真的再現(xiàn)復(fù)雜的物理場(chǎng)景。

[ 1 ]MELL W, JENKINS M A, GOULD J. A physics-based approach to modeling grassland fires[J]. International Journal of Wildland Fire, 2007, 16(1): 1–22.

[ 2 ]柳有權(quán), 劉學(xué)慧, 吳恩華. 基于 GPU 帶有復(fù)雜邊界的三維實(shí)時(shí)流體模擬[J]. Journal of Software, 2006, 17(3): 568–576. LIU You-quan, LIU Xue-hui, WU En-hua. Real-time 3D fluid simulation on GPU with complex obstacles[J]. Journal of Software, 2006. 3(17): p. 568–576.

[ 3 ]王長波, 張卓鵬, 張強(qiáng). 基于LBM的自由表面流體真實(shí)感繪制[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2011, 23(1): 104–110. WANG Chang-bo, ZHANG Zhuo-peng, ZHANG Qiang. Realistic rendering of free surface fluid based on lattice boltzmann method[J]. Journal of Computer-Aided Design & Computer Graphics, 2011. 1(23): p. 104–110.

[ 4 ]王順利, 康鳳舉. 基于物理模型的艦艏浪三維可視化[J]. 計(jì)算機(jī)工程與應(yīng)用, 2014, 2014(3): 18–21. WANG Shun-li, KANG Feng-ju. 3D visualization of bow wave based on physical model[J]. Computer Engineering and Applications, 2014. 3(2014): p. 18–21.

[ 5 ]JIANG Chen fan-fu. The material point method for the Physics-Based simulation of solids and fluid[D]. Los Angeles: UNIVERSITY OF CALIFORNIA Los Angeles, 2015: 1–88.

[ 6 ]HRVOJE J, VUKO V, DOMINIK C. Rapid free surface simulation for Steady-State hull resistance with FVM using openfoam[C]//30th Symposium on Naval Hydrodynamics. Hobart, Tasmania, Australia, 2014: 1–7.

[ 7 ]TSUKAMOTO T, ANZAI K, NIYAMA E. A tracing technique for free surface of fluid flow by FVM method using an irregular grid[J]. Nihon Kikai Gakkai Ronbunshu B Hen/transaction, 1997, 63(6): 77–82.

[ 8 ]HERMELINE F, LAYOUNI S, OMNES P. A finite volume method for the approximation of Maxwell’s equations in two space[J]. Journal of Computational Physics, 2008, 227(1): 9365–9388.

[ 9 ]ZHANG M-l, WM W. A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid[J]. Applied Ocean Research, 2011, 33(4): 297–308.

[10]BORIS A, MOSTAFA B, RICARDO R. ANALYSIS OF A FINITE VOLUME METHOD FOR A CROSS-DIFFUSION MODEL IN POPULATION DYNAMICS[J]. Mathematical Models & Methods in Applied Sciences, 2011, 21(2): 307–344.

[11]POOCHINAPAN K. Numerical implementations for 2D Lid-Driven cavity flow in stream function formulation[J]. ISRN Applied Mathematics, 2012, 2012(4): 1–17.

Application and research of numerical calculation of partial differential equations in virtual reality

ZOU Chang-jun, YIN Yong, LI Hai-jiang, TANG Huang
(Marine Dynamic Simulation and Control Lab, Dalian Maritime University, Dalian 116021, China)

With the development of virtual reality hardware and software, VR technology is entering people’s daily life. However, there is no denying that the core issue of VR-physical reality does not have a good solution. Mostly it is because of the complexity of control equation which takes a lot of time to solve; hence it forces people to adopt alternative method of sacrificing precision to meet the real-time requirements. Whereas many of the physical models of the solution are attributed to solving partial differential equations, so how to accurately and efficiently solve partial differential equation (PDE) and PDEs has a crucial role in enhancing VR-physical sense in the virtual reality system. Based on the solution of typical PDE, the corresponding PDE solver is established and the calculation results are verified. Finally, the method is applied in the field of computational fluid dynamics, the solution of shallow water wave equation of the fluid mechanics in the NS equations and hydraulics, the velocity field in the area or height field are calculated and validated. The results show that the method has high reliability.

virtual reality；PDE；numerical calculation；physical reality

TP391

A

1672 – 7649(2017)08 – 0164 – 06

10.3404/j.issn.1672 – 7649.2017.08.035

2016 – 07 – 13；

2016 – 10 – 18

863課題(2015AA016404)；海洋公益性行業(yè)科研專項(xiàng)(201505017-4)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助(3132016310)

鄒長軍(1987–)，男，博士研究生，從事航海動(dòng)態(tài)仿真、交通系統(tǒng)虛擬現(xiàn)實(shí)技術(shù)等研究。