复杂淹没区域风暴潮洪水演进三维动态全过程仿真方法

摘要

本发明属于洪水演进可视化仿真领域。为实现风暴潮洪水演进三维动态全过程仿真，实现海堤安全管理、风暴潮灾害的快速科学评估及科学制定抢险救灾实施方案，本发明采取的技术方案是，复杂淹没区域风暴潮洪水演进三维动态全过程仿真方法，包括如下步骤：(1)复杂淹没区域风暴潮洪水演进三维数学模型；(2)复杂淹没区域风暴潮洪水淹没水深图的绘制；(3)复杂淹没区域风暴潮洪水演进液相流的动态模拟；(4)复杂淹没区域风暴潮洪水演进动力学解算；(5)复杂淹没区域风暴潮洪水演进虚实耦合建模；(6)三维动态全过程仿真后期处理与可视化输出。本发明主要应用于风暴潮洪水演进三维动态全过程仿真。

说明书

技术领域

本发明属于复杂淹没区域洪水演进可视化仿真领域，涉及复杂淹没区域风暴潮洪水演进 的三维数值模拟及基于数值模拟的复杂淹没区域洪水演进三维动态全过程仿真方法。具体讲， 涉及基于数值模拟的复杂淹没区域风暴潮洪水演进三维动态全过程仿真方法。

背景技术

中国海洋灾害损失约占全部自然灾害总损失的10％，而风暴潮灾害又为海洋灾害之首。 随着CFD(Computational Fluid Dynamics，计算流体动力学)技术的广泛应用，风暴潮洪水演 进数值模拟已成为风暴潮灾害风险评估的重要手段之一。基于复杂淹没区域风暴潮洪水演进 数值模拟的洪水演进三维动态全过程仿真，在三维虚拟世界中实现对洪水的仿真体验，便于 对复杂淹没区域中洪水演进的现象和规律进行全方位的观察、操作和分析，更好地了解复杂 淹没区域洪水演进的全过程，对复杂淹没区域的风暴潮洪水演进到不同地区的时间、水深和 流速提供快速、准确和直观的预报。

国外关于洪水演进数值模拟的相关技术主要是利用一些以GIS(Geographic Information  System，地理信息系统)技术为基础开发出的商用软件包进行初步的洪水淹没与演进模拟， 相关的软件主要有DAMBRK、MIKE-11、HEC-2等等。这些技术在洪水演进与GIS集成方 面做了不少的工作，为洪水演进仿真技术的深层次开发提供了基础，但都不能够动态模拟复 杂淹没区域中流域洪水的三维淹没演进过程，因而不能直观、高效和实时动态地提供流域周 边经济带和生活区的防洪对策[1]。

在国内，董文锋[2]等应用了一维水动力学模型结合OpenGL(Open Graphics Library)和 GIS技术研制了清江流域“洪水演进模拟仿真系统”模拟流域洪水的淹没过程；韩敏[3-4]等 以给定水位的方式采用VC和OpenGL技术建立了洪水演进模拟系统，实现了洪水演进的可 视化；李云[5]等采用平面一维和二维水流数学模型结合GIS和VR(Virtual Reality，虚拟现 实)技术实现了淮河临淮岗大型行蓄洪区洪水演进的可视化；张成才[6]等以GIS技术为基础， 基于二维浅水方程，采用有限体积法，建立洪水演进数学模型，并将其制成动态链接库，通 过编程调用该库实现洪水淹没分析功能，计算出给定水位下的洪水淹没面积和水量，采用面 向对象的编程语言C#，结合ArcEngine组件库，开发了黄河东平湖洪水淹没场景三维可视化 系统，并成功应用于东平湖蓄滞洪区的洪水计算中；冶运涛[7]等采用网格逼近的方法生成了 三维河道地形，与大尺度流域三维场景形成嵌套结构，在与一维水动力学模型集成的基础上， 采用动态纹理技术和纹理坐标求解实现了大范围水面动态模拟。

综上所述，国内外针对洪水演进数值模拟的研究大多数都集中于一维或二维的水流数值 模拟计算，不能真实反映流场的细节变化，且均未涉及复杂淹没区域地形建模的洪水演进全 过程动态仿真，目前基于三维数值模拟的复杂淹没区域风暴潮溃堤洪水演进全过程仿真还未 见报道，缺乏对含有复杂地貌的大范围洪水演进的数值模拟。

参考文献：

[1]袁艳斌，王乘，杜迎泽，等.洪水演进模拟仿真系统研制的技术和目标分析[J].水电能源 科学，2001，19(3)：30-33.

[2]董文锋，袁艳斌，杜迎泽，等.流域三维仿真地形仿真及洪水演进动态模拟[J].水电能源 科学，2001，19(3)：37-39.

[3]韩敏，陈明.洪水演进模拟的三维可视化研究[J].计算机应用研究，2005，25(8)：1906-1907.

[4]韩敏，孟宇，顾峰峰，等.VR技术在扎龙湿地三维洪水仿真系统中的应用[J].计算机工程 与应用，2006，42(5)：185-188.

[5]李云，范子武，吴时强，等.大型行蓄洪区洪水演进数值模拟与三维可视化技术[J].水利 学报，2005，36(10)：1158-1164.

[6]张成才，常静，孙喜梅，等.基于GIS的洪水淹没场景三维可视化研究[J].北京师范大学 学报(自然科学版)，2010，46(3)：329-332.

[7]冶运涛，张尚弘，王兴奎.三峡库区洪水演进三维可视化仿真研究[J].系统仿真学报，2009， 21(14)：4379-4382.

发明内容

