动力学系统时域响应计算的六种方法Matlab源程序.docx
动力学系统时域响应计算的六种方法Matlab源程序Newmark法Matlab源程序function acc,vel,dsp=Newmark_2(kk,cc,mm,fd,nt,dt,q0,dq0)%输入参数% kk-刚度矩阵% mm-质量矩阵% cc-阻尼矩阵% q0-初始位移% dq0-初始速度% dt-时间步长% nt-总的计算步数，等于结束时间除以dt%返回值% dsp-位移% vel-速度% acc-加速度sdof,n2=size(kk);dsp=zeros(sdof,nt); % displacement matrixvel=zeros(sdof,nt); % velocity matrixacc=zeros(sdof,nt); % acceleration matrixdsp(:,1)=q0; % initial displacementvel(:,1)=dq0; % initial velocityalpha=0.5; beta=0.5; % select the parametersacc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1); % compute the initial acceleration (t=0)ekk=kk+mm/(alpha*dt2)+cc*beta/(alpha*dt);for it=1:nt % loop for each time stepcfm=dsp(:,it)/(alpha*dt2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);cfc=dsp(:,it)*beta/(alpha*dt)+vel(:,it)*(beta/alpha-1).+acc(:,it)*(0.5*beta/alpha-1)*dt;efd=fd(:,it)+mm*cfm+cc*cfc; % compute the effective force vectordsp(:,it+1)=inv(ekk)*efd; % find the displacement at time t+dtacc(:,it+1)=(dsp(:,it+1)-dsp(:,it)/(alpha*dt2)-vel(:,it)/(alpha*dt).-acc(:,it)*(0.5/alpha-1); % find the acceleration at time t+dtvel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt;% find the velocity at time t+dtend中心差分法Matlab源程序function acc,vel,dsp=central_diff(kk,cc,mm,fd,nt,dt,q0,dq0)%输入参数% kk-刚度矩阵% mm-质量矩阵% cc-阻尼矩阵% q0-初始位移% dq0-初始速度% dt-时间步长% nt-总的计算步数，等于结束时间除以dt%返回值% dsp-位移% vel-速度% acc-加速度sdof,n2=size(kk);dsp=zeros(sdof,nt); % displacement matrixvel=zeros(sdof,nt); % velocity matrixacc=zeros(sdof,nt); % acceleration matrixdsp(:,1)=q0; % initial displacementvel(:,1)=dq0; % initial velocityacc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1);dsp0=dsp(:,1)-vel(:,1)*dt+0.5*acc(:,1)*dt2;ekk=mm/dt2+cc/(2*dt);efd=fd(:,1)-(kk-2*mm/dt2)*dsp(:,1)-(mm/dt2-cc/(2*dt)*dsp0;dsp(:,1+1)=inv(ekk)*efd;for it=2:nt % loop for each time step after first stepefd=fd(:,it)-(kk-2*mm/dt2)*dsp(:,it)-(mm/dt2-cc/(2*dt)*dsp(:,it-1);% compute the effective force vectordsp(:,it+1)=inv(ekk)*efd; % find the dsp at t+dtacc(:,it)=(dsp(:,it+1)-2*dsp(:,it)+dsp(:,it-1)/dt2; % find the acc at tvel(:,it)=(dsp(:,it+1)-dsp(:,it-1)/(2*dt); % find the vel at tacc(:,it+1)=acc(:,it); vel(:,it+1)=vel(:,it);endHoubolt法Matlab源程序function acc,vel,dsp=Houbolt(kk,cc,mm,fd,nt,dt,q0,dq0)%输入参数% kk-刚度矩阵% mm-质量矩阵% cc-阻尼矩阵% q0-初始位移% dq0-初始速度% dt-时间步长% nt-总的计算步数，等于结束时间除以dt%返回值% dsp-位移% vel-速度% acc-加速度sdof,n2=size(kk);dsp=zeros(sdof,nt); % displacement matrixvel=zeros(sdof,nt); % velocity matrixacc=zeros(sdof,nt); % acceleration matrixdsp(:,1)=q0; % initial displacementvel(:,1)=dq0; % initial velocityacc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1);dsp0=dsp(:,1)-vel(:,1)*dt+0.5*acc(:,1)*dt2;ekk=mm/dt2+cc/(2*dt); % compute the effective stiffness matrixit=1;efd=fd(:,it)-(kk-2*mm/dt2)*dsp(:,it)-(mm/dt2-cc/(2*dt)*dsp0;dsp(:,it+1)=inv(ekk)*efd; % find the dsp at t+dtfor it=2:3 % loop for two steps after first stepefd=fd(:,it)-(kk-2*mm/dt2)*dsp(:,it)-(mm/dt2-cc/(2*dt)*dsp(:,it-1);dsp(:,it+1)=inv(ekk)*efd; % find the dsp at t+dtacc(:,it)=(dsp(:,it+1)-2*dsp(:,it)+dsp(:,it-1)/dt2; % find the acc at tvel(:,it)=(dsp(:,it+1)-dsp(:,it-1)/(2*dt); % find the vel at tendekk=kk+cc*11/(6*dt)+mm*2/(dt*dt);for it=3:nt-1 % loop for each time step after initial stepcfm=dsp(:,it)*5/dt2-dsp(:,it-1)*4/dt2+dsp(:,it-2)*1/dt2;cfc=dsp(:,it)*3/dt-dsp(:,it-1)*3/(2*dt)+dsp(:,it-2)*1/(3*dt);efd=fd(:,it+1)+mm*cfm+cc*cfc;dsp(:,it+1)=inv(ekk)*efd; % find the dsp at t+dtacc(:,it)=(2*dsp(:,it+1)-5*dsp(:,it)+4*dsp(:,it-1)+dsp(:,it-2)/dt2;vel(:,it)=(11*dsp(:,it+1)-18*dsp(:,it)+9&dsp(:,it-1)-2*dsp(:,it-2)/(6*dt);acc(:,it+1)=acc(:,it); vel(:,it+1)=vel(:,it);end
