CN114818550B - 一种飞机振动试验中时变气动载荷地面等效模拟方法 - Google Patents

一种飞机振动试验中时变气动载荷地面等效模拟方法 Download PDF

Info

Publication number
CN114818550B
CN114818550B CN202210759016.5A CN202210759016A CN114818550B CN 114818550 B CN114818550 B CN 114818550B CN 202210759016 A CN202210759016 A CN 202210759016A CN 114818550 B CN114818550 B CN 114818550B
Authority
CN
China
Prior art keywords
time
point
function
model
order
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
CN202210759016.5A
Other languages
English (en)
Other versions
CN114818550A (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.)
AVIC Aircraft Strength Research Institute
Original Assignee
AVIC Aircraft Strength Research Institute
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 AVIC Aircraft Strength Research Institute filed Critical AVIC Aircraft Strength Research Institute
Priority to CN202210759016.5A priority Critical patent/CN114818550B/zh
Publication of CN114818550A publication Critical patent/CN114818550A/zh
Application granted granted Critical
Publication of CN114818550B publication Critical patent/CN114818550B/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
    • 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
    • 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
    • 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

Landscapes

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

Abstract

本发明提供了一种飞机振动试验中时变气动载荷地面等效模拟方法，属于飞机振动试验技术领域。方法包括以下步骤：S1、基于时间冻结假设得到适用于飞机振动试验整个时变历程的非定常气动力快速计算模型；S2、对非定常气动力快速计算模型进行仿真验证，得到仿真验证后的非定常气动力快速计算模型；S3、采用仿真验证后的非定常气动力快速计算模型进行飞机振动试验中时变气动载荷地面等效模拟。本发明解决了当前飞机振动试验模拟只能应用于定常系统的问题，具有有效完善当前飞机振动试验模拟技术，推动飞机振动试验模拟工程应用进展的优点。

Description

