昨天发表的《一文搞懂Matlab时频傅里叶转换》一文中存在严重的理论错误,感谢CAE中学生的热心群友提醒,错误修正的地方将在下文详细说明。
众所周知,振动中会用到响应谱,响应谱分析的主要用途是代替瞬态动力学的时程分析,如下图,将一个复杂的时域信号使用拟合的方法,分离成众多的谐波信号,即正弦或余弦信号,由于频谱图省略了各谐波的相位信息,只记录它们的频率与幅值,所以是近似且快速的计算手段。这个内容在教主的动力学课程中说得很详细了。
严格意义上,响应谱分析中的边界条件是一种响应,因为输出的反应谱是位移(速度、加速度或力)与频率之间的关系,这种关系称为频谱。
很多教程说响应谱不需要工程师自己生成,可以直接查询相关标准。在一些标准试验中的确是这样,比如舰船的海浪冲击可以直接查询GJB150相应文件、地震对建筑物的冲击可以查询国标或欧标的建筑规范标准。
但是大多数时候,我们要计算的频谱并不是标准手册上有的,而是根据我们自身产品在冲击或振动试验中测得的。一般使用传感器测试产品某处的时间历程响应,即时间与位移(或速度、加速度)的时域信号。如果是在固定约束处,也可直接使用振动台的时域信号,而无需传感器。将这个时域信号转换为频域信号,导入有限元软件中对产品进行仿真计算。如果振动台控制器自带频谱输出,那么都不需要时频转换了。
工作中,甲方也不会给你配一个人专门进行信号处理,往往甲方自己也不清楚边界或加载情况,更有甚者试验都让仿真者参与设计,所以不要想着有人给你个频谱,你拿进ansys/abaqus中咔咔一顿操作交差完事,这么轻松为啥甲方自己不搞呢。
在有限元教材中不会讲Matlab的时频转换,在Matlab的教材中又不会讲振动,这个最重要的环节似乎被各位大佬故意遗忘了,我们本文中要探讨的正是这个中间被忽略的部分。
上面划线这部分存在理论错误。先看周老师的说明:
从严格意义上来说,响应谱分析中的边界条件不是传统意义上的载荷,而是一种响应,因为输入的反应谱是位移、速度、加速度或力与频率之间的关系,这种关系被称为频谱。例如:在振动试验台上安装4个单自由度的弹簧质量系统,频率分别为f1、f2、f3、f4,且f1<f2<f3<f4。当振动台以频率f1激振时,可得4个系统的位移响应图,再叠加f3频率激振效果,可得4个系统的位移响应图,最后叠加所有频率的综合激振效果,得到的曲线为频谱,如图 4-1-2 所示。
个人理解,振动台只是起到扫频作用,比如从0~1000HZ扫频,在某些频率下,产品会发生共振而导致振幅(位移、速度或加速度)变大,记录下0~1000HZ下产品所有的振幅,形成响应谱图。
所以可见,响应谱是产品在连续频率的输入载荷下表现出的自身响应特性。同样0~1000HZ扫频,产品A与产品B的响应谱可能完全不同。
在Matlab中,时频转换是通过快速傅里叶变换FFT完成了,我们不需要详细研究FFT内部转换过程,只需要直接调用就行了,以下结合实例介绍如何转换。
本文使用的是16版本的Matlab,按照Matlab桀骜不驯的性格,每年升级都会悄悄咪 咪地改动几个命令,所以不能保证适用于其他版本的Matlab。
Step1 数据分析
如图,打开ZXS.csv文件,可以看到这个信号共1s,每个采样数据间隔0.01s,即采样率100。这里一共有100个数据,注意数据的数量必须是偶数个,数据最上方不能有表头,否则需要修改后文的代码。
Step2 Matlab转换。
首先将工作目录设置到ZXS.csv所在的文件夹。
输入如下代码,%之后的文字为注释。
需要注意的是,需要核对输入文件名称ZXS.csv是否与我们保存的数据文件名相同,核对采样率是否正确。
clc;clear all;close all %清理屏幕与缓存
data = xlsread('ZXS.csv'); %19版本读取命令为data = readmatrix('ZXS.csv')
t = data(:,1); %对应第一列
y1 = data(:,2); %对应第二列
figure
y=y1; % 读取时域数据
Fs=100; % 采集频率,1秒采集了多少个点(需要自行根据数据更改)
T=1/Fs; % 采集时间间隔
N=length(y); % 采集信号的长度
t=(0:1:N-1)*T; % 定义整个采集时间点
t=t'; % 转置成列向量
% fft变换
Y=fft(y); % Y为fft变换结果,复数向量
Y=Y(1:N/2+1); % 只看变换结果的一半即可
A=abs(Y); % 复数的幅值(模)
f=(0:1:N/2)*Fs/N; % 生成频率范围
f=f'; % 转置成列向量
% 幅值修正
A_adj=zeros(N/2+1,1);
A_adj(1)=A(1)/N; % 频率为0的位置
A_adj(end)=A(end)/N; % 频率为Fs/2的位置
A_adj(2:end-1)=2*A(2:end-1)/N;
%写出数据
dataa(:,1)=f; %频率
dataa(:,2)=A_adj; %幅值
str = ['DXS.csv']; %写出文件名称
xlswrite(str,dataa); %19版本为 writematrix(str,dataa)
% 绘制频率幅值图和频谱图
subplot(2,1,1)
plot(t,y,'LineWidth',1)
xlabel('时间')
ylabel('信号值')
title('时域信号')
subplot(2,1,2)
plot(f,A_adj,'LineWidth',3)
xlabel('频率(Hz)')
ylabel('幅值(修正后)')
title('FFT变换幅值图')
输入后,点击上方运行工具,将自动导出转换的频域信号DXS.csv,并生成图像,上面显示的是输入的时域图,下面显示的是输出的频域图。可以看到,这个型号是由三个信号组成,分别是频率=1幅值=5的信号1,频率=8幅值=1的信号2,频率=10幅值=0.5的信号3。
到输出的DXS.csv文件中,插入散点图能更清晰看到组成。
可以看到,频域信号底部前后各扩展了1hz,形成了三角形,这是因为频谱分辨率导致,分辨率=1/总时间,我们的信号总时间是1s,所以分辨率是1Hz。
现在,你可以将自己的时频信号导入测试了。
动力学的仿真实操其实很简单,与静力学无异,难就难在理论的解读,我也感觉到自身理论的短板,也请大家多大指教批评。
本文特别感谢卧虎藏龙的CAE中学生的初一(2)班的jiff、黄荣等大佬的热心解答。
参考书籍:《ANSYS Workbench 有限元分析实例详解 动力学》—周炬、苏金英
参考视频:《ANSYS高级视频教程》——张晔