PROGRAM DYNA1 C *** TIME INTEGRATION IMPLICIT ALGORITHM C * PLANE STRAIN (CONSTANT STRAIN 3 AND 4 NODED ELEMENTS) C * ORDER OF STRAINS & STRESSES: X, Y, XY, Z C * SOIL DOES NOT BEAR TENSILE STRESS C MINIMUM PRINCIPAL STRESS & SIGUMA-Z C * STRESS-STRAIN RELATIONSHIP C ACCORDING TO D. R. J. OWEN & E. HINTON C * NON-ASSOCIATED FLOW RULE C * IN CALCULATION, TENSILE STRESS IS POSITIVE, C BUT IN OUTPUT, COMPRESSIVE STRESS IS POSITIVE. C * GROUND-WATER LEVEL (STATIC PORE-WATER PRESSURE) C CAN BE CONSIDERED. C * WHEN USING INTERFACE-ELEMENTS IN ACTUAL SIZE MODEL, C DIVERGENCE MAY OCCUR DUE TO VIBRATION BY GRAVITY FORCE, C BECAUSE THE ELEMENTS BEAR NO TENSILE STRESS. C * IN ACTUAL SIZE MODEL, CONSIDER GRAVITY FORCE. C TIME LENGTH OF A FEW SECONDS IS REQUIRED FOR STABILIZING C THE EFEECT OF GRAVITY FORCE. C * YIELD FUNCTION FOR INTERFACE ELEMENTS IS ADDED C IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA11/IREQD(10),IREQS(10),NOUTD(6550) COMMON/LA12/ACCEH(6550),ACCEV(6550) COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA18/ACCEI(2220),ACCEL(2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) COMMON/LA20/IYIEL(2150),IUNLO(2150) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) C OPEN(5,FILE='DADYNA-1',STATUS='UNKNOWN') OPEN(6,FILE='PRDYNA-1',STATUS='UNKNOWN') OPEN(8,FILE='PRDISP',STATUS='UNKNOWN') OPEN(9,FILE='PRVELO',STATUS='UNKNOWN') OPEN(10,FILE='PRACCE',STATUS='UNKNOWN') OPEN(12,FILE='M-AVS2.INP',STATUS='UNKNOWN') OPEN(13,FILE='M-AVS3.INP',STATUS='UNKNOWN') C CALL INPUTD CALL LINKIN CALL INTIME CALL PREVOS CALL LOADPL CALL LUMASS CALL GSTIFF C DO 10 ISTEP=1,NSTEP DO 20 IITER=1,MITER CALL IMPEXP(IITER,ISTEP) CALL RESEPL(IITER,ISTEP) CALL ITRATE(IITER) IF(NCHEK.EQ.1) GO TO 15 20 CONTINUE 15 CALL OUTDYN(IITER,ISTEP) IOUT=NOUTD(ISTEP) IF(IOUT.LE.0) GO TO 10 CALL MAVS2 CALL MAVS3 10 CONTINUE C CLOSE(5) CLOSE(6) CLOSE(8) CLOSE(9) CLOSE(10) CLOSE(12) CLOSE(13) STOP END C * * * * * * * * * * * * SUBROUTINE INPUTD C *** INPUT ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) DIMENSION TITLE(20),TY(2150),ZZ(1998) DIMENSION IDA(10),IDAQ(10),NET(10),IWA(2150) C C NN: MAXIMUM NUMBER OF DISPLACEMENT VARIABLES NN=2220 DO 2 I=1,10 2 IDAQ(I)=0 DO 4 I=1,6 4 NET(I)=0 C TITLE(I): PROJECT NAME READ(5,100) (TITLE(I),I=1,15) WRITE(6,100) (TITLE(I),I=1,15) C NPOIN: NUMBER OF NODES, NELEM:NUMBER OF ELEMENTS C NMATS: NUMBER OF MATERIALS C NIST : NUMBER OF INITIAL STRESS GROUPS C NPREV=1: CONSIDER INITIAL STRESSES IN EACH ELEMENT C 0: NOT CONSIDER C IAS=0: ASSOCIATED FLOW RULE, 1:NON-ASSOCIATED C IPR=0: NOT PRINT NODAL & ELEMENT DATA, 1:PRINT C IPW=0: NOT CONSIDER WATER PRESSURE C 1: CALCULATE PRE-EFFECTIVE STRESSES C PERFORM FE-ANALYSIS USING UNDER-WATER UNIT-WEIGHT C 2: CALCULATE & CONSIDER WATER PRESSURE C READ(5,101) NPOIN,NELEM,NMATS,NIST,NPREV,IAS,IPR,IPW WRITE(6,200) NPOIN,NELEM,NMATS,NIST,NPREV,IAS,IPR,IPW C *** NODAL DATA C K:NODE NO., XX(I):X-COORDINATE, YY(I):Y-COORDINATE C IX(I)=0: X-DIRECTION FREE, 1=FIXED C IY(I)=0: Y-DIRECTION FREE, 1=FIXED C IQ(I)=0: ROTATION FREE, 1=FIXED C WRITE(6,201) DO 10 IP=1,NPOIN READ(5,102) K,XX(IP),YY(IP),IX(IP),IY(IP),IQ(IP) IF(IPR.EQ.0) GO TO 10 WRITE(6,102) K,XX(IP),YY(IP),IX(IP),IY(IP),IQ(IP) 10 CONTINUE C *** ELEMENT DATA C K: ELEMENT NO., IJK(L,I): NODE NO. C K2(IE)=ELEMENT TYPE: TRUSS=1, BEAM=2, TEXTILE=3, C PLANE-STRAIN (ELASTIC=4, ELASTIC-PLASTIC=5), C INTERFACE=6 (I-J:SHEAR DIRECTION, ANTI-CLOCKWISE) C MAT(IE): NO. OF MATERIAL C IST(IE): NO. OF SET OF ACTUAL INITIAL STRESSES C IWA(IE): 0=UPPER WATER-TABLE, 1=UNDER WATER-TABLE C *** IF DATA(I)=0, PRECEDING VALUE IS EMPLOYED C DO 12 IE=1,NELEM READ(5,103) K,(IJK(IE,I),I=1,4),(IDA(I),I=1,4) DO 14 I=1,4 IF(IDA(I).EQ.0) IDA(I)=IDAQ(I) IF(IDA(I).LT.0) IDA(I)=0 IDAQ(I)=IDA(I) 14 CONTINUE K2(IE)= IDA(1) MAT(IE)=IDA(2) IST(IE)=IDA(3) IWA(IE)=IDA(4) N=K2(IE) NET(N)=NET(N)+1 12 CONTINUE WRITE(6,202) (NET(I),I=1,6) WRITE(6,203) DO 20 IE=1,NELEM IF(IPR.EQ.0) GO TO 20 WRITE(6,204) IE,(IJK(IE,I),I=1,4),K2(IE),MAT(IE),IST(IE),IWA(IE) 20 CONTINUE C *** MATERIAL PROPERTIES C TRUSS 1:YOUNG'S MODULUS, 2:CROSS AREA, 3:0, 4:RHO, 5:0, 6:0, 7:0 C 8:ALFA, 9:BETA C BEAM 1;YOUNG'S MODULUS, 2:CROSS AREA, 3:MOMENT OF INERTIA, 4:RHO C 5:0, 6:0, 7:0, 8:ALFA, 9:BETA, 10:0 C PLANE-STRAIN 1:YOUNG'S MODULUS, 2:POISSON'S RATIO, 3:1.0, 4:DENSITY C 5:COHESION C, 6:ANGLE OF SHEAR RESISTANCE, 7:ANGLE OF DILATANCY C 8:ALFA, 9:BETA (DAMPING PARAMETERS), 10:0 C 10:RATIO OF REDUCING COMPRESSIVE STRENGTH TO TENSILE STRENGTH C INTERFACE 1:SHEAR RIGIDITY, 2:YOUNG'S MODULUS, 3:POISSON'S RATIO, C 4:DENSITY, 5:COHESION C, 6:ANGLE OF SHEAR RESISTANCE C 7:ANGLE OF DILATANCY,8:ALFA, 9:BETA, C CALCULATE DAMPING PARAMETERS, ACCORDING TO SPECIFIC PERIOD. C WRITE(6,205) DO 30 MT=1,NMATS READ(5,104) MT1,(PRP(MT,I),I=1,10) IF(IAS.LE.0) PRP(MT,7)=PRP(MT,6) WRITE(6,206) MT1,(PRP(MT,I),I=1,10) COHES=PRP(MT,5) COMST=2.D0*COHES TENST=COMST*PRP(MT,10) C HEREAFTER PRP(MT,10) MEANS TENSILE STRENGTH PRP(MT,10)=TENST 30 CONTINUE C C *** CALCULATE PORE-WATER PRESSURE WRITE(6,207) DO 40 IE=1,NELEM DO 40 I=1,4 40 PW(IE,I)=0.D0 IF(IPW.LT.2) GO TO 52 C USE STRIN(IE,I) & STRAG(IE,I) TEMPORALILY DO 42 I=1,2 42 READ(5,105) (STRIN(IE,I),IE=1,NELEM) DO 44 I=1,2 44 READ(5,105) (STRAG(IE,I),IE=1,NELEM) DO 46 IE=1,NELEM IF(IWA(IE).LE.0) GO TO 46 DO 48 I=1,2 48 PW(IE,I)=-(STRIN(IE,I)-STRAG(IE,I)) PW(IE,4)=PW(IE,1) 46 CONTINUE DO 50 I=1,4 50 WRITE(6,105) (PW(IE,I),IE=1,NELEM) 52 CONTINUE C 100 FORMAT(15A4) 101 FORMAT(15I5) 102 FORMAT(I5,2F10.3,3I5) 103 FORMAT(10I5) 104 FORMAT(I2,10E10.3) 105 FORMAT(10F8.3) 200 FORMAT(/'NODES=',I5,' ELEMENTS=',I5,' MATERIALS=',I3 */' INITIAL STRESS GROUPS=',I3 */' GIVE INITIAL STRESSES IN EACH ELEMENT (YES=1)=',I2 */' FLOW-RULE=',I2,' PRINT=',I2 */' PORE-WATER PRESSURE=',I2) 201 FORMAT(/'NODAL DATA'/' NO.',9X,'X',9X,'Y',' X-F Y-F R-F') 202 FORMAT(/'TRUSS=',I3,' BEAM=',I3,' TEXTILE=',I3,' P-STRAIN=',I3 */' P-STRAIN(NL)=',I3,' INTERFACE=',I3) 203 FORMAT(/'ELEMENT DATA'/' NO. NODES TYPE MAT IST WATER') 204 FORMAT(10I5) 205 FORMAT(/'MATERIAL PARAMETERS' * /'NO. E NYU THICK RHO C FAI DELTA ALFA BETA TENSION'/) 206 FORMAT(I2,10E10.3) 207 FORMAT(/'PORE-WATER PRESSURE') RETURN END C * * * * * * * * * * SUBROUTINE LINKIN C *** LINKS WITH PROFILE SOLVER IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) DIMENSION LX(6000),NDF(1998) C DO 10 IP=1,NPOIN 10 NDF(IP)=2 DO 12 IE=1,NELEM IF(K2(IE).NE.2) GO TO 12 I=IJK(IE,1) J=IJK(IE,2) NDF(I)=3 NDF(J)=3 12 CONTINUE DO 20 IP=1,NPOIN LY(3*IP-2)=IX(IP) LY(3*IP-1)=IY(IP) IF(NDF(IP).EQ.2) IQ(IP)=1 20 LY(3*IP)=IQ(IP) NX=3*NPOIN LOC=0 DO 22 L=1,NX IF(LY(L).EQ.0) GO TO 24 LY(L)=NN GO TO 22 24 LOC=LOC+1 LX(LOC)=L LY(L)=LOC 22 CONTINUE NSIZE=LOC WRITE(6,200) NSIZE 200 FORMAT(/'NUMBER OF VARIABLES',I5) IF(NSIZE.GE.NN) STOP RETURN END C * * * * * * * * * * * * SUBROUTINE INTIME C *** INITIAL VALUES AND TIME INTEGRATION DATA IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA11/IREQD(10),IREQS(10),NOUTD(6550) COMMON/LA12/ACCEH(6550),ACCEV(6550) DIMENSION LOUTD(100) C C *** READ TIME STEPPING AND SELECTIVE OUTPUT PARAMETERS C NSTEP: NUMBER OF TOTAL TIME-STEPS C NOUTP: NUMBER OF SELECTED TIME-STEPS FOR OUTPUT C NREQD: NUMBER OF SELECTED NODES FOR OUTPUT C NREQS: NUMBER OF SELECTED ELEMENTS FOR OUTPUT C NACCE: NUMBER OF STEPS FOR ACCERALATION INPUT C IFUNC=0: INPUT ACCERALATION TIME RECORD C =1: HEAVISIDE FUNCTION C =2: FUNCTION=AZERO+BZERO*SIN(ARGUM) C ARGUM=OMEGA*JSTEP*DTIME C IFIXD: DIRECTION OF ACCERALATION INPUT C =0: HORIZONTAL & VERTICAL, 1:VERTICAL, 2:HORIZONTAL C MITER: MAXIMUM NUMBER OF ITERATIONS (20 IS RECOMENNDED) C IPRED: TIME-STEPS FOR GRAVITY FORCE ONLY C READ(5,100) NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC, * IFIXD,MITER,IPRED WRITE(6,200) NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC, * IFIXD,MITER,IPRED C DTIME: TIME-STEP INCREMENT C DTEND: END-TIME C DTREC: TIME-STEP AT ACCERALATION RECORD C DELTA, GAMMA: CONSTANTS IN NEWMARK METHOD C DELTA=0.25, GAMMA=0.5 C AZERO, BZERO, OMEGA C FUNCTS=AZERO+BZERO*SIN(OMEGA*STEP*DTIME) C TOLER: CONVERGENCE CONSTANT (TOLER=0.0001) C READ(5,101) DTIME,DTEND,DTREC,DELTA,GAAMA,AZERO,BZERO, * OMEGA,TOLER WRITE(6,201) DTIME,DTEND,DTREC,DELTA,GAAMA,AZERO,BZERO, * OMEGA,TOLER C C *** SELECTED TIME-STEPS FOR OUTPUT READ(5,100) (LOUTD(N),N=1,NOUTP) WRITE(6,208) (LOUTD(N),N=1,NOUTP) DO 2 N=1,NSTEP 2 NOUTD(N)=0 DO 4 L=1,NOUTP N=LOUTD(L) 4 NOUTD(N)=1 C *** SELECTED NODES AND ELEMENTS FOR OUTPUT WRITE(6,202) READ(5,100) (IREQD(I),I=1,NREQD) WRITE(6,203) (IREQD(I),I=1,NREQD) READ(5,100) (IREQS(I),I=1,NREQS) WRITE(6,204) (IREQS(I),I=1,NREQS) C IF(IFUNC.NE.0) GO TO 50 C *** READ ACCELEROGRAM DATA C IFIXD=0: H & V ACC, 1: V, 2: H C AFACT=DTREC/DTIME IF(IFIXD-1) 40,42,44 40 WRITE(6,209) DTREC READ(5,103) (ACCEH(N),N=1,NACCE) WRITE(6,103) (ACCEH(N),N=1,NACCE) WRITE(6,210) DTREC READ(5,103) (ACCEV(N),N=1,NACCE) WRITE(6,103) (ACCEV(N),N=1,NACCE) GO TO 50 42 WRITE(6,210) DTREC READ(5,103) (ACCEV(N),N=1,NACCE) WRITE(6,103) (ACCEV(N),N=1,NACCE) GO TO 50 44 WRITE(6,209) DTREC READ(5,103) (ACCEH(N),N=1,NACCE) WRITE(6,103) (ACCEH(N),N=1,NACCE) 50 CONTINUE C 100 FORMAT(16I5) 101 FORMAT(9F10.4) 103 FORMAT(10F8.4) 200 FORMAT(/'TIME STEPPING PARAMETERS'/ * /2X,'NSTEP=',I5,12X,'NOUTP=',I5,12X,'NREQD=',I5 * /2X,'NREQS=',I5,12X,'NACCE=',I5,12X,'IFUNC=',I5 * /2X,'IFIXD=',I5,12X,'MITER=',I5,12X,'IPRED=',I5) 201 FORMAT(/2X,'DTIME=',G12.4,5X,'DTEND=',G12.4,5X,'DTREC=',G12.4 * /2X,'DELTA=',G12.4,5X,'GAAMA=',G12.4,5X,'AZERO=',G12.4 * /2X,'BZERO=',G12.4,5X,'OMEGA=',G12.4,5X,'TOLER=',G12.4) 202 FORMAT(/'SELECTIVE OUTPUT' ) 203 FORMAT(/2X,' NODE ',10I5) 204 FORMAT(/2X,' ELEM ',10I5) 205 FORMAT(/'SELECTIVE TIME-STEP FOR MICRO-AVS=',I5) 208 FORMAT(/'SELECTIVE TIME-STEPS',30(/10I5)) 209 FORMAT(/'HORIZONTAL ACCELERATION ORDINATES AT',F9.4,2X,'SEC'/) 210 FORMAT(/'VERTICAL ACCELERATION ORDINATES AT',F9.4,2X,'SEC'/) RETURN END C * * * * * * * * * * * * SUBROUTINE PREVOS C *** INITIAL STATE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA20/IYIEL(2150),IUNLO(2150) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) DIMENSION STI(20,6) C C KTENS=0 DO 2 IE=1,NELEM EFFST(IE)=0.D0 EPSTN(IE)=0.D0 IYIEL(IE)=0 IUNLO(IE)=0 C ITENS(IE)=0 DO 4 I=1,6 STRSG(IE,I)=0.D0 STRAG(IE,I)=0.D0 STRIN(IE,I)=0.D0 4 CONTINUE 2 CONTINUE DO 6 IT=1,NSIZE DISPI(IT)=0.D0 DISPT(IT)=0.D0 DISPL(IT)=0.D0 VELOI(IT)=0.D0 6 CONTINUE C *** INPUT INITIAL STRESS GROUPS WRITE(6,200) IF(NIST.LE.0) GO TO 20 DO 10 IS=1,NIST READ(5,100) IS1,(STI(IS,I),I=1,6) WRITE(6,100) IS1,(STI(IS,I),I=1,6) 10 CONTINUE DO 12 IE=1,NELEM IS=IST(IE) DO 14 I=1,4 14 STRIN(IE,I)=STI(IS,I) 12 CONTINUE 20 CONTINUE C *** INPUT INITIAL STRESSES IN EACH ELEMENT IF(NPREV.EQ.0) RETURN OPEN(7,FILE='DASTRE',STATUS='UNKNOWN') DO 30 I=1,4 READ(7,101) (STRIN(IE,I),IE=1,NELEM) DO 32 IE=1,NELEM 32 STRIN(IE,I)=-STRIN(IE,I) WRITE(6,101) (STRIN(IE,I),IE=1,NELEM) 30 CONTINUE CLOSE(7) 100 FORMAT(I2,6F8.3) 101 FORMAT(10F8.3) 200 FORMAT(/'INITIAL STRESSES'/) RETURN END C * * * * * * * * * * * * SUBROUTINE LOADPL C *** STANDARD LOAD ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) DIMENSION ELOAD(8) C DO 2 IV=1,NSIZE 2 FORCE(IV)=0.D0 C *** LOAD DATA C NFL: NUMBER OF NODAL POINT LOADS C IGRAV=0: NOT CONSIDER GRAVITY LOADS, 1:CONSIDER C READ (5,100) NFL,IGRAV WRITE(6,200) NFL,IGRAV C *** READ NODAL POINT LOADS C NOD: NODE NO. C IXY: DIRECTION OF LOAD (1=HORIZONTAL, 2=VERTICAL, 3=MOMENT) C FLO: LOAD C IF(NFL.EQ.0) GO TO 12 WRITE(6,201) DO 10 IL=1,NFL READ(5,101) NOD,IXY,FLO WRITE(6,101) NOD,IXY,FLO IT=(NOD-1)*3+IXY IT=LY(IT) IF(IT.EQ.NN) GO TO 10 FORCE(IT)=FORCE(IT)+FLO 10 CONTINUE 12 CONTINUE C *** GIVE GRAVITY LOADS IF(IGRAV.EQ.0) GO TO 22 GRAVZ=-9.81D0 DO 20 IE=1,NELEM KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 CALL TRUSS(IE) NNODE=2 NEVAB=6 GO TO 30 32 CALL BEAM(IE) NNODE=2 NEVAB=6 GO TO 30 34 CALL PLANES(IE,KOL2) NNODE=4 NEVAB=8 MM=IJK(IE,4) IF(MM.LE.0) NNODE=3 IF(MM.LE.0) NEVAB=6 30 CONTINUE C MT=MAT(IE) RHO= PRP(MT,4) IF(RHO.LE.0.D0) RHO=0.D0 DMASS=AES(IE)*RHO*GRAVZ/DFLOAT(NNODE) NEVA2=NEVAB/2 IF(KOL2.EQ.2) NEVA2=2 DO 36 I=1,NEVA2 I2=I*2 IF(KOL2.EQ.2) I2=3*I-1 IT=LLL(IE,I2) IF(IT.EQ.NN) GO TO 36 FORCE(IT)=FORCE(IT)+DMASS 36 CONTINUE 20 CONTINUE 22 CONTINUE C WRITE(6,202) C WRITE(6,203) (FORCE(IV),IV=1,NSIZE) C 100 FORMAT(10I5) 101 FORMAT(2I5,F10.2) 200 FORMAT(/'NODAL LOADS=',I5,' GRAVITY LOADS=',I5/) 201 FORMAT(/'NODAL LOAD'/' NODE X-Y LOAD') 202 FORMAT(/'FORCE') 203 FORMAT(10E8.2) RETURN END C * * * * * * * * * * * * SUBROUTINE LUMASS C *** CALCULATES LUMPED MASS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) DIMENSION EMASS(8,8) C DO 5 IV=1,NSIZE DO 5 JV=1,NSIZE XMASS(IV,JV)=0.D0 5 DAMPI(IV,JV)=0.D0 C DO 100 IE=1,NELEM KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 CALL TRUSS(IE) NNODE=2 NEVAB=4 GO TO 30 32 CALL BEAM(IE) NNODE=2 NEVAB=6 GO TO 30 34 CALL PLANES(IE,KOL2) NNODE=4 NEVAB=8 MM=IJK(IE,4) IF(MM.LE.0) NNODE=3 IF(MM.LE.0) NEVAB=6 30 CONTINUE C MT=MAT(IE) RHO= PRP(MT,4) IF(RHO.LE.0.D0) RHO=0.001D0 ALFA=PRP(MT,8) DO 12 I=1,NEVAB DO 12 J=1,NEVAB 12 EMASS(I,J)=0.D0 C DMASS=AES(IE)*RHO/DFLOAT(NNODE) DO 14 I=1,NEVAB IF(KOL2.EQ.2.AND.I.EQ.3) DMASS=0.001D0 IF(KOL2.EQ.2.AND.I.EQ.6) DMASS=0.001D0 EMASS(I,I)=DMASS 14 CONTINUE C DO 20 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 20 DO 22 J=1,NEVAB JT=LLL(IE,J) IF(JT.EQ.NN) GO TO 22 XMASS(IT,JT)=XMASS(IT,JT)+EMASS(I,J) DAMPI(IT,JT)=DAMPI(IT,JT)+EMASS(I,J)*ALFA 22 CONTINUE 20 CONTINUE 100 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE GSTIFF C *** EVALUATES STIFFNESS MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) C DO 4 IV=1,NSIZE DO 4 JV=1,NSIZE 4 STIFF(IV,JV)=0.D0 C DO 100 IE=1,NELEM MT=MAT(IE) BETA=PRP(MT,9) KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 NEVAB=4 GO TO 30 32 NEVAB=6 GO TO 30 34 NEVAB=8 IF(IJK(IE,4).EQ.0) NEVAB=6 30 CONTINUE DO 20 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 20 DO 22 J=1,NEVAB JT=LLL(IE,J) IF(JT.EQ.NN) GO TO 22 STIFF(IT,JT)=STIFF(IT,JT)+EKK(IE,I,J) DAMPI(IT,JT)=DAMPI(IT,JT)+EKK(IE,I,J)*BETA 22 CONTINUE 20 CONTINUE 100 CONTINUE C *** CALCULATES K-STAR MATRICES CONSC=1.D0/(DTIME*DTIME*DELTA) CONSD=GAAMA/(DTIME*DELTA) DO 36 IV=1,NSIZE DO 36 JV=1,NSIZE 36 STIFS(IV,JV)=XMASS(IV,JV)*CONSC+DAMPI(IV,JV)*CONSD+STIFF(IV,JV) C CALL MATINV(STIFS,NN,NSIZE,0,DET,IND) IF(IND.EQ.1) GO TO 80 WRITE(6,200) IND 200 FORMAT(//'INDEX (STIFS)=',I3) STOP 80 RETURN END C * * * * * * * * * * * * SUBROUTINE IMPEXP(IITER,ISTEP) C *** GENERATES PARTIAL EFFECTIVE LOAD VECTOR IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA12/ACCEH(6550),ACCEV(6550) COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA18/ACCEI(2220),ACCEL(2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) DIMENSION TACCE(2220),ACCEJ(2220),ACCEK(2220) C IF(ISTEP.GT.1.OR.IITER.GT.1) GO TO 100 CONSA=DTIME*DTIME*(0.5D0-DELTA) CONSB=DTIME*(1.D0-GAAMA) CONSC=1.D0/(DTIME*DTIME*DELTA) CONSD=GAAMA/(DTIME*DELTA) C DO 10 IP=1,NPOIN IT=(IP-1)*3+1 IT=LY(IT) IF(IT.EQ.NN) GO TO 12 ACCEI(IT)=1.D0 ACCEL(IT)=0.D0 12 CONTINUE IT=(IP-1)*3+2 IT=LY(IT) ACCEI(IT)=0.D0 ACCEL(IT)=1.D0 10 CONTINUE C *** CALCULATES VECTORS FOR HORIZONTAL AND VERTICAL EXCITATION DO 30 IV=1,NSIZE ACCEK(IV)=0.D0 ACCEJ(IV)=0.D0 DISPL(IV)=0.D0 VELOL(IV)=0.D0 DO 32 JV=1,NSIZE ACCEK(IV)=ACCEK(IV)+XMASS(IV,JV)*ACCEL(JV) ACCEJ(IV)=ACCEJ(IV)+XMASS(IV,JV)*ACCEI(JV) DISPL(IV)=DISPL(IV)+STIFF(IV,JV)*DISPI(JV) VELOL(IV)=VELOL(IV)+DAMPI(IV,JV)*VELOI(JV) 32 CONTINUE 30 CONTINUE C *** CALCULATES INITIAL ACCELERATION DO 14 IV=1,NSIZE 14 TACCE(IV)=FORCE(IV)-DISPL(IV)-VELOL(IV) DO 16 IV=1,NSIZE DO 16 JV=1,NSIZE 16 STIFF(IV,JV)=XMASS(IV,JV) C CALL MATINV(STIFF,NN,NSIZE,0,DET,IND) IF(IND.EQ.1) GO TO 18 WRITE(6,202) IND 202 FORMAT(//'INDEX (XMASS)=',I5) STOP C 18 CONTINUE DO 34 IV=1,NSIZE ACCEI(IV)=0.D0 DO 34 JV=1,NSIZE 34 ACCEI(IV)=ACCEI(IV)+STIFF(IV,JV)*TACCE(JV) WRITE(6,200) WRITE(6,201) (ACCEI(IV),IV=1,NSIZE) C 100 CONTINUE IF(IITER.GT.1) GO TO 110 C *** CALCULATES PREDICTED DISPLACEMENT AND VELOCITY VECTOR DO 20 IV=1,NSIZE DISPI(IV)=DISPI(IV)+DTIME*VELOI(IV)+CONSA*ACCEI(IV) VELOI(IV)=VELOI(IV)+CONSB*ACCEI(IV) DISPT(IV)=DISPI(IV) VELOT(IV)=VELOI(IV) ACCEI(IV)=CONSC*(DISPT(IV)-DISPI(IV)) 20 CONTINUE C *** CALCULATES LOAD VECTORS FACTS=FUNCTS(ISTEP) FACTH=FUNCTA(ACCEH,AFACT,ISTEP) FACTV=FUNCTA(ACCEV,AFACT,ISTEP) C 110 CONTINUE IF(ISTEP.EQ.1) GO TO 120 C *** CALCULATES PARTIAL EFFECTIVE LOAD VECTOR DO 36 IV=1,NSIZE VELOL(IV)=0.D0 DO 36 JV=1,NSIZE 36 VELOL(IV)=VELOL(IV)+DAMPI(IV,JV)*VELOT(JV) 120 DO 40 IV=1,NSIZE IF(IFUNC.NE.0) GO TO 42 IF(IFIXD.EQ.2) * DISPL(IV)=-VELOL(IV)-FACTH*ACCEJ(IV)+FORCE(IV) IF(IFIXD.EQ.1) * DISPL(IV)=-VELOL(IV)-FACTV*ACCEK(IV)+FORCE(IV) IF(IFIXD.EQ.0) * DISPL(IV)=-VELOL(IV)-FACTH*ACCEJ(IV)-FACTV*ACCEK(IV)+FORCE(IV) IF(IFUNC.EQ.0) GO TO 40 42 CONTINUE DISPL(IV)=-VELOL(IV)+FORCE(IV)*FACTS 40 CONTINUE 200 FORMAT(/'INITIAL ACCELERATION '/) 201 FORMAT(10E8.2) RETURN END C * * * * * * * * * * * * SUBROUTINE RESEPL(IITER,ISTEP) C *** EVALUATES RESIDUAL FORCES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA18/ACCEI(2220),ACCEL(2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) COMMON/LA20/IYIEL(2150),IUNLO(2150) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) DIMENSION AVECT(4),DVECT(4),DEVIA(4),ELOAD(8),TU(2220), * STRAN(6),STRES(6),SIGMA(6),SGTOT(6),DESIG(6), * BVECT(4),EVECT(4),EU(6) C DO 2 IV=1,NSIZE RESID(IV)=0.D0 TU(IV)=DISPT(IV) 2 CONTINUE TU(NN)=0.D0 IYIEL(IE)=0 PI1(IE)=0.D0 PI2(IE)=0.D0 C DO 10 IE=1,NELEM KOL2=K2(IE) MT=MAT(IE) YOUNG=PRP(MT,1) POISS=PRP(MT,2) THICK=PRP(MT,3) COHES=PRP(MT,5) FRICT=PRP(MT,6)*0.017453292D0 DAILT=PRP(MT,7)*0.017453292D0 UNIAX=COHES*DCOS(FRICT) HARDS=0.D0 PREYS=UNIAX C DO 8 I=1,6 8 STRES(I)=0.D0 GO TO (11,12,11,14,14,14),KOL2 11 NEVAB=4 NSTRE=1 GO TO 20 12 NEVAB=6 NSTRE=6 GO TO 20 14 NEVAB=8 NSTRE=4 IF(IJK(IE,4).EQ.0) NEVAB=6 20 CONTINUE C CALL LINGNL(IE,KOL2,TU,YOUNG,POISS,STRAN,STRES) DO 22 I=1,NSTRE 22 STRAG(IE,I)=STRAG(IE,I)+STRAN(I) IF(ISTEP.GT.1.OR.IITER.GT.1) GO TO 23 DO 24 I=1,NSTRE 24 STRES(I)=STRES(I)+STRIN(IE,I) 23 CONTINUE DO 25 I=1,NSTRE DESIG(I)=STRES(I) SIGMA(I)=STRSG(IE,I)+STRES(I) 25 CONTINUE IF(KOL2.LE.4) GO TO 30 C IF(ISTEP.LE.IPRED) GO TO 27 DO 28 I=1,NSTRE 28 SIGMA(I)=SIGMA(I)-PW(IE,I) 27 CONTINUE C CALL INVAR(DEVIA,MT,SINT3,STEFF,SIGMA,THETA,VARJ2,YIELD,KOL2) ESPRE=EFFST(IE)-PREYS C ESPRE.GE.0: THE ELEMENT YIELDED AT PRECEDING STEP IF(ESPRE.GE.0.D0) GO TO 32 ESCUR=YIELD-PREYS C ESCUR.GT.0: THE ELEMENT HAS YIELDED AT PRESENT STEP IF(ESCUR.LE.0.D0) IYIEL(IE)=0 IF(ESCUR.LE.0.D0) GO TO 38 RFACT=ESCUR/(YIELD-EFFST(IE)) GO TO 34 32 CONTINUE ESCUR=YIELD-EFFST(IE) C ESCUR.LE.0: THE ELEMENT DOES NOT YIELD AT PRESENT STEP IF(ESCUR.LE.0.D0) IYIEL(IE)=0 IF(ESCUR.LE.0.D0) GO TO 38 RFACT=1.D0 34 CONTINUE IYIEL(IE)=1 C SX=-SIGMA(1) SY=-SIGMA(2) TA=-SIGMA(3) SZ=-SIGMA(4) CALL MOHRCO(IE,SX,SY,TA,F1,S,P2,B0) PI1(IE)=F1 IF(F1.GE.0.D0) IYIEL(IE)=1 CALL PRINCE(SX,SY,TA,S1,S3,T) F2=(S1-SZ)-(2.D0*COHES+(S1+SZ)*DSIN(FRICT)) PI2(IE)=F2 IF(F2.GE.0.D0) IYIEL(IE)=2 C MSTEP=ESCUR*8.D0/UNIAX+1.D0 IF(MSTEP.GT.10) MSTEP=10 ASTEP=MSTEP C REDUC=1.D0-RFACT DO 36 I=1,NSTRE SGTOT(I)=STRSG(IE,I)+REDUC*STRES(I) IF(ISTEP.GT.IPRED) SGTOT(I)=SGTOT(I)-PW(IE,I) STRES(I)=RFACT*STRES(I)/ASTEP 36 CONTINUE C DO 40 JSTEP=1,MSTEP CALL INVAR (DEVIA,MT,SINT3,STEFF,SGTOT,THETA,VARJ2,YIELD,KOL2) CALL YIELDF(AVECT,BVECT,DEVIA,SGTOT,SINT3,STEFF,THETA, * VARJ2,FRICT,DAILT,KOL2) CALL FLOWPL(AVECT,BVECT,ABETA,DVECT,EVECT,YOUNG,POISS,THICK,KOL2) AGASH=0.D0 DO 42 I=1,NSTRE 42 AGASH=AGASH+AVECT(I)*STRES(I) DLAMD=AGASH*ABETA IF(DLAMD.LT.0.D0) DLAMD=0.D0 BGASH=0.D0 DO 44 I=1,NSTRE BGASH=BGASH+AVECT(I)*SGTOT(I) 44 SGTOT(I)=SGTOT(I)+STRES(I)-DLAMD*EVECT(I) EPSTN(IE)=EPSTN(IE)+DLAMD*BGASH/YIELD 40 CONTINUE C CALL INVAR(DEVIA,MT,SINT3,STEFF,SGTOT,THETA,VARJ2,YIELD,KOL2) CURYS=UNIAX BRING=1.D0 IF(YIELD.GT.CURYS) BRING=CURYS/YIELD DO 46 I=1,NSTRE 46 STRSG(IE,I)=BRING*SGTOT(I) EFFST(IE)=BRING*YIELD IF(ISTEP.LE.IPRED) GO TO 47 DO 48 I=1,NSTRE 48 STRSG(IE,I)=STRSG(IE,I)+PW(IE,I) 47 CONTINUE C *** ALTERNATIVE LOCATION OF STRESS REDUCTION LOOP TERMINATION CARD GO TO 50 38 CONTINUE DO 54 I=1,NSTRE 54 STRSG(IE,I)=SIGMA(I) IF(ISTEP.LE.IPRED) GO TO 53 DO 55 I=1,NSTRE 55 STRSG(IE,I)=STRSG(IE,I)+PW(IE,I) 53 CONTINUE GO TO 52 C 30 CONTINUE DO 56 I=1,NSTRE 56 STRSG(IE,I)=STRSG(IE,I)+DESIG(I) 52 CONTINUE EFFST(IE)=YIELD 50 CONTINUE C *** NO-TENSION CALL NOTENS(IE,KOL2) C C *** CALCULATE THE EQUIVALENT NODAL FORCES DVOLU=AES(IE) NSTR1=NSTRE IF(KOL2.EQ.2) GO TO 60 IF(KOL2.GE.4) NSTR1=3 DO 62 I=1,NEVAB ELOAD(I)=0.D0 DO 62 J=1,NSTR1 62 ELOAD(I)=ELOAD(I)+BMX(IE,J,I)*STRSG(IE,J)*DVOLU GO TO 64 60 CONTINUE DO 66 I=1,NEVAB IT=LLL(IE,I) 66 EU(I)=TU(IT) DO 68 I=1,NEVAB ELOAD(I)=0.D0 DO 68 J=1,NEVAB 68 ELOAD(I)=ELOAD(I)+EKK(IE,I,J)*EU(J) 64 CONTINUE DO 70 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 70 RESID(IT)=RESID(IT)+ELOAD(I) 70 CONTINUE 10 CONTINUE RETURN END C * * * * * * * * * * SUBROUTINE ITRATE(IITER) C *** CALCULATES INCREMENT IN DISPLACEMENT AND APPLIES CONVERGENCE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA13/FORCE(2220),EPSTN(2150),EFFST(2150) COMMON/LA14/XMASS(2220,2220),STIFF(2220,2220),STIFS(2220,2220) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA18/ACCEI(2220),ACCEL(2220) COMMON/LA19/DAMPI(2220,2220),RESID(2220) DIMENSION TVECT(2220),DISPD(2220) C NCHEK=0 DO 2 IV=1,NSIZE ACCEL(IV)=0.D0 DO 2 JV=1,NSIZE 2 ACCEL(IV)=ACCEL(IV)+XMASS(IV,JV)*ACCEI(JV) C *** CALCULATES TOTAL EFFECTIVE LOAD VECTOR DO 10 IV=1,NSIZE 10 TVECT(IV)=DISPL(IV)-ACCEL(IV)-RESID(IV) C *** CALCULATES DELTA DISPLACEMENT DO 11 IV=1,NSIZE DISPD(IV)=0.D0 DO 11 JV=1,NSIZE 11 DISPD(IV)=DISPD(IV)+STIFS(IV,JV)*TVECT(JV) C *** APPLIES CONVERGENCE SUMPP=0.D0 SUMPQ=0.D0 DO 12 IV=1,NSIZE DISPP=DISPD(IV) DISPQ=DISPT(IV)+DISPP DISPT(IV)=DISPQ SUMPP=SUMPP+DISPP*DISPP SUMPQ=SUMPQ+DISPQ*DISPQ 12 CONTINUE CONSD=DTIME*(1.D0-GAAMA) CONSF=1.D0/(DTIME*DTIME*DELTA) DO 14 IV=1,NSIZE ACCEI(IV)=CONSF*(DISPT(IV)-DISPI(IV)) 14 VELOT(IV)=VELOI(IV)+CONSD*ACCEI(IV) SUMPP=DSQRT(SUMPP/SUMPQ) WRITE(6,200) IITER,SUMPP,SUMPQ IF(SUMPP.GT.TOLER) GO TO 30 NCHEK=1 GO TO 20 30 IF(IITER.LT.MITER) GO TO 32 20 DO 22 IV=1,NSIZE VELOI(IV)=VELOT(IV) 22 DISPI(IV)=DISPT(IV) 32 CONTINUE 200 FORMAT(I5,2E10.3) RETURN END C * * * * * * * * * * * * SUBROUTINE LINGNL(IE,KOL2,TU,YOUNG,POISS,STRAN,STRES) C *** ELASTIC STRAINS AND STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) DIMENSION EU(8),TU(2220),STRAN(6),STRES(6) C DO 5 I=1,6 STRAN(I)=0.D0 5 STRES(I)=0.D0 GO TO (11,12,11,14,14,14),KOL2 C *** TRUSS & TEXTILE 11 CONTINUE DO 20 I=1,4 IT=LLL(IE,I) 20 EU(I)=TU(IT) STRAN(1)=0.D0 DO 22 I=1,4 22 STRAN(1)=STRAN(1)+BMX(IE,1,I)*EU(I) STRAN(1)=STRAN(1)-STRAG(IE,1) STRES(1)=STRAN(1)*YOUNG GO TO 10 C *** BEAM 12 CONTINUE DO 30 I=1,6 IT=LLL(IE,I) 30 EU(I)=TU(IT) DO 32 I=1,6 STRES(I)=0.D0 DO 32 J=1,6 32 STRES(I)=STRES(I)+DBM(IE,I,J)*EU(J) DO 34 I=1,6 STRAN(I)=0.D0 34 STRES(I)=STRES(I)-STRSG(IE,I) C *** STRES(I) MEAN PARALLEL & NORMAL NODAL FORCES C AND BENDING MOMENT AT EACH NODE GO TO 10 C *** PLANE-STRAIN AND INTERFACE 14 CONTINUE NEVAB=8 IF(IJK(IE,4).LE.0) NEVAB=6 DO 40 I=1,NEVAB IT=LLL(IE,I) 40 EU(I)=TU(IT) DO 42 I=1,4 STRAN(I)=0.D0 DO 42 J=1,NEVAB 42 STRAN(I)=STRAN(I)+BMX(IE,I,J)*EU(J) STRAN(4)=0.D0 DO 44 I=1,4 44 STRAN(I)=STRAN(I)-STRAG(IE,I) DO 46 I=1,4 STRES(I)=0.D0 DO 46 J=1,4 46 STRES(I)=STRES(I)+DMX(IE,I,J)*STRAN(J) IF(STRES(4).GT.0.D0) STRES(4)=0.D0 10 RETURN END C * * * * * * * * * * * * SUBROUTINE INVAR(DEVIA,MT,SINT3,STEFF,STEMP,THETA,VARJ2, * YIELD,KOL2) C *** STRESS INVARIANTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) DIMENSION DEVIA(4),STEMP(4) C PHIRA=PRP(MT,6)*0.017453292 SNPHI=DSIN(PHIRA) COPHI=DCOS(PHIRA) GO TO (10,10,10,10,15,16),KOL2 15 CONTINUE C *** INVARIANTS ROOT3=1.73205080757 SMEAN=(STEMP(1)+STEMP(2)+STEMP(4))/3.D0 DEVIA(1)=STEMP(1)-SMEAN DEVIA(2)=STEMP(2)-SMEAN DEVIA(3)=STEMP(3) DEVIA(4)=STEMP(4)-SMEAN VARJ2=DEVIA(3)*DEVIA(3)+0.5D0*(DEVIA(1)*DEVIA(1)+ * DEVIA(2)*DEVIA(2)+DEVIA(4)*DEVIA(4)) VARJ3=DEVIA(4)*(DEVIA(4)*DEVIA(4)-VARJ2) STEFF=DSQRT(VARJ2) IF(VARJ2.EQ.0.D0.OR.STEFF.EQ.0.D0) GO TO 5 SINT3=-2.5980762113*VARJ3/(VARJ2*STEFF) GO TO 6 5 SINT3=0.D0 6 CONTINUE IF(SINT3.LT.-1.) SINT3=-1.D0 IF(SINT3.GT. 1.) SINT3=1.D0 THETA=DASIN(SINT3)/3.D0 C *** MOHR-COULOMB YIELD=SMEAN*SNPHI+STEFF*(DCOS(THETA)-DSIN(THETA)*SNPHI/ROOT3) GO TO 10 C *** COULOMB 16 CONTINUE TAU=STEMP(3) TAU=DABS(TAU) SIG=-STEMP(2) IF(SIG.LT.0.D0) SIG=0.D0 YIELD=TAU*COPHI-SIG*SNPHI 10 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE YIELDF(AVECT,BVECT,DEVIA,STEMP,SINT3,STEFF, * THETA,VARJ2,FRICT,DAILT,KOL2) C *** SELECTS YIELD FUNCTION AND CALCULATES VECTOR 'AVECT' IMPLICIT REAL*8(A-H,O-Z) DIMENSION AVECT(4),BVECT(4),DEVIA(4),STEMP(4), * VECA1(4),VECA2(4),VECA3(4) C GO TO (50,50,50,50,55,56),KOL2 55 CONTINUE IF(STEFF.EQ.0.D0) RETURN TANTH=DTAN(THETA) SINTH=DSIN(THETA) COSTH=DCOS(THETA) COST3=DCOS(3.D0*THETA) ROOT3=1.73205080757 C *** VECTOR A1 VECA1(1)=1.D0 VECA1(2)=1.D0 VECA1(3)=0.D0 VECA1(4)=1.D0 C *** VECTOR A2 DO 2 ISTR1=1,4 2 VECA2(ISTR1)=DEVIA(ISTR1)/(2.D0*STEFF) VECA2(3)=DEVIA(3)/STEFF C *** VECTOR A3 VECA3(1)=DEVIA(2)*DEVIA(4)+VARJ2/3.D0 VECA3(2)=DEVIA(1)*DEVIA(4)+VARJ2/3.D0 VECA3(3)=-2.D0*DEVIA(3)*DEVIA(4) VECA3(4)=DEVIA(1)*DEVIA(2)-DEVIA(3)**2+VARJ2/3.D0 C *** MOHR-COULOMB CONS1=DSIN(FRICT)/3.D0 ABTHE=DABS(THETA*57.29577951308) IF(ABTHE.LT.29.) GO TO 4 CONS3=0.D0 PLUMI=1.D0 IF(THETA.GT.0.D0) PLUMI=-1.D0 CONS2=0.5D0*(ROOT3+PLUMI*CONS1*ROOT3) GO TO 6 4 TANT3=DTAN(3.D0*THETA) CONS2=COSTH*((1.D0+TANTH*TANT3)+CONS1*(TANT3-TANTH)*ROOT3) CONS3=(ROOT3*SINTH+3.D0*CONS1*COSTH)/(2.D0*VARJ2*COST3) 6 CONTINUE DO 8 ISTR1=1,4 8 AVECT(ISTR1)=CONS1*VECA1(ISTR1)+CONS2*VECA2(ISTR1)+ * CONS3*VECA3(ISTR1) C *** NON-ASSOCIATED FLOW RULE CONS1=DSIN(DAILT)/3.D0 IF(ABTHE.LT.29.) GO TO 10 CONS3=0.D0 PLUMI=1.D0 IF(THETA.GT.0.D0) PLUMI=-1.D0 CONS2=0.5D0*(ROOT3+PLUMI*CONS1*ROOT3) GO TO 12 10 TANT3=DTAN(3.*THETA) CONS2=COSTH*((1.D0+TANTH*TANT3)+CONS1*(TANT3-TANTH)*ROOT3) CONS3=(ROOT3*SINTH+3.D0*CONS1*COSTH)/(2.D0*VARJ2*COST3) 12 CONTINUE DO 14 I=1,4 14 BVECT(I)=CONS1*VECA1(I)+CONS2*VECA2(I)+CONS3*VECA3(I) GO TO 50 C *** COULOMB 56 CONTINUE TAU=STEMP(3) SIG=STEMP(2) C1=1.D0 IF(TAU.LT.0.D0) C1=-1.D0 C2=-1.D0 IF(SIG.GT.0.D0) C2=0.D0 AVECT(1)=0.D0 AVECT(2)=C2*(-DSIN(FRICT)) AVECT(3)=C1*DCOS(FRICT) AVECT(4)=0.D0 BVECT(1)=0.D0 BVECT(2)=C2*(-DSIN(DAILT)) BVECT(3)=C1*DCOS(DAILT) BVECT(4)=0.D0 50 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE FLOWPL(AVECT,BVECT,ABETA,DVECT,EVECT, * YOUNG,POISS,THICK,KOL2) C *** EVALUATES THE PLASTIC D VECTOR IMPLICIT REAL*8(A-H,O-Z) DIMENSION AVECT(4),DVECT(4),BVECT(4),EVECT(4) C GO TO (10,10,10,10,15,16),KOL2 C *** COULOMB 16 CONTINUE G=YOUNG YOUNG=POISS POISS=THICK GO TO 20 C *** MOHR-COULOMB 15 CONTINUE G=0.5D0*YOUNG/(1.D0+POISS) 20 CONTINUE FMUL1=YOUNG/(1.D0+POISS) FMUL2=YOUNG*POISS*(AVECT(1)+AVECT(2)+AVECT(4))/ * ((1.D0+POISS)*(1.D0-2.D0*POISS)) DVECT(1)=FMUL1*AVECT(1)+FMUL2 DVECT(2)=FMUL1*AVECT(2)+FMUL2 DVECT(3)=AVECT(3)*G DVECT(4)=FMUL1*AVECT(4)+FMUL2 FMUL4=YOUNG*POISS*(BVECT(1)+BVECT(2)+BVECT(4))/ * ((1.D0+POISS)*(1.D0-2.D0*POISS)) EVECT(1)=FMUL1*BVECT(1)+FMUL4 EVECT(2)=FMUL1*BVECT(2)+FMUL4 EVECT(3)=BVECT(3)*G EVECT(4)=FMUL1*BVECT(4)+FMUL4 DENOM=0.D0 DO 22 I=1,4 22 DENOM=DENOM+AVECT(I)*EVECT(I) ABETA=1.D0/DENOM 10 CONTINUE RETURN END C * * * * * * * * * * SUBROUTINE MOHRCO(L,SX,SY,TA,F,S,P2,B0) IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) C M=MAT(L) CC=PRP(M,5) FAI=PRP(M,6)*3.14159D0/180.D0 SFI=DSIN(FAI) CFI=DCOS(FAI) P1=SX+SY P2=SX-SY A1=P1*SFI+2.D0*CC*CFI B0=P2*P2+4.D0*TA*TA IF(B0.LT.0.0001D0) B0=0.0001D0 B1=DSQRT(B0) F=B1-A1 S=A1/B1 RETURN END C * * * * * * * * * * SUBROUTINE NOTENS(IE,KOL2) C *** NO-TENSION ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) DIMENSION SIG(6) C DO 2 I=1,4 2 SIG(I)=-STRSG(IE,I) JTENS(IE)=0 MT=MAT(IE) TENST=-PRP(MT,10) GO TO (10,10,11,10,15,16),KOL2 11 CONTINUE SG=SIG(1) IF(SG.LE.0.D0) GO TO 18 SIG(1)=0.D0 KTENS=KTENS+1 ITENS(KTENS)=IE JTENS(IE)=KTENS GO TO 18 15 CONTINUE SX=SIG(1) SY=SIG(2) TA=SIG(3) CALL PRINCE(SX,SY,TA,S1,S3,TH) IF(S3.GE.TENST) GO TO 18 KTENS=KTENS+1 ITENS(KTENS)=IE JTENS(IE)=KTENS C=1.D0 IF(TA.LT.0.D0) C=-1.D0 C IF(S1.GT.TENST) GO TO 20 SIG(1)=TENST SIG(2)=TENST SIG(3)=0.D0 GO TO 18 20 CONTINUE IF(SY.GE.SX) GO TO 30 C SX.GT.SY IF(SX) 50,50,40 40 CONTINUE SY=TENST+(S1-SX) SIG(2)=SY SIG(1)=SX GO TO 35 30 CONTINUE C SY.GE.SX IF(SY) 52,52,42 42 CONTINUE SX=TENST+(S1-SY) SIG(1)=SX SIG(2)=SY 35 CONTINUE A=((S1-TENST)*0.5D0)**2-((SX-SY)*0.5D0)**2 IF(A.LT.0.000001D0) A=0.000001D0 SIG(3)=C*DSQRT(A) GO TO 18 50 C1=1.D0 GO TO 55 52 C1=-1.D0 55 CONTINUE R=(S1-S3)*0.5D0 B2=R*R-TA*TA IF(B2.GT.0.D0) GO TO 57 T=0.5D0*3.14159D0 GO TO 58 57 CONTINUE B=DSQRT(B2) AT=TA/B T=DATAN(AT) 58 CONTINUE S=(S1-TENST)*0.5D0*DCOS(T) SIG(1)=(S1+TENST)*0.5D0+C1*S SIG(2)=(S1+TENST)*0.5D0-C1*S SIG(3)=(S1-TENST)*0.5D0*DSIN(T)*C GO TO 18 16 CONTINUE SG=SIG(2) IF(SG.GE.0.D0) GO TO 18 KTENS=KTENS+1 ITENS(KTENS)=IE JTENS(IE)=KTENS SIG(2)=0.001D0 TA=SIG(3) C=1.D0 IF(TA.LT.0.D0) C=-1.D0 MT=MAT(IE) SIG(3)=C*PRP(MT,5) 18 CONTINUE IF(SIG(4).LT.0.D0) SIG(4)=0.D0 DO 4 I=1,4 4 STRSG(IE,I)=-SIG(I) 10 CONTINUE RETURN END C * * * * * * * * * * FUNCTION FUNCTS(JSTEP) C *** HEAVISIDE AND HARMONIC TIME FUNCTION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT C FUNCTS=0.D0 IF(IFUNC.EQ.0) RETURN IF(JSTEP.EQ.0.OR.DFLOAT(JSTEP)*DTIME.GT.DTEND) RETURN IF(IFUNC.EQ.1) FUNCTS=1.D0 IF(IFUNC.EQ.2) ARGUM=OMEGA*JSTEP*DTIME IF(IFUNC.EQ.2) FUNCTS=AZERO+BZERO*DSIN(ARGUM) RETURN END C * * * * * * * * * * * * FUNCTION FUNCTA(ACC,AFACT,IST) C *** ACCELEROGRAM INTERPOLATION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA DIMENSION ACC(6550) C FUNCTA=0.D0 IF(IFUNC.NE.0) RETURN TE=DFLOAT(IST)*DTIME IF(IST.EQ.0.OR.TE.GT.DTEND) RETURN X=(DFLOAT(IST)-1.D0)/AFACT+1.D0 M=X N=M+1 X=X-DFLOAT(M) FUNCTA=ACC(M)*(1.D0-X)+X*ACC(N) RETURN END C * * * * * * * * * * * * SUBROUTINE TRUSS(L) C *** TRUSS ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) DIMENSION LL(4),BM(4),EK(4,4),T(4,4),TEK(4,4),TET(4,4) C II=IJK(L,1) JJ=IJK(L,2) C MT=MAT(L) E=PRP(MT,1) A=PRP(MT,2) XI=XX(II) YI=YY(II) XJ=XX(JJ) YJ=YY(JJ) DX=XJ-XI DY=YJ-YI EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL BM(1)=-1.D0/EL BM(2)=0.D0 BM(3)=1.D0/EL BM(4)=0.D0 DO 10 I=1,4 DO 10 J=1,4 10 EK(I,J)=BM(I)*BM(J)*E*EL*A C DO 20 I=1,4 DO 20 J=1,4 20 T(I,J)=0.D0 T(1,1)=C T(1,2)=S T(2,1)=-S T(2,2)=C T(3,3)=C T(3,4)=S T(4,3)=-S T(4,4)=C CALL WULT(T,EK,TEK,4,4,4) CALL MULT(TEK,T,TET,4,4,4) DO 22 I=1,4 DO 22 J=1,4 22 EKK(L,I,J)=TET(I,J) AES(L)=A*EL DO 24 J=1,4 BMX(L,1,J)=0.D0 DO 24 K=1,4 24 BMX(L,1,J)=BMX(L,1,J)+BM(K)*T(K,J) LL(4)=3*JJ-1 LL(3)=LL(4)-1 LL(2)=3*II-1 LL(1)=LL(2)-1 DO 30 I=1,4 IT=LL(I) 30 LLL(L,I)=LY(IT) RETURN END C * * * * * * * * * * * * SUBROUTINE BEAM(L) C *** BEAM ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) DIMENSION LL(10),EK(6,6),T(6,6),TEK(6,6),TET(6,6),EKT(6,6) C II=IJK(L,1) JJ=IJK(L,2) C MT=MAT(L) E= PRP(MT,1) A= PRP(MT,2) AI=PRP(MT,3) XI=XX(II) YI=YY(II) XJ=XX(JJ) YJ=YY(JJ) DX=XJ-XI DY=YJ-YI EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL G=E*A/EL G5=2.D0*E*AI/EL G4=2.D0*G5 G3=3.D0*G5/EL G2=2.D0*G3/EL C DO 20 I=1,6 DO 20 J=1,6 T(I,J)=0.D0 EK(I,J)=0.D0 20 CONTINUE EK(1,1)=G EK(2,2)=G2 EK(3,3)=G4 EK(4,4)=G EK(5,5)=G2 EK(6,6)=G4 EK(1,4)=-G EK(2,3)=G3 EK(2,5)=-G2 EK(2,6)=G3 EK(3,5)=-G3 EK(3,6)=G5 EK(5,6)=-G3 DO 22 I=1,5 IP1=I+1 DO 24 J=IP1,6 24 EK(J,I)=EK(I,J) 22 CONTINUE DO 26 K=1,4,3 T(K,K)=C T(K,K+1)=S T(K+1,K)=-S T(K+1,K+1)=C T(K+2,K+2)=1.D0 26 CONTINUE CALL WULT(T,EK,TEK,6,6,6) CALL MULT(TEK,T,TET,6,6,6) CALL MULT(EK,T,EKT,6,6,6) DO 28 I=1,6 DO 28 J=1,6 DBM(L,I,J)=EKT(I,J) EKK(L,I,J)=TET(I,J) 28 CONTINUE AES(L)=A*EL LL(6)=3*JJ LL(5)=LL(6)-1 LL(4)=LL(5)-1 LL(3)=3*II LL(2)=LL(3)-1 LL(1)=LL(2)-1 DO 30 I=1,6 IT=LL(I) 30 LLL(L,I)=LY(IT) RETURN END C * * * * * * * * * * * * SUBROUTINE PLANES(L,KOL2) C *** PLANE STRAIN & INTERFACE ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB5/DMX(2150,4,4),BMX(2150,4,8),DBM(2150,6,8) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) DIMENSION LL(10),D(4,4),ES(8,8),IOE(6),B(4,6),DB(4,6), * EK(6,6),BM(4,8),BMT(8),T(8,8),TEK(8,8),XXL(4),YYL(4), * COR(8),COD(8) C II=IJK(L,1) JJ=IJK(L,2) KK=IJK(L,3) MM=IJK(L,4) MT=MAT(L) C DO 2 I=1,4 DO 2 J=1,4 2 D(I,J)=0.D0 C E=PRP(MT,1) P=PRP(MT,2) H=PRP(MT,3) GO TO (80,80,80,14,14,16),KOL2 14 CONTINUE GP=1.D0+P GM=1.D0-P GN=1.D0-2.D0*P G=E/(GP*GN) G1=GM*G G2=P*G G3=0.5D0*E/GP D(1,1)=G1 D(2,2)=G1 D(3,3)=G3 D(4,4)=G1 D(1,2)=G2 D(2,1)=G2 D(1,4)=G2 D(2,4)=G2 D(4,1)=G2 D(4,2)=G2 GO TO 18 16 CONTINUE GP=1.D0+H GM=1.D0-H GN=1.D0-2.D0*H G=P/(GP*GN) G1=GM*G G2=H*G D(1,1)=G1 D(2,2)=G1 D(3,3)=E D(4,4)=G1 D(1,2)=G2 D(2,1)=G2 D(1,4)=G2 D(2,4)=G2 D(4,1)=G2 D(4,2)=G2 DX=XX(JJ)-XX(II) DY=YY(JJ)-YY(II) EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL DO 20 I=1,8 DO 20 J=1,8 20 T(I,J)=0.D0 DO 22 I=1,7,2 T(I,I)=C T(I,I+1)=S T(I+1,I)=-S T(I+1,I+1)=C 22 CONTINUE COR(1)=XX(II) COR(2)=YY(II) COR(3)=XX(JJ) COR(4)=YY(JJ) COR(5)=XX(KK) COR(6)=YY(KK) COR(7)=XX(MM) COR(8)=YY(MM) DO 24 I=1,8 COD(I)=0.D0 DO 24 J=1,8 24 COD(I)=COD(I)+T(I,J)*COR(J) DO 26 I=1,4 IE=2*I IO=IE-1 XXL(I)=COD(IO) YYL(I)=COD(IE) 26 CONTINUE C 18 CONTINUE DO 30 I=1,4 DO 30 J=1,6 30 B(I,J)=0.D0 DO 32 I=1,8 DO 32 J=1,8 32 ES(I,J)=0.D0 DO 34 J=1,8 34 BMT(J)=0.D0 AE=0.D0 I1=0 I2=1 I3=2 N4=4 NV=8 IF(MM.EQ.0) N4=1 IF(MM.EQ.0) NV=6 DO 40 LI=1,N4 I1=I1+1 I2=I2+1 IF(I2.EQ.5) I2=1 I3=I3+1 IF(I3.EQ.5) I3=1 I1E=2*I1 I1O=I1E-1 I2E=2*I2 I2O=I2E-1 I3E=2*I3 I3O=I3E-1 IOE(1)=I1O IOE(2)=I1E IOE(3)=I2O IOE(4)=I2E IOE(5)=I3O IOE(6)=I3E GO TO (80,80,80,44,44,46),KOL2 44 CONTINUE IIT=IJK(L,I1) JJT=IJK(L,I2) KKT=IJK(L,I3) XI=XX(IIT) YI=YY(IIT) XJ=XX(JJT) YJ=YY(JJT) XK=XX(KKT) YK=YY(KKT) GO TO 48 46 CONTINUE XI=XXL(I1) YI=YYL(I1) XJ=XXL(I2) YJ=YYL(I2) XK=XXL(I3) YK=YYL(I3) 48 CONTINUE AI=XJ*YK-XK*YJ AJ=XK*YI-XI*YK AK=XI*YJ-XJ*YI BI=YJ-YK BJ=YK-YI BK=YI-YJ CI=XK-XJ CJ=XI-XK CK=XJ-XI DA=(AI+AJ+AK)*0.5D0 AE=AE+DABS(DA)*0.5D0 D2=0.125D0/DA IF(MM.EQ.0) D2=0.5D0/DA BMT(I1E)=BMT(I1E)+BI*D2 BMT(I1O)=BMT(I1O)+CI*D2 BMT(I2E)=BMT(I2E)+BJ*D2 BMT(I2O)=BMT(I2O)+CJ*D2 BMT(I3E)=BMT(I3E)+BK*D2 BMT(I3O)=BMT(I3O)+CK*D2 DEL2=1.D0/(AI+AJ+AK) B(3,1)=CI*DEL2 B(3,2)=BI*DEL2 B(3,3)=CJ*DEL2 B(3,4)=BJ*DEL2 B(3,5)=CK*DEL2 B(3,6)=BK*DEL2 DO 50 I=1,3 IE=2*I IO=IE-1 B(1,IO)=B(3,IE) B(2,IE)=B(3,IO) 50 CONTINUE CALL MULT(D,B,DB,4,4,6) CALL WULT(B,DB,EK,6,4,6) DO 52 I=1,6 IE=IOE(I) DO 54 J=1,6 JE=IOE(J) ES(IE,JE)=ES(IE,JE)+EK(I,J) 54 CONTINUE 52 CONTINUE 40 CONTINUE C IF(MM.EQ.0) AE=DABS(DA) AES(L)=AE GO TO (80,80,80,64,64,66),KOL2 64 CONTINUE C1=0.25D0 IF(MM.EQ.0) C1=1.D0 DO 70 I=1,NV DO 70 J=1,NV 70 EKK(L,I,J)=ES(I,J)*AE*C1 DO 72 I=1,4 DO 72 J=1,NV 72 BMX(L,I,J)=0.D0 DO 74 J=1,NV 74 BMX(L,3,J)=BMT(J) N1=NV/2 DO 76 J=1,N1 JE=2*J JO=JE-1 BMX(L,1,JO)=BMX(L,3,JE) 76 BMX(L,2,JE)=BMX(L,3,JO) GO TO 68 66 CONTINUE CALL WULT(T,ES,TEK,8,8,8) CALL MULT(TEK,T,ES,8,8,8) DO 78 I=1,8 DO 78 J=1,8 78 EKK(L,I,J)=ES(I,J)*AES(L)*0.25D0 DO 90 I=1,4 DO 90 J=1,8 90 BM(I,J)=0.D0 DO 92 J=1,8 92 BM(3,J)=BMT(J) DO 94 J=1,4 JE=2*J JO=JE-1 BM(1,JO)=BM(3,JE) BM(2,JE)=BM(3,JO) 94 CONTINUE DO 96 I=1,4 DO 96 J=1,8 BMX(L,I,J)=0.D0 DO 96 K=1,8 96 BMX(L,I,J)=BMX(L,I,J)+BM(I,K)*T(K,J) 68 CONTINUE LL(8)=3*MM-1 LL(7)=LL(8)-1 LL(6)=3*KK-1 LL(5)=LL(6)-1 LL(4)=3*JJ-1 LL(3)=LL(4)-1 LL(2)=3*II-1 LL(1)=LL(2)-1 DO 98 I=1,NV IT=LL(I) 98 LLL(L,I)=LY(IT) DO 100 I=1,4 DO 100 J=1,NV DBM(L,I,J)=0.D0 DO 100 K=1,4 100 DBM(L,I,J)=DBM(L,I,J)+D(I,K)*BMX(L,K,J) DO 102 I=1,4 DO 102 J=1,4 102 DMX(L,I,J)=D(I,J) 80 RETURN END C * * * * * * * * * * * * SUBROUTINE MATINV(AA,N0,N1,N2,DET,IND) C *** CALCULATES INVERSE MATRIX IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(N0,N0),IPERM(2220),X(2220) N=N1 M=N+N2 NMX=N0+1 IF(0.GE.N.OR.N2.LT.0.OR.N.GE.NMX.OR.M.GE.NMX) GO TO 80 IND=1 DO 1 I=1,N 1 IPERM(I)=I EPS=0.D0 DO 2 K=1,N RMAX=0.D0 DO 3 J=K,N V=DABS(AA(K,J)) IF(RMAX-V) 4,3,3 4 RMAX=V L=J 3 CONTINUE IF(EPS-RMAX) 5,6,6 6 IF(EPS*0.01D0-RMAX) 7,8,8 8 DET=0.D0 IND=3 DO 9 I=1,N DO 9 J=1,N 9 AA(I,J)=1.0038D0 WRITE(6,200) GO TO 11 7 IND=2 5 PIVOT=AA(K,L) PIVI=1.D0/PIVOT IF(L-K) 12,13,12 12 IW=IPERM(K) IPERM(K)=IPERM(L) IPERM(L)=IW DO 14 I=1,N W=AA(I,K) AA(I,K)=AA(I,L) AA(I,L)=W 14 CONTINUE 13 CONTINUE DO 15 J=1,M X(J)=AA(K,J)*PIVI AA(K,J)=X(J) 15 CONTINUE DO 16 I=1,N IF(I-K) 17,16,17 17 W=AA(I,K) IF(W) 18,16,18 18 DO 19 J=1,M IF(J-K) 20,19,20 20 AA(I,J)=-W*X(J)+AA(I,J) 19 CONTINUE AA(I,K)=-W*PIVI 16 CONTINUE AA(K,K)=PIVI EPS=DMAX1(RMAX*1.E-33,EPS) 2 CONTINUE DO 21 I=1,N 22 K=IPERM(I) IF(K-I) 23,21,23 23 IW=IPERM(K) IPERM(K)=IPERM(I) IPERM(I)=IW DO 24 J=1,M W=AA(I,J) AA(I,J)=AA(K,J) AA(K,J)=W 24 CONTINUE GO TO 22 21 CONTINUE 11 RETURN 80 CONTINUE WRITE(6,201) N,N2 IND=4 GO TO 11 201 FORMAT(/'N1=',I5,' N2=',I5,' MEMORY-OVER') 200 FORMAT(/'THE GIVEN MATRIX IS SINGULAR') END C * * * * * * * * * * * * FUNCTION VFUNC(C) VFUNC=0.D0 IF(C.GT.0.D0) VFUNC=C RETURN END C * * * * * * * * * * * * SUBROUTINE PRINCE(SX,SY,TA,S1,S3,T) C *** CALCULATES PRINCIPAL STRESSES IMPLICIT REAL*8(A-H,O-Z) C=0.5D0*(SX+SY) A1=(SY-SX)*0.5D0 A2=A1*A1+TA*TA A=0.D0 IF(A2.GT.0.D0) A=DSQRT(A2) S1=C+A S3=C-A IF(SY.EQ.S3) GO TO 1 T1=TA/(SY-S3) T=DATAN(T1) GO TO 2 1 T=3.14159D0*0.5D0 2 RETURN END C * * * * * * * * * * * * SUBROUTINE MULT(A,B,C,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(L,M),B(M,N),C(L,N) DO 1 I=1,L DO 1 J=1,N C(I,J)=0.D0 DO 1 K=1,M 1 C(I,J)=C(I,J)+A(I,K)*B(K,J) RETURN END C * * * * * * * * * * * * SUBROUTINE WULT(A,B,C,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M,L),B(M,N),C(L,N) DO 1 I=1,L DO 1 J=1,N C(I,J)=0.D0 DO 1 K=1,M 1 C(I,J)=C(I,J)+A(K,I)*B(K,J) RETURN END C * * * * * * * * * * * * SUBROUTINE VULT(A,U,V,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M,N),U(N),V(M) DO 1 I=1,M V(I)=0.D0 DO 1 J=1,N 1 V(I)=V(I)+A(I,J)*U(J) RETURN END C * * * * * * * * * * SUBROUTINE OUTDYN(IITER,ISTEP) C *** OUTPUT ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB4/MAT(2150),IST(2150),IX(1998),IY(1998),IQ(1998) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA11/IREQD(10),IREQS(10),NOUTD(6550) COMMON/LA12/ACCEH(6550),ACCEV(6550) COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(2220),DISPL(2220),DISPT(2220) COMMON/LA17/VELOI(2220),VELOL(2220),VELOT(2220) COMMON/LA18/ACCEI(2220),ACCEL(2220) COMMON/LA20/IYIEL(2150),IUNLO(2150) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) DIMENSION STRSP(2150,8),UX(1998),UY(1998),UXH(20), * VX(1998),AX(1998),VXH(20),AXH(20),AXA(20),TT(2150) DIMENSION IYITE(2150) C KSTEP=ISTEP IOUT=NOUTD(KSTEP) XTIME=DFLOAT(KSTEP)*DTIME WRITE(6,200) KSTEP,XTIME C *** REARRANGE DISPLACEMENT VECTOR DO 10 IP=1,NPOIN I1=(IP-1)*3+1 IT=LY(I1) UX(IP)=DISPI(IT) VX(IP)=VELOI(IT) AX(IP)=ACCEI(IT) I2=(IP-1)*3+2 IT=LY(I2) UY(IP)=DISPI(IT) 10 CONTINUE IF(IOUT.LE.0) GO TO 11 WRITE(6,201) DO 12 IP=1,NPOIN 12 WRITE(6,202) IP,UX(IP),UY(IP) 11 CONTINUE DO 52 ID=1,NREQD IP=IREQD(ID) WRITE(6,202) IP,UX(IP),UY(IP) UXH(ID)=UX(IP) VXH(ID)=VX(IP) AXH(ID)=AX(IP) AXA(ID)=AXH(ID)+ACCEH(KSTEP) 52 CONTINUE WRITE(8,205) XTIME,(UXH(ID),ID=1,NREQD) WRITE(9,205) XTIME,(VXH(ID),ID=1,NREQD) WRITE(10,205) XTIME,(AXH(ID),AXA(ID),ID=1,NREQD) C *** WRITES STRESSES WRITE(6,203) DO 14 IE=1,NELEM KOL2=K2(IE) MT=MAT(IE) TT(IE)=0.D0 GO TO (21,21,21,24,24,21),KOL2 21 CONTINUE DO 28 I=1,6 28 STRSP(IE,I)=STRSG(IE,I) STRSP(IE,7)=0.D0 A=1.D0 IF(KOL2.EQ.1.OR.KOL2.EQ.3) A=PRP(MT,2) TT(IE)=STRSG(IE,1)*A GO TO 20 24 CONTINUE SX=-STRSG(IE,1) SY=-STRSG(IE,2) TA=-STRSG(IE,3) SZ=-STRSG(IE,4) CALL PRINCE(SX,SY,TA,S1,S3,AG) STRSP(IE,1)=SX STRSP(IE,2)=SY STRSP(IE,3)=TA STRSP(IE,4)=SZ STRSP(IE,5)=S1 STRSP(IE,6)=S3 STRSP(IE,7)=PI1(IE) STRSP(IE,8)=PI2(IE) 20 CONTINUE IF(IOUT.LE.0) GO TO 14 WRITE(6,204) IE,IYIEL(IE),(STRSP(IE,I),I=1,8) 14 CONTINUE DO 30 ID=1,NREQS IE=IREQS(ID) 30 WRITE(6,204) IE,IYIEL(IE),(STRSP(IE,I),I=1,8) DO 32 IE=1,NELEM IYITE(IE)=0 IF(JTENS(IE).GE.1) IYITE(IE)=2 IF(IYIEL(IE).GE.1) IYITE(IE)=1 32 CONTINUE WRITE(6,207) (IYITE(IE),IE=1,NELEM) C *** OUTPUT FOR INITIAL STRESSES OPEN(7,FILE='DASTRE0',STATUS='UNKNOWN') DO 34 I=1,4 34 WRITE(7,206) (STRSP(IE,I),IE=1,NELEM) CLOSE(7) 80 CONTINUE C 200 FORMAT(/'*** DISPLACEMENTS AT TIME STEP',I5,3X,'TIME=',E12.4) 201 FORMAT(/1X,'NODE',6X,'X-DISP',6X,'Y-DISP') 202 FORMAT(I5,3X,E10.3,3X,E10.3) 203 FORMAT(/'STRESSES'/2X,'EL YIELD',5X,'X-ST',5X,'Y-ST', * 5X,'XY-ST',5X,'Z-ST',5X,'MAX-ST',5X,'MIN-ST',5X,'F-1 F-2') 204 FORMAT(I5,I3,2X,8E10.3) 205 FORMAT(F12.4,6E12.4) 206 FORMAT(10E10.3) 207 FORMAT(24I2) 221 FORMAT(I5,3E12.4) RETURN END C * * * * * * * * * * * * SUBROUTINE MAVS2 C *** MICRO-AVS10 OUTPUT (PART 2: DISPLACEMENTS & YIELD ELEMENTS) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 AC(900) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA20/IYIEL(2150),IUNLO(2150) COMMON/LA21/NTENS,KTENS,ITENS(2150),JTENS(2150) DIMENSION UX(1280),UY(1280),UM(1280),WST(1260,6),PI9(1260) C DO 2 IP=1,NPOIN I1=(IP-1)*3+1 IT=LY(I1) UX(IP)=DISPI(IT) I2=(IP-1)*3+2 IT=LY(I2) UY(IP)=DISPI(IT) 2 CONTINUE WRITE(12,310) 1 WRITE(12,311) WRITE(12,300) 1 C WRITE(12,301) NPOIN,NELEM DO 10 I=1,NPOIN 10 WRITE(12,302) I,XX(I),YY(I),0 DO 12 L=1,NELEM AC(L)='quad' IF(IJK(L,4).EQ.0) AC(L)='tri' IF(IJK(L,3).EQ.0) AC(L)='line' WRITE(12,303) L,K2(L),AC(L),(IJK(L,I),I=1,4) 12 CONTINUE WRITE(12,304) WRITE(12,305) WRITE(12,306) WRITE(12,307) WRITE(12,308) DO 14 I=1,NPOIN 14 WRITE(12,309) I,UX(I),UY(I),0 WRITE(12,313) WRITE(12,314) DO 20 L=1,NELEM PI9(L)=0.1D0 IF(JTENS(L).GE.1) PI9(L)=1.D0 IF(IYIEL(L).GE.1) PI9(L)=2.D0 20 CONTINUE DO 22 L=1,NELEM 22 WRITE(12,315) L,PI9(L) C 300 FORMAT('step',I1) 301 FORMAT(2I5) 302 FORMAT(I5,3F10.4) 303 FORMAT(2I5,A6,20I5) 304 FORMAT(4X,'3',4X,'1') 305 FORMAT(4X,'3',4X,'1',4X,'1',4X,'1') 306 FORMAT(1X,'disp_x, m') 307 FORMAT(1X,'disp_y, m') 308 FORMAT(1X,'disp_z, m') 309 FORMAT(I5,3E14.4) 310 FORMAT(I5) 311 FORMAT('data_geom') 313 FORMAT(4X,'1',4X,'1') 314 FORMAT(1X,'no') 315 FORMAT(I5,E12.3) RETURN END C * * * * * * * * * * * * SUBROUTINE MAVS3 C *** MICRO-AVS10 OUTPUT (PART 3: PRINCIPAL STRESSES) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 AC(900) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,NSIZE,NN COMMON/LAB2/LY(6100),LLL(2150,8),IJK(2150,4),K2(2150) COMMON/LAB3/EKK(2150,8,8),AES(2150),PW(2150,4),PRP(20,12) COMMON/LAB6/XX(1998),YY(1998),PI1(2150),PI2(2150) COMMON/LAB7/NSTEP,NOUTP,NREQD,NREQS,NACCE,IFUNC COMMON/LA15/STRIN(2150,6),STRAG(2150,6),STRSG(2150,6) COMMON/LA20/IYIEL(2150),IUNLO(2150) DIMENSION WST(6),AX(1260),AY(1260),XXX(1280),YYY(1280) DIMENSION J(1690) C DO 2 L=1,NELEM J(L)=0 AX(L)=0.D0 AY(L)=0.D0 2 CONTINUE C DO 10 L=1,NELEM KOL2=K2(L) IF(KOL2.LE.3) GO TO 10 DO 12 I=1,3 12 WST(I)=STRSG(L,I) SX=WST(1) SY=WST(2) TA=WST(3) CALL PRINCE(SX,SY,TA,S1,S3,TH) AX(L)=S1*DSIN(TH) AY(L)=S1*DCOS(TH) 10 CONTINUE C NNE1=0 DO 6 L=1,NELEM IF(IJK(L,3).EQ.0) GO TO 6 NNE1=NNE1+1 6 CONTINUE WRITE(13,301) NPOIN,NNE1 C DO 30 I=1,NPOIN XXX(I)=XX(I)*100.D0 YYY(I)=YY(I)*100.D0 WRITE(13,302) I,XXX(I),YYY(I),0 30 CONTINUE DO 32 L=1,NELEM AC(L)='quad' IF(IJK(L,4).EQ.0) AC(L)='tri' IF(IJK(L,3).EQ.0) GO TO 32 WRITE(13,303) L,J(L),AC(L),(IJK(L,I),I=1,4) 32 CONTINUE WRITE(13,304) WRITE(13,305) WRITE(13,306) DO 34 L=1,NELEM 34 WRITE(13,307) L,AX(L),AY(L) C 301 FORMAT(2I5,4X,'0',4X,'2',4X,'0') 302 FORMAT(I5,2F12.3,F9.3) 303 FORMAT(2I5,A8,I4,20I5) 304 FORMAT(3X,'2',3X,'1',4X,'1') 305 FORMAT('strX , (kN/m2)') 306 FORMAT('strY , (kN/m2)') 307 FORMAT(I5,3E13.4) RETURN END