垂直轴潮流能水轮机直叶片流固耦合分析

李 霆，陈 兵

（大连理工大学 水利工程学院，辽宁 大连 124221）

垂直轴潮流能水轮机直叶片流固耦合分析

李 霆，陈 兵*

（大连理工大学 水利工程学院，辽宁 大连 124221）

在开放海域运行的潮流能水轮机需要应对复杂多变的海洋环境，要有较好的刚度、强度和抗疲劳性能。这对水轮机结构，尤其是为其运行提供动力的叶片结构的设计提出了较高要求。文中建立了垂直轴潮流能水轮机叶片及流场的三维数值计算模型，基于ANSYS商业软件提供的CFX与瞬时结构模块的双向流固耦合分析功能，分析了叶片在不同流速下的流体力学性能，分析了流场的三维效应及其引起叶片升力损失的原因，研究了其结构的刚度和强度特性。还对比分析了使用铝合金和钢材作为水轮机叶片材料的优劣。

数值模拟；潮流能水轮机；叶片；流固耦合；三维效应；结构分析

垂直轴潮流能水轮机具有结构简单、可以利用各方向来流无需换向、发电机密封成本低、噪音较小等优点[1-2]。但垂直轴式水轮机的运行方式导致其叶片会受到来自潮流作用的周期性荷载，易发生疲劳破坏[3]。因此分析垂直轴潮流能水轮机叶片结构的应力响应对于叶片的设计有一定价值。

水轮机叶片截面一般选用风力发电机的翼型，或在其基础上进行改进[4]。实践中各国研究人员开发了多种叶片结构形式和材料。中国海洋大学王树杰等[5]研究了柔性材料叶片的应用，显示柔性材料可以有效提高水轮机获能利用率。李冬等[6]提出了采用三维CAD、反求技术和快速原型技术快速制造潮流能水轮机叶片原型的方法。Gaurier等[7]研究了DCPD树脂叶片在水流单独作用和浪流同时作用下的应力应变规律，并指出：浪流同时作用下叶片应力有显著提高，在疲劳分析中需充分考虑其作用。应用较为广泛的是源自风力发电机叶片设计的蒙皮、主梁加腹板空心结构，该结构构造简单，力的传导清晰明确，便于设计和计算；其材料可以是铝、钢材、复合材料等。Harper等[8]提出了一种研究复合材料水轮机叶片的裂纹生长和层合板失效规律的方法。Ng等[9]分析了使用高强度碳纤维制作潮流能水轮机叶片的可能性。

随叶轮旋转，水轮机叶片承受的水流作用力呈周期性变化。采用流固耦合方法可以研究叶片受水流作用下的瞬时响应。荆丰梅等[10]使用单向流固耦合计算方法，将二维CFX流体计算结果作为水轮机受力，将流体与结构中相对应的节点压力施加到结构有限元模型上，不考虑结构变形对流场的扰动，研究了垂直轴潮流能水轮机的力学性能，其结果与实验吻合较好。刘雪峰等[11]使用单向流固耦合计算方法分析了潮流能水轮机叶片的强度，得到了叶片的各阶模态振型、频率及最大变形量。Singh等[12]使用流固耦合方法分析了韩国某1 MW水轮机叶片结构的强度，并对比分析了钢和铝两种材料叶片受力的差异。在每一个循环中，叶片承受的最大作用力与水流的攻角和流速有关。对于采用NACA0018翼型的叶片，当雷诺数在0.25×106～2×106范围内时，在水流攻角为12°附近叶片会受到较大升力，研究此时叶片结构的应力响应对叶片的设计有一定价值。本文选用蒙皮、主梁、腹板结构建立了垂直轴潮流能水轮机叶片模型，使用流固耦合方法分析其受12°攻角的不同流速水流冲击下的应力响应。

1 流固耦合计算原理

1.1 流场控制方程及湍流模型

本文中，CFX使用不可压的雷诺时均N-S方程为流场控制方程，其表达式为：

