基于SPH方法的立面二维涌浪数值模拟

缪吉伦^1,2，陈景秋¹，张永祥¹

(1. 重庆大学 资源环境学院，重庆，400015；

2. 重庆交通大学 西南水运工程研究所，重庆，400016)

摘要：基于无网格SPH法模拟块体下滑所引起的涌浪流动过程，并将模拟得到的自由表面速度、高度变化结果与解析解以及VOF紊流模型计算结果进行比较。在此基础上，计算滑块下滑速度分别为0.2，0.6，1.0 m/s时的自由面变化。结果表明：滑块下滑引起水槽自由液面变化并形成卷碎波，并有水气混掺现象。自由面孤立波波形、水平速度与解析解以及VOF紊流模型计算结果吻合较好。下滑速度较慢时波形平坦，无卷碎波产生；滑速大时水气完全分离；涌浪形成的波峰峰值随滑块入水速度增大而增加。SPH方法可进一步用于模拟水体中由大型崩岸产生的高速远程滑坡涌浪或复杂的海啸模拟。

关键词：SPH；无网格法；涌浪；数值模拟

中图分类号：TV145          文献标志码：A         文章编号：1672-7207(2012)08-3244-06

Two-dimension numerical simulation of impulsive waves by SPH method

MIAO Ji-lun^1,2, CHEN Jin-qiu¹, ZHANG Yong-xiang¹

(1. College of Resource & Environment Engineering, Chongqing University, Chongqing 400015, China;

2. South Westen Research Institute of Waterway Transport Engineering,

Chongqing Jiaotong University, Chongqing 400016, China)

Abstract: The meshless method of smoothed particle hydrodynamics (SPH) was used to simulate impulsive wave generated by a rectangular block falling inside water. The velocity and height change of free surface were contrasted with the analytical solution and VOF turbulent flow model results. The free surface transition was calculated at slide falling speed of 0.2, 0.6 and 1.0 m/s, respectively. The results show that sliding block results in free surface transformation, scroll wave generation and mixture of liquid with air. The solitary wave profile and horizontal velocity agree well with the analytical solution and VOF results. The wave profile is flat with no scroll waves at slowly sliding speed; the water is separate from the air at high speed. The peak wave crest increases with the increase of the sliding speed. The SPH model can serve as useful model of Large-scale rock avalanches from steep cliffs near substantial bodies of water, or of complex tsunami prediction.

Key words: smoothed particle hydrodynamics (SPH); meshless method; impulsive wave; numerical simulation

光滑粒子流体动力学法(SPH法)是一种无网格的纯Lagrange方法，1977年由Gingold等^[1]提出，最初主要是用于解决天体物理中在流体质团无边界情况下的三维空间任意流动的计算问题。随着计算理论日趋成熟，近年来逐渐在流体计算领域引起了关注^[2-9]。水库库岸或河岸由于土体特性或自然条件变化常会发生崩岸滑坡。高速崩岸滑坡往往引起自由水面的剧烈变化，大型滑坡往往产生巨大的涌浪。滑坡涌浪不仅对其波及的水域造成冲击，威胁航行船只安全，随着涌浪的传播和叠加，还可能会造成溃坝等重大事故。滑坡还会堵塞河道，如我国长江著名的云阳鸡扒子及新滩滑坡等，使长江航运一度中断，造成了巨大的经济损失。目前滑坡涌浪的研究手段主要采用物理模型和数值模拟方法。Iwasaki^[10]在固定水深的渠道中对滑坡的水平运动进行了试验研究；Monaghan等^[9]在Scott Russell试验的基础上，对块体快速垂直落入不同水深的水池的过程进行了物理试验。宋新远等^[11]基于FLUENT软件，结合RNG κ-ε湍流模型并采用VOF方法，模拟了二维滑坡涌浪。袁晶等^[12]采用动网格技术，建立了可变网格下的平面二维水库滑坡涌浪FVM数学模型。

1  SPH基本原理

与传统的基于网格的方法如有限差分法FDM和有限元法FEM相比，SPH法的主要优点在于不需要使用任何提前定义的提供结点连接信息的网格，它用一系列任意分布的粒子质点来代表整个连续介质流体并估计相应的偏微分方程。

SPH基于以密度、速度、能量等为变量的偏微分方程组，将描述场的函数用“核函数逼近”近似表达为任意函数和核函数的乘积的积分^[13]：

        (1)

式中，x与x′为计算域内两点坐标；áf(x)?表示坐标x处的核估计值；E为整个求解区域； f(x′)为坐标x′点处的场量值；h为光滑长度；W(x-x′，h)称为核函数，它有2个自变量：粒子间距离|x-x′|和光滑长度h。

