海上风力发电结构动力响应分析方法与流程

本发明属于海上风力发电，涉及风、浪等荷载作用下海上风机整体结构动力响应的计算分析方法，尤其涉及一种海上风力发电结构动力响应分析方法。

背景技术：

1、永不枯竭的风能是大自然的恩赐，相比陆上风电，海上风力发电具有占地少、风速大、风向稳等优点，正成为新能源领域开发的重点。我国海上风能资源丰富、开发效益高，是经济可行的绿色环保能源。目前海上风力发电多采用水平轴风机系统，风电机组主要由叶片、轮毂、机舱、塔架和水下基础等部分组成，设计、建造出具有足够刚度和强度的风机基础是维持风机在使用期限内稳定运转、创造经济价值的必要条件。

2、综上，本发明建立风机气动力-水动力-结构动力响应-土抗力的一体化耦合计算模型，提出海上风机整体结构动力特性的研究分析方法，能够为海上风机结构设计、施工及运行维护提供技术支撑(海上风机基础结构的研究现状与发展前景，张济，何海华，张兆德)。

技术实现思路

1、针对现有技术存在的问题，本发明提供一种海上风机在风、浪联合作用下整体结构动力响应分析方法，该分析方法能够处理流体-海上风机结构-海床地基相互作用问题，在大幅度提高计算效率的基础上，保证风机整体结构动力响应计算结果的准确性。

2、为了达到上述目的，本发明采用的技术方案为：

3、一种海上风力发电结构动力响应分析方法，是一种海上风机在风、浪荷载下的整体结构动力响应预测分析方法，该分析方法基于海上风机整体结构动力响应耦合计算模型实现，其计算模型包括四个模块：空气动力学模块、水动力学模块、桩土相互作用模块、海上风机整体结构动力响应计算模块，其中，水动力学模块也为波浪荷载计算模块。从图1中可以看出，整体耦合计算模型中的三个荷载计算子模块(空气动力学模块、水动力学模块、桩土相互作用模块)之间通过海上风机整体结构动力响应计算模块进行联系。首先，分别利用三个荷载计算子模块计算波浪荷载、风荷载、土荷载，并转化为海上风机整体结构动力响应计算模块中的等效节点荷载，随后利用海上风机整体结构动力响应计算模块求解风机整体结构动力响应得到各节点处的位移、速度和加速度，并分别将其传递给对应的荷载计算模块。因此，本发明建立的模型在时域上是相互耦合的，能够充分考虑各部分外荷载与海上风机结构的相互作用。具体包括以下步骤：

4、第一步，通过建立高阶谱数值波浪水槽模型模拟非线性极端波浪场

5、真实的海洋波浪具有随机性、非线性，数值分析中通常利用不规则波浪进行模拟。虽然在线性理论中，不规则波可表示为不同波浪成分的线性叠加；但当波浪非线性较强时，线性理论模拟的结果不能满足精度要求，本发明需要建立完全非线性高阶谱数值波浪水槽模型模拟非线性极端波浪。

6、所述的高阶谱数值波浪水槽模型首先将水槽内流体的总速度势分解为周期速度势和附加速度势，利用周期速度势处理周期性边界条件，并利用附加速度势处理非周期性造波边界条件，从而实现数值水槽的造波功能；其次，为了提高计算效率，采用拟谱方法求解分解后的两个速度势，将二者分别展开为满足控制方程和相应边界条件的傅里叶级数形式；再次，每个计算时刻根据水槽造波边界条件计算附加速度势；最后，将求得的附加速度势代入到非线性自由水面边界条件，并进行时间积分得到下一时刻的水槽波面及波面处的周期速度势，从而构建出完整的高阶谱数值波浪水槽模型，周而复始，直至计算结束。本发明构建的高阶谱数值波浪水槽模型可以将理论计算的目标波浪经非线性传播演化成更加复杂和真实的非线性波浪，即能够用于高效、精确地模拟随机性、强非线性的真实海上波浪，也就是实现高阶谱数值波浪水槽模型。具体步骤如下：

7、1)构建高阶谱数值波浪水槽模型，将其水槽域定义为计算域，如图2所示，将计算域内流体的总速度势分解为两部分，分别为周期速度势和附加速度势，具体的：

