垂直轴水轮机CFD数值模拟有效性验证

， ， ，，

(1.中海油研究总院，北京 100027;

2.中国石油大学博士后流动站，北京 102200;

3.中国石油大学机械与储运工程学院，北京 102200;

4. 哈尔滨工程大学海洋可再生能源研究所，哈尔滨 150001)

0 引 言

潮流能作为一种新兴的清洁、可再生能源正在为世界各种争相开发，而垂直轴水轮机作为潮流能转换的核心装置之一，已成为潮流能开发的主要途径。水轮机技术的研究主要来自实验验证与计算机的数值模拟。垂直轴水轮机和垂直轴风力机基本原理是相同的，在数值模拟方面可以相互借鉴。2005年，Horiuchi等人[1]使用STAR-CD 模拟二维假设下的垂直轴风力机速度场，并与试验值进行比较。2007年，Guerri等人[2]使用STAR-CD模拟二维垂直轴风力机叶片的受力。2008年孙科提出应用滑移网格理论研究H型叶轮的动边界CFD数值模拟问题[3-4]。2009年，Hwang等人[5]应用STAR-CD，采用k-ε湍流模型对摆线式垂直水轮机进行了研究。

CFD数值模拟的方法在水轮机技术研究过程中得到了广泛的应用，但是模拟结果的精度及可信度一直是所有研究者共同关注的焦点。尽管垂直轴水轮机或风力机CFD数值预报方面取得了很大的发展，但是有限的精度仍旧迫切需要得到改进。这其中如物理模型、网格划分、边界条件、湍流模型、CFD求解器等因素引起的计算结果中的误差，都是不能被忽略的。因此，本文针对CFD方法在垂直轴水轮机数值模拟中的有效性进行了详细的研究，包括湍流模型、网格尺度，以及时间步长对计算精度的影响。

1 CFD理论基础

计算流体力学(Computational Fluid Dynamics，简称CFD)是流体力学研究领域中的重要技术之一，其基本思路为使用数值方法在计算机中对流体力学的控制方程进行求解，从而达到模拟预测流场的目的。下面介绍本文中涉及的湍流模型及其特点，公式的推导过程可以参考Ansys CFX用户使用手册[6]。

1.1 k-ε模型

(1)标准的k-ε模型：

标准k-ε模型：是从实验现象总结出来的半经验公式，也是最简单的两方程模型，主要基于湍流动能k和湍流耗散率ε。k方程是个精确方程，ε方程是个由经验公式导出的方程。该模型假定流场完全是湍流，分子之间的粘性可以忽略。因此，只对完全是湍流的流场有效。

(2)RNGk-ε模型：

1.2 k-ω模型

(1)标准的k-ω模型：

两方程模型，包含了低雷诺数影响、可压缩性影响和剪切流扩散，因此适用于尾迹流动计算、混合层计算、射流计算以及受壁面限制的流动计算和自由剪切流计算。

(2)SSTk-ω模型：

它综合了k-ω模型在近壁区计算的优点和k-ε模型在远场计算的优点，将k-ω模型和标准k-ε模型都乘以一个混合函数后再相加就得到这个模型。与标准k-ω相比，SSTk-ω模型中增加了横向耗散导数项，同时在湍流粘性定义中考虑了湍流剪切应力的输运过程，模型中使用的湍流常数也有所不同。这些特点使得 SSTk-ω模型的适用范围更广，比如可以用于带逆压梯度的流动计算、翼型计算、跨音速激波流动计算等等。

1.3 雷诺应力模型(RSM)

雷诺应力模型中没有涡粘度的各向同性假设，因此从理论上来说湍流模式理论要更精确。从理论上说，RSM 应该比一方程和两方程模型的计算精度更高，但实际上雷诺应力模型的精度受限于模型的封闭形式，因此雷诺应力明显具有各向异性的特点时才必须使用 RSM，比如龙卷风、燃烧室内流动等带有强烈旋转的流动问题。

2 垂直轴水轮机的CFD有效性验证

2.1 计算算例

1983年Strikland等人[5]等人做的一系列研究垂直轴风力机动态失速试验中，测量和记录的数据包括叶片的瞬时切向和法向受力、叶片表面压力分布以及尾流流场等，试验结果作为一个可靠的参考依据。 Deglaire等人[7]、Hwang等人、Tchon和Paraschivoiu[8]和汪鲁兵[9]的模拟研究中都是与其比较得出结论的。本文研究同样采用Strikland等人的试验模型为验证算例。模型的具体参数如表1所示。

表1 Strickland试验轮机的主要参数

2.2 网格划分

