ABSTRACT

The motion of surging breaking wave is simulated here with the help of SPH (smoothed particle hydrodynamics) method. According to the mesh-less property, the wave shape is described by the profiles of particle distributions, and the evolution of the tiny breaking waves is described by particle velocity vectors. Comparing with the experimental photographs, it is found that the SPH method performs well in wave breaking problem. Then the breaking position is located by the evolution of wave peak. And the parameter at the breaking position is very approximate to empirical data.

Keywords:Breaking Wave, Surging Breaking Wave, SPH

基于SPH的激破波破碎过程的 数值模拟

吴宗铎¹，许斐²，严谨^1*，张建³

¹广东海洋大学海洋工程学院，广东 湛江

²上海船舶研究设计院，上海

³广东省航道局，广东 广州

收稿日期：2018年11月13日；录用日期：2018年11月29日；发布日期：2018年12月6日

摘 要

本文利用SPH方法模拟了激破波的运动过程。通过该方法的无网格特性，可利用粒子分布图来描绘波浪的形态，并利用粒子的速度矢量图来描绘碎浪的变化。和实验的图片进行对比，发现SPH方法在波浪破碎问题上表现出良好的计算性能。随后利用波峰位置的演化来确定发生破碎的位置，得到的破碎位置处的参数也和经验公式比较接近。

关键词 :波浪破碎，激破波，SPH

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

波浪破碎是波浪传播到近岸区域后，受到地形的限制而产生的变形、折射、绕射、反射等作用后产生的剧烈形态变化。波浪破碎对海岸带、沿海城市、交通枢纽、国防建设等一系列的基础设施造成了威胁，对沿岸人民的生活有较大的影响，所以对于近岸地区波浪传播与破碎运动的研究十分之重要 [1] 。

根据破碎后的形态，波浪破碎可分为崩破波、卷破波、激破波三种类型 [2] 。三种类型中，激破波的破碎形式最复杂，其现象为：在波峰前沿根部开始出现破碎现象，随后波峰前沿部分呈非常杂乱的破碎状态，并沿斜坡上爬。波峰前沿的杂乱部分，由于水质点存在极大的不规则特性，则成为了数值模拟的大难题。河海大学诸裕良 [3] 利用波浪水槽实验研究了复合坡度下的破碎变化形态，虽然地形符合激破波的条件，但实验过程中却没有观测到明显的激破现象。郄禄文 [4] [5] 则用SPH方法模拟了三种破碎形式下的速度变化和波峰压力，但对激破波的形态并未作特别的关注。天津大学赵京津 [6] 做了激破波的相关水槽实验，并用光滑粒子方法模拟了波浪的爬坡过程，但未针对激破波的现象做细致的研究。

由于在波浪破碎的过程中，波浪的自由表面发生扭曲，随后大量的水质点发生分离，整个运动过程中自由表面的形态产生了较大的变化。如用常规网格来计算，处理难度较大。面对形态复杂的激破波，SPH方法是一种比较实用的数值方法。由于流体的自由表面在粒子形态下被离散化，流体模拟技术能逼真地模拟各种复杂流体现象，能直观地呈现各种自由表面的复杂变形 [7] 。目前，在波浪破碎的数值模拟中，SPH已大量得到应用 [8] [9] [10] 。

本文利用SPH方法对激破波展开研究，首先将通过试验结果来验证SPH方法的准确性，随后将对激破波展开细致的研究，结合数值模拟图形来分析激破波波前散碎波浪的特性。由于波浪形态用光滑的离散粒子来表示，得到的结果具备很好的直观性。对激破波发生波浪破碎的时间和位置，本文将做重点的关注。

2. 激破波形式及发生条件

2.1. 激破波简介

激破波是破碎形态中形态变化比较剧烈一种破碎现象。在波浪的运动推进的过程中，在近岸的斜坡上发生破碎时，首先前沿面会首先变得陡立，然后卷曲拱形。但是，当坡度比较陡的时候，波浪无法形成卷曲形态进行翻卷，而是会转化成大量的散碎的水花向前运动。此时，波峰前后逐渐变的非常不对称，波浪的形态已经基本消失。

而激破波此时所表现出的现象为波峰前沿根部开始出现破碎现象，随后波峰前沿部分呈非常杂乱的破碎状态，并沿斜坡上爬。但是，由于波浪前端的形态已经转化为杂乱的水质点，这些水质点无法形成状态规则的自由表面。等到水质点重新落回水面以后，形成流动状态比较复杂的一股水流。整个激破波的变化形态如图1所示，波峰前存在着从陡立形态到激破形态的几个阶段。

图1. 激破波形态变化

