岸式振荡水柱波能转换装置的数值模拟

宁德志，石进，滕斌，赵海涛

(1.大连理工大学海岸和近海工程国家重点实验室，辽宁大连116024;2.国家海洋局第二海洋研究所工程海洋学重点实验室，浙江杭州310012)

海洋面积占地球表面积的71%，它不仅为人 类提供航运、水产和丰富的矿藏，而且蕴含着大量的能量。海洋能主要包括波浪能，潮汐能，海洋温差能，海洋盐差能和海流能等，其中波浪能作为一种重要的可再生清洁能源具有分布广泛、能量密度大、采集结构简单等特点。振荡水柱式波能转换装置(oscillating water column，OWC)是目前各国最为重视、公认的最有前途、投入研究力量最大、建造使用最多的一种装置。该装置通过气室将波浪能转换为空气气流的能量，再通过空气透平将气流能量转换为电机转轴的轴功，最后通过发电机将转轴轴功转换为电能。岸式的振荡水柱转换装置是振荡水柱波能转换装置的一种，它具有结构简单，易于安装和维护等特点。关于 OWC的研究有很多。Evans等［1-6］在 OWC的二维和三维频域势流理论分析与数值计算中做了很多工作。时域势流理论方面，Koo等［7-8］采用考虑完全非线性自由水面边界条件的常数边界元方法并在气室内引入压强模型计算了一种带底坡的岸式的振荡水柱波能转换装置。Zhang等［9-11］应用基于N-S方程的粘性流模型对振荡水柱式波能转换装置进行了数值模拟。此外，Sarmento等［12-14］实验研究了振荡水柱式波能转换装置。

本文应用时域高阶边界元方法与压强模型建立了岸式振荡水柱波能转换装置的二维数值模型，研究了前墙尺寸，气室宽度及入射波要素对岸式振荡水柱波能转换装置能量转换效率的影响。

1 OWC时域数值模型

1.1 控制方程和边界条件

本文研究的岸式振荡水柱波能转换装置如图1所示，H代表水深，A、B、C、D分别代表气孔宽度、气室宽度、前墙壁厚和入水深度。在流体无粘、无旋、不可压缩的假定下，计算域内流体的速度可用速度势的梯度表示:

图1 岸式振荡水柱波能转换装置示意图Fig.1 Schematic of land-based OWC

本模型采用源造波方法造波，控制方程为泊松方程:

式中:q*(xs，z，t)=2vδ(x-xs)为源造波强度，v为流体质点的水平速度，本文给定为二阶Stokes波速度解析解，波浪产生位置x=xs。

在自由水面上，满足完全非线性运动学和动力学边界条件:

式中:η为自由水面的铅垂位移，P为自由水面上的压强，g为重力加速度。本模型采用混合欧拉-拉格朗日法更新自由水面边界条件，利用物质导数公式d/dt=∂/∂t+v·∇，并在源造波上游区域的自由水面布置人工阻尼层防止波浪反射，自由水面边界条件转化成以下形式:

其中:

式中:X0=(x0，0)为水质点初始位置，造波位置为x=0，α为阻尼层的阻尼系数，βλ为阻尼层的长度，λ为入射波波长，β为阻尼层宽度系数。

在岸式振荡水柱波能转换装置表面、水槽侧壁和水底等固定物体边界上，流体的法向速度为0，即

本模型在时域内求解，在静水面上初始条件为

1.2 气室内的压强模型

振荡水柱式波能转换装置气室内的气体压强与液面的运动互相影响，存在复杂的气液耦合现象，需要引入压强模型加以处理。本模型假定气室内的压强与孔口处的气体流速呈线性关系［7］:

或二次关系［7］:

式中:Ud(t)表示气孔处的气体流速，Cdm和Ddm分别为线性阻尼系数和二次阻尼系数，它们取决于气孔处透平的物理特性，例如威尔斯透平两端的气体压强差与流速就基本满足线性关系，式(8)中的绝对值符号是为了保证压强符号和孔口处气体流速方向的关联性。在孔口气体流速与气室内压强很小的情况下可以假设空气不可压缩，由此可得

