水沙环境下海流能机组翼型水动性能

2023-09-01高艳婧倪艺铭刘宏伟林勇刚

中南大学学报（自然科学版） 2023年7期

高艳婧，倪艺铭，刘宏伟，林勇刚

(浙江大学流体动力与机电系统国家重点实验室，浙江杭州，310027)

我国近海岸海流因水土流失导致含沙量较高[1-3]。泥沙中含有大量的固体颗粒和其他杂质，会对水力机械水动力性能造成严重影响，加速水力机械的磨损，严重影响其工作效率和使用寿命[4]。浙江省杭州湾位于河流入海口处，属于近海沿岸海域，水域较浅，底部泥沙在水流的曳力作用下悬浮移动，导致海流能机组在含沙的工况下运行[5]。颗粒与流体之间的相互作用是否会改变翼型截面的流场以及影响叶轮的捕获功率尚且未知，因此，研究两相流环境下的叶轮捕获功率可以为后续海流能机组选址以及叶片设计提供参考。

国内外大部分的研究均针对海流能机组叶片表面因沙粒冲击、空化以及海水腐蚀引起的粗糙度增加而造成的功率损失[6-8]展开，对于水沙环境直接对海流能机组捕获功率的影响的研究甚少。ORME等[9]利用风洞研究了藤壶对机翼升力系数和阻力系数的影响，发现机翼上覆盖不同尺寸和密度的理想藤壶时，翼型升阻比降低。BATTEN等[10]基于叶素动量理论(BEM)研究了污垢对叶片粗糙度增加的潜在影响，发现在叶片被污染的情况下阻力系数增加，而升力系数没有改变，基于叶素动量理论预测在高的叶尖速下功率系数下降6%～8%。WALKER 等[11]通过风洞实验测试了NACA63618 翼型在清洁和粗糙条件下的升阻力系数，而且通过水池实验测试了模型叶轮在不同表面条件下的捕获性能，发现升力最大下降幅度的平均值为11%，阻力在重度粗糙下能上升150%，捕获功率最大降低19%。而对于沙粒与流体的相互作用对海流能叶片的捕获功率的影响尚不明确。不同海域的颗粒直径、粒径分布以及颗粒浓度均有差异[12]，研究这些因素对于叶轮捕获功率的影响，可为海流能机组选址和叶片设计提供参考。

在本研究中，以120 kW 海流能叶片翼型NACA63-为研究对象，采用CFD-DPM 模型(考虑颗粒与流体相互作用)模拟不同颗粒属性下翼型的升力系数和阻力系数。根据沿同一颗粒编号轨迹运动的颗粒速度和有/无颗粒的流体速度，分析颗粒对流体速度造成的影响，再依据颗粒偏离流线的程度揭示在流体速度变化的情况下的升阻力系数的变化。

1 理论基础

翼型是海流能叶片的重要组成部分，翼型参数直接影响海流能叶片的性能。NACA63-翼型是海流能叶片最常用的翼型。水流经过翼型前缘附近的压力最高点即驻点后，分成上、下两股，分别沿翼型上、下表面流过，在后缘重新汇合，翼型上下表面分别形成变化的压力。根据伯努利方程，当水流流经翼型凸面时，平均流速增大、压力减小，此面被称为吸力面。当水流流经翼型凹面时，平均流速减小、压力增大，此面被称为压力面。叶片功率与翼型升力系数呈正相关，与翼型阻力系数呈负相关。

翼型升力源于吸/压力面压差，吸力面侧流体速度先增大后减小，压力面侧流体速度先减小再增大。同一流线位置的流体速度变化大小可以反映升力系数大小，即同一流线位置流体速度变化越大，升力系数越大。在相同的雷诺数和攻角下，流线距离翼型表面越远，沿流线流体速度变化越小，如图1所示。当升力系数变化值相同时，流线距离翼型吸/压力面越远，沿同一流线流体速度变化值越小，如图2所示。由图2可以推断，沿同一流线流体速度变化值相同时，流线距离翼型吸/压力面越远，升力系数变化值越大。翼型阻力系数也与雷诺数和攻角相关，随着流体速度增大，翼型阻力系数减小。

图1 不同流线位置的流体速度Fig. 1 Fluid velocities at different streamline positions

图2 升力系数变化值相同(攻角从6°到8°)时不同流线位置的流体速度变化值Fig. 2 Change values of fluid velocity at different streamline positions when the lift coefficient changes by the same value(angle of attack from 6° to 8°)

考虑到随着颗粒直径或颗粒形状因子的变化，颗粒偏离流线的程度不同，故在研究不同直径、不同浓度和不同颗粒形状因子下的颗粒速度和无颗粒流体速度时，以同一颗粒编号的颗粒轨迹为研究对象。由远离吸力侧指向远离压力侧方向，颗粒编号逐渐增大，颗粒编号为30000 和编号为36000的颗粒轨迹位置如图3所示。

图3 颗粒编号为30000和编号为36000的颗粒轨迹位置Fig. 3 Particle track positions when particle IDs are 30000 and 36000

2 计算模型

2.1 计算域和边界条件

以叶尖NACA63-翼型为研究对象，基于叶素理论[13-14]，将120 kW 海流能叶片分成为若干翼型叶素。在ICEM软件中对翼型和流域进行几何建模和网格划分。图4所示为计算域和边界条件，以翼型尾部为圆心，上游与翼型的距离为弦长c的12.5倍，以减小边界对翼型附近流体流动的影响；下游与翼型的距离为16 倍弦长，以保证流体能够充分流动；展向取1倍弦长。颗粒喷射面与翼型的距离为弦长的6倍，颗粒喷射面是在FLUENT中创建的虚拟面，通过设置节点数来保证颗粒喷射面网格均匀[15]，颗粒喷射面长度为10 倍弦长，保证升阻力系数不再随颗粒喷射面长度的增加而变化。

进口处采用速度入口边界，出口处采用压力出口边界，将翼型定义为无滑移壁面，颗粒在颗粒喷射面以与流体相同的速度进入流场。颗粒在翼型壁面时设置为“reflect”，在出口处设置为“escaped”。湍流模型采用SSTk-ω模型，边界层湍流和自由剪切流均具有较理想的精度与稳定性。速度与压力的耦合选择SIMPLE算法，动量、湍流动能和耗散率均采用二阶迎风格式，并采用基于压力耦合的隐式求解器[16]。

图4 计算域和边界条件Fig. 4 Calculation domain and boundary conditions

2.2 网格划分

流场网格划分采用结构化网格方法，并用C-Block 方法对翼型附近的网格进行加密[17]。基于SSTk-ω模型，首层网格量纲一高度y+要小于1[18]，首层网格高度为0.02 mm。

为保证求解的准确性，进行网格无关性验证[19]。随着网格数量的增加(从网格1 到网格6)，翼型下游8倍弦长位置处的速度以及湍流动能基本上不发生改变(图5)。考虑到计算精度和时间成本，选择网格3，翼型表面分布340 个网格节点，翼型壁面法向分布180个网格节点，网格增长率为1.2，网格单元总数为32 万个左右。图6 所示为结构化网格划分后的流体域网格和翼型表面网格。

图5 不同网格分辨率翼型下游8倍弦长位置处的速度和湍流动能变化曲线Fig. 5 Velocity and turbulent kinetic energy curves at 8 chord lengths downstream of airfoils with different grid resolutions

图6 计算域网格Fig. 6 Computational meshes used in the simulation

2.3 模拟方法验证

2.3.1 SSTk-ω模型验证

图7所示为无颗粒海水条件下的翼型升阻力系数仿真结果与通过风洞实验得到的升阻力系数[20](相同雷诺数)的对比，可见升阻力系数仿真结果与实验结果[20]具有较高的吻合性。

图7 升阻力系数仿真和实验结果对比Fig. 7 Comparison of the simulated and experimental results of lift and drag coefficients

2.3.2 DPM模型验证

