本发明涉及一种基于cfd技术的风力机气动计算模型建模方法，属于风力发电机空气动力分析技术领域。
背景技术：
目前中国风能迅速发展，但同时弃风限电现象时有发生。进一步降低风力发电成本势在必行，其关键在于对风电场的流场分析，图1是风力机流场示意图。通过风电场流场分析，可以得到风力机的性能指标数据和风电场的流场特性，为功率预测、风机性能分析、风机微观选址等多方面风场优化运行提供有效基础理论。因此建立一种有效的风力机气动模型成为风电场流场分析的重要部分。
流场分析始终离不开风力机本身的建模工作，目前使用风力机真实模型所造成的计算量不可能在现有计算机上完成，而多种高精度模拟手段并不能满足工程所需的时间要求。而上世纪70年代丹麦等风电强国的风机模型主要依靠工程经验所设计，已不能满足目前风电行业所需的分析要求。风机建模工作，从维度上可分为二维、三维建模，从气动理论上可分为经验推理、涡方法、叶素动量方法等。总的来说，目前流场分析所需的计算模型在精度上能够提供理想量的有效数据，同时能够满足工程可接受的计算速度。
技术实现要素：
发明目的
本发明的目的在于克服现有工程风机模型精度上的不足，提供一种工程可接受速度并能够提供详细有效数据的风机气动计算模型。本模型建立方法立足于：依靠传统气动理论，并提出多种修正，得出适用于雷诺平均方法的计算模型，并在实际风场数据得到有效验证。
技术方案
一种基于cfd技术的风力机气动计算模型建模方法，其特征在于：所述模型将整个建模过程分为三个阶段依次进行：第一阶段，采用基于bem气动理论的建模方式将风力机简化为计算盘体，保留风力机基本的物理信息；第二阶段，采用基于dns方法的cfd技术将第一阶段的气动模型转化为体积力源项，加入预先设定好的控制方程；第三阶段，加入风力机模型的边界条件，针对网格部分进行有效设计，求解控制方程，并与实际数据收敛，调整相关参数，校正模型准确性。
1、叶片气动模型(第一阶段)
(1)代表风力机叶片作用的体积力基于叶素假设，将叶片的微元段(即叶素)看成二维翼型截面，根据当地的流场信息与已知的翼型性能数据、叶片几何外形信息计算所得。图2就是某叶素的翼型截面图，其中x方向为来流方向。叶片翼型截面的当地速度为：
其中，ω为风轮转动的角速度，r为叶片翼型截面展向位置，vx和vθ分别为叶片翼型截面的轴向速度和切向速度。
叶片翼型截面的当地速度与风轮平面之间的夹角，也就是入流角为：
叶片翼型截面的当地攻角为入流角减去当地桨距角：
α＝φ-γ(3)
其中γ为当地桨距角。
在确定了叶片翼型截面当地速度与当地攻角之后，即可以得到叶片各个截面单位展长升力与阻力：
其中为nb叶片数目，l,d分别为单位展长的升力与阻力，ρ为空气密度，cl＝cl(a,re)和cd＝cd(a,re)分别为升力系数和阻力系数，re为当地弦长c的雷诺数，el和ed为升力和阻力方向的单位向量。
这就是作为体积力动量源项添加到navier-stokes方程中的气动力。
(2)叶尖损失修正
一个真实的风轮，叶片的数目是有限的，叶尖损失也必须被考虑。为了使bem理论更加准确，经典的prandtl叶尖损失修正被引入到该方法中，对动量方程中的诱导速度进行修正。与bem理论相类似，广义致动盘方法也是将风轮简化为一个圆盘，也是无穷个叶片；但是不同的是bem理论是基于一维假设的，是完全不考虑叶片展向流动的，而广义致动盘方法在叶尖处仍是存在一定的叶尖绕流现象，能够通过数值计算模拟出来。然而广义致动盘方法在计算气动力时仍然是基于二维翼型性能数据，但是叶尖附近三维流动现象显著，因此必须对叶尖处的二维翼型性能数据进行修正，以提高广义致动盘方法在叶尖处的计算精度。
本模型提出修正因子：
参数g表达式如下：
g＝exp(-c1(nbωr/u∞-c2))+c3
其中u∞为来流风速，经验参数c1、c2、c3建议值为0.125、21、0.1
于是，修正后的叶片气动模型为：
(3)一维高斯分布
必须采取合适的措施去防止分布体积力周围的空间振荡。由于风轮对流场的作用力并不是单独的存在于一个平面中，而是有一个光滑连续的过渡过程，因此选择一维高斯分布来近似模拟风轮作用力在流场中的分布。它能够保证数值模拟的稳定性并且增加计算的收敛速度。具体的实现方式如图所示，对于风轮平面上的每一个点，体积力在风轮平面法线方向上呈一维高斯分布。
首先将单位展长的升力与阻力在360°方位角上进行平均分布：
再沿风轮平面法线方向进行一维高斯分布的单位体积力为：
其中l为网格单元中心点与风轮平面之间的距离，
με(l)为密度函数，有
其中ε为控制分布密度的分布因子。
2、控制方程的建立(第二阶段)
本模型方案方程采用雷诺平均方法，经过推导，以下给出最终方程形式和相关参数的设定
连续性方程：
其中：uj为xj方向上的速度，j＝1,2,3.
动量方程：
其中：μ为动力黏性系数；sij为平均应变张量；
p为正应力；μ′i为脉动速度；
sd为一阶段所计算的动量源项。
湍流模型方程：
其中：
dk＝β*ρkω
α＝f1α1+(1-f1)α2
β＝f1β1+(1-f1)β2
(y为离最近壁面距离)
封闭常数为：
β*＝0.09,a1＝0.31,α1＝5/9,α2＝0.44,β1＝3/40
β2＝0.0828,σk1＝0.085,σk2＝1,σw1＝0.5,σw2＝0.856
3、求解控制及收敛(第三阶段)
求解边界条件为中性大气层，风速控制在10m/s。如图3整个计算域为长方体，上游距离风轮3d(d为风轮直径)的距离，下游延伸至8d的距离，横截面为边长7d的正方形。为验证模型正确性，关注的重点为轮盘附近的流场信息，因此在轮盘附近划分了一个圆柱形内区进行网格加密(如图4)，保证计算的精度。圆柱形内区上游与下游分别距离轮盘d和2d，横截面圆形直径为3d。收敛判断将使用美国国家可再生能源实验室的phazevi实验数据，部分模拟结果如图5所示。
本专利的优点与积极效果如下：
本发明提供一种基于cfd技术的风力机气动模型建模方法,本发明整体为三维建模，其中扇面采用所基于bem气动理的二维物理建模，轴向采用基于高斯分布的一维数学建模，考虑了叶尖损失效应并做了优化处理；
所使用的cfd技术为雷诺平均方法，建立应用于不可压缩湍流问题的三维纳维斯托克斯方程，其中动量方程使用雷诺平均方程，利用一阶段计算所得的体力源项，湍流模型使用sstk-ω模型；
求解采用基于simple算法的一种解法，在空间上采用有限体积法进行离散，时间离散采用隐式时间积分方法，采用格心格式储存数据，封闭的控制方程进行求解，并与实体数据对比验证模型准确性。
本模型优势在于克服现有工程风机模型精度上的不足，提供一种工程可接受速度并能够提供详细有效数据的风机气动计算模型。本模型建立方法立足于：依靠传统气动理论，并提出多种修正，得出适用于雷诺平均方法的计算模型，并在实际风场数据得到有效验证。
通过本模型，能够有效对风电场流场分析，可以得到风力机的性能指标数据和风电场的流场特性，为功率预测、风机性能分析、风机微观选址等多方面风场优化运行提供了有效基础理论。
附图说明:
图1为计算流程第一模块(基于bem的动量源项)；
图2为计算流程第二模块(基于rans的流场收敛)；
图3为风力机流场示意图；
图4为叶素动量理论下的翼型截面受力示意图；
图5为风力机cfd计算域设计图；
图6为风力机cfd网格划分图；
图7为体积力计算结果图；
具体实施方式
本计算模型的关键在于气动力的准确计算以及合理分布的体积力动量源项加载到计算流场中，整个流程如图1、2所示。基本可以分为两大模块:
在第一模块中，对风力机整体进行数据计算：
第一步，获得风轮各截面的当地攻角。本文采用有限体积法在空问上对控制方程进行空问离散，控制体定义以及流场参数存储使用格心格式，因此将网格单元的中心点坐标作为网格的坐标。通过判断网格单元的坐标可以确定风轮所在的网格单元编号，以及这些网格单元对应的风轮不同的截面。读取这些网格单元的速度信息，并且考虑风轮的旋转速度，获得风轮各截面的当地速度和入流角，当地入流角再减去当地桨距角，即是当地攻角。
第二步，计算风轮各截面的气动力。通过风轮各截面的展向位置，插值得出各白的弦长，根据当地攻角获得相应的升、阻力系数，再加上当地速度以及叶尖修正因了，即可获得当地的升力与阻力。通过坐标转换就可以得到笛卡尔坐标系下风轮各截面的气动力。
第三步，将风轮各截面的气动力通过一维高斯分布合理加载到计算流场中。根据网格单元的坐标，确定气动力的分布位置，由一维高斯分布计算出相应网格单元的气动力，即作为体积力动量源项加载到控制方程中。
通过第一个模块获得加载到控制方程中的体积力动量源项后，进入第二模块。利用图2中所示计算流程来完成整体过程的数值模拟，整体过程基于simple算法，主要划分为四个步骤:
第一步，初始化流场。按照给定的边界条件，对计算域内的单元赋初值；
第二步，计算体积力动量源项。读取流场中的相关信息，以及给定的叶片儿何外形信息、翼型性能数据，按照流程图中的计算流程获得加载到控制方程中的体积力动量源项；
第三步，求解控制方程和湍流模型方程。在空问上采用有限体积法进行离散，时问离散采用隐式时问积分方法，采用格心格式储存数据，对附带体积力动量源项的控制方程和湍流模型方程进行求解；
第四步，判断收敛性。根据残差来判断计算是否收敛，如果未收敛，则返回第二步重新计算体积力动量源项，再求解控制方程，直至收敛。计算收敛后输出相关结果。
所述模型将整个建模过程分为三个阶段依次进行：第一阶段，采用基于bem气动理论的建模方式将风力机简化为计算盘体，保留风力机基本的物理信息；第二阶段，采用基于dns方法的cfd技术将第一阶段的气动模型转化为体积力源项，加入预先设定好的控制方程；第三阶段，加入风力机模型的边界条件，针对网格部分进行有效设计，求解控制方程，并与实际数据收敛，调整相关参数，校正模型准确性。
1、叶片气动模型(第一阶段)
(1)代表风力机叶片作用的体积力基于叶素假设，将叶片的微元段(即叶素)看成二维翼型截面，根据当地的流场信息与已知的翼型性能数据、叶片几何外形信息计算所得。图2就是某叶素的翼型截面图，其中x方向为来流方向。叶片翼型截面的当地速度为：
其中，ω为风轮转动的角速度，r为叶片翼型截面展向位置，vx和vθ分别为叶片翼型截面的轴向速度和切向速度。
叶片翼型截面的当地速度与风轮平面之间的夹角，也就是入流角为：
叶片翼型截面的当地攻角为入流角减去当地桨距角：
α＝φ-γ(3)
其中γ为当地桨距角。
在确定了叶片翼型截面当地速度与当地攻角之后，即可以得到叶片各个截面单位展长升力与阻力：
其中为nb叶片数目，l,d分别为单位展长的升力与阻力，ρ为空气密度，cl＝cl(a,re)和cd＝cd(a,re)分别为升力系数和阻力系数，re为当地弦长c的雷诺数，el和ed为升力和阻力方向的单位向量。
这就是作为体积力动量源项添加到navier-stokes方程中的气动力。
(2)叶尖损失修正
一个真实的风轮，叶片的数目是有限的，叶尖损失也必须被考虑。为了使bem理论更加准确，经典的prandtl叶尖损失修正被引入到该方法中，对动量方程中的诱导速度进行修正。与bem理论相类似，广义致动盘方法也是将风轮简化为一个圆盘，也是无穷个叶片；但是不同的是bem理论是基于一维假设的，是完全不考虑叶片展向流动的，而广义致动盘方法在叶尖处仍是存在一定的叶尖绕流现象，能够通过数值计算模拟出来。然而广义致动盘方法在计算气动力时仍然是基于二维翼型性能数据，但是叶尖附近三维流动现象显著，因此必须对叶尖处的二维翼型性能数据进行修正，以提高广义致动盘方法在叶尖处的计算精度。
本模型提出修正因子：
参数g表达式如下：
g＝exp(-c1(nbωr/u∞-c2))+c3
其中u∞为来流风速，经验参数c1、c2、c3建议值为0.125、21、0.1
于是，修正后的叶片气动模型为：
(3)一维高斯分布
必须采取合适的措施去防止分布体积力周围的空间振荡。由于风轮对流场的作用力并不是单独的存在于一个平面中，而是有一个光滑连续的过渡过程，因此选择一维高斯分布来近似模拟风轮作用力在流场中的分布。它能够保证数值模拟的稳定性并且增加计算的收敛速度。具体的实现方式如图所示，对于风轮平面上的每一个点，体积力在风轮平面法线方向上呈一维高斯分布。
首先将单位展长的升力与阻力在360°方位角上进行平均分布：
再沿风轮平面法线方向进行一维高斯分布的单位体积力为：
其中l为网格单元中心点与风轮平面之间的距离，
με(l)为密度函数，有
其中ε为控制分布密度的分布因子。
2、控制方程的建立(第二阶段)
本模型方案方程采用雷诺平均方法，经过推导，以下给出最终方程形式和相关参数的设定
连续性方程：
其中：uj为xj方向上的速度，j＝1,2,3.
动量方程：
其中：μ为动力黏性系数；sij为平均应变张量；
p为正应力；μ′i为脉动速度；
sd为一阶段所计算的动量源项。
湍流模型方程：
其中：
dk＝β*ρkω
α＝f1α1+(1-f1)α2
β＝f1β1+(1-f1)β2
(y为离最近壁面距离)
封闭常数为：
β*＝0.09,a1＝0.31,α1＝5/9,α2＝0.44,β1＝3/40
β2＝0.0828,σk1＝0.085,σk2＝1,σw1＝0.5,σw2＝0.856
3、求解控制及收敛(第三阶段)
求解边界条件为中性大气层，风速控制在10m/s。如图3整个计算域为长方体，上游距离风轮3d(d为风轮直径)的距离，下游延伸至8d的距离，横截面为边长7d的正方形。为验证模型正确性，关注的重点为轮盘附近的流场信息，因此在轮盘附近划分了一个圆柱形内区进行网格加密(如图4)，保证计算的精度。圆柱形内区上游与下游分别距离轮盘d和2d，横截面圆形直径为3d。收敛判断将使用美国国家可再生能源实验室的phazevi实验数据，部分模拟结果如图5所示。
本发明整体为三维建模，其中扇面采用所基于bem气动理的二维物理建模，轴向采用基于高斯分布的一维数学建模，考虑了叶尖损失效应并做了优化处理；
所使用的cfd技术为雷诺平均方法，建立应用于不可压缩湍流问题的三维纳维斯托克斯方程，其中动量方程使用雷诺平均方程，利用一阶段计算所得的体力源项，湍流模型使用sstk-ω模型；
求解采用基于simple算法的一种解法，在空间上采用有限体积法进行离散，时间离散采用隐式时间积分方法，采用格心格式储存数据，封闭的控制方程进行求解，并与实体数据对比验证模型准确性。