柔性飞翼飞机的模态分析
本算例给出了柔性翼飞机弯曲模态的计算。机翼的振动响应沿其跨的多点采集。数据用于识别系统的动态模型。从辨识出的模型中提取模态参数。将模态参数数据与传感器位置信息相结合，以可视化机翼的各种弯曲模态。这个例子需要信号处理工具箱™。
柔性翼飞机
在本例中，我们研究了从明尼苏达大学无人飞行器实验室(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 / fs;%样品时间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对象是在System Identification Toolbox™中打包时域数据的标准方法。输入信号被视为带限信号。
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 = Merge（数据{：}）
Data =包含10个实验的时域数据集。实验样本样本时间Exp1 600001 0.0005 Exp2 600001 0.0005 Exp3 600001 0.0005 Exp4 600001 0.0005 Exp5 600001 0.0005 Exp7 600001 0.0005 Exp8 600001 0.0005 Exp9 600001 0.0005 Exp10 600001 0.0005 output Unit(如果指定)y1 y2 Inputs Unit(如果指定)u1
模型识别
我们想要确定一个动态模型，其频率响应与实际飞机尽可能接近。动态模型将系统的输入和输出之间的数学关系封装为微分或差分方程。动态模型的例子有传递函数和状态空间模型。在系统识别工具箱中，动态模型由诸如idtf(传递函数),idpoly(适用于AR、ARMA模型)和中的难点(状态方程模型)。动态模型可以通过运行评估命令来创建，例如特遣部队和党卫军对时域或频域测量数据的命令。
对于本例，我们首先将实测的时域数据转换为频响数据，采用经验传递函数估计etfe命令。然后使用估计的FRF来识别飞机的振动响应的状态空间模型。可以直接使用时域数据进行动态模型识别。然而，FRF的数据形式允许将大型数据集压缩成更少的样本以及更容易调整相关频率范围的估计行为。FRFS被封装在一起idfrd对象。
估计每个数据实验的双输出/单输入频响函数(FRF)。使用没有窗口。使用24,000个频率点进行响应计算。
G = cell(1,10);N = 24000;为k = 1:10%使用ETFE命令将时域数据转换为频响Data_k = getexp(Data, k); / /数据G{k} = etfe(Data_k， []， N);％g {k}是@idfrd对象结束
将所有频响连接成一个20输出/ 1输入频响。
G = cat(1, G{:});%连接所有估计频响的输出g.outputname =“y”；% name输出'y1'， 'y2'，…,“日元”G.Intersample =.“提单”；
为了获得估计的频率响应的感觉，绘制输出1和输出15的振幅(任意选择)。放大到感兴趣的频率范围(4-35 Hz)。
选择= bodeoptions;％绘图选项opt.FreqUnits =“赫兹”；%表示频率，单位为Hzopt.PhaseVisible =“关闭”；％不显示阶段输出:[1 15];% pick输出1和15用于绘图clf bodeplot(G(OutputNum，:)， opt)％绘图频率响应网格xlim (35 [4])在
频响显示至少9个共振频率。为了进行分析，我们希望将重点放在6- 35hz的频率跨度上，这是飞机最关键的柔性弯曲模态所在的位置。因此将频响减小到这个频率区域。
f = g.frequency / 2 / pi;%提取频率矢量，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)' Input channels: 'u1' Status: Estimated model
Gs从而包含20个测量位置的频率响应测量。接下来，确定一个适合的状态空间模型Gs．子空间估计n4sid提供快速的估计。状态空间模型结构配置如下：
我们估计了一个18阶连续时间模型。对不同的顺序进行了一些试验，并检查了模型与频响函数的拟合结果。
该模型包含馈通术语（D矩阵是非零）。这是因为我们限制了我们的分析 35赫兹，而机翼的带宽明显大于35赫兹(响应显著在35赫兹)。
为了加快计算速度，我们抑制了参数协方差的计算。
FRF幅度在频率范围内变化显着变化。为了确保低幅度接收到相同的重点作为较高的幅度，我们应用自定义权重，其变化与第11响应的平方根相同。第11次输出的选择有点任意但有效，因为所有20个FRF都有类似的档案。
我们设置了评估选项n4sid使用n4sidOptions．
FRF =挤压（GS.Responsedata）;加权=平均值（1./sqrt（abs（frf）））。';n4opt = n4sidoptions（“EstimateCovariance”假的,......“WeightingFilter”，加权，......'输出'、眼睛(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搜索方法“显示”，“上”，......'输出'、眼睛(20));syso2 = syso2 (Gs, syso2, sopt);％估计状态 - 空间模型（这需要几分钟）适合= sys2.Report.Fit.FitPercent '
适合=1×2089.7225 89.585 89.7260 90.4986 88.5522 88.875 75.9215 83.5975 75.9215 83.1763 91.9215 89.7310 90.6844 89.8444 89.8444 89.8985 89.84467 87.8532 88.8532 88.03858 89898 89.8532 89.8532 88.8532 88.8989888
和之前一样，我们把最坏的和最好的结合起来。我们还可视化了模型频率响应的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 / pi;%数据频率(Hz)Fn = modalfit(sys2, f, 12);%固有频率(Hz)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.提取模态参数。对FRF的检测表示大约10个模态频率，大致围绕频率[5 7 10 15 17 23 27 30] Hz。我们可以使用该评估更多具体modalsd将模态参数的稳定性作为底层模型的顺序进行检查的命令。这个操作需要很长时间。因此，所得到的曲线直接作为图像插入。在下面的注释块内执行代码以重现该图。
frf = permute（gs.responsedata，[3 1 2]）;f = gs.frequency / 2 / pi;% {数字pf =模态（frf，f，fs，'maxmodes'，20，'fitmethod'，'lsrf'）;%}
检查地块和pfValues建议一个真实的自然频率的精炼列表:
Freq = [7.8 9.4 15.3 19.3 23.0 27.3 29.2 31.7];
我们使用弗里克向量作为从模型中选择最主要模式的指南SYS2.．这是用modalfit命令。
[Fn，DR，MS] = ModalFit（Sys2，F，长度（频率），“PhysFreq”、频率);
fn为固有频率(Hz)，博士相应的阻尼系数和女士归一化残差作为列向量(每一列对应一个固有频率)。在提取这些模态参数的过程中，只使用模型的稳定的欠阻尼极点。这女士列只包含具有正虚部的极点的数据。
模式形状可视化
为了可视化各种弯曲模态，我们需要上面提取的模态参数。此外，我们还需要关于测量点位置的信息。这些位置(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科技论坛。（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 (1);yAft = AccelPos(船尾,2);xFwd = AccelPos(录象,1);yFwd = AccelPos(前轮驱动,2);wn = fn (ModeNum) * 2 *π;%模式频率，单位为rad/秒t = 1 / fn（modenum）;%模态振荡周期Np = 2.25;%要模拟的周期数tmax = 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])视图(5,55)title(sprintf())% d’模式。FN = %s Hz'、ModeNum num2str (fn (ModeNum), 4)));％通过更新z-cocordate来动画。持有在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 = 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)结束结束
你也可以从以下列表中选择一个网站:
