极端服役环境下的风电机组塔架结构参数优化研究

戴巨川 刘 旋 杨书仪 文泽军 沈祥兵

1.湖南科技大学，湘潭，411201 2.湘电风能有限公司，湘潭，411000

0 引言

风能是一种可再生的清洁能源，是各国新能源开发和利用的战略选择，近年来对它的开发利用发展迅速。塔架作为风电机组的支撑结构，其稳定性与可靠性直接影响整个机组的运行安全，塔架成本昂贵，制造成本占整机成本的15%～20%。因此，近年来越来越多的学者对塔架设计优化进行了研究：Negm等［1］以横截面面积、壁厚和高度为变量，最小塔架质量、最大刚度和最大自然频率为目标进行了塔架优化设计；Horvath等［2］以塔架等应力分布为目标，采用有限元方法（FEM）和计算流体力学（CFD）方法探讨了塔架高度和横截面参数取值范围；Uys等［3］同时考虑材料成本与制造成本，以综合成本最小为目标，获得了塔架壁厚、加强筋数目和直径等的优化值；Clifton－smith等［4］运用有限元法研究了塔架的屈曲稳定性，对屈曲应力计算公式进行了修正；Hassan等［5］以塔架固有频率最大为目标，塔架质量为约束，对塔架的壁厚变化方式进行了优化设计；Yildirim等［6］以塔架质量最小为目标，采用遗传算法进行了优化求解；赵吉文等［7－8］以固有频率最大为目标，采用支持向量机建立了塔架结构参数与固有频率之间的快速计算模型，引入粒子群优化、混沌搜索等方法进行优化求解。尽管如此，针对风电机组塔架的优化设计仍有许多问题需要解决，如获取塔架直径、壁厚等结构参数与应力、固有频率及寿命等的关联特性，寻求合理的替代模型解决优化求解时需迭代计算带来的耗时长、效率低等问题。

本文采用一种壁厚分段式线性变化塔架结构模型，借助有限元数值仿真分析结构参数与应力、固有频率间的关联特性，采用支持向量机构建结构参数与应力、固有频率间的快速计算模型。

1 塔架几何模型

大型风电机组塔架大多采用按一定锥度变化的变截面圆筒形结构。为方便运输，通常将塔架分为几节，安装时用法兰将其连接成整体。考虑到塔架应力分布与圆筒半径、壁厚等为非线性关系，为尽可能地改善塔架应力分布，本文采用分段式线性变化塔架结构，分为三节，如图1所示。塔架主要结构参数包括底端直径、顶端直径、塔架总高度、各节塔架高度及厚度。由于风电机组在设计时一般以轮毂高度处的风速为设计依据，而塔架总高度主要由风电机组容量、风场风速和所处地理位置等决定，故进行塔架设计优化时该值视为常量；同时考虑到塔架顶端与机舱配合，顶端直径值也可视为常量。

图1 塔架结构示意图

2 风电机组载荷分析

对风电机组进行载荷分析，获得准确的塔架载荷是进行塔架受力分析、优化设计的基础。风电机组承受的载荷特别是气动载荷非常复杂，由于来流的时变性，作用在风轮上的气动载荷在风电机组运行过程中是时刻变化的，是影响风电机组塔架载荷变化的主要因素。为便于分析，本文采用风电机组坐标系如图2所示，叶片上某一截面处的载荷分解为坐标系xeyeze中沿xe轴方向分力Fxe和绕xe轴的力矩Mxe（挥舞力矩），沿ye轴方向分力Fye和绕ye轴的力矩Mye（摆振力矩），沿ze轴方向分力Fze和绕ze轴的力矩Mze（变桨距力矩）。

图2 风电机组坐标系［9］

图2 中，γ为风轮锥角，τ为风能仰角。坐标系xgygzg为惯性坐标系，固定于塔架底部；坐标系xnynzn固定在机舱上，其yn轴在风电机组主轴轴线上；坐标系xsyszs固定在主轴上，与主轴同步旋转，转速为n；xpypzp为假定无风能仰角时，固定于旋转主轴上的坐标系；坐标系xeyeze固定在旋转的叶片1上，ze轴沿叶片变桨轴线方向，ye轴垂直于ze轴指向塔架（上风向风电机组），xe轴垂直于叶片变桨轴线和主轴轴线。

