柔性飞翼飞机的模态分析
本算例给出了柔性翼飞机弯曲模态的计算。机翼的振动响应沿其跨的多点采集。数据用于识别系统的动态模型。从辨识出的模型中提取模态参数。将模态参数数据与传感器位置信息相结合，以可视化机翼的各种弯曲模态。这个例子需要信号处理工具箱™。
柔性翼飞机
在本例中，我们研究了从明尼苏达大学无人飞行器实验室(Uninhabited Aerial Vehicle Laboratories, University of Minnesota[1])建造的小型柔性飞行翼飞机收集的数据。飞机的几何形状如下所示。
飞机机翼可以在负载发生大的变形。柔性模式频率比与更刚性的翼共同飞机低。这种灵活的设计降低了材料成本，增加灵活性和飞机的飞行范围。但是，如果不加以控制，灵活的方式可能导致灾难性的气动弹性不稳定（扑）。设计有效的控制规律为抑制这些不稳定性，需要机翼的各种弯曲模式中的准确测定。
实验装置
实验的目的是收集飞行器在不同位置对外界激励的振动响应。飞机在一个木框架上悬挂，在其重心上使用一个弹簧。弹簧足够灵活，因此弹簧-质量振荡的固有频率不会干扰飞机的基频。输入力通过靠近飞机中心的Unholtz-Dickie模型20电动激振器施加。
沿翼幅放置20个PCB-353B16加速度计，采集振动响应，如下图所示。
激振器输入命令被指定为该形式的恒定振幅啁啾输入 ．啁啾频率随时间线性变化，即: ．输入信号所覆盖的频率范围为3 - 35hz。数据由两个加速度计(一个是前缘加速度计，一个是后缘加速度计)收集x-location)每次。为此，进行了10次实验，采集了20个加速度计的全部响应。加速度计和力传感器的测量都在2000赫兹采样。
数据准备
数据由10组双输出/单输入信号表示，每组包含600K样本。这些数据可在MathWorks支持文件站点获得。金宝app看到免责声明．下载数据文件，并将数据加载到MATLAB®工作空间。
url ='//www.tatmou.com/金宝appsupportfiles/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:length(u)-1)' * t;图subplot(211) plot(t,y) ylabel(“输出(m / s ^ 2)”)传说(“前沿”，“后缘”) subplot(212) plot(t,u) ylabel('输入'）xlabel（的时间(秒)）
为了准备用于模型识别的数据，数据被打包成iddata对象。的iddata对象在系统辨识工具箱包装时域数据™的标准方式。作为带限输入信号进行处理。
ExpNames =字段名(MeasuredData);Data = cell(1,10);为k = 1:10 y = MeasuredData.(ExpNames{k}).y;k u = MeasuredData。(ExpNames {}) .u;Data{k} = iddata(y, u, t, t)“InterSample”，“提单”）;结束
将数据集对象合并为一个多实验数据对象。
Data =合并(数据{:})
含有10个实验数据=时域数据集。实验样品采样时间EXP1 600001 0.0005 EXP2 600001 0.0005 EXP3 600001 0.0005 EXP4 600001 0.0005 Exp5 600001 0.0005 Exp6 600001 0.0005 Exp7 600001 0.0005 Exp8 600001 0.0005 Exp9 600001 0.0005 Exp10 600001 0.0005输出单元（如果指定）Y1 Y2输入单元（如果指定）U1
模式识别
我们希望，以确定其频率响应，实际飞机尽可能地匹配的动态模型。动态模型封装输入和系统输出为差分或差分方程之间的数学关系。动态模型的例子是传输函数和状态空间模型。在系统辨识工具箱，动态模型是由物体，如封装idtf(传递函数),idpoly（AR为，ARMA模型）和中的难点(状态方程模型)。动态模型可以通过运行评估命令来创建，例如特遣部队和ssest对时域或频域测量数据的命令。
对于本例，我们首先将实测的时域数据转换为频响数据，采用经验传递函数估计ETFE命令。然后利用估计的频响函数来识别飞机振动响应的状态空间模型。可以直接使用时域数据进行动态模型识别。然而，FRF形式的数据允许将大数据集压缩成更少的样本，并更容易调整估计行为到相关的频率范围。frf封装为idfrd对象。
估计每个数据实验的双输出/单输入频响函数(FRF)。使用没有窗口。使用24,000个频率点进行响应计算。
G = cell(1,10);N = 24000;为k = 1:10%使用ETFE命令将时域数据转换为频响Data_k = getexp（数据，K）;ģ{K} = ETFE（Data_k，[]，N）;% G{k}是一个@idfrd对象结束
串连所有的FRF成单个20输出/单输入FRF。
G = cat(1, G{:});%连接所有估计频响的输出G.OutputName =“y”；% name输出'y1'， 'y2'，…,“日元”G.InterSample =“提单”；
为了获得估计的频率响应的感觉，绘制输出1和输出15的振幅(任意选择)。放大到感兴趣的频率范围(4-35 Hz)。
选择= bodeoptions;%绘图选项opt.FreqUnits =“赫兹”；％显示在赫兹频率opt.PhaseVisible =“关闭”；%不显示相位OutputNum = [1 15];％挑输出1和15，用于绘制clf bodeplot(G(OutputNum，:)， opt)图频率响应百分比网格xlim (35 [4])在
的FRF显示至少9的谐振频率。为了分析，我们希望把重点放在6-35赫兹频率范围在飞机谎言的最关键的灵活弯曲模式。因此减少了FRF该频率区域。
f = G.Frequency / 2 /π;%提取频率矢量，Hz (G存储频率，rad/s)Gs = fselect(G, f>6 & f<=32)% "fselect"选择要求范围内的频响(6.5 - 35hz)
Gs = IDFRD模型。包含20个输出和1个输入的频率响应数据。响应数据可在624个频率点获得，范围从37.96 rad/s到201.1 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提供快速的非迭代估计。状态空间模型结构配置如下:
我们估计的18阶连续时间模型。该命令是一些试验用的各种命令和检查所产生的模型到FRF的配合后找到。
该模型包含一个馈通项(D矩阵非零)。这是因为我们的分析局限于 35赫兹而机翼的带宽比显著更大（响应是显著在35赫兹）。
为了加快计算速度，我们抑制了参数协方差的计算。
频响幅度在频率范围内变化显著。为了确保低振幅和高振幅得到同样的重视，我们采用了自定义权重，其变化与第11个响应的平方根成反比。第11个输出的选择有些随意，但由于所有20个frf都有类似的配置文件，所以是可行的。
我们设置了评估选项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）测量频率响应。最贫穷和最适合的未来绘制。
[〜，伊明] =分钟（FIT）;[〜，IMAX] = MAX（FIT）;CLF bodeplot度（Gs（[伊明，IMAX]，:)，SYS1（[伊明，IMAX]，:)，优化）;XLIM（[6 32]）标题（“数据和模型之间最糟糕(上)和最好(下)吻合”网格)在传奇(“数据”，'模型'）
与模型的匹配sys1通过对模型参数的迭代非线性最小二乘精化，可以进一步改进。可以使用ssest命令。我们设置了评估选项ssest使用ssestOptions命令。这一次也估计了参数协方差。
ssOpt = ssestOptions (“EstimateCovariance”，真的，．..“WeightingFilter”, n4Opt。WeightingFilter,．..“SearchMethod”，“lm”，．..%使用Levenberg-Marquardt搜索方法“显示”，“上”，．..“OutputWeight”、眼睛(20));系统2 = ssest（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(Fit);[~， Imax] = max(Fit);clf h = bodeplot (Gs (Imin, Imax,:), sys2 ((Imin, Imax):),选择);showConfidence(h, 1) xlim([6.5 35]) title(“最糟糕(上)和最好(下)的数据与优化模型相吻合”网格)在传奇(“数据”，“模型(sys2)”）
该精炼状态空间模型sys2在7-20赫兹范围内吻合得很好。不确定性边界在大多数共振位置附近都很紧。我们估计了一个24阶模型，这意味着系统中最多有12种振动模式sys2．使用modalfit命令获取这些模式的以Hz为单位的固有频率。
f = Gs.Frequency / 2 /π;％的数据的频率（Hz）的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]Hz左右。我们可以用modalsd命令，在底层模型的顺序更改时检查模态参数的稳定性。这个操作需要很长时间。因此，得到的图作为图像直接插入。执行下面注释块中的代码来重新生成图。
频=排列(Gs。ResponseData,中[3 1 2]);f = Gs.Frequency / 2 /π;% {数字润扬悬索桥pf = modalsd (f, f, MaxModes, 20日“FitMethod”,“lsrf”);％}
检查地块和PF.Values建议一个真实的自然频率的精炼列表:
Freq = [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 Hz的模式。
AnimateOneMode(3, fn, ms, AccelPos);
结论
这个例子展示了一种基于参数模型识别的模态参数估计方法。利用24阶状态空间模型，在6-32 Hz频段提取了8种稳定的振荡模式。将模态信息与加速度计位置知识相结合，以可视化各种弯曲模态。其中一些模式如下图所示。
参考文献
Gupta, Abhineet, Peter J. Seiler, Brian P. Danowsky。柔性飞行翼飞机地面振动试验AIAA大气飞行力学大会，AIAA科技论坛。张仁(2016 - 1753)。
功能AnimateOneMode（modenum所，FN，ModeShapes，AccelPos）%动画一个模式。% ModeNum:模式的索引。重新排列传感器位置，以便绘制连续的，绕机身绘制非相交曲线。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（AFT，1）;yAft = AccelPos（AFT，2）;xFwd = AccelPos（FWD，1）;yFwd = AccelPos（FWD，2）;WN = FN（modenum所）* 2 * PI;％模式频率以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,“- - -”，X，Y，0 * Z0，“。”， 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])视图(5,55)title(sprintf())% d’模式。FN = %s Hz'、ModeNum num2str (fn (ModeNum), 4)));%通过更新z坐标进行动画。持有在Col2 = [0.87 0.5 0];h1 = plot3 (x, y, 0 * z0,“- - -”，X，Y，0 * Z0，“。”， 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 = real(z0)*Rot1 - imag(z0)*Rot2;zAft_t = real(zAft)*Rot1 - imag(zAft)*Rot2;zFwd_t = real(zFwd)*Rot1 - imag(zFwd)*Rot2;[[zAft_t, zFwd_t]'; / /修改zAft_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)结束结束
您还可以从以下列表中选择一个网站：
