CN114896830A - 导弹非线性非定常气动力微分方程模型辨识方法 - Google Patents

导弹非线性非定常气动力微分方程模型辨识方法 Download PDF

Info

Publication number
CN114896830A
CN114896830A CN202210825423.1A CN202210825423A CN114896830A CN 114896830 A CN114896830 A CN 114896830A CN 202210825423 A CN202210825423 A CN 202210825423A CN 114896830 A CN114896830 A CN 114896830A
Authority
CN
China
Prior art keywords
model
aerodynamic
differential equation
aerodynamic force
missile
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.)
Granted
Application number
CN202210825423.1A
Other languages
English (en)
Other versions
CN114896830B (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202210825423.1A priority Critical patent/CN114896830B/zh
Publication of CN114896830A publication Critical patent/CN114896830A/zh
Application granted granted Critical
Publication of CN114896830B publication Critical patent/CN114896830B/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/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • 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)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明提供了一种导弹非线性非定常气动力微分方程模型辨识方法，包括：S1：数据准备：利用风洞试验或CFD计算得到导弹静态气动力和力矩系数、大振幅俯仰振荡气动力和力矩系数时间历程的动态数据表，并经过数据处理后生成气动建模的输入数据文件；S2：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，采用一阶微分方程描述非定常增量，构建气动力微分方程模型；S3：将气动力微分方程模型辨识问题转化为动态系统的参数辨识问题；S4：利用所述气动数据，基于最小二乘准则，采用Gauss‑Newton优化算法辨识获取模型参数的估计值。本发明适用于
Figure 100004_DEST_PATH_IMAGE001
全攻角范围，模型泛化性能强。

Description