目前，进行风电机组气动载荷分析有多种方法，如叶素—动量（BEM）理论、基于动态入流理论的GDW方法、CFD等。在实际工程中，由于BEM理论简便可信，故是常用的方法。以下仅对基于BEM理论的气动载荷计算方法进行简约介绍。在BEM理论中，取半径r处长度为dr的叶素为对象，气流速度与空气动力分量如图3所示［9－10］。

图3 叶素上气流速度与空气动力分量

图3中，A为前缘点，C为压力中心，B为后缘点。v0为轴向风速，v1为气流相对于叶素的速度；ω为风轮角频率；a为轴向诱导速度系数，a′为切向诱导速度系数；α为攻角，β为节距角，φ为来流角度；dR为叶素上的气动力，dL为dR的升力分量，dD为dR的阻力分量，dQ为气动力切向分量（与旋转面相切），dN为气动力轴向分量；ve－Tip和ve－Top分别是塔架振动引起的x方向和y方向速度，ve－Bip和ve－Bop分别是叶片振动引起的x方向和y方向速度。对于叶片上不同半径r处的叶素而言，ve－Tip和ve－Top不变，而ve－Bip和ve－Bop随当前叶素半径r的变化而变化。

作用在元素上的升力和阻力分别为

式中，c为半径r处的叶片弦长；Cl、Cd分别为升力系数和阻力系数；ρ为空气密度。

dL和dD的合力dR可以分解为轴向力dN和切向力dQ：

产生的摆振方向弯矩为

产生的挥舞方向弯矩为

式中，Ct、Cn分别为切向力系数和法向力系数。

设v为风轮中心高度h0处的来流风速，如果风速分布服从赫尔曼Holman规律［11］，则有

式中，r′为风速廓线指数，与地面粗糙度有关；ψ为叶片方位角。

沿叶片方向将整个积分空间［0，R］分成若干个子空间［ri，ri＋1］，i＝1，2，…，N，叶片产生的摆振方向总弯矩可写为

式中，ωBx、ωBy分别为x、y方向叶片的振动角速度。

同样，叶片产生的挥舞方向总弯矩可以写为

叶片摆振方向剪力为

叶片挥舞方向剪力为

同一空间向量，在不同坐标系中坐标值不同，通过变换矩阵可以实现相互转换，变换矩阵的获取可采用方向余弦或欧拉方法。这样就可以把风轮上的气动载荷转化至塔架上。

除了风轮上气动载荷以外，塔架自身上风载和风电机组运行载荷（变桨距）、惯性载荷等对塔架也有影响，本文对此未做介绍，相关理论可参考有关文献。在具体确定设计载荷时，《风力发电机组规范》（中国船级社，2008）规定以具体的装配、吊装、维修、运行状态或设计工况同外部条件的组合为依据，应考虑到以合理概率出现的相关载荷，如极端外部条件与故障工况相关，两者应组合为一种载荷情况。对于某2MW风电机组，作用在塔架基础上的极限载荷按表1给定（由相关厂家提供）。

表1 塔架极限载荷

3 结构参数与应力、固有频率关联特性分析

塔架结构参数对应力、固有频率、屈曲稳定性、疲劳寿命和动态响应特性等都有直接影响，获取结构参数与这些属性的关联特性是塔架优化设计的前提。一般来说，塔架在极端服役环境下发生破坏的可能性较大，因此本文选取塔架结构参数与应力、固有频率间的关系进行分析，不考虑其疲劳寿命和动态响应的影响。

3.1 结构参数与应力关联特性

