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)')