导弹非线性非定常气动力微分方程模型辨识方法
技术领域
本发明涉及空气动力建模技术领域，具体涉及一种导弹非线性非定常气动力微分方程模型辨识方法。
背景技术
大攻角机动能力是现代战术导弹的一项重要设计指标。特别是对于空射导弹，利用大攻角机动实现越肩发射，可以极大地缩短作战响应时间，提高载机生存概率和导弹作战效能。不论是传统的越肩发射还是新型自翻转越肩发射，导弹都会经历大攻角过失速状态，气动力呈现高非定常、强非线性特性。此时，基于准定常假设的气动导数模型已不适用，需要建立非线性非定常气动力模型。
关于导弹非线性非定常气动力建模问题，目前并未给出相应的解决方法；现有技术都只针对飞机或机翼，主要有：飞机大攻角非定常气动力建模方法（汪清等，飞机大迎角空间机动气动力建模研究，航空学报，2004，25(5)：447-451）。将飞机大攻角气动力分解为三部分，即静态气动力分量、由定常旋转和下洗迟滞产生的气动力增量以及由涡破裂和恢复迟滞引起的气动力增量；采用一阶微分方程描述由非定常效应产生的附加气动力，右端项中的非线性函数
Figure 100002_DEST_PATH_IMAGE001
关于攻角
Figure 468654DEST_PATH_IMAGE002
和攻角速率
Figure 100002_DEST_PATH_IMAGE003
作多项式展开。由于非线性函数
Figure 241831DEST_PATH_IMAGE001
采用多项式展开，气动力模型适用的攻角范围受到较大限制，通常仅适用于
Figure 177426DEST_PATH_IMAGE004
攻角范围，对于
Figure 100002_DEST_PATH_IMAGE005
全攻角范围，由于要求的多项式阶次过高，将导致模型参数辨识困难以及模型泛化性能严重下降等问题。
飞机非定常气动力微分方程模型（Abramov N B, et al. Aircraft dynamics athigh incidence flight with account of unsteady aerodynamic effects. AIAA2004-5274, 2004.）。将飞机大攻角气动力分解为附着流气动力、附着流动导数增量、动态气动力增量三部分；采用一阶微分方程描述动态气动力增量；利用动导数风洞试验数据，采用两步线性回归辨识线性气动力模型；利用大振幅动态风洞试验数据辨识非线性模型中的附加项。文献未明确说明模型参数化方法，也未明确说明其适用的攻角范围，应用实例的攻角范围为
Figure 901800DEST_PATH_IMAGE006
。
发明内容
为解决上述问题，本发明提供了一种导弹非线性非定常气动力微分方程模型辨识方法，该方法应用叠加原理，将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量三个部分，采用一阶微分方程描述非定常气动力增量，并采用二次样条函数近似模型中的非线性函数，从而构建气动力微分方程模型，基于最小二乘准则，采用Gauss-Newton优化算法辨识获取模型参数的估计值。
本发明提供了一种导弹非线性非定常气动力微分方程模型辨识方法，具体技术方案如下：
S1：数据准备：利用风洞试验或CFD计算得到导弹静态气动力和力矩系数、大振幅俯仰振荡气动力和力矩系数时间历程的动态数据表，并经过数据处理后生成气动建模的输入数据文件；
S2：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，采用一阶微分方程描述非定常增量，构建气动力微分方程模型；
S3：将气动力微分方程模型辨识问题转化为动态系统的参数辨识问题；
S4：利用气动数据，基于最小二乘准则，采用Gauss-Newton优化算法辨识获取模型参数的估计值。
进一步的，在步骤S1中，所述数据处理包括如下子步骤：
a.将气动力矩系数的参考点统一变换到导弹的质心处；
b.计算动态数据表的攻角
Figure 100002_DEST_PATH_IMAGE007
、攻角速率
Figure 354778DEST_PATH_IMAGE008
和俯仰角速率
Figure DEST_PATH_IMAGE009
等参数；c.计算无因次时间
Figure 871210DEST_PATH_IMAGE010
、无因次攻角速率
Figure 100002_DEST_PATH_IMAGE011
和无因次俯仰角速率
Figure 321914DEST_PATH_IMAGE012
：
Figure 100002_DEST_PATH_IMAGE013
其中，
Figure 94436DEST_PATH_IMAGE014
为时间；
Figure 882263DEST_PATH_IMAGE016
为弹长；
Figure 100002_DEST_PATH_IMAGE017
为来流速度。
进一步的，步骤S2包括如下子步骤：
S201：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，表示为：
Figure 925305DEST_PATH_IMAGE018
其中，
Figure 100002_DEST_PATH_IMAGE019
为气动力和力矩系数，
Figure 671544DEST_PATH_IMAGE020
；
Figure 100002_DEST_PATH_IMAGE021
为静态气动力系数；
Figure 167248DEST_PATH_IMAGE022
为非定常增量；
Figure 100002_DEST_PATH_IMAGE023
为附着流俯仰阻尼和洗流时差组合导数；
Figure 165291DEST_PATH_IMAGE007
为攻角，
Figure 656315DEST_PATH_IMAGE024
为轴向力系数，
Figure 100002_DEST_PATH_IMAGE025
为法向力系数，
Figure 950287DEST_PATH_IMAGE026
为俯仰力矩系数；
S202：采用一阶微分方程描述非定常增量，表示为：
Figure 100002_DEST_PATH_IMAGE027
其中，
Figure 198865DEST_PATH_IMAGE028
为反映流动迟滞特性的无因次时间常数；
Figure 328495DEST_PATH_IMAGE011
为无因次攻角速率；
Figure 100002_DEST_PATH_IMAGE029
、
Figure 80551DEST_PATH_IMAGE030
、
Figure 100002_DEST_PATH_IMAGE031
表示非线性函数；
S203：对微分方程模型中的非线性函数采用二次样条函数进行近似，表示为：
Figure 434172DEST_PATH_IMAGE032
其中，
Figure 100002_DEST_PATH_IMAGE033
；
Figure 107729DEST_PATH_IMAGE034
为样条系数
Figure 100002_DEST_PATH_IMAGE035
；
Figure 742848DEST_PATH_IMAGE036
为样条节点
Figure 100002_DEST_PATH_IMAGE037
；
Figure 677306DEST_PATH_IMAGE038
为节点数。
进一步的，步骤S3包括如下子步骤：
S301：根据气动力模型，构建气动力微分方程模型辨识问题的动态系统模型，所述动态系统模型包括状态方程和观测方程：
状态方程为：
Figure 100002_DEST_PATH_IMAGE039
观测方程为：
Figure 139511DEST_PATH_IMAGE040
其中，下标“m”表示测量值；
Figure 100002_DEST_PATH_IMAGE041
为测量噪声，
Figure 769207DEST_PATH_IMAGE042
；N为测量数据总点数；
S302：构建待辨识的模型参数
Figure 100002_DEST_PATH_IMAGE043
，表示为：
Figure 302956DEST_PATH_IMAGE044
其中，
Figure 232866DEST_PATH_IMAGE023
为附着流俯仰阻尼和洗流时差组合导数；
Figure 193869DEST_PATH_IMAGE028
为反映流动迟滞特性的无因次时间常数；
Figure 169915DEST_PATH_IMAGE034
为样条系数
Figure 606889DEST_PATH_IMAGE035
。
进一步的，步骤S4包括如下子步骤：
S401：对待辨识的模型参数
Figure 781518DEST_PATH_IMAGE043
赋初值
Figure 100002_DEST_PATH_IMAGE045
；
S402：采用四阶龙格-库塔法对状态方程进行数值积分，得到
Figure 523209DEST_PATH_IMAGE046
，
Figure 100002_DEST_PATH_IMAGE047
；
S403：计算最小二乘准则函数：
Figure 189814DEST_PATH_IMAGE048
S404：采用Gauss-Newton法修正待辨识的模型参数
Figure 268628DEST_PATH_IMAGE043
：
Figure DEST_PATH_IMAGE049
其中，上标“（l）”表示第l步迭代；
Figure 969868DEST_PATH_IMAGE050
为第l步迭代得到的模型参数估计值；
Figure DEST_PATH_IMAGE051
为模型参数
Figure 741515DEST_PATH_IMAGE043
的修正量；
S405：重复执行步骤S402～S404直至收敛（即准则函数
Figure 128372DEST_PATH_IMAGE052
达到最小值），获得模型参数
Figure 745298DEST_PATH_IMAGE043
的估计值。
进一步的，在步骤S4之后，还进行模型预测误差计算，定义均方根误差和相对均方根误差，具体计算如下：
Figure DEST_PATH_IMAGE053
其中，N为建模样本点数；
Figure 301044DEST_PATH_IMAGE054
表示气动力风洞试验或CFD结果；
Figure DEST_PATH_IMAGE055
表示气动力模型预测结果。
本发明的有益效果如下：
1、本发明将导弹大攻角气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量三个部分，采用一阶微分方程描述非定常气动力增量，从而建立导弹非线性非定常气动力微分方程模型，适用于导弹非线性非定常气动力建模，针对导弹大攻角机动场景。
2、本发明基于建立的非线性非定常气动力微分方程模型，采用二次样条函数近似模型中关于攻角的非线性函数，使气动力模型适用的攻角范围更宽，能够适用于
Figure 243592DEST_PATH_IMAGE056
全攻角范围，同时也避免了多项式需要高阶次展开造成模型参数辨识困难以及模型泛化性能严重下降等问题。
3、本发明的模型参数辨识方法更鲁棒，能够保证辨识结果具有物理意义，避免两步线性回归辨识方法在动导数随减缩频率变化规律不好或动导数/动态气动数据存在较大误差的情况下，可能导致模型中的时间常数出现小于0的情况，进而避免模型预测结果发散。
附图说明
图1为本发明的导弹非线性非定常气动力微分方程模型辨识方法的流程图；
图2为实施例中的模型参数辨识过程的前部分流程图；
图3为实施例中的模型参数辨识过程的后部分流程图。
具体实施方式
在下面的描述中对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例
本发明的实施例公开了一种导弹非线性非定常气动力微分方程模型辨识方法，本实施例中，以导弹布局气动力CFD计算数据建立其非线性非定常气动力微分方程模型，辨识模型参数，此作为具体实例，进行如下说明。
如图1所示，导弹非线性非定常气动力微分方程模型辨识方法具体步骤流程如下：
S1：数据准备：利用风洞试验或CFD计算得到导弹静态气动力和力矩系数、大振幅俯仰振荡气动力和力矩系数时间历程的动态数据表，并经过数据处理后生成气动建模的输入数据文件；
所述计算数据包括静态气动力和力矩系数以及大振幅俯仰振荡气动力和力矩系数时间历程；
本实施例中，采用导弹布局气动力CFD计算数据，计算状态如下表所示：
表1：导弹布局气动力CFD计算状态
Figure DEST_PATH_IMAGE057
对所述气动数据依次进行如下处理：
a.将气动力矩系数的参考点统一变换到导弹的质心处；b.计算动态数据表的攻角
Figure 88052DEST_PATH_IMAGE007
、攻角速率
Figure 305406DEST_PATH_IMAGE008
和俯仰角速率
Figure 184501DEST_PATH_IMAGE009
参数；c.计算无因次时间
Figure 829109DEST_PATH_IMAGE010
、无因次攻角速率
Figure 488760DEST_PATH_IMAGE011
和无因次俯仰角速率
Figure 886637DEST_PATH_IMAGE012
；
Figure 744871DEST_PATH_IMAGE013
其中，
Figure 763643DEST_PATH_IMAGE014
为时间；
Figure 582694DEST_PATH_IMAGE016
为弹长；
Figure 141852DEST_PATH_IMAGE017
为来流速度；
经过上述处理后，截取一个周期数据，生成气动建模的输入数据文件。纵向气动力和力矩系数
Figure 57855DEST_PATH_IMAGE058
呈现较明显的非线性和非定常特性，迟滞环随攻角、频率等参数变化的规律性较好。
S2：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，采用一阶微分方程描述非定常增量，构建气动力微分方程模型，具体过程如下：
S201：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量；
由于在攻角平面，导弹气动力和力矩是马赫数、攻角、俯仰角速率、俯仰舵偏角的非线性函数；采用叠加原理，将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量三个部分，如下所示：
Figure DEST_PATH_IMAGE059
其中，
Figure 716370DEST_PATH_IMAGE019
为气动力和力矩系数，
Figure 350613DEST_PATH_IMAGE060
；M为马赫数；
Figure 588828DEST_PATH_IMAGE007
为攻角；
Figure 156075DEST_PATH_IMAGE009
为俯仰角速率；
Figure DEST_PATH_IMAGE061
为俯仰舵偏角；
Figure 687288DEST_PATH_IMAGE062
为无因次俯仰角速率，
Figure DEST_PATH_IMAGE063
；
Figure 74407DEST_PATH_IMAGE016
为弹长；
Figure 116313DEST_PATH_IMAGE017
为来流速度，
Figure 803646DEST_PATH_IMAGE064
为轴向力系数，
Figure DEST_PATH_IMAGE065
为法向力系数，
Figure 741646DEST_PATH_IMAGE026
为俯仰力矩系数；
Figure 147220DEST_PATH_IMAGE066
为静态气动力系数，是马赫数、攻角、舵偏角的非线性函数；
Figure 789554DEST_PATH_IMAGE023
为附着流俯仰阻尼和洗流时差组合导数，是马赫数的函数，通常不考虑攻角和舵偏角的影响；
Figure DEST_PATH_IMAGE067
为非定常增量，在这里为简便起见，省略表示增量的D，是运动历程的非线性泛函，这里用[]表示泛函，以区别于通常的函数，其中舵偏角的影响也可以忽略；
本实施例中，在CFD计算过程中，设定马赫数固定不变，舵偏角
Figure 3497DEST_PATH_IMAGE061
保持为0，则气动力和力矩系数
Figure 705874DEST_PATH_IMAGE019
表示如下：
Figure 709996DEST_PATH_IMAGE068
S202：采用一阶微分方程描述非定常增量；
大振幅振荡风洞试验和CFD数值计算表明，在大攻角下，导弹动态运动的流动状态对运动历程存在“记忆”，从而
Figure 218337DEST_PATH_IMAGE022
与瞬时攻角
Figure 83525DEST_PATH_IMAGE007
及其变化历程有关；并且，气动力的过渡过程近似为指数律收敛过程；因此，本实施例中，采用下列一阶微分方程描述由非定常效应产生的附加气动力
Figure 894486DEST_PATH_IMAGE022
：
Figure DEST_PATH_IMAGE069
保证了在
Figure 477914DEST_PATH_IMAGE070
的定常情况下，非定常效应
Figure DEST_PATH_IMAGE071
；
其中，
Figure 196472DEST_PATH_IMAGE010
为无因次时间，
Figure 588270DEST_PATH_IMAGE072
；
Figure 632449DEST_PATH_IMAGE011
为无因次攻角速率，
Figure DEST_PATH_IMAGE073
；
S203：对微分方程模型中的非线性函数采用二次样条函数进行近似；
对微分方程模型中的非线性函数进行参数化，本实施例中，为保证气动模型适用于
Figure 968753DEST_PATH_IMAGE074
全攻角范围，采用二次样条函数来近似，如下所示：
Figure 927219DEST_PATH_IMAGE032
其中，
Figure DEST_PATH_IMAGE075
；
Figure 298158DEST_PATH_IMAGE034
为样条系数
Figure 185342DEST_PATH_IMAGE035
；
Figure 274521DEST_PATH_IMAGE036
为样条节点
Figure 538143DEST_PATH_IMAGE076
；
Figure 232430DEST_PATH_IMAGE038
为节点数。
S3：将气动力微分方程模型辨识问题转化为动态系统的参数辨识问题，具体过程如下：
S301：根据所述气动力模型，构建气动力微分方程模型辨识问题的动态系统模型；
根据所述气动力模型，构建气动力微分方程模型辨识问题的动态系统模型，所述动态系统方程包括状态方程和观测方程；
所述状态方程为：
Figure DEST_PATH_IMAGE077
所述观测方程为：
Figure 87253DEST_PATH_IMAGE040
其中，下标“m”表示测量值；
Figure 601411DEST_PATH_IMAGE041
为测量噪声，
Figure 793358DEST_PATH_IMAGE042
；N为测量数据总点数；
S302：构建待辨识的模型参数
Figure 515720DEST_PATH_IMAGE043
；
所述模型参数
Figure 869341DEST_PATH_IMAGE043
由气动力模型中的所有未知参数构成，在静态气动力
Figure 870795DEST_PATH_IMAGE078
已知的情况下，待辨识参数为：
Figure DEST_PATH_IMAGE079
其中，
Figure 476220DEST_PATH_IMAGE023
为附着流俯仰阻尼和洗流时差组合导数；
Figure 676257DEST_PATH_IMAGE028
为反映流动迟滞特性的无因次时间常数；
Figure 341725DEST_PATH_IMAGE034
为样条系数
Figure 830475DEST_PATH_IMAGE035
。
S4：利用所述气动数据，基于最小二乘准则，采用Gauss-Newton优化算法辨识获取模型参数的估计值；
图2-图3给出了模型参数辨识过程的流程图，具体过程如下：
S401：对待辨识的模型参数
Figure 364224DEST_PATH_IMAGE043
赋初值
Figure 559713DEST_PATH_IMAGE045
；
本发明对模型参数初值
Figure 458399DEST_PATH_IMAGE045
无严格要求，只需符合物理意义，以下本实施例的取值仅供参考，在本实施例中，轴向力系数
Figure 231183DEST_PATH_IMAGE080
模型参数初值取为：
Figure DEST_PATH_IMAGE081
法向力系数
Figure 676946DEST_PATH_IMAGE082
模型参数初值取为：
Figure DEST_PATH_IMAGE083
俯仰力矩系数
Figure 258100DEST_PATH_IMAGE084
模型参数初值取为：
Figure DEST_PATH_IMAGE085
S402：采用四阶龙格-库塔法对状态方程进行数值积分，得到
Figure 124425DEST_PATH_IMAGE046
，
Figure 259871DEST_PATH_IMAGE086
；
S403：计算最小二乘准则函数；
对于单观测量系统，最大似然准则等价于最小二乘准则，在本实施例中，模型参数辨识采用最小二乘准则：
Figure 338685DEST_PATH_IMAGE048
待辨识参数的估计值为：
Figure DEST_PATH_IMAGE087
S404：采用Gauss-Newton法修正待辨识的模型参数
Figure 774346DEST_PATH_IMAGE043
；
采用Gauss-Newton优化算法从初值
Figure 280414DEST_PATH_IMAGE045
出发，进行迭代求解模型参数的最小二乘估计
Figure 762210DEST_PATH_IMAGE088
，所述待辨识参数的修正公式如下：
Figure DEST_PATH_IMAGE089
其中，上标“
Figure 564424DEST_PATH_IMAGE090
”表示第l步迭代；
Figure 448066DEST_PATH_IMAGE050
为第l步迭代得到的模型参数估计值；
Figure 62719DEST_PATH_IMAGE051
为模型参数
Figure 235074DEST_PATH_IMAGE043
的修正量；
本实施例中，观测灵敏度
Figure DEST_PATH_IMAGE091
采用解析法，由灵敏度方程计算得到，具体过程如下；
观测灵敏度方程：
Figure 124532DEST_PATH_IMAGE092
状态灵敏度方程：
Figure DEST_PATH_IMAGE093
Figure 534785DEST_PATH_IMAGE094
S405：重复执行步骤S402～S404直至收敛（即准则函数
Figure 818874DEST_PATH_IMAGE052
达到最小值），获得模型参数
Figure 275263DEST_PATH_IMAGE043
的估计值。
在步骤S4之后，还进行模型预测误差计算，定义均方根误差和相对均方根误差，具体计算如下：
Figure 233992DEST_PATH_IMAGE053
其中，N为建模样本点数；
Figure 967592DEST_PATH_IMAGE054
表示气动力风洞试验或CFD结果；
Figure 783102DEST_PATH_IMAGE055
表示气动力模型预测结果。
以下通过具体示例对本发明的导弹非线性非定常气动力微分方程模型辨识方法的步骤S1～S4的有益效果进行说明。
示例一：建模方法验证
将所述气动数据分为两组，一组作为建模数据，另一组作为校验数据，以检验气动力微分方程建模方法的可行性。为此，将下列状态的气动数据作为校验数据，其余状态的气动数据作为建模数据：中心攻角
Figure DEST_PATH_IMAGE095
、振幅
Figure 336574DEST_PATH_IMAGE096
、频率
Figure DEST_PATH_IMAGE097
；中心攻角
Figure 567835DEST_PATH_IMAGE098
、振幅
Figure DEST_PATH_IMAGE099
、频率
Figure 280576DEST_PATH_IMAGE100
；中心攻角
Figure DEST_PATH_IMAGE101
、振幅
Figure 378239DEST_PATH_IMAGE096
、频率
Figure 809220DEST_PATH_IMAGE102
；中心攻角
Figure 109751DEST_PATH_IMAGE098
、振幅
Figure DEST_PATH_IMAGE103
、频率
Figure 614682DEST_PATH_IMAGE104
；
利用建模数据，采用本发明的导弹非线性非定常气动力微分方程模型辨识方法，建立
Figure DEST_PATH_IMAGE105
、
Figure 647360DEST_PATH_IMAGE106
、
Figure DEST_PATH_IMAGE107
的微分方程模型，其中攻角节点
Figure 706583DEST_PATH_IMAGE108
取为10、20、30、40、50、60、70、80、90、100、110、120、130、140、150、160、170°；然后，采用所建立的微分方程模型预测校验状态的气动力响应。
建模-校验误差如下表所示：
表2：微分方程模型建模-校验误差
Figure DEST_PATH_IMAGE109
对于校验状态，
Figure 607543DEST_PATH_IMAGE106
和
Figure 668777DEST_PATH_IMAGE107
模型预测结果与CFD数据符合较好，表明微分方程模型具有较好的泛化性能；
Figure 200253DEST_PATH_IMAGE105
模型预测结果与CFD数据差别稍大，主要原因是，俯仰振荡
Figure 605826DEST_PATH_IMAGE105
的CFD数据存在较大的迟滞环，但迟滞环大小随振荡频率变化不明显，换言之，从静态到低频动态运动，
Figure 248160DEST_PATH_IMAGE105
存在“跳跃”现象，模型对这一现象的描述存在一定困难。
示例二：气动力建模
利用所述全部静态和大振幅俯仰振荡CFD计算数据，采用本发明的导弹非线性非定常气动力微分方程模型辨识方法，建立
Figure 930946DEST_PATH_IMAGE105
、
Figure 430060DEST_PATH_IMAGE106
、
Figure 932717DEST_PATH_IMAGE107
的数学模型，其中攻角节点
Figure 378742DEST_PATH_IMAGE110
取为10、20、30、40、50、60、70、80、90、100、110、120、130、140、150、160、170°。模型结构如下所示：
轴向力系数模型
Figure DEST_PATH_IMAGE111
法向力系数模型
Figure 447192DEST_PATH_IMAGE112
俯仰力矩系数模型
Figure DEST_PATH_IMAGE113
微分方程建模误差如下表所示：
表3：微分方程建模误差
Figure 117207DEST_PATH_IMAGE114
其中，
Figure DEST_PATH_IMAGE115
和
Figure 343046DEST_PATH_IMAGE116
模型预测结果与CFD数据符合很好，建模误差很小；
Figure 655078DEST_PATH_IMAGE117
模型预测结果与CFD数据符合稍差，建模误差略大，这是因为CFD数据的
Figure DEST_PATH_IMAGE118
迟滞环随振荡频率变化很小，而模型预测结果的
Figure 515718DEST_PATH_IMAGE119
迟滞环随振荡频率变化较明显。总体来看，所建立的微分方程模型能够较好地描述导弹气动特性随攻角变化的非线性和非定常特性。本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (6)

