摘 要：光滑粒子流体动力学（Smoothed Particle Hydrodynamics，简称SPH）方法已经广泛应用于流体计算中，但在模拟三维问题时常花费大量的时间，为了提高计算速度与精度，对SPH模拟三维问题程序进行并行化是非常必要的。在详细阐述SPH方法在并行化过程中的数据通信方案和搜索策略的基础上，对三维溃坝算例进行数值模拟和分析，合理再现了溃坝的过程，验证了并行方案的正确性。
关键词：SPH方法；溃坝模拟；并行化
中图分类号：TP301 文献标识码：A 文章编号：16727800（2012）011002402
________________________________________
作者简介：潘中建 （1986-），男，硕士，江苏省南通商贸高等职业学校管理系助教，研究方向为无网格数值模拟。0 引言
SPH方法是近30年来发展起来的一种纯Lagrange粒子方法。SPH方法的基本思想是利用核函数插值来计算空间各离散点上物理量的近似值，物理量对空间的导数则被转化为对核函数求导数，从而避免了高维Lagrange网格方法中的网格生成、缠结和扭曲等令人头痛的问题。SPH方法具有抗畸变能力强，能处理大变形和高速碰撞等问题，避免了大量的单元划分并克服了有限元方法中由于场函数近似引起的误差，能够计算复杂介质流动、自由表面流动等问题。
然而，SPH方法也遇到了一些缺陷，最严重的是其执行速度问题。在二维情况下，每个粒子有50个相邻粒子，而在三维模拟中，相邻粒子达到了200个，在整个模拟中，计算时间往往达数小时，因此对SPH代码的加速已经成为该方法发展的一部分。
本文在详细阐述并行化过程中数据通信的基础上，利用SPH方法对三维溃坝进行了数值模拟，合理再现了溃坝的过程。
1 SPH方法的控制方程
SPH是一种完全无网格的Lagrange方法，流体的流动以跟随流体的粒子运动来描述。在SPH方法中，某一粒子上的任一宏观变量都可以表示为计算区域内除此粒子之外所有粒子上该变量的积分插值。
描述流体运动的质量守恒方程和动量守恒方程的SPH形式为：
人工粘性项的表达式为：
2 并行方案
2.1 数据交换
本文利用区域分解方法将问题域分解成各个子区域，每个子区域被分配给并行系统的各个处理器，每个处理器只对分配到该处理器的子区域粒子进行计算。当一个粒子从一个处理器中的子区域移动到另外一个处理器中的子区域时，该粒子的相关属性也被转移到新的处理器，子区域之间是通过边界间的数据交换进行通信的。
将子区域分解成网格单元，这样对处于网格单元中的计算粒子与其附近粒子间的搜索就转变为网格单元之间的搜索。以二维情况举例（如图1所示），假定一个处理器处理该计算区域中的ncx×ncz网格单元，每个处理器能够对所有的内网格单元（图中实线内的白色区域）进行计算，但是在计算边界网格单元（图中灰色区域）中的粒子时，处理器还必须知道子区域边界另一边的网格单元中的粒子数据，这些数据必须从相邻处理器中取得并存储在图中的扩展网格单元中（图中虚线区域）。同时，边界网格单元中的数据还要发送到相邻子区域中的扩展网格单元中。因此在每步积分之前，相邻的处理器要通过数据交换进行通信（如图2所示）。
图1 区域分解
图2 局部边界通信
2.2 搜索策略
在SPH方法中，每个时间步都要进行相互作用粒子的搜索，这也是SPH程序在运行过程中主要的CPU时间。在程序中，用链表算法进行粒子间相互粒子作用的搜索。计算区域被分解成边长为2h的网格单元，对于给定的粒子，其相邻粒子只能在同一网格单元内，或者在紧密相邻的单元内。所以对于三维问题需要搜索27个网格单元。
但是当我们在进行计算时发现，如果搜索某个粒子的全部相邻网格单元，将会大大增加计算量，所以采用相互粒子作用对的方法（如图3所示）。以二维为例，对于网格单元中的粒子在进行相邻粒子搜索时，ik网格单元从计算背景网格的左下开始进行，那么只需要搜索E、N、NW和NE 4个方向的网格单元及自身一共5个网格单元，其它4个方向的网格单元通过粒子间的相互作用在前续的计算中已经被计算过。同样在三维的问题中，相邻粒子只需要搜索14个网格单元。
3 溃坝模拟
在三维溃坝模型中，初始时刻水体位于长为1.7m、宽度为0.72m、高度为0.4m的长方体容器固壁中，数值模拟中共使用25 057个粒子，其中布置固定边界粒子12 083个。图4分别为0.25s、0.6s、1s和1.2s 4个时刻流体的速度等值线图。从图中可以看出，当水体到达容器顶端时，水体与容器壁发生碰撞并升高，然后在重力作用下水体下落并形成卷波。
表1是在HP ProLiant DL140 G1/ DL140 G2服务器上运行得到的数值模拟的加速比和并行效率。从表1中可以看出并行效率得到了一定的提高，但随着处理器个数的增加，并行效率下降速度很快，这是因为模拟中粒子个数较少，因而导致细粒度并行。在计算规模不变的情况下，随着处理器个数增加，各处理器计算粒子的个数减少，处理器间的数据通信量增大，使并行效率不断降低。
4 结语
本文开发了基于SPH方法的并行程序，较为详细地阐述了在并行化过程中的数据通信方案和搜索策略，对三维溃坝算例进行了数值模拟和分析，模拟结果比较合理地展现了水体的流动过程，验证了并行方案的正确性。
参考文献：
[1] MONAGHAN J J.Simulating free surface flows with SPH[J].Journal of Computational Physics， 1994.
[2] TAKEDA H，MIYAMA SM，SEKIYA M.Numerical simulation of viscous flow by smoothed particle hydrodynamics[J].Prog Theor Physics ，1994.
[3] 张学莹，潘中建，马小飞.复杂边界溃坝流与结构物作用过程的数值模拟[J].水动力学研究与进展：A辑，2010（3）.
[4] 张林波，迟学斌，莫则尧，等.并行计算导论[M].北京：清华大学出版社，2006.
[5] MO ZEYAO， ZHANG JINGLIN.Parallelization and optimization for a two dimensional molecular dynamics code（MDP）[J].Chinese J Comput Phys，2000（1）.
（责任编辑：孙娟）