CN117108445B - 一种串列式双风轮风电机组数字孪生仿真方法 - Google Patents
一种串列式双风轮风电机组数字孪生仿真方法 Download PDFInfo
- Publication number
- CN117108445B CN117108445B CN202310917715.2A CN202310917715A CN117108445B CN 117108445 B CN117108445 B CN 117108445B CN 202310917715 A CN202310917715 A CN 202310917715A CN 117108445 B CN117108445 B CN 117108445B
- Authority
- CN
- China
- Prior art keywords
- wind
- wind wheel
- tower
- wheel
- load
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 82
- 238000004088 simulation Methods 0.000 title claims abstract description 49
- 238000004458 analytical method Methods 0.000 claims abstract description 39
- 238000004364 calculation method Methods 0.000 claims abstract description 39
- 230000003993 interaction Effects 0.000 claims abstract description 36
- 238000013528 artificial neural network Methods 0.000 claims abstract description 17
- 230000008878 coupling Effects 0.000 claims abstract description 12
- 238000010168 coupling process Methods 0.000 claims abstract description 12
- 238000005859 coupling reaction Methods 0.000 claims abstract description 12
- 238000012544 monitoring process Methods 0.000 claims abstract description 12
- 238000013461 design Methods 0.000 claims abstract description 8
- 230000010354 integration Effects 0.000 claims abstract description 7
- 230000005540 biological transmission Effects 0.000 claims description 39
- 239000013598 vector Substances 0.000 claims description 38
- 230000033001 locomotion Effects 0.000 claims description 35
- 230000004044 response Effects 0.000 claims description 29
- 238000005452 bending Methods 0.000 claims description 24
- 238000006073 displacement reaction Methods 0.000 claims description 16
- 230000006870 function Effects 0.000 claims description 13
- 238000010801 machine learning Methods 0.000 claims description 13
- 238000010008 shearing Methods 0.000 claims description 13
- 230000008859 change Effects 0.000 claims description 11
- 230000009471 action Effects 0.000 claims description 9
- 230000001133 acceleration Effects 0.000 claims description 8
- 230000009467 reduction Effects 0.000 claims description 8
- 238000013179 statistical model Methods 0.000 claims description 6
- 238000013016 damping Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 238000011161 development Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 230000004048 modification Effects 0.000 claims description 3
- 238000012986 modification Methods 0.000 claims description 3
- 230000000306 recurrent effect Effects 0.000 claims description 3
- 238000011410 subtraction method Methods 0.000 claims 2
- 238000012360 testing method Methods 0.000 abstract description 2
- 238000004904 shortening Methods 0.000 abstract 1
- 230000008569 process Effects 0.000 description 6
- 238000010248 power generation Methods 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 230000003595 spectral effect Effects 0.000 description 4
- KFLQGJQSLUYUBF-WOJBJXKFSA-N Phyllanthin Chemical compound C([C@H](COC)[C@@H](COC)CC=1C=C(OC)C(OC)=CC=1)C1=CC=C(OC)C(OC)=C1 KFLQGJQSLUYUBF-WOJBJXKFSA-N 0.000 description 3
- KFLQGJQSLUYUBF-PMACEKPBSA-N Phyllanthin Natural products C([C@@H](COC)[C@H](COC)CC=1C=C(OC)C(OC)=CC=1)C1=CC=C(OC)C(OC)=C1 KFLQGJQSLUYUBF-PMACEKPBSA-N 0.000 description 3
- 229910052799 carbon Inorganic materials 0.000 description 3
- 125000004122 cyclic group Chemical group 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 239000000853 adhesive Substances 0.000 description 1
- 230000001070 adhesive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000004215 lattice model Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F03—MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
- F03D—WIND MOTORS
- F03D7/00—Controlling wind motors
- F03D7/02—Controlling wind motors the wind motors having rotation axis substantially parallel to the air flow entering the rotor
- F03D7/04—Automatic control; Regulation
- F03D7/042—Automatic control; Regulation by means of an electrical or electronic controller
- F03D7/043—Automatic control; Regulation by means of an electrical or electronic controller characterised by the type of control logic
- F03D7/046—Automatic control; Regulation by means of an electrical or electronic controller characterised by the type of control logic with learning or adaptive control, e.g. self-tuning, fuzzy logic or neural network
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F03—MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
- F03D—WIND MOTORS
- F03D17/00—Monitoring or testing of wind motors, e.g. diagnostics
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F03—MACHINES OR ENGINES FOR LIQUIDS; WIND, SPRING, OR WEIGHT MOTORS; PRODUCING MECHANICAL POWER OR A REACTIVE PROPULSIVE THRUST, NOT OTHERWISE PROVIDED FOR
- F03D—WIND MOTORS
- F03D7/00—Controlling wind motors
- F03D7/02—Controlling wind motors the wind motors having rotation axis substantially parallel to the air flow entering the rotor
- F03D7/04—Automatic control; Regulation
- F03D7/042—Automatic control; Regulation by means of an electrical or electronic controller
- F03D7/043—Automatic control; Regulation by means of an electrical or electronic controller characterised by the type of control logic
- F03D7/045—Automatic control; Regulation by means of an electrical or electronic controller characterised by the type of control logic with model-based controls
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/06—Wind turbines or wind farms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E10/00—Energy generation through renewable energy sources
- Y02E10/70—Wind energy
- Y02E10/72—Wind turbines with rotation axis in wind direction
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Chemical & Material Sciences (AREA)
- Software Systems (AREA)
- Sustainable Energy (AREA)
- Combustion & Propulsion (AREA)
- Mechanical Engineering (AREA)
- Sustainable Development (AREA)
- Life Sciences & Earth Sciences (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Artificial Intelligence (AREA)
- Algebra (AREA)
- Fuzzy Systems (AREA)
- Computational Mathematics (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Computer Graphics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Medical Informatics (AREA)
- Wind Motors (AREA)
Abstract
本发明涉及风电机组仿真测试技术领域，具体是指一种串列式双风轮风电机组数字孪生仿真方法。包括如下步骤：S1：多模块集成建立数字孪生仿真架构；S2：前后风轮串列连接反向旋转，耦合安装在同一个塔筒上；S3：由入流风模块实现前后风轮风速模拟与计算；S4：由空气动力模块实现前后风轮气动性能分析；S5：由结构动力学模块实现载荷计算和耦合；S6：由空气动力学与结构动力学模块结合实现前后风轮气动载荷干涉；S7：由控制模块实现协调控制，多控制器间可进行信息交互；S8：使用神经网络动态补偿得到双风轮机组的数字孪生体。本发明提供一种能降低成本并缩短设计周期，实时监测关键部件状态，准确描述双风轮间气动载荷干涉特性的数字孪生仿真方法。
Description
技术领域
本发明涉及风电机组仿真测试技术领域，具体是指一种串列式双风轮风电机组数字孪生仿真方法。
背景技术
随着双碳战略的提出，减少碳排放以应对气候变化逐步成为全球共识。风力发电中的传统风力机的迎风面积小，输出功率低，一般采取扩大扫掠面积的方法提高其输出功率，但过长的叶片增加了额外的成本。新型双风轮风电机组可通过风能梯次利用提高风能捕获效率并降低风电度电成本，有望为双碳目标的实现做出重大贡献。
由于风电物理样机的研制周期长，成本高，因此采用数字化仿真技术可以降低成本并缩短设计周期，但如何描述双风轮之间的干涉特性，实现高精度的非线性气动-弹性耦合，是一个关键问题，同时，由于传感器的高制造和安装成本，如何实现对风机关键部件的高精度监测也是一个关键问题；
现有技术一：
建立双风轮风机的简化模型，简化或忽略某些部件(如塔筒、前后风轮的传动系统等)，基于所建立的简化模型计算风机载荷，不能准确描述发电过程中的特性；由于传感器较高的安装和制造成本，在某些关键部件，无法实现运行状态的实时监测。
现有技术二：
将前后风轮作为单独运行的个体，或仅考虑前风轮影响后风轮的输入风速，而未考虑前后风轮的干涉特性，不能计算双风轮风电机组运行过程中前后风轮对塔筒或另一风轮的影响。
发明内容
本发明要解决的技术问题是提供一种能够降低成本并缩短设计周期，实时监测关键部件状态，准确描述发电过程中前后风轮对塔筒或另一风轮的影响，描述双风轮之间气动载荷干涉特性的一种串列式双风轮风电机组数字孪生仿真方法。
本发明所采用的技术方案为：一种串列式双风轮风电机组数字孪生仿真方法，包括如下步骤：
S1：一种串列式双风轮风电机组数字孪生仿真方法，包含前后风轮风速模拟与计算、前后风轮气动载荷计算、塔筒结构载荷计算和耦合、前后风轮气动载荷干涉、前后风轮协调控制、整机关键部件状态监测、误差动态补偿功能，各功能由入流风、空气动力学、结构动力学和控制模块有机结合实现，多模块集成建立双风轮风电机组的数字孪生仿真架构；
S2：所述串列式双风轮风电机组由塔筒、前后风轮、前后传动轴、前后发电机、轮毂组成，前后风轮串列连接反向旋转，经各自的传动轴连接至发电机，耦合安装在同一个塔筒上；
S3：所述前后风轮风速模拟与计算由入流风模块实现，其中前风轮输入风速使用统计模型模拟得到矢量时间序列，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速使用尾迹计算模型求解方法处理后得到；
S4：所述气动载荷计算由空气动力学模块实现，采用气动特性精细化建模方法进行前后风轮的气动性能分析；
S5：所述结构载荷计算和耦合由结构动力学模块实现，将风机结构认为刚体与柔体结合的多体系统，采用动力学特性分析和谐响应分析(如模态叠加法、完全法、缩减法)对风机系统进行建模；
S6：所述前后风轮气动载荷干涉结合空气动力学与结构动力学模块实现，前后风轮的气动载荷影响其传动系统和塔筒，塔筒的结构动力学载荷和响应传递给另一风轮以影响其气动载荷；
S7：所述协调控制功能由控制模块实现，控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；
S8：所述整机关键部件状态监测和误差动态补偿由结构动力模块实现，使用虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿，得到双风轮风电机组的数字孪生体。
作为本发明进一步的方案：所述步骤S1中，多模块集成建立了双风轮风电机组的仿真架构；调用入流风、空气动力学、控制和结构动力学进行多领域、全工况的仿真；通过入流风模块模拟或计算得到前后风轮的输入风速，传递给空气动力学模块计算得到前后风轮和塔筒的气动载荷，计算得到的风轮气动推力叠加在塔筒上，通过结构动力学模块计算塔筒的载荷和结构运动信息，塔筒载荷反馈给前风轮改变其气动载荷，塔筒结构运动信息传递给后风轮改变其输入风速，进一步改变其气动力和气动载荷，控制模块可随时读取机组的运行状态，并进行前后风轮控制器之间的信息交互，通过虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿，实现了双风轮风电机组的数字孪生。
作为本发明进一步的方案：所述步骤S2中，双风轮风电机组主要包括以下部件：前后风轮、塔筒、前后传动轴、前后发电机、轮毂；串列式双风轮风电机组的物理结构为：前后风轮水平排列、串列连接、对向旋转，并共用一个塔架结构，气动转矩经各自的传动轴传递给发电机产生电能；塔筒载荷同时受到前后风轮作用，风轮运动时均受经塔筒传递的另一风轮的干涉影响，前后传动轴和发电机也受到经风轮传递的影响。
作为本发明进一步的方案：所述步骤S3中，前后风轮风速模拟与计算由入流风模块实现；前风轮输入风速由统计模型模拟得到，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速由尾迹计算模型求解得到；
使用流场计算方法中的涡流理论方法(如自由涡尾迹FVW、预定尾迹方法、刚性尾迹方法)得到前风轮的尾流；
首先将叶片中弧面沿弦向和展向分别划分网格，在每个网格的1/4弦向宽度处沿展向放置附着涡线(boundvortexline)，该涡线的强度代表了对应网格展向分布涡片的涡量总和；在附着涡线两端，沿弦向往下游拖出半无限长的自由涡线(trailingvortexline)，其强度与附着涡强度相同；该自由涡线可进一步分为两个部分：叶片部分和尾迹部分；
叶片上拖出的自由涡线经过一段时间的发展，将分别向叶根和叶尖方向卷起而形成叶根和叶尖涡，叶根处拖出的尾涡很快消散，但叶尖涡的影响一直存在；基于这种情况，同时也为了降低计算量，将上述尾迹部分的自由涡进一步分成两个部分：近场尾迹和远场叶尖涡线，近场尾流的长度常为一到两倍的风轮直径，故后风轮位于近场尾流区中，选取前风轮尾流中的近场尾流速度v1设置为后风轮的初始输入风速，塔筒振动造成的后风轮输入风速变化量为vt，故后风轮的最终输入风速Vrear为：
Vrear＝V1+Vt
在相同或相邻的两个时间步长内，前风轮首先计算其尾流，计算完成后处理并通过glue-code传递给后风轮，作为其输入风速。
作为本发明进一步的方案：所述步骤S4中，采用气动特性精细化建模方法进行前后风轮的气动性能分析；前后风轮叶片的气动载荷均采用叶素-动量理论(BEM)计算；风轮是能量转换的起始段，捕获来流风的动能并转化为机械能，不同段上的风机翼形存在差异，需要根据设计要求对风机叶片进行分段设计；综合考量叶片的分布式特性参数(如弦长、相对厚度、密度、刚度和气动扭转角)，沿径向从叶根到叶尖将叶片划分为诸多叶素，在某一风速下计算各叶素所受升力和阻力，计算得到叶片所受整体力矩。
作为本发明进一步的方案：所述S5中，采用动力学特性分析和谐响应分析(如模态叠加法、完全法、缩减法)对风机系统进行建模；
前后风轮叶片的塔筒和叶片的结构动力学响应用谐响应分析(如模态叠加法、完全法、缩减法)进行描述，与动力学特性分析方法结合建立其动力学方程；叶片和塔筒简化为两个自由度的柔性悬臂梁，叶片与刚性的轮毂连接，塔筒与大地连接；基于正交单位向量建立了坐标系统以定义刚性结构参考系，各个部件坐标系之间可以由坐标转移矩阵实现转化，根据塔筒的分布式特性参数(如直径、密度和刚度)，沿塔架的中心从下到上用n个广义坐标描述，风力机上任何一点的受力情况和位移都可以由公式来进行表示；
前后风轮的传动链被模拟成连接轮毂和发电机的等效扭转弹簧系统，有线性的扭转弹性和阻尼，同时对以下部件做了简化近似处理：发电机、机舱。
作为本发明进一步的方案：所述步骤S6中，风轮气动载荷干涉结合空气动力学与结构动力学模块实现：
首先，根据前后风轮入流风速计算各风轮的初步气动推力Ff和Fr，计算入流风作用而产生的塔筒初步气动载荷Lt1(包括受入流风作用产生的塔顶剪切力Ftt、弯矩Btt)，将前后风轮的初步气动推力传递至塔筒。
其次，塔筒受前后风轮气动推力影响产生塔筒次生载荷Lt2(包括次生塔顶剪切力Ftb、次生弯矩Btb)，在初步气动载荷和次生载荷的作用下产生结构响应，由结构动力学模块计算以下结构运动信息：整体塔顶剪切力Ft、弯矩Bt以及位移xt、速度vt和加速度at；
最后，分别计算前后风轮受塔筒结构运动影响后的整体气动载荷Lfore和Lrear。
式中，前风轮的整体载荷Lfore受到自身的初步气动载荷、后风轮的气动载荷以及塔筒载荷影响；Fr *为塔筒运动导致后风轮输入风速变化后的气动力，塔筒运动与其初步载荷和次生载荷均有关，即后风轮也同时受到其入流风造成的初步气动载荷、前风轮气动载荷和塔筒载荷影响。
作为本发明进一步的方案：所述步骤S7中，所述控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；在单风轮机组中，DLL控制器通过glue-code中一个固定数组与机组进行信息交互，该数组的内容被提前定义，用于表征机组的各项信息；本平台中，前后风轮分别用两个DLL控制器进行控制，使用协调控制模块进行信息交互；
协调控制模块中对前后风轮的机组信息(如以下信息：实时的转矩、转速、功率、桨距角)进行耦合，改造机组的信息交互数组，给出了另一风轮的机组信息定义，这些信息通过glue-code进行传递；基于这种方法，前后风轮的DLL控制器可以随时读取另一风轮的状态信息，做到整机状态的实时监测；
协调控制模块加装了DLL控制器的交互通道，通过对信息交互数组的进一步改造，开辟了用于存放控制器交互数据的数组空间，按照机组信息传递的方式也可交换控制器信息；前后风轮的控制器既可分别读取自身风轮的信息并给出控制指令，也可以读取另一风轮的信息并给出控制指令，还可以进行信息交互；前风轮的控制器可以输出数据给同一个时间步长内的后风轮控制器，后风轮控制器也可以输出数据给下一个时间步的前风轮控制器。
作为本发明进一步的方案：所述步骤S8中，使用虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿；
所述虚拟传感器在结构动力学模块中实现，使用动力学特性分析和谐响应分析可得到机组任意一点处的动力学响应和载荷；所述关键部件包括但不限于：叶根、塔顶、塔底、传动轴，关键部件状态包括但不限于如下较难测量或测量成本较高的载荷：叶根弯矩、叶片净空距离、塔顶剪切力、塔顶位移、塔底弯矩；
使用机器学习方法(如递归神经网络、循环神经网络)建立多输入多输出神经网络补偿模型，实现对风机模型的全工况补偿；定义状态误差补偿向量ex(k)和输出误差补偿向量ey(k)为输出向量，误差补偿回归向量E(k)为输入向量；ey(k)＝y(k)-yr(k)，ex(k)＝x(k)-xr(k)，xr和yr为实际值；
E(k)＝[ex T(k-1),ex T(k-2),…,ex T(k-nex),ey T(k-1),ey T(k-2),…,ey T(k-ney),yT(k-1),yT(k-2),…,yT(k-ny),uT(k-1),uT(k-2),…,uT(k-nu)]T，nex、ney、ny、nu为各对应向量的自回归阶次；
使用得到的状态误差补偿向量和输出误差补偿向量，对机组状态和载荷进行动态补偿，实现了该机组的数字孪生：
其中，x(k)是机组状态向量，包括前后风轮转速、前后风轮气动推力、高速轴弯矩、低速轴弯矩；u(k)为控制量，包括前后风轮桨距角和前后发电机电磁转矩；w(k)为扰动项，包括前后风轮输入风速；y(k)为关键部件载荷；fx(k)和fy(k)为多输入多输出非线性函数，表示所建立的高精度复杂模型。
与现有技术相比，本发明的有益效果：
前后风轮的控制器既可分别读取自身风轮的信息并给出控制指令，也可以读取另一风轮的信息并给出控制指令，还可以进行信息交互塔筒同时受入流风和前后风轮的载荷作用振动变形产生位移、旋转或加速度结构运动信息，结构动力学模块输出塔筒的整体载荷和结构运动信息在结构动力学模块中输出的风轮转速、风轮转矩、桨距角及其他载荷信息均考虑了前后风轮的干涉影响；
整体来看，本发明能够降低成本并缩短设计周期，实时监测关键部件状态，准确描述发电过程中前后风轮对塔筒或另一风轮的影响，描述双风轮之间气动载荷干涉特性。
附图说明
图1为本发明一种串列式双风轮风电机组数字孪生仿真方法的串列式双风轮风电机组模型图。
图2为本发明一种串列式双风轮风电机组数字孪生仿真方法的仿真架构示意图。
具体实施方式
下面结合附图对本发明作进一步说明。
为了串列式双风轮风电机组数字孪生仿真方法能够降低成本并缩短设计周期，实时监测关键部件状态，准确描述发电过程中前后风轮对塔筒或另一风轮的影响，描述双风轮之间气动载荷干涉特性，本发明提供一种如图所示的一种串列式双风轮风电机组数字孪生仿真方法，包括如下步骤：
S1：一种串列式双风轮风电机组数字孪生仿真方法，包含前后风轮风速模拟与计算、前后风轮气动载荷计算、塔筒结构载荷计算和耦合、前后风轮气动载荷干涉、前后风轮协调控制、整机关键部件状态监测、误差动态补偿功能，各功能由入流风、空气动力学、结构动力学和控制模块有机结合实现，多模块集成建立双风轮风电机组的数字孪生仿真架构；
S2：所述串列式双风轮风电机组由塔筒、前后风轮、前后传动轴、前后发电机、轮毂部件组成，前后风轮串列连接反向旋转，经各自的传动轴连接至发电机，耦合安装在同一个塔筒上；
S3：所述前后风轮风速模拟与计算由入流风模块实现，其中前风轮输入风速使用统计模型模拟得到矢量时间序列，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速使用尾迹计算模型求解方法处理后得到；
S4：所述气动载荷计算由空气动力学模块实现，采用气动特性精细化建模方法进行前后风轮的气动性能分析；
S5：所述结构载荷计算和耦合由结构动力学模块实现，将风机结构认为刚体与柔体结合的多体系统，采用动力学特性分析和谐响应分析(如模态叠加法、完全法、缩减法)对风机系统进行建模；
S6：所述前后风轮气动载荷干涉结合空气动力学与结构动力学模块实现，前后风轮的气动载荷影响其传动系统和塔筒，塔筒的结构动力学载荷和响应传递给另一风轮以影响其气动载荷；
S7：所述协调控制功能由控制模块实现，控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；
S8：所述整机关键部件状态监测和误差动态补偿由结构动力模块实现，使用虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿，得到双风轮风电机组的数字孪生体。
本发明中所述步骤S1中，多模块集成建立了双风轮风电机组的仿真架构；调用入流风、空气动力学、控制和结构动力学多个模块进行多领域、全工况的仿真；通过入流风模块模拟或计算得到前后风轮的输入风速，传递给空气动力学模块计算得到前后风轮和塔筒的气动载荷，计算得到的风轮气动推力叠加在塔筒上，通过结构动力学模块计算塔筒的载荷和结构运动信息，塔筒载荷反馈给前风轮改变其气动载荷，塔筒结构运动信息传递给后风轮改变其输入风速，进一步改变其气动力和气动载荷，控制模块可随时读取机组的运行状态，并进行前后风轮控制器之间的信息交互，通过虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿，实现了双风轮风电机组的数字孪生。
本发明中所述步骤S2中，双风轮风电机组主要包括前后风轮、塔筒、前后传动轴、前后发电机、轮毂部件；串列式双风轮风电机组的物理结构为：前后风轮水平排列、串列连接、对向旋转，并共用一个塔架结构，气动转矩经各自的传动轴传递给发电机产生电能；塔筒载荷同时受到前后风轮作用，风轮运动时均受经塔筒传递的另一风轮的干涉影响，前后传动轴和发电机也受到经风轮传递的影响。
本发明中所述步骤S3中，前后风轮风速模拟与计算由入流风模块实现；
前风轮输入风速由统计模型模拟得到，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速由尾迹计算模型求解得到；
定常风模型是最简单的一种风速模型，风速和风向都不随时间的变化而变化；定常风模型可以由以下参数来表示：水平风速、垂直风速、风方向、水平剪切、垂直剪切和剪切指数；
阶跃风模型的风向不随时间发生变化，风速在特定时间发生阶跃；阶跃风模型由以下参数来表示：阶跃时间、水平风速、垂直风速、风方向、水平剪切、垂直剪切和剪切指数；湍流风模型可模拟全流场随机风，其统计特征首先定义在频域下的风功率谱中，然后经过快速傅里叶逆变换转化时域风，根据空间相干函数来生成空间多个点位的三维湍流风，从而生成空间中各点的风速相关参数文件，采用功率谱模型(如Kaimal模型和Von Karman模型)对湍流风进行模拟；
使用流场计算方法中的涡流理论方法(如自由涡尾迹FVW、预定尾迹方法、刚性尾迹方法)得到前风轮的尾流；
首先将叶片中弧面沿弦向和展向分别划分多个网格，在每个网格的1/4弦向宽度处沿展向放置附着涡线(bound vortex line)，该涡线的强度代表了对应网格展向分布涡片的涡量总和；在附着涡线两端，沿弦向往下游拖出半无限长的自由涡线(trailingvortex line)，其强度与附着涡强度相同；该自由涡线可进一步分为两个部分：叶片部分和尾迹部分；
叶片上拖出的自由涡线经过一段时间的发展，将分别向叶根和叶尖方向卷起而形成叶根和叶尖涡，叶根处拖出的尾涡很快消散，但叶尖涡的影响一直存在；基于这种情况，同时也为了降低计算量，将上述尾迹部分的自由涡进一步分成两个部分：近场尾迹和远场叶尖涡线，近场尾流的长度常为一到两倍的风轮直径，故后风轮位于近场尾流区中，选取前风轮尾流中的近场尾流速度v1设置为后风轮的初始输入风速，塔筒振动造成的后风轮输入风速变化量为vt，故后风轮的最终输入风速Vrear为：
Vrear＝V1+Vt
在相同或相邻的两个时间步长内，前风轮首先计算其尾流，计算完成后处理并通过glue-code传递给后风轮，作为其输入风速。
本发明中所述步骤S4中，采用气动特性精细化建模方法进行前后风轮的气动性能分析；前后风轮叶片的气动载荷均采用叶素-动量理论(BEM)计算；风轮是能量转换的起始段，捕获来流风的动能并转化为机械能，不同段上的风机翼形存在差异，需要根据设计要求对风机叶片进行分段设计；综合考量叶片的分布式特性参数(如弦长、相对厚度、密度、刚度和气动扭转角)，沿径向从叶根到叶尖将叶片划分为诸多叶素，在某一风速下计算各叶素所受升力和阻力，计算得到叶片所受整体力矩。
本发明中所述S5中，采用动力学特性分析和谐响应分析(如模态叠加法、完全法、缩减法)对风机系统进行建模；
前后风轮叶片的塔筒和叶片的结构动力学响应用谐响应分析(如模态叠加法、完全法、缩减法)进行描述，与动力学特性分析方法结合建立其动力学方程；叶片和塔筒简化为两个自由度的柔性悬臂梁，叶片与刚性的轮毂连接，塔筒与大地连接；基于正交单位向量建立了坐标系统以定义刚性结构参考系，各个部件坐标系之间可以由坐标转移矩阵实现转化，根据塔筒的分布式特性参数(如直径、密度和刚度)，沿塔架的中心从下到上用n个广义坐标描述，风力机上任何一点的受力情况和位移都可以由公式来进行表示；
前后风轮的传动链被模拟成连接轮毂和发电机的等效扭转弹簧系统，有线性的扭转弹性和阻尼，同时对发电机、机舱部件做了简化近似处理。
本发明中所述步骤S6中，风轮气动载荷干涉结合空气动力学与结构动力学模块实现：
首先，根据前后风轮入流风速计算各风轮的初步气动推力Ff和Fr，计算入流风作用而产生的塔筒初步气动载荷Lt1(包括受入流风作用产生的塔顶剪切力Ftt、弯矩Btt)，将前后风轮的初步气动推力传递至塔筒。
其次，塔筒受前后风轮气动推力影响产生塔筒次生载荷Lt2(包括次生塔顶剪切力Ftb、次生弯矩Btb)，在初步气动载荷和次生载荷的作用下产生结构响应，由结构动力学模块计算整体塔顶剪切力Ft、弯矩Bt以及位移xt、速度vt和加速度at结构运动信息
最后，分别计算前后风轮受塔筒结构运动影响后的整体气动载荷Lfore和Lrear。
式中，前风轮的整体载荷Lfore受到自身的初步气动载荷、后风轮的气动载荷以及塔筒载荷影响；Fr *为塔筒运动导致后风轮输入风速变化后的气动力，塔筒运动与其初步载荷和次生载荷均有关，即后风轮也同时受到其入流风造成的初步气动载荷、前风轮气动载荷和塔筒载荷影响。
本发明中所述步骤S7中，所述控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；在单风轮机组中，DLL控制器通过glue-code中一个固定数组与机组进行信息交互，该数组的内容被提前定义，用于表征机组的各项信息；本平台中，前后风轮分别用两个DLL控制器进行控制，使用协调控制模块进行信息交互；
协调控制模块中对前后风轮的机组信息(如实时的转矩、转速、功率、桨距角信息)进行耦合，改造机组的信息交互数组，给出了另一风轮的机组信息定义，这些信息通过glue-code进行传递；基于这种方法，前后风轮的DLL控制器可以随时读取另一风轮的状态信息，做到整机状态的实时监测；
协调控制模块加装了DLL控制器的交互通道，通过对信息交互数组的进一步改造，开辟了用于存放控制器交互数据的数组空间，按照机组信息传递的方式也可交换控制器信息；前后风轮的控制器既可分别读取自身风轮的信息并给出控制指令，也可以读取另一风轮的信息并给出控制指令，还可以进行信息交互；前风轮的控制器可以输出数据给同一个时间步长内的后风轮控制器，后风轮控制器也可以输出数据给下一个时间步的前风轮控制器。
所示步骤S8中，使用虚拟传感器实时计算关键部件状态，通过机器学习方法对计算结果进行神经网络动态补偿。
所述虚拟传感器在结构动力学模块中实现，使用动力学特性分析和谐响应分析可得到机组任意一点处的动力学响应和载荷。所述关键部件包括但不限于：叶根、塔顶、塔底、传动轴，关键部件状态包括但不限于如下较难测量或测量成本较高的载荷：叶根弯矩、叶片净空距离、塔顶剪切力、塔顶位移、塔底弯矩。
使用机器学习方法(如递归神经网络、循环神经网络)建立多输入多输出神经网络补偿模型，实现对风机模型的全工况补偿。定义状态误差补偿向量ex(k)和输出误差补偿向量ey(k)为输出向量，误差补偿回归向量E(k)为输入向量。ey(k)＝y(k)-yr(k)，ex(k)＝x(k)-xr(k)，xr和yr为实际值；
E(k)＝[ex T(k-1),ex T(k-2),…,ex T(k-nex),ey T(k-1),ey T(k-2),…,ey T(k-ney),yT(k-1),yT(k-2),…,yT(k-ny),uT(k-1),uT(k-2),…,uT(k-nu)]T，nex、ney、ny、nu为各对应向量的自回归阶次。
使用得到的状态误差补偿向量和输出误差补偿向量，对机组状态和载荷进行动态补偿，实现了该机组的数字孪生：
其中，x(k)是机组状态向量，包括前后风轮转速、前后风轮气动推力、高速轴弯矩、低速轴弯矩；u(k)为控制量，包括前后风轮桨距角和前后发电机电磁转矩；w(k)为扰动项，包括前后风轮输入风速；y(k)为关键部件载荷；fx(k)和fy(k)为多输入多输出非线性函数，表示所建立的高精度复杂模型。
本技术主要分为入流风模块、控制模块、空气动力学模块和结构动力学模块部分。其中，入流风模块负责前风轮风速的输入；控制模块中，可针对前后风轮构建可交互的外部控制器，对各风轮桨距角和发电机的电磁转矩进行控制；空气动力学模块负责计算前后风轮和塔架的气动载荷；结构动力学模块定义了结构的建模选项和几何形状，包括塔架、前后传动轴和前后风轮部位，可得到叶片和塔架上任意一点的形变。
模块间的耦合是通过模块化接口和耦合器(即glud-code)实现。载荷和响应通过粘合代码在结构动力学模块、空气动力学模块和控制模块之间传递，可以实现任意时刻的气动-弹性-伺服的交互。
如图2所示，建立仿真技术架构，首先由入流风速计前后风轮的初步气动载荷。入流风模块将输入风速传递给前后风轮和塔架，通过空气动力学模块计算其初步气动载荷，其中前风轮的尾流经计算后提取近场尾迹作为后风轮的输入风速。
其次计算塔筒受风轮传递载荷与自身气动载荷叠加而产生的整体载荷与结构运动。前后风轮上的气动载荷传递给结构动力学模块中双风轮的叶片，经各自的传动系统传递给塔筒，塔筒同时受入流风和前后风轮的载荷作用振动变形产生位移、旋转或加速度结构运动信息，结构动力学模块输出塔筒的整体载荷和结构运动信息；
然后计算风轮结合塔筒结构运动而产生的整体气动载荷。此时后风轮受到塔筒结构运动影响而产生振动与位移，改变其入流风产生的初步气动载荷，即前后风轮同时受到其入流风造成的气动力、塔筒气动力和另一风轮的气动力影响，在结构动力学模块中输出的风轮转速、风轮转矩、桨距角及其他载荷信息均考虑了前后风轮的干涉影响。
最后计算传动系统受各结构干涉影响的整体载荷。前后风轮的传动系统均受塔筒振动位移、前后风轮整体载荷和发电机电磁转矩的影响。即传动系统输入载荷考虑了前后风轮干涉的影响，将后风轮的运动状态纳入到前风轮输出转矩的计算过程中，后风轮同理；同时考虑了塔筒振动位移对齿轮箱产生的影响，改变了传动系统的输出载荷。传动系统将整体载荷传递给各自的发电机转子，并受到发电机电磁转矩的影响，转子转动切割磁感线产生电能，经机器学习方法进行动态误差补偿之后，得到数字孪生体。
本实例搭建2.7MW串列式双风轮风电机组的机组模型，其中前风轮容量为1.5MW，后风轮容量为1.2MW。两个风轮串列连接、对向旋转，并共用同一个塔架结构，组成串列式双风轮风电机组的物理结构。
如图1所示，首先搭建整体输入文件，说明全局仿真时间、仿真步长、所需调用的前后风轮fst配置文件、前后风轮间距信息；在前后风轮的fst配置文件中，分别对使用功能、输入文件和输出的配置做了说明，其中输入文件主要包括入流风、空气动力学、结构动力学、控制四个模块的配置文件；入流风模块中使用InflowWind.dat文件规定了前风轮输入风速的形式和风速网格，可使用BTS或WND风速文件输入风速，传递给空气动力学模块用于前风轮气动力的计算，后风轮的输入风速通过前风轮的尾迹计算得到；空气动力学模块中使用AeroDyn15.dat文件建立了塔架的精细化模型，调用Airfoil.dat与AeroDyn15_blade.dat规定了叶片的翼型和叶素的结构信息；在结构动力学模块中使用ElastoDyn.dat文件规定了叶片、传动、发电机、塔架和基础的自由度，并设置了这些部位的结构、质量和惯性，调用ElastoDyn_blade.dat和ElastoDyn_tower.dat对叶片和塔筒的结构、质量和惯性做了详细设置；在控制模块中调用ServoDyn.dat规定了所用控制器的形式并给出了Control.dll外部控制器的文件路径，其中外部控制器可使用基于Visual Studio编译生成的X64版本bladed-style DLL。
风速模型用于描述风速在空间尺度与时间尺度上的波动，所用Kaimal功率谱模型的功率谱密度表达式为：
其中，下标k代表速度分量的方向(k＝u,v,w)；Sk为功率谱密度；f为湍流风的频率；σk为风速分量标准差；Lk为积分尺度参数；vhub为轮毂高度出10分钟之内的平均风速。IEC61400-1标准规定Lk的取值为：
其中，Λk为湍流尺度参数，设风机轮毂高度为ht，则Λk可表示为：
Von Karman功率谱模型中风速u方向的功率谱密度表达式为：
湍流风在v,w方向的功率谱密度表达式为：
Von Karman功率谱模型中积分尺度参数L为：
L＝3.5Λu(1.20)
风速模型中的参考系主要采取惯性参考坐标系，z轴与塔筒中心线重合，正向竖直向上，x轴为正顺风方向，y轴为xz平面的垂直方向。
在势流理论中，常使用一系列基本解(如涡、源、偶极子)的叠加来模拟实际中具有复杂边界条件的流场。OLAF模型是基于拉格朗日方法，其中涡轮尾迹离散成拉格朗日标记。使用一种混合晶格/细丝方法表示尾迹。该方法在叶片参考面(如中弧面)上沿展向和弦向布置涡格或涡线以代替真实叶片，其中展向涡为叶片提供气动升力，弦向涡模拟了叶片上流动沿径向的变化；同时，从叶片尾缘向下游拖出自由涡系，模拟真实流动中的尾迹结构。拖出的自由涡从叶片尾缘一直延伸到下游无穷远处，将上述自由涡进一步分成两个部分：近场尾迹和远场叶尖涡线，使用近场尾迹作为后风轮初步输入风速，拖出的叶尖涡用一根自由涡线来表示。
用自由尾迹模型来描述叶尖涡线在空间中的运动轨迹。叶尖涡线被离散为一系列首尾相接的直线段，用控制点的位置rw描述涡线形状。涡线在运动中不受力的作用，是随当地流场运动的，那么叶尖涡控制点的运动方程可以写为：
其中，Vloc为当地流场的速度矢量，ψ为叶片方位角，ξ为尾迹寿命角。将叶尖涡的卷起点设在近场尾迹终止位置处，因此叶尖涡控制点的位置矢量在寿命角方向上存在边界条件，即寿命角ξ的取值范围为：
ξ∈[ξ0，∞] (1.22)
叶尖涡自由尾迹问题就是在给定的边界条件下，求解微分控制方程，最后得到控制点的位置矢量用以描述自由尾迹的几何结构。基于预测得到的尾迹形状，结合升力面涡格模型，即通过求解方程式可得到下一时刻的叶片附着涡环量预测值，同时更新近场尾迹涡和叶尖涡可得到相应的环量预测值。得到叶片环量后，结合来流速度、叶片旋转速度，即可得到叶片的有效近场尾流速度Vfvw。
风机的风轮是能量转换的起始段，捕获来凤的动能并转化为机械能。目前风电机组逐渐大型化，风机的叶片尺寸也逐渐增大，对风轮直接进行受力分析较为困难，故常采用BEM模型对风机叶片做精细化建模。
综合考量叶片的弦长、相对厚度和气动扭转角翼型参数，沿径向从叶根到叶尖将叶片划分为诸多长度为dr的微元，假设各微元之间的空气流动互不干扰，并且作用在所有微元上的力和力矩可以无损失地合成为叶片承受的合力和总力矩，则叶片所受的整体力矩可由叶片所有微元受力叠加得到。将这些微元称为叶素，前风轮划分得到40个叶素，后风轮划分得到31个叶素。
在前后风轮叶片的n个叶素中随机取一个距离叶根r米的叶素，如图2，设该叶素受到了相对速度为v的来流风的作用，产生了一个垂直于弦的气动力dF，将该气动力沿来流风的方向和与之垂直的方向分解，得到该叶素在来流风下的升力dL和阻力dD：
其中，ρ为空气密度，Cl、Cd分别为该叶素的升力系数和阻力系数，通过其翼型信息进行气动性能分析得到，l为该叶素的弦长，dr为该叶素的长度。
将气动力dF沿平行和垂直于风轮旋转平面的方向分解，得到局部风轮推力dFa和转矩dFu，用升力dL和阻力Dd来表示：
其中，φ为来流风与风轮旋转平面的夹角。
通过空气动力学方法计算出各微段的受力情况，则风轮叶片整体所受的气动力矩可由各微段受力叠加求得，其表达式为
式中，Fr为风轮总推力，Qr为风轮总转矩；ρ为空气密度；R为风轮半径，Rhub为轮毂半径；为该处气流倾角。
N个自由度的风机的结构动力学响应用假设模态法进行描述，使用Kane动力学方程建立其动力学方程：
F(r)+FR(r)＝0 (1.6)
其中，F(r)与FR(r)为风机系统所受的广义主动力和广义惯性力，广义惯性力FR:
R(r)＝FT(r)+FN(r)+FH(r)+FB(r)，r＝1，2，…，N (1.7)
其中，FT(r)、FN(r)、FH(r)、FB(r)分别为塔架、机舱、轮毂和叶片产生的广义惯性力，塔架和叶片的惯性力F根据下述公式得到：
其中ρ(l)、a(l)分别为沿塔架中心线和叶片轴向的线密度和线性加速度矢量。机舱和轮毂的惯性力F根据下述公式得到：
F＝v(-ma)+w(-H) (1.9)
其中，v、w和a为机舱和轮毂在惯性系中的偏线速度、偏角速度和线性加速度，m为机舱和轮毂的质量，H为机舱和轮毂在惯性系中绕质心的角动量的微分。
风机系统上广义主动力F(r)由作用在风机上的气动力Faero(r)、重力FG(r)、弹性力Felastic(r)组成：
F(r)＝Faero(r)+FG(r)+Felastic(r) (1.10)
其中，Faero(r)由下式计算得到：
F＝vF+wM (1.11)
其中，F和M分别为空气动力学模块中作用在叶片上的气动力和俯仰力矩。
叶片和塔筒简化为两个自由度的柔性悬臂梁，叶片与刚性的轮毂连接，塔筒与大地连接。基于正交单位向量建立了坐标系统以定义刚性结构参考系，各个部件坐标系之间可以由坐标转移矩阵实现转化，根据塔筒的分布式特性参数(如直径、密度和刚度)，沿塔架的中心从下到上用n个广义坐标描述，任意时刻悬臂梁上任意一点的变形U(z,t)可以写成：
其中，z为沿悬臂方向的位置，t为时刻，φi(z)是悬臂的振型，qi(t)为塔顶或叶根的广义坐标；mij、kij分别代表广义质量和广义刚度；c(t)是关于模态振型φi(z)的广义坐标，风力机上任何一点的受力情况和位移都可以由上述公式来进行表示。
传动链被模拟成连接轮毂和发电机的等效扭转弹簧系统，有线性的扭转弹性和阻尼。传动链的方程可写成如下方程：
Tshaft＝Kd(q1-q2)+Cd(q3-q4) (1.14)
式中Tshaft是低速轴上电机端的扭矩，Kd是传动链扭转刚度系数，Cd是传动链扭转阻尼系数。q1q2和q3q4分别为风轮和发电机的方位角及其一阶微分。
首先，根据前后风轮入流风速计算各风轮的初步气动推力Ff和Fr，计算入流风作用而产生的塔筒初步气动载荷Lt1(包括受入流风作用产生的塔顶剪切力Ftt、弯矩Btt)，在glue-code中将前后风轮的初步气动推力经各自的传动系统传递至塔筒。
其次，塔筒受前后风轮气动推力影响产生塔筒次生载荷Lt2(包括次生塔顶剪切力Ftb、次生弯矩Btb)，在初步气动载荷和次生载荷的作用下产生结构响应，由结构动力学模块计算整体塔顶剪切力Ft、弯矩Bt以及位移xt、速度vt和加速度at结构运动信息
最后，分别计算前后风轮受塔筒结构运动影响后的整体气动载荷Lfore和Lrear。
式中，前风轮的整体载荷Lfore受到自身的初步气动载荷、后风轮的气动载荷以及塔筒载荷影响；Fr *为塔筒运动导致后风轮输入风速变化后的气动力，塔筒运动与其初步载荷和次生载荷均有关，即后风轮也同时受到其入流风造成的初步气动载荷、前风轮气动载荷和塔筒载荷影响。
控制器之间的信息交互通过glud-code实现，在单风轮机组中，DLL控制器通过一个固定数组与机组进行信息交互，该数组的内容被提前定义，用于表征机组的各项信息。
本平台中两个风轮有自身的glud-code，同时有公共的缓存区以进行信息交互。为了使其中一个DLL控制器可以读取另一风轮的实时信息，在两个风轮的数组中提前定义另一风轮机组信息，该数组保存在glud-code的公共缓存区，其中的信息每台风轮都可以读取，数据可以在两台风轮中自由交换。为了使某一DLL控制器给另一DLL控制器传递信息，在上述数组中开辟了独立空间，不再对其进行提前定义，由用户进行自定义，可传递包括机组信息、控制指令、时间步在内的诸多信息。
以上对本发明及其实施方式进行了描述，这种描述没有限制性，附图中所示的也只是本发明的实施方式之一。总而言之如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。
Claims (4)
1.一种串列式双风轮风电机组数字孪生仿真方法，其特征在于：包括如下步骤：
S1：一种串列式双风轮风电机组数字孪生仿真方法，包含以下功能：前后风轮风速模拟与计算、前后风轮气动载荷计算、塔筒结构载荷计算和耦合、前后风轮气动载荷干涉、前后风轮协调控制、整机关键部件状态监测、误差动态补偿，各功能由入流风、空气动力学、结构动力学和控制模块有机结合实现，多模块集成建立双风轮风电机组的数字孪生仿真架构；
多模块集成建立了双风轮风电机组的仿真架构；调用入流风、空气动力学、控制和结构动力学进行多领域、全工况的仿真；通过入流风模块模拟或计算得到前后风轮的输入风速，传递给空气动力学模块计算得到前后风轮和塔筒的气动载荷，计算得到的风轮气动推力叠加在塔筒上，通过结构动力学模块计算塔筒的载荷和结构运动信息，塔筒载荷反馈给前风轮改变其气动载荷，塔筒结构运动信息传递给后风轮改变其输入风速，进一步改变其气动力和气动载荷，控制模块可随时读取机组的运行状态，并进行前后风轮控制器之间的信息交互；
S2：所述串列式双风轮风电机组由塔筒、前后风轮、前后传动轴、前后发电机、轮毂组成，前后风轮串列连接反向旋转，经各自的传动轴连接至发电机，耦合安装在同一个塔筒上；
S3：所述前后风轮风速模拟与计算由入流风模块实现，其中前风轮输入风速使用统计模型模拟得到矢量时间序列，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速使用尾迹计算模型求解方法处理后得到；
前后风轮风速模拟与计算由入流风模块实现；
前风轮输入风速由统计模型模拟得到，包括全风速范围内稳态风、阶跃风、湍流风，后风轮的输入风速由尾迹计算模型求解得到；
使用流场计算方法中的涡流理论方法得到前风轮的尾流，使用的涡流理论方法包括如自由涡尾迹FVW、预定尾迹方法、刚性尾迹方法；
首先将叶片中弧面沿弦向和展向分别划分为多个网格，在每个网格的1/4弦向宽度处沿展向放置附着涡线，该涡线的强度代表了对应网格展向分布涡片的涡量总和；在附着涡线两端，沿弦向往下游拖出半无限长的自由涡线，其强度与附着涡强度相同；该自由涡线可进一步分为两个部分：叶片部分和尾迹部分；
叶片上拖出的自由涡线经过一段时间的发展，将分别向叶根和叶尖方向卷起而形成叶根和叶尖涡，叶根处拖出的尾涡很快消散，但叶尖涡的影响一直存在；基于这种情况，同时也为了降低计算量，将上述尾迹部分的自由涡进一步分成两个部分：近场尾迹和远场叶尖涡线，近场尾流的长度常为一到两倍的风轮直径，故后风轮位于近场尾流区中，选取前风轮尾流中的近场尾流速度v1设置为后风轮的初始输入风速，塔筒振动造成的后风轮输入风速变化量为vt，故后风轮的最终输入风速Vrear为：
Vrear＝V1+Vt
在相同或相邻的两个时间步长内，前风轮首先计算其尾流，计算完成后处理并通过glue-code传递给后风轮，作为其输入风速；
S4：所述气动载荷计算由空气动力学模块实现，采用气动特性精细化建模方法进行前后风轮的气动性能分析；
S5：所述结构载荷计算和耦合由结构动力学模块实现，将风机结构认为刚体与柔体结合的多体系统，采用动力学特性分析和谐响应分析对风机系统进行建模，动力学特性分析和谐响应分析包括模态叠加法、完全法、缩减法；
采用动力学特性分析和谐响应分析对风机系统进行建模，动力学特性分析和谐响应分析包括模态叠加法、完全法、缩减法；
前后风轮叶片的塔筒和叶片的结构动力学响应用谐响应分析进行描述，结构动力学响应用谐响应分析包括模态叠加法、完全法、缩减法，与动力学特性分析方法结合建立其动力学方程；叶片和塔筒简化为两个自由度的柔性悬臂梁，叶片与刚性的轮毂连接，塔筒与大地连接；基于正交单位向量建立了坐标系统以定义刚性结构参考系，各个部件坐标系之间可以由坐标转移矩阵实现转化，根据塔筒的分布式特性参数，塔筒的分布式特性参数包括直径、密度和刚度，沿塔架的中心从下到上用n个广义坐标描述，风力机上任何一点的受力情况和位移都可以由公式来进行表示；
前后风轮的传动链被模拟成连接轮毂和发电机的等效扭转弹簧系统，有线性的扭转弹性和阻尼，同时对发电机、机舱做了简化近似处理；
S6：所述前后风轮气动载荷干涉结合空气动力学与结构动力学模块实现，前后风轮的气动载荷影响其传动系统和塔筒，塔筒的结构动力学载荷和响应传递给另一风轮以影响其气动载荷；风轮气动载荷干涉结合空气动力学与结构动力学模块实现：
首先，根据前后风轮入流风速计算各风轮的初步气动推力Ff和Fr，计算入流风作用而产生的塔筒初步气动载荷Lt1(包括受入流风作用产生的塔顶剪切力Ftt、弯矩Btt)，将前后风轮的初步气动推力传递至塔筒；
其次，塔筒受前后风轮气动推力影响产生塔筒次生载荷Lt2，塔筒次生载荷Lt2包括次生塔顶剪切力Ftb、次生弯矩Btb，在初步气动载荷和次生载荷的作用下产生结构响应，由结构动力学模块计算以下结构运动信息：整体塔顶剪切力Ft、弯矩Bt以及位移xt、速度vt和加速度at；
最后，分别计算前后风轮受塔筒结构运动影响后的整体气动载荷Lfore和Lrear；
式中，前风轮的整体载荷Lfore受到自身的初步气动载荷、后风轮的气动载荷以及塔筒载荷影响；Fr *为塔筒运动导致后风轮输入风速变化后的气动力，塔筒运动与其初步载荷和次生载荷均有关，即后风轮也同时受到其入流风造成的初步气动载荷、前风轮气动载荷和塔筒载荷影响
S7：所述协调控制功能由控制模块实现，控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；
S8：所述整机关键部件状态监测和误差动态补偿由结构动力模块实现，使用虚拟传感器实时计算关键部件状态，使用机器学习方法对计算结果进行神经网络动态补偿，得到双风轮风电机组的数字孪生体，使用虚拟传感器实时计算关键部件状态，使用机器学习方法对计算结果进行神经网络动态补偿；
所述虚拟传感器在结构动力学模块中实现，使用动力学特性分析和谐响应分析可得到机组任意一点处的动力学响应和载荷；所述关键部件包括但不限于：叶根、塔顶、塔底、传动轴，关键部件状态包括但不限于如下较难测量或测量成本较高的载荷：叶根弯矩、叶片净空距离、塔顶剪切力、塔顶位移、塔底弯矩；
使用机器学习方法建立多输入多输出神经网络补偿模型，机器学习方法包括递归神经网络、循环神经网络，实现对风机模型的全工况补偿；定义状态误差补偿向量ex(k)和输出误差补偿向量ey(k)为输出向量，误差补偿回归向量E(k)为输入向量；ey(k)＝y(k)-yr(k)，ex(k)＝x(k)-xr(k)，xr和yr为实际值；E(k)＝[ex T(k-1),ex T(k-2),…,ex T(k-nex),ey T(k-1),ey T(k-2),…,ey T(k-ney),yT(k-1),yT(k-2),…,yT(k-ny),uT(k-1),uT(k-2),…,uT(k-nu)]T，nex、ney、ny、nu为各对应向量的自回归阶次；
使用得到的状态误差补偿向量和输出误差补偿向量，对机组状态和载荷进行动态补偿，实现了该机组的数字孪生：
其中，x(k)是机组状态向量，包括前后风轮转速、前后风轮气动推力、高速轴弯矩、低速轴弯矩；u(k)为控制量，包括前后风轮桨距角和前后发电机电磁转矩；w(k)为扰动项，包括前后风轮输入风速；y(k)为关键部件载荷；fx(k)和fy(k)为多输入多输出非线性函数，表示所建立的高精度复杂模型。
2.根据权利要求1所述的一种串列式双风轮风电机组数字孪生仿真方法，其特征在于：所述步骤S2中，双风轮风电机组主要包括前后风轮、塔筒、前后传动轴、前后发电机、轮毂；串列式双风轮风电机组的物理结构为：前后风轮水平排列、串列连接、对向旋转，并共用一个塔架结构，气动转矩经各自的传动轴传递给发电机产生电能；塔筒载荷同时受到前后风轮作用，风轮运动时均受经塔筒传递的另一风轮的干涉影响，前后传动轴和发电机也受到经风轮传递的影响。
3.根据权利要求1所述的一种串列式双风轮风电机组数字孪生仿真方法，其特征在于：所述步骤S4中，采用气动特性精细化建模方法进行前后风轮的气动性能分析；前后风轮叶片的气动载荷均采用叶素-动量理论计算；风轮是能量转换的起始段，捕获来流风的动能并转化为机械能，不同段上的风机翼形存在差异，需要根据设计要求对风机叶片进行分段设计；综合考量叶片的分布式特性参数，所述分布式特性参数包括弦长、相对厚度、密度、刚度和气动扭转角，沿径向从叶根到叶尖将叶片划分为诸多叶素，在某一风速下计算各叶素所受升力和阻力，计算得到叶片所受整体力矩。
4.根据权利要求1所述的一种串列式双风轮风电机组数字孪生仿真方法，其特征在于：所述步骤S7中，所述控制模块包括前后风轮控制器，每个控制器可以读取整机的实时运行状态，控制器之间可进行信息交互；在单风轮机组中，DLL控制器通过glue-code中一个固定数组与机组进行信息交互，该数组的内容被提前定义，用于表征机组的各项信息；本平台中，前后风轮分别用两个DLL控制器进行控制，使用协调控制模块进行信息交互；
所述协调控制模块中对前后风轮的机组信息进行耦合，前后风轮的机组信息包括以下信息：实时的转矩、转速、功率、桨距角，改造机组的信息交互数组，给出了另一风轮的机组信息定义，这些信息通过glue-code进行传递；基于这种方法，前后风轮的DLL控制器可以随时读取另一风轮的状态信息，做到整机状态的实时监测；
所述协调控制模块加装了DLL控制器的交互通道，通过对信息交互数组的进一步改造，开辟了用于存放控制器交互数据的数组空间，按照机组信息传递的方式也可交换控制器信息；前后风轮的控制器既可分别读取自身风轮的信息并给出控制指令，也可以读取另一风轮的信息并给出控制指令，还可以进行信息交互；前风轮的控制器可以输出数据给同一个时间步长内的后风轮控制器，后风轮控制器也可以输出数据给下一个时间步的前风轮控制器。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310917715.2A CN117108445B (zh) | 2023-07-25 | 2023-07-25 | 一种串列式双风轮风电机组数字孪生仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310917715.2A CN117108445B (zh) | 2023-07-25 | 2023-07-25 | 一种串列式双风轮风电机组数字孪生仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117108445A CN117108445A (zh) | 2023-11-24 |
| CN117108445B true CN117108445B (zh) | 2024-05-03 |
Family
ID=88797440
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310917715.2A Active CN117108445B (zh) | 2023-07-25 | 2023-07-25 | 一种串列式双风轮风电机组数字孪生仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117108445B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117494618B (zh) * | 2024-01-03 | 2024-04-09 | 北京亚能电气设备有限公司 | 一种基于实时网络的风电场智能图像监控系统 |
| CN117993303B (zh) * | 2024-04-02 | 2024-06-11 | 上海勘测设计研究院有限公司 | 一种海上风电叶片变形的气动载荷校核方法 |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2020100825A4 (en) * | 2020-05-22 | 2020-07-02 | China University Of Mining And Technology, Yinchuan College | Method for evaluating aerodynamic performance of wind wheel of wind turbine based on three-factor fitting integral method |
| CN113236491A (zh) * | 2021-05-27 | 2021-08-10 | 华北电力大学 | 一种风力发电数字孪生系统 |
| CN113742862A (zh) * | 2021-08-30 | 2021-12-03 | 重庆大学 | 一种基于重心平衡的双风轮风电机组齿轮传动系统设计方法 |
| CN113931793A (zh) * | 2021-11-12 | 2022-01-14 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组仿真的变桨环节等值方法及系统 |
| CN114065514A (zh) * | 2021-11-12 | 2022-02-18 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组仿真的bem数据映射方法及系统 |
| CN114169197A (zh) * | 2021-11-30 | 2022-03-11 | 中国电力科学研究院有限公司 | 一种风电机组在环仿真系统及方法 |
| WO2022134639A1 (zh) * | 2020-12-22 | 2022-06-30 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组协同控制方法 |
| CN115755639A (zh) * | 2022-10-10 | 2023-03-07 | 中车株洲电力机车研究所有限公司 | 一种风电机组数字孪生仿真方法及系统 |
| CN116292095A (zh) * | 2023-03-27 | 2023-06-23 | 华北电力大学 | 一种串列式双风轮风电机组多场景协调运行方法 |
-
2023
- 2023-07-25 CN CN202310917715.2A patent/CN117108445B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2020100825A4 (en) * | 2020-05-22 | 2020-07-02 | China University Of Mining And Technology, Yinchuan College | Method for evaluating aerodynamic performance of wind wheel of wind turbine based on three-factor fitting integral method |
| WO2022134639A1 (zh) * | 2020-12-22 | 2022-06-30 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组协同控制方法 |
| CN113236491A (zh) * | 2021-05-27 | 2021-08-10 | 华北电力大学 | 一种风力发电数字孪生系统 |
| CN113742862A (zh) * | 2021-08-30 | 2021-12-03 | 重庆大学 | 一种基于重心平衡的双风轮风电机组齿轮传动系统设计方法 |
| CN113931793A (zh) * | 2021-11-12 | 2022-01-14 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组仿真的变桨环节等值方法及系统 |
| CN114065514A (zh) * | 2021-11-12 | 2022-02-18 | 中国华能集团清洁能源技术研究院有限公司 | 一种串列式双风轮风电机组仿真的bem数据映射方法及系统 |
| CN114169197A (zh) * | 2021-11-30 | 2022-03-11 | 中国电力科学研究院有限公司 | 一种风电机组在环仿真系统及方法 |
| CN115755639A (zh) * | 2022-10-10 | 2023-03-07 | 中车株洲电力机车研究所有限公司 | 一种风电机组数字孪生仿真方法及系统 |
| CN116292095A (zh) * | 2023-03-27 | 2023-06-23 | 华北电力大学 | 一种串列式双风轮风电机组多场景协调运行方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117108445A (zh) | 2023-11-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN117108445B (zh) | 一种串列式双风轮风电机组数字孪生仿真方法 | |
| Pao et al. | Control co-design of 13 MW downwind two-bladed rotors to achieve 25% reduction in levelized cost of wind energy | |
| Bortolotti et al. | Combined preliminary–detailed design of wind turbines | |
| Pechlivanoglou | Passive and active flow control solutions for wind turbine blades | |
| Lee et al. | Multi‐flexible‐body dynamic analysis of horizontal axis wind turbines | |
| Barlas et al. | State of the art and prospectives of smart rotor control for wind turbines | |
| CN108694277A (zh) | 基于多平台联合仿真的风力发电机动力响应分析方法 | |
| Rodriguez et al. | Strongly-coupled aeroelastic free-vortex wake framework for floating offshore wind turbine rotors. Part 1: Numerical framework | |
| CN105626378A (zh) | 基于rbf神经网络的大型风电机组独立变桨控制方法 | |
| Selvam | Individual pitch control for large scale wind turbines | |
| Fasel et al. | Aeroservoelastic optimization of morphing airborne wind energy wings | |
| Dai et al. | Load and dynamic characteristic analysis of wind turbine flexible blades | |
| CN110714881B (zh) | 一种风力机性能预测控制方法及装置 | |
| CN116484671A (zh) | 风力机叶片全耦合动力响应精细化分析方法 | |
| Boulamatsis et al. | Active control of wind turbines through varying blade tip sweep | |
| Licitra et al. | Viability assessment of a rigid wing airborne wind energy pumping system | |
| CN111120222B (zh) | 一种具有风剪切和塔影效应的真实风况风力发电模拟装置 | |
| Montinari et al. | Ultimate and fatigue load mitigation by an inertial-driven passive flap, using a geometrically exact multibody formulation | |
| CN111639428B (zh) | 带襟翼的风力发电机的动力学仿真方法、设备及介质 | |
| Saito et al. | Application of an active device for helicopter noise reduction in JAXA | |
| Cherubini | Kite dynamics and wind energy harvesting | |
| Resor et al. | Impact of Higher Fidelity Models on Active Aerodynamic Load Control for Fatigue Damage Reduction | |
| Scott | Optimal Aeroelastic Tailoring of Wind Turbines | |
| Chen et al. | Impact of blade flexibility on wind turbine loads and pitch settings | |
| Wendell | Simplified aeroelastic modeling of horizontal axis wind turbines |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |