风能作为清洁的可再生能源，是进入21世纪以来世界上发展最快、具有大规模商业应用前景的新能源^[1]. 开展应用风能的相关研究，着力解决风力机运行过程中面临的技术难题，对于风能产业的发展具有重要意义.

现有的关于风力机问题的研究大多是在恒定入流风速条件下运行的. 然而，在实际的运行环境中，入流风速随时间不断变化^[2]，大气边界层中的风速沿高度方向不均匀分布，存在风切变. 研究风力机在动态入流条件下的运行特性，对于指导风力机运行控制、优化风力机功率输出意义重大. 国内外已有一些关于风力机动态入流的初步研究. Knudsen等^[3]在动量-叶素理论的基础上发展了简化的入流模型来描述和评估风力机的动态入流，但是没有对数值模拟的精度进行有效检验. Odgaaard等^[4]在Knudsen等^[3]的入流模型基础上，借助FAST软件仿真，研究动态入流条件下功率与载荷的波动特性. Kamel等^[5]通过Matlab仿真研究风电微网在动态入流条件下的功率响应，将比例积分（proportional-integral，PI）控制方法应用于变桨距控制以实现微网功率随入流风速的平滑波动. 刘雄等^[6]发展了考虑动态入流效应的动态气动载荷分析模型，对风力机在风速变化和快速变桨时的载荷波动特性进行了仿真. 陈严等^[7]开发了适用于柔性叶片翼型的动态入流模型，对比动态入流模型与其他模型的仿真结果，验证了动态入流模型的可行性. 但是，现有的研究大多只是在动量-叶素理论基础上进行理论分析或依托于Matlab、FAST等软件进行风力机运行仿真，利用计算流体力学（computational fluid dynamics, CFD）开展动态入流条件下风力机的运行控制的研究还非常少.

目前很难通过实验测量研究大型商用风力机，数值模拟成为研究大型风力机的重要方法^[8]. 模拟风力机流场时，大涡模拟比雷诺平均方法的精度更高^[9]. 特别是在模拟非定常湍流大气流动时，大涡模拟在流场瞬时值的模拟精度上具有显著的优势^[10]. 由于商用风力机尺寸较大，采用大涡模拟完整解析风力机几何结构需要庞大的网格量，计算成本高昂^[11]. 为此，Sorensen等^[12]提出致动线模型，在流场中引入体积力代替叶片的作用，不直接对风力机叶片几何外形进行建模. 研究表明，基于致动线模型的大涡模拟方法精度较高^[13-14]，可广泛应用于风力机流场的数值模拟. 在流动模拟的基础上，越来越多的研究者使用致动线模型分析风力机载荷特性及功率输出特性^[11，15-16].

以我国典型的陆上1.5 MW风力发电机组为研究对象，采用基于致动线模型的大涡模拟方法对其在动态入流条件下的运行特性进行数值模拟研究. 重点研究在变桨距控制和转矩控制作用下，风力机转速、转矩及功率对入流风速的响应特性，对比验证模拟得到的功率曲线，为所提方法在实际风电场模拟研究中的应用打下基础.

1 数学模型与控制方法 1.1 致动线模型

致动线模型的基本思想是使用虚拟的、承受体积力的致动线代替风力机叶片，模型不求解叶片表面边界层，大大降低了计算成本^[17]. 基于叶素理论，致动线模型将每条致动线沿展向分割为一系列叶素^[18]，每个叶素受到的升力L、阻力D大小为

${\rm{d}}L = \frac{1}{2}{C_{\rm{L}}}(\alpha )\rho U_{{\rm{rel}}}^2c{\rm{d}}r,$

(1)

${\rm{d}}D = \frac{1}{2}{C_{\rm{D}}}(\alpha )\rho U_{{\rm{rel}}}^2c{\rm{d}}r. $

(2)

式中： ${C_{\rm{L}}}(\alpha )$ 、 ${C_{\rm{D}}}(\alpha )$ 分别为升力和阻力系数，α为攻角，ρ为密度，U_rel为空气与风力机叶片的相对速度，c为弦长，dr为叶素径向宽度. U_rel的大小为

${U_{{\rm{rel}}}} = \sqrt {U_{{a}}^{\rm{2}} + {{({\varOmega r} + {U_{{t}}})}^2}}. $

(3)

