偏微分方程數值計算在虛擬現實中應用與研究

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

(大連海事大學 航海動態仿真和控制實驗室，遼寧 大連 116021)

偏微分方程數值計算在虛擬現實中應用與研究

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

(大連海事大學 航海動態仿真和控制實驗室，遼寧 大連 116021)

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

虛擬現實；PDE；數值計算；物理真實感

0 引 言

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

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

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

從已有的研究中可看出，虛擬現實技術的物理真實感一直都是一項具有挑戰性的課題。

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

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

1 有限體積方法（FVM）和PDE的數值計算中的應用

1.1 PDE方程的分類

偏微分方程的求解方法取決于偏微分方程的類型。因此，不同類型的PDE對應著不同的數值解法，求解時所需要的初值條件和邊值條件是不盡相同的。

以一個2階PDE為例：

引入記號

1.2 PDE邊界條件

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

在進行數值計算時，必須根據特定問題進行邊界條件或初始條件的指定，然后才能開始數值計算。否則可能無法得到準確的結果。

1.3 有限體積法

目前數值離散方法主要有：有限差分方法，有限元方法和有限體積法。每種方法都會把一個簡單的PDE方程轉換成相應的數值近似方程。其中有限差分是基于泰勒級數來建立方程的數據庫，把變量的導數變為空間和時間上不同點之間的差分方法。

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

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

1.4 有限體積法在PDE求解中的應用

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

2 典型PDE的數值計算及驗證

選取具有代表性的拋物線型PDE進行解算和驗證，其余類型方程求解器的設計方法類似。

算例驗證與分析：

對于熱擴散方程

3 FVM方法在計算流體力學中的應用

為了對不可壓縮流體求解器進行驗證，文章選用流體力學的典型問題方腔驅動問題和潰壩問題進行驗證。通過采用PISO算法實現流體力學運動控制方程的求解。

3.1 PISO算法介紹

PISO表示有分裂算子的壓力隱式算法，由Issa于1986年提出，它最初作為對非定常可壓縮流提出的一種非迭代的壓力-速度計算方法，已經成功應用于求解定常狀態的迭代問題。PISO包含1個預測過程和2個校正過程，可以把它看作是Simple的擴展，即對Simple算法再多一個修正過程以達到改善他的目的。PISO算法對壓力修正方程進行2次求解，為了計算壓力修正方程的源項，需增加存儲量。

3.2 NS方程數值計算驗證

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

對于不可壓縮流體有控制方程：

針對該控制方程設計了專門的求解器，其中非穩態項采用隱式Crank-Nicolson格式，梯度項、散度項采用高斯離散格式。

求解器形式如下式所示：

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

3.3 淺水波方程驗證

淺水方程是水力學中的典型問題，有淺水波運動控制方程：

相應的求解器如下式：

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

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

1）理想潰壩模擬

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

理想潰壩問題是典型的非線性流的例子，經常用來檢驗數值格式對強間斷的捕捉能力，并且具有解析解，便于數值結果的對比。

2）二維部分潰壩數值模擬

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

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

4 結 語

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

同時在后續的研究中希望能夠將數值計算的結果引入到虛擬現實系統中，將該成果應用于虛擬現實領域系統中，充分發揮出計算機在基于物理數值計算方面的價值。既滿足虛擬現實系統的實時性的要求，又具有較高的準確性，更加逼真的再現復雜的物理場景。

[ 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 ]柳有權, 劉學慧, 吳恩華. 基于 GPU 帶有復雜邊界的三維實時流體模擬[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 ]王長波, 張卓鵬, 張強. 基于LBM的自由表面流體真實感繪制[J]. 計算機輔助設計與圖形學學報, 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]. 計算機工程與應用, 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)；海洋公益性行業科研專項(201505017-4)；中央高校基本科研業務費專項資金資助(3132016310)

鄒長軍(1987–)，男，博士研究生，從事航海動態仿真、交通系統虛擬現實技術等研究。