一种飞机振动试验中时变气动载荷地面等效模拟方法
技术领域
本发明涉及飞机振动试验技术领域，具体是涉及一种飞机振动试验中时变气动载荷地面等效模拟方法。
背景技术
飞机振动试验模拟技术作为近年来发展的一种以真实结构作为验证对象的半物理仿真颤振验证技术，具有风险小、周期短、精度高的优势。该技术利用传感器测量结构的振动响应信号，将其代入预先计算得到的非定常气动力模型获得结构的等效气动力，再采用激振器将等效气动力实时加载在试验结构上，形成闭环系统在地面实现对结构颤振现象的模拟。通过对试验结构施加初始扰动并观察其响应，即可开展针对真实结构颤振边界测试验证。
在上述飞机振动试验模拟的步骤中，考虑到常规结构的颤振发生频率，试验系统完成一次循环计算的周期应不高于0.002秒。为了满足飞机振动试验模拟系统对计算效率的需求，飞机振动试验模拟设计时依据结构动力学特性对频域气动力影响系数矩阵进行降阶缩聚处理，进而采用有理函数近似以及拉氏变换得到由数个低阶矩阵组成的时域非定常气动力模型。
上述建模方法能够成立的一个大前提是结构的动力学特性不发生改变，而在实际飞行器结构中无论是舵面偏转影响或气动热效应导致的结构热模态时变特性，均会使现有非定常气动力计算方法不再完全适用。
本发明提出一种基于Kriging代理模型的一种飞机振动试验中时变气动载荷地面等效模拟方法。本发明依据时间冻结假设建立多个离散时间点的非定常气动力模型后，采用Kriging代理模型对离散时间点的非定常气动力模型进行拟合建模，使得新得到的气动力模型在保证系统运行效率的同时，针对动力学特性时变的结构具有更高的计算精度，有效扩充了地面颤振模拟试验的应用范围。
发明内容
针对当前地面颤振模拟试验只能应用于定常系统的缺陷，本发明给出了可应用于动力学特性时变结构的非定常气动力快速计算方法，有效完善当前飞机振动试验模拟技术，推动飞机振动试验模拟工程应用进展。
为解决上述问题，本发明的技术方案如下：
本发明提供一种飞机振动试验中时变气动载荷地面等效模拟方法，包括以下步骤：
S1、基于时间冻结假设，通过将多个离散时间点的非定常气动力快速计算模型进行Kriging代理模型拟合，得到适用于飞机振动试验整个时变历程的非定常气动力快速计算模型，具体包括以下步骤：
S1-1、建立试验件控制点承受非定常气动力的计算公式，
S1-2、建立离散时间点的试验件控制点应承受非定常气动力计算模型，
S1-3、结合Kriging代理模型建立对应随机函数
Kriging代理模型设定
Figure DEST_PATH_IMAGE002
为对应随机函数
Figure DEST_PATH_IMAGE004
的均方误差最小无偏估计结果，
Figure DEST_PATH_IMAGE005
由采样时间点数据建立的全局近似模型以及随机过程两部分组成，
Figure 756649DEST_PATH_IMAGE004
的表达公式如下：
Figure DEST_PATH_IMAGE007
（3）
上式中，
Figure DEST_PATH_IMAGE008
为对应随机函数，
Figure DEST_PATH_IMAGE010
为已知的基函数，
Figure DEST_PATH_IMAGE012
为基函数系数，
Figure DEST_PATH_IMAGE014
为设计空间内由采样点求解得到的全局近似模型，
Figure DEST_PATH_IMAGE016
是均值为0、方差为
Figure DEST_PATH_IMAGE018
的静态随机过程函数，
S1-4、建立静态随机过程函数在任意两个采样时间点间的协方差计算模型
在飞机地面颤振试验定义的时间范围内，得到的静态随机过程函数在任意两个采样时间点间的协方差计算模型如下：
Figure DEST_PATH_IMAGE020
（4）
上式中，
Figure DEST_PATH_IMAGE021
为方差，
Figure DEST_PATH_IMAGE023
为第一采样时间点的静态随机过程函数，
Figure DEST_PATH_IMAGE025
为第二采样时间点的静态随机过程函数，
Figure DEST_PATH_IMAGE027
为第一采样时间点，
Figure DEST_PATH_IMAGE029
为第二采样时间点，
Figure DEST_PATH_IMAGE031
为相关函数，相关函数的形式为高斯函数或三次样条函数，
S1-5、建立非定常气动力快速计算模型，即Kriging代理模型在任意时刻
Figure DEST_PATH_IMAGE033
的预测值计算模型
通过公式（4）构建形如式（3）的样本点数据列向量定义公式如下：
Figure DEST_PATH_IMAGE035
（5）
上式中，
Figure DEST_PATH_IMAGE036
为样本点数据列向量，
Figure 722725DEST_PATH_IMAGE012
为基函数系数，
Figure DEST_PATH_IMAGE038
为n个采样时间点的样本点数据，
Figure DEST_PATH_IMAGE040
为基函数列向量，
Figure DEST_PATH_IMAGE042
，
Figure DEST_PATH_IMAGE044
为随机函数列向量，
Figure DEST_PATH_IMAGE046
，
根据公式（4）中的相关函数
Figure DEST_PATH_IMAGE047
定义
Figure DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵
Figure DEST_PATH_IMAGE051
，以及表示任意时间点与所有采样时间点相关性的相关向量
Figure DEST_PATH_IMAGE053
，相关矩阵
Figure 627096DEST_PATH_IMAGE051
的计算公式如下：
Figure DEST_PATH_IMAGE055
（6）
上式中，
Figure 553463DEST_PATH_IMAGE051
为
Figure 367835DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 100002_DEST_PATH_IMAGE057
为
Figure 393560DEST_PATH_IMAGE049
个时间采样点中第i个采样时间点，
Figure 100002_DEST_PATH_IMAGE059
为
Figure DEST_PATH_IMAGE060
个时间采样点中第j个采样时间点，
Figure DEST_PATH_IMAGE062
为第i个采样时间点和第j个采样时间点之间的相关函数，
Figure DEST_PATH_IMAGE063
表示大小为
Figure DEST_PATH_IMAGE064
的矩阵，
采样时间点之外的时间节点处的矩阵函数
Figure DEST_PATH_IMAGE066
由已知采样点函数值的线性加权得到：
Figure DEST_PATH_IMAGE068
（8）
上式中，
Figure DEST_PATH_IMAGE069
为加权行向量，
Figure 581965DEST_PATH_IMAGE066
为采样时间点外的时间节点函数，
Figure DEST_PATH_IMAGE070
为样本点数据列向量，
为满足Kriging代理模型无偏估计
Figure DEST_PATH_IMAGE072
的要求，根据公式（5）与公式（8）可得基函数与基函数列向量关系的计算公式：
Figure DEST_PATH_IMAGE074
（9）
上式中，
Figure DEST_PATH_IMAGE075
为基函数列向量，
Figure DEST_PATH_IMAGE076
为基函数行向量，
Figure 100002_DEST_PATH_IMAGE077
为加权行向量，
Figure 660428DEST_PATH_IMAGE012
为基函数系数，
Figure DEST_PATH_IMAGE079
为数学期望，
Figure 594886DEST_PATH_IMAGE066
为采样时间点外的时间节点函数，
Figure DEST_PATH_IMAGE080
为样本点数据列向量，
根据公式（3）、公式（8）与公式（9）得到Kriging代理模型的预测值误差计算公式如下：
Figure DEST_PATH_IMAGE082
（10）
上式中，
Figure DEST_PATH_IMAGE083
为采样时间点外的时间节点函数值，
Figure DEST_PATH_IMAGE085
为
Figure DEST_PATH_IMAGE086
通过Kriging代理模型得到的对应随机函数，
Figure DEST_PATH_IMAGE087
为任一时变气动力影响系数矩阵中的元素，
Figure DEST_PATH_IMAGE088
为随机函数列向量，
Figure DEST_PATH_IMAGE089
是均值为0、方差为
Figure 447304DEST_PATH_IMAGE021
的静态随机过程函数，
Figure DEST_PATH_IMAGE091
为加权行向量，
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，即为在
Figure DEST_PATH_IMAGE093
的前提下，寻找满足均方误差最小的加权列向量
Figure DEST_PATH_IMAGE095
，为此引入拉格朗日乘数
Figure DEST_PATH_IMAGE097
构建拉格朗日函数
Figure DEST_PATH_IMAGE099
如下：
Figure DEST_PATH_IMAGE101
（12）
上式中，
Figure DEST_PATH_IMAGE102
为拉格朗日函数，
Figure DEST_PATH_IMAGE103
为方差，
Figure 404896DEST_PATH_IMAGE095
为加权列向量，
Figure 594437DEST_PATH_IMAGE091
为加权行向量，
Figure 100002_DEST_PATH_IMAGE104
为相关向量，
Figure 117823DEST_PATH_IMAGE051
为
Figure DEST_PATH_IMAGE105
个时间采样点之间相关性的相关矩阵，
Figure DEST_PATH_IMAGE107
为已知的基函数，
Figure DEST_PATH_IMAGE109
为拉格朗日乘数，
Figure DEST_PATH_IMAGE111
为基函数行向量，
分别对公式（12）中加权列向量
Figure 219771DEST_PATH_IMAGE095
和拉格朗日乘数
Figure 461396DEST_PATH_IMAGE109
求偏导数，得到满足要求的加权列向量
Figure 100002_DEST_PATH_IMAGE112
计算公式为：
Figure DEST_PATH_IMAGE114
（13）
上式中，
Figure 189050DEST_PATH_IMAGE051
为
Figure 832521DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 574212DEST_PATH_IMAGE104
为相关向量，
Figure DEST_PATH_IMAGE115
为已知的基函数，
Figure 303133DEST_PATH_IMAGE075
为基函数列向量，
Figure 381948DEST_PATH_IMAGE111
为基函数行向量，
将公式（13）带入公式（8）得到非定常气动力快速计算模型，即Kriging代理模型在任意时刻
Figure DEST_PATH_IMAGE116
的预测值计算模型如下：
Figure DEST_PATH_IMAGE118
（14）
上式中，
Figure 145504DEST_PATH_IMAGE066
为采样时间点外的时间节点函数，
Figure 575873DEST_PATH_IMAGE080
为样本点数据列向量，
Figure DEST_PATH_IMAGE119
为
Figure 526512DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 612279DEST_PATH_IMAGE104
为相关向量，
Figure 964763DEST_PATH_IMAGE115
为已知的基函数，
Figure 907312DEST_PATH_IMAGE075
为基函数列向量，
Figure 876405DEST_PATH_IMAGE111
为基函数行向量；
S2、对步骤S1建立的非定常气动力快速计算模型进行仿真验证，得到仿真验证后的非定常气动力快速计算模型，若非定常气动力快速计算模型的准确性满足要求，得到仿真验证后的非定常气动力快速计算模型，若不满足要求，加密采样时间点后再次进行仿真验证，直至非定常气动力快速计算模型的准确性满足要求；
S3、采用步骤S2仿真验证后的非定常气动力快速计算模型进行飞机振动试验中时变气动载荷地面等效模拟。
由公式（5）与公式（13）可知，公式（14）中除
Figure 297022DEST_PATH_IMAGE095
表达式中的
Figure 769591DEST_PATH_IMAGE104
外其他参数只与采样时间点有关，因此可以提前计算得到，而
Figure 883041DEST_PATH_IMAGE104
需要将当前时刻
Figure DEST_PATH_IMAGE121
代入相关函数公式（7）求解，模型的输入参数为时间，因此基于公式（14）可以获得任意时刻下与结构特性相匹配的气动力。
进一步地，步骤S1-1具体包括以下内容：
根据飞机地面颤振试验中试验件的结构振动响应计算试验件控制点承受的非定常气动力，试验件控制点承受非定常气动力的计算公式如下：
Figure DEST_PATH_IMAGE123
（1）
上式中，
Figure DEST_PATH_IMAGE125
为试验件控制点应承受的非定常气动力，
Figure DEST_PATH_IMAGE127
为来流动压，
Figure DEST_PATH_IMAGE129
为降阶后的气动力影响系数矩阵，
Figure DEST_PATH_IMAGE130
是通过与马赫数相匹配的非定常气动力理论计算后并依据结构动力学特性进行降阶处理得到的，
Figure DEST_PATH_IMAGE132
表示控制点位移，
Figure DEST_PATH_IMAGE134
为减缩频率。
进一步地，减缩频率
Figure 932905DEST_PATH_IMAGE134
的计算公式如下：
Figure DEST_PATH_IMAGE136
上式中，
Figure DEST_PATH_IMAGE137
为减缩频率，
Figure DEST_PATH_IMAGE139
为简谐运动的圆频率，
Figure DEST_PATH_IMAGE141
为参考长度，对于飞机的翼面类结构，参考长度取根弦长度的一半，
Figure DEST_PATH_IMAGE143
为来流速度。
进一步地，步骤S1-2具体包括以下内容：
由于实际飞机地面颤振试验需要试验件控制点应承受非定常气动力的时域关系式，将公式（1）通过最小状态法有理函数近似延拓至拉氏域，再通过拉氏反变换得到离散时间点的试验件控制点应承受非定常气动力计算模型如下：
Figure DEST_PATH_IMAGE145
（2）
上式中，
Figure DEST_PATH_IMAGE147
为与时间相关的试验件控制点应承受非定常气动力，
Figure DEST_PATH_IMAGE149
为位移气动力影响系数矩阵，
Figure DEST_PATH_IMAGE151
为速度气动力影响系数矩阵，
Figure DEST_PATH_IMAGE153
为加速度气动力影响系数矩阵，
Figure DEST_PATH_IMAGE155
为第一辅助气动力影响系数矩阵，
Figure DEST_PATH_IMAGE157
为第二辅助气动力影响系数矩阵，
Figure DEST_PATH_IMAGE159
为第三辅助气动力影响系数矩阵，
Figure DEST_PATH_IMAGE161
为控制点位移，
Figure DEST_PATH_IMAGE163
为控制点速度，
Figure DEST_PATH_IMAGE165
为控制点加速度，
Figure DEST_PATH_IMAGE166
为来流动压，
Figure DEST_PATH_IMAGE167
为来流速度，
Figure 59343DEST_PATH_IMAGE141
为参考长度，对于飞机的翼面类结构，参考长度取根弦长度的一半。
由于公式（2）的建模方法只能适用于时不变系统，而在部分实际飞机地面颤振试验各个时间节点的气动力影响系数矩阵数值均会发生变化，因此，需在公式（2）的基础上引入Kriging代理模型，将公式（1）和公式（2）的求解过程替换为简单矩阵运算，实现针对时变动力学特性结构的非定常气动力快速计算。
进一步地，步骤S1-5中：
相关向量
Figure DEST_PATH_IMAGE168
的计算公式如下：
Figure DEST_PATH_IMAGE170
（7）
上式中，
Figure 386419DEST_PATH_IMAGE121
为
Figure DEST_PATH_IMAGE171
个时间采样点中任意一个样本点，
Figure DEST_PATH_IMAGE173
，
Figure DEST_PATH_IMAGE175
为相关向量，
Figure DEST_PATH_IMAGE177
为
Figure DEST_PATH_IMAGE178
个时间采样点中任意一个样本点
Figure 529825DEST_PATH_IMAGE121
与第1个采样时间点的相关函数，
Figure DEST_PATH_IMAGE180
为
Figure DEST_PATH_IMAGE181
个时间采样点中任意一个样本点
Figure 942352DEST_PATH_IMAGE121
与第2个采样时间点的相关函数，
Figure DEST_PATH_IMAGE183
为
Figure 501509DEST_PATH_IMAGE181
个时间采样点中任意一个样本点
Figure 683092DEST_PATH_IMAGE121
与第
Figure DEST_PATH_IMAGE184
个采样时间点的相关函数。
更进一步地，步骤S1-5还包括：
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，通过公式（10）得到Kriging代理模型的均方误差计算公式如下：
Figure DEST_PATH_IMAGE186
（11）
上式中，
Figure DEST_PATH_IMAGE188
为Kriging代理模型的预均方误差，
Figure 544868DEST_PATH_IMAGE083
为采样时间点外的时间节点函数，
Figure 444691DEST_PATH_IMAGE085
为
Figure DEST_PATH_IMAGE189
通过Kriging代理模型得到的对应随机函数，
Figure DEST_PATH_IMAGE190
为任一时变气动力影响系数矩阵中的元素，
Figure DEST_PATH_IMAGE191
为加权行向量，
Figure DEST_PATH_IMAGE192
为数学期望，
Figure 869856DEST_PATH_IMAGE018
为方差，
Figure DEST_PATH_IMAGE194
为加权列向量，
Figure 171525DEST_PATH_IMAGE104
为相关向量，
Figure DEST_PATH_IMAGE195
为
Figure 797678DEST_PATH_IMAGE105
个时间采样点之间相关性的相关矩阵。
优选地，步骤S2包括以下步骤：
S2-1、在时变历程选取数个采样时间点；
S2-2、将数个采样时间点数据导入步骤S1，得到步骤S1-2计算的每个采样时间点的非定常气动力模型、步骤S1-5构建的非定常气动力快速计算模型；
S2-3、采用matlab.simulink开展非定常气动力模型和非定常气动力快速计算模型变温地面颤振系统仿真，仿真对象为高超声速舵面结构，时变历程为高超声速舵面结构由初始温度场升至稳态温度场的过程，得到每个采样时间点的非定常气动力模型的定常颤振仿真以及整个时变历程的仿真、非定常气动力快速计算模型对每个采样时间点的定常颤振仿真以及整个时变历程的仿真，并分别与每个采样时间点的颤振频域数值计算结果对比；
S2-4、步骤S2-3的对比结果为不一致时，加密采样时间点后，回到步骤S2-2，直至步骤S2-3的对比结果为一致。
优选地，初始温度场的温度范围为45~55℃，稳态温度场的温度范围为150~490℃。
优选地，步骤S2-3还包括：
分别向每个采样时间点的非定常气动力模型和非定常气动力快速计算模型中引入颤振临界速度追踪模块，颤振临界速度追踪模块根据仿真对象的响应变化趋势自动调整所述非定常气动力模型和所述非定常气动力快速计算模型中的风速，使仿真对象保持在颤振临界状态，实现对时变颤振速度边界的实时辨识与记录。
进一步优选地，步骤S3具体包括以下步骤：
S3-1、对飞机地面颤振试验颤振计算中的气动平面和结构间的插值点进行缩聚，通过优化方法选取最优插值点，最优插值点为：激振点和拾振点，将最优插值点作为试验件控制点；
S3-2、在试验件上方布设加热器，在试验件底部固定激振器，与激振器电性连接的功率放大器，在试验件上试验件控制点布置激光位移传感器和加速度传感器，功率放大器、激光位移传感器、加速度传感器与加载控制器电性连接；
S3-3、通过激光位移传感器和加速度传感器将获取的试验件控制点响应信号传入加载控制器，加载控制器将试验件控制点响应信号代入步骤S2仿真验证后的非定常气动力快速计算模型，非定常气动力快速计算模型计算试验件的等效气动力；
S3-4、加载控制器通过功率放大器控制激振器将等效气动力实时加载在试验件上，加载完成后回到步骤S3-3，直至飞机振动试验结束。
本发明的有益效果是：
（1）本发明将Kriging代理模型引入地面颤振模拟试验，依据时间冻结假设获取采样时间点数据并提出了根据采样时间点数据的完整建模方法，获得了可应用于动力学特性时变结构的非定常气动力快速计算方法，增强了此类试验系统对于时变颤振系统的处理能力；
（2）本发明所需工作模块只是在现有时不变系统所采用非定常气动力模块的基础上增加数个低阶矩阵运算，对地面颤振模拟试验运算效率影响极小，能够满足地面颤振模拟试验对于气动力运算实时性的要求，此外，由于Kriging模型在应对非线性系统时具有明显优势，因此在计算非定常气动力时拟合精度较高，能够保证在试验件时变历程中都保持较高的非定常气动力计算精度。
附图说明
图1是实施例1一种飞机振动试验中时变气动载荷地面等效模拟方法流程图；
图2是实施例1中离散时间点的试验件控制点应承受非定常气动力计算模型结构图；
图3是实施例1步骤S2-1中仿真对象的温度分布图；
图4是实施例1步骤S2-3中仿真及数值计算结果对比图；
图5是实施例1中步骤S1流程图；
图6是实施例3最优插值点位置示意图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
在本发明实施例中使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本发明。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义，“多种”一般包含至少两种。
实施例1
本实施例为一种飞机振动试验中时变气动载荷地面等效模拟方法，如图1所示，包括以下步骤：
S1、基于时间冻结假设，通过将多个离散时间点的非定常气动力快速计算模型进行Kriging代理模型拟合，得到适用于飞机振动试验整个时变历程的非定常气动力快速计算模型，如图5所示，具体包括以下步骤：
S1-1、建立试验件控制点承受非定常气动力的计算公式
根据飞机地面颤振试验中试验件的结构振动响应计算试验件控制点承受的非定常气动力，试验件控制点承受非定常气动力的计算公式如下：
Figure 184797DEST_PATH_IMAGE123
（1）
上式中，
Figure DEST_PATH_IMAGE196
为试验件控制点应承受的非定常气动力，
Figure 7128DEST_PATH_IMAGE127
为来流动压，
Figure DEST_PATH_IMAGE197
为降阶后的气动力影响系数矩阵，
Figure DEST_PATH_IMAGE198
是通过与马赫数相匹配的非定常气动力理论计算后并依据结构动力学特性进行降阶处理得到的，
Figure DEST_PATH_IMAGE199
表示控制点位移，
Figure 304249DEST_PATH_IMAGE134
为减缩频率，减缩频率
Figure 101303DEST_PATH_IMAGE134
的计算公式如下：
Figure DEST_PATH_IMAGE200
上式中，
Figure 975719DEST_PATH_IMAGE137
为减缩频率，
Figure 149211DEST_PATH_IMAGE139
为简谐运动的圆频率，
Figure 818614DEST_PATH_IMAGE141
为参考长度，对于飞机的翼面类结构，参考长度取根弦长度的一半，
Figure 786570DEST_PATH_IMAGE167
为来流速度，
S1-2、建立离散时间点的试验件控制点应承受非定常气动力计算模型
由于实际飞机地面颤振试验需要试验件控制点应承受非定常气动力的时域关系式，将公式（1）通过最小状态法有理函数近似延拓至拉氏域，再通过拉氏反变换得到离散时间点的试验件控制点应承受非定常气动力计算模型如下：
Figure 882702DEST_PATH_IMAGE145
（2）
上式中，
Figure DEST_PATH_IMAGE201
为与时间相关的试验件控制点应承受非定常气动力，
Figure 797569DEST_PATH_IMAGE149
为位移气动力影响系数矩阵，
Figure 662756DEST_PATH_IMAGE151
为速度气动力影响系数矩阵，
Figure 801614DEST_PATH_IMAGE153
为加速度气动力影响系数矩阵，
Figure 916200DEST_PATH_IMAGE155
为第一辅助气动力影响系数矩阵，
Figure 165916DEST_PATH_IMAGE157
为第二辅助气动力影响系数矩阵，
Figure DEST_PATH_IMAGE202
为第三辅助气动力影响系数矩阵，
Figure 541403DEST_PATH_IMAGE161
为控制点位移，
Figure 851161DEST_PATH_IMAGE163
为控制点速度，
Figure 921885DEST_PATH_IMAGE165
为控制点加速度，
Figure 178554DEST_PATH_IMAGE166
为来流动压，
Figure 18334DEST_PATH_IMAGE167
为来流速度，
Figure 498994DEST_PATH_IMAGE141
为参考长度，对于飞机的翼面类结构，参考长度取根弦长度的一半，离散时间点的试验件控制点应承受非定常气动力计算模型结构如图2所示，
S1-3、结合Kriging代理模型建立对应随机函数
Kriging代理模型设定
Figure 791435DEST_PATH_IMAGE002
为对应随机函数
Figure 179691DEST_PATH_IMAGE004
的均方误差最小无偏估计结果，
Figure 873978DEST_PATH_IMAGE005
由采样时间点数据建立的全局近似模型以及随机过程两部分组成，
Figure 525539DEST_PATH_IMAGE004
的表达公式如下：
Figure 961068DEST_PATH_IMAGE007
（3）
上式中，
Figure 356278DEST_PATH_IMAGE008
为对应随机函数，
Figure 905071DEST_PATH_IMAGE010
为已知的基函数，
Figure 524271DEST_PATH_IMAGE012
为基函数系数，
Figure 791304DEST_PATH_IMAGE014
为设计空间内由采样点求解得到的全局近似模型，
Figure 990204DEST_PATH_IMAGE016
是均值为0、方差为
Figure 393504DEST_PATH_IMAGE018
的静态随机过程函数，
S1-4、建立静态随机过程函数在任意两个采样时间点间的协方差计算模型
在飞机地面颤振试验定义的时间范围内，得到的静态随机过程函数在任意两个采样时间点间的协方差计算模型如下：
Figure 855709DEST_PATH_IMAGE020
（4）
上式中，
Figure 610039DEST_PATH_IMAGE021
为方差，
Figure 347050DEST_PATH_IMAGE023
为第一采样时间点的静态随机过程函数，
Figure 526228DEST_PATH_IMAGE025
为第二采样时间点的静态随机过程函数，
Figure 956072DEST_PATH_IMAGE027
为第一采样时间点，
Figure 932118DEST_PATH_IMAGE029
为第二采样时间点，
Figure 738400DEST_PATH_IMAGE031
为相关函数，相关函数的形式为高斯函数或三次样条函数，
S1-5、建立非定常气动力快速计算模型，即Kriging代理模型在任意时刻
Figure 381871DEST_PATH_IMAGE033
的预测值计算模型
通过公式（4）构建形如式（3）的样本点数据列向量定义公式如下：
Figure DEST_PATH_IMAGE203
（5）
上式中，
Figure 717038DEST_PATH_IMAGE036
为样本点数据列向量，
Figure 852484DEST_PATH_IMAGE012
为基函数系数，
Figure 931298DEST_PATH_IMAGE038
为n个采样时间点的样本点数据，
Figure 429276DEST_PATH_IMAGE040
为基函数列向量，
Figure 384944DEST_PATH_IMAGE042
，
Figure 601162DEST_PATH_IMAGE044
为随机函数列向量，
Figure 483667DEST_PATH_IMAGE046
，
根据公式（4）中的相关函数
Figure 836151DEST_PATH_IMAGE047
定义
Figure 44278DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵
Figure 482213DEST_PATH_IMAGE051
，以及表示任意时间点与所有采样时间点相关性的相关向量
Figure 902830DEST_PATH_IMAGE168
，相关矩阵
Figure 781924DEST_PATH_IMAGE051
的计算公式如下：
Figure 629794DEST_PATH_IMAGE055
（6）
上式中，
Figure 555025DEST_PATH_IMAGE051
为
Figure 231863DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 558939DEST_PATH_IMAGE057
为
Figure 843290DEST_PATH_IMAGE049
个时间采样点中第i个采样时间点，
Figure 255817DEST_PATH_IMAGE059
为
Figure 549395DEST_PATH_IMAGE060
个时间采样点中第j个采样时间点，
Figure 730977DEST_PATH_IMAGE062
为第i个采样时间点和第j个采样时间点之间的相关函数，
Figure 186230DEST_PATH_IMAGE063
表示大小为
Figure DEST_PATH_IMAGE204
的矩阵，
相关向量
Figure 226998DEST_PATH_IMAGE168
的计算公式如下：
Figure 793108DEST_PATH_IMAGE170
（7）
上式中，
Figure 829198DEST_PATH_IMAGE121
为
Figure 986509DEST_PATH_IMAGE171
个时间采样点中任意一个样本点，
Figure 373628DEST_PATH_IMAGE173
，
Figure 477851DEST_PATH_IMAGE175
为相关向量，
Figure DEST_PATH_IMAGE205
为
Figure 86555DEST_PATH_IMAGE178
个时间采样点中任意一个样本点
Figure 883610DEST_PATH_IMAGE121
与第1个采样时间点的相关函数，
Figure 492446DEST_PATH_IMAGE180
为
Figure 72463DEST_PATH_IMAGE181
个时间采样点中任意一个样本点
Figure 83144DEST_PATH_IMAGE121
与第2个采样时间点的相关函数，
Figure 51100DEST_PATH_IMAGE183
为
Figure 412812DEST_PATH_IMAGE181
个时间采样点中任意一个样本点
Figure 858836DEST_PATH_IMAGE121
与第
Figure 989603DEST_PATH_IMAGE184
个采样时间点的相关函数，
采样时间点之外的时间节点处的矩阵函数
Figure 580991DEST_PATH_IMAGE066
由已知采样点函数值的线性加权得到：
Figure 164419DEST_PATH_IMAGE068
（8）
上式中，
Figure DEST_PATH_IMAGE206
为加权行向量，
Figure 414135DEST_PATH_IMAGE066
为采样时间点外的时间节点函数，
Figure 805933DEST_PATH_IMAGE070
为样本点数据列向量，
为满足Kriging代理模型无偏估计
Figure 115691DEST_PATH_IMAGE072
的要求，根据公式（5）与公式（8）可得基函数与基函数列向量关系的计算公式：
Figure DEST_PATH_IMAGE207
（9）
上式中，
Figure 451995DEST_PATH_IMAGE075
为基函数列向量，
Figure 505401DEST_PATH_IMAGE076
为基函数行向量，
Figure 345181DEST_PATH_IMAGE077
为加权行向量，
Figure 825841DEST_PATH_IMAGE012
为基函数系数，
Figure 777004DEST_PATH_IMAGE079
为数学期望，
Figure DEST_PATH_IMAGE208
为采样时间点外的时间节点函数，
Figure 634102DEST_PATH_IMAGE080
为样本点数据列向量，
根据公式（3）、公式（8）与公式（9）得到Kriging代理模型的预测值误差计算公式如下：
Figure 531651DEST_PATH_IMAGE082
（10）
上式中，
Figure DEST_PATH_IMAGE209
为采样时间点外的时间节点函数值，
Figure DEST_PATH_IMAGE210
为
Figure DEST_PATH_IMAGE211
通过Kriging代理模型得到的对应随机函数，
Figure DEST_PATH_IMAGE212
为任一时变气动力影响系数矩阵中的元素，
Figure 714370DEST_PATH_IMAGE088
为随机函数列向量，
Figure DEST_PATH_IMAGE213
是均值为0、方差为
Figure DEST_PATH_IMAGE214
的静态随机过程函数，
Figure 149900DEST_PATH_IMAGE091
为加权行向量，
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，通过公式（10）得到Kriging代理模型的均方误差计算公式如下：
Figure 545109DEST_PATH_IMAGE186
（11）
上式中，
Figure DEST_PATH_IMAGE215
为Kriging代理模型的预均方误差，
Figure 297164DEST_PATH_IMAGE209
为采样时间点外的时间节点函数，
Figure 385206DEST_PATH_IMAGE210
为
Figure DEST_PATH_IMAGE216
通过Kriging代理模型得到的对应随机函数，
Figure DEST_PATH_IMAGE217
为任一时变气动力影响系数矩阵中的元素，
Figure 183398DEST_PATH_IMAGE191
为加权行向量，
Figure 116719DEST_PATH_IMAGE192
为数学期望，
Figure 785598DEST_PATH_IMAGE103
为方差，
Figure DEST_PATH_IMAGE218
为加权列向量，
Figure 434754DEST_PATH_IMAGE104
为相关向量，
Figure 923504DEST_PATH_IMAGE051
为
Figure 457253DEST_PATH_IMAGE105
个时间采样点之间相关性的相关矩阵，
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，即为在
Figure DEST_PATH_IMAGE219
的前提下，寻找满足均方误差最小的加权列向量
Figure 980638DEST_PATH_IMAGE095
，为此引入拉格朗日乘数
Figure 348166DEST_PATH_IMAGE097
构建拉格朗日函数
Figure DEST_PATH_IMAGE220
如下：
Figure DEST_PATH_IMAGE221
（12）
上式中，
Figure 589791DEST_PATH_IMAGE220
为拉格朗日函数，
Figure 786286DEST_PATH_IMAGE103
为方差，
Figure 164178DEST_PATH_IMAGE095
为加权列向量，
Figure 764924DEST_PATH_IMAGE091
为加权行向量，
Figure 493845DEST_PATH_IMAGE104
为相关向量，
Figure 838239DEST_PATH_IMAGE051
为
Figure 336216DEST_PATH_IMAGE105
个时间采样点之间相关性的相关矩阵，
Figure 842284DEST_PATH_IMAGE115
为已知的基函数，
Figure 465026DEST_PATH_IMAGE109
为拉格朗日乘数，
Figure 347532DEST_PATH_IMAGE111
为基函数行向量，
分别对公式（12）中加权列向量
Figure 700016DEST_PATH_IMAGE095
和拉格朗日乘数
Figure 826585DEST_PATH_IMAGE109
求偏导数，得到满足要求的加权列向量
Figure 264520DEST_PATH_IMAGE112
计算公式为：
Figure DEST_PATH_IMAGE222
（13）
上式中，
Figure 216295DEST_PATH_IMAGE119
为
Figure 423285DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 536735DEST_PATH_IMAGE104
为相关向量，
Figure 461966DEST_PATH_IMAGE115
为已知的基函数，
Figure 92798DEST_PATH_IMAGE075
为基函数列向量，
Figure 419874DEST_PATH_IMAGE111
为基函数行向量，
将公式（13）带入公式（8）得到非定常气动力快速计算模型，即Kriging代理模型在任意时刻
Figure 704225DEST_PATH_IMAGE033
的预测值计算模型如下：
Figure DEST_PATH_IMAGE223
（14）
上式中，
Figure 303703DEST_PATH_IMAGE208
为采样时间点外的时间节点函数，
Figure 331702DEST_PATH_IMAGE080
为样本点数据列向量，
Figure 513284DEST_PATH_IMAGE195
为
Figure 499695DEST_PATH_IMAGE049
个时间采样点之间相关性的相关矩阵，
Figure 399518DEST_PATH_IMAGE104
为相关向量，
Figure 965628DEST_PATH_IMAGE115
为已知的基函数，
Figure 142663DEST_PATH_IMAGE075
为基函数列向量，
Figure 768816DEST_PATH_IMAGE111
为基函数行向量；
S2、对步骤S1建立的非定常气动力快速计算模型进行仿真验证，得到仿真验证后的非定常气动力快速计算模型，若非定常气动力快速计算模型的准确性满足要求，得到仿真验证后的非定常气动力快速计算模型，若不满足要求，加密采样时间点后再次进行仿真验证，直至非定常气动力快速计算模型的准确性满足要求，具体包括以下步骤：
S2-1、在时变历程选取12个采样时间点，
S2-2、将12个采样时间点数据导入步骤S1，得到步骤S1-2计算的每个采样时间点的非定常气动力模型、步骤S1-5构建的非定常气动力快速计算模型，
S2-3、采用matlab.simulink开展变温地面颤振系统仿真，仿真对象为高超声速舵面结构，仿真对象温度分布如图3所示，时变历程为高超声速舵面结构由初始温度场升至稳态温度场的过程，初始温度场的温度为45℃，稳态温度场的温度为150℃，得到每个采样时间点的非定常气动力模型的定常颤振仿真以及整个时变历程的仿真、非定常气动力快速计算模型对每个采样时间点的定常颤振仿真以及整个时变历程的仿真，并分别与每个采样时间点的颤振频域数值计算结果对比，并分别向每个采样时间点的非定常气动力模型和非定常气动力快速计算模型中引入颤振临界速度追踪模块，颤振临界速度追踪模块根据仿真对象的响应变化趋势自动调整所述非定常气动力模型和所述非定常气动力快速计算模型中的风速，使仿真对象保持在颤振临界状态，实现对时变颤振速度边界的实时辨识与记录，
S2-4、结果如图4所示，可见三类结果一致性较好，证明了本方法的建模精度；
S3、采用步骤S2仿真验证后的非定常气动力快速计算模型进行飞机振动试验中时变气动载荷地面等效模拟，具体包括以下步骤：
S3-1、对飞机地面颤振试验颤振计算中的气动平面和结构间的插值点进行缩聚，通过优化方法选取最优插值点，最优插值点为：激振点和拾振点，将最优插值点作为试验件控制点，
S3-2、在试验件上方布设加热器，在试验件底部固定激振器，与激振器电性连接的功率放大器，在试验件上试验件控制点布置激光位移传感器和加速度传感器，功率放大器、激光位移传感器、加速度传感器与加载控制器电性连接，
S3-3、通过激光位移传感器和加速度传感器将获取的试验件控制点响应信号传入加载控制器，加载控制器将试验件控制点响应信号代入步骤S2仿真验证后的非定常气动力快速计算模型，非定常气动力快速计算模型计算试验件的等效气动力，
S3-4、将非定常气动力快速计算模型得出试验件的等效气动力并将等效气动力传入半实物仿真系统NI，半实物仿真系统NI的由通过NI VeriStand仿真配置软件对接Matlab.simulink的上位机和嵌入式控制器下位机构成，半实物仿真系统NI的上位机通过NI VeriStand仿真配置软件对接Matlab.simulink，将试验件的等效气动力转化为激振器控制信号，加载控制器通过功率放大器控制激振器将等效气动力对应的激振器控制信号实时加载在试验件上，加载完成后回到步骤S3-3，直至飞机振动试验结束。
实施例2
本实施例与实施例1的区别在于：
步骤S2-3中，初始温度场的温度为55℃，稳态温度场的温度为490℃。
实施例3
本实施例与实施例1的区别在于：
S2、通过实际飞机地面颤振试验对步骤S1建立的非定常气动力快速计算模型进行试验验证，具体包括以下步骤：
S2-1、如图6所示，对飞机地面颤振试验颤振计算中的气动平面和结构间的插值点进行缩聚，通过优化方法选取最优插值点，最优插值点为：激振点和拾振点，将最优插值点作为试验件控制点，再通过步骤S1建立适用于试验件控制点的非定常气动力快速计算模型；
S2-2、在试验件上方布设加热器，在试验件底部固定激振器，与激振器电性连接的功率放大器，在试验件上试验件控制点布置激光位移传感器和加速度传感器，功率放大器、激光位移传感器、加速度传感器与加载控制器电性连接；
S2-3、通过激光位移传感器和加速度传感器将获取的试验件控制点响应信号传入加载控制器搭载的非定常气动力快速计算模型，非定常气动力快速计算模型得出试验件的等效气动力并将等效气动力传入半实物仿真系统NI，半实物仿真系统NI的由通过NIVeriStand仿真配置软件对接Matlab.simulink的上位机和嵌入式控制器下位机构成，半实物仿真系统NI的上位机通过NI VeriStand仿真配置软件对接Matlab.simulink，将试验件的等效气动力转化为激振器控制信号，加载控制器通过激振器控制信号控制激振器将等效气动力实时加载在试验件上，加载完成后回到步骤S2-2，直至飞机振动试验等效模拟过程结束，在此过程中，加载控制器记录飞机振动试验等效模拟过程中的颤振临界速度，
S2-4、开展实际飞机地面颤振试验，并采用加载控制器的颤振临界速度追踪模块记录试验件在飞机地面颤振试时变历程中的颤振临界速度，与步骤S2-3中加载控制器记录飞机振动试验等效模拟过程中的颤振临界速度进行比对，判定非定常气动力快速计算模型的准确性。

