CN104991999A - 一种基于二维sph的溃坝洪水演进模拟方法 - Google Patents

Abstract

一种基于二维SPH的溃坝洪水演进模拟方法，属于地学过程模拟与地理信息系统技术领域。该方法先获得水库及下游地理空间数据信息，再根据获得的数据，建立基于浅水方程的洪水演进二维模型，对仿真过程进行初始化，包括粒子属性的初始配置，建立粒子的邻域粒子列表以及参数设置；然后，进入仿真循环，根据粒子位置判断是否设置边界虚粒子，更新每个水粒子的邻域粒子列表；计算每个水粒子所受的力后，计算每个水粒子在下一个时间步的位置和速度；进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，根据粒子分布提取流体表面，对提取的流体表面进行渲染。本发明所述方法可以对较大规模的溃坝洪水演进情况进行比较真实的仿真模拟。

Description

一种基于二维SPH的溃坝洪水演进模拟方法

技术领域

本发明属于地学过程模拟与地理信息系统技术领域，涉及一种仿真建模方法，特别涉及一种基于二维SPH的溃坝洪水演进模拟方法。

背景技术

频发的地质灾害给社会带来了巨大损失，对于诸如溃坝洪水等自然或人为灾难的过程仿真与模拟，以其涉及地学、仿真计算、物理学、计算机图形学等多个学科领域，已成为目前研究的热点和难点。

在地理信息系统领域，传统的仿真模拟是利用地学模型或者数学模型等经验公式模型来实现，即通过分析溃坝洪水演进过程的主要影响因子，例如地形坡度、汇流量等，利用空间分析方法(如汇流分析、元胞自动机分析、种子淹没法分析等)建立地学模型。使用这类方法实现的模拟速度较快，适用于大规模的快速地进行洪水过程模拟。但同时因为忽略了很多物理学与动力学因素，其模拟的可靠性受到模型的参数选取与空间分析的约束，因此其计算结果的精度往往较低，仿真效果简单，在某些应用中不能满足高精度高仿真的需求。

目前，主流的仿真模拟是使用基于物理模型的方法，以计算流体动力学中著名的Navier-Stokes(NS)方程为理论基础，包括基于欧拉网格的模拟和基于拉格朗日粒子的模拟。但是，基于网格的方法在处理大变形的流体时网格易发生畸变而造成计算复杂度的上升以及精度的下降，而以光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)为代表的基于拉格朗日粒子法在求解控制方程时，往往追求模拟效果的真实性而增加粒子数量，并且需要进行流体表面的重构(如Marching Cubes方法)，使得计算量剧增，增加了对计算机硬件的要求。在一定程度上，这些方法都不满足例如溃坝洪水这种大范围场景的流 体实时仿真的需求。

发明内容

为了同时满足诸如溃坝洪水这种大范围场景的流体仿真对于较高精度以及较低计算复杂度的需求，本发明提出一种基于二维SPH的溃坝洪水演进模拟方法。该方法摒弃三维的NS方程，使用二维的浅水方程为流体的控制方程，利用光滑粒子流体动力学离散化流体，建立地形地物与水粒子的受力模型和状态更新方法，并提供快捷的边界检测方法，以及使用双线性差值的方法重构流体表面进行渲染。

本发明提出的一种基于二维SPH的溃坝洪水演进模拟方法，包括如下几个步骤：

一、获得水库及下游地区的地理空间数据信息，主要包括水库水位、水库库区以及下游地区的地形高程数据。

二、根据步骤一中获得的数据，建立基于浅水方程的洪水演进二维SPH模型，并对仿真过程进行初始化，包括粒子属性的初始配置，建立粒子的邻域粒子列表以及参数设置。首先根据步骤一中获取的水库范围以及水位信息，将水体粒子化，即将水体用粒子的形式表达，同时完成参数设置以及粒子属性的初始化：

1)根据水库范围以及设定的水粒子间距，计算得到水粒子的数量；

2)建立水粒子集合的位置数组和速度数组；

3)初始化水粒子的位置坐标，并将粒子的速度设置为零，粒子的高度设置为水库的初始水位，计算得到每个粒子的体积，初始化粒子的光滑核半径；

4)设置光滑核函数、粘滞系数以及阻尼系数等参数，并定义时间步长。

