基于物理模型的海浪真实感模拟研究

朱天惠，任鸿翔，赵玉然

(大连海事大学航海动态仿真和控制交通行业重点实验室，辽宁 大连116026)

1 引言

航海模拟器的高速发展为航海院校及科研院所培养输送大批高素质海员提供了坚实的基础。航海模拟器的优劣主要体现在功能或称为适任能力、物理真实感、行为真实感和操作环境真实感等几个方面。视景系统是近年来航海模拟器发展最快的部分之一，是模拟器研制成功与否的关键技术，其主要体现在提升模拟器的操作环境真实感[1]。

航海模拟器中海浪约占整个视景的一半区域，其模拟的逼真与否对环境真实感影响较大。现阶段航海模拟器中有关海浪的模拟主要采用的是基于构造的方法，该方法主要包括基于几何模型的方法和基于统计和谱的方法。基于构造的方法的优点是实时性较强，缺点是真实感较差，难以展现流体细节，而且在较高视点处观察海浪，磁砖效应比较明显[2]。

随着计算机硬件设备的高速发展，图形处理器的计算速度大大加快，基于物理模型的方法逐渐步入视野。基于物理模型的海浪绘制方法现阶段主要分为欧拉网格法(Grid-basedEulerian Approach)以及拉格朗日粒子法(Particle-based Lagrangian Approach)欧拉法通过固定网格单元来模拟流体在单元上的移动，其缺点是易受时间步长条件限制，且当网格分辨率增加，真实感较强时，限制效果会更明显。拉格朗日法则通过分析流体质点在任意时刻的速度、运动轨迹、密度、压力等物理量，计算机图形学中使用粒子来表示质点，通过粒子来离散控制方程。拉格朗日法的优势在于模拟真实感较强，且可以在三维空间模拟例如水花、泡沫等更加丰富的流体细节、处理更加复杂的流体表面[3]。常用的拉格朗日方法有SPH(SmoothedParticleHydrodynamics)方法以及PBF(PositionBasedFluids)方法等。因为PBF方法相比较传统SPH方法，具有允许大时间步长，实时性较好等特点，基于位置的流体绘制方法成为海浪仿真研究的热点。Müller等[5]最先提出PBD(PositionBasedDynamics)框架，并用来模拟布料。Macklin[5]等首次将PBD框架用于流体模拟，提出一个融合了密度迭代求解器的PBF框架，实现了与当前SPH求解器接近的不可压缩性和收敛性，稳定性强，时间步长大，适合实时渲染[3]。Kang等[6]通过构建无散速度场，加快收敛速度。KÖster等[7]提出一种自适应迭代的APBF(Adaptive PBF)方法，适用于更大场景流体实时模拟。

结合上述分析，综合考虑物理海浪流体建模方法，提出一种基于位置动力学的海浪真实感模拟方法。主要通过构建PBF流体模型，利用涡度约束函数进行Curl运算求解粒子的旋度场，再根据所求得的旋度场与位置向量进行计算，得到涡度力场，最后根据力场所得旋向加速度实时更新流体粒子的速度，从而实现实时的真实感海浪模拟。

2 流体模型构建

PBF方法是一种基于PBD框架下的无网格、拉格朗日粒子法。PBD框架的本质是构建约束条件，然后通过对约束条件进行投影，更新流体粒子的位置信息，并借此更新速度值。PBF作为基于PBD框架的流体模拟方法，通过构建密度约束，迭代求解流体粒子的位移偏量，根据位移偏量，实时更新流体粒子的速度，再将流体粒子绘制出来。

2.1 基于位置的动力学框架

传统方法通过求解力的值，再根据力的值进行数值积分，从而更新速度值。与传统的基于力的求解方法不同，基于位置的动力学框架采用的是约束的手段。在PBD方法中，当检测到两个物体发生穿透时，直接根据所采用的约束条件修正物体的位置，再更新物体的速度信息。PBD框架[4]算法如表1所示：

表1 基于位置的动力学框架算法

