陈立家 刘锭坤 田延飞 黄 巍

(武汉理工大学航运学院1) 武汉 430063) (内河航运技术湖北省重点实验室2) 武汉 430063) (长江引航中心江阴引航站3) 江阴 214400)

0 引 言

水流模拟的方法包括基于纹理贴图、几何形态、粒子系统、流体动力学、海浪波谱、柏林噪声,以及元胞自动机等仿真方法[1-3].对于不同类型水体的仿真，由于对物理真实感和渲染真实感以及实时性、交互性的要求不同，选用的模拟方法各有不同，使用过程中需结合实际情况进行必要的设置或改进.

在河流仿真与流体绘制的研究中，冶运涛等[4]采用基于纹理的流场动态可视化方法来模拟河道水流，综合了几何形状映射和颜色映射方法的优点.刘晨龙[5]基于Perlin噪声和海浪谱的方法绘制水面波浪效果，结合纹理坐标变换和图像变换的纹理动画技术实现河流流动，但该方法生成的河流物理特性较弱.余伟等[6]利用断面水位线结合DEM地形构建多层水面叠加模型，采取流速驱动纹理块自适应于河道移动的河流渲染方法，该方法能良好表达水体的自然特性，但在支流处无法智能实现断面处理.翟小明等[7]基于真实河道岸界数据求解河流速度场，通过瓦片式向量图实现河流表面的实时绘制，该方法计算简便，绘制过程基于GLSL着色器并在GPU中执行，体现出真实的河流流动方向.

在河流水体的三维模拟研究中，康苏海[8]基于非静压假设下的三维水动力数学模型，结合立体成像原理实现了三维流场可视化，并将其应用于长江复杂水域的三维水流求解，计算结果较为精确，但未注重水体渲染；冯小香等[9]基于平面曲线坐标系和垂向坐标系建立了三维水流数学模型，对弯曲河段和船闸引航道口门区的水流进行建模，模拟精度能满足通航水流条件研究的要求；扬程等[10]用非结构化网格对黄河沙坡头的连续弯道处水流运动进行研究，垂直方向上进行等距分层，采用有限体积法进行离散，较好的模拟了天然河道的水流运动.目前，在大规模场景的河流绘制研究上并未有较好的方法既体现水体的真实物理特性又兼顾实时、细致的渲染效果.

本文主要研究船舶操纵模拟器中对内河特殊航段的三维水流建模与仿真.为增强河流绘制的物理特性，选用三维水动力模型进行流场解算，以葛洲坝下游水域为研究对象，水平方向上采用结构化网格，垂向上采用σ坐标系进行非等距分层，运用多重纹理映射技术进行河流流体渲染，结合流场的三维特点对水面高度图及法向贴图、流动纹理贴图进行制作，并通过折射、反射纹理实现水面光照效应，整个过程在OpenGL平台中利用GLSL着色器进行流体渲染，利用GPU并行算法提高绘制质量和效率.

1 水动力三维数学模型

1.1 控制方程

控制方程基于水平边界拟合曲线坐标系和σ坐标系，转换关系为

(1)

笛卡尔坐标系下河流的三维紊流控制方程包括动量方程、连续方程、温盐方程和体积分数方程.

(2)

(3)

(4)

(5)

(6)

(7)

式中：x，y，z分别为笛卡尔直角坐标系的东、北和垂直方向的坐标；u，v，w分别为笛卡尔坐标下x，y，z方向上的流速分量；ρ为水温T和盐度S下的水密度；ρ0为水密度平均值，取999.7 kg/m3；g为重力加速度；ζ为自由表面水位；h为水深；H为实际水深，H=h+ζ；pa为大气压力；fc=2ωsinφ，为科里奥利系数；K，vT为水平和垂直方向上的涡流黏性系数；ωs为单位时间从外界流入的水量；c为水体积分数；w0为浮力引起的垂直方向上的速度；Ts，Ss，cs分别为水源的温度、盐度和体积分数；λ为非守恒系数；Kc，vc为水平和垂直方向上的扩散系数.

为节省运算时间，本次模型不考虑温度、盐度和浓度的影响，采用如下所示的忽略二阶水平对流项的标准k-ε模型求得垂直涡流黏性系数，其中：

(8)

(9)

(10)

(11)

(12)

式中：k为紊动动能；ε为紊动耗散率；Pr为普朗特数，Pr=0.96；σε=1.3;cμ，c1ε，c2ε，c3ε为经验常数，cμ=0.09，c1ε=1.44，c2ε=1.92，c3ε在稳定条件下，vTN>0，c3ε=0.001，在非稳定条件下，vTN<0，c3ε=2.

模型采用有限差分方法进行数值离散，流场在水平方向上采用矩形网格进行划分，垂直方向采用σ坐标系进行非等比例分层.

1.2 边界条件

1) 自由表面与底部边界条件[11].

在自由表面z=ζ(x,y,t)处：

(13)

(14)

式中：τs(x),τs(y)为表面风应力在x和y方向上的分量；ρa为空气密度，取1.225 kg/m3；W=(wx,wy)为水面上10 m高度处的风速；Cd为风应力拖曳系数，Cd=(0.63+0.066|W|)10-3.

