风力机柔塔系统固有振动特性快速建模及仿真方法

1.本发明涉及多体系统动力学技术领域，特别涉及一种风力机柔塔系统。


背景技术：

2.随着风能产业的快速发展，为了更有效地捕获风能，大型化的兆瓦级风力机正逐步成为主流。近年来，兆瓦级风力机主要应用高柔性塔筒支撑。柔性塔可支撑风力机轮毂到110～150m的高度，成本较低。塔架一般为薄壁钢空心圆筒，具有柔性大、阻尼小、重量轻的特点。但是风力机柔塔引起的流体
‑
结构相互作用(fsi)振动更加突出。风力机的生命周期主要受其动力学响应影响，为了避免风力机在流速范围内产生共振，必须知道其固有频率。而风力机的低阶固有频率主要受塔筒影响，因此预测风力机整机振动特性对研究风力机塔筒结构稳定性具有重要意义，为进一步研究塔筒减振提供技术支持。另外，风力机塔筒承受着复合荷载，主要包括风荷载、自身重量及其引起的偏心矩。这些载荷可导致机电设备的运行失效，也会显著增加塔架上的变形和应力。柔塔作为风力机的重要承载部件，其变形和振动会降低结构强度和稳定。在严重情况下，可能会导致疲劳失效，甚至塔架倒塌。因此，风力机柔塔的振动特性预估对于确定结构可靠性具有重要意义，同时也需要一种风力柔塔系统建模和振动特性快速计算方法。
3.风力机可归结为由多个刚体(机舱、轮毂等)和柔性体(叶片和塔架)组成的多刚柔体系统。有限元法(fem)是描述柔性体变形的常用方法。已有文献对风力机塔筒进行了有限元建模，但在其模型中，刚体往往被简化为集中质量。有限元法虽然具有较强的通用性，但不能单独完成多刚柔体系统的动力学建模。此外，构造高阶连续场函数难度大，计算效率低。多体系统传递矩阵法(mstmm)是芮筱亭院士提出的一种新的多体动力学方法，已应用于多种多刚柔体系统的动力学建模中。该方法可用于计算多体系统的振动特性和动力响应。mstmm具有矩阵阶低、公式简单、计算成本低等优点。mstmm在风力机柔塔中的应用，既能考虑柔性塔的几何非线性，又能考虑柔体变形与刚体运动的耦合。因此，该方法适用于风力机柔塔系统的建模和振动特性分析。


技术实现要素：