1.一种导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，包括：
S1：数据准备：利用风洞试验或CFD计算得到导弹静态气动力和力矩系数、大振幅俯仰振荡气动力和力矩系数时间历程的动态数据表，并经过数据处理后生成气动建模的输入数据文件；
S2：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，采用一阶微分方程描述非定常增量，构建气动力微分方程模型；
S3：将气动力微分方程模型辨识问题转化为动态系统的参数辨识问题；
S4：利用气动数据，基于最小二乘准则，采用Gauss-Newton优化算法辨识获取模型参数的估计值。
2.根据权利要求1所述的导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，在步骤S1中，所述数据处理包括如下子步骤：
a.将气动力矩系数的参考点统一变换到导弹的质心处；
b.计算动态数据表的攻角
Figure DEST_PATH_IMAGE001
、攻角速率
Figure 637033DEST_PATH_IMAGE002
和俯仰角速率
Figure DEST_PATH_IMAGE003
参数；
c.计算无因次时间
Figure 275956DEST_PATH_IMAGE004
、无因次攻角速率
Figure DEST_PATH_IMAGE005
和无因次俯仰角速率
Figure 985286DEST_PATH_IMAGE006
：
Figure DEST_PATH_IMAGE007
其中，
Figure 467083DEST_PATH_IMAGE008
为时间；
Figure 21692DEST_PATH_IMAGE010
为弹长；
Figure DEST_PATH_IMAGE011
为来流速度。
3.根据权利要求1所述的导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，步骤S2包括如下子步骤：
S201：将气动力分解为静态气动力、俯仰阻尼和下洗迟滞增量、非定常增量，表示为：
Figure 639755DEST_PATH_IMAGE012
其中，
Figure DEST_PATH_IMAGE013
为气动力和力矩系数，
Figure 275312DEST_PATH_IMAGE014
；
Figure DEST_PATH_IMAGE015
为静态气动力系数；
Figure 854192DEST_PATH_IMAGE016
为非定常增量；
Figure DEST_PATH_IMAGE017
为附着流俯仰阻尼和洗流时差组合导数；
Figure 478071DEST_PATH_IMAGE001
为攻角；
S202：采用一阶微分方程描述非定常增量，表示为：
Figure 481799DEST_PATH_IMAGE018
其中，
Figure DEST_PATH_IMAGE019
为反映流动迟滞特性的无因次时间常数；
Figure 267352DEST_PATH_IMAGE005
为无因次攻角速率；
Figure 927004DEST_PATH_IMAGE020
、
Figure DEST_PATH_IMAGE021
、
Figure 853109DEST_PATH_IMAGE022
表示非线性函数；
S203：对微分方程模型中的非线性函数采用二次样条函数进行近似，表示为：
Figure DEST_PATH_IMAGE023
其中，
Figure 321131DEST_PATH_IMAGE024
；
Figure DEST_PATH_IMAGE025
为样条系数
Figure 808744DEST_PATH_IMAGE026
；
Figure DEST_PATH_IMAGE027
为样条节点
Figure 752429DEST_PATH_IMAGE028
；
Figure DEST_PATH_IMAGE029
为节点数。
4.根据权利要求3所述的导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，步骤S3包括如下子步骤：
S301：根据气动力模型，构建气动力微分方程模型辨识问题的动态系统模型，所述动态系统模型包括状态方程和观测方程：
状态方程为：
Figure 655794DEST_PATH_IMAGE030
观测方程为：
Figure DEST_PATH_IMAGE031
其中，下标“m”表示测量值；
Figure 542104DEST_PATH_IMAGE032
为测量噪声，
Figure DEST_PATH_IMAGE033
；N为测量数据总点数；
S302：构建待辨识的模型参数
Figure 262936DEST_PATH_IMAGE034
，表示为：
Figure DEST_PATH_IMAGE035
其中，
Figure 366021DEST_PATH_IMAGE017
为附着流俯仰阻尼和洗流时差组合导数；
Figure 666552DEST_PATH_IMAGE019
为反映流动迟滞特性的无因次时间常数；
Figure 109166DEST_PATH_IMAGE025
为样条系数
Figure 532057DEST_PATH_IMAGE036
。
5.根据权利要求4所述的导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，步骤S4包括如下子步骤：
S401：对待辨识的模型参数
Figure 325701DEST_PATH_IMAGE034
赋初值
Figure DEST_PATH_IMAGE037
；
S402：采用四阶龙格-库塔法对状态方程进行数值积分，得到
Figure 961081DEST_PATH_IMAGE038
，
Figure DEST_PATH_IMAGE039
；
S403：计算最小二乘准则函数：
Figure 756737DEST_PATH_IMAGE040
S404：采用Gauss-Newton法修正待辨识的模型参数
Figure 350529DEST_PATH_IMAGE034
：
Figure DEST_PATH_IMAGE041
其中，上标“（l）”表示第l步迭代；
Figure 100310DEST_PATH_IMAGE042
为第l步迭代得到的模型参数估计值；
Figure DEST_PATH_IMAGE043
为模型参数
Figure 477065DEST_PATH_IMAGE034
的修正量；
S405：重复执行步骤S402～S404直至收敛（即准则函数
Figure 222167DEST_PATH_IMAGE044
达到最小值），获得模型参数
Figure 721282DEST_PATH_IMAGE034
的估计值。
6.根据权利要求5所述的导弹非线性非定常气动力微分方程模型辨识方法，其特征在于，在步骤S4之后，还进行模型预测误差计算，定义均方根误差和相对均方根误差，具体计算如下：
Figure DEST_PATH_IMAGE045
其中，N为建模样本点数；
Figure 20676DEST_PATH_IMAGE046
表示气动力风洞试验或CFD结果；
Figure DEST_PATH_IMAGE047
表示气动力模型预测结果。
CN202210825423.1A 2022-07-14 2022-07-14 导弹非线性非定常气动力微分方程模型辨识方法 Active CN114896830B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210825423.1A CN114896830B (zh) 2022-07-14 2022-07-14 导弹非线性非定常气动力微分方程模型辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210825423.1A CN114896830B (zh) 2022-07-14 2022-07-14 导弹非线性非定常气动力微分方程模型辨识方法

