海上风机基础超单元计算方法对比研究

周 润，李 昕，王文华

(大连理工大学 水利工程学院 海岸与近海工程国家重点实验室， 辽宁 大连 116024)

目前海上风机动力响应分析主要采用半整体方法和整体耦合分析方法。半整体方法的研究思路是选取基础结构与风机塔筒连接位置为过渡点，分别建立叶片-轮毂-机舱-塔筒-基础超单元气弹计算模型和基础结构有限元模型，并针对上述模型分别施加风浪荷载。此外，为充分考虑气动荷载对于基础结构影响，半整体方法将气弹模型所得风机荷载作用于基础结构，作用位置为前述所选取法兰点位置。

整体耦合分析方法则是直接建立叶片-轮毂-机舱-塔筒-基础机构的整体模型，同时施加风、浪、流等外部荷载，并考虑相应控制策略，开展整体结构动力反应计算。

海上风机半整体计算方法的核心为海上风机基础超单元计算，即采用相关矩阵缩聚方法在保证基础结构动力特性的前提下减小基础结构特征矩阵维度。目前，结构工程领域常见的Guyan方法[1]、SEREP方法[2]和C-B方法[3]等均已被应用于海上风机基础结构超单元计算，比如Valk等[4]基于2.3MW风机模型开展了C-B缩聚方法与RuBin缩聚方法对比研究，并进一步提出了适用于风机结构缩聚的Dual C-B方法。Voormeeren等[5]则重点验证了C-B方法在海上风机基础结构超单元计算中的适用性，Jason等[6]将C-B方法应用于FAST以建立固定式海上风机整体结构耦合模型，通过C-B方法实现了基础结构顶部节点自由度与塔筒底部节点自由度的耦合。

由上述介绍可知半整体方法虽然能够保证较高的计算精度，不过其缺点仍显著，比如无法充分考虑气动荷载、转子系统与基础结构之间的相互作用，以及气动阻尼对于波浪作用下海上风机结构反应的影响。

随着海上风机研究技术的不断发展和提高，整体耦合模型在海上风机整体结构动力反应分析及安全评价方面得到了广泛应用，特别是针对浮式海上风机结构，比如Luxcey等[7]建立了叶片-轮毂-机舱-塔筒-基础整体耦合模型，施加机电控制策略，系统研究了基础结构运动响应与转子系统之间的耦合效应。方通通[8]分别建立了叶片-轮毂-机舱-塔筒-基础结构的整体耦合模型和叶片-轮毂-机舱-塔筒-超单元的半整体模型，开展了随机风场作用下固定式海上风机基础结构疲劳分析，验证得出环境荷载、控制系统与结构反应之间的耦合效应对于基础结构疲劳具有显著影响，有必要直接基于整体耦合模型开展海上风机结构疲劳计算。

结合海上风机整体耦合模型研究进展可知，环境荷载、控制系统及支撑系统之间的耦合效应不可忽视，对于海上风机结构动力反应计算有必要建立整体耦合计算模型。不过需指出，现阶段海上风机整体耦合模型中所包含基础结构有限元模型均采用线性梁单元建立[9]，所以对于复杂基础结构仍无法直接建立相应整体结构模型。

因此，对于复杂基础结构海上风机仍需采用半整体方法开展基础结构设计及承载力校核。此外，相对于整体方法，在保证计算精度的前提下，采用半整体方法将极大地提高计算效率。海上风机基础超单元作为半整体方法计算核心，基础超单元选取是否合理与气弹模型计算精度密切相关，并将进一步影响基础结构设计。故本文拟选取不同基础结构类型展开海上风机整体耦合与半整体分析方法对比，并进一步研究不同凝聚方法对于海上风机基础超单元及气弹模型计算精度的影响。

1 凝聚方法理论

1.1 C-B凝聚方法

随机荷载作用下通用基础运动方程可写为

(1)

基于C-B方法开展基础结构特征矩阵缩聚，首先将结构节点区分边界节点和内节点[3]，对应的节点位移向量分别表示为{δB}和{δI}。进一步，整理结构运动方程可得

(2)

式中：上标B和I分别表示边界节点和内节点；{FB}和{FI}分别为作用于边界节点和内节点的外荷载向量；[KBB]和[KII]分别为边界节点刚度矩阵和内节点刚度矩阵；[KBI]和[KIB]分别为位于重排列刚度矩阵副对角线的刚度矩阵。[MBB]和[MII]分别为边界节点和内节点的质量矩阵。

基于文献[3]可知，按照式(2)形式改写运动方程之后，内节点位移将主要由两部分组成，分别为：

(1) 由边界节点位移所造成的内部节点位移，将式(2)中内部节点外荷载向量假定为零，则可得由边界点位移所造成的内部节点位移为

(3)

(2) 由相邻节点位移所造成的位移反应为

