多重分形的海上漂浮式风力机系泊状态特性分析

2023-01-11李蜀军肖俊青金江涛王江波

哈尔滨工业大学学报 2023年1期

闻麒，李蜀军，孙康，肖俊青，金江涛，李春,2，王江波，陈泳

(1．上海理工大学能源与动力工程学院，上海 200093；2．上海市动力工程多相流动与传热重点实验室(上海理工大学)，上海 200093；3．江苏东华测试技术股份有限公司，江苏靖江 214500)

中国浅海风资源开发趋于饱和[1]，深海风资源较为丰富，因此海上风电正逐步向深海发展[2]。漂浮式风力机作为深海风能主要开发设备，在长期受海水腐蚀及复杂海洋环境载荷作用下，系泊极易失效进而影响整个漂浮式风力机的安全稳定运行。因此，采用合理有效的故障诊断方法，在系泊蠕变早期阶段实现预警，对保证漂浮式风力机安全运行及延长其使用寿命具有重要意义。

文献[3]提出采用变分模态分解(variational mode decomposition, VMD)对故障信号进行预处理，可有效改善上述问题，并具有良好的鲁棒性和运算效率；文献[4]采用VMD方法开展了轴承的故障诊断，研究发现VMD可提取多个特征；文献[5]将快速傅里叶变换、包络分析与VMD结合，实现了更有效的智能故障诊断；文献[6]研究了VMD方法在故障特征提取中的应用，结果表明VMD与EMD和EEMD相比更有效。但现有振动分析法在强噪声环境下难以有效识别故障状态[7]。因此，须引入新方法对非线性特征进行表征与反映。

分形学为表征振动信号的非线性特征提供了一种分析方法。因故障信号具有明显自相似性[8]，故可利用分形对其进行分析。分形维数是描述信号非线性程度的特征量，包含Hausdorff位数、盒维数和关联维数等。其中，盒维数抗噪性好，可有效表现信号的非线性特征，在故障诊断领域已取得广泛应用[9]。文献[10]将分形盒维数应用于轴承信号的特征提取和故障诊断，证明了其在故障状态识别上的有效性；文献[11]结合盒维数与EMD方法，在传统单一分形维数判断的基础上，提高了特征提取的效率和准确性；文献[12]采用多重分形去趋势波动分析(multifractal detrended fluctuation analysis, MF-DFA)方法对机械故障非线性特征进行分析，以多重分形的多测度刻画分形特征优势弥补了单分形维数存在的缺陷。

目前国内外学者对风力机故障诊断的研究主要集中于发电机与齿轮箱方面，对海上风力机系泊故障诊断研究鲜见。本文以搭载NREL 5 MW风力机ITI Energy Barge平台的系泊系统为研究对象，分析不同位置系泊蠕变及其失效后对漂浮式风力机动态响应的影响；通过VMD对不同位置系泊蠕变下的平台响应数据进行解耦，结合分形盒维数对分量进行优选，提取系泊故障信号中更具代表性的非线性特征；并采用MF-DFA方法解析系泊故障信号多重分形特征，为漂浮式风力机系泊健康状况诊断提供了理论基础和技术支持。

1 研究对象

1.1 漂浮式风力机模型

以ITI Energy Barge平台的NERL 5 MW漂浮式风力机为研究对象，其整机系统见图1，主要参数见表1，漂浮式平台主要参数见表2[13]。

图1 Barge平台漂浮式风力机

表1 NREL 5 MW风力机参数

表2 ITI Energy Barge平台参数

1.2 系泊系统

海上漂浮式风力机平台需安装系泊系统，减小平台位移并确保平台相对于海面的位置，实现在深海较恶劣环境下运行[14]。常见系泊方式有悬链线式系泊和张紧式系泊2种。

系泊系统中各根系泊能通过辐射状布置提供各方向的回复力。目前系泊缆通常为合成纤维缆，其为工程提供了便利，但相较传统锚链及钢丝绳具有更加复杂的非线性特征，为相关研究带来更大的难度。系泊发生蠕变时特征为刚度与系泊长度的小幅变化，在数值模拟上常通过系泊刚度与长度按百分比变化实现。

Barge平台利用8根从四角上导缆孔引出的系泊与海底锚点相连，采用悬链线式系泊，平台俯视图见图2，系泊参数见表3[15]。

图2 系泊示意

表3 系泊参数

1.3 锚泊力分析