1.3 能量转换效率的计算

入射波浪会使振荡水柱波能转换装置气室内的水体做上下振动，从而推动气室内的气体在气孔两侧往复运动，推动透平做功。波浪推动气体做功的速率可由下式求得

仿照波能流的推导，在一个周期内做平均可得

式中:z(t)表示气室内平均水面的速度。假定气室内的空气不可压缩且孔口处的气体流速呈正弦变化:

将式(8)、(12)代入式(11)可以推得

气室内的能量转换效率κ可由下式求得

式中:分母为入射波波能流，Ai为入射波波幅，Cg为群速度。

1.4 数值求解

在流体计算域内对速度势应用第二格林公式，可得到如下边界积分方程［15］:

式中:p=(x0，z0)为源点，q=(x，z)为场点，C(p)为固角系数，Γ为包括自由水面边界和固体边界的流体计算边界，Ω代表整个流体计算域。为了消除水底的积分，减小计算量，选取考虑水底镜像的格林函数:

式中:r1为p和q两点距离，r2为p和q关于水底的镜像之间距离。

本文用三节点的二次单元将计算域边界离散成一些曲线单元，每个单元可以通过数学变换变换成参数坐标下的等参元，在等参元内引入形函数可以插值得到单元内任一点的几何坐标和速度势。将源点取在各个节点上可以将离散后的积分方程化为线性方程组的形式，从而求解未知量。每一步的计算都认为当前时刻自由水面上的速度势和物面上的速度势法向导数是已知的，通过积分方程计算当前时刻自由水面上的速度势法向导数和物面上的速度势。在气室内应用压强模型后，下一时刻的自由水面上水质点的位置和速度势可通过四阶Runga-Kutta法计算，然后重新划分网格，继续应用积分方程计算自由水面上的速度势法向导数和物面上的速度势。这样计算周而复始，直到计算结束［16］。

2 OWC数值模型验证

气室内的压强模型已经成功被用于气动浮式防波堤的研究中，为了验证本文压强模型的准确性，对比计算了文献［7］研究过的一种气动防波堤的相应工况。计算水深30 m，气室宽度30 m，前、后墙壁厚10 m，入水深度10 m，气孔宽度1 m。入射波周期9～12 s，入射波波幅0.5 m。本研究中计算域长度取10倍波长距离，气室箱体布置在计算域中间，水槽两端均布置长度为1.5λ的阻尼层。通过数值收敛性试验，自由水面上每个波长长度布置20个网格，沿水深方向布置10个网格，时间步长Δt=T/60。图2(a)、(b)分别给出了入射波周期为12 s和10 s时，透射波高增大率γ随线性阻尼系数Cdm与二次阻尼系数Ddm变化的对比曲线。从图中可以看出，本文结果与文献［7］的数值结果吻合很好。

图2 不同压强模型的对比Fig.2 Comparison of the pressure models

图3是线性压强模型下，Cdm=150，入射波周期T=10 s时防波堤气室内压强的时间历程曲线。可以看到气室内压强波动周期与入射波周期相同。本模型气室内的压强在一定时间后趋于稳定，说明本模型可以消除二次反射的影响进行长时间模拟。

图3 气动浮式防波堤气室内压强时间历程，Cdm=150，T=10 sFig.3 Time history of the pressure in the air chamber of pneumatic floating breakwater，Cdm=150，T=10 s

为了进一步验证数值模型的准确性，将本文结果与文献［3，9，14］解析解、数值解和实验值进行了对比。振荡水柱式波能转换装置如图1，具体参数见表1。气孔宽度A=0.005 m，计算水深H=0.92 m，入射波幅0.04 m。计算域长度取5倍波长距离，通过开展数值收敛性验证，单位波长长度内划分30个网格，时间步长Δt=T/80。本文采用线性压强模型，经过试算，阻尼系数Cdm选定为 3.5。

