CN103902810A - 一种用于风力机的涡面/涡环混合自由涡尾迹方法 - Google Patents

Abstract

本发明公开一种用于风力机的涡面/涡环混合自由涡尾迹方法，包括以下几个步骤：步骤一：给定风轮和叶片基本参数和初始入流条件；步骤二：计算初始尾迹形状；步骤三：根据叶素理论计算叶片各叶素的入流角和迎角，确定叶素的附着环量、尾随涡强度和脱体涡强度；步骤四：计算近尾迹节点及涡环控制点的诱导速度，并确定新的尾迹形状，包括近尾迹涡面和远尾迹涡环；步骤五：由新尾迹形状计算叶素控制点的诱导速度，计算方法同步骤四；步骤六：计算几何残差，当几何残差小于设定的阈值时，尾迹形状收敛，进行步骤七；否则返回步骤三；步骤七：计算叶片的气动载荷及气动性能。该涡尾迹方法能缩短计算时间。

Description

一种用于风力机的涡面/涡环混合自由涡尾迹方法

技术领域

本发明涉及一种用于风力机的涡面/涡环混合自由涡尾迹方法，属于风力机空气动力学应用领域。

背景技术

风力机气动特性计算方法主要有三种：叶素动量理论、涡尾迹方法、计算流体动力学方法。综合考虑方法的准确性和计算成本，涡尾迹方法更具优势，而且其本质上具有旋涡特性，能更准确地计及尾流场尾涡间相互诱导作用，是模拟风力机气动特性较为灵活的数值工具。

自由涡尾迹方法不需要流场的先验数据，直接通过解涡线控制方程得到尾迹结构，有较强的理论根据，是涡尾迹方法中应用最广泛的。自由涡尾迹方法中，涡线从叶片尾缘拖出流向下游，涡线的控制节点随当地流速自由移动，通过求解涡线控制方程得到最终的尾迹几何形状。目前，许多学者采用自由涡尾迹建模计算风力机气动特性，尾迹形状的描述均是采用上述从叶片尾缘拖出自由移动至远场的方式。这些模型中尾迹节点的当地流速需要叠加所有涡线对该节点的诱导速度，故计算量与节点数目是平方关系，虽然计算时间远小于计算流体动力学方法，但是在对效率要求很高的工程应用中仍不被看重，因此建立一种准确快捷的自由涡尾迹方法意义重大。

发明内容

本发明的目的，在于提供一种用于风力机的涡面/涡环混合自由涡尾迹方法，可以缩短计算时间。

为了达成上述目的，本发明的解决方案是：

一种用于风力机的涡面/涡环混合自由涡尾迹方法，包括以下几个步骤：

步骤一：采集风力机中风轮和叶片的以下参数：风轮半径R、风轮转速Ω、来流速度

叶素段数N_E、方位角一圈离散段数N_T和翼型数据；

步骤二：计算初始尾迹形状，包括近尾迹涡面中节点在直角坐标系中三个方向的坐标、远尾迹的叶根涡环中节点在直角坐标系三个方向的坐标和远尾迹的叶尖涡环中节点在直角坐标系三个方向的坐标；

步骤三：根据叶素理论计算叶片各叶素的入流角和迎角，确定叶素的附着环量、尾随涡强度和脱体涡强度；

步骤四：计算近尾迹节点及涡环控制点的诱导速度，并确定新的尾迹形状，包括近尾迹涡面和远尾迹涡环；

步骤五：根据步骤四确定的新的尾迹形状，计算叶素控制点的诱导速度；

步骤六：计算前述新的尾迹形状的几何残差，该几何残差用新的尾迹形状与初始尾迹形状中各节点坐标的均方差表示，当几何残差小于设定的阈值时，该新的尾迹形状收敛，进行步骤七；否则返回步骤三；

步骤七：计算叶片的气动载荷及气动性能。

进一步的，所述步骤二中，近尾迹涡面中节点的计算方法为：将叶片分成N_E段叶素，则有h个控制点和g个边界点，各边界点距风轮中心距离表示为：

