波浪对直立墙作用力数值模拟方法研究

摘要： 本文基于GPU-SPH方法建立了数值波浪水槽。通过加入Shepard密度过滤器和Hugehes & Graham修正算法，改进数值模型压力测定方法。计算立波对直立墙的作用力之后，将数值结果与理论计算结果进行对比、分析，前者与后者符合较好。

关键词：GPU-SPH；数值波浪水槽；波浪对直立墙作用压力；

中图分类号：S969.1 文献标识码：A

引言

近年来在计算流体力学领域发展起来一种新型无网格拉格朗日数值计算方法——光滑粒子流体动力学算法（SPH）。由于原SPH方法为防止粒子的边界穿透现象而让边界对粒子的作用力偏大，在边界附近粒子压力较大，影响计算精度和压力测量结果。为达到边界压力计算准确的目的，本研究在模型计算中加入Shepard【[Panizzo A (2004) Physical and Numerical Modelling of Subaerial Landslide Generated Waves. PhD thesis, Universita degli Studi di L'Aquila. ]】密度过滤器过滤因压缩而密度过大的粒子，同时使用Hugehes & Graham【[Hughes, J. and Graham, D., Comparison of incompressible and eakly-compressible SPH models fro free-surface water flows, J. Hyd. Res., in press, 2010.]】算法对边界进行密度进行修正来减少计算运行后的压力波动现象，并且选择了适用于波浪对结构物作用模拟的参数和波浪作用力测量方法。

基于本文模型，建立弱可压缩SPH-GPU算法建立数值波浪水槽，采用Wendland 5次光滑函数，研究波浪与直立墙相互作用，并将模拟结果与理论值和前人实验数据进行比较，分析了规则波作用下直立墙前压力场特性。

SPH方法数值模拟

SPH原理

SPH方法中，核插值近似法将场变量的积分形式表达如下：

(1)

其中：r为粒子的矢径；h是光滑长度，积分区域的大小由其控制；W为核函数，本文采用Wendland提出的五次型函数，其表达式为

(2)

对于二维的计算情况，式中为，。

控制方程

对于水槽中的流体，流体运动的N-S方程：

(3) du/dt=-1/ρ ∇P+g+θ 

(4) dρ/dt=-ρ(∂v^β)/(∂x^β )

其中：为水的密度，为重力加速度；为压力；为动量方程中的扩散项，针对本模型，使用SPH方法中最常用的人工粘度法计算扩散项 。

边界条件及压力计算方法调整

本文采用动态边界条件来处理边界粒子，包括水槽固壁边界、造波板边界、以及结构物边界。边界粒子以一种错列的形式安排如图1所示：

图1 边界粒子排列形式图

当一个流体粒子接近边界时，边界粒子的密度会增加以至于压力也随其增加。当距离流体粒子和边界粒子之间的距离小于2h时，施加在入射流体粒子上的压力作为排斥力随着距离减小而增加如图2所示。

图2 边界密度距离、压力距离曲线图

图中h为粒子的光滑长度，也叫做影响域或者光滑域，它的大小控制着每个粒子多其周围其他粒子的作用力大小。h的计算公式如下：

(5)   h=α*√(dx*dx+dy*dy+dz*dz)  

                  

本文中设置参数，2维坐标轴的微变量dx、dz都为0.01m。因此计算得为。而一般情况下，了防止粒子穿透边界层现象的发生，一般边界作用力是远大于现实的水压力的。因此，如果将水流对结构物作用压力测点直接设置在其边界粒子上得到的压力值往往远大于理论值。在距离边界h处边界粒子对此处粒子的直接作用力便可忽略不计。但是考虑到水粒子与边界粒子在距离h之间时，其排斥力影响水粒子后由此粒子将仍会将较大的压力传递至距离边界h之外的粒子上，选择将测点设置在距离边界粒子2h与1h之间的某处位置。这样能更加准确地测出直立墙附近波浪压力。本文选取此距离为0.002m。

计算参数选取

本文中所有计算参数根据汪立军【[汪立军，郑永来，基于SPH的数值波流水槽系统及验证，实验室研究与探索[J]，324-328，2013.]】的研究成果设置，具体参数如下表：

表1 参数设置表

数 参数值

重力加速度（N/kg）

粒子直径（m）

CLF值

声速率

γ

rhop0值

eps压缩率值 9.8

0.01

0.2

10

0.92

7

1000

0.3

Shepard密度过滤器

当SPH的动态模拟越来越真实的同时，粒子压力值出现很大的压力波动。前人致力于可以通过粒子修正或发展一种不可压缩解等方法来解决。解决这个问题的方法包括粒子修正或发展一种不可压缩解等。其中一种最直接且计算成本最低的方法就是提供一个粒子密度过滤器，同时重置每个粒子的密度（Colagrossi 和 Landrini）。Shepard密度修正法，也称作零阶修正法，是一种快捷简单的密度过滤器，以下计算过程每30个时间步长进行一次：

在已经被修正过的粒子上使用一个零阶修正：

波浪对直立墙作用数值模拟

波浪对直立墙作用力的数值模拟

《海港水文规范》【[海洋水文规范 JTJ 213-98[M]. 北京：人民交通出版社.1998]】中规定的方法是目前计算波浪对直立堤作用力的一种比较普遍的方法，本文通过SPH模拟规则波正向入射对直立墙的作用，并且将结果与计算方法结果及前人试验结果进行比较。

水槽布置

水槽布置如图3所示，水槽长11.00m,模拟水深0.33m，波高为0.05m，波长为1.04m，周期为0.8s。直立墙设置在距离左边界10m处，造波板位于水槽最左侧。左侧为造波边界，底部、右侧和直立墙为固边界。粒子间距取0.01m，整个水槽共布置140460个粒子。根据本文上部分的研究成果，在距离直立墙0.02m处设置波压力测点布置线，沿布置线每间隔0.0025m设置一压力测点，共321个测点，检测波浪对直立墙作用力随时间及深度的变化情况。

图3 立波与直立墙作用状况模拟水槽布置图

造波边界

本文造波板运动周期即为波浪周期。需要波高根据多次改变造波板运动位移，凑谱得到。

模拟结果及分析

取波峰作用时刻进行分析。将波峰对直立墙作用压力计算结果与《海港水文规范》计算结果进行比较。根据杨成渝[杨成渝，波浪正向入射对直立堤作用研究，第十四届中国海洋（岸）工程学术讨论会论文集，609-613，2009.]的试验结果拟合曲线计算最大相对波压力值为431pa，模拟值为550pa，误差率为27%，根据规范计算最大相对波压力为500pa，误差率为10%。图4给出了SPH模拟立波波峰作用下墙面上的压力分布与规范计算值的比较。综上所述，本文数值模拟与计算值和试验吻合较好，较为精确地反映了波浪对直立墙作用的力学特性。

图4 海港水文规范计算值与数值模拟值波峰作用压强深度曲线比较

结语

基于SPH方法的基本原理，加入Shepard密度过滤器和Hugehes & Graham算法对压力求解进行修正改进。同时通过理论分析，提出了一种新的边界压力测量计算方法，解决了模拟边界压力远大于真实值的问题。并且通过静水压力测定对新方法进行了初步验证。模拟了立波对直立墙的作用，数值模拟结果与规范计算结果和前人试验结果进行对比分析，吻合较好，验证该模型具有较为准确地模拟水动力学问题的能力。