CN117574723A - 一种面向风电机组塔筒在线运行状态监测的数字孪生方法 - Google Patents

Abstract

本发明公开了一种面向风电机组塔筒在线运行状态监测的数字孪生方法，包括以下步骤：1)构建风电机组塔筒的简化模型，并使用有限元法对简化模型进行离散得到有限元模型；2)使用本征正交分解对塔筒的有限元模型进行降阶，分析不同阶数下降阶模型的精度，选取当精度符合实际工程需求前提下，降阶阶数最小的降阶模型作为最终的降阶模型；3)在计算机中编写上位机软件，将降阶模型部署到上位机软件中，同时搭建塔筒的物理实体，使用降阶模型对物理实体的应力与应变进行在线监测；使用本方法构建的降阶模型能有效的对物理实体的物理量进行实时的计算，且达到较高的精度，进而对塔筒的运行状态进行在线监测。

Description

一种面向风电机组塔筒在线运行状态监测的数字孪生方法

技术领域

本发明涉及风电机组的状态评估技术领域，具体涉及一种面向风电机组塔筒在线运行状态监测的数字孪生方法。

背景技术

塔筒是风电机组的关键支撑部件，其健康状态影响着整台风电机组的运行安全。因此，需要对风电机组塔筒进行在线状态监测。

数字孪生是一种以数字化的方式将装备物理实体全空间尺度、全生命周期地映射到数字空间中构建数字孪生模型的方法；通过数字孪生模型与物理实体的虚实交互反馈，感知装备运行情况并对其进行优化控制，进而使得装备发挥出最大效能，数字孪生在装备的在线运行状态监测方面有着广阔的应用前景。然而，复杂物理实体的数字孪生涉及求解装备实体的物理场微分方程，并且对微分方程的求解速度有着较高要求，而传统数值方法(如有限元法、有限体积法、有限差分法)在解算物理场微分方程时，因离散后的方程中自由度数量过高导致计算成本过大，较难直接用于构建风电机组塔筒的数字孪生模型。

模型降阶方法对高维的有限元模型通过数值拟合、物理简化、数学投影等方法获得低维的降阶模型，是一种有效降低计算成本的方法。使用模型降阶方法构建塔筒的数字孪生模型可实现对塔筒的物理量(如应力，应变等)进行实时的计算分析，进而实现塔筒运行状态的在线感知，有益于保障风电机组的安全高效运行。

发明内容

本发明的目的在于：提供一种面向风电机组塔筒在线运行状态监测的数字孪生方法，以及整体的实施过程，技术实施方案如图1所示。

本发明的技术方案如下：

一种面向风电机组塔筒在线运行状态监测的数字孪生方法，包括如下步骤：

步骤1：根据风电机组塔筒的几何特征获得其物理简化模型，分析塔筒实际受力情况，构建塔筒有限元模型；

步骤2：使用本征正交分解(Proper Orthogonal Decomposition，POD)对步骤1所构建的有限元模型进行降阶分析，对比不同阶数下降阶模型的精度，并确定最终的降阶阶数，最终降阶阶数应满足在该阶数下的降阶模型符合工程需求，本发明要求最终降阶模型的预测值与有限元模型的计算值之间的相对误差小于1％；

步骤3：在计算机中设计并构建上位机软件，并将塔筒的降阶模型部署到上位机软件中；搭建塔筒的物理实体并安装应变传感器与倾角传感器，使用单片机读取传感器数据并将数据上传到上位机中转换为位移边界条件，再将位移边界条件输入到降阶模型中进行计算；通过计算应变传感器的应变测量值与上位机中降阶模型的应变计算值之间的误差来判断用于塔筒运行状态监测的数字孪生方法的准确度。

进一步的，步骤1所述流程为：简化风电机组塔筒间的法兰连接结构，认为法兰连接紧固，将塔筒简化为梁结构，进而构建出塔筒的物理简化模型。塔筒简化模型的弯曲变形示意图如图2所示。图中，坐标系XYZ的原点O为塔筒弯曲变形前的轴线与塔筒下端面的交点；n₀，n₁分别为塔筒弯曲变形前与弯曲变形后塔筒顶端端面的法向量，两向量的夹角为塔筒顶端处端面的倾角θ。塔筒顶端处的挠度w与端面倾角θ的表达式为：

