CN106934192B

CN106934192B - 一种参数优化的浅水方程模型水体建模方法

Info

Publication number: CN106934192B
Application number: CN201511022387.1A
Authority: CN
Inventors: 艾明晶; 徐涵; 杜安定
Original assignee: Beijing University of Aeronautics and Astronautics
Current assignee: Beijing University of Aeronautics and Astronautics
Priority date: 2015-12-30
Filing date: 2015-12-30
Publication date: 2019-11-19
Anticipated expiration: 2035-12-30
Also published as: CN106934192A

Abstract

本发明是一种参数优化的浅水方程模型水体建模方法。本方法首先构建三维动态坐标系统，并推导出参数优化的浅水方程模型。在求解模型时，采用反向路径追踪和特征线法相结合的方法求解方程中的平流项，采用向量投影和傅立叶变换相结合的方法求解外作用力对水体的作用，以此来更新速度场和高度场。接着进行边界检查，对超出边界的网格点采用边界截断法进行处理。最后绘制出高度场，并循环处理下一时间片。本发明与传统的浅水方程相比，当水体与固体接触面过于陡峭时，能更加逼真准确地构建出水体模型，并且在模拟水体建模过程中能够达到平均24.9帧/秒的速率，与传统方法相比速度提高了29.016％，基本上能够实现实时模拟。

Description

一种参数优化的浅水方程模型水体建模方法

技术领域

本发明涉及计算机图形学领域，具体涉及一种参数优化的浅水方程模型水体建模方法。

背景技术

复杂流体模拟技术一直以来都是计算机图形学中一个热门的话题。在计算机流体模拟的众多方法中，基于物理的流体模拟使用流体力学方程来描述水体内部速度、压力、密度等物理量随时间的变化，从而能够真实呈现水体波动以及其他细节。

在描述流体动力学方程方面有著名的纳维-斯托克斯方程(参考文件1：M.Griebel,T.Dornsheifer,T.Neunhoeffer.Numerical Simulation in Fluid Dynamics:A Practical Introduction[M].SIAM:Society for Industrial and AppliedMathematics,1997.)，简称N-S方程。但N-S方程是一个非线性偏微分方程，在求解过程中需要大量的计算节点，导致计算复杂度非常高，因此使用传统的N-S方程来对水体进行建模，从计算机应用的角度出发是不可行的。但是可以通过数值推导来简化N-S方程的复杂度，从而进行相应的求解。从这个角度出发，将水体占据的计算空间按一定规则分为网格，然后分析空间中每一个网格位置上的水体速度、压强、密度等参数随时间和空间的变化，从而衍生出基于欧拉网格的浅水方程方法(Shallow Water Equation)，简称SWE方法。