式中：Ω为风轮转速，U_a、U_t分别为轴向和切向速度，r为叶素距离中心轴的径向长度、速度的矢量关系如图1所示. 由图中的速度三角形可求得U_rel与旋转平面的入流角ϕ，根据入流角ϕ和桨距角θ，计算得到攻角α = ϕ－θ. 确定相对速度和攻角后，计算单位长度叶素的受力^[19]：

${{f}} = \frac{1}{2}\rho U_{{\rm{rel}}}^2c({C_{\rm{L}}}(\alpha){{{e}}_{\rm{L}}} + {C_{\rm{D}}}(\alpha){{{e}}_{\rm{D}}}). $

(4)

式中：e_L和e_D分别为沿升力和阻力方向的方向向量.

每个致动元产生的体积力通过高斯映射反作用于计算域，高斯分布函数的表达式为^[19]

${\eta _\varepsilon } = \frac{1}{{{\varepsilon ^3}{{\text{π }}^{\!\!3/2}}}}\exp \;\left[ { - {{\left( {\frac{\;d\;}{\;\varepsilon\; }} \right)}^2}} \right]. $

(5)

式中：d为流场中网格点与致动点间的距离，ε为高斯光顺参数，Troldborg等^[19]认为ε=2Δx (Δx为致动点附近网格尺度). 流场中任意点(x, y, z)处的体积力为

$f(x,y,z,t) = \sum\limits_{i = 1}^N {f({x_i},{y_i},{z_i},t)} \frac{1}{{{\varepsilon ^3}{{\text{π }}^{\!\!3/2}}}}\exp\; \left[ { - {{\left( {\frac{{{\,d_i\,}}}{\,\varepsilon\, }} \right)}^2}} \right]. $

(6)

式中：(x_i, y_i, z_i)为第i个致动点，t为计算的物理时刻，N为致动点总个数.

1.2 风力机变桨距控制

风力机变桨距控制采用比例积分控制器，使用转速偏差量的比例信号和积分信号的线性组合作为变桨控制的输出量，桨距角变化量为^[20]

$\Delta \theta = {K_{\rm{P}}}{N_{{\rm{Gear}}}}{\rm{\Delta }}{\varOmega } + {K_{\rm{I}}}\int_0^t {{N_{{\rm{Gear}}}}} {\rm{\Delta }}{\varOmega }{\rm{d}}t. $

(7)

式中：K_P、K_I分别为变桨距控制的比例、积分增益，N_Gear为齿轮转速比，ΔΩ为转速偏差. 比例、积分增益的计算方法为^[20]

${K_{\rm{P}}} = \frac{{2{I_{{\rm{Dr}}}}{{\varOmega }_0}\zeta \omega }}{{{N_{{\rm{Gear}}}}( - \partial P/\partial \theta )}},$

(8)

${K_{\rm{I}}} = \frac{{{I_{{\rm{Dr}}}}{\varOmega _0}\omega _{}^{\rm{2}}}}{{{N_{{\rm{Gear}}}}( - \partial P/\partial \theta )}}. $

(9)

式中：I_Dr为传动系统转动惯量，取7.98×10⁶ kg∙m²；Ω₀为风轮额定转速；ω为二阶系统自然频率，ζ为阻尼系数，根据文献[20]，ω和ζ分别取0.6 rad/s和0.7； $\partial P/\partial \theta $ 为气动功率对桨距角的敏感性. $\partial P/\partial \theta $ 是风速、转速、桨距角的函数，为了简化模型，在计算中采用额定风速、额定转速、零桨距角下的 $\partial P/\partial \theta $ 值. 通过引入无量纲修正系数来修正简化带来的误差. 修正系数为

${\rm{GK}}(\theta ) = \frac{1}{{1 + \theta /{\theta _{\rm{k}}}}}. $

(10)

式中：θ_k为桨距角敏感性为额定运行状态下的2倍时对应的桨距角.

修正后的K_P、K_I的计算方法为^[16]

${K^\prime_{{{\rm{P}} }}}(\theta ) = - \frac{{2{I_{{\rm{Dr}}}}{\varOmega _0}\zeta \omega }}{{{N_{{\rm{Gear}}}}(\partial P/\partial\theta )}}{\rm{GK}}(\theta ),$

(11)

${K{^\prime}_{{{\rm{I}}}}}(\theta ) = - \frac{{{I_{{\rm{Dr}}}}{\varOmega _0}\omega _{}^2}}{{{N_{{\rm{Gear}}}}\left( {\partial P/\partial \theta } \right)}}{\rm{GK}}(\theta ). $