式中：U，p为时均的速度和压力值；ρ为密度；u为速度脉动项；SM为变形率张量；τij为分子粘性应力张量；uiuj为雷诺应力项，下标i,j=x,y,z表示空间坐标方向。为解决方程的封闭问题，引入湍流模型以构造脉动项与时均项的关系式。基于平均方程和湍流模型的研究方法可以有效模拟小尺度涡的湍流运动。

目前研究人员已开发了多种湍流模型，其中基于Boussinesq涡粘性假设的二方程模型运用最为广泛。k-ωSST模型是一种混合模型，它在边界层内使用k-ω模型计算，在外部流场使用k-ε模型计算。同时通过对粘度系数的限幅，k-ωSST模型修正了k-ω模型未考虑剪切压力输送的问题，可以更精确地预测涡的分离，不可压缩流k-ωSST模型的湍动能和比耗散率方程为：

其中，Pk是湍动能的产生项：

式中：参数 φ3（φ 代指 α，β,σk，σω）是 k-ω 模型和k-ε模型中参数的线性组合：

式中：F在边界层内取值为1，边界层外取值为0；式中参数见下表：

β' α1 β1 σk1 σω1 α2 β2 σk2 σω2 0.09 5/9 0.075 2 2 0.44 0.082 8 1 1/0.856

k-ωSST模型修正后的粘度系数表达式为：

其中，vt=μt/ρ；S 为不变测度的应变率。

1.2 结构应力应变控制方程

本文中使用的材料均为均质线性材料，其应力应变关系可由式(7)表示：

式中：{σ}=[σxσyσzσxyσyzσxz]T为应力向量，[D]为弹性刚度矩阵，表达式如下：

式中：Ei为杨氏模量；vij为泊松比；Gij为剪切模量；下标i,j=x,y,z表示空间坐标方向。

{εel}={ε}-{εth}为弹性应变向量，其中{εth}为温度应变向量，{ε}为应力应变向量，本文中未考虑温度的影响，所以应变向量表达式如下：

2 ANSYS流固耦合计算方法

2.1 耦合计算流程

ANSYS提供了多场耦合计算（MFX）功能，该功能可以实现CFX和瞬时结构模块的双向流固耦合计算。有别于单向流固耦合方法仅考虑流场对结构的作用力，双向流固耦合方法还考虑了结构变形对流场的干扰作用，可以更真实地模拟流场和结构的瞬时变化。一个典型的流固耦合模型包含流场和结构场两个场，两个场有至少一个耦合面。耦合计算中，流场、结构场和耦合面处都有独立的收敛目标。流场计算时间步为Δtf，结构场计算时间步为Δts，两个场间隔时间Δtm耦合一次，它们之间满足Δtf=Δtm，Δtm=nΔts，n=1,2,3…。一次耦合中对结构场的多时间步求解可以帮助复杂模型收敛，提高计算的精度和稳定性。

耦合计算的累积时间步反映了真实的模拟时间t，即tn=nΔtm。tn时刻开始的一个耦合时间步的计算流程如下：

（1）求解流场，收敛后传递耦合面处各向力给结构场；

（2）求解结构场，收敛后返回耦合面处位移给流场，耦合面处力和位移的传递通过节点插值实现；

（3）检查流场、结构场和耦合面处是否均已收敛，若收敛则开始下一个耦合时间步的计算（tn+1=tn+Δtm），若不收敛则返回tn时刻重新计算。

可见，在一个耦合计算时间步内流场和结构场会反复求解数次，导致流固耦合计算有占用资源高、计算时间长、不易收敛的问题。为减少计算量，加速收敛，在耦合计算初期仅计算流场，耦合面处不发生力和位移的传递，待流场流态稳定后再加入结构场计算。为避免流场计算初期扰动的影响，使用静止流场开始计算，控制入口流速从零开始缓慢加速至稳定流速。

2.2 耦合计算模型

叶片截面选用NACA0018翼型，材料为均质钢材，弦长为0.5 m，展弦比为4。叶片蒙皮厚2 mm，由等间距的三根梁和三块腹板支撑。叶片在第二腹板距前缘1/4处使用4根直径5 mm的螺栓固定在旋臂上，计算中假设旋臂为刚体，固定处的约束采用圆柱面约束（Cylindrical Support），禁止轴向和径向的位移，允许切向旋转。由于叶片在展向的几何对称性，基于对称性原理，假设叶片内部应力分布及流场流态关于叶片中部横截面对称，为节省计算时间使用叶片展长的一半建模计算，如图1所示。叶片剖面处的边界条件为无摩擦支撑条件：允许叶片在剖面所在平面内平移和旋转，限制叶片剖面脱离原平面的各种运动。

水动力计算使用长方体计算域，入口边界与侧方边界各距叶片10倍叶片弦长，出口边界距叶片20倍叶片弦长，顶部边界距叶片2倍叶片弦长。为保证边界层处网格保持原有形状随叶片变形运动，动网格参数为各节点参照前一时间步的位置移动和增加小体积网格的刚度。叶片表面是流场和结构场的耦合面，其边界条件为无滑移边界条件。入口为速度边界条件，出口边界压力为0。其余边界均为对称边界条件。

图1 叶片结构模型图

耦合计算为瞬态计算，CFX使用二阶隐式欧拉法求解，对流项离散格式为高分辨率格式。ANSYS结构求解方式选择直接求解，迭代算法采用完全Newton-Raphson算法，该算法在每一次平衡迭代后都会修改刚度矩阵。

3 网格独立性、时间步无关性验证

本文使用结构化网格划分流场，按照叶片表面首层网格y+值为5、15、30分别生成了规模为50万、35万、20万的3种流场网格。由于叶片内部结构较为复杂，在其固定位置和梁板连接的倒角处难以使用较高质量的结构化网格模拟，所以本文使用非结构化网格以更好地模拟叶片的结构细节，防止因模型失真可能引发的应力集中。通过在叶片固定处、叶片前缘处及梁、板、蒙皮的连接处加密，得到了规模分别为13万、35万和50万的3种结构场网格。

图2 流场网格

图3 叶片附近流场网格

图4 叶片表面网格

图5 开始耦合后叶片最大应变随时间变化曲线

选取水流流速为2 m/s，分别计算叶片的升力系数和叶片受水流作用产生的最大变形和最大应力，以比选合适的网格数量，结果见表1。耦合计算中叶片最大应力、应变出现在叶片固定位置附近蒙皮与腹板的连接处，叶片最大应变随时间的变化见图5。耦合计算开始时，叶片突然承受流场施加的荷载会产生较大应力，同时伴随着较大的位移，随时间发展初始扰动的能量被逐渐消耗，叶片最终在平衡位置处作周期性振动。

表1 网格试算中升、阻力系数及最大应力统计表

表1中，各项误差率的计算均以流场网格数55万、叶片网格数50万的计算结果为参照。叶片关键部位网格加密后，叶片的应力、应变和变形的计算精度都有了大幅度提高，但叶片的升力系数和阻力系数并无明显变化。这是因为叶片刚度较高，在水流作用下变形很小，未对流场产生明显扰动。边界层处网格加密后，叶片的升力系数和阻力系数的计算精度明显提高，叶片的应力、应变和变形的计算精度也有一定提升。在流场和叶片的网格数均为35万时，各计算参数的误差率均在1%以内，可以满足计算精度的要求。

为选取合适的时间步长，使用从0.01～0.001 s的4种不同时间步试算，计算结果并未发生较大改变，选择使用0.005 s作为计算的时间步长。

4 流速变化下叶片的水动力性能

4.1 叶片升力系数和阻力系数变化规律

相关实验和计算显示，雷诺数的变化会显著影响叶片的力学性能，使叶片升力、阻力系数随之改变。本模型中，入口流速从0.5 m/s增加至4 m/s时，流场雷诺数从0.25×106提高到2×106，本文选择其中间距为0.5 m/s的8种流速计算。本文设置了对比组来验证计算的可靠性，建立叶片与流场等高的模型，使流场处于准二维状态，其计算结果应接近二维实测值[13]。结果显示，准二维流场中叶片各截面受到的升力和阻力相同，各流速下叶片升力系数和阻力系数计算结果与二维实测值吻合较好。叶片升力系数和阻力系数随雷诺数变化曲线见图6～图7。

图6 升力系数随雷诺数变化的规律

图7 阻力系数随雷诺数变化的规律

图6显示了叶片升力系数随雷诺数的变化规律。随雷诺数增加，叶片升力系数随之变大，与二维实测值变化趋势相同。但是，三维流场中叶片升力系数增加的幅值远小于二维实测值的增加值。与准二维流场不同，三维流场的叶片高度小于流场高度，水流绕过叶片端部后会流向叶片背流面低压区域，流场中有不可忽视的三维流动。叶片端部的三维流动导致叶片背流面压力升高，叶片两侧压差减小，使叶片的升力系数小于二维实测值。随流速增加，流场的三维效应变强，三维流动影响了更大的范围，使叶片升力系数衰减比例增大。雷诺数为0.25×106时（流速0.5 m/s）升力系数比二维实测值衰减约11.8%，随雷诺数增加，升力系数的衰减比例提高至30%。

图7显示了叶片阻力系数随雷诺数的变化规律。随雷诺数增加，叶片阻力系数稍稍变大，与二维实测值变化趋势相反，可见在该范围内，流速越大叶片受到的阻力越高。计算所得阻力系数远大于实测值，随雷诺数增大，它与二维实测值的比值从2增长到3.4。

4.2 三维效应的分析

由于叶片的展弦比数值较小，叶片周围流体有较强的三维效应，使叶片升力系数不达预期。观察流场发现，叶片端部有涡产生，叶片近端部的流场受涡影响流态发生改变，导致叶片近端部的升力有较大损失。

图8为流速为4 m/s时叶片迎流面和背流面的压力分布图。从叶中剖面处到端部，叶片迎流面正压基本保持不变，直至到达端部时断崖式下跌至负压；背流面负压的值在叶中剖面处达到最大值，在靠近端部过程中持续加速下降，在叶片端部达到最小值，这也使整个叶片升力系数下降、阻力系数提高。

随流速增加，流场的三维效应变强，叶片升力系数衰减比例增大。为研究流场三维流动影响的规律，提取从页中剖面处至叶片上端不同高度处截面的升力和阻力，计算出各截面处升力系数在不同流速时的损失比例，绘制了图9，其横轴为截面高度与展长之比β，竖轴为该截面处升力相对于二维实测值的损失比例。当流速为0.5 m/s时，在远离端部的叶中剖面处，叶片的升力系数计算值与实测值基本吻合，升力损失比例仅为1%，这里的流场流态与理想二维平面流动基本一致；从叶中剖面到端部，升力损失比例迅速增大，叶端的升力损失达到了50%。流速增大后，三维流动扩散到了更大范围的流场，叶中剖面处也逐渐受到三维流动的影响，这里的升力损失比例逐渐增大到4 m/s时的20%；叶片端部的升力损失比例随流速增大也有一定增加，其值从0.5 m/s时的50%增加至4 m/s时的60%，可见在低流速时叶片端部附近流场的三维流动已充分发育，流速增加后并未使端部升力的损失显著上升。观察图9发现，流速增大后升力衰减曲线斜率变小，最大值缓慢增加，最小值迅速增加，曲线趋向平坦，显示三维流动更加平均地影响到了更大范围的流场。所以，流速增大后，叶片升力损失提高的原因主要是三维流动影响范围的扩大，同时伴随着现有影响范围内升力的进一步损失。

图8 流速为4 m/s时的压力分布图

图9 不同流速下叶片升力衰减曲线

5 流速变化下叶片结构力学性能

5.1 钢制叶片的力学性能

表2列出了各流速时叶片的最大变形、最大应力和最大应变。可见，随着流速增大，叶片承受的作用力不断提高，叶片内部最大应力、应变随之变大，叶片变形也越来越严重。对表2数据分析发现，叶片的最大变形、应力、应变都近似与流速的二次方成正比。叶片的升力（阻力）公式中包含流速的平方项和升力系数（阻力系数）的一次方项，即叶片的升力（阻力）与速度的二次方成正比，与升力系数（阻力系数）成正比。本文中流速提高了8倍的同时升力系数和阻力系数却并未发生较大改变，所以叶片的变形、应力和应变与流速的二次方近似成正比关系。

