%   Newmark-Beta法求解動力平衡方程   文獻本P31
%目的：用于檢查ANSYS軟件算出的Y、YY和YYY是否正確。P113算例
%************************************************
clc;clear;
%--------------------------------
%  輸入結構基本數據,單位：KN,m
%--------------------------------
NN=31;   %節點數
NE=30;   %單元數
L=4;  %單元長度
ND=2*NN;    %節點位移總數
%--------------------------------
%     單元剛度矩陣ki,單元質量矩陣mi
%--------------------------------
ki=8.04e9*EleSiff(L);    %     i=1,2,,NE=30
mi=9878*EleMass(L);    %     i=1,2,,NE=30
%--------------------------------
%    荷載Pt
%--------------------------------
t=0:0.01:1.99;
randn('state',2)
Pt=-1.12e5*randn(size(t));                    %隨機干擾
Pt=[zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    Pt
    zeros(size(t))
    Pt
    zeros(size(t))
    Pt
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))
    zeros(size(t))];
%--------------------------------
%    K,M,C
%--------------------------------
K=zeros(ND);
K=StiffAssem(K,ki,1,2);    %   62*62
K=StiffAssem(K,ki,2,3);
K=StiffAssem(K,ki,3,4);
K=StiffAssem(K,ki,4,5);
K=StiffAssem(K,ki,5,6);
K=StiffAssem(K,ki,6,7);
K=StiffAssem(K,ki,7,8);
K=StiffAssem(K,ki,8,9);
K=StiffAssem(K,ki,9,10);
K=StiffAssem(K,ki,10,11);
K=StiffAssem(K,ki,11,12);
K=StiffAssem(K,ki,12,13);
K=StiffAssem(K,ki,13,14);
K=StiffAssem(K,ki,14,15);
K=StiffAssem(K,ki,15,16);
K=StiffAssem(K,ki,16,17);
K=StiffAssem(K,ki,17,18);
K=StiffAssem(K,ki,18,19);
K=StiffAssem(K,ki,19,20);
K=StiffAssem(K,ki,20,21);
K=StiffAssem(K,ki,21,22);
K=StiffAssem(K,ki,22,23);
K=StiffAssem(K,ki,23,24);
K=StiffAssem(K,ki,24,25);
K=StiffAssem(K,ki,25,26);
K=StiffAssem(K,ki,26,27);
K=StiffAssem(K,ki,27,28);
K=StiffAssem(K,ki,28,29);
K=StiffAssem(K,ki,29,30);
K=StiffAssem(K,ki,30,31);
K=[K(2:20,2:20)  K(2:20,22:40)  K(2:20,42:60) K(2:20,62)%引入邊界條件，縮減自由度數,58*58
   K(22:40,2:20)  K(22:40,22:40) K(22:40,42:60) K(22:40,62)
   K(42:60,2:20)  K(42:60,22:40) K(42:60,42:60) K(42:60,62)
   K(62,2:20)  K(62,22:40) K(62,42:60) K(62,62)];
M=zeros(ND);
M=MassAssem(M,mi,1,2);    %   62*62
M=MassAssem(M,mi,2,3);
M=MassAssem(M,mi,3,4);
M=MassAssem(M,mi,4,5);
M=MassAssem(M,mi,5,6);
M=MassAssem(M,mi,6,7);
M=MassAssem(M,mi,7,8);
M=MassAssem(M,mi,8,9);
M=MassAssem(M,mi,9,10);
M=MassAssem(M,mi,10,11);
M=MassAssem(M,mi,11,12);
M=MassAssem(M,mi,12,13);
M=MassAssem(M,mi,13,14);
M=MassAssem(M,mi,14,15);
M=MassAssem(M,mi,15,16);
M=MassAssem(M,mi,16,17);
M=MassAssem(M,mi,17,18);
M=MassAssem(M,mi,18,19);
M=MassAssem(M,mi,19,20);
M=MassAssem(M,mi,20,21);
M=MassAssem(M,mi,21,22);
M=MassAssem(M,mi,22,23);
M=MassAssem(M,mi,23,24);
M=MassAssem(M,mi,24,25);
M=MassAssem(M,mi,25,26);
M=MassAssem(M,mi,26,27);
M=MassAssem(M,mi,27,28);
M=MassAssem(M,mi,28,29);
M=MassAssem(M,mi,29,30);
M=MassAssem(M,mi,30,31);
M=[M(2:20,2:20)  M(2:20,22:40)   M(2:20,42:60)  M(2:20,62)%引入邊界條件，縮減自由度數,58*58
   M(22:40,2:20)  M(22:40,22:40) M(22:40,42:60) M(22:40,62)
   M(42:60,2:20)  M(42:60,22:40) M(42:60,42:60) M(42:60,62)
   M(62,2:20)     M(62,22:40)     M(62,42:60)    M(62,62)];
C=0.6247*M+0.0039*K;
%--------------------------------
%    Newmark-Beta法
%--------------------------------
DeltT=0.01;
Beta=0.25;
Gama=0.5;
a0=1/(Beta*DeltT^2);
a1=Gama/(Beta*DeltT);
a2=1/(Beta*DeltT);
a3=1/(2*Beta)-1;
a4=Gama/Beta-1;
a5=(Gama/(2*Beta)-1)*DeltT;
KK=K+a0*M+a1*C;                  %等效剛度矩陣
Npoint=200;
Y=zeros(58,Npoint);YY=zeros(58,Npoint); YYY=zeros(58,Npoint);   %賦初值
Y(:,1)=0;YY(:,1)=0;YYY(:,1)=M\Pt(:,1); %初始值
for i=1:199
    PPt(:,i+1)=Pt(:,i+1)+M*(a0*Y(:,i)+a2*YY(:,i)+a3*YYY(:,i))+C*(a1*Y(:,i)+a4*YY(:,i)+a5*YYY(:,i));       %t+DeltaT 時刻的等效荷載 
    Y(:,i+1)=KK\PPt(:,i+1);  %位移
    YYY(:,i+1)=a0*(Y(:,i+1)-Y(:,i)-DeltT*YY(:,i)-(0.5-Beta)*DeltT*DeltT*YYY(:,i)); %加速度
    YY(:,i+1)=YY(:,i)+(1-Gama)*DeltT*YYY(:,i)+Gama*DeltT*YYY(:,i+1); %速度
end
for i=1 :200     %檢查結果是否正確，若xxx=yyy，則正確。經檢查證明，本程序正確。
xxx(:,i)=M*YYY(:,i)+C*YY(:,i)+K*Y(:,i) ;
yyy(:,i)=Pt(:,i);
end