${\left(r_{\mathrm{bp}}\right)}_g=\left\{\begin{array}{cc}{R}_t& g=1\\ R_t+\frac{\frac{2}{\mathrm{\π}}\arccos (1-\frac{g-1}{N_E})-0.1}{0.9}(R-R_t)& g=\mathrm{2,3},...,N_E+1\end{array}\right.$

式中，R_t为叶根与风轮中心的距离，控制点距风轮中心距离为：

${\left(r_{\mathrm{cp}}\right)}_h=\frac{1}{2}[{\left(r_{\mathrm{bp}}\right)}_h+{\left(r_{\mathrm{bp}}\right)}_{h+1}]h=\mathrm{1,2},...,N_E$

叶片在方位角为ψ_j时，第g个叶素边界拖出的尾迹中寿命角为ζ_k的节点在直角坐标系中三个方向的坐标为：

$\left\{\begin{array}{c}{x}_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_n(g,j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_k}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.$

式中，V₀为来流速度

的数值，Δζ为时间差Δt对应的寿命角步长，k表示尾迹离散后的节点序号，从叶片尾缘拖出的第一个涡线节点开始计数，k的最大值等于尾迹总的寿命角除以步长，步长等于360度除以N_T；；远尾迹的叶根涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{root}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.$

式中，ζ_control为涡环控制点的寿命角；叶尖涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{tip}}(j,k)=R\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{tip}}(j,k)=R\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{tip}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right..$

进一步的，所述步骤三中，入流角的计算方法是：

${\mathrm{\φ}}_i=\arctan \frac{V_0(1-a_i)}{\mathrm{\Ω}{\left(r_{\mathrm{cp}}\right)}_i(1+a_i^{\′})}$

式中，a_i为第i个叶素处的轴向诱导因子，a′_i为第i个叶素处的切向诱导因子，假设初始诱导因子均为零，第i个叶素的迎角的计算方法是：

α_i=φ_i-θ_i

式中，θ_i为第i个叶素的桨距角，确定第i个叶素的附着环量为：

${\left({\mathrm{\Γ}}_b\right)}_i=\frac{1}{2}{W}_i{C}_L{c}_i$

式中，W_i为第i个叶素的合速度，c_i为第i个叶素的弦长，C_L为翼型升力系数；尾随涡强度定义为相邻叶素附着环量之差，因此，确定第g个叶素边界拖出的尾随涡强度为：

${\left({\mathrm{\Γ}}_t\right)}_{g,j}\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{1,j}& g=1\\ {\left({\mathrm{\Γ}}_b\right)}_{g,j}-{\left({\mathrm{\Γ}}_b\right)}_{g-1,j}& g=\mathrm{2,3},...,N_E\\ {\left({\mathrm{\Γ}}_b\right)}_{N_E,j}& g=N_E+1\end{array}\right.$

式中，j=1,2,...,N_T，(Γ_t)_g,j中的下标j表示对应第j个方位角；

脱体涡强度定义为相邻方位角上叶素附着环量之差，因此，确定第j个方位角下的脱体涡强度为：

${\left({\mathrm{\Γ}}_s\right)}_{i,j}=\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{i,1}-{\left({\mathrm{\Γ}}_b\right)}_{i,N_T}& j=1\\ {\left({\mathrm{\Γ}}_b\right)}_{i,j}-{\left({\mathrm{\Γ}}_b\right)}_{i,j-1}& j=\mathrm{2,3},...,N_T\end{array}\right.$

式中，i=1,2,...,N_E。

进一步的，所述步骤四中每段离散的涡线段对空间某点的诱导速度可以由Biot-Sarvart定律求得：

${\stackrel{\→}{V}}_{\mathrm{ind}}=\frac{\mathrm{\Γ}}{4\mathrm{\π}}{\∫}_A^B\frac{d\stackrel{\→}{l}\×\stackrel{\→}{r}}{r^3}$

式中，A、B分别表示涡线段的两个端点，

是涡线段的向量，且方向由A指向B；

是空间点的向量；环量Γ是指每个涡线段的涡元强度：在计算叶片附着涡对空间节点的诱导速度时采用各个叶素的附着环量Γ_b代入；在计算尾随涡对空间节点的诱导速度时采用各叶素边界拖出尾随涡的强度Γ_t代入；在计算脱体涡对空间节点的诱导速度时采用各叶素尾缘拖出的脱体涡的强度Γ_s代入；将空间所有涡线段对空间某点的诱导速度叠加，得到该点的诱导速度；

