线弹性时程分析法子程序.docVIP

线弹性时程分析法（wilson-θ）子程序 ********************************************************** * * * ====== STEP-BY-STEP INTEGRAL PROGRAM ====== * * * ********************************************************** SUBROUTINE ESSIP(NI,NZ,ID,XI1,XI2,W1,W2,IC) W1，W2是第一、二阶频率，XI1,XI2是第一、二阻尼比，ID是地震烈度，IC是阻尼算法 REAL WK(30),TT(30),CM(30),PM(30,2),S(30) REAL*8 A(330,38),W(330),EK(30,30),FK(30,30) INTEGER NX(30) CHARACTER AA*20 COMMON /C5/A/C6/P(330)/C7/W/C17/JDW(150,3)/N0/FK COMMON /N1/CD(30,30)/N2/DX(30),VX(30),AX(30) COMMON /N3/YD(30),YV(30),YA(30)/N4/GG(1000) READ(12,*)IX,(CM(I),I=1,IX),(NX(I),I=1,IX)‘输入：振动位移个数IX，质量矩阵（向量）CM，振动位移对应的节点编号NX CLOSE (12) TSTEP=.02 ‘时间步长 THETA=1.4 IF(W2.EQ.0) IC=1‘如果W2＝0,则自由一个频率，取第一种阻尼算法 WRITE(*,*)Name of earthquake wave[EL-CENTRO]? READ(*,*) read(*,(A\))AA IF(AA.EQ. ) AA=EL-CENTRO OPEN(5,FILE=AA) READ (5,*)NMAX,(GG(I),I=1,NMAX) ‘输入地面运动加速度点数，加速度值 CLOSE(5) DO 4 I=1,IX S(I)=CM(I) “[M]”质量向量幅值给S（） DO 4 J=1,IX 4 CD(I,J)=0 “[C]”阻尼矩阵清零 OPEN (11,FILE=DSP.MID,STATUS=NEW)打开文件保存位移 OPEN (12,FILE=VEL.MID,STATUS=NEW) 打开文件保存速度 7 OPEN (13,FILE=ACC.MID,STATUS=NEW) 打开文件保存加速度 DO 10 I=1,IX 10 NX(I)=JDW(NX(I),1)利用节点定位向量，把振动位移自由度编号取出来 IF (ID.EQ.7) C=.22*.0981根据烈度确定加速度最大值 IF (ID.EQ.8) C=.4*.0981 IF (ID.EQ.9) C=.62*.0981 DO 15 I=1,NMAX 15 GG(I)=GG(I)*C（调整规格化地震波，从数据文件输入的地震加速度最大峰值都是100gal） TAU=THETA*TSTEP (τ=θΔt) READ(9)((A(I,J),J=1,NZ),I=1,NI) （总刚，NZ半带宽，NI总自由度数目） DO 60 I=1,IX DO 61 J=1,NI 61 P(J)=0 P(NX(I))=1.0对振动自由度上施加单位荷载 CALL JFC(NI,NZ,I) DO 65 J=1,IX 65 FK(J,I)=W(NX(J))保存柔度系数的一列,逐步形成柔度矩阵 60 CONTINUE CALL CVK(IX,FK) （求抗侧力总刚[K]E=[F]-1） IF(IC.EQ.1) THEN计算阻尼矩阵([C]=α[M]) T1=2*XI1*W1