用OpenFOAM实现数值水池造波和消波

� � 文章编号: 1005�9865( 2011) 03�0001�12 用OpenFOAM 实现数值水池造波和消波 查晶晶, 万德成 (上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室, 上海 � 200240) 摘 � 要:基于 OpenFOAM 求解器 interDyMFoam,开发实现数值水池造波(包括推板和摇板造波)和阻尼消波。所编写的造波边界 条件,可依据实验造波理论将各种造波模式植入其中, 从而实现各种波的数值造波。首先进行了线性波的数值造波实验, 通 过结合阻尼消波段的应用,造波水池可提供稳定的线性波。还进行了瞬时极限波和有限振幅的数值造波实验, 与实验数据或 同类数值结果吻合很好,进一步验证了的数值造波和消波方式的可靠性。 关键词: OpenFOAM; 动边界; 数值造波; 阻尼消波 中图分类号:TV139. 2� � � 文献标识码: A Numerical wave generation and absorption based on OpenFOAM CHA Jing�jing, WAN De�cheng ( State Key Laboratory of Ocean Engineering , School of Naval Architecture, Ocean and Civil Engineering , Shanghai Jiao Tong University, Shanghai 200240, China) Abstract: Based on the interDyMFoam package in an open source code software�OpenFOAM, piston�type and flap�type numerical wave�makers are set up to generate linear wave and finite�amplitudewaves, and theway of damping wave absorption is proposed. In order to absorb wave at the end of the wave tank, a damping zone funcuion is added to the solver. Several numerical experiments of the generation and absorption of linear wave and transient extreme wave in a numerical wave flume are performed to testify the effectiveness of the presented approach. Key words: OpenFOAM; moving boundary; numerical wave maker; damping wave absorption 收稿日期: 2010�12�09 基金项目:国家自然科学基金资助项目( 50739004, 11072154) ;海洋工程国家重点实验室研究课 ( GKZD010053�11) 和上海高校东方学者 特聘教授岗位计划项目资助 作者简介:查晶晶( 1983�) ,女,湖北人,硕士生,主要从事计算流体力学方向的研究。 通讯作者:万德成。E�mail: dcwan@ sjtu. edu. cn 随着高性能计算机的快速发展和数值模拟技术的改善提高, 船舶与海洋工程数值水池的开发和利用在 近几十年受到了越来越大的关注。船舶与海洋工程数值水池的开发和利用需要解决一些关键问题: 时域内 对非线性自由面运动的模拟; 数值造波和消波;时域内对非线性浮体运动的模拟等。与此相关的理论和数值 技术在过去三十年里得到了很大地发展, 其中数值造波和消波是建立数值水池所必需的最基本功能。 数值造波方法一般可分为物理造波和人工造波两类[ 1]。所谓物理造波,就是模拟实验室造波,即根据实 验造波理论来设定数值造波板的运动造出波[ 2�3] ;而人工造波, 是数值模拟中所特有的方法,实际物理水池中 无法实现, 比如边界输入法,是依据所造波的解析解或数值解,在计算域的边界设定波面和水质点速度[ 4�5] ; 或者在入口边界上, 设定一个模拟柔性造波板运动的速度分布[ 6�7] ;而源造波法,同样利用波浪理论的解 析解,在动量方程或连续方程中添加源项[ 8]。人工造波法避开了要处理造波板运动界面所带来的问题, 在技 第 29 卷第 3期 2011 年8 月 海 洋 工 程THE OCEAN ENGINEERING Vol. 29 No. 3Aug. 2011 术上容易实现且计算开销较小,因此许多学者采用此种方法造波;物理造波法在基于边界元法建立的完全非 线性势流模型( fully nonlinear potential flow model)中较容易实现, 而在基于 Navier�Stokes方程建立的粘性流水 池中, 往往都需要引入动边界处理技术。由于采用物理造波法可以充分利用实验造波理论和技术造出各种 波形, 模拟结果也方便与实验情况比较, 便于检验数值模型的可靠性,因此采用物理造波的数值水池更有望 接近实验水池, 从长远来看,有益于数值模拟与实验模拟两种研究手段更为密切的结合。 波浪在数值水槽一端产生并传播,到达水池另一端时会发生反射,反射回的波浪干扰入射波场, 使数值 实验无法继续下去。因此数值水池如同实验水池一样,需要采取消波措施或设置能让波浪完全透过出流边 界,使水槽末端不会发生波反射。出流边界(也称开边界)是数值计算中采用的人工截断边界,其作用是要让 流体自然的流出计算域边界而对计算域内的流场没有任何影响,这种理想效果实际难以达到。在提出的众 多类型中, Sommerfeld辐射边界的效果最好[ 9] ,但 Clement[ 10]研究指出, 由于该条件在时间空间上是局部的, 因此仅限于频率已知的规则入射波和长波情况,而不适用于完全非定常入射波情况。应用得更多的数值消 波方法:一种是设置阻尼消波区;另一种是模拟实验室主动吸收式消波。阻尼消波是由 Baker等人引入的, 其思想是在需要消波的区域对自由面边界条件或动量方程添加一阻尼耗散项, 研究发现消波区随着入射波 波长的增大需要相应增长,因此该方法对于长波的消波效率低下。主动吸收式消波, 是在水池末端设置推 板,让推板运动制造一个与反射波(若推板静止,波到达推板处会产生反射)波幅相同相位相反的波, 抵消掉 反射波。对于水池中的入射波,推板像透明的一样仿佛能让波自然地流出水池。Clement将阻尼消波与主动 吸收式两种方法耦合起来[ 10] , 阻尼消波区的长度保持一定主要是针对高频短波,消波板则针对低频长波,耦 合式方法不依赖于入射波频谱适应很宽的频带。而在粘性数值水池的消波处理上, Troch [ 11]等人在他们开发 的VOFbreak中也采用了主动吸收式消波技术,适用于规则和不规则波,而对于入射波到达结构物后引起的 反射波,通过数字过滤并与测得的速度信号叠加, 使其与入射波分离开来。另外还有一种简便易行的方法是 在水池后段设置逐渐稀疏的网格,利用数值粘性将波完全衰减掉,同时也可减少计算量, 不过此种方法很难 保证消波段的数值粘性不会影响到前面的波,因此网格稀疏比例不容易设置。实际采用何种消波方法, 要看 消波技术的可行性和消波效果。 基于 OpenFOAM求解器, 来开发和实现数值水池的物理造波和阻尼消波。OpenFOAM 是一款在 CFD领 域脱颖而出、可与商业软件媲美的开源 CFD软件。其开源性、源于 C+ + 面向对象所获得的良好扩展 性,稳定强大的底层类库,丰富的前后处理接口,并行化计算等诸多优点,有助于研究者在求解器、边界条件、 物理模型库、算法、几何网格等不同层面和深度上集中精力进行自己的研究开发并提高研究效率增进工作的 可延续性。另外,开源软件的开放模式大大促进了 OpenFOAM自身的壮大, 日益丰富的算法、模型库和求解 器使得研究者在解决问题时可有更多的选择, 并从中选出最优。OpenFOAM 的这些优势无疑给数值水 池的建立带来了极大的方便。因此,选择了OpenFOAM 作为数值水池的开发平台。由于OpenFOAM�1. 5中已 经植入了动网格技术[ 12�13]的两相流求解器 interDyMFOAM,因此可基于 OpenFOAM提供的动网格边界类设置 动边界的运动, 并调用网格运动求解库,模拟推板或摇板造波。 1 � 数值造波模型 数值造波模型是基于 OpenFOAM 中的基本两相流模型 interDyMFoam 求解器,结合编写的造波边界和阻 尼消波。即如何通过 C+ + 类继承的方式创建造波所需的边界条件, 如何在顶层求解器代码中为动量方程 添加阻尼消波项,以及如何结合修改过的 interDyMFoam求解器和造波边界建立二维数值造波水池。 1. 1 � 基本两相流模型 1. 1. 1 � 控制方程 在OpenFOAM 中, interDyMFOAM是求解一类不可压缩、不可渗透、等温的两相流问题的求解器。它基于 N�S方程, 使用VOF法捕捉界面,并结合了动网格技术。对于不可压缩粘性流体流, 若速度场用 u 表示,由 质量守恒有: � � u = 0 ( 1) � � 由动量守恒有: 2 海 � � 洋 � � 工 � � 程 第 29卷 ��u� t + � � ( �uu ) - � � �� u - �g = - �p - f � ( 2) 式中: �表示整个流场的密度, �表示动力粘性系数, g 表示重力加速度, p 表示压力场, f � 表示自由面上的 张力,且由下面式子给出: f � = ��( x ) n (3) n = �� | ��| (4) �( x ) = � � u ( 5) � � 由于液体在自由面上受到表面张力的影响, 这里采用了连续表面力模型。其中, �是张力系数, 表示自 由面每增加单位面积所需做的功, �表示界面平均曲率, n 表示界面的单位法向量, 用 f �表示表面张力沿界 面法向的一个分量, 此力的作用是平衡界面两边的压力差, 它只在界面处起作用,在非界面处其值为零。 在对自由面的处理上,根据VOF 法, 引入体积分数函数 �: �( x , t ) = 1, � 流体 1 0, � 流体 2 0 < �< 1, � 两种流体的交界面 (6) � � 这里将流体的过渡区域厚度记为 �, 体积分数函数满足: D� D t = ���t + u � ��= 0 (7) � � 在求解时, 将两种流体当作一种流体来处理, 它满足方程( 1)、( 2) ,其中该流体的物理属性密度 �和粘性 系数 �由两种流体各自的密度( �1, �2)、粘性( �1, �2)和体积分数函数来表示: �= ��1+ (1 - �) �2 � = ��1+ (1- �) �2 (8) 1. 1. 2 � 数值方法 在二阶精度有限体积法离散框架下,所有物理量安排在单元体中心上( colocated) , 对控制方程在每个单 元上做时间空间积分,各积分项的数值估算需要通过插值。OpenFOAM中提供了多种插值格式,这里时间积 分采用隐式欧拉格式,体积分中对流项用 Guass limitedlinearV 1, 粘性扩散项用 linear corrected(对非正交网格 做显式校正) ,其余用线性插值。 为了求解体积分数函数 �的迁移问题,Weller引入了额外的人工压缩项, 将式( 3)改为: � �� t + � � ( �u) + � � ( �( 1- �) u�) = 0 (9) 式中: u�是一个适于压缩界面的速度场,这一项只对界面区域产生影响[ 14]。体积分数方程中的对流项采用 Gauss vanleer格式,为构造界面而添加的一项采用 Gauss interfaceCompression格式。网格运动的求解基于拉普 拉斯方程;速度压力场采用分离式( segregated)算法 PISO;离散后的线性方程组采用预处理共轭梯度法求解。 1. 2 � 造波边界及消波阻尼 1. 2. 1 � 自定义造波边界 OpenFOAM在组织架构上充分利用了 C+ + 的抽象分层思想, 它的边界条件类型都是由相关基类( fv� PatchField, pointPatchField)派生出来的(表 1)。OpenFOAM 中的网格运动边界条件基类是pointPatchField, 流场 边界条件基类是 fvPatchField,所有的动网格边界条件类型和流场边界条件类型分别是由这两个基类派生出 来的, 因此在设定边界条件类型时, 要从它们的派生类中选取。这里的数值造波是基于实验造波理论,通过 设定边界的运动来模拟推板和摇板造波。基于OpenFOAM 提供的能使网格边界作往复运动的边界条件类型 ( angularOscillatingVelocity, angularOscillatingDisplacement, oscillatingDisplacement, oscillating Velocity) , 植入各种 造波公式(如公式( 16)、( 17) ) ,构造了新的推板造波边界 pistonMaker Displacement类(摇板造波边界的构造方 法也与此类似) , 如表 2所示, 编译及使用可参考相关资料。 3第 3期 查晶晶, 等:用 OpenFOAM 实现数值水池造波和消波 表 1� OpenFOAM中的边界条件 Tab. 1 � The boundary conditions in OpenFOAM 边界条件类型 流场边界 网格边界 基类 一级派生类 次级派生类 fvPatchField< Type> fixedValue fixedGradient mixed � advective movingWallVelocity slip pressureInletOutletVelocity � pointPatchField< Type> fixedValue basicSymmtry mixed � oscillatingVelocity oscillatingDisplacement angularOscillatingDisplacement angularOscillatingVelocity � 表 2 � 推板造波边界 Tab. 2� The pistonmaker boundary pistonMakerDisplacement类的构造 class pistonMakerDisplacementPointPatchVectorField: public fixedValuePointPatchField< vector> { private: vector stroke ; scalar omega ; scalar depth ; scalar height ; scalar k ; scalar w ; word waveType ; const static scalar g = 9. 81; public: TypeName(� pistonMakerDisplacement�) ; pistonMakerDisplacementPointPatchVectorField ( const pointPatch&, const DimensionedField< vector, pointMesh> &, const dictionary& ) ; scalar k( ) ; scalar w( ) ; virtual void updateCoeffs( ) ; } 推板造波边界继承自 fixedValuePatchField 数据成员: 推板冲程; 推板运动频率; 水深; 目标波高; 波数(满足线性频散关系) ; 目标波高与推板冲程之比; 波的类型; 成员函数: 将此边界条件类型取名为: pistonMakerDisplacement 主要的构造函数, 对象初始化时,部分信息通过字典 从算例中读取 迭代求解 �2= kgtanhkd ,由 �, g, d 求得 k ; 按公式( 14) 求得 w ; 每一时间步内更新边界处的值,各种造波公式(如公 式( 13)、( 17) )在此植入, 具体的造波公式在下一章 给出。 1. 2. 2 � 添加阻尼消波项 为了避免水池右端墙面处的波浪反射,采用阻尼消波技术, 在动量方程中加一阻尼项 ��u : 4 海 � � 洋 � � 工 � � 程 第 29卷 ��u� t + � � ( �uu) - � � ��u + ��u = �g - �p - f � ( 10) � � 其中定义 �: �= x - x 0 x 1- x 0 �1, � x 0 < x < x 1 0, � � � � � 0 < x < x 0 ( 11) 式中: x 0是消波段起始点, x 1是消波段终点即为水池末端, �1 是消波系数。 具体在OpenFOAM 中实现时,在 interDyMFOAM的求解器程序文件 UEqu. H 里, 增加以下语句并对速度 方程进行修改, 见表 3。 表 3 � 增加阻尼消波功能 Tab. 3 � The function for wave damping 对 interDyMFoam 中UEqu.H 增加的语句 说明 dimensionedScalar theta0 ( � theta0� , dimensionSet( 0, 0,�1, 0, 0, 0, 0) , 0 ) ; dimensionedScalar theta1 ( environmentalProperties. lookup(� dampCoffe�) ) ; �Info< < � damp cofficient is�< < theta1. value( ) < < endl; dimensionedScalar � x1= max( mesh. Cf( ) . component( vector: : X) ) - min( mesh. Cf( ) . component( vector: : X) ) ; �Info< < � the length of the wave tank is � < < x1< < endl; dimensionedScalar x0 ( environmentalProperties. lookup(� dampZoneStartPoint� ) ) ; �Info< < � the damp zone start point � � � � is� < < x0. value( ) < < endl; volScalarField thetaField ( � thetaField� , min( max( theta1* ( mesh. C( ) . component( vector: : X)�x0) / ( x1�x0) , theta0) , theta1) ) ; �Info< < thetaField. dimensions( ) < < endl; fvVectorMatrix UEqn ( � � � � fvm: : ddt( rho, U ) � � � � + fvm: : div( rhoPhi, U) � � � � - fvm: : laplacian( muf, U ) � � � � - ( fvc: : grad( U ) & fvc: : grad( muf) ) � � � � + U* rho* thetaField ) ; 定义一个有量纲的标量 �0,取值为 0,单位为 s- 1; 定义消波因子�1, 其量纲和取值是通过字典对象 en� vironmental Properties读入的, 因此在建立算例时需要 在 waveTank/ constant/ environmentalProperties 文件中设 置该量的量纲和取值, 其量纲跟 �0 一样, 取值可根 据计算结果来调整 ; 计算二维水池的长度; 定义阻尼消波起始点, 量纲和取值通过字典对象 en� vironmental Properties读入; 按照式( 11)定义 �; 检查 �的量纲是否设置正确; 在动量方程左端增加 ��u 项; 5第 3期 查晶晶, 等:用 OpenFOAM 实现数值水池造波和消波 1. 3 � 数值造波水池的建立 1. 3. 1 � 网格划分 本工作的数值水池模型采用结构化网格,它的生成直接利用 OpenFOAM中自带的 blockMesh。为了提高 波面的捕捉精度,对波面处的网格做垂向加密,而远离波面网格逐渐变稀。造波水池如图 1所示。 1. 3. 2 � 边界条件的设置 在建立算例时,需要同时设定网格边界和流场边界。网格边界在waveTank/ 0/ pointDisplacement文件中 设定。例如要使方形水池的左边界模拟推板运动, 可在 pointDisplacement文件中将左边界设为 pistonMak� erDisplacement类型, 并设置造波类型和必需的造波参数。整个网格的上下边界要设置成 slip,即左端的动边 图 1� 造波水池示意 Fig. 1� Definition sketch of the wave tank 界可在上下边界之间滑移。网格右边界要设成 fixed� Value类型,且速度为 0, 若设为 zeroGradient, 会导致整 个水池动起来。在设定流场边界时,在算例文件 wave� Tank/ 0/U中,模拟造波板运动的动边界处的流场边界 设为 movingWallVelocity 类型, 它能保证动边界处的流 体质点跟着动边界一起运动。水池底部一律设为 fixedValue,且速度为零,即无滑移条件。 1. 3. 3 � 动网格库的调用 在waveTank/ constant/ dynamicMeshDict文件中需要设置动网格信息,这里选取网格不发生拓补变换的 dy� namicFvMesh, 并进一步选取适合网格小幅度运动的 dynamicMotion SolverFvMesh, 然后调用网格运动库 libfvMo� tionSolvers. so,并设置相关参数。 1. 3. 4 � 波面的提取 后处理中, 为得到波面位置,对体积分数值沿水池垂向作积分,然后可以简单的换算出波面高度。此方 法的优点是所获得的波面位置不受界面重构质量的影响。具体的程序改写自 patchAverage( OpenFOAM中已 有的用于提取边界处平均值的程序)。 2 � 数值收敛测试 依据线性造波理论[ 15] ,用推板造线性波时,若要产生一个正弦波面 �= - H 2 sin( kx - �t ) ( 12) 推板的位移应设为 �( t ) = S 2 cos�t ( 13) � � 并且造波板的冲程 S 与目标波波高H 有如下关系: H S = 2( cosh2kh - 1) sinh2kh + 2kh ( 14) 其中, h为水深, k 为波数。 推板的线性造波运动式( 13)以及关系式( 12)、( 14)已植入上一节介绍的造波边界 pistonMakerDisplacement 中。下面是调用此边界造线性波的效果分析。 由于实际求解过程中,网格和时间步长的选取是直接影响离散误差的因素,数值解的精度往往表现出对 它们的依赖性, 因此需要测试出这些因素对结果的影响程度或影响规律。数值收敛测试可以为我们合理选 取网格和时间步长提供指导, 有助于我们用尽可能少的计算获得足够的精度。在数值收敛测试中,建立一个 长30 m,高 1. 2 m,水深为 0. 6 m的二维矩形水池模型。利用 pistonMakerDisplacement造波边界在此水池中产 生一个波高为 0. 02 m, 周期为 1 s的线性波。在水池 20~ 30 m的区域设置消波阻尼,消波系数取 1. 7。对时 间步长,选取了 0. 008 s、0. 005 s、0. 002 s、0. 001 s 四种情况; 对水池水平方向单元数目, 选取了 600、800、 1 000、1 500四种情况;对水面区域在垂向波高范围内的单元数,选取了 10、20、30三种情况。整个水池的水 平方向网格均匀分布,垂向网格在波面波高范围内均匀分布而在非波面区域网格远离波面逐渐变稀疏。由 于在 40 s时刻,水池左端最初产生的波已通过消波区的有效衰减, 到达水池右端时归于平静,因此选取 40 s 6 海 � � 洋 � � 工 � � 程 第 29卷 时刻整个水池的波面曲线作为比较对象,观察它在不同的时间步长, 不同的纵向网格疏密度和不同的垂向网 格疏密度下的变化趋势。这里需要说明的是,对于所有波面(空间和时间序列)位置的获取,采取的方法是对 体积分数沿水池垂向作积分, 因此这里可以排除界面重构的质量对所获波面的影响。 图2显示的是用四种时间步长算得的 40 s时刻水池波面曲线,其中水池水平方向都取600个单元, 波面 处垂向都取 20个单元, 可以看到算得的波面随着时间步长的增大有降低的趋势。说明时间步长取得越大, 所产生的数值粘性造成的波面衰减也变大,因此在接下来的计算中时间步长要选得稍小一些。 图 2� 四种时间步长下得到的 40 s时刻水池中的波面曲线 Fig. 2� Wave profile along wave tank at the moment of 40 s, horizontal cells: 600, vertical cells around the surface: 20, four types of time step: 0. 008, 0. 005, 0. 002, 0. 001 图3显示的是在水平方向采用四种不同单元数目算得的40 s时的水池波面曲线。其中计算时间步长都 取0. 002, 波面区域垂向单元数目都取 20,可以看到,算得的波面其幅度随着纵向网格的加密而增大, 当纵向 单元数增加到 1 500时,波面幅度达到 1 cm,此时与目标波幅最吻合, 说明此时纵向网格的稀疏度引起的误 差已经足够小了。由此看来, 纵向网格越密, 所获波的波幅越接近目标波幅。 图 3� 四种不同单元网格数目下算得的 40 s时刻水池波面曲线 Fig. 3� Wave profile along wave tank at the moment of 40s, time step: 0. 002, vertical cells around the surface: 20, four types of horizontal cells: 600, 800, 1 000, 1 500 图4显示的是垂向波面区域采用三种不同单元数目算得的 40 s时刻水池波面曲线。其中水平方向单元 数目都取 800,计算时间步长都取 0. 002,可以看到波面区域垂向网格越密造出的波波幅反而越偏小, 反之在 波面区域垂向取十个单元的情况下,波幅最接近目标波幅, 此时由垂向网格疏密度引起的误差已足够小。 由此看来,其他条件一定,波面垂向网格越稀,所获波的波幅越接近目标波幅,这一点值得关注。 从图 3到图 4的比较中发现, 纵向网格、垂向网格和时间步长在取值 1 500、20、0. 002和 800、10、0. 002这 两种情况下的结果都很好,因此很自然的将两个结果作一比较。正如图 5中所示它们非常接近,即用前面一 种较密的网格与用后面一种较粗的网格算得的结果非常接近。进一步还可以发现,尽管两套网格疏密程度 差异较大,但是网格单元的长宽比例很接近(大概为20比 1)。结合前面图 3和图 4各自的分析:存在纵向网 格越密垂向网格越稀获得波的波幅越接近目标波幅的趋势。由此可以得出结论:计算结果的精度主要受波 面区域矩形单元的长宽比影响, 网格单元越细长, 引入的数值粘性越大, 由此造出的波与目标波波幅相差 太大。 从上面的数值收敛测试可以看出以下几点: 7第 3期 查晶晶, 等:用 OpenFOAM 实现数值水池造波和消波 1)网格对造波精度影响较大。通过分析结果随着纵向和垂向网格疏密度的变化趋势, 得知由波高区域 单元长宽比引起的数值粘性会导致造波波幅偏差。通过以上多组计算发现,控制网格长宽比是提高结果精 度的关键,在此规律下,实际计算中利用稀疏的网格不但可降低计算开销而且有可能获得同样较好的结果; 2)对比最好的结果(图 5)和最差的结果(图 2)可以看到,网格引起的误差是非常可观的,尤其是针对波 幅较小的波,应该特别注意由此引起的波面非物理衰减。 图 4� 三种不同垂向单元数目下算得的 40 s时刻水池中的波面曲线 Fig. 4� Wave profile along wave tank at the moment of 40 s, time step: 0. 002, horizontal cells: 800, three types of vertical cells: 10, 20, 30 图 5 � 两组不同参数组合情况下计算结果的比较 Fig. 5� Comparison of results of two conditions 对于本套数值求解模式,在时间步长和网格疏密度的选取上, 提供如下的建议:实际计算中可采用变步 长模式,限定柯朗数小于 0. 2,求解程序会根据柯朗数的限定动态调整时间步长;在划分网格时,针对不同的 造波情况,最好先通过试算考察波面单元长宽比对造波精度的影响程度,以获得合理的波面区域网格单元长 宽比。比如本例造波测试中发现波面长宽比控制在 20�1以下,波高范围内垂向单元取 10个左右,波长范围 内纵向单元取 40个左右,可获得较好的造波效果。本数值模式中造波精度对网格的依赖性较大, 推测是源 于在离散求解VOF 法的体积分数方程( 9)的过程中引入的数值粘性过大导致的,或许可在该方程的离散格 式上想办法加以改进。 3 � 数值造波实验 3. 1 � 线性波 为了检验线性波的造波效果,选取了上一节数值收敛验证中的一组计算结果( 800/ 10/ 0. 002)来做分析, 较为理想的造波效果是要能在数值水池里持续不断的产生与理论解吻合较好的波。为了观察水池各段波面 运动的稳定性, 在距造波边界一个波长,两个波长,四个波长,六个波长,八个波长,十个波长, 十二个波长,十 三个波长(已接近消波段)远处设置探测点,提取波面历时曲线(图 6)。为了观察波在水池中的传播性态、消 波段的衰减效果和长时间内波在水池中的可重复性, 每隔五个周期依次提取了 30 s、35 s、40 s、45 s各个时 刻水池中的波形曲线(图 7)。 图6的历时曲线显示了水池中各个位置处的水面从初始平静状态然后开始振荡到进入稳定波动状态的 过程, 可以看到各探测点处的波面能够长时间稳定波动,另外仔细观察可以发现这一波陡(波高波长比)为 0. 013, 波高水深比为 0. 033的数值波(实线)与线性解析解(虚线)相比,波峰重合,波谷略有抬升, 说明这种 波实际应该还是有微弱的非线性效应存在,总的来说吻合得较好。 8 海 � � 洋 � � 工 � � 程 第 29卷 图7中可以观察到水池前十几米内的波段与理论解(虚线)吻合得很好,波到达阻尼消波段 20 m 处后经 过五六个波长距离的逐渐衰减之后消失殆尽, 还可以看到水池中的波面形态在 35 s、40 s、45 s时刻几乎一 样。另外,为了更清楚显示消波效果,将不加消波段的情况与之作了比较, 图 8给出未经消波段衰减的波经 水池右壁反射后的情况。综合图 6至图 8来看, 所造线性波可以在水池中长时间稳定的生成和传播。 图 6 � 距离造波边界不同距离处波面的时间序列 Fig. 6� Time history of wave elevation at the distance of eight/ ten/ twelve/ thirteen wave lengths from the wave maker 图 7 � 不同时刻水池内的波面曲线 Fig. 7� Wave profile along wave tank at moments of 30, 35, 40, 45 seconds 9第 3期 查晶晶, 等:用 OpenFOAM 实现数值水池造波和消波 图 8 � 有无消波段 45 s 时刻水池波面 Fig. 8� The effect of wave at 45 s without damping zone and damping zone 3. 2 � 瞬时极限波 为检验造波效果,对 Daniel T Cox[ 16]为研究甲板上浪现象在实验水池中造出的一段瞬时极限波进行数值 模拟。根据Daniel实验中的摇板驱动信号(如图 9)设置动边界的运动: �= A 1sin( �1( t - �1) ) , � t 1 � t < t 2 A 2sin( �2( t - �2) ) , � t 2 � t < t 3 ( 15) 式中: �表示摇板的角位移, A 1, A 2 是前后两段做正弦摆动的幅度, �1, �2是相应的角频率。此种造波运动 是依据( 15)对边界条件类 angularOscillat ingDisplacement中的 updateCoeffs( )函数里的角位移公式进行修改后实 现的。由于不清楚具体的实验细节,比如摇板转轴的位置,模拟中将摇板转轴设在池底处, 而实验中摇板的 冲程也不可知, 因此在实际模拟中根据电压信号图保证在前后两个时间段摇板在水面处的最大位移为 1�2, 而摇板冲程则视模拟结果进行调整。模拟中发现若摇板的运动从 1. 25 s开始, 造出的波与实验相比有一点 相位差;若将 t1 调整为 1. 5 s, 模拟结果与实验吻合的较好。分析原因可能是实验中摇板的运动较电压信号 稍稍滞后,此时 t 2取值 3. 5 s, t 3取值 7. 25 s, �1和 �2 取 0. 5。 从图 10对比曲线来看,这段瞬时极限波开始在 4. 5 m 处出现两个较大的波峰,传到 7 m左右变换成一 个更大的波峰, 到 10 m左右又分解成两个波峰。模拟结果与实验吻合的比较好,能反映出这段波在传播过 程中的演变特点,不过仍然可以发现模拟出来的波始终达不到实验中的陡度, 即表现为小幅度波的波幅偏 大,大幅度波则始终略低于实验中的波峰高度。 图 9 � 实验中的造波板电压信号 Fig. 9� The input electric signal for wave�maker of experiment � 图 10� 瞬时极限波的数值模拟结果与实验的对比 Fig . 10� The comparison of time histories of wave elevation of presented results with the experimental ones 3. 3 � 有限振幅波 根据Madsen [ 17]的造波实验和理论,利用推板线性造波理论造有限振幅波时,造出的波含有二阶自由波 成分, 其相速与主波成分不一样,在传播过程中会导致波形不稳定。Madsen对推板运动公式进行修正后,经 实验验证, 在 Ur< 8�2 3 的范围内, 修正后的推板运动能有效消除二阶自由项成分的干扰, 从而产生稳定的波 形。若让推板仅作正弦运动, 推板的位移公式: �( t ) = - �0cos�t ( 16) � � Madsen修正后的二阶造波板运动位移公式: 10 海 � � 洋 � � 工 � � 程 第 29卷 �= �1 + �2 = - �0[ cos�t + a2h0n1( 34sinh2 kh0 - n1 2 ) sin2�t ] ( 17) 其中, n1 = 1 2 (1+ 2kh0 sinh2kh0 ) , �0 = an1tanhkh0 ( 18) 式中: �0 表示推板冲程; h0表示水深; k 表示波数,满足线性波频散关系; a 表示波幅。由此产生的是 stokes 二阶波,在传播过程中能保证稳定的波形: �= H2 cos( kx - �t ) + �8 H 2 L coshkh sinh 3 kh ( cosh2kh + 2) cos[ 2( kx - �t ) ] ( 19) � � 二阶造波式( 17)已植入 pistonMakerDisplacement边界条件中。调用时指明造波类型为二阶波并设置必需 的造波参数即可。为了跟Madsen的实验结果以及Huang C J的数值模拟结果[ 18]进行比较,水深、推板运动周 期和冲程等按实验中的设置: h = 38 cm, T= 2. 75 s, �0= 6. 1 cm, Ur= 27. 2。分别利用一阶线性造波公式和修 正后的二阶造波公式造一个有限振幅波,取 21. 8 s时刻的水池波面进行观察。图 11中上方是推板按照线性 造波理论运动公式造出的波, 此时造波板的运动幅度比之前造微幅线性波时要大,板子运动产生的波从水池 左端传到右端, 可以看到在水池30 m 至 35 m 波前处的一个波形还是看起来波峰陡波谷平坦的规则样子,但 是此刻远离波前的波形已经明显发现了变化,要么波谷变成陡三角( 15 m至 20 m处) , 要么波谷处出现二次 波峰( 10 m至 15 m处) ,可以推知此时的波已经不是单频成分,并且次级成分的波已经造成整个波在空间上 不再保持稳定波形。再看图 11中下面一波形曲线,此时波形在空间上保持了稳定,说明给造波板运动增加 一个次级的倍频成分,可以有效消除次级波对主波的影响, 从而能够产生稳定的有限振幅波。 图 11� 推板线性造波运动和二阶造波运动得到的 21. 8 s 时刻的波面曲线 Fig. 11� Wave profiles at t= 21. 8 s by piston linear and second order wave�makers 对造出的有限振幅波取其稳定后一个周期内的波面运动情形进行观察,图 12是分别于水池两个探测点 处探测到的波面历时曲线,实线表示此处模拟结果, 虚线是 Huang C J的数值模拟结果[ 18] ,空心圈是 Madsen 的实验结果,实心圈是解析解(式( 19) ) ,可以看到此处模拟结果与Huang C J的数值结果非常接近。 图 12� 距造波边界不同距离处的波面在稳定后一个周期内的历时曲线 Fig. 12� The time histories of wave elevation at different x locations 11第 3期 查晶晶, 等:用 OpenFOAM 实现数值水池造波和消波 4 � 结 � 语 基于 OpenFOAM平台,通过OpenFOAM的边界条件类和求解器层面的代码改写,实现了二维数值水池的 数值造波和消波。在线性波的数值收敛测试中发现,在此套数值模式下,要想获得好的计算结果需要特别注 意网格的合理划分, 波面区域网格的长宽比不能过大, 否则会引起严重的数值粘性,从而使造出的波出现衰 减的假象, 这种因数值粘性造成波的衰减往往大于物理粘性造成的衰减。在合理的选取网格、时间步长、阻 尼消波区长度和消波系数后, 所建立的数值水池可以长时间稳定的产生线性波。随后对实验室的一段瞬时 极限波的模拟结果显示与实验情况吻合得比较好,造出的有限振幅波与 Madsen 的实验结果以及Huang C J 的数值模拟结果非常接近,证明基于OpenFOAM 的数值造波的效果非常好。 通过工作, 还表明在使用开源代码工具 OpenFOAM时,有利于对所研究问题作集中深入的开发, 容易将 自己所需的功能融入到整个OpenFOAM系统构架中, 可以大大提高研究工作的延续性和合作共享性。在目 前二维造波工作的基础上,将进一步开展三维数值水池造波消波和应用方面的工作。 参考文献: [ 1] � Tanizawa K. The state of the art on numerical wave tank[ C] �Proceedings of 4th Osaka Colloquium on Seakeeping Performance of Ships. 2000: 95�114. [ 2] � Grilli S T , Horrillo J. Numerical generation and absorption of fully nonlinear periodic waves[ J] . Journal of Engineering Mechanics, 1997, 123( 10) : 1060�1069. [ 3] � Huang C�J, Zhang E�C, Lee J�F. Numerical simulation of nonlinear viscous wave fields generated by a piston�type wavemaker[ J] . Journal of Engineering Mechanics, 1998, 124( 10) : 1110�1120. [ 4] � Lin PengZhi, Liu P L�F. A numerical study of breaking waves in the surf zone[ J] . Journal of Fluid Mechanics, 1998, 359: 239�264. [ 5] � Baudic S F, Williams A N, Kareem A. A two�dimensional numerical wave flume�part 1: nonlinear wave generation, propagation, and absorption[ C] �The 19th International Symposium and Exhibit on Offshore Mechanics and Arctic Engineering. 2001. [ 6] � 王永学. 无反射造波数值波浪水槽[ J] . 水动力学研究与进展: A 辑, 1994, 9( 2) : 205�214. [ 7] � 齐 � 鹏, 王永学. 三维数值波浪水池技术与应用[ J] . 大连理工大学学报, 2003, 43( 6) : 825�830. [ 8] � Ge Wei, James T K , Amar S. Generation of waves in Boussinesq models using a source function method[ J] . Coastal Engineering, 1999, 36: 271�299. [ 9] � Hatayama S. Comparison of effectiveness of four open boundary conditions for incompressible unbounded flows[ R] . Techanical Report of National Aerospace Laboratory, 1998. [ 10] � Clement A. Coupling of two absorbing boundary conditions for 2D time�domain simulations of free surface gravity waves[ J] . J. Com� put. Phys. , 1996, 126: 139�151. [ 11] � Troch P, Rouck J D. An active wave generating�absorbing boundary condition for VOF type numerical model[ J] . Coastal Engineering, 1999, 38: 223�247. [ 12] � Jasak H, Tukovic Z. Automatic mesh motion for the unstructured finite volume method[ J] . Transactions of FAMENA, 2007, 30: 1� 18. [ 13] � Jasak H, Tukovic Z. Dynamic mesh handling in OpenFOAM applied to fluid�structure interaction simulations[ C] �Presented at the V European Conference on Computational Fluid Dynamics. 2010. [ 14] � Rusche H. Computational Fluid Dynamics of Dispersed Two�Phase Flows at High Phase Fractions[ D] . London: Imperial Colledge of London, 2002. [ 15] � Dean R G, Dalrymple R A. Water Wave Mechanics for Engineers and Scientists[ M ] . New Jersey: Prentice�Hall Inc. , 1984. [ 16] � Cox Daniel T , Ortega J A . Laboratory observations of green water overtopping a fixed deck[ J] . Ocean Engineering, 2002, 29: 1827� 1840. [ 17] � Madsen O S. On the generation of long waves[ J] . Journal of Geophysical Research, 1971, 76: 8672�8683. [ 18] � Dong CM , Huang C J. On a 2�D numerical wave tank in viscous fluid[ C] �Proceedings of the Eleventh ( 2001) International Offshore and Polar Engineering Conference. 2001: 146�155. 12 海 � � 洋 � � 工 � � 程 第 29卷