0
  • 聊天消息
  • 系统消息
  • 评论与回复
登录后你可以
  • 下载海量资料
  • 学习在线课程
  • 观看技术视频
  • 写文章/发帖/加入社区
会员中心
创作中心

完善资料让更多小伙伴认识你,还能领取20积分哦,立即完善>

3天内不再提示

MATLAB中的振动分析与信号处理分析

西门子EDA 来源:MATLAB 作者:刘海伟 2021-08-16 10:09 次阅读

模态分析主要研究频率域内系统动态特性。

通过模态分析方法搞清楚了结构物在某一易受影响的频率范围内的各阶主要模态的特性,就可以预言结构在此频段内在外部或内部各种振源作用下产生的实际振动响应。

模态分析主要分为解析模态分析和试验模态分析

解析模态分析

通常我们针对一个线性定常系统进行动力学描述可以得到方程组:

56576ce0-fdc9-11eb-9bcf-12bb97331649.png

其中,[M] 是质量矩阵,[C] 是阻尼矩阵,[K] 是刚度矩阵,{x(t)} 是位移向量, {F(t)} 是力矩阵。 我们的目标是求解这个线性定常系统振动微分方程组得到 {x(t)},也就是系统上各点随时间的位移。 对于单自由度的系统是方便求解的,所以对于这种多自由度系统我们首先想到的是通过坐标变换,用一组新的正交基(模态空间里的基)去描述 {x(t)},例如 [P⁻¹]x(t),来实现对方程组 (1) 解耦,从而将问题转化为互相独立的子方程(组),用更少自由度甚至单自由度的方程进行求解。 其中[P⁻¹]就是模态矩阵,其每列为模态振型,它描述的是在新的解耦后的坐标系中的各维坐标与原来坐标系中各个维度的线性关系。

5684d310-fdc9-11eb-9bcf-12bb97331649.png

例如对于一个简化的多自由度的弹簧系统,这个线性定常系统由 3 个相同重量的模块组成,m₁=m₂=m₃=m,他们中间用弹簧链接, 为了简化问题,我们设弹簧的劲度系数 k₀=k₁=k₂=k₃=k,阻尼系数 c₀=c₁=c₂=c₃=0。 定义 x₁,x₂,x₃‍为每个质量块的位移,另外 k₀ ,k₄边缘两端的位移。由于两端固定,都为0。每个质量块的运动方程分别为:

56918e7a-fdc9-11eb-9bcf-12bb97331649.png

将上述方程写为矩阵形式,上述运动学方程组可以简化为:

56a2bbfa-fdc9-11eb-9bcf-12bb97331649.png

其中

56aea1ae-fdc9-11eb-9bcf-12bb97331649.png

对于方程组 (2),如何进行坐标解耦呢? 计算时运动学方程组(2) 右侧项可以忽略, 只保留质量矩阵项 [M] 和刚度矩阵 [K] 项,即

56b9f824-fdc9-11eb-9bcf-12bb97331649.png

对于方程组(3) 通常的做法是转换为:

56c65b14-fdc9-11eb-9bcf-12bb97331649.png

本质对方程组(4) 解耦实际上就是求解质量矩阵关于刚度矩阵的广义特征值的问题。也就是计算得到特征值,

56d28628-fdc9-11eb-9bcf-12bb97331649.png

和特征向量,

56dba8ac-fdc9-11eb-9bcf-12bb97331649.png

使得[M][P]=[K][P][D] (5) 其中特征向量 [P] 即为模态向量(空间基向量),为对应的特征值对角阵对应解耦后固有频率的平方,即

56f7d342-fdc9-11eb-9bcf-12bb97331649.png

具体此处不做推导,但可以简单的从必要性上进行解释: 假设已经通过空间变换矩阵得到新的坐标

5706efb2-fdc9-11eb-9bcf-12bb97331649.png