近尾迹涡面中涡线节点满足涡线控制方程：

$\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ψ}}+\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ξ}}=\frac{1}{\mathrm{\Ω}}[{\stackrel{\→}{V}}_0+{\stackrel{\→}{V}}_{\mathrm{ind}}\left(\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})\right)]$

式中，为涡线节点的位置矢量，ψ为方位角，ζ为尾迹寿命角，Ω为风轮旋转速度，

为来流速度，

为流场中所有涡线和涡环对该节点的诱导速度总和；

远尾迹涡环控制点的周向位置与所替代的实际涡线中点一致，以叶片旋转一圈为周期，第一个涡环的径向位置和轴向位置由近尾迹末端节点随当地流速移动半个周期确定，第二个涡环的径向位置和轴向位置由第一个涡环控制点随当地流速移动一个周期确定，依次类推，即

$\left\{\begin{array}{c}{r}_r^1=r_r^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^2=r_r^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^3=r_r^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_r^n=r_r^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.$

$\left\{\begin{array}{c}{r}_z^1=r_z^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^2=r_z^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^3=r_z^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_z^n=r_z^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.$

式中，

和

分别为近尾迹末端节点的径向位置和轴向位置，

和

分别为由前向后第n个涡环控制点的径向位置和轴向位置，其中，n为自然数，T为时间周期，Δt为近尾迹相邻节点的时间差，和分别为近尾迹末端节点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，

和

分别为由前向后第n-1个涡环控制点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，其中，n为自然数，当地流速为入流速度与当地诱导速度之和。

进一步的，所述步骤六中的几何残差用新的尾迹形状与初始尾迹形状中各节点坐标的均方差表示：

$R_{\mathrm{MS}}=\sqrt{\frac{\underset{j=1}{\overset{j_{\max }}{\mathrm{\Σ}}}\underset{k=1}{\overset{k_{\max }}{\mathrm{\Σ}}}{(r_{j,k}^{\mathrm{new}}-r_{j,k}^{\mathrm{old}})}^2}{j_{\max }{k}_{\max }}}$

式中，

表示执行步骤四之后的新尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；

表示执行步骤四之前的初始尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；j_max=N_T，k_max等于k的最大值，当几何残差小于设定的阈值时，尾迹形状收敛，进行步骤七；否则返回步骤三。

进一步的，根据步骤五中叶素控制点的诱导速度得到轴向诱导因子和切向诱导因子，根据步骤三中方法计算各叶素的迎角，再通过迎角和翼型数据得到各叶素的气动载荷，再通过对各叶素的气动载荷的积分得到整个叶片的气动性能。

采用上述方案后，本发明所提供的一种用于风力机的涡面/涡环混合自由涡面/涡环混合自由涡尾迹方法，将远尾迹涡线简化成涡环进行尾迹迭代，涡环形状和位置由涡环的控制点确定，大大减少了尾迹迭代时间，以减少自由涡尾迹方法的计算时间。

附图说明

图1为本发明的涡环控制点示意图。

图2为本发明的涡面/涡环混合尾迹三维图。

图3为本发明的涡面/涡环混合尾迹轴向视图。

图4为本发明的涡面/涡环混合尾迹侧向视图。

图5为本发明计算时间与未用本发明计算时间比较图。

图6为本发明的低速轴扭矩与实验值和未用本发明计算的比较图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本发明提出一种用于风力机的涡面/涡环混合自由涡尾迹方法，尤其用于风力机气动性能计算的涡面/涡环混合自由涡面/涡环混合自由涡尾迹方法。风力机的叶片旋转时，尾缘会有一系列尾涡产生，尾涡从尾缘拖出后一般会自由卷起成集中涡，包括叶尖涡和叶根涡。涡线会对风轮平面产生诱导，改变叶片的入流特性。在风轮下游一定距离的远尾迹区，涡线对风轮平面的诱导作用越来越弱。基于这些特征，将远尾迹区的涡线简化叶尖涡环和叶根涡环进行自由涡尾迹方法建模，整个尾迹由近尾迹涡面和远尾迹涡环组成。涡环与所替代的实际涡线中点相交，该交点作为涡环的控制点，涡环的半径及位置由控制点确定。

参照图1，为本发明实施例的涡环控制点示意图，远尾迹涡线2由涡环3代替，涡环3与所替代的远尾迹涡线2中点相交，该交点作为涡环的控制点1，涡环的半径及位置由控制点1确定。

图2至图4为利用本发明实施例风力机风轮尾迹结果，分别为混合尾迹三维图，轴向视图和侧向视图。风轮尾迹由近尾迹涡面4和远尾迹涡环组成，远尾迹涡环分为叶尖涡环5和叶根涡环6，近尾迹涡面4遵循涡线控制方程，远尾迹涡环遵循涡环控制方程。综合考虑模型计算效率和计算准确度，上述近尾迹涡面的尾迹寿命角ζ范围为120度至150度。以一个叶尖涡环和叶根涡环为一对涡环，涡环的对数等于风轮后两倍直径尾迹区内实际尾迹的圈数减1。该实施例中，所述近尾迹涡面4的尾迹寿命角ζ范围取120度。该实施例风力机风轮直径10.058米，风轮转速72转/分，计算风速10米/秒，涡环的对数为3。

近尾迹涡面中涡线节点满足涡线控制方程：

$\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ψ}}+\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ξ}}=\frac{1}{\mathrm{\Ω}}[{\stackrel{\→}{V}}_0+{\stackrel{\→}{V}}_{\mathrm{ind}}\left(\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})\right)]---\left(1\right)$