本发明旨在解决克服现有技术的不足，实现复杂淹没区域风暴潮洪水演进三维动态全过 程仿真，为海堤安全管理、风暴潮灾害的快速科学评估及科学制定抢险救灾实施方案提供理 论与技术支持，对提高风暴潮灾害应急处置能力，确保沿海地区防潮安全具有非常重要的理 论意义和实践价值。为达到上述目的，本发明采用的技术方案是：复杂淹没区域风暴潮洪水 演进三维动态全过程仿真方法，具体包括如下内容：

(1)复杂淹没区域风暴潮洪水演进三维数学模型；

(2)复杂淹没区域风暴潮洪水淹没水深图的绘制；

(3)复杂淹没区域风暴潮洪水演进液相流的动态模拟；

(4)复杂淹没区域风暴潮洪水演进动力学解算；

(5)复杂淹没区域风暴潮洪水演进虚实耦合建模；

(6)三维动态全过程仿真后期处理与可视化输出。

复杂淹没区域风暴潮洪水演进三维数学模型是采用耦合VOF法的三维κ-ε复杂淹没 区域风暴潮洪水演进数学模型：

在计算域内，VOF法，即Volume of Fluid法，流体体积函数法，对模型中的每一相引入 一个变量-相体积分数F，F是空间和时间的函数，即F＝F(x，y，z，t)，在离散网格内，F值 是网格内各相流体的体积与能够被流体通过的空间体积的比值，水气界面的追踪可通过求解 F的输运方程完成：

$\frac{\partial F}{\partial t}+u\frac{\partial F}{\partial x}+v\frac{\partial F}{\partial y}+w\frac{\partial F}{\partial z}=0$F ∈[0，1]

其中，F为相体积分数，无因次；t为时间，单位s；x，y，z为直角坐标系下的坐标分量， 单位m；u，v，w分别为直角坐标系下x，y，z方向的速度分量，单位m/s；模型的基本微分方程 包括连续性方程、动量方程、紊动能k方程和紊动能耗散率ε方程，分别表示如下：

连续性方程：

$\frac{\partial \rho }{\partial t}+\frac{{\partial \mathrm{\rho u}}_i}{{\partial x}_i}=0$

雷诺时均动量方程：

$\frac{\partial \left({\mathrm{\rho u}}_i\right)}{\partial t}+\frac{\partial \left({\mathrm{\rho u}}_i{u}_j\right)}{{\partial x}_j}=\frac{\partial P}{{\partial x}_i}+\frac{\partial }{{\partial x}_j}\left[(\mu +{\mu }_t)(\frac{{\partial u}_i}{{\partial x}_j}+\frac{{\partial u}_j}{{\partial x}_i})\right]$

紊动能k方程：

$\frac{\partial \left(\mathrm{\rho k}\right)}{\partial t}+\frac{\partial }{{\partial x}_j}\left({\mathrm{\rho u}}_{i}k\right)=\frac{\partial }{{\partial x}_i}\left[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{{\partial x}_i}\right]+G-\mathrm{\rho ϵ}$

紊动能耗散率ε方程：

$\frac{\partial \left(\mathrm{\rho ϵ}\right)}{\partial t}+\frac{\partial \left({\mathrm{\rho u}}_iϵ\right)}{{\partial x}_i}=\frac{\partial }{{\partial x}_i}\left[(\mu +\frac{\left({\mu }_t\right)}{{\sigma }_{ϵ}})\frac{\partial ϵ}{{\partial x}_i}\right]+C_{1ϵ}\frac{ϵ}{k}G{-C}_{2ϵ}\rho \frac{{ϵ}^2}{k}$

其中，μ_t为紊动粘性系数，C_μ为经验常数，取0.09；G为紊动能生成率， u_i、u_j为直角坐标系下的速度分量，单位m/s；x_i、x_j为直角坐标分 量，单位m；ρ为流体密度，单位kg/m³；μ为分子动力粘性系数，单位N·m/s；k为紊动 动能，单位m²/s²；ε为紊动能耗散率，单位m²/s³；P为修正压力，单位Pa；σ_k、σ_ε分 别为k、ε的紊流普朗特数，无因次；C_1ε、C_2ε为经验常数，无因次；流体密度ρ和分子动 力粘性系数μ是体积分数F的函数；

复杂研究区域模型边界条件的确定：

(1)进口边界条件：如果进口条件已知，给定速度如下式，湍流参数k和ε根据经验公 式计算得到：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}$

$k_{\mathrm{in}}={0.03V}_{\mathrm{in}}^2$

${ϵ}_{\mathrm{in}}=C_{\mu }{k}_{\mathrm{in}}^{1.5}/l_{\mathrm{in}}$

其中，为速度矢量，单位m/s；为进口速度标量，单位m/s；k_in为进口紊动动能， 单位m²/s²；ε_in为进口紊动能耗散率，单位m²/s³；C_μ为经验常数，取0.09；l_in为进口的 混合长度，单位m；设进口当量直径为D，则l_in＝0.5D；

如果不知道准确的来流条件，只知道来流的流量条件，那么可以使用均匀来流条件，则 进口速度如下式：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}=\frac{Q}{A}\overrightarrow{n}$

其中，为速度矢量，单位m/s；为进口速度矢量，单位m/s；Q为进口流量，单位 m³/s；A为进口面积，单位m²；为进口断面法向单位矢量；

实际上的进口断面法向并不一定就是流速的方向，可能相差还比较大，这时速度方向可 以给定在网格方向上，条件是网格为结构化网格，并设为网格的方向矢量，则进口速度给 定如下：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}=\frac{Q·\overrightarrow{m}}{(\overrightarrow{n}·\overrightarrow{m})A}$

(2)出口边界条件：通常在计算域的出口，各速度分量(u，v，w)以及k和ε均取为第二类齐次 边界条件，代表出口断面的法向，即：

$\frac{\partial \phi }{\partial n}=0,$φ＝u，v，w，k，ε

