贾进生

(山西大同大学数学与统计学院，山西大同 037009)

在海洋中，地震、海啸引起的波浪通常是孤立波，其波长较长且波高较大，对海岸附近的人类生产生活和建筑物安全有着重大影响。因此，针对孤立波开展研究，有着重要的实际意义。

目前，对于波浪的常用研究手段包括物理模型实验和数值模拟。物理模型实验是在实验室中按照一定的比尺对实际工程进行还原，用造波机产生目标波浪，通过实验研究波浪的相关问题。然而，物理模型实验对实验条件要求较高，并且实验成本较贵。数值模拟的成本相对较低，同时较物理模型实验更为灵活，因此近年来已成为相关领域的研究热点。

在实验室中，孤立波通常由造波机用推板造波[1-4]的方法产生，而数值模型中通常采用速度入口法或质量源方法造波[2-4]。基于动网格方法，采用类似推板造波的方式在数值模型中产生孤立波，并对网格尺寸对模拟结果的影响进行研究。

1 数值模型

流体的基本控制方程包括质量守恒方程和动量守恒方程。在研究波浪问题时，通常将水视为不可压缩的流体，即水的密度为常数，此时质量守恒方程可以简化为：

其中，U是水体的速度矢量，该方程通常称为连续性方程。

动量方程表示为：

其中Pd是动水压力，g表示重力加速度，X是位置矢量。

在数值模型中，水体和气体的交界面采用VOF方法进行计算，其控制方程为：

其中，α 为体积分数，代表水体占网格体积的百分数，当α=0 时代表气体，α=1 时代表水体，当网格在界面上时α介于0和1之间。

在空间上，对整个计算域用有限体积法进行离散。为了求解上述控制方程，将方程中的每一项进行离散，转化为代数方程，从而采用PISO 算法耦合求解速度和压力。用库朗数来限制时间步长，从而保证数值计算的精度，库朗数的表达式为：

屈哨兵：这个问题问得好。这也是我最近一直在思考的一个问题，就是好教育的质量观。这个问题我此前还没有做很多思考，最近稍微梳理了一下，大概有几个基本观点可以先说一下。

其中，Δt表示时间步长，δx是网格的长度，计算中的最大库朗数设为0.5，时间步长在计算中由于库朗数的限制不断变化。

2 动网格造波

采用动网格技术，可以在数值计算过程中，实现网格大小和位置的不断变化。对于造波边界处的网格，可模拟造波机的往复运动，从而用类似于推板造波的方式产生目标波浪。

对于推板造波，在造波边界上应当满足条件：

其中，u代表水质点的水平速度，在浅水条件下，由连续性方程可得水质点水平速度的关系如下：

其中，C表示波速，对于孤立波，有：

其中，H为孤立4波的波高，d为水深，将式(6)、(7)代入(5)，积分可得造波边界的位移公式：

造波边界的冲程为S，造波边界运动时间为T，

在数值模型中，当时间t＜T时，在每一时间步开始时迭代求解式(9)，求得的X(t)即当前时刻造波边界的位置，当造波结束后，造波边界在X(t)=S处静止不动。

3 数值计算验证

基于上述数值模型，建立数值波浪水槽进行验证，水槽长40 m，高1.7 m，如图1 所示，整个水槽均划分为长方形网格进行计算，为了保证数值计算精度，在水面附近对网格加密。水槽长度方向（x方向）的网格尺寸为Δx=0.25 m，在水面附近的网格加密区内，水槽高度方向（y方向）的网格尺寸为Δy=0.025 m。造波边界位于水槽左侧，模拟的孤立波波高为0.25 m，水深1.2 m。经计算，在水槽中x=5 m，10 m和15 m 处的波面历时曲线如图2 所示，其中x是与造波边界的距离。

图1 数值水槽示意图

由图2 可以看出，在距离造波边界5 m 处，模拟的孤立波波高与目标波高基本一致，说明动网格造波方法是有效的。但在孤立波继续传播的过程中，波高有所衰减，尤其是到达距离造波边界15 m 处时，衰减较为明显。

图2 距造波边界5 m，10 m和15 m处的波面历时曲线

将x方向上的网格加密，网格尺寸减小为Δx=0.08 m，y方向的网格仍为Δy=0.025 m，再次进行计算。此外，为了探讨y方向网格尺寸的影响，取网格尺寸Δx=0.08 m，Δy=0.065 m作为对照。由图2中的结果可以看出，在x=5 m和10 m处，由于距离造波边界较近，其沿程衰减并不明显，因此仅对比x=15 m处的计算结果。如图3所示，红色和蓝色虚线分别为网格尺寸Δx=0.08 m，Δy=0.025 m和Δx=0.08 m，Δy=0.065 m时，x=15 m 处的波面历时曲线，黑色线是该位置处的理论波面，用公式(7)计算。

图3 距造波边界5 m，10 m和15 m处的波面历时曲线

可以看出，当Δx=0.08 m，Δy=0.025 m 时，计算结果与理论波面几乎重合。与图2 对比可以看出，x方向的网格尺寸对孤立波的沿程衰减程度有一定影响，当x方向的网格足够密时，沿程衰减几乎可以忽略不计。通过与Δx=0.08 m，Δy=0.065 m 时的计算结果对比可以看出，当y方向的网格较疏时，不仅会造成孤立波的沿程衰减，还会造成孤立波相位的差异，这可能是由于此时水面处的速度计算存在误差，因此造成了孤立波的相位差。

需要指出的是，模拟孤立波时，网格尺寸的选择还与波浪参数有关。当波高较小时，网格尺寸也需要相应减小。例如当水深0.55 m 时，模拟波高为0.04 m 的孤立波，应当选取尺寸较小的网格，Δx=0.05 m，Δy=0.008 m，计算结果如图4 所示。可以看出，模拟的孤立波波高与目标波高一致，并且在传播至距离造波边界20 m 处时仍没有明显的衰减。此外，评判孤立波的造波效果时通常会观察是否存在尾波，即孤立波传播过后，是否存在后续的波动。在以上多个算例中均没有明显的尾波，说明孤立波的数值模拟方法是行之有效的。

图4 距造波边界10 m，15 m和20 m处的波面历时曲线

4 结语

基于动网格方法，用模拟推板造波的方式在数值模型中产生了孤立波。在数值计算中，沿水槽长度和高度方向的网格较疏时，会造成孤立波波高的沿程衰减，而水槽高度方向的网格较疏还会造成孤立波的相位变化。整体而言，产生的孤立波与目标波高一致，并且没有明显的尾波，说明文中的方法是有效的。