对粒子进行初始化配置后，为每个粒子开辟邻域粒子列表存储空间，即为了对每一个粒子寻找其支持域内的所有粒子，将支持域内的邻近粒子的粒子索引值存入邻近粒子列表中。

三、进入仿真循环步骤。首先根据粒子的位置设置虚粒子，其中，虚粒子的设置是为了提高边界处粒子由于邻近粒子的数量不足造成的计算精度偏差；

四、计算粒子的光滑核长度，并更新每个粒子的邻近粒子列表。由于此步骤在仿真循环中最为耗时，因此本发明也涉及一种基于动态网格的近邻粒子搜索方法，该方法能够降低更新邻近粒子列表时的计算复杂度。具体步骤为：

1)根据粒子的二维坐标值，将所有粒子分别按照x坐标和y坐标进行升序排序；

2)对待搜索邻近粒子的粒子，先在x坐标的正反两个方向上搜索距该粒子为光滑核半径范围内的粒子，并标记；

3)再在y坐标的正反两个方向上距粒子为光滑核半径范围内的粒子进行遍历，如被标记，则计算粒子间距，如小于粒子的光滑核长度，则将该粒子存入邻近粒子索引列表中，同时将此距离存入邻近粒子距离列表中。

五、计算每个粒子所受的力，具体包括：

1)计算粒子所受的压力f_i ^press；

2)计算粒子所受的地形梯度力f_i ^terrain；

3)计算粒子所受的人工粘滞力f_i ^visc；

4)计算粒子所受的合力F_i＝f_i ^press+f_i ^terrain+f_i ^visc。

六、计算每个水粒子在下一个时间步的位置和速度。具体包括：

1)计算粒子在合力F_i下的加速度a_i；

2)计算粒子的速度增量Δv_i＝a_iΔt；

3)更新粒子的下一刻速度

4)更新粒子的下一刻位置

七、进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，否则进入下一步。

八、根据粒子分布提取流体表面。由于粒子携带高度信息，因此步骤可归纳如下：

1)根据粒子坐标范围，建立一个规则网格高程模型；

2)将每个粒子按照坐标值映射到该网格中，并在网格中存储粒子高程值，如果网格被多个粒子映射，则此网格的高度值取所有粒子高度的平均值；

3)对于没有被映射的网格，进行双线性插值以补全数据，从而得到表示流体表面的网格数据。

九、对提取的流体表面进行渲染：采用OpenSceneGraph三维渲染引擎进行图形渲染，将步骤八中提取的网格数据传输到显卡中进行水面的绘制和显示，并添加光照、阴影等增强真实感的效果。

十、判断仿真是否结束。如果仿真时间到达设定的时长或者人为关闭仿真进程，则停止运行，否则重复执行步骤三至步骤九。

本发明的一种基于二维SPH的溃坝洪水演进模拟方法，能够针对溃坝洪水这一大范围的流体行为，利用现有的地理空间数据信息，基于浅水方程的二维SPH求解，进行较高精度的仿真模拟，并提供较为真实逼真的模拟结果，弥补了利用现有的地学模型与空间分析算法在溃坝洪水演进过程模拟上的精度不足，以及基于求解NS方程造成的计算量巨大的缺陷，为溃坝洪水演进模拟的进一步分析与应用提供有效的技术手段。

附图说明

图1为本发明实施中一种基于二维SPH的溃坝洪水演进模拟方法的总体流程图。

图2为使用二维SPH方法离散空间水粒子的示意图。

其中，水粒子i在水平平面内运动，u_i为粒子的速度矢量，h_i为粒子的水深属性值，粒子所在位置的地形高度为H(x,y)。在本方法中，粒子可以直观地理解为具有一定高度的水柱。

图3为固壁边界处虚粒子设置及边界惩罚力示意图。

其中，l为光滑核长度，d为粒子距边界的距离。

图4为水粒子的存储结构示意图。

其中，N表示水粒子的数量，M为每个时间步内虚粒子的实际数量，且有M<N。

图5为基于动态网格的邻近粒子搜索方法示意图。

其中，P_i表示待查找邻域的粒子，其坐标为(x_i,y_i),l为光滑核长度。

具体实施方式

为更加清晰明了地展示本发明的目的、技术方案以及优点，下面结合附图详细叙述本发明的具体实施例。

如图1，为本发明实施中一种基于二维SPH的溃坝洪水模拟方法的程序流程图。具体步骤如下：

一、获得水库以及下游地理空间数据信息。

