CN113806871B - 一种考虑结构非线性的柔性飞行动力学建模方法 - Google Patents

一种考虑结构非线性的柔性飞行动力学建模方法 Download PDF

Info

Publication number
CN113806871B
CN113806871B CN202111365361.2A CN202111365361A CN113806871B CN 113806871 B CN113806871 B CN 113806871B CN 202111365361 A CN202111365361 A CN 202111365361A CN 113806871 B CN113806871 B CN 113806871B
Authority
CN
China
Prior art keywords
generalized
aircraft
aerodynamic force
elastic
freedom
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
Application number
CN202111365361.2A
Other languages
English (en)
Other versions
CN113806871A (zh
Inventor
李道春
邵浩原
阚梓
申童
姚卓尔
赵仕伟
向锦武
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202111365361.2A priority Critical patent/CN113806871B/zh
Publication of CN113806871A publication Critical patent/CN113806871A/zh
Application granted granted Critical
Publication of CN113806871B publication Critical patent/CN113806871B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft 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/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
    • 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
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明公开了一种考虑结构非线性的柔性飞行动力学建模方法，利用非定常气动力建立包含非线性气动弹性机翼的飞行器飞行动力学方程；采用有理式拟合迟滞非线性广义力与位移之间的关系，以此建立了同时含有非线性刚度和非线性阻尼的飞行器飞行动力学方程；最终利用数值方法或谐波平衡法求解非线性飞行动力学模型的响应特征和稳定性。本发明建立的描述结构非线性的统一表达形式，可以解决控制面铰链以及连接件之间的安装间隙、部件老化松动等情况下飞行器间隙非线性动力学建模及后续动力学分析问题，从而为相关控制器设计提供指导思想，具有一定的工程应用价值。

Description

一种考虑结构非线性的柔性飞行动力学建模方法
技术领域
本发明涉及飞行器动力学技术领域，具体涉及一种考虑结构非线性的柔性飞行动力学建模方法，可以高效准确的解决高速飞行器和大展弦比无人机飞行动力学建模中的结构非线性问题及动力学分析，为飞行器后续的控制设计提供了模型依据。
背景技术
随着航空科学技术的不断进步，现代飞行器性能要求的不断提高，经典线性气动弹性理论已经不能完全满足飞行器设计的需要，非线性气动弹性问题逐渐引起了人们的重视。气动弹性问题中的非线性因素可以分为气动非线性和结构非线性两种。其中，结构非线性又可以分为分布非线性和集中非线性，分布非线性一般是由弹性变形引起，它影响整个结构；集中非线性作用于局部，在控制装置、连接件、吊舱等部位最为常见。
结构非线性因素可能会大大降低气动弹性系统的临界颤振速度，目前针对结构非线性问题常采用立方非线性模型、间隙非线性模型和迟滞非线性模型进行动力学建模，但是上述三种典型的结构非线性模型难以准确地描述气动弹性中的非线性关系，而且间隙和双线性迟滞模型的分段特性也会给研究带来不便。因此需要建立合理的考虑结构非线性的动力学模型，进而实现非线性因素对气动弹性系统稳定性、动态响应特征等的影响分析，以使飞行器在设计过程中尽可能地避免非线性因素带来的不利影响。
申请号为CN201911211566.8的中国发明专利申请公开了一种高精度亚音速固定翼飞行器非线性动力学建模方法，通过结构辨识和真实飞行科目中的飞行数据的参数辨识，建立飞行动力学模型，该方法是用于降低研发时间成本，避免使用CFD和风洞实验的飞行动力学建模方法，难以解决具量化的结构非线性因素评估和评价问题；申请号为CN202110247995.1的中国发明专利申请公开了一种基于非线性气动力降阶的弹性飞行动力学建模方法，通过非线性气动力降阶模型和刚弹耦合动力学方程对弹性飞机进行飞行动力学仿真。上述方法均对于结构非线性因素，尤其是舵面非线性因素的考虑具有一定的局限性，难以准确地描述气动弹性中的非线性关系。
发明内容
有鉴于此，本发明针对柔性飞行动力学建模中的结构非线性问题，提供了一种考虑结构非线性的柔性飞行动力学建模方法，包括以下步骤：
一种考虑结构非线性的柔性飞行动力学建模方法，包括以下步骤：
（1）确定飞行器运动学参数关系，建立飞行器运动学与动力学方程；
（2）求解机体刚体运动和弹性变形产生的广义气动力；
（3）建立舵面偏转产生的广义气动力模型，为飞行器动力学模型增加控制输入；
（4）对于已建立的飞行器动力学模型进行求解仿真，利用Runge-Kutta数值方法或谐波平衡法求解飞行器动力学模型的响应特征和稳定性。
特别地，所述步骤（1）具体为：
考虑未变形状态，与刚体飞行器固联的体轴系在惯性系的位置向量
Figure 298158DEST_PATH_IMAGE001
，用欧拉角表达从惯性系到体轴系的姿态变换
Figure 99892DEST_PATH_IMAGE002
，体轴系的速度和转动速度在体轴系的表达为
Figure 179843DEST_PATH_IMAGE003
，
Figure 533464DEST_PATH_IMAGE004
，描述飞行器质点运动和姿态运动的运动学方程为：
Figure 439978DEST_PATH_IMAGE005
Figure 904458DEST_PATH_IMAGE006
其中，
Figure 104495DEST_PATH_IMAGE007
、
Figure 238804DEST_PATH_IMAGE008
、
Figure 524292DEST_PATH_IMAGE009
分别为位置向量在惯性系三个轴上的投影；
Figure 792462DEST_PATH_IMAGE010
、
Figure 456793DEST_PATH_IMAGE011
、
Figure 152216DEST_PATH_IMAGE012
分别为滚转角、俯仰角和偏航角；
Figure 659421DEST_PATH_IMAGE013
、
Figure 606648DEST_PATH_IMAGE014
、
Figure 515699DEST_PATH_IMAGE015
分别为滚转角速度、俯仰角速度和偏航角速度；
Figure 44942DEST_PATH_IMAGE016
、
Figure 773863DEST_PATH_IMAGE017
、
Figure 383836DEST_PATH_IMAGE018
分别为飞行器速度在体轴系三个轴上的投影；
转换矩阵：
Figure 757180DEST_PATH_IMAGE019
，
Figure 59985DEST_PATH_IMAGE020
；
集中质量的转动速度
Figure 276203DEST_PATH_IMAGE021
由体轴系的转动速度
Figure 565233DEST_PATH_IMAGE022
和变形带来的转动速度
Figure 183296DEST_PATH_IMAGE023
共同组成：
Figure 532369DEST_PATH_IMAGE024
节点的弹性位移包括平移和转动，假设弹性变形的广义坐标满足正交条件的自由-自由边界条件下的模态表达，设弹性变形的广义坐标为
Figure 235883DEST_PATH_IMAGE025
，节点弹性位移表达为
Figure 453237DEST_PATH_IMAGE026
其中，
Figure 299708DEST_PATH_IMAGE027
为节点的运动变形，
Figure 678737DEST_PATH_IMAGE028
为节点的扭转变形，
Figure 135126DEST_PATH_IMAGE029
和
Figure 500380DEST_PATH_IMAGE030
分别为广义坐标
Figure 358614DEST_PATH_IMAGE031
中平动和转动自由度对应的广义坐标模态；
采用准坐标系下表达的拉格朗日方程，建立动力学方程：
Figure 642965DEST_PATH_IMAGE032
其中，拉格朗日算子
Figure 196437DEST_PATH_IMAGE033
，
Figure 755594DEST_PATH_IMAGE034
、
Figure 202756DEST_PATH_IMAGE035
分别为动能和势能，
Figure 64533DEST_PATH_IMAGE036
，
Figure 964356DEST_PATH_IMAGE037
，
Figure 796046DEST_PATH_IMAGE038
分别为平动自由度、转动自由度、弹性自由度对应的非保守力；
对于动能
Figure 2774DEST_PATH_IMAGE034
，由每个集中质量
Figure 894507DEST_PATH_IMAGE039
的平动动能
Figure 281626DEST_PATH_IMAGE040
和转动动能
Figure 57952DEST_PATH_IMAGE041
构成，有：
Figure 479706DEST_PATH_IMAGE042
其中，
Figure 683285DEST_PATH_IMAGE043
为集中质量的位置矢量，
Figure 557700DEST_PATH_IMAGE044
为集中质量的转动惯量；
对于势能
Figure 731193DEST_PATH_IMAGE035
，势能包含重力势能和弹性势能，利用线弹性假设，全机势能表达式为
Figure 148399DEST_PATH_IMAGE045
其中，
Figure 647513DEST_PATH_IMAGE046
，
Figure 9224DEST_PATH_IMAGE047
，
Figure 360309DEST_PATH_IMAGE048
为当地重力加速度，
Figure 491076DEST_PATH_IMAGE049
为刚度阵，
Figure 161092DEST_PATH_IMAGE050
为弹性变形；
全机质量表达为
Figure 885465DEST_PATH_IMAGE051
，对拉格朗日变量求微分，结合非保守力
Figure 931919DEST_PATH_IMAGE052
得到平动自由度方程为：
Figure 917192DEST_PATH_IMAGE053
拉格朗日变量对角度和角速度求导，得到转动自由度方程：
Figure 633475DEST_PATH_IMAGE054
拉格朗日变量对广义弹性坐标求导，得到飞行器弹性方程：
Figure 969779DEST_PATH_IMAGE055
对参考平衡状态的飞行器，动力学方程进行小扰动线性化处理，其中参考平衡状态的平动运动参数、转动运动参数分别为常数，弹性运动参数定义为0，保留一阶扰动项，整合动力学模型并进行简化，有：
Figure 429710DEST_PATH_IMAGE056
其中，
Figure 535069DEST_PATH_IMAGE057
为飞行器转动惯量，
Figure 281309DEST_PATH_IMAGE058
为飞行器动量矩，
Figure 744389DEST_PATH_IMAGE059
为广义质量矩阵，
Figure 867065DEST_PATH_IMAGE060
为广义阻尼矩阵，
Figure 826931DEST_PATH_IMAGE061
为广义刚度矩阵，非保守力项
Figure 885017DEST_PATH_IMAGE062
，其中
Figure 195913DEST_PATH_IMAGE063
表示气动力分量，
Figure 856701DEST_PATH_IMAGE064
表示推力分量；
对于非保守力项中的气动力分量进一步分为：
Figure 812019DEST_PATH_IMAGE065
其中
Figure 165640DEST_PATH_IMAGE066
为机体刚体运动和弹性变形产生的广义气动力，
Figure 698252DEST_PATH_IMAGE067
为舵面偏转产生的广义气动力。
特别地，所述步骤（2）具体为：
给定t时刻运动变形，利用非定常气动力求解算法或者CFD中动网格方法，求解气动力分布矩阵P；结合弹性模态阵，获得结构所有模态阵；用
Figure 303677DEST_PATH_IMAGE068
和
Figure 238135DEST_PATH_IMAGE069
表示刚体运动模态振型，
Figure 497078DEST_PATH_IMAGE070
表示弹性变形模态振型；结构所有模态阵表示为
Figure 156467DEST_PATH_IMAGE071
；通过插值函数将结构模态振型插值到气动网格上，得到对应气动分布矩阵自由度的广义气动模态振型
Figure 159058DEST_PATH_IMAGE072
；则t时刻广义气动力表达为
Figure 948023DEST_PATH_IMAGE073
。
特别地，所述步骤（2）具体为：
给定t时刻运动变形，用
Figure 784392DEST_PATH_IMAGE074
和
Figure 557176DEST_PATH_IMAGE075
表示刚体运动模态和弹性变形模态振型表达，结构所有模态阵表示为
Figure 973245DEST_PATH_IMAGE076
；对t时刻刚体运动及弹性变形具体值参数导致的广义气动力
Figure 882295DEST_PATH_IMAGE077
、
Figure 748620DEST_PATH_IMAGE078
求解，结构所有模态阵
Figure 618487DEST_PATH_IMAGE079
与气动力模型的模态对应，将模态阵代入给定非定常气动力模型，则t时刻广义气动力表达为：
Figure 962880DEST_PATH_IMAGE080
。
特别地，所述步骤（3）具体为：
考虑双线性迟滞模型，将位移上升段和位移下降段分开描述；舵面广义气动力
Figure 992016DEST_PATH_IMAGE081
，其中
Figure 403144DEST_PATH_IMAGE082
为舵面偏转广义自由度矩阵；
针对每一个舵面偏转广义自由度
Figure 884941DEST_PATH_IMAGE083
，通过实验测得双线性迟滞非线性广义力与位移之间的关系，分别用上行有理式
Figure 33025DEST_PATH_IMAGE084
和下行有理式
Figure 792034DEST_PATH_IMAGE085
拟合，表达为：
Figure 161DEST_PATH_IMAGE086
采用非线性有理式进行函数曲线拟合，以统一的表达形式描述模型中的结构非线性，得到广义位移与广义力的表达式。
特别地，所述步骤（4）中，利用谐波平衡法进行模型气动弹性求解分析，推导出以颤振频率为未知数的高次代数方程，由高次代数方程解出颤振频率，求出两个方向的颤振振幅。
本发明相对于现有技术的有益效果：
（1）本发明的建模周期短，仿真计算效率高，可将机体运动和控制面运动进行单独仿真计算并引入飞行器飞行动力学仿真中。
（2）本发明通过建立飞行器动力学模型中结构非线性因素的统一表达形式，有效解决了飞行器动力学建模中的结构非线性建模问题，进而提高模型动力学分析中的合理性及准确性。
（3）本发明考虑柔性飞行动力学建模过程中的结构非线性因素，在完整的弹性飞机气动-结构耦合分析基础上，可以解决控制面铰链以及连接件之间的安装间隙、部件老化松动等情况下飞行器结构非线性动力学建模及后续动力学分析问题，从而为相关控制器设计提供指导思想，具有一定的工程应用价值。
附图说明
图1为动力学建模方法的流程示意图；
图2为实施例极限环震荡时间历程；
图3为含有间隙非线性机翼飞行器俯仰位移的时间响应历程；
图4为含有间隙非线性机翼飞行器无量纲沉浮位移的时间响应历程；
图5为飞行器无量纲临界颤振速度与初始角位移关系。
具体实施方式
下面结合附图和实施例对本发明的一种考虑结构非线性的柔性飞行动力学建模方法进一步详细说明。
一种考虑结构非线性的柔性飞行动力学建模方法，包括以下步骤：
(1)确定飞行器运动学参数关系，建立飞行器运动学与动力学方程。考虑未变形状态，与刚体飞行器固联的体轴系在惯性系的位置向量
Figure 844621DEST_PATH_IMAGE087
，用欧拉角表达从惯性系到体轴系的姿态变换
Figure 530817DEST_PATH_IMAGE088
，飞行器速度和转动速度在体轴系的表达为
Figure 534545DEST_PATH_IMAGE089
，
Figure 788940DEST_PATH_IMAGE090
，描述飞行器质点运动和姿态运动的运动学方程为：
Figure 979750DEST_PATH_IMAGE091
Figure 735216DEST_PATH_IMAGE092
其中，
Figure 436194DEST_PATH_IMAGE093
、
Figure 251703DEST_PATH_IMAGE094
、
Figure 929809DEST_PATH_IMAGE095
分别为滚转角、俯仰角和偏航角；
Figure 98753DEST_PATH_IMAGE096
、
Figure 811494DEST_PATH_IMAGE097
、
Figure 407692DEST_PATH_IMAGE098
分别为滚转角速度、俯仰角速度和偏航角速度；
Figure 573094DEST_PATH_IMAGE099
、
Figure 404784DEST_PATH_IMAGE100
、
Figure 112977DEST_PATH_IMAGE101
分别为飞行器速度在体轴系三个轴上的投影。
转换矩阵：
Figure 739130DEST_PATH_IMAGE019
，
Figure 657408DEST_PATH_IMAGE020
；
集中质量的转动速度
Figure 666690DEST_PATH_IMAGE102
由体轴系的转动速度
Figure 88444DEST_PATH_IMAGE103
和变形带来的转动速度
Figure 151078DEST_PATH_IMAGE104
共同组成：
Figure 900859DEST_PATH_IMAGE105
节点的弹性位移可以表达为平移和转动，根据假设条件，可以满足正交条件的自由-自由边界条件下的模态表达，设弹性变形的广义坐标为
Figure 605510DEST_PATH_IMAGE025
，节点弹性位移可以表达为
Figure 881770DEST_PATH_IMAGE026
其中，
Figure 256251DEST_PATH_IMAGE027
为节点的运动变形，
Figure 617962DEST_PATH_IMAGE028
为节点的扭转变形，
Figure 736091DEST_PATH_IMAGE029
和
Figure 866858DEST_PATH_IMAGE030
为广义坐标
Figure 271294DEST_PATH_IMAGE031
对应的广义坐标模态。
采用准坐标系下表达的拉格朗日方程，建立动力学方程：
Figure 759782DEST_PATH_IMAGE106
其中拉格朗日算子
Figure 275077DEST_PATH_IMAGE033
，
Figure 791509DEST_PATH_IMAGE034
、
Figure 242213DEST_PATH_IMAGE035
分别为动能和势能，
Figure 844096DEST_PATH_IMAGE036
，
Figure 163082DEST_PATH_IMAGE037
，
Figure 409387DEST_PATH_IMAGE038
为非保守力。
对于动能项，由每个集中质量
Figure 155626DEST_PATH_IMAGE039
的平动动能
Figure 713646DEST_PATH_IMAGE040
和转动动能
Figure 977268DEST_PATH_IMAGE041
构成，有：
Figure 937134DEST_PATH_IMAGE042
其中，
Figure 854274DEST_PATH_IMAGE043
为集中质量的位置矢量，
Figure 804651DEST_PATH_IMAGE044
为集中质量的转动惯量。
对于势能项，势能包含重力势能和弹性势能，利用线弹性假设，全机势能表达式为
Figure 465439DEST_PATH_IMAGE045
其中，
Figure 420757DEST_PATH_IMAGE046
，
Figure 774378DEST_PATH_IMAGE047
，
Figure 306990DEST_PATH_IMAGE048
为当地重力加速度，
Figure 912415DEST_PATH_IMAGE049
为刚度阵，
Figure 846873DEST_PATH_IMAGE050
为弹性变形。
利用上述表达式，全机质量表达为
Figure 371395DEST_PATH_IMAGE051
，对拉格朗日变量求微分，结合广义力
Figure 266670DEST_PATH_IMAGE052
可以得到平动自由度方程为：
Figure 534840DEST_PATH_IMAGE053
拉格朗日变量对角度和角速度求导，得到转动自由度方程：
Figure 323805DEST_PATH_IMAGE054
其中，
Figure 658709DEST_PATH_IMAGE058
为飞行器动量。
拉格朗日变量对广义弹性坐标求导，得到飞行器弹性方程：
Figure 900334DEST_PATH_IMAGE055
对参考平衡状态的飞行器，动力学方程进行小扰动线性化处理，其中参考平衡状态的平动运动参数、转动运动参数分别为常数，弹性运动参数定义为0，保留一阶扰动项。整合动力学模型并进行简化，有：
Figure 706616DEST_PATH_IMAGE056
其中，
Figure 756612DEST_PATH_IMAGE059
为广义质量矩阵，
Figure 357358DEST_PATH_IMAGE060
为广义阻尼矩阵，
Figure 492804DEST_PATH_IMAGE061
为广义刚度矩阵，非保守力项
Figure 837198DEST_PATH_IMAGE107
，其中
Figure 600754DEST_PATH_IMAGE063
表示气动力分量，
Figure 513347DEST_PATH_IMAGE064
表示推力分量。
对于非保守力项中的气动力分量可以进一步分为：
Figure 995143DEST_PATH_IMAGE065
其中
Figure 143228DEST_PATH_IMAGE066
为机体刚体运动和弹性变形产生的广义气动力，
Figure 400772DEST_PATH_IMAGE067
为舵面偏转产生的广义气动力。
(2) 求解机体刚体运动和弹性变形产生的广义气动力，求解基于非定常广义气动力的过程可以描述为：
给定t时刻运动变形，利用非定常气动力求解算法或者CFD中动网格方法，求解气动力分布矩阵P。结合Nastran等结构计算软件导出的弹性模态阵，获得结构所有模态阵。用
Figure 608899DEST_PATH_IMAGE068
和
Figure 577992DEST_PATH_IMAGE108
表示刚体运动模态振型，
Figure 139555DEST_PATH_IMAGE109
表示弹性变形模态振型。结构所有模态阵可以表示为
Figure 877704DEST_PATH_IMAGE110
。通过插值函数将结构模态振型插值到气动网格上，得到对应气动分布矩阵自由度的广义气动模态振型
Figure 256732DEST_PATH_IMAGE072
。进而可以把t时刻模型输出的广义气动力表达为
Figure 322908DEST_PATH_IMAGE073
。
特别的，为简便计算采用Theodorsen等简单非定常气动力计算方法时，结构所有模态阵往往与气动力分布矩阵P对应，因而可以不用通过插值函数将结构模态振型插值到气动网格，即
Figure 78375DEST_PATH_IMAGE111
。
(3) 建立舵面偏转产生的广义气动力模型，为飞行器动力学模型增加控制输入。本发明中考虑双线性迟滞模型，将位移上升段和位移下降段分开描述。舵面广义气动力
Figure 671030DEST_PATH_IMAGE112
，其中
Figure 361906DEST_PATH_IMAGE113
为舵面偏转广义自由度矩阵。
针对每一个舵面偏转广义自由度
Figure 40012DEST_PATH_IMAGE083
，可以通过实验测得双线性迟滞非线性广义力与位移之间的关系，分别用上行有理式
Figure 68011DEST_PATH_IMAGE084
和下行有理式
Figure 154653DEST_PATH_IMAGE085
拟合，表达为：
Figure 875484DEST_PATH_IMAGE086
采用非线性有理式进行函数曲线拟合，以统一的表达形式来描述模型中的结构非线性，可以给出具体的广义位移与广义力的表达式。
(4)对于已建立的飞行器动力学仿真模型，利用Runge-Kutta数值方法对其进行求解。基于飞行器动力学模型可对时间步内弹性飞行稳定性进行分析，在此基础上，加入气动力模型以对气动力响应进行分析，再加入舵面气动力模型以进行气动伺服弹性及飞行仿真分析。
特别的，利用谐波平衡法可以进行模型气动弹性求解分析，首先推导出以颤振频率为未知数的高次代数方程，由此方程解出颤振频率，然后即可求出两个方向的颤振振幅。该方法与数值积分方法相比，可以避免积分过程的累积误差，计算量较小，但是推导过程较复杂，而且不能考察初始条件对系统响应的影响。
具体实施例一
如图1所示，本实施例针对只考虑纵向平面的飞行器模型（图2），对其进行考虑结构非线性的柔性飞行动力学建模。
第一步，对飞行器进行坐标系定义，体轴系原点位于中心处，初始位置与惯性系重合，确定飞行器简化模型自由度分别为：沉浮自由度
Figure 40886DEST_PATH_IMAGE114
，向下为正；俯仰自由度
Figure 13522DEST_PATH_IMAGE115
，飞行器抬头为正；控制面偏转自由度
Figure 315190DEST_PATH_IMAGE116
，向下为正。飞行器弹性轴到重心的距离为
Figure 347868DEST_PATH_IMAGE117
，弦向中点到弹性轴的距离为
Figure 566DEST_PATH_IMAGE118
；操纵面铰链到机翼中点的距离为
Figure 370368DEST_PATH_IMAGE119
，操纵面重心到操纵面铰链的距离为
Figure 933067DEST_PATH_IMAGE120
，其中
Figure 261280DEST_PATH_IMAGE121
为机翼半弦长，
Figure 135696DEST_PATH_IMAGE122
为机翼关于弹性轴的回转半径。
Figure 214248DEST_PATH_IMAGE123
其中，
Figure 490508DEST_PATH_IMAGE124
，
Figure 724043DEST_PATH_IMAGE125
为广义气动力矩阵，质量矩阵
Figure 226700DEST_PATH_IMAGE126
、刚度矩阵
Figure 938304DEST_PATH_IMAGE127
分别为
Figure 334650DEST_PATH_IMAGE128
其中，m t为飞行器的总质量，S α 、S β 分别为俯仰方向和副翼偏转的静矩，I α 、I β 分别为相应的惯性矩。k α 为俯仰自由度刚度，k h 为沉浮自由度刚度，k β 为线性操纵面刚度。将方程无量纲化，并引入阻尼，则质量矩阵、阻尼矩阵和刚度矩阵分别为：
Figure 880032DEST_PATH_IMAGE129
其中，r α 为I α 的无量纲量，r β 为I β 的无量纲量，x α 为S α 的无量纲量，x β 为S β 的无量纲量，m为参考质量一般取为1个单位展长的机翼的质量，
Figure 729040DEST_PATH_IMAGE130
为俯仰自由度的自然频率比，
Figure 385280DEST_PATH_IMAGE131
为舵面偏转自由度的自然频率比，
Figure 636133DEST_PATH_IMAGE132
为沉浮自由度的自然频率比，M α 和M β 分别为与系统中俯仰自由度与舵面偏转相关的非线性因素组成的非线性函数，
Figure 211471DEST_PATH_IMAGE133
为方程
Figure 187255DEST_PATH_IMAGE134
的特征根
Figure 506241DEST_PATH_IMAGE135
的特征向量矩阵，
Figure 877179DEST_PATH_IMAGE136
为
Figure 498784DEST_PATH_IMAGE137
的线性化形式。
Figure 322384DEST_PATH_IMAGE138
为自然耦合频率，
Figure 445061DEST_PATH_IMAGE139
为矩阵
Figure 280293DEST_PATH_IMAGE140
的对角元素，
Figure 463012DEST_PATH_IMAGE141
为测定阻尼比。
第二步，利用Theodorsen理论求解气动力，无量纲形式的非定常气动力与力矩向量
Figure 773908DEST_PATH_IMAGE142
，其中各分量具体表达式为：
Figure 310063DEST_PATH_IMAGE143
其中，U为无量纲来流速度，
Figure 390014DEST_PATH_IMAGE144
为空气密度，
Figure 743635DEST_PATH_IMAGE145
（i=1,2,3,…,14）为Theodorsen常数，具体为
Figure 939166DEST_PATH_IMAGE146
建模得到完整的以矩阵形式表达的三自由度飞行器运动方程为
Figure 138066DEST_PATH_IMAGE147
其中，
Figure 72524DEST_PATH_IMAGE148
为新定义变量，矩阵
Figure 737991DEST_PATH_IMAGE149
、
Figure 492321DEST_PATH_IMAGE150
、
Figure 901437DEST_PATH_IMAGE151
分别为：
Figure 690401DEST_PATH_IMAGE152
其它向量具体为：
Figure 120245DEST_PATH_IMAGE153
其中，
Figure 768395DEST_PATH_IMAGE154
其中，
Figure 574677DEST_PATH_IMAGE155
为机翼面积，
Figure 483728DEST_PATH_IMAGE156
为空气密度，
Figure 723954DEST_PATH_IMAGE157
和
Figure 718455DEST_PATH_IMAGE158
分别为单位功角对应的升力系数和力矩系数，
Figure 328428DEST_PATH_IMAGE159
为机翼质量,
Figure 232930DEST_PATH_IMAGE160
为机翼和支架的总质量；
Figure 4577DEST_PATH_IMAGE161
为机翼惯性矩；
Figure 627319DEST_PATH_IMAGE162
为沉浮方向的结构阻尼系数与刚度。
第三步，建立舵面偏转产生的广义气动力模型，增加控制输入
Figure 775403DEST_PATH_IMAGE163
，定义向量
Figure 393467DEST_PATH_IMAGE164
，
Figure 476960DEST_PATH_IMAGE165
，有
Figure 180474DEST_PATH_IMAGE166
其中，
Figure 132250DEST_PATH_IMAGE167
矩阵中各元素定义为：
Figure 244300DEST_PATH_IMAGE168
所用到的参数为：
Figure 888908DEST_PATH_IMAGE169
，初始条件为
Figure 79718DEST_PATH_IMAGE170
。
加入迟滞非线性模型，其中
Figure 444971DEST_PATH_IMAGE171
，通过曲线拟合，有理式近似可以分为上行和下行两个有理式，分别为：
Figure 303206DEST_PATH_IMAGE172
实际迟滞非线性与拟合结果如图2所示。
第四步，利用数值方法求解飞行器动力学模型。在无量纲来流速度U=1时，俯仰方向含有迟滞非线性飞行器俯仰和沉浮位移的时间响应历程分别如图3和图4所示。进行飞行器颤振分析，以初始角位移为例，在
Figure 853136DEST_PATH_IMAGE173
的范围内，以0.01弧度为间隔逐渐改变初始角位移的大小，保持其他初始条件不变，可以做出系统临界颤振速度随初始角位移变化的颤振边界，如图5所示。
以上所述仅为本发明的具体实施方式，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种考虑结构非线性的柔性飞行动力学建模方法，其特征在于，包括以下步骤：
（1）确定飞行器运动学参数关系，建立飞行器运动学与动力学方程；
（2）求解机体刚体运动和弹性变形产生的广义气动力；
（3）建立舵面偏转产生的广义气动力模型，为飞行器动力学模型增加控制输入；
（4）对于已建立的飞行器动力学模型进行求解仿真，利用Runge-Kutta数值方法或谐波平衡法求解飞行器动力学模型的响应特征和稳定性；
所述步骤（1）具体为：
考虑未变形状态，与刚体飞行器固联的体轴系在惯性系的位置向量
Figure 340850DEST_PATH_IMAGE001
，用欧拉角表达从惯性系到体轴系的姿态变换
Figure 251037DEST_PATH_IMAGE002
，体轴系的速度和转动速度在体轴系的表达为
Figure 989317DEST_PATH_IMAGE003
，
Figure 541521DEST_PATH_IMAGE004
，描述飞行器质点运动和姿态运动的运动学方程为：
Figure 929251DEST_PATH_IMAGE005
Figure 275919DEST_PATH_IMAGE006
其中，
Figure 767074DEST_PATH_IMAGE007
、
Figure 857390DEST_PATH_IMAGE008
、
Figure 350557DEST_PATH_IMAGE009
分别为位置向量在惯性系三个轴上的投影；
Figure 602547DEST_PATH_IMAGE010
、
Figure 315419DEST_PATH_IMAGE011
、
Figure 209425DEST_PATH_IMAGE012
分别为滚转角、俯仰角和偏航角；
Figure 825608DEST_PATH_IMAGE013
、
Figure 248499DEST_PATH_IMAGE014
、
Figure 448667DEST_PATH_IMAGE015
分别为飞行器速度在体轴系三个轴上的投影；
Figure 880786DEST_PATH_IMAGE016
、
Figure 974644DEST_PATH_IMAGE017
、
Figure 83283DEST_PATH_IMAGE018
分别为滚转角速度、俯仰角速度和偏航角速度；
转换矩阵:
Figure 754436DEST_PATH_IMAGE019
，
Figure 6557DEST_PATH_IMAGE020
；
集中质量的转动速度
Figure 813976DEST_PATH_IMAGE021
由体轴系的转动速度
Figure 362025DEST_PATH_IMAGE022
和变形带来的转动速度
Figure 520474DEST_PATH_IMAGE023
共同组成：
Figure 45127DEST_PATH_IMAGE024
节点的弹性位移包括平移和转动，假设弹性变形的广义坐标满足正交条件的自由-自由边界条件下的模态表达，设弹性变形的广义坐标为
Figure 972632DEST_PATH_IMAGE025
，节点弹性位移表达为
Figure 688653DEST_PATH_IMAGE026
其中，
Figure 334398DEST_PATH_IMAGE027
为节点的运动变形，
Figure 397163DEST_PATH_IMAGE028
为节点的扭转变形，
Figure 710333DEST_PATH_IMAGE029
和
Figure 334605DEST_PATH_IMAGE030
分别为广义坐标
Figure 467646DEST_PATH_IMAGE031
中平动和转动自由度对应的广义坐标模态；
采用准坐标系下表达的拉格朗日方程，建立动力学方程：
Figure 724315DEST_PATH_IMAGE032
其中，拉格朗日算子
Figure 377145DEST_PATH_IMAGE033
，
Figure 920122DEST_PATH_IMAGE034
、
Figure 789726DEST_PATH_IMAGE035
分别为动能和势能，
Figure 709141DEST_PATH_IMAGE036
，
Figure 216477DEST_PATH_IMAGE037
，
Figure 195934DEST_PATH_IMAGE038
分别为平动自由度、转动自由度、弹性自由度对应的非保守力；
对于动能
Figure 544046DEST_PATH_IMAGE034
，由每个集中质量
Figure 1572DEST_PATH_IMAGE039
的平动动能
Figure 628993DEST_PATH_IMAGE040
和转动动能
Figure 779352DEST_PATH_IMAGE041
构成，有：
Figure 357970DEST_PATH_IMAGE042
其中，
Figure 619187DEST_PATH_IMAGE043
为集中质量的位置矢量，
Figure 225749DEST_PATH_IMAGE044
为集中质量的转动惯量；
对于势能
Figure 297741DEST_PATH_IMAGE035
，势能包含重力势能和弹性势能，利用线弹性假设，全机势能表达式为
Figure 114387DEST_PATH_IMAGE045
其中，
Figure 431493DEST_PATH_IMAGE046
，
Figure 17195DEST_PATH_IMAGE047
，
Figure 260088DEST_PATH_IMAGE048
为当地重力加速度，
Figure 564031DEST_PATH_IMAGE049
为刚度阵，
Figure 416318DEST_PATH_IMAGE050
为弹性变形；
全机质量表达为
Figure 122106DEST_PATH_IMAGE051
，对拉格朗日变量求微分，结合非保守力
Figure 535901DEST_PATH_IMAGE052
得到平动自由度方程为：
Figure 327139DEST_PATH_IMAGE053
拉格朗日变量对角度和角速度求导，得到转动自由度方程：
Figure 720468DEST_PATH_IMAGE054
拉格朗日变量对广义弹性坐标求导，得到飞行器弹性方程：
Figure 421708DEST_PATH_IMAGE055
对参考平衡状态的飞行器，动力学方程进行小扰动线性化处理，其中参考平衡状态的平动运动参数、转动运动参数分别为常数，弹性运动参数定义为0，保留一阶扰动项，整合动力学模型并进行简化，有：
Figure 990092DEST_PATH_IMAGE056
其中，
Figure 19359DEST_PATH_IMAGE057
为飞行器转动惯量，
Figure 229761DEST_PATH_IMAGE058
为飞行器动量矩，
Figure 893829DEST_PATH_IMAGE059
为广义质量矩阵，
Figure 898694DEST_PATH_IMAGE060
为广义阻尼矩阵，
Figure 149678DEST_PATH_IMAGE061
为广义刚度矩阵，非保守力项
Figure 632612DEST_PATH_IMAGE062
，其中
Figure 685275DEST_PATH_IMAGE063
表示气动力分量，
Figure 595462DEST_PATH_IMAGE064
表示推力分量；
对于非保守力项中的气动力分量进一步分为：
Figure 333742DEST_PATH_IMAGE065
其中
Figure 885946DEST_PATH_IMAGE066
为机体刚体运动和弹性变形产生的广义气动力，
Figure 524607DEST_PATH_IMAGE067
为舵面偏转产生的广义气动力。
2.根据权利要求1所述的考虑结构非线性的柔性飞行动力学建模方法，其特征在于，所述步骤（2）具体为：
给定t时刻运动变形，利用非定常气动力求解算法或者CFD中动网格方法，求解气动力分布矩阵P；结合弹性模态阵，获得结构所有模态阵；用
Figure 871275DEST_PATH_IMAGE068
和
Figure 96851DEST_PATH_IMAGE069
表示刚体运动模态振型，
Figure 62533DEST_PATH_IMAGE070
表示弹性变形模态振型；结构所有模态阵表示为
Figure 306432DEST_PATH_IMAGE071
；通过插值函数将结构模态振型插值到气动网格上，得到对应气动分布矩阵自由度的广义气动模态振型
Figure 87917DEST_PATH_IMAGE072
；则t时刻广义气动力表达为
Figure 315636DEST_PATH_IMAGE073
。
3.根据权利要求1所述的考虑结构非线性的柔性飞行动力学建模方法，其特征在于，所述步骤（2）具体为：
给定t时刻运动变形，用
Figure 694796DEST_PATH_IMAGE074
和
Figure 793202DEST_PATH_IMAGE075
表示刚体运动模态和弹性变形模态振型表达，结构所有模态阵表示为
Figure 996519DEST_PATH_IMAGE076
；对t时刻刚体运动及弹性变形具体值参数导致的广义气动力
Figure 445955DEST_PATH_IMAGE077
、
Figure 628806DEST_PATH_IMAGE078
求解，结构所有模态阵
Figure 847298DEST_PATH_IMAGE079
与气动力模型的模态对应，将模态阵代入给定非定常气动力模型，则t时刻广义气动力表达为：
Figure 958866DEST_PATH_IMAGE080
。
4.根据权利要求1所述的考虑结构非线性的柔性飞行动力学建模方法，其特征在于，所述步骤（3）具体为：
考虑双线性迟滞模型，将位移上升段和位移下降段分开描述；舵面广义气动力
Figure 630019DEST_PATH_IMAGE081
，其中
Figure 6774DEST_PATH_IMAGE082
为舵面偏转广义自由度矩阵；
针对每一个舵面偏转广义自由度
Figure 830505DEST_PATH_IMAGE083
，通过实验测得双线性迟滞非线性广义力与位移之间的关系，分别用上行有理式
Figure 860777DEST_PATH_IMAGE084
和下行有理式
Figure 534073DEST_PATH_IMAGE085
拟合，表达为：
Figure 42415DEST_PATH_IMAGE086
采用非线性有理式进行函数曲线拟合，以统一的表达形式描述模型中的结构非线性，得到广义位移与广义力的表达式。
5.根据权利要求1所述的考虑结构非线性的柔性飞行动力学建模方法，其特征在于，所述步骤（4）中，利用谐波平衡法进行模型气动弹性求解分析，推导出以颤振频率为未知数的高次代数方程，由高次代数方程解出颤振频率，求出两个方向的颤振振幅。
CN202111365361.2A 2021-11-18 2021-11-18 一种考虑结构非线性的柔性飞行动力学建模方法 Active CN113806871B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111365361.2A CN113806871B (zh) 2021-11-18 2021-11-18 一种考虑结构非线性的柔性飞行动力学建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111365361.2A CN113806871B (zh) 2021-11-18 2021-11-18 一种考虑结构非线性的柔性飞行动力学建模方法

Publications (2)

Publication Number Publication Date
CN113806871A CN113806871A (zh) 2021-12-17
CN113806871B true CN113806871B (zh) 2022-02-08

Family

ID=78938314

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111365361.2A Active CN113806871B (zh) 2021-11-18 2021-11-18 一种考虑结构非线性的柔性飞行动力学建模方法

Country Status (1)

Country Link
CN (1) CN113806871B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114707370A (zh) * 2022-01-28 2022-07-05 北京航空航天大学 一种适用于弹性飞机的飞行仿真方法
CN114476123B (zh) * 2022-03-21 2023-04-14 中国空气动力研究与发展中心高速空气动力研究所 一种用于软管式空中加油对接模拟试验装置的设计方法
CN114942595B (zh) * 2022-07-25 2022-11-18 西安爱生技术集团有限公司 一种考虑降雨影响的无人机阵风响应建模和分析方法
CN115390512B (zh) * 2022-09-08 2023-05-30 上海交通大学 柔性龙门双驱系统及其机电联合解耦运动控制方法
CN115657729B (zh) * 2022-12-27 2023-03-10 北京航空航天大学 一种考虑探测任务约束的无人机边界保护控制方法
CN117077296B (zh) * 2023-10-17 2024-01-09 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器气动运动结构控制耦合仿真方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260492A (zh) * 2015-04-27 2016-01-20 中国直升机设计研究所 一种旋翼与机体耦合动力学模态综合建模方法
CN105843076A (zh) * 2016-03-31 2016-08-10 北京理工大学 一种柔性飞行器气动弹性建模与控制方法
CN112634704A (zh) * 2020-12-23 2021-04-09 上海科梁信息工程股份有限公司 一种飞行仿真系统、方法及存储介质

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6764047B2 (en) * 2001-01-26 2004-07-20 Todd Scott Miller Model airplane hinge construction
CN103593524B (zh) * 2013-11-13 2014-10-01 北京航空航天大学 一种变体飞行器的动力学建模与分析方法
US20150203215A1 (en) * 2014-01-17 2015-07-23 Eric T. Falangas Early performance evaluation of conceptual flight and space vehicles
CN105676672B (zh) * 2016-01-06 2017-05-31 西北工业大学 扑翼飞行器复合飞行策略仿真建模方法
CN106484969B (zh) * 2016-09-23 2019-08-09 中国运载火箭技术研究院 一种大包线强机动飞行器动力学高精度仿真方法
CN108052787B (zh) * 2018-02-01 2020-06-23 南京航空航天大学 基于飞行动态的高超声速飞行器机翼颤振损伤估计方法
CN108427322B (zh) * 2018-03-16 2020-04-17 南京航空航天大学 一种大柔性飞行器基于在线辨识的建模方法
CN110929339B (zh) * 2019-11-29 2024-05-10 海丰通航科技有限公司 高精度亚音速固定翼飞行器非线性动力学建模方法
CN112580240B (zh) * 2020-12-15 2022-08-05 北京航空航天大学 一种适用于复杂大柔性飞机建模的非线性子结构方法
CN113111430B (zh) * 2021-03-06 2022-09-30 北京航空航天大学 基于非线性气动力降阶的弹性飞机飞行动力学建模方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260492A (zh) * 2015-04-27 2016-01-20 中国直升机设计研究所 一种旋翼与机体耦合动力学模态综合建模方法
CN105843076A (zh) * 2016-03-31 2016-08-10 北京理工大学 一种柔性飞行器气动弹性建模与控制方法
CN112634704A (zh) * 2020-12-23 2021-04-09 上海科梁信息工程股份有限公司 一种飞行仿真系统、方法及存储介质

Also Published As

Publication number Publication date
CN113806871A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
CN113806871B (zh) 一种考虑结构非线性的柔性飞行动力学建模方法
CN113111430B (zh) 基于非线性气动力降阶的弹性飞机飞行动力学建模方法
CN108427322B (zh) 一种大柔性飞行器基于在线辨识的建模方法
Patil et al. On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ratio wings
CN113868771B (zh) 一种考虑结构和气动非线性的飞行动力学建模方法
Frulla et al. Critical behaviour of slender wing configurations
CN115659523B (zh) 一种大展弦比无人机刚柔耦合建模分析方法
Sun et al. Time-domain modeling and control of a wing-section stall flutter
Castellani et al. Flight loads prediction of high aspect ratio wing aircraft using multibody dynamics
Takarics et al. Flight control oriented bottom-up nonlinear modeling of aeroelastic vehicles
Xie et al. Geometrically nonlinear aeroelastic stability analysis and wind tunnel test validation of a very flexible wing
Suleman et al. Non-linear aeroelastic analysis in the time domain of high-aspect-ratio wings: Effect of chord and taper-ratio variation
Shafaghat et al. Nonlinear aeroelastic analysis of a HALE aircraft with flexible components
Koreanschi et al. Flutter analysis of a morphing wing technology demonstrator: numerical simulation and wind tunnel testing
Maza et al. A cost-effective algorithm for the co-simulation of unsteady and non-linear aeroelastic phenomena
Bath et al. Reduced order aerodynamics and flight mechanics coupling for tiltrotor aircraft stability analysis
Ricci et al. Active control of three-surface aeroelastic model
Castrichini Parametric Assessment of a Folding Wing-Tip Device for Aircraft Loads Alleviation
Kier et al. Development of aircraft flight loads analysis models with uncertainties for pre-design studies
Rizzi et al. Virtual aircraft design and control of transcruiser-A canard configuration
Wang et al. Static aeroelastic analysis of flexible aircraft with large deformations
Pontillo et al. Low-order Aeroelastic Modelling of a High Aspect Ratio Wing Aircraft Under Constrained Motion
Altememe et al. Preliminary design and computational fluid dynamic analysis of flapping wing of micro aerial vehicle for low Reynolds numbers regime
Kotikalpudi et al. Sensitivity of Robust Flutter Boundary to Model Uncertainties in Aeroservoelastic Systems-Invited
Seber et al. Nonlinear flutter calculations using finite elements in a direct Eulerian-Lagrangian formulation

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