数值模拟推板式波浪水槽的关键参数取值探讨

关 超 吴静萍* 詹成胜 刘博宇

(武汉理工大学交通学院1) 武汉 430063) (中交机电工程局有限公司船舶及贸易事业部2) 北京 100088)

0 引 言

随着计算机技术和计算流体动力学(computational fluid dynamics，CFD)的发展，从20世纪60年代开始至今，计算流体动力学取得了丰硕的成果，出现了众多的CFD商业软件，如PHOENICS，CFX，FLUENT，STAR CCM+，FLOW-3D，FIDIP等，使得数值模拟方法以其成本低、周期短、精度高等优势在流体力学学科的科学研究和工程项目中已经得到认可和被推广使用[1-2].数值模拟方法在波浪运动、波浪与结构物相互作用问题中的应用，因自由面捕捉技术、造波技术、消波技术、动网格技术等关键技术起步相对较晚，与物理波浪水槽并称为数值波浪水槽，在海岸工程、港口与航道工程、船舶与海洋工程专业的波浪问题研究中起着重要作用.

数值波浪水槽造波有3种方法：仿物理(推板式或摇板式)造波法[3-4]、边界造波法[5]和源造波法[6-7].对于水波问题，属于非定常流动.推板式造波适用于浅水波及中等水深波，而摇板式造波则适用于深水波.数值模拟非定常的水波问题，计算相当耗时，除了必须采用高性能计算机和并行技术之外，其网格尺度和时间步长对减少计算时间和保证计算精度起着关键作用.

网格尺度、时间步长的取值问题一直受到研究者重视.郭晓宇[8]对源造波的数值波浪水槽进行验证时采用的波参数为波长L=2 m，波高H=0.05 m，水平方向网格尺度取L/20和L/40两种，垂向波高范围内网格尺度取H/8，结果显示两者与解析波形均吻合良好；此外又通过一系列具有不同非线性程度的波浪输入条件进行验证，波长分别为L=30，60和90 m，波高分别为H=1，3和4 m，水平方向网格尺度为L/120，L/180，L/240，发现随着波浪非线性程度的增大，单位波长需要较多网格节点才可以满足足够的精度.查晶晶等[9]用OpenFOAM实现数值波浪水槽造波，采用推板式造波法模拟波高H=0.02 m，周期T=1 s的线性波，所选网格尺度在30 m水池长度范围内的单元数为600，800，1 000，1 500，在垂向波高范围内的单元数为10，20，30，时间步长选取0.008，0.005，0.002，0.001 s进行探讨.通过与解析解比较，最终确定纵向网格、垂向网格和时间步长为1 500，20，0.002和800，10，0.002为最优参数，并得出结论，计算结果的精度主要受波面区域矩形单元的长宽比影响，网格单元越细长，引入的数值黏性越大，由此造出的波与目标波波幅相差太大.Finnegan等[10]采用CFX构建二维摇板式数值波浪水槽，模拟线性深水波与有限水深波，讨论网格依赖性，推荐了最小网格尺度计算式；推荐的时间步长为T/50.管陈等[11]采用FLUENT软件实现摇板式造波，模拟波长L=2 m，波高分别为H=0.04，0.06和0.1 m的规则波，选取RNG雷诺时均模型，时间步长约为0.001 s，自由液面附近网格最小尺度为波高的1/10～1/20，水平方向网格尺度沿池长方向逐渐加大，并将数值结果与解析解进行比较，发现同一波长下，波高越小，其误差越小，模拟精度越高.

