模拟梁的固有频率和阻尼

对于单输入/锤激发

一系列的锤击激发了3DOF系统，传感器记录产生的位移。系统是按比例阻尼的，这样阻尼矩阵是质量矩阵和刚度矩阵的线性组合。

导入两组测量数据，包括激励信号、响应信号、时间信号和地面真实频率响应函数。第一组响应信号，Y1，测量第一个质量块和第二个质量块的位移，Y2，测量第二质量。每个激励信号由十个锤击串联而成，每个响应信号包含相应的位移。每个撞击信号的持续时间为2.53秒。激励和响应信号中都存在加性噪声。可视化第一次测量的第一次激励和响应通道。

[t, fs, X1, X2, Y1, Y2, f0, H0) = helperImportModalData ();X0 = X1 (: 1);Y0 = y₁(:1);helperPlotModalAnalysisExample ([t ' X0 Y0]);

根据动态柔度计算并绘制第一激励和响应通道的FRF，动态柔度是力的位移测量[1].默认情况下，FRF是通过对加窗段的频谱进行平均来计算的。由于每个锤子激励在下一次激励之前基本上会衰减，因此可以使用矩形窗口。将传感器指定为位移。

winLen = 2.5275 * fs;%样本中的窗口长度winLen modalfrf (X0, Y0, fs,“传感器”，“说”）

使用默认值估计的FRF“标题”估计器，在被测频带中包含三个突出的峰值，对应于三种灵活的振动模式。相干性接近这些峰值，在反共振区域低，响应测量的信噪比低。接近1的相干表示高质量的估计。的“标题”当噪声只存在于输出测量时，估计是最优的，而“氢气”当仅在输入[2]上存在加性噪声时，估计器是最优的。计算“标题”和“氢气”此FRF的估计值。

(频响1,f1) = modalfrf (X0, Y0, fs, winLen“传感器”，“说”);%计算FRF（H1）(误差2,f2) = modalfrf (X0, Y0, fs, winLen“传感器”，“说”，“估计员”，“氢气”);

当存在显著的测量噪声或激励较差时，参数化方法可以提供从数据中准确提取频响的额外选项。“子空间”方法首先对数据[3]拟合一个状态空间模型，然后计算它的频响函数。状态空间模型的顺序(等于极点的数量)和是否存在馈通可以指定来配置状态空间估计。

(频响3 f3) = modalfrf (X0, Y0, fs, winLen“传感器”，“说”，“估计员”，“子空间”，“引线”,真正的);

这里，FRF3是通过拟合一个包含馈通项的状态空间模型来估计的，该状态空间模型的阶数为1:10。比较估计的频响“标题”，“氢气”和“子空间”方法的理论频响。

helperPlotModalAnalysisExample（f1、FRF1、f2、FRF2、f3、FRF3、f0、H0）；

在响应峰附近，估计器的性能相当“氢气”估计者高估了反共振的响应。相干性不受估计器选择的影响。

接下来，使用峰值拾取算法估计每个模式的固有频率。峰值拾取算法是识别FRF中峰值的简单而快速的程序。它是一种局部方法，因为每个估计都是从单个频率响应函数生成的。它也是一个单自由度（SDOF）方法，因为每个模式的峰值是独立考虑的。因此，为每个FRF生成一组模态参数。根据之前的曲线图，指定200到1600 Hz的频率范围，其中包含三个峰值。

fn = modalfit(频响1 f1 fs 3“FitMethod”，“页”，“频率范围”1600年[200])

Fn = 1.0e+03 * 0.3727 0.8525 1.3707

固有频率大约是373、853和1371赫兹。绘制一个重建频响，并将其与测量数据进行比较modalfit. 使用从频响函数矩阵估计的模态参数重构FRF，频1呼叫modalfit同样没有输出参数，以生成包含重构的FRF的图。

modalfit（FRF1、f1、fs、3、，“FitMethod”，“页”，“频率范围”1600年[200])

重建的频响与第一激励和响应通道的实测频响基本一致。在下一节中，将考虑另外两个激励位置。

粗纱锤激

使用默认值计算并绘制所有三个传感器响应的频响函数“标题”估计量。指定测量类型为“流动输入”因为我们有粗纱锤激励。

winLen modalfrf (X1, Y1, fs,“传感器”，“说”，“测量”，“流动输入”）

在前一节中，从单个频响计算了一组模态参数。现在，使用最小二乘复指数(LSCE)算法估计模态参数。lce和LSRF算法通过同时分析多个响应信号生成一组模态参数。这些都是全局的多自由度(mof)方法，因为所有模态的参数都是从多个频响函数中同时估计的。

LSCE算法生成结构中不存在的计算模式。使用稳定图，通过检查极点随模式数量增加的稳定性来识别物理模式。物理模式的固有频率和阻尼比往往保持在同一位置，或“稳定”。创建一个并输出频率稳定的极点的固有频率。

[FRF，f]=modalfrf（X1，Y1，fs，winLen，“传感器”，“说”，“测量”，“流动输入”); fn=modalsd（FRF，f，fs，“MaxModes”, 20岁,“FitMethod”，“lsce”);识别物理模式

默认情况下，如果极点的固有频率变化小于1%，则极点被归类为频率稳定。在频率上稳定的极点，如果阻尼比的变化小于5%，就进一步归类为阻尼上稳定的极点。两个标准都可以调整为不同的值。根据稳定极点的位置，选择固有频率373、852.5和1371 Hz。这些频率包含在输出中fn的modalsd，以及其他稳频极点的固有频率。通常需要比实际存在的模态数量更高的模型阶数，才能使用lce算法产生良好的模态参数估计。在这种情况下，四个模态的模型阶数表示三个稳定极点。感兴趣的频率出现在第4行的前三列fn．

physFreq = fn(4，[1 2 3]);

估计固有频率和阻尼，绘制重建和测量的频响。从稳定性图中指定四种模式和物理频率，“PhysFreq”．modalfit仅返回指定模式的模态参数。

modalfit（FRF，f，fs，4，“PhysFreq”，physFreq）

(fn1,根据dr1) = modalfit(润扬悬索桥,f, f 4“PhysFreq”，physFreq）

Fn1 = 1.0e+03 * 0.3727 0.8525 1.3706 dr1 = 0.0008 0.0018 0.0029

接下来，计算FRF并绘制第二组锤击的稳定图，传感器位于不同位置。将频率稳定性标准更改为0.1%，阻尼稳定性标准更改为2.5%。

润扬悬索桥[f] = modalfrf (X2, Y2、fs winLen,“传感器”，“说”，“测量”，“流动输入”); fn=modalsd（FRF，f，fs，“MaxModes”, 20岁,“SCriteria”[0.001 - 0.025]);

在更严格的标准下，大多数极点被归类为频率不稳定。频率和阻尼稳定的极点与平均频响紧密对齐，表明它们存在于测量数据中。

physFreq = fn(4，[1 2 3]);

为这组测量提取模态参数，并与第一组测量的模态参数进行比较。指定驱动点频响的指标，对应于激励和响应测量重合的位置。固有频率一致在1%以内，阻尼比一致在4%以内，表明模态参数从测量到测量是一致的。

[fn2, dr2] = modalfit(润扬悬索桥,f, f 4“PhysFreq”，physFreq，“DriveIndex”[1, 2])

fn2=1.0e+03*0.3727 0.8525 1.3705 dr2=0.0008 0.0018 0.0029

Tdiff2=表格（（fn1-fn2）。/fn1，（dr1-dr2）。/dr1，“VariableNames”, {“diffFrequency”，“扩散阻尼”})