式中，

为涡线节点7的位置矢量，ψ为方位角，ζ为尾迹寿命角，Ω为风轮旋转速度，为来流速度，

为流场中所有涡线和涡环对该节点的诱导速度总和。

远尾迹涡环控制点1的周向位置与所替代的实际涡线中点一致。以叶片旋转一圈为周期，第一个涡环的径向位置和轴向位置由近尾迹末端节点8随当地流速移动半个周期确定，第二个涡环的径向位置和轴向位置由第一个涡环控制点随当地流速移动一个周期确定，依次类推，即

$\left\{\begin{array}{c}{r}_r^1=r_r^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^2=r_r^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^3=r_r^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_r^n=r_r^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.---\left(2\right)$

$\left\{\begin{array}{c}{r}_z^1=r_z^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^2=r_z^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^3=r_z^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_z^n=r_z^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.---\left(3\right)$

式中，

和

分别为近尾迹末端节点的径向位置和轴向位置，

和

分别为由前向后第n个涡环控制点的径向位置和轴向位置，其中，n为自然数，T为时间周期，Δt为近尾迹相邻节点的时间差，

和

分别为近尾迹末端节点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，

和

分别为由前向后第n-1个涡环控制点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，其中，n为自然数，当地流速为入流速度与当地诱导速度之和。

参照图5，为本发明实施例计算时间与未用本发明计算时间比较图。可以看出，本实施例在未用本发明时，尾迹迭代时间为531.3秒，利用本发明后迭代时间为7.1秒，缩短了75倍，计算效率得到大幅提高。

参照图6，为本发明实施例的低速轴扭矩与实验值和未用本发明计算的比较图可以看出，利用本发明后的计算结果与未用本发明计算结果相差不大，且与实验值比较接近，说明本发明在提高模型计算效率的同时，又能够保证模型计算的准确度。

本发明的一种用于风力机的涡面/涡环混合自由涡尾迹方法，包括以下几个步骤：

步骤一：给定风轮和叶片基本参数和初始入流条件，如：风轮半径R、风轮转速Ω、来流速度

、空气密度ρ、叶片数目N_b、叶素段数N_E、方位角一圈离散段数NT和翼型数据；

步骤二：计算初始尾迹形状，近尾迹为等螺距涡面，涡面中节点计算方法为：将叶片分成N_E段叶素，则有h个控制点和g个边界点，各边界点距风轮中心距离表示为：

${\left(r_{\mathrm{bp}}\right)}_g=\left\{\begin{array}{cc}{R}_t& g=1\\ R_t+\frac{\frac{2}{\mathrm{\π}}\arccos (1-\frac{g-1}{N_E})-0.1}{0.9}(R-R_t)& g=\mathrm{2,3},...,N_E+1\end{array}\right.---\left(4\right)$

式中，R_t为叶根与风轮中心的距离，控制点距风轮中心距离为：

${\left(r_{\mathrm{cp}}\right)}_h=\frac{1}{2}[{\left(r_{\mathrm{bp}}\right)}_h+{\left(r_{\mathrm{bp}}\right)}_{h+1}]h=\mathrm{1,2},...,N_E---\left(5\right)$

叶片在方位角为ψ_j时，第g个叶素边界拖出的尾迹中寿命角为ζ_k的节点在直角坐标系中三个方向的坐标为：

$\left\{\begin{array}{c}{x}_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_n(g,j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_k}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.---\left(6\right)$

式中，V₀为来流速度

的数值，Δζ为时间差Δt对应的寿命角步长，k表示尾迹离散后的节点序号，从叶片尾缘拖出的第一个涡线节点开始计数，k的最大值等于尾迹总的寿命角除以步长，步长等于360度除以N_T；远尾迹的叶根涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{root}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.---\left(7\right)$

式中，ζ_control为涡环控制点的寿命角；叶尖涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{tip}}(j,k)=R\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{tip}}(j,k)=R\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{tip}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.---\left(8\right)$

步骤三：根据叶素理论计算叶片各叶素的入流角，入流角为：

式中，a_i为第_i个叶素处的轴向诱导因子，a′_i为第i个叶素处的切向诱导因子，假设初始诱导因子均为零。第i个叶素的迎角为：

α_i=φ_i-θ_i   （10）

式中，θ_i为第i个叶素的桨距角。根据迎角和翼型数据得到叶素的载荷，确定第i个叶素的附着环量为：

${\left({\mathrm{\Γ}}_b\right)}_i=\frac{1}{2}{W}_i{C}_L{c}_i---\left(11\right)$

式中，W_i为第i个叶素的合速度，c_i为第i个叶素的弦长，C_L为翼型升力系数；尾随涡强度定义为相邻叶素附着环量之差，因此，确定第g个叶素边界拖出的尾随涡强度为：

${\left({\mathrm{\Γ}}_t\right)}_{g,j}\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{1,j}& g=1\\ {\left({\mathrm{\Γ}}_b\right)}_{g,j}-{\left({\mathrm{\Γ}}_b\right)}_{g-1,j}& g=\mathrm{2,3},...,N_E\\ {\left({\mathrm{\Γ}}_b\right)}_{N_E,j}& g=N_E+1\end{array}\right.---\left(12\right)$

式中，j=1,2,...,N_T，(Γ_t)_g,j中的下标j表示对应第j个方位角；

脱体涡强度定义为相邻方位角上叶素附着环量之差，因此，确定第j个方位角下的脱体涡强度为：

${\left({\mathrm{\Γ}}_s\right)}_{i,j}=\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{i,1}-{\left({\mathrm{\Γ}}_b\right)}_{i,N_T}& j=1\\ {\left({\mathrm{\Γ}}_b\right)}_{i,j}-{\left({\mathrm{\Γ}}_b\right)}_{i,j-1}& j=\mathrm{2,3},...,N_T\end{array}\right.---\left(13\right)$

式中，i=1,2,...,N_E；

步骤四：计算近尾迹节点及涡环控制点的诱导速度，每段离散的涡线对空间某点的诱导速度可以由Biot-Sarvart定律求得：

${\stackrel{\→}{V}}_{\mathrm{ind}}=\frac{\mathrm{\Γ}}{4\mathrm{\π}}{\∫}_A^B\frac{d\stackrel{\→}{l}\×\stackrel{\→}{r}}{r^3}---\left(14\right)$

式中，A、B分别表示涡线段的两个端点，

是涡线段的向量，且方向由A指向B；

是空间点的向量；环量Γ是指每个涡线段的涡元强度：在计算叶片附着涡对空间节点的诱导速度时采用各个叶素的附着环量Γ_b代入；在计算尾随涡对空间节点的诱导速度时采用各叶素边界拖出尾随涡的强度Γ_t代入；在计算脱体涡对空间节点的诱导速度时采用各叶素尾缘拖出的脱体涡的强度Γ_s代入；将空间所有涡线段对空间某点的诱导速度叠加，得到该点的诱导速度；根据尾迹控制方程（1）、（2）和（3）确定新尾迹形状，包括近尾迹涡面和远尾迹涡环；

步骤五：由新尾迹形状计算叶素控制点的诱导速度，计算方法同步骤四；

步骤六：计算几何残差，几何残差用新的尾迹形状与初始尾迹形状中各节点坐标的均方差表示：

$R_{\mathrm{MS}}=\sqrt{\frac{\underset{j=1}{\overset{j_{\max }}{\mathrm{\Σ}}}\underset{k=1}{\overset{k_{\max }}{\mathrm{\Σ}}}{(r_{j,k}^{\mathrm{new}}-r_{j,k}^{\mathrm{old}})}^2}{j_{\max }{k}_{\max }}}---\left(15\right)$

式中，

表示执行步骤四之后的新尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；表示执行步骤四之前的初始尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；j_max=N_T，k_max等于k的最大值（见步骤二）。当几何残差小于设定的阈值（在本实施例中，该阈值为1×10^-4）时，尾迹形状收敛，进行步骤七；否则返回步骤三；

步骤七：计算叶片的气动载荷及气动性能。根据步骤五中叶素控制点的诱导速度得到轴向诱导因子a和切向诱导因子a′，根据步骤三中方法计算各叶素的迎角，再通过迎角和翼型数据得到各叶素的气动载荷，再通过对各叶素的气动载荷的积分得到整个叶片的气动性能。

本发明所提供的一种用于风力机的涡面/涡环混合自由涡尾迹方法，将远尾迹涡线简化成涡环进行尾迹迭代，涡环形状和位置由涡环的控制点确定，大大减少了尾迹迭代时间。并且，由于远尾迹涡线对风轮平面的诱导影响已经变弱，远尾迹涡线的简化对风力机气动性能的计算影响不大，因此涡面/涡环混合自由涡尾迹方法在提高计算效率的同时，又能保证计算的准确性。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (6)

1.一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：包括以下几个步骤：

步骤一：采集风力机中风轮和叶片的以下参数：风轮半径R、风轮转速Ω、来流速度

叶素段数N_E、方位角一圈离散段数N_T和翼型数据；

步骤二：计算初始尾迹形状，包括近尾迹涡面中节点在直角坐标系中三个方向的坐标、远尾迹的叶根涡环中节点在直角坐标系三个方向的坐标和远尾迹的叶尖涡环中节点在直角坐标系三个方向的坐标；

步骤三：根据叶素理论计算叶片各叶素的入流角和迎角，确定叶素的附着环量、尾随涡强度和脱体涡强度；

步骤四：计算近尾迹节点及涡环控制点的诱导速度，并确定新的尾迹形状，包括近尾迹涡面和远尾迹涡环；

步骤五：根据步骤四确定的新的尾迹形状，计算叶素控制点的诱导速度；

步骤六：计算前述新的尾迹形状的几何残差，该几何残差用新的尾迹形状与初始尾迹形状中各节点坐标的均方差表示，当几何残差小于设定的阈值时，该新的尾迹形状收敛，进行步骤七；否则返回步骤三；

步骤七：计算叶片的气动载荷及气动性能。

2.如权利要求1所述的一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：所述步骤二中，近尾迹涡面中节点的计算方法为：将叶片分成N_E段叶素，则有h个控制点和g个边界点，各边界点距风轮中心距离表示为：

${\left(r_{\mathrm{bp}}\right)}_g=\left\{\begin{array}{cc}{R}_t& g=1\\ R_t+\frac{\frac{2}{\mathrm{\π}}\arccos (1-\frac{g-1}{N_E})}{0.9}& g=\mathrm{2,3},...,N_E+1\end{array}\right.$

式中，R_t为叶根与风轮中心的距离，控制点距风轮中心距离为：

${\left(r_{\mathrm{cp}}\right)}_h=\frac{1}{2}[{\left(r_{\mathrm{bp}}\right)}_h+{\left(r_{\mathrm{bp}}\right)}_{h+1}]h=\mathrm{1,2},...,N_E$

叶片在方位角为ψ_j时，第g个叶素边界拖出的尾迹中寿命角为ζ_k的节点在直角坐标系中三个方向的坐标为：

$\left\{\begin{array}{c}{x}_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_n(g,j,k)={\left(r_{\mathrm{bp}}\right)}_g\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_n(g,j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_k}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.$

式中，V₀为来流速度

的数值，Δζ为时间差Δt对应的寿命角步长，k表示尾迹离散后的节点序号，从叶片尾缘拖出的第一个涡线节点开始计数，k的最大值等于尾迹总的寿命角除以步长，步长等于360度除以N_T；远尾迹的叶根涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{root}}(j,k)={\left(r_{\mathrm{bp}}\right)}_1\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{root}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.$

式中，ζ_control为涡环控制点的寿命角；叶尖涡环中节点在直角坐标系三个方向的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{tip}}(j,k)=R\cos ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ y_{\mathrm{tip}}(j,k)=R\sin ({\mathrm{\ψ}}_j-{\mathrm{\ζ}}_k)\\ z_{\mathrm{tip}}(j,k)=\mathrm{\Δt}\frac{{\mathrm{\ζ}}_{\mathrm{control}}}{\mathrm{\Δ\ξ}}{V}_0\end{array}\right.$

3.如权利要求2所述的一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：所述步骤三中，入流角的计算方法是：

${\mathrm{\φ}}_i=\arctan \frac{V_0(1-a_i)}{\mathrm{\Ω}{\left(r_{\mathrm{cp}}\right)}_i(1+a_i^{\′})}$

式中，a_i为第i个叶素处的轴向诱导因子，a′_i为第i个叶素处的切向诱导因子，假设初始诱导因子均为零，第i个叶素的迎角的计算方法是：

α_i=φ_i-θ_i

式中，θ_i为第i个叶素的桨距角，确定第i个叶素的附着环量为：

${\left({\mathrm{\Γ}}_b\right)}_i=\frac{1}{2}{W}_i{C}_L{c}_i$

式中，W_i为第i个叶素的合速度，c_i为第i个叶素的弦长，C_L为翼型升力系数；尾随涡强度定义为相邻叶素附着环量之差，因此，确定第g个叶素边界拖出的尾随涡强度为：

${\left({\mathrm{\Γ}}_t\right)}_{g,j}\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{1,j}& g=1\\ {\left({\mathrm{\Γ}}_b\right)}_{g,j}-{\left({\mathrm{\Γ}}_b\right)}_{g-1,j}& g=\mathrm{2,3},...,N_E\\ {\left({\mathrm{\Γ}}_b\right)}_{N_E,j}& g=N_E+1\end{array}\right.$

式中，j=1,2,...,N_T，(Γ_t)_g,j中的下标j表示对应第j个方位角；

脱体涡强度定义为相邻方位角上叶素附着环量之差，因此，确定第j个方位角下的脱体涡强度为：

${\left({\mathrm{\Γ}}_s\right)}_{i,j}=\left\{\begin{array}{cc}{\left({\mathrm{\Γ}}_b\right)}_{i,1}-{\left({\mathrm{\Γ}}_b\right)}_{i,N_T}& j=1\\ {\left({\mathrm{\Γ}}_b\right)}_{i,j}-{\left({\mathrm{\Γ}}_b\right)}_{i,j-1}& j=\mathrm{2,3},...,N_T\end{array}\right.$

式中，i=1,2,...,N_E。

4.如权利要求3所述的一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：所述步骤四中每段离散的涡线段对空间某点的诱导速度可以由Biot-Sarvart定律求得：

${\stackrel{\→}{V}}_{\mathrm{ind}}=\frac{\mathrm{\Γ}}{4\mathrm{\π}}{\∫}_A^B\frac{d\stackrel{\→}{l}\×\stackrel{\→}{r}}{r^3}$

式中，A、B分别表示涡线段的两个端点，

是涡线段的向量，且方向由A指向B；

是空间点的向量；环量Γ是指每个涡线段的涡元强度：在计算叶片附着涡对空间节点的诱导速度时采用各个叶素的附着环量Γ_b代入；在计算尾随涡对空间节点的诱导速度时采用各叶素边界拖出尾随涡的强度Γ_t代入；在计算脱体涡对空间节点的诱导速度时采用各叶素尾缘拖出的脱体涡的强度Γ_s代入；将空间所有涡线段对空间某点的诱导速度叠加，得到该点的诱导速度；

近尾迹涡面中涡线节点满足涡线控制方程：

$\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ψ}}+\frac{\∂\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})}{\∂\mathrm{\ξ}}=\frac{1}{\mathrm{\Ω}}[{\stackrel{\→}{V}}_0+{\stackrel{\→}{V}}_{\mathrm{ind}}\left(\stackrel{\→}{r}(\mathrm{\ψ},\mathrm{\ξ})\right)]$

