SPH方法在风沙流研究中的应用进展

赵杰，金阿芳，楚花明

(新疆大学机械工程学院,新疆 乌鲁木齐 830047)

光滑粒子流体动力学(Smoothed Particle Hydrodynamics)又被称为SPH方法，是在1977年最早由Gingold和Monaghan等人提出来并加以运用的[1,2]。该方法最初是应用在天体物理中以解决流体质团在无穷大三维空间中任意流动的问题。随着其他研究人员和学者的应用以及改进，SPH方法快速发展，并开始推广应用于流体力学方面的研究。SPH作为一种数值模拟方法同时具有无网格性质和粒子性质，特别适合于处理大变形问题、多相流问题以及追踪单个粒子运动[3]，在溃坝、水下爆炸、飞沫方面的问题研究中都取得了成果。

另一方面，为了研究风沙流的运动机理，防治风沙灾害，对风沙运动的研究在20世纪30年代初便引起了当时西方学者的关注。英国学者拜格诺(Bagnold)在1941发表了自己第一篇关于风沙地貌学的著作，标志着风沙物理学的开端[4]。此后，关于风沙运动的研究进一步发展，逐渐发展到现有的3种研究方法：野外观测、风洞实验和数值模拟。而对风沙流的研究也着眼于宏观和微观两个方面：宏观上注重对风沙流整体结构的研究，微观上则对单颗沙粒的运动规律进行研究，其中，跃移沙粒占整个输移沙粒总量的75%，成为连接宏观与微观研究之间的桥梁[5]。在现有的研究中，尚未对风沙运动的内在机理形成统一认识，对风沙流的定量描述也多是基于观测数据的半经验公式。所以，使用SPH方法对风沙运动进行数值模拟方法将有助于研究风沙流的运动规律并帮助完善其基本理论，同时也拓宽了SPH方法的应用领域，促进方法本身的成熟。

1 SPH模拟风沙流

风沙流是一种典型的气固两相流，其中包括固体和气体两种相互作用的物质。SPH方法作为一种流体拟颗粒模型，在将沙粒相作为离散相来处理同时，将气体也离散为“颗粒”性质的流体微团，通过两种不同物质颗粒之间的相互作用，来实现风沙流物理特性的模拟[6]。由于在模拟过程中不需要建立网格，同时又能携带单个粒子的信息，所以特别适合处理风沙流这种大变形的物理运动。

金阿芳等人从气固两相流的力学观点出发，通过沙粒相和气体相之间的碰撞建立了两相之间的耦合关系，建立了基于光滑粒子流体动力学方法的风沙运动控制方程[5]。在他们的工作中，编程实现了SPH算法对风沙流运动的模拟，提出了使用SPH方法对风沙流模拟时的一些关键参数的选取规则，在处理边界条件时，他们使用虚粒子来阻止粒子的穿透，为今后在对风沙流动进行数值模拟时，确定边界条件给出了一定的参考。同时，他们通过仿真模拟了跃移沙粒运动状态，确定了其运动参数和几何参数，其中运动参数包括起跳速度、沉降速度、碰撞速度等。几何参数包括跃移长度和高度、起跳角和碰撞角等，并对计算结果进行了分析。 牟阳阳等人使用SPH方法建立了倾斜沙床面下风沙运动的物理模型，他们模拟分析了在不同时间段沙粒在空间中的分布特点；统计了在沙床斜面上不同位置沙粒的起跳高度和距离，得出沙粒的平均跃移距离是单个沙波纹长度的1.16倍[7]。在逯博等人的研究中，对气体相使用N-S方程，对沙粒相使用牛顿定律进行处理，建立了两相的控制方程，他们使用SPH方法对沙粒的跃移进行了模拟，得到了沙粒的典型跃移轨迹，并在模拟中加入了风的加载和卸载模块，统计了在风的加载和卸载情况下沙粒的轨迹和法向速度的变化[8]。为了对沙粒的跃移参数进行总体性地描述，金阿芳等人在模拟中随机选取了1 000个沙粒进行分析，得到了它们的起跳角和起跳速度的概率密度分布。他们发现沙粒的跃移起跳角度的概率密度函数符合指数分布，起跳速度的概率密度函数符合对数正态分布；并在模拟中，他们对沙粒相和气体相分别选取了不同的光滑长度，得到了更好的耦合效果[9]。