在底部z=-h(x,y)处：

(15)

(16)

fb=[k/ln(zb/z0)]2

(17)

式中：τb(x),τb(y)为底部应力在x和y方向上的分量；fb为底部拖曳系数；κ为卡曼常数,κ=0.4；zb为无量纲层厚度；zo为底部粗糙度参数.

2) 固边界采用无滑移边界条件，无质量和热量交换，固边界处法向和切向速度为零.

3) 开边界条件包括水位、水量边界条件和流速边界条件.采用实测水位、水量和流速资料作为开边界条件.

4) 出水浅滩采用动边界处理，即对计算区域内滩地干湿过程采用水位判别法处理，设置判别水深为干水深δ1和湿水深δ2.当某单元水深小于δ1时，该处流速为零，滩地干出，不参与计算；单元水深在δ1和δ2之间时，只进行连续方程的计算；单元水深大于δ2时参与所有控制方程的计算.

2 水面渲染方案

本文基于GLSL着色器的多重纹理映射技术进行水面渲染，最终形成具有真实感的河流流体，渲染过程在GPU中执行，使用可编程渲染管线处理数据，以减小CPU运算量，提高计算效率，增强渲染效果.渲染流程见图1，首先将水动力模型计算出的水位信息结合三维水波算法[12]制作水流高度图，并利用凹凸纹理映射技术将水流高度图转化为法向贴图，经着色器渲染生成水波波纹；添加流动纹理贴图控制波纹的偏移，从而控制河流流向和大小，实现河流流动效果；最后结合菲涅尔效应生成的混合纹理实现水面光照效应.

图1 基于多重纹理映射技术的流动水面渲染

2.1 法向贴图

使用凹凸映射技术将高度图中储存的高度值进行转化，生成法向贴图，将法向贴图任一点周围的四个法向量混合叠加形成混合法向量，再经过着色器渲染即可形成水面波纹.

首先对水动力方程计算出的水位值插值得出水面上任一点高度值，生成初始水面高度；再根据三维水波算法，生成具备起伏变换效果的水波，并将水波高度与初始水面高度叠加，最终得出可实时变换的动态水面高度图.在顶点着色器中通过纹理坐标读取高度信息，并将其置换到水面模型中，可获得凹凸不平的波浪效果.

法向贴图中各点处光线的法向量值可以将相邻点高度值作差后叉乘得到，空间坐标系中法线由X，Y，Z轴上三个分量组成，法线信息储存在贴图的R，G，B三个通道中.

将高度图转化后的法线贴图进行法向量叠加，为

1) 已知点pa，将其坐标分别在x、y和沿45°角方向上进行位移，得到pb，pc，pd点的坐标.Normalpa、Normalpb、Normalpc、Normalpd分别为坐标点对应的法向量.

2) 将pa点与pb点法向量在x方向上进行混合生成新的法向量Normalpapb，将pc点与pd点法向量在x方向上混合生成法向量Normalpcpd，再将Normalpapb与Normalpcpd在y方向上混合得到混合法向量Normal，经过着色器渲染即可形成水面波纹.

3) 引入二维偏移变量，该偏移值大小为该任一点处x,y方向上的速度-时间函数，用该值去影响pa,pb,pc,pd坐标值，可得到随时间变化的波纹.

因高度图信息基于流场的水位信息，偏移变量与水面实际速度相关，因此水面波纹的效果与实际波动状态相符.

2.2 流动纹理

要实现水面的流动效果，还需要采用动态纹理技术，通过流动纹理贴图控制纹理的移动体现出实际水流的方向和流速.流动纹理的R,G,B,A各通道用于保存流体信息，R通道保存x方向的速度信息，G通道保存y方向的速度信息，通过流动纹理的R值和G值对法向纹理坐标的方向进行旋转，改变纹理映射的方向，从而改变波纹的流动方向；B通道保存河流的速度大小，通过B的值对纹理坐标进行缩放，可改变法向贴图在着色器中的映射数量，从而控制波纹大小；A通道为各点处混合法向量参数，通过A的值控制河道外混合法向量为零向量，对河道外部分不进行渲染，以控制河道边界.

2.3 菲涅尔效应混合纹理

在现实中，对水面进行观察时，观察者水面角度越小，反射效果越明显；角度越大，折射效果越明显，这种现象称为菲涅尔效应.为提高水面模拟的真实感，考虑到菲涅尔效应的影响，分别获取水面反射纹理与折射纹理，并将两者进行混合，使用顶点着色器和片元着色器对水面进行渲染，最后得到与实际水面光照效应接近的渲染效果.采用以下方式进行简单计算.

FresnelRatio=F+(1-F)×

pow((1-dot(eveDir,normal)),5)

(18)

finalColor=waterColor+reflectionColor×

FresnelRatio+refractionColor×

(1-FresnelRatio)

(19)

式中：F=0.200 37，为空气到水面的反射系数；eyeDir为视线方向；normal为光线与水面交点处的法线方向；finalColor为菲涅尔效应渲染的水面颜色；waterColor为初始水的顶点颜色；reflectionColor为水面的反射颜色；reflractionColor为水面的折射颜色.