Claims (10)

1.一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，包括以下步骤：
S1、基于时间冻结假设，通过将多个离散时间点的非定常气动力快速计算模型进行Kriging代理模型拟合，得到适用于飞机振动试验整个时变历程的非定常气动力快速计算模型，具体包括以下步骤：
S1-1、建立试验件控制点承受非定常气动力的计算公式，
S1-2、建立离散时间点的试验件控制点应承受非定常气动力计算模型，
S1-3、结合Kriging代理模型建立对应随机函数
Kriging代理模型设定
Figure DEST_PATH_IMAGE001
为对应随机函数
Figure 377279DEST_PATH_IMAGE002
的均方误差最小无偏估计结果，
Figure DEST_PATH_IMAGE003
由 采样时间点数据建立的全局近似模型以及随机过程两部分组成，
Figure 904338DEST_PATH_IMAGE002
的表达公式如下：
Figure 326092DEST_PATH_IMAGE004
（3）
上式中，
Figure 388726DEST_PATH_IMAGE005
为对应随机函数，
Figure DEST_PATH_IMAGE006
为已知的基函数，
Figure 404086DEST_PATH_IMAGE007
为基函数系数，
Figure 577579DEST_PATH_IMAGE008
为设计空 间内由采样点求解得到的全局近似模型，
Figure DEST_PATH_IMAGE009
是均值为0、方差为
Figure 244052DEST_PATH_IMAGE010
的静态随机过程函数，
S1-4、建立静态随机过程函数在任意两个采样时间点间的协方差计算模型
在飞机地面颤振试验定义的时间范围内，得到的静态随机过程函数在任意两个采样时间点间的协方差计算模型如下：
Figure DEST_PATH_IMAGE011
（4）
上式中，
Figure 743167DEST_PATH_IMAGE012
为方差，
Figure DEST_PATH_IMAGE013
为第一采样时间点的静态随机过程函数，
Figure 245823DEST_PATH_IMAGE014
为第二采样时间点 的静态随机过程函数，
Figure DEST_PATH_IMAGE015
为第一采样时间点，
Figure 223006DEST_PATH_IMAGE016
为第二采样时间点，
Figure DEST_PATH_IMAGE017
为相关函数，相关 函数的形式为高斯函数或三次样条函数，
S1-5、建立非定常气动力快速计算模型，即Kriging代理模型在任意时刻
Figure 884932DEST_PATH_IMAGE018
的预测值计算 模型
通过公式（4）构建形如式（3）的样本点数据列向量定义公式如下：
Figure DEST_PATH_IMAGE019
（5）
上式中，
Figure 443696DEST_PATH_IMAGE020
为样本点数据列向量，
Figure 27124DEST_PATH_IMAGE007
为基函数系数，
Figure 683364DEST_PATH_IMAGE021
为n个采 样时间点的样本点数据，
Figure DEST_PATH_IMAGE022
为基函数列向量，
Figure 199796DEST_PATH_IMAGE023
，
Figure DEST_PATH_IMAGE024
为随机函数列向 量，
Figure 165347DEST_PATH_IMAGE025
，
根据公式（4）中的相关函数
Figure DEST_PATH_IMAGE026
定义
Figure 501651DEST_PATH_IMAGE027
个时间采样点之间相关性的相关矩阵
Figure DEST_PATH_IMAGE028
，以及 表示任意时间点与所有采样时间点相关性的相关向量
Figure 961582DEST_PATH_IMAGE029
，相关矩阵
Figure 66941DEST_PATH_IMAGE028
的计算公式如下：
Figure DEST_PATH_IMAGE030
（6）
上式中，
Figure 78760DEST_PATH_IMAGE028
为
Figure 636780DEST_PATH_IMAGE031
个时间采样点之间相关性的相关矩阵，
Figure DEST_PATH_IMAGE032
为
Figure 385555DEST_PATH_IMAGE027
个时间采样点中第i个采样 时间点，
Figure 611000DEST_PATH_IMAGE033
为
Figure 528141DEST_PATH_IMAGE027
个时间采样点中第j个采样时间点，
Figure DEST_PATH_IMAGE034
为第i个采样时间点和第j个采样 时间点之间的相关函数，
Figure 714402DEST_PATH_IMAGE035
表示大小为
Figure 375191DEST_PATH_IMAGE036
的矩阵，
采样时间点之外的时间节点处的矩阵函数
Figure DEST_PATH_IMAGE037
由已知采样点函数值的线性加权得 到：
Figure 845355DEST_PATH_IMAGE038
（8）
上式中，
Figure DEST_PATH_IMAGE039
为加权行向量，
Figure 198976DEST_PATH_IMAGE037
为采样时间点外的时间节点函数，
Figure 731589DEST_PATH_IMAGE040
为样本点数 据列向量，
为满足Kriging代理模型无偏估计
Figure DEST_PATH_IMAGE041
的要求，根据公式（5）与公式（8）可得基函数 与基函数列向量关系的计算公式：
Figure 337014DEST_PATH_IMAGE042
（9）
上式中，
Figure 271472DEST_PATH_IMAGE022
为基函数列向量，
Figure DEST_PATH_IMAGE043
为基函数行向量，
Figure 530415DEST_PATH_IMAGE044
为加权行向量，
Figure 439071DEST_PATH_IMAGE007
为基函数系数，
Figure DEST_PATH_IMAGE045
为数学期望，
Figure 707242DEST_PATH_IMAGE046
为采样时间点外的时间节点函数，
Figure 496206DEST_PATH_IMAGE047
为样本点数据列向量，
根据公式（3）、公式（8）与公式（9）得到Kriging代理模型的预测值误差计算公式如下：
Figure DEST_PATH_IMAGE048
（10）
上式中，
Figure 66996DEST_PATH_IMAGE049
为采样时间点外的时间节点函数值，
Figure DEST_PATH_IMAGE050
为
Figure 839780DEST_PATH_IMAGE051
通过Kriging代理模型 得到的对应随机函数，
Figure DEST_PATH_IMAGE052
为任一时变气动力影响系数矩阵中的元素，
Figure 770696DEST_PATH_IMAGE053
为随机函数列向 量，
Figure DEST_PATH_IMAGE054
是均值为0、方差为
Figure 679746DEST_PATH_IMAGE012
的静态随机过程函数，
Figure 546071DEST_PATH_IMAGE055
为加权行向量，
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，即为在
Figure 150358DEST_PATH_IMAGE056
的前提下，寻找满足均方误差最小的加权列向量
Figure DEST_PATH_IMAGE057
，为此引入拉格朗日乘数
Figure 25911DEST_PATH_IMAGE058
构建拉格朗日 函数
Figure DEST_PATH_IMAGE059
如下：
Figure 415566DEST_PATH_IMAGE060
（12）
上式中，
Figure 452792DEST_PATH_IMAGE059
为拉格朗日函数，
Figure 669010DEST_PATH_IMAGE010
为方差，
Figure 817094DEST_PATH_IMAGE057
为加权列向量，
Figure 310524DEST_PATH_IMAGE055
为加权行向量，
Figure DEST_PATH_IMAGE061
为相 关向量，
Figure 784230DEST_PATH_IMAGE062
为
Figure 487744DEST_PATH_IMAGE031
个时间采样点之间相关性的相关矩阵，
Figure 173940DEST_PATH_IMAGE063
为已知的基函数，
Figure 771144DEST_PATH_IMAGE064
为拉格朗日乘 数，
Figure DEST_PATH_IMAGE065
为基函数行向量，
分别对公式（12）中加权列向量
Figure 415752DEST_PATH_IMAGE066
和拉格朗日乘数
Figure 606562DEST_PATH_IMAGE064
求偏导数，得到满足要求的加权列 向量
Figure DEST_PATH_IMAGE067
计算公式为：
Figure 237394DEST_PATH_IMAGE068
（13）
上式中，
Figure 830050DEST_PATH_IMAGE028
为
Figure 114401DEST_PATH_IMAGE031
个时间采样点之间相关性的相关矩阵，
Figure 415675DEST_PATH_IMAGE061
为相关向量，
Figure 709254DEST_PATH_IMAGE063
为已知的基函 数，
Figure 156415DEST_PATH_IMAGE022
为基函数列向量，
Figure 877247DEST_PATH_IMAGE065
为基函数行向量，
将公式（13）带入公式（8）得到非定常气动力快速计算模型，即Kriging代理模型在任意 时刻
Figure 511490DEST_PATH_IMAGE069
的预测值计算模型如下：
Figure 218546DEST_PATH_IMAGE070
（14）
上式中，
Figure 785794DEST_PATH_IMAGE037
为采样时间点外的时间节点函数，
Figure DEST_PATH_IMAGE071
为样本点数据列向量，
Figure 943106DEST_PATH_IMAGE028
为
Figure 454859DEST_PATH_IMAGE027
个 时间采样点之间相关性的相关矩阵，
Figure 90239DEST_PATH_IMAGE061
为相关向量，
Figure 511993DEST_PATH_IMAGE063
为已知的基函数，
Figure 574627DEST_PATH_IMAGE022
为基函数列向 量，
Figure 324409DEST_PATH_IMAGE065
为基函数行向量；
S2、对步骤S1建立的非定常气动力快速计算模型进行仿真验证，得到仿真验证后的非定常气动力快速计算模型，若非定常气动力快速计算模型的准确性满足要求，得到仿真验证后的非定常气动力快速计算模型，若不满足要求，加密采样时间点后再次进行仿真验证，直至非定常气动力快速计算模型的准确性满足要求；
S3、采用步骤S2仿真验证后的非定常气动力快速计算模型进行飞机振动试验中时变气动载荷地面等效模拟。
2.如权利要求1所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S1-1具体包括以下内容：
根据飞机地面颤振试验中试验件的结构振动响应计算试验件控制点承受的非定常气动力，试验件控制点承受非定常气动力的计算公式如下：
Figure 497901DEST_PATH_IMAGE072
（1）
上式中，
Figure DEST_PATH_IMAGE073
为试验件控制点应承受的非定常气动力，
Figure 39741DEST_PATH_IMAGE074
为来流动压，
Figure 7697DEST_PATH_IMAGE075
为降阶后的气 动力影响系数矩阵，
Figure 995507DEST_PATH_IMAGE076
是通过与马赫数相匹配的非定常气动力理论计算后并依据结构动 力学特性进行降阶处理得到的，
Figure DEST_PATH_IMAGE077
表示控制点位移，
Figure 238269DEST_PATH_IMAGE078
为减缩频率。
3.如权利要求2所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征 在于，所述减缩频率
Figure 369036DEST_PATH_IMAGE078
的计算公式如下：
Figure 648839DEST_PATH_IMAGE079
上式中，
Figure 497846DEST_PATH_IMAGE080
为减缩频率，
Figure DEST_PATH_IMAGE081
为简谐运动的圆频率，
Figure 544300DEST_PATH_IMAGE082
为参考长度，对于飞机的翼面类结构， 参考长度取根弦长度的一半，
Figure 654207DEST_PATH_IMAGE083
为来流速度。
4.如权利要求2所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S1-2具体包括以下内容：
由于实际飞机地面颤振试验需要试验件控制点应承受非定常气动力的时域关系式，将公式（1）通过最小状态法有理函数近似延拓至拉氏域，再通过拉氏反变换得到离散时间点的试验件控制点应承受非定常气动力计算模型如下：
Figure DEST_PATH_IMAGE084
（2）
上式中，
Figure 229545DEST_PATH_IMAGE085
为与时间相关的试验件控制点应承受非定常气动力，
Figure 565848DEST_PATH_IMAGE086
为位移气动力影 响系数矩阵，
Figure 884834DEST_PATH_IMAGE087
为速度气动力影响系数矩阵，
Figure 865560DEST_PATH_IMAGE088
为加速度气动力影响系数矩阵，
Figure 611799DEST_PATH_IMAGE089
为第一 辅助气动力影响系数矩阵，
Figure DEST_PATH_IMAGE090
为第二辅助气动力影响系数矩阵，
Figure 700977DEST_PATH_IMAGE091
为第三辅助气动力影响 系数矩阵，
Figure DEST_PATH_IMAGE092
为控制点位移，
Figure 446823DEST_PATH_IMAGE093
为控制点速度，
Figure DEST_PATH_IMAGE094
为控制点加速度，
Figure 672268DEST_PATH_IMAGE095
为来流动压，
Figure DEST_PATH_IMAGE096
为来流 速度，
Figure 730354DEST_PATH_IMAGE082
为参考长度，对于飞机的翼面类结构，参考长度取根弦长度的一半。
5.如权利要求1所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S1-5中：
相关向量
Figure 775670DEST_PATH_IMAGE029
的计算公式如下：
Figure 436459DEST_PATH_IMAGE097
（7）
上式中，
Figure DEST_PATH_IMAGE098
为
Figure 906623DEST_PATH_IMAGE099
个时间采样点中任意一个样本点，
Figure DEST_PATH_IMAGE100
，
Figure 260244DEST_PATH_IMAGE101
为相关向量，
Figure 792857DEST_PATH_IMAGE102
为
Figure 132702DEST_PATH_IMAGE103
个时间采样点中任意一个样本点
Figure 67160DEST_PATH_IMAGE098
与第1个采样时间点的相关函数，
Figure DEST_PATH_IMAGE104
为
Figure 591683DEST_PATH_IMAGE105
个时间采样 点中任意一个样本点
Figure 237690DEST_PATH_IMAGE098
与第2个采样时间点的相关函数，
Figure DEST_PATH_IMAGE106
为
Figure 505860DEST_PATH_IMAGE105
个时间采样点中任意一个 样本点
Figure 294825DEST_PATH_IMAGE098
与第
Figure 865614DEST_PATH_IMAGE105
个采样时间点的相关函数。
6.如权利要求1所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S1-5还包括：
根据Kriging代理模型预测值满足无偏估计且均方误差最小的要求，通过公式（10）得到Kriging代理模型的均方误差计算公式如下：
Figure 107240DEST_PATH_IMAGE107
（11）
上式中，
Figure DEST_PATH_IMAGE108
为Kriging代理模型的预均方误差，
Figure 179101DEST_PATH_IMAGE046
为采样时间点外的时间节 点函数，
Figure 822572DEST_PATH_IMAGE050
为
Figure 547951DEST_PATH_IMAGE109
通过Kriging代理模型得到的对应随机函数，
Figure 542452DEST_PATH_IMAGE110
为任一时变气动力 影响系数矩阵中的元素，
Figure 886846DEST_PATH_IMAGE111
为加权行向量，
Figure DEST_PATH_IMAGE112
为数学期望，
Figure 791348DEST_PATH_IMAGE010
为方差，
Figure 828574DEST_PATH_IMAGE113
为加权列向量，
Figure 44792DEST_PATH_IMAGE061
为 相关向量，
Figure 192876DEST_PATH_IMAGE062
为
Figure 434108DEST_PATH_IMAGE114
个时间采样点之间相关性的相关矩阵。
7.如权利要求1所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S2包括以下步骤：
S2-1、在时变历程选取数个采样时间点；
S2-2、将数个采样时间点数据导入步骤S1，得到步骤S1-2计算的每个采样时间点的非定常气动力模型、步骤S1-5构建的非定常气动力快速计算模型；
S2-3、采用matlab.simulink开展非定常气动力模型和非定常气动力快速计算模型变温地面颤振系统仿真，仿真对象为高超声速舵面结构，时变历程为高超声速舵面结构由初始温度场升至稳态温度场的过程，得到每个采样时间点的非定常气动力模型的定常颤振仿真以及整个时变历程的仿真、非定常气动力快速计算模型对每个采样时间点的定常颤振仿真以及整个时变历程的仿真，并分别与每个采样时间点的颤振频域数值计算结果对比；
S2-4、步骤S2-3的对比结果为不一致时，加密采样时间点后，回到步骤S2-2，直至步骤S2-3的对比结果为一致。
8.如权利要求7所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述初始温度场的温度范围为45~55℃，稳态温度场的温度范围为150~490℃。
9.如权利要求7所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S2-3还包括：
分别向每个采样时间点的非定常气动力模型和非定常气动力快速计算模型中引入颤振临界速度追踪模块，颤振临界速度追踪模块根据仿真对象的响应变化趋势自动调整所述非定常气动力模型和所述非定常气动力快速计算模型中的风速，使仿真对象保持在颤振临界状态，实现对时变颤振速度边界的实时辨识与记录。
10.如权利要求1所述的一种飞机振动试验中时变气动载荷地面等效模拟方法，其特征在于，所述步骤S3具体包括以下步骤：
S3-1、对飞机地面颤振试验颤振计算中的气动平面和结构间的插值点进行缩聚，通过优化方法选取最优插值点，最优插值点为：激振点和拾振点，将最优插值点作为试验件控制点；
S3-2、在试验件上方布设加热器，在试验件底部固定激振器，与激振器电性连接的功率放大器，在试验件上试验件控制点布置激光位移传感器和加速度传感器，功率放大器、激光位移传感器、加速度传感器与加载控制器电性连接；
S3-3、通过激光位移传感器和加速度传感器将获取的试验件控制点响应信号传入加载控制器，加载控制器将试验件控制点响应信号代入步骤S2仿真验证后的非定常气动力快速计算模型，非定常气动力快速计算模型计算试验件的等效气动力；
S3-4、加载控制器通过功率放大器控制激振器将等效气动力实时加载在试验件上，加载完成后回到步骤S3-3，直至飞机振动试验结束。
CN202210759016.5A 2022-06-30 2022-06-30 一种飞机振动试验中时变气动载荷地面等效模拟方法 Active CN114818550B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210759016.5A CN114818550B (zh) 2022-06-30 2022-06-30 一种飞机振动试验中时变气动载荷地面等效模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210759016.5A CN114818550B (zh) 2022-06-30 2022-06-30 一种飞机振动试验中时变气动载荷地面等效模拟方法