2  流体SPH数值方法

2.1  Navier-Stokes方程的SPH粒子近似式

流体力学中拉格朗日型的控制方程^[14]为

质量守恒方程：

             (2)

动量守恒方程：

        (3)

式中：ρ为密度；t为时间；v为粒子速度；p为压力；g为重力加速度；υ₀为水的动力黏性系数(1.0×10^-6 m²/s)。按照SPH 的粒子近似过程, 可以将以上Navier-Stokes方程离散成以下的粒子近似方程^[15]：

           (4)

        (5)

式中：v_ij=v_i-v_j；N为粒子i的支持域中粒子的总数；m_j为粒子j的质量；W_ij为粒子j对粒子i产生影响的光滑函数，r_ij为粒子i和粒子j的距离。

状态方程：通过引入人工压缩性，把一般的不可压缩流体看作可压缩流体，在水的状态方程中引入压缩率可模拟具有自由表面的流动。这样不用求解压力泊松方程，编程简单。

              (6)

式中：ρ₀为粒子的初始密度；γ为系数，对于水取γ=7；B用于限制密度的最大改变量，一般作为初始压力,本文取，C为人工音速，选取总体流动最大流

速的10倍作为人工音速，即，H为水面高度。

粒子的位移变化为：

                (7)

方程(2)，(3)，(6)和(7)组成了基于SPH的离散化N-S方程组，加入定解条件就组成了完整的迭代方程组。

2.2  光滑核函数的选取

核函数除满足归一性和δ函数条件外，Monaghan认为还必须具有良好的紧支性，目前常用的核函数有高斯核函数、三次样条函数、五次样条函数等。基于计算精度和效率考虑，本文采用B-样条光滑函数^[4]：

        (8)

式中：q=(x-x′)/h；α_d为归一化常数，在二维情况下。

2.3  边界处理方法

入口边界、出口边界看成特殊自由边界处理，界面处粒子压力赋值为0或一定的外压力。在SPH 数值模拟中, 将固壁边界离散成为“边界粒子”, 固壁粒子参与控制方程计算，壁面采用滑移边界，但它们在模拟中位置固定，即设置位移为0 m。为阻止流体粒子穿越固壁边界，假定边界粒子对靠近它的流体粒子施加一个大小适当的中心排斥力作用,且排斥力只在近距离上起作用。

2.4  数值解法

先通过连续方程计算粒子的密度变化，捕捉自由表面粒子并进行自由表面粒子的密度校正，然后通过状态方程计算粒子的压力，同时计算外力引起的加速度，再通过动量方程计算粒子总的加速度。通过迭代方法实现粒子密度、速度、位置的更新。粒子搜索采用关联链表搜索技术。

对离散的SPH方程组进行时间积分一般采用显式格式, 如四阶Runge-Kutta法、标准蛙跳(Leapfrog)法等, 还可用中心差分、Verlet差分格式等。本文采用预测—校正方法^[16]。为保证解的收敛性, 时间步长Δt应满足CFL条件。

3  应用算例

3.1  数值水槽模型

为验证数值模型的有效性，本文选取Monaghan & Kos(2000)刚性滑块垂直滑入恒定水深的水槽中产生孤立波的物理模型试验作比较，计算模型布置如图l所示。

图1  试验模型剖面示意图

Fig.1  Cross-section initial condition of model

在数值模型中，水槽的长度为2.0 m，高度为0.78 m，水槽水深为0.21 m。滑块尺寸(长×宽×高) 为0.3 m×0.4 m×0.4 m，初始时刻让滑块底部没入水中0.5 cm，并距水槽左侧壁2.5 cm(图1)，这样可避免快速运动的滑块和水体接触瞬间造成的水体溅出而影响到试验环境。滑块的竖直下落速度为：

           (9)

式中：v在t时刻垂直下落的速度；D为水槽水深；z为在t时刻滑块底部距水槽底部的高度。初始速度及滑块到达槽底的速度均为0 m/s。

3.2  计算参数设置

为防止发生粒子穿透现象引起计算崩溃，水槽边界及滑块粒子均布置为3层且交错排列。模型共布置了22 732个粒子，流体粒子间距dx=dz=5 mm，初始时间步长dt= 0.1 ms。粒子布置如图2(a)所示。

3.2  计算结果分析

3.2.1  涌浪的产生及演变

