基于改进LS-STAG浸没边界法的风力机转子自旋数值模拟

摘要

利用改进的LS-STAG水平集/切割单元浸入边界法，开发了一套用于风力机转子自旋数值模拟的软件。水平集函数用于浸入式边界描述。提出了基于Bézier曲线的复杂翼型水平集函数构造算法。此外，还描述了在不重构Bézier曲线的情况下，对每一个新的转子位置重新计算水平集函数的算法。提出了一种用于气动力矩滤波的二阶巴特沃斯低通滤波器，并采用粗网格进行仿真。为了验证改进的LS-STAG方法，对两叶片自旋萨伏尼乌斯旋翼的流动进行了模拟 ．

1.介绍

转子是风力涡轮机功能元件链中的第一个元件。因此，其气动和动力特性在许多方面对整个系统都有决定性的影响。设计人员面临的问题是找出转子的实际形状，例如，转子叶片的数量或其叶片的翼型与其气动性能之间的关系。

为了模拟转子的动力学[1- - - - - -6特别是它的自旋，这就需要解决耦合的气动弹性问题。这类问题的数值求解比较复杂，因为需要考虑流动与移动的浸入体之间的干扰。因此，本研究的目的是建立一种有效的数值模拟风力机转子自旋的方法。

对于质量较大的物体，当其平均密度明显高于流动密度时，可采用“分步”弱耦合数值算法求解气动弹性耦合问题。首先模拟物体绕流，物体按已知参数运动，然后计算物体在已知水动力载荷下的动力学。本研究就考虑了这种情况。

浸入边界法[7]适用于耦合气动弹性问题的数值模拟，因为它不需要单元边与计算域边界重合，允许在域形状不规则或由于气动弹性体运动在模拟过程中发生显著变化时求解。这些方法的主要优点是网格不需要在每一个时间步进行重建。

在本研究中，LS-STAG切割细胞浸入边界法[8]用于旋翼的自旋模拟。这种方法允许在笛卡尔网格上求解问题。浸入边界用水平集函数表示[9］．由LS-STAG离散化的Navier-Stokes方程或reynolds -平均Navier-Stokes方程得到的线性系统采用BiCGStab方法求解[10]与ILU及Multigrid [11，12预处理。求解器代价系数估计的一种原始算法[13，用于多电网预调节器的最优参数选择。

2.控制方程

该问题是在二维非定常情况下考虑的，假设绕刚性翼型的流动是粘性和不可压缩的。连续性方程和动量方程如下: 在这里是速度矢量，的压力,是时间,是流密度，和为流动运动粘度。计算域外边界的边界条件为: 在这里速度矢量在入口边界和是单位外法向量。翼型表面线上的边界条件为无滑移条件: 在这里为浸入边界的速度。

为了模拟风力机转子的旋转，需要求解以下动力学方程: 在这里为转子的旋转角度，为转子的极惯性矩，粘性摩擦系数在轴上，和为气动力矩。考虑二维问题，所以 (图1）.

3.数值方法

带有单元格的笛卡尔网格 在矩形计算领域中引入。人们认为是…的脸细胞和 是它的中心。未知的组件和的速度矢量在单元面流体部分的中间计算。这些点是控制卷的中心 和 与脸和 ，分别(图2）.这种数值方法的主要思想在我们以前的论文中有描述，例如，在[14］．

浸入式边界相交的单元格称为切割单元格(图3.和4）.这些细胞包含固体部分和液体部分。水平集函数用于浸没边界描述。的边界是由切割单元格上的线段表示的 ．这段端点的位置由变量的线性插值定义 ．

关于水平集函数的构造应该提到一些注意事项。对于复杂形状的旋翼，即Darrieus旋翼，对于圆形翼型和简单的旋翼形状，例如萨伏尼乌斯旋翼，不能解析地构造水平集函数(图)5和6）.因此，需要用曲线来近似转子表面线，这样既可以模拟光滑部分的边界，也可以模拟尖锐的边缘。此外，希望能很容易地计算出从任意点到边界的距离。对于任意形状翼型的水平集函数构造的一种有效方法，它符合上述要求，在[15］．在此基础上，提出了一种基于Bézier曲线的翼型表面线近似算法。

为了避免对翼型面线在每个时间步进行Bézier曲线重构，可以采用以下方法。在计算之初，需要建立动叶翼型在位置处的Bézier曲线及其导数 ．然后是水平集函数在点可以计算为 为转子按角度旋转逆时针方向。在这里是否为位置的转子表面线建立了水平集函数 通过使用Bézier曲线和是点的像吗后顺时针旋转上角 ．因此，水平集函数可以在任何时候重新计算，而无需为每个新的转子位置重建Bézier曲线。

细胞面分数比和介绍了。它们在区间内取值代表了东面和北面的流体部分 ，分别。一维线性插值 上段和 上段用于计算网格面比:

根据LS-STAG方法的概念，在基底网格上取样法向雷诺应力分量(类似于压力离散)，在基底网格单元的右上角取样剪切应力分量。

用半隐式欧拉格式对由空间连续性和动量方程采样得到的微分代数系统进行时间积分[8］．预测器步进导致下一时间层速度预测的亥姆霍兹方程的离散模拟。校正步骤导致了压力校正的泊松方程的离散模拟。通过计算流动变量，得到翼型运动的动力学方程(4)应该被解决。

转子角速度为 ．所以不同的类比(4)可以下列形式书写: 在-第n个时间步长可以计算如下: 在这里中心的坐标是细胞,是机翼旋转的点的坐标，和和阻力和升力是否作用在切割单元的实心部分在th时间步骤: 在这里 ， ， ， ,是切应力的平方，切应力取决于切槽[8］．

转子下一个时间步长的角速度值由(6）.用该值可以定义新的转子位置和浸入边界速度: 这种方法允许重建水平集函数和所有矩阵以及下一个时间步骤所需的源项。

在研究一些复杂的物理现象时，对所考虑的施工行为进行初步的定性估计是非常重要的。它们允许寻找网格细化的域，预测结构的动态，并估计CFL数。但气动载荷的较大波动会发生在粗网格上。所以作用在翼型上的扭矩值应该被低通滤波器过滤(图)7）.为此，二阶巴特沃斯低通滤波器[16)设计。

有必要解释选择这个过滤器的原因。从计算的角度来看，具有无限脉冲响应的滤波器比具有有限脉冲响应的滤波器更便宜。巴特沃斯滤波器的频率响应在通频带内是最大平坦的(即没有波纹)，在阻频带内滚向零[16］．比切比雪夫型I / II型滤波器或椭圆滤波器,巴特沃斯滤波器具有慢转出,因此需要一个高阶来实现一个特定的阻带规范,但巴特沃斯滤波器具有线性相位响应在通频带比切比雪夫型I / II型和椭圆滤波器可以实现(16］．

巴特沃斯二阶滤波器的传递函数为: 在这里 ， ， ,是通频带中的失真水平。实践表明，高阶滤波器会导致滤波信号出现数值不稳定性。对于这个函数, ．因此，归一化巴特沃斯二阶滤波器对应的传递函数为: 为了控制滤波器截止频率，需要使用以下传递函数:

转矩的采样频率为 在数值模拟(为时间离散化步骤)。必须使滤波器截止频率等于 Hz，则截止频率上的抑制电平等于 dB。因此，下面的条件施加在所需的滤波器频率响应函数上(图8)必须满足所需的滤波器传递函数 ：

解决方案(13)会导致以下结果: 因此，所期望的滤波器传递函数为:

为了得到相应的数字滤波系数，需要使用双线性变换[16]: 在这里 ．众所周知 对于具有无限脉冲响应的二阶低通滤波器。因此，由式(16)和(17)，则所设计的数字滤波器系数为: 所以，过滤后的扭矩作用在机翼上从当时的气流 是

4.数值实验

我们考虑两个叶片的萨伏尼乌斯旋翼的自旋(图)5) ．旋转半径从旋转轴到桶的外缘测量。turbine-swept区域等于 ．雷诺数是基于转子直径。仿真中使用了以下参数:

需要注意的是，上述水平集重构算法可以很容易地应用于其他形状的转子，即Darrieus转子。

无量纲平均转子角速度为: 平均在16个无量纲时间单位进行。这段时间足以使转子以平均转子角速度最小的值完成全转 ．在选定的参数值下，叶尖速比等于 转矩系数值是通过非平稳依赖项的平均得到的随着时间的推移, 同样的，功率系数是通过非平稳依赖项的平均得到的随着时间的推移,

在非均匀网格上进行了大量的计算 ．采用空间离散步骤的均匀网格块体 是在转子附近使用的。时间离散步长为 ．计算是在基于Intel C610平台的服务器上进行的，使用Intel Xeon E5-1620 V3 4核处理器(3.5 GHz)，支持超线程(8个逻辑核)。服务器配备16gb ECC DDR4-2133 RAM和2块硬盘(2tb)，组成RAID1盘卷。该服务器为Windows server 2012 R2操作系统。模拟一个无量纲时间单位，大约需要24小时。

功率系数的计算值与实验数据对比[17，18]在[18］．从图中可以看出9，计算结果与实验结果吻合较好。一个模拟非平稳依赖关系的例子并计算出相应的功率系数值在 和 如图所示10．

5.结论

我们有以下结论(我)利用改进的LS-STAG水平集/切割单元浸入边界法，开发了一套用于风力机转子自旋数值模拟的软件。该软件包可用于其他形状的转子(萨沃尼乌斯转子，达瑞乌斯转子，等)。(2)描述了水平集函数构造和重算的算法。(3)设计了一种二阶巴特沃斯低通滤波器用于粗网格仿真时的气动力矩滤波。(iv)本文对两叶片自旋萨伏尼乌斯转子的流动进行了数值模拟 ．(v)计算结果与实验数据定性一致。

的利益冲突

作者声明本文的发表不存在利益冲突。

致谢

这项研究得到了俄罗斯教育和科学部的部分资助(项目编号:020920917)。9.2422.2017/PP)和俄罗斯联邦总统青年博士资助项目(项目编号:92422.2017 /PP)。MK-7431.2016.8)。

参考文献

1. 中国。Krogstad和J. a . Lund，“模型涡轮机性能的实验和数值研究”，风能，第15卷，第5期。3，页443-457,2012。视图:出版商的网站|谷歌学术搜索
2. N. N. Sorensen, J. A. Michelsen和S. Schreck，“Navier-stokes预测的NREL第VI阶段转子在NASA Ames 80 × 120风洞，”在2002年ASME风能研讨会论文集，第94-105页，美国内华达州里诺，2002年1月。视图:谷歌学术搜索
3. K. Rogowski和R. Maroński，“萨伏尼乌斯转子的CFD计算”，理论与应用力学学报，第53卷，第53期1, pp. 37-45, 2015。视图:出版商的网站|谷歌学术搜索
4. R. Gupta和K. K. Sharma，“使用计算流体动力学(CFD)的Darrieus-Savonius组合转子的流动物理”，Int。研究j·英格。科学、技术和创新，第1卷，第1 - 13页，2012。视图:谷歌学术搜索
5. 刘志强，“风力机转子的CFD模拟与实验比较”，载可再生能源研究会议论文集，rrc 2014，第111-116页，2014年6月。视图:出版商的网站|谷歌学术搜索
6. M. H. Mohamed, a . M. Ali, a . a . Hafiz，“作为低速风能转换器的h -转子Darrieus涡轮机的CFD分析”，《工程科学与技术》，国际期刊第18卷第2期1, pp. 1 - 13, 2015。视图:出版商的网站|谷歌学术搜索
7. R.米塔尔和G. Iaccarino的《浸入式边界法》流体力学年度回顾。37卷。，第37卷为基础。启流体机械。，第239-261页，年度回顾，帕洛阿尔托，加州，2005。视图:出版商的网站|谷歌学术搜索|MathSciNet
8. Y. Cheny和O. Botella，“LS-STAG方法:计算具有良好守恒特性的复杂运动几何中的不可压缩粘性流的一种新的浸入边界/水平集方法”，计算物理学杂志第229卷，第229期4, pp. 1043-1076, 2010。视图:出版商的网站|谷歌学术搜索|MathSciNet
9. S. Osher和R. Fedkiw，水平集方法与动态隐式曲面， vol. 153，施普林格，纽约，美国，纽约，2003。视图:出版商的网站|MathSciNet
10. H. a . van der Vorst，“BI-CGSTAB:用于求解非对称线性系统的BI-CG的快速和平滑收敛变体”，科学与统计计算杂志，第13卷，第2期2，第631-644页，1992。视图:出版商的网站|谷歌学术搜索|MathSciNet
11. p·韦塞尔,多重网格方法简介， John Wiley & Sons，奇切斯特，英国，1992。视图:MathSciNet
12. J. van Kan, C. Vuik, P. Wesseling，“基于时间的二维和三维不可压缩流动的快速压力计算”，数值线性代数及其应用，第7卷，第5期6，页429-447,2000。视图:出版商的网站|谷歌学术搜索|MathSciNet
13. I. Marchevsky和V. Puzikova，“线性系统求解的OpenFOAM迭代方法效率分析”，RAS系统编程研究所论文集， vol. 24, pp. 71-86, 2013。视图:出版商的网站|谷歌学术搜索
14. V. V. Puzikova和I. K. Marchevsky，“LS-STAG浸入边界法在耦合气动弹性问题数值模拟中的应用”，发表于第十一届世界计算力学大会论文集，WCCM 2014，第五届欧洲计算力学大会论文集，ECCM 2014，第六届欧洲计算流体力学大会论文集，ECFD 2014，第1995-2006页，特别是2014年7月。视图:谷歌学术搜索
15. V. V. Puzikova，“当使用LS-STAG方法建模经过它的气流时，任意拓扑翼型的水平集函数的构建”，在莫斯科国立技术大学鲍曼先驱报。自然科学，页163-173,2012(俄语)。视图:谷歌学术搜索
16. r·g·里昂理解数字信号处理培生教育，2010。
17. R. E. Sheldahl, B. F. Blackwell, L. V. Feltz，“两桶和三桶萨伏尼乌斯转子的风洞性能数据”，J能量，第2卷，第2期3，页160-164,1978。视图:出版商的网站|谷歌学术搜索
18. p . n . Shankar国家航空实验室， AE-TM-3-76，国家航空实验室，班加罗尔，印度，1976。

版权

版权所有©2017 Ilia K. Marchevsky和Valeria V. Puzikova。这是一篇发布在知识共享署名许可协议，允许在任何媒介上不受限制地使用、传播和复制，但必须正确引用原作。