Publications (2)

Publication Number Publication Date
CN114818550A CN114818550A (zh) 2022-07-29
CN114818550B true CN114818550B (zh) 2022-09-13

Family

ID=82523243

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210759016.5A Active CN114818550B (zh) 2022-06-30 2022-06-30 一种飞机振动试验中时变气动载荷地面等效模拟方法

Country Status (1)

Country Link
CN (1) CN114818550B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116067602B (zh) * 2023-03-06 2023-07-14 西安航天动力研究所 一种随机振动信号生成方法、装置及电子设备

Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8393206B1 (en) * 2010-02-09 2013-03-12 Ping-Chih Chen Dry wind tunnel system
CN103310060A (zh) * 2013-06-19 2013-09-18 西北工业大学 一种跨音速极限环颤振分析方法
CN104899365A (zh) * 2015-05-27 2015-09-09 南京航空航天大学 一种可减小气动干扰不利影响的直升机气动布局优化方法
CN105631125A (zh) * 2015-12-28 2016-06-01 北京理工大学 基于降阶模型的气动-热-结构耦合分析方法
CN105843978A (zh) * 2014-12-18 2016-08-10 达索系统西姆利亚公司 基于数据的交互式3d体验
CN107341297A (zh) * 2017-06-21 2017-11-10 东南大学 一种基于kl展开的分布随机动载荷识别方法
CN107505842A (zh) * 2017-09-04 2017-12-22 重庆邮电大学 一种数控机床广义空间切削稳定性预测与优化方法
CN108416141A (zh) * 2017-08-31 2018-08-17 北京理工大学 一种线性时变结构模态振型辨识方法
CN109459206A (zh) * 2018-12-17 2019-03-12 西北工业大学 地面试验非定常气动力加载方法
CN109738144A (zh) * 2018-10-31 2019-05-10 中国飞机强度研究所 一种突风响应地面模拟试验方法
CN111767685A (zh) * 2020-07-08 2020-10-13 中国飞机强度研究所 一种热颤振特性试验方法
CN113392545A (zh) * 2021-05-31 2021-09-14 东南大学 一种考虑空间相关性的时域分布随机动载荷识别方法
CN113685314A (zh) * 2021-08-24 2021-11-23 浙江大学 变桨控制方法、系统和可读存储介质
CN114282310A (zh) * 2021-12-31 2022-04-05 北京航空航天大学 一种基于自适应加点代理模型的气弹结构耦合寻优方法

