线弹性时程分析法子程序

1、 线弹性时程分析法（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/

2、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 如果W20,则自由一个频率，取第一种阻尼算法 WRITE(*,*)'Name of earthquake waveEL-CENTRO?' RE

3、AD(*,*) 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,IX4 CD(I,J)=0 “C”阻尼矩阵清零 OPEN (11,FILE='DSP.MID',STATUS='NEW')打开文件保存位移 OPEN (12,FILE=

4、9;VEL.MID',STATUS='NEW') 打开文件保存速度7 OPEN (13,FILE='ACC.MID',STATUS='NEW') 打开文件保存加速度 DO 10 I=1,IX10 NX(I)=JDW(NX(I),1)利用节点定位向量，把振动位移自由度编号取出来根据烈度确定加速度最大值 DO 15 I=1,NMAX15 GG(I)=GG(I)*C（调整规格化地震波，从数据文件输入的地震加速度最大峰值都是100gal） TAU=THETA*TSTEP (=t) READ(9)(A(I,J),J=1,NZ),I=1,NI) （总

5、刚，NZ半带宽，NI总自由度数目） DO 60 I=1,IX DO 61 J=1,NI61 P(J)=0 P(NX(I)=1.0对振动自由度上施加单位荷载 CALL JFC(NI,NZ,I) DO 65 J=1,IX65 FK(J,I)=W(NX(J)保存柔度系数的一列,逐步形成柔度矩阵60 CONTINUE CALL CVK(IX,FK) （求抗侧力总刚KE=F-1）计算阻尼矩阵(C=M) T1=2*XI1*W1 DO 5 I=1,IX 5 CD(I,I)=T1*CM(I) ENDIF (C=M+K) T1=2.*W1*W2*(Xi1*W2-Xi2*W1)/(W2*W2-W1*W1) T2=

6、2*(Xi2*W2-Xi1*W1)/(W2*W2-W1*W1) DO 6 I=1,IX CD(I,I)=T1*CM(I) DO 6 J=1,IX 6 CD(I,J)=CD(I,J)+T2*FK(I,J) ENDIFC Introduce the Initial Conditions DO 35 I=1,IX35 PM(I,1)=-CM(I)*GG(1)时间步长初始的惯性力，即荷载 DO 70 I=1,IX YD(I)=W(NX(I)初始位移（如竖向荷载作用下的位移） YV(I)=070 YA(I)=-GG(1)取初始加速度为地面运动加速度 WRITE(11,*) 'Step 1'

7、; WRITE(11,'(6F12.4)') (YD(I),I=1,IX) WRITE(12,*) 'Step 1' WRITE(12,'(6F12.4)') (YV(I),I=1,IX)C - Start of Step-by-Step Integral - DO 200 JJ=1,NMAX-1 WRITE(*,*)' Step ',JJ DO 64 J=1,NI64 P(J)=0.C Form the Equivalent Stiffness Matrix EK DO 20 I=1,IX DO 21 J=1,IX21 EK(I

8、,J)=FK(I,J)+3./TAU*CD(I,J)计算等效刚度矩阵20 EK(I,I)=EK(I,I)+CM(I)*6./TAU/TAU DO 80 I=1,IX 80 PM(I,2)=-CM(I)*GG(JJ+1)时间步长末的惯性力，荷载C Form the Equivalent Loading Vector DO 90 I=1,IX “C” 90 WK(I)=3.*YV(I)+TAU/2.*YA(I) DO 92 I=1,IX TT(I)=0 DO 92 K=1,IX92 TT(I)=TT(I)+CD(I,K)*WK(K) DO 95 I=1,IX C=CM(I)*(6./TAU*YV(

9、I)+3.*YA(I) “M” 95 AX(I)=TT(I)+C+(PM(I,2)-PM(I,1)*THETA等效荷载()C Solve for the Incremental Displacement Vector CALL CCHOL(EK,AX,IX) “解方程”C Determine the Incremental Acceleration Vector DO 110 I=1,IX110 AX(I)=6.*AX(I)/TAU/TAU-6.*YV(I)/TAU-3.*YA(I) DO 120 I=1,IX120 AX(I)=AX(I)/THETA正常步长t时加速度增量C Determin

10、e the Displacement,Velocity Vectors at Time t(j+1) DO 130 I=1,IX DX(I)=YV(I)*TSTEP+YA(I)*TSTEP*2/2.+AX(I)*TSTEP*2/6. YD(I)=DX(I)+YD(I) VX(I)=YA(I)*TSTEP+AX(I)*TSTEP/2.130 YV(I)=YV(I)+VX(I)C Determine the Acceleration Vectors at Time t(j+1) DO 132 I=1,IX TT(I)=0 DO 132 K=1,IX132 TT(I)=TT(I)+FK(I,K)*YD(K)+CD(I,K)*YV(K) DO 160 I=1,IXC160 YA(I)=YA(I)+AX(I)不采用近似公式计算t(j+1)时刻的加速度160 YA(I)=(PM(I,2)-TT(I)/CM(I)代入运动方程求解t(j+1)时刻的加速度 WRITE(11,*)' Step ',JJ+1 WRITE(11,'(6F12.4)') (YD(I),I=1,IX) WRITE(12,*)' Step ',JJ+1 WRITE(12,'(6F12.4)') (YV(I),I=1,IX) WRITE(13,*)&