首先确定待模拟的水库，在谷歌地球上确定经纬度坐标，根据下游地形情况确定下游可能被淹没的范围。并使用多边形工具选取研究区域(此区域需要把水库以及可能被淹没的下游区域包括在内)。之后根据此区域的经纬度范围在 免费获得的数字高程图上截取待研究范围的高程数据。

同时，查阅水库信息确定水库正常的水位高度。

二、建立基于浅水方程的洪水演进二维模型，并完成一系列初始化过程。

1)浅水方程通过2D高度场的变化来描述流体表面行为，其Lagrangian形式的方程为：

$\frac{\mathrm{Dh}}{\mathrm{Dt}}=-h\&dtri;\&CenterDot;u---\left(1\right)$

$\frac{\mathrm{Du}}{\mathrm{Dt}}=-g\&dtri;(h+H)---\left(2\right)$

其中，h表示地面上水体的高度，u＝(u,v)表示流体的水平速度，g表示重力加速度，H表示地形高程，表示梯度算子。

使用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.)所述的方法得到

$h_i=\underset{j}{\overset{N}{\mathrm{\&Sigma;}}}{V}_{j}W(r_i-r_j.l)---\left(3\right)$

$\frac{{\mathrm{Du}}_i}{\mathrm{Dt}}=-g\underset{j}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j\&dtri;W(r_i-r_j,l)-g\&dtri;H---\left(4\right)$

同时，为提高计算精度，引入人工粘滞项。另外，为提高边界处的粒子行为真实性，在边界处使用惩罚力与虚粒子结合的方法，如此，速度更新公式修正为：

$\frac{{\mathrm{Du}}_i}{\mathrm{Dt}}=-g\underset{j}{\mathrm{\&Sigma;}}{V}_j\&dtri;W(r_i-r_j,l)-g\&dtri;H+\underset{j}{\mathrm{\&Sigma;}}{V}_j{\mathrm{\&Pi;}}_{\mathrm{ij}}\&dtri;W(r_i-r_j,l)+{f_i}^{\mathrm{pb}}---\left(5\right)$

其中，

${\mathrm{\&Pi;}}_{\mathrm{ij}}=\left\{\begin{array}{cc}\frac{-{\mathrm{\&alpha;}}_{\mathrm{\&Pi;}}{\stackrel{\&OverBar;}{c}}_{\mathrm{ij}}{\mathrm{\&phi;}}_{\mathrm{ij}}+{\mathrm{\&beta;}}_{\mathrm{\&Pi;}}{\mathrm{\&phi;}}_{\mathrm{ij}}^2}{h_{\mathrm{ij}}},& v_{\mathrm{ij}}\&CenterDot;r_{\mathrm{ij}}<0\\ 0,& v_{\mathrm{ij}}\&CenterDot;r_{\mathrm{ij}}\&GreaterEqual;0\end{array}\right.---\left(6\right)$

式中，标准常数α_Π＝0.1， β_Π＝0.01，因子c和v分别表示粒子的声速以及速度矢量。

$f^{\mathrm{pb}}=\left\{\begin{array}{cc}(k_{s}d-k_d(u\&CenterDot;n))n,& d\&le;d_s\\ 0,& d>d_s\end{array}\right.---\left(7\right)$

光滑核函数取：

$W(R,l)={\mathrm{\&alpha;}}_d\left\{\begin{array}{cc}(\frac{2}{3}-\frac{9}{8}{R}^2+\frac{19}{24}{R}^3-\frac{5}{32}{R}^4),& 0\&le;R\&le;2\\ 0,& R>2\end{array}\right.---\left(8\right)$

式中，α_d在二维空间中取值为

2)初始化粒子的属性，以及SPH的相关参数。 

在使用二维SPH方法离散浅水方程时，水粒子将表示为带有一定高度的水柱，粒子在水平面内运动，因此其速度为一个二维矢量，这些粒子的特点在附图2中已经展现。为此，每个粒子所具备的属性有：位置(二维矢量)、速度(二维矢量)、水高、所在地形处的地形高、光滑核长度、所受合力(二维矢量)、粒子索引值。

初始化时，首先根据步骤一中的水库范围以及水位高度设置水粒子。根据指定的粒子间距，设置粒子的位置、水高，最终得到粒子的总数目。之后，遍历粒子，将粒子的速度、合力设为零值，根据粒子的位置检索得到处于地形处的高度，将此数值设置为粒子所在地形处的地形高，并设置光滑核长度，并按照遍历顺序设置粒子的索引值。

3)初始化粒子邻近粒子列表

为每个粒子开辟一个链表空间，以存储邻近粒子索引。

三、进入仿真循环步骤中。先根据粒子的位置判断是否需要设置边界粒子。附图3中给出了虚粒子设置的方案：对每个粒子，判断其距固壁边界的距离，如果小于光滑核长度的二倍，即在边界外的对称处设置一虚粒子，该虚粒子的 速度与实粒子相反，其他属性相同。

虚粒子存储在实粒子之后，如附图4所示。首先在步骤二中，确定了实粒子的数量，在此后继续开辟与实粒子数量相同的空间以存储虚粒子，并在每个循环步中，通过上述的边界条件，可确定虚粒子的数目M。因此实际上，每次循环中M的大小是在变化的，其上限为N(因为一个实粒子只能存在一个虚粒子)。

四、计算粒子的光滑核长度，并更新每个水粒子的邻域粒子列表。

对于大变形问题，为保证全局精度的一致性，需对光滑核长度进行修正，根据粒子自身的密度更新其独自的光滑长度，以保证邻近粒子数量的相对稳定。计算方法为：

$l=l_0{\left(\frac{h_0}{h}\right)}^{\frac{1}{2}}---\left(9\right)$

其中，l₀为初始时刻的光滑长度，h₀分别表示初始时刻的粒子水深。

在每次仿真循环中，每个粒子的光滑核长都会按照公式(9)计算更新，因此，需要动态地调整邻域粒子的搜索范围。本发明涉及的一种基于动态网格的邻域粒子搜索方法如下：

对于N个粒子P_i(x_i,y_i),i＝1,2,…,N，首先分别根据粒子的x,y坐标值进行排序，并分别查找两个坐标方向上与粒子距离为l的粒子，则在两个坐标方向上都被查找到的粒子组成了邻近粒子备选集合，然后在此集合内进行进一步查找即可得到邻域粒子。

如图5为邻近粒子搜索示意图。粒子的邻域搜索过程表述如下：

1)根据N个粒子的坐标值，分别按照x,y坐标值大小进行升序排序；