8、φ＝φp+φadd   (1)

9、所述的周期速度势为φp，其需要满足自由水面边界条件、水底边界条件、左右两端不透水边界条件：

10、

11、式中，x表示水平坐标，表示对x求偏导数；z表示垂向坐标，表示对z求偏导数；lx为水槽长度，d为水槽水深。

12、所述的附加速度势为φadd，其需要满足水槽左端的造波边界条件、水底边界条件、右端不透水侧壁边界条件：

13、

14、式中，u为在水槽入口处输入的目标波浪的水平速度，所述u值根据目标波浪参数利用二阶stokes波浪理论直接计算得到，具体方法参见：水波理论及其应用.邹志利.科学出版社.2005年，第二章。所述的目标波浪参数包括波浪周期、波高以及水深，由上述参数可求得目标波浪的水平速度。

15、2)求解附加速度势φadd及附加速度势的时间导数

16、2.1)确定附加速度势φadd

17、φadd在扩展的计算域dadd中求解，如图3所示，扩展域的总宽度lz＝hadd+d，其中hadd表示延拓计算域上边界垂向坐标，d表示水深；输入边界的水平速度关于中心线zc＝(hadd-d)/2反对称分布，因此本发明选择hadd＝3d，中心线高度zc＝d，保证中心线zc始终高于自由水面。

18、为了精确模拟更陡峭的非线性波，本发明提出了从水槽底部z＝-d到瞬时波面z＝η输入目标波浪水平速度的造波方法。具体的：

19、水槽入口处波面η(0,t)可由目标波浪参数根据二阶stokes波浪理论直接计算得到(水波理论及其应用.邹志利.科学出版社.2005年，第二章)。在扩展域dadd中根据第一步公式(3)中所述的φadd满足的边界条件，将其按照特征函数展开为如下形式：

20、

21、式中，bn(t)是附加速度势φadd第n个模态的振幅；κn＝(2n-1)π/(hadd+d)是扩展域dadd垂直方向上的第n个本征模态，hadd表示延拓计算域上边界处的垂向坐标，nz表示垂直方向的总模态数量；lx表示水槽长度。则根据附加速度势φadd的特征展开式可直接求取其x方向偏导数

22、

23、将偏导数代入水槽入射造波边界条件：

24、在x＝0   (6)

25、式中，u(z)表示在水槽入口处输入的目标波浪的水平速度，其根据目标波浪参数利用二阶stokes波浪理论直接计算得到(水波理论及其应用.邹志利.科学出版社，2005，第二章)，利用垂向特征函数cos[κn(z+d)]的正交性，根据造波边界处水平速度u(z)的垂向分布，即可确定系数bn(t)，进而确定φadd。

26、2.2)确定附加速度势的时间导数

27、所述的附加速度势的时间导数满足以下边界条件：

28、在x＝0   (7)

29、式中，ut(z)表示在水槽入口处输入的目标波浪的水平加速度，其根据目标波浪参数利用二阶stokes波浪理论直接计算得到(水波理论及其应用.邹志利.科学出版社，2005，第二章)。

30、则附加速度势时间导数按照特征函数展开为如下形式：

31、

32、则附加速度势时间导数的x方向偏导数可解析表达成如下形式：

33、

34、式中，cn(t)是第n个模态的振幅；将偏导数代入式(7)，利用垂向特征函数cos[κn(z+d)]的正交性，根据造波边界处水平速度时间导数ut(z)的垂向分布，确定系数cn(t)，进而确定

35、为了避免初始效应，u(z)和ut(z)需要乘以缓冲函数rm(t)，表达式如下：

36、

37、式中，tm为缓冲时长，本发明取为两倍的波浪特征周期。

38、3)求解周期速度势φp

39、3.1)将水槽的自由水面处的周期性速度势函数φp表示成仅依赖于水平坐标的形式：

40、

41、式中，t为当下计算时刻，η为t时刻的波面高度。

42、3.2)则水槽的完全非线性自由水面边界条件表达式如下：

43、在z＝η上   (12)

44、在z＝η上   (13)

45、式中，g表示重力加速度；w表示自由水面处水质点垂向速度；表示对时间求偏导数；ν表示在水槽末端加入的消波区空间阻尼参数，加入消波区后能有效防止波浪反射，本发明选取三次多项式形式的空间阻尼分布函数：

