目前我国的潮流能和波浪能技术发展较为成熟^[1]，出现了各种各样形式的海洋能发电装置。如基于双棘轮式的新型振荡浮子式海洋波浪能发电装置^[2]，将振荡浮子和振荡气室结构相结合的振荡水柱式波浪能发电装置^[3]，以及振荡水翼式潮流能发电装置等。其中振荡水翼式潮流能发电装置是一种对海洋环境适应性强以及破坏小的新型海洋能发电装置。水翼的结构参数直接关系到震荡水翼的运动状态，以及系统能否运转和能量转化效率的大小，是本领域研究的重点和难点。目前国内外的学者都以二维水翼作为研究对象，并建立其运动的水动力学方程。其中ZHU等 ^{[4, 5, 6]}首先提出了二维振荡水翼模型及其水动力学方程，得出转动中心的位置和回转运动弹簧刚度系数是影响系统运动状态主要参数的结论。随后邓见等^[7]基于该方程，在考虑了质量和转动惯量的影响下，采用数值模拟方法得出了同样的结果。 根据PENG和ZHU^[4]的模型,振荡水翼吸收的能量主要来自垂直来流方向的升力作用，可是文献^[8]的结论表明，平直翼水平方向的诱导阻力比梯形翼要大。而诱导阻力变大会降低升力的作用，影响振荡水翼的能量转化效率。所以本文根据叶素理论，以梯形翼作为三维水翼模型的代表，在建立了梯形振荡水翼的水动力学模型和运动方程之后，着重研究了梯形比(表征梯形翼的主要参数)对振荡水翼的运动状态和能量转化效率的影响。

根据PENG和ZHU^[4]等的研究成果建立的两自由度的NACA0012翼型振荡水翼模型，其几何模型如图1所示。图中U为来流速度，a为水翼弦长，k_α为转动中心的扭转弹簧刚度，b为转动中心与水翼中心的距离，c为升沉方向作为能量输出的阻尼器阻尼系数。考虑质量及转动惯量后，可得动力学方程如下：

式中：h为升沉方向的位移，α为俯仰角，L为水翼所受的水动力学升力，M为水翼所受的俯仰力矩，m、I分别为系统质量和转动惯量。

图1(Fig.1)

图1 两自由度的振荡水翼二维模型 Fig.1 2D oscillating hydrofoil model of two degrees of freedom

根据PENG和ZHU^[4]的结论，二维振荡水翼在k_α-b平面内被分成4种运动状态：(1)水翼静止不动；(2)水翼呈现周期性俯仰振荡与升沉振荡的组合运动；(3)水翼沿着两个平衡点作周期性俯仰振动，而升沉运动则无规律；(4)水翼作周期性俯仰振动，而升沉运动则向一边偏移。

要建立三维梯形振荡水翼模型，首先要引用叶素理论。叶素理论基本思想是：将叶片沿着半径方向分成若干微段，称这些微段为叶素；假设各段不互相影响，即叶素可看成二维翼型，这样就可以对每个叶素进行分析，将作用在每个叶素上的作用力叠加，就可以得到作用在水轮上的力和力矩等参数^[9]。叶素理论应用基于以下假设：

(1)不考虑相邻叶素之间的相互干扰；

(2)作用在每个叶素上的力由翼型的动力特性决定。

因此结合梯形翼结构特点，选择NACA0012作为水翼剖面翼型，其余参数选取b=0，k_α=0.2，建立了如图2所示的三维振荡水翼模型，图(a)为水翼参数模型，图(b)为水翼的俯仰-升沉运动模型。

图2(Fig.2)

图2 两自由度的梯形振荡水翼三维模型 ar为水翼翼根处的弦长，at为翼尖处的弦长，二者之比称为梯形比(k=at/ar),kα为水翼转动中心的扭转弹簧刚度，c为升沉方向的阻尼器阻尼系数。 Fig.2 3D trapezoid oscillating hydrofoil model of two degree of freedom

根据叶素理论，水翼整体的受力可以表示为单个叶素截面上的受力之和，可建立三维水翼的动力学方程如下所示：

式中：L_j、M_j分别代表每个截面所受的二维水翼升力和转矩。

在振荡水翼系统中，根据文献^[5]的结论，其输出功率为P=ch.²，周期性振荡下的平均输出功率为：

式中：T为振荡周期。经过水翼的来流动能为：

式中：Y_p为水翼扫过的面积(取前缘或后缘扫过的面积较大者)，d为翼展的宽度，因此得到系统能量转化效率为：

采用有限体积法求解三维非定常黏性不可压缩Navier-Stokes方程，雷诺数为Re=1×10⁵，流动为湍流，湍流模型选用RNG k-ε模型。因为梯形水翼的三维结构复杂，采用非结构网格，所以对于水翼的流固耦合问题适合采用动网格技术；另外水翼运动属于大变形问题，网格重构采用弹簧光顺法和网格局部重构法耦合^[10]。则其控制方程为：

俯仰运动需要用到随体运动的旋转坐标系。现假设(x,y,z)为惯性坐标系，(x₀,y₀,z₀)为旋转坐标系，考虑到俯仰运动只以角速度ω绕z轴旋转，其余绕轴运动不作考虑，得到了如图3所示的旋转坐标系与惯性坐标系之间的坐标变换关系图。

图3(Fig.3)

图3 坐标变换关系图 Fig.3 The relation of transformation of coordinates