(12)

1.3 发电机转矩控制

风力机风轮转子的转速是由风轮的气动转矩和发电机转矩决定的，风轮转速的变化率为

$\frac{{{\rm{d}}{\varOmega }}}{{{\rm{d}}t}} = \frac{{{T_{{\rm{Aero}}}} {\eta _{{\rm_{GB}}}} - {T_{{\rm{Gen}}}} {N_{{\rm{Gear}}}}}}{{{I_{{\rm{Dr}}}}}}. $

(13)

式中：T_Aero为气动转矩，η_GB为齿轮箱效率，取0.98，T_Gen为发电机转矩，N_Gear取104.

发电机转矩是电机转速的函数，根据转速的大小将转矩控制分为5个区域：1、1-1/2、2、2-1/2、3^[20]. 5个区域电机转矩对转速的响应曲线如图2所示，其中 $n$ 为电机转速. 区域1对应低于切入风速的风况，电机转速小于最小工作转速，发电机转矩为0；区域2对应优化运行转速区间，电机转矩与转速的平方成正比，风力发电机组在最优叶尖速比时运行；在区域3，电机功率保持为额定功率，转矩与转速成反比；区域1-1/2为风机启动加速区域，是区域1、2的线性过渡区，通过设定电机转速的下限来限定风力机的工作转速范围；区域2-1/2是区域2、3的线性过渡区，用于限制达到额定运行功率时的转速，减弱额定运行状态时风力机的噪声^[21].

如图2所示，虚线是预先通过FAST软件模拟得到的风力机理想的转矩-转速关系曲线，风力机功率系数在曲线上达到最优. 但是风力机在实际运行中存在最小启动转速，电机转矩也不能无限增大，不可能完全按照理想曲线进行转矩控制，图中实线是实际的电机转矩-转速响应曲线. 区域1与区域1-1/2的分界线为电机启动转速（1 009 r/min）；根据文献[20]，区域2的起始转速取值高于启动转速30%（1 312 r/min）；区域2-1/2的起始转速为额定转速的90%（1 620 r/min）；区域3的起始转速为额定转速（1 800 r/min）.

2 数值方法与计算设定 2.1 控制方程

对流场的流动细节及湍流参数的精确模拟是获得精确的风力机运行参数的重要基础. 基于开源软件OpenFOAM，采用对于瞬态N-S方程求解精度更高的大涡模拟方法模拟动态入流条件下的风力机运行特性. 经过滤波处理的不可压流体的流动控制方程为

$\frac{{\partial {{\bar u}_i}}}{{\partial {x_i}}} = 0,$

(14)

$\begin{split}\displaystyle\frac{{\partial {{\bar u}_i}}}{{\partial t}} + \frac{{\partial ({{\bar u}_i}{{\bar u}_j})}}{{\partial {x_j}}} = - \displaystyle\frac{1}{\rho }\displaystyle\frac{{\partial \bar p}}{{\partial {x_i}}} + \nu \displaystyle\frac{{{\partial ^2}{{\bar u}_i}}}{{\partial {x_j}\partial {x_j}}} - \displaystyle\frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}} + {f_i}.\end{split}$

(15)

式中： $\bar u$ 为滤波后的速度（i, j=1、2、3，分别代表x、y、z方向的速度）， $\bar p$ 为滤波后的压力，ν为运动黏度，τ_ij为亚格子应力项，f_i为体积力，由式（6）计算得出. 大涡模拟亚格子模型采用Smagorinsky模型，动量方程的求解采用PISO（pressure-implicit with splitting of operators）算法，时间项的离散采用Backward格式，空间离散采用高斯线型格式.

2.2 计算域网格划分

模拟的对象为内蒙古某风电场实际运行的额定功率为1.5 MW的风力机，风力机直径D=77 m，中心高度为70 m. 风力机切入风速、切出风速和额定风速分别为3、25、11 m/s，额定转速为17.3 r/min. 选用典型的长方体计算域，长度为1 232 m（16D）、宽度为616 m（8D）、高度为312 m（约4D），风力机布置在距离入口边界面6D处的位置，计算域示意图如图3所示. 为了优化网格，划分网格时采用文献[22]的方法，从外到内进行双重加密. 文献[22]中最小网格尺度为3 m，为了满足网格无关性要求，将最小网格尺度细化到2 m. 如图3所示，外侧D1区域，网格分辨率为8 m，D2区域经过一重加密，网格分辨率为4 m，最内侧D3区域，再次加密后网格的分辨率为2 m，整个计算域的总网格数约为188万.