将受沙粒冲击的120 kW 海流能机组叶片的磨损数据作为实验数据，120 kW 海流能机组叶片参数如表1 所示，包含30 个翼型叶素。选择考虑颗粒形状因子和冲击角的Ahlert侵蚀模型[21]，仿真模型中颗粒属性和叶片壁面属性如表2所示。详细求解过程如下：首先，通过求解雷诺平均Navier-Stokes方程得到流体流动的解；然后，考虑流体对颗粒的阻力获得颗粒碰撞信息；最后，依据颗粒碰撞信息通过Ahlert侵蚀模型预测磨损特性[22]。

表1 120 kW海流能机组叶片参数Table 1 Parameters of the 120 kW tidal current turbine blade

表2 材料属性Table 2 Material properties

图8(a)所示为实际海况下120 kW 海流能叶片的沙粒冲蚀图，图8(b)所示为考虑粒径分布(表3)的120 kW 海流能叶片仿真冲蚀云图。图9 所示为最大侵蚀位置(以翼型前缘点为原点，沿弦长方向沙粒对壁面造成的最大侵蚀距离)的仿真和实验结果对比，可见最大侵蚀位置从叶根到叶尖逐渐向前缘移动，与实验结果较吻合。

图8 120 kW海流能叶片的颗粒冲击冲蚀图Fig. 8 Erosion diagrams of the 120 kW tidal current turbine blade

表3 舟山海况颗粒粒径分布Table 3 Particle size distribution in Zhoushan sea condition

图9 最大侵蚀位置仿真和实验结果对比Fig. 9 Comparison of the maximum erosion location of the simulation and the experiment results

本研究利用CFD-DPM 模型计算120 kW 海流能机组叶片各翼型的升力系数和阻力系数，再根据叶素动量理论计算叶轮的输出转矩和功率。不同流速下叶轮功率的仿真结果与实验结果的对比见图10。由图10 可以看出，模拟的叶轮功率与实验结果较吻合，验证了DPM模型的准确性。

3 结果及分析

仿真算例的运行条件如表4所示，研究不同颗粒属性下翼型升阻力系数。在本算例中，颗粒直径dP为500 μm，质量浓度CP为1.5 g/L，形状因子f为1，弦长L为1 m，攻角α为6°，相对流体速度U为16 m/s。

表4 仿真算例运行条件Table 4 Operating conditions of the simulation cases

3.1 颗粒直径

3.1.1 颗粒直径对翼型升力系数的影响

图11所示为翼型升力系数随颗粒直径的变化。由图11可知：当颗粒直径小于100 μm时，翼型升力系数大于无颗粒时的翼型升力系数；当颗粒直径大于100 μm 时，翼型升力系数小于无颗粒时的翼型升力系数；随着颗粒直径增加，翼型升力系数在先减小后增大的趋势中波动。

图11 翼型升力系数随颗粒直径的变化Fig. 11 Change of airfoil lift coefficient with particle diameter

图12所示为无颗粒时的流线，图13所示为不同颗粒直径下的颗粒轨迹。由图12 和图13 可知：当颗粒直径为20 μm 时，颗粒轨迹与流线基本重合；随着颗粒直径的增加，颗粒轨迹偏离流线的程度增大。

图12 无颗粒流线Fig. 12 The streamlines without particles

图13 不同颗粒直径下的颗粒轨迹Fig. 13 Particle trajectories with different particle diameters

图14 所示为当颗粒直径为20 μm 时，沿同一颗粒编号轨迹运动的颗粒和流体速度曲线。由于颗粒直径较小、惯性较小，颗粒轨迹主要受流体阻力控制[23-25]。由图14 可知：在颗粒喷射面颗粒速度急剧下降，随流体速度变化而变化，等效于流体密度增加，从而导致翼型升力增加。

图14 沿5号颗粒轨迹运动的颗粒速度和流体速度(颗粒直径为20 μm)Fig. 14 Velocities of particle and fluid along the particle trajectory of the No.5 particle(particle diameter is 20 μm)

图15所示为当颗粒直径大于100 μm时，沿相同颗粒编号轨迹运动的颗粒速度和无颗粒流体速度。由图15 可知：受翼型影响，压力面流体速度下降，流体对颗粒的阻力同样导致颗粒速度下降。但流体速度的减小程度大于颗粒速度的减小程度，即颗粒速度大于流体速度。考虑颗粒对流体的影响，压力面流体速度减小的程度降低，导致升力系数减小。