陈福振等人则将SPH方法和传统的有限体积法(FVM)结合起来，建立了SPH-FVM耦合方法来模拟风沙运动[10]。在他们的研究中，对连续气体相采用FVM，将沙粒相视为拟流体采用SPH方法处理，气体相和沙粒相之间通过拖曳力、压力梯度和浓度等参数进行耦合。同时，他们对SPH方法进行了改造，SPH粒子代表颗粒群，包含沙粒群诸如密度、质量、体积、粒径分布和体积分数等物理信息，由此减小了问题规模，提高了计算效率。陈福振等人在SPH边界条件的设置上使用了罚函数方法施加边界力，对气体相使用无滑移边界条件。在对风沙流的模拟中，他们得到了沙粒的运动轨迹、沙粒在气体中的水平速度分布以及气体在沙粒的反作用下速度的变化，验证了方法的有效性。

2 模型和参数选取

在已有的基于SPH方法建立的风沙运动的物理模型中，因为编程实现难度以及研究内容的侧重点不同，都对模型做了一定程度的简化。同时，因为现在尚未有对风沙运动的内在机理形成统一的共识，在处理物理模型时，使用的方法也不同。例如，在现有基于SPH方法的风沙流研究中都未考虑沙粒的带电模型。在构建沙床面时，虽然大多数沙床面由混合粒径的沙粒构成，但不同粒径的沙粒在沙床面的随机分布并不完全符合不同其在实际沙床上的分布规律。在研究沙粒的运动过程中，大多数模型都忽略了沙粒自身的旋转和沙粒的蠕移，这将影响到沙粒的跃移高度以及当沙粒落回床面时，其击溅起沙粒的起跳角度，同时，对沙床面的形态也会造成影响。对来流风的处理上，大多简化为均匀流或对数流，对紊流的影响并未加以考虑。

当前，限于建立的物理模型规模，研究的主要内容仍是单颗沙粒运动的微观特征，对大尺度的风沙运动特性如风沙流的输移特征、沙丘的形成过程等还有待进一步的研究。在使用SPH方法建模时，对方法中一些重要参数的选取，如：核函数、光滑长度、时间步长等大都根据已有的经验，或考虑到计算效率予以选择，这些都将影响到模拟结果的精度[11]。

3 边界条件

在SPH中，对边界条件的处理仍然没有一致的结论。在风沙运动的模拟中，金阿芳等人在他们的工作中通常在沙床面下方设定一个固定边界，防止粒子向下穿透，气体上方设定自由边界或固定边界，计算域的左右两端往往是周期边界，即运动出计算域的粒子将携带已有的物理信息从另一侧重新进入[5]。自从Monaghan于1994年首次在自由表面流动中使用了排斥边界后[12]，后来的学者又提出了一系列方法，这些方法可以归为3类：虚粒子法、排斥函数和边界积分。其中，因为虚粒子法和排斥函数的应用条件相对较为简单，其使用较为普遍。虚粒子法是在边界外排布一组虚粒子，参与问题域内粒子的粒子近似计算中，使得原问题域内粒子的支持域不被截断，提高了近似精度。这种方法虽然在处理边界条件时十分有效，但同时它无法完全阻止粒子穿透边界，或会产生非物理分离。镜像粒子是另外一种虚粒子方法，它通过将粒子沿着边界生成的镜像来防止粒子穿透，但这种方法只适用于简单几何形状的边界。排斥函数的方法则用另外一种方式使用虚粒子，它在边界处布置一排虚粒子，对接近虚粒子的粒子产生排斥力，排斥力大小和两种粒子之间的距离有关。Monaghan和Kajtar随后对这种方法进行了改进，解决了当粒子平行于边界运动时受力不稳定的问题[13]。但是这种方法的缺陷在于，靠近于边界的粒子，其支持域被截断，无法得到很精确的结果。

