模态分析的一个灵活的飞行翼飞机
这个例子展示了计算弯曲模式的柔性翼飞机。机翼的振动响应是收集在多个点沿跨度。数据被用来确定一个系统的动态模型。提取模态参数识别模型。模态参数数据结合传感器位置信息可视化的各种弯曲模式。这个例子需要信号处理工具箱™。
柔性翼飞机
在这个例子中,我们研究收集的数据从一个小灵活的飞翼飞机建造的无人飞行器实验室、明尼苏达大学[1]。飞机的几何图形如下所示。
飞机机翼载荷作用下可以发生大的变形。灵活的模式频率低于普通飞机的翅膀更加僵化。这种灵活的设计降低了材料成本,提高敏捷性和飞机的飞行范围。然而,如果不控制,灵活的模式可能会导致灾难性的气动弹性不稳定(颤振)。设计有效的抑制控制律这些不稳定需要准确确定机翼弯曲的各种模式。
实验装置
实验的目的是收集飞机在不同位置的振动响应来响应外部激励。飞机是暂停一个木制框架使用一个春天的重心。春天是足够灵活,这样的弹簧-质量振动固有频率不干扰飞机的基本频率。一个输入的力通过Unholtz-Dickie应用模型飞机的中心附近20电动振动器。
二十pcb - 353 b16转椅感应器放置沿着翼展收集第二图所示的振动响应。
瓶输入命令被指定为一个恒定的振幅唧唧喳喳的输入形式 。短促声波频率随时间线性变化,也就是说, 。输入信号的频率范围覆盖3-35赫兹。收集的数据是两个加速度计(领导和后缘加速计x一次位置)。因此10实验进行收集所有20个加速度计的响应。加速度计和测力传感器测量采样在2000赫兹。
数据准备
数据是由10套两个输出/一个输入信号,每个都包含600 k样本。的数据可用MathWorks支持文件的网站。万博1manbetx看到免责声明。下载数据文件并将数据装载到MATLAB®工作区。
url =“//www.tianjin-qmedu.com/万博1manbetxsupportfiles/ident/flexible-wing-GVT-data/FlexibleWingData.mat”;websave (“FlexibleWingData.mat”url);负载FlexibleWingData.matMeasuredData
的变量MeasuredData是一个结构域Exp1,Exp2、……Exp10。每个字段是一个结构字段y包含两个反应和u和相应的数据输入的力。第一个实验的数据。
fs = 2000;%数据采样频率Ts = 1 / f;%样品时间y = MeasuredData.Exp1.y;%输出数据(2列,一个用于每个加速度计)u = MeasuredData.Exp1.u;%输入力数据t =(0:长度(u) 1)“* Ts;图次要情节(211)情节(t、y) ylabel (“输出(m / s ^ 2)”)传说(“前沿”,“后缘”)次要情节(212)情节(t, u) ylabel (“输入”)包含(的时间(秒))
为了准备数据模型识别、数据打包成iddata对象。的iddata对象包装时域数据的标准方法在系统辨识工具箱™。输入信号被视为带宽受限。
ExpNames =字段名(MeasuredData);Data =单元(10);为k = 1:10 y = MeasuredData。(ExpNames {k}) .y;k u = MeasuredData。(ExpNames {}) .u;数据{k} = iddata (y、u, Ts,“InterSample”,“提单”);结束
合并成一个multi-experiment数据对象的数据集对象。
Data =合并(数据{:})
Data =时域数据集包含10个实验。实验样品样本Exp2 Exp1 600001 0.0005 600001 0.0005 Exp3 600001 0.0005 Exp4 600001 0.0005 Exp5 600001 0.0005 600001年Exp6 0.0005 Exp7 600001 0.0005 Exp8 0.0005 Exp10 600001 0.0005 600001 0.0005 Exp9 600001输出单位(如果指定了的话)y1 y2输入单位(如果指定了的话)u1
模式识别
我们想要确定一个动态模型,其频率响应匹配实际的飞机尽可能。一个动态模型封装了一个数学系统的输入和输出之间的关系作为一个微分或差分方程。动态模型传递函数和状态空间模型的例子。在系统辨识工具箱,动态模型封装等对象idtf(传递函数),idpoly(AR, ARMA模型)中的难点(状态方程模型)。动态模型可以通过运行评估创建命令等特遣部队和党卫军命令在时域或频域测量数据。
对于这个示例,我们第一次测量的时域数据转换成频率响应数据的经验传递函数估计使用etfe命令。用于识别估计结构的飞机状态空间模型的振动响应。可以直接使用时域数据进行动态模型识别。然而,结构形式的数据允许大型数据集的压缩到更少的样品以及更容易调整估计行为相关的频率范围。频被封装idfrd对象。
估计两个输出/一个输入频率响应函数(降维)为每个数据实验。使用没有窗口。使用24000个频率点进行响应计算。
细胞G = (10);N = 24000;为k = 1:10%时域数据转换成频使用ETFE命令Data_k = getexp(数据、k);G {k} = etfe (Data_k [] N);% G {k}是一个@idfrd对象结束
所有结构合并到一个20-output /一个输入误差。
G =猫(1 G {:});%连接输出的估计误差G。OutputName =“y”;%名称输出y1, y2,…,“日元”G。InterSample =“提单”;
为估计频率响应,感受情节振幅输出1和15(任意选择)。放大感兴趣的频率范围(4-35赫兹)。
选择= bodeoptions;%绘图选项opt.FreqUnits =“赫兹”;%显示赫兹的频率opt.PhaseVisible =“关闭”;%不展示阶段OutputNum = 15 [1];%选择输出1和15策划clf bodeplot (G (OutputNum:),选择)%绘制频率响应网格xlim (35 [4])在
频显示至少9共振频率。分析我们要关注6-35赫兹的频率跨度最关键的灵活弯曲模式的飞机谎言。因此减少了结构频率区域。
f = G.Frequency / 2 /π;%提取频率向量赫兹(G频率存储在rad / s)Gs = fselect (G、f > 6 & f < = 32)%”fselect”选择了结构在请求范围内(6.5 - 35赫兹)
Gs = IDFRD模型。包含频率响应数据输出20 (s)和1输入(s)。频率响应数据存在于624点,从37.96到201.1 (rad / s rad / s。样品时间:0.0005秒输出通道:‘y (1)’,‘y (2)’,‘y (3)’,‘y (4)’,‘y (5)’,‘y (6)’,‘y (7)’,‘y (8)’,‘y (9)’,‘y (10)’,‘y (11)’,‘y (12)’,‘y (13)’,‘y (14)’,‘y (15)’,‘y (16)’,‘y (17)’,‘y (18)’,‘y (19)’,‘y(20)输入通道:“u1”状态:估计模型
Gs因此包含20个测量的频率响应测量的位置。接下来,确定一个适合状态空间模型Gs。子空间估计量n4sid提供了一种快速noniterative估计。状态空间模型的结构配置如下:
我们估计一个18阶连续时间模型。订单后发现一些试验各种订单,检查结果符合模型的降维。
模型包含一个直通的术语(D矩阵非零)。这是因为我们限制我们的分析 35赫兹在机翼的带宽明显大于35赫兹)(反应是重要的。
加速计算,我们抑制参数协方差的计算。
频级整个频率范围的差异很大。为了确保低振幅更高的振幅得到平等的强调,我们应用一个自定义权重的平方根成反比11日回应。11日的选择输出有些武断,但因为所有工作20误差也有类似的资料。
我们建立的评估选项n4sid使用n4sidOptions。
频=挤压(Gs.ResponseData);权重=平均(1. /√(abs(降维)))。”;n4Opt = n4sidOptions (“EstimateCovariance”假的,…“WeightingFilter”权重,…“OutputWeight”、眼睛(20));sys1 = n4sid (Gs, 24岁,“t”0,“引线”,真的,n4Opt);适合= sys1.Report.Fit.FitPercent '
适合=1×2057.0200 57.9879 85.0160 86.3815 80.4879 80.4430 58.2216 45.2692 61.5057 76.7612 84.7305 86.2600 86.4266 85.0251 76.9208 82.1191 74.7982 79.6837 67.9078 76.7249
“适合”数字显示之间的配合比例数据(Gs)和模型(sys1)频率响应使用归一化均方根误差(NRMSE)拟合优度的措施。最贫穷和最适合绘制。
[~,Imin] = min(配合);[~,Imax] = max(配合);clf bodeplot (Gs (Imin, Imax,:), sys1 ((Imin, Imax):),选择);xlim(32[6])标题(“最好最差(上)和(下)适合数据与模型之间的网格)在传奇(“数据”,“模型”)
实现符合模型sys1可以提高更多的迭代非线性最小二乘优化模型的参数。这可以通过使用党卫军命令。我们建立的评估选项党卫军使用ssestOptions命令。这一次参数协方差估计。
ssOpt = ssestOptions (“EstimateCovariance”,真的,…“WeightingFilter”n4Opt.WeightingFilter,…“SearchMethod”,“lm”,…%使用Levenberg-Marquardt搜索方法“显示”,“上”,…“OutputWeight”、眼睛(20));sys2 = ss (Gs、sys1 ssOpt);%状态空间模型的估计(这需要几分钟)适合= sys2.Report.Fit.FitPercent '
适合=1×2089.7225 89.5185 89.7260 90.4986 88.5522 88.8727 81.3225 83.5975 75.9215 83.1763 91.1358 89.7310 90.6844 89.8444 89.6685 89.1467 87.8532 88.0385 84.2898 83.3578
像以前我们最糟糕和最适合的阴谋。我们也可视化1-sd置信区域模型的频率响应。
[~,Imin] = min(配合);[~,Imax] = max(配合);clf h = bodeplot (Gs (Imin, Imax,:), sys2 ((Imin, Imax):),选择);showConfidence (h, 1) xlim(35[6.5])标题(“最好最差(上)和(下)适合数据和改进模型之间的网格)在传奇(“数据”,“模型(sys2)”)
状态空间模型的精制sys27-20赫兹地区符合误差也很好。大多数共振位置周围的不确定性边界紧张。我们估计这意味着有24阶模型系统中最多12个振荡模式sys2。使用modalfit命令来获取固有频率在赫兹对这些模式。
f = Gs.Frequency / 2 /π;%数据频率(赫兹)fn = modalfit (sys2 f (12);%固有频率(赫兹)disp (fn)
7.7721 7.7953 9.3147 9.4009 9.4910 15.3463 19.3291 23.0219 27.4164 28.7256 31.7014 63.3034
中的值fn建议两个频率非常接近7.8赫兹和三个接近9.4赫兹。检查这些频率附近的频率响应显示在输出峰值位置移一点。这些差异可能会被更好地控制实验过程,执行直接时域识别输入带宽限制在小范围集中在这些频率,或安装一个振荡模式在这些频率频率响应。在本例中不探讨这些方案。
模态参数计算
我们现在可以使用模型sys2提取模态参数。降维的检查表明大约10模态频率,频率大约在[5 7 10 15 17日23日27日30]赫兹。我们可以做这个评估更具体的使用modalsd命令检查模态参数的稳定性作为底层模型的顺序发生了变化。这个操作需要很长时间。因此产生的直接插入情节作为一个形象。执行下面的代码注释块内繁殖。
频=排列(Gs。ResponseData,中[3 1 2]);f = Gs.Frequency / 2 /π;% {图润扬悬索桥pf = modalsd (f, f, MaxModes, 20日“FitMethod”,“lsrf”);%}
情节和检验pf值表明真正的自然频率的精制列表:
频率= (7.8 9.4 15.3 - 19.3 23.0 - 27.3 29.2 - 31.7);
我们使用的值频率向量作为挑选最具优势的指导模式的模型sys2。这是使用来完成的modalfit命令。
(fn,博士、女士)= modalfit (sys2, f,长度(频率),“PhysFreq”、频率);
fn是固有频率(赫兹),博士相应的阻尼系数女士标准化残差为列向量(一列每个固有频率)。这些模态参数提取过程中,只有稳定、欠阻尼模型的极点。的女士列包含的数据只有两极积极的虚部。
模式形状可视化
可视化的各种弯曲模式,我们需要上面的模态参数提取。另外我们需要测量的位置点信息。这些位置(x - y值)记录了每个矩阵的加速度计AccePos:
AccelPos = […%见图216.63 - 18.48;%的中心16.63 - 24.48;27.90 - 22.22;27.90 - 28.22;37.17 - 25.97;37.17 - 31.97;46.44 - 29.71;46.44 - 35.71;55.71 - 33.46;55.71 - 39.46;%最远的权利-16.63 - 18.48;%最近的中心-16.63 - 24.18;-27.90 - 22.22;-27.90 - 28.22;-37.17 - 25.97;-37.17 - 31.97;-46.44 - 29.71;-46.44 - 35.71;-55.71 - 33.46;-55.71 - 39.46);%最远的离开
中包含的模式形状矩阵女士每一列对应的形状为一个频率。动画模式通过叠加模式形状振幅在传感器坐标和不同模式的固有频率的振幅。动画显示了弯曲无阻尼。作为一个例子,考虑在15.3赫兹的模式。
AnimateOneMode (3 fn ms, AccelPos);
结论
这个例子显示了一个基于参数模型识别模态参数估计方法。使用24日阶状态空间模型,8稳定振荡模式提取6-32赫兹频率地区。模态信息结合加速度计职位的知识可视化的各种弯曲模式。这些模式所示下面的数据。
引用
[1]Gupta, Abhineet,彼得·j·西勒和布莱恩·Danowsky。“地面振动测试在一个灵活的飞翼飞机。”张仁大气飞行力学会议,张仁科技论坛。张仁(2016 - 1753)。
函数AnimateOneMode (ModeNum fn、ModeShapes AccelPos)%动画一个模式。% ModeNum:指数模式。%为策划调整传感器的位置,这样一个连续的,% non-intersecting曲线跟踪飞机的身体。PlotOrder = [19: 2:11 1:2:9 10: 2:2, 12:2:20, 19);Fwd = PlotOrder (1:10);尾= PlotOrder (20: 1:11);x = AccelPos (PlotOrder, 1);y = AccelPos (PlotOrder 2);船尾xAft = AccelPos (1);yAft = AccelPos(船尾,2);xFwd = AccelPos(录象,1);yFwd = AccelPos(前轮驱动,2);wn = fn (ModeNum) * 2 *π;%模式频率在rad /秒T = 1 / fn (ModeNum);%的模态振动Np = 2.25;%的时间来模拟达峰时间= Np * T;%模拟Np时期元= 100;%的时间步数的动画t = linspace(0,达峰时间,Nt);ThisMode = ModeShapes (:, ModeNum) / max (abs (ModeShapes (:, ModeNum)));%规范化的阴谋z0 = ThisMode (PlotOrder);zFwd = ThisMode (Fwd);zAft = ThisMode(尾);clf col1 = (1 1 1) * 0.5;xx =重塑([[xAft xFwd] ';南(2,10)],[2]20);yy =重塑([[yAft yFwd] ';南(2,10)],[2]20);plot3 (x, y, 0 * z0,“- - -”0 * z0, x, y,,“。”xx (:), yy(:), 0(40岁,1),“——”,…“颜色”col1,“线宽”,1“MarkerSize”10“MarkerEdgeColor”,col1);包含(“x”)ylabel (“y”)zlabel (“振幅”)ht = max (abs (z0));轴([-100 100 10 40 ht ht)视图(55)标题(sprintf (% d’模式。FN = % s赫兹'、ModeNum num2str (fn (ModeNum), 4)));%更新z坐标的动画。持有在col2 = (0.87 - 0.5 0);h1 = plot3 (x, y, 0 * z0,“- - -”0 * z0, x, y,,“。”xx (:), yy(:), 0(40岁,1),“——”,…“颜色”、col2“线宽”,1“MarkerSize”10“MarkerEdgeColor”、col2);h2 = fill3 (x, y, 0 * z0 col2,“FaceAlpha”,0.6);持有从为k = 1: Nt Rot1 = cos (wn * t (k));Rot2 =罪(wn * t (k));z_t =实际(z0) * Rot1 -图像放大(z0) * Rot2;zAft_t =实际(zAft) * Rot1 -图像放大(zAft) * Rot2;zFwd_t =实际(zFwd) * Rot1 -图像放大(zFwd) * Rot2;zz =重塑([[zAft_t zFwd_t] ';南(2,10)],[2]20);集(h1 (1),“ZData”z_t)组(h1 (2),“ZData”z_t)组(h1 (3),“ZData”zz (:)) h2.Vertices (:, 3) = z_t;暂停(0.1)结束结束