在有限元分析软件ANSYS中采用壳单元进行塔架建模，考虑到叶片、轮毂和机舱等部件质量影响塔架固有频率，将叶片、轮毂和机舱简化为一集中质量点，即mass21单元节点，用mpc184无质量刚性梁单元连接mass21单元与塔架顶端节点，材料属性的弹性模量定义为205GPa，泊松比为0.3，材料密度为7.85×103kg／m3，塔架壁厚值大小通过设置壳单元实常数实现，塔架底部节点全约束作为边界条件。塔架基础上的载荷是塔架顶部存在的推力产生的，塔架顶部推力大小可由基础弯矩除以塔架高度得到，在mass21节点上施加等效载荷后可求解塔架的最大应力。考虑设计的具体情况，本文选取上段壁厚δ1变化范围为12～18mm，中段壁厚δ2变化范围为18～26mm，下段壁厚δ3变化范围为26～45mm，下段塔架长度H1变化范围为15～20m，中段塔架长度H2变化范围为20～35m，塔架底端直径D变化范围为3.8～4.4m，由于不同结构参数的数量级不同，为便于比较各结构参数与塔架应力的关系，对结构参数进行归一化处理，得到如图4所示的塔架结构参数与最大应力的关系曲线。为观察不同塔架结构参数组合下最大应力点位置，选取4组结构参数求解其最大应力，最大应力及其节点位置如表2所示。

图4 塔架最大应力与结构参数的变化关系

表2 最大应力值与节点位置

从图4可知，不同结构参数对塔架的影响程度不同，在图4a中，随着δ1增大，最大应力值显著减小，但是当δ1增加到一定程度以后，对应力的影响逐渐放缓；δ2和δ3的变化对最大应力值没有显著影响，总体趋势比较平缓。从图4b可以看出，随着塔架单段长度、塔架底端直径增大，最大应力值总体趋势均是逐渐减小。在表2中，4组不同塔架结构参数对应的最大应力值分别为153MPa、195MPa、177MPa和209MPa，分别位于距塔架底端48m、38m、18.5m和6m处。不难看出，塔架结构参数组合不但影响塔架的最大应力，而且影响了最大应力点的位置。

3.2 结构参数与固有频率关联特性

风电机组在运行过程中，叶轮、塔架等相互激励，有可能发生共振引发机组倒塌事故，因此在设计时，应充分考虑各子系统固有频率，避免共振的发生。接近风轮旋转频率的塔架一阶固有频率最易导致共振，对其采用有限元模态分析法进行求解。需要指出的是，在进行模态分析时，叶片、轮毂和机舱等部件质量不能忽略，采用集中质量的方式简化，建模时在集中质量重心处设置一个质量节点单元mass21，用MPC184无质量刚性梁单元连接质点单元与塔顶节点，以底部节点全约束为边界条件求解塔架固有频率。为便于比较各结构参数与塔架固有频率的关系，对结构参数进行归一化处理。图5所示为结构参数与塔架一阶固有频率的关系曲线，表3所示为4组结构参数不同时塔架的固有频率。

图5 塔架一阶固有频率与结构参数关系

表3 塔架固有频率

从图5可知，随着各结构参数的增大，固有频率均呈增大趋势；不同结构参数对塔架固有频率的影响程度不同，在所分析的结构参数中上段壁厚δ1和下段直径D的变化对塔架的固有频率影响较大，下段壁厚δ3和中段塔架长度H2影响较小。从表3中可以看出，不同结构参数组合下，塔架固有频率值差异比较明显，表中4组不同结构参数下的塔架一阶固有频率变化范围为0.372 09～0.532 96Hz，而风电机组运行时风轮最大旋转频率为0.375Hz左右，第3组、4组参数下塔架易发生共振。

4 结构参数与应力、固有频率快速计算模型

塔架结构参数与应力、固有频率之间是一种复杂非线性关系，近似的数学表达式计算误差大，应该采用有限元计算结果。由于在优化设计时需要反复多次计算，并且优化算法存在不确定性，一次优化结果一般不能作为最终结果，需要多次优化后进行结果比较确定，直接采用有限元分析软件计算效率较低，寻求合理的替代模型可以有效解决这一问题。在各种人工智能方法中，支持向量机在解决小样本、非线性及高维模式识别中表现出许多特有的优势，并能推广到函数拟合等其他机器学习问题中。本文采用支持向量机建模方法建立塔架结构参数与应力、固有频率间的快速计算替代模型。

4.1 样本数据获取

为了构建基于支持向量机的结构参数与应力、固有频率的快速计算模型，首先需要获得足够的样本数据。采用正交试验设计与随机组合方法相结合的混合试验设计方法确定参数组合，通过有限元数值仿真试验获得样本值。根据塔架设计中可变的结构参数确定正交试验的6个因素，即因素A（塔架上段壁厚δ1）、因素B（塔架中段壁厚δ2）、因素C（塔架下段壁厚δ3）、因素D（塔架下段长度H1）、因素E（塔架中段长度H2）、因素F（塔架底端直径D）；根据塔架结构参数与最大应力、固有频率的变化关系确定每个因素的5个水平值。试验因素水平表见表4，按照正交试验方法需要进行L25（56）＝25组试验。除正交试验获得的25个样本数据外，再从整个排列组合156 25种试验方案中随机选取200组进行试验，共获得225个样本数据。