2)针对每一个待查找邻域的粒子P_i(x_i,y_i)，首先沿x坐标的正反两个方向搜索区间[x_i-l,x_i+l]的粒子，并进行标记；

3)然后沿y坐标的正反两个方向，搜索区间[y_i-l,y_i+l]内的粒子，如果粒子已被标记，则将粒子索引存入邻近粒子备选集合；

4)在邻近粒子备选集合内计算每个粒子与粒子P_i(x_i,y_i)的距离，如小于l，则将该粒子索引与距离值存入粒子P_i(x_i,y_i)的邻近粒子索引链表和邻近粒子距离链表中。

在每次循环中都要重新搜索每个粒子邻域粒子，因此图2中以阴影标示的网格大小是动态改变的。该方法中，两次排序的时间复杂度为O(NlgN)(采用快速排序实现)，遍历搜索区间粒子的复杂度为O(M)(M为搜索区间内粒子的数量,且N>M)，因此总的时间复杂度为O(NlgN)，此方法操作简洁，空间复杂度较小。

五、计算每个水粒子所受的力。

1)计算粒子所受的压力 ${f_i}^{\mathrm{press}}=-g\underset{j}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j\&dtri;W(r_i-r_j,l);$

2)计算粒子所受的地形梯度力

3)计算粒子所受的人工粘滞力

4)计算粒子所受的合力F_i＝f_i ^press+f_i ^terrain+f_i ^visc。

六、计算每个水粒子在下一个时间步的位置和速度。具体包括：

1)计算粒子在合力F_i下的加速度a_i；

2)计算粒子的速度增量Δv_i＝a_iΔt；

3)更新粒子的下一刻速度

4)更新粒子的下一刻位置

七、进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，否则进入下一步。

根据图3所示的方法进行边界判断，如遇到边界，施加一边界力，力的计算方法如公式(7)，并根据粒子现在所处的位置以及速度进行修正。即分别对于位置和速度施加修正量，已达到边界条件的处理。

八、根据粒子分布提取流体表面。由于粒子携带高度信息，因此步骤可归 纳如下：

1)根据粒子坐标范围，建立一个规则网格高程模型；