集中质量法常应用于求解动态系泊模型[16]，见图3。自底部锚固点开始，将系泊缆切分为以N+1个节点连接的N个大小相等的系泊段，每个系泊线段i+1/2具有相同的未拉伸长度L、体积等效直径d、弹性模量E和惯性阻尼系数Cint等特性。每个节点i位置由位置向量ri定义，包含节点位置的三维坐标。

图3 动态系泊模型示意

在集中质量法中，通过将每个节点分配为相邻2个线段的组合质量的一半，将系泊缆的质量离散为每个节点的点质量。线段i+1/2的内部刚度和阻尼力分别表示为Ti+1/2和Ci+1/2，在节点i集中的线段的质量记为Wi，Bi为节点i的净浮力。

节点i的3×3质量矩阵表示为

(1)

式中：ρ为系泊缆的密度,I为单位矩阵。

节点i的附加质量矩阵表示为

(2)

节点处切线方向可表示为

(3)

式中ri为节点i的位置向量。

进而，各节点i的系泊缆运动方程为

Ci+1/2-Ci-1/2+Wi+Bi+Dpi+Dqi)

(4)

式中Dpi和Dqi分别为横向和切向上的阻力。

2 环境载荷与参数

复杂的海洋环境下，漂浮式风力机会受到风浪流等多种载荷作用，其中上部风轮受到的风载荷及下部漂浮式平台受到的海流与波浪载荷影响最大[17]。

2.1 风载荷

基于叶素动量理论，通过动态入流理论对风轮平面的诱导速度进行求解，进而解得沿叶片展向各位置攻角[18]。

其中，诱导速度可表示为

(5)

进而，叶片所受推力与扭矩可表示为

(6)

式中：ρ为空气密度，c为距轮毂中心r处风力机叶片弦长，Ω为风轮转速，b为切向诱导因子，a为轴向诱导有因子，θ为来流攻角，U∞为来流速度，Cl为升力系数，Cd为阻力系数。

2.2 海流载荷

一般情况下海流速度较缓，对海工结构的主要作用力为拖曳力，可表示为

(7)

式中：A为平台沿海流方向等效面积，Uo为流速，Cd为阻力系数，ρs为海水密度。

2.3 波浪载荷

由于Barge平台体型较大，其对波浪场的影响不可忽略，但须考虑平台在海水中的绕射效应，故本文采用辐射/绕射理论求解Barge平台在海水中的波浪载荷，速度势函数可表示为[19]

(8)

式中：φi为入射势，φd为绕射势，φr为辐射势。

其中，速度势函数需满足拉普拉斯方程与边界条件：

(9)

式中：g为重力加速度，d为海深，η为入射波波面函数，n为浮体湿表面外法向量。

漂浮式风力机平台受到的波浪力与波浪力矩可表示为：

FW=∬SB-pndSB

(10)

MW=∬SB-p(r×n)dSB

(11)

式中：p为压强，FW为波浪力，SB为浮体湿表面积；MW为波浪力矩，r为结构物表面到基点矢量。

2.4 环境参数

基于TurbSim[20]，参照IEC 61400-3，选用Kaimal风谱，生成11.4 m/s湍流风场，见图4(a)。海流为均匀流，流速0.8 m/s。波浪基于P-M谱，谱峰周期10.1 s，有义波高6 m，水深为150 m，不规则波高时域曲线见图4(b)。

(a)湍流风场

(b)波高

3 非线性特征提取

3.1 变分模态分解

变分模态分解算法基于外差解调、Hilbert变换和维纳滤波，将原始信号分解为K个模态分量，可有效抑制EMD及EEMD等方法存在的模态混淆现象[10]。

建立约束变分模型表达式如下：

(12)

式中：K为模态分量数，j为虚数单位，σ(t)为单位脉冲函数，∂t为偏导运算，*为卷积运算，{uk}={u1,u2,u3,…,uK}、{ωk}={ω1,ω2,ω3,…,ωK}分别为经变分模态分解得到的K个模态分量与其中心频率。

同时，利用惩罚系数α和拉格朗日乘子λ(t)求解变分约束问题。其增广拉格朗日表达式为

L({uk},{ωk},λ)=

(13)

其迭代表达式为：

(14)

(15)

(16)

式中：τ为保真系数，^为傅里叶变换，n为迭代次数。

求解过程中各模态分量的带宽和中心频率亦不停迭代，直至达到更新停止条件，其判别表达式为

(17)

式中β为判别精度。