2.3 边界条件

为了模拟动态入流条件下风力机的运行特性，采用用户自定义的方法编译动态速度入口边界条件，引入正弦函数定义速度随时间的周期波动. 为了使入口速度更贴近真实的大气入流，使用对数风廓线定义风速沿垂直方向上的分布. 入口速度是时间t和垂直高度H的函数：

${u_{{\rm{in}}}}(z,t) = {U_{{\rm{ref}}}}\frac{{\ln\; (z/{z_0})}}{{\ln\; ({H_{{\rm{ref}}}}/{z_0})}}({c_v} + \left| {\sin\; ({\rm{\pi }}t/T)} \right|). $

(16)

式中：U_ref为参考速度，取风力机额定风速11 m/s；H_ref为参考高度，取风力机轮毂中心高度70 m；z为垂直高度；z₀为粗糙度长度，取0.03 m；c_v为最小速度修正常数，为了保证风力机轮毂中心高度处速度高于风力机切入风速（3 m/s），c_v取0.3；T为速度波动周期，取300 s. 根据以上定义得到不同时刻的入口速度风廓线，如图4所示. 使用的边界条件如表1所示.

2.4 时间步长

依据CFL(Courant-Friedrichs-Lewy)准则选定时间步长，保证流体在单位时间步长内的位移不超过1个网格间距，时间步长满足

$\frac{{U\Delta t}}{{\Delta s}} \leqslant 1. $

(17)

式中： $U$ 为流体的速度， $\Delta t$ 为时间步长， $\Delta s$ 为网格尺寸.

基于致动线模型的风力机流动模拟的最大流动速度出现在叶尖. 时间步长的选择主要受限于叶尖速度，叶尖处的流动需要满足CFL准则^[23]. 将叶尖速度和叶尖处网格尺寸代入式（17）中，得出时间步长应当满足Δt≤0.029 s，时间步长取0.02 s. 计算物理时长为3倍的动态入流周期，即900 s. 模拟采用64核（处理器型号：Intel Xeon E5-2660）并行计算，总计算时长约为16 h.

3 结果与分析 3.1 转速的响应

根据2.3节中定义的入流条件，风力机轮毂中心高度处风速随时间的波动曲线如图5所示，图中水平虚线表示风力机额定风速U_rated，根据入流条件计算出模拟时间段内达到额定风速的临界时刻分别为t=74、226、374、526、674、826 s.

在给定的动态入流条件下，风力机风轮转速随入流速度的变化而波动，当风力机转速超过额定转速时，控制变桨距以抑制转速的进一步增长. 如图6所示为动态入流条件下风轮转速与桨距角响应曲线模拟结果，图中水平虚线表示风轮额定转速Ω_rated的位置，垂直虚线表示变桨距控制作用的起始与终止时刻. 结果显示，第2个入流周期与第3个入流周期风轮转速及桨距角随时间的变化基本一致，表明在当前入流条件下，流动在第2个入流周期已经基本稳定，不再受初始条件影响，模拟算到第3个入流周期后终止. 模拟结果显示，风轮转速分别在t=74、374、674 s时达到额定值，与图5中入流风速达到额定风速的时刻吻合. 当风轮转速达到额定转速后，由于入流风速还在继续增加，为了控制风轮转速，变桨距控制开始起作用，桨距角随之逐渐增加. 在变桨距控制的作用下，风轮转速在短暂超过额定转速后，很快又回落到额定值，风轮转速相对于额定转速的最大偏移量约为0.35 r/min. 当风轮转速回落到额定值后，桨距角不再增加，逐渐减小到0，在此过程中，风轮转速比额定值有微小的减小. 模拟结果表明，在3个入流周期中，风轮转速分别从t=236、536、836 s时开始明显减小，同时变桨距控制结束，相比于图5中转速回落到额定转速的时间各有10 s的延迟.

3.2 转矩的响应

在给定的动态入流条件下，风力机气动转矩与电机转矩的响应如图7所示，为了使两者具有可比性，根据式（13），对比气动转矩与齿轮箱效率的乘积和发电机转矩与齿轮转速比的乘积. 结果表明，气动转矩和电机转矩在t=74、374、674 s时达到峰值，正是入流风速达到额定风速即变桨距控制开始作用的时刻. 下面以第3个入流周期为例分析1个周期内转矩的变化规律.