此外，虽然一些学者研究数值波浪水槽时未对相关数值造波技术参数进行依赖性的讨论，但在文献中给出了其取值.严汝建等[12]使用FLUENT软件基于摇板式造波，数值模拟波参数为L=8 m，H=0.2 m的规则波，水平方向和自由液面附近垂向网格尺度分别为L/20和H/10，并选取RNGk-ε湍流模型.李宏伟[13]使用FLUENT建立二维数值水槽，采用摇板式造波法模拟参数为L=8 m，H=0.4 m，T=2.267 8 s的规则波，水槽长度方向上网格尺度为L/100，自由液面附近进行加密，时间步长为0.001 s，采用RNGk-ε湍流模型.黄华等[14]利用推板式造波法模拟二维二阶Stokes波，波浪参数为L=5.783 2 m，H=0.129 m，T=2 s，给出的网格尺度水平和竖直方向分别约为L/111和H/21，时间步长为T/400.廉静静等[15]采用FLUENT软件，通过摇板式造波模拟参数为L=8 m，T=2.284 8 s的规则波时，时间步长选为0.001 s，采用k-ε湍流模型.

本文以某流体力学实验室的小型推板式波浪水槽为原型，基于商业软件FLUENT建立推板式造波的二维数值波浪水槽.选取波浪水槽有代表性的波浪参数造波.考虑水波理论多数基于势流理论，而真实的流体具有黏性，所以假设本文模拟的波浪运动为层流流动.首先讨论了沿波长方向不同网格尺度、自由面附近沿波高方向上不同网格尺度以及不同时间步长对计算精度和计算耗时的影响，通过将数值结果与解析解的比较，综合考虑波高误差、单位距离波高相对衰减量和计算耗时，推荐综合性能较优的网格尺度、时间步长；然后将基于推荐的网格尺度、时间步长的数值结果与试验结果进行比较，两种吻合程度相当令人满意.

1 数值方法的基本控制方程

本文中二维波浪流动属于黏性、不可压缩和非定常类型，其波浪传播方向为x轴正向，波浪高度方向向上为z轴正向，则流体连续性方程为

(1)

动量方程为

(2)

(3)

式中：u为x方向速度；w为z方向速度；t为时间；ρ为流体密度；p为压强；ν为运动黏度；g为重力加速度.

采用VOF方法捕捉自由液面，控制方程为

(4)

式中：F为体积分数的函数，F=1为单元体全为水，F=0为单元体全为空气，0

2 数值计算方法

本文建立二维数值波浪水槽.水槽长为19 m(后4 m为消波段)，高为0.8 m，造波水深为0.5 m.生成目标波形的参数为L=1.5 m，T=0.995 1 s，H=0.06 m(为波峰致波谷之间的高度差)，属于二阶Stokes波.

采用建模软件Gambit进行数值波浪水槽的模型建立及网格划分，网格均为结构化网格.水槽底部设置边界层网格，自由液面附近进行网格加密，加密总高度为波高的2倍.数值波浪水槽与物理波浪水槽一样，需采要取一定的措施消除水槽末端池壁反射波的影响.数值消波技术主要采用人工阻尼区[16]、辐射边界条件[17]和主动吸收[18]，或者在数值波浪水槽末端使用稀疏网格进行数值耗散[19].本文在消波段采用渐变的稀疏网格耗散波能，从而消除水槽末端反射波的影响.

数值波浪水槽上方设为压力出口边界条件；两端及底部设为壁面边界条件，其中左边模拟物理推板造波，赋予推板一个UDF函数，采用layer的动网格方法，使其按照给定的方式进行运动.

数值波浪水槽模型和边界条件见图1，网格划分见图2.

图1 数值波浪水槽模型和边界条件

图2 数值波浪水槽网格划分

3 数值波浪水槽模拟计算

3.1 时间步长变化

时间步长的选取会影响数值计算结果的准确度和计算耗时.本文仅考虑波浪传播方向及自由液面附近网格尺度，首先取计算网格尺度Δx=L/70，Δz=H/10，考虑黏性，采用层流模型，取Δt=T/2 000，T/1 000，T/500三种时间步长进行数值模拟.

图3为t=19 s时沿池长方向的波形空间分布，z为波面抬高.表1为x=5 m与x=8 m处时间步长Δt变化时数值波高与解析波高比较的相对误差，以及单位距离波高相对衰减量.

图3 时间步长变化对数值模拟波浪传播的影响