Patent Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8393206B1 (en) * 2010-02-09 2013-03-12 Ping-Chih Chen Dry wind tunnel system
CN103310060A (zh) * 2013-06-19 2013-09-18 西北工业大学 一种跨音速极限环颤振分析方法
CN105843978A (zh) * 2014-12-18 2016-08-10 达索系统西姆利亚公司 基于数据的交互式3d体验
CN104899365A (zh) * 2015-05-27 2015-09-09 南京航空航天大学 一种可减小气动干扰不利影响的直升机气动布局优化方法
CN105631125A (zh) * 2015-12-28 2016-06-01 北京理工大学 基于降阶模型的气动-热-结构耦合分析方法
CN107341297A (zh) * 2017-06-21 2017-11-10 东南大学 一种基于kl展开的分布随机动载荷识别方法
CN108416141A (zh) * 2017-08-31 2018-08-17 北京理工大学 一种线性时变结构模态振型辨识方法
CN107505842A (zh) * 2017-09-04 2017-12-22 重庆邮电大学 一种数控机床广义空间切削稳定性预测与优化方法
CN109738144A (zh) * 2018-10-31 2019-05-10 中国飞机强度研究所 一种突风响应地面模拟试验方法
CN109459206A (zh) * 2018-12-17 2019-03-12 西北工业大学 地面试验非定常气动力加载方法
CN111767685A (zh) * 2020-07-08 2020-10-13 中国飞机强度研究所 一种热颤振特性试验方法
CN113392545A (zh) * 2021-05-31 2021-09-14 东南大学 一种考虑空间相关性的时域分布随机动载荷识别方法
CN113685314A (zh) * 2021-08-24 2021-11-23 浙江大学 变桨控制方法、系统和可读存储介质
CN114282310A (zh) * 2021-12-31 2022-04-05 北京航空航天大学 一种基于自适应加点代理模型的气弹结构耦合寻优方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A comprehensive preference-based optimization framework with application to high-lift aerodynamic design;Carrese,R;《ENGINEERING OPTIMIZATION》;20121231;全文 *
Ground Flutter Simulation Test Based on Reduced Order Modeling of Aerodynamics by CFD/CSD Coupling Method;Binwen Wang 等;《International journal of applied mechanics》;20190221;全文 *
变厚度翼型的非定常气动力模型降阶及影响分析;李文成等;《应用力学学报》;20180319(第02期);全文 *
水平轴风力机翼型非定常气动特性研究;刘鹏寅;《中国博士学位论文全文库 工程科技Ⅱ辑》;20200315;全文 *
航空发动机叶盘系统的多场耦合振动特性及优化设计研究;杨文军;《中国博士学位论文全文库 》;20200615;全文 *
高超声速飞行器气动—热—结构建模及模型降阶研究;陈鑫;《中国博士学位论文全文库 工程科技Ⅱ辑》;20160715;全文 *