计算得到的不同时刻水槽内液面变化见图2，模拟过程中假定滑块在下滑过程中为刚体运动。由图2可见，由于受到竖直下落滑块的推挤，滑块下部的水体向右快速运动，t=0.15 s时在滑块右侧附近形成第1次卷碎波。从0.35 s开始形成第2次卷碎波，持续到0.6 s左右，由于波在传播过程中的能量耗散，第2次卷碎波的强度明显弱于第1次。初始时刻，滑块距数值水槽左壁2.5 cm，滑块下落后在水槽壁和滑块之间的部分水体将高速向上运动，漫过滑体顶部后将向右移动。由于滑体的入水速度大于波速，在大约0.2 s时，紧靠滑块右侧的水体中有气泡产生，t=0.6 s时气泡完全被周围水体包围，波浪夹杂着气泡向右侧运动。随着滑体运动到池底，气泡随之变小破碎。由图2还可以看出：在t＜0.65 s时间内，右侧水体液面基本不变，涌浪对水槽右侧影响较小。

3.2.2  自由面波形与水平速度

根据Boussinesq方程得到的孤立波波形的解析解公式为

        (10)

式中： a为波幅，取a=0.092 m；c为波速，。

为了验证模拟结果的准确性，将本文计算得到的t=0.65s时自由面波形与解析解和采用VOF法计算的FLUENT模拟结果进行了比较(见图3)，自由面波形吻合较好。除与解析解略有差异外(小于0.02 m)，SPH法的计算结果与FLUENT解基本一致。

自由面波形的水平速度u的解析解表达式为：

              (11)

本文将t=0.65 s时自由面水平速度与解析解(式11)作对比，对比结果如图4所示。从整体来看，本文的计算结果和FLUENT模拟解符合较好，能够较准确地模拟滑块涌浪速度的变化。在滑块右侧附近(0.8 m≤x≤1.2 m)的数值结果与解析解存在一定的误差，是因为在滑块右侧附近的水体发生高度紊流运动，误差在0.1 m/s左右。

图2  不同时刻滑块附近自由液面变化

Fig.2  Change of free-surface at various time

图3  t=0.65 s时自由面波形比较

Fig.3  Comparison of waves of free-surface when t=0.65 s

图4  t=0.65 s时自由面水平速度比较

Fig.4  Comparison of horizontal velocity underneath wave profile when t=0.65 s

3.2.3  滑块下滑速度对涌浪特征的影响

为研究滑块下滑速度对涌浪特征的影响，本文给滑块取3种不同的速度，使其匀速下滑，滑速分别为v=0.2，0.6和1.0 m/s。图5所示为滑块刚刚滑入到水槽底部时刻的自由面波形比较图。由图5可见，滑块下滑速度较慢时波谷较平坦，整个水面平缓波动，无卷碎波产生，但对右侧水体影响范围大；当滑速v=0.6 m/s时可以明显地看到滑块右下侧有部分空气被卷入水体中；当滑块以v=1.0 m/s下滑时，由于滑速非常快，右侧水体被迅速排开，水面短时间内被抬升很高，被抬高的水体随之以很大的速度冲向滑块，在回落的过程中产生卷碎波，波浪尚未及向右传播。此时由于滑块的快速下滑，滑块右侧水气完全分离。研究还表明，涌浪形成的波峰峰值随滑块入水速度增大而增加，见表1。

图5  不同下滑速度时自由面波形比较

Fig.5  Comparison of wave profile at various falling speeds

表1  自由面涌浪高度最大值

Table 1  Maximum height of impulsive wave

4  结论

(1) 基于SPH法建立的滑坡涌浪数值模型，运用预测—校正法(分数步长法)求解不可压缩流体的SPH方程组,有效模拟了因滑体下滑激发的水体自由面变形、涌浪高度和速度场的变化，模拟结果与解析解、试验结果及VOF紊流模型计算结果一致。表明用SPH法来追踪自由液面、研究工程中的自由液面波动问题是可行的。

(2) 滑块下滑引起水槽自由液面变化。下滑速度较慢时波形平坦，无卷碎波产生；滑速大时形成卷碎波，并有水气混掺现象；涌浪形成的波峰峰值随滑块入水速度增大而增加。该数值模型对大体积高速崩岸滑坡所引起的涌浪灾害预防有一定的参考价值，也可用于复杂的海啸模拟。

(3) 由于SPH方法无需生成网格, 避免了大量的单元划分,解决了通常拉氏方法中的网格缠结和扭曲以及网格重划的问题，能够适应扭曲变形，特别是在大变形问题中自适应性强，而且克服了有限元方法中局部近似所引起的误差。SPH法通过核函数插值使场函数及其梯度函数在求解域内连续，控制方程更便于求解，在流体研究领域前景广阔。

