﻿ 偏微分方程数值计算在虚拟现实中应用与研究
 舰船科学技术  2017, Vol. 39 Issue (8): 164-169 PDF

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
Abstract: 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.
Key words: virtual reality     PDE     numerical calculation     physical reality
0 引　言

1 有限体积方法（FVM）和PDE的数值计算中的应用 1.1 PDE方程的分类

 $\begin{split}{a_{11}}\frac{{{\partial ^2}\varPhi }}{{\partial {x^2}}} + & 2{a_{12}}\frac{{{\partial ^2}\varPhi }}{{\partial x\partial y}} + {a_{22}}\frac{{{\partial ^2}\varPhi }}{{\partial {y^2}}} + \\& {a_1}\frac{{\partial \varPhi }}{{\partial x}} + {b_1}\frac{{\partial \varPhi }}{{\partial y}} + c\varPhi + f = 0\text{，}\end{split}$ (1)

 $\varDelta = - \left| {\begin{array}{*{20}{c}}{{a_{11}}}&{{a_{12}}}\\{{a_{21}}}&{{a_{22}}}\end{array}} \right| = a_{12}^2 - {a_{11}}{a_{12}}\text{，}$ (2)

1）若在点 $(x,y)$ 处， $\varDelta > 0$ ，成方程在点 $(x,y)$ 为双曲线型。

2）若在点 $(x,y)$ 处， $\varDelta = 0$ ，成方程在点 $(x,y)$ 为抛物线型。

3）若在点 $(x,y)$ 处， $\varDelta < 0$ ，成方程在点 $(x,y)$ 为椭圆型。

1.2 PDE边界条件

1）Dirichlet边界条件也称为第一类边界条件问题。在这种边界条件下，已知边界条件上变量的值是常数。这就允许进行简单代替来给边界定值。比如当给定边界的为恒定入流条件，这时就可以设置入流边界的速度为常数U，即 $u{\left. \begin{array}{l}\\\end{array} \right|_\varGamma } = g$

2）Neumann边界条件也称为第二类边界条件问题。这种边界条件下，边界给定的导数值是已知的，这时就给出了一个附加的方程，可以用这个方程得到边界上的值。如当给定的条件为出口压力恒定条件，这是压力的导数为0，即 $\alpha \displaystyle\frac{{\partial u}}{{\partial n}}{\left. \begin{array}{l}\\\end{array} \right|_\varGamma } = g$ ，从而推出边界出口前后网格节点的压力值相等。

3）Bobin边界条件也称为第三类边界问题也称为混合边界值问题。在边界 $\varGamma$ 上给定函数值及其方向偏导数值的组合，形式如式（3）所示：

 $(\alpha \frac{{\partial u}}{{\partial n}} + \beta u){\left. \begin{array}{l}\\\end{array} \right|_\varGamma } = g\text{。}$ (3)

1.3 有限体积法

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的数值计算及验证

 $\frac{{\partial T}}{{\partial t}} = \alpha (\frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{{\partial ^2}T}}{{\partial {y^2}}})\text{，}$ (4)

 ${\rm ddt}(T) = \alpha *{\rm laplacian}(T)\text{。}$ (5)

 图 1 数值计算结果与解析解对比 Fig. 1 Comparison between numerical results and analytical solutions
3 FVM方法在计算流体力学中的应用

3.1 PISO算法介绍

PISO表示有分裂算子的压力隐式算法，由Issa于1986年提出，它最初作为对非定常可压缩流提出的一种非迭代的压力-速度计算方法，已经成功应用于求解定常状态的迭代问题。PISO包含1个预测过程和2个校正过程，可以把它看作是Simple的扩展，即对Simple算法再多一个修正过程以达到改善他的目的。PISO算法对压力修正方程进行2次求解，为了计算压力修正方程的源项，需增加存储量。

3.2 NS方程数值计算验证

 图 2 方腔驱动模型设置 Fig. 2 Cavity drive model setup

 \left. \begin{aligned}\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho V) = 0\\\frac{{\partial \rho u}}{{\partial t}} \!+\! \nabla \cdot (\rho uV)\! =\! \! -\! \frac{{\partial p}}{{\partial x}} + \frac{{\partial {\tau _{xx}}}}{{\partial x}} \!+\! \frac{{\partial {\tau _{yx}}}}{{\partial y}} \!+\! \frac{{\partial {\tau _{zx}}}}{{\partial z}} \!+\! \rho {f_x}\\\frac{{\partial \rho v}}{{\partial t}} \!+\! \nabla \cdot (\rho vV)\! =\! \! -\! \frac{{\partial p}}{{\partial y}} \!+\! \frac{{\partial {\tau _{xy}}}}{{\partial x}} \!+\! \frac{{\partial {\tau _{yy}}}}{{\partial y}} \!+\! \frac{{\partial {\tau _{xy}}}}{{\partial z}} \!+\! \rho {f_y}\\\frac{{\partial \rho w}}{{\partial t}} \!+\! \nabla \cdot (\rho wV)\! =\! \! -\! \frac{{\partial p}}{{\partial z}} \!+\! \frac{{\partial {\tau _{xz}}}}{{\partial x}} \!+\! \frac{{\partial {\tau _{yz}}}}{{\partial y}} \!+\! \frac{{\partial {\tau _{zx}}}}{{\partial z}} \!+\! \rho {f_z}\end{aligned} \right\}\text{。} (6)

 ${\rm ddt}(U) + div(phi,U) - {\rm laplacian}(nu,U) = - grad(p)\text{。}$ (7)

 图 3 雷诺数为1 000计算结果对比 Fig. 3 Comparison of calculation results at reynolds numberr 1 000
3.3 浅水波方程验证

 \left. \begin{aligned}\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}} + w\frac{{\partial u}}{{\partial z}} - fv = - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} + \frac{1}{\rho }\frac{{\partial {\tau _x}}}{{\partial z}}\\\frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}} + w\frac{{\partial v}}{{\partial z}} + fu = - \frac{1}{\rho }\frac{{\partial p}}{{\partial y}} + \frac{1}{\rho }\frac{{\partial {\tau _y}}}{{\partial z}}\\0 = - \frac{1}{\rho }\frac{{\partial p}}{{\partial z}} - g\end{aligned} \right\}\text{。} (8)

 ${\rm ddt}(hU) + div(phi,hU) = - grad(h + {h_0}) + src(f)\text{。}$ (9)

1）理想溃坝模拟

 图 4 文献[11]计算结果和解析解对比 Fig. 4 Results comparison from paper [11] and analytical solutions

 图 5 不同数值模式结果 Fig. 5 Results from different numerical scheme

2）二维部分溃坝数值模拟

 图 6 部分溃坝模型设置 Fig. 6 Partial dam-break model setup

 图 7 溃坝数值计算结果 Fig. 7 Numerical results of dam-break

 图 8 文献[9]计算结果 Fig. 8 Result from pater [9]
4 结　语

 [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. DOI: 10.1071/WF06002 [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. DOI: 10.1016/j.apor.2011.07.004 [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.