program1.m

program about base isolation structure

% program about base isolation structure
% 行末注有@表明这行参数在不同的算例中有可能会变化
XZ=0;                                                % 陡线中点横坐标值
XS=0.02;                                             % 隔震层屈服位移                @
cs1=0.05;                                            % 上部结构第一振型阻尼比        @
cs2=0.05;                                            % 上部结构第二振型阻尼比        @
k1=1.5;                                              % 隔震层弹性刚度                @
k2=0.6;                                              % 隔震层塑性刚度                @
ug=load('d:\data\input\wave.txt','-ascii');          % 输入地震波
ug=ug';                                              % 转置
ug=ug(:);                                            % 转为单下标标识
n1=length(ug);                                       % 输入地震波记录点数
n1=2000;                                              % 计算点数
ugm=341.7/max(abs(ug));                              % 调整加速度峰值
ug=0.01*ug*ugm;                                      % 单位转化（1厘米/秒*秒＝0.01米/秒*秒）
m0=[41 53.2 61.3 54.5 54.5 54.5 32.4]*1e+3;          % 输入质量阵 %%%@  黄建文
n=length(m0);                                        % 质点数
M=diag(m0);                                          % 化为对角阵
P=zeros(n,1);                                        % 为各层的恢复力列向量(方程右边项）
p(1)=0;                                              % 为隔震层刚度判别因子
k0=[k1 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;        % 输入刚度阵                   @
[K]=matrix(k0,n);                                    % 调用子程序求刚度阵
[x,d]=eig(K,M);                                      % 求结构的自振特性
d=sqrt(d);%
w=sort(diag(d));                                     % 结构的自振频率
T=2*3.1415926./w;                                    % 结构的自振周期
c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);  % 上部结构瑞雷阻尼系数
c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);            % 上部结构瑞雷阻尼系数
%C=c1*M+c2*K;                                        % 阻尼阵
c0=[43.5 171.5 171.5 111.4 111.4 111.4 111.4]*1e+3;  % 输入各层阻尼值
[C]=matrix(c0,n)                                     % 调用子程序求阻尼阵
dt=0.005;                                             % 时间步长                      @
ug0=[ug(1)*ones(n,1)];%
F0=-M*ug0;                                           % 计算初始力
X0=zeros(n,1);                                       % 计算初始位移
X1=zeros(n,1);                                       % 计算初始速度
X2=inv(M)*(F0-C*X1-K*X0);                            % 计算初始加速度
XX0(1:n,1)=X0;%
XX1(1:n,1)=X1;%
XX2(1:n,1)=X2;%

%计算wilson法的几个系数
sita=1.37;%
a0=6/(sita*dt)^2;%
a1=3/(sita*dt);%
a2=2*a1;                  
a3=sita*dt/2;%
a4=a0/sita;%
a5=-a2/sita;%
a6=1-3/sita;%
a7=dt/2;%
a8=dt*dt/6;%
Ki=K+a0*M+a1*C;                                      % 形成有效刚度矩阵



% 对每一时间步长作如下计算
i=2;
while(i<n1+1);
    ug0=[ug(i)*ones(n,1)];                           % 地震加速度
    F_t=-M*ug0-P;%
    Ft0=F0+sita*(F_t-F0)+M*(a0*X0+a2*X1+2*X2)+C*(a1*X0+2*X1+a3*X2);%
    Xt0=inv(Ki)*Ft0;                                 % (t+sita*dt)的位移
    
    
    % 计算t+dt的加速度，速度，位移
    X_2=X2+(a0*Xt0-a0*X0-2*a1*X1-3*X2)/sita;          %计算t+dt的加速度
    X_1=X1+a7*(X_2+X2);                               %速度
    X_0=X0+dt*X1+a8*(X_2+2*X2);                       %位移
    X_2=-inv(M)*C*X_1-inv(M)*K*X_0-ug(i)-inv(M)*P;    %直接由动力方程计算加速度，消除累积误差
    
    

       %%%%%%%%%% 求出刚度变化判断因子p %%%%%%%%%%
    if (X_0(1)-XZ)>=XS&X_1(1)>0;%
        p(i)=1;                                       % 上缓线
        P(1)=(k1-k2)*(1e+6)*XS;                       % 隔震层的恢复力列向量(方程右边项）    @
        Q(i)=k2*(1e+6)*X_0(1)+P(1);                   % 隔震层的恢复力列向量                 @
    elseif (X_0(1)-XZ)<=-XS&X_1(1)<0%
        p(i)=-1;                                      % 下缓线
        P(1)=(k2-k1)*(1e+6)*XS;                       % 隔震层的恢复力列向量(方程右边项）    @
        Q(i)=k2*(1e+6)*X_0(1)+P(1) ;                  % 隔震层的恢复力列向量                 @
    else p(i)=0;                                      % 陡线
        P(1)=(k2-k1)*(1e+6)*XZ;%                      % 隔震层的恢复力列向量(方程右边项）    @
        Q(i)=k1*(1e+6)*X_0(1)+P(1);                   % 隔震层的恢复力列向量                 @
      end
   
      
    
     %%%%%%%%%%   隔震层第二刚度赋值   %%%%%%%%%%
    if (p(i)==1)|(p(i)==-1);                                %
        XZ=X_0(1)-p(i)*XS;                                  %
        k0=[k2 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;      % 输入刚度阵            @
        [K]=matrix(k0,n);
        [x,d]=eig(K,M);%
        d=sqrt(d);%
        w=sort(diag(d));
        c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);  % 上部结构瑞雷阻尼系数
        c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);            % 上部结构瑞雷阻尼系数
        %C=c1*M+c2*K;                                        % 阻尼阵
        Ki=K+a0*M+a1*C;                                      % 形成有效刚度矩阵
    end
    
     if (p(i)==0);                                  % 
        k0=[k1 115.3 115.3 74.9 74.9 74.9 74.6]*1e+6;         % 输入刚度阵            @
        [K]=matrix(k0,n);
        [x,d]=eig(K,M);                                        %
        d=sqrt(d);                                             %
        c1=2*w(1)*w(2)*(cs1*w(2)-cs2*w(1))/(w(2)^2-w(1)^2);    % 上部结构瑞雷阻尼系数
        c2=2*(cs2*w(2)-cs1*w(1))/(w(2)^2-w(1)^2);              % 上部结构瑞雷阻尼系数
        %C=c1*M+c2*K;                                          % 阻尼阵
        Ki=K+a0*M+a1*C;                                        % 形成有效刚度矩阵
    end;
 %变量传递
    X0=X_0;%
    X1=X_1;%
    X2=X_2;%
    F0=-M*ug0-P;%
    XX0(1:n,i)=X0;%
    XX1(1:n,i)=X1;%
    XX2(1:n,i)=X2;%
    XX3(1:n,i)=X2+ug(i);%
    i=i+1;%
end