溃坝流动是一种典型的自由表面流动，具有运动边界和复杂的大变形问题，自由表面流动模拟一直是计算流体力学领域的一大难题。19世纪70年代，圣维南提出一维非恒定流基本方程组，奠定了溃坝水流问题研究的理论基础^[1]；1892年Ritter通过特征理论获得了著名的Ritter瞬间溃坝问题解^[2]；Hunt B对溃坝波的近似解和摄动解做了一系列的研究工作^[3]；Abdolmaleki等采用FVM法结合VOF格式追踪自由表面模拟了二维溃坝问题^[4]。杨小亭采用MAC法求解不可压缩非定常Navier-Stokes方程^[5]，张永祥等建立基于CE/SE格式的溃坝洪水波计算模型^[6]，其他一些学者通过传统的网格方法(如有限差分法、有限体积法)对溃坝流动过程进行了研究^[7-8]，但基于网格的方法在模拟水流飞溅、高速冲击等大变形问题时非常复杂，甚至会造成求解错误。文中采用无网格方法—光滑粒子流体动力学(SPH)法研究二维溃坝问题。

1 SPH基本原理

SPH(smoothed particle hydrodynamics)是一种无网格粒子法，由Lucy & .Gingold & .Monaghan于1977年提出，最先用于天体物理现象的模拟，现已被广泛地应用于连续固体力学和流体力学中。SPH法是一种纯拉格朗日无网格法，相对于传统的基于网格的数值方法它具备一些独特的优点，如SPH法的自适应性、问题域用粒子表示和将粒子作为场变量近似的框架。由于SPH法自适应性，公式的构造不受粒子随意分布的影响，因此，能够解决一些水下爆炸、波浪破碎、高速水流、高速冲击碰撞等极大变形的问题，SPH法日益成熟并且在诸多领域快速发展^[9]。

SPH方程的构造分2个关键步骤，第一步为积分表示法，又称核近似法；第二步为粒子近似法。核近似法是对任意函数和光滑函数进行逐步积分。函数积分表达式再通过最近相邻粒子的值进行累加求和近似。在SPH方法中，函数积分表达式定义为

$ f\left( x \right) = \int\limits_\Omega {f\left( {x'} \right)W\left( {x - x',h} \right){\rm{d}}x'} , $

(1)

式中：Ω为整个求解域；f为三维坐标向量x的函数；W(x-x′, h)为核函数; |x-x′|为粒子间距；h为光滑长度。任意一点i的函数值都可以表示为SPH的离散形式：

$ f\left( {{x_i}} \right) = \sum\limits_{j = 1}^N {{m_j}\frac{{{f_j}}}{{{\rho _j}}}W\left( {\left| {{x_i} - {x_j}} \right|,h} \right),} $

(2)

式中：N为支持域内相邻粒子数；m_j为粒子j的质量；ρ_j为粒子j的密度。|x_i-x_j|为粒子i和粒子j的间距，h为光滑长度。

2 溃坝流体SPH数值方法 2.1 用SPH法解Navier-Stokes方程组

拉格朗日描述下的流体控制方程为^[10-12]

连续方程：

$ \frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} = - \rho \nabla \cdot v, $

(3)

动量方程：

$ \frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - \frac{1}{\rho }\nabla P + g + \Theta , $

(4)

$ \Theta = {\upsilon _0}{\nabla ^2}v + \frac{1}{\rho }\nabla \tau , $

(5)

式中:P为粒子压力；ρ为流体粒子密度；t为时间；g为重力加速度；v为粒子运动速度；Θ为紊动项；υ₀为水的运动粘性系数；τ为应力张量。

模拟流体动力学问题时，须对算法进行一些特别的处理，避免求解域内产生非物理振荡，研究可知，加上人工粘度可减少非物理振荡且防止粒子接近时的非物理穿透。文中采用Monaghan(1989)提出的公式，具体表达式为^[13]