网格划分是CFD模拟过程中较为耗时的环节，也是直接影响模拟精度和效率的关键因素之一。对于所有的CFD数值模拟，结构网格模拟精度和速度要明显优于非结构网格情况，因此本文系采用结构化网格对流场进行划分。根据垂直轴水轮机的结构特点，叶片在展长方向上具有相同的横截面，为了节省计算时间，本文简化为二维模型。

水蛭交配一般发生在早上或夜间，交配后一个月左右，雌体生殖器分泌出稀薄的黏液，中包被卵带，形如“蚕茧”，排出体外，在湿泥中孵化。产茧一般在土壤（含30%水分）中进行，产茧时间一般为4月中下旬温度在20-25℃以上。

划分网格之前，还需要对模型中的流场进行基本的预估。过疏或过密的网格都会对计算有极大影响。过疏的网格往往会得到不精确甚至完全错误的结果；过密的网格会使计算量增大，使计算难以收敛。文中的计算网格如图1和图2所示。

图1 Turbine 网格分布图

图2 Turbine 叶片周围网格分布图

2.3 边界条件的设置

流体域内是以大气压为参考气压的；在进口边界处设置为速度进口条件，给定均匀来流速度V及湍流参数；出口边界采用开敞边界；流体计算域的侧面采用滑动墙面(slip wall)设定，这种墙面设定可以有效地提高计算速度，叶片表面为不可滑移壁面；旋转域给定旋转角速度，旋转域和静止域通过交界面连接；使用入口处流动参数来初始化流场；采用瞬态求解器。

图3 边界条件设定

3 结果对比

3.1 不同湍流模型计算结果对比

图4给出了标准、RNG 、SST和Omega雷诺应力四种湍流模型计算得到的水轮机在速比为5.0时，叶片法向力系数和切向力系数随位置角的变化历程与试验测量结果的比较。从总体来看，不同湍流模型计算结果和试验结果趋势基本一致，法向力系数要比切向力系数吻合的好。但是在数值上，不同湍流模型结果还是存在一些差异。

从图4(a)法向力系数曲线可以看出，结果的主要差别出现在谷值和峰值上，四种湍流模型模拟结果谷值均大于试验结果，雷诺应力模型结果偏差更大一些，其他三种模型结果相当；模拟结果峰值均小于试验值，SST模型模拟结果与试验值更为接近，其次是雷诺应力模型，标准模型和RNG 模型结果相当。从4(b)切向力系数曲线可以看出，四种湍流模型模拟结果在上游盘面()与试验结果符合程度好于下游盘面时的情况，这主要是因为下游盘面经过叶片扰动，流场变得更加复杂，使模拟精度下降。在上游盘面谷值处四种湍流模型模拟结果很接近，均略大于试验结果，同时与试验结果相位也有偏差；上游盘面峰值处，RNG 模型结果高于试验结果，标准模型结果低于试验结果，SST模型和雷诺应力模型结果与试验结果都比较接近；下游盘面谷值处四种湍流模型模拟结果很接近，但是都小于试验值，同样也存在相位偏差；下游盘面峰值情况，四种湍流模型模拟结果均出现了平峰，这一点与试验值差别较大，SST模型模型结果更接近试验值。

总的来说，采用SST模型所计算的结果与试验结果吻合更好，这也进一步体现了SST模型对逆向压力梯度与流动分离的流场模拟的优势，虽然采用雷诺应力模型计算结果与SST模型计算的准确性相差不大，但是雷诺应力模型比SST模型计算耗时更长，因此在本文以后的计算中选用SST模型。

3.2 网格尺度对计算结果的影响

网格的疏密程度与计算迭代的收敛速度以及最终的计算精度紧密相关。网格节点数量大则意味着计算结果精度较高，但这是以增加计算时间为代价的。通常情况下，网格的生成过程中人的主观因素对网格的质量将产生很大的影响。 CFD 用户在生成网格时仅仅凭个人经验，在平衡计算时间和精度的前提下确定网格分布的拓扑结构、密度等。文中需要对垂直轴水轮机非定常绕流流场做大量计算，希望网格在保证一定精度的前提下，网格数尽可能少。因此在开展数值计算前，有必要对网格进行测试，测试具有不同疏密程度的网格系统对垂直轴水轮机计算结果的影响。

垂直轴水轮机数值模拟中，实质是叶片的绕流问题，因此叶片周围网格最为重要，为了进行网格测试，叶片周围生成4套密度不同的“⊙”形结构网格，周向节点布置相同，仅变化径向的网格密度，相关信息见表2。

表2 不同网格结构相关信息