基于位置的动力学框架核心在于第7)～15)行。第7)行通过显式欧拉积分方法估算物体顶点的新位置，9)～11)行是一个迭代解算器，通过采用高斯-赛德尔迭代的方法反复进行约束投影来控制估算位置从而保证位移满足约束条件。通过公式可以表示为：

(1)

在13)、14)行，物体顶点移动到了最佳估算位置，速度也随之更新，速度的控制通过第5)、6)及16)行实现。无法转换为位置约束的外力通过第5)行进行速度更新，系统只考虑重力对物体的影响，所以合外力与质量倒数的乘积等于重力加速度。第6)行用于特定情况下添加速度阻尼，第16)行根据摩擦系数和补偿系数调整碰撞顶点的速度。

2.2 PBF流体模型

通过与传统SPH插值方法结合，Macklin等人[5]提出一种基于PBD框架下的流体模拟方法PBF。不同于传统PBD框架用于解决固体问题而采用的位置约束，PBF方法通过计算流体质点的密度，采用一种密度约束方法来更新流体质点的位置。PBF粒子系统建模流程如图1所示。

图1 PBF粒子系统建模流程

为保证流体的不可压缩性，针对每一个流体粒子，都需要满足密度约束条件

(2)

式中，ρ0是静态密度，ρi通过标准SPH密度估算器给出

(3)

在流体模型中，认为所有的流体粒子拥有均等的质量，所以后续公式中舍掉质量项mj。式中W为光滑核函数，根据Muller[8]首次将SPH方法用于模拟流体时采用的核函数，本文选用Poly6核函数用于密度估算，选用Spiky核函数用于梯度计算，核函数如下

(4)

(5)

Monaghan等人。给出了在粒子层面上的梯度公式，将其用于密度约束，可以得到对于任意流体粒子k的约束条件的梯度为

(6)

梯度公式随着流体粒子的类型不同而产生变化，流体粒子k可分为对象粒子以及邻域粒子

(7)

将其带入到式(1)，可以求得式中拉格朗日算子λ的值为

(8)

每个流体粒子都具有拉格朗日算子λ，算子的值对于同一流体模型中的所有流体粒子都是相同的。

因为式(2)中的约束条件是非线性的，且在光滑核函数的边界梯度逐渐消失。这会导致所求得的拉格朗日算子在分母无穷小时的不稳定。根据Smith提出的混合约束力理论，可以在密度约束阶段添加一个大小为ξλ的约束力，从而使约束更加合理，ξ是人为定义的松弛参数常量，考虑影响的拉格朗日算子为

(9)

由此，可以得到考虑邻域粒子密度约束修正情况下，且保证流体不可压缩性条件的流体粒子的位移值为

(10)

3 旋涡限制力

为实现海浪表面的逼真模拟，并体现海浪粒子运动实时性的特点，需要为海浪表面粒子添加旋涡限制力。本文在添加旋涡限制力阶段采用的方法是涡度约束的方法[11]。不同于文献[4]利用涡度约束函数来抵消因为方法本身产生的附加阻尼导致的原有涡旋的耗散，本文通过讨论选用合理化调节系数，通过涡度约束函数产生旋涡控制力调节流体粒子系统速度场，从而达到逼真模拟海浪的目地。

图2 一个简单向量的旋度场

3.1 Curl运算

在进行流体粒子涡度计算的时候，采用了Curl运算。向量分析中，旋度(Curl)是一个向量算子，表示在三维欧几里得空间中的向量场的无穷小量旋转，在向量场每个点上，点的旋度表示为一个向量，称为旋度向量。这个向量的特征刻画了在这个点上的旋转。

在不同的坐标系下，向量场的旋度有不同的表达方式。PBF流体模型中采用的是流体粒子三维直角坐标系，所以三维直角坐标系中展开Curl运算。

设向量场A为

A(x，y，z)=Ax(x，y，z)i+Ay(x，y，z)j+Az(x，y，z)k

(11)

(12)

采用向量积的行列式记号形式可表示为：

(13)

3.2 涡度约束

涡度约束[12]的方法需要首先对于每一个流体粒子i采用Curl运算计算其涡度ωi：