图15 不同颗粒直径下沿5号颗粒轨迹运动的颗粒速度和无颗粒流体的速度Fig. 15 Particle velocity and the velocity of particle-free fluid along the trajectory of No.5 particle with different particle diameters

从图15 还可以看出：随着颗粒直径增加，Up-fn增加(Up-fn为沿相同颗粒编号轨迹运动的颗粒速度与无颗粒流体的速度之差)，单个颗粒对流体的影响增大，Ufp-fn增加(Ufp-fn为沿相同颗粒编号轨迹运动的下流体速度与无颗粒流体的速度之差)；随着颗粒直径继续增大，颗粒质量浓度不变，与流体作用的颗粒数减少，导致Ufp-fn减小；随着颗粒直径增大，Ufp-fn先增大后减小，沿同一颗粒轨迹的流体速度变化值先增大后减小。

随着颗粒直径的增大，颗粒的惯性增大，具有相同颗粒编号的颗粒远离压力面，如图13所示。从图13 可见：当沿同一颗粒编号轨迹运动的流体速度变化值相同时，大直径颗粒会使升力系数(相对于无颗粒时的升力系数)减小更多，即升力系数更小；由于同一编号颗粒的位置不同导致流体速度变化值与升力系数的变化值不具有统一性，因此，随着直径增大，升力系数在先减小后增大的趋势中波动；当颗粒直径为2 500 μm 时，升力系数变化最大，和无颗粒相比升力系数减小了0.141%。

3.1.2 颗粒直径对翼型阻力系数的影响

图16所示为翼型阻力系数随颗粒直径的变化。由图16 可知：随着颗粒直径增大，翼型阻力系数先增大后减小。

图16 翼型阻力系数随颗粒直径的变化Fig. 16 Change of airfoil drag coefficient with particle diameter

图17 所示为当颗粒直径为1 500 μm 时，沿23000 号和45000 号颗粒轨迹运动的颗粒流体速度和无颗粒流体速度，其中，路径长度以颗粒喷射面作为起点。由图17 可知：翼型阻力系数随速度的减小而增大，与无颗粒流体速度相比，等效流体速度的减小程度越大，翼型阻力系数越大；随着粒径的增大，Ufp-fn先增大后减小，等效流体速度减小程度先增大后减小，阻力系数先增大后减小；当颗粒直径为1 500 μm时阻力系数变化最大，与无颗粒流体相比阻力系数增加了1.872%。

图17 沿45000号和23000号颗粒轨迹运动的无颗粒流体速度与颗粒作用下的流体速度Fig. 17 Velocities of fluid with or without particles influence along the trajectories of No.23000 and No.45000 particles

3.2 颗粒质量浓度的影响

3.2.1 颗粒浓度对翼型升力系数的影响

翼型升力系数随颗粒质量浓度的变化如图18所示。由图18 可知：翼型升力系数与颗粒质量浓度呈线性关系。图19 所示为不同质量浓度下沿5号颗粒轨迹运动的颗粒速度和无颗粒流体的速度。由图19 可知：随着颗粒质量浓度的增加，颗粒与无颗粒流体速度差保持不变，即单个颗粒对流体速度变化的影响不变。图20 所示为不同颗粒质量浓度下颗粒轨迹。从图20 可以看出：随着颗粒质量浓度增加，颗粒偏离流线的程度不变，单个颗粒对翼型升力系数变化的影响保持不变；与流场相互作用的颗粒数与颗粒质量浓度成正比，即对于相同的粒径，当颗粒质量浓度增加值相同时，与流体相互作用的颗粒数增加值相同。因此，翼型升力系数与颗粒质量浓度呈线性关系。

图18 翼型升力系数随颗粒质量浓度的变化Fig. 18 Change of airfoil lift coefficient versus particle mass concentration

当颗粒质量浓度为2 g/L，粒径为2 500 μm时，升力系数降低了0.188%。考虑到舟山海况中的最大颗粒质量浓度为5.522 g/L，当粒径为2 500 μm时，升力系数降低了0.478%。