表1 岸式振荡水柱波能转换装置结构尺寸Table 1 Geometrical parameters of the land-based OWC

图4给出了3种前墙尺寸的波能转换效率κ随无量纲入射波周期T*变化的对比结果，T*=从图中可以看出在高频区(T*＜6)，3种不同尺寸的模型计算中，本文结果和文献［9，14］吻合的都较好。在低频区(T*＞7)，本文结果与文献［9］的数值结果吻合较好，但与文献［14］的实验结果相比数值偏小，可能的原因是在实验中的造波板与结构距离为37.5 m，在入射波波长较长的情况下二次反射很快会影响到结构处的测量结果，造成能量转换效率被高估。本文结果与文献［3］解析解的趋势和共振频率基本一致，但数值上有一定差别，这是由于解析解对模型做了很多理想化假定，波能的最大转化效率可以达到100%，这在实际中是不可能达到的。总体来说，相对于文献［9］的两相流模型，本模型无论是可计算的频率宽度还是与实验结果的对比都较好，具有较高的计算精度。

图4 不同前墙厚度下，能量转换效率随无量纲入射周期T*的变化Fig.4 Energy conversion efficiency versus T*with different thickness of the front wall

3 结构尺寸对能量转换效率的影响

在水深一定的情况下，岸式振荡水柱波能转换装置的前墙入水深度、前墙厚度和气室宽度是影响其能量转换效率的主要因素。本文计算水深统一为0.92 m，入射波幅0.04 m。图5给出了在气室宽0.64 m，前墙厚度0.04 m时，结构在不同前墙入水深度下的能量转换效率对比曲线。从图中可以看到，整个结构存在一个明显的共振频率，在此频率上的能量转换效率最高，可以达到70%以上。随着入射波频率远离共振频率，振荡水柱波能转换装置的能量转换效率逐渐变小。在当前计算范围内，前墙入水深度的增大会导致结构的共振频率区间向低频区域移动，同时高频区(T*＜6)的能量转换效率变小，共振频率上的能量转换效率增大。前墙入水深度的变化对低频区(T*＞7)的能量转换效率影响不大，这是因为在低频区前墙入水深度相对于波长很小。

图6给出了气室宽度0.64 m，前墙入水深度0.15 m情况下，3种不同前墙厚度下的能量转换效率曲线。同样，在低频区(T*＞7)，前墙厚度的变化对能量转换效率的影响很小。随着前墙厚度的增大结构的共振区间向低频区的方向移动，在高频区和共振频率上，前墙厚度的增大会导致能量转换效率的减小。

图7给出了在前墙入水深度0.15 m，前墙厚度0.04 m，不同气室宽度下的能量转换效率对比曲线。改变气室的宽度会同时改变气室在低频区与高频区的能量转换效率，随着气室宽度的增大结构的共振频率向低频方向移动，且能量转换效率最大值会随气室宽度的增大而变小。气室宽度的增大会改善结构低频区的能量转换效率，降低结构在高频区的能量转换效率。

图5 不同前墙入水深度下，能量转换效率随T*的变化Fig.5 Energy conversion efficiency versus T*with different immersion depth of the front wall

图6 不同前墙厚度下，能量转换效率随T*的变化Fig.6 Energy conversion efficiency versus T*with different thickness of front wall

图7 不同气室宽度下，能量转换效率随T*的变化Fig.7 Energy conversion efficiency versus T*with different width of the chamber

4 结论