46、

47、式中，x0表示消波区起点，消波区域长度为lx-x0，则无量纲长度系数μ＝(x-x0)/(lx-x0)；α为阻尼强度；

48、3.3)周期速度势φp对应的水质点垂向速度通过下式递推计算得到：

49、

50、式中，m是非线性保留总阶数，w(m)表示第m阶分量，j为整数；

51、根据第一步公式(2)中所述的周期速度势φp满足的边界条件，将其展开成傅里叶级数，并按下式递推计算得到：

52、

53、式中，表示第m阶周期速度势第n个模态的振幅，kn＝nπ/lx为波浪水槽的第n个自然本征模态，nx为水平方向的总模态数量。

54、最终构建得到高阶谱数值波浪水槽模型，即高阶谱数值波浪水槽模型由公式(1)、(4)、(8)、(12)、(13)及(16)构成。

55、4)将构建得到的高阶谱数值波浪水槽模型用于高效、精确地模拟随机性、强非线性的真实海上波浪，如图4所示，包括以下步骤：

56、4.1)：通过构建的数值波浪水槽模型，划分水槽在水平向和垂向的网格，并给定水槽的初始条件。设置数值波浪水槽模型在初始时刻各处的自由表面势和自由表面高程η如下：

57、η＝0，在t＝0   (17)

58、4.2)：在当前t时刻，根据水槽造波边界条件公式(6)、(7)求解得到附加速度势φadd和具体包含如下子步骤：

59、4.2.1)根据目标波浪参数，利用二阶stokes波浪理论计算目标波浪在水槽入口处的波面η、水质点速度u和加速度ut；

60、4.2.2)根据水槽造波边界条件式(6)、(7)，在扩展域中计算附加速度势φadd、附加速度势时间导数

61、4.2.3)在得到附加速度势φadd特征展开系数的基础上，采用解析的方法求解自由水面处的附加速度势φadd的空间导数

62、4.3)：计算周期速度势φp及其空间导数、以及波面η的空间导数，并计算流场中水质点速度与加速度，具体包含如下子步骤：

63、4.3.1)已知当前时刻的根据递推关系式(16)进行递推计算最终进行累加求和可求得该时刻的

64、4.3.2)分别计算波面处周期速度势和波面的水平导数以及周期速度势φp的垂向导数w；

65、4.3.3)在计算得到周期速度势φp、附加速度势φadd的基础上，利用解析方法进一步计算二者的空间、时间导数并求和，即可得到流场内部任意位置处的水质点的速度、加速度。

66、4.4)：将步骤4.2)、4.3)计算得到的物理量作为强迫项代入到完全非线性自由水面条件式(12)、(13)中进行时间积分，得到下一时刻的自由表面势和自由表面高程η。判断计算时间是否结束，如未结束周而复始重复步骤4.2)～4.4)，进入下一时刻t+δt的计算；

67、4.5)：当到达波浪模拟结束时间，则高阶谱数值波浪水槽模型计算结束，即可实现数值波浪水槽对非线性波浪生成与传播的全过程数值模拟。

68、本发明采用高阶谱数值波浪水槽模型生成聚焦波模拟海上的非线性极端波浪，如图2所示，非线性数值波浪水槽的入口速度根据模拟的目标波浪参数采用二阶stokes波浪理论计算，设置聚焦波各个波浪成分均在时刻t＝t0、空间位置x＝x0处波面高程达到最大。各个波浪成分在聚焦波的频率范围内按照等间隔分布，波幅根据newwave模型由波浪频谱确定：

69、

70、式中，ai为第i个波浪成分对应的线性聚焦波幅；n为波浪成分总个数；a为总聚焦波幅；s(fi)为第i个波浪成分对应的频谱值；s(fn)为第n个波浪成分对应的频谱值。其中，s(f)为波浪的频谱，本发明采用jonswap谱：

71、

72、式中，h1/3为有效波高；tp为谱峰周期；fp为谱峰频率；γ为频谱的谱峰升高因子，

73、

74、