3.2 盒维数

考虑分形具有自相似性、无标度性及自仿射性，而分形维数是一种刻画分形不规则性的有效度量方式，其中盒维数因其原理简单及抗噪性好的特点，广泛用于信号处理中，可度量信号的不规则度。其定义为

(18)

设离散信号y(i)⊂Y，Y为n维欧氏空间Rn上的闭集。利用网格划分Rn，网格ε尽量小，Nε为划分网格数。仅定义难以求得，故采用近似法求解。扩大网格ε至kε，k∈+。

min{yk(i-1)+1,yk(i-1)+2,…,yk(i-1)+k+1}|

(19)

Nkε=P(kε)/(kε)+1

(20)

式中：i=1,2,…,N/k；N为采样点数；k=1,2,…,M,M

自lg(kε)-lgNkε选取线性程度良好的范围，令其两端点分别为k1与k2，则：

lgNkε=alg(kε)+bk1≤k≤k2

(21)

通过数据拟合范围内曲线斜率：

(22)

可计算得盒维数DB为

(23)

4 多重分形

系泊蠕变下平台响应信号较为复杂，单分形仅从单一测度无法准确描述信号的分形特性[21]。多重分形是在分形结构上具有多个标度指数的奇异测度组成的集合，可对分形测度或物理量在支集上的分布情况进行精细刻画[22]，弥补单测度分析时易导致判别位置混淆的缺陷，且保留了单分形维数抗噪性优势。因此，本文采用多重分形去趋势波动分析法，研究不同位置系泊蠕变下平台响应信号的多重分形特征。

4.1 广义Hurst指数

现有振动信号分析法在强噪声环境下难以有效识别故障状态，为研究非平稳信号的动力学特性，采用MF-DFA进行分析，具体步骤为[23]：

1)计算振动信号{x(k),k=1,2,…,N}的均值累积离差Y(i)：

(24)

2)等分划分区间Y(i)，m为子区间个数，s为子区间长度。因m可能为非整数，反向重复划分，共得2m个子区间。

3)子区间均方差F2(s,v)常以最小二乘多项式拟合：

(25)

式中v=1,2,…,m。

区间Y(i)无法单向等分划分时：

(26)

式中：v=m+1,m+2,…,2m；yv(i)为子区间v的拟合多项式。

4)通过二元多项式拟合计算得到测度波动函数Fq(s)为

(27)

当q等于2时，式(27)则为去趋势波动分析。

5)改变子区间长度s，重复上述步骤。若x(k)存在长程相关性，可知q和hq、Fq(s)和s关系如下：

Fq(s)∝shq

(28)

若广义Hurst指数hq不随q变化而发生变化，可知x(k)具有单分形特征；当广义Hurst指数hq随q变化而变化，则x(k)具有多重分形特性。

4.2 多重分形谱

通过配分函数得关系式为

τq=qhq-1

(29)

通过勒让德变换得奇异指数α和多重分形谱f(α)为

αq=dτq/dq

(30)

f(α)=qαq-τq

(31)

利用式(29)～(31)计算，得奇异指数α、谱函数f(α)及广义Hurst指数hq间的关系为：

αq=hq+qdhq/dq

(32)

f(α)=q(αq-hq)+1

(33)

多重分形谱f(α)反映了x(k)的分形情况。若f(α)图像表现为单峰状曲线，则为多重分形；若f(α)图像为水平的直线，则仅为单重分形。

5 结果与分析

5.1 动态响应分析

设置风浪流均以垂直风力机平面方向入射，因Barge平台8根系泊对称分布，故仅对其中4根进行分析。漂浮式风力机平台于系泊正常、发生蠕变及失效3种情况下六自由度上响应曲线与1～4号系泊在各阶段最大响应幅值对比见图5。

图5 系泊正常、蠕变及失效下平台六自由度响应

如图5所示，0～3 000 s，3 000～6 000 s，6 000～9 000 s分别为系泊正常、蠕变及失效情况，图中上部各曲线分别表示各系泊在3个阶段的漂浮式风力机平台动态响应，系泊1～4在3个阶段的最大响应幅值对应图中下部柱状图。通过对比图5中平台于系泊不同状况下在六个自由度上的响应曲线及响应幅值可知，系泊发生蠕变与失效后平台响应产生不同程度地增大。其中，纵荡、横荡、横摇及艏摇受到的影响最大，对纵摇和横荡影响较小；系泊蠕变后，平台响应幅值增加极小，系泊失效后，响应幅值急剧增大；不同位置系泊状态变化对平台动态响应影响不一，靠近迎风浪侧的1～2号系泊对平台影响较大，远离迎风浪侧的3～4号系泊影响较小。