(3)固壁边界条件：在壁面上采用无滑移条件，即：直角坐标系下的各速度分量 u＝v＝w＝0；为计算近壁区的紊流，采用壁面函数法，设近壁点p到壁面的距离为y_p，则 近壁点p的速度u_p，紊动动能k_p和紊动能耗散率e_p的值由下式确定：

$\frac{u_p}{u_{\tau }}=\frac{1}{\kappa }\ln \left({\mathrm{Ey}}_p^{+}\right)$

$k_p=\frac{u_{\tau }^2}{\sqrt{C_{\mu }}},{ϵ}_p=\frac{u_{\tau }^3}{{\mathrm{\kappa y}}_p}$

式中，y_p为近壁点p到壁面的距离，单位m；u_p为近壁点p的流速，单位m/s；k_p和ε_p分别为近壁点p的湍流参数紊动动能和紊动能耗散率；E为经验系数，取9.0；κ为卡门常 数，取0.419；为距壁面的距离，无因次，$y_p^{+}=\frac{{\mathrm{\rho u}}_{\tau }{y}_p}{\mu }=\frac{{\mathrm{\rho C}}_{\mu }^{1/4}{k}_p^{1/4}{y}_p}{\mu };$u_τ为壁面摩阻流 速，τ_w为壁面的剪应力，ρ为流体密度，单位kg/m³；C_μ为经验常数，取 0.09；μ为分子动力粘性系数，单位N·m/s；

(4)压力边界条件：在计算域的边界上，压力应满足Neumann条件，即第二类边界条 件，为保证计算的稳定性，在规定的某一内点上，压力为给定值，且定义为参考压力P_ref；

(5)复杂下垫面边界条件：风暴潮洪水淹没区域的下垫面错综复杂，有城市建筑物、街 道、河道、村庄、农田；当风暴潮洪水向城市中心演进时，将遇到密集的建筑物，通常将密 度高的多个建筑物作为一个整体即称建筑群来模拟，单个建筑物可视为建筑群的特殊情况。

复杂淹没区域风暴潮洪水淹没水深图的绘制：基于复杂淹没区域风暴潮洪水演进三维精 细化数值模拟结果，按照划分网格的顶点得到不同时刻的空间坐标为(x，y，z)的数据类型， 通过数据结构存储构建水面网格的顶点坐标以及各顶点之间的拓扑关系，构成空间三角网； 首先通过将水流空间三角网与三维数字场景中的地形、建筑物等模型间的相交融合来解决遮 挡关系与混叠问题，实现水流与三维场景的精准、有效结合，其次将三维数值模拟结果通过 标量场的动态可视化方法来模拟包括水位、淹没面积的相关物理量的变化；标量场的动态模 拟，是利用颜色的变化来表示物理量随时间的变化过程：将物理量的变化范围与所选定颜色 的变化范围建立一一对应的函数关系，分别生成不同时刻的水位高程-颜色渐变曲面，通过 不同时刻水深图的连续更迭显示，实现复杂淹没区域风暴潮洪水演进过程中水流的包括淹没 范围、水位的物理量随时间的变化过程。

复杂淹没区域风暴潮洪水演进液相流的动态模拟：

基于复杂淹没区域风暴潮洪水演进三维水流精细化数值模拟结果，划分网格的顶点得到 不同时刻的空间坐标，为(x，y，z)的数据类型，每相邻的4个点组成一个四边形，再将所 有的四边形联接成网，形成反映高低起伏变化的三维交错网格；经过以上处理后，初始显示 流场质点，进行质点流速场的绘制：表示各质点流速矢量的示踪质可以用球体、线段或箭头 表示，选择适当的时间间隔，逐时生成一系列不同时刻流体质点所在的位置及其速度矢量图， 就得到了动态的流场；复杂淹没区域风暴潮洪水演进过程中液相流运动模型采用GPGPU即 General-Purpose Computing on Graphics Processing Unit，通用计算图形处理单元技术，实现水 流示踪质即粒子间的交互。

复杂淹没区域风暴潮洪水演进动力学解算：

根据复杂淹没区域风暴潮洪水演进三维仿真的具体特点，应用Autodesk公司3ds max软件 和以NURBS为理论基础的3D造模软件Rhinoceros制作场景模型，利用程序接口将制作成型的场 景模型导入到NextLimit技术公司的Realflow中，设置模型属性模拟真实场景，设置粒子发 射器和网格模拟水流，设置包括重力、风力、表面张力、摩擦力的辅助器模拟真实场景中的 各种物理力场，进而对选定的动水场景进行动力学解算；将Realflow中最终解算所得到的水 流效果即网格或粒子及场景模型参数导回到制作场景模型的3ds max软件中，并在此基础上 进行灯光照明、赋予材质和渲染等后期处理工作，从而顺利地实现逼真的动态水流效果模拟。

复杂淹没区域风暴潮洪水演进虚实耦合建模：

基于数值模拟结果分析复杂淹没区域风暴潮洪水演进的空间位置，以确定研究区域，进 而提取相应的三维模型数据实现复杂淹没区域的三维场景中的地形与建筑物、绿化设施之间 的匹配，将数值模拟结果所得到的水情信息，在三维场景中动态的表现出来；基于复杂淹没 区域的三维场景中风暴潮洪水演进过程的模拟效果，选取具有典型意义的标志性建筑物，并 前往实地捕获大量现场图片以得到详细、全面的信息数据，保证视景仿真的精度与可靠性， 进而在现场图片的基础上耦合动态的水流，实现动静、虚实的结合，逼真的模拟复杂淹没区 域风暴潮的洪水演进过程；