式中，F为塔筒顶端施加的集中载荷，单位为N；E为塔筒材料的弹性模量，单位为Pa；I为塔筒截面惯性矩，单位为m⁴；L为塔筒长度，本发明实例，L设为1m。

在三维空间描述塔筒弯曲变形不仅需要挠度w还需要确定弯曲的方位角方位角的定义如图2所示，是用于表达塔筒的弯曲方位。因此塔筒有限元模型的位移边界条件由w与/>这两个参数确定，w与/>的定义如图2所示，根据线性叠加原理与力的分解确定塔筒位移边界条件，得到：

式中，q_bx＝[1,0,1,0,...,1,0]^T，q_by＝[0,1,0,1,...,0,1]^T，向量长度由具体位移边界条件决定，在本实例中，向量长度为24；q_b为塔筒的位移边界条件决定的节点位移向量，即塔筒有位移约束的节点对应的节点位移向量；w_x为总挠度在X方向的分量；w_y为总挠度在Y方向的分量。

将塔筒顶端处的挠度w转换为位移边界条件，并构建出塔筒的有限元模型，位移边界条件下塔筒有限元模型的平衡方程如式(3)所示：

式中，K_aa为塔筒无位移约束节点对应的刚度矩阵；K_ab、K_ba为塔筒位移约束节点与无位移约束节点之间的耦合矩阵；K_bb为塔筒有位移约束节点对应的刚度矩阵；q_a为塔筒无位移约束节点的节点位移向量；g_a、g_b分别为塔筒无位移约束、有位移约束节点对应的节点重力载荷向量；R_r为塔筒有位移约束节点对应的节点支反力向量。

进一步的，步骤2所述过程为：将本征正交分解与有限元方法结合构建塔筒位移场降阶模型，包括：

求解步骤1的有限元模型，得到节点的位移解集合，从集合中选取m个样本组成矩阵X＝[q₁,q₂,...,q_m]；其中q_i＝[u₁,u₂,...,u_z]^T是节点的位移解，z是有限元模型的自由度数量；设{ζ₁,ζ₂,...,ζ_n}为X的一组数量为n的标准正交基，将其写为矩阵形式记为Φ；对于塔筒的有限元模型而言，任意时刻的节点位移值q_t可由{ζ₁,ζ₂,...,ζ_n}表示：

式中，

选择{ζ₁,ζ₂,...,ζ_n}中前k个标准正交基向量来表示qt有：

式中，Φ_k＝[ζ₁,ζ₂,...,ζ_k]，k＜n；b_t(k)＝[b_t1,b_t1,...,b_tk]^T。

令误差函数为：

为获得误差函数在约束条件{ζ₁,ζ₂,...,ζ_n}为标准正交基下的最优解，引入Lagrange因子u_ij(i,j＝k+1,k+2,…,n)，构造Lagrange函数，寻找最优解对应的标准正交基。

式中，δ_ij为克罗内克函数。

将式(7)两端对ζ_j求偏导，得到：

式中，Φ_n-k＝[ζ_k+1,ζ_k+2,...,ζ_n]^T，u_j＝[u_k+1,j,u_k+2,j,...,u_n,j]^T。将上式写为矩阵形式得到下式：

式中，U_n-k＝[u_k+1,u_k+2,...,u_n]^T。为求最优解，令式(9)等于0，并对其左乘由于矩阵Φ_n内部向量标准正交，因此得到/>易知，U_n-k为半正定矩阵，故存在正交矩阵P可将U_n-k对角化，得到/>将上式两端左乘Φ_n-kP，得到下式：

式中，Λ为U_n-k的对角矩阵。将Λ的对角元素降序排列，从式(10)可知Λ的第i个对角元素为XX^T的特征值λ_i，而矩阵Φ_n-kP的第i列是XX^T对应于λ_i的特征向量。考虑到正交矩阵在Frobenius范数下是保范的，因此式(6)可重写为：

取X前n项奇异值对应的左奇异值向量组成变换矩阵Φ_n，此时Φ_n为误差函数在约束条件{ζ₁,ζ₂,...,ζ_n}为标准正交基下的最优解，误差函数有最小值为矩阵XX^T最后n-k个特征值的和。定义下式：

式中，λ_i为XX^T按照降序排列的特征值。若I(k)≥d％则称k维基向量保留有原始样本d％的特征信息。构建不同阶数下的降阶模型，并分别计算不同降阶阶数下降阶模型的精度，确定合适的降阶阶数并得到降阶模型，最终确定的降阶阶数应满足该阶数下的降阶模型符合实际工程需求。

进一步的，对于步骤3，搭建塔筒的物理实体，构建塔筒缩小的简化模型，在塔筒物理实体上安装应变传感器与倾角传感器；编写单片机程序，实现读取传感器测量数据并向计算机上传数据；在计算机上设计并构建上位机软件，将降阶模型部署到上位机软件中；上位机实现读取单片机上传的传感器数据，将数据转换为位移边界条件，并输入到降阶模型中实现快速计算，得到当前塔筒物理实体的应变与应力计算值，以及节点的位移解。降阶模型完成计算后，在上位机中对计算结果进行处理，即基于节点位移解快速地量化塔筒物理实体变形并在上位机中进行显示；通过降阶模型计算得到的应变值与传感器测量的应变值得到降阶模型的应变计算结果与物理实体真实应变之间的误差；结合节点的位移解与应力计算值实现应力云图的可视化，直观并清楚地表达出塔筒物理实体的应力分布情况。塔筒的物理实体与降阶模型之间的虚实交互分析表明，使用基于本征正交分解面向风电机组塔筒在线运行状态监测的数字孪生方法能有效地对物理实体的物理量(如应力、应变与位移等)进行实时的计算，且达到了较高的精度，进而实现对塔筒的运行状态进行在线监测。

本发明的有益效果如下：

本发明提供一种面向风电机组塔筒在线运行状态监测的数字孪生方法，相较于现有技术，该方法能快速地获得风电机组塔筒当前状态并以应力云图、应变曲线的方式进行了可视化，分析结果具有较高精度，达到在线监测塔筒运行状态的目的。

附图说明

图1是本发明技术实施方案路线图；

图2是风电机组塔筒变形示意图；

图3是奇异值对数分布曲线图；

图4各阶变换矩阵未保留原始样本的特征信息图；

图5是5阶降阶模型在各自由度上相对误差对数值图；

图6是简化塔筒物理实体示意图、物理实体图、数字孪生模型图；

图7是上位机软件界面图；

图8是简化塔筒最大von-mises应力在线监测计算图；

图9是简化塔筒应变在线监测计算图。

具体实施方式

以下结合说明书附图以及实施例，对本发明作进一步描述。

一种面向风电机组塔筒在线运行状态监测的数字孪生方法，具体步骤如下：

步骤1：建立风电机组塔筒的物理简化模型，将塔筒近似成一个一端有法兰结构的细长梁，梁截面为圆环形状，材料为聚氯乙烯(Polyvinyl Chloride，PVC)；使用3D8节点实体结构单元对简化模型划分网格，得到一个含有7080个节点与3780个单元的网格模型。

为计算塔筒无位移约束节点的节点位移向量q_a，对式(3)变形后取下式：

使用拉丁超立方采样方法对塔筒顶端处的挠度w与塔筒弯曲方位角进行采样，确定挠度的采样范围为[0,0.15]，方位角的采样范围为[0,2π)，采样点数量为20，结果如表1所示：

表1：塔筒顶端处挠度、塔筒弯曲方位角使用拉丁超立方采样方法获得的采样结果

将表1数据代入到式(2)获得位移边界条件，然后再代入式(13)进行计算，计算得到节点的位移解集合并将其用于步骤2用于构建降阶模型。本发明最终在线监测的物理量为应力与应变，可以使用节点的位移解参照式(14)计算任一单元节点处的应力与应变：

式中，σ(x₀,y₀,z₀)为坐标(x₀,y₀,z₀)处单元节点的应力；D为材料弹性系数矩阵；ε(x₀,y₀,z₀)为坐标(x₀,y₀,z₀)处单元节点的应变；为几何方程的算子矩阵；N(x₀,y₀,z₀)为坐标(x₀,y₀,z₀)处形状函数矩阵；σ_xx、σ_yy、σ_zz、σ_xy、σ_yz、σ_zx分别为单元节点的应力分量，ε_xx、ε_yy、ε_zz、ε_xy、ε_yz、ε_zx分别为单元节点的应变分量。