4.针对现有技术中存在的不足，本发明提供了一种风力机柔塔系统固有振动特性快速建模及仿真方法，实现了风力机柔塔振动特性的快速建模和仿真，方便了风力机柔塔系统的振动特性分析。
5.本发明的目的是这样实现的：一种风力机柔塔系统固有振动特性快速建模方法，包括以下过程：
6.步骤一：对风力机柔塔系统进行简化，将系统离散成各多体动力学元件；
7.步骤二：基于多体系统传递矩阵法推导考虑重力作用的塔筒梁传递矩阵；
8.步骤三：确定其他多体动力学元件的传递矩阵；
9.步骤四：确定边界条件，建立柔塔系统总传递方程；
10.作为本发明的进一步限定，步骤一具体包括：忽略影响塔筒振动特性的次要部分，对风力机柔塔结构模型做出如下假设和简化：
11.(1)忽略塔筒模型扭转运动；
12.(2)假设风力机叶片、轮毂和机舱为一个存在质量偏心的空间刚体，忽略其中复杂部件的影响；
13.(3)忽略塔筒自身的结构阻尼；
14.(4)根据塔架的高质量比特性，忽略塔筒内部空气质量；
15.(5)假设塔筒内部法兰为集中质量，忽略螺栓、螺丝和螺孔造成的影响；
16.(6)假设塔筒截面始终保持圆形，不考虑壳体截面变形；
17.(7)假设塔筒表面平滑，不存在焊缝凹陷；
18.(8)假设塔筒为梁结构，假设不存在剪切效应和局部屈曲效应；
19.(9)假设塔筒与机舱为固定连接，忽略机舱相对于塔筒的扭转；
20.(10)忽略风力机内部爬梯、电缆架、平台等对塔架整体强度影响较小的复杂部件；
21.根据风力机柔塔系统的结构，将塔筒离散成多段考虑顶部构件压力作用、不同结构参数的euler
‑
bernoulli梁；将风力机塔筒内部的法兰、螺栓和支撑其的平台处理为集中质量；将作用在塔筒顶部的风力机叶片、轮毂和机舱处理为一个存在质量偏心的空间刚体；
22.根据模型中所有元件的运动和相互作用，状态矢量可统一为z＝[y,z,θ
y
,θ
z
,m
y
,m
z
,q
y
,q
z
]
t
，其中，z为沿坐标轴z位移对应的模态坐标列阵，y为沿坐标轴y位移对应的模态坐标列阵，θ
z
为该点处相对于平衡位置相对于在z轴的角位移对应的模态坐标列阵，θ
y
为该点处相对于平衡位置相对于y轴的角位移对应的模态坐标列阵，m
z
为沿坐标轴z内力矩对应的模态坐标列阵，m
x
为沿坐标轴x内力矩对应的模态坐标列阵，q
z
为沿坐标轴z内力对应的模态坐标列阵，q
y
为沿坐标轴y内力对应的模态坐标列阵。
[0023]
作为本发明的进一步限定，步骤二中所述推导塔筒的梁传递矩阵包括：
[0024]
建立考虑重力的双自由度euler
‑
bernoulli梁的振动微分方程：
[0025][0026][0027]
式中，ei为梁的弯曲刚度，为单位长度质量，t为时间，z为塔筒流向位移，y为塔筒横向位移，x表示在x轴上的位置，g为梁轴向所受的压力；
[0028]
由模态坐标转换，令
[0029]
y(x,t)＝y(x)e
iωt
,z(x,t)＝z(x)e
iωt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0030]
式中，ω为角速度，y(x)、z(x)为y(x,t)、z(x,t)通过模态坐标转换得到，i代表虚数；
[0031]
将式(3)其代入式(1)、(2)，使偏微分方程变为4阶常微分方程得：
[0032]
[0033][0034]
式(4
‑
5)可计算得通解：
[0035][0036]
式中a1～a8为常数，并且
[0037][0038][0039]
根据euler
‑
bernoulli梁特性，可得到模态坐标系下的角位移θ
y
、θ
z
，力矩m
y
、m
z
和内力q
y
、q
z
的表达式：
[0040][0041]
将式(6)和式(9)写成矩阵的形式可得：
[0042][0043]
即z(x)＝b(x)
·
a,a＝[a1,a2,a3,a4,a5,a6,a7,a8]
t
，其中z(x)为状态矢量，即方程(10)等式左边项，b(x)代表方程(10)等式右边的第一个矩阵；令x＝0代入方程(10)可得到z(0)＝b(0)
·
a；由于b
‑1(0)
·
z(0)＝a，因此有z(x)＝b(x)
·
b
‑1(0)
·
z(0)＝u(x)
·
z(0)；
[0044]
所以该段euler
‑
bernoulli梁从输入点到x长位置处的传递矩阵为：
[0045]
u(x)＝b(x)
·
b
‑1(0)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0046]
此处，
[0047][0048]
根据式(11)，可写出长度为x的双自由度euler
‑
bernoulli梁传递矩阵形式为：
[0049][0050]
其中，
[0051][0052][0053][0054][0055][0056][0057][0058]
[0059][0060]
作为本发明的进一步限定，步骤三中确定其他多体动力学元件的传递矩阵包括：
[0061]
法兰、螺栓和支撑其的平台被简化为集中质量，作用在塔筒内部；此类常规多体动力学元件，可以直接写出其柔塔系统坐标系中的传递矩阵形式：
[0062][0063]
式中，m
i
为该元件的质量，ω同样为角频率；
[0064]
塔筒顶部支撑的风力机部件被简化为一空间刚体；作为常规多体动力学元件，根据状态矢量和坐标系，同样可直接写出质量为m的空间刚体传递矩阵形式：
[0065][0066]
式中，
[0067][0068]
其中j为刚体相对于输入点的惯量矩阵，下标代表转动轴。
[0069]
作为本发明的进一步限定，步骤四中，所述建立柔塔系统总传递方程包括：
[0070]
确认系统两端边界条件，塔筒底端为固定边界，位移和角位移都为0，则状态矢量为z
1,0
＝[0,0,0,0,m
y
,m
z
,q
y
,q
z
]
t
；顶部构件为自由振动，力和力矩都为0，则状态矢量为z
n+1,n
＝[y,z,θ
y
,θ
z
,0,0,0,0]
t
；根据mstmm对链式多刚柔体系统的定义，按传递方向拼接各元件传递矩阵，写出系统总传递方程：
[0071]
z
n+1,n
＝u
n
u
n
‑1...u
i
...u2u1z
1,0
＝u
all
z
1,0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0072]
式中，u
i
下标为元件编号为i的元件传递矩阵，u
all
＝u
n
u
n
‑1...u
i
...u2u1为总传递矩
阵。
[0073]
一种风力机柔塔系统固有振动特性快速仿真方法，对所述建模方法得到的柔塔系统总传递方程求解，得到柔塔系统的振动特性。
[0074]
作为本发明的进一步限定，所述求解总传递方程，得到柔塔系统的振动特性包括：假设u
all
为：
[0075][0076]
将边界条件带入式(18)，可得到：
[0077][0078]
根据式(19)中的线性齐次项，柔塔系统的特征方程可写成：
[0079][0080]
对方程(20)求解，即可求出柔塔系统固有频率；求解对应于固有频率的系统边界点状态矢量z
1,0
、z
n+1,n
，进而通过元件传递方程得到系统全部联接点的状态矢量，即为系统的振型。
[0081]
与现有技术相比，本发明的有益效果在于：
[0082]
(1)本发明提出一种风力机柔塔系统固有振动特性快速建模及仿真方法，推导过程简单且无需划分大量网格单元，省去了每改变一次参数就要对塔筒进行重新建模的过程；同时也能给工程上类似的细长柔性结构的动力学建模与快速仿真提供参考；
[0083]
(2)本发明采用多体系统传递矩阵法，无需划分大量网格单元，具有矩阵阶次低、计算量小、计算效率高的优点，更能满足工程上快速计算的需求；
[0084]
(3)本发明在风力机柔塔模型的建模中，考虑了顶部质量，考虑柔性塔截面非线性变化的特性，使模型更接近实际；同时，还考虑了柔塔与机舱等其它风力机部件之间的刚柔
耦合效应。
附图说明
[0085]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
[0086]
图1是本发明方法流程示意图。
[0087]
图2是风力机柔塔系统简化物理模型示意图。
[0088]
图3是柔塔的梁模型示意图。
[0089]
图4是空间刚体模型。
[0090]
图5是基于mstmm计算出的系统圆频率计算结果图。
[0091]
图6是计算出的柔塔第一阶振型(a)是z方向；(b)是y方向。
[0092]
图7是计算出的柔塔第二阶振型(a)是z方向；(b)是y方向。
[0093]
图8是计算出的柔塔第三阶振型(a)是z方向；(b)是y方向。
具体实施方式
[0094]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0095]
现有技术中，多体系统传递矩阵法(mstmm)由芮筱亭院士及其团队建立，用于多体动力学分析。该方法能实现一般多刚体、多刚柔体系统动力学研究，其高计算效率方便运用于工程中。运用于本发明中对风力机柔塔系统的多体动力学建模，可以实现振动特性的快速计算，本发明可为工程上类似的细长柔性结构的动力学建模与快速仿真提供参考。
[0096]
本发明的一种风力机柔塔系统固有振动特性快速建模及仿真方法，参见图1所示。具体包括以下步骤：
[0097]
步骤一：对风力机柔塔系统进行合理的简化，将系统离散成各多体动力学元件；
[0098]
现有的三叶片水平轴风力机可归结为由多个刚体(机舱、轮毂等)和柔体(叶片和塔架)组成的多刚柔体系统。风力机塔筒底部靠桩土固定，塔筒顶部连接机舱，而机舱通过传动轴与外部的轮毂和叶片相连。为方便生产和运输，塔筒通常需要分成多段，在安装时通过法兰和螺栓分段连接。
[0099]
考虑到风力机系统一个复杂的多刚柔体系统，忽略影响塔筒振动特性的次要部分，对风力机柔塔结构模型做出如下假设和简化：(1)忽略塔筒模型扭转运动；(2)假设风力机叶片、轮毂和机舱为一个存在质量偏心的空间刚体，忽略其中复杂部件的影响；(3)忽略塔筒自身的结构阻尼；(4)根据塔架的高质量比特性，忽略塔筒内部空气质量；(5)假设塔筒内部法兰为集中质量，忽略螺栓、螺丝和螺孔造成的影响；(6)假设塔筒截面始终保持圆形，不考虑壳体截面变形；(7)假设塔筒表面平滑，不存在焊缝凹陷；(8)假设塔筒为梁结构，假设不存在剪切效应和局部屈曲效应；(9)假设塔筒与机舱为固定连接，忽略机舱相对于塔筒
的扭转；(10)忽略风力机内部爬梯、电缆架、平台等对塔架整体强度影响较小的复杂部件。
[0100]
考虑柔塔系统各部件之间的相互作用，基于mstmm将整个系统离散成各多体动力学元件，如图2所示。根据风力机柔塔系统的结构，将塔筒离散成多段考虑顶部构件压力作用、不同结构参数的euler
‑
bernoulli梁；将风力机塔筒内部的法兰、螺栓和支撑其的平台处理为集中质量；将作用在塔筒顶部的风力机叶片、轮毂和机舱处理为一个存在质量偏心的空间刚体。建立坐标系在塔筒底部，以来流方向为z轴，以垂直来流方向、平行于地面的横向为y轴，以塔筒中心轴线为x轴。
[0101]
根据mstmm，定义传递方向为塔筒底部向顶端刚体方向，如图2中箭头指向，按传递方向对系统中的元件进行编号。根据模型中所有元件的运动和相互作用，状态矢量可统一为z＝[y,z,θ
y
,θ
z
,m
y
,m
z
,q
y
,q
z
]
t
。z为沿坐标轴z位移对应的模态坐标列阵，y为沿坐标轴y位移对应的模态坐标列阵，θ
z
为该点处相对于平衡位置相对于在z轴的角位移对应的模态坐标列阵，θ
y
为该点处相对于平衡位置相对于y轴的角位移对应的模态坐标列阵，m
z
为沿坐标轴z内力矩对应的模态坐标列阵，m
x
为沿坐标轴x内力矩对应的模态坐标列阵，q
z
为沿坐标轴z内力对应的模态坐标列阵，q
y
为沿坐标轴y内力对应的模态坐标列阵。
[0102]
步骤二：基于mstmm推导考虑重力作用的塔筒梁传递矩阵；
[0103]
如图3所示为简化为euler
‑
bernoulli梁元件的某段塔筒，该段梁长为l。图中g为受到的梁上方元件的压力，z
i
为输入端状态矢量，z
o
为输出端状态矢量，u(x)代表该段梁的传递矩阵。根据以上条件，建立考虑重力的双自由度euler
‑
bernoulli梁的振动微分方程：
[0104][0105][0106]
式中，ei为梁的弯曲刚度，为单位长度质量，t为时间，z为塔筒流向位移，y为塔筒横向位移，x表示在x轴上的位置。
[0107]
由模态坐标转换，令
[0108]
y(x,t)＝y(x)e
iωt
,z(x,t)＝z(x)e
iωt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0109]
式中，ω为角速度，y(x)、z(x)为y(x,t)、z(x,t)通过模态坐标转换得到。
[0110]
将式(3)其代入式(1
‑
2)，使偏微分方程变为4阶常微分方程得：
[0111][0112][0113]
式(4
‑
5)可计算得通解：
[0114][0115]
式中a1～a8为常数，并且
[0116][0117][0118]
根据euler
‑
bernoulli梁特性，可得到模态坐标系下的角位移θ
y
、θ
z
，力矩m
y
、m
z
和内力q
y
、q
z
的表达式：
[0119][0120]
将式(6)和式(9)写成矩阵的形式可得：
[0121][0122]
即z(x)＝b(x)
·
a,a＝[a1,a2,a3,a4,a5,a6,a7,a8]
t
，其中z(x)为状态矢量，即方程(10)等式左边项，b(x)代表方程(10)等式右边的第一个矩阵。令x＝0代入方程(10)可得到z(0)＝b(0)
·
a。由于b
‑1(0)
·
z(0)＝a，因此有z(x)＝b(x)
·
b
‑1(0)
·
z(0)＝u(x)
·
z(0)。
[0123]
所以该段euler
‑
bernoulli梁从输入点到x长位置处的传递矩阵为：
[0124]
u(x)＝b(x)
·
b
‑1(0)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0125]
此处，
[0126][0127]
根据式(11)，可写出长度为x的双自由度euler
‑
bernoulli梁传递矩阵形式为：
[0128][0129]
其中，
[0130][0131][0132][0133][0134][0135][0136][0137][0138][0139]
本发明中，由塔筒离散成的梁元件传递矩阵形式都为式(12)中的u
beam
，但各个梁元件结构参数不同，需根据各梁元件的参数填入矩阵中。
[0140]
步骤三：确定其他多体动力学元件的传递矩阵；
[0141]
法兰、螺栓和支撑其的平台被简化为集中质量，作用在塔筒内部。此类常规多体动力学元件，可以直接写出其柔塔系统坐标系中的传递矩阵形式：
[0142][0143]
式中，m
i
为该元件的质量，ω同样为角频率。
[0144]
塔筒顶部支撑的风力机部件(机舱、轮毂、叶片等)被简化为一空间刚体，模型如图4所示。i(0,0,0)代表输入点坐标，o(b1,b2,b3)代表输出点坐标，c(c1,c2,c3)代表质心坐标。作为常规多体动力学元件，根据状态矢量和坐标系，同样可直接写出质量为m的空间刚体传递矩阵形式：
[0145][0146]
式中，
[0147][0148]
其中j为刚体相对于输入点的惯量矩阵，下标代表转动轴。
[0149]
步骤四：确定边界条件，建立柔塔系统总传递方程；
[0150]
确认系统两端边界条件，塔筒底端为固定边界，位移和角位移都为0，则状态矢量为z
1,0
＝[0,0,0,0,m
y
,m
z
,q
y
,q
z
]
t
；顶部构件为自由振动，力和力矩都为0，则状态矢量为z
n+1,n
＝[y,z,θ
y
,θ
z
,0,0,0,0]
t
。根据mstmm对链式多刚柔体系统的定义，按传递方向拼接各元件传递矩阵，写出系统总传递方程：
[0151]
z
n+1,n
＝u
n
u
n
‑1...u
i
...u2u1z
1,0
＝u
all
z
1,0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0152]
式中，u
i
下标为元件编号为i的元件传递矩阵，u
all
＝u
n
u
n
‑1...u
i
...u2u1为总传递矩阵。
[0153]
步骤五：求解总传递方程，得到柔塔系统的振动特性；
[0154]
假设u
all
为：
[0155][0156]
将边界条件带入式(18)，可得到：
[0157][0158]
根据式(19)中的线性齐次项，柔塔系统的特征方程可写成：
[0159][0160]
方程(20)只与系统的结构参数和固有频率有关，对于实际的线性振动系统，式(20)必有非零解，即可求出柔塔系统固有频率。求解对应于固有频率的系统边界点状态矢量z
1,0
、z
n+1,n
，进而通过元件传递方程得到系统全部联接点的状态矢量，即为系统的振型。
[0161]
通过该方法可以快速计算系统每个部件的参数对风力机柔塔振动特性，即振型和频率的影响。
[0162]
实施例
[0163]
针对型号为nordex s70/1500的三叶片水平轴风力机柔塔进行多体动力学建模并振动特性分析。得到的柔塔系统的固有圆频率计算结果如图5所示。图中横坐标为圆频率，纵坐标表示δ值的大小，当δ值接近于0时即可以求出圆频率。图中显示的竖线即搜根过程，每条竖线下对应的就是该阶模态的圆频率。图6为计算出的风力机柔塔的第一阶振型，(a)是z方向；(b)是y方向；图7为风力机柔塔的第二阶振型，(a)是z方向；(b)是y方向；图8为风力机柔塔的第三阶振型，(a)是z方向；(b)是y方向。
[0164]
同样，本方法对其他风力机柔塔建立多体动力学模型并实现振动特性分析，只需修改简化多体动力学模型和其中具体参数即可。
[0165]
上实施例的说明只是用于帮助理解本发明的方法及其核心思想。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以对本发明进行若
干改进和修饰，这些改进和修饰也落入本发明权利要求的保护范围内。