复杂淹没区域风暴潮洪水演进虚实耦合建模过程中，基于复杂淹没区域的三维模型场景 采用纹理变换的方法，即通过动态改变水流的纹理信息模拟水流的运动，实现动态水流的三 维可视化，模拟出水流沿地势高低向前演进的过程，反映水位的起伏变化。

三维动态全过程仿真后期处理与可视化输出：

计算机动画是指用程序或工具生成一系列的静态画面，然后通过画面的连续播放来反映 对象的连续变化过程，动画中的运动不仅指物体的运动，还包括视点的移动、光照的变换、 纹理和色彩的变化；

在制作动画之前，先根据要求制作脚本，要考虑好动画的主要内容是什么，如何布置重 要的场景，动作如何描述，设置音响效果；进而对动画场景赋予适当的材质和贴图，使场景 具有真实感，必要时还可以通过材质参数的关键帧设置实现动画效果；调整好材质和贴图后， 设置场景光源类型、强度、色彩及其相应的位置和角度，以烘托良好的环境氛围；然后根据 精度要求，设置好需要的分辨率，渲染出图片和动画文件；在制作过程中，考虑到渲染时间 较长，故应该将整个动画分为几块来进行渲染。

本发明的技术特点及效果：

(1)与已有技术相比，本发明建立了耦合VOF法的三维κ-ε复杂淹没区域风暴潮洪 水演进数学模型，实现复杂淹没区域风暴潮洪水演进的精细化数值模拟，结果更加符合真实 场景下洪水波的流动规律，为复杂淹没区域风暴潮洪水演进提供基础水情信息；

(2)基于复杂淹没区域风暴潮洪水演进三维数值模拟结果，实现了不同频率、不同历时 下的复杂淹没区域风暴潮洪水淹没水深图的绘制；

(3)基于液相流动态模拟、虚实耦合建模等技术，实现了兼顾真实性和科学性的复杂淹 没区域的水流动态模拟，使得复杂淹没区域风暴潮洪水演进模拟的制作流程更简单，图形显 示技术更先进，可视化结果更逼真。

(4)利用洪水演进动力学解算，制作出了复杂淹没区域风暴潮洪水演进过程中海浪拍打 堤岸、风暴潮来袭、巨浪翻过海堤向近海港口城市涌动、涌入城市的海浪撞击建筑物、质量 较轻的物体漂浮在水面上并随水流一起流动等全部动水效果。

附图说明

图1基于数值模拟的复杂淹没区域风暴潮洪水演进三维动态全过程仿真方法的总体技 术流程。

图2复杂淹没区域的淹没水深图。

图3GPGPU中的纹理数据。

图4复杂淹没区域三维洪水演进模拟解算流程。

图5虚实耦合建模技术数据流分析。

图6复杂淹没区域风暴潮洪水演进动态全过程仿真的三维演示原理流程。

具体实施方式

本发明基于耦合VOF法的三维κ-ε紊流方程对复杂淹没区域风暴潮洪水演进过程进 行了精细化三维数值模拟，以所建立的研究区域的三维数字模型为载体，结合淹没水深图绘 制、液相流动态模拟、洪水演进动力学解算、虚实耦合建模等技术，实现复杂淹没区域风暴 潮洪水演进的三维全过程动态仿真，对复杂淹没区域风暴潮洪水演进到不同地区的时间、水 深和流速提供快速、准确和直观的预报。

该方法具体包括：

(1)复杂淹没区域风暴潮洪水演进三维数学模型；

(2)复杂淹没区域风暴潮洪水淹没水深图的绘制；

(3)复杂淹没区域风暴潮洪水演进液相流的动态模拟；

(4)复杂淹没区域风暴潮洪水演进动力学解算；

(5)复杂淹没区域风暴潮洪水演进虚实耦合建模；

(6)三维动态全过程仿真后期处理与可视化输出。

下面结合附图和具体实施实例进一步详细说明本发明。

(一)基于数值模拟的复杂淹没区域风暴潮洪水演进三维动态全过程仿真方法的总体技术流 程

总体技术流程见图1，具体包括如下步骤：

(1)建立了耦合VOF法的三维κ-ε复杂淹没区域风暴潮洪水演进数学模型，对复杂 淹没区域的风暴潮洪水水气两相流进行精细化模拟，根据实际情况确定进口、出口、固壁、 压力、复杂淹没区域下垫面等边界条件。

(2)利用已获得的监测数据、历史数据、实时数据等构建复杂淹没区域风暴潮洪水演进 信息数据库。

(3)通过数值模拟软件，在对研究区域进行网格划分的基础上，结合复杂淹没区域风暴 潮洪水演进三维数学模型，调用复杂淹没区域风暴潮洪水演进信息数据库，完成复杂淹没区 域风暴潮洪水演进的三维数值模拟。

(4)利用工程设计CAD图和文档资料、纸质地图、电子地图、卫星遥感图像等资料， 建立复杂研究区域地形、防潮工程以及周边环境，如建筑群、农田、水域等的三维数字模型， 以此构建出复杂淹没区域风暴潮洪水演进三维数字场景。

(5)结合复杂淹没区域风暴潮洪水演进数值模拟成果与复杂淹没区域的三维可视化仿真 场景，通过复杂淹没区域的淹没水深图绘制、液相流动态模拟、洪水演进动力学解算、虚实 耦合建模等技术，实现复杂淹没区域风暴潮洪水演进三维动态全过程仿真。

(二)复杂淹没区域风暴潮洪水演进三维数学模型

采用耦合VOF法的三维κ-ε复杂淹没区域风暴潮洪水演进数学模型。

在计算域内，VOF法，即Volume of Fluid法，流体体积函数法，对模型中的每一相引入 一个变量-相体积分数F，F是空间和时间的函数，即F＝F(x，y，z，t)，在离散网格内，F值 是网格内各相流体的体积与能够被流体通过的空间体积的比值。水气界面的追踪可通过求解 F的输运方程完成：

