三维数值波浪水池的构建和粘性影响研究

赵 艳，朱仁庆，刘 珍

(江苏科技大学，江苏 镇江 212003）

三维数值波浪水池的构建和粘性影响研究

赵 艳，朱仁庆，刘 珍

(江苏科技大学，江苏 镇江 212003）

以商业软件Ansys14.0中的CFD模块Fluent为平台，通过程序接口加载UDF程序进行二次开发，采用流体体积法VOF(Volume of Fraction）和几何重构法(Geo-Reconstruct）捕捉流体分界面，利用动网格铺层技术Layering实现造波板的运动，同时在动量方程中添加源项实现阻尼消波，建立具有良好造波、消波功能的三维推板式数值波浪水池。利用本数值水池，结合理想流体Inviscid模型和湍流RNG k－ε模型，通过模拟具有不同波长、波高的余弦波研究粘性在数值波浪水池中的影响。结果分析表明，粘性力对流体质点的运动存在阻滞作用，从而使波高沿水池长度方向发生了明显的衰减，衰减幅度与波陡呈正比，衰减形式符合指数规律。

粘性力;衰减规律;数值波浪水池;VOF

0 引言

近几十年来，随着计算机性能的不断提升和数值模拟技术的快速发展，数值波浪水池在水动力研究领域的应用越来越广泛。与传统的物理实验波浪水池相比，其成本更低，实验周期更短，对流体质点的监测和观察更方便易行，而且不需要耗费人力物力进行日常维护，因此被越来越多的学者所关注，并取得一定的研究成果。

谷汉斌等［1］研究了基于完全非线性Boussinesq方程的源函数数值造波。王大国等［2］用有限元求解拉普拉斯方程，以半拉格朗日法跟踪流体自由表面，建立了三维完全非线性数值波浪水池。宁德志等［3］基于无旋、不可压的势流理论，利用高阶边界元法建立了一种可应用于无限水深的无粘三维完全非线性数值波浪水槽，继而在此基础上建立了波流混合作用的水槽模型。William Finnegan等［4］以 CFD软件Ansys CFX为平台，建立了基于N-S方程的摇板式数值水池，并分析了流体质点的速度和水池的适用范围。Guillaume Ducrozet等［5］基于一种改进的高阶谱方法，建立了三维数值波浪水池，并对不规则波和聚焦波进行了模拟研究。

本文以商业软件Ansys14.0中的CFD模块Fluent为平台，通过程序接口加载UDF程序进行二次开发，建立具有良好造波、消波功能的三维推板式数值波浪水池。并结合理想流体Inviscid模型和湍流RNG k－ε模型，通过模拟具有不同波长、波高的余弦波研究了粘性在数值波浪水池中的影响。

1 数值波浪水池的数学模型

1.1 控制方程

本文数值波浪水池的控制方程有2种形式:当考虑粘性力时，控制方程由连续性方程(1）和N-S方程(2）构成;当忽略粘性力影响时，控制方程由连续性方程(1）和 Euler方程(3）构成［6］。

式中:t为时间;ρ为流体密度;v为速度矢量;μ为动力粘性系数;Fb为质量力;p为压强;S为应变率张量。式(2）左端代表单位体积流体的惯性力;右端第1项代表单位体积的质量力;第2项代表作用于单位体积流体的压强梯度力;第3项代表粘性变形应力;第4项代表粘性体膨胀应力。

由于所建立的数值水池中只有空气和水2种密度不相容的流体，所以适合采用流体体积法VOF(Volume of Fraction）来追踪2层流体的交界面。这种方法不仅计算时间短、存储量少，而且可以跟踪复杂变形的自由面，如在自由面上发生的翻转、运并、飞溅等强非线性现象。其基本原理是通过网格单元的流体体积分数αq来确定流体界面，即αq定义为单元内第q相流体所占有体积与该单元的体积之比。若αq=1，表示该单元内全部为第q相流体;若αq=0，表示该单元内没有第q相流体;若0＜αq＜1，表示该单元为交界面单元［7］。

各项流体的体积分数满足:

式中:u为x方向速度分量;v为y方向速度分量;w为z方向速度分量。

每个单元中流体的各项属性都由单元中的所有分相决定，例如单元中流体的密度:

如果仅采用VOF法捕捉自由面，得到的是一个不光滑的模糊界面，因此，还需要结合几何重构方法 (Geo-Reconstruct）进行界面重构，以获得精度较高的光滑自由液面。

1.2 Inviscid模型和RNG k－ε模型

Inviscid模型是一种描述理想流体运动的数值模型。它假设所有流体没有粘性，在流动时各层之间没有相互作用的切应力，对切向变形没有任何抗拒能力。应该强调指出，真正的理想流体在客观实际中不存在，它只是实际流体在某种条件下的一种近似模型。Inviscid模型中无粘性流体的动量守恒方程即Euler方程(3）。

RNG k － ε 湍流模型最早由 Yakhot及 Orzag［8］提出，是一种对瞬时的N-S方程用重整化群的方法推导出来的湍流粘性模型，目前在船舶CFD领域中应用较为广泛。其k方程与ε方程为:

与标准k－ε湍流模型相比，RNG k－ε湍流模型可以更好地处理高应变率以及流线弯曲程度较大的流动。RNG k－ε模型仍是针对充分发展的湍流有效，对于近壁区内的流动及雷诺数较低的流动，必须采用壁面函数法进行处理，即采用一组半经验公式将近壁区域的物理量与湍流充分发展区的物理量联系起来。

1.3 造波和消波

数值波浪水池最关键的技术是造波和消波。目前，常用的造波方法主要有源造波、设置速度入口造波和仿物理造波。仿物理造波较前两者要复杂，但是造波性能好、质量高且稳定性好。仿物理造波又被称为动边界造波，主要有摇板式造波和推板式造波2种。摇板造波动边界区域需要划分非结构网格，动网格方法过于复杂。推板造波只需划分结构网格，配合使用Layering铺层动网格方法即可实现，因此推板造波在有限水深数值水池中使用更为普遍。

根据推板造波理论，推板造波的传递函数［9］为:对于平衡位置在原点的推板式造波机，其往复式推板在水池中作正弦运动，运动轨迹及速度方程为:

式中:A为波幅;k为波数;ω为角频率;h为水深。

为防止波浪到达水池后方边界后产生反射波，对水池工作区域的正常使用产生干扰，还需要在水池后方进行消波处理。目前，数值波浪水池中常用的消波方法有设置辐射边界条件法、主动消波法和设置阻尼区消波法。本文采用阻尼消波法，具体实现方式为在水池后端设置消波区域，通过在消波区的动量方程中加载自编的DEFINE_SOURCE(source，c，t，dS，eqn）宏源项来实现尾端消波。消波区内动量方程为:

式中:C(x）为单调递增的消波系数，在消波区起点其值为0，以消除波浪传播到阻尼区时因阻尼存在产生提前反射的现象，保证波浪在平稳前进的过程中实现消波;υ为运动粘性系数。为增强固定长度的消波区的适用性，本文提出一种和波长、密度都有关的计算公式

式中:λ为波长;x0，xL分别为消波区左、右边界的x坐标值;ρ为流体密度。

2 数值计算模型

2.1 数值水池

三维数值水池总长40 m，宽4 m，高3 m，水深2 m，左端推板造波区长0.5 m，右端消波区长10 m，中间部分为工作区，水池最左侧放置推板造波机，笛卡尔正交坐标系的原点设置在水池中纵剖面上推板初始位置与水面的交界处，整个造波水池如图1所示。

图1 三维数值水池示意图Fig.1 Schematic diagram of the 3D numericalwave tank

网格全部采用结构化正六面体网格，总网格数为1 294 440。网格质量的好坏对数值结果有直接影响，为了保证精度，需参照目标波浪的波高和波长进行网格划分，至少保证1个波高范围内有5个网格，1个波长范围内有20个网格。在推板和自由液面附近对网格进行加密，造波区网格长度为0.01 m，自由液面附近网格高度为0.02 m，远离自由液面处网格以一比例系数向上下两侧逐渐稀疏，工作区网格长度为0.08 m，尾端消波区网格则沿x轴正方向逐渐稀疏，以利用数值耗散加强消波效果，y方向网格宽度统一为0.1 m。oxz平面的网格划分如图2所示。

图2 oxz平面网格划分示意图Fig.2 Schematic diagram of the grid on oxz plane

数值波浪水池初始条件设置:将水池以oxy平面为分界面分为两流体域，上方1 m水体积分数定义为0，即全部填充空气，密度为1.225 kg/m3;下方2 m水体积分数定义为1，即全部填充水，密度为998.2 kg/m3。流体初始速度均为0，参考压力值为101 325 Pa，通过 Define Custom Field Function函数定义静水压强。

2.2 数值计算

本文算例以Fluent瞬态求解器进行数值计算。推板的往复式运动通过动网格铺层技术Layering实现，其中分离因子Split Factor和合并因子Collapse Factor均设置为0.4。控制方程采用有限体积法进行离散，压力速度耦合方法为PISO(Pressure Implicit with Splitting of Operator），压力插值采用Body Force Weight体积力方式。为避免动网格运动剧烈出现负体积，时间步长不易过大，文中算例均为0.005 s，最大迭代次数为20。为了研究粘性对造波的影响，文中每个算例均采用Inviscid理想模型和RNG k－ε湍流粘性模型各算1次，其余设置保持一致。

通过Iso-Surface在水池y=0 m和y=2.0 m两纵剖面上平行设置波高监测仪，监测仪沿x方向的坐标位置分别为x=1 m，5 m，10 m，15 m，20 m，25 m，30 m，33 m，36 m，39 m。

本文要模拟的余弦波工况如表1所示。

表1 余弦波计算工况Tab.1 Cases of the cosine wave

3 数值模拟结果与分析

3.1 造波和消波结果的验证

图3～图5均为由Inviscid模型得到的Case1的数值计算结果。

图3 Case1，λ=4.0 m，x=10 m，20 m处的波高监测值Fig.3 Case1，λ =4.0 m，wave elevations at different positions:x=10 m，20 m

图3为分布在水池中4个位置 (x=10 m，y=2m;x=10m，y=0m;x=20 m，y=2 m;x=20 m，y=0 m）的波高仪监测到的数值结果与理论值的对比图。由图3可知:本文提出的平滑函数效果良好，使流体由静态到动态得到了平滑过渡，有效降低了静水初始效应;利用数值波浪水池推板造波得到的波形与理论波形吻合良好，证实了所构建的三维数值波浪水池具有良好的造波功能。

图4 Case1，λ=4.0 m，t=30 s时刻，数值波浪水池Fig.4 Case1，λ =4.0 m，t=30 s，numerical wave tank

图4真实展现了t=30 s时的数值波浪水池情况。由图可见:水池工作区内波形规则平稳，完全可以放置结构物进行波浪与结构物相互作用方面的研究;消波区内则因阻尼耗散作用使得波浪逐渐衰减，避免了波浪反射;整个水池内沿宽度方向波形保持一致，避免了固壁影响。由图3也可以看出，y=2 m与y=0 m两平面上的数值结果完全一致，这主要是因为对以y轴为法向量的前后两侧面进行了symmetry边界条件设置，保证了两侧面上无动量交换，从而避免了以往固壁边界对流体的不良影响。

图5 Case1，λ =4.0 m，A=0.15 m，t=30 s时刻，数值水池中纵剖面x方向速度、z方向速度云图Fig.5 Case1，λ =4.0 m，A=0.15 m，t=30 s，x-velocity and z-velocity contour at the NWT y-mid plan

图5所示的速度分量云图清晰显示了数值水池内平稳有序的流体运动，由图亦可以看出，本文提出的消波方法效果良好，进入消波区后，x方向速度和z方向速度因受阻尼耗散和数值耗散的双重作用而逐渐减小，进而趋于0，成功避免了波浪到达后端固壁后产生反射波。

经上述验证，本文所建立的推板造波数值水池不仅具有良好的造波功能，而且具有良好的消波功能，完全可以进行有关波浪的数值模拟研究。

3.2 粘性影响分析

图6 Case5，λ =8.0 m，A=0.15 m，t=50 s时刻，数值波浪水池中纵剖面波形图Fig.6 Case5，λ =8.0 m，A=0.15 m，t=50 s，shape of waves at the NWT y-mid plan

以Case5为例，从图6可以看出，在数值波浪水池的工作区 (0 m≤x≤30 m），由Invicid理想流体模型得到的数值模拟结果和余弦波的线性理论值吻合良好。但由RNG k－ε湍流粘性模型获得的波高值在水池长度方向上有一定的衰减，这主要是由于RNG k－ε模型考虑了流体粘性作用而导致的 (x＞30 m时的衰减主要为消波所致）。

图7 Case5，λ =8.0 m，A=0.15 m，t=50 s，不同横剖面(x=5 m，x=15 m，x=25 m）的速度曲线对比图Fig.7 Case5，λ =8.0 m，A=0.15 m，t=50 s，comparison of velocity at different transverse sections(x=5 m，x=15 m，x=25 m）

进一步分析RNG k－ε模型计算得到的Case5数值结果，将水池x=5 m，x=15 m，x=25 m三个横剖面上的x方向速度分量和z方向速度分量与由公式(13）和式(14）得到的理论值进行绘图对比 (图例中后缀S.表示模拟值，L.表示理论值），发现自由液面附近x方向和y方向上的速度分量绝对值均小于理论速度绝对值，从而说明粘性力对流体质点的运动存在阻滞作用。

图8 波高沿程衰减曲线Fig.8 Attenuation curve of wave height

图8为波高沿程衰减曲线，当横坐标x＞30 m时，波高的急剧衰减是由所采取的强制阻尼消波造成的，同时进一步验证了文中所构建的数值波浪水池具有良好的消波效果。除此以外，波高的沿程衰减趋势基本符合指数形式，衰减程度与波陡正相关。由Case4，Case5，Case6对比可知，同一波长下，波陡越大，衰减幅度越大，波能耗散先快后慢的现象越明显;由Case1，Case4或Case2，Case3，Case6对比可知，即使波长不同，只要波陡相等，波浪稳定后的数值水池两端的波高衰减幅度相等，只是波长大的波能耗散速度相对更为均匀。

3.3 波高衰减拟合研究

为了更准确地预测波高衰减，从而在数值波浪水池中特定区域获得预定波高，进一步开展波浪与海洋结构物相互作用方面的研究，本文利用Matlab软件进行回归分析编程，对所构建的数值水池的衰减规律作了数值拟合研究。

波高沿程衰减形式符合指数规律:

应用上述拟合公式绘制的波高衰减曲线如图9虚线所示。通过与波高仪监测数据(标记点所示）进行对比可以看到，利用本文提出的公式(17）描述数值水池中波浪因粘性作用导致的沿程衰减规律是合适，拟合效果令人满意，较之常用的线性拟合(图8所示）更为平滑准确。因此，其他从事粘性数值波浪水池研究的学者亦可以采用式(17）的形式对波高进行拟合预测，只要根据实际情况将拟合系数a，b，c稍加修改即可。

图9 波高衰减拟合曲线Fig.9 Fitted attenuation curve ofwave height

4 结语

本文以商业软件Ansys14.0中的CFD模块Fluent为平台，通过程序接口加载UDF程序进行二次开发，采用VOF法结合几何重构对自由液面进行追踪，利用动网格铺层技术Layering实现造波板的运动，同时在动量方程中添加源项实现阻尼消波，成功建立了具有良好造波、消波功能的三维推板式数值波浪水池。

利用此数值水池，结合理想流体Inviscid模型和湍流RNG k－ε模型，通过模拟具有不同波长、波高的余弦波研究了粘性在数值波浪水池中的影响。结果分析表明，粘性力对流体质点的运动存在阻滞作用，从而使波高沿水池长度方向发生了明显的衰减，衰减幅度与波陡呈正比，衰减形式符合指数规律。

利用Matlab软件进行回归分析编程，对所构建的数值水池中波浪的衰减规律作了数值拟合研究，提出了一种拟合效果十分令人满意的衰减公式形式，从而为更准确地预测波高衰减，进一步开展波浪与海洋结构物相互作用方面的研究提供参考。

［1］谷汉斌，李炎保，李绍武.基于完全非线性Boussinesq方程的源函数数值造波［J］.海洋技术，2004，23(2）:80-85.

GU Han-bin，LIYan-bao，LIShao-wu.The source function wave-generating and wave data analysis method based on fully nonlinear Boussinesq equations［J］. Ocean Technology，2004，23(2）:80-85.

［2］WANG Da-guo，ZOU Zhi-li，LIU Xia.Fully nonlinear 3-D numerical wave tank simulation［J］.Journal of Ship Mechanics，2010，14(6）:577-586.

［3］宁德志，陈丽芬，田宏光.波流混合作用的完全非线性数值水槽模型［J］.哈尔滨工程大学学报，2010，31(11）:1450-1455.

NING De-zhi，CHEN Li-fen，TIAN Hong-guang.A fully nonlinear numerical flume model for wave-current interactions［J］.Journal of Harbin Engineering University，2010，31(11）:1450-1455.

［4］FINNEGANW，GOGGINS J.Numerical simulation of linear water waves and wave-structure interaction［J］.Ocean Engineering，2012，43:23-31.

［5］DUCROZET G，BONNEFOY F，TOUZÉ D L，et al.A modified high-order spectral method for wavemaker modeling in a numericalwave tank［J］.European Journal of Mechanics B/Fluids，2012，34:19-34.

［6］周光埛，严宗毅，许世雄，等.流体力学［M］.北京:高等教育出版社，2011:114-123.

ZHOU Guang-jiong，YAN Zong-yi，XU Shi-xiong，et al.Fluid mechanics［M］.Beijing:Higher Education Press，2011:114-123.

［7］杨锦凌，孙大鹏.基于FLUENT二次开发的数值波浪水槽［J］.中国水运，2012，12(5）:59-61.

YANG Jin-ling，SUN Da-peng.Numerical wave tank based on FLUENT secondary development［J］.China Water Transport，2012，12(5）:59-61.

［8］YAKHOT V，ORZAG S A.Renormalization group analysis of turbulence:basic theory［J］.JScient Comput，1986，1:3-11.

［9］俞聿修，柳淑学.随机波浪及其工程应用［M］.大连:大连理工大学出版社，2011:149-153，241-248.

YU Yu-xiu，LIU Shu-xue.Random wave and its applications to engineering［M］.Dalian:Dalian University of Technology Press，2011:149-153，241-248.

Simulation of 3D numerical wave tank and viscosity research

ZHAO Yan，ZHU Ren-qing，LIU Zhen
(Jiangsu University of Science and Technology，Zhenjiang 212003，China）

In this paper，based on CFD software Fluent，one module of the Ansys14.0，a threedimensional numericalwave tank with piston-like wavemaker is established by adding UDF programs to secondary development.The volume of fluid(VOF）method and Geo-Reconstructmethod is used to track the free surface.The motion of wave maker is achieved by using of dynamic mesh method Layering.Meanwhile，source terms are added in themomentum equation in order tomake wave damp.The relevant simulations verify that this numericalwave tank has effective wave-making and wave-absorbing function.On this basis，using the Inviscid model and RNG k－εturbulencemodel，series of cosine waveswith different length or height are simulated to research the influence of viscosity in NWT.The simulation result shows that viscosity can reduce velocity of fluid particle，so that results in wave height attenuation along the length of NWT obviously.The attenuation amplitude is in direct proportion to wave steepness and attenuation form corresponds to exponential law.

viscosity;attenuation law;numerical wave tank;VOF

TV139.2

A

1672-7649(2014）05-0042-07

10.3404/j.issn.1672-7649.2014.05.009

2013-10-21;

2013-11-04

国家自然科学基金资助项目(51179077，51209107），江苏高校优势学科建设工程资助项目

赵艳(1987-），女，硕士研究生，从事波浪与海洋结构物相互作用研究。