线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书
此类帖子之前有人发过，不过我拿到matlab里调试不出来，于是又整理一个，以便交流学习。改进之处在于应用纽马克beta法，将beta改为1/4，即为纽马克beta法 ，程序如下
function maker4(m,c,k,x0,v0,time,delt)
%线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书
%将beta改为1/4，即为纽马克beta法
%m-为质量，c-阻尼，k-刚度，x0-初位移，v0-初速度；time-仿真时间，delt-时间步长,num-自由度数
% MX''+CX'+KX=F
% X(t+delt)=^(-1){F(t+delt)-C-K}
% X(t+delt)'=X(t)'+delt*/2
% X(t+delt)=X(t)+delt*X(t)'+(delt^2)*X(t)''/3+(delt^2)*X(t+delt)''/6
n=time/delt;
disp=zeros(3,n);%设定存储位移矩阵
m_inv=inv(m+delt*c/2+delt^2*k/6);
=eig(inv(m)*k);
diag(sqrt(fre));%固有频率
i=1;
beta=1/4;
for t=0:delt:time
f=';%外扰力
if t==0
xdd0=inv(m)*(f-k*x0-c*v0);%初始加速度
else
xdd=m_inv*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-beta)*delt^2*xdd0));%计算加速度
x=m_inv*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);%计算位移
xd=v0+delt*(xdd0+xdd)/2;%计算速度
xdd0=xdd;
v0=xd;
x0=x;
disp(:,i)=x0;
i=i+1;
end
end
t=1:n;
figure('numbertitle','off','name','自由度1的位移','pos',);
plot(t,disp(1,:)),grid,xlabel('时间(s/10)'),title('自由度1的时程曲线');
figure('numbertitle','off','name','自由度2的位移','pos',);
plot(t,disp(2,:)),grid,xlabel('时间(s/10)'),title('自由度2的时程曲线');
figure('numbertitle','off','name','自由度3的位移','pos',);
plot(t,disp(3,:)),grid,xlabel('时间(s/10)'),title('自由度3的时程曲线');
%end
图上的初始数据也是按照论坛上前辈所给数据
实例：
m=2*;
c=1.5*;
k=50*;
x0=';
v0=';
delt=0.1;
time=100;
作用力由于用到较多的数据，直接写在原程序中。
页:
[1]