$\frac{\partial F}{\partial t}+u\frac{\partial F}{\partial x}+v\frac{\partial F}{\partial y}+w\frac{\partial F}{\partial z}=0$F∈[0.1]

其中，F为相体积分数，无因次；t为时间，单位s；x，y，z为直角坐标系下的坐标分量，单位 m；u，v，w分别为直角坐标系下x，y，z方向的速度分量，单位m/s。

模型的基本微分方程包括连续性方程、动量方程、紊动能k方程和紊动能耗散率ε方程， 分别表示如下。

连续性方程：

$\frac{\partial \rho }{\partial t}+\frac{{\partial \mathrm{\rho u}}_i}{{\partial x}_i}=0$

雷诺时均动量方程：

$\frac{\partial \left({\mathrm{\rho u}}_i\right)}{\partial t}+\frac{\partial \left({\mathrm{\rho u}}_i{u}_j\right)}{{\partial x}_j}=\frac{\partial P}{{\partial x}_i}+\frac{\partial }{{\partial x}_j}\left[(\mu +{\mu }_t)(\frac{{\partial u}_i}{{\partial x}_j}+\frac{{\partial u}_j}{{\partial x}_i})\right]$

紊动能k方程：

$\frac{\partial \left(\mathrm{\rho k}\right)}{\partial t}+\frac{\partial }{{\partial x}_j}\left({\mathrm{\rho u}}_{i}k\right)=\frac{\partial }{{\partial x}_i}\left[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{{\partial x}_i}\right]+G-\mathrm{\rho ϵ}$

紊动能耗散率ε方程：

$\frac{\partial \left(\mathrm{\rho ϵ}\right)}{\partial t}+\frac{\partial \left({\mathrm{\rho u}}_iϵ\right)}{{\partial x}_i}=\frac{\partial }{{\partial x}_i}\left[(\mu +\frac{\left({\mu }_t\right)}{{\sigma }_{ϵ}})\frac{\partial ϵ}{{\partial x}_i}\right]+C_{1ϵ}\frac{ϵ}{k}G{-C}_{2ϵ}\rho \frac{{ϵ}^2}{k}$

其中，μ_t为紊动粘性系数，C_μ为经验常数，取0.09；G为紊动能生成率， u_i、u_j为直角坐标系下的速度分量，单位m/s；x_i、x_j为直角坐标分 量，单位m；ρ为流体密度，单位kg/m³；μ为分子动力粘性系数，单位N·m/s；k为紊动 动能，单位m²/s²；ε为紊动能耗散率，单位m²/s³；P为修正压力，单位Pa；σ_k、σ_ε分 别为k、ε的紊流普朗特数，无因次；C_1ε、C_2ε为经验常数，无因次；流体密度ρ和分子动 力粘性系数μ是体积分数F的函数。

复杂研究区域模型边界条件的确定：

(1)进口边界条件：如果进口条件已知，给定速度如下式，湍流参数k和ε根据经验公式 计算得到：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}$

$k_{\mathrm{in}}={0.03V}_{\mathrm{in}}^2$

${ϵ}_{\mathrm{in}}=C_{\mu }{k}_{\mathrm{in}}^{1.5}/l_{\mathrm{in}}$

其中，为速度矢量，单位m/s；为进口速度矢量，单位m/s；k_in为进口紊动动能， 单位m²/s²；ε_in为进口紊动能耗散率，单位m²/s³；C_μ为经验常数，取0.09；l_in为进口的 混合长度，单位m；设进口当量直径为D，则l_in＝0.5D。

如果不知道准确的来流条件，只知道来流的流量条件，那么可以使用均匀来流条件，则 进口速度如下式：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}=\frac{Q}{A}\overrightarrow{n}$

其中，为速度矢量，单位m/s；为进口速度矢量，单位m/s；Q为进口流量，单位 m³/s；A为进口面积，单位m²；为进口断面法向单位矢量。

实际上的进口断面法向并不一定就是流速的方向，可能相差还比较大，这时速度方向可 以给定在网格方向上，条件是网格为结构化网格，并设为网格的方向矢量，则进口速度给 定如下：

$\overrightarrow{V}={\overrightarrow{V}}_{\mathrm{in}}=\frac{Q·\overrightarrow{m}}{(\overrightarrow{n}·\overrightarrow{m})A}$

(2)出口边界条件：通常在计算域的出口，各速度分量(u，v，w)以及k和ε均取为第二 类齐次边界条件，代表出口断面的法向，即：

$\frac{\partial \phi }{\partial n}=0,$φ＝u，v，w，k，ε

(3)固壁边界条件：在壁面上采用无滑移条件，即：直角坐标系下的各速度分量 u＝v＝w＝0；为计算近壁区的紊流，采用壁面函数法，设近壁点p到壁面的距离为y_p，则 近壁点p的速度u_p，紊动动能k_p和紊动能耗散率e_p的值由下式确定：

$\frac{u_p}{u_{\tau }}=\frac{1}{\kappa }\ln \left({\mathrm{Ey}}_p^{+}\right)$

$k_p=\frac{u_{\tau }^2}{\sqrt{C_{\mu }}},$${ϵ}_p=\frac{u_{\tau }^3}{{\mathrm{\kappa y}}_p}$

式中，y_p为近壁点p到壁面的距离，单位m；u_p为近壁点p的流速，单位m/s；k_p和e_p分别 为近壁点p的湍流参数紊动动能和紊动能耗散率；E为经验系数，取9.0；κ为卡门常数， 取0.419；为距壁面的距离，无因次，$y_p^{+}=\frac{{\mathrm{\rho u}}_{\tau }{y}_p}{\mu }=\frac{{\mathrm{\rho C}}_{\mu }^{1/4}{k}_p^{1/4}{y}_p}{\mu };$u_τ为壁面摩阻流速， τ_w为壁面的剪应力，ρ为流体密度，单位kg/m³；C_μ为经验常数，取0.09； μ为分子动力粘性系数，单位N·m/s。