(4)

设凝聚转换矩阵为[α]，可得转换公式为

(5)

根据式(3)、式(4)，可得凝聚转换矩阵为

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

1.2 SEREP凝聚法

将结构位移向量{x}用广义模态坐标{qr}∈{Rr}表示，

{x}=[φ]{qr}

(14)

式中：[φ]为结构模态矩阵。

选取结构位移向量{x}的主、副自由度，并同时选取模态矩阵[φ]各阶模态主自由度，进一步可将式(14)改写为，

(15)

式中：{xm}和[φm]分别为结构主自由度位移及其对应的模态矩阵；{xs}和[φs]分别为结构副自由度位移及其对应的模态矩阵。

由式(15)可得，结构主自由度位移向量{xm}与广义模态坐标{qr}的换算关系为，

{xm}=[φm]{qr}

(16)

进一步得到qr为，

(17)

由此，可建立结构位移向量与所选取主自由度位移向量的换算关系，

(18)

将SEREP方法所得缩聚矩阵及其转置分别右乘和左乘结构初始刚度和特征矩阵，可得凝聚后的结构刚度和质量矩阵为，

(19)

(20)

需指出，基于式(15)开展缩聚矩阵计算时，需涉及到将结构整体模态矩阵运算。当结构自由度数量较为庞大时，计算量将十分繁巨。因此，在指定基础结构主副自由度之前，可优先选取结构模态矩阵主模态。基于结构主模态矩阵展开结构特征矩阵缩聚，在保证计算精度的前提下，可以显著降低计算量[2]。

2 样本风机计算模型

2.1 半整体模型运动方程

2.1.1 上部结构运动方程

基于结构动力学[10-11]基本理论可知，空气动力荷载作用下，法兰盘以上海上风机结构运动方程可表示为：

(21)

2.1.2 基础结构运动方程

水动力荷载作用下基础结构运动方程可表示为：

[Kb]*{ub}={fTP}+{fhydro}

(22)

2.2 整体耦合模型运动方程

基于结构动力学[10-11]、空气动力学[12]和水动力学[13]基本理论可得，风和波浪荷载联合作用下海上风力机整体结构运动方程为：

(23)

2.3 样本风机基本参数

2.3.1 基础结构基本参数

本文拟选取我国某海域5 MW单桩基础和五桩基础结构海上风机开展整体与半整体模型及半整体模型超单元计算方法对比，所选取单桩基础和五桩基础结构基本参数如表1所示，基于ANSYS[14]所建立的基础结构有限元模型如图1所示。

表1 基础结构及有限元模型基本参数

图1单桩和五桩基础结构有限元模型

2.3.2 风机基本参数

采用NREL 5MW基准风机[15]作为本文的样本风机，NREL 5MW风机叶片数目为3，轮毂中心高度为90 m。风机塔筒高度为78.6 m，塔筒顶部和底部直径分别6 m和3.76 m，NREL 5MW风机详细参数可参考文献[15]。

3 基础结构超单元计算

本文拟分别基于C-B方法和SEREP方法生成单桩基础和五桩基础结构超单元矩阵，主要计算流程如图2所示。

图2超单元生成过程流程图

在基于上述方法开展基础结构超单元计算时，需注意C-B方法和SEREP方法的如下不同：

(1) SEREP凝聚法在凝聚时涉及到主自由度的选择，本文在基于该方法进行基础结构凝聚时将塔筒顶部节点x向自由度、离塔筒顶部最近的5个节点的x向自由度共同选取为单桩基础和五桩基础结构自由度。

(2) C-B凝聚法需区分内外节点。本文中将单桩基础和五桩基础的外节点设置为基础结构边界节点和过渡节点，剩余节点为内节点。

4 样本风机动力特性对比

4.1 基础结构特征频率

基于单桩基础和五桩基础结构总刚度和质量矩阵，分别采用C-B方法和SEREP方法缩聚得到的基础结构特征矩阵，并进行动力特性计算和对比，相关对比如表2、表3所示。

表2 单桩基础结构前六阶频率的对比

表3 五桩基础结构前六阶频率对比

从基础结构的频率对比来看，基于SEREP法凝聚后计算的频率前6阶与结构频率完全一致，C-B法计算的频率与整体稍有误差，最大相对误差为1.03%，仍满足计算精度要求。由理论介绍可知，造成上述计算误差的原因在于SEREP基于基础结构全部自由度开展动力特性计算，而C-B法则只考虑了内部节点自由度影响。

4.2 转子-机舱-基础超单元特征频率

基于FAST建立转子-机舱-塔筒-基础超单元的上部结构有限元模型，假定风机处于停机状态，在塔筒顶部施加3 m初始位移，开展自由衰减测试，进一步验证采用基础超单元代替基础结构有限模型建立上部结构有限元模型开展气动分析的可行性，对比结果如表4所示。