表4 塔架结构参数因素水平表

4.2 快速计算模型的构建

样本空间为 ｛（x1，z1），（x2，z2），…，（x225，z225）｝，其中，xi＝（δ1i，δ2i，δ3i，H1i，H2i，Di），为塔架结构参数；zi＝（σi，fi），为塔架应力和固有频率；i＝1，2，…，225。支持向量机中采用的方法是将输入向量xi通过某种非线性关系映射到高维空间，则在高维空间回归函数为

式中，φ为一个从输入空间到高维特征空间的映射函数；b为偏差值；w为权向量。

式（12）即为待求的塔架结构参数与应力或固有频率替代模型。设拟合精度为ε，求解最佳回归函数转换为固定间隔下寻找最小的‖w‖问题。采用ε－SVR型支持向量，考虑样本的容错性，给拟合精度ε这一硬性值加一松弛变量ξ，则有［12］：

式中，C为惩罚因子为松弛变量。

使用Lagrange乘子法解决上述问题，得到原问题的对偶问题为

对偶问题完全根据样本数据表达，由样本数据求得最优的Lagrange乘子后即可得到问题的最优解。则式（12）可写为

将通过数值试验获得的225组样本数据分为两部分，200组用于训练，其余用于测试。支持向量机类型选用ε－SVR类型，核函数设置类型选用径向基核函数（RBF）。在MATLAB中分别用网格搜索法、粒子群方法和遗传算法函数求解核函数参数，选取平均平方误差（MSE）最小、平方相关系数最接近1的核参数，最终应力模型中取核参数c＝68.3341，g＝0.083 923，固有频率模型中取核参数c＝24.5844，g＝0.003 814 7。

基于上述方法得到的应力回归模型由179个支持向量组成，偏差值为－0.5034，固有频率回归模型由153个支持向量组成，偏差值为0.0192。表5给出了5组模型输出与有限元仿真结果，可以看出构建的模型最大应力输出值与有限元仿真结果误差为－0.004%～4.7%，一阶固有频率输出值与有限元仿真结果误差为－0.005%～0.3%。结果表明，构建的塔架结构参数与最大应力、一阶固有频率的替代模型可用于工程实际。

表5 模型输出值与有限元仿真结果

5 塔架结构参数优化

5.1 参数优化数学模型

选取塔架的结构参数为设计变量为

以塔架总质量最小为目标，即

式中，mtower为塔架总质量，mi为每段塔架的质量（可由塔架结构参数求得）。

约束条件包括：

（1）强度约束。为保证塔架的安全，外载荷作用下的塔架最大应力σ应小于材料的屈服应力，即：

式中，nst为安全系数；σ为材料屈服应力。

（2）固有频率约束。风电机组运行时，如果风轮旋转频率或叶片通过频率接近塔架的固有频率，则会出现共振现象，引发严重的安全问题。某型风电机组运行时风轮的最大转速为22.5r／min，风轮旋转频率和叶片通过频率分别为0.375Hz和1.125Hz，根据GB18451.1－2001规定，塔架工作状态下的弯曲固有频率f0，n和激振频率fR、fR，m之间应有适当的间隔，按下述公式计算：

式中，fR为正常运行范围内风轮的最大旋转频率；f0，1为塔架（整机状态下）的一阶固有频率；fR，m为m个风轮叶片的通过频率；f0，n为塔架（整机状态下）的n阶固有频率。

由式（19）得到塔架一阶固有频率约束为

（3）边界约束。根据工程具体情况，塔架的结构参数在一定范围内取值，这里选取变量的范围为

本文采用遗传算法进行参数优化求解，遗传算法是一种基于生物遗传和进化机制的自适应概率优化算法，具有很好的全局寻优能力。对于一些非线性、多模型、多目标的函数优化问题用其他优化方法较难求解，而遗传算法可以方便地得到较好的结果。种群中的个体数目设为30，采用浮点编码方式，交叉概率为0.95，变异概率为0.08，最大迭代次数为250次。