Δt网格单元长宽比H1/m(x=5 m)波高相对误差|H1-H|H×100%H2/m(x=8 m)波高相对误差|H2-H|H×100%单位距离波高相对衰减量/%T/2 0003.570.054 88.670.053 710.500.61T/1 0003.570.056 46.000.056 26.330.11T/5003.570.057 14.830.058 03.330.50

由表1可知，随着时间步长的增大，数值结果与解析解的误差反而更小.由图3的波形及单位距离波高相对衰减量来看，当Δt=T/1 000时，波形的稳定性相对较好.综合考虑计算精度和波形稳定性，Δt=T/1 000是合适的时间步长选择.

3.2 网格尺度变化

网格尺度会对计算结果的精度产生影响，因此需对计算网格进行依赖性测试.首先讨论波长方向上网格尺度对数值造波的影响.取自由面附近波浪高度方向上网格尺度为Δz=H/10不变，波浪传播方向(水平方向)上网格尺度分别取Δx=L/50，L/60，L/70，L/80，L/90，L/100，L/110.

图4为在t=19 s时波形的空间分布.表2为x=5 m与x=8 m处网格尺度Δx变化时数值波高与解析波高比较的相对误差，以及单位距离波高相对衰减量.

图4 水平网格尺度变化对数值模拟波浪传播的影响

Δx网格单元长宽比H1/m(x=5 m)波高相对误差|H1-H|H×100%H2/m(x=8 m)波高相对误差|H2-H|H×100%单位距离波高相对衰减量/%L/505.000.054 49.330.053 910.170.28L/604.170.055 87.000.055 57.500.17L/703.570.056 46.000.056 26.330.11L/803.130.057 05.000.056 95.170.06L/902.780.057 24.670.057 14.830.05L/1002.500.057 44.330.057 34.500.06L/1102.270.057 54.170.057 44.330.05

由图4与表2可知，当波浪传播方向上网格尺度Δx=L/80时，数值结果与解析解的误差已比较小.

保持Δx=L/80不变，取波浪高度方向上、在自由面附近加密的网格尺度分别为Δz=H/10,H/20,H/30,H/40.图5为在t=19 s时波形的空间分布.表3为x=5 m与x=8 m处网格尺度Δz变化时数值波高与解析波高比较的相对误差，以及单位距离波高相对衰减量.

图5 竖直网格尺度变化对数值模拟波浪传播的影响

Δz网格单元长宽比H1/m(x=5 m)波高相对误差|H1-H|H×100%H2/m(x=8 m)波高相对误差|H2-H|H×100%单位距离波高相对衰减量/%H/103.130.05705.000.056 95.170.06H/206.250.056 85.330.056 46.000.22H/309.380.056 85.330.056 46.000.22H/4012.500.056 95.170.056 46.000.28

由图5和表3可知，波浪高度方向上网格尺度为Δz=H/10时数值计算结果基本趋近解析解.一味减小Δz并不能增加计算精度，要考虑合适的网格单元长宽比.

综合上述工作，本文波浪数值模拟的计算网格尺度为Δx=L/80，Δz=H/10已经能够达到较好的精度.

4 与物理试验结果比较

依据数值模拟推板式波浪水槽波浪传播的关键参数取值分析，所推荐的网格尺度Δx=L/80，Δz=H/10和时间步长Δt=T/1 000，选取层流模型，数值模拟波长为1.5 m、波高为0.06 m的波，将此数值计算结果与物理试验比较，进一步验证该数值算法的可靠性.物理试验水槽为武汉理工大学流体力学实验室小型推板式波浪水槽，长为18 m、宽为0.6 m、高为0.8 m、水深为0.5 m.试验采用伺服电机控制推板造波，在水槽x=5 m与x=10 m处各布置一个浪高仪，用以测量波高随时间的变化.

在数值计算中，同样在x=5 m与10 m处设置虚拟浪高仪，得出随时间变化的波形图.

图6为数值模拟结果与试验结果的比较.表4为x=5 m与10 m处数值平均波高与试验平均波高的比较及相对误差(图6方框内平均波高).

