某型波浪发电装置3D打印样机飞逸转速的风洞测试与数值仿真

郑贵林，梁灿棉

(武汉大学动力与机械学院，武汉 430072)

0 引 言

凭借存在范围广、储量大、能量形式为机械能等特点，波浪能被视作可再生自然能源开发利用的重要一员，各种形式的波浪能采集装置也不断涌现。而海洋具有表面波浪起伏，海面20～30 m以下平稳的特点[1,2]，于是差动垂直波浪能提取发电装置应运而生。

该发电装置的特点是利用导流叶片使内部旋转叶轮在海水相对其从上向下流还是从下向上流都保持一个旋转方向[3,4]。

目前，对于因受到外界流体的作用而被动转动的旋转叶轮例如风力发电机和水轮发电机组等在某一特定负载力矩下的旋转速度的预测的研究尚少，一般都是通过具体的样机试验测定，而叶轮的旋转速度又是决定变速机构的设计的一个重要参数，如果能够使用CFD软件准确地模拟计算出叶轮在不同流体速度和不同负载力矩情况下的旋转速度，将对变速机构的设计带来极大的方便，并且方便对叶轮的效率分析以及进一步的优化设计工作。近年来，国内外有部分学者使用CFD软件来进行水轮机飞逸状态(负载力矩为零的状态)的三维数值模拟并取得可观的成效。周大庆等[5]对轴流转桨式水轮机在不同桨叶启闭规律下的飞逸状态进行计算，得出飞逸转速在不同桨叶启闭规律下的变化情况以及尾水管负压值、压力脉动等的变化。李金伟等[6]对某型混流式水轮机在不同开度下的飞逸状态进行三维非定常计算，得出飞逸状态下单位转速、单位流量以及水力矩的变化情况。罗兴锜等[7]对贯流式水轮机的飞逸状态进行数值计算，获得飞逸转速、流量、力矩、轴向力等外特性参数以及水轮机内部流场的动态特性。

如今对于利用海洋能量来进行发电的被动型旋转叶轮的飞逸状态的三维数值模拟的研究并不多见，差动垂直波浪能提取发电装置其中的能量获取机构便是一种旋转叶轮，对该旋转叶轮的转速预测可以为变速机构的设计、效率分析以及优化设计提供参考指导作用。由于对大尺寸的发电装置进行海试实验的成本高、操作难度大、改进的研发周期长等，而Fluent等一些计算流体力学软件进行三维数值仿真具有成本低、试验方便、改进研发周期短的特点，因此可以使用Fluent为差动垂直波浪能提取发电装置的优化改进提供一定的辅助作用。另外随着3D打印技术日渐成熟，能够打印出具有较为复杂表面的实体，因而便借助3D打印来制作小尺寸的发电装置样机来进行实物实验来验证装置模型设计的正确性和数值计算方法的可行性等[8]。而要进行实物实验除了需要发电装置样机还需要实验环境，由于发电装置的工作过程与风力发电机的工作过程比较类似，即是利用流体与装置之间的相对运动的能量，另外波浪的垂直方向速度曲线近似梯形波时存在速度比较平稳的阶段，因而决定使用相对较简单的风洞实验来取替水槽实验来模拟起伏运动中速度比较平稳的阶段并近似为单向恒速状态。为此，本文通过Fluent软件建立差动垂直波浪能提取发电装置中旋转叶轮在单向恒定流速空气流作用下的飞逸转速的三维湍流数值计算方法，并与样机的风洞测试实验进行对比，以验证数值计算方法中被动型动网格、湍流模型等参数设置的正确性，为后续双向变流速的水流作用下以及带负载作用下的仿真和优化提供参考指导。

1 差动垂直波浪能提取发电装置的结构、工作原理以及3D打印的样机

1.1 结构组成及工作原理

差动垂直波浪能提取发电装置的总体结构如图1所示，其中包括海面漂浮平台、漂浮平台上安置的发电机、传动转轴和能量转换波轮机，而波轮机由导流圆管、导流叶片和转子等组成。