步骤2：按照表1确定位移边界条件并求解有限元模型，将计算得到的节点的位移解对应采样点序号排序组成矩阵V＝[v₁,v₁,...,v₂₀]，对矩阵V使用奇异值分解，得到奇异值矩阵Σ和左奇异值向量矩阵W。

将奇异值矩阵Σ对角线上的前20个奇异值绘制在图3中，因数值大小差异较大，故对数值取以10为底的对数后进行绘制，并按照奇异值在矩阵Σ中所在的位置给定序号，给定对角线上第一个奇异值序号为1，以此类推；分析图3中所得的结果，可以知道序号6到20的奇异值大小较为相近。选取前k个奇异值对应的左奇异值向量构建变换矩阵，使用变换矩阵构建降阶模型，并使用下式计算变换矩阵未保留原始样本的特征信息U(k)：

U(k)＝log₁₀((1-I(k))×100％) (15)

当k小于10时，变换矩阵未保留原始样本的特征信息如图4所示。由图4可知，当k取值大于等于5后，变换矩阵的未保留原始样本的特征信息变化幅度较小；通过计算可知降阶阶数为5时，变换矩阵未保留原始样本0.0001％的特征信息，因此初步确定降阶阶数为5阶。

选取左奇异值向量矩阵W的前五列构建变化矩阵Φ₅，使用该变化矩阵对有限元模型进行降阶得到5阶降阶模型，将式(5)代入式(13)，并对式(13)两端左乘得到下式：

式中，b为系数向量，长度为降阶阶数，按照式(5)重构可得位移向量。

在相同的位移边界条件下，分别使用降阶模型与有限元模型进行计算得到节点的位移解，计算降阶模型在节点的位移解中各自由度上的相对误差，由于大部分相对误差的数值大小较为接近，为了清楚地表示降阶模型的相对误差，采用相对误差对数值，定义如下：

式中，err_i为降阶模型在第i个自由度上的相对误差对数值；q_i为有限元模型计算得到的节点的位移解中第i个自由度对应的位移值，q_i(5)为5阶降阶模型计算得到的节点的位移解中第i个自由度对应的位移值。对节点的位移解中所有自由度进行计算，并将计算结果绘制于图5，从图5中可以知道相对误差对数值集中在-10附近，相对误差对数值最大不超过-7。通过计算可知最大的相对误差为0.82×10^-7，满足工程实际应用的需求，因此将该塔筒的降阶阶数确定为5阶。

步骤3：塔筒结构示意图如图6(a)所示物理实体总长度为1m，在塔筒的顶端处安装一倾角传感器WT61C，倾角传感器可以检测塔筒顶端端面的倾角；在塔身上距离地面0.5m位置上安装一应变传感器BF350-3EB，应变传感器可测量塔筒安装位置处的平均应变并输出电压；使用STM32单片机采集应变传感器与倾角传感器的数据，通过单片机自带的串行通信端口将采集的数据上传到计算机中，计算机处理数据并计算塔筒的应力与应变。最终搭建的塔筒物理实体如图6(b)所示，数字孪生模型如图6(c)所示。

使用倾角传感器检测塔筒顶端端面的倾角，结合式(1)与式(2)得到位移边界条件；使用模拟数字转换器(Analog-to-digital Converter，ADC)采样读取应变传感器的电压值，根据式(18)可以计算物理实体在应变传感器安装位置处的平均应变。

式中，e₀为应变传感器的输出电压，单位为V；E为应变传感器的供电电压，单位为V；R₁、R₂、R₃、R₄为应变传感器的电桥电阻，单位为Ω；K_s为应变传感器的灵敏系数，与应变传感器使用的材料相关，BF350-3EB应变传感器的灵敏系数为2.1；ε为物理实体在应变传感器安装位置的平均应变。

基于下式可以通过单元节点处的六个应力分量计算出该节点处的von-mises应力：

式中，σ_xx、σ_yy、σ_zz、σ_xy、σ_yz、σ_zx分别为单元节点处应力的六个应力分量，其相关的定义已在式(14)给出。