式中，

为涡线节点的位置矢量，ψ为方位角，ζ为尾迹寿命角，Ω为风轮旋转速度，

为来流速度，

为流场中所有涡线和涡环对该节点的诱导速度总和；

远尾迹涡环控制点的周向位置与所替代的实际涡线中点一致，以叶片旋转一圈为周期，第一个涡环的径向位置和轴向位置由近尾迹末端节点随当地流速移动半个周期确定，第二个涡环的径向位置和轴向位置由第一个涡环控制点随当地流速移动一个周期确定，依次类推，即

$\left\{\begin{array}{c}{r}_r^1=r_r^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^2=r_r^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_r^3=r_r^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_r^n=r_r^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{r,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.$

$\left\{\begin{array}{c}{r}_z^1=r_z^0+\underset{j=1}{\overset{\frac{T}{2\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^0(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^2=r_z^1+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^1(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ r_z^3=r_z^2+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^2(r,\mathrm{\ψ},z)\mathrm{\Δt}\\ \·\\ \·\\ \·\\ r_z^n=r_z^{n-1}+\underset{j=1}{\overset{\frac{T}{\mathrm{\Δt}}}{\mathrm{\Σ}}}{v}_{z,j}^{n-1}(r,\mathrm{\ψ},z)\mathrm{\Δt}\end{array}\right.$

式中，和

分别为近尾迹末端节点的径向位置和轴向位置，和

分别为由前向后第n个涡环控制点的径向位置和轴向位置，其中，n为自然数，T为时间周期，Δt为近尾迹相邻节点的时间差，

和分别为近尾迹末端节点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，

和

分别为由前向后第n-1个涡环控制点在经过j个方位角步长的时间段后当地流速的径向分量和轴向分量，其中，n为自然数，当地流速为入流速度与当地诱导速度之和。

5.如权利要求4所述的一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：所述步骤六中的几何残差用新的尾迹形状与初始尾迹形状中各节点坐标的均方差表示：

$R_{\mathrm{MS}}=\sqrt{\frac{\underset{j=1}{\overset{j_{\max }}{\mathrm{\Σ}}}\underset{k=1}{\overset{k_{\max }}{\mathrm{\Σ}}}{(r_{j,k}^{\mathrm{new}}-r_{j,k}^{\mathrm{old}})}^2}{j_{\max }{k}_{\max }}}$

式中，

表示执行步骤四之后的新尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；

表示执行步骤四之前的初始尾迹的第j个方位角第k个尾迹离散后的节点序号的节点坐标；j_max=N_T，k_max等于k的最大值，当几何残差小于设定的阈值时，尾迹形状收敛，进行步骤七；否则返回步骤三。

6.如权利要求5所述的一种用于风力机的涡面/涡环混合自由涡尾迹方法，其特征在于：根据步骤五中叶素控制点的诱导速度得到轴向诱导因子和切向诱导因子，根据步骤三中方法计算各叶素的迎角，再通过迎角和翼型数据得到各叶素的气动载荷，再通过对各叶素的气动载荷的积分得到整个叶片的气动性能。

Images