从发生波浪破碎位置，到岸边的一段区域，称为破碎带。在破碎带中，水的深度从破碎点到岸边逐渐减小，波浪将一直延续破碎状态。水体的湍流和涡度非常强，是波浪能大量的消耗，且水质点的运动形态非线性较强，波形，能量消耗，湍流和漩涡的剧烈变化非常强烈。

2.2. 激破波发生条件

根据Galvin的实验 [2] ，波浪破碎形态可根据破波相似参数来进行区分，这里定义深水处的相似参数为 ${\xi }_0$ ，破碎处的相似参数为 ${\xi }_b$ 。它们的计算方式为：

${\xi }_0=\frac{\tan \beta }{\sqrt{H_0/L_0}},{\xi }_b=\frac{\tan \beta }{\sqrt{H_b/L_0}}$ (1)

这里，H₀和H_b分别为深水处和破碎处的波高，L₀为深水处波长，b为坡度的倾斜角。波浪破碎后的类型，区分依据如表1：

表1. 破碎形态分类

从表1中可以看出，波浪类型与坡面的斜度关系较为密切。当坡度较为缓和时，波浪能维持正常形态。随着水底的坡度逐渐倾斜，水深不断变浅，此时波长逐渐缩短，波高逐渐增大。当波高波长比增大到一定程度时，波浪维持正常形态的难度会变大。若坡度较缓和，则波峰顶部开始出现白色浪花，此时表现为崩破波。坡度倾斜度加大时，则破碎位置从顶部向下移，破碎位置在波峰中部是，则发生翻卷，表现为卷破。当坡度倾斜度很大时，则波峰根部就发生破碎。波峰根部以上的自由表面全部被零散而杂乱的浪花所取代，如图2所示。

图2. 不同坡度下的波浪形态

3. SPH基本理论

SPH方法是一种用运动的粒子来代替流体的无网格、纯拉格朗日粒子方法。物理参数的表达式主要是通过周边区域内粒子物理参数的积分插值。基本理论如下 [11] ：

3.1. 基本关系式

对于物理参数A，可以表示成关于矢量坐标 $r$ 的函数，即 $A\left(r\right)$ 。函数 $A\left(r\right)$ 可以通过周边区域内相关质点的数值积分来得到：

$A\left(r\right)={\displaystyle \int A\left(r^{\prime }\right)W\left(r-r^{\prime },h\right)}\text{d}{r}^{\prime }$ (2)

这里，h表示光滑长度，W称为核函数， $W\left(r-r^{\prime },h\right)$ 它表示在光滑长度h为特征长度的圆圈内，众多粒子的相关参数的加权函数。在实际的数值计算中，函数 $A\left(r\right)$ 往往用离散的关系式来表示：

$A\left(r\right)={\displaystyle \underset{b}{\sum }{m}_b\frac{A_b}{{\rho }_b}}{W}_{ab}$ (3)

这里，a表示需要进行拟合函数的质点，b代表拟合函数的周边相关质点。W_ab的物理意义为：a点的物理参数，需要将周边多个质点b的数值通过加权函数W_ab来进行拟合。

而W_ab的形式可以有很多种，这里列出几种常见的：

高斯型： $W\left(r,h\right)={\alpha }_D\exp \left(-q^2\right)$

二次样条型： $W\left(r,h\right)={\alpha }_D\left[\frac{3}{16}{q}^2-\frac{3}{4}{q}^2+\frac{3}{4}\right]$

