风沙流结构的SPH模拟研究

楚花明，金阿芳，闻腾腾，梁岚博

（新疆大学机械工程学院，新疆 乌鲁木齐 830047）

1 引言

受人类活动的影响，自然原有环境遭到了极大破坏，土地沙漠化日益加剧，沙尘灾害成为世界面临的严重自然灾害之一，极大的影响了人类生存环境和生活质量。在风沙运动的研究中，输沙量是风沙流结构的一个重要参数，它为制定有效的防沙治沙措施提供科学指导[1]，是各国学者研究的重点。目前对输沙量沿高程上的分布规律不同学者有各自的研究结论还未形成共识。根据文献[2]的研究，输沙量与高度表现指数关系，即输沙量随高度的增加呈指数减小变化。他的结论被文献[3]通过实验所证实并加以引用。然而在文献[4]的研究中，他们发现输沙量在高度上的分布并不严格遵守指数分布。当前对风沙流结构的研究主要通过野外观测、风洞试验和数值模拟两种手段进行。目前国内外关于输沙量有三种表达式，国内的主要结论还是指数函数。但受地表属性的影响，近地层5cm内可能发生明显变化，但上层还是没有脱离指数函数的变化规律。由于风沙输运过程主要集中在近地层表面，沙粒浓度较高，通过风洞实验观察测量很难得到精准的结果。文献[5]的风洞实验结果显示，在近地层附近，输沙量比指数分布数值要小。因此，输沙量沿高度服从指数分布是在一定条件下得到的结果，并不是一个普适性结论。以往的研究结果表明，输沙量在高程的分布结构应该与沙粒起跳时的速度分布存在关联。因此，在研究风沙流结构时对沙粒起跳速度的分布规律的研究也必不可少。

现有对风沙流运动的研究方法中主要包括：野外实地观测，风洞实验和数值模拟这三种方式[6]。随着计算机硬件和软件以及计算方法的发展，数值模拟方法因为可以避免实地观测的随机性和不可复制性，以及相对于风洞实验有成本较低和便于实现的优点而在风沙流运动的研究中越来越普及。光滑粒子流体动力学（SPH）方法作为一种无网格拉格朗日粒子法，将待处理的问题区域离散为一系列粒子，这些粒子既可以作为计算插值点，其身也携带有相关物理信息，特别适合处理水动力学、爆炸冲击以及天体运动等大变形问题，应用于风沙流运动的研究上具有独特的优势［7-8］。运用SPH方法将气体与沙粒均离散为SPH粒子，并建立气固两相之间的耦合关系，模拟了风沙流运动过程并统计分析了沙粒的起跳速度分布，以及风沙流达到动态稳定后的输沙量沿高度的分布特点等风沙流结构的主要研究内容。

2 研究内容与方法

2.1 SPH方法介绍

SPH方法最早是由文献［9］提出并应用于模拟天体运动上的。它是将连续的物体离散为一系列相互联系的SPH粒子，这些粒子携带有物体的诸如密度、硬度等物理信息，通过SPH控制方程来描述物体的宏观性质。其求解过程一般包括核近似和粒子近似两个步骤，基本方程分别如下［10］。

2.1.1 核近似

对任意函数使用积分来近似表示，形如f(x)的函数，其积分表示式为：

式中：Ω—包含x的积分体积；x′—x影响域内粒子；δ(x-x′)—Dirac函数，有以下性质：

其中，使用核函数代替w(x-x′，h)，即狄拉克函数完成了核近似过程。如下式：

式中：w(x-x′，h)—核函数；h—光滑半径。

2.1.2 粒子近似

通过对相邻粒子的值进行累加求和对函数的核近似表达式进行近似，即产生目标粒子处函数的粒子近似式。假设粒子的体积为ΔVj，则粒子的质量为：mj=ΔVjρj（式中：ρj—粒子j的密度）。用ΔVj取代前述方程中粒子j处的无穷小体元dx′，则f(x)的连续SPH积分表示式为：

最终，粒子i处函数的粒子近似式可写为：

式中：WIJ=W(xi-xj，h)；N—粒子i的支持域内的粒子总数。

2.2 风沙运动模型

