基于BESO算法的高桩承台式水平轴风力机支撑结构优化

占玲玉, 何文君, 周 岱,2, 韩兆龙,2, 朱宏博, 张 凯, 涂佳黄

(1. 上海交通大学 船舶海洋与建筑工程学院,上海 200240; 2. 上海交通大学 海洋工程国家重点实验室,上海 200240; 3. 湘潭大学 土木工程与力学学院,湖南 湘潭 411105)

风力机大型化可有效提高风力发电量和风能利用率,是海上风力发电的世界发展趋势[1].研发建造包括10、15和20 MW在内的大型风力机是我国海洋绿色能源开发利用的战略重点.然而随着风力机转子、塔筒尺寸的不断增加[2],大型风力机结构系统柔性增加,抗风能力明显降低[3],风力机风致动力响应更加复杂,系统风致动力灾变风险更加凸显.

研究并优化大型水平轴风力机支撑结构对保证系统运行安全十分重要.风力机大型化导致塔筒顶部大型转子和发电机质量增加、支撑塔筒高度增大.常用的圆锥塔筒无法满足风力机结构的强度和刚度要求.迄今,学者们对大型水平轴风力机支撑结构优化已有一定研究.Gentils等[4]基于遗传算法,考虑多准则约束,对美国国家可再生能源实验室(NREL)提出的5 MW单桩式风力机(NREL 5 MW)塔架和基础进行优化分析并计算风力机支撑结构优化尺寸参数;Oest等[5]基于疲劳损伤荷载、准静态分析和动态分析,对海上风力机支撑结构进行梯度优化设计;葛旭等[6]结合有限元分析与进化策略,对海上三桩单立柱风力机支撑结构主尺度参数进行多参数同步优化设计;Tian等[7]对导管架式海上风力机的导管架基础进行拓扑优化,结合尺寸和形状优化实现结构重分布,得到质量更轻的导管架结构.上述研究仅为针对特定风力机结构形式的参数优化,未改变其支撑结构形式,虽对风力机结构稳定性有一定提高,但设计方法具有局限性.相较于参数优化,拓扑优化方法具有更多的设计自由度,是优化大型水平轴支撑结构的新途径.

拓扑优化方法应用于包括结构优化在内的广泛领域,其双向渐进结构优化(BESO)算法优化效果较好.王章骏等[8]采用BESO算法对吸能管进行拓扑优化,提高了带隔板薄壁方管的耐撞性能, 实现了吸能管轻量化的设计目标;刘昆等[9]基于BESO算法对船体梁腹板开孔进行设计优化,提出适用于船体结构孔型的BESO算法;BESO还在罗马小体育宫[10]、卡塔尔国家会议中心[11]等建筑工程设计中得到应用.目前,大型水平轴风力机塔筒支撑结构拓扑优化还鲜有研究.

基础是风力机系统的关键结构之一,高桩承台基础的结构安全性能高,适合我国海床地质条件,施工工艺成熟,在我国近海风电场建设中多有应用[12].但不考虑桩-土相互作用会导致风力机位移响应被严重低估[13-15],是结构安全设计的重要隐患.

考虑桩土相互作用的风力机支撑结构性能,开展支撑结构优化分析.基于反比例型删除率的BESO算法,提出大型高桩承台式水平轴风力机支撑结构拓扑的优化算法,运用有限元软件ANSYS开展大型水平轴风力机支撑结构寻优分析,获得一种风力机优化支撑结构,显著降低了风力机风致动力响应.

1 风力机结构建模与分析

1.1 风力机支撑结构建模

以应用较广的NREL 5 MW[16]风力机为背景,进行考虑桩土作用的近海坐底式高桩承台大型水平轴风力机塔筒结构优化研究.高桩承台式水平轴风力机初始模型如图1(a)所示,其支撑系统包括上部塔筒结构、承台和群桩.水平轴风力机塔筒为圆筒式支撑结构,如图1(b)所示,塔筒底部外径D=6 m,塔筒高度H0=90 m[16].塔筒选用Q345钢材[17],密度ρ0=7.85 g/m3,弹性模量E0=211 N/m2,泊松比ν0=0.3,风力机额定风速为11.4 m/s,具体参数如表1所示.

表1 5 MW水平轴风力机塔架参数Tab.1 5 MW horizontal axis wind turbine parameters

图1 水平轴风力机塔筒结构初始模型Fig.1 Horizontal axis wind turbine initial model