(4)压力边界条件：在计算域的边界上，压力应满足Neumann条件，即第二类边界条 件。为保证计算的稳定性，在规定的某一内点上，压力为给定值，且定义为参考压力P_ref；

(5)复杂下垫面边界条件：风暴潮洪水淹没区域的下垫面错综复杂，有城市建筑物、街 道、河道、村庄、农田等。当风暴潮洪水向城市中心演进时，将遇到密集的建筑物，通常将 密度高的多个建筑物作为一个整体即称建筑群来模拟，单个建筑物可视为建筑群的特殊情况。

(三)复杂淹没区域风暴潮洪水淹没水深图的绘制

风暴潮洪水在复杂淹没区域的演进过程中，水位的升降、淹没范围的变化以及流速场的 迁移等现象，在本质上都可视为研究区域内水流模型三维网格节点处相关物理量的函数值随 时间变化的过程。为此，基于复杂淹没区域风暴潮洪水演进三维数值模拟结果，按照数值模 拟计算中所划分网格的顶点得到其在不同时刻的空间坐标，为(x，y，z)的数据类型，通过 数据结构可以存储构建水流网格的顶点坐标以及各顶点之间的拓扑关系，构成空间三角网。 首先通过将水流空间三角网与三维数字场景中的地形、建筑物等模型间的相交融合来解决遮 挡关系与混叠问题，实现水流与三维场景的精准、有效结合，其次将三维数值模拟结果通过 标量场的动态可视化方法来模拟水位、淹没面积等相关物理量的变化。标量场的动态模拟， 通常利用颜色的变化来表示物理量随时间的变化过程，主要是将物理量的变化范围与所选定 颜色的变化范围建立一一对应的函数关系。以水位图的标量场的实现为例，介绍复杂淹没区 域风暴潮洪水演进过程中淹没范围内的水位高程与颜色的映射关系。首先确定水深的变化范 围为Z₀-Z_n，设定起始颜色三元组(startRed，startGreen，startBlue)和终止颜色三元组(endRed， endGreen，endBlue)，将Z₀与RGB(startRed，startGreen，startBlue)相对应，Z_n与RGB(endRed， endGreen，endBlue)相对应。再根据实际需要将起始颜色和终止颜色等分成n份，则三元颜色 递增量分别为redInc、greenInc、blueInc，即：

redInc＝(endRed-startRed)/n；

greenInc＝(endGreen-startGreen)/n；

blueInc＝(endBlue-startBlue)/n。

同样水深高程在Z₀-Z_n范围内也划分成n份，将每段高程按顺序依次对应每组颜色三元组， 可得到水深高程与从起始颜色到终止颜色的渐变色带间的映射关系。为了使显示结果更为精 细美观，可分别给定起始颜色和终止颜色为一组色系，在两组色系间实现颜色的渐变，这样 所得到的渐变色带能够满足更复杂、更精准的显示要求。对于所建立的复杂淹没区域的水流 模型空间三角网，该网格仅能提供网格节点上的水位高程，而对于网格内部面域和边界线的 各点高程可以通过插值方法得到。首先基于建立的水位高程值和渐变色带的映射函数关系， 逐个三角网格作为研究对象，将各网格顶点的数据为基准，通过差值方法得到水流模型上每 一点的水位高程及其相对应的颜色像素值，并在网格内部实现颜色的渐变，从而得到某一时 刻完整的水位高程与颜色间的映射曲面，即风暴潮洪水淹没范围的水深图，如图2所示。依照 此方法及相应的高程颜色对应函数关系，分别生成不同时刻的水位高程-颜色渐变曲面，通 过不同时刻水深图的连续更迭显示，可实现复杂淹没区域风暴潮洪水演进过程中水流的淹没 范围、水位等物理量随时间的变化过程。

(四)复杂淹没区域风暴潮洪水演进液相流的动态模拟

基于复杂淹没区域风暴潮洪水演进三维数值模拟结果，按照数值模拟计算过程中所划分 网格的顶点得到不同时刻的空间坐标，为(x，y，z)的数据类型，每相邻的4个点组成一个 四边形，再将所有的四边形联接成网，建立反映水流高低起伏变化的三维交错网格。对于复 杂淹没区域风暴潮洪水演进过程中流速场的可视化模拟，各质点的流速矢量可以通过球体、 线段或箭头等示踪质实现直观、形象的表达。以箭头表示为例进行详细说明，可将箭头的方 向、长短以及颜色等属性与速度的相应矢量建立对应的映射关系，箭头的方向表示为速度在 三维空间上的运动方向，箭头的长短表示为流速值的大小，箭头颜色的深浅同样可以表示流 速值的大小。在进行复杂淹没区域风暴潮洪水演进液相流动态模拟的过程中，基于示踪质逐 时生成一系列不同时刻的矢量图，其能逼真形象的展示流体质点的空间位置及其速度，并依 时间顺序实现各时刻流速矢量图的更迭显示得到了动态的流场。动态的流场是时间与速度矢 量二者的函数表达，能显示各质点的空间属性随时间的变化过程，跟踪各流体质点的运动轨 迹。

复杂淹没区域风暴潮洪水演进过程中液相流运动模型的复杂度很高，通过采用GPGPU (General-Purpose Computing on Graphics Processing Unit，通用计算图形处理单元)技术实现 水流示踪质(粒子)间的交互，可以极大地加快液相流运动的计算模拟速度，进而缩短模型 的运算时间、提高液相流运动过程模拟的处理性能与效率。基于GPGPU加速图像处理的过程 中，数据是以纹理的形式存储的，针对水流粒子的速度和空间位置等属性分别建立了相应的 纹理，如图3所示。