$ \prod\nolimits_{ij} { = \left\{ \begin{array}{l} \frac{{ - {\alpha _\Pi }{{\bar c}_{ij}}{\varphi _{ij}} + {\beta _\Pi }\varphi _{ij}^2}}{{{{\bar \rho }_{ij}}}},{v_{ij}} \cdot {x_{ij}} < 0;\\ 0,{v_{ij}} \cdot {x_{ij}} \ge 0. \end{array} \right.} $

(6)

式中: $\bar c{_{ij}} = \frac{1}{2}{\rm{ }}({c_i} + {c_j})$，${\varphi _{ij}} = \frac{{{h_{ij}}{v_{ij}}\cdot{x_{ij}}}}{{\left| {{x_{ij}}} \right|{^2} + {\varphi ^2}}}$，${v_{ij}} = {v_i} - {v_j}$，$\bar \rho {_{ij}} = \frac{1}{2}({\rho _i} + {\rho _j})$，${h_{ij}} = \frac{1}{2}({h_i} + {h_j})$，${x_{ij}} = {x_i} - {x_j}$。

在以上方程中α_Π和β_Π为标准常数，一般取值在1.0左右，c为声速，v_ij为粒子i对粒子j的相对速度，x_ij为粒子i对粒子j的相对位移，h_ij为粒子i和粒子j光滑长度之和的平均值，因子φ=0.1h_ij可以避免粒子靠近时产生的数值发散。式中与α_Π相关的项可得到体积粘度，与β_Π相关的项可以避免在高马赫数时发生粒子间的穿透。

根据式(2)得出控制方程的SPH离散格式为

$ \frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = {\rho _i}\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}\frac{{\partial {W_{ij}}}}{{\partial {x_i}}},} $

(7)

$ \frac{{{\rm{d}}{v_i}}}{{{\rm{d}}t}} = - \sum\limits_{j = 1}^N {{m_j}\left( {\frac{{{p_i}}}{{\rho _i^2}} + \frac{{{p_j}}}{{\rho _j^2}} + {\Pi _{ij}}} \right)\frac{{\partial {W_{ij}}}}{{\partial {x_i}}} + g} , $

(8)

式中：v_ij为粒子i对粒子j的相对速度；Π为人工粘性项；p为粒子压力。

2.2 选取光滑核函数

光滑函数具有紧支性、归一性、对称性及光滑性等。目前，常用的核函数有Gaussian型、三次样条核函数、四次样条核函数及五次样条核函数等。文中选用Quintic核函数进行模拟^[14-15]。

