cmi码型变换matlab程序_【转载】动⼒时程分析⼊门[含
Matlab源码]
Duhamel积分、Fourier变换法、中⼼差分⽅法、Newmark法、功率谱法……傻傻分不清楚，相信有不少⼩伙伴像笔者⼀样好奇：这些⽅法精度如何，结果有什么特点？以下笔者将以单⾃由度弹性体系为例，对⽐这⼏种⽅法的精度。为了便于感兴趣的⼩伙伴重现，同时附上了Matlab代码。
算例为《结构动⼒学》(刘晶波主编)习题5.2
1  Duhamel积分⽅法
Duhamel积分是⼀种时域分析⽅法，它将荷载分解为⼀系列脉冲，获得每⼀个脉冲作⽤下结构的反应，然后叠加每⼀脉冲作⽤下的反应得到结构总的反应。Duhamel积分⽅法以积分的⽅式给出了体系运动的解析表达式(忽略离散采样带来的误差的情况下)，但从实际应⽤上看，当采⽤数值积分时，其计算效率不⾼，因为对于计算任⼀个时间点t的反应，积分都要从0积到t，⽽实际要计算⼀时间点系列，可能要⼏百到⼏千个点，计算量很⼤。因为使⽤了叠加原理，仅适⽤于线弹性分析。
Matlab程序如下：
function y=duhamel(dt)
%duhamel
m=17.5; %质量
k=875.5; %刚度
zeta=0.14138;%阻尼⽐
%dt=0.01; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
w0=sqrt(k/m);
w1=w0*sqrt(1-zeta^2);
t=t0:dt:t2;
y=t;
for i=1:(length(t))
x=linspace(t(1),t(i));
px=(100*x).*(x>=0&x<=0.4)+(80-100*x).*(x>0.4&x<0.8);
a=px.*exp(zeta*w0*x).*cos(w1*x);
A=trapz(x,a);
b=px.*exp(zeta*w0*x).*sin(w1*x);
B=trapz(x,b);
y(i)=exp(-zeta*w0*t(i))*(A*sin(w1*t(i))-B*cos(w1*t(i)))/(m*w1);
end
ymax=max(y)
figure
plot(t,y);
2 Fourier变换法
Fourier变换法是⼀种频域分析⽅法，其基本计算步骤是：
1)对外荷载做Fourier变换，得到外荷载的Fourier谱；
2)利⽤复频反应函数得到反应的频域解；
3)应⽤Fourier逆变换，得到反应的时域解。
在⽤频域法分析中涉及到两次Fourier变换，均为⽆穷域积分，特别是Fourier逆变换，被积函数是复数，有时涉及复杂的围道积分。
在⽤频域法分析中涉及到两次Fourier变换，均为⽆穷域积分，特别是Fourier逆变换，被积函数是复数，有时涉及复杂的围道积分。当外荷载是复杂的时间函数(如地震动)时，⽤解析型的Fourier变换⼏
乎是不可能的，实际计算中⼤量采⽤的是离散Fourier变换。因为使⽤了叠加原理，仅适⽤于线弹性分析。
Matlab程序如下：
function [un,tf]=qiaoFFT(dt,N)
% Fourier
%% 执⾏FFT点数为64
% 构建原信号
% dt=0.1;
% N=64;
t=[0:N-1]*dt;  % 时间序列
xn=(100*t).*(t>=0&t<=0.4)+(80-100*t).*(t>0.4&t<0.8);
subplot(2,2,1)
plot(t,xn)  % 绘出原始信号
xlabel('时间/s'),title('原始信号')
axis([0 6.4 0 50])  % 调整坐标范围
% FFT分析
NN=N;  % 执⾏64点FFT
XN=fft(xn,NN);  %
f0=1/(dt*NN);  % 基频
f=[0:NN-1]*f0;  % 频率序列
A=real(XN);  % 幅值序列
subplot(2,2,2);
stem(f,A),xlabel('频率/Hz')  % 绘制频谱
axis([0 10 -200 200])  % 调整坐标范围
title('荷载傅⾥叶谱');
%% H(iw)
k=875.5;
zeta=0.14138;
fn=1.126288;
HN=ones(1,NN)./(1-(f/fn).*(f/fn)+2i*zeta*(f/fn))/k;
for i=1:NN/2
HN(NN+1-i)=conj(HN(i));
end
UN=HN.*XN;
subplot(2,2,3);
subplot(2,2,3);
stem(f,abs(HN)),xlabel('频率/Hz')  % 绘制频谱
axis([0 10 -0.01 0.01])  % 调整坐标范围
title('复频反应函数');
un=ifft(UN,NN);
subplot(2,2,4)
plot(t,un)  %
xlabel('时间/s'),title('ut')
axis([0 6.4 -0.1 0.1])  % 调整坐标范围
tf=t
ymax=abs(max(un))
荷载Fourier谱、复频反应函数及时域位移反应如下图所⽰。
3 中⼼差分法
中⼼差分⽅法⽤有限差分代替位移对时间的求导(即速度和加速度)，在计算i+1时刻的运动时，需要已知i和i-1两个时刻的运动，属于两步法，它具有2阶精度，是有条件稳定的，稳定条件为, 是显式积分⽅法，不需要对刚度矩阵求逆，具有较⾼的计算效率。
Matlab程序如下：
function u=central(dt)
% central
m=17.5; %质量
k=875.5; %刚度
c=35;%阻尼⽐
%dt=0.1; %时间步长
t0=0; %起始时间
t2=6.4; %结束时间
t=t0:dt:t2;
u=t;u(1)=0;u(2)=0;
k1=m/dt/dt+c/2/dt;
b=m/dt/dt-c/2/dt;
c=2*m/dt/dt;
for i=2:(length(t)-1)
x=t(i);
pi=(100*x)*(x>=0&x<=0.4)+(80-100*x)*(x>0.4&x<0.8);
fs=k*u(i);
pi1=pi-fs+c*u(i)-b*u(i-1)
u(i+1)=pi1/k1;
end
ymax=max(u)
matlab求导
4 Newmark⽅法
Newmark⽅法同样将时间离散化，运动⽅程仅要求在离散的时间点上满⾜。与中⼼差分法不同的是，它不是⽤差分对i时刻的运动⽅程展开，得到外推计算i+1时刻位移的公式，⽽是通过对加速度的假设，以i时刻的运动量为初始值，通过积分得到计算i+1时刻的运动公式，计算过程中需要对刚度矩阵求逆，是隐式⽅法。当γ=1/2, β=1/4时，是⽆条件稳定的，就是常加速度法；当γ=1/2,β=1/6时，就是线性加速度法。
Matlab程序如下：
function u=newmark(dt,beta)
% newmark
m=17.5; %质量
k=875.5; %刚度
c=35;%阻尼⽐