CN106407577B - 一种模拟风力机尾流的改进致动面模型建立方法 - Google Patents

Abstract

本发明公开了一种模拟风力机尾流的改进制动面模型建立方法，其特征是，包括如下步骤：1)建立致动模型；2)计算得到沿叶片展向单位长度上的体积力源项，其中，体积力分布采取分段线性分布；3)将体积力源项添加到制动面模型所定义的无厚度平面上；4)将叶片所在平面作为制动面所在平面，根据待确定点与各个叶片的向量积对待确定点的位置进行判断，实现网格点自动识别；5)计算模拟流场。本发明所达到的有益效果：本发明提出的一种模拟风力机尾流的改进制动面模型，可以减少网格数量和计算时间，适应用于风力机尾流场的计算。可对近海风电场微观选址有一定的指导意义，在工程中有很好的应用前景。

Description

一种模拟风力机尾流的改进致动面模型建立方法

技术领域

本发明涉及一种模拟风力机尾流的改进致动面 模型建立方法，属于致动 模型计算风力机尾流技术领域。

背景技术

风电是一种可再生的清洁能源，中国风能资源丰富，风电是一种可再生的清洁能源，中国风能资源丰富，可开发利用的风资源14亿kW，其中陆上6亿 kW、海上8亿kW。近年来，计算流体力学方法在风力机流场及其气动性能研究中得到越来越多的运用，该方法能够准确的描述出风力机及其周围的复杂流场。常规的风机CFD计算，为了保证风力机叶片及尾流区域的计算准确性，需要建立相对复杂的实物模型，先将固体的实体模型用专业的三维建模软件进行建模，建模完成后再对整个流场区域进行网格划分。三维模型的复杂性使得在整个过程中网格的划分成为了一个难点，而且由于在风轮下游远尾流区域内网格还要保持一定的精度，网格的数量大量增加，从而必然会导致计算量的增加。

研究人员通过将BEM理论与常规的计算流体力学方法相结合，开发出了致动模型方法，即先用BEM理论求解风轮叶片的气动力，再将气动力作为体积力源项添加到N-S方程中求解，模拟叶片与流场的作用力。因为模型中没有真实的叶片的固壁边界，所以能够大大减少模型复杂程度以及网格数量，进而节约大量的计算资源。致动模型主要有致动盘模型、致动线模型和致动面 模型。

目前致动模型的主要计算研究方向是风力机叶片的载荷特性，转轮区域流场特性，目前现有致动面 模型还不能准确计算远尾流处流场特性。而风力机远尾流区域的计算对于风场开发、微观选址有着重要的参考意义。

发明内容

为解决现有技术的不足，本发明的目的在于提供一种模拟风力机尾流的改进致动面 模型建立方法，运用改进致动面 模型对风力机的尾流场进行数值模拟计算，并与致动线模型的数值模拟结果进行对比，验证改进致动面 模型方法运用在风力机远尾流区域的计算方面以及在风场选址上的可行性。

为了实现上述目标，本发明采用如下的技术方案：

一种模拟风力机尾流的改进致动面 模型建立方法，其特征是，包括如下步骤：

1)将风力机叶片简化成无厚度的平面，在平面上施加不连续的压力来模拟叶片对气流的作用，形成致动模型；

2)通过BEM理对致动模型计算各截面位置的入流角φ和攻角α大小，利用迭代法计算诱导因子，再根据二维翼型气动数据查表获得每个翼型升力系数和阻力系数；计算得到叶片叶展方向单位长度上的体积力，从而得到沿叶片展向单位长度上的体积力源项，其中，体积力分布采取分段线性分布；

3)将体积力源项添加到致动面 模型所定义的无厚度平面上，并添加一个体积力修正系，得到修正后的体积力；

4)将叶片所在平面作为致动面 所在平面，根据待确定点与各个叶片的向量积对待确定点的位置进行判断，实现网格点自动识别；并将修正后的体积力线性分段分布加载到致动面 上的点中；

5)计算模拟流场：叶片旋转后，在每个时间步长上对致动面 网格重复以上步骤进行识别。

进一步地，所述步骤2)中对体积力分布采取分段线性分布时，以翼型的 1/4c弦长位置为分界点，且保持1/4c弦长位置处翼型的俯仰力矩为零。

进一步地，所述步骤2)中具体计算过程如下：

相对于叶片的空气流相对速度式中，Ω表示风轮转速； V_z表示相对径向速度；r为风轮半径；V_θ表示切向速度；V_rel表示合速度；旋转平面的入流角φ可用下式求得

攻角α＝φ-γ，其中γ为安装角；叶片叶展方向单位长度上的体积力由下式确定式中：C_L＝C_L(α,Re)表示升力系数；

C_D＝C_D(α,Re)表示阻力系数，是以攻角α和雷诺数Re为变量的函数；e_L为升力方向向量；e_D为阻力方向向量；雷诺数Re由c弦长和来流速度确定。