3 案例研究

3.1 计算范围

模型计算区域为葛洲坝下游水域，包括主航道与三江航道，及其交汇水域，该处河段水流状况复杂，三江下引航道口门与大江来流夹角过大，存在的较大横流与回流使得船舶通航难度加大.因此，有必要对该类水文条件复杂的通航水域进行三维流场模拟，并开展水域内的船舶操纵模拟研究，从而提高内河航行安全.研究所用资料为长江勘测规划设计研究院在葛洲坝河段的河道水文数据，计算域范围为30°40′29.9N″～30°43′31.9″N，111°15′26.6″E～111°18′34.2″E，水域面积约6 km2，零水位处水深分布见图2.

图2 模型计算区域水深

3.2 计算条件

根据《三峡(正常运行期)—葛洲坝水利枢纽梯级调度规程》，葛洲坝大江1号船闸最大通航流量为35 000 m3/s，葛洲坝上游运行水位为66.0 m.模拟方案控制下游大江航道下泄流量为35 000 m3/s，同时二号船闸充水，控制三江航道流量为1 600 m3/s.试验采用Windows7操作系统，以OSG为视景渲染引擎，使用GLSL着色器语言进行河流流体渲染，OpenGL为4.3版.

计算网格平面上采用矩形网格，分别采用方案一30 m网格与方案二10 m网格尺度计算，见图3，其中CS1、CS2分别为两种网格在同一位置处的河流横截面，纬度为30°42′15.3″N.

图3 交汇河段计算网格

垂向网格划分为6层，其中表层与底层厚度较小，垂向网格划分方法见图4.

图4 垂向网格划分方法

3.3 参数设置及调整

参照河段地貌与水流阻力特征，河段模型糙率系数取值为全区域一致0.025；垂直涡流黏度系数最小值取分子黏度1.2×10-6；水平涡流黏度系数为参数化亚尺度湍流涡旋中能量耗散的经验系数，为有效防止数值震荡，该值大小为方案一1.1，方案二1.0.

考虑到对流项对水动力计算过程中的影响，时间步长取值受到条件限制为

(20)

式中：Δt为时间步长；Δx为网格步数；U为水流流速.根据该河段实际水流特征，时间步长取值为方案一1.5 s，方案二0.5 s.

3.4 结果分析

图5～6为主河段和三江河段交汇处的表层、中层、底层流场及回流区表层流场分布状况.垂向流速从表层到底层逐渐减小，河口交汇处存在斜流与回流.三江河段与主河段流速、流量相差较大，结合图7截面处流速的横向分布，可明显看出长江主流相对于三江航道中心线而言的斜流流态.沿三江流向的水流在交汇区域靠左岸处流速降低，并在河口处形成回流，回流以逆时针方向旋转.受网格精度影响，30 m网格中无明显回流流态，10 m网格中回流明显，汇流口流线清晰，汇流角明显.从实验结果来看，三维流场模拟精度与实际吻合度较高.

图5 河流交汇处流场分布(30 m网格)

图6 河流交汇处流场分布(10 m网格)

图7 CS2截面处横向流速分布

图8为交汇河口水面流动效果图，方案二湍流特性较方案一好，该处部分区域流速较大，水面波动状况复杂，方案一受网格精度影响，渲染无回流流态.由于回流区流速相较于主流流速太小，因而回流流动并不明显，回流区渲染效果体现为回流流速慢，流态相对平静.

图8 河流交汇处水面流动效果

本文所用方法的河流流体绘制过程在GPU中执行，10 m网格和30 m网格的河流流体渲染帧率平均分别在42帧和45帧左右，可实现场景的流畅运行.总体上水流渲染结合流场的物理特性，水体流动变换自然，体现出河流的真实流向，河口汇流处波动现象与实际符合度较好.

4 结 束 语

本文以葛洲坝下游水域为例，建立了内河河流三维数值流场，对坝后弯曲河段及主、支流交汇河口的复杂流态进行模拟.使用GLSL着色器进行流体渲染，以计算出的水位信息为基础生成水流高度图，以此制作了符合水面实际光影效果的法向贴图，并根据计算出的速度场信息制作了流动纹理贴图，最后结合菲涅尔效应混合纹理实现了河流流体绘制.

实验表明该方法计算出的水流仿真精度较高，能还原出斜流、回流等复杂流态，结构化网格的数据结构简单，模型整体的计算效率高.流体渲染通过GPU并行计算完成，可编程渲染管线能灵活的进行纹理映射，渲染过程结合流场的物理特性，生成的河面起伏波动效果较好，通过河流的速度场信息驱动纹理移动，能体现出实际的流速流向.总体来说，基于GLSL的多重纹理映射技术能较为真实的绘制出河流水体，生成与实际相符的水面光照效应和流动效果.

本文方法建立的三维水流模型有效提高了内河航道水流仿真精度，渲染方案能够实现水流模拟的真实感.对特定航段水流仿真的成果可应用于航道水流条件分析、内河船舶操纵模拟器系统建设、船舶安全操纵决策等方面.