線彈性時程分析法子程序

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)輸入：振動位移個數(shù)IX，質(zhì)量矩陣（向量）CM，振動位移對應(yīng)的節(jié)點編號NX CLOSE (12) TSTEP=.02 時間步長 THETA=1.4 IF(W2.EQ.0) IC=1如果W20,則自由一個頻率，取第一種阻尼算法 WRITE(*,*)'Name of earthquake waveEL

3、-CENTRO?' READ(*,*) read(*,'(A)')AA IF(AA.EQ.' ') AA='EL-CENTRO' OPEN(5,FILE=AA) READ (5,*)NMAX,(GG(I),I=1,NMAX) 輸入地面運動加速度點數(shù)，加速度值 CLOSE(5) DO 4 I=1,IX S(I)=CM(I) “M”質(zhì)量向量幅值給S（） DO 4 J=1,IX4 CD(I,J)=0 “C”阻尼矩陣清零 OPEN (11,FILE='DSP.MID',STATUS='NEW')打開文件保存位移 O

4、PEN (12,FILE='VEL.MID',STATUS='NEW') 打開文件保存速度7 OPEN (13,FILE='ACC.MID',STATUS='NEW') 打開文件保存加速度 DO 10 I=1,IX10 NX(I)=JDW(NX(I),1)利用節(jié)點定位向量，把振動位移自由度編號取出來 IF (ID.EQ.7) C=.22*.0981根據(jù)烈度確定加速度最大值 IF (ID.EQ.8) C=.4*.0981 IF (ID.EQ.9) C=.62*.0981 DO 15 I=1,NMAX15 GG(I)=GG(I)*C（

5、調(diào)整規(guī)格化地震波，從數(shù)據(jù)文件輸入的地震加速度最大峰值都是100gal） TAU=THETA*TSTEP (=t) READ(9)(A(I,J),J=1,NZ),I=1,NI) （總剛，NZ半帶寬，NI總自由度數(shù)目） 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)保存柔度系數(shù)的一列,逐步形成柔度矩陣60 CONTINUE CALL CVK(IX,FK) （求抗側(cè)力總剛KE=F-1） IF(IC.EQ.1) THEN計算阻尼矩陣(

6、C=M) T1=2*XI1*W1 DO 5 I=1,IX 5 CD(I,I)=T1*CM(I) ENDIF IF(IC.EQ.2) THEN(C=M+K) T1=2.*W1*W2*(Xi1*W2-Xi2*W1)/(W2*W2-W1*W1) T2=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)

7、時間步長初始的慣性力，即荷載 DO 70 I=1,IX YD(I)=W(NX(I)初始位移（如豎向荷載作用下的位移） YV(I)=070 YA(I)=-GG(1)取初始加速度為地面運動加速度 WRITE(11,*) 'Step 1' 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 WRI

8、TE(*,*)' 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,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

9、 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(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 V

10、ector 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 Determine 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)=

11、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(