参照上海东海大桥海上大型风力发电工程项目采用高桩承台式基础[18],混凝土承台直径D1=20 m,承台厚度S=5 m,设置4根直径d1=3 m、壁厚t=20 mm的全直钢管基桩,桩基长度H1=52 m,入土深度h=40 m.承台材料选用C45标号混凝土,弹性模量E=33.5 GPa,密度为2 500 kg/m3,泊松比ν=0.2[19],基础模型如图2(a)所示.基于Jung等[20]和Chen等[21]确定土体模型和力学参数:土体浮重度γ′、弹性模量E′、泊松比ν′、内摩擦角φ′和不排水抗剪强度su、黏土泊松比νu、50%的最大有效应力下的轴向应变ε50、土层深度z、膨胀角ψ′、黏聚力c′、地基反力模量k,如图2(b)所示.

图2 基础与土体模型Fig.2 Foundation geometric model and soil profiles

1.2 风力机支撑结构风荷载

对高桩承台大型水平轴风力机结构系统,基于计算流体动力学(CFD)方法并运用Star-CCM+计算软件[22]数值模拟获得风力机叶片气动荷载,风力机风轮直径为d,选取长方体风场计算域,计算域长度、宽度和高度分别为20d、8d和8d,如图3(a)和3(b)所示.计算域划分共 1 700 万个网格.计算域中,风力机叶片圆柱体旋转域在竖平面的直径为1.5d,旋转域高度即顺风向尺度为0.75d;旋转域之外的计算域等其他部分为静止域.对旋转域和静止域设置交界面,采用滑移网格处理风轮运行旋转问题,如图3(c)和3(d)所示.叶片表面采用棱柱层网格,延伸因子为1.3,总厚度为0.08 m,无量纲化壁面距离为y+<1.

图3 风场计算域和边界条件 Fig.3 Wind field computation domain and boundary condition

对大型风力机支撑结构采用锥筒结构[23],塔筒高度为90 m.将塔筒简化为受弯曲变形、轴向压缩和扭转变形的悬臂梁.风力机运行时,需考虑塔筒沿高度的风压变化.沿塔筒高度l处的风荷载[24]为

(1)

式中:Cd为风阻力系数,取值1.2;ρ为空气密度;Dl为塔筒高度l处的塔筒截面直径;V(l)为沿塔架高度l处的瞬时风速.

1.3 桩土相互作用

美国石油协会(API)规范推荐的土体p-y曲线,即非线性弹簧模型广泛应用于海洋工程桩基设计.引入用于分析风力机支撑结构的桩土相互作用[25],该模型考虑了土体非线性、土层变化和土体种类等因素的影响,软土中各层p-y曲线为

(2)

式中:p为深度X处单位桩长的极限水平土抗力标准值;yc为实验室中对未扰动工试样做不排水压缩试验时,其应力达到最大应力一半时的变形;y为泥面以下深度X处桩的侧向水平变形;XR为泥面以下到土抗力减少区域底部的深度;pu为深度X处单位桩长的极限水平土抗力标准值.

对于砂土情况,p-y曲线为

(3)

式中:K为地基反力初始模量;A为计入静力荷载和循环荷载条件的参数.

基于API规范的p-y曲线法,计算循环荷载作用情况下的桩侧土抗力,使用ANSYS软件的Combine39单元模拟p-y曲线,该单元是具有非线性功能的单向拉压单元,常用于桩土相互作用模拟.在搭建有限元模型时,根据计算得到的土体p-y曲线参数,每间隔1 m设置一个Combine39土弹簧单元.在桩土相互作用模拟过程中,为简化计算未考虑桩土循环弱化作用及桩端设置竖向简支约束.

1.4 风力机支撑结构系统振动响应

采用 Newmark-β法[26]中的平均常加速度法对风力机支撑结构系统风致振动进行瞬态分析,其中控制Newmark-β法精度和稳定性的参数取值为γ=0.5,β=0.25.采用瑞利阻尼模拟风力机结构阻尼,刚度阻尼系数取值为0.03[27].风力机结构系统动力激励下的振动方程[28]为

FI(t)+FC(t)+FK(t)

(4)

2 大型水平轴风力机支撑拓扑优化

2.1 BESO算法的数学模型

针对大型水平轴风力机支撑结构优化问题,建立以结构静刚度最大为目标、以结构体积为约束条件的数学力学模型[29],表示为

(5)

(6)