设水翼表面上的一点在惯性坐标系中的坐标为P=(x,y,z)，在旋转坐标系上的坐标为(P₀=(x₀,y₀,z₀)，振荡水翼绕z轴转动角度为θ(t)，则P点在旋转坐标系中的速度为：

式中：v_r为相对速度(旋转系中观察的速度)；v为绝对速度(惯性系中观察的速度)；u_r为牵连速度(由于旋转系所导致的速度)。通过式(10)可以实现振荡水翼速度在旋转坐标系和惯性坐标系之间的转化，以便于数值仿真中对参数的计算。

求解程序的基本算法：空间离散上，对流项采用二阶迎风格式，拉普拉斯项采用中心差分格式，时间上采用一阶隐式，不可压缩流动模拟中压力与速度的耦合采用PISO算法。动力学方程与流体方程采用时间上的松耦合^[7]。

计算区域与网格如图4所示。进口边界距离水翼前缘15倍弦长，出口边界距离后缘15倍弦长，上下、左右边界距离水翼同为10倍弦长。边界条件设置：入口设置为速度入口，出口为充分发展的流动，水翼表面为无滑移边界条件，其余为壁面条件，将进口速度作为初始化条件。

图4(Fig.4)

图4 计算区域与网格 Fig.4 Calculation region and mesh grid

常微分方程采用欧拉法求解。在均匀来流U=1 m/s，雷诺数条件下Re=1×10⁵，采用RNG k-ε模型，对基于动网格法的梯形振荡水翼三维模型进行了数值模拟，得到以下的结果。

水翼的位移和转动角度随着时间的变化关系如图5所示，水翼的运动状态呈现周期性俯仰振荡与升沉振荡的组合运动，此时水翼能够稳定地获得能量输出。与文献^{[4, 7]}对二维振荡水翼模型的研究结果相比，运动状态相同。这一方面表面k_α-b相平面内，振荡水翼无论是二维还是三维的模型，运动状态不发生变化；另一方面表明采用动网格技术对由叶素理论建立的梯形振荡水翼水动力学模型进行数值模拟是可靠的。

图5(Fig.5)

图5 水翼的位移和转动角度随着时间的变化曲线(κ=0.1) Fig.5 The curve of hydrofoil’s displacement and rotation angle changing with time(κ=0.1)

不同梯形比的振荡水翼的位移随时间变化曲线，以及转动角度随时间变化曲线分别如图6所示。从图中可以看出，随着梯形比的升高水翼的振动周期明显变长，但是俯仰运动的角度幅值只有细微的增加；升沉运动的振动周期与俯仰运动的变化相同，但是其位移幅值随着梯形比的降低而降低。这表明梯形比不会影响水翼的工作状态，但是会改变水翼的振动周期和升沉运动的幅值，相关文献^{[4, 7]}得出的k_α-b是影响水翼工作状态的决定性因素的结论不变。

图6(Fig.6)

图6 水翼的位移和转动角度分别随着时间的变化曲线 Fig.6 The curve of hydrofoil’s displacement and rotation angle changing with time(κ=0.1)

根据图6(a)，可以假设水翼振荡为余弦振动，其运动方程表达式为：

式中：h_a表示位移的幅值，φ₀为相位角。则根据式(6)、(7)，可得到水翼作周期性运动的能量转化效率为：

由式(12)计算得到各个梯形比下的水翼能量转化效率，如图7所示。从图中可以看出，水翼的转化效率随着梯形比的降低先升高后降低。根据《空气动力学》^[8]的结论，当梯形翼的梯形比接近0.5时，梯形翼的升力分布与椭圆翼越相似，而椭圆翼的诱导阻力最低，即梯形比为0.5时梯形翼的诱导阻力最低；而随着梯形比过大或者过小，都会导致诱导阻力上升。由此可以得出诱导阻力和振荡水翼的能量转化效率随着梯形比的变化趋势是相同的。需要指出的是当梯形比接近0时，在翼尖附近的剖面升力系数会急剧地加大，翼尖处水翼上下表面的环流相交，环量为0，即升力系数为0，致使此处的升力发生突变。另外考虑到振荡水翼的周期性振荡，在往复力的作用下，水翼翼尖处极易发生疲劳破坏。因此小梯形比的梯形翼是不适合用在振荡水翼上的。综上，选取梯形比为0.5的梯形翼较为合适，这也是符合空气动力学的结论的。

图7(Fig.7)

图7 水翼能量转化效率 Fig.7 The energy conversion efficiency of hydrofoil

基于PENG和ZHU^[4]的两自由度二维振荡水翼式海流能采集系统，本文提出了三维梯形振荡水翼模型，并对其进行了数值模拟研究。通过结果对比分析，得出三维水翼的运动状态与二维水翼的相同，符合叶素理论的假设，所以本文采用的模型与数值仿真方法可靠。在此基础上本文分别对梯形比分别为0.25、0.50、0.75、1.00的梯形水翼进行了模拟仿真，得到了以下的结论：1)k_α-b是影响水翼的工作状态的决定性因素，梯形比只能够影响系统的振动频率，而不会改变系统的稳定振动状态；2)采用梯形翼能够提高振荡水翼的能量转化效率，这是因为随着梯形比接近0.5，水翼的诱导阻力变小；同时当梯形比进一步减小时，诱导阻力反而会升高，系统的能量转化效率因此开始快速降低，所以采用梯形比为0.5的梯形翼是相对较优的选择(在参数k_α=0.2，b=0条件下)。

振荡水翼在k_α-b相面内处于稳定振荡状态的区域较大，本文只选取了其中较为典型的一组参数进行了数值模拟分析，其能量转化效率也不是最高的组合。因此梯形翼的几何参数设计还有待进一步优化。但是采用梯形翼设计可以作为一个着重研究的方向。