动力学系统时域响应计算的六种方法Matlab源程序(Newmark,Houbolt法,中心差分法)
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
vel(:,1)=dq0; % initial velocity
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1));
dsp0=dsp(:,1)-vel(:,1)*dt+0.5*acc(:,1)*dt^2;
dsp(:,1)=q0; % initial displacement
vel(:,1)=dq0; % initial velocity
dsp(:,it+1)=inv(ekk)*efd; % find the displacement at time t+dt
acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...
alpha=0.5; beta=0.5; % select the parameters
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1)); % compute the initial acceleration (t=0)
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 vector
-acc(:,it)*(0.5/alpha-1); % find the acceleration at time t+dt
vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt;% find the velocity at time t+dt
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 matrix
vel=zeros(sdof,nt); % velocity matrix
acc=zeros(sdof,nt); % acceleration matrix
ekk=mm/dt^2+cc/(2*dt);
efd=fd(:,1)-(kk-2*mm/dt^2)*dsp(:,1)-(mm/dt^2-cc/(2*dt))*dsp0;
% mm------质量矩阵
% cc------阻尼矩阵
% q0------初始位移
% dq0------初始速度
% dt------时间步长
% nt------总的计算步数，等于结束时间除以dt
%返回值
% dsp------位移
% vel------速度
% acc------加速度
acc=zeros(sdof,nt); % acceleration matrix
dsp(:,1)=q0; % initial displacement
end
中心差分法Matlab源程序
function [acc,vel,dsp]=central_diff(kk,cc,mm,fd,nt,dt,q0,dq0)
%输入参数
% kk------刚度矩阵
[sdof,n2]=size(kk);
dsp=zeros(sdof,nt); % displacement matrix
vel=zeros(sdof,nt); % velocity matrix
ekk=kk+mm/(alpha*dt^2)+cc*beta% loop for each time step
cfm=dsp(:,it)/(alpha*dt^2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1));
dsp0=dsp(:,1)-vel(:,1)*dt+0.5*acc(:,1)*dt^2;
dsp(:,1)=q0; % initial displacement
vel(:,1)=dq0; % initial velocity
dsp(:,it+1)=inv(ekk)*efd; % find the displacement at time t+dt
acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...
alpha=0.5; beta=0.5; % select the parameters
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1)); % compute the initial acceleration (t=0)
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 vector
-acc(:,it)*(0.5/alpha-1); % find the acceleration at time t+dt
vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt;% find the velocity at time t+dt
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 matrix
vel=zeros(sdof,nt); % velocity matrix
acc=zeros(sdof,nt); % acceleration matrix
ekk=mm/dt^2+cc/(2*dt);
efd=fd(:,1)-(kk-2*mm/dt^2)*dsp(:,1)-(mm/dt^2-cc/(2*dt))*dsp0;
% mm------质量矩阵
% cc------阻尼矩阵
% q0------初始位移
% dq0------初始速度
% dt------时间步长
% nt------总的计算步数，等于结束时间除以dt
%返回值
% dsp------位移
% vel------速度
% acc------加速度
acc=zeros(sdof,nt); % acceleration matrix
dsp(:,1)=q0; % initial displacement
end
中心差分法Matlab源程序
function [acc,vel,dsp]=central_diff(kk,cc,mm,fd,nt,dt,q0,dq0)
%输入参数
% kk------刚度矩阵
[sdof,n2]=size(kk);
dsp=zeros(sdof,nt); % displacement matrix
vel=zeros(sdof,nt); % velocity matrix
ekk=kk+mm/(alpha*dt^2)+cc*beta% loop for each time step
cfm=dsp(:,it)/(alpha*dt^2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);