5.2 优化结果

基于上述方法，得到塔架结构优化参数如表6所示，同时给出了优化前后塔架最大应力和固有频率的有限元数值计算结果，优化后的塔架总质量为100.4t，相比优化前的104t，用钢量减少3.5%（3.6t），一阶塔架固有频率从0.4788Hz变为0.4793Hz，最大应力从124MPa降为120MPa。图6给出了优化前后塔架应力分布曲线，优化后应力变化更平缓，塔架的安全性有所提高。

表6 塔架结构参数优化结果

图6 塔架应力分布曲线

6 结论

（1）通过试验设计和有限元数值仿真分析了各参数变化对应力、固有频率的影响趋势，获得了塔架结构参数与应力、固有频率的关联特性，在塔架设计时可为各结构参数的分布提供一定的参考。

（2）采用支持向量机理论分别构建了塔架结构参数与应力、塔架结构参数与固有频率间的快速计算模型，计算结果显示，替代模型输出结果能与有限元仿真结果较好吻合，表明本文构建的快速模型是可行的。

（3）构建了某2MW风电机组以塔架质量最小为目标的塔架结构参数优化设计模型，采用遗传算法进行优化求解，优化后的塔架质量减少3.6t，减少用钢量3.5%，为风电机组塔架实际工程设计优化提供了一种参考方法。

［1］Negm H M，Maalawi K Y.Structural Design Optimization of Wind Turbine Towers［J］.Computers ＆Structures，2000，74（6）：649－666.

［2］Horvath G，Toth L.New Methods in Wind Turbine Tower Design［J］.Wind Engineering，2001，25（3）：171－178.

［3］Uys P，Farkas J，Jarmai K，et al.Optimisation of a Steel Tower for a Wind Turbine Structure［J］.Engineering Structures，2007，29（7）：1337－1342.

［4］Clifton－smith M，Wood D.Optimisation of Selfsupporting Towers for Small Wind Turbines［J］.Wind Engineering，2010，34（5）：561－578.

［5］Hassan A，Maalawi K，Negm H.Frequency Optimization of Conical Shells Under Mass Equality Constraint［J］.JASR，2006.

［6］Yildirim S，Özkol I·.Wind Turbine Tower Optimization under Various Requirements by Using Genetic Algorithm［J］.Engineering，2010，2（8）：641－647.

［7］赵吉文，谢芳，刘永斌，等.600kW风力机塔架结构参数粒子群优化设计［J］.中国机械工程，2007，18（15）：1776－1779.Zhao Jiwen，Xie Fang，Liu Yongbin，et al.Optimization Design of Tower Structure of 600kW－Wind Turbine Based on Particls Swarm［J］.China Mechanical Engineering，2007，18（15）：1776－1779.

［8］赵吉文，刘永斌，陈军宁，等.基于混沌搜索的600kW风力机塔架结构优化设计［J］.系统仿真学报，2007，19（11）：2425－2428.Zhao Jiwen，Liu Yongbin，Chen Junning，et al.Optimization Design of Tower Structure of 600kWWind Turbine Based on Chaotic Search［J］.Journal of System Simulation，2007，19（11）：2425－2428.

［9］Dai Juchuan，Hu Yanping，Liu Deshun，et al.Aerodynamic Loads Calculation and Analysis for Large Scale Wind Turbine Based on Combining BEM Modified Theory with Dynamic Stall Model［J］.Renewable Energy，2011，36（3）：1095－1104.

［10］戴巨川，胡燕平，刘德顺，等.MW 级变桨距风电机组叶片转矩计算与特性分析［J］.太阳能学报，2010，31（8）：1030－1036.Dai Juchuan，Hu Yanping，Liu Deshun，et al.Torque Calculation and Characteristics Analysis of Blade for MW Adjustable－ptich Wind Turbine［J］.Acta Energiae Solaris Sinica，2010，31（8）：1030－1036，2（11）：821－829.

［11］Gourieres D L.Wind Power Plants：Theory and Design［M］.Oxford：Pergamon Press，1982.

［12］Vapnik V.Statistical Learning Theory［M］.New York：Wiley，1998.