声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1577|回复: 1

[编程技巧] 线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书

[复制链接]
发表于 2012-9-3 10:34 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
此类帖子之前有人发过,不过我拿到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)=[M+delt*C/2+delt^2*K/6]^(-1){F(t+delt)-C[X(t)+delt*X(t)'']-K[X(t)+delt*X(t)'+delt^2*X(t)''/3]}
% X(t+delt)'=X(t)'+delt*[X(t)''+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);
[mod,fre]=eig(inv(m)*k);
diag(sqrt(fre));%固有频率
i=1;
beta=1/4;
for t=0:delt:time
    f=[2.0*sin(3.5*t) -2.0*cos(2*t) 1.0*sin(3*t)]';%外扰力
    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',[450 180 400 420]);
plot(t,disp(1,:)),grid,xlabel('时间(s/10)'),title('自由度1的时程曲线');
figure('numbertitle','off','name','自由度2的位移','pos',[450 180 400 420]);
plot(t,disp(2,:)),grid,xlabel('时间(s/10)'),title('自由度2的时程曲线');
figure('numbertitle','off','name','自由度3的位移','pos',[450 180 400 420]);
plot(t,disp(3,:)),grid,xlabel('时间(s/10)'),title('自由度3的时程曲线');
%end

三自由度简谐受迫振动

三自由度简谐受迫振动
回复
分享到:

使用道具 举报

 楼主| 发表于 2012-9-3 10:38 | 显示全部楼层
图上的初始数据也是按照论坛上前辈所给数据
实例:
m=2*[1 0 0;0 1 0;0 0 1];
c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
k=50*[2 -1 0;-1 2 -2;0 -1 2];
x0=[1 1 1]';
v0=[1 1 1]';
delt=0.1;
time=100;
作用力由于用到较多的数据,直接写在原程序中。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-28 23:27 , Processed in 0.071346 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表