〔求助〕如何用MATLAB计算结构响应？

ME!

用几种方法一起验证
%function [d,v,a] = Wilson( K, M, C, f, d1, v1, dt, tend )
% 利用Wilson-theta (Newmark) 解析解法计算结构的动力响应
  %K=sysK;% K ----- 调用刚度矩阵
%M=sysM; % M ----- 调用质量矩阵

% C=0.03*M+0.02*K;% C ----- 阻尼矩阵
K=[2,-1;-1,3];
  M=[4,0;0 8];
  C=0;
dt=0.1;% dt ----- 时间步长
tend=40;% tend --- 结束时间
w=600;
[n,n] = size( K ) ;

%Wilson算法
theta = 1.37 ;
tao = theta*dt ;
alpha0 = 6/tao^2 ;
alpha1 = 3/tao ;
alpha2 = 2*alpha1 ;
alpha3 = tao/2 ;
alpha4 = alpha0/theta ;
alpha5 = -alpha2/theta ;
alpha6 = 1-3/theta ;
alpha7 = dt/2 ;
alpha8 = dt^2/6 ;
K1 = K + alpha0*M + alpha1*C ;
d = zeros( n, floor(tend/dt) + 1 ) ;
v = zeros( n, floor(tend/dt) + 1 ) ;
a = zeros( n, floor(tend/dt) + 1 ) ;
f = zeros( n, floor(tend/dt) + 1 ) ;
d(1,1) = 1.2 ; %初始位移
d(2,1)=0;
v(:,1) = 0 ;%初始速度
f(:,1) =0 ;%初始载荷
a(:,1) = inv(M)*(f(:,1)-K*d(:,1)-C*v(:,1)) ;  %初始加速度
t=0:dt:tend;
for i=2:1:length(t)
%t(i) = (i-1)*dt ;
f(:,i)=0*100*sin(w*t(i));
ftheta = floor(theta) ;
%fq = f(i-1+ftheta-1)+ (theta-ftheta)*( f(i+ftheta-1) - f(i+ftheta-2) ) ;
fq=f(:,i-1)+ftheta*(f(:,i)-f(:,i-1));
f1 = fq + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+2*a(:,i-1))+ C*(alpha1*d(:,i-1)+2*v(:,i-1)+alpha3*a(:,i-1)) ;
dq= inv(K1)*f1 ;
a(:,i) = alpha4*(dq-d(:,i-1)) + alpha5*v(:,i-1) + alpha6*a(:,i-1) ;
v(:,i) = v(:,i-1) + alpha7 * ( a(:,i) + a(:,i-1) ) ;
d(:,i) = d(:,i-1) + dt*v(:,i-1) + alpha8 * ( a(:,i)+2*a(:,i-1) ) ;
end
%解析解

for i=1:length(t)
x(i)=0.4*cos(0.5*t(i))+0.8*cos(1.581*0.5*t(i));
end

%Newmark算法
gama = 0.5 ;
beta = 0.25 ;
[n,n] = size( K ) ;
Nalpha0 = 1/beta/dt^2 ;
Nalpha1 = gama/beta/dt ;
Nalpha2 = 1/beta/dt ;
Nalpha3 = 1/2/beta - 1 ;
Nalpha4 = gama/beta - 1 ;
Nalpha5 = dt/2*(gama/beta-2) ;
Nalpha6 = dt*(1-gama) ;
Nalpha7 = gama*dt ;
NK1 = K + Nalpha0*M + Nalpha1*C ;
Nd = zeros( n, floor(tend/dt) + 1 ) ;
Nv = zeros( n, floor(tend/dt) + 1 ) ;
Na = zeros( n, floor(tend/dt) + 1 ) ;
Nd(1,1) = 1.2 ; %初始位
Nd(2,1)=0;
f(:,1) =0 ;%初始载荷
Na(:,1) = inv(M)*(f(:,1)-K*Nd(:,1)-C*Nv(:,1)) ;  %初始加速度
t=0:dt:tend;
for i=2:1:length(t)
f(:,i)=0*100*sin(w*t(i));
f2 = f(:,i) + M*(Nalpha0*Nd(:,i-1)+Nalpha2*Nv(:,i-1)+Nalpha3*Na(:,i-1))+ C*(Nalpha1*Nd(:,i-1)+Nalpha4*Nv(:,i-1)+Nalpha5*Na(:,i-1)) ;
Nd(:,i) = inv(NK1)*f2 ;
Na(:,i) = Nalpha0*(Nd(:,i)-Nd(:,i-1)) - Nalpha2*Nv(:,i-1) - Nalpha3*Na(:,i-1) ;
Nv(:,i) = Nv(:,i-1) + Nalpha6*Na(:,i-1) + Nalpha7*Na(:,i) ;
end

plot(t,d(1,:),'-b^',t,Nd(1,:),'g+',t,x,'r')
xlabel('t');
ylabel('u(t)');
title('位移与时间的关系');

听风声解花语

feifeifool 发表于 2006-8-14 09:15
clc;clear;close;
M=[1 1 1];M=diag(M);
C=diag([0.1,0.0,0.0]);

如果是个移动荷载怎么处理啊