塔筒物理实体与计算机上的降阶模型交互分析时，首先使用STM32单片机读取倾角传感器与应变传感器的数据；获取数据后通过STM32单片机将将数据上传到计算机中；在已经编写完成的计算机上位机软件中对倾角数据进行处理转换为位移边界条件，将位移边界条件代入已在上位机部署完成的5阶降阶模型中求解得到节点的位移解，上位机软件界面如图7所示，界面分为三个区块，分别为塔筒三维模型显示区、实时数据显示区、通信以及数据显示控制功能区；通过节点的位移解计算各节点处的von-mises应力，使用三维云图将结果绘制在塔筒三维模型显示区中，并在数据显示区中实时显示最大的节点von-mises应力值，结果如图8所示；使用节点的位移解计算塔筒物理实体在应变传感器安装处的平均应变，并使用节点的位移解在塔筒三维模型显示区中绘制变形后的塔筒模型，将STM32单片机上传的应变传感器测量值和降阶模型计算值绘制在实时数据显示区的同一数据可视化控件中，同时计算降阶模型的相对误差，将相对误差显示在实时数据显示区的另一数据可视化控件，简化塔筒应变在线监测结果如图9所示。

从图8可以看出，最大von-mises应力曲线平缓变化且数据无间断、跳变，符合物理实际；从图9可以看出，降阶模型的应变计算值虽然与应变传感器的应变测量值有一定误差，但应变计算值的变化规律与应变测量值的变化规律一致。

综上所述，基于本征正交分解的数字孪生方法能有效的对风电机组塔筒的运行状态进行在线监测。

Claims (4)

1.一种面向风电机组塔筒在线运行状态监测的数字孪生方法，其特征在于，包括如下步骤：

步骤1：根据风电机组塔筒的几何特征获得其物理简化模型，分析塔筒实际受力情况，构建塔筒有限元模型；

步骤2：使用本征正交分解POD对步骤1所构建的有限元模型进行降阶分析，对比不同阶数下降阶模型的精度，并确定最终的降阶阶数，最终降阶阶数应满足在该阶数下的降阶模型符合工程需求；

步骤3：在计算机中设计并构建上位机软件，并将塔筒的降阶模型部署到上位机软件中；搭建塔筒的物理实体并安装应变传感器与倾角传感器，使用单片机读取传感器数据并将数据上传到上位机中转换为位移边界条件，再将位移边界条件输入到降阶模型中进行计算；通过计算应变传感器的应变测量值与上位机中降阶模型的应变计算值之间的误差来判断用于塔筒运行状态监测的数字孪生方法的准确度。

2.根据权利要求1所述的一种面向风电机组塔筒在线运行状态监测的数字孪生方法，其特征在于，

步骤1所述流程为：简化风电机组塔筒间的法兰连接结构，将塔筒简化为梁结构，进而构建出塔筒的物理简化模型；

塔筒弯曲变形前的轴线与塔筒下端面的交点为坐标系XYZ的原点O；塔筒弯曲变形前与弯曲变形后塔筒顶端端面的法向量分别为n₀，n₁，两向量的夹角为塔筒顶端处端面的倾角θ，塔筒顶端处的挠度w与端面倾角θ的表达式为：

式中，F为塔筒顶端施加的集中载荷，单位为N；E为塔筒材料的弹性模量，单位为Pa；I为塔筒截面惯性矩，单位为m⁴；L为塔筒长度，单位为m；

在三维空间描述塔筒弯曲变形需要确定挠度w和方位角其中方位角/>是用于表达塔筒的弯曲方位；因此塔筒有限元模型的位移边界条件由w与/>这两个参数确定，根据线性叠加原理与力的分解确定塔筒位移边界条件，得到：

式中，q_bx＝[1,0,1,0,...,1,0]^T，q_by＝[0,1,0,1,...,0,1]^T，向量长度由具体位移边界条件决定；q_b为塔筒的位移边界条件决定的节点位移向量，即塔筒有位移约束的节点对应的节点位移向量；w_x为总挠度在X方向的分量；w_y为总挠度在Y方向的分量；

将塔筒顶端处的挠度w转换为位移边界条件，并构建出塔筒的有限元模型，位移边界条件下塔筒有限元模型的平衡方程如式(3)所示：