3.2.2 颗粒质量浓度对翼型阻力系数的影响

图21 所示为翼型阻力系数随颗粒质量浓度的变化，由图21 可知：随着颗粒质量浓度增加，翼型阻力系数呈线性增大；当粒径为1 500 μm，颗粒质量浓度分别为2.000 g/L 和5.522 g/L 时，与无颗粒流体相比阻力系数分别增加了2.418% 和6.676%。

图21 翼型阻力系数随颗粒质量浓度的变化Fig. 21 Change of airfoil drag coefficient with particle mass concentration

3.3 颗粒形状因子的影响

3.3.1 颗粒形状因子对翼型升力系数的影响

颗粒形状是造成海流能机组磨损严重的重要原因之一，下面对颗粒形状因子对翼型升阻系数的影响进行研究。颗粒的形状因子f指与颗粒体积相同的球体表面积与颗粒实际表面积的比值，其定义如下：

式中：a为与非球形颗粒体积相同的球形颗粒表面积；As为非球形颗粒的表面积。

图22 所示为不同颗粒形状因子下翼型升力系数随粒径的变化。由图22 可知：随着颗粒形状因子的增大，使颗粒翼型升力系数大于无颗粒时流体翼型升力系数的临界颗粒直径减小。小形状因子的颗粒具有较大的表面积，受到流体更大的阻力，促使颗粒跟随流体移动。不同形状因子下的颗粒轨迹如图23所示。由图23可知：随着颗粒形状因子增大，颗粒偏离流线的程度增大。颗粒能够与流体保持相对静止运动的临界颗粒直径减小，因此，能够使颗粒翼型升力系数大于无颗粒时流体翼型升力系数的临界颗粒直径减小。

图22 不同颗粒形状因子下翼型升力系数随粒径的变化Fig. 22 Change of airfoil lift coefficient with particle diameter under different particle shape factors

图23 不同颗粒形状因子下颗粒轨迹Fig. 23 Particle trajectories under different particle shape factors

图24 所示为翼型升力系数随颗粒形状因子的变化，由图24 可知：随着颗粒形状因子增加，翼型升力系数减小。图25 所示为不同颗粒形状因子下沿相同编号颗粒轨迹运动的颗粒速度与无颗粒流体速度。由图25 可知：随颗粒形状因子增加，颗粒偏离流体流线的程度增大，Up-fn增大。图26所示为不同颗粒形状因子下沿相同颗粒编号轨迹运动的颗粒速度和无颗粒流体速度之差。由图26可知：颗粒对流体的影响增强，沿同一颗粒轨迹的流体速度变化程度减小，Ufp-fn增大。

图24 翼型升力系数随颗粒形状因子的变化Fig. 24 Change of airfoil lift coefficient with particle shape factor

图25 不同颗粒形状因子下沿相同编号颗粒轨迹运动的颗粒速度与无颗粒流体速度Fig. 25 Particle velocity and the particle-free fluid velocity along the same particle ID trajectory under different particle shape factors

经分析可得，随着颗粒形状因子增加，相同颗粒编号的颗粒远离压力表面。在沿同一编号颗粒轨迹运动的流体速度变化值相同的情况下，形状因子大的颗粒升力系数(与无颗粒升力系数相比)降低程度更大，升力系数较小。而且随着颗粒形状因子增加，Ufp-fn增大，沿同一编号颗粒轨迹运动的流体速度变化值增大。故随颗粒形状因子增大，升力系数减小。当形状因子为0.9、颗粒直径为2 500 μm 时，升力系数减小幅度最大，与无颗粒时相比减少了0.18%。

3.3.2 颗粒形状因子对翼型阻力系数的影响

图27所示为翼型阻力系数随形状因子的变化。由图27 可知：翼型阻力系数随颗粒形状因子的增大而增大；随着颗粒形状因子的增加，Ufp-fn增大，等效于相对流体速度的减少程度增大，导致阻力系数增大；当形状因子为0.9，颗粒直径为1 500 μm时，阻力系数变化最大，与无颗粒时相比增加了1.315%。