由表4可知，对于单桩基础，基于转子-机舱-塔筒-超单元模型得到的整体结构基频与整体耦合模型计算结果基本一致；而对于五桩基础风机，相对于SEREP方法，采用C-B方法得到的结构基频更为接近整体耦合模型计算结果，此时相对误差仅为3.18%。

表4 海上风机整体与半整体模型动力特性对比

综上，对于单桩和多桩基础固定式海上风机，采用等效基础超单元代替基础结构有限元模型，即将基础结构等效为弹性边界条件可以得到较为准确的整体结构基频，同时所对比的两种超单元计算方法均取得了较高的计算精度。

5 耦合分析响应对比

基于FAST分别建立海上风机转子-机舱-塔筒-基础结构整体耦合模型和转子-机舱-塔筒-基础超单元半整体计算模型开展脉动风速作用下海上风机整体和半整体计算模型对比，相应脉动风速时程如图3所示。

图38m/s脉动风风速时程曲线

5.1 单桩基础风机计算

脉动风场作用下，分别基于整体模型、半整体模型(SEREP)和半整体模型(C-B)得到的单桩基础风机塔筒顶部YBFX、底部风机荷载TBMY和塔筒顶部位移YBDX时程统计值对比如表5所示，塔筒顶部位移YBDX动力响应局部时程图及频域反应图如图5、图6所示。

图4 风机结构坐标系示意图

图5 单桩基础结构塔筒顶部位移

图6单桩基础塔筒顶部位移频域反应

由风机荷载和塔筒位移时频域反应对比可知，基于上述整体和半整体模型得到的海上风机结构反应基本一致。时程反应均值、标准差和协方差对比如表5所示，统计值相关计算公式如式(26)—式(28)所示。

(26)

(27)

(28)

式中：xi为半整体法计算的每时刻响应结果；yi为整体法计算的每时刻响应结果；n为样本数量。

表5 单桩基础风机动力响应统计值

由海上风机结构时程反应均值可得，稳态风场作用下整体模型与半整体模型的计算结果基本一致。进一步综合协方差对比可得，相对于SEREP方法，基于C-B方法所建立的半整体模型得到的海上风机结构反应更为精确，如基于半整体模型(C-B)和半整体模型(SEREP)所得风机荷载TBMY的协方差分别约为1 306 kN·m和1 405 kN·m。

5.2 五桩基础风机计算

基于平均风速为8 m/s脉动风场，如图3所示，开展五桩基础海上风机整体模型与半整体模型及半整体模型超单元计算方法对比，该基础结构相对于单桩基础的区别为显著增加基础结构自由度数目，约为1 224。随机风场作用下基于五桩基础海上风机整体和半整体结构动力响应时频域反应对比如图7和图8所示，统计值对比分别表6所示。

图7 五桩基础结构塔筒顶部位移

图8五桩基础塔筒顶部位移频域反应

由塔筒顶部位移时程反应对比可知，随着基础结构自由度数量的增加，半整体计算模型仍能保持较高的计算精度，采用分别C-B方法和SEREP方法建立的半整体模型计算结果基本一致。同时基于上述整体和半整体模型所得塔筒顶部位移频域反应控制频率基本相同，如图8所示。

由表6所列风机结构时程反应统计值对比可知，分别基于半整体(C-B)和半整体(SEREP)模型得到的塔筒顶部位移DX和塔筒顶部风机荷载FX时程反应统计值基本一致，比如基于上述半整体模型分别得到的塔筒顶部位移极大值分别为0.1553 m和0.1554 m，与整体耦合模型计算误差分别为2.88%和2.81%。

表6 五桩基础风机动力响应统计值

进一步对比可得，基于半整体模型(C-B)所得时程反应协方差明显小于半整体模型(SEREP)计算的结果，比如基于上述两个计算模型所得塔筒底部弯矩协方差分别为1 720.6 kN·m和3 435.7 kN·m。上述显著差异说明，基于C-B方法所得基础超单元建立的半整体模型计算结果更为符合海上风机整体耦合模型得到的结构动力响应。

6 结 论

本文基于单桩基础和五桩基础海上风机系统开展了整体耦合与半整体计算模型及半整体模型超单元计算方法对比研究。基于所述对比可得如下结论：

(1) 对于固定式海上风机，在保证基础超单元选取合理的前提下，基于转子-机舱-塔筒-基础超单元所得结构响应与整体耦合模型计算结果基本一致。

(2) 基于C-B方法和SEREP方法计算固定式基础结构超单元，在保证基础结构低阶动力特性相似的前提下，均能够准确模拟随机风荷载作用下转子系统和塔筒结构动力响应。

(3) 由时程反应协方差对比可知，相对于SEREP方法，基于C-B方法所建立的半整体计算模型计算结果更为接近整体耦合模型计算结果。