由于液相流动态模拟的过程中任意时刻各粒子的属性值为上一时刻的数据更新替换而 成，采用Ping-Pong传递技术实现当前的输入纹理属性与上一时刻输出纹理属性的不断迭代， 从而实现水流示踪质属性信息的更新。在每一次的迭代过程后，都需要进行水流示踪质的空 间定位并进行排序，首先实现对液相流模型的空间剖分，再根据示踪质此时刻所处的空间位 置将其映射到三维空间模型网格之中，然后按照网格的编号对示踪质进行排序。基于上述方 法实现了复杂淹没区域风暴潮洪水演进液相流的动态模拟。实践证明，GPGPU技术将液相流 运动模型的计算模拟速度大大提高，且达到了良好的视觉效果。

(五)复杂淹没区域风暴潮洪水演进动力学解算

根据复杂淹没区域风暴潮洪水演进三维仿真的具体特点，采用Realflow(由NextLimit 技术公司开发的3D可视化的流体模拟和物理动态模拟软件)来快速真实地实现动水模拟效 果。Realflow软件是一种建立在流体动态计算技术上的物理粒子系统，是目前解决流体动力 学模拟的一套比较完善的解决方案。其软件内部的重力、风力、表面张力等十几种辅助器能 够完美再现真实物理世界中的力场，使模拟的效果非常真实。但是Realflow软件本身没有照 明和渲染模块，材质模块也相对简陋。因此为了制作出精彩的视觉效果，需要在3ds max(由 Autodesk公司开发的基于PC系统的三维动画渲染和制作软件)中进行后期的灯光照明、赋 予材质和渲染工作。通过程序接口就可以方便快捷地实现Realflow与3ds max等三维制作软 件之间的衔接。从而可以顺利地实现逼真的动态水流效果模拟。复杂淹没区域风暴潮洪水演 进三维仿真采用的动水效果可视化模拟流程如图4所示，技术路线如下：

(1)选择动水场景

在制作动水效果之前，首先应确定需要表现动水模拟效果的场景。复杂淹没区域风暴潮 洪水演进过程中需要进行模拟的动水效果主要有海浪拍打堤岸、风暴潮来袭、巨浪翻过海堤 向近海港口城市涌动、涌入城市的海浪撞击建筑物、质量较轻的物体漂浮在水面上随水流一 起流动等。

(2)制作场景模型

应用3ds max，Rhinoceros(以NURBS为理论基础的3D造模软件)等软件制作场景模 型，不同的场景其模型的精细程度略有不同，这都是由所需表达的动水效果决定的，例如， 当所需表达的动水效果为海浪涌入城市的大场景时，其场景模型的精细程度可以降低；当所 需表达的动水效果是海浪拍打海堤或是海水撞击建筑物等诸如此类的小场景时，其场景模型 的精细程度就需要提高。

(3)导入Realflow，设置参数

应用程序接口将制作成型的场景模型导入到Realflow软件中，设置模型属性模拟真实场 景，设置粒子发射器和网格模拟水流，设置重力、风力、表面张力、摩擦力等辅助器模拟真 实场景中的各种物理力场，有时还需要设置约束器模拟各种环境约束和边界条件，初设各种 参数，准备解算。

(4)调整参数，多次解算

这一步是应用Realflow进行动水模拟的最重要也是最耗费时间的一个步骤。根据初设参 数，进行初次解算，分析解算成果，调整相关参数，再进行解算，反复解算调整多次，直到 解算效果满足要求时为止。初期解算时可以把一些参数(例如网格数量，解算精度等)设置 的精度稍粗略一些，这样可以大大地节省解算过程中所耗费的时间，进而进行参数调整，当 解算成果大体上满足要求时，再把部分参数设置的精细些，进行最后的解算和成果输出。

(5)导入3ds max，后期处理

将Realflow中最终解算所得到的水流效果(网格或粒子)及场景模型参数导回到制作场 景模型的3ds max软件中，为了使动水模拟达到更加逼真的效果可将Realflow中的部分网格 在3ds max软件中进行替换，转化为具有雾化、水花等特效的粒子系统，并在此基础上进行 灯光照明、赋予材质和渲染等后期处理工作，使场景及动水效果模拟得更加逼真。

(六)复杂淹没区域风暴潮洪水演进虚实耦合建模

针对复杂淹没区域风暴潮洪水演进三维数值模拟的计算结果得到的洪水演进淹没面积、 水深、流速等水情信息，所采取的常规的分析手段，如水位变化过程线、流速(流量)变化 过程线、瞬时等深线和瞬时流场图等已不能满足当前数学模型发展和水利信息化的要求，有 必要利用一些新的技术和手段开发出基于数值模拟结果的复杂淹没区域风暴潮洪水演进三维 可视化系统。应用虚实耦合建模技术手段，可以达到较好的三维可视化仿真的效果。基于数 值模拟结果分析复杂淹没区域风暴潮洪水演进的空间位置，以确定研究区域，进而提取相应 的三维模型数据实现复杂淹没区域的三维场景中的地形与建筑物、绿化设施等之间的匹配， 将数值模拟结果所得到的水情信息，在三维场景中动态的表现出来。基于复杂淹没区域的三 维场景中风暴潮洪水演进过程的模拟效果，选取具有典型意义的标志性建筑物，并前往实地 捕获大量现场图片以得到详细、全面的信息数据，保证视景仿真的精度与可靠性，进而在现 场图片的基础上耦合动态的水流，实现动静、虚实的结合，逼真的模拟复杂淹没区域风暴潮 的洪水演进过程。虚实耦合建模技术数据流分析如图5所示。

