[0068] 2)初始化粒子的属性，以及SPH的相关参数。
[0069] 在使用二维SPH方法离散浅水方程时，水粒子将表示为带有一定高度的水柱，粒 子在水平面内运动，因此其速度为一个二维矢量，这些粒子的特点在附图2中已经展现。为 此，每个粒子所具备的属性有：位置（二维矢量）、速度（二维矢量）、水高、所在地形处的地 形高、光滑核长度、所受合力（二维矢量）、粒子索引值。
[0070] 初始化时，首先根据步骤一中的水库范围以及水位高度设置水粒子。根据指定的 粒子间距，设置粒子的位置、水高，最终得到粒子的总数目。之后，遍历粒子，将粒子的速度、 合力设为零值，根据粒子的位置检索得到处于地形处的高度，将此数值设置为粒子所在地 形处的地形高，并设置光滑核长度，并按照遍历顺序设置粒子的索引值。
[0071] 3)初始化粒子邻近粒子列表
[0072] 为每个粒子开辟一个链表空间，以存储邻近粒子索引。
[0073] 三、进入仿真循环步骤中。先根据粒子的位置判断是否需要设置边界粒子。附图 3中给出了虚粒子设置的方案：对每个粒子，判断其距固壁边界的距离，如果小于光滑核长 度的二倍，即在边界外的对称处设置一虚粒子，该虚粒子的速度与实粒子相反，其他属性 相同。
[0074] 虚粒子存储在实粒子之后，如附图4所示。首先在步骤二中，确定了实粒子的数 量，在此后继续开辟与实粒子数量相同的空间以存储虚粒子，并在每个循环步中，通过上述 的边界条件，可确定虚粒子的数目M。因此实际上，每次循环中M的大小是在变化的，其上限 为N(因为一个实粒子只能存在一个虚粒子）。
[0075] 四、计算粒子的光滑核长度，并更新每个水粒子的邻域粒子列表。
[0076] 对于大变形问题，为保证全局精度的一致性，需对光滑核长度进行修正，根据粒子 自身的密度更新其独自的光滑长度，以保证邻近粒子数量的相对稳定。计算方法为：
[0077]
[0078] 其中，L为初始时刻的光滑长度，h ^分别表示初始时刻的粒子水深。
[0079] 在每次仿真循环中，每个粒子的光滑核长都会按照公式（9)计算更新，因此，需要 动态地调整邻域粒子的搜索范围。本发明涉及的一种基于动态网格的邻域粒子搜索方法如 下：
[0080] 对于N个粒子Pi(Xi, yj, i = 1，2,…，N，首先分别根据粒子的X，y坐标值进行排 序，并分别查找两个坐标方向上与粒子距离为1的粒子，则在两个坐标方向上都被查找到 的粒子组成了邻近粒子备选集合，然后在此集合内进行进一步查找即可得到邻域粒子。
[0081] 如图5为邻近粒子搜索示意图。粒子的邻域搜索过程表述如下：
[0082] 1)根据N个粒子的坐标值，分别按照x，y坐标值大小进行升序排序；
[0083]2)针对每一个待查找邻域的粒子Pi(Xi，yi)，首先沿x坐标的正反两个方向搜索区 间的粒子，并进行标记；
[0084] 3)然后沿y坐标的正反两个方向，搜索区间[yi_l，yi+l]内的粒子，如果粒子已被 标记，则将粒子索引存入邻近粒子备选集合；
[0085] 4)在邻近粒子备选集合内计算每个粒子与粒子Pi(Xi，yi)的距离，如小于1，则将 该粒子索引与距离值存入粒子Pi(Xi，yj的邻近粒子索引链表和邻近粒子距离链表中。
[0086] 在每次循环中都要重新搜索每个粒子邻域粒子，因此图2中以阴影标示的网格大 小是动态改变的。该方法中，两次排序的时间复杂度为〇 (NlgN)(采用快速排序实现），遍历 搜索区间粒子的复杂度为〇(M) (M为搜索区间内粒子的数量，且N>M)，因此总的时间复杂度 为O(NlgN)，此方法操作简洁，空间复杂度较小。
[0087] 五、计算每个水粒子所受的力。
[0091] 4)计算粒子所受的合力Fi= f f+fitewn+f,'
[0092] 六、计算每个水粒子在下一个时间步的位置和速度。具体包括：
[0093] 1)计算粒子在合力Fi下的加速度a i;
[0094] 2)计算粒子的速度增量Avi= a^t ;
[0095] 3)更新粒子的下一刻速度= v,+ + ;
[0096] 4)更新粒子的下一刻位置= x,. + 。
[0097] 七、进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，否则进入 下一步。
[0098] 根据图3所示的方法进行边界判断，如遇到边界，施加一边界力，力的计算方法如 公式（7)，并根据粒子现在所处的位置以及速度进行修正。即分别对于位置和速度施加修正 量，已达到边界条件的处理。
[0099] 八、根据粒子分布提取流体表面。由于粒子携带高度信息，因此步骤可归纳如 下：
[0100] 1)根据粒子坐标范围，建立一个规则网格高程模型；
[0101] 2)将每个粒子按照坐标映射到该网格中，并存储粒子高程值；
[0102] 3)对于没有被映射的网格，进行双线性插值补全数据，从而得到表示流体表面的 网格数据。
[0103] 流体是用存储了高度信息的粒子来表示的，因此可以将所有粒子的坐标值以及高 度值视为高程图，这样方便了渲染，同时也避免了在三维情况下需要进行流体表面重构的 复杂操作，进一步简化了计算复杂度。
[0104] 九、对提取的流体表面进行植染：采用OpenSceneGraph三维植染引擎进行图形植 染，将步骤八中提取的网格数据传输到显卡中进行水面的绘制和显示。
[0105] 十、判断仿真是否结束。如仿真时间到达设定的时长或者人为关闭仿真进程，则停 止运行，否则重复执行步骤三至步骤九。
[0106] 使用本发明涉及的方法，对一真实水库进行溃坝洪水的演进仿真，得到特定时刻 的洪水演进效果。
[0107] 以上所述是本发明的优选实施方式，应当指出，对于本技术领域的相关人员，在不 脱离本发明原理的前提下，还可以进行若干改进，但这些改进也应视为本发明的保护范围。
【主权项】
1. 一种基于二维SPH的溃坝洪水演进模拟方法，其特征在于包括以下步骤： (一） 获得水库以及下游地理空间数据信息：首先确定待模拟的水库，确定经炜度坐标 范围，之后根据此经炜度范围在获得的数字高程图上截取待研宄范围的高程数据； (二） 建立基于浅水方程的洪水演进二维模型，并完成一系列初始化过程： 1) 通过二维浅水方程建立溃坝洪水演进过程的物理模型，并添加边界处理条件； 2) 初始化SPH的相关参数以及粒子的属性，包括位置、速度、水高、所在地形处的地形 高、光滑核长度、所受合力、粒子索引值，根据水库范围设置水粒子的位置和数目，根据水位 高度设置水高，并将粒子的速度、合力设为零值； 3) 初始化粒子邻近粒子列表，为每个粒子开辟一个链表空间，以存储邻近粒子索引； (三） 进入仿真循环步骤，根据粒子位置判断是否设置边界虚粒子，根据粒子距边界的 间距判断是否设置其镜像虚粒子，新设置的虚粒子存储在实粒子空间之后； (四） 根据粒子的水深，计算修正粒子自身的光滑核长度，并进行邻近粒子搜索算法更 新每个水粒子的邻近粒子列表； (五） 计算每个水粒子所受的力； 1) 计算粒子所受的压力2) 计算粒子所受的地形梯度力=-公▽"; 3) 计算粒子所受的人工粘滞力4) 计算粒子所受的合力巧=+_/T_+_/T; (六） 计算每个水粒子在下一个时间步的位置和速度；具体包括： 1) 计算粒子在合力Fi下的加速度ai; 2) 计算粒子的速度增量AVi=aiAt; 3) 更新粒子的下一刻速度'^广=力+夂； 4) 更新粒子的下一刻位置=A+v广; (七） 进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度； (八） 根据粒子分布提取流体表面； (九） 对提取的流体表面重构：采用OpenSceneGraph三维植染引擎进行图形植染； (十）判断仿真是否结束。2. 根据权利要求1所述的溃坝洪水演进模拟方法，其特征在于，所述的邻近粒子搜索 方法完成粒子的搜索，具体步骤为： 1) 根据粒子的坐标值，分别按坐标值x，y大小进行升序的计数排序； 2) 针对每一个待查找邻域的粒子Pi(Xi，yi)，首先沿x坐标的正反两个方向搜索区间 [Xi-hXi+l]的粒子，并进行标记； 3) 然后沿y坐标的正反两个方向，搜索区间[yi_l，yi+l]内的粒子，如果粒子已被标记， 则将粒子索引存入邻近粒子备选集合； 4) 在邻近粒子备选集合内计算每个粒子与粒子Pi(Xi，yi)的距离，如小于1，则将该粒 子索引与距离值存入粒子的邻近粒子索引链表和邻近粒子距离链表中。3.如权利要求1或2所述的溃坝洪水演进模拟方法，其特征在于，所述的流体表面重构 过程步骤为： 1) 根据粒子坐标范围，建立一个规则网格高程模型； 2) 将每个粒子按照坐标映射到该网格中，并存储粒子高程值； 3) 对于没有被映射的网格，进行双线性插值补全数据，从而得到表示流体表面的网格 数据。
【专利摘要】一种基于二维SPH的溃坝洪水演进模拟方法，属于地学过程模拟与地理信息系统技术领域。该方法先获得水库及下游地理空间数据信息，再根据获得的数据，建立基于浅水方程的洪水演进二维模型，对仿真过程进行初始化，包括粒子属性的初始配置，建立粒子的邻域粒子列表以及参数设置；然后，进入仿真循环，根据粒子位置判断是否设置边界虚粒子，更新每个水粒子的邻域粒子列表；计算每个水粒子所受的力后，计算每个水粒子在下一个时间步的位置和速度；进行碰撞检测，如遇边界，施加边界惩罚力以改变粒子位置与速度，根据粒子分布提取流体表面，对提取的流体表面进行渲染。本发明所述方法可以对较大规模的溃坝洪水演进情况进行比较真实的仿真模拟。
【IPC分类】G06F17/50
【公开号】CN104991999
【申请号】CN201510337225
【发明人】韩敏, 边茂松
【申请人】大连理工大学
【公开日】2015年10月21日
【申请日】2015年6月17日