我们计算一下特征向量矩阵是否将原始方程组坐标由 {X} 变换为后 {Q} 得到单自由度的二阶常微分方程组。 我们将{X}=[P]{Q} 带入方程(3)得到

571353c4-fdc9-11eb-9bcf-12bb97331649.png

同时我们将(5) 带入(6) 可以得到

57202482-fdc9-11eb-9bcf-12bb97331649.png

在 [K][P] 均可逆的条件下,我们得到方程

57392f2c-fdc9-11eb-9bcf-12bb97331649.png

即:

57466160-fdc9-11eb-9bcf-12bb97331649.png

也就是完全解耦的单自由度的二阶常微分方程,接下来可以单独求解 q₁(t), q₂(t), q₃(t), 最后只需要再做一次 [P] 变换将模态空间坐标变换到物理坐标即可。

‍‍‍‍‍‍‍‍‍

575239fe-fdc9-11eb-9bcf-12bb97331649.png

% 'M' 是质量矩阵,'K' 是刚度矩阵. 'm' 质量块质量,单位 Kgs

% 'k' 刚度系数,单位N/m. 'P' 和'D' 是特征向量和特征值.

m = 0.003;

k = 180000;

M = m*eye(3);

K = k*[2 -1 0 ;

-1 2 -1 ;

0 -1 2 ];

% P为特征向量:变换矩阵,D为特征值:固有频率的平方,w_nat 包含自然响应频率.

[P,D] = eig(K,M)

w_nat = sqrt(D)

% 我们查看低阶模态振型,也就是低频振型,可以尝试设置number

% 查看三阶模态振型'EIGS' 函数可以用来计算特征值和特征向量

number = 3;

[P,D]=eigs(K,M,number,'smallestabs');

% 因为系统两端固定,模态振型坐标在这两端为0

vect_aug1 = [0 0 0;P;0 0 0];

c = ['m','b','r'];

figure(1)

for i=1:size(vect_aug1,2)

plot(vect_aug1(:,i),c(i))

hold on;

end

575c8e04-fdc9-11eb-9bcf-12bb97331649.png

最终

57a82e68-fdc9-11eb-9bcf-12bb97331649.png

57d2c0f6-fdc9-11eb-9bcf-12bb97331649.png

的求解以及绘制都可以用通过 MATLAB 脚本实现。大家可以自己尝试。 当然对于我们的系统不可能是这种简单的系统,通常要结合有限元的手段进行计算。 MATLAB 中的 Partial Differential EquationToolbox 对于满足我们一些基础的需求可以提供求解方案。 我们看一个基于 MATLAB 有限元计算模态的示例。 本示例是对 KinovaGen3 机械臂肩部关节部件进行模态和频率响应分析。机械臂通过多个关节链接,一端固定。这些链接结构强度要够大以避免电机带动负载运动时产生振动。 机械臂终端的负载会让每个链接处产生压力。压力的方向取决于负载的方向。此示例主要展示如何分析 Kinova Gen3 超轻量级机械臂的肩部关节连接部件在一定压力下可能的形变。

模态分析MATLAB 中有限元求解流程

57fe67d8-fdc9-11eb-9bcf-12bb97331649.png

model = createpde('structural','modal-solid');

importGeometry(model,'Gen3Shoulder.stl');

generateMesh(model);

structuralProperties(model,'YoungsModulus',E, ...

'PoissonsRatio',nu, ...

'MassDensity',rho);

将 facelabel 可视化出来方便设置边界约束和负载。

5822fe9a-fdc9-11eb-9bcf-12bb97331649.png

face3 是固定的,face4 是活动的。设置约束

structuralBC(model,'Face',3,'Constraint','fixed');

在关心的频率范围进行模型求解。

RF = solve(model,'FrequencyRange',[-1,10000]*2*pi);

modeID = 1:numel(RF.NaturalFrequencies);

通过对结果除以2pi,得到Hz单位的频率值:

tmodalResults = table(modeID.',RF.NaturalFrequencies/2/pi);

tmodalResults.Properties.VariableNames = {'Mode','Frequency'};

disp(tmodalResults);

5832f610-fdc9-11eb-9bcf-12bb97331649.png

同样我们需要可视化模态振型。最好的方式是以各阶模态频率的简谐振动动画显示出来,此处显示前六阶模态振型:

频率响应分析模拟机械臂关节在压力载荷下的动力学,假设附加其上的连杆对各半面分别施加大小相等方向相反的压力。分析面上某点的频率响应和形变。

同样跟上述流程一样,创建结构,导入几何形状,生成网格。

其他过程跟模态分析相同,区别在于加一个力,使用 pressFcnFR 函数在面 4 上施加边界载荷。 这个函数作用一个推力和一个扭转压力信号。推压分量是均匀的。扭力对左侧面施加正压力,对右侧施加负压力。

structuralBoundaryLoad(fmodel,'Face',4,'Pressure',@(region,state)pressFcnFR(region, state),'Vectorized','on');

这个压力函数跟频率无关,作用于所有频率。

同样进行求解

R = solve(fmodel,flist,'ModalResults',RF)

我们可以看面 4对应的负向负荷最大的点(0.003; 0.0436; 0.1307)对应的位移

queryPoint= [0.003; 0.0436; 0.1307];

queryPointDisp =R.interpolateDisplacement(queryPoint);

58e286f2-fdc9-11eb-9bcf-12bb97331649.png

响应的峰值出现在 2662Hz 附近,与二阶模态接近。在接近 1947Hz 的一阶模态上也会出现小的响应。

通过使用 max 函数找到峰值响应频率对应的峰值和索引

[M, I] = max(abs(queryPointDisp.uy))

M = 1.1256e-04

I = 152 绘制峰值响应频率处的变形。

可以看到所施加的载荷主要激发了部件的开口模态和弯曲模态。 通过上面两个示例,我们针对计算模态的场景可以在 MATLAB 中通过相应的内置函数做探索研究。

试验模态分析

试验模态分析主要是通过实测实验数据进行频率响应估计并计算相应模态参数的一种方法。

我们通过 MATLAB 自带的一个例子来介绍试验模态分析。

详见:https://ww2.mathworks.cn/help/releases/R2021a/ident/ug/modal-analysis-of-a-flexible-flying-wing-aircraft.html

这是明尼苏达大学无人飞行器实验室的小型柔性飞翼飞机的示例。这个例子展示了柔性机翼飞机弯曲模态的计算过程。

通过沿翼型进行机翼的振动响应的多点采集获得数据,用于辨识系统的动态模型。

从辨识出的模型中提取模态参数。

将模态参数数据与传感器位置信息相结合,可视化机翼的各种弯曲模态。

实验设置

实验的目的是收集飞机在外部激励下不同位置的振动响应。

这架飞机悬挂在一个木制框架上,其重心由一根弹簧支撑。该弹簧具有足够的弹性保证弹簧-质量振荡的固有频率不会干扰飞机的基频。通过一个电动激振器在飞机中心附近施加输入力。

沿着翼展放置 20 个加速度计来收集振动响应,如下图所示

激振器输入命令指定为一个恒定振幅的 chirp 信号 Asin(ω(t)t), chirp 信号的频率随时间线性增加,ω(t)=c0+c1t, 频率范围为 3 - 35hz。试验数据由两个加速度计(在 x 方向的前缘和后缘)一次收集。总共进行了 10 组实验来收集所有 20 个加速度计的响应。加速度计和力传感器的测量频率都是 2000hz。

数据准备

数据由 10 组单输入/双输出信号表示,每组包含 600K 个样本。

变量 MeasuredData 是一个结构体。结构体的每个域还是一个结构体,包含两个响应 y 和对应输入 u。随机绘制第一次实验的数据。

fs = 2000; % data sampling frequency

Ts = 1/fs; % sample time

y = MeasuredData.Exp1.y; % 加速度计的输出值,每个包含两列

u = MeasuredData.Exp1.u; % input force data

t = (0:length(u)-1)' * Ts;

5982e0c0-fdc9-11eb-9bcf-12bb97331649.png

为了用于系统辨识,将数据封装到 iddata 对象中,并将 10 次试验合并。

Data = merge(Data{:})

59d0b3e0-fdc9-11eb-9bcf-12bb97331649.png

系统辨识

我们想识别一个频率响应与实际飞机尽可能接近的动态模型。

动态模型将系统的输入和输出之间的数学关系转化为微分方程或差分方程。例如传递函数和状态空间模型都是动态模型。

动态模型可以通过在时域或频域对试验测量数据运行 fest 和 sest 之类的模型估计命令来创建。

这个例子中,我们先用 etfe 命令进行传递函数估计,从而将测量数据从时域转换为频率响应。然后利用估计的频响函数用于辨识飞机振动响应的状态空间模型。

当然直接利用时域数据进行状态空间模型辨识也是可以的。但频响函数的形式可以将大数据集压缩成更少的样本,并且更方便的在相关的频率范围调整估计过程。

针对每组实验数据(两输出/单输入)进行频响函数(FRF)估算。使用 24000 个频率点进行无窗响应计算。

G = cell(1, 10);

N = 24000;

for k = 1:10

% Convert time-domain data into a FRF using ETFE command

Data_k = getexp(Data, k);

G{k} = etfe(Data_k, [], N); % G{k} is an @idfrd object

end

G = cat(1,G{:}); % 将所有的频响函数合并为一个(20输出/一个输入)的频响

G.OutputName = 'y'; % name outputs 'y1', 'y2', ..., 'y20'

G.InterSample = 'bl';

我们随便挑选第 1 个和第 15 个输出幅值进行绘制来看一下频率响应函数的估计结果。我们关注的频率范围(4 - 35hz)。

59dea266-fdc9-11eb-9bcf-12bb97331649.png

频响函数显示至少 9 个谐振频率(峰值)。我们最关心飞机的机翼弯曲模态,这些模态主要集中在 6- 35hz 的频率范围,因此我们只选择这个范围的频响。

f =G.Frequency/2/pi; % 单位变换

Gs = fselect(G,f>6 & f<=32)    % 选择频响频率范围 (6.5 - 35 Hz)

接下来,计算一个状态空间模型来逼近 Gs。子空间估计器 n4sid 提供了一个快速的非迭代估计。状态空间模型参数配置会影响精度:

1. 模型阶数的选择。通常要多次尝试。

2. 模型结构。例如是否包含馈通项(状态空间模型中的 D 矩阵是否为零),等等。

3. 设置权重项,确保低振幅和高振幅对结果有相同的影响。

FRF =squeeze(Gs.ResponseData);

Weighting = mean(1./sqrt(abs(FRF))).';

n4Opt =n4sidOptions('EstimateCovariance',false,...

'WeightingFilter',Weighting,...

'OutputWeight',eye(20));

sys1 = n4sid(Gs,24,'Ts',0,'Feedthrough',true,n4Opt);

Fit = sys1.Report.Fit.FitPercent'

通过查看 Fit 的效果,显示这 20 个响应中最好的和最差的

59f0835a-fdc9-11eb-9bcf-12bb97331649.png

可以看到 sys1 结果仍然有待提高。通过对模型参数进行非线性最小二乘优化迭代,可以进一步改善模型的拟合效果。这可以使用 sest 命令来实现。这一次也估计了参数协方差。

ssOpt = ssestOptions('EstimateCovariance',true,...

'WeightingFilter',n4Opt.WeightingFilter,...

'SearchMethod','lm',... % use Levenberg-Marquardt search method

'Display','on',...

'OutputWeight',eye(20));

sys2 = ssest(Gs, sys1, ssOpt); % estimate state-space model (this takesseveral minutes)

Fit = sys2.Report.Fit.FitPercent'

我们同样看看拟合效果:最好的和最差的幅值进行可视化。同时将 1-sd 置信区间可视化出来。

5a00c92c-fdc9-11eb-9bcf-12bb97331649.png

优化后的状态空间模型 sys2 在 7 - 20hz 区域的频响拟合效果很好。多数共振位置的不确定性区间都很窄。我们最开始设置的是一个 24 阶模型,这意味着在系统 sys2 中最多有 12 个谐振模态。使用 modalfit 命令获取这些模态的固有频率。

f = Gs.Frequency/2/pi; % data frequencies (Hz)

fn = modalfit(sys2, f, 12); % naturalfrequencies (Hz)

5a377d0a-fdc9-11eb-9bcf-12bb97331649.png

fn 中的值表明两个频率非常接近 7.8 Hz,三个接近 9.4 Hz。查看这些频率附近的各点频率响应,峰值位置在输出中确实发生了一些偏移。

这些差异可以通过更好地控制实验过程,直接利用频率为中心的狭窄范围内的时域数据进行直接辨识,或将这些频率附近的频率响应拟合为单个振动模态。本例中不讨论这些替代方案。

模态参数计算

现在我们可以使用模型 sys2 来提取模态参数。通过查看频响结果找到 10 个模态频率,大约在频率 [5 7 10 15 17 23 27 30]Hz 附近左右。通过使用 modalsd 函数可让估计更加准确,该函数尝试不同模型的阶数来检查模态参数的稳定性。

通过稳定图可以看到一组更好的自然响应频率

Freq = [7.8 9.4 15.3 19.3 23.0 27.3 29.231.7];

我们使用 Freq 向量的值作为从模型 sys2 中选择主要模态的基准。使用 modalfit 函数实现

[fn,dr,ms] = modalfit(sys2,f,length(Freq),'PhysFreq',Freq);

fn 是固有频率 (Hz), dr 是相应的阻尼系数,ms 是模态振型向量 (每个固有频率一列)。

模态振型可视化

我们使用上述模态参数可视化各种弯曲模态。此外,我们需要各测量点位置的信息。

模态振型包含在矩阵 ms 中,其中每一列对应于一个频率的振型。通过在加速度传感器位置坐标上叠加模态振型的振幅并以模态固有频率改变振幅来做动画展示。

结论

这个例子展示了一种基于参数模型辨识的模态参数估计方法。利用 24 阶状态空间模型,在 6 ~ 32 Hz 频率范围内提取了 8 个稳定的振动模态。将模态信息与位置信息相结合,可视化各种弯曲模态。

当然,您也可以对其他设备例如风力发电机的叶片振型进行计算,了解风力机叶片的动态特性从而优化运行效率和预测叶片失效,可以按同样的思路实现。

声明:本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人,不代表电子发烧友网立场。文章及其配图仅供工程师学习之用,如有内容侵权或者其他违规问题,请联系本站处理。 举报投诉
  • matlab
    +关注

    关注

    186

    文章

    2983

    浏览量

    231281
  • 仿真分析
    +关注

    关注

    3

    文章

    105

    浏览量

    33717

原文标题:MATLAB 中的振动分析与信号处理 —— 模态分析

文章出处:【微信号:Mentor明导,微信公众号:西门子EDA】欢迎添加关注!文章转载请注明出处。

收藏 人收藏

    相关推荐

    函数信号分析仪的原理和应用场景

    )、脑电图(EEG)等。通过分析这些信号的频谱特性,可以揭示生物体的生理状态和病理变化。 振动分析:在机械工程,函数
    发表于 01-20 14:13

    卡尔曼滤波在信号处理的应用分析

    卡尔曼滤波在信号处理的应用十分广泛,其强大的滤波和预测能力使其成为信号处理领域的一种重要工具。以下是对卡尔曼滤波在
    的头像 发表于 12-16 09:14 1689次阅读

    Simulink与 MATLAB 的结合使用 Simulink信号处理方法

    在工程和科学研究信号处理是一个重要的领域,涉及到信号的采集、分析处理和生成。
    的头像 发表于 12-12 09:25 566次阅读

    提高设备健康状态 KMWIS无线振动分析仪解决化工集团风机振动异常问题!

    风机振动过大这一难题怎么解决?需要用到专业的仪器来进行振动分析!通过对风机振动信号的采集、处理
    的头像 发表于 12-04 10:29 132次阅读
    提高设备健康状态  KMWIS无线<b class='flag-5'>振动</b><b class='flag-5'>分析</b>仪解决化工集团风机<b class='flag-5'>振动</b>异常问题!

    分析滤波器在信号处理应用

    滤波器在信号处理的应用十分广泛,其主要功能是从信号中去除不需要的频率成分,保留所需的频率成分,从而实现对信号的有效
    的头像 发表于 11-27 15:56 1144次阅读

    信号分析的目的意义是什么

    分析的定义、基本原理、应用领域以及信号分析的重要性。 一、信号分析的定义 信号
    的头像 发表于 06-03 10:31 1539次阅读

    信号分析的基本思想是什么

    信号进行数学处理,提取信号的有用信息,以便对信号进行识别、分类、估计和预测。信号
    的头像 发表于 06-03 10:28 951次阅读

    信号分析处理信号与系统的区别

    信号分析处理信号与系统是电子工程和信息科学领域中的两个重要概念。尽管它们在某些方面有相似之处,但它们之间存在明显的区别。本文将详细探讨这两个概念的定义、特点、应用以及它们之间的联系
    的头像 发表于 06-03 10:15 3019次阅读

    MATLAB信号处理常用函数详解

    MATLAB是一款功能强大的数学软件,尤其在信号处理领域,它提供了众多的函数和工具箱,使得信号分析
    的头像 发表于 05-17 14:31 2687次阅读

    基于MATLAB信号处理系统与分析

    基于MATLAB信号处理系统与分析,包括信号的导入、预处理
    的头像 发表于 05-17 14:24 1277次阅读

    信号分析信号处理必须遵循的原则

    在信息技术的快速发展信号分析信号处理作为信息科学的重要组成部分,扮演着至关重要的角色。无论是通信、控制、图像
    的头像 发表于 05-17 14:19 1288次阅读

    信号分析信号处理的基本方法有哪些

    在电子工程、通信、生物医学工程、地球物理学等众多领域中,信号分析处理扮演着至关重要的角色。信号分析是指对
    的头像 发表于 05-16 17:25 3285次阅读

    信号分析信号处理的区别

    在通信、电子工程、生物医学工程、地球物理学等众多领域中,信号分析信号处理是两个至关重要的概念。它们都是对信号进行
    的头像 发表于 05-16 17:16 1059次阅读

    工程监测振弦采集仪在振动监测的应用与数据处理技术

    工程监测振弦采集仪在振动监测的应用与数据处理技术 振弦采集仪是一种用于振动监测和分析的仪器设备。它采用振弦传感器作为
    的头像 发表于 04-07 13:59 615次阅读
    工程监测振弦采集仪在<b class='flag-5'>振动</b>监测<b class='flag-5'>中</b>的应用与数据<b class='flag-5'>处理</b>技术

    振动分析仪如何工作?

    振动分析仪基本上是一台通过一个或多个加速度计记录振动的计算机。加速度计内部的振动运动被转换成与加速度成比例的电流。该信号在计算机
    的头像 发表于 02-21 18:08 910次阅读
    <b class='flag-5'>振动</b><b class='flag-5'>分析</b>仪如何工作?