75、传统不规则波二阶stokes波浪理论(水波理论及其应用.邹志利.科学出版社.2005年，第二章)仅能计算出聚焦波的二阶和频与差频成分，而本发明建立的高阶谱数值水槽模型可以模拟强非线性聚焦波中的更高阶的成分，计算结果更精确，能够模拟海上真实的极端波浪。

76、第二步，构建波浪荷载计算模块，每个时刻根据rainey细长体理论计算第一步模拟的极端波浪作用在风机结构上的波浪荷载，其分为两部分：线波浪荷载frainey(t,z)贡献部分、以及水面集中力(fx)surface贡献部分。具体如下：

77、(1)计算线波浪荷载frainey(t,z)贡献部分

78、本发明采用rainey细长体理论直接计算极端波浪作用在桩基中心线任意高度处的波浪力，即风机结构单位长度上的线波浪荷载frainey(t,z)：

79、

80、式中，ρ为海水密度；u、分别为入射波浪场中桩柱中心轴线z高度处的水质点水平速度、加速度；wz表示水质点垂向速度在z方向上的偏导数。x、分别表示风机结构在z高度处的位移、速度、加速度；对于聚焦波入射情形，本发明利用高阶谱数值波浪水槽模型模拟的聚焦波代表极端波浪，在每个时间步还原水槽内波浪场并计算流场中风机结构处的水质点速度、加速度。cm为附加惯性力系数或质量系数；cd为阻力系数或者速度力系数；a＝πd2/4为截面面积，d为桩柱直径。值得注意的是，与经典的morison公式相比，本发明采用的rainey细长体理论额外考虑了轴向散度项

81、所述的线波浪荷载frainey(t,z)贡献部分的波浪荷载需将公式(22)计算得到的frainey(t,z)沿风机结构从海底z＝-d到瞬时波面z＝η(t)进行垂向积分，表达式如下：

82、

83、式中，z表示垂向坐标；d表示水深；η(t)表示当下计算时刻风机轴线处的波面；fx表示frainey(t,z)波浪作用在风机结构上的水平力；my表示波浪作用在风机结构上的绕泥面处的倾覆力矩。

84、(2)水面集中力贡献部分

85、所述的由水面集中力(fx)surface贡献部分的波浪荷载计算公式如下：

86、

87、式中，(my)surface表示水面集中力(fx)surface作用在风机结构上绕泥面处的倾覆力矩；

88、将式(23)、(24)计算得到的两部分荷载累加求和即可得到每个时刻波浪作用在风机结构上的总波浪荷载。

89、第三步，构建空气动力荷载计算模块，每个时刻根据叶素动量理论计算风机结构上的风荷载，具体如下：

90、根据风机上部结构的叶片参数，本发明利用经典的叶素动量理论计算风机气动荷载。该理论根据动量理论与叶素理论计算的推力和转矩相等，迭代求解得到轴向和切向诱导因子a和a'，如图5所示，风机叶片整体所受的水平推力t通过沿叶片径向积分得到：

91、

92、式中，v0为风速，ρair为空气密度，rhub为轮毂半径，r为风机叶片半径，f为损失因子。

93、第四步，构建桩土相互作用模块，每个时刻利用p-y曲线法计算风机结构上的土荷载，具体如下：

94、海上风机结构插入砂土地基一定深度进行固定，在风、浪荷载作用下风机结构会产生侧向挠曲变形，一部分侧向荷载由桩身承担，另一部分荷载通过桩基传递给土体。

95、本发明利用p-y曲线法计算土荷载，对于砂土地基建立的单位桩长上土体抗力p与桩身侧向位移y的函数关系如下：

96、

97、式中，as为与荷载类型相关的系数，由下式确定：

98、

99、k为初始地基反力模量，与土体内摩擦角有关；z为土体深度，由泥面往下为正；pu表示任意深度z处土体的极限承载力，按下式计算：

100、

101、式中，d为单桩基础直径；γ为土的有效容重；zr＝((c3-c2)d)/c1为临界深度；系数c1、c2和c3为土体内摩擦角φ'的函数，由下式确定：

102、

103、图6给出了土体抗力计算示意图，对地基中的土抗力沿垂向积分可得到风机结构上的总抗力的表达式如下：

104、

105、式中，lp为桩基入土部分的长度，d为水深。