三次样条型： $W\left(r,h\right)=\{\begin{array}{ll}{\alpha }_D\left[1-\frac{3}{2}{q}^2\text{+}\frac{3}{4}{q}^3\right]\hfill & 0\le q\le 1\hfill \\ \frac{{\alpha }_D}{4}{\left(2-q\right)}^3\hfill & 1\le q\le 2\hfill \\ 0\hfill & q\ge 2\hfill \end{array}$

五次样条型： $W\left(r,h\right)={\alpha }_D{\left(1-\frac{q}{2}\right)}^4\left(2q+1\right),0\le q\le 2$

3.2. 控制方程

光滑粒子的运动特性及相关参数的变化，可通过以下方程来进行确定：

• 连续性方程：

流体密度的变化用连续性方程表示：

$\frac{\text{d}{\rho }_a}{\text{d}t}={\displaystyle \underset{b}{\sum }{m}_b{v}_{ab}{\nabla }_a{W}_{ab}}$ (4)

这里符号 ${\nabla }_a$ 表示在a点的梯度值

• 动量方程：

动量守恒的关系可以表示为：

$\frac{Dv}{Dt}=\frac{1}{\rho }\nabla p+g+\Theta$ (5)

这里，为 $\Theta$ 对流项。常见的对流项包括人工粘性项、层流粘性和完整粘性(层流粘性 + 亚粒子湍流粘性)。

• 状态方程

这里采用形式简单的状态方程来描述水和气体：

$p=B\left[{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }-1\right]$ (6)

对于水，参数 $$ ， $B={\rho }_0{c}_0^2$ ，分别为水的密度和声速。对于理想气体， $\gamma =\text{1}.\text{4}$ ，B为标准大气压。

• 粒子运动方程

粒子运动轨迹遵循如下方程：

$\frac{\text{d}{r}_a}{\text{d}t}=v_a+\epsilon {\displaystyle \underset{b}{\sum }\frac{m_b}{{\bar{\rho }}_{ab}}{v}_{ab}{W}_{ab}}$ (7)

这里 $\epsilon =0.\text{5}$ ， ${\bar{\rho }}_{ab}=\left({\rho }_a+{\rho }_b\right)/2$ ，其物理意义为：当a点的粒子在运动过程中，它的运动轨迹与周边粒子的加权关系。

• 能量方程

内能e_a可以用如下方程来表示：

$\frac{\text{d}{e}_a}{\text{d}t}=\frac{1}{2}{\displaystyle \underset{b}{\sum }{m}_b\left(\frac{p_a}{{\rho }_a^2}+\frac{p_b}{{\rho }_b^2}+{\psi }_{ab}\right)v_{ab}{\nabla }_a{W}_{ab}}$ (8)

这里 ${\psi }_{ab}$ 代表粘性项。

3.3. 边界条件

SPH中常见的边界条件有动力边界和排斥力边界。

动力边界条件的主体思想为，边界粒子和水粒子一样，满足连续方程(4)、动量方程(5)、状态方程(6)和能量方程(8)，但是由于运动范围被固定，粒子运动方程(7)不满足。而排斥力边界条件则设定为边界上粒子为虚拟粒子，虚拟粒子不能被正常运动的粒子穿透，且边界粒子会给正常粒子施加排斥力。

本文所研究的对象为波浪破碎问题，该问题中，水底和斜坡的边界均与正常运动的粒子无明显的动力作用，因此这里选用排斥力边界条件。排斥力边界条件中，假设边界粒子与正常流体粒子之间的间距为r，那么排斥力 $f$ 的计算方式为：

$f=nR\left(\psi \right)P\left(\xi \right)\epsilon \left(z,u_{\perp }\right)$ (9)

这里， $\stackrel{\to }{n}$ 代表固体边界法向，代表R(ψ)为排斥力函数，其表达式为：

$R\left(\psi \right)=A\frac{1}{\sqrt{q}}\left(1-q\right)$ (10)

其中， $q=\psi /2h$ ， $A=0.01c^2/h$ 。c为声速。函数P(ξ)代表平行于墙运动时受到的持续的作用力：

$P\left(\xi \right)=\frac{1}{2}\left[1+\cos \left(\frac{2\text{π}\xi }{\Delta b}\right)\right]$ (11)

Δb是两个相邻边界粒子之间的间距。最后，函数 $\epsilon \left(z,u_{\perp }\right)$ 则是根据水深z和速度 $u_{\perp }$ 的修正。速度 $u_{\perp }$ 为水质点垂直于墙的速度分量。

4. 数值算例

4.1. 概况

本文计算时的参数设定与赵津京在文献 [6] 中的水槽实验保持一致。波浪参数为周期T为2 s，水深h为0.26 m。造波方式为推板式造波，造波板的推程E为0.12 m，斜坡的倾斜角度为10˚。如图3所示。

图3. 激破波运动初始状态

由于实验水槽的水深有限，这里用有限水深下色散关系来考虑波高H与波长L，具体表达式为：

${\omega }^2=gk\tanh \left(kh\right),{\omega }^2=g{k}_0$ (12)

这里ω为圆频率，h为水深，k为波数， $k=\text{2π}/L$ ， $k_0={\text{2π}/L}_0$ 。根据色散关系(12)可以算出k = 0.99，即有限水深L = 6.24 m。而对应在深水下的波数k₀ = 0.25，相应深水波长L₀ = 24.97m。而波高H根据造波板运动幅值E与波幅A的传递函数来计算：

$\frac{A}{E}=\frac{2{\sinh }^2\left(kh\right)}{2kh+\sinh \left(2kh\right)}$ (13)

其中波高H = 2A，根据计算可得H约为0.066 m。深水下的波高H₀与H相差不大，可根据查曲线图4来得到 [1] 。此前的波长比L/L₀也可以通过图4来确定。

图4. H/H₀和L/L₀的变化曲线

由于破碎点的位置不明确，这里利用(1)计算 ${\xi }_0$ 来判断破碎的判断条件，计算得到的 ${\xi }_0$ 约为3.43，符合激破波的发生条件。

4.2. 计算结果与实验照片的对比

本文用SPH方法模拟了激破波在爬坡过程中的破碎过程，并将计算的结果同文献 [6] 的实验结果进行了比较，对比的结果如图5所示。

图5. 不同时刻与实验的对比

从图5上可以观察到整个激破波发生破碎的全过程。在4.7 s时刻，在波浪爬坡的最高点位置，前端已开始出现较薄的一层碎浪。4.8 s开始，碎浪的范围慢慢加大，在水深很浅的前端区域，几乎全部为碎浪，到4.9 s时刻，前端的零散的碎浪依旧比较明显。但是从5.0 s开始，之前产生的碎浪重新落回水面，激破波的现象慢慢消逝。

由于部分波高转化为碎浪落回水面，此时波高已经减小，这样。从5.1 s到5.4 s，波浪的破碎形态略有转变。由于水深变浅，破碎位置更靠近岸边，此时破碎波高H_b进一步加大，参数 ${\xi }_b$ 则减小，因此后续的运动过程和卷破波有点类似。经过一次破碎过程后的波峰，顶端延伸出去，随后直接运动至爬坡运动的最高点。此时，波峰前端已不再是零散的碎浪，而是层次分明的翻卷。

4.3. 速度矢量场

为了更好的观察破碎情况，这里计算几个时刻的矢量场。从图6的速度场可以观察到波峰前的碎浪情况。由于碎浪的具有无规律性的特点，因此在矢量图6中可以观察到碎浪速度的方向比较散乱。在4.7s时刻，波浪前端的散碎的浪花较多，且分布范围也比较广，在x = 2.1 m~2.2 m的区间内几乎都被碎浪覆盖。到了4.8 s时刻，波峰继续向岸边运动并挤压了碎浪的运动空间，碎浪的分布范围开始缩小，但是在x = 2.3 m附近仍然有少量碎浪。到5.1 s，激破波的碎浪已经不太明显，此时水质点的运动方向基本斜向上方。

图6. 不同时刻的速度矢量场

4.4. 破碎位置判定

按照根据破碎的定义，最大波高受地形限制达到极限波陡时，波浪就将进行破碎。波陡定义为H/L。根据Miche提供的经验公式 [1] ，波浪在有限水深下的极限波陡有：

$\frac{H_b}{L_b}=0.142\tanh \frac{2\text{π}{h}_b}{L_b}\approx 0.142\frac{2\text{π}{h}_b}{L_b}$

进行化简后，破碎波高与破碎处水深的关系大致为：H_b/L_b = 0.89。

根据图5中各个时刻激破波形态，大约在5.0 s时刻达到最大值。该时刻下，波峰前方的碎浪区域处于最活跃状态，并开始回落至水面。破碎位置可大致认定为图7中x = 1.7 m处。而破碎处的波高则可以从图7中虚线所示的水位差中大致读出。破碎位置处，从图7上可大致估算H_b = 0.105 m，而水深h_b = 0.114 m。此时H_b/L_b = 0.92，与经验公式大致吻合。

图7. 破碎时刻的波高与水深

5. 结论

利用SPH方法，完成了激破波爬坡的数值模拟。通过对数值模拟的结果进行观察和分析，可以得出了以下结论：

• SPH方法能较好的模拟激破波爬坡过程。爬坡过程中，波峰前的碎浪在粒子分布图上不明显，但是可以通过速度矢量图来观察碎浪的运动情况。

• 激破波的破碎过程中，波峰前出现碎浪以后的一段时候，碎浪落回水面。随后波峰会继续向前运动，波高进一步增大，在满足条件的情况下，可能出现卷破波。

• 当激破波破碎完成时，前方碎浪达到最活跃状态。此时波峰位置可视为破碎位置，破碎位置的波高可通过水位线的差值来得到。SPH结果中，破碎波高与破碎水深的比值与经验公式的比值比较接近。

虽然SPH方法能较好地模拟实验水槽的激破波运动，但是在工程实际中，激破波的运动容易受到水底地形及淤积泥沙的影响。如水底坡度变化剧烈，可将地形考虑成不同坡角组成的复合坡度 [3] ；如有较多的泥沙，可把泥沙考虑成圆形的颗粒物 [12] ，再进行数值计算或模型试验。

基金项目

国家自然科学基金资助项目(11702066)；广东省自然科学基金(2017A030313275)；广东海洋大学“创新强校工程”省财政资金支持项目(GDOU2016050258, GDOU2017052503)。

文章引用

吴宗铎,许 斐,严 谨,张 建. 基于SPH的激破波破碎过程的数值模拟
Numerical Simulation of the Breaking Process of Surging Breaking Wave Based on SPH[J]. 流体动力学, 2018, 06(04): 100-108. https://doi.org/10.12677/IJFD.2018.64013

参考文献