Hagen等人(参考文件2：T.R.Hagen,J.M.Hjelmervik,K.A.Lieet al.VisualSimulation of Shallow-Water Waves[J].Simulation Modelling Practice andTheory,2005,13(8):716-726.)第一次将SWE引入到计算机图形学中，利用有限元方法在GPU(Graphics Processing Unit，图形处理器)上求解SWE，进而模拟河流、海浪冲刷不规则地形的效果。利用高度场的方法，无法处理破碎的波浪效果，因此，Thury和Müller等人(参考文件3：N.Thurey,M.Muller-Fischer,S.Schirmet al.Real-Time Breaking Waves forShallow Water Simulations[A].Computer Graphics and Applications,2007.PG'07.15th Pacific Conference on[C],2007:39-46)提出了一种能够自动生成细化三角网格补丁的方法来模拟海水冲击岸边后产生的破碎效果，但是该算法不能很好地进行并行计算，从而难以达到实时模拟。Chentanez和Müller(参考文件4：N.Chentanez,M.Müller.Real-Time Simulation of Large Bodies of Water with Small Scale Details[A].Proceedings of the 2010ACM SIGGRAPH/Eurographics Symposium on ComputerAnimation[C],2010:197-206)调研整合了之前的SWE算法，并描述了两种水固交互以及干湿跟踪的方法，这两种方法也在SWE的网格模型的效率和稳定性上进行了进一步的加强，其中还详细介绍了如何使用SWE构建一个复杂水体场景，包括对海浪拍打岛屿岸边、浪花飞溅、泡沫、漩涡等效果，并在其中实现了瀑布效果。但是，上述基于传统浅水方程的水体建模都是建立在高度域的坐标轴方向固定为重力方向的基础上的，因此当与水体交互的固体表面过于陡峭时便会出现失真的现象，无法准确描述水体特征。其次，传统的浅水方程并没有描述外作用力之间是如何共同影响水流变化的，包括水流表面的张力和人为的控制力，同时对于任意的外作用力并没有给出隐式的描述，这样会导致不能逼真地对水体进行建模。

发明内容

本发明的目的是为了解决上述现有传统浅水方程模型水体建模方法的无法逼真模拟水体，以及计算过程中高数值耗散带来的计算速度慢等缺点。本发明提出一种参数优化的浅水方程模型水体建模方法，旨在提高水体建模过程中的计算效率和水体模拟效果的逼真度。

本发明提出的一种参数优化的浅水方程模型水体建模方法，具体包括如下步骤：

步骤一：对N-S方程进行数值推导和条件限制，推导出参数优化的浅水方程模型。

浅水方程中的坐标系包括，高度域空间坐标轴以及平面空间坐标轴。首先，沿着水体与固体接触面的法向量方向建立高度域的空间坐标轴，即h方向；然后，以垂直于该坐标轴的平面建立平面的空间坐标轴，即位置x方向，从而构建浅水方程中使用的坐标系统。参数优化的浅水方程由式(1)和式(2)确定：

其中，u(x,t)是t时刻水体速度在平面x处的水平分量，ρ为水体的密度，p为外作用力(包括空气压力、表面张力以及重力沿h轴的分量)产生的压强沿h轴的分量，a为外作用力在x平面方向产生的水平加速度，h(x,t)为t时刻位置x处的高度域。d＝h-b，其中d为水体的深度，b为地势的高度。u和v为速度u(x,t)在x平面处的分量。符号为空间偏导，在二维曲面中运算符“·”表示向量相乘。式中速度域u和高度域h为待求参数。

上述浅水方程模型的参数优化主要体现在对N-S方程依据动量守恒和质量守恒进行数值推导时的过程简化和形式变换。根据动量守恒定理优化水体速度域偏微分方程时，将式(1)等号右边部分优化为三个组成部分：水体在x平面方向的平流项外作用力在h轴方向对水体的影响以及外作用力在x平面方向对水体的影响。其中，外作用力在h轴方向对水体的影响由项决定，外作用力在x平面方向对水体的影响由x平面方向的水流加速度a决定。

根据质量守恒定理优化水体高度域偏微方程时，将式(2)等号右边部分优化为两个部分：水体在h轴方向的平流项和水流速度对高度域的影响。其中，水流速度对高度域的影响由(u+v)d项决定。在水体运动过程中，物理量如速度u和高度h是随着时间和空间变化的，即速度和高度是时间t和空间x的函数，在式(1)和式(2)中分别表现为和统称为平流项。平流项的存在，使得浅水方程是非线性的，因此在求解过程中比较复杂。本发明对传统浅水方程进行上述参数优化有利于建立求解模型，提高方程的求解效率。

步骤二：为步骤一构建的参数优化的浅水方程模型建立求解模型。

步骤2.1，离散化步骤一中的待解参数。设浅水方程的空间计算网格为△τ，网格大小为△x×△x，时间步长为△t，在网格点(i,j)处，高度域为h_i,j，水平速度分量分别为u_i,j和v_i,j，作如下形式离散化：

步骤2.2，初始化浅水方程中的各项参数，包括初始高度场h，速度场u＝(u,v)以及地形高度b等。本发明中各参数的初始化如下：u＝0，v＝0，h＝2，b＝1。这样，水体的初始深度为：d＝h-b＝1。实际操作过程中可以根据需要修改各项参数的初始值。

步骤2.3，计算浅水方程中的平流项。采用反向路径追踪和特征线法来计算平流项和反向路径跟踪的思想是，要计算t时刻位置x处的粒子在新的时间间隔△t后的速度，首先定义一条由t时刻到t+△t时刻的路径p(x,s)，其中s为时间间隔[t,t+△t]中的一个时间点，新的速度项w_new由△t时间前旧的速度项w_old通过计算得来：

w_new(x)＝w_old(p(x,-△t)) (6)

特征线法是一种解偏微分方程比较成熟的方法，其基本思想是通过把偏微分方程转化为两组常微分方程，再对常微分方程进行求解。两组常微分方程中的一组用于定义特征线，另一组用以描述解沿给定特征线的变化。这样，结合反向路径跟踪和特征线法便可计算出平流项，详见具体实施方式中步骤2.3。

步骤2.4，计算方程中外作用力对流体的作用，更新流体的速度场和高度场。设初始速度项为w₀，计算外作用力对流体作用后的速度项，首先计算外作用力的直接作用效果：

w₁＝w₀(x)+△tF(x,t) (7)

其中，F(x,t)为外作用力在平面x处的分量，由加速度项a决定；w₁为外作用力作用于水体后水体的速度。

然后对求得的w₁进行傅立叶变换：

对变换后的项进行投影：

对投影后的项进行傅立叶反变换，得到最后的结果：

上述计算所得结果w₂即为新的速度域，将该新的速度域代入式(2)中并结合式(3)可以求出新的高度域。最后，用新的速度域和高度域更新旧的速度域和高度域，准备下一轮的迭代求解。

步骤2.5，进行边界条件检查，越过阈值处采用边界截断法。在计算平流项时，反向路径追踪到的网格点可能超过定义的网格阈值范围产生越界。在处理这种情况时，当粒子距离边界的距离小于半个网格，甚至超过了网格，设置其位置为距离边界半个网格距离的位置。并根据上述步骤中计算出的速度域和高度域的值，调用绘图程序绘制在当前时间t和空间位置x处的水体。在模拟时间内跳转到步骤2.3，处理下一时间片。

步骤三：利用步骤二所建立的参数优化的浅水方程模型构建水体场景，输出结果并显示。步骤二中建立的水体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列网格和粒子构成，并不能逼真地呈现现实中的水体场景，因此还需要对上述的模型进行相应的渲染操作，才能真实呈现水体模型。本发明中采用OpenGL提供的库函数进行编程，对水体进行渲染，从而展现逼真的水体效果。

本发明方法的优点和积极效果在于：本发明在建立参数优化的浅水方程模型过程中，采用动态变化的坐标系统，即h坐标的方向不是固定为重力方向，而是根据水体与固体接触面的法向量来决定，提高了模型的稳定性，减轻了由于接触面过于陡峭带来的破碎效果；然后在求解方程的过程中采用向量投影与傅立叶变换相结合的方式计算外作用力的影响，以及采用反向路径跟踪和特征线法相结合的方式计算平流项，提高了方程的求解速率，同时更加逼真地模拟了水体场景。

附图说明

图1是本发明参数优化的浅水方程模型水体建模方法的整体步骤流程图；

图2是本发明三维动态坐标系统图；

图3是本发明物理模型示意图；

图4是本发明反向路径跟踪示意图。

具体实施方式

下面将结合附图和实施例对本发明的技术方案作进一步的详细说明。

本发明提出一种参数优化的浅水方程模型水体建模方法，通过在建立浅水方程模型过程中，采用动态坐标系统建模，提高了模型的稳定性；然后在求解方程的过程中，对于平流项的计算采用反向路径跟踪和特征线法相结合的方式，对于外作用力的计算采用向量投影与傅立叶变换相结合的方式，提高了模型的计算效率和性能。

本发明提出一种参数优化的浅水方程模型水体建模方法，旨在提高传统浅水方程模型对水体进行建模的性能，下面具体对实现步骤进行说明。

步骤1.1，沿着水体与固体表面接触面的法向量方向建立一维高度域坐标系h，再以该高度域坐标系为基础建立二维平面坐标系x，从而构建本模型的三维动态变化的坐标系统。

步骤1.2，对水体建模过程中的所有场元素，包括高度场、速度场以及外作用力等进行物理建模。图2三维动态坐标系统图中展示了各个场元素的物理模型。为了避免水体之间的相互碰撞，首先，在一个低分辨率接触面处计算其平均曲面法线的值；然后，定义地势b(x)为接触面原始曲面法线与平均曲面法线的差值。这样根据水体的高度域h(x,t)函数就可以计算出水体深度为d＝h-b。定义u(x,t)为水平速度(这里的水平并不是传统意义上的水平面，而是平面x方向)，垂直速度由隐式给出。最后对外作用力进行物理建模，外作用力包括重力、大气压力、表面张力以及人为控制力等。将外作用力分解为h轴方向的一维作用力p和x平面方向的二维作用力F(x,t)，其所产生的水平加速度为a。

步骤1.3，根据步骤1.2中所建立的物理模型，将三维的N-S方程限制到二维曲面，来构建本发明的参数优化的浅水方程，如上述所构建的式(1)和式(2)。其中式(1)由动量守恒推导而来，外作用力对水体的作用以加速度的形式体现，通过式(1)来更新水体的速度u；式(2)由质量守恒推导而来，通过式(2)来更新水体的高度域h。

本发明步骤一中的三维动态坐标系统与传统的浅水方程坐标系统相比，可以避免因水体与固体接触面过于陡峭带来的水体破碎等不稳定的弊端，具有更好的稳定性；同时，对水体的物理建模以及从三维N-S方程推导至二维的浅水方程，与传统浅水方程相比具有更高的求解效率。

步骤2.1，离散化步骤一中建立的物理模型，转化为易于求解的数学模型。设该物理模型所使用的空间网格为△τ，每个网格的大小为△x×△x，每一个时间片的大小为△t，在网格在网格点(i,j)处，高度域为h_i,j，水平速度u_i,j分量分别为u_i,j和v_i,j。在空间上进行离散，将连续变化的流体离散为由一个个小网格构成，对任意场的空间微分，只需求得与相邻值的中央差分即可。在时间上进行离散，将连续的的时间离散为由一个个时间片构成，只需在每个时间片上计算流体运动即可。因此，离散结果由式(3)～(5)给出。

步骤2.2，初始化方程中的各项参数，包括初始高度场，速度场以及地形高度等。本发明中各参数的初始化如下：u＝0，v＝0，h＝2，b＝1。这样，水体的初始深度为：d＝h-b＝1。实际操作过程中可以根据需要修改各项参数的初始值。

步骤2.3，计算方程中的平流项。采用反向路径追踪和特征线法来计算平流项和图3反向路径跟踪示意图表明了这个过程。特征线法可以求解类似如下形式的平流问题：

其中为标量，初始条件为定义向量u的特征线n(x₀,t)是由位置x₀处的粒子从时间t＝0开始平移而来，则有：

其中，初始条件为n(x₀,0)＝x₀。

然后，定义矢量是由特征线从初始状态x₀，t＝0平移而来。这样可以通过下面的微分方程进行计算：

其中特别的有所以在计算t时刻位置x处的值时，首先沿着特征线反向追踪位置x处先前的位置，直到位置x₀为止，然后初始化该位置的值：

根据上述特征线法可以计算本模型中的平流项，如上述式(6)。

步骤2.4，计算方程中外作用力对流体的作用，更新流体的速度场和高度场。该步骤主要针对式(2)动量守恒形式。

著名的Helmholtz-Hodege分解可以将一个任意的向量η分解成如下形式：

其中，向量ξ具有零散度，即a是标量。这个分解表明任何一个向量可以由一个质量守恒的域和一个标量域的和组成。因此，可以定义一个算子K将任意向量η投影到一个零散度的向量ξ＝Kη，于是有：

将上式代入到式(2)的等式两边可得：

其中Ku＝u，并且根据该分解可以建立本发明中的投影模型。

将本模型中的边界条件限制在周期边界中，利用如下傅立叶变换进行相应的计算：

其中k＝|k|。

因此，在计算外作用力对水体的作用以及更新速度场和高度域时，首先计算力的直接作用效果如式(7)，然后对求得的项进行傅立叶变换如式(8)，接着采用Helmholtz-Hodege分解投影如式(9)，最后对投影后的项进行傅立叶反变换求得最终结果，并进行相应的域更新。

步骤二建立的计算模型通过将基于Helmholtz-Hodege分解的向量投影和傅立叶变换相结合来求解外作用力对水体的作用，同时将反向路径跟踪和特征线法相结合来求解平流项，能够提高浅水方程的求解速度，从而提高水体建模的效率。

表1展示了传统的浅水方程模型水体建模方法与本发明提出的参数优化的浅水方程模型水体建模方法的模拟速率对比表。可以看到，本发明提出的方法平均模拟速率为24.9帧/秒，相对于传统浅水方程方法平均模拟速率提高了29.016％。

表1模拟速率比较

本发明所有测试都是在联想Y460笔记本电脑上进行的，CPU型号是英特尔Corei3M 350@2.27GHz，显示卡型号是ATI Mobility Radeon HD 5650(1024MB)。平均模拟速率是通过记录多次模拟实验的结果，求平均值得来的，不同硬件环境所得结果可能有所不同。

Claims

1.一种参数优化的浅水方程模型水体建模方法，其特征在于，实现步骤如下：

步骤一：对纳维-斯托克斯方程(N-S方程)进行数值推导和条件限制，推导出参数优化的浅水方程模型；

首先，构建参数优化的浅水方程的三维动态坐标系统；其次，对水体建模过程中所有涉及到的域进行物理建模；然后，将三维N-S限制到二维曲面，结合数学推导简化来构建参数优化的浅水方程：

其中，u(x,t)是t时刻速度在平面x处的水平分量，ρ为水体的密度，p为外作用力产生的压强沿h轴的分量，a为外作用力在x平面方向产生的水平加速度，h(x,t)为t时刻位置x处的高度域；d＝h-b，其中d为水体的深度，b为地势的高度；u和v为速度u(x,t)在x平面处的分量；符号为空间偏导，在二维曲面中运算符“·”表示向量相乘，式中速度域u和高度域h为待求参数；

步骤二：为步骤一构建的参数优化的浅水方程模型建立求解模型；

步骤2.1，离散化步骤一中的待解参数，设方程的空间计算网格为Δτ，网格大小为Δx×Δx，时间步长为Δt，在网格点(i,j)处，高度域为h_i,j，水平速度分量分别为u_i,j和v_i,j，作如下形式离散化：

步骤2.2，初始化方程中的各项参数，包括初始高度场，速度场以及地形高度，各参数的初始化如下：u＝0，v＝0，h＝2，b＝1，这样，水体的初始深度为：d＝h-b＝1，实际操作过程中可以根据需要修改各项参数的初始值；

步骤2.3，计算方程中的平流项，采用反向路径追踪和特征线法来计算平流项和具体是，要计算t时刻位置x处的粒子在新的时间间隔Δt后的速度，首先定义一条由t时刻到t+Δt时刻的一条路径p(x,s)，新的速度项由Δt时间前的速度项来决定w_new(x)＝w_old(p(x,-Δt))；

步骤2.4，计算方程中外作用力对水体的作用，更新流体的速度场和高度场，采用向量投影和傅立叶变换相结合的方法来计算外作用力对水体的作用，设初始速度项为w₀，计算外作用力对水体作用后的速度项，首先计算外作用力的直接作用效果w₁＝w₀(x)+ΔtF(x,t)，其中F(x,t)为外作用力平面x处的分量，由公式中的加速度项a决定；然后对求得的项进行傅立叶变换接着对变换后的项进行投影最后对投影后的项进行傅立叶反变换，得到最后的结果计算所得w₂即为新的速度域，将该新的速度域代入参数优化的浅水方程中可以求出新的高度域；最后，用新的速度域和高度域更新旧的速度域和高度域，准备下一轮的迭代求解；

步骤2.5，进行边界条件检查，越过阈值处采用边界截断法，在计算平流项时，反向路径追踪到的网格点可能超过定义的网格阈值范围产生越界，在处理这种情况时，当粒子距离边界的距离小于半个网格，或者超过了网格，设置其位置为距离边界半个网格距离的位置，并根据上述步骤中计算出的速度域和高度域的值，调用绘图程序绘制在当前时间t和空间位置x处的水体，在模拟时间内跳转到步骤2.3，处理下一时间片；

步骤三：利用步骤二所建立的参数优化的浅水方程模型构建水体场景，输出结果并显示，步骤二中建立的水体模型，仅仅是基于数学物理求解得出的结果，从视觉上描述，是由一系列网格和粒子构成，并不能逼真地呈现现实中的水体场景，因此还需要对上述的模型进行相应的渲染操作，才能真实呈现水体模型，采用OpenGL提供的库函数进行编程，对水体进行渲染，从而展现逼真的水体效果。

2.根据权利要求1所述的一种参数优化的浅水方程模型水体建模方法，其特征在于，步骤一中所建立的参数优化的浅水方程，具体是：以水体与固体接触面的曲面法向量方向为高度域的方向，动态建立三维坐标系统，然后对各个场元素进行相应的物理建模，最后通过对N-S方程进行相应的推导，建立最终的参数优化的浅水方程模型。