在风沙运动中，气流拖动床面沙粒使其由静止进入运动状态，而跃起的沙粒进入空气后获得能量，在回落到床面时又击溅起其他沙粒起跳，从而形成稳定的运动状态［11］。在这个过程中，风沙流运动的一切能量都是由气体提供的，而沙粒对气体的运动具有反作用，改变气体原先的运动状态，从而形成了风沙流这一复杂的、非线性且具有自组织性的两相流系统［12］。

2.2.1 气体控制方程

对于气体，使用Navier-Stokes方程来描述其运动规律，其主要包括质量守恒方程、动量守恒方程公式。

式中：uβ—风速在β方向上的分量；ρ—气体密度。

式中：ν=μ/ρ—运动粘度；μ—气体的动力粘度；P和ρ—压力和密度；u—气体速度；α—方向分量。

式中：Rα—摩尔气体常数，对于干燥气体，Rα=287J/(kg·K)。

2.2.2 沙粒的运动方程

沙粒在运动过程中的受力情况，如图1所示。从图中可以看出沙粒在空中的受力情况较为复杂，所以将在沙粒运动过程中作用不明显和目前尚无法量化的力，诸如萨夫曼力FS、Basset力FB以及电磁力忽略掉，则可以得到沙粒的运动方程为［13］：

图1 沙粒受力分析Fig.1 Force Analysis of Sand

式中：FD—气体拖曳力；FL—气流升力；Fg—沙粒自身重力；FM—马格纳斯力。

2.2.3 沙粒的碰撞

沙粒之间的碰撞是风沙运动过程中的一个重要机制，沙粒在落回床面以及在起跳进入空中运动后有一定概率与其他沙粒发生碰撞。将发生碰撞的沙粒分别命名为i，j，k，i是一颗发生跃移后从空中降落的沙粒，如图2所示。图中：j—正在沙床上做蠕移运动的沙粒；k—沙床上的一颗随机沙粒。假设，沙粒i以水平夹角为α的速度Vi与以速度Vj运动的沙粒j发生碰撞，而沙粒j又与沙床上静止的沙粒k发生碰撞。碰撞点如图中标识为1，2。沙粒i与沙粒j碰撞时两沙粒的质心连线和水平方向的夹角记为β，沙粒j与沙粒k的质心连线与竖直方向的夹角记为γ。

图2 沙粒与沙床以及与空中沙粒的碰撞Fig.2 The Collision of Sand-Bed and Sand-Sand

不同于其他模拟方法单独设定碰撞概率来模拟碰撞，SPH方法本身就会在每一时间步自动搜索目标粒子支持域范围内的沙粒，对支持域内判定为沙粒的粒子会根部光滑长度和两粒子之间的距离来判断是否发生碰撞，对符合碰撞条件的使用动量定理来计算出碰撞后两粒子的速度。

2.3 计算模型及参数设置

主要是基于沙粒与气体粒子之间的受力机制来模拟沙粒的运动过程并分析风沙流结构特点，分别模拟了在摩阻风速为0.1（m/s）和0.5（m/s）时，0.1mm、0.25mm两种粒径沙粒的运动过程。模拟过程中涉及的各类参数，如表1所示。

表1 算例中气体与沙粒的参数列表Tab.1 Parameter List of Gas and S and in Numerical Example

在模拟区域的上下边界各设置一层虚粒子防止粒子穿透，模拟区域的左右边界采取周期性边界条件模拟粒子的运动过程，如图3所示。初始来流风速遵循式（11）所示的对数分布[14]，如图4所示。

图3 计算区域初始状态示意图Fig.3 Initial State of Calculated Area

图4 风速廓线图Fig.4 Wind Speed Profile

式中：uy—高度为y处的风速；k—冯·卡曼系数；Y0—沙床面粗糙度。

3 结果与分析

3.1 模型验证