进一步地，所述步骤4)中判断方法为：记O点为叶片旋转中心，P点是风轮旋转平面内任意一点，任意选择一个叶片，Q为该叶片弦线上一点；

若则P点在叶片弦线右边；

若则P点在叶片弦线左边；式中，表示Z轴正方向单位向量，

与来流风速方向相反；

若P点在弦线右边，且则P点在该叶片上；否则，判定P点不在该叶片上，则再依次与其他叶片进行匹配，确定点P在平面上的位置，从而实现计算中对平面上属于致动面 范围内的网格点自动识别。

本发明所达到的有益效果：本发明提出的一种模拟风力机尾流的改进致动面 模型，可以减少网格数量和计算时间，适应用于风力机尾流场的计算。可对近海风电场微观选址有一定的指导意义，在工程中有很好的应用前景。

附图说明

图1是改进致动面 模型流程图；

图2是叶素受力分析；

图3是体积力分布方式；

图4是致动面 识别示意图；

图5是计算域划分示意图，其中(a)为主视图，(b)为左视图；

图6是体积力修正系数对比图；

图7是风轮后2.5D、6D和7.5D处水平线上致动模型计算风速与实验数据对比图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图1-图7所示，本发明所建立的一种模拟风力机尾流的改进致动面 模型，建立的步骤如下：

步骤1) 致动面 模型实质上是将风力机叶片简化成无厚度的平面，在平面上施加不连续的压力来模拟叶片对气流的作用，形成致动面 模型，并对致动面 区域进行识别。

步骤2)叶片上各点压力用二维翼型气动数据，根据BEM理论计算确定。将分布的体积力在流场区域内用三维N-S方程求解，不可压缩N-S方程：式中：为速度；ρ为流体密度；p为压强；μ为粘度系数；f表示作用在旋转叶片上的体积力。

然后，如图2所示，相对于叶片的空气流相对速度可以从上述速度三角形中得出式中：Ω表示风轮转速；V_z表示相对径向速度；r 为风轮半径；V_θ表示切向速度；V_rel表示合速度；旋转平面的入流角φ可用下式求得当地攻角α＝φ-γ，其中γ为安装角。

叶片叶展方向单位长度上的体积力由下式确定：

式中：C_L＝C_L(α,Re)表示升力系数； C_D＝C_D(α,Re)表示阻力系数，是以攻角α和雷诺数Re为变量的函数；e_L为升力方向向量；e_D为阻力方向向量；雷诺数Re由c弦长和来流速度确定。

进一步，通过上一步体积力的计算，得到沿叶片展向单位长度的体积力源项，考虑到叶片弦长对叶片上体积力分布的影响，模型中体积力分布采取分段线性分布方式，如图3所示，以翼型的1/4c弦长位置为分界点，且保持1/4c弦长位置处翼型的俯仰力矩为零。这样的分布方式更符合真实情况下的叶片固壁边界效应，能有效的改进近尾流区域的计算准确度。

步骤3)致动模型计算叶片体积力通过BEM理论，计算各截面位置的入流角φ和攻角α大小，然后利用迭代法计算诱导因子，再根据二维翼型气动数据查表获得每个翼型升力系数和阻力系数最后由上述体积力公式得出。

因此在数值模拟过程中，叶片的三维特性不能够得到充分表现，且迭代后拟合过程也存在一定误差，从而影响计算结果准确度。

综上考虑，添加一个体积力修正系数C_f，那么体积力可表示为： f_2D'＝C_ff_2D。

步骤4)如图4，叶片所在平面即为致动面 所在平面，O点为叶片旋转中心， P点是风轮旋转平面内任意一点，任意选择一个叶片，Q为该叶片弦线上一点。

若则P点在叶片弦线右边；

若则P点在叶片弦线左边。

式中，表示Z轴正方向单位向量(与来流风速方向相反)。

若P点在弦线右边，且(c为叶片弦长)，则P点在该叶片上；若P点不在该叶片上，则再依次与另外两个叶片进行匹配，最终可确定点P在平面上的位置，从而实现计算中对平面上属于致动面 范围内的网格点自动识别，并将步骤2中的体积力f_2D'线性分段分布加载到致动面 上。

步骤5)计算模拟流场，叶片旋转后，随着时间t和角速度ω变化，在每个时间步长上对致动面 网格重复以上步骤进行识别。

下面通过丹麦Nibe A型水平轴风力机实际数据对模型进行验证。风机的轮毂高度为45m，风轮直径为40m，叶片设计采用气动翼型NACA 44系列，额定风速为13m·s^-1。以NibeA型风力机为模型在来流风速为8.5m·s^-1，转速为3.5rad ·s^-1的工况下，用致动面 方法对风机尾流场进行了数值模拟，并将结果与实验测量值进行比较。整个计算域设计为一个规则圆柱体，半径140m，总长1203m，具体划分如图5(其中a图中2、3是计算模型前端流场，4、5、6是计算模型后端流场，b图7是致动面 所在位置，其余是流场分块部分)。