2)将每个粒子按照坐标映射到该网格中，并存储粒子高程值；

3)对于没有被映射的网格，进行双线性插值补全数据，从而得到表示流体表面的网格数据。

流体是用存储了高度信息的粒子来表示的，因此可以将所有粒子的坐标值以及高度值视为高程图，这样方便了渲染，同时也避免了在三维情况下需要进行流体表面重构的复杂操作，进一步简化了计算复杂度。

九、对提取的流体表面进行渲染：采用OpenSceneGraph三维渲染引擎进行图形渲染，将步骤八中提取的网格数据传输到显卡中进行水面的绘制和显示。

十、判断仿真是否结束。如仿真时间到达设定的时长或者人为关闭仿真进程，则停止运行，否则重复执行步骤三至步骤九。

使用本发明涉及的方法，对一真实水库进行溃坝洪水的演进仿真，得到特定时刻的洪水演进效果。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的相关人员，在不脱离本发明原理的前提下，还可以进行若干改进，但这些改进也应视为本发明的保护范围。

Claims (3)

1.一种基于二维SPH的溃坝洪水演进模拟方法，其特征在于包括以下步骤：

(一)获得水库以及下游地理空间数据信息：首先确定待模拟的水库，确定经纬度坐标范围，之后根据此经纬度范围在获得的数字高程图上截取待研究范围的高程数据；

(二)建立基于浅水方程的洪水演进二维模型，并完成一系列初始化过程：

1)通过二维浅水方程建立溃坝洪水演进过程的物理模型，并添加边界处理条件；

2)初始化SPH的相关参数以及粒子的属性，包括位置、速度、水高、所在地形处的地形高、光滑核长度、所受合力、粒子索引值，根据水库范围设置水粒子的位置和数目，根据水位高度设置水高，并将粒子的速度、合力设为零值；

3)初始化粒子邻近粒子列表，为每个粒子开辟一个链表空间，以存储邻近粒子索引；

(三)进入仿真循环步骤，根据粒子位置判断是否设置边界虚粒子，根据粒子距边界的间距判断是否设置其镜像虚粒子，新设置的虚粒子存储在实粒子空间之后；

(四)根据粒子的水深，计算修正粒子自身的光滑核长度，并进行邻近粒子搜索算法更新每个水粒子的邻近粒子列表；

(五)计算每个水粒子所受的力；

1)计算粒子所受的压力 ${f_i}^{\mathrm{press}}=-g\underset{j}{\overset{N}{\mathrm{\&Sigma;}}}{V}_j\&dtri;W(r_i-r_j,l);$

2)计算粒子所受的地形梯度力

3)计算粒子所受的人工粘滞力

4)计算粒子所受的合力

(六)计算每个水粒子在下一个时间步的位置和速度；具体包括：

1)计算粒子在合力F_i下的加速度a_i；

2)计算粒子的速度增量Δv_i＝a_iΔt；

3)更新粒子的下一刻速度

4)更新粒子的下一刻位置

(七)进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度；

(八)根据粒子分布提取流体表面；

(九)对提取的流体表面重构：采用OpenSceneGraph三维渲染引擎进行图形渲染；

(十)判断仿真是否结束。

2.根据权利要求1所述的溃坝洪水演进模拟方法，其特征在于，所述的邻近粒子搜索方法完成粒子的搜索，具体步骤为：

1)根据粒子的坐标值，分别按坐标值x，y大小进行升序的计数排序；

2)针对每一个待查找邻域的粒子P_i(x_i,y_i)，首先沿x坐标的正反两个方向搜索区间[x_i-l,x_i+l]的粒子，并进行标记；

3)然后沿y坐标的正反两个方向，搜索区间[y_i-l,y_i+l]内的粒子，如果粒子已被标记，则将粒子索引存入邻近粒子备选集合；

4)在邻近粒子备选集合内计算每个粒子与粒子P_i(x_i,y_i)的距离，如小于l，则将该粒子索引与距离值存入粒子的邻近粒子索引链表和邻近粒子距离链表中。

3.如权利要求1或2所述的溃坝洪水演进模拟方法，其特征在于，所述的流体表面重构过程步骤为：

1)根据粒子坐标范围，建立一个规则网格高程模型；

2)将每个粒子按照坐标映射到该网格中，并存储粒子高程值；

3)对于没有被映射的网格，进行双线性插值补全数据，从而得到表示流体表面的网格数据。