Publications (2)

Publication Number Publication Date
CN114896830A true CN114896830A (zh) 2022-08-12
CN114896830B CN114896830B (zh) 2022-11-08

Family

ID=82729993

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210825423.1A Active CN114896830B (zh) 2022-07-14 2022-07-14 导弹非线性非定常气动力微分方程模型辨识方法

Country Status (1)

Country Link
CN (1) CN114896830B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116956471A (zh) * 2023-09-19 2023-10-27 中国空气动力研究与发展中心计算空气动力研究所 大型运输机的气动力预测方法、装置、设备及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101587328A (zh) * 2008-05-23 2009-11-25 朱豫才 非线性过程动态模型辨识装置
US20140136586A1 (en) * 2012-11-13 2014-05-15 Korea Institute Of Geoscience And Mineral Resources Method of seeking semianalytical solutions to multispecies transport equations coupled with sequential first-order reactions
CN112800543A (zh) * 2021-01-27 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 一种基于改进Goman模型的非线性非定常气动力建模方法
CN112904898A (zh) * 2021-01-28 2021-06-04 上海机电工程研究所 旋转弹箭非定常气动响应特性评估方法和系统
CN113919194A (zh) * 2021-09-07 2022-01-11 中国空气动力研究与发展中心计算空气动力研究所 一种基于滤波误差法的非线性飞行动力学系统辨识方法
WO2022033183A1 (zh) * 2020-08-13 2022-02-17 重庆邮电大学 动静态数据混合驱动的hammerstein非线性工业系统简约灰箱子空间辨识方法
CN114611416A (zh) * 2022-05-12 2022-06-10 中国空气动力研究与发展中心计算空气动力研究所 导弹非线性非定常气动特性ls-svm建模方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101587328A (zh) * 2008-05-23 2009-11-25 朱豫才 非线性过程动态模型辨识装置
US20140136586A1 (en) * 2012-11-13 2014-05-15 Korea Institute Of Geoscience And Mineral Resources Method of seeking semianalytical solutions to multispecies transport equations coupled with sequential first-order reactions
WO2022033183A1 (zh) * 2020-08-13 2022-02-17 重庆邮电大学 动静态数据混合驱动的hammerstein非线性工业系统简约灰箱子空间辨识方法
CN112800543A (zh) * 2021-01-27 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 一种基于改进Goman模型的非线性非定常气动力建模方法
CN112904898A (zh) * 2021-01-28 2021-06-04 上海机电工程研究所 旋转弹箭非定常气动响应特性评估方法和系统
CN113919194A (zh) * 2021-09-07 2022-01-11 中国空气动力研究与发展中心计算空气动力研究所 一种基于滤波误差法的非线性飞行动力学系统辨识方法
CN114611416A (zh) * 2022-05-12 2022-06-10 中国空气动力研究与发展中心计算空气动力研究所 导弹非线性非定常气动特性ls-svm建模方法

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
N.ABRAMOV 等: "Aircraft Dynamics at High Incidence Flight with Account of Unsteady Aerodynamic Effects", 《AIAA ATMOSPHERIC FLIGHT MECHANICS CONFERENCE AND EXHIBIT》 *
QING WANG 等: "Unsteady aerodynamics modeling for flight dynamics application", 《ACTA MECHANICA SINICA》 *
张天姣 等: "基于最大似然法的风洞自由飞试验气动力参数辨识技术研究", 《实验流体力学》 *
曹德一 等: "基于气动力数据集的气动偏差建模与辨识", 《战术导弹技术》 *
杨文 等: "面向复杂构型飞机的非定常气动力建模与辨识", 《航空学报》 *
汪清 等: "飞机大攻角空间机动气动力建模研究", 《航空学报》 *
汪清 等: "飞机大攻角非定常气动力建模与辨识", 《航空学报》 *
蔡金狮: "飞行器气动参数辨识进展", 《力学进展》 *
龚正 等: "非定常气动力非线性微分方程建模方法", 《航空学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116956471A (zh) * 2023-09-19 2023-10-27 中国空气动力研究与发展中心计算空气动力研究所 大型运输机的气动力预测方法、装置、设备及介质
CN116956471B (zh) * 2023-09-19 2024-01-12 中国空气动力研究与发展中心计算空气动力研究所 大型运输机的气动力预测方法、装置、设备及介质