表2 不同流速下叶片变形及最大应力、应变统计表

随着流速增加，叶片变形形状及叶片较大应力、应变分布位置并未发生改变。图10～图11分别给出了流速为2 m/s时叶片的整体变形示意图和叶片表面应力分布示意图。

图10 叶片在2 m/s水流作用下整体变形图

图11 叶片在2 m/s水流作用下应力分布图

图10是叶片整体变形图，该图中线框为叶片未变形位置，实体为叶片变形后形状。叶片变形主要集中在蒙皮中间未被支撑的位置和叶片端部。叶片的最大形变发生在叶片前缘背流面两肋之间的蒙皮处，这是因为该位置迎流面和背流面之间的压差很大，是叶片大部分升力的来源。叶片端部和剖面处都产生了向背流方向的位移，其中端部位移相对较大，这是因为叶片端部是悬臂结构，更易产生变形。总之，由主梁和腹板组成的框架结构为叶片提供了很好的骨架，整体而言，叶片刚度足够，结构未发生明显变形，未因形变影响叶片的获能能力。但在叶片设计中仍需对叶片前缘蒙皮作足够的支撑，对承受鱼类、漂浮物撞击等意外作用作好预判，保障叶片长期有效运行。

图11是叶片应力分布图，可见较大的应力多集中在蒙皮与梁、腹板的连接处，其中最大应力出现在叶片前缘背流面的蒙皮与腹板连接处。在流速为4 m/s时，叶片内部最大应力为16.8 MPa，并未超出钢材的线弹性极限，即叶片只发生了线弹性变形，叶片强度满足要求且还有一定的强度储备。但考虑到蒙皮与腹板、梁连接处会存在焊缝，水轮机在其25 a的设计寿命中会经历高达108次应力循环，这些位置在高周循环应力的作用下仍有疲劳失效的可能，在叶片设计、制作中需着重关注。

5.2 铝合金叶片与钢叶片的对比

铝合金的密度较低，仅为钢密度的1/3左右，可以使叶片重量减轻50%以上，节省材料；铝合金的强度也比较高，优质铝合金的强度可以接近或超过钢材；铝合金还有较好的塑性，可以加工成各种形状，可以满足叶片加工的需求；另外，不同于钢材极易锈蚀的缺点，铝合金的抗腐蚀性优良，对于海洋环境中的应用有先天优势。

改变叶片材料为铝合金，参照上文设置重新计算，得到铝合金叶片的升力系数、阻力系数和最大变形、应力、应变，见表3。

分析表3数据发现，材料的变化并未使叶片的升力系数和阻力系数有明显改变。由于来自水流的升力和阻力不变，铝合金叶片的最大应力也没有变化，叶片仍在线弹性范围内。与钢叶片相比，铝合金叶片的变形和应变都有所增大，其数值大约是钢叶片数值的3倍。在流速为4 m/s时，铝合金叶片的最大变形约为0.2 mm，叶片变形虽有所增大，但强度仍然满足使用要求。

表3 铝合金叶片升、阻力系数和最大变形、应力、应变统计表

6 总结

本文使用双向流固耦合计算方法研究了垂直轴潮流能水轮机叶片的流体力学性能和叶片结构的力学性能，升力系数、阻力系数计算结果与实测值吻合较好，显示计算有较高可靠性；发现叶片腹板、主梁与蒙皮的连接处是叶片的应力集中区域，对工程中的应用有一定借鉴意义。本文研究了流场的三维效应导致叶片升力损失的规律，发现流速增大后，叶片升力损失提高的原因主要是三维流动影响范围的扩大，同时伴随着现有影响范围内升力的进一步损失。本文还对比分析了铝合金和钢材在水轮机叶片中的应用，铝合金材料可以在保证强度的前提下有效减轻叶片重量、节约材料，并提高叶片抗腐蚀性能，但其变形和应变是钢材的3倍。未来，可以基于双向流固耦合计算方法进一步研究垂直轴潮流能水轮机运行中水轮机与水流的瞬时相互作用，分析叶片结构的应力循环疲劳特性，在设计中预测其使用寿命。