仍然以水轮机在速比5.0时的试验结果与数值结果进行对比，以此判断所选网格密度是否合理。图5给出了四种不同网格密度(Mesh 1～4)计算得到叶片法向力系数和切向力系数随位置角的变化历程与试验测量结果的比较。从图5(a)法向力系数的对比曲线看出，网格疏密对法向力的影响不大，四种不同网格密度的计算结果十分接近。而网格疏密对切向力的影响比较大，从图5(b)切向力系数的对比曲线看出，随着网格密度的增加而呈越来越接近试验结果的趋势，但是Mesh2和Mesh1的结果差别不大，基本满足网格无关性的要求。

图5 叶片法向力和切向力系数随位置角的变化(λ=5.0)

这也和文献推荐值很吻合，当时可以充分利用低雷诺数模型预测分离流的优势。鉴于速比5.0是一个比较高的速比，分离流动还不是很明显，因此文中后面计算全部采用网格密度类型Mesh2。

3.3 时间步长对计算结果的影响

在非稳态问题求解时，时间步长影响模拟的收敛性及CPU计算时间。正如存在网格无关解一样，非稳态数值模拟应当存在时间步长无关解：数值解在工程允许的偏差范围内不再随时间步长的减小而发生变化。

但是，时间步长的减小导致了计算时间的迅速增加。由于描述计算对象的微分方程的刚性较强，较小的时间步长有助于得到较为准确的结果，且能提高时间步内的收敛速度；但是若计算的总时间一定，这并不足以抵消时间步数增加而引起的计算量增加。时间步长往往需要在计算精度和收敛速度之间进行折衷。下面采用简单的增加时间步长的方法来分析时间步长对垂直轴水轮机数值模拟结果的影响，时间步长主要根据叶轮每步转过角度来设定，采用上节中介绍的Mesh2网格模型，时间步长设置见表3。

表3 Strickland实验轮机的主要参数

图6给出了四个不同时间步长计算得到的水轮机在速比为5.0时，叶片法向力系数和切向力系数随位置角的变化历程与试验测量结果的比较。从图中可以看出，时间步长对计算结果的影响和网格影响很相似，都是对法向力影响比较小，对切向力影响比较大，湍流模型的影响也有类似结论，究其原因：一是因为叶片切向力系数的大小比法向力系数小一个量级；二是切向力系数的变化规律比法向力系数复杂，在一个周期内，切向力系数在上游盘面和下游盘面各出现一次峰值和谷值，即其变化频率比法向力系数大。从图6(b)叶片切向力系数曲线中可以看出，时间步长越大，计算结果与试验结果的偏差越大，但是Tstep2与Tstep1的计算结果相差无几，可以说当叶轮每个时间步转过角度小于3°时，计算结果基本达到了时间步长无关性。

图6 叶片法向力和切向力系数随位置角的变化(λ=5.0)

4 结 论

通过对比以上计算结果可以得到以下结论:

(1)通过与试验结果对比，验证了CFD方法在垂直轴水轮机性能预报方面应用的可靠性；

(2)湍流模型选择方面，和其他湍流模型相比，SST湍流模型计算精度更高；

(3)当满足时基本满足网格无关性要求；

(4)时间步长设定时，当叶轮每个时间步转过角度小于3°时，计算结果基本达到了时间步长无关性要求。

[1] Horiuchi, K., Ushiyama, I., and Seki, K., Straight Wing Vertical Axis Wind Turbine: A Flow Analysis, Wind Engineering, Vol. 29, No.3, pp.243-252, May 2005.

[2] Guerri, O., Sakout, A. and Bouhadef, K., Simulations of the Fluid Flow around a rotating Vertical Axis Wind Turbine, Wind Engineering Vol.31, no. 3, pp 149-163, 2007.

[3] 孙 科, 竖轴H型叶轮及导流罩流体动力性能数值模拟[J].哈尔滨工程大学,2008.

[4] 姜子豪，李浩，阚江明.减振降噪在机械设计中的应用[J].森林工程，2013，29(5)：55-59.

[5] Hwang, S., Lee, Y. H., and Kim, S. J., Optimization of Cycloidal Water Turbine and the Performance Improvement by Individual Blade Control, Applied Energy, Vol. 86, Issue 9, 2009.

[6] Ansys Inc: “ANSYS CFX User Manual”, 2005

[7] Deglaire, P., Engbom, S., Agren, O., and Bernhoff, H., Analytical Solutions for a Single Blade in Vertical Axis Turbine Motion in Two-Dimensions, European J. of Mechanics B/Fluids, Volume 28, Issue 4, pp. 506-520, July-August 2009.

[8] Tchon, Ko-Foa and Paraschivoiu, I., Navier-Stokes Simulation of the Flow Around an Airfoil in Darrieus Motion, Journal of Fluids Engineering, Vol. 116, pp. 870-876., ASME, 1994.

[9] 汪鲁兵.竖轴潮流水轮机水动力性能理论与实验研究[J]. 哈尔滨工程大学,2006.