Also Published As

Publication number Publication date
CN114896830B (zh) 2022-11-08

Similar Documents

Publication Publication Date Title
CN107220403B (zh) 飞行器弹性模态的控制关联建模方法
Ghoreyshi et al. Transonic aerodynamic load modeling of X-31 aircraft pitching motions
CN114611416B (zh) 导弹非线性非定常气动特性ls-svm建模方法
CN115238397B (zh) 高超声速飞行器热环境计算方法、装置和计算机设备
CN114896830B (zh) 导弹非线性非定常气动力微分方程模型辨识方法
CN113609600B (zh) 一种适用于飞行器多体分离相容性度量与表征方法
Hallissy et al. High-fidelity aeroelastic analysis of very flexible aircraft
CN110287505B (zh) 飞行器稳定性分析方法
CN107976908B (zh) 一种飞行器耦合动稳定性特征分析方法
CN112068444B (zh) 一种采用非线性自适应滑模的飞行器攻角控制方法
Yang et al. An improved nonlinear reduced-order modeling for transonic aeroelastic systems
Chauhan et al. Review of aerodynamic parameter estimation techniques
CN113392599A (zh) 一种弹性飞行器动响应的确定方法
Grauer Real-time data-compatibility analysis using output-error parameter estimation
CN113221237B (zh) 一种基于降阶建模的大迎角颤振分析方法
CN115793449B (zh) 直升机飞行姿态控制器设计方法及装置、设备、存储介质
Artola et al. Proof of Concept for a Hardware-in-the-Loop Nonlinear Control Framework for Very Flexible Aircraft
CN115840992A (zh) 一种弹性飞行器飞行仿真方法、系统、计算机存储介质及终端
Liu et al. Unsteady vibration aerodynamic modeling and evaluation of dynamic derivatives using computational fluid dynamics
Lee et al. Predicting aerodynamic rotor-fuselage interactions by using unstructured meshes
Zhang et al. Nonlinear aerodynamics and nonlinear structures interaction for F-16 limit cycle oscillation prediction
Guner et al. An Approximate Finite State Dynamic Wake Model for Predictions of Inflow below the Rotor
Aksu Aerodynamic parameter estimation of a missile without wind angle measurements
CN117521561B (zh) 巡航飞行器的气动力与推力在线预示方法
Reveles et al. Application of a Novel Stability Method to CFD-Based Whirl Flutter Analysis

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