式中:CT为结构平均柔顺度;V*为目标体积比;vi为第i个单元的体积;n为迭代次数;xi为设计变量;xmin和1分别代表某一单元在设计域内缺失或存在.

2.2 单元灵敏度过滤及更新

为避免产生奇异单元刚度矩阵,算法以xmin代替0表示空单元,在该软删单元的BESO算法[30]中所有单元的材料弹性模量可根据材料插值模型设置为

(7)

式中:q为惩罚因子.

单元灵敏度表示增删单元对结构平均柔顺度的影响程度,单元灵敏度大小为单元灵敏度的相对排序,表示为

(8)

式中:Ki和ui分别为单元刚度矩阵和单元位移向量.

为抑制优化过程中出现实体单元与空洞单元交替出现的棋盘格等数值不稳定现象[31],通过权重系数引入节点灵敏度节点概念:

(9)

在迭代过程中,根据节点灵敏度对单元灵敏度进行过滤并更新单元灵敏度,更新过程为

(10)

(11)

2.3 基于BESO算法的结构拓扑优化

鉴于水平轴风力机结构系统对称性,将支撑结构三维优化问题转化为二维平面优化,以风力机塔筒支撑形式和支撑点位置为优化对象.将风力机塔筒主轴高度方向和塔筒半径方向围成的平面作为初始设计域,采用ANSYS软件中的4节点Plane182单元对计算模型进行网格划分.设置底部为固定端约束,选取最不利气动荷载作为外部风荷载.采用反比例型删除率优化模型[32]对风力机支撑结构进行拓扑优化:

(12)

式中:ω为结构材料的删除率;ωmax为给定的最大删除率,取0.057;ωmin为给定的最小删除率,取0.012;Vi为优化过程中每一步迭代过程的体积比;V*值取0.4.

在迭代过程中,根据删除率对体积比进行更新

Vk+1=Vk(1±ω),k=1, 2, 3, …

(13)

式中:Vk为第k次迭代的体积比:Vk+1为第k+1次迭代的体积比.

当迭代结果满足根据体积约束和收敛准则[31]时,终止迭代.为提高迭代效率,设置收敛准则,表示当前迭代步之前10次迭代中平均柔顺度的相对变化量,当目标函数变化随着迭代变化足够小时,则判定迭代结束:

(14)

式中:Ck为第k次迭代的平均柔顺度;τ为收敛因子,取值为0.001.

根据上述计算过程实现基于BESO算法的水平轴风力机支撑结构优化,流程如图4所示.

图4 拓扑优化流程图Fig.4 Flow chart of topology optimization

3 优化前后模型受力特性

3.1 基于BESO算法的支撑结构优化

运用BESO双向渐进结构优化算法,对5 MW大型水平轴风力机塔筒结构(高度为90 m)进行拓扑优化,如图5所示.为防止风力机塔筒顶部发生过大位移,通常对风力机塔筒受弯矩作用较大的底部区域进行加固,所选取的底部拓扑区域为宽度a=10 m、长度b=50 m所围成的平面,以此作为初始设计域Ω,采用4节点Plane182单元进行网格划分.在设计域区域施加外部荷载,包括风力机机舱、转子质量等纵向荷载Fz,风力机叶片运动产生的推力、沿塔筒高度分布的梯度风荷载等在设计域顶部产生的荷载Fxy和弯矩Mxy.底部边界条件为固定端约束,拓扑优化的结构目标体积分数为40%.

图5 塔架结构底部受力及拓扑设计区域Fig.5 Tower structure force diagram and topology design area

分析拓扑过程中体积比即单元删除量变化和平均柔顺度相对迭代次数的变化历程,如图6所示.拓扑迭代初期,拓扑单元较多,单元删除量较大,拓扑迭代速率较快.随着迭代进行,单元体积比逐渐减小,迭代20次后,结构单元体积比开始趋近于目标体积约束条件,拓扑形态开始稳定,最终得到反比例型删除率的优化模型柔顺度量值为1 557.1 N·m.

图6 迭代过程结构平均柔顺度和体积分数变化Fig.6 Average flexibility and volume fraction of optimized structure

对初始设计域,通过上述二维拓扑找形得到带有侧向支撑的塔筒形状如图7(a)所示,其侧向支撑高度为48 m、宽度为8 m.基于支撑结构的对称性,考虑基础为四桩布局,将二维拓扑结果扩展到三维空间,获得的四支撑结构优化形式如图7(b)所示.根据拓扑优化结果,构建高桩承台式水平轴风力机支撑结构模型.

