I=0.05;
A=0.6;
L=2;
d=2500;
damage=zeros(10,1);
%%
n_number=11;
e_number=10;
idof=2;
T_dof=22;
gBC1=[1 1 0
11 1 0];
gNode=[0 0
2 0
4 0
6 0
8 0
10 0
12 0
14 0
16 0
18 0
20 0];
%%
%形成平面纯弯梁的单元刚度矩阵，单元质量矩阵
k=(EI/L^3)[12 6L -12 6L
6L 4L^2 -6L 2L^2
-12 -6L 12 -6L
6L 2L^2 -6L 4L^2];
m=(dAL/420)[156 22L 54 -13L^2
22L 4L^2 13L -3L^2
54 13L 156 -22L
-13L -3L^2 -22L 4*L^2];
K=zeros(T_dof);
M=zeros(T_dof);
%集成整体刚度和质量矩阵
for i=1:e_number
K=Beam2D2Node_Assembly(K,(1-damage(i))*k,i,i+1);
M=Beam2D2Node_massAssembly(M,m,i,i+1);
i=i+1
end
%% NewMark-beta法
% 1:频率 100HZ，步长 1/100，采样时间1s；
timestep=100;
gDeltaT=1/timestep;
tend=1;
%
Acce=zeros(n_numberidof,timestep);
Velo=zeros(n_numberidof,timestep);
Disp=zeros(n_number*idof,timestep);
Disp(:,1)=zeros(n_numberidof,1);
Velo(:,1)=ones(n_numberidof,1);
% 3:选择相关参数并计算下列有关常数
gama = 0.5 ;
beta = 0.25 ;
C = zeros( size( K ) ) ;
[N,N] = size( K ) ;
alpha0 = 1/beta/gDeltaT^2 ;
alpha1 = gama/beta/gDeltaT ;
alpha2 = 1/beta/gDeltaT ;
alpha3 = 1/2/beta - 1 ;
alpha4 = gama/beta - 1 ;
alpha5 = gDeltaT/2*(gama/beta-2) ;
alpha6 = gDeltaT*(1-gama) ;
alpha7 = gamagDeltaT ;
%计算有效刚度矩阵
K1 = K + alpha0M + alpha1*C;
%% 对K1进行边界条件处理
K1(1,:)=zeros(1,22);
K1(:,1)=zeros(22,1);
K1(21,:)=zeros(1,22);
K1(:,21)=zeros(22,1);
% step3.3 计算初始加速度
Acce(:,1) = M(-KDisp(:,1)-CVelo(:,1)) ;
% step3.4 对每一个时间步计算
for i=2：1：timestep
%计算t+delta_t时刻的有效荷载向量
f1 = M*(alpha0Disp(:,i-1)+alpha2Velo(:,i-1)+alpha3Acce(:,i-1)) ...
+ C(alpha1Disp(:,i-1)+alpha4Velo(:,i-1)+alpha5Acce(:,i-1)) ;
% 对f1进行边界条件处理
f1(1,1)=0;
f1(21,1)=0;
y = K1\f1 ;
%计算t+delta_t时刻的位移，加速度和速度
Disp(:,i) = K1\y ;
Acce(:,i) = alpha0(Disp(:,i)-Disp(:,i-1)) - alpha2Velo(:,i-1) - alpha3Acce(:,i-1) ;
Velo(:,i) = Velo(:,i-1) + alpha6Acce(:,i-1) + alpha7Acce(:,i) ;
end