为了修正上述在边界处存在的缺陷、提高边界精度，不断有学者提出对SPH方法的改进形式。如Dilts等人基于伽辽金近似构造的移动最小二乘光滑粒子法(MLSPH)[14]， Chen建立的正则化形式的修正光滑粒子法 ( CSPH)[15]， Liu使用泰勒展开在非连续区域两端分段构造的非连续光滑粒子法 ( DSPH)[16]，Liu对核函数采用校正函数进行修正的再生核粒子方法 ( RKPM)[16]。总体而言，当前在发展SPH新的边界处理方法方面取得了一定进展，改善了边界出的计算精度。但是在边界处理上仍存有许多未被发现的问题，有待进一步解决。

4 计算效率

在现有的基于SPH方法对风沙流的模拟中，大多数为二维模型。虽然因为SPH程序的特点，由二维环境扩展到三维并不困难，但粒子数的增加在提高结果精度的同时也将大幅增加计算量。为兼顾精度和计算效率，目前通常忽略风沙流在平行于地面、垂直于气流方向上的变化，将风沙运动模型简化为二维情况。在SPH程序的计算用时中，绝大部分时间是用于其粒子搜索过程。目前在风沙运动模拟中较常使用的全配对粒子搜索法最为简单，其基本原理是：在SPH方法中若要得到粒子i的物理信息，首先要得到其支持域内所有粒子的信息，对其叠加求和后加权平均得到粒子i的近似解。所以在每个时间步的开始都需要在粒子i的支持域内进行粒子搜索，在更新了粒子i的信息后，再进行下一个时间步。当粒子的个数为n时，则每个粒子需要搜索n-1次，则每一步最终的计算量为 O( n2)。因此，在风沙流这种粒子数众多的问题中，计算速度较慢。此外，其他的粒子搜索法还有链表搜索法，但当问题域在运算过程中发生大的变形时，这种方法仍然效率较低[17]。

近来，基于图形处理单元(GPU)的开源代码快速发展，逐渐应用到数值计算当中。因为SPH非常适合使用GPU进行数据并行处理，已经发展出来大量基于GPU的开源代码，例如Hérault等人提出的GPUSPH[18]，Crespo等人的DualSPHysics[19]。它们可以在个人计算机上只用数小时便模拟上百万个粒子，这对之前所需的几个月时间是一个巨大的飞跃，也为SPH方法应用于工程实践提供了可能。

5 总结

相较于传统的网格方法，SPH以其无网格、具有粒子性质的特点使得其在应用于诸如风沙流等大变形问题上具有优势。目前，SPH方法在计算稳定性和可行性上都取得了一定的进展。但若要将基于SPH的风沙流模拟应用到工程实践中，仍需要解决以下问题：

5.1 对风沙运动的模拟还有赖于风沙物理学的进一步发展。当前风沙物理学仍处于发展过程，尚未形成完整成熟的理论。在使用控制方程对风沙两相流进行数学描述时，还无法将影响因素考虑完全，如还未考虑能量耗散等问题。由此，使用SPH方法对风沙运动进行模拟时还需要风沙物理学理论的发展和完善。

5.2 使用SPH方法计算效率还需要进一步提高。由于SPH方法的计算量十分巨大，用SPH方法模拟风沙流还仅限于二维情况，对三维风沙运动的模拟还未出现。当前已有国内外学者针对GPU并行算法展开了研究并取得了巨大成果，未来将基于GPU的并行运算引入到风沙运动的模拟中，将进一步扩大模拟规模，提高模拟精度。

5.3 SPH方法还需要进一步完善。当前，SPH方法在某些领域的模拟中取得了较好的结果，但其还应该继续完善以下4个方面：数值稳定性、收敛性、边界条件以及自适应性。随着SPH方法不断地发展和完善，使用SPH方法对风沙运动进行模拟的精度将更加可靠。