5.2 信号处理

由上述分析可知，平台艏摇响应在各自由度中最为敏感，故下文对各系泊蠕变状态下的平台艏摇响应数据进行分析。正常系泊平台艏摇时域信号与上节中系泊1～4发生蠕变后平台艏摇时域响应见图6。由图6可知，各系泊蠕变状态下平台艏摇响应微弱，时域图差异甚微，难以判断系泊蠕变位置。

采用VMD方法对5种系泊状态下平台艏摇响应信号进行分解，每种响应信号均具有3个模态分量。因盒维数可表征信号间的自相似性，故以其作为筛选特征分量的依据。结果见图7。

因盒维数越小表明信号非线性越弱，越规则[24]。故以其值最小的为特征信号，计算原始信号与采用本文方法所选取信号的盒维数，结果见表4。

图6 5种系泊状态下平台艏摇响应时域

图7 各模态分量分形盒维数

表4 5种系泊状态下分形盒维数对比

由表4可知，各系泊状况下平台艏摇响应原信号分形盒维数较大，自相似性弱；经本文方法处理后所提取信号盒维数小，自相似性强。表明经VMD方法可较大程度滤除噪声，信号中非线性特征可得到有效提取。

5.3 多重分形特性分析

非线性特征经VMD捕捉后采用MF-DFA方法进行多重分形分析，求取各多重分形参数与波动函数阶数q的变化关系曲线。5种系泊状态下平台艏摇响应信号的τq与波动函数阶数q的变化曲线见图8。

图8 平台响应信号标度指数变化曲线

由图8可知，5种系泊状态下，τq在取值范围内，其值均随q值增大而呈非线性增加。q<0时，系泊4状态下τq最大，系泊1最小，系泊3、正常系泊和系泊2的τq依次减小；q>0时，系泊1状态下τq最大，系泊4最小，系泊2、系泊3及正常系泊状态下τq顺序减小；当波动函数阶数=0时，各信号标度指数均为0.6，进一步说明5种系泊状态下平台响应信号具有多重分形特征。

图9为5种系泊状态下平台响应信号Hurst数与q的关系曲线。

图9 平台响应信号Hurst指数变化曲线

由图9可知，5种系泊状态下hq与q的变化曲线呈非线性变化。在全取值范围内，随q增加，hq持续减小，由此进一步表明各信号均为多重分形序列。当q>0时，系泊1状态下hq随q降速最大呈弱非线性，系泊3状态下降速最小呈强非线性；当q>0时，系泊1状态下hq降速最小呈强非线性，系泊4状态下降速最大呈弱非线性。由此表明，不同系泊状态下的广义Hurst指数存在较大差异，通过hq与q的斜率关系可判断各系泊状态下平台响应信号具备不同的非线性特征。

图10为5种系泊状况下平台响应信号的奇异谱函数f(α)与奇异指数α的关系曲线。

图10 平台响应信号多重分形谱变化曲线

由图10可知，5种系泊状态下平台响应信号多重分形谱均呈几字形，曲线峰值处大致相同，各系泊状态下分形谱变化存在差异，说明各系泊状态下平台响应信号具备多重分形特征。

α0为变化关系曲线峰值的横坐标，其可说明信号在q=0下的不均匀程度[25]。文献[26]研究表明信号非线性状况可通过分析α0得到有效判断。因此本文对5种系泊状况下平台响应信号各随机选取30组数据，计算并分析其α0值，解析各系泊状态下平台响应信号的非线性特征。

原信号与经VMD处理后5种系泊状态下各30组平台响应信号α0值变化趋势见图11。

由图11可知，各系泊状态下平台响应信号的α0值各异。系泊1状态下α0最大，说明该信号包含复杂的数据模式。而系泊2状态下α0最小，可见此数据模式单调，整体系泊系统受非线性作用较小。但α0在噪声环境下无法准确描述系统非线性状况。因此，当系泊处于不同状况时，α0值在图11(a)中出现了交叉现象，系泊1～4的α0值无法准确描述信号不规则程度。但采用VMD方法后，提取信号的非线性特征更具代表性，能够精准刻画信号规则状况，消除曲线交织。由此可依靠数据复杂度判断系泊状态。