本发明涉及海岸工程技术领域，尤其涉及一种海啸波越浪过程的二维数值模拟方法。
背景技术：
海堤作为重要海岸防护建筑物之一，在波浪潮流等动力荷载的作用，时常被破坏，尤其是在极端的动力条件下，如海啸、风暴潮等。溃堤导致沿海地区直接承受波流的影响，给沿海地区的人类生活和生产带来了巨大的灾难。海啸，是一种具有强大破坏力的海浪，以海底地震、滑坡变动处形成向四周传播的表面波，在向近岸传播的过程中几乎无能量损失，随地形的变化发生浅水变形，对近岸地区构成巨大的威胁。
根据破坏位置的不同，波浪对海堤的破坏可以分为以下三个类型：(1)堤脚冲刷和前坡失稳；(2)堤顶冲刷破坏；(3)后坡冲刷侵蚀破坏。越浪流是造成海堤破坏的重要因素之一，而越浪流对海堤后坡的破坏体现在以下三个方面：(1)对堤顶及后坡的冲刷破坏；(2)对后坡的冲刷导致海堤发生滑坡；(3)使海堤土体含水量增大导致海堤失稳。
目前，关于越浪流对海堤的影响的研究主要是针对平均越浪量的研究。越浪流的研究方法主要有三种：(1)理论研究；(2)物理实验；(3)数值模拟，而物理实验是最为常用的方法，但是物理实验的研究成本相对较高。相比于之下，数值模拟更加的经济适用。受自由表面处理技术的影响，研究越浪流的数值模型还不是很成熟，需要进一步地去研究和发展。
技术实现要素：
有鉴于现有技术的上述缺陷，本发明所要解决的技术问题是提供一种海啸波越浪过程的二维数值模拟方法，可以模拟完整的海啸波的越浪过程以及流速场的变化过程、计算越浪量和堤顶越浪流流速、得到堤身受到的最大冲击压力点及其位置。
本发明的一种海啸波越浪过程的二维数值模拟方法，包括以下步骤：
步骤s1、根据海堤结构轮廓和水深条件，确定海堤模型尺寸、数值水槽尺寸以及所有粒子的初始空间坐标；
步骤s2、根据目标海啸波波高，迭代生成采用推波板造波法所需的初始波高值；
步骤s3、通过直接离散的拉格朗日形式的连续方程和动量方程，控制流体的运动，模拟海啸波的越浪过程，得到波浪的运动形态、水流结构的速度场和压力场；
步骤s4、对输出的结果进行计算，得到最大冲击压力值以及其发生的位置。
其中，上述步骤s1具体为：
(1-1)根据实际海堤结构形态，提取海堤轮廓线，确定数值水槽中的海堤
结构形式；
(1-2)选取计算水深和海啸波的波浪要素，即目标波高度h，并结合海堤
结构形式，确定水槽长宽尺度，即确定数值模型的计算域范围；
(1-3)根据计算域大小，确定计算粒子直径d；
(1-4)根据海堤轮廓线，确定固定墙粒子的布置区域；
(1-5)根据初始粒子数密度相等条件，通过牛顿迭代法计算，确定墙粒子
之间的初始距离，确定墙粒子的初始空间坐标，粒子数密度相等条件表示为
n0|i'＝n0|i
式中：n0为初始粒子数密度；i、i′为任意的墙粒子；并且，
式中：i为目标粒子；j为粒子i的影响域内的粒子；ni为粒子数密度；w为核函数；re为影响域半径；rij为粒子之间的距离；
(1-6)在水槽的最左端设置可移动的墙粒子，构成推波板；
(1-7)根据初始的水深h，以每个固定的墙粒子为中心，从上至下依次布
置流体粒子，粒子间的垂直间距以粒子直径为标准，确定流体粒子的初始
空间坐标。
上述粒子直径d选为0.005m。
上述步骤2具体包括：
(2-1)构建一个平坦水槽，输入一个任意的波浪高度h0，通常取目标海
啸波波高h；
(2-2)根据输入的目标海啸波波高h，控制推波板运动，生成入射波；
(2-3)计算生成的波浪在设定位置处的波高h，并计算该波浪高度值与目标海啸波波高h的相对误差re＝|h/h-1|，若re<ε，ε为允许误差，则该波浪高度值即为控制推波板运动的初始波高h0，反之，则对波浪高度h0进行迭代修正：
式中：h为目标海啸波波高值；h′为设定位置处生成波波高值；h0为迭代计算初始波高值；γ为松弛系数，大于1；
迭代该式，直至计算结果满足误差要求。
上述步骤3具体包括：
(3-1)输入步骤s1中确定的粒子初始空间坐标和计算域的范围；
(3-2)输入推波板运动的初始波高值h0；
(3-3)输入模型基本参数，模型基本参数包括流体密度ρ和运动粘滞系
数ν；
(3-4)根据cfl常数计算时间步长；
式中：δt为时间步长；δtmax为模型计算容许的最大时间步长；δtc为cfl常数条件允许的最大时间步长；δtν为运动粘滞项允许的最大时间步长；|ux|max为最大水平速度值；|uy|max为最大垂向速度值；l0为粒子直径值；
(3-5)采用活塞式推波板造波边界，根据推波板运动方程，控制可移动墙粒子的位移，推波板运动方程为
式中：xw为推波板位移；s为推波板的最大冲程；t为推波板运动时长；h0为初始波高；c为波浪传播速度；h为初始水深；n为计算参数；
(3-6)利用连续方程和动量方程控制流体粒子的运动，拉格朗日形式的控制方程可表示为
式中：ρ为流体密度；u为粒子的速度矢量；p为压力值；g为重力加速度；ν为运动粘滞系数；
采用映射法对控制方程进行求解，每步计算过程分为两步：第一步先计算动量方程中的重力项和粘滞力项，得到临时的速度场和位移场
式中：ui*为临时速度矢量；ri*为临时位移矢量；uik为k步长时刻的速度矢量；rik为k步长时刻的位移矢量；δt为时间步长；
第二步将临时的速度场和位移场带入连续方程中得到一个压力泊松方程，通过求解泊松方程得到压力场，再将压力梯度项
式中：uik+1为k+1步长时刻的速度矢量；rik+1为k+1步长时刻的位移矢量；ρ0为流体初始密度；
(3-7)更新计算时刻；
(3-8)每固定时间间隔输出模拟结果，模拟结果包括所有粒子的空间位置、速度和压强。
上述步骤4具体包括：
(4-1)绘制完整的海啸波越浪运动过程图，得到海啸波在海堤上的运动形态；
(4-2)绘制海啸波越浪过程的流场图、速度场图和压力场图；
(4-3)确定越浪量计算断面，根据堤顶上方自由表面与堤顶间的垂直距离，计算该断面处的堤顶越浪流厚度，并沿断面方向平均布置流速计算点，采用加权平均法计算堤顶越浪流断面平均流速，通过堤顶越浪流厚度乘以越浪流断面平均流速得到越浪量；
式中：q为堤顶越浪量；vc(t)为t时刻时的堤顶越浪流平均流速；hc(t)为t时刻时的堤顶越浪流厚度；t1为越浪量计算的起始时刻；t2为越浪量计算的终止时刻；
(4-4)根据压力场的输出结果，以每个墙粒子为中心，计算其1.5d空间范围内的平均压力值提取墙粒子的压力值，结合海堤的位置，得到海堤所受的压力，从而确定最大冲击压力值以及其发生的位置。
本发明的有益效果是：
本发明提供的模拟海啸波越浪的二维数值计算方法，利用粒子法的拉格朗日运动系统的特点，很好地重现海啸波的完整越浪过程，并且有效地计算了海啸波的越浪量、越浪流速以及海堤上可能出现的冲击压破坏位置，为越浪流研究提供了新的方法，同时该发明本身也促进了粒子法数值模型在波浪越浪领域的发展。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明的实施例的流程示意图；
图2是本发明实施例的数值水槽初始粒子分布图，不同颜色的粒子代表不同的粒子类型；
图3是本发明实施例中计算得到的海啸波越浪过程图；
图4(a)是堤顶越浪流层流水深hc的历时曲线图；
图4(b)是越浪流断面平均流速vc的历时曲线图；
图4(c)是堤顶越浪量q的历时曲线图。
具体实施方式
如图1所示，一种海啸波越浪过程的二维数值模拟方法，包括以下步骤：
步骤s1、根据具体的海堤结构轮廓和水深条件，确定海堤模型尺寸、数值水槽尺寸以及所有粒子的初始空间坐标；
步骤s2、根据目标海啸波波高，迭代生成推波板造波所需的初始波高值；
步骤s3、通过直接离散的拉格朗日形式的连续方程和动量方程，控制流体的运动，模拟海啸波的越浪过程，得到波浪的运动形态、水流结构的速度场和压力场；
步骤s4、对输出的结果进行后处理。
本实施例中，步骤s1具体为：
(1-1)根据实际海堤轮廓，确定模型中海堤结构形式，海堤堤身高0.2m，
前坡坡度1:6，后坡坡度1:3；
(1-2)根据海堤结构大小，确定海啸波的波浪要素，目标波高h＝0.072m，确定水槽长宽尺度4.85m×0.6m，即确定计算域范围；
(1-3)根据计算域尺度，粒子直径d选为0.005m；
(1-4)提取海堤轮廓线，确定固定墙粒子布置区域；
(1-5)根据初始粒子数密度相等条件，通过迭代计算，确定墙粒子之间的
初始距离，粒子数密度相等条件表示为
n0|i'＝n0|i
并且，
式中：n——粒子数密度；
w——核函数，一种权函数；
re——影响域半径；
rij——粒子之间的距离。
根据粒子的空间分布，利用核函数计算得到每一时刻的粒子数密度。
(1-6)在水槽的最左端设置可移动的墙粒子，构成推波板；
(1-7)根据初始的水深h＝0.18m，以每个固定的墙粒子为中心，从上至下
依次布置流体粒子，粒子间的垂直间距以粒子直径为标准。
本实施例中粒子的初始分布如图1所示。本实施例中，所述步骤2具体包括：
(2-1)构建一个平坦水槽，输入一个任意的波浪高度h0，通常取目标海
啸波波高h，即h0＝0.072m；
(2-2)根据输入的波高值，控制推波板运动，生成入射波；
(2-3)计算生成波浪在指定位置处的波面高度h，并计算该波高値与
目标波高的相对误差re＝|h/h-1|，若re<ε，ε＝0.01，则该波高值即
为控制推波板运动的初始波高h0，反之，则对波高h0进行迭代修正
式中：h——目标海啸波波高值；
h′——设定位置处生成波波高值；
h0——迭代计算初始波高值；
γ——松弛系数，γ＝1.05。
通过实际生成的波浪高度h，与目标波高h的反比，对初始设定的波浪高度进行修正，利用迭代收敛得到控制推波板运动的初始波高值h0。
本实施例中，所述步骤3具体包括：
(3-1)输入步骤s1中获得的粒子初始空间坐标和计算域的范围；
(3-2)输入推波板运动的初始波高值h0；
(3-3)输入模型基本参数，主要包括：流体初始密度ρ0＝1000kg/m3和运动粘滞系数ν＝10-6m2/s；
(3-4)计算时间步长；
式中：δt——时间步长；
δtmax——模型计算容许的最大时间步长，δtmax＝0.001s；
δtc——cfl常数条件允许的最大时间步长；
δtν——运动粘滞项允许的最大时间步长；
|ux|max——最大水平速度值；
|uy|max——最大垂向速度值；
l0——粒子直径值，l0＝d＝0.005m。
计算出cfl常数条件允许的最大时间步长、运动粘滞项允许的最大时间步长以及数值模型允许的固定最大时间步长，取其最小值作为实际选用的时间步长。
(3-5)利用推波板运动方程，控制可移动墙粒子的运动，在造波边界处生成波浪，推波板的运动方程为
式中：xw——推波板位移；
s——推波板的最大冲程；
t——推波板运动时长；
h0——初始波高；
c——波浪传播速度；
h——初始水深，h＝0.18m；
n——计算参数。
通过控制每一时刻推波板的位置，得到推波板的运动速度，进而推动水体运动，生成波浪。
(3-6)利用连续方程和动量方程控制流体粒子的运动，拉格朗日形式的控制方程可表示为
式中：ρ——流体密度；
u——粒子的速度矢量；
p——压力值；
g——重力加速度，g＝9.81m/s2；
ν——运动粘滞系数，ν＝10-6m2/s。
整个过程分为两步：第一步先计算动量方程的重力项和粘滞力项得到临时的速度场和位移场
式中：ui*——临时速度矢量；
ri*——临时位移矢量；
uik——k步长时刻的速度矢量；
rik——k步长时刻的位移矢量；
δt——时间步长。
第二步将临时的速度场和位移场带入连续方程中可得到一个压力泊松方程，通过求解泊松方程得到压力场，再将压力梯度项带入动量方程中，更新得到下一时步的速度场和位移场
式中：uik+1——k+1步长时刻的速度矢量；
rik+1——k+1步长时刻的位移矢量；
ρ0——流体初始密度，ρ0＝1000kg/m3。
(3-7)更新计算时刻；
(3-8)每固定时间间隔输出模拟结果，主要包括所有粒子的空间位置、
速度和压强。
本实施例中，所述步骤4具体包括：
(4-1)绘制完整的海啸波越浪过程图，及其速度场图和压力场图；
(4-2)计算堤顶越浪流厚度和堤顶越浪流流速，通过堤顶越浪流厚度乘以
越浪流断面平均流速得到越浪量；
式中：q——堤顶越浪量；
vc(t)——t时刻时的堤顶越浪流平均流速；
hc(t)——t时刻时的堤顶越浪流厚度；
t1——越浪量计算的起始时刻；
t2——越浪量计算的终止时刻。
(4-3)计算海堤受到的压力值，得到最大冲击压力值以及其发生的位置。
图2为本发明实施例的数值水槽初始粒子分布图，不同颜色的粒子代表不同的粒子类型。
图3给出了本实施例得到的海啸波越浪过程，从图3可以确定该海啸波的越浪运动形态图，波浪回落的过程中发生了水跃现象；图4(a)至(c)给出了本实施例得到的堤顶越浪流层流水深hc、越浪流断面平均流速vc和堤顶越浪量q的历时曲线图，从图4可以得到该海啸波的越浪量q＝0.021m3。
具体地，本发明所述的基于粒子法的海啸波越浪的二维数值方法包括：
根据实际海堤结构轮廓，提取固壁边界的形状，确定海堤在数值模型里的结构形式；
根据水深、波浪条件以及海堤大小，选择合理的粒子直径d和粒子的影响半径re；
根据初始粒子数密度相等条件，确定墙粒子之间的间距，确定固定墙粒子的初始位置，并且以墙粒子中心线为基准，向外扩展布置虚粒子层，弥补核函数的截断误差；
采用活塞式推波板造波法生成目标波浪，在水槽的最左端生成可移动墙边界；
根据初始水深条件，以墙粒子为基点，从下到上依次布置流体粒子，粒子间距取与粒子直径相等；
根据推波板造波原理和海堤的位置点，在平坦水槽的相同位置设置波高计算点，选用合适的位移方程控制推波板的运动，通过迭代尝试多组波浪参数，确定初始波浪要素，即推波板运动方程中的初始波高。选用一阶精度的孤立波波面方程代表海啸波，其推波板的控制方程选择显示的运动方程，推波板位移xw表示为
式中，s为推波板的最大冲程，t为推波板的运动时长，两者分别表示为
式中：h0——初始波高；
c——波浪传播速度；
h——初始水深；
n——计算参数。
根据推波板每一时刻的位置，推动推波板运动，从而生成海啸波。
计算在指定位置处的生成波浪的波高h’是否满足目标波高h的要求，即计算该波高値与目标波高的相对误差re＝|h/h-1|是否满足re<ε，如果满足，确定推波板运动方程中的波浪高度h0，否则更新初设波高h0，重新生成波浪，直到满足波高要求,波高的迭代公式选取为：
基于移动粒子半隐式法(mps法)建立一个二维的数值波浪水槽，计算每个时刻的波浪场的水动力条件，具体包括：
采用新的自由表面判别方法，以粒子数密度ni和目标粒子的影响域内粒子的填充率ai双重标准来识别出自由表面粒子，自由表面粒子必须满足条件：
ai＜αa0&ni＜βn0(α，β＜1.0)
式中：n0——初始粒子数密度；
a0——初始粒子间线性距离和；
ni——粒子数密度；
ai——粒子间线性距离和；
w——核函数；
α，β——松弛系数，α＝0.9，β＝0.97；
该公式利用自由表面粒子的粒子数密度ni和粒子间线性距离和ai均小于其各自的初始值的特点来识别自由表面粒子。
采用拉格朗日形式的连续方程和纳维-斯托克斯(navier-stokes)方程作为控制方程，追踪每个流体粒子的运动过程，其控制方程形式如下：
根据粒子间的相互作用，得到各个算子的离散形式，其梯度算子模型、散度算子模型以及laplace算子模型的表示分别如下所示：
式中：n0——初始粒子数密度；
φ——任意物理标量；
ds——空间维度，二维空间取2；
w——核函数；
λ——计算参数，表示为：
n0是初始粒子数密度，物理意义等价于流体密度，n取决于核函数的形式，通常表示为：
核函数选用kondo和koshizuka于2011年提出的标准核函数，其形式为
为了提高原始的mps方法的计算精度，采用一个新的压力梯度算子模型，根据泰勒级数展开，新的压力梯度模型可以表示为
式中：pj——粒子j的压力值；
ci——计算参数，表示如下：
通过压力梯度模型，可以得到由压力场生成的加速度场。
将各个算子模型带入到控制方程中，离散后的控制方程表示为
式中：uik+1——k+1步长时刻的速度矢量；
uik——k步长时刻的速度矢量；
ρ0——流体初始密度，ρ0＝1000kg/m3。
离散后的方程式为最终控制模型运行的公式。
方程组中的第二等式又被称为压力泊松方程。为了提高计算结果的精度，采用kondo和koshizuka于2011年提出的改进的多项源项的压力泊松方程
式中，b和γ是经验系数，通常取值500和50,000。
采用分步投影法，第一步忽略动量方程中的压力项，只计算粘滞项和重力项，得到一个临时的速度场，并更新所有流体粒子的位移，
式中：ui*——临时速度矢量：
ri*——临时位移矢量
第二步，将得到的临时的位移场带入压力泊松方程，采用双共轭梯度迭代法求解得到下一时刻的压力场，将新的压力值带入动量方程中的压力梯度项，更新第一步中得到的临时速度场，得到下一时步的速度场和位移场，具体表示如下
式中：uik+1——k+1步长时刻的速度矢量；
rik+1——k+1步长时刻的位移矢量；
ρ0——流体初始密度。
该方程利用压力梯度项对临时速度场和位移场进行校正，可以得到下一时刻的准确的速度场和位移场。根据最新的速度场和柯朗-弗里德里希斯-列维条件(cfl条件)，计算下一时步的时间步长δt
δt＝min(δtc,δtν,δtmax)
其中，
在计算得到的所有时间步长值中选取其最小值，从而满足所有的数值稳定条件的要求。
在堤顶最前端xs处设置判别函数f，若存在流体粒子的x坐标值大于xs，则f＝1，表明存在越浪流，若无流体粒子的x坐标值大于xs，则f＝0，表明无越浪发生。针对f＝1情况，采用schüttrumpf越浪量计算方法，在海堤前坡与堤顶交界处、海堤堤顶中部、海堤后坡与堤顶交界处等三处位置设置水位、流速计算点，每间隔固定时间，计算堤顶越浪流水深hc和断面平均流速vc。越浪量q的计算公式可以表示为
式中：vc(t)——t时刻时的堤顶越浪流平均流速；
hc(t)——t时刻时的堤顶越浪流厚度；
t1——越浪量计算的起始时刻；
t2——越浪量计算的终止时刻。
每间隔固定时间，提取构成海堤的墙粒子的压力值，取每个墙粒子在x轴上的投影为横坐标，压力值为垂向坐标，时间为纵坐标，绘制海堤断面压力图，计算可能出现的波浪冲击压的作用位置和最大冲击压大小。
本发明与现有技术相比，其显著优点是：本发明提供的模拟海啸波越浪的二维数值计算方法，利用粒子法的拉格朗日运动系统的特点，很好地重现海啸波的完整的越浪过程，并且有效地计算了海啸波的越浪量、越浪流速以及海堤上可能出现的冲击压破坏位置，为越浪流研究提供了新的方法，同时该发明本身也促进了粒子法数值模型在波浪越浪领域的发展。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。