C SIMULATION PROGRAM FOR THREE-DIMENSIONAL MECHANICAL SYSTEMS C *********************************************************** C C ----STAGE II MAINLINE---- C INTEGER*4 STMASS,STNVDR,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,NAME COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STNVDR,STDAMP,STSPRG,STFORC, 2 KOUNT C<<<<< NETWORK SIZE SPECIFICATION FOLLOWS>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DIMENSION M(7,10),V(10,3),X(10,3),F(10,3),SIZE(10),SPLENG(10) DIMENSION VZERO(10,3),XZERO(10,3),RK(10,3,4),FINALK(10,3) DIMENSION NOD(7),PRIMEK(10,3) C C ----READ IN VARIOUS DATA----- CALL READER(V,X,SPLENG,NCOL,NSPRNG,STEP) C<<<< CALCULATE INDEXES FOR DIFFERENT ELEMENTS NT=NN-1 STMASS=1 STNVDR=NMAS+1 STDAMP=STNVDR+NVDR STSPRG=STDAMP+NDAM STFORC=STSPRG+NSPR ENDMAS=NMAS ENDVDR=ENDMAS+NVDR ENDAMP=ENDVDR+NDAM ENDSPR=ENDAMP+NSPR ENDFOR=ENDSPR+NFOR C NROW=NN NCOL=NB C-----SET START TIME TO ZERO ----- T=0.0 NSPRNG=NSPR NUMASS=NMAS C -----FORM CUTSET MATRIX----- CALL TREE(M,SIZE,NOD,NAME,NROW,NCOL) C CALL CLEAR(.TRUE.) CALL ROTER CALL AXIS C ----TOP OF MAIN SIMULATION LOOP----- H=STEP 100 CONTINUE IF((NVDR+NFOR).EQ.0) GO TO 201 CALL DRIVER(V,X,F,NCOL) 201 CONTINUE CALL VELPOS(V,X,M,NROW,NCOL) IF(NDAM.EQ.0) GO TO 203 CALL DAMPER(V,X,F,SIZE,NROW,NCOL) 203 CONTINUE IF(NSPR.EQ.0) GO TO 204 CALL SPRING(X,F,SIZE,SPLENG,NCOL,NSPRNG) 204 CONTINUE CALL FORCE1(F,M,NROW,NCOL) C C ----CALL STAGE III FOR GRAPHIC DISPLAY----- C CALL GRAFIX(NT,NN,NB,X,F) C C C -----INTEGRATE AND RETURN TO GRAPHICS SECTION---- 888 CONTINUE CALL RUNGE(M,V,X,F,SIZE,SPLENG,VZERO,XZERO,RK,NROW,NCOL,NSPRNG, 1 NUMASS,FINALK,PRIMEK) IF(T.GT.TLIMIT) GO TO 8 GO TO 100 C 8 CONTINUE STOP 1 END C SUBROUTINE READER(V,X,SPLENG,NCOL,NSPRNG,STEP) C INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT COMMON/LIMITS/XMIN,XMAX,YMIN,YMAX COMMON/EYCEP/EYE(3),CEN(3),UP(3) COMMON/DRVDAT/BASE(10,3),AMPL(10,3),OMEGA(10,3), 9PHASE(10,3),DRTIME(10,2) COMMON/MODDAT/IAX,DIST,SMODEL(7,17,2),MODNO(10) DIMENSION V(10,3),X(10,3),SPLENG(10),AN(2),BN(2) DATA AN(1),AN(2)/5HDYNDA,4HTSRC/ DATA BN(1),BN(2)/5HMODEL,4HSSRC/ C C ----READ IN MODEL SYMBOLS---- CALL OPEN(1,'MODELSSRC') DO 15 K=1,7 READ(1,18)MODNO(K) IND=MODNO(K) DO 16 IL=1,IND 16 READ(1,19)(SMODEL(K,IL,IE),IE=1,2) 15 CONTINUE 18 FORMAT(1X,I3) 19 FORMAT(2F10.4) CALL CLOSE(1) C ------READ IN ELEMENT QUANTITY---- READ(1,2)IAX,NN,NB,NMAS,NVDR,NDAM,NSPR,NFOR C ----READ MASS INITIAL COND.------ DO 10 J=1,NMAS 10 READ(1,4)X(J,1),X(J,2),X(J,3),V(J,1),V(J,2),V(J,3) C -----READ DRIVER FUNCTION DATA----- K=NVDR+NFOR IF(K.EQ.0)GO TO 20 DO 11 J=1,K READ(1,4)BASE(J,1),BASE(J,2),BASE(J,3),AMPL(J,1),AMPL(J,2), 9AMPL(J,3) READ(1,4)OMEGA(J,1),OMEGA(J,2),OMEGA(J,3),PHASE(J,1), 9PHASE(J,2),PHASE(J,3) READ(1,4)DRTIME(J,1),DRTIME(J,2) 11 CONTINUE C ------READ UNSTRETCHED SPRING------ 20 IF(NSPR.EQ.0)GO TO 30 DO 12 J=1,NSPR READ(1,4)SPLENG(J) 12 CONTINUE C -----INTEGRATION PARAMETERS----- 30 READ(1,4)STEP,TLIMIT C -----LIMITS----- READ(1,4)XMIN,XMAX,YMIN,YMAX READ(1,4)EYE(1),EYE(2),EYE(3) READ(1,4)CEN(1),CEN(2),CEN(3) C ----CALCULATE UP AND DIST---- UP(1)=CEN(1) UP(2)=CEN(2)+10. UP(3)=CEN(3) DIST=EYE(3)-CEN(3) 2 FORMAT(8I5) 4 FORMAT(6F10.4) C RETURN END C <<<<<<<<<<<<<<<<<<<<<<< SUBROUTINE TREE(M,SIZE,NOD,NAME,NROW,NCOL) C <<<<<<<<<<<<<<<<<<<<<<< C C*****TREE GENERATION SUBROUTINE INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,NAME,NAME1,NAME2 COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT COMMON/GENDAT/MCON(20,2),NTYPE(20),BSCALE(20), 9KOLOR(20),IDENT(20) DIMENSION M(7,10),NOD(7),AN(2),SIZE(10) DATA AN(1),AN(2)/5HDYNDA,4HTSRC/ C*****SET VALUES OF MATRICES AND COUNTERS DO 8 I=1,NN DO 8 J=1,NB M(I,J)=0 8 CONTINUE DO 1 I=1,NN NOD(I)=I 1 CONTINUE IT=0 IC=NN C C*****BEGIN TO INPUT BRANCHES DO 4 I=1,NB READ(1,30)NSTART,NEND,IDENT(I),NTYPE(I),KOLOR(I), 9BRMAG,BSCALE(I) 30 FORMAT(5I5,2F10.4) MCON(I,1)=NSTART MCON(I,2)=NEND NUMERR=1 C C*****SORT INTO TREE/COTREE IF(IT.EQ.NT) GO TO 2 IF(NOD(NSTART).EQ.NOD(NEND)) GO TO 2 C C*****INSERT VALUES IN TREE SIZE(I)=BRMAG M(NSTART,I)=1 M(NEND,I)=2 NSS=NOD(NSTART) NNN=NOD(NEND) DO 3 I2=1,NN IF(NOD(I2).EQ.NSS) NOD(I2)=NNN 3 CONTINUE IT=IT+1 GO TO 4 C C*****INSERT VALUES IN COTREE 2 SIZE(IC)=BRMAG M(NSTART,IC)=1 M(NEND,IC)=2 IC=IC+1 4 CONTINUE C C*****PRINT OUT NAMES OF ELEMENTS IN TREE & COTREE C C******BEGIN GENERATING CUTSET MATRIX DO 11 I=1,NT IF(M(I,I).NE.0.0) GO TO 24 C C*****CHECK COLUMN ENTRIES BELOW DIAGONAL I1=I+1 IF(I1.GT.NT) GO TO 1000 DO 12 K3=I1,NT IF(M(K3,I).GT.0) GO TO 23 12 CONTINUE NUMERR=2 GO TO 1000 C C*****ADD NON-ZERO ROW TO ROW(I) 23 DO 13 K4=1,NB M(I,K4)=M(I,K4)+M(K3,K4) 13 CONTINUE C C*****SET ALL COLUMN VALUES TO ZERO EXCEPT DIAGONAL 24 DO 14 K1=1,NT IF(K1.EQ.I) GO TO 14 21 IF(M(K1,I).EQ.0) GO TO 14 25 DO 15 J=1,NB M(K1,J)=M(K1,J)+M(I,J) IF(M(K1,J).GT.2) M(K1,J)=M(K1,J)-3 15 CONTINUE GO TO 21 14 CONTINUE 11 CONTINUE C C*****REPLACE TWOS BY +1 OR -1 NUMERR=3 26 DO 16 I=1,NT IF(M(I,I).EQ.0) GO TO 1000 IF(M(I,I).EQ.1) GO TO 28 M(I,I)=1 27 DO 17 K=NN,NB IF(M(I,K).EQ.0) GO TO 17 IF(M(I,K).EQ.2) GO TO22 M(I,K)=-1 GO TO 17 22 M(I,K)=1 17 CONTINUE GO TO 16 28 DO 18 K=NN,NB IF(M(I,K).GT.1) M(I,K)=-1 18 CONTINUE 16 CONTINUE C C*****PRINT OUT CUTSET MATRIX RETURN C C*****ERROR MESSAGES ORIGINATE HERE 1000 CONTINUE STOP 3 END C SUBROUTINE SPECAL(V,X,F,M) C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT DIMENSION V(10,3),X(10,3),F(10,3) C C THIS SUBROUTINE CONTAINS THE TERMINAL EQUATIONS FOR THE SPECIFIED C FUNCTIONS.THERE SHOULD BE ONE EQUATION FOR EACH VECTOR COMPONENT O C FUNCTIONS. IF THERE ARE (NVDR) VELOCITY AND POSITION DRIVERS AND C FORCE DRIVERS THEN THERE SHOULD BE (6*NVDR+3*NFOR) EQUATIONS FOLLO RETURN END C SUBROUTINE DRIVER(V,X,F,NCOL) C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT COMMON/DRVDAT/BASE(10,3),AMPL(10,3),OMEGA(10,3), 9PHASE(10,3),DRTIME(10,2) COMMON/GENDAT/MCON(20,2),NTYPE(20),BSCALE(20), 9KOLOR(20),IDENT(20) DIMENSION V(10,3),X(10,3),F(10,3) C* LOOP FOR CALCULATING POSITION AND VELOCITY DRIVERS ISPEC=0 IF(NVDR.EQ.0) GO TO 200 DO 100 J=1,NVDR I=J+NMAS TSTART=DRTIME(J,1) TEND =DRTIME(J,2) EFTIME=T-TSTART ID=IDENT(I) IF(ID.EQ.5)GO TO 50 X(I,1)=BASE(J,1) X(I,2)=BASE(J,2) X(I,3)=BASE(J,3) C ---SET EFFECTIVE TIME AND BRANCH TO TYPE--- GO TO (10,20,30,40),ID C ---STEP DRIVER TYPE---- 20 IF(T.GT.TEND)GO TO 10 IF(EFTIME.LE.0.0)GO TO 10 C DO 21 K=1,3 21 X(I,K)=X(I,K)+AMPL(J,K) GO TO 10 C ----RAMP TYPE---- C 30 IF(EFTIME.LE.0.0)GOTO 10 IF(T.GT.TEND)EFTIME=TEND-TSTART DO 31 K=1,3 31 X(I,K)=X(I,K)+EFTIME*AMPL(J,K) IF(T.GT.TEND)GO TO 10 DO 32 K=1,3 32 V(I,K)=AMPL(J,K) GO TO 100 C ----SINUSOID DRIVER------ 40 ANGX=PHASE(J,1) ANGY=PHASE(J,2) ANGZ=PHASE(J,3) IF(T.GT.TEND)GO TO 41 IF(EFTIME.LE.0.0)GO TO41 ANGX=OMEGA(J,1)*EFTIME+ANGX ANGY=OMEGA(J,2)*EFTIME+ANGY ANGZ=OMEGA(J,3)*EFTIME+ANGZ GO TO 42 41 EFTIME=0.0 42 X(I,1)=X(I,1)+AMPL(J,1)*SIN(ANGX) X(I,2)=X(I,2)+AMPL(J,2)*SIN(ANGY) X(I,3)=X(I,3)+AMPL(J,3)*SIN(ANGZ) C BRANCH---- IF(EFTIME.EQ.0.0)GO TO 10 C ---SET VELOCITY--- V(I,1)=AMPL(J,1)*OMEGA(J,1)*COS(ANGX) V(I,2)=AMPL(J,2)*OMEGA(J,2)*COS(ANGY) V(I,3)=AMPL(J,3)*OMEGA(J,3)*COS(ANGZ) GO TO 100 C ---ST C ------SET VELOCITY TO ZERO---- 10 V(I,1)=0.0 V(I,2)=0.0 V(I,3)=0.0 GO TO 100 C* SPECIAL DRIVER - SPECIAL SUBROUTINE WILL BE CALLED BEFORE RETURN 50 CONTINUE 100 CONTINUE C* C* LOOP FOR CALCULATING FORCE DRIVERS 200 IF(NFOR.EQ.0) GO TO 500 DO 300 J2=1,NFOR J=J2+NVDR I=STFORC+J2-1 TSTART=DRTIME(J,1) TEND =DRTIME(J,2) ID=IDENT(I) EFTIME=T-TSTART C* CONSTANT FORCE - VALUE OF FO SET INITIALLY 210 DO 16 K=1,3 F(I,K)=BASE(J,K) 16 CONTINUE GO TO (300,220,230,240,250),ID C* STEP FORCE INPUT 220 IF(EFTIME.LE.0.0)GO TO 300 IF(T.GT.TEND)GO TO 300 DO 221 K=1,3 221 F(I,K)=BASE(J,K) + AMPL(J,K) GO TO 300 C* RAMP FORCE INPUT 230 IF(EFTIME.LE.0.0)GO TO 300 IF(T.GT.TEND) EFTIME=TEND-TSTART DO 231 K=1,3 231 F(I,K)=BASE(J,K) + EFTIME*AMPL(J,K) GO TO 300 C* SINUSOIDAL FORCE INPUT 240 IF(EFTIME.LE.0.0)GO TO 300 IF(T.GT.TEND) GO TO 300 DO 241 K=1,3 W=OMEGA(J,K) ANG=W*EFTIME + PHASE(J,K) 241 F(I,K)=BASE(J,K) + AMPL(J,K)*SIN(ANG) GO TO 300 C* SPECIAL DRIVER - SPECIAL SUBROUTINE WILL BE CALLED BEFORE RETURN 250 CONTINUE 300 CONTINUE 500 IF(ISPEC.EQ.0) GO TO 1000 CALL SPECAL(V,X,F,M) 1000 RETURN END C ------------------------------- SUBROUTINE VELPOS(V,X,M,NROW,NCOL) INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT DIMENSION V(10,3),X(10,3),M(7,10) ITRANS=1 N1=1 N2=NT DO 100 I=STDAMP,NB IROW=I CALL SYSUM (V,M,IROW,N1,N2,ITRANS,NROW,NCOL) V(I,1)=-SUM(1) V(I,2)=-SUM(2) V(I,3)=-SUM(3) CALL SYSUM (X,M,IROW,N1,N2,ITRANS,NROW,NCOL) X(I,1)=-SUM(1) X(I,2)=-SUM(2) X(I,3)=-SUM(3) 100 CONTINUE RETURN END C SUBROUTINE DAMPER (V,X,F,SIZE,NROW,NCOL) INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT DIMENSION UR(3),F(10,3),V(10,3),X(10,3),SIZE(10) DO 700 I=STDAMP,ENDAMP C<<< CALCULATE UNIT VECTOR RMAG=X(I,1)*X(I,1)+X(I,2)*X(I,2)+X(I,3)*X(I,3) R=SQRT(RMAG) IF(R.EQ.0.0) GO TO 702 UR(1)=X(I,1)/R UR(2)=X(I,2)/R UR(3)=X(I,3)/R C<<<< CALCULATE DAMPER FORCE D=SIZE(I)*(V(I,1)*UR(1)+V(I,2)*UR(2)+V(I,3)*UR(3)) D=-D F(I,1)=D*UR(1) F(I,2)=D*UR(2) F(I,3)=D*UR(3) GO TO 700 702 CONTINUE C<<< PROCEDURE FOR ZERO DAMPER LENGTH F(I,1)=SIZE(I)*V(I,1) F(I,2)=SIZE(I)*V(I,2) F(I,3)=SIZE(I)*V(I,3) 700 CONTINUE RETURN END C SUBROUTINE SPRING(X,F,SIZE,SPLENG,NCOL,NSPRNG) INTEGER*4 STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR COMMON T,TLIMIT,H,SUM(3) COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS, 1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC, 2 KOUNT DIMENSION X(10,3),F(10,3),SIZE(10),SPLENG(10),UR(3) N=0 DO 800 I=STSPRG,ENDSPR N=N+1 C<<