柔性飞翼飞机的模态分析
本算例给出了柔性翼飞机弯曲模态的计算。机翼的振动响应沿其跨的多点采集。数据用于识别系统的动态模型。从辨识出的模型中提取模态参数。将模态参数数据与传感器位置信息相结合，以可视化机翼的各种弯曲模态。这个例子需要信号处理工具箱™。
柔性翼飞机
在本例中，我们研究了从明尼苏达大学无人飞行器实验室(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(“输入”)包含(的时间(秒)）
为了准备用于模型识别的数据，数据被打包成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 =合并(数据{:})
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封装为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 /π;%提取频率矢量，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阶连续时间模型。对不同的顺序进行了一些试验，并检查了模型与频响函数的拟合结果。
该模型包含一个馈通项(D矩阵非零)。这是因为我们的分析局限于 35赫兹，而机翼的带宽明显大于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)拟合优度测量。最不适合和最好适合的是接下来绘制的。
[~, 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));syso2 = syso2 (Gs, syso2, sopt);%估计状态空间模型(这需要几分钟)适合= 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);%固有频率(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提取模态参数。对频响的检查表明大约有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”);%}
检查地块和pfValues建议一个真实的自然频率的精炼列表:
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为固有频率(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科技论坛。张仁(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;%要模拟的周期数达峰时间= 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坐标进行动画。持有在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)结束结束
你也可以从以下列表中选择一个网站:
