SUBROUTINE ABQMAIN C C PROGRAM LDYN C C----PROGRAM TO GENERATE "EXACT" SOLUTION C C READ AMPLITUDE DATA ON UNIT 1 C C FOR DIRECT INTEGRATION ANALYSIS: C WRITE TO UNIT 9 C REQUIRES UNIT 8 C INCLUDE 'aba_param.inc' C CHARACTER*80 FNAME DIMENSION OMEGA(20),BL(20),AMASS(20),GAMMA(20),DAMPR(20),D(20), 1 V(20),A(20),ACC(2500),F(8),FREQD(20),S1(20),S2(20),S3(20),S4(20) 1 ,S5(20),S6(20) DIMENSION ARRAY(513),JRRAY(NPRECD,513),ABA(2,3) DIMENSION LRUNIT(2,1) EQUIVALENCE(ARRAY(1),JRRAY(1,1)) DATA BE/30000000./, BI/.6666667/, BM/0.001456/, BLL/300./ DATA NMODES/10/ DATA DAMP/0.0/ DATA NACC/1001/ C FOR MODAL DYNAMIC C DATA NACC/2500/ C DT=0.01 C STORAGE ALLOCATED FOR A MAXIMUM OF 20 MODES C AND 2500 TIMEPOINTS IN THE EARTHQUAKE HISTORY C FNAME='3cantilever_restart' C-----OPEN FILE TO READ EARTHQUAKE RECORD C-----THIS INPUT ASSUMES THAT CANTILEVER_QUAKEDATA.INP C-----HAS BEEN COPIED TO QUAKE.AMP OPEN(UNIT=1,STATUS='OLD',FILE='QUAKE.AMP') C C-----ASSIGN UNITS AND OPEN FILES TO WRITE ASCII DATA OF EXACT RESULTS OPEN(UNIT=50,FILE='exactdisp') OPEN(UNIT=51,FILE='exactvelo') OPEN(UNIT=52,FILE='exactaccl') C-----INITIALIZE RESULTS FILE (.fil) FOR VALUES CALCULATED IN THIS CODE LRUNIT(1,1)=8 LRUNIT(2,1)=2 LOUTF=2 CALL INITPF(FNAME,1,LRUNIT,LOUTF) JUNIT=8 CALL DBRNU(JUNIT) 10 FORMAT(8F10.0) C-----WRITE MODEL DATA DEFINED AS DATA TO SCREEN WRITE(6,20) BE,BI,BM,BLL 20 FORMAT(////,5X,'YOUNG''S MODULUS =',5X,1PG12.4, 1 /,5X,'MOMENT OF INERTIA =',3X,1PG12.4, 1 /,5X,'MASS PER UNIT LENGTH =',1PG12.4, 3 /,5X,'LENGTH OF CANTILEVER =',1PG12.4) WRITE(6,40) NMODES 40 FORMAT(/,5X,'NUMBER OF MODES',7X,I5) WRITE(6,42) 42 FORMAT(/) C-----START CALCULATION C-----FIRST CALCULATE VARIABLES REQUIRED FOR EVALUATION OF EXACT SOLUTION ONE=1. TWO=2. PI=TWO*ASIN(ONE) A1=SQRT((BE*BI)/(BM*BLL**4)) DO 60 KM=1,NMODES BL(KM)=(KM-0.5)*PI DO 45 KITER=1,100 X=BL(KM) IF(X.GT.75.) GO TO 48 DX=(ONE+COSH(X)*COS(X))/(COSH(X)*SIN(X)-SINH(X)*COS(X)) BL(KM)=BL(KM)+DX IF(ABS(DX).LT.1.E-10) GO TO 48 45 CONTINUE 48 CONTINUE OMEGA(KM)=BL(KM)**2*A1 A22=OMEGA(KM)/(TWO*PI) N1=201 AM=0. DX=BLL/FLOAT(N1-1) X=0. AREA1=DX*0.66666666667 AG=0. DO 55 K1=1,N1 AREA=AREA1 IF(K1.EQ.1.OR.K1.EQ.N1) AREA=0.5*AREA IF(((K1/2)*2).EQ.K1) AREA=TWO*AREA BX= BL(KM)*X/BLL DZ=2.*(SIN(BL(KM))*COSH(BL(KM))-COS(BL(KM))*SINH(BL(KM))) A2=((SINH(BL(KM))+SIN(BL(KM)))*(COSH(BX)-COS(BX))- 1 (COSH(BL(KM))+COS(BL(KM)))*(SINH(BX)-SIN(BX)))/DZ AM=AM+A2*A2*AREA AG=AG+A2*AREA X=X+DX 55 CONTINUE AMASS(KM)=AM*BM GAMMA(KM)=AG*BM C-----WRITE VARIOUS PRELIMINARY CALCULATION RESULTS TO THE SCREEN WRITE(6,50) KM,A22,AMASS(KM),GAMMA(KM) 50 FORMAT(5X,'MODE',I10,2X,'FREQUENCY',1PG12.5,2X, 1 'GENERALIZED MASS',1PG12.5,2X,'PARTICIPATION FACTOR',1PG12.5) 60 CONTINUE WRITE(6,65) DAMP 65 FORMAT(/,5X,'DAMPING RATIO =',1PG12.4) DO 70 KM=1,NMODES DAMPR(KM)=DAMP*OMEGA(KM) FREQD(KM)=SQRT(OMEGA(KM)**2-DAMPR(KM)**2) 70 CONTINUE WRITE(6,75) NACC 75 FORMAT(/,5X,'NUMBER OF ACCELERATION POINTS =',I10) FACT=-386.09 I1=1 DO 80 K1=1,99999 READ(1,10) (F(I2),I2=1,8) ACC(I1)=F(2)*FACT ACC(I1+1)=F(4)*FACT ACC(I1+2)=F(6)*FACT ACC(I1+3)=F(8)*FACT I1=I1+4 IF(I1.GT.NACC) GO TO 82 80 CONTINUE 82 CONTINUE DO 90 KM=1,NMODES D(KM)=0. V(KM)=0. A(KM)=0. 90 CONTINUE Y1=ACC(1) IACC=2 NACC1=NACC-1 T=0. DO 100 KM=1,NMODES BT=DAMPR(KM)*DT OM2=OMEGA(KM)**2 EBT=EXP(-BT) S1(KM)=EBT*COS(FREQD(KM)*DT) S2(KM)=EBT*SIN(FREQD(KM)*DT)/FREQD(KM) S3(KM)=GAMMA(KM)/(AMASS(KM)*OMEGA(KM)**2) S4(KM)=TWO*DAMPR(KM)/(OMEGA(KM)**2) S5(KM)=DAMPR(KM)*GAMMA(KM)/(AMASS(KM)*OMEGA(KM)**2) S6(KM)=DAMPR(KM)**2-FREQD(KM)**2 100 CONTINUE WRITE(6,110) 110 FORMAT(//,8X,'TIME',15X,'MODAL ANALYSIS',22X,'ABAQUS *DYNAMIC', 1 /,22X,'D',11X,'V',11X,'A',11X,'D',11X,'V',11X,'A',5X, 2 'STRAIN ENERGY',2X,'KIN. ENERGY',2X,'TOTAL ENERGY',/) C----LOOP TO CALCULATE EXACT RESULTS AND EVALUATE RELATIVE RESULTS FROM ABAQUS ANALYSIS cprotect keyprv=0 cprotect DO 200 K1=1,NACC1 SLOPE=(ACC(IACC)-Y1)/DT C-----INCREMENT TIME (T) AND RESET EXACT RESULTS VARIABLES T=T+DT TD=0. TV=0. TA=0. 111 CONTINUE C-----ACCESS ABAQUS RESULTS FILE (cantilever_restart.fil) DO 112 K3=1,99999 CALL DBFILE(0,JRRAY,JRCD) cprotect if(jrcd.ne.0 .and. keyprv.eq.2001) then c normal end of file reached go to 960 else if(jrcd.ne.0) then go to 900 end if keyprv=jrray(1,2) cprotect CALL DBFILW(1,JRRAY,JRCD) IF(JRRAY(1,2).EQ.2000) GO TO 113 112 CONTINUE GO TO 900 113 CONTINUE TABA=ARRAY(3) IF(TABA.LT.0.1*DT) GO TO 111 DO 120 K3=1,99999 CALL DBFILE(0,JRRAY,JRCD) cprotect if(jrcd.ne.0 .and. keyprv.eq.2001) then c normal end of file reached go to 960 else if(jrcd.ne.0) then go to 900 end if keyprv=jrray(1,2) cprotect C-----DETERMINE RECORD KEY OF RESULTS READ KEY=JRRAY(1,2) IF(KEY.EQ.2001) GO TO 122 IF(KEY.EQ.1999) GO TO 118 C-----IF RECORD READ IS DISPLACEMENT, ACCELERATION OR VELOCITY CONTINUE IF(KEY.GE.101.AND.KEY.LE.103) GO TO 115 CALL DBFILW(1,JRRAY,JRCD) GO TO 120 115 CONTINUE I1=0 C-----DETERMINE WHETHER RECORD READ IS FOR NODE 1 OR NODE 11 IN MODEL IF(JRRAY(1,3).EQ.11) I1=1 IF(JRRAY(1,3).EQ.1) I1=2 IF(I1.EQ.0) GO TO 120 I2=KEY-100 ABA(I1,I2)=ARRAY(4) GO TO 120 118 CONTINUE ESTR=ARRAY(5) EKIN=ARRAY(4) ETOT=EKIN+ESTR 120 CONTINUE 122 CONTINUE C-----CALCULATE RELATIVE DISPLACEMENT, VELOCITY AND ACCELERATION FOR ABAQUS RESULTS C-----RELATIVE DISPLACEMENT=(U1 AT NODE 1) - (U1 AT NODE 11) C-----SIMILAR FOR VELOCITY AND ACCELERATION ABA(2,1)=ABA(2,1)-ABA(1,1) ABA(2,2)=ABA(2,2)-ABA(1,2) ABA(2,3)=ABA(2,3)-ABA(1,3) DO 150 KM=1,NMODES D1=D(KM) V1=V(KM) A1=A(KM) OM2=OMEGA(KM)**2 G1=GAMMA(KM)*Y1/AMASS(KM) C-----CALCULATE EXACT SOLUTION FOR CURRENT MODE FOR DISPLACEMENT (D), VELOCITY (V) C-----AND ACCELERATION (A) D(KM)=S1(KM)*(D1+S3(KM)*(S4(KM)*SLOPE-Y1)) 1 +S2(KM)*(D1*DAMPR(KM)+V1+S3(KM)*(S4(KM)*DAMPR(KM)*SLOPE 2 -Y1*DAMPR(KM)-SLOPE)) 3 +S3(KM)*(Y1+SLOPE*DT-S4(KM)*SLOPE) V(KM)=S1(KM)*(V1-S3(KM)*SLOPE) 1 +S2(KM)*(-D1*OM2-DAMPR(KM)*V1-S5(KM)*SLOPE+G1)+S3(KM)*SLOPE A(KM)=S1(KM)*(-TWO*DAMPR(KM)*V1-OM2*D1+G1) 1 +S2(KM)* (S6(KM)*V1+GAMMA(KM)*SLOPE/AMASS(KM) 2 +DAMPR(KM)*(D1*OM2-G1)) C-----RUNNING SUM FOR CURRENT MODE TD=TD+D(KM) TV=TV+V(KM) TA=TA+A(KM) 150 CONTINUE C-----WRITE RESULTS FOR CURRENT MODE TO SCREEN WRITE(6,160) T,TD,TV,TA,ABA(2,1),ABA(2,2),ABA(2,3),ESTR,EKIN,ETOT 160 FORMAT(5X,1P10G12.4) C C-----WRITE ASCII DATA OF EXACT SOLUTION TO FILES WRITE(50,170) T,TD WRITE(51,170) T,TV WRITE(52,170) T,TA 170 FORMAT (G12.4,5X,G12.4) C-----SET RECORD KEYS AND WRITE RELATIVE DISPLACEMENT, VELOCITY AND ACCELERATION C-----OF ABAQUS RESULTS (cantilever_restart) TO CURRENT RESULTS FILE C-----ALSO WRITE EXACT SOLUTION (TD, TV, TA) TO CURRENT RESULTS FILE C-----DISPLACEMENT RESULTS, RECORD KEY 101 JRRAY(1,1)=5 JRRAY(1,2)=101 JRRAY(1,3)=1 ARRAY(4)=ABA(2,1) ARRAY(5)=TD CALL DBFILW(1,JRRAY,JRCD) C-----VELOCITY RESULTS, RECORD KEY 102 JRRAY(1,2)=102 ARRAY(4)=ABA(2,2) ARRAY(5)=TV CALL DBFILW(1,JRRAY,JRCD) C-----ACCELERATION RESULTS, RECORD KEY 103 JRRAY(1,2)=103 ARRAY(4)=ABA(2,3) ARRAY(5)=TA C-----WRITE END OF CURRENT INCREMENT RECORD CALL DBFILW(1,JRRAY,JRCD) JRRAY(1,1)=2 JRRAY(1,2)=2001 CALL DBFILW(10,JRRAY,JRCD) Y1= ACC(IACC) IACC=IACC+1 200 CONTINUE GO TO 990 900 CONTINUE WRITE(6,950) 950 FORMAT(' UNEXPECTED ERROR ') GO TO 990 960 CONTINUE WRITE(6,970) 970 FORMAT(' NORMAL END OF FILE REACHED ') 990 CONTINUE STOP END