Tdiff2 = 3x2 table diffFrequency diffDamping _____________ ___________ 2.9972 -06 -0.031648 -5.9335e-06 -0.0099076 1.965e-05 0.0001186

当频响中存在测量噪声或频响具有较高的模态密度时，模态参数估计的参数化方法可以提供一个有用的选峰方法和lce方法。最小二乘有理函数(LSRF)方法将共享分母传递函数拟合到多输入、多输出频响中，从而得到模态参数[4]的单一全局估计。使用LSRF方法的过程与使用lce方法的过程类似。您可以使用稳定图来识别静止模态，并提取与识别的物理频率相对应的模态参数。

[FRF，f]=modalfrf（X1，Y1，fs，winLen，“传感器”，“说”，“测量”，“流动输入”); fn=modalsd（FRF，f，fs，“MaxModes”, 20岁,“FitMethod”，“lsrf”);%使用lsfr识别物理模式physFreq = fn(4，[1 2 3]);[fn3, dr3] = modalfit(润扬悬索桥,f, f 4“PhysFreq”，physFreq，“DriveIndex”(1 - 2),“FitMethod”，“lsrf”）

fn3=372.6832 372.9275 852.4986 dr3=0.0008 0.0003 0.0018

Tdiff3=表格（（fn1-fn3）。/fn1，（dr1-dr3）。/dr1，“VariableNames”, {“diffFrequency”，“扩散阻尼”})

Tdiff3 = 3x2 table diffFrequency diffDamping _____________ ___________ -7.8599e-06 0.007982 0.56255 0.83086 0.37799 0.37626

关于参数化方法的最后一个注意:频响估计方法(“子空间”)和模态参数估计方法(“lsrf”)类似于“系统辨识工具箱”中用来拟合时域信号或频响函数的动态模型。如果您有这个工具箱，您可以使用诸如tfest和党卫军. 您可以使用比较和渣油命令。一旦验证了模型的质量，就可以使用它们来提取模态参数。这是用状态空间估计器简单地展示的。

Ts=1/fs；%样品时间%创建用于模型估计的数据对象。EstimationData = iddata(Y0(1:1000)， X0(1:1000)， 1/fs);%创建用于模型验证的数据对象ValidationData = iddata(Y0(1001:2000)， X0(1001:2000)， 1/fs);