在入流风速达到额定风速之前（t<674 s），气动转矩大于电机转矩，根据式（13），风轮转子处于加速阶段. 当t>674 s时，变桨距控制开始，气动转矩随着桨距角的快速增大骤降至略小于电机转矩，该过程对应于风轮转子先加速到超过额定转速然后逐渐减小到接近额定转速. 当674 s<t≤836 s时，转子转速基本维持在额定转速，电机转矩受转矩控制作用基本维持在额定转矩，风轮转矩在入流风速变化和变桨距控制的作用下，在额定转矩附近存在微小波动. 当t>836 s时，入流风速降至额定风速以下，气动转矩迅速减小，转子减速，电机转矩随转速的减小而减小，在该过程中电机转矩大于气动转矩，转子一直处于减速状态. 从t=674 s时的转矩对比结果可以看出，当入流风速处于额定转速、桨距角为0时，气动转矩大于电机转矩. 当入流风速降到额定风速时，为了保证转子转速维持在额定转速附近，桨距角没有马上降为0，而是继续起作用直到气动转矩与电机转矩相等，这是图6中存在着10 s延迟时间的原因.

3.3 功率的响应

在给定的动态入流条件下，风力发电机输出功率响应曲线如图8所示. 为了更清晰地观察输出功率与入流风速之间的关系，将风力机轮毂中心高度的入流风速也呈现在图8中. 图中的2条水平虚线分别表示风力机额定功率P_rated和额定风速U_rated，垂直虚线分别代表模拟时间段内达到额定功率或额定风速的临界时刻. 结果表明，在入流风速增长到额定风速的时刻（t=74、374、674 s），输出功率也同时达到额定值；此后，虽然入流风速超过额定风速，但是在变桨距控制和转矩控制的共同作用下，风力机输出功率维持在额定功率附近；与转速和转矩的变化一样，当入流风速开始降到额定风速以下时（t=226、526、826 s），输出功率并没有马上下降，同样存在10 s的延迟，之后功率才随着入流风速的下降而减小. 提取第3个入流周期内的加速阶段（t=600~750 s）的功率与入流风速数据，得到风力机功率曲线，与风力机实际的功率曲线进行对比，结果如图9所示. 可以看出，当入流风速小于9 m/s时，模拟的功率曲线与实际功率曲线存在一定的误差，最大误差约为75 kW（额定功率的5%）；当入流风速大于9 m/s时，模拟的功率曲线与实际功率曲线基本吻合. 实际上，当风速较低时，相应的转速也低，风速小于9 m/s的阶段对应着转矩控制的1-1/2区和2区的前段，偏差主要是在1-1/2区段产生的，区段转矩与转速采用线性控制方法，精度相对不高. 在将来的研究中，可以通过优化1-1/2区段的转矩控制方式进一步提高模拟的精度.

4 结　论

（1）采用提出的变桨距控制和转矩控制方法，风力机转速能够较好地响应动态入流风速的变化，转速在额定入流风速下达到额定值. 当入流风速超过额定风速时，变桨距控制能够及时作用，使风轮转速维持在额定转速附近，转速相对于额定转速的最大偏移量约为0.35 r/min.

（2）风轮气动转矩和电机转矩对入流风速变化的响应较好. 转矩控制使电机转矩在风力机达到额定转速后维持在额定转矩，变桨距控制在风力机达到额定转速后迅速启动，气动转矩减小，风力机转速维持在额定转速附近. 当入流风速减到额定风速时，由于气动转矩与电机转矩在零桨距角运行时存在转矩差，变桨距控制继续作用10 s后才停止，相应的转速、转矩、功率的下降也存在10 s的延迟.

（3）风力机输出功率能够较好地响应转速的变化，当风速超过额定风速时，在转矩控制与变桨距控制的作用下，风力机输出功率维持在额定功率. 模拟得到的功率曲线与实际功率曲线整体上较吻合，但是在风速小于9 m/s的区间段还存在一定的偏差，需要进一步优化. 功率曲线的对比结果也表明数值模拟框架对于动态入流的模拟具有较高的精度. 今后将会在真实风电场的模拟研究中进一步验证和优化提出的模拟方法.