式中，K_aa为塔筒无位移约束节点对应的刚度矩阵；K_ab、K_ba为塔筒位移约束节点与无位移约束节点之间的耦合矩阵；K_bb为塔筒有位移约束节点对应的刚度矩阵；q_a为塔筒无位移约束节点的节点位移向量；g_a、g_b分别为塔筒无位移约束、有位移约束节点对应的节点重力载荷向量；R_r为塔筒有位移约束节点对应的节点支反力向量。

3.根据权利要求1所述的一种面向风电机组塔筒在线运行状态监测的数字孪生方法，其特征在于，

步骤2所述过程为：将本征正交分解与有限元方法结合构建塔筒位移场降阶模型，包括：

求解步骤1的有限元模型，得到节点的位移解集合，从集合中选取m个样本组成矩阵X＝[q₁,q₂,...,q_m]；其中q_i＝[u₁,u₂,...,u_z]^T是节点的位移解，z是有限元模型的自由度数量；设{ζ₁,ζ₂,...,ζ_n}为X的一组数量为n的标准正交基，将其写为矩阵形式记为Φ；对于塔筒的有限元模型而言，任意时刻的节点位移值q_t由{ζ₁,ζ₂,...,ζ_n}表示：

式中，

选择{ζ₁,ζ₂,...,ζ_n}中前k个标准正交基向量来表示qt有：

式中，Φ_k＝[ζ₁,ζ₂,...,ζ_k]，k＜n；b_t(k)＝[b_t1,b_t1,...,b_tk]^T；

令误差函数为：

为获得误差函数在约束条件{ζ₁,ζ₂,...,ζ_n}为标准正交基下的最优解，引入Lagrange因子u_ij(i,j＝k+1,k+2,…,n)，构造Lagrange函数，寻找最优解对应的标准正交基；

式中，δ_ij为克罗内克函数；

将式(7)两端对ζ_j求偏导，得到：

式中，Φ_n-k＝[ζ_k+1,ζ_k+2,...,ζ_n]^T，u_j＝[u_k+1,j,u_k+2,j,...,u_n,j]^T；

将上式写为矩阵形式得到下式：

式中，U_n-k＝[u_k+1,u_k+2,...,u_n]^T；

为求最优解，令式(9)等于0，并对其左乘由于矩阵Φ_n内部向量标准正交，因此得到/>易知，U_n-k为半正定矩阵，故存在正交矩阵P可将U_n-k对角化，得到/>将上式两端左乘Φ_n-kP，得到下式：

式中，Λ为U_n-k的对角矩阵；

将Λ的对角元素降序排列，从式(10)可知Λ的第i个对角元素为XX^T的特征值λ_i，而矩阵Φ_n-kP的第i列是XX^T对应于λ_i的特征向量；

考虑到正交矩阵在Frobenius范数下是保范的，因此式(6)可重写为：

取X前n项奇异值对应的左奇异值向量组成变换矩阵Φ_n，此时Φ_n为误差函数在约束条件{ζ₁,ζ₂,...,ζ_n}为标准正交基下的最优解，误差函数有最小值为矩阵XX^T最后n-k个特征值的和；

定义下式：

式中，λ_i为XX^T按照降序排列的特征值；

若I(k)≥d％则称k维基向量保留有原始样本d％的特征信息；构建不同阶数下的降阶模型，并分别计算不同降阶阶数下降阶模型的精度，确定合适的降阶阶数并得到降阶模型，最终确定的降阶阶数应满足该阶数下的降阶模型符合实际工程需求。

4.根据权利要求1所述的一种面向风电机组塔筒在线运行状态监测的数字孪生方法，其特征在于，

步骤3中，上位机实现读取单片机上传的传感器数据，将数据转换为位移边界条件，并输入到降阶模型中实现快速计算，得到当前塔筒物理实体的应变与应力计算值，以及节点的位移解；降阶模型完成计算后，在上位机中对计算结果进行处理，即基于节点位移解快速地量化塔筒物理实体变形并在上位机中进行显示；通过降阶模型计算得到的应变值与传感器测量的应变值得到降阶模型的应变计算结果与物理实体真实应变之间的误差；结合节点的位移解与应力计算值实现应力云图的可视化，直观并清楚地表达出塔筒物理实体的应力分布情况。