(4) 由于滑坡涌浪过程较为复杂，实际应用中对于滑体的变形特性、水体阻力以及连续下滑对水体自由面变形的影响，速度场、压力波对构筑物的冲击破坏等的影响还有待进一步研究。

参考文献：

[1] Gingold R A, Monaghan J J, Smoothed Particle Hydrodynamics: theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375-389.

[2] 缪吉伦, 陈景秋, 张永祥. SPH方法在自由表面流体研究中的应用[J]. 水利水电科技进展, 2011, 31(3): 20-23.
MIAO Ji-lun, CHEN Jing-qiu, ZHANG Yong-xiang. Application of SPH method to studies on free surface flows[J]. Advance in Science and Technology of Water Resources, 2011, 31(3): 20-23.

[3] 郑兴, 段文洋. 光滑质点流体动力学(SPH)及其算法特性[J]. 船舶力学, 2008(8): 550-559.
ZHENG Xing, DUAN Wen-Yang. Smoothed Particle Hydrodynamics (SPH) and its numerical behavior[J]. Journal of Ship Mechanics, 2008(8): 550-559.

[4] 孙晓艳, 王军. SPH方法的理论及应用[J]. 水利水电技术, 2007(3): 44-46.
SUN Xiao-yan, WANG Jun. Theories and application on Smoothed Particle method[J]. Water Resources and Hydropower Engineering, 2007, 138(3): 44-46.

[5] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics[J]. Journal of Computational Physics, 2003, 191(2): 448-475.

[6] Mutsuda H, Shinkura Y. An eulerian scheme with lagrangian particles for solving impact pressure caused by wave breaking[J]. Proceedings of the Eighteenth International Offshore and Polar Engineering Conference, Canada, 2008: 162-165

[7] Oger G. Two-dimensional SPH simulation of wedge water entries[J]. Journal of Computational Physics, 2006, 213: 803-822.

[8] Souto A, Delorme L, Pérez-Rojas L, et al. Liquid moment amplitude assessment in sloshing type problems with Smooth Paritcle Hydrodynamics[J]. Ocean Engineering, 2006, 33: 1462-1484.

[9] Monaghan J J, Kos A. Scott Russell’s wave generator[J]. Physics of fluids, 2000, 12(3): 622-630.

[10] Iwasaki S I. A semi-analytical and semi-empirical method for the evaluotion of tsunami generating efficiency due to landslides plunging into the water bodies[C]//Preceding of the Fourth Pacific Congress on Marine Science and Technology. Tokyo, 1990: 162-166.

[11] 宋新远, 邢爱国, 陈龙珠. 基于FLUENT的二维滑坡涌浪数值模拟[J]. 水文地质工程地质, 2009(3): 90-94.
SONG Xin-yuan, XING Ai-guo, CHEN Long-zhu. Numerical simulation off two-dimensional waters due to landslides based on FLUENT[J]. Journal of Hydrology and Engineering Geology, 2009(3): 90-94.

[12] 袁晶, 张小峰, 张为. 可变网格下的水库滑坡涌浪数值模拟研究[J]. 水科学进展, 2008(4): 546-550.
YUAN Jing, ZHANG Xiao-feng, ZHANG Wei. Horizontal 2-D flow model with variable grid for simulating surges due to landslide in reservoirs[J]. Advances in Water Science, 2008(4): 546-550.

[13] LIU Gui-rong, LIU Mou-bing. Smoothed particle hydrodynamics—A meshfree particle method[M]. Singapore: World Scientific Publishing Company, 2003: 26-32.

[14] Roubustsova V, Kahavita R. The SPH technique applied to free surface flows[J]. Computers & Fluids, 2006, 35 (10): 1359-1371.

[15] Shao S D, Gotoh H. Pressure analysis of dam-break and wave-breaking by SPH model[J]. Annual Journal of Hydraulic Engineering, JSCE, 2003, 47: 403-408.

[16] Valizadeh A, Shafieefar M, Monaghan J, Modeling two-phase flows using SPH method[J]. Journal of Applied Sciences, 2008, 21(8): 3817-3826.

(编辑 赵俊)

收稿日期：2011-09-19；修回日期：2011-12-12

基金项目：国家自然科学基金资助项目(50879097)；教育部博士点基金资助项目(2011019111036)

通信作者：缪吉伦(1971-)，男，四川阆中人，副研究员，博士研究生，从事计算流体力学研究；电话：13012368881；E-mail：jimcqcu@163.com