M=[2,0;0 1 ];                      %质量矩阵

K=[6 -2;-2 4];                   %刚度矩阵

a=0;b=0;

C=a*K+b*M;

dt=0.28;

t=0:dt:2.8;

ft0=zeros(length(K),length(t));

for i=1:length(t)

ft0(1,i)=10;  %在节点4的竖直方向加大小为200N的阶跃力

end

dsp=zeros(length(K),length(t));                                         % 位移

vel=zeros(length(K),length(t));                                             % 速度

acc=zeros(length(K),length(t));                                          % 加速度

%--------------------------------------------------------------------------

%  (2) Newmark

alpha=0.3; beta=0.6;                                        % 稳定条件

acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1));   % 计算初始加速度 (t=0)

ekk=K+M/(alpha*dt^2)+C*beta/(alpha*dt);  % 计算有效刚度矩阵

for it=2:length(t)-1                                              % 时间步循环

cfm=dsp(:,it)/(alpha*dt^2)+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=ft0(:,it)+M*cfm+C*cfc;                  %  计算有效力矢量

dsp(:,it+1)=inv(ekk)*efd;                        %  t+dt时刻的位移

acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...

-acc(:,it)*(0.5/alpha-1);              %  t+dt时刻的加速度

vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt; %  t+dt时刻的速度

end

plot(t, dsp(2,:),'-b*')

hold on

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')

%  (2) Wilson

%--------------------------------------------------------------------------

theta=1.4;                                          % 稳定条件参数

acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1)); % 计算初始加速度 (t=0)

ekk=K+M*(6/(theta*dt)^2)+C*(3/(theta*dt));  % 计算有效刚度矩阵

for it=2:length(t)-1

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=ft0(:,it)+theta*(ft0(:,it+1)-ft0(:,it))+M*cfm+C*cfc; %计算有效力矢量

dtheta=inv(ekk)*efd;                         %  t+ dt时刻的位移

acc(:,it+1)=(dtheta-dsp(:,it))*(6/(theta^3*dt^2))...

-vel(:,it)*(6/(theta^2*dt))+acc(:,it)*(1-3/theta);  %  t+dt时刻的加速度

vel(:,it+1)=vel(:,it)+acc(:,it+1)*dt/2+acc(:,it)*dt/2; %  t+dt时刻的速度

dsp(:,it+1)=dsp(:,it)+vel(:,it)*dt...

+(acc(:,it+1)+2*acc(:,it))*(dt^2/6);         %  t+dt时刻的位移

end

plot(t, dsp(2,:),'-g*')

hold on

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')

%--------------------------------------------------------------------------

%  (2) 振型叠加法

%施加不平衡激振力

t=t';

[V,D]=eig(K,M);                       % 计算特征值和特征向量

[lambda,ki]=sort(diag(D));                     % 给特征值和特征向量排序

omega=sqrt(lambda);                                       % 角频率w.

omega1=sqrt(lambda)/(2*pi);                           % 固有频率 Hz.

V1=V(:,ki);

h=V1'*M*V1;

Factor=diag(V1'*M*V1);                             %模态质量

Vnorm=V1*inv(sqrt(diag(Factor)));                  % 正则化振型向量

VnormK=diag(Vnorm'*K*Vnorm);                       %模态刚度

omega0=diag(sqrt(Vnorm'*K*Vnorm)) ;                % 正则化模态刚度

VnormM=diag(Vnorm'*M*Vnorm);                       %正则化模态质量

Fnorm=Vnorm'*ft0 ;                                     % 模态力矢量

Modamp=Vnorm'*(a*M+b*K)*Vnorm;                     % 模态阻尼矩阵

zeta=diag((1/2)*Modamp*inv(diag(omega0)));           % 阻尼比

%--------------------------------------------------------------------------

%  %  (4) 模态坐标的响应

%--------------------------------------------------------------------------

q0=zeros(length(K),1);

dq0=zeros(length(K),1);                         % 初始化位移和速度

eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0;     %2.2-53% 初始条件模态坐标的位移和速度

eta=zeros(length(t),length(K));                %  %初始化位移2.2-52

for i=1:length(K)                                 % t(i)时刻的响应

phase0=omega0(i)*t;                                %w*t omega0为无阻尼固有频率

omegad=omega(i)*sqrt(1-zeta(i)^2);

phase=omegad*t;                            %wd*t其中(omegad为有阻尼固有频率)

Exx=exp(-zeta(i)*omega(i)*t);                %中间变量

C1=eta0(i);                                        %初始位移

C2=(deta0(i)+eta0(i)*zeta(i)*omega0(i))/omegad;          %中间变量

D1=zeta(i)*omega(i)/omegad;

II=ones(length(t),1);

XX=Fnorm(i)/(omegad^2+zeta(i)^2*omega(i)^2);

eta(:,i)=C1*Exx.*cos(phase)+C2*Exx.*sin(phase)+ XX*(II-Exx.*cos(phase)-D1*Exx.*sin(phase));%有阻尼系统

%  eta(:,i)=Fnorm(i)/((omega(i))^2)*(1-cos(omega(i)*t)); %无阻尼系统

end

%--------------------------------------------------------------------------

%   (6) 将模态坐标转化到物理坐标系

%--------------------------------------------------------------------------

eta=eta';            %模态坐标

y=Vnorm*eta;           %物理坐标

plot(t, y(2,:),'-r*')

xlabel('Time(seconds)')

ylabel(' Vertical displ. (m)')