图6 数值结果与试验结果的比较

H1/m(x=5 m)H2/m(x=10 m)试验结果 0.057 20.055 9数值结果 0.056 50.055 3相对误差/%1.221.07

由图6与表4可知，数值模拟的波形与试验波形吻合程度很好，其变化趋势表现出良好的一致性，验证了该数值算法的有效性.

5 结 论

1) 本文数值波浪水槽中推荐采用：波浪传播方向网格尺度取Δx=L/80，自由液面附近波高方向网格加密尺度取Δz=H/10，时间步长取Δt=T/1 000；网格尺度确定还要满足一定的单元长宽比.

2) 数值模拟的结果与物理试验结果吻合程度很好，验证了本文数值算法的有效性.

在展开新型浮式防波堤衰减长周期波性能研究时遭遇物理波浪水槽长度尺度限制的问题.本文工作为发挥数值波浪水槽的优势，利用数值模拟手段辅助物理试验进行研究奠定了基础.也为读者应用数值波浪水槽进行研究提供有益参考.

参考文献

[1] 王福军.计算流体动力学分析:CFD软件原理与应用[M].北京:清华大学出版社，2004.

[2] 韩占忠.FLUENT：流体工程仿真计算实例与应用[M].北京:北京理工大学出版社，2010.

[4] PRASAD D D, AHMED M R, LEE Y H, et al. Validation of a piston type wave-maker using Numerical Wave Tank [J]. Ocean Engineering, 2017,131:57-67.

[5] 刘霞，谭国焕，王大国.基于边界造波法的二阶Stokes波的数值生成[J].辽宁工程技术大学学报(自然科学版)，2010，29(1):107-111.

[6] LIN P, LIU P, L F. Internal wave-maker for navier-stokes equations models[J]. Journal of Waterway, Port, Coastal and Ocean Engineering, 1999,125(4):207-217.

[7] WEI G, KIRBY J T, SINHA A. Generation of waves in boussinesq models using a source function method[J]. Coastal Engineering, 1999,36(4):271-299.

[8] 郭晓宇.数值波浪水槽及其应用研究[D].上海：上海交通大学，2011.

[9] 查晶晶，万德成.用OpenFOAM实现数值水池造波和消波[J].海洋工程,2011,29(3):1-12.

[10] FINNEGAN W, GOGGINS J. Numerical simulation of linear water waves and wave-structure interaction [J]. Ocean Engineering, 2012,43(4):23-31.

[11] 管陈，董国祥，金允龙.三维数值波浪水池造波技术研究[J].上海船舶运输科学研究所学报,2013,36(2):11-15.

[12] 严汝建，庞永杰，李宏伟，等.深水池造波系统数值造波仿真研究[J].哈尔滨工程大学学报,2010,31(1):32-41.

[13] 李宏伟.数值水池造波方法研究[D].哈尔滨：哈尔滨工程大学，2011.

[14] 黄华，邓冰，陈昱松，等.数值波浪水槽构建与二阶Stokes波仿真[J].系统仿真学报,2012,24(1):227-231.

[15] 廉静静，尹勇，杨晓，等.基于黏性流船舶数值波浪水池造波和消波方法研究[J].船舶力学,2013(1):56-62.

[16] 赵西增，孙昭晨，梁书秀.非线性波浪水槽内波浪的产生和传播[J].船舶力学,2010,14(3):203-216.

[17] ORLANSKI I. A simple boundary condition for unbounded hyperbolic flows[J]. J Comput Phys, 1976,21(3):251-269.

[18] DEEPAK D P,MOHAMMED R A,YOUNG H L,et al. Validation of a piston type wave-maker using Numerical Wave Tank [J]. Ocean Engineering, 2017,131:57-67.

[19] PARK J C, UNO Y, SATO T, et al. Numerical reproduction of fully nonlinear multi-directional waves by a viscous 3D numerical wave tank[J]. Ocean Engineering, 2004,31(11):1549-1565.