差动垂直波浪能提取发电装置的运行原理为：如图1所示，漂浮平台位于海洋表面波浪区，波轮机置于稳定区。漂浮平台随波浪上下运动并通过传动转轴带动波轮机在垂直方向上下运动，因此波轮机与其周围稳定区中的水流存在垂直方向上的相对往复运动，把参考系放在波轮机上，便可以看作稳定区中的水流在垂直方向上下往复地流经波轮机。水流在上下往复的过程中，经由上部和下部导流叶片的引导始终单方向作用于转子叶片从而使转子单方向地旋转，再经由传动转轴带动发电机旋转发电。

图1 差动垂直波浪能提取发电装置示意图Fig.1 The differential vertical wave energy converter

1.2 3D打印的样机

由于使用风洞实验取替水槽实验，直接有空气流体与波轮机作相对运动，不需要漂浮平台来带动波轮机做上下起伏的运动，故只使用3D打印机打印差动垂直波浪能提取发电装置中的波轮机部分。打印的波轮机主要参数如下：导流圆管最大直径：D1=100 mm；导流圆管最小直径：D2=80 mm；转子轮毂直径：D=20 mm；导流叶片数量：Z1=12或16；转子叶片数量：Z2=12或16；转子叶片截面形状：圆弧形和角形。打印完成的转子零件以及进行组装所得的波轮机样机如图2，图3所示。

图2 圆弧形叶片转子以及角形叶片转子Fig.2 Rotor with circular arc blades and rotor with folded angle blades

图3 完整组装的波轮机样机Fig.3 The whole propotype

2 数值计算方法与模型

2.1 运动方程及转速计算方程

在波轮机转子飞逸过渡过程中，其外负载为0，转子被水流冲击而转速迅速上升，其运动方程为：

(1)

式中：Mt为主动力矩，N·m；Mg为负载力矩，N·m；J为机组总转动惯量，kg·m2；ω为转子旋转角速度，rad/s。

每一步的转速ωi+1可以通过下式来求得：

(2)

式中：Δt是计算时间步长，s。

2.2 湍流模型与边界条件

本文采用商业CFD软件Fluent16.0来模拟波轮机在风洞测试仪中的飞逸过程，计算中采用非定常(即瞬态)的标准k-ε湍流模型。流体介质选择为空气，密度为1.225 kg/m3，动力黏度为1.789 4×10-5kg/s，进口采用速度边界条件，流体速度设置为5 m/s，方向垂直于进口截面。出口采用自由出流边界条件。由于转子转速是根据式(2)来确定，属于被动转动，因此采用6DOF被动型动网格模型来设置转子所在的流体域，其他部分的流体域则为静流体域，不同的流体域之间的重合面设置为interface交界面。通过编写用户自定义函数(UDF)来输入机组的总转动惯量、质量、运动约束条件等参数。计算过程中以静止状态为初始值，时间步长设为0.001 s。

由于本文的目的是通过对波轮机样机进行单向恒速空气流风洞测试与仿真计算来对比以验证被动型动网格、湍流模型等参数设置的正确性，因此对于实际环境中的双向变流速水流的边界条件未作考虑。对于垂直方向速度曲线近似为梯形波的波浪，由于存在速度比较平稳的阶段，此阶段内的力矩特性可以用单向恒速流体作用下的力矩特性来近似，假若速度曲线梯形波的斜线段时间很短即速度曲线接近矩形波的时候，整个周期内的力矩特性都能够用单向恒速流体作用下的力矩特性来近似。工作介质为空气时，计算结果所得的力矩值小于0.001 N·m，现有的测量设备难以测量如此小的力矩值，因此风洞测试的力矩值没有进行测量。计算工作介质改为水时的飞逸转速有约10%的提高，力矩值提升约三个数量级，与水的密度和空气的密度比值相接近。

2.3 流体场网格

按照风洞测试仪的尺寸来建立直径为25 cm，长度为53 cm的圆柱体风洞流体场，波轮机轴线与该圆柱体轴线重合，波轮机前端距离风洞流体场入口截面16 cm，后端距离风洞流体场出口截面30 cm，对于圆柱体风洞流体场可以通过切割分块并采用六面体网格划分，如下图4所示。由于转子叶片和导流叶片的几何外形比较复杂，故这些叶片附近的流体域采用非结构化四面体网格来划分，而距离叶片较远的区域则切割分块并采用六面体网格划分以减少网格数量，如图5所示。

