水下流场模拟可视化渲染绘制方法研究

1、水下流场模拟可视化渲染绘制方法研究摘要:在3D水下场景渲染中融入流场的模拟可视化是增强现实领域近年来研究工作的重点之一。水下流场模拟可视化 框架,包括流场模拟和可视化渲染绘制两部分。流场模拟基于N-S方程，温度为唯一的输入参数。为了更准确描述水体 温差、密度不均匀导致的对流现象，浮力作为流上下浮动的外力项；同时，通过波浪扰动项描述流和水面波动的一致性。 可视化渲染绘制部分，采用基于纹理的直接体渲染(Direct Volume Rendering，DVR)和一维的转移函数(transfer-function) 描述流场特征，展现时变流场数据(温度，速度，压强等)。框架通过GPU引擎对上述算法进行

2、加速，结合上帝之光实时模 拟可视化时变流场的扰动效果,满足场景实时交互可视化渲染绘制需求。关键词:水下场景；流场模拟；体渲染;可视化;增强现实穿过水面的光线会与水中粒子发生碰撞形成上帝之光和 焦散等奇妙的光学现象。同时温度、密度和压强等的变化会引 起水体内对流和扩散等现象。隐藏在这些现象下的本质就是 流体动力学。流模拟需要昂贵的计算开销，而且渲染绘制流体 不能依赖传统的基于光栅化渲染平台实现，因此实时流动态模 拟可视化渲染绘制仍然是一个具有挑战性的课题。本文构建水下流场模拟可视化框架，实现水下场景中流扰 动过程的模拟和可视化渲染绘制。模拟阳光穿透水面和水流 发生相互作用时产生的复杂的光学现象,

3、结合体渲染和后处理 体积光生成算法增强流特征。文章从下面几部分介绍相关工作：流场模拟可视化框架， 包括流场模拟N-S方程求解、基于纹理的直接体渲染,实验结 果分析、未来的改进方向等。1研究背景流场模拟技术可归纳为两个范畴:基于粒子的拉格朗日方 法和基于网格的欧拉方法。本文流场模拟采用基于网格的欧 拉方法。Stam是第一位采用网格欧拉方法模拟流现象的研究者。 欧拉流场模拟方法遵循克林形式(Chorin-style)，采用操作符分 割方法和对流-投影方法来对流施加不可压缩的条件限制。传 统的对流-投影框架采用半拉格朗日对流方法,该方法可提供 无条件的稳定性。近十年来许多基于该框架的方法被提出用 以

4、处理不同的流场模拟问题。Kang et al. 2引进了一个GFM (GhostFluidMethod)方法模拟多种状态的不可压缩流体，该方 法被Hong and Kim3采用对流体的表面变化进行追踪,采用半 拉格朗日对流和粒子集(PLS)扩展了该方法。Kang et al.|4扩展 了 GFM,模拟可融合和不可融合的流之间的相互作用。基于网 格的流场模拟方法也被用来模拟高粘度的流现象。Carlson et al.151在MAC(Mark And Cell)框架上采用一个隐式的粘度方程来 模拟高粘度的流和融化效果。Goktekin et al61采用黏性和塑性 系数模拟粘弹性的流现象。最近，L

5、arionov et al71在一个压力- 粘度求解器中采用基于网格的离散化方法来模拟高粘度的牛 顿流体。为兼顾性能和渲染质量，本文采用基于屏空间的体渲染方 法来渲染时变3D流体数据。2流场模拟可视化框架本文框架包含大量的纹理渲染技术，由可编程着色器实 现,N-S方程求解的流体动力学运算采用CUDA加速。为渲染 一个真实的水下场景以增强流体渲染的视觉冲击力，框架还包 括水面的建模和渲染、水下焦散和上帝之光的模拟。框架结构 如图1所示：图1流场模拟可视化框架2.1流体动力学方程求解本文采用基于网格的流场模拟方法对N-S方程进行求解。为了展现水下流和水面一致的波动，在N-S方程中加入波浪扰 动和浮

6、力两个外力项。波浪扰动项描述风对水体的推动而形 成的加速度,浮力项模拟水体中温差引起的对流现象。框架实 现的流场模拟被限制在一个3D的立方体区域内。假设流的体 积在水中的特定深度保持不变，水中流的运动可以由不可压缩 流体的N-S方程来描述：三=- r- =(1)采用对流-投影(advection-projection)方法对N-S方程进 行求解。将速度场u在欧拉网格上离散化，公式拆成对流、 扩散、外力三个部分。u(x-1 + 8t)= u (x - u (x-1)&)(2)式(2)采用半拉格朗日离散化方法计算公式中的第一项， 对流项(advection)。不仅作用于速度场,还可更新温度场计算

7、浮力项(buoyance)。通过亥姆霍兹分解(Helmholtz -Hodge de - composition)指定的速度场u分为与其最相近的无源场和一个 保守矢量场(旋度为零)：u = w + p(3)将一个矢量场分为两个矢量场的过程称为投影，通过泊松 方程求解：(4)(5)(4)(5)=2V p = u 式(1)中的第二项是扩散项(diffusion)普=yv2 u(6)根据半拉格朗日离散化方法，偏微分方程(6)写成下面的 形式：u (t + 8t) - y&V2u (t + 8t) = u(t)(7)式(7)是一个隐式的泊松方程。为了书写简单，采用b代替 u(t)，采用u(t+ &)代