106、第五步，构建风机整体结构动力响应有限元求解模型，计算每个时刻海上风机整体结构动力响应，具体如下：

107、5.1)建立海上风机整体结构动力响应有限元求解模型

108、海上风机在静力状态下结构竖直，长度为l、桩基的线密度为风机材料的弹性模量为e，横截面惯性矩是i，截面抗弯刚度为ei，将风机叶片与机舱等上部结构视为集中质量。由于风机结构横纵比小，采用有限元方法应用euler-bernoulli梁单元(有限单元法[m].王勖成.北京:清华大学出版社,2003.)离散风机结构，建立风机结构线弹性有限元模型。分布参数体系满足如下的运动微分方程：

109、

110、式中，x、分别表示风机结构在z高度处的位移、速度、加速度，c为结构阻尼，f(t,z)为风机结构单位长度所受的波浪荷载。海上风机顶部和底部均无约束，满足的边界条件如下：

111、

112、海上风机结构采用标准的四自由度梁单元进行离散，如图7所示，单元中每个节点有两个自由度(节点位移、转角)。应用galerkin有限单元法(王勖成.有限单元法基本原理和数值方法[m].北京:清华大学出版社,1997.)，将权函数取为插值函数，建立加权余量的积分方程并求解。对每个单元进行积分可以得到单元质量阵、单元刚度阵：

113、

114、将风机单位长度所受的波浪荷载f(t,z)转换为单元等效节点荷载，转化过程如图8所示，并将单元矩阵组装到整体矩阵，即可获得风机系统整体结构动力响应的时域微分方程如下：

115、

116、式中，x、为有限元模型中节点位移、速度、加速度组成的列向量；mg、cg和kg分别是整体质量、阻尼和刚度矩阵，结构阻尼矩阵cg采用瑞利阻尼矩阵形式，即阻尼阵由刚度阵和质量阵线性叠加得到：

117、cg＝αmg+βkg   (35)

118、式中，阻尼系数α＝2ξω1ω2/(ω1+ω2)，β＝2ξ/(ω1+ω2)，ξ是阻尼比，ω1和ω2是风机结构前两阶振动频率。ftotal为每个时刻将前述第二步、第三步、第四步计算得到的波浪荷载、空气动力荷载以及土荷载，将其转化为风机有限元模型的等效节点荷载，转化过程如图8所示，并进行组装形成的荷载列向量。

119、5.2)通过逐步积分法数值求解海上风机整体结构动力响应

120、逐步积分法的基本思想是将风机整体结构动力响应时域微分方程式(34)在时域离散，转化为对时间的差分格式，根据前一时刻t－δt的响应量(位移、速度和加速度)，以及t时刻的结构荷载(注：由于风机结构响应相对水质点运动、风速较小，因此波浪荷载、风荷载的计算和风机整体结构动力响应的求解采用弱耦合形式，即波浪荷载计算公式中的风机结构动力响应量仍用上一时刻t－δt的响应结果。)，采用高精度的数值积分方法逐步求出t时刻的响应量。本发明采用经典的newmark方法(newmark n m.a method ofcomputation for structural dynamics[j].journal of engineering mechanics-asce,1959,85(3):67-94.)完成时间积分，即可求解t时刻的各节点处的位移、速度和加速度，并将其传递给对应的荷载计算模块，周而复始，直到计算结束。

121、5.3)输出每个时刻的风机整体结构动力响应，并对结果进行分析找到极值响应及其对应的时间，为海上风机实际工程设计提供高效、精确的分析手段。

122、本发明的创新点分析：

123、本发明提出并建立海上风机气动力-水动力-结构动力响应-土抗力的一体化耦合计算模型，实现海上风机整体结构动力响应的数值分析，从而为单桩式海上风机的实际工程设计提供了数值计算方法。

124、本发明的有益效果：

125、(1)本发明提出的算法可以充分考虑风、浪等联合作用下海上风机塔架振动与气动作用的耦合效应，可以准确计算风机整体结构动力响应；

126、(2)本发明开发的海上风机动力响应全耦合计算模型具有计算高效、稳定的特点，能够充分考虑非线性“波-波”相互作用的极端波浪与单桩式风机相互作用，对于海上风机强非线性、瞬态性的结构高频共振现象提供了有效的分析手段。