[1]张亮，李志川，张学伟，等.垂直轴潮流能水轮机研究与利用现状[J].应用能源技术，2011，9：1-7.

[2]马勇，张亮，马良，等.垂直轴水轮机式潮流能发电装置开发现状与发展趋势[J].科技导报，2012，30（12）：71-75.

[3]Lago L，Ponta L，Chen L.Advances and trends in hydrokinetic turbine systems[J].Energy for Sustainable Development，2010，14：287-296.

[4]BahajA，Moland F，Chaplin R.Hydrodynamicsofmarine current turbines[J].RenewableEnergy，2006，31（2）：249-256.

[5]王树杰.柔性叶片潮流能水轮机水动力学性能研究［D］.青岛：中国海洋大学，2009.

[6]李冬，王树杰，袁鹏，等.潮流能水轮机多种翼型系列叶片优化设计中的实验模型制作方法研究[C]//全国海洋能学术研讨会，2009：2541-2544.

[7]Gaurier B，Davies P，DeuffA.Flume tank characterization ofmarine current turbine blade behaviour undercurrent and wave loading[J].Renewable Energy，2013（59）：1-12.

[8]Harper P，Hallett S.Advanced numerical modelling techniques for the structural design of composite tidal turbine blades[J].Ocean Engineering，2015（96）：272-283.

[9]Ng K，Lam W，Pichiah S.A review on potential applications of carbonnanotubesin marine current turbines[J].Renewable and Sustainable EnergyReviews，2013（28）:331-339.

[10] 荆丰梅，肖钢，熊志民.潮流能水轮机单向流固耦合计算方法[J].振动与冲击，2013，32（8）：91-95，105.

[11]刘雪峰，王晋宝，田美灵.基于单向流固耦合的潮流能水轮机叶片强度及模态分析[J].制造业信息化，2014，11：94-98.

[12]Singh P，Choi Y.Shape design and numerical analysis on a 1 MW tidal current turbinefor the south-western coast of Korea[J].Renewable Energy，2014（68）：485-193.

[13]訚耀保.海洋波浪能综合利用：发电原理与装置[M].上海：上海科学技术出版社，2013：342-345.

Analysis on a Two-Way Fluid-Solid Interaction of VATT Straight Blades

LI Ting1,CHEN Bing2
1.School of Hydraulic Engineering,Dalian University of Technology,Dalian 116024,Liaoning Province,China;
2.School of Ocean Science and Technology,Dalian University of Technology,Dalian 116024,Liaoning Province,China

Designed to be situated in the open sea,tidal turbines must cope with the complex and volatile marine environment,and this poses a strict requirement for its structural stiffness,strength and anti-fatigue property.This requires high-standard design of the turbine structure,especially for the blades that power the turbine.Based on the two-way fluid-solid coupling analysis function of CFX and transient structure module provided by ANSYS commercial software,a three-dimensional numerical simulation model is established for the vertical axis tidal turbine blades and the fluid domain.The hydrodynamic performance of the blades is analyzed at different flow rates.The three-dimensional effect of the fluid domain,the loss of lift force,as well as the stiffness and strength characteristics of the blade structure are also studied in this paper.It also analyzes the use of aluminum alloy and steel as blade materials,and illustrates both their advantages and disadvantages.

CFD;VATT;Blade;two-way FSI;Three-dimensional effect;structural analysis

P743.1

A

1003-2029（2017）05-0112-08

10.3969/j.issn.1003-2029.2017.05.018

2017-04-10

国家自然科学基金资助项目（51379036）

李霆（1993-），男，硕士，主要研究方向为海洋可再生能源。E-mail：thunderbolt_lee@outlook.com

陈兵（1970-），男，博士，副教授，主要从事海洋工程及新能源技术的研究工作。E-mail：chenbing@dlut.edu.cn