8、替u,式(7)可以转换成能够采用雅各比迭代 求解的形式：uk +1=kkkkkkui + l.J. k + ui-l.J. k +k + ui,.i-1,k + k + 1 + u0k-1 + ba3(8) 在这里a= (dx)2/y&，3= a+ 6，k表示迭代次数。式(5)也 同样可以修改成其离散形式：pk +1= kkkkkkp + p + p + p + p + p bap i + 1-1-kp I - 1,1-kpt. j + 1, kp i-i - 1- kp I.I-k + 1p I.I-k - 1 a3(9) 这里b表示速度场的散度值，Mu, a和13分别等于-(dx )2 和

9、4。本文模拟的不是高粘性流体，所以将粘度系数设定 Y =0.001 o计算了 50次雅各比迭代对式(8)和式(9)进行求解. 公式的最后一项是外力项F。对外力F的计算分为两步 进行： TOC o 1-5 h z F =+ fbuoyance(10)其中 Me = bPh(k)W彳(11)Ph(k )和W分别为菲利普频谱和风速，用来生成FFT波浪：exp ( -1/( kL )2) -2Ph(k)=A(12)kL = w2/g是在风速w下可能的最大波浪，k是波浪向量的 模。对于流场模拟域内的每个体素(x,y,z)，体素的波浪向量 为(K K： )=(2nn/s, 2nm/s)，-N/2Wn WN

10、/2,-M/2 Wm WM/2，N和 M是三维网格的第一个维度和第二个维度的大小。在一个立 方体网格内离散求解流体动力学方程,M=N。K是单位化的波 浪向量，b是一个缩放系数。通过在式(11)的结尾乘以一个水 平方向的向量X=(i,j)，可以给每个体素内的速度值一个水平面 内的加速度。给每一个体素一个初始的速度值T,式(10)中的 浮力项可以由下面的式子计算：T0是水体的平均温度值,a是一个物理常量，在实验中需 要反复调整来取得比较好的视觉效果。j是一个归一化的垂直 方向向量，给予体素中的速度值一个垂直向上的加速度。由于 温度是流场模拟的唯一初始输入,本文的流场模拟对温度值较 为敏感。下面的结

11、果是用不同的温度输入值模拟10秒的结 果。图2中第一排的温度初始化采用的高斯随机函数均值为 10,第二排均值为20,网格分辨率为64 X 64 X 64。(a1,a2,a3)均值为10的高斯随机数，温度场扩散模拟结果(a1,a2,a3)均值为10的高斯随机数，温度场扩散模拟结果(a1- a2- a3) simulation of dye dispersion in a temperature field by Gauss-random function with initial mean 10(b1,b2,b3)均值为20的高斯随机数，温度场扩散模拟结果(b1- b2- b3) simulat

12、ion of dye dispersion in a temperature field by Gauss-random function with initial mean 20图2不同温度场的染料扩散情况流体数据集(速度和温度)存储在3D数组中，数据集的每 一帧采用公式(1)得到的速度场进行更新,采用DVR技术可视 化渲染时变流体数据。公式(14)计算沿视线方向累积光强和模 糊度，其中G是样本点颜色，&.是模糊度,用来模拟光线到达样 本点j的衰减状况：n i - Jc=Gn(i-&)(14)i=1I = 1Ai=1-11 (1 - AI)I = 1基于视点的体渲染将体数据的约束体采用代理几

13、何 (proxy geometries )进行切分，并按照式(14)进行采样和累积光 强，每个代理几何和视线方向垂直，如图3所示：图3代理几何法每一步代理几何片光栅化后的结果都存在一个帧缓存 中，整个体渲染过程总共采用两个帧缓存按照式(15)沿视线反 方向迭代的累积颜色和模糊值：Ci = (1 - A) 乂 Ci + G-1(15)A =( 1 - Ai - )x Ai + Ai -!架对渲染海面和水下光照采用的法线贴图方式实现参数化控 制。计算量最大的部分是流体动力学方程的求解。其中FFT 法线贴图的计算在Cuda上进行。在一张64 64的2D网格上 生成FFT波浪，在一张64 64 64的网格上模拟流和采用 DVR可视化流特征。为了兼顾效果和性能，在低分辨率的网格 上生成波浪高度场，采用高分辨率的法线贴图渲染水面。当法 线贴图的分辨率增加后，海面会呈现更多高频的波浪,水下也 可以看到更多波动的细腻光束，如图8所示。为了求解流体动 力学方程，本文采用30次雅各比迭代来求解扩散项，采用50次 迭代进行投影操作。根据实验结果，在计算扩散项时，采用20 次以上的迭代可以得到比较小的误差。在本文的模拟中，dx= dy = dz=1/64=0.0156,dt