图7 塔筒支撑结构拓扑优化Fig.7 Topology optimization of tower support structures

3.2 优化后风力机支撑结构特性

基于CFD数值方法计算优化后风力机支撑结构的风致水平推力和风力机功率.为考虑最不利荷载工况,在进行CFD数值计算时,设速度入口风速等于切出风速为25 m/s,与边界元法(BEM)、致动盘(AD)[33]及有限体积法(FVM)[34]中的结果进行对比分析,如图8所示.可知,本文风力机功率、风力机旋转时叶片对塔筒的水平推力(T=750 kN)均与相关文献结果接近.

图8 本文结果及其分析对比Fig.8 Power comparison between CFD model and reference

分别计算优化后支撑结构和原塔筒结构的静力特性,注意到风力机塔筒经受周期性脉动风荷载作用,可将风力机旋转多圈并稳定后叶片对塔筒的水平推力用于计算最不利工况下的塔筒位移, 运用式(1)计算沿风机塔筒高度的风荷载.风力机支撑结构顺风向(x向)、横风向(y向)和竖向(z向)最大动位移均发生于塔筒顶部,其中顺风向位移明显大于横风向位移.针对5 MW大型水平轴风力机,分别计算初始支撑结构和优化后支撑结构在额定风速下的塔顶位移.初始结构塔顶顺风向位移为550 mm,竖向位移为-2.5 mm,而优化后结构塔顶的顺风向位移为250 mm,竖向位移为-1.7 mm.优化后支撑结构的位移明显减少,顺风向和竖向位移相对初始支撑结构减少了54.5%和32%,如图9所示.

图9 初始结构与优化结构位移响应云图Fig.9 Comparison of displacement of initial structure and optimized structure

3.3 桩土作用优化后支撑结构风致动力响应分析

考虑桩土作用,分析5 MW风力机初始和优化后结构在顺风向和竖向的动力响应.设来流风速为25 m/s,风力机叶片转速为0.664 rad/s,风力机运动周期为9.463 s.风力机旋转6～10圈后气动性能基本稳定,相应的风力机旋转时长为50～200 s.初始和优化后结构的塔筒顶端顺风向位移(Ux)时程曲线如图10所示.计算发现在t=68 s时刻,风力机结构振幅较大.因此分析该时刻优化前后支撑结构的顺风向和竖向位移,得到初始和优化后塔顶顺风向最大位移分别为 1 130 和300 mm,减少73.45%;竖向位移分别为 -49～-26 mm和-18～-8 mm,最大竖向位移减少了63.27%.由此可见,在保持与初始结构体积相同情况下,优化后风力机支撑结构位移明显减少,结构刚度大于初始结构.

图10 初始结构和优化结构的塔顶顺风向位移时程曲线Fig.10 Comparison of displacement time history curve of two models at the top of tower

考虑桩土作用的初始和优化后结构在50～200 s时段顺风向的速度(vx)响应时程和加速度(ax)响应时程如图11所示.可见优化后结构的速度和加速度响应小于初始结构,风振特性得到显著改善,在风力机旋转多圈如100 s进入稳态后改善效果更明显,如塔顶加速度由初始结构的5.02 m/s2减小为优化后结构的3.8 m/s2,减幅达24%.考虑桩土作用的支撑结构振动效应明显减小,具有较好的减振效果.

图11 初始和优化结构的塔顶顺风向速度时程和加速度时程Fig.11 Comparison of velocity and acceleration time history curve at the top of tower

4 结论

针对大型高桩承台式水平轴风力机支撑结构系统,运用CFD方法计算风致效应,考虑桩土作用效应,运用基于反比例函数的双向渐进结构优化算法对风力机塔筒支撑结构进行拓扑优化,开展风力机结构系统的静动力响应特性分析,得到结论如下:

(1) 运用基于反比例型删除率的BESO算法并考虑桩土作用影响,开展水平轴风力机塔筒支撑结构拓扑优化,得到风力机支撑结构的新形式.

(2) 优化后支撑结构风致效应明显减小,在保持结构体积不变情况下,优化后支撑结构整体刚度明显提高,塔顶顺风向和竖向静力位移分别减小54.5%和32%,动力分析下塔顶在顺风向和竖向最大位移分别减少73.45%和63.27%.

运用BESO算法进行水平轴风力机支撑结构的二维平面优化,后续将对风力机支撑结构开展三维优化研究.