2024中国可再生能源大会：风机关键部件动力学分析的精细积分法研究

风机关键部件动力学分析的精细积分法研究风机关键部件动力学分析的精细积分法研究2研究背景与现状风机基础结构的基本计算模型线性与非线性动力学问题的精细积分法350050新增累计400装机容量(GW)装机容量(GW)300•我国风电装机容量新增、累计呈快速发展趋势装机容量(GW)装机容量(GW)300200•我国海上风电发展迅猛，近几年装机规模庞大20010000201420152016201720182019202020212022202320142015201620172018201920202021202220232014201520162017201820192020202120222023•风电关键结构尺寸规模逐年增大，建造成本变高•风电机组荷载会对海上风机的基础结构产生影响•海上风电机组关键结构常用梁模型进行计算量缩减4风机基础结构的基本计算以及算法改进基于SubDyn展开单桩三角架单桩三角架SubDyn兼容两类梁模型：Euler-Bernoulli梁单元以及Timoshenko梁单元，在通过单元刚度形成全局刚度的时候，局部与全局坐标系的示意图如上图；最终形成的动力学方程可以形成如下的线性形式5!•Newmark法!•Newmark法•Wilson法•广义u法•精细积分方法风机柔性叶片基本计算模型风力发电机计算仿真•国内外开发了Bladed、OpenFAST、ENFAST、HAWC2等的各类风力发电机的仿真软件•OpenFAST中有专门进行叶片动力学仿真的动力学模块BeamDyn•叶片在空间中的形态具有较大的几何非线性。对于柔性叶片而言，其模型具有大位移、大转角甚至大应变的特点•Euler-Bernoulli梁•Timoshenko梁•共旋坐标法•绝对坐标法•几何精确梁模型6风机柔性叶片基本计算模型风机柔性叶片基本计算模型•将风机叶片视为梁模型•荷载以外力荷载形式施加在叶片上•叶片通过锚点固结在中心转子上ω锚点锚点!!转子•在锚点处建立局部参考系，与锚点运动状态一致，局部参考系中定义了梁的参考构型•梁在局部坐标系内柔性变形，随局部坐标系在整体坐标系中旋转转子随叶片一起运动的坐标系ω锚点随叶片一起运动的坐标系整体坐标系随叶片一起运动的坐标系7•每个单元定义了一个共旋构型•每个单元定义了一个共旋构型•单元变形与受力分析是基于这个共旋构型来计算的•共旋构型与参考构型之间的变换是刚体旋转变换•适用于大旋转大变形问题，但是单元还是小应变的2EI0--2EI0--M1M2L02EI--M1M2L02EI--内力通过旋转变换至整体坐标系中，有了整体坐标下的KM=TTKM0T•位移变换一部分是全局位移转换到单元坐标系后的值，另一部分是由于单元坐标系发生旋转所产生的位移•整体坐标中的刚度阵满足与内力的切线关系∂pK=∂u8设置时间步长度η,可以得到按时间步更新的设置时间步长度η,可以得到按时间步更新的迭代格式，在第k时间步有：Vtk+η=eAηVt+∫eAη一ξg(tk+ξ)dξ若在时间步内，外荷载可以表现为线性的则非齐次项的积分可以积分得出，在第k时间步的迭代格式如下：r1线性与非线性动力学问题的精细积分法精细积分方法（PIM）上式是一个离散后的动力学方程，以如下方式引入状态矢量V，表示为矩阵形式通过这种方法，可以将原方程化成如上的等价的一阶微分方程形式V(t)=eAtV0+∫eA(t一ξ)g(ξ)dξ如此形式的微分方程，其解可以表示为：齐次方程的通解+非齐次项的特解。∫eA(t一ξ)f(ξ)dξ是Duhamel积分的形式9在时间步开始阶段，在时间步开始阶段，−fi+kiu的值为。1.首先假定时间步内−fi+kiu项不发生改变，即认为时间步内的−fi+kiu为。2.计算r0和r1，通过迭代格式从vi计算vi+13.而后通过vi+1修正ki，再计算r0和r1，从vi重新计算vi+14.如此往复直到计算获得的vi+1和上一循环的vi+1的差别在精度范围内。如此完成非线性问题的预估-矫正格式线性与非线性动力学问题的精细积分法非线性动力学问题的精细积分法若在精细积分的计算中，于时间步中考虑全局刚度与质量的变化，考虑其非线性效应：内力项fi在时间步内并不一定时刻满足fi=kiu，全局刚度与内力仅满足切线关系ki=∂fi/∂ui。则可以参考这样的思路采用类似的预估-矫正方法构造求解方法：回到最初的动力学方程，可以构造出时间步内的线性系统Mi+ciuǗ+kiu=Qt−fiu+kiu即在前文推导的基础上，外荷载取得g(t)=tt−fiu+kiu)。M−1(Q()()通过多次积分的预估-矫正方法来对f进行矫正，来构造非线性的迭代格式。基于SubDyn软件，修改并发展SubDyn中的时程积分方法。计算左图所示的基础结构：•单桩底部与海底固结•顶部接口处输入为零•考虑单桩自重的影响分别采用龙格库塔与精细积分方法在不同的时间步长下进行计算，对于不同时间步长下各类方法的结果比较，以4阶龙格库塔（RK4）方法在0.001s时的计算结果作为参考解使用MATLAB进行数据的可视化•在时间步长为0.001s时，PIM的计算•在时间步长为0.0015s时，这是RK4所能计算的最大时间步长，P•在时间步长为0.003s时，PIM的计算三脚架模型以及离散计算简图基于SubDyn软件，修改并发展SubDyn中的时程积分方法。计算左图所示的基础结构：•底部固结•顶部输入x=0.1m•考虑单桩自重影响分别采用龙格库塔与精细积分方法在不同的时间步长下进行计算，并对比结果；对于不同时间步长下各类方法的结果比较，以4阶龙格库塔（RK4）方法在0.005s时的计算结果作为参考解•在时间步长为0.005s时，PIM的计•在时间步长为0.1s时，PIM的计三脚架模型以及离散计算简图基于SubDyn软件，修改并发展SubDyn中的时程积分方法。计算左图所示的基础结构：•底部桩腿固结•顶部输入外力为5sin(10t)•并考虑单桩自重分别采用龙格库塔与精细积分方法在不同的时间步长下进行计算，并对比结果；对于不同时间步长下各类方法的结果比较，以4阶龙格库塔（RK4）方法在0.005s时的计算结果作为参考解•在时间步长为0.005s时，PIM的计算结果与参考解一致，验证了精细积分法的准确性•在时间步长为0.05s时，为RK4能计算的最大时间步长，RK4的计算结果略微偏差，PIM计算结果准确•在时间步长为0.1s时，PIM的计算结果仍与参考解一致，说明了精细积分法的稳定性•基于SubDyn的精细积分法能计算更大时间步长的数值积分，在求解大规模多自由度问题有显著优势若先假设u(t)的解析形式，则可得出外力Qt的对应形式假设u(t)具有如下解析形式u(t)=Lsinwt代入动力学微分方程，可得Q(t)的形式如下=—MLw2sin(wt)+K1Lsinwt+K2L3sin3wtQt=MLw2sin(wt)+K1Lsinwt+K2L3sin3wt在后续具体计算模拟过程中取得参数设置Y=设置Y=0.5,β=0.25，产生了较大的数值阻尼其计算的角速度较解析解明显变慢，同时振幅也很快衰竭使用广义u法对此问题进行计算使用与Newmark法中相同的Y,β参数，并取广义α参数P=0，其计算结果与Newmark法类似使用广义u法对此问题进行计算使用与Newmark使用广义u法对此问题进行计算使用与Newmark法中相同的y,β参数，并取广义α参数P=0.5，其计算结果有所改善效果较Newmark法与广义α法好没有出现角速度与振幅的偏移共旋坐标梁的静力学算例设置共旋坐标梁的静力学算例设置•设置梁长度为10m，采用φ20Q235圆钢的材料参数；•使梁的最左侧节点固结于空间，约束其全部的自由度；•施加弯矩于最右侧节点，弯矩大小EIz，为其弯成圆周的理论弯矩大小；•设置30荷载步，通过NR迭代方法求解共旋坐标梁的静力在静力学与动力学求解过程中，由于模型的几何非线性，在形成全局坐标下的刚度阵时，需要记录单元的旋转历史。•形成10节点的有限元梁模型：采用共旋欧拉梁的单元刚度，取6×6的阻尼阵C为对角阵，于对角线取得阻尼值•在最右侧节点施加大小为弯矩大小EIz，为梁右端节点旋转角度为90°的理论弯矩；•时间步长为0.01s，总时间步长为2000步；外力矩在500时间步内线性加载至约定的外力矩M。πMM计算结构简图M0弯矩施加a)梁变形后的构型图b)精细积分法的计算过程a)梁变形后的构型图b)精细积分法的计算过程•通过MATLAB进行编程，实现共旋坐标梁模型，并采用精细积