(14)

式中，mj为流体粒子的质量，ρj为流体粒子的密度，vi流体粒子的速度，vj为粒子i的邻域粒子的速度。采用粒子i与邻域粒子的速度差与梯度核函数进行叉乘，得到Curl运算的结果。

得到流体粒子的涡度ωi后，需计算流体粒子的位置矢量N：

(15)

(16)

将计算所得位置矢量与已得涡度进行叉乘，得到流体粒子的旋涡控制力F：

(17)

最后，根据由旋涡控制力所求得的加速度，在时间步长内更新流体粒子的速度。在调节系数ε的控制下，旋涡控制力大小的不一致产生不同程度的流体粒子表面的起伏，从而实现逼真地模拟海浪。

4 结果与讨论

4.1 实验结果

实验平台是Intel(R)Core(TM)i5-7300HQ CPU 2.50 GHZ，主存16GB，GPU为NVIDIAGeForceGTX 1050 Ti，显卡内存为4GB，开发环境VisualStudio 2015，采用OpenGL作为图形的编程接口。

图3显示了粒子系统的初始形态，流体粒子数为12000，调节系数选定为0，即不添加旋涡控制力时流体粒子的溃堤模型。从图中可以看出，粒子模型表面平整，无涡旋痕迹，较为光滑，是一种静水中水面的模拟。但切合航海实际，大多数情况下海面不会是平静的，在风流、浪涌的作用下，海面会产生浪花等不规则运动。所以添加旋涡限制力来模拟海面是必要地。

图3 无添加旋涡限制力的流体粒子模型

图4显示了当调节系数选定为1.0时的流体粒子模型。此时随着旋涡控制力的加入，流体模型表层粒子产生涡旋效果，涡旋与涡旋的边界相叠加从而产生类似于波峰波谷的变化。

图4 调节系数1.0下的流体粒子模型

随着调节系数的不断加大，流体粒子的运动剧烈程度逐渐加大。图5为调节系数2.0下的流体粒子模型。与调节系数1.0下的粒子模型不同的是，当系数加大至2.0时，可以清楚地看到在模型的上方有部分粒子飞溅，同时整个模型的运动程度相比较系数1.0下更为剧烈。从航海实际来看，剧烈的海浪运动常常带有翻卷、破碎飞溅等细节特征。当调节系数选取合适时，飞溅的粒子能恰好的模拟出海浪运动产生浪花飞溅等特效。

图5 调节系数2.0下的流体粒子模型

4.2 比较分析

从实验结果可以看出，流体粒子系统的运动剧烈程度随着调节系数的增大而变得更加剧烈，但是海浪模拟的真实感不会随着系数的增大而变得更为真实。当系数增大到一定程度时，流体粒子的旋向加速度会使其旋向速度过大，从而导致流体粒子发生横向漂移。此时，流体粒子的运动已不符合自然物理规律，无法真实地模拟海浪运动。

表2给出了在不同调节系数下海浪运动的模拟效果。当调节系数选取在0～0.9之间时，流体粒子模型的涡旋程度只能模拟出原有流体粒子涡旋的耗散；当系数选取在1.0～2.9之间时，伴随着所选取系数的增大，海面的波动程度不断加深；系数选取在3.0～3.9之间时，整个流体粒子系统演变为一个大涡旋的旋转运动，此时已无法真实地模拟海浪运动；当系数的选取大于4.0时，流体粒子系统发散。

表2 不同调节系数下的模拟效果

通过综合比较分析，当调节系数取值在1.0～3.0之间，可以较好地模拟海浪运动。

5 结语

本文主要讨论海浪真实感模拟问题。利用基于位置的动力学框架下的PBF方法进行流体粒子建模，再利用涡度约束的方法丰富流体表面细节。实验结果表明，本文提出的合理化系数涡度约束方法能有效的模拟海浪表面的翻卷与破碎，提高了海面模拟的真实感。但对于流体表面重建与绘制以及调节系数的最优化选择还需要做进一步的研究，从而继续提升模拟的真实感。