PROGRAM lfemi4lpcmech C C23456789012345678901234567890123456789012345678901234567890123456789012 C 1 2 3 4 5 6 7 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * M A I N P R O G R A M * C * NAME: ( LFEMI ) * C * L INEAR F INITE E LEMENT M ODEL I N-CORE * C * * C * CORE STORAGE MINIMIZED * C * 1. STORES ELEMENT-NODE DATA ON SCRATCH FILES (TAPE 1-4,7-9). * C * 2. DIMENSIONED VARIABLES SHARE HIGH SPEED CORE STORAGE. * C * 3. USES THE SKYLINE STORAGE SCHEME, DEVELOPED BY BATHE AND * C * WILSON ( NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS), TO * C * STORE THE STIFFNESS ARRAY AS A SINGLE DIMENSIONED ARRAY * C * SUCH THAT A COLUMN ORIENTED GAUSS ELIMINATION TECHNIQUE IS * C * USED IN SUBROUTINE COLSOL TO SOLVE THE SET OF LINEAR EQNS. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL NU12,NU13,NU23 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E11,E22,E33,G12,G13,G23,NU12,NU13,NU23 COMMON/RDK3/AL1,AL2,AL3,BT1,BT2,BT3 C DIMENSION TITLE(20,5) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * THE FOLLOWING TWO CARDS DETERMINE THE MAXIMUM HIGH SPEED CORE * C * STORAGE SHARED BY DIMENSIONED VARIABLES IN VARIOUS SUBROUTINES. * C * TO CHANGE STORAGE, CHANGE BOTH MTOT AND SHARE(MTOT). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON SHARE(119197) MTOT=119197 C C DEFINE READ AND WRITE UNIT NUMBERS C LR=5 LW=6 C C ASSIGN TAPE NO.S 5 AND 6 TO INPUT AND OUTPUT FILES RESPECTIVELY C OPEN(5,FILE='lfemi4lpc_mech.DAT',STATUS='OLD',ERR=888) OPEN(6,FILE='lfemi4lpc_mech.out',STATUS='UNKNOWN',ERR=888) OPEN(10,FILE='lfemi4lpc_mechstress.out',STATUS='UNKNOWN',ERR=888) C C H E A D I N G C WRITE(LW,999) 999 FORMAT(/,80(1H*),/) READ(LR,10)((TITLE(I,J),I=1,20),J=1,5) 10 FORMAT(20A4) WRITE(LW,11)((TITLE(I,J),I=1,20),J=1,5) 11 FORMAT(1X,20A4) WRITE(LW,999) C C G R I D P A R A M E T E R S C READ(LR,20)NE,NDS,NLOAD,NDISP,KEY 20 FORMAT(5I6) C WRITE(LW,22)NE,NDS,NLOAD,NDISP,KEY 22 FORMAT(//,5X,'G R I D P A R A M E T E R S',/, 1' NUMBER OF ELEMENTS .........................',I4,/, 2' NUMBER OF NODES ............................',I4,/, 3' NUMBER OF NODE LOADS .......................',I4,/, 4' NO. OF PRESCRIBED NODE DISPLACEMENTS .......',I4,/, 5' OUTPUT FILE PRINT KEY ......................',I4) C C M A T E R I A L P R O P E R I T I E S C WRITE(LW,30) 30 FORMAT(//,10X,'M A T E R I A L P R O P E R T I E S',/) C READ(LR,40)E11,E22,E33 40 FORMAT(3(6X,E10.3)) WRITE(LW,41)E11,E22,E33 41 FORMAT(/,15X,'YOUNGS MODULI',/,' E11=',E10.3,' E22=',E10.3, *' E33=',E10.3) READ(LR,40)G12,G13,G23 WRITE(LW,42)G12,G13,G23 42 FORMAT(/,15X,'SHEAR MODULI',/,' G12=',E10.3,' G13=',E10.3, *' G23=',E10.3) READ(LR,40)NU12,NU13,NU23 WRITE(LW,43)NU12,NU13,NU23 43 FORMAT(/,15X,'POISSONS RATIO',/,' NU12=',E10.3,' NU13=',E10.3, *' NU23=',E10.3,//) C READ(LR,50)AL1,AL2,AL3 READ(LR,50)BT1,BT2,BT3 50 FORMAT(3(8X,E10.3)) C WRITE(LW,51)AL1,AL2,AL3 51 FORMAT(/,15X,'THERMAL EXPANSION COEFFICIENTS',/, 1' ALPHA-1=',E10.3,' ALPHA-2=',E10.3,' ALPHA-3=',E10.3,/) WRITE(LW,52)BT1,BT2,BT3 52 FORMAT(/,15X,'MOISTURE EXPANSION COEFFECIENTS',/, 2' BETA-1=',E10.3,' BETA-2=',E10.3,' BETA-3=',E10.3,//) C C L O A D I N F O R M A T I O N C O T H E R T H E N L O A D S A T N O D E P O I N T S C WRITE(LW,60) 60 FORMAT(1H1,10X,'L O A D I N F O R M A T I O N',/,2X, *'O T H E R T H E N L O A D S A T N O D E P O I N T S', *//) C READ(LR,61)LTYPE 61 FORMAT(I6) READ(LR,62)SLOAD 62 FORMAT(F12.0) READ(LR,63)TLOAD,WLOAD 63 FORMAT(2F12.0) C WRITE(LW,64)LTYPE 64 FORMAT(/,10X,'LOAD TYPE = ',I1,' ( = 0, NODE LOADS ONLY)', */,23X,' ( = 1, STRAIN APPLIED NORMAL TO ELEMENTS ONLY)', */,23X,' ( = 2, CHANGE IN THERMAL-MOISTURE CONDITIONS ONLY)',//) WRITE(LW,65)LTYPE 65 FORMAT(5X,'LOAD TYPE =',I2,///) IF(LTYPE.LE.1)WRITE(LW,66)SLOAD 66 FORMAT(' STRAIN APPLIED NORMAL TO H-Z PLANE OF ELEMENTS',/, 1' SLOAD=',E11.4,//) IF(LTYPE.EQ.2)WRITE(LW,67)TLOAD,WLOAD 67 FORMAT(' STRAINS NORMAL TO PLANE OF ELEMENTS',/, 4' TO BE CALCULATED',/, 5' FOR NONZERO TEMPERATURE CHANGE, TLOAD=',E11.4,/, 6' MOISTURE CHANGE, WLOAD=',E11.4,/) C C N O D E - E L E M E N T I N F O R M A T I O N C C ADDRESSES N1=1 N2=N1+3*NDS N3=N2+NDS C C NODE BOUNDARY CONDITION CODES ID(NDS,3) C NODE HORIZONTAL CO-ORDINATES H(NDS) C NODE VERTICAL CO-ORDINATES Z(NDS) C C ID(NDS,3) H(NDS) Z(NDS) CALL NEDATA(SHARE(N1),SHARE(N2),SHARE(N3),NDS,NE) C C C A L C U L A T E M A T E R I A L S T I F F N E S S E S C F O R E A C H E L E M E N T C CALL CIJ(NE) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * C * F. E. M. S O L U T I O N * C * * * * * * * * * * * * * * * * * * * * * * * * * * * C C A S S I G N E Q U A T I O N N U M B E R S TO C E A C H D E G R E E O F F R E E D O M C STORED IN ARRAY ID(NDS,3) C C ID(NDS,3) CALL EQNUM(SHARE(N1),NDS,NEQ,NEQF,LTYPE) C C C A L C U L A T E A D D R E S S E S U S E D T O S T O R E C T H E G L O B A L S T I F F N E S S A R R A Y C C GLOBAL ARRAY SIZE, NCOFS (NO. OF COEFFICIENTS) C GLOBAL ADDRESSES OF DIAGONAL TERMS, MAXA(NEQ1) C GLOBAL HALF BANDWIDTH INCLUDING DIAGONAL, MBAND C ELEMENT ARRAY SIZE, NDOFE (DEGREES OF FREEDOM PER ELEMENT) C THERMAL-MOISTURE LOAD EQN. NO., NEQF C C CALCULATE NO. EQN.S PLUS ONE, NEQ1 USED TO C STORE DIAGONAL TERMS OF GLOBAL STIFFNESS ARRAY NEQ1=NEQ+1 C C ID(NDS,3) MAXA(NEQ1) CALL ADDRS(SHARE(N1),SHARE(N2),NE,NDS,NEQ1,LTYPE,NEQF,MBAND,NCOFS, * NDOFE) C C C A L C U L A T E T O T A L S I Z E O F A R R A Y S C T O B E S T O R E D I N C O M M O N C C ARRAYS USED IN SUBROUTINE ELMNT AND SHARED BY MAIN PROG. COMMON C C ARRAY NAMES COMMON ADDRESSES C ID( NDS,3 ) = SHARE(N1) N1=1 C MAXA( NEQ1 ) = SHARE(N2) N2=N1+3*NDS C A( NCOFS ) = SHARE(N3) N3=N2+NEQ1 C V( NEQ ) = SHARE(N4) N4=N3+NCOFS C NTOT1=3*NDS+NEQ1+NCOFS+NEQ C C ARRAYS USED IN SUBROUTINE STR AND SHARED BY MAIN PROG. COMMON C C ARRAY NAMES COMMON ADDRESSES C ID( NDS,3 ) = SHARE(N1) N1=1 SPACE ALLOCATED FOR ID ARRAY C PRIOR TO SUBPROG. STR C U( NDS,3 ) = SHARE(N2) N2=N1+NDS*3 C EN( NE ) = SHARE(N3) N3=N2+NDS*3 C EH( NE ) = SHARE(N4) N4=N3+NE C EZ( NE ) = SHARE(N5) N5=N4+NE C EHZ( NE ) = SHARE(N6) N6=N5+NE C ENZ( NE ) = SHARE(N7) N7=N6+NE C ENH( NE ) = SHARE(N8) N8=N7+NE C NTOT2=6*NDS+6*NE C C CHOOSE LARGEST SHARED ARRAY SIZE C NTOT=NTOT1 IF(NTOT2.GT.NTOT1)NTOT=NTOT2 C WRITE(LW,999) WRITE(LW,70)MTOT,NTOT1,NTOT2 70 FORMAT(' TOTAL SIZE OF ARRAYS STORED IN COMMON',/, 1' AVAILBLE(',I6,') REQUIRED BY SUBROUTINE ELMNT (',I6,')',/, 233X,'SUBROUTINE STR (',I6,')') WRITE(LW,999) C IF(NTOT.GT.MTOT)STOP C C C A L C U L A T E V A R I O U S E L E M E N T A R R A Y S C CONNECTIVITY STIFFNESS FORCE C LM( NDOFE ) AK( NDOFE,NDOFE ) F( NDOFE ) C C ---------- AND ---------- C A S S E M B L E E L E M E N T A R R A Y S I N T O C G L O B A L A R R A Y S C STIFFNESS LOAD C A( NCOFS ) V( NEQ ) C C ADDRESSES N3=N2+NEQ1 N4=N3+NCOFS C C ID(NDS,3) MAXA(NEQ1) A(NCOFS) V(NEQ) CALL ELMNT(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),NE,NDS, * LTYPE,SLOAD,TLOAD,WLOAD,NCOFS,NEQ,NEQ1,NEQF,NDOFE) C C I F A V A I L A B L E C S U P E R P O S E N O D E F O R C E S C C ID(NDS,3) V(NEQ) IF(NLOAD.GT.0)CALL NDLOAD(SHARE(N1),SHARE(N4),NDS,NEQ,NLOAD) C C I F A V A I L A B L E I N C O R P O R A T E C P R E S C R I B E D N O D E D I S P L A C E M E N T S C C ID(NDS,3) MAXA(NEQ1) A(NCOFS) V(NEQ) IF(NDISP.GT.0)CALL NDDISP(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4), * NDS,NEQ1,NCOFS,NEQ,NDISP) C C********** S O L V E S Y S T E M O F E Q N. S ************* C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CALL COLSOL TWICE; 1ST TO TRIANGULARIZE A(NCOFS) C 2ND TO SOLVE FOR DISPLACEMENTS C LOADS STORED IN V(NEQ) REPLACED BY DISPLACEMENTS C C MAXA(NEQ1) A(NCOFS) V(NEQ) CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,1) CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,2) C C************* S Y S T E M O F E Q N. S S O L V E D ************ C C I F D I S P L A C E M E N T S W H E R E P R E S C R I B E D C R E C O V E R C O R R E S P O N D I N G N O D E F O R C E S C C ID(NDS,3) MAXA(NEQ1) A(NCOFS) V(NEQ) IF(NDISP.GT.0)CALL NDFORC(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4), * NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) C C R E C O V E R D I S P L A C E M E N T S U(NDS,3) C F R O M COLSOL S O L T I O N S S T O R E D I N A R R A Y V(NEQ) C C ID(NDS,3) U(NDS,3) V(NEQ) CALL DISP(SHARE(N1),SHARE(N2),SHARE(N4),NDS,NEQ,NEQF,LTYPE,UKEEP) C C C A L C U L A T E S T R E S S E S A N D S T R A I N S C C ADDRESSES N2=N1+3*NDS N3=N2+3*NDS N4=N3+NE N5=N4+NE N6=N5+NE N7=N6+NE N8=N7+NE C C SYMMETRIC STRAIN TENSOR COMPONENTS C W.R.T. LAMINATE AXIS (N,H,Z) C ( EN , EH , EZ , EHZ , ENZ , ENH ) C C U(NDS,3) EN(NE) EH(NE) EZ(NE) EHZ(NE) CALL STR(SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5),SHARE(N6), * SHARE(N7),SHARE(N8),NE,NDS,LTYPE,SLOAD,UKEEP) C ENZ(NE) ENH(NE) C REWIND 1 REWIND 2 REWIND 3 REWIND 4 CLOSE(1) CLOSE(2) CLOSE(3) CLOSE(4) 888 STOP END C*********************************************************************** SUBROUTINE NEDATA(ID,H,Z,NDS,NE) C CALLING PROGRAM C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * READS NODE - ELEMENT DATA FROM TAPE UNIT LR. * C * CALCULATES ELEMENT AREAS AND GEOMETRIC PARAMETERS * C * STORES PARAMENTERS FOR REFERENCE ON TAPES (1,2,3) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),H(NDS),Z(NDS) C C R E A D N O D E D A T A C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,10) 10 FORMAT(1H1,5X,'N O D E I N F O R M A T I O N',//, 16X,'(FIXED=1',/,8X,'FREE=0)',/,' NODE MOTION CO-ORDINATES',/, 2' NO. N H Z ',6X,'H',8X,'Z') C C L O O P T H R U A L L N O D E S C DO 40 I=1,NDS C C READ NODE DATA FROM INPUT FILE C READ(LR,20)ID(I,1),ID(I,2),ID(I,3),H(I),Z(I) 20 FORMAT(6X,3I3,2F12.0) C C WRITE NODE DATA ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,30)I,ID(I,1),ID(I,2),ID(I,3), * H(I),Z(I) 30 FORMAT(1X,I6,3(1X,I1,1X),1X,2F9.5) C 40 CONTINUE C C E N D N O D E L O O P C C C R E A D - S T O R E E L E M E N T D A T A C C C A L C U L A T E : RESCALED NODE CO-ORDINATES C ASSOCIATED WITH EACH ELEMENT C C C A L C U L A T E - S T O R E 1. ELEMENT FIBER ANGLE ORIENTATION C 2. ELEMENT GEOMETRIC PARAMETERS C AND AREAS C C INITIALIZE TOTAL GRID AREA C TAREA=0 C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,50) 50 FORMAT(1H1,5X,'E L E M E N T I N F O R M A T I O N',//, 1' ELMT------NODES----- PLY',/, 2' NO. 1ST 2ND 3RD ANGLE AREA') C C CALCULATE PI (RADIANS TO DEGREES CONVERSION) C PI=ACOS(-1.) C C L O O P T H R U A L L E L E M E N T S C DO 90 K=1,NE C C READ ELEMENT DATA FROM INPUT FILE C READ(LR,60)ND1,ND2,ND3,THETA 60 FORMAT(6X,3I6,F8.2) C C STORE NODES ASSOCIATED WITH EACH ELEMENT ON TAPE 1 C WRITE(1)ND1,ND2,ND3 C C STORE FIBER ANGLE (IN RADIANS) FOR EACH ELEMENT ON TAPE 2 C ANGLE ORIENTATION (CLOCKWISE ROTATION, POSITIVE) C WITH RESPECT TO LAMINATE NORMAL N-AXIS C THRAD=-THETA*PI/180. WRITE(2)THRAD C C NODE CO-ORDINATES ASSOCIATED WITH EACH ELEMENT RESCALED C H1=H(ND1) H2=H(ND2) H3=H(ND3) Z1=Z(ND1) Z2=Z(ND2) Z3=Z(ND3) C C ELEMENT AREA CALCULATED C TERM1=-(-H1+H2)*(Z1+Z2)/2.0 TERM2=-(-H2+H3)*(Z2+Z3)/2.0 TERM3=-(-H3+H1)*(Z3+Z1)/2.0 AREA=TERM1+TERM2+TERM3 C C ELEMENT GEOMETRIC PARAMETERS CALCULATED C AA=(Z2-Z3)/2.0 BB=(H3-H2)/2.0 CC=(Z3-Z1)/2.0 DD=(H1-H3)/2.0 EE=(Z1-Z2)/2.0 GG=(H2-H1)/2.0 C C STORE AREA AND GEOMETRIC PARAMETERS ON TAPE 3 C WRITE(3)AA,BB,CC,DD,EE,GG,AREA C C WRITE ELEMENT DATA ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,70)K,ND1,ND2,ND3,THETA,AREA 70 FORMAT(4(1X,I6),1X,F6.2,E11.4) C C CALCULATE TOTAL GRID AREA C TAREA=TAREA+AREA C 90 CONTINUE C C E N D E L E M E N T L O O P C C WRITE GRID AREA ON OUTPUT FILE C WRITE(LW,100)TAREA 100 FORMAT(//,5X,'TOTAL GRID AREA=',E15.8) C RETURN END C*********************************************************************** SUBROUTINE EQNUM(ID,NDS,NEQ,NEQF,LTYPE) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * GENERATES EQUATION NUMBERS BY REPLACING THE FREE DISPL * C * BOUNDARY CONDITION CODE (ID=0) WITH AN EQUATION NUMBER. * C * EQUATION NUMBERS ASSOCIATED WITH EACH FREE NODE DEGREE * C * OF FREEDOM (DOF) ARE STORED IN ARRAY ID(NDS,3) WHICH IS * C * SHARED WITH MAIN PROGRAM CONNOM SHARE(MTOT) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLE SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3) C C INITIALIZE VARIABLES DEFINING EQUATION NUMBERS C IEQN=0 NEQF=0 C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,10) 10 FORMAT(1H1,13X,'EQUATION NUMBERS',//,5(12X,'NODE MOTION '),/, * 5(' # NODE NORM HORIZ VERT')) C C SEQUENTIALLY TEST THE DISPLACEMENT BOUNDARY CONDITION C CODES( FIXED=1 OR FREE=0 ) FOR EACH NODE. THERE ARE C THREE POSSIBLE DEGREES OF FREEDOM ASSOCIATED WITH EACH NODE. C ID(I,1)=0, I-TH NODE IS FREE TO MOVE NORMAL TO THE H-Z PLANE C ID(I,2)=0, I-TH NODE IS FREE TO MOVE IN THE HORIZ. H-DIRECTION C ID(I,3)=0, I-TH NODE IS FREE TO MOVE IN THE VERT. Z-DIRECTION C IF A NODE IS FREE TO MOVE THEN AN EQUATION NUMBER (IEQN) IS C ASSOCIATED WITH THAT NODE DEGREE OF FREEDOM. C DO 40 I=1,NDS DO 20 J=1,3 IF(ID(I,J))6,5,6 5 IEQN=IEQN+1 ID(I,J)=IEQN GO TO 20 6 ID(I,J)=0 20 CONTINUE 40 CONTINUE C C WRITE EQN. NO.S ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,3)(I,(ID(I,J),J=1,3),I=1,NDS) 3 FORMAT(5(' #',I5,3I6)) C C LAST VALUE OF IEQN IS EQUAL TO THE TOTAL NO. OF EQUATIONS C NEQ=IEQN C C WRITE TOTAL NUMBER OF EQUATIONS DUE TO NODAL D.O.F. S C ON OUTPUT FILE C WRITE(LW,45)NEQ 45 FORMAT(//,' TOTAL NO. OF EQN.S DUE TO NODAL DEGREES OF FREEDOM', *I5,//) C C IF LOAD TYPE IS MECHANICAL LOADS ONLY C THEN RETURN C IF(LTYPE.LE.1)RETURN C C WHEN LOAD TYPE (LTYPE=2) IS EXCLUSIVELY THERMAL AND/OR MOISTURE C ADD ONE ADDITIONAL EQN TO THE SET OF SIMULTANEOUS EQUATIONS C NEQ=NEQ+1 C C NEQ = TOTAL NO. OF EQN.S TO BE SOLVED C NEQF = EQN. NO. WHICH RELATES UNKNOWN THERMAL OR MOISTURE C STRAIN (UKEEP) TO THE THERMAL OR MOISTURE LOAD B(NEQF) C C C ADDITIONAL EQN. NO. NEQF IS LOCATED AS THE LAST EQN. C OF THE ASSEMBLED (GLOBAL) STIFFNESS ARRAY SUCH THAT THE C COLUMN ORIENTED STORAGE OF THE NO. OF COEFF.S (NCOFS) IS MINIMIZED C NEQF=NEQ C C WRITE ADDITIONAL EQUATION NUMBER ON OUTPUT FILE C WRITE(LW,50)NEQF 50 FORMAT(' ADDITIONAL EQN. FOR THERMAL-MOISTURE LOAD, NEQF=',I6) C RETURN END C*********************************************************************** SUBROUTINE CIJ(NE) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATE AND STORE MATERIAL STIFFNESSES FOR * C * EACH ELEMENT OF THE GRID * C * * * * * * * * * * * * * * * * * * * * C REAL NU12,NU13,NU23,NU21,NU31,NU32 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E11,E22,E33,G12,G13,G23,NU12,NU13,NU23 C C CALCULATE MATERIAL STIFFNESSES C NU21=NU12*E22/E11 NU31=NU13*E33/E11 NU32=NU23*E33/E22 C DEL=1-NU12*NU21-NU23*NU32-NU31*NU13-2.0*NU21*NU32*NU13 C C11=(1.0-NU23*NU32)*E11/DEL C12=(NU21+NU31*NU23)*E11/DEL C13=(NU31+NU21*NU32)*E11/DEL C22=(1.0-NU13*NU31)*E22/DEL C23=(NU32+NU12*NU31)*E22/DEL C33=(1.0-NU12*NU21)*E33/DEL C44=G23 C55=G13 C66=G12 C C WRITE MATERIAL STIFFNESSES ON OUTPUT FILE C IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(LW,20)C11,C12,C13,C22,C23,C33, * C44,C55,C66 20 FORMAT(1H1,/,15X,'MATERIAL STIFFNESSES',//, 1' C11=',E11.4,' C12=',E11.4,' C13=',E11.4,/, 216X,' C22=',E11.4,' C23=',E11.4,/,32X,' C33=',E11.4,/, 348X,' C44=',E11.4,/,64X,' C55=',E11.4,/,80X,' C66=',E11.4) C C STIFFNESS ARRAY TRANSFORMATIONS USING THRAD=THETA(RADIANS) C FOR THE K-TH ELEMENT C REWIND 2 C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(LW,25) 25 FORMAT(//,15X,'ELEMENT MATERIAL STIFFNESSES') C DO 10 K=1,NE READ(2)THRAD C=COS(THRAD) SN=SIN(THRAD) C2=C*C C4=C2*C2 S2=SN*SN S4=S2*S2 D11=C11*C4+2.*C2*S2*(C12+2.*C66)+C22*S4 D12=C2*S2*(C11+C22-4.*C66)+C12*(S4+C4) D13=C13*C2+C23*S2 D16=-C*SN*(C11*C2-C22*S2-(C12+2.*C66)*(C2-S2)) D22=C11*S4+2.*C2*S2*(C12+2.*C66)+C22*C4 D23=C13*S2+C23*C2 D26=-C*SN*(C11*S2-C22*C2+(C12+2.*C66)*(C2-S2)) D33=C33 D36=C*SN*(C23-C13) D44=C44*C2+C55*S2 D45=C*SN*(C44-C55) D55=C55*C2+C44*S2 D66=(C11+C22-2.*C12)*C2*S2+C66*(C2-S2)**2 C C STORE ELEMENT STIFFNESSES ON TAPE 4 C WRITE(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66 C C WRITE ELEMENT STIFFNESSES ON OUTPUT FILE C IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(LW,35)K,D11,D12,D13,D16,D22,D23, *D26,D33,D36,D44,D45,D55,D66 35 FORMAT(' ELEMENT NO.=',I4,/,' D11=',E11.4,' D12=',E11.4, 1' D13=',E11.4,33X,'D16=',E11.4,/,16X,' D22=',E11.4,' D23=',E11.4, 2,32X,' D26=',E11.4,/,32X,' D33=',E11.4,32X,' D36=',E11.4,/, 348X,' D44=',E11.4,' D45=',E11.4,/,64X,' D55=',E11.4,/,80X, 4' D66=',E11.4) C 10 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE ADDRS(ID,MAXA,NE,NDS,NEQ1,LTYPE,NEQF,MBAND,NCOFS,NDOFE) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. DETERMINES COLUMN HEIGHTS, STORED TEMPORARILY * C * IN ARRAY MAXA(NEQ). * C * 2. DETERMINES ADDRESSES = ( I-TH SUBSCRIPT IN A(I) ) FOR * C * THE DIAGONAL TERMS IN THE GLOBAL ARRAY A(I) * C * ADDRESSES STORED IN MAXA(NEQ+1). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES ID AND MAXA SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),MAXA(NEQ1),LM(10) C C CLEAR LM AND MAXA ARRAYS C DO 5 IEQN=1,NEQ1 5 MAXA(IEQN)=0 DO 10 I=1,10 10 LM(I)=0 C C USING EQN. NO.S STORED IN ID(NDS,3) ARRAY TOGETHER WITH NODE C NO.S (ND1,ND2,ND3) STORED ON TAPE 1 CALCULATE CONNECTIVITY C LM ARRAYS FOR EACH ELEMENT. C REWIND 1 C C L O O P T H R U A L L E L E M E N T S C DO 40 K=1,NE C C READ NODE NO.S ASSOCIATED WITH K-TH ELEMENT C READ(1)ND1,ND2,ND3 C C ASSIGN EQUATION NUMBERS TO ELEMENT CONNECTIVITY ARRAY C C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 1ST NODE)=ID(ND1,1)=LM(1) C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 2ND NODE)=ID(ND2,1)=LM(2) C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 3RD NODE)=ID(ND3,1)=LM(3) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 1ST NODE)=ID(ND1,2)=LM(4) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 2ND NODE)=ID(ND2,2)=LM(5) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 3RD NODE)=ID(ND3,2)=LM(6) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 1ST NODE)=ID(ND1,3)=LM(7) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 2ND NODE)=ID(ND2,3)=LM(8) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 3RD NODE)=ID(ND3,3)=LM(9) C LM(1)=ID(ND1,1) LM(2)=ID(ND2,1) LM(3)=ID(ND3,1) LM(4)=ID(ND1,2) LM(5)=ID(ND2,2) LM(6)=ID(ND3,2) LM(7)=ID(ND1,3) LM(8)=ID(ND2,3) LM(9)=ID(ND3,3) C C ******* EXCLUSIVELY ******* C FOR THERMAL - MOISTURE LOADS ONLY INCLUDE EXTRA DOF C BY ASSINGING EQN NO. NEQF TO LM(10) FOR EACH ELEMENT C IF(LTYPE.EQ.2)LM(10)=NEQF C C CHOOSE NO. OF DEGREES OF FREEDOM PER ELEMENT, NDOFE C C DEFINE NDOFE=9 1. FOR A STRAIN APPLIED NORMAL TO THE C ELEMENT H-Z PLANE, SLOAD C 2. FOR NODE LOADS OR PRESCRIBED DISPL.S NDOFE=9 C C EXCLUSIVELY DEFINE C NDOFE=10 1. FOR A CHANGE IN TEMPERATURE, TLOAD C 2. FOR A CHANGE IN ABSORBED MOISTURE, WLOAD C IF(LTYPE.EQ.2)NDOFE=10 C C CALCULATE THE COLUMN HEIGHTS ABOVE THE DIAGONAL C OF THE ASSEMBLED GLOBAL STIFFNESS ARRAY C (COLUMN HEIGHTS TEMPORARILY STORED IN MAXA ARRAY) C C THIS IS DONE BY: C COMPARING VALUES IN LM ARRAYS FOR ELEMENTS C WHICH SHARE A DEGREE OF FREEDOM(DOF) C COLUMN HEIGHT = MAXA( SHARED DOF ) - LOWEST DOF IN SHARED LM ARRAY C DO 30 I=1,NDOFE L=LM(I) IF(L.EQ.0)GO TO 30 DO 20 J=1,NDOFE M=LM(J) IF(M.EQ.0)GO TO 20 IF(M.GT.L)GO TO 20 N=L-M+1 IF(N.GT.MAXA(L))MAXA(L)=N 20 CONTINUE 30 CONTINUE 40 CONTINUE C C E N D E L E M E N T L O O P C C CALCUALTE WIDTH OF HALFBAND PLUS DIAGONAL, MBAND C WRITE COLUMN HEIGHTS AND MBAND ON OUTPUT FILE C MBAND=0 NEQ=NEQ1-1 DO 55 I=1,NEQ IF(MAXA(I).GT.MBAND)MBAND=MAXA(I) 55 CONTINUE IF((KEY.EQ.4).OR.(KEY.EQ.7))WRITE(LW,50)(I,MAXA(I),I=1,NEQ1) 50 FORMAT(4(' # COLUMN=',I6,' HEIGHT=',I6,' #')) WRITE(LW,56)MBAND 56 FORMAT(' WIDTH OF HALFBAND PLUS DIAGONAL, MBAND=',I4) C C USING COLUMN HEIGHTS STORED IN MAXA( NEQ ) C DETERMINE ADDRESSES = ( I-TH SUBSCRIPT IN A(I) ) C FOR DIAGONAL TERMS IN THE GLOBAL ARRAY A(I) C AND STORE VALUES IN MAXA( NEQ + 1 ). C C FIRST COLUMN DIAGONAL ADDRESS BY DEFINITION IS ONE C NN=1 MAXA(1)=1 MAXM1=1 C C ADD COLUMN HEIGHTS TO PREVIOUS DIAGONAL ADDRESS C DO 60 I=2,NEQ1 NN=NN+MAXM1 MAXM1=MAXA(I) MAXA(I)=NN 60 CONTINUE C C WRITE GLOBAL ARRAY SIZE AND DIAGONAL ADDRESSES ON THE OUTPUT FILE C NCOFS=MAXA(NEQ1) WRITE(LW,70)NCOFS 70 FORMAT(/,' NO. OF COEFF.S IN GLOBAL STIFFNESS ARRAY, NCOFS=',I6,/) C IF((KEY.EQ.6).OR.(KEY.EQ.7))WRITE(LW,75)(I,MAXA(I),I=1,NEQ) 75 FORMAT(3(' # COLUMN=',I6,' DIAGONAL ADDRESS=',I6,' #')) C RETURN END C*********************************************************************** SUBROUTINE ELMNT(ID,MAXA,A,V,NE,NDS,LTYPE,SLOAD,TLOAD,WLOAD, * NCOFS,NEQ,NEQ1,NEQF,NDOFE) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATES THERMAL-MOISTURE STRAINS. * C * CALCULATES ELEMENT ARRAYS LM(I),AK(I,J),AND F(I). * C * ASSEMBLES GLOBAL ARRAYS A(I) AND V(I) FROM ELEMENT ARRAYS. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY COMMON/RDK3/AL1,AL2,AL3,BT1,BT2,BT3 C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ) C C DIMENSION ELEMENT ARRAYS C DIMENSION AK(10,10),F(10),LM(10) C C CLEAR A(NCOFS) AND V(NEQ) ARRAYS C DO 1 I=1,NCOFS 1 A(I)=0.0 DO 2 I=1,NEQ 2 V(I)=0.0 C C C A L C U L A T E T H E R M A L - M O I S T U R E S T R A I N S C C STRN1, STRAIN PARALLEL TO PLY FIBER DIRECTION C STRN2, STRAIN NORMAL TO PLY FIBER DIRECTION C STRN12, SHEAR STRAIN IN THE THIN PLY 1-2 PLANE C STRN3, STRAIN NORMAL TO THE PLY 1-2 PLANE C IF(LTYPE.LE.1)GO TO 5 C STRN1=AL1*TLOAD+BT1*WLOAD STRN2=AL2*TLOAD+BT2*WLOAD STRN12=TLOAD*(AL2-AL1)+WLOAD*(BT2-BT1) STRN3=AL3*TLOAD+BT3*WLOAD C C R E W I N D S T O R E D E L E M E N T D A T A C C TAPE 1 ND1,ND2,ND3 ....... NODE NO.S ASSIGNED TO EACH ELEMENT 5 REWIND 1 C C TAPE 2 THRAD=THETA(RADIANS) .. ELEMENT FIBER AXIS ORIENTATION REWIND 2 C C TAPE 3 AA,BB,...,GG,ARR ...... ELEMENT GEOMETRIC PARAMETERS REWIND 3 C C TAPE 4 D11,D12,...,D66 ....... ELEMENT STIFFNESSES REWIND 4 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * L O O P T H R U A L L E L E M E N T S * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DO 70 K=1,NE C C F O R M A T I O N O F C O N N E C T I V I T Y A R R A Y C LM ( NDOFE ) C C READ NODE NO.S ASSOCIATED WITH THE K-TH ELEMENT C READ(1)ND1,ND2,ND3 C C ASSIGN EQUATION NUMBERS TO ELEMENT CONNECTIVITY ARRAY C C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 1ST NODE)=ID(ND1,1)=LM(1) C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 2ND NODE)=ID(ND2,1)=LM(2) C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 3RD NODE)=ID(ND3,1)=LM(3) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 1ST NODE)=ID(ND1,2)=LM(4) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 2ND NODE)=ID(ND2,2)=LM(5) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 3RD NODE)=ID(ND3,2)=LM(6) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 1ST NODE)=ID(ND1,3)=LM(7) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 2ND NODE)=ID(ND2,3)=LM(8) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 3RD NODE)=ID(ND3,3)=LM(9) C LM(1)=ID(ND1,1) LM(2)=ID(ND2,1) LM(3)=ID(ND3,1) LM(4)=ID(ND1,2) LM(5)=ID(ND2,2) LM(6)=ID(ND3,2) LM(7)=ID(ND1,3) LM(8)=ID(ND2,3) LM(9)=ID(ND3,3) C C ******* EXCLUSIVELY ******* C FOR A THERMAL - MOISTURE LOADS ONLY INCLUDE EXTRA DOF C BY ASSIGNING EQN. NO. NEQF TO LM(10) FOR EACH ELEMENT C IF(LTYPE.EQ.2)LM(10)=NEQF C C WRITE LM ARRAY FOR EACH ELEMENT ON OUTPUT FILE C IF((KEY.EQ.3).OR.(KEY.EQ.7))GO TO 9 GO TO 15 C 9 IF(NDOFE.EQ.9)WRITE(LW,10)K,NDOFE,(LM(I),I=1,NDOFE) IF(NDOFE.EQ.10)WRITE(LW,11)K,NDOFE,(LM(I),I=1,NDOFE) 10 FORMAT(1H1,' ELEMENT NO.',I4,' CONNECTIVITY ARRAY LM(',I2,')', *9(I4,1X)) 11 FORMAT(1H1,' ELEMENT NO.',I4,' CONNECTIVITY ARRAY LM(',I2,')', *10(I4,1X)) C C F O R M A T I O N O F C E L E M E N T A L S T I F F N E S S C A R R A Y AK ( NDOFE,NDOFE ) C ------- AND ------- C E L E M E N T A L R I G H T - H A N D - S I D E C F O R C E V E C T O R F ( NDOFE ) C C READ ELEMENT DATA C 15 READ(3)AA,BB,CC,DD,EE,GG,ARR READ(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66 C C CALCULATE STIFFNESS ARRAY FOR EACH ELEMENT C AK(1,1)=(D55*BB*BB+D66*AA*AA)/ARR AK(1,2)=(D55*BB*DD+D66*AA*CC)/ARR AK(1,3)=(D55*BB*GG+D66*AA*EE)/ARR AK(1,4)=(D26*AA*AA+D45*BB*BB)/ARR AK(1,5)=(D26*CC*AA+D45*BB*DD)/ARR AK(1,6)=(D26*EE*AA+D45*BB*GG)/ARR AK(1,7)=(D36*BB*AA+D45*BB*AA)/ARR AK(1,8)=(D36*DD*AA+D45*BB*CC)/ARR AK(1,9)=(D36*GG*AA+D45*BB*EE)/ARR AK(2,2)=(D55*DD*DD+D66*CC*CC)/ARR AK(2,3)=(D55*DD*GG+D66*CC*EE)/ARR AK(2,4)=(D26*AA*CC+D45*DD*BB)/ARR AK(2,5)=(D26*CC*CC+D45*DD*DD)/ARR AK(2,6)=(D26*EE*CC+D45*DD*GG)/ARR AK(2,7)=(D36*BB*CC+D45*DD*AA)/ARR AK(2,8)=(D36*DD*CC+D45*DD*CC)/ARR AK(2,9)=(D36*GG*CC+D45*DD*EE)/ARR AK(3,3)=(D55*GG*GG+D66*EE*EE)/ARR AK(3,4)=(D26*AA*EE+D45*GG*BB)/ARR AK(3,5)=(D26*CC*EE+D45*GG*DD)/ARR AK(3,6)=(D26*EE*EE+D45*GG*GG)/ARR AK(3,7)=(D36*BB*EE+D45*GG*AA)/ARR AK(3,8)=(D36*DD*EE+D45*GG*CC)/ARR AK(3,9)=(D36*GG*EE+D45*GG*EE)/ARR AK(4,4)=(D22*AA*AA+D44*BB*BB)/ARR AK(4,5)=(D22*AA*CC+D44*BB*DD)/ARR AK(4,6)=(D22*AA*EE+D44*BB*GG)/ARR AK(4,7)=(D44*BB*AA+D23*AA*BB)/ARR AK(4,8)=(D44*BB*CC+D23*AA*DD)/ARR AK(4,9)=(D44*BB*EE+D23*AA*GG)/ARR AK(5,5)=(D22*CC*CC+D44*DD*DD)/ARR AK(5,6)=(D22*CC*EE+D44*DD*GG)/ARR AK(5,7)=(D44*DD*AA+D23*CC*BB)/ARR AK(5,8)=(D44*DD*CC+D23*CC*DD)/ARR AK(5,9)=(D44*DD*EE+D23*CC*GG)/ARR AK(6,6)=(D22*EE*EE+D44*GG*GG)/ARR AK(6,7)=(D44*GG*AA+D23*EE*BB)/ARR AK(6,8)=(D44*GG*CC+D23*EE*DD)/ARR AK(6,9)=(D44*GG*EE+D23*EE*GG)/ARR AK(7,7)=(D33*BB*BB+D44*AA*AA)/ARR AK(7,8)=(D33*BB*DD+D44*AA*CC)/ARR AK(7,9)=(D33*BB*GG+D44*AA*EE)/ARR AK(8,8)=(D33*DD*DD+D44*CC*CC)/ARR AK(8,9)=(D33*DD*GG+D44*CC*EE)/ARR AK(9,9)=(D33*GG*GG+D44*EE*EE)/ARR C C CALCULATE THE RIGHT-HAND-SIDE FORCE VECTOR DUE TO A C UNIFORMLY APPLIED STRAIN (SLOAD) NORMAL TO THE H-Z PLANE C F(1)=-D16*AA*SLOAD F(2)=-D16*CC*SLOAD F(3)=-D16*EE*SLOAD F(4)=-D12*AA*SLOAD F(5)=-D12*CC*SLOAD F(6)=-D12*EE*SLOAD F(7)=-D13*BB*SLOAD F(8)=-D13*DD*SLOAD F(9)=-D13*GG*SLOAD C C IF LOAD TYPE IS A MECHANICALLY APPLIED STRAIN NORMAL TO C THE ELEMENT PLANE THEN LOOP AROUND ADDITIONAL EQUATIONS C IF(LTYPE.LE.1)GO TO 20 C C FOR THERMAL-MOISTURE LOADING (EXCLUSIVELY) C 1. INCLUDE TERMS IN THE TENTH ROW/COLUMN OF STIFFNESS ARRAY C AK(1,10)=D16*AA AK(2,10)=D16*CC AK(3,10)=D16*EE AK(4,10)=D12*AA AK(5,10)=D12*CC AK(6,10)=D12*EE AK(7,10)=D13*BB AK(8,10)=D13*DD AK(9,10)=D13*GG AK(10,10)=D11*ARR C C 2. CALCULATE RIGHT-HAND-SIDE FORCE VECTORS C C 2.1 TRANSFORM C THERMAL-MOISTURE STRAINS C FROM MATERIAL AXIS(1,2,3) TO LAMINATE AXIS(N,H,Z) C READ(2)THRAD C C=COS(THRAD) S=SIN(THRAD) ENO=C*C*STRN1+S*S*STRN2 EHO=S*S*STRN1+C*C*STRN2 ENHO=2.0*S*C*STRN12 EZO=STRN3 C C 2.2 STORE THERMAL-MOISTURE STRAINS ON TAPE 7 C WRITE(7)ENO,EHO,ENHO,EZO C C 2.3 CALCULATE LOADS DUE TO THERMAL-MOISTURE STRAINS C F(1)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*AA F(2)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*CC F(3)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*EE F(4)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*AA F(5)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*CC F(6)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*EE F(7)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*BB F(8)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*DD F(9)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*GG F(10)=(D11*ENO+D12*EHO+D13*EZO+D16*ENHO)*ARR C C ASSIGN VALUES TO SYMMETRIC PORTION OF ELEMENTAL C STIFFNESS ARRAY AK(I,J)=AK(J,I) C 20 JJ=1 DO 40 II=1,NDOFE DO 30 J=JJ,NDOFE 30 AK(J,II)=AK(II,J) 40 JJ=JJ+1 C C WRITE STIFFNESS ARRAY AND FORCE VECTOR ON OUTPUT FILE C IF((KEY.EQ.3).OR.(KEY.EQ.7))GO TO 45 GO TO 60 C 45 IF(NDOFE.EQ.10)GO TO 55 C WRITE(LW,50)K 50 FORMAT(//,15X,'ELEMENT NO.',I4,/,12X,'STIFFNESS ARRAY AK(I,J)', */,' ROW/COL 1',9X,'2',9X,'3',9X,'4',9X,'5',9X,'6',9X,'7',9X,'8' *,9X,'9',5X,'FORCE') DO 51 I=1,NDOFE 51 WRITE(LW,52)I,(AK(I,J),J=1,NDOFE),F(I) 52 FORMAT(1X,I4,3X,10(E10.3)) GO TO 60 C 55 WRITE(LW,56)K 56 FORMAT(//,15X,'ELEMENT NO.',I4,/,12X,'STIFFNESS ARRAY AK(I,J)', */,' ROW/COL 1',9X,'2',9X,'3',9X,'4',9X,'5',9X,'6',9X,'7',9X,'8' ,9X,'9',9X,'10',5X,'FORCE') DO 57 I=1,NDOFE 57 WRITE(LW,58)I,(AK(I,J),J=1,NDOFE),F(I) 58 FORMAT(1X,I4,3X,11(E10.3)) C C A S S E M B L E C ELEMENTAL ARRAYS ( AK , F ) INTO GLOBAL ARRAYS ( A , V ) C USING THE CONNECTIVITY AND DIAGONAL ADDRESS ARRAYS ( LM , MAXA ) C 60 CALL ASSEMB(NDOFE,AK,F,LM,MAXA,A,V,NEQ1,NCOFS,NEQ) C 70 CONTINUE C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * E N D E L E M E N T L O O P * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C WRITE ASSEMBLED 1. GLOBAL STIFFNESS ARRAY ON OUTPUT FILE C 2. RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE C IF((KEY.EQ.4).OR.(KEY.EQ.7))GO TO 74 RETURN C C WRITE HEADER AND FIRST DIAGONAL TERM ON OUTPUT FILE C 74 WRITE(LW,75)A(1) 75 FORMAT(1H1,//,' GLOBAL STIFFNESS ARRAY',//,' COLUMN NO. 1', */,7X,' A( 1)=',E11.4,' 1 RST DIAGONAL POSITION') C C START LOOP AT 2 ND DIAGONAL TERM C IDIAG=2 C C LOOP THRU A(I) TERMS C DO 80 I=2,NCOFS IF(I.EQ.MAXA(IDIAG))GO TO 77 WRITE(LW,76)I,A(I) 76 FORMAT(8X,'A(',I6,')=',E14.7) GO TO 80 77 IF(I.NE.NCOFS)WRITE(LW,78)IDIAG,I,A(I),IDIAG 78 FORMAT(/,' COLUMN NO.',I6,/,8X,'A(',I6,')=',E14.7, *' DIAGONAL NO.',I6) C IF(I.EQ.NCOFS)WRITE(LW,79)I,A(I) 79 FORMAT(//,8X,'A(',I6,')=',E14.7,' LAST STORED VALUE IN A(I)') C IDIAG=IDIAG+1 C 80 CONTINUE C C WRITE RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE C WRITE(LW,81) 81 FORMAT(1H1,//,' RIGHT-HAND-SIDE LOAD VECTOR',//) C DO 82 I=1,NEQ 82 WRITE(LW,83)I,V(I) 83 FORMAT(7X,'V(',I6,')=',E14.7) C RETURN END C*********************************************************************** SUBROUTINE ASSEMB(NDOFE,AK,F,LM,MAXA,A,V,NEQ1,NCOFS,NEQ) C CALLING PROGRAM - ELMNT C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * ASSEMBLES : 1. EACH ELEMENT STIFFNESS ARRAY AK(I,J) * C * ONTO THE GLOBAL ARRAY A(I). * C * 2. EACH ELEMENT FORCE VECTOR F(I) ONTO THE * C * RIGHT-HAND-SIDE LOAD VECTOR V(I). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DIMENSION AK(10,10),F(10),LM(10) DIMENSION MAXA(NEQ1),A(NCOFS),V(NEQ) C DO 20 I=1,NDOFE L=LM(I) IF(L.EQ.0)GO TO 20 NN=MAXA(L) DO 10 J=I,NDOFE M=LM(J) IF(M.EQ.0)GO TO 10 MM=MAXA(M) KK=M-L MM=MM+KK IF(KK.LT.0)MM=NN-KK C A(MM)=A(MM)+AK(I,J) 10 CONTINUE C V(L)=V(L)+F(I) 20 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE NDLOAD(ID,V,NDS,NEQ,NLOAD) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * ADDS NODAL LOADS TO THE EXISTING ASSEMBLED GLOBAL * C * RIGHT-HAND-SIDE LOAD VECTOR V(NEQ) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),V(NEQ) C C USING EQN NO.S STORED IN ID(NDS,2) ARRAY C RECOVER EQN NO.(IEQN) WHICH CORRESPONDS TO NODAL LOAD C DIRECTION (IDIRN) AND ADD NODE LOAD (FLOAD) TO THE C RIGHT-HAND-SIDE LOAD VECTOR (V(IEQN)). C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NLOAD 10 FORMAT(1H1,//,' NUMBER OF APPLIED LOADS=',I4,/) C C L O O P T H R U L O A D S C DO 60 ILOAD=1,NLOAD C C READ 1. NOD , NODE NO. WHERE I-TH LOAD IS APPLIED C 2. IDIRN, DIRECTION OF I-TH LOAD C 3. FLOAD, MAGNITUDE OF I-TH LOAD C READ(LR,20)NOD,IDIRN,FLOAD 20 FORMAT(2I6,F12.0) C C RECOVER EQN NO. C IEQN=ID(NOD,IDIRN) C WRITE(LW,30)ILOAD,NOD,IDIRN,IEQN,FLOAD,IEQN,V(IEQN) 30 FORMAT(/,5X,'ILOAD=',I4,' NODE NO.=',I4,' LOAD DIRECTION=', *I3,' EQN NO.=',I5,/,10X,'LOAD MAGNITUDE=',E10.3,/, *15X,'OLD RIGHT-HAND-SIDE LOAD VECTOR, V(',I4,')=',E11.4) C C ADD NODAL LOAD (FLOAD) TO EXISTING LOAD VECTOR C V(IEQN)=V(IEQN)+FLOAD C WRITE(LW,50)IEQN,V(IEQN) 50 FORMAT(15X,'NEW RIGHT-HAND-SIDE LOAD VECTOR,V(',I4,')=',E11.4,/) C 60 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE NDDISP(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M C * ADDS A HUGE SPRING STIFFNESS (HUGEK) TO THE ASSEMBLED * C * STIFFNESS (A(MAXA(IEQN))) AND LOAD (V(IEQN)) COMPONENTS * C * SO AS TO SIMULATE A PRESCRIBED NODAL DISPLACEMENT. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ) C C USING EQN NO.S STORED IN ID(NDS,3) ARRAY C RECOVER EQN NO. (IEQN) CORRESPONDING TO NODE DISPLACEMENT (DISP). C THEN 1. ADD LARGE SPRING STIFFNESS (HUGEK) TO DIAGONAL C TERN ( A(MAXA(IEQN)) ) OF THE GLOBAL STIFFNESS ARRAY. C 2. ADD LARGE SPRING FORCE ( HUGEK*DISP ) TO THE C RIGHT-HAND-SIDE LOAD VECTOR ( V(IEQN) ). C HUGEK=1.0E+20 C C PRIOR TO ADDITION OF HUGE SPRING STIFFNESS C STORE GLOBAL STIFFNESS ARRAY ON TAPE 8 C WRITE(8)(A(I),I=1,NCOFS) C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NDISP 10 FORMAT(1H1,//,' NUMBER OF PRESCRIBED DISPL.S=',I4,/) C C L O O P T H R U P R E S C R I B E D D I S P L. S C DO 60 IDISP=1,NDISP C C READ 1. NOD , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED C 2. IDIRN, DIRECTION OF I-TH DISPL. C 3. DISP , MAGNITUDE OF I-TH DISPL. C READ(LR,20)NOD,IDIRN,DISP 20 FORMAT(2I6,F12.5) C C STORE PRIESCRIBED DISPL. DATA ON TAPE 9 C WRITE(9)NOD,IDIRN,DISP C C RECOVER EQN NO. (IEQN) AND DIAGONAL ADDRESSES MAXA(IEQN) C IEQN=ID(NOD,IDIRN) II=MAXA(IEQN) C WRITE(LW,30)IDISP,NOD,IDIRN,IEQN,DISP,II,A(II),IEQN,V(IEQN) 30 FORMAT(/,5X,' IDISP=',I6,' NODE NO.=',I6,' DISPL. DIRECTION=', *I3,/,6X,'EQN NO.=',I6,14X,'DISPL. MAGNITUDE=',E10.3,//, *5X,'OLD DIAGONAL A(',I6,')=',E11.4, *' OLD RIGHT-HAND-SIDE VECTOR V(',I6,')=',E11.4) C C ADD SPRING STIFFNESS AND SPRING FORCE C A(II)=A(II)+HUGEK V(IEQN)=V(IEQN)+HUGEK*DISP C WRITE(LW,50)II,A(II),IEQN,V(IEQN) 50 FORMAT(5X,'NEW DIAGONAL A(',I6,')=',E11.4, *' NEW RIGHT-HAND-SIDE VECTOR V(',I6,')=',E11.4,/) C 60 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE COLSOL(MAXA,A,V,NN,NWK,NNM,KKK) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * C * THIS SUBROUTINE IS CALLED TWICE * C * * C * THE FIRST TIME * C * * TO PERFORM L D L TRANSPOSE OPERATION * C * ON THE GLOBAL STIFFNESS MATRIX * C * * C * THE SECOND TIME * C * TO SOLVE THE SYSTEM OF EQUATIONS * C * * C * PROGRAM * C * TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN * C * CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME * C * - - INPUT VARIABLES - - * C * A(NWK) = STIFFNESS MATRIX STORED IN COMPACTED FORM * C * V(NN) = RIGHT-HAND-SIDE LOAD VECTOR * C * MAXA(NNM) = VECTOR CONTAINING ADDRESSES OF DIAGONAL * C * ELEMENTS OF STIFFNESS MATRIX IN A * C * NN = NEQ = NUMBER OF EQUATIONS * C * NWK (NOT USED)= NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX * C * NNM (NOT USED)= NN + 1 * C * KKK = INPUT FLAG * C * EQ. 1 TRIANLGULARIZATION OF STIFFNESS MATRIX * C * EQ. 2 REDUCTION AND BACK-SUBSTITUTION OF LOAD VECTOR * C * IOUT = NUMBER OF OUTPUT DEVICE * C * - - OUTPUT - - * C * A(NWK) = D AND L - FACTORS OF STIFFNESS MATRIX * C * V(NN) = DISPLACEMENT VECTOR * C * * * * * * * * * * * * * * * * * * * * * * * * C C IMPLICIT REAL*8(A-H,0-Z) C C * * * * * * * * * * * * * * * * * * * * * * * * C * THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON * C * CDC EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM * C * OR UNIVAC MACHINES. ACTIVATE, DEACTIVATE OR ADJUST ABOVE * C * CARD FOR SINGLE OR DOUBLE PRECISION ARITHMETIC. * C * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/IOUT,LR,KEY C DIMENSION MAXA(NNM),A(NWK),V(NN) C C *** PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX C IF (KKK-2) 40,150,150 40 DO 140 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL IF (KH) 110,90,50 50 K=N-KH IC=0 KLT=KU DO 80 J=1,KH IC=IC+1 KLT=KLT-1 KI=MAXA(K) ND=MAXA(K+1)-KI-1 IF (ND) 80,80,60 60 KK=IC IF (ND.LT.IC)KK=ND C=0.0 DO 70 L=1,KK 70 C=C+A(KI+L)*A(KLT+L) A(KLT)=A(KLT)-C 80 K=K+1 90 K=N B=0.0 DO 100 KK=KL,KU K=K-1 KI=MAXA(K) C=A(KK)/A(KI) B=B+C*A(KK) 100 A(KK)=C A(KN)=A(KN)-B 110 IF (A(KN)) 120,120,140 120 WRITE (IOUT,2000) N,A(KN) STOP 140 CONTINUE RETURN C C *** REDUCE RIGHT-HAND-SIDE LOAD VECTOR C 150 DO 180 N=1,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF (KU-KL) 180,160,160 160 K=N C=0.0 DO 170 KK=KL,KU K=K-1 170 C=C+A(KK)*V(K) V(N)=V(N)-C 180 CONTINUE C C *** BACK SUBSTITUTE C DO 200 N=1,NN K=MAXA(N) 200 V(N)=V(N)/A(K) IF (NN.EQ.1) RETURN N=NN DO 230 L=2,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF (KU-KL) 230,210,210 210 K=N DO 220 KK=KL,KU K=K-1 220 V(K)=V(K)-A(KK)*V(N) 230 N=N-1 RETURN C 2000 FORMAT (//48H STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE ,/ */ 32H NONPOSITIVE PIVOT FOR EQUATION ,I4// * 10H PIVOT = ,E20.12) END C*********************************************************************** SUBROUTINE NDFORC(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * RECOVER FORCES AT NODES CORRESPONDING TO PRESCRIBED * C * NODAL DISPLACEMENTS. FORCES CALCULATED FROM COMPONENTS * C * OF GLOBAL ARRAY A(I) STORED ON TAPE 8 PRIOR TO SUPERPOSING * C * A LARGE SPRING STIFFNESS (HUGEK) ON DIAGONAL TERMS. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ) C C CALCULATE BANDWITH (MB) NOT INCLUDING DIAGONAL TERM C CALCULATE (LIMIT) EQN NO. WHERE TRUNCATION OF TERMS TO C RIGHT OF THE DIAGONAL OCCURS FOR (IEQN.GT.LIMIT) C MB=MBAND-1 LIMIT=NEQ1-MBAND C C REPLACE GLOBAL STIFFNESS ARRAY TERMS A(NCOFS) USED IN COLSOL WITH C VALUES STORED ON TAPE 8 PRIOR TO THE ADDITION OF HUGEK STIFFNESS. C REWIND 8 READ(8)(A(I),I=1,NCOFS) C C USING EQN NO.S STORED IN ID(NDS,3) AND MAXA(NEQ+1) ARRAYS C RECOVER EQN NO. (IEQN) CORRESPONDING TO PRESCRIBED NODAL DISPL. C C TO CALCULATE THE FORCE IN THE DIRECTION OF THE PRESCRIBED C DISPLACEMENT, SUM UP THE PRODUCT OF KNOWN DISPLACEMENTS STORED C IN ARRAY V(I) WITH THE ROW OF STIFFNESSES A(I) TERMS STORED C ON TAPE 8, WHICH COMPRISES EQN NO. (IEQN). C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NDISP 10 FORMAT(1H1,//,' NUMBER OF FORCES CORRESPONDING TO PRESCRIBED NODA *L DISPLACEMENTS=',I4) WRITE(LW,999) 999 FORMAT(/,80(1H*),/) C REWIND 9 C C L O O P T H R U P R E S C R I B E D D I S P L. S C DO 70 IDISP=1,NDISP C C RECOVER PRESCRIBED DISPL. DATA STORED ON TAPE 9 C NOD , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED C IDIRN, DIRECTION OF I-TH DISPL. C DISP , MAGNITUDE OF I-TH DISPL. C READ(9)NOD,IDIRN,DISP C IEQN=ID(NOD,IDIRN) II=MAXA(IEQN) C C INITIALIZE VARIABLES USED TO POSITION SUMMATION OF TERMS C AT LEFT SIDE OF EQN NO. (IEQN) C C N , NO. OF COMPONENTS TO THE LEFT AND INCLUDING DIAGONAL C IN, SUBSCRIPT OF FIRST DISPL. IN EQN NO. (IEQN) C IM, SUBSCRIPT OF FIRST STIFFNESS COMPONENT IN EQN NO. (IEQN) C N=MAXA(IEQN+1)-II IN=IEQN-N+1 IM=MAXA(IEQN+1)-1 C C SUM UP PRODUCTS OF DISPLACEMENTS V(I) AND STIFFNESSES A(I) C FOR TERMS TO THE LEFT AND INCLUDING DIAGONAL C FORCE=0.0 C DO 30 K=1,N FORCE=FORCE+V(IN)*A(IM) C IF(KEY.EQ.7)WRITE(LW,20)IN,V(IN),IM,A(IM),FORCE 20 FORMAT(/,' DISPLACEMENT, V(',I4,')=',E11.4, *' STIFFNESS, A(',I6,')=',E11.4,' FORCE=',E11.4) C IM=IM-1 IN=IN+1 30 CONTINUE C C SUM UP REMAINING PRODUCT OF TERMS TO RIGHT OF DIAGONAL. C BUT FIRST RECALCULATE TRUNCATED BANDWIDTH (MB) TO THE RIGHT C OF THE DIAGONAL IF THE EQN NO. (IEQN) IS GREATER THAN (LIMIT) C IF(IEQN.GT.LIMIT)MB=MBAND-1-IEQN+LIMIT IF(MB.EQ.0)GO TO 55 C DO 50 IB=1,MB IM=MAXA(IN)+IB IF(IM.GE.MAXA(IN+1))GO TO 40 FORCE=FORCE+V(IN)*A(IM) IF(KEY.EQ.7)WRITE(LW,20)IN,V(IN),IM,A(IM),FORCE 40 IN=IN+1 50 CONTINUE C 55 WRITE(LW,60)NOD,IDIRN,IEQN,DISP,FORCE 60 FORMAT(///,' NODE NO.=',I6,' DIRECTION=',I2,' EQN. NO.=',I6,/, *10X,'PRESCRIBED DISPLACEMENT=',E11.4,' RESULTING FORCE=',E11.4) WRITE(LW,999) 70 CONTINUE C C E N D P R E S C R I B E D D I S P L. S L O O P C RETURN END C*********************************************************************** SUBROUTINE DISP(ID,U,V,NDS,NEQ,NEQF,LTYPE,UKEEP) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. RECOVERS NODAL DISPL.S STORED IN V(NEQ) BY COLSOL. * C * AND ASSIGNS NODE DISPL.S TO ARRAY U(NDS,3). * C * 2. RECOVERS THERMAL-MOISTURE STRAIN (UKEEP). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,3),U(NDS,3),V(NEQ) C C RECOVER TOTAL STRAIN NORMAL TO ELEMENT H-Z PLANE C DUE TO A CHANGE IN THERMAL-MOISTURE CONDITION C IF(LTYPE.EQ.2)UKEEP=V(NEQF) IF(LTYPE.EQ.2)WRITE(LW,1)UKEEP 1 FORMAT(1H1,' FOR CHANGES IN THERMAL-MOISTURE CONDITIONS',/, *' TOTAL STRAIN NORMAL TO H-Z ELEMENT PLANE=',E14.7) C C RECOVER NODAL DISPLACEMENTS C IEQN=1 DO 10 I=1,NDS DO 10 J=1,3 U(I,J)=0.0 IF(IEQN.EQ.NEQF)IEQN=IEQN+1 IF(ID(I,J).EQ.0)GO TO 10 U(I,J)=V(IEQN) IEQN=IEQN+1 10 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE STR(U,EN,EH,EZ,EHZ,ENZ,ENH,NE,NDS,LTYPE,SLOAD,UKEEP) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATES 1. STR AINS FROM DISPLACEMENTS * C * 2. STR ESSESES FROM STRAINS * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION U(NDS,3),EN(NE),EH(NE),EZ(NE),EHZ(NE),ENZ(NE),ENH(NE) C C R E W I N D S T O R E D E L E M E N T D A T A C C TAPE 1 ND1,ND2,ND3 ....... NODE NO.S ASSIGNED TO EACH ELEMENT REWIND 1 C C TAPE 3 AA,BB,...GG,ARR ... ELEMENT GEOMETRIC PARAMETERS REWIND 3 C C TAPE 4 D11,D12,...,D66 ... ELEMENT STIFFNESSES REWIND 4 C C TAPE 7 ENO,EHO,ENHO,EZO .. ELEMENT THERMAL-MOISTURE STRAINS REWIND 7 C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10) 10 FORMAT(1H1,//,T10,14H** STRESSES **,//,T2,8HELE. NO.,T15,4HSIGN, *T30,4HSIGH,T45,4HSIGZ,T60,5HTAUHZ,T75,5HTAUNZ,T90,5HTAUNH,//) C C L O O P T H R U E L E M E N T S C DO 40 K=1,NE C C READ ELEMENT DATA READ(1)ND1,ND2,ND3 READ(3)AA,BB,CC,DD,EE,GG,ARR READ(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66 C C CALCULATE TOTAL STRAINS FROM DISPLACEMENTS C EN(K)=SLOAD IF(LTYPE.EQ.2)EN(K)=UKEEP EH(K)=(AA*U(ND1,2)+CC*U(ND2,2)+EE*U(ND3,2))/ARR EZ(K)=(BB*U(ND1,3)+DD*U(ND2,3)+GG*U(ND3,3))/ARR EHZ(K)=(BB*U(ND1,2)+DD*U(ND2,2)+GG*U(ND3,2)+AA*U(ND1,3)+ 1CC*U(ND2,3)+EE*U(ND3,3))/ARR ENZ(K)=(BB*U(ND1,1)+DD*U(ND2,1)+GG*U(ND3,1))/ARR ENH(K)=(AA*U(ND1,1)+CC*U(ND2,1)+EE*U(ND3,1))/ARR C C IF LOAD CASE IS A CHANGE IN THERMAL-MOISTURE CONDITION C THEN RECOVER ELASTIC STRAINS FROM TOTAL STRAINS C OTHERWISE LOOP TO STATEMENT 20 C IF(LTYPE.LE.1)GO TO 20 C C READ THERMAL-MOISTURE STRAINS C READ(7)ENO,EHO,ENHO,EZO C C RECOVER MECHANICAL (ELASTIC) STRAINS FROM TOTAL STRAINS C EN(K)=EN(K)-ENO EH(K)=EH(K)-EHO ENH(K)=ENH(K)-ENHO EZ(K)=EZ(K)-EZO C C CALCULATE STRESSES FROM ELASTIC STRAINS C 20 SIGN=D11*EN(K)+D12*EH(K)+D13*EZ(K)+D16*ENH(K) SIGH=D12*EN(K)+D22*EH(K)+D23*EZ(K)+D26*ENH(K) SIGZ=D13*EN(K)+D23*EH(K)+D33*EZ(K)+D36*ENH(K) SIGHZ=D44*EHZ(K)+D45*ENZ(K) SIGNZ=D45*EHZ(K)+D55*ENZ(K) SIGNH=D16*EN(K)+D26*EH(K)+D36*EZ(K)+D66*ENH(K) C C WRITE STRESSES ON OUTPUT FILE C WRITE(LW,30)K,SIGN,SIGH,SIGZ,SIGHZ,SIGNZ,SIGNH 30 FORMAT(T2,I5,6(1X,E14.7)) write(10,31)K,SIGN,SIGH,SIGZ,SIGHZ,SIGNZ,SIGNH 31 format(I5,6(1x,e12.5)) C 40 CONTINUE C C E N D E L E M E N T L O O P C C WRITE STRAINS AND DISPLACEMENTS ON OUTPUT FILE C IF((KEY.EQ.5).OR.(KEY.EQ.7))GO TO 45 RETURN C C STRAINS C 45 WRITE(LW,50) 50 FORMAT(1H1,//,T10,13H** STRAINS **,//,T2,8HELE. NO.,T15,2HEN,T30, *2HEH,T45,2HEZ,T60,3HEHZ,T75,3HENZ,T90,3HENH,//) C DO 60 K=1,NE 60 WRITE(LW,70)K,EN(K),EH(K),EZ(K),EHZ(K),ENZ(K),ENH(K) 70 FORMAT(T2,I5,6(1X,E14.7)) C C DISPLACEMENTS C WRITE(LW,80) 80 FORMAT(1H1,//,16X,29H NODAL DISPLACEMENTS,/,10X,65HNODE * NO. U-DISPLACEMENT V-DISPLACEMENT W-DISPLACEMENT) C DO 90 I=1,NDS 90 WRITE(LW,100)I,U(I,1),U(I,2),U(I,3) 100 FORMAT(13X,I5,3(7X,E12.5)) RETURN END