我们选取了6颗具有典型跃移沙粒运动特征的沙粒，所选粒子位置，如图5所示。把它们在各个时间步的位置依次连接得到它们的运动轨迹，其中1、4、6号粒子为直径0.1mm的沙粒，2、3、5号粒子为直径0.25mm的沙粒，如图6所示。从图中可以看出，沙粒的跃移轨迹具有非常明显的抛物线特性，这验证了我们建立的沙粒受力模型是正确的。沙粒的跃移高度和长度会因沙粒的起跳速度和角度不同而不同：图中1号和5号沙粒在跃移结束后落回床面不再运动；2号沙粒在落回床面后因与床面上的沙粒发生碰撞而反弹，反弹后会做非常短距离的跃移；3号、4号和6号沙粒在落回床面后以较大的反弹速度和角度重新弹回气流中，其中6号沙粒的反弹速度很大；可见，大多数的跃移沙粒在落回床面时都会因与其他沙粒的碰撞而反弹，但反弹的速度和角度根据每个沙粒的初始速度和位置不同有较大差异。

图5 所选粒子位置示意图Fig.5 Position of Selected Particle

图6 沙粒的跃移轨迹Fig.6 Jump Trajectory of S and Particles

3.2 沙粒的起跳速度分布

在风沙流结构的研究中，沙粒的起跳速度分布将影响到输沙量沿高度的分布。我们模拟了沙粒在沙床表面被风吹起，在空中运动后落回床面并击起其它沙粒起跳的整个过程。随机选取了100颗起跳沙粒，分别对它们起跳时的水平速度和垂直速度进行统计分析，如图7、图8所示。从图7（a）和图7（b）可以看出，沙粒的水平起跳速度分布具有明显的单峰特性，大部分的水平起跳速度分布在（0.1～0.4）m/s之间；同等粒径的沙粒在来流风速较大时，其水平起跳速度较大；0.1mm和0.25mm两种粒径沙粒的起跳速度分布区别不大。图8（a）和图8（b）显示了跃移沙粒垂直起跳速度的概率密度分布。可见，两种粒径的沙粒其跃移的垂直起跳速度呈指数分布，随着垂直速度的增大起跳数量呈指数下降。相对而言，小粒径沙粒的垂直起跳速度要大于同等摩阻风速条件下的大粒径沙粒，同一粒径的沙粒在两种摩阻风速下的垂直起跳速度分布基本一致，说明沙粒粒径对垂直起跳的速度影响更大。

图7 沙粒起跳水平速度概率分布Fig.7 Distribution of Horizontal of Take-Off Sand

图8 沙粒起跳垂直速度概率分布Fig.8 Distribution of Vertical Velocity of Take-Off Sand

不同粒径沙粒在不同摩阻风速下的起跳速度概率分布，如图9所示。可以看出，随着起跳速度的增加，沙粒的起跳数量经历了一个先增加达到峰值后减小的过程，单峰向左偏移，这符合gamma分布。这与Anderson和Hallet等人的研究相吻合[13]。同时，在同一粒径条件下，摩阻风速为0.5（ms-1）时，其起跳速度分布函数概率密度相对于摩阻风速为0.1m/s时要向右分散。

图9 跃移沙粒起跳速度的分布Fig.9 Take-off Velocity Distribution of Leaping S and Particle

3.3 单宽输沙率

当沙床面上起跳的沙粒数与落回床面的沙粒数相差不大时，风沙运动进入了稳态。两种粒径沙粒在不同摩阻风速下，单宽输沙率随时间的变化关系，如图10所示。

图10 单宽输沙率随时间步发展的过程Fig.10 The Process of Single Width S and Transport with Time

从图中可知，在相同粒径下，摩阻风速为0.5（ms-1）的单宽输沙率要显著大于摩阻风速为0.1（ms-1）时，而在相同的摩阻风速下，粒径较小沙粒的单宽输沙率要大于粒径较大的沙粒；不同粒径沙粒的单宽输沙率都随着时间步的增加而先增加后稍有减小并趋于稳定，在一定时间步后，单宽输沙率稳定在一数值上下。当风沙运动处于动态平衡时，0.1mm粒径沙粒在摩阻风速0.1（ms-1）和0.5（ms-1）下的单宽输沙率分别为：0.015（kgm-1s-1）、0.0461（kgm-1s-1），0.25mm粒径沙粒分别为0.0074（kgm-1s-1），0.038（kgm-1s-1）。由图中可以看到，0.1mm粒径沙粒的单宽输沙率在两种摩阻风速下的增长速度都高于0.25mm粒径沙粒，在0.5（ms-1）摩阻风速下更快达到动态平衡。但在0.1（ms-1）摩阻风速下，0.25mm粒径沙粒要先达到动态平衡，这与风速较小且0.25mm粒径沙粒起跳并不充分有关，其在起跳沙粒数量较少的情况下达到动态平衡。