Also Published As

Publication number Publication date
CN114818550A (zh) 2022-07-29

Similar Documents

Publication Publication Date Title
Liu et al. Robust modeling method for thermal error of CNC machine tools based on ridge regression algorithm
Houtzager et al. Recursive predictor-based subspace identification with application to the real-time closed-loop tracking of flutter
Miao et al. Study on the effects of changes in temperature-sensitive points on thermal error compensation model for CNC machine tool
Hazra et al. Aerodynamic shape optimization using simultaneous pseudo-timestepping
CN108121856B (zh) 一种全飞行域飞行器动稳定性分析方法
CN112528411B (zh) 一种基于模态减缩的几何非线性结构噪声振动响应计算方法
CN111368466B (zh) 一种基于频响函数参数修正的机械振动预测方法
CN109902404A (zh) 不同阻尼形式的结构时程响应积分的统一递推计算方法
CN114818550B (zh) 一种飞机振动试验中时变气动载荷地面等效模拟方法
CN115587512A (zh) 基于ANSYS TwinBuilder的锂电池热电耦合数字孪生模型构建方法
CN114347029B (zh) 一种用于气动软体机器人快速模拟的模型降阶方法
Zhou et al. Data-driven real-time prediction of pouch cell temperature field under minimal sensing
Tóth et al. Order and structural dependence selection of LPV-ARX models revisited
CN112580855A (zh) 基于自适应变异pso-bp神经网络的电缆群稳态温升预测方法
Hazra An efficient method for aerodynamic shape optimization
CN111241728A (zh) 一种欧拉方程的间断伽辽金有限元数值求解方法
Zhu et al. Two degree-of-freedom µ synthesis control with Kalman filter for flight environment simulation volume with sensors uncertainty
CN114970402A (zh) 飞机地面振动试验中气动载荷模拟计算系统及方法
CN113221064B (zh) 模型处理方法、装置、计算机设备和存储介质
CN111639463B (zh) 一种基于XGBoost算法的电力系统扰动后频率特征预测方法
CN104795812A (zh) 一种考虑定子阻尼效应并采用变量代换的慢同调分区法
CN114491730A (zh) 一种高速铁路路基结构动力安定分析迭代方法及装置
CN110309617B (zh) 传感器上温度、压力、振动三种载荷直接耦合的仿真方法
CN106773697A (zh) 一种时间维度拓展极限学习机模型的工业过程软测量建模方法
Winter et al. Application of Unsteady Aerodynamic Reduced-Order Modeling Techniques to a Complex Configuration

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