- 积分
- 0
- 注册时间
- 2011-9-9
- 仿真币
-
- 最后登录
- 1970-1-1
|
%--------------------------------------------------------
% Ex.2.2 plane 9-bar truss
%--------------------------------------------------------
% To compute the dynamic response of node 4 on the vertical direction
%
%
E=2.0e11;
A=2.5e-3;
density=7860;
node_number=6;
elment_number=9;
nc=[0 0;4 0;4 3;8 0;8 3;12 0]; %node_coordinate
en=[1 2;1 3;2 3;2 4;3 4;3 5;4 5;4 6;5 6]; % element_node
ed(1:node_number,1:2)=1; % elment_displacement
dof=0;
for loopi=1:node_number
for loopj=1:2
dof=dof+1;
ed(loopi,loopj)=dof;
end
end
ek=E*A*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; % elment_stiffness_matrix
em=(density*A)/6*[2 0 1 0; ... %elment_mass_matrix
0 2 0 1; ...
1 0 2 0; ...
0 1 0 2];
k(1:dof,1:dof)=0; %tructural_stiffness_matrix
m=k; % structural_mass_matrix, same size with k
theta(1:9)=0;
el(1:9)=0;
e2s(1:4)=0; %index of transform the elment displament number to structural
for loopi=1:elment_number
for zi=1:2
e2s((zi-1)*2+1)=ed(en(loopi,zi),1);
e2s((zi-1)*2+2)=ed(en(loopi,zi),2);
end
el(loopi)=sqrt((nc(en(loopi,1),1)-nc(en(loopi,2),1))^2+(nc(en(loopi,1),2)-nc(en(loopi,2),2))^2);
theta(loopi)=asin((nc(en(loopi,1),2)-nc(en(loopi,2),2))/el(loopi));
lmd=[cos(theta(loopi)) sin(theta(loopi)); -sin(theta(loopi)) cos(theta(loopi))];
t=[lmd zeros(2); zeros(2) lmd];
dk=t'*ek*t/el(loopi);
dm=t'*em*t*el(loopi);
for jx=1:4
for jy=1:4
k(e2s(jx),e2s(jy))=k(e2s(jx),e2s(jy))+dk(jx,jy);
m(e2s(jx),e2s(jy))=m(e2s(jx),e2s(jy))+dm(jx,jy);
end
end
end
c=0*k+0*m;
nt=1500;dt=0.0001;
time=0:dt:nt*dt;
q0=zeros(dof,1);
dq0=zeros(dof,1);
bcdof=zeros(dof,1);
fd=zeros(dof,nt);
for i=1:nt
fd(6,i)=200;
end
[acc,vel,dsp]=TransResp4(k,c,m,fd,bcdof,nt,dt,q0,dq0);
plot(time, dsp(8,:))
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
——————————————————————————————————————————————
用到的函数 wilson
——————————————————————————————————————————————
function [acc,vel,dsp]=TransResp4(kk,cc,mm,fd,bcdof,nt,dt,q0,dq0)
%--------------------------------------------------------------------------
% Purpose:
% The function subroutine TransResp4.m calculates transient response of
% a structural system using Wilson integration scheme.
% Synopsis:
% [acc, vel, dsp]=TransResp4(kk,cc,mm,fd,bcdof,nt,dt,q0,dq0)
% Variable Description:
% Input parameters
% kk, cc, mm - stiffness, damping and mass matrices
% fd - Input or forcing influence matrix
% bcdof - Boundary condition dofs vector
% nt - Number of time steps
% dt - Time step size
% q0, dq0 - Initial condition vectors
% Output parameters
% acc - Acceleration response
% vel - Velocity response
% dsp - Displacement response
%--------------------------------------------------------------------------
% (1) initial condition
%--------------------------------------------------------------------------
[sdof,n2]=size(kk);
dsp=zeros(sdof,nt); % displacement matrix
vel=zeros(sdof,nt); % velocity matrix
acc=zeros(sdof,nt); % acceleration matrix
dsp(:,1)=q0; % initial displacement
vel(:,1)=dq0; % initial velocity
theta=1.4; % select the parameters
%--------------------------------------------------------------------------
% (2) Wilson integration scheme for time integration
%--------------------------------------------------------------------------
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1));
% compute the initial acceleration (t=0)
ekk=kk+mm*(6/(theta*dt)^2)+cc*(3/(theta*dt));
% compute the effective stiffness matrix
for i=1:sdof % assign zero to dsp, vel, acc of the dofs associated with bc
if bcdof(i)==1
dsp(i,1)=0;
vel(i,1)=0;
acc(i,1)=0;
end
end
for it=1:nt % loop for each time step
cfm=dsp(:,it)*(6/(theta*dt)^2)+vel(:,it)*(6/(theta*dt))+2*acc(:,it);
cfc=dsp(:,it)*(3/(theta*dt))+2*vel(:,it)+acc(:,it)*(theta*dt/2);
efd=fd(:,it)+theta*(fd(:,it+1)-fd(:,it))+mm*cfm+cc*cfc;
% compute the effective force vector
dtheta=inv(ekk)*efd; % find the displacement at time t+ dt
acc(:,it+1)=(dtheta-dsp(:,it))*(6/(theta^3*dt^2))...
-vel(:,it)*(6/(theta^2*dt))+acc(:,it)*(1-3/theta);
% find the acceleration at time t+dt
vel(:,it+1)=vel(:,it)+acc(:,it+1)*dt/2+acc(:,it)*dt/2;
% find the velocity at time t+dt
dsp(:,it+1)=dsp(:,it)+vel(:,it)*dt...
+(acc(:,it+1)+2*acc(:,it))*(dt^2/6);
% find the displacement at time t+dt
for i=1:sdof % assign zero to acc, vel, dsp of the dofs associated with bc
if bcdof(i)==1
dsp(i,it+1)=0;
vel(i,it+1)=0;
acc(i,it+1)=0;
end
end
end
if cc(1,1)==0
disp('The transient response results of undamping system')
else
disp('The transient response results of damping system')
end
disp('The method is Wilson integration')
%--------------------------------------------------------------------------
% The end
%--------------------------------------------------------------------------
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
|