复杂淹没区域风暴潮洪水演进虚实耦合建模过程中，基于复杂淹没区域的三维模型场景 采用纹理变换的方法进行动态水流的三维可视化，实现水流沿地势高低向前演进的模拟过程， 以反映水位起伏变化，进而逼真展示淹没建筑物的过程等视觉效果以及水流的流速流态等物 理特性。纹理变换就是通过动态改变水流的纹理信息而模拟水流的运动，其分两种情况，纹 理贴图位置的移动及多个纹理图片的替换。

(1)纹理移动方法

纹理移动方法是指对于某三维场景研究区域内的水流模型，其由多个相同的面叠加组成， 将每一个面进行同一图片的纹理贴图，而对于纹理贴图的位置应设定一定距离的间隔，以保 证水面的动态效果，并依据纹理贴图位置的不同对各个面进行依次排序，通过循环显示实现 水流动态的模拟。针对复杂淹没区域风暴潮洪水演进过程中水流的动态模拟，首先基于精细 化数值模拟结果得到的水情信息将每一场景下的水流模型分成若干区域，将每一区域水流模 型都依据三维场景可视化模拟的时长复制n份，并叠加显示，然后将n个水流模型赋予同一纹 理图片但纹理贴图的位置不同，若此区域水流模型的长度为l，则每两个面之间的纹理贴图位 置差距为l/n，进而依据纹理贴图位置的顺序实现n个水面模型的循环显示，耦合各个区域水 流模型的动态模拟结果即可完成复杂三维场景研究区域内水流流动的效果。

(2)多纹理替换方法

多纹理替换是指对于某三维场景研究区域内的水流模型，设定多张纹理图片用于赋予水 流模型，在三维场景的显示时长内依次将各张纹理图片赋予水流模型，并按顺序显示以实现 水流的动态模拟。对于复杂淹没区域风暴潮洪水演进过程中水流的动态模拟，首先将每一场 景下的水流模型分别赋予n张纹理图片，若此区域水流模型的显示时长为t，则每两张纹理贴 图赋予水流模型的时间间隔为t/n，进而依据纹理贴图的时间顺序实现水流模型的渲染输出显 示，以完成复杂三维场景研究区域内水流模型流动的效果。在实现水流模型的动态渲染输出 时，应指定各向同性，以保证每两张纹理图片在更迭替换过程中既满足水流流动的动态效果， 又能确保水流模型在经多次纹理贴图后水流的运动形态无跳变。

(七)三维动态全过程仿真后期处理与可视化输出

计算机动画是指用程序或工具生成一系列的静态画面，然后通过画面的连续播放来反映 对象的连续变化过程。动画中的运动不仅指物体的运动，还包括视点的移动、光照的变换、 纹理和色彩的变化等。3ds max软件提供了强大的三维图像处理功能，在动画实现的过程中， 只需要控制空间物体的坐标点变化和关键帧等，其余功能由软件实现，使仿真实现过程大大 简化。复杂淹没区域风暴潮洪水演进动态全过程仿真的三维演示流程如图6所示。

3ds max的动画效果都是基于其参数的变化，在关键帧中通过位置参数、尺寸参数等参数 的变化来实现运动，形状的变化等动画效果。以下是进行动画制作的几个要点。

(1)制作脚本。在动画制作之前，根据要求制作脚本，是非常重要的第一步。要清楚动 画的主要内容是什么，如何布置重要的场景，动作如何描述，设置音响效果等。由于渲染的 时间一般较长，因此，在最终生成动画之前，需要通过预览、部分渲染等方法首先生成低分 辨率的动画，以观察其效果是否达到要求。

(2)材质。材质编辑是3ds max动画制作的重要环节之一。材质描述物体的表面特性， 影响物体表面的颜色、光亮程度和不透明度等。正确运用材质和贴图，才能使物体具有真实 感。还可以通过材质参数的关键帧设置实现动画效果，如水的波动和流动效果。

(3)灯光。调好材质后，必须将场景中配上光才能看出真实效果。材质和灯光是相互联 系的。设置好场景光源类型、强度、色彩及其相应的位置和角度，才能创造出良好的环境氛 围。有时需要对某些物体单独加光的效果，比如在制作水流时，为了更真实，要对制作水的 粒子系统单独加入较强的白灯光，如果这种灯光作用于其他物体必定反光很强烈，因此需要 让灯单独作用于水，3ds max自带的灯光INCLUDE/EXCLUDE(包含/排除)功能可以有效地 解决这个问题。

(4)对场景设置背景。好的背景能烘托出良好的氛围，表现出实际的工程效果。背景图 的获取有很多种方法，可以直接使用3ds max自带的图片，也可以在Photoshop中制作自己需 要的效果图作为背景图。

(5)动画运动的设置。3ds max中设置运动的方法有关键帧法、变形法和关节法，较常 用的是使用关键帧法。通过动画约束，可以使摄像机始终对着某一物体，可以使物体只能沿 着某条路径或者某个表面运动，这可以使摄像机沿着某条光滑线运动，摄像机的运动比较平 滑，使动画不会出现抖动的效果。

(6)渲染。建立模型后，下一步的工作就是完成平面效果图以及仿真动画的制作。根据 精度要求，设置好需要的分辨率，渲染出图片和动画文件。在制作过程中，考虑到渲染时间 较长，故应该将整个动画分为几块来进行渲染。由于渲染过程耗时较长，制作过程中应首先 应用低分辨率，在完成整个部分的敷色及调整后再应用高分辨率渲染最终效果；在渲染的过 程中应同时把背景图片作为渲染场景，这对调整视角及着色效果的微调很有帮助。