流场中各个分块网格数量如下表：

表1计算域尺寸及网格参数

采用的非结构性网格，在各区域单独加密，同时考虑到网格无关性，具体的网格大小和数量如表1所示，单个致动面 弦长方向上网格数10个左右，保证了模拟结果的准确性和网格的无关系。

进一步，如图6所示，本发明选取1.0，1.1，1.2，1.3四个体积力系数分别作计算与对比，选取风轮后2.5D及6D处截面的计算风速进行对比。比较后可见C_f＝1.2时风速曲线与实验数据更吻合，因此选定的体积力修正系数大小为1.2 时准确可靠。

进一步的，如图7，在流场中轮毂高度截取水平面，分别导出在转轮后侧 2.5D、6D和7.5D三处水平线上速度进行分析，同时将用致动线方法模拟所得结果一并列出，与实验数据进行比较。

从图7中可以看出，三个位置处致动线方法和致动面 方法数值模拟的结果整体接近，与实验数据(EXP)相比，两侧风速基本一致，最小风速都出现在中心位置。±1D范围内尾流速度曲线坡度明显减小。

在2.5D处，致动面 方法与致动线方法的模拟结果和实验数据都很吻合，整体上都呈V型，致动线方法中心位置速度略有波动；在6D和7.5D处，致动线方法计算结果和实验数据偏差较大，轴线位置的实验最小风速为0.7U左右，而致动线模型最小风速只有0.5U。致动线模型截面上风速曲线的趋势和实验数据有较大出入，依然保持V型速度曲线，尤其风轮位置(r/D＝±0.5D)风速明显偏小，误差超过20％，而实验数据中速度变化曲线呈浅弧形(U型)。这说明在致动线模型的数值模拟中，远尾流区域动能耗散与恢复和真实情况误差较大，而致动面 方法的模拟结果就相对准确很多，曲线走势基本一致，尤其是在远尾流区域依然保持了较高的吻合度。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (3)

1.一种模拟风力机尾流的改进致动面模型建立方法，其特征是，包括如下步骤：

1)将风力机叶片简化成无厚度的平面，在平面上施加不连续的压力来模拟叶片对气流的作用，形成致动模型；

2)通过BEM理论对致动模型计算各截面位置的入流角φ和攻角α大小，利用迭代法计算诱导因子，再根据二维翼型气动数据查表获得每个翼型升力系数和阻力系数；计算得到叶片叶展方向单位长度上的体积力，从而得到沿叶片展向单位长度上的体积力源项，其中，体积力分布采取分段线性分布；

3)将体积力源项添加到致动面模型所定义的无厚度平面上，并添加一个体积力修正系数，得到修正后的体积力；

4)将叶片所在平面作为致动面所在平面，根据待确定点与各个叶片的向量积对待确定点的位置进行判断，实现网格点自动识别；并将修正后的体积力线性分段分布加载到致动面上的点中，判断方法为：记O点为叶片旋转中心，P点是风轮旋转平面内任意一点，任意选择一个叶片，Q为该叶片弦线上一点；

若则P点在叶片弦线右边；

若则P点在叶片弦线左边；式中，表示Z轴正方向单位向量，与来流风速方向相反；

若P点在弦线右边，且则P点在该叶片上；否则，判定P点不在该叶片上，则再依次与其他叶片进行匹配，确定点P在平面上的位置，从而实现计算中对平面上属于致动面范围内的网格点自动识别；

5)计算模拟流场：叶片旋转后，在每个时间步长上对致动面网格重复以上步骤进行识别。

2.根据权利要求1所述的一种模拟风力机尾流的改进致动面模型建立方法，其特征是，所述步骤2)中对体积力分布采取分段线性分布时，以翼型的1/4c弦长位置为分界点，且保持1/4c弦长位置处翼型的俯仰力矩为零。

3.根据权利要求1所述的一种模拟风力机尾流的改进致动面模型建立方法，其特征是，所述步骤2)中具体计算过程如下：

相对于叶片的空气流相对速度式中，Ω表示风轮转速；V_z表示相对径向速度；r为风轮半径；V_θ表示切向速度；V_rel表示合速度；旋转平面的入流角φ可用下式求得

攻角α＝φ-γ，其中γ为安装角；叶片叶展方向单位长度上的体积力由下式确定式中：C_L＝C_L(α,Re)表示升力系数；C_D＝C_D(α,Re)表示阻力系数，是以攻角α和雷诺数Re为变量的函数；e_L为升力方向向量；e_D为阻力方向向量；雷诺数Re由c弦长和来流速度确定。