CC--PROGRAM TO INTEGRATE EQUATION OF MOTION SUBJECTED TO
CC ACCELERATION INPUT. EQUATION IS INTEGRATED FOR A GIVEN
CC RANGE OF FREQUENCIES AND DAMPING PARAMETER.
CC MAXIMUM RESPONSE WILL DEFINE RESPONSE SPECTRUM WHICH
CC CAN BE USED IN MODAL DYNAMICS OPTION TO ESTIMATE MAXIMUM
CC RESPONSE OF THE CONSIDERED STRUCTURE
CC PROGRAM RESPON
SUBROUTINE ABQMAIN
C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(2,2),B(2,2),Q(2500),QP(2500),QPP(2500),
1ACC(2500),TAB(8),FR(1000),QM(1000),QPM(1000),
1QPPM(1000)
parameter(one=1.d0,two=2.d0,zero=0.d0,dt=1.d-2,fact=-386.09d0)
DATA NACC/1000/
C DATA NACC/2500/
C DATA FREQMIN/0.1d0/
C DATA FREQMAX/30.d0/
DATA INT/501/
C STORAGE ALLOCATED FOR A MAXIMUM OF 2500 TIMEPOINTS
C IN THE EARTHQUAKE HISTORY
CC READ ACCELERATION HISTORY FROM FILE QUAKE.AMP
CC (COPIED FROM CANTILEVER_QUAKEDATA.INP).
CC DISPLACEMENT SPECTRUM (FOR ABAQUS OPTION \*RESPONSE SPECTRUM)
CC WILL BE WRITTEN TO FILE QUAKE.DIS, VELOCITY SPECTRUM TO FILE QUAKE.VEL.
C
C**** TIME INTEGRATION FOR LINEAR ACCELERATION (EXACT SOLUTION)
C DATA DAMP DENOTES DAMPING AS PERCENTAGE OF CRITICAL DAMPING
C DATA FREQMIN AND FREQMAX DEFINE FREQUENCY RANGE
C DATA INT DEFINES NUMBER OF POINTS IN FREQUENCY RANGE
C** THIS INPUT ASSUMES THAT CANTILEVER_QUAKEDATA.INP HAS BEEN COPIED
C** TO QUAKE.AMP
OPEN(UNIT=100,STATUS='OLD',FILE='QUAKE.AMP')
C** OPEN the output files required, to generate figures 10,11 and 12
OPEN(UNIT=116,STATUS='NEW',FILE='QUAKE0.DIS')
OPEN(UNIT=117,STATUS='NEW',FILE='QUAKE2.DIS')
OPEN(UNIT=118,STATUS='NEW',FILE='QUAKE4.DIS')
OPEN(UNIT=119,STATUS='NEW',FILE='QUAKER0.DIS')
OPEN(UNIT=120,STATUS='NEW',FILE='QUAKER2.DIS')
OPEN(UNIT=121,STATUS='NEW',FILE='QUAKER4.DIS')
OPEN(UNIT=122,STATUS='NEW',FILE='QUAKE0.VEL')
OPEN(UNIT=123,STATUS='NEW',FILE='QUAKE2.VEL')
OPEN(UNIT=124,STATUS='NEW',FILE='QUAKE4.VEL')
OPEN(UNIT=125,STATUS='NEW',FILE='QUAKER0.VEL')
OPEN(UNIT=126,STATUS='NEW',FILE='QUAKER2.VEL')
OPEN(UNIT=127,STATUS='NEW',FILE='QUAKER4.VEL')
FRAC=ONE/DBLE(INT-1)
C*** INITIATE AMAT,BMAT BEFORE TIME INTEGRATION
DO 10 I=1,2
DO 10 J=1,2
A(I,J)=ZERO
B(I,J)=ZERO
10 CONTINUE
C** READ AMPLITUDE DATA AND STORE ON ACC(2500)
DO 30 I=1,625
READ(100,*) (TAB(LL),LL=1,8)
DO 40 K=1,4
ACC(K+4*(I-1))=TAB(2*K)*FACT
40 CONTINUE
30 CONTINUE
C** DO 990 is loop to generate data for frequency ranges of 0.1-30Hz
C** and 0.01-5Hz, for data shown in figures 10,11 and 12
DO 990 IFRQ=1,2
IF (IFRQ.EQ.1) THEN
FREQMIN=0.1
FREQMAX=30
IA=14
IB=20
ELSE IF (IFRQ.EQ.2) THEN
FREQMIN=0.01
FREQMAX=5
IA=16
IB=22
END IF
DAMP=ZERO
DFREQ=FREQMAX-FREQMIN
C** CHOOSE DAMPING.
C** DAMPING MUST BE LESS THAN CRITICAL (BETWEEN 0.0 AND 1.0).
DO 300 IKSI=1,3
IF(DAMP.GT.ONE) WRITE(6,11)
11 FORMAT(/,3X,50HTHIS PROGRAM IS WRITTEN FOR UNDERDAMPED CASES ONLY)
C** CHOOSE FREQUENCY FROM THE RANGE (FREQMIN,FREQMAX)
DO 200 IFREQ=1,INT
FREQIN=FREQMIN+FRAC*DFREQ*DBLE(IFREQ-1)
C** FREQMIN MUST BE GREATER THAN ZERO
PI=TWO*ASIN(ONE)
VXSI=DAMP
FREQ=FREQIN*TWO*PI
DEN=SQRT(ONE-VXSI*VXSI)
VRATIO=ONE/DEN
FREQE=DEN*FREQ
XSIWT=VXSI*FREQ*DT
ETAU=EXP(-XSIWT)
SINWT=SIN(FREQE*DT)
COSWT=COS(FREQE*DT)
A(1,1)=ETAU*(VXSI*VRATIO*SINWT+COSWT)
A(1,2)=ETAU*SINWT/FREQE
A(2,1)=-ETAU*FREQ*VRATIO*SINWT
A(2,2)=ETAU*(COSWT-VXSI*VRATIO*SINWT)
C
TXSI=(TWO*VXSI*VXSI-ONE)/FREQ/FREQ/FREQE/DT
XSIF=VXSI/FREQ/FREQE
FREQI=ONE/FREQ/FREQ
C
B(1,1)=ETAU*(-(XSIF+TXSI)*SINWT-
1 (FREQI+TWO*VXSI*FREQI/FREQ/DT)*COSWT)+
2 TWO*VXSI*FREQI/FREQ/DT
B(1,2)=ETAU*(TXSI*SINWT+TWO*VXSI*FREQI/FREQ/DT*COSWT)+
1 FREQI-TWO*VXSI*FREQI/FREQ/DT
B(2,1)=ETAU*(-(FREQE*COSWT-VXSI*FREQ*SINWT)*(TXSI+XSIF)+
1 (FREQE*SINWT+VXSI*FREQ*COSWT)*(FREQI+TWO*VXSI*FREQI/
1 FREQ/DT))-FREQI/DT
B(2,2)=ETAU*((FREQE*COSWT-VXSI*FREQ*SINWT)*TXSI-
1 (FREQE*SINWT+VXSI*FREQ*COSWT)*TWO*VXSI*FREQI/
2 FREQ/DT)+FREQI/DT
DO 100 IT=1,NACC
IF(IT.EQ.1)THEN
C** INITIAL CONDITIONS
T=0.d0
Q(1)=0.d0
QP(1)=0.d0
QPP(1)=0.d0
ELSE
T=T+DT
Q(IT)=A(1,1)*Q(IT-1)+A(1,2)*QP(IT-1)+B(1,1)*ACC(IT-1)+B(1,2)
1 *ACC(IT)
QP(IT)=A(2,1)*Q(IT-1)+A(2,2)*QP(IT-1)+B(2,1)*ACC(IT-1)+
1 B(2,2)*ACC(IT)
C absolute acceleration record
QPP(IT)=-QP(IT)*TWO*VXSI*FREQ-Q(IT)*FREQ*FREQ
ENDIF
100 CONTINUE
QMAX=0.d0
QPMAX=0.d0
QPPMAX=0.d0
DO 110 II=1,NACC
QABS=ABS(Q(II))
QPABS=ABS(QP(II))
QPPABS=ABS(QPP(II))
QMAX=MAX(QMAX,QABS)
QPMAX=MAX(QPMAX,QPABS)
QPPMAX=MAX(QPPMAX,QPPABS)
110 CONTINUE
QM(IFREQ)=QMAX
QPM(IFREQ)=QPMAX
QPPM(IFREQ)=QPPMAX
FR(IFREQ)=FREQIN
200 CONTINUE
DO 210 LI=1,INT
WRITE(100+IFRQ+IKSI+IA,211) QM(LI),FR(LI),DAMP
WRITE(100+IFRQ+IKSI+IB,211) QPM(LI),FR(LI),DAMP
ID=1
210 CONTINUE
211 FORMAT(1X,E12.5,1H,,E12.5,1H,,E12.5)
DAMP=DAMP+0.02d0
300 CONTINUE
990 CONTINUE
C
CLOSE(100)
CLOSE(116)
CLOSE(117)
CLOSE(118)
CLOSE(119)
CLOSE(120)
CLOSE(121)
CLOSE(122)
CLOSE(123)
CLOSE(124)
CLOSE(125)
CLOSE(126)
CLOSE(127)
C
STOP
END