一种基于二维sph的溃坝洪水演进模拟方法
【技术领域】
[0001] 本发明属于地学过程模拟与地理信息系统技术领域，涉及一种仿真建模方法，特 别涉及一种基于二维SPH的溃坝洪水演进模拟方法。
【背景技术】
[0002] 频发的地质灾害给社会带来了巨大损失，对于诸如溃坝洪水等自然或人为灾难的 过程仿真与模拟，以其涉及地学、仿真计算、物理学、计算机图形学等多个学科领域，已成为 目前研宄的热点和难点。
[0003] 在地理信息系统领域，传统的仿真模拟是利用地学模型或者数学模型等经验公式 模型来实现，即通过分析溃坝洪水演进过程的主要影响因子，例如地形坡度、汇流量等，利 用空间分析方法（如汇流分析、元胞自动机分析、种子淹没法分析等）建立地学模型。使用 这类方法实现的模拟速度较快，适用于大规模的快速地进行洪水过程模拟。但同时因为忽 略了很多物理学与动力学因素，其模拟的可靠性受到模型的参数选取与空间分析的约束， 因此其计算结果的精度往往较低，仿真效果简单，在某些应用中不能满足高精度高仿真的 需求。
[0004] 目前，主流的仿真模拟是使用基于物理模型的方法，以计算流体动力学中著名的 Navier-St〇ke S(NS)方程为理论基础，包括基于欧拉网格的模拟和基于拉格朗日粒子的模 拟。但是，基于网格的方法在处理大变形的流体时网格易发生畸变而造成计算复杂度的上 升以及精度的下降，而以光滑粒子流体动力学（Smoothed Particle Hydrodynamics，SPH) 为代表的基于拉格朗日粒子法在求解控制方程时，往往追求模拟效果的真实性而增加粒子 数量，并且需要进行流体表面的重构（如Marching Cubes方法），使得计算量剧增，增加了 对计算机硬件的要求。在一定程度上，这些方法都不满足例如溃坝洪水这种大范围场景的 流体实时仿真的需求。
【发明内容】
[0005] 为了同时满足诸如溃坝洪水这种大范围场景的流体仿真对于较高精度以及较低 计算复杂度的需求，本发明提出一种基于二维SPH的溃坝洪水演进模拟方法。该方法摒弃 三维的NS方程，使用二维的浅水方程为流体的控制方程，利用光滑粒子流体动力学离散化 流体，建立地形地物与水粒子的受力模型和状态更新方法，并提供快捷的边界检测方法，以 及使用双线性差值的方法重构流体表面进行渲染。
[0006] 本发明提出的一种基于二维SPH的溃坝洪水演进模拟方法，包括如下几个步骤：
[0007] -、获得水库及下游地区的地理空间数据信息，主要包括水库水位、水库库区以及 下游地区的地形高程数据。
[0008] 二、根据步骤一中获得的数据，建立基于浅水方程的洪水演进二维SPH模型，并对 仿真过程进行初始化，包括粒子属性的初始配置，建立粒子的邻域粒子列表以及参数设置。 首先根据步骤一中获取的水库范围以及水位信息，将水体粒子化，即将水体用粒子的形式 表达，同时完成参数设置以及粒子属性的初始化：
[0009] 1)根据水库范围以及设定的水粒子间距，计算得到水粒子的数量；
[0010] 2)建立水粒子集合的位置数组和速度数组；
[0011] 3)初始化水粒子的位置坐标，并将粒子的速度设置为零，粒子的高度设置为水库 的初始水位，计算得到每个粒子的体积，初始化粒子的光滑核半径；
[0012] 4)设置光滑核函数、粘滞系数以及阻尼系数等参数，并定义时间步长。
[0013] 对粒子进行初始化配置后，为每个粒子开辟邻域粒子列表存储空间，即为了对每 一个粒子寻找其支持域内的所有粒子，将支持域内的邻近粒子的粒子索引值存入邻近粒子 列表中。
[0014] 三、进入仿真循环步骤。首先根据粒子的位置设置虚粒子，其中，虚粒子的设置是 为了提高边界处粒子由于邻近粒子的数量不足造成的计算精度偏差；
[0015] 四、计算粒子的光滑核长度，并更新每个粒子的邻近粒子列表。由于此步骤在仿真 循环中最为耗时，因此本发明也涉及一种基于动态网格的近邻粒子搜索方法，该方法能够 降低更新邻近粒子列表时的计算复杂度。具体步骤为：
[0016] 1)根据粒子的二维坐标值，将所有粒子分别按照x坐标和y坐标进行升序排序；
[0017] 2)对待搜索邻近粒子的粒子，先在x坐标的正反两个方向上搜索距该粒子为光滑 核半径范围内的粒子，并标记；
[0018] 3)再在y坐标的正反两个方向上距粒子为光滑核半径范围内的粒子进行遍历， 如被标记，则计算粒子间距，如小于粒子的光滑核长度，则将该粒子存入邻近粒子索引列表 中，同时将此距离存入邻近粒子距离列表中。
[0019] 五、计算每个粒子所受的力，具体包括：
[0020] 1)计算粒子所受的压力c
[0021] 2)计算粒子所受的地形梯度力
[0022] 3)计算粒子所受的人工粘滞力f/、
[0023]4)计算粒子所受的合力Fi= f f+f^^n+fV^o
[0024] 六、计算每个水粒子在下一个时间步的位置和速度。具体包括：
[0025] 1)计算粒子在合力Fi下的加速度a i;
[0026] 2)计算粒子的速度增量AVi=aiA t;
[0027] 3)更新粒子的下一刻速度=v, +Av;;
[0028] 4)更新粒子的下一刻位置工广=x,. +v广Af。
[0029] 七、进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，否则进入 下一步。
[0030] 八、根据粒子分布提取流体表面。由于粒子携带高度信息，因此步骤可归纳如下：
[0031]1)根据粒子坐标范围，建立一个规则网格高程模型；
[0032] 2)将每个粒子按照坐标值映射到该网格中，并在网格中存储粒子高程值，如果网 格被多个粒子映射，则此网格的高度值取所有粒子高度的平均值；
[0033] 3)对于没有被映射的网格，进行双线性插值以补全数据，从而得到表示流体表面 的网格数据。
[0034] 九、对提取的流体表面进行植染：采用OpenSceneGraph三维植染引擎进行图形植 染，将步骤八中提取的网格数据传输到显卡中进行水面的绘制和显示，并添加光照、阴影等 增强真实感的效果。
[0035] 十、判断仿真是否结束。如果仿真时间到达设定的时长或者人为关闭仿真进程，则 停止运行，否则重复执行步骤三至步骤九。
[0036] 本发明的一种基于二维SPH的溃坝洪水演进模拟方法，能够针对溃坝洪水这一大 范围的流体行为，利用现有的地理空间数据信息，基于浅水方程的二维SPH求解，进行较高 精度的仿真模拟，并提供较为真实逼真的模拟结果，弥补了利用现有的地学模型与空间分 析算法在溃坝洪水演进过程模拟上的精度不足，以及基于求解NS方程造成的计算量巨大 的缺陷，为溃坝洪水演进模拟的进一步分析与应用提供有效的技术手段。
【附图说明】
[0037] 图1为本发明实施中一种基于二维SPH的溃坝洪水演进模拟方法的总体流程图。
[0038] 图2为使用二维SPH方法离散空间水粒子的示意图。
[0039] 其中，水粒子i在水平平面内运动，Ui为粒子的速度矢量，h 粒子的水深属性值， 粒子所在位置的地形高度为H(x，y)。在本方法中，粒子可以直观地理解为具有一定高度的 水柱。
[0040] 图3为固壁边界处虚粒子设置及边界惩罚力示意图。
[0041] 其中，1为光滑核长度，d为粒子距边界的距离。
[0042] 图4为水粒子的存储结构示意图。
[0043] 其中，N表示水粒子的数量，M为每个时间步内虚粒子的实际数量，且有M〈N。
[0044] 图5为基于动态网格的邻近粒子搜索方法示意图。
[0045] 其中，？1表示待查找邻域的粒子，其坐标为（x i，yi)，1为光滑核长度。
【具体实施方式】
[0046] 为更加清晰明了地展示本发明的目的、技术方案以及优点，下面结合附图详细叙 述本发明的具体实施例。
[0047] 如图1，为本发明实施中一种基于二维SPH的溃坝洪水模拟方法的程序流程图。具 体步骤如下：
[0048] 一、获得水库以及下游地理空间数据信息。
[0049] 首先确定待模拟的水库，在谷歌地球上确定经炜度坐标，根据下游地形情况确定 下游可能被淹没的范围。并使用多边形工具选取研宄区域（此区域需要把水库以及可能被 淹没的下游区域包括在内）。之后根据此区域的经炜度范围在免费获得的数字高程图上截 取待研宄范围的高程数据。
[0050] 同时，查阅水库信息确定水库正常的水位高度。
[0051] 二、建立基于浅水方程的洪水演进二维模型，并完成一系列初始化过程。
[0052] 1)浅水方程通过2D高度场的变化来描述流体表面行为，其Lagrangian形式的方 程为：
[0055] 其中，h表示地面上水体的高度，u = (u，v)表示流体的水平速度，g表示重力加速 度，H表示地形高程，V:表示梯度算子。
[0056] 使用SPH方法离散化上述两个方程，并根据Lee等（Lee H, Han S. Solving the Shallow Water equations using 2D SPH particles for interactive applications[J]. The Visual Computer, 2010, 26 (6-8): 865-872.)所述的方法得到
[0057]
[0058]
[0059] 同时，为提高计算精度，引入人工粘滞项。另外，为提高边界处的粒子行为真实性， 在边界处使用惩罚力与虚粒子结合的方法，如此，速度更新公式修正为：
) 浪准常数％ = 0. 1，0 n= 0. 01，因子供=G_%_，c和v分别表示粒子的声速以及速度矢量。
[0064]
[0065]
[0066]
[0067]