第１４卷第５期 ２０１０年５月 船舶力学 Ｊｏｕｒｎａｌ ｏｆ Ｓｈｉｐ Ｍｅｃｈａｎｉｃｓ Ｖｏ１．１４ Ｎｏ．５
Ｍａｙ ２０１０
文章编号：１００７—７２９４（２０１０）０５—０４８１—０６
不规则波浪的数值模拟
梁修锋，杨建民，李 欣，李 俊
（上海交通大学海洋工程国家重点实验室，上海２０００３０）
摘要：优良的波浪环境条件是研究非线性波浪和浮式结构物相互作用问题的基础。文章在数值波浪水池中对不
规则波浪进行了数值模拟，造波通过摇板的运动实现，消波通过在水池后段设置消波区实现，自由液面由ＶＯＦ
方法捕捉。对计算所得的波浪时历进行谱分析并将所得谱与理论谱进行了比较，吻合良好。
关键词：数值波浪水池；动网格；ＶＯＦ方法；ＪＯＮＳＷＡＰ谱
中图分类号：Ｕ６６１．７ 文献标识码：Ａ
Ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｉｒｒｅｇｕｌａｒ ｗａｖｅ
ＬＩＡＮＧ Ｘｉｕ－ｆｅｎｇ，ＹＡＮＧ Ｊｉａｎ－ｍｉｎ，ＬＩｘｉｎ，Ｌｌｊｕｎ
（Ｓｔａｔｅ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆ Ｏｃｅａｎ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｈａｎｇｈ￣Ｊｉａｏ Ｔｏｎｇ Ｕｎｉｖｅｒｓｉｔｙ，Ｓｈａｎｇｈａｉ ２０００３０，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｑｕａｌｉｔｙ ｏｆ ｗａｖｅ ｅｎｖｉｒｏｎｍｅｎｔ ｉｓ ｔｈｅ ｂａｓｉｓ ｆｏｒ ｓｔｕｄｙｉｎｇ ｔｈｅ ｎｏｎｌｉｎｅａｒ ｉｎｔｅｒａｃｔｉｏｎｓ ｂｅｔｗｅｅｎ ｗａｔｅｒ
ｗａｖｅｓ ａｎｄ ｆｌｏａｔｉｎｇ ｓｔｒｕｃｔｕｒｅｓ．Ｉｎ ｔｈｉｓ ｐａｐｅｒ，ａ ｎｕｍｅｒｉｃａｌ ｗａｖｅ ｔａｎｋ ｗａｓ ｕｓｅｄ ｔｏ ｓｉｍｕｌａｔｅ ｔｈｅ ｉｒｒｅｇｕｌａｒ ｗａｖｅ．
Ｆｌａｐ ｔｙｐｅ ｗａｖｅ——ｍａｋｅｒ ｗａｓ ａｄｏｐｔｅｄ ｆｏｒ ｇｅｎｅｒａｔｉｎｇ ｗａｖｅｓ ａｎｄ ｗａｖｅ ｒｅｆｌｅｃｔｉｏｎ ｗａｓ ｈａｎｄｌｅｄ ｐｒｏｐｅｒｌｙ ｔｈｒｏｕｇｈ
ｓｅｔｔｉｎｇ ｕｐ ｗａｖｅ ａｂｓｏｒｐｔｉｏｎ ｚｏｎｅ．ＶＯＦ ｍｅｔｈｏｄ ｗａｓ ｅｍｐｌｏｙｅｄ ｆｏｒ ｔｒａｃｋｉｎｇ ｔｈｅ ｆｒｅｅ ｓｕｒｆａｃｅ．Ｓｐｅｃｔｒｕｍ ａｎａｌｙ—
ｓｉｓ ｗａｓ ｐｅｒｆｏｒｍｅｄ ｏｎ ｔｉｍｅ ｔｒａｃｅ ｏｆ ｗａｖｅｓ ａｎｄ ｔｈｅ ｃａｌｃｕｌａｔｅｄ ｓｐｅｃｔｒｕｍ ａｇｒｅｅｄ ｗｅｌｌ ｗｉｔｈ ｔｈｅ ｔｈｅｏｒｅｔｉｃａｌ
ｓｐｅｃｔｒｕｍ．
Ｋｅｙ ｗｏｒｄｓ：ｎｕｍｅｒｉｃａｌ ｗａｖｅ ｔａｎｋ；ｄｙｎａｍｉｃ ｍｅｓｈ；ＶＯＦ ｍｅｔｈｏｄ；ＪＯＮＳＷＡＰ ｓｐｅｃｔｒｕｍ
１引 言
随着计算机技术的发展，建立数值水池应用于数值模拟的条件越来越成熟。与物模实验相比数值
水池具有经济性高、无触点测量流场、比尺效应小、消除物模中由传感器尺寸及模型变形等因素对流
场的影响、可获得较详细的流场信息等优点。但是目前数值水池的端部反射问题，尤其是造波板的二
次反射限制了数值水池的应用。通常的处理方法是在二次反射到达结构物之前停止实验，或模型宽度
设计得比水池的宽度ｔＪｑｌ￣多，对数值模拟意味着需要大量增加计算区域，这受到计算机内存及计算机
速度的限制。王永学【】１基于线性造波机理论在水池的一端设置了可吸收式造波边界条件，造波板的运
动除了产生行进波外，同时还产生一个抵消反射波的局部波动。Ｂｒｏｒｓｅｎ和Ｌａｒｓｅｎｔ２１提出了适合于边界
积分方程方法（ＢＩＥＭ）的源造波方法（ｓｏｕｒｃｅ ｇｅｎｅｒａｔｉｏｎ），即在计算域内设置一造波源，在源项两边同
时产生方向相反的两列波，源项处可透过波浪遇建筑物形成的反射波。
收稿日期：２ｏ０９—１２—２４
基金项目：深水油气开发工程试验技术（２００６ＡＡ０９Ａ１０７）；自然科学基金（５０７０９０７９）
作者简介：梁修锋（１９８１－），男，博士研究生，Ｅ－ｍａｉｌ：ｌｉａｎｇｘｉｕｆｅｎｇ＠ｓｊｔｕ．ｅｄｕ．ｃａ；
杨建民（联系人），男，教授，博士生导师，Ｅ—ｍａｉｌ：ｊｍｙａｎｇ＠ｓｊｔｕ．ｅｄｕ．ｃｎ。
４８２ 船舶力学 第１４卷第５期
为了消除边界处的反射，Ｌａｒｓｅｎ和Ｉ）ａｎｃｙｔ３￣提出了用海绵层方法吸收传递到边界处的波能，通过海
绵层消除波能的大部分，同时在出流边界处利用Ｓｏｍｍｅｒｆｌｅｄ条件，使未能衰减的部分波浪透过边界传
到域外。孙大鹏和李玉成［４１利用上述思想采用时域内对波面运动位置追踪的边界元方法，建立了一种
非线性波浪变形计算的三维数值模式，开发了一个三维非线性波的数值造波水池，进而对水池内的
Ｓｔｏｋｅｓ波进行了波浪变形计算。高学平等翻利用ＭＡＣ（ｍａｒｋｅｒ ａｎｄ ｃｅｌ１）法直接数值求解连续方程和Ｎ—
Ｓ方程，为模拟不规则波和长时段内连续造波及消除波浪遇结构物形成的反射，采用了源造波方法，在
开敞边界处采用了海绵层阻尼消波和Ｓｏｍｍｅｒｆｌｅｄ条件相结合的处理方式。
本文首先推导出用于不规则波浪模拟的摇板信号，结合ＦＬＵＥＮＴ软件中的用户自定义函数
（ＵＤＦ）功能和动网格（Ｄｙｎａｍｉｃ Ｍｅｓｈ）功能进行造波，同时利用ＦＬＵＥＮＴ软件中可以加入源项（Ｓｏｕｒｃｅ
Ｔｅｒｍｓ）的特征，在水池后段加入了消波源项，使得水质点速度在水池的后段逐步衰减为零，从而消除
了反射波对水池内波形的影响，为长时问的波浪模拟提供了保证，自由液面由两相流中的ＶＯＦ方法
捕捉。基于以上对不规则波浪的数值模拟进行了研究，得到了高品质的数值波浪。借助于ＦＬＵＥＮＴ软
件可以提供复杂网格的特性，该水池可以为复杂浮式结构物水动力特性的模拟提供数值波浪环境，具
有广阔的应用前景。
２数学模型
控制方程包拈连续性万程：
＋警＋ ＝０ （ ）
Ｎ—Ｓ方程
誓＋ｕ鲁 等一古誓＋ ａｘ＋ ＋ （２）
＋ 鲁 等一古等＋ －＋ Ｊ ｖ一 ＋ ｃ３
上式中，“、 分别为ｘ，ｙ方向的速度，Ｐ为流体的密度，／．ｔ为动力学粘性系数， ＝０， ＝＿ｇ， 和 为附
加源项。
本文采用ＶＯＦ方法捕捉自由液面，因此还需要引入另外一个控制方程，即流体体积分数的输运
方程
＋旦！． ＋ ）＿：０ ０ｔ ０ｘ
： （４）
，（ ，２）
文中ｑ＝ｌ表示空气相，ｑ＝２表示水相。 ，是一个标量，表示ｇ相流体在网格内占的体积分数
—．网格内ｑ相流体的体积 一 一网格的体积
本文所采用的Ｇｅｏｍｅｔｒｉｃ Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ方法是由Ｙｏｕｎｇｓｔ６１的方法发展而来的。该方法可以用于非
结构网格，计算精度较高。几何重构的第一步是计算交界面在网格内相对于网格中心的位置，第二步
由速度矢量和交界面的位置确定通过网格面的流体流量，第三步由前两步的结果计算每一个网格的
体积分数。它在每一时问步内根据网格内的流体体积分数将网格标记为流体网格（ ＝１），空网格（ ＝
０）或自由液面网格（０＜ ＜１）。并根据自由液面网格内的体积分数％重构自由液面的几何形状。然后利
用离散后的传输方程确定体积分数得到下一个时间步的自由液面形状，基于交错网格或同位网格求 第５期 梁修锋等：不规则波浪的数值模拟 ４８３
解不可压流体的Ｎ—Ｓ方程，得到下一个时间步流域内各点的速度值和压力值，如此循环而求出流体在
各个时刻的运动情形。本文默认体积分数 ＝０．５处为自由液面。
３几何模型
采用ＧＡＭＢＩＴ建立计算所需的几何模型，见图１。坐标原点０位于水池左端的初始自由液面处，
水池长４０ｍ，水深２ｍ，消波区长６ｍ，初始自由液面以上的高度为０．４ｍ，造波板绕Ｒ点做定轴转动。在
ｘ＝１５ｍ处设置了虚拟浪高仪器以监测波浪升高时历。流场网格划分如图２所示，在流场的大部分区
域都采用了结构网格，这样既可以减少网格数量，同时又可以得到光顺的自由液面，在自由液面附近
加密了网格以提高计算精度，见图２，自由液面附近网格大小为Ｏ．０１５ｍ，纵向为０．０３ｍ，网格总数约为
７５ ０００。 单位：ｍ
图１计算区域几何简图
Ｆｉｇ．１ Ｇｅｏｍｅｔｒｉｃ ｓｋｅｔｃｈ ｏｆ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎａｌ ｄｏｍａｉｎ
４摇板的运动规律 图２流场网格划分
Ｆｉｇ．２ Ｍｅｓｈ ｏｆｔｈｅｆｌｏｗｆｉｅｌｄ
本文的目标波波谱为ＪＯＮＳＷＡＰ谱川，水深为１２８ｍ，有义波高为１４．３ｍ，谱峰周期为１５．３ｓ，谱型参
数 ＝２．４，缩尺比为６４：１，计算时采用了模型值，最后将计算结果换算成了实型值。
ｓ ： Ｉ ５ ］ ㈣ ｝．．．．． ，，．．．Ｌ
『．． ， ４８４ 船舶力学 第１４卷第５期
式中ａ＝０．０６２ ４／［０．０２３ ０＋０．３３６３／一０．１８５（１．９＋ｙ） 】；ｔｏ．ｔｏＰ时ｏ－＝０．０７，ｔＯ＞ｔＯＰ时ｏ－＝０．０９。
零阶谱距（谱密度曲线下的面积）为：
％： （ ＝ 唧 者
因 舶＝４（ｍ。），所以 ＝ 。
为了得到不规则波需要对规则波进行波浪叠加阴，划分频率区间的方法有等分频率法和等分能量
法两种，本文采用等分能量法。定义累积谱为
Ｅ（０９）：ｆ Ｓ（ｔＯ） （８）
如按等分能量法分成 份，则分界频率ｔＯ 可按式
Ｅ（ｔＯ ）＝ｉＥ（∞）ＩＭ＝ｉｍｏ／Ｍ （９）
的解来决定。对于可以积分的谱，通过积分可得Ｅ（ｔＯ）和Ｅ（。。），然后将Ｅ（ｃ【，）和Ｅ（ｏ。）代人（９）式得ｔＯ 。
当谱形无法积分时，可通过数值积分计算累积谱ＪＺ（ｔＯ）并按等分能量法确定分界频率ｔＯ 。把选定的波
浪谱等分为Ｍ份后，求解得到每－－＋Ｎ率对应的波幅值ｑ＝、／ 珊（ｔＯｉ）ＡｔＯ ，式中ｓ聊（ ）为ｔＯ 处的波
谱值，当ｉ＝１时△ 产（ｔＯｍ—ｔＯ￡）／２，当ｉ＝Ｍ时△ （ｔＯｉ－ｔＯＨ）／２，当ｉ∈（２，３，… １）时△∞ （ｔＯ州— ｆ－１）／２；ｔＯｉ
为第ｉ个组成波的圆频率，本文中圆频率的区间为（０．２５－１．５），对于极高频和极低频的部分都不予考
虑，然后按等分能量法将频率区间分成Ｍ＝８５份，分解后得到的组成波的波幅见图４，这样对应的重复
周期约为８分钟，相应的波浪时历的时长也必须大于８分钟，不然所得波浪时历将不能完全反映目标
谱的统计特征。
为了给定造波板的运动规律，需要知道水力
传递函数，即不同周期时板前的波浪振幅％与造
波板振幅ｅ的关系。
为了适用于一般情况，假定造波板的运动为
推摇混合，距离池底距离ｄ ，ｄ２处板的振幅分别为
ｅ ，ｅ，，水面处振幅为ｅ，见图３，由势流理论可得［８１。 图３造波板布置图
Ｆｉｇ．３ Ｓｋｅｔｃｈ ｏｆ ｗａｖｅ—ｍａｋｅｒ ｓｅｔｕｐ
圆频率（ｒａｄ／ｓ）
图４组成波的波幅 图５组成波的随机初相位
Ｆｉｇ．４ Ｗａｖｅ ａｍｐｌｉｔｕｄｅ ｏｆ ｗａｖｅ ｃｏｍｐｏｎｅｎｔｓ Ｆｉｇ．５ Ｒａｎｄｏｍ ｉｎｉｔｉａｌ ｐｈａｓｅ ｏｆ ｗａｖｅ ｃｏｍｐｏｎｅｎｔｓ
ｒ￣（ｔＯ，ｅ ＝罢＝ （ ｓｉｎｈｋｄ￣ ｓｉｎ
＋ 壶壶 ｃｏｓ （１０
对于摇板，ｅ１＝ｅ，ｅ２＝０，ｄｌ＝ｄ，
＝０，则