确定一个包含馈通项的六阶连续时间状态空间模型。

sys=ssest（估计数据，6，“引线”,真正的)

sys=连续时间识别状态空间模型：dx/dt=ax（t）+bu（t）+ket（t）y（t）=cx（t）+du（t）+e（t）4.6-4.6 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 5 5 5 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 4×4×5×6y1 3.135e-05 2.511e-06 8.634e-06-1.416e-05 2.218e-06-6.271e-06 D=u1 y1 7.564e-09 K=y1 x1 3.513e+07X2-3.244e+06X3-3.598e+07X4-1.059e+07X51.724e+08X6 7.521e+06参数化：自由形式（A、B、C中的所有系数均为自由形式）。馈通：是干扰分量：估计自由系数数：55使用“idssdata”、“getpvec”、“getcov”表示参数及其不确定性。状态：使用SSEST对时域数据“EstimationData”进行估计。拟合估计数据：99.3%（预测焦点）FPE:1.235e-16，MSE:1.189e-16

通过检查模型与验证数据的吻合程度来评估模型的质量。

clf比较(ValidationData sys)%这张图很适合

使用该模型sys计算模态参数。

[fn4，dr4]=modalfit（sys，f，3）；

风力机叶片的模态形状矢量

了解风力机叶片的动态特性对于优化运行效率和预测叶片失效具有重要意义。本节对风力机叶片的实验模态分析数据进行分析，并对叶片的模态振型进行可视化。一个锤在20个位置激励涡轮叶片，一个参考加速度计测量位置18的响应。在叶片的底部安装一个铝块，叶片在垂直于叶片的平坦部分的挥动方向上被激励。为每个位置收集一个FRF。FRF数据由马萨诸塞大学洛厄尔分校的结构动力学和声学系统实验室提供。首先，可视化测量位置的空间布局。

负荷和绘图风力涡轮机叶片频响估计位置18和20。放大前几个峰值。

[FRF，f，fs]=helperImportModalData（）；helperPlotModalAnalysisExample（FRF，f，[18-20]）；

前两种模式在37hz和111 Hz左右出现峰值。绘制稳定图来估计固有频率。14阶模型的前两个返回值在频率和阻尼比上是稳定的。

润扬悬索桥fn = modalsd (f, f,“MaxModes”, 20);physFreq = fn(14，[1 2]);

接下来，使用下面的方法提取前两个模式的模态形状modalfit．根据前面的图，限制适合的频率范围从0到250hz。

[~, ~,女士]= modalfit(润扬悬索桥,f, f, 14日“PhysFreq”，physFreq，“频率范围”250年[0]);

振型量化了结构在每个位置的每个模态的运动幅度。为了估计模态振型矢量，需要频率响应函数矩阵的一行或一列。实际上，这意味着要么在结构的每个测量位置都需要一个激励(在这种情况下，是一个粗纱锤)，要么在每个位置都需要一个响应测量。模态振型可以通过检查频响函数的虚部来显示。为叶片一侧的位置绘制频响矩阵的虚部的瀑布图。将频率范围限制在150赫兹以检查前两种模式。图的顶点代表振型。

Measlocs = [3 6 9 11 13 15 17 19 20];%测量叶片边缘位置helperPlotModalAnalysisExample（FRF，f，measlocs，150）；

图中由峰值等高线表示的形状代表了叶片的第一和第二弯矩。接下来，绘制相同测量位置的模态振型矢量的大小。

helperPlotModalAnalysisExample (ms, measlocs);

虽然振幅的缩放方式不同（振型向量缩放为单位模态A），但振型轮廓形状一致。第一个振型的形状有一个较大的尖端位移和两个节点，其中振幅为零。第二个振型也有一个较大的尖端位移和三个节点。

总结

本实例对粗纱锤激励下的三自由度系统的模态分析仿真数据集进行了分析和比较。它估计的固有频率和阻尼使用稳定图和lce和LSRF算法。模态参数在两组测量中是一致的。在一个单独的用例中，风力涡轮机叶片的模态振型使用频响矩阵的虚部和模态振型向量来可视化。

承认

感谢马萨诸塞大学洛厄尔分校结构动力学和声学系统实验室的Peter Avitabile博士，为风力涡轮机叶片实验数据的收集提供了便利。

参考文献

[1] 布兰特，安德斯。噪声和振动分析:信号分析和实验程序．英国奇切斯特:约翰·威利父子公司，2011。

[2] Vold, Håvard, John Crowley和G. Thomas Rocklin。"估计频响函数的新方法"声音和振动第18卷，1984年11月，第34-38页。

彼得·范·奥弗切和巴特·德摩尔。《N4SID:确定-随机联合系统识别的子空间算法》。自动化．第30卷，1994年1月，第75-93页。

[4] Ozdemir, A. A.和S. Gumussoy。“基于向量拟合的系统辨识工具箱中的传递函数估计”。国际自动控制联合会第20届世界大会论文集．2017年7月，法国图卢兹。