3.4 输沙量

当沙床上沙粒的吹蚀和堆积数量相当时，风沙运动处于动态平衡，此时输沙量的大小可由该高度处沙粒的数量来表示，因此，输沙量沿高度的变化即为在单位时间内单位面积上沙粒浓度沿高度的变化。由上文得出沙粒的起跳速度服从gamma分布，可见，沙粒以不同的速度进入到空中，其所能达到的最大高度也不同。文中统计了在风沙流稳定后，输沙量沿高度的变化。随着高度的增加，输沙量先增加而后减小呈指数衰减，在接近沙床面的上方附近存在一个最大值，如图10所示。

这与文献[14]所做的模拟结果趋势上相同。由图11（a）和图11（b）可以看出，在沙粒粒径为0.1mm和0.25mm两种情况下，摩阻风速为0.5（ms-1）时，其输沙量达到最大值的高度要高于摩阻风速为0.1（ms-1）时，输沙量也要远大于摩阻风速为0.1（ms-1）时，这是因为随着风速的增加沙粒起跳时携带的能量也增加，所以输沙量的最大值会随着风速的增大而增大，即“象鼻效应”。

图11 输沙量沿高度的分布Fig.11 Distribution of S and Transport Along Height

4 讨论

（1）影响沙粒起跳速度的因素有很多，文中只考虑了0.1mm和0.25mm两种粒径沙粒在0.1（ms-1）和0.5（ms-1）两种摩阻风速下的运动过程。由起跳速度的分布情况可知，在文中条件下沙粒的粒径对起跳的速度影响更大，这是因为在沙粒粒径为0.1mm时，给出的两种摩阻风速下都能够充分起跳；沙粒粒径为0.25mm时，两种摩阻风速下的起跳都不够充分。

（2）输沙量沿高度的分布和沙粒的起跳速度分布具有很大的相似性，沙粒的起跳速度具有一单峰值，并偏移在低速度区域内，输沙量则沿着高度先增加而后指数减少，也存在一个输沙量峰值，同样偏移在距离地面较近的高度，该高度值与沙粒的起跳速度分布有关，在沙粒的起跳速度分布中存在一个速度区间，在该速度区间内，起跳的沙粒数量最多，以该速度起跳的沙粒所能达到的高度即为输沙量达到最大值时的高度。

（3）不同于其他参数化建模的模拟方式，文中基于SPH方法建立的风沙流模型是从沙粒及其气体粒子的受力情况入手，将沙粒的受力情况量化，选取沙粒在运动过程中一直起主要作用的力，忽略掉一些次要的作用力，建立风沙两项间的耦合机制以及沙粒之间的碰撞机制，自然得到沙粒的起跳速度及其它一系列运动情况，更容易探究沙粒的受力情况，符合物理本质。

5 结论

（1）不同粒径的沙粒，其风沙运动演变过程基本一致，其输沙量都随着摩阻风速的增加而增加，随着时间的增加而先呈指数形式的增长而后略有减小，最后趋于稳定。粒径大的沙粒其输沙量达到稳定状态需要的时间略长。

（2）通过SPH方法对沙粒运动过程进行模拟，发现其起跳速度呈gamma分布，输沙量沿高度先增加而后呈指数减小。摩阻风速较大时，输沙量的达到的峰值以及该峰值对应的高度都大于摩阻风速较低时的相应值。

（3）输沙量沿高度的分布和沙粒起跳的速度分布有内在关联，输沙量其饱和层的高度和沙粒起跳速度分布的峰值所对应的速度有关，以该速度起跳所达到的高度即为饱和层高度。

（4）SPH方法基于自身特点非常适合风沙两相流的模拟。课题组首次运用SPH方法，建立了风沙两相流的运动模型，模拟结果与理论吻合较好，在风沙运动的机理研究中具有很大应用前景。