图4 风洞流体场网格划分示意图Fig.4 Flow field grid division

图5 波轮机内部流体域网格划分示意图Fig.5 Internal flow field grid division of the converter

2.4 网格无关性验证

为减小网格数量对仿真计算结果的影响，分别选取网格数为442万(具体是4 423 871)、513万(具体是5 134 559)、550万(具体是5 503 763)、606万(具体是6 061 361) 这4种网格模型对进行了计算，以波轮机转子飞逸转速为目标变量来进行网格无关性验证，验证结果如图4所示。由图6可知，网格数量达到513万以上波轮机转子飞逸转速几乎无变化，考虑计算精度和计算经济性，确定模型网格数量为513万。

图6 网格无关性测试结果Fig.6 Gird independence test results

3 数值计算结果与风洞测试结果

3.1 飞逸转速对比

对定子叶片数量与转子叶片数量均为12和定子叶片数量与转子叶片数量均为16，转子叶片形状分别为圆弧形和折角形这四种样机在风速5 m/s下的飞逸转速的数值计算结果与风洞测试结果如表1所示。从表1中可以看出数值计算所得的飞逸转速与风洞测试所得的飞逸转速相对误差在10%以内，说明数值计算飞逸转速的可行性，以及圆弧形叶片的转子的飞逸转速比折角形叶片的转子的飞逸转速更高。

表1 两种形状叶片转子的数值计算飞逸转速与风洞测试飞逸转速的对比Tab.1 Comparison of the runway speeds of the numerical simulations and the wind tunnel tests(rotors with circular arc blades and folded angle blades)

对于数值计算所得的飞逸转速与风洞测试所得的飞逸转速的误差分析，作者认为有主要两个部分：第一部分，由于用于风洞测试的样机内部安装有轴承，虽然是测试飞逸状态也就是不带外负载力矩，但是轴承的滚动摩擦力矩依然会对样机的飞逸转速产生轻微的影响，而在数值计算中没有考虑对轴承滚动摩擦力矩的设置；第二部分，3D打印所得的样机模型与在三维绘图软件中的原始设计模型仍然存在一些外形和尺寸误差，而在数值计算中使用的是三维绘图软件中的原始设计模型。

3.2 转动惯量与飞逸转速关系分析

对定子叶片数量与转子叶片数量均为12和定子叶片数量与转子叶片数量均为16，转子叶片形状为圆弧形，但转子转动惯量不一样的4种样机在风速5 m/s下的飞逸转速的数值计算结果与风洞测试结果如表2所示。

从表2中可以看出，转子转动惯量的改变对飞逸转速的影响几乎没有(在不考虑因为转动惯量变化带来机械摩擦变化的情况下)。经过仿真得到该装置的力矩-转速曲线如下图7所示。结合力矩-转速曲线以及式(2)显而易见，当Mt=Mg=0以后，ωi+1=ωi，即转子最终稳定的转速总是达到力矩为0时对应的转速值也就是飞逸转速值，而转动惯量主要影响了转子从初始状态到飞逸状态之间变化的时间，该时间长度T=iΔt，i为式(2)达到飞逸状态迭代计算的所使用的步数。

表2 不同转动惯量转子的数值计算飞逸转速与风洞测试飞逸转速的对比Tab.2 Comparison of the runway speeds of the numerical simulations and the wind tunnel tests(rotors with different moments of inertia but the same circular arc blades)

图7 仿真计算的力矩-转速曲线Fig.7 Curves of torque against rotate speed

4 结 语

通过对差动垂直波浪能提取发电装置中的能量转换机构波轮机样机的飞逸转速进行单向恒速风洞测试以及使用Fluent软件进行三维数值仿真，得出以下结论。

(1)样机在风洞测试中的飞逸转速值与Fluent仿真计算所得的飞逸转速值接近，误差在10%以内，说明数值计算的可行性。

(2)在其他条件相同的情况下，圆弧形叶片的转子与折角形叶片的转子相比具有更高的飞逸转速。

(3)保持其他条件不变的情况下，转子转动惯量的变化对转子的飞逸转速的影响几乎没有。