本研究利用时域高阶边界元理论与气体压强模型建立了岸式振荡水柱波能转换装置的完全非线性数值模型。从研究中可以发现岸式振荡水柱波能转换装置具有一个明显的共振区间，适当减小前墙入水深度、厚度与气室宽度都可以增大共振频率、提高高频区的能量转换效率;前墙尺寸的变化对低频区能量转换效率的影响十分微小;气室宽度的增加会提高低频区的能量转换效率，降低高频区的能量转换效率。和前人的势流模型相比，本模型采用高阶单元与源造波技术，精度更高，所需单元更少且消波问题很容易处理。和粘性流模型相比，本模型的计算效率更高，在合理的阻尼系数下同样可以与实验结果吻合良好。由于无法计算由于流动分离、涡旋等引起的粘性耗散并且尚未考虑空气的可压缩性，本模型存在一定的局限性，下一步工作需要建立更加贴近实际的压强模型并考虑粘性耗散的影响，同时对模型进行三维拓展，以达到更好的模拟效果。

［1］EVANS D.The oscillating water column wave-energy device［J］.IMA Journal of Applied Mathematics，1978，22(4):423-433.

［2］游亚戈.复杂地形下岸式波能装置的水动力学计算［J］.海洋技术，1993，12(1):20.YOU Yage.Hydrodynamic calculation of shore-mounted wave-power device at complicated topography［J］.Ocean Tecnology，1993，12(1):20.

［3］EVANS D，PORTER R.Hydrodynamic characteristics of an oscillating water column device［J］.Applied Ocean Research，1995，17(3):155-164.

［4］DELAUR Y，LEWIS A.3D hydrodynamic modelling of fixed oscillating water column wave power plant by a boundary element methods［J］.Ocean Engineering，2003，30(3):309-330.

［5］WANG D，KATORY M，LI Y.Analytical and experimental investigation on the hydrodynamic performance of onshore wave-power devices［J］.Ocean Engineering，2002，29(8):871-885.

［6］HONG D，HONG S Y，HONG S W.Numerical study of the motions and drift force of a floating OWC device［J］.Ocean Engineering，2004，31(2):139-164.

［7］KOO W，KIM M，LEE D.Nonlinear time-domain simulation of pneumatic floating breakwater［J］.International Journal of Offshore and Polar Engineering，2006，16(1):25-32.

［8］KOO W，KIM M-H.Nonlinear time-domain simulation of a land-based oscillating water column［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2010，136(5):276-285.

［9］ZHANG Y，ZOU Q P，GREAVES D.Air-water two-phase flow modelling of hydrodynamic performance of an oscillating water column device［J］.Renewable Energy，2011，41:159-170.

［10］刘臻.岸式振荡水柱波能发电装置的试验及数值模拟研究［D］.青岛:中国海洋大学，2008:70-107.LIU Zhen.Experimental and numerical investigation of oscillating water column wave energy converter［D］.Qingdao:Ocean University of China，2008:70-107.

［11］纪君娜，刘臻，纪立强.振荡水柱波能发电装置气室的三维数值模拟研究［J］.海岸工程，2011，30(2):7-13.JI Junna，LIU Zhen，JI Liqiang.3D numerical simulation for oscillating water column chamber in wave-power converter［J］.Coastal Engineering，2011，30(2):7-13.

［12］SARMENTO A，FALCAO A F O.Wave generation by an oscillating surface-pressure and its application in wave-energy extraction［J］.Journal of Fluid Mechanics，1985，150:467-485.

［13］TSENG R S，WU R H，HUANG C C.Model study of a shoreline wave-power system［J］.Ocean Engineering，2000，27(8):801-821.

［14］MORRIS-THOMAS M T，IRVIN R J，THIAGARAJAN K P.An investigation into the hydrodynamic efficiency of an oscillating water column［J］.Journal of Offshore Mechanics and Arctic Engineering，2007，129(4):273-278.

［15］NING D，TENG B，EATOCK T R，et al.Numerical simulation of non-linear regular and focused waves in an infinite water-depth［J］.Ocean Engineering，2008，35(8):887-899.

［16］NING D，TENG B.Numerical simulation of fully nonlinear irregular wave tank in three dimension［J］.International Journal for Numerical Methods in Fluids，2007，53(12):1847-1862.