二阶入射边界条件在完全非线性数值波浪水槽中的应用

收稿日期：2009-10- 基金项目：国家自然科学基金（50709005和50921001），海洋工程国家重点实验室开放基金（0901）和中央高校基本科研业务费专项资金（DUT10JN03）资助.作者简介：周斌珍（1984-），女，福建省宁德市人，博士生，主要从事非线性波浪与海洋结构物作用的研究。Email: zhoubinzhen4827. 完全非线性深水波的数值模拟周斌珍1, 宁德志1 2, 滕斌1, 宋伟华1(1. 大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024；2. 上海交通大学 海洋工程国家重点实验室，上海， 200030)摘 要：基于势流理论，并结合深水波质点运动从水面向下呈e指数衰减的特性，建立了完全非线性数值变深水槽模型，通过实时模拟活塞式造波机运动来产生波浪。采用时域高阶边界元法进行模拟，利用混合欧拉-拉格朗日方法和四阶Runge-Kutta方法追踪流体瞬时水面，应用镜像格林函数消除了水槽两个侧面的积分，在水槽末端布置人工阻尼层消除反射波浪。利用GMRES加速算法，提高了计算效率。利用所建模型对深水波进行了模拟研究，得到了稳定的波形，在造波板运动幅值较小时，与线性解析结果吻合良好，当造波板运动幅值较大时，体现出波浪的非线性特性。同时利用所建模型研究了造波板所在位置（上部台阶）水深对生成波浪高度的影响，由此可以选择合适的造波板所在位置水深及造波板运动幅值来得到所需要的深水波浪。关键词：造波板运动；数值变深水槽；完全非线性；高阶边界元；深水波中图文分类号： U661 文献标识码：A 文章编号：0253-4193（2009）00-0000-001 引 言近年来，由于中等和浅水深的开发技术已经成熟，以及当今世界大部分地区浅海油气资源已日渐减少，因而深海资源开发已经成为当今世界能源开发的热点和必然趋势。许多专家和学者对深水问题开展了各种理论研究工作，Newman1通过水平台阶理论分析了波浪从浅水到无限水深和无限水深到浅水的传播变化；Fenton2 通过解析研究得到了五阶深水波浪速度势和波面的解析解；Ferront3通过频域和时域方法研究了深水中非线性波浪辐射问题。然而，由于理论分析适用范围的局限性，开展深水问题的数值模型试验是很重要的。国内外专家开展了许多相关研究工作，且都希望建立一个与物理波浪水槽功能相当的数值模型。譬如Lee等4通过截取有限长度计算域来研究深水中运动物体的荷载；Fonseca 和 Soares5利用满足自由水面条件的格林函数对无限水深中船舶做大振幅运动问题进行了数值模拟研究；宁德志等6, 7利用源造波技术对无限水深的完全非线性波浪进行了数值模拟研究。然而，大多数数值模型中都是通过在入射边界处给定特定的理论波浪特性（如速度或波高）来实现造波，尽管这样可以令入射边界条件固定而避免网格的重新划分，使问题求解相对简单，但在入射边界给定一个与实际波浪非线性相匹配的理论速度或波高通常是很难的8，特别是对于造波板做大振幅运动的情况。而通过模拟造波板实时运动产生波浪则可以解决入射边界不匹配问题，最终真实模拟物理波浪水槽中的自然现象。如Bai和Eatock Taylor9-11 采用高阶边界元方法，在数值波浪水槽中模拟造波板实时运动产生波浪，分别对完全非线性规则波、聚焦波与垂直圆柱相互作用问题及完全非线性波浪对固定和漂浮结构物作用问题进行了模拟；周斌珍等12和Ning等13采用类似的方法建立了完全非线性数值波浪水槽模型，并对规则波、不规则波进行了模拟。以往的波浪水槽大多以平底水槽为主，在水槽一端设置造波机，通过给定造波机的运动形式产生所需要的波浪。对于深水波而言，质点运动从水面向下呈e指数衰减，因此只需要把造波机设置在水面附近一定深度，只在表面造波，由此建立台阶式变深水槽。Zhou等14等利用特征函数展开法对台阶式变深水槽中推板式造波机造波问题进行了解析研究。 本文利用时域高阶边界元方法建立了一个完全非线性数值变深水槽，通过实时模拟造波板在上部台阶运动来产生深水波浪，且自由水面满足完全非线性边界条件。在每一时间步内，利用四阶Runge-Kutta法和混合欧拉-拉格朗日法对瞬时水面和入射边界面进行更新；为了避免由于网格运动变形而引起数值不稳定问题，在一定时间步内对网格进行重新划分，通过采用二次形状函数计算新网格的物理量和几何量；将镜像格林函数运用到所建立的边界积分方程中，进而消除水槽两侧面的积分；通过在入射边界面和台阶立面上网格沿水深方向逐渐变粗的方法，减少了计算所需要的网格；利用GMRES加速算法，提高了计算效率。利用所建模型对深水波进行了模拟，并与线性解析结果进行了对比验证，同时利用所建模型研究了造波板所在位置（上部台阶）水深对生成波浪高度的影响。2 控制方程及初边值条件 考虑如图1所示的台阶式变深水槽，在静水面上建立笛卡尔坐标系，z轴向上为正，并于台阶立面重合，x轴向右为正。波浪沿x轴正方向传播。在水槽的左端装有造波机，水槽右端为吸收边界。这里，用t代表时间，h表示波面高度，自由水面、固体边界和水槽右端散射边界分别用Sf 、SN（包括上部台阶水底Sd1、下部台阶水底Sd2、水槽两个侧面Sc、入射边界面S1、台阶立面S2）和Sr表 示 。 假定流体均匀、不可压缩、无黏和无旋，且忽略表面张力，则流速可以定义为速度势的梯度f，速度势f(x, y, z, t)满足拉普拉斯方程： (1)式中，是三维拉普拉斯算子。图1 变深水槽示意图 在瞬时自由水面Sf上，满足完全非线性运动和动力边界条件。 本文采用拉格朗日法更新自由水面，故自由水面边界条件可以写成以下拉格朗日形式15： (2)式中，g是重力加速度，X=(x, y, z)为瞬时自由表面上任意流体质点的位置矢量， 为物质导数算子。自由表面采用拉格朗日法就是自由表面位置由自由表面流体质点运动速度对时间的积分确定。 在固体物面SN上，速度势满足边界条件： (3)式中，n为物面单元单位法向导数，指出物体为正，Un是预先知道的物面速度，在水槽的两个侧壁Sc和水槽底Sd1、Sd2上，Un =0；在水槽左侧为造波板S1上对应点的运动速度。假定自由水面初始时是静止的，即 (4) 在水槽的右端Sr上，波浪满足向外传播的辐射条件。本文采用在下游水面区加一人工阻尼层来吸收向右传播的波浪，通过在自由面运动学和动力学边界条件(2)中加入阻尼项来实现，即 (5)式中，阻尼系数 (6)X0=(x0, y0, 0)是指水质点静止时的位置，x0为阻尼层起始坐标，a为阻尼系数，bl为阻尼层长度，l为波长，b为岸滩宽度系数，本文取a=1，b=1。3 数值模拟方法3.1 边界积分方程 在整个流域内应用格林第二定理，则上述边值问题可转化为如下的边界积分方程： (7)式中，=(x, y, z)为源点，q=(x0, y0, z0)为场点，a(p)为固角系数。边界S包括自由水面Sf和固体表面SN，G是简单格林函数。 为减少计算量，本文将镜像格林函数应用到水槽侧壁，消除了两个侧面的积分，同时，自由水面、入射边界面、台阶立面、上部台阶底部和下部台阶底部可以采用不同密度的网格，不受侧壁的约束，上部台阶和下部台阶可采用较疏的网格，减小了整体的网格数。格林函数通过两水槽侧壁无限次反射的像叠加得到。为确保格林函数收敛， 每一项均减去因子。此格林函数可以写成如下形式： (8) 式中，X=x-x0，Y=y-y0，Z=z-z0，B 是水槽宽度。式(8)中的叠加求和各项分别代表点源在两水槽侧壁中所形成的镜像。对于式中这个趋于无穷的计算式子，按照Newman 方法将其转化为下式快速精确求解16。 (9)式中: K0 是修正汉科尔(Hankel) 函数，g= 0. 577778 是常数，当水槽宽度大于1 时，m 最大项数截取到六项时，式(9) 就可以达到7 位小数的精度。3.2 积分方程的离散 本文用高阶边界元离散计算域成一些曲面单元，对每个单元通过数学变换，将其转换成参数坐标(x, )下的等参单元，采用二次形状函数插值方法保证单元内物理量分布的连续性。 单元内任一点的几何坐标和速度势等物理量可以写成如下形式： ， ， (10)式中，和 分别是节点k的坐标，速度势，速度势法向导数和形状函数； K 是单元节点个数。 将式(8) 、(10)代入式(7) 中，积分方程可以离散为如下形式： (11)式中，J(x, )是联系广义坐标和局部坐标的雅可比矩阵，Ne1 和 Ne2分别是自由水面和固体边界面SN（入射边界面S1、台阶立面S2、上部台阶底部Sd1和下部台阶底部Sd2，计算深水波时Sd2的影响可忽略，不用划分网格）上的单元个数。当源点p在自由水面上时，速度势为已知量，这时方程(11)左端第一项将被移到方程右端；当源点在固体边界面上时，速度势为未知量，这时方程(11)左端第一项将保留在方程左端，这样方程左端均为未知量，右端均为已知量。由此通过配点法，将源点p分别取在各个节点上，可建立如下线性方程组： (12) 其中， (13) (14) (15) (16) 由于积分边界是不断地随着时间变化的，在每一计算时刻都要重新建立系数矩阵，并且在每一计算时刻都要对方程进行求解。计算中认为当前时刻物面SN上的速度势法向导数和自由水面Sf上的速度势是已知的，根据积分方程(11)计算当前时刻物面SN上的速度势和自由水面Sf上的法向速度，然后应用四阶Runge-Kutta法，根据自由水面条件(2)计算下一时刻的水面高度和自由水面Sf上的速度势，根据物体运动方程和物面条件，计算下一时刻物体的运动位置、物体的运动速度和物面上的法向速度，再对物面和自由水面重新划分网格，重新应用积分方程计算下一时刻物面上的速度势和自由水面上的法向速度。这样计算周而复始，直到计算结束。3.3 网格划分 网格划分包括直立面（造波板、台阶立面）和水平面（自由水面、台阶底面）。由于自由水面和造波板位置的变化，计算域的大小每一时间步都在改变，网格相应地跟着改变，为了防止由于造波板或结构物附近网格大小不均匀引起数值不稳定问题，需要定期进行网格重新划分来调整这种网格不均匀问题，当造波板运动幅度较大时每一时间步都要进行网格重分。网格的重新划分过程详见参考文献12。 对于深水波而言，质点运动从水面向下呈e指数衰减，因此网格沿水深方向可逐渐变粗，这样可以减少网格的数量，如下式所示： (17)Z1,j指水面坐标或台阶立面顶高程的坐标，Mz指造波板或台阶立面上z方向划分的单元个数，g表示网格逼近顶面的疏密程度，单元中间点的坐标等于相邻两节点的平均值。以下举例说明g大小对网格以及计算效率的影响。应用所建立模型模拟深水波，考虑造波板所在位置水深d1=5.0m，水槽工作区水深d2=25.0m，造波板与台阶立面距离B1=3.0m，造波板运动周期T=4s，造波板运动幅值为a=0.7m。定义水槽长度和宽度分别为L=4.0l（波长l24.96m）和B=0.125l，为了吸收出流波浪，在水槽末端布置1.5l长的阻尼层。取空间步长Dx=l/16，Dy=B/2，时间步长Dt=T/60，其余网格步长如表1所示（其中g1表示造波板上网格逼近自由水面的程度，g2表示台阶立面网格逼近台阶顶面高程的程度）。图2是入射边界面取6个网格，台阶立面取10个网格， g取不同的值时，水槽网格布置图。从图中可以看出，当g趋近于0时，网格接近于均分，随着g的增大，网格从上往下逐渐增大。取表1中较密的网格一作为参照标准，当其它情况的结果与网格一结果基本接近时，认为网格收敛。从表1可以看出，当g2很小，网格接近均分时，台阶立面网格需取6个，结果才能收敛，当g2=0.2时，台阶立面网格只需3个就能使结果收敛，因此可看出选择合适的g值可以减小网格数，对计算效率有着重要的作用，尤其对于强非线性波浪的工况。 (a) g1=g2=0.001 (b) g1=g2=0.1 (c) g1=g2=0.3 (d) g1=g2=0.5图2 变深水槽网格布置图表1 网格分布与收敛情况网格入射边界面单元数台阶立面单元数g1g2收敛性网格一6100.10.1收敛网格二230.0010.001不收敛网格三230.0010.2收敛网格四260.0010.001收敛3.4 数值计算方法传统边界元法求解线性方程组时，一般采用高斯消去法或是由高斯消去法派生的相关方法(如LU分解法)，解方程所需的计算次数与方程秩的三次方成比例。因此当方程组规模变大时，利用高斯消去法为基础的直接法所需的求解时间迅速增大，对此类问题的求解，迭代方法成为了有力的工具。而GMRES方法17是在Arnodi算法的基础上提出的解决具有稠密系数矩阵的线性方程组的有力的数值解法，其计算次数明显减少，因此使用GMRES算法求解大型边界元法系数方程组是十分有效的。本文以高斯消去法和GMRES法对比来说明GMRES算法的快速性。 图3为对整个区域划分不同单元数目时，GMRES算法和高斯消去法（Gauss）求解每一时间步时所用时间的比较图。N表示未知量个数。从图中可以看出，当未知量个数较少时，采用高斯消去法和GMRES法用时相近，而随着未知量的增加，GMRES方法则更为省时高效。图3 计算量对比图4 数值计算及结果 作为算例，考虑推板式造波机做正弦运动时产生波浪的情况（通过调整初相位j也可得到造波机余弦运动情况），造波板运动方程写为： (16)进而可以得到造波板瞬时位置的运动速度： (17)式中，a是造波板的运动幅值，w是造波板运动圆频率。应用所建立模型模拟深水波，考虑造波板所在位置水深d1=5.0m，水槽工作区水深d2=25.0m，造波板与台阶立面距离B1=3.0m，造波板运动周期T=4s，相位角j=0。定义水槽长度和宽度分别为L=4.0l（波长l24.96m）和B=0.125l，为了吸收出流波浪，在水槽末端布置1.5l长的阻尼层。通过对时间步长和空间步长开展收敛性数值试验，确定空间步长Dx=l/16，Dy=B/2，时间步长Dt=T/60，上部台阶底部取3个网格，g1=0.001，台阶立面取4个网格，g2=0.2。图4和图5分别是造波板运动幅值a=0.05m和a=0.7m时，由本文方法计算的不同位置处波面时间历程与线性解析结果14在当前边界条件下计算结果的对比。从图4可以看出，在造波板运动幅值很小，即波浪非线性很小的情况下，两种方法得到的结果吻合良好，验证了本文方法的正确性。从图5可以看出，当造波板运动幅值较大时，由本文方法计算得到的波形的非线性特性更加明显，波峰更加陡峭，波谷更加趋于平缓。为了进一步分析这种非线性特性，运用傅里叶变换的方法分别对图5中波面结果进行处理，得到如图6所示的波浪频谱变化关系，从图中可以看出本文结果和线性解析结果在主频（f0=1/T）上吻合很好，进一步验证了本文方法的正确性，而本文结果除了主频的贡献外，还存在零频、二倍频甚至三倍频的贡献。图7和图8分别是两种波况下t=9T和t=10T沿水槽中线处波面分布图，从图中可以看出，两个时刻的波面完全重合，说明数值结果已经达到稳定。图9为t=10T，a=0.7m时三维波面变化图，从图中可以看出，波面在Y方向上保持稳定。 (a) x/d1=-0.3 (b) x/d1=0.0 (c) x/d1=5.0图4 不同位置波面升高时间历程（a=0.05m） (a) x/d1=-0.3 (b) x/d1=0.0 (c) x/d1=5.0图5 不同位置波面升高时间历程（a=0.7m） (a) x/d1= -0.3 (b) x/d1=0.0(c) x/d1=5.0图6 不同位置波面频谱变化关系图（a=0.7m） 图7 水槽中线处波面变化图(a=0.05m) 图8 水槽中线处波面变化图(a=0.7m)图9 t=10T时三维波面变化图（a=0.7m）5 特性分析 为了分析台阶式变深水槽中工作区生成波浪高度的影响因素，本文研究了造波板所在位置（上部台阶）水深对水槽中工作区波高的变化，以下算例采用与上例相同的时间步长和空间步长。图10为水槽工作区水深d2=25.0m，造波板与台阶立面距离B=3.0m，造波板运动周期T=4s，造波板运动幅值a=0.7m，造波板所在位置水深分别为d1 =4.0m，5.0 m，6.0 m时，x=l处波面升高时间历程图。图11是对图10中的波面进行傅里叶分析得到的频谱变化关系图。从图中可以看出，随着造波板所在位置水深的增加，水槽工作区某一固定水质点的波高增大，且波浪的非线性更强。 图10 x=l处波面升高时间历程 图11 x=l处波面频谱变化关系从图4中的线性解析解可以看出，通过增加造波板的运动幅值，可以加大水槽工作区的波高，从图10和11的分析分析中可以看出，通过增加造波板所在位置水深，也可以加大水槽工作区的波高。当其他条件一定，只增加造波板运动幅值时，试验中会因为造波板推程的限制无法得到所需要的波高，数值中会出现数值不稳定问题而导致计算无法进行，因此通过调节造波板所在位置水深得到所需要的波高是非常重要的。6 结论为了避免入射边界条件和完全非线性自由水面条件的不匹配，并结合深水波传播特性，建立了一个与真实物理波浪水槽功能相当的台阶式变深水槽数值模型。本文通过实时模拟造波板的运动来产生波浪，并通过混合欧拉-拉格朗日方法和四阶RK技术更新自由水面和造波板的瞬时位置。数值模型通过时域高阶边界元方法进行求解，其中镜像格林函数被应用到水槽两个侧壁，使得水槽侧面的积分被消除，同时其余固体边界可以采用不同密度的网格，不受侧壁的约束，减小了整体的网格数，节省了大量的计算量。利用GMRES加速算法求解线性方程组，提高了计算效率。利用所建模型对深水规则波进行了数值模拟并对波浪的非线性进行了分析，在造波板运动幅值较小时，与线性解析结果吻合很好，当造波板运动幅值较大时，波浪的非线性特性显现出来。同时利用所建模型研究了造波板所在位置（上部台阶）水深对生成波浪高度的影响，研究发现随着造波板所在位置水深的增加，生成波浪高度增大，且波浪的非线性更强，由此可以选择合适的造波板所在位置水深及造波板运动幅值来得到所需要的波浪高度。通过以上研究，验证了本文模型可以准确模拟非线性规则波。本文模型为研究深水波浪与海洋结构物相互作用问题提供了重要的参考，并可进一步应用到深水波浪与海洋结构物相互作用的模拟。参考文献:1 NEWMAN J N. Propagation of water waves over an infinite stepJ. Fluid Mech, 1965, 23(2):399-4152 FENTON J D. A fifth-order Stokes theory for steady wavesJ. Journal of Waterway, Port, Coastal and Ocean Engineering, 1985,111(2) :216-234.3 FERRANT P. A coupled time and frequency approach for nonlinear wave radiationC. Proceedings of the 18th symposium on Naval and Hydrodynamics, Michigan, USA, 1990.4 LEE C C, LIU Y H, KIM C H. Simulation of nonlinear waves and forces due to transient and steady motion of submerged sphereJ. International Journal of Offshore and Polar Engineering, 1994, 4 (3): 174182.5 FONSECA N, SOARES C. Time-domain analysis of large-amplitude vertical ship motions and wave loadJ. Journal of Ship Research, 1998, 42(2):139-153.6 宁德志，滕斌，刘珍，等. 无限水深聚焦波的完全非线性数值模拟J. 海洋工程，2008，26(4):27-34.7 NING D Z, TENG B, EATOCK TAYLOR R, et al, Numerical simulation of nonlinear regular and focused waves in an infinite water-depthJ. Ocean Engineering, 2008, 35:887-899.8 RYU S, KIM M H, LYNETT P. Fully nonlinear wave-current interactions and kinematics by a BEM-based numerical wave tank J. Computat. Mech. 2003, 32:336-346.9 BAI W, EATOCK TAYLOR R. Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders J. Applied Ocean Research, 2006, 28:247265.10 BAI W, EATOCK TAYLOR R. Numerical simulation of fully nonlinear regular and focused wave diffraction around a vertical cylinder using domain decomposition J.Applied Ocean Research, 2007, 29:55-71.11 BAI W, EATOCK TAYLOR R. Fully nonlinear simulation of wave interaction with fixed and floating flares structures J.Ocean engineering, 2009, 36:223-236.12 周斌珍，宁德志，滕斌. 造波板运动造波的实时模拟J. 水动力学研究与进展，A辑，2009，24(4)：406-416.13 NING D Z, ZhOU B Z, TENG B, et al. Numerical simulation of nonlinear regular and focused waves generated by a piston wave makerC. The 5th international conference on Asian and Pacific Coasts. Singapore, 2009.14 ZHOU B Z, NING D Z, TENG B, et al. Analytical study on wave making in a deep wave flume in step-typeJ, Wave Motion , 2010, 47(1):1-11.15 LIN W M, NEWMAN J N, YUE D K. Nonlinear forced motion of floating bodiesC. Proc.15th Symp. On Naval Hydrodynamics, Hamburg, Germany, 1984, 33-47.16 NEWMAN J N. Approximation of free2surface Green functions A . In : P. A. Martin & G. R. Whickham (eds) , Wave AsympoticsM . Cambridge University Press , 1992. 107-142.17 SAAD Y and SCHULTZ M H. A generalized GMRES algorithm for solving nonsymmetric linear systemsJ. SIAM J. Sci. Statist. Comput, 1986, 7:856-869Numerical simulation of fully nonlinear deep-water wavesZHOU Bin-zhen1, NING De-zhi1 2, TENG Bin1, SONG Wei-hua1(1. State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China; 2. Stae Key Laboratory of Ocean Engineering, Shanghai Jiaotong University, Shanghai, 200030, China.)Abstract: Based on potential theory and the particle velocity decreasing exponentially with the increasing of water depth for the deep-water waves, a fully nonlinear numerical variable deep wave flume model was developed. Waves generated by a piston wave maker were real-time simulated. The model is developed by using a time domain higher-order boundary element method (HOBEM). A mixed Eulerian-Lagrangian technique and a 4th-order Runge-Kutta scheme are utilized to track the free surface. Image Green function is used in the whole fluid domain so that the integration on the lateral surfaces are excluded. An artificial damping layer is distributed at the end of the flume to eliminate wave reflection. The GMRES accelerated algorithm is utilized to improve calculation efficiency. Numerical experiments are carried out to model the deep water waves. Steady wave profiles are obtained and good agreements between numerical solutions and analytical solutions are obtained for the small motion amplitude of the wave maker. The wave nonlinear features are shown for the large motion amplitude. Numerical experiments are also carried out to study the effect on the generated wave height due to the water depth at the up step. According to the proposed numerical simulation, the suitable water depth at the up step and the motion amplitude of the wave maker can be derived to get the required wave in deep water. Key words： Wave maker motion; numerical wave tank; fully nonlinear; HOBEM; deep water waves14