線彈性時(shí)程分析法子程序

1、 線彈性時(shí)程分析法（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)輸入：振動(dòng)位移個(gè)數(shù)IX，質(zhì)量矩陣（向量）CM，振動(dòng)位移對(duì)應(yīng)的節(jié)點(diǎn)編號(hào)NX CLOSE (12) TSTEP=.02 時(shí)間步長(zhǎng) THETA=1.4 如果W20,則自由一個(gè)頻率，取第一種阻尼算法 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) 輸入地面運(yùn)動(dòng)加速度點(diǎn)數(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')打開文件保存位移 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)利用節(jié)點(diǎn)定位向量，把振動(dòng)位移自由度編號(hào)取出來根據(jù)烈度確定加速度最大值 DO 15 I=1,NMAX15 GG(I)=GG(I)*C（調(diào)整規(guī)格化地震波，從數(shù)據(jù)文件輸入的地震加速度最大峰值都是100gal） TAU=THETA*TSTEP (=t) READ(9)(A(I,J),J=1,NZ),I=1,NI) （總

5、剛，NZ半帶寬，NI總自由度數(shù)目） DO 60 I=1,IX DO 61 J=1,NI61 P(J)=0 P(NX(I)=1.0對(duì)振動(dòng)自由度上施加單位荷載 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）計(jì)算阻尼矩陣(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)時(shí)間步長(zhǎng)初始的慣性力，即荷載 DO 70 I=1,IX YD(I)=W(NX(I)初始位移（如豎向荷載作用下的位移） YV(I)=070 YA(I)=-GG(1)取初始加速度為地面運(yùn)動(dòng)加速度 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)計(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)時(shí)間步長(zhǎng)末的慣性力，荷載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正常步長(zhǎng)t時(shí)加速度增量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)不采用近似公式計(jì)算t(j+1)時(shí)刻的加速度160 YA(I)=(PM(I,2)-TT(I)/CM(I)代入運(yùn)動(dòng)方程求解t(j+1)時(shí)刻的加速度 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,*)&