$ W\left( {x - x',h} \right) = {\alpha _D}{\left( {1 - \frac{q}{2}} \right)^4}\left( {2q + 1} \right),0 \le q \le 2, $

(9)

$ q = \left| {\frac{{x - x'}}{h}} \right|, $

(10)

式中，α_D为归一化常数，二维问题α_D=7/(4πh²)。

2.3 边界处理方法

SPH数值模拟中，在固定边界上分布1组虚粒子，这些虚粒子对邻近边界的粒子中心作用强排斥力，排斥力只在近距离上起作用，从而阻止邻近边界的粒子非物理穿透边界^[16]。固壁粒子参与控制方程的计算，壁面采用滑移边界，但设置位移为0，即在模拟中位置固定。

2.4 状态方程

在求解可压缩流体问题的SPH方法中，粒子的运动由压力梯度产生，而粒子的压力是由状态方程中粒子自身的密度和内能计算得到。理论上不可压缩流体事实上是可压缩的，因此，需要引入人工压缩率，目的是求解压力对时间的导数^[17]。文中采用水的状态方程模拟溃坝问题:

$ P = B\left[ {{{\left( {\frac{\rho }{{{\rho _0}}}} \right)}^r} - 1} \right], $

(11)

$ B = \frac{{{\rho _0}{c^2}}}{\gamma }, $

(12)

式中：P为粒子压力；γ为常数，一般情况下取γ=7；ρ₀为参照密度；B用于限制密度的最大改变量，一般作为初始压力。C为人工音速, $ \mathit{c}=\sqrt{200gH}$，H为水面高度，粒子的位移变化：

$ \frac{{{\rm{d}}x}}{{{\rm{d}}t}} = v。$

(13)

方程(7)(8)(11)(13)组成SPH离散化的N-S方程组，加入定解条件得到完整的迭代方程组。

2.5 数值解法

先由连续方程计算得到粒子密度变化，捕捉自由表面粒子且对表面粒子进行密度校正，再由状态方程计算粒子的压力，同时，计算得到外力引起的加速度，最后，通过动量方程计算粒子的总加速度。采用迭代方法对粒子的密度、位置、速度进行更新。文中采用关联列表搜索技术搜索最近相邻粒子。

离散的SPH方程组对时间的积分一般采用显式格式，如标准蛙跳(Leapfrog)法、四阶Runge-Kutta法、Verlet差分格式及中心差分法等。文中采用Beeman预测校正法。时间步长Δt应满足CFL条件以保证解的收敛。

3 算例模拟 3.1 模型有效性验证

为验证SPH法模拟溃坝这种自由表面流体大变形问题的有效性，通过建立上下游不同水深的水槽模型，将模拟计算结果和理论解进行对比^[18-19]。

如图 1所示，模型为单宽矩形槽, 长40 m，距左侧5 m处设一闸门，上游初始水位为10 m，下游初始水位为5 m, 水体密度为1 000 kg/m³。水体在初始时刻处于静力平衡状态，通过迅速提升闸门来模拟溃坝问题。

坝址下游有水的情况下，Stoker推导得到了理论解，考虑的初始条件为：矩形河道、平底、无阻力、下游有水、瞬间全溃、溃坝前上下游为静水。如图 3所示，将t=30 s, 60 s时液面高度的数值解和Stoker理论解进行对比，从图 3中可见，数值解和理论解吻合较好，误差在允许的范围内，表明建立的溃坝数值模型是有效的。导致误差的主要原因为，理论解假定了河底和边界无摩擦力及SPH方法数值计算精度受核函数形式、边界条件处理、时间步长选择等因素的影响。

3.2 2种消能坎的溃坝流动模拟

参考Martin等的水槽模型，在距水槽左侧位置x=0.6 m处布置矩形障碍物，分别模拟了无孔障碍物和有孔障碍物2种情形下的自由表面流动现象。通过比较2种模型中水流撞击障碍物及右壁时的水位爬高和水流压强，分析比较2种障碍物的消能效果。

如图 4~图 5所示，水槽模型为单宽矩形槽, 长1.2 m，距左侧0.4 m处设一闸门，上游初始水位为0.5 m，水体密度为1 000 kg/m³。水体在初始时刻处于静力平衡状态，通过迅速提升闸门来模拟水柱溃坝时的情形。距左侧位置x=0.6 m处布置障碍物，高0.2 m。

对整个计算域离散化，初始粒子布置分别如图 6和图 7所示，时间步长Δt根据CFL条件给出。粒子初始压力由静水压力给出。数值模拟计算的模型参数如表 1所示。

图 8显示数值模拟计算后溃坝在几个时刻中粒子分布情况的对比，从图中可以看出，溃坝后水位下降，水体向右高速运动，水体在翻越障碍物前有序流动。

在t=0.25 s时，水体撞击障碍物后向上飞溅，从图中可以看出水体撞击有孔隙障碍物时的飞溅高度要低于撞击无孔障碍物时的高度；之后在重力的作用下快速跌落翻卷，在t=0.7 s时水体下落撞击底部后，无孔障碍物的右边有封闭的气穴，而有孔障碍物右边没有；在t=0.9 s时，水体撞击右壁面后水位迅速爬高。

为了比较无孔障碍物和有孔障碍物的在溃坝后的消能效果，如图 9所示，文中计算了水体在撞击右壁面的压强对比，从图中可看出在有孔障碍物的阻碍下，水体撞击右壁面时的压强要小于无孔障碍物时的情形，所以在一定范围内，有孔障碍物有更好的消能效果。

4 结论

根据SPH方法的基本原理建立了溃坝数值模型，通过将数值解和Stoker理论解进行对比分析，验证了模型的有效性，模拟中采用预测校正法求解弱可压缩流体方程组，较好地模拟了溃坝水体自由面变形、飞溅、融合等复杂自由表面现象。对障碍物有孔和无孔2种情形进行模拟，通过计算及对比分析，表明有孔障碍物在一定程度下有更好的消能效果。