C C TUBEPARM C C Fortran program for parametric studies of unidirectional C and laminated tubes subjected to axial, torsional, C thermal or pressure loading C The theory behind the program can be found in Chapter 10 C of Mechanics of Fibrous Composites by Carl T. Herakovich C (John Wiley and Sons, Inc, 1998, ISBN 0-471-10636-4) C Print command corrected Feb. 15, 1998 by CTH C C**THIS PROGRAM WORKS ON PC'S AS IS WITH FILE INPUT=5 AND C OUTPUT=6 C ELASTICITY SOLUTION FOR AN ANGLE-PLY COMPOSITE TUBE W/THERMAL C AND MECHANICAL LOADS C C Modified in late '96 or '97 to calculate effective C Poisson's Ratio nu,xtheta (or change in radius) C Modified by CTH, Jan. 1996, to simplify naming of output files C Revisions by CTH to calculate tube moduli and coupling ratios C Average stresses added jan. 12, 1995 C New Year's Day January 1, 1995 C Revision by CTH to permit parametric study c r'cvd from Deidre by Herakovich on March 3, 1990 c Modified for SUN by CTH 3/4/90 C original program by CARL ROUSSEAU AND MIKE HYER C MODIFIED BY DEIDRE HIRSCHFELD TO INCLUDE ISOTROPIC LAYERS C OR TRANSVERSE ISOTROPIC SYMMETRY C (CB22=CB33 & CB12=CB13) C C C IMPLICIT REAL*8 (A-H,M-Z) COMMON /IN/ THETA(20), T(20), CB11(20), CB12(20), CB13(20), $ CB16(20),CB22(20), CB23(20), CB26(20), CB33(20), CB36(20), $ CB66(20), RO(20), RI(20), ALPAX(20), ALPA0(20), ALPAR(20), $ ALPAX0(20), IFLAG(20),DT, P, J ,INCR ,IPRNT COMMON /MATL/ E1(20), E2(20), E3(20), G12(20), NU23(20), NU32(20), $ NU13(20), NU31(20), NU12(20), NU21(20), ALPA1(20), ALPA2(20) $ , ALPA3(20), C11(20), C12(20), C13(20), C22(20), C23(20), $ C33(20), C66(20),TW,PIN,POUT COMMON /BOUND/ BC(42,42), ET(42), V(20), V1(20), V2(20), $ Z1(20), Z2(20), Z3(20),SIGB(20) COMMON/OUT/ R(200),EPSX(11,20),EPS0(11,20), $ EPSR(11,20),GAMX0(11,20), $ SIGX(11,20),SIG0(11,20),SIGR(11,20),TAUX0(11,20), $ EPS1(11,20),EPS2(11,20),EPS3(11,20),GAM12(11,20), $ SIG1(11,20),SIG2(11,20),SIG3(11,20),TAU12(11,20), $ RR(11,20),TS(20),F1S1(20), $ F2S2(20),F2S3(20),F11S1(20), $ F22S2(20),F22S3(20),F66S6(20),F12S12(20),F12S13(20), $ F23S23(20) DIMENSION V3(20), V4(20), ZZ(20), Z4(20), Z5(20), $ Z6(20),ANGLE(20) C Basic commands for modifying a fortran program for the SUN C CHARACTER*4 ITITLE(18) CHARACTER*22 IFNAME,output,tubeplt,sigplt,epsplt CHARACTER*8 PROBLM C Read the names of the input and standard lead-in for C output files from the screen. C WRITE (*,3) 3 FORMAT(/'Enter your input data file name (default=input.data):') READ (*,5) IFNAME IF(IFNAME.EQ.' ') IFNAME='input.data' WRITE (*,4) 4 FORMAT(/'Enter lead-in,8char max,output files(default=PROBLM):') READ (*,206) PROBLM 206 FORMAT (A8) do 210 i=8,2,-1 if (PROBLM(i:i).gt.' ') goto 7 210 continue 7 IF(PROBLM.EQ.' ') PROBLM='PROBLM' output = PROBLM(:i)//'.output' sigplt = PROBLM(:i)//'.sigplt' tubeplt = PROBLM(:i)//'.tubeplt' epsplt = PROBLM(:i)//'.epsplt' C displt = PROBLM(:i)//'.displt' 5 FORMAT(A22) C C Open and rewind the input and output data files C OPEN (UNIT=15,FILE=IFNAME,STATUS='OLD') REWIND 15 OPEN (UNIT=16,FILE=output,STATUS='UNKNOWN') REWIND 16 OPEN (UNIT=17,FILE=tubeplt,STATUS='UNKNOWN') REWIND 17 OPEN (UNIT=18,FILE=sigplt,STATUS='UNKNOWN') REWIND 18 OPEN (UNIT=19,FILE=epsplt,STATUS='UNKNOWN') REWIND 19 C OPEN (UNIT=20,FILE=displt,STATUS='UNKNOWN') C REWIND 20 C unit 16 stores all output for reading C unit 18 stores stresses for plotting radially C unit 19 stores strains for plotting radially C future - unit 20 stores displacements for plotting radially C unit 17 stores tube properties for plotting in phi 1111 READ (15,900) ITITLE 900 FORMAT (18A4) WRITE (16,888) 888 FORMAT(////,10x,'***************************************',/, a10x,'*************************************************',/, 110X,'TUBEPARM vH98.0 - LAMINATED COMPOSITE TUBE PROGRAM',/, 210X,'CARL T. HERAKOVICH - UNIVERSITY OF VIRGINIA',/, 310X,'EMAIL: HERAK@VIRGINIA.EDU',/, 410X,'For theory see MECHANICS OF FIBROUS COMPOSITES by ',/ 510X,'CARL T. HERAKOVICH, John Wiley & Sons, Inc., 1998 ',/ 610X,'******************************************************',//) WRITE (16,901) ITITLE C WRITE (20,901) ITITLE 901 FORMAT (//,18A4,/) WRITE(17,999) ITITLE 999 FORMAT (2x,18A4,/2x,'Angle',6x,'eps-0',9x,'gam-0',9x,'Zeta',9x, $'Ex',8x,'Gxt',8x,'W(1,1)/ Ri'/,'*DATA',) WRITE(18,998) ITITLE 998 FORMAT(x,18A4,/,x,'Layr',x,'r/RI',6x,'sigx',9x,'sigh',12x,'sigr', $12x,'tauxh',9x,'sig1',9x,'sig2',9x,'sig3',9x,'sig12',/,'*DATA',) WRITE(19,997) ITITLE 997 FORMAT(x,18A4,/,x,'Layr',x,'r/RI',6x,'epsx',9x,'epsh',12x,'epsr', $12x,'gamxh',9x,'eps1',9x,'eps2',9x,'eps3',9x,'gam12',/,'*DATA',) c PI = 3.1415927 PI = dacos(-1.0d0) READ(15,*) IPRNT, INCR, AMAX, AMIN, AINC, ITYPE WRITE (16,902) IPRNT, ITYPE, INCR, AMAX, AMIN, AINC IF (INCR .GT. 11 ) WRITE (16,930) 930 FORMAT (/' Number of Increments Exceeds Allowable',/) 902 FORMAT (/,5x,'IPRNT = ',I2,4x,'ITYPE = ',I2,/, 15x,'INCR = ',I2,2x,'AMAX = ', 2F6.1,2x,'AMIN = ',F6.1,2x,'AINC = ',F6.1/) C AMAX, AMIN, AINC are maximum, minimum and increment of angle C for parametric results on fiber orientation, ITYPE = 0 for C no parametric study, ITYPE > 0 for parametric study on phi C ITYPE = 2 for unidirectional off-axis helical wound tube C C MATERIAL PROPERTIES READ(15,*) J, RI(1), DT, P,TW,PIN,POUT WRITE(16,60) J, RI(1),DT,P,TW,PIN,POUT 60 FORMAT(/5x,'NUMBER LAYERS = ',I3,5x,'INSIDE RADIUS = ',D10.4, $/5x,'DELTA-TEMPERATURE LOAD =',F5.0,' F'/5x,'AXIAL LOAD=', $D10.4,' LB/'/5x,'TORQUE =',D10.4,'in-lb',/ $5x,'INTERNAL PRESS= ',D10.4,' psi',5x, $' XTERNAL PRESSURE=' ,D10.4,' psi'/) DO 240 I=1,J READ(15,*) E1(I), E2(I), E3(I), G12(I), NU23(I), NU13(I), $ NU12(I), ALPA1(I), ALPA2(I), ALPA3(I), T(I), $ ANGLE(I) IF (I.EQ.1) GO TO 240 RI(I) = RO(I-1) C RI is the inside radius of the layer C RO is the outside radius of the layer C RO(J) is the outside radius of the tube C RI(1) is the inside radius of the tube 240 RO(I) = RI(I) + T(I) C PI = 3.14159265.. PI = dacos(-1.0d0) C avgsigx, avgsigt, avgtauxt are average stresses over C tube thickness, Ny is force per unit length for Pi C loading, Nx and Nxy are CLT loads per unit length avgsigx = P/(PI*(RO(j)**2 - RI(1)**2)) thick = RO(J)-RI(1) avgsigt = PIN*RI(1)/thick Ny = PIN*RI(1) sigtN = Ny/thick C sigtN is normalization factor for sig,theta stresses Nx = avgsigx * thick avgtauxt = 3.0 * TW/(2.*PI*(RO(j)**3 - RI(1)**3)) Nxy = avgtauxt * thick WRITE(16,910) avgsigx,avgsigt,avgtauxt,Nx,Ny,Nxy 910 FORMAT (5x,'avgsigx =',D10.4,' avgsigt =',D10.4, $'avgtauxt =',D10.4,/3x,'Nx =',D10.4,2x,' Ny =',D10.4, $2x,'Nxy =',D10.4/) IF ( IPRNT . EQ . 1 ) WRITE(16,15) 15 FORMAT(////'LAMINA MATERIAL PROPERTIES'////) NANG = ((AMAX - AMIN)/AINC) + 1 IF (ITYPE .EQ. 0 ) NANG = 1 C NANG is the number of angles for parametric study C C THIS PROGRAM AUTOMATICALLY CHECKS FOR E22=E33 AND NU12=NU13 C THEN USES THE APPROPIATE EQUATIONS TO HANDLE TRANSVERSE C ISOTROPY AND ISOTROPIC MATERIALS C C TO SOLVE SIMULTANEOUS EQNS - GASJON BY REDDY IS USED C C END EFFECTS ARE IGNORED C C C ALL INPUT IS FREE FORMAT THUS DATA CAN BE SEPERATED BY COMMAS C OR AT LEAST 1 BLANK SPACE C ***************** C * FIRST CARD * IPRNT - =1 PRINT C AND CBAR MATRICES OF EACH C * OR * LAYER C * FIRST SEVERAL* =0 DO NOT PRINT C AND CBAR C * AS IN SAMPLE * C * INPUT * INCR - NUMBER OF LOCATIONS IN WHICH THE C ***************** STRAINS AND STRESSES ARE CALCULATED C WITHIN EACH LAYER C C J = NUMBER OF LAYERS - MUST BE AT LEAST TWO TO C SATISFY NO. EQNS - IF YOU HAVE ONE MATERIAL C SPLIT IT INTO TWO LAYERS C RI(1) = INSIDE RADIUS OF TUBE C DT = CHANGE IN TEMPERATURE C P = AXIAL LOAD CAN BE EITHER POS(TENSION) OR NEG C TW = TORQUE APPLIED TO TUBE END C PI = INTERNAL PRESSURE C PO = EXTERNAL PRESSURE C NOW THE MAT'L PROPERTIES FOR EACH LAYER BEGINNING AT THE C INSIDE RADIUS, also layer thickness and fiber orientation C E1,E2,E3,G12,NU23,NU13,NU12,ALPHA1,ALPHA2, C ALPHA3, THICKNESS OF LAYER, FIBER ANGLE (IF ISOTROPIC MAT'L C = 0 ) C C REPEAT FOR EACH LAYER " C C PROGRAM ALSO HAS AN OPTION TO DETERMINE TSAI-HILL OR C TSAI-WU FAILURE CRITERIA FOR A ONE MATERIAL TUBE C INPUT CARD ONE - FLAG = 0 - DO NOT DET. FAILURE CRIT. C = 1 - TSAI-HILL INPUT XF,YF,SF C = 2 - TSAI-WU INPUT F1,F2,F11,F22,F66 C C DO 105 II = 1, NANG IF (ITYPE .EQ. 0 ) GO TO 300 C ITYPE = 0, NO PARAMETRIC STUDY ANGLE(1) = AMIN + AINC * (II-1) ANGLE(2) = -AMIN - AINC * (II-1) C ITYPE = 2 for uni off-axis tube (2 layers req'd) IF (ITYPE .EQ. 2) ANGLE(2) = ANGLE(1) 300 CALL INPUT (ANGLE) MCONST=2*J+2 DO 6 K=1,MCONST DO 6 KK=1,MCONST BC(K,KK)=0. 6 CONTINUE C DO 10 I=1,J C IF(IFLAG(I).EQ.1) THEN SIGB(I) =((CB23(I)-CB22(I))*(ALPA0(I)-ALPAR(I))+ & (CB36(I)-CB26(I))*ALPAX0(I))*DT V(I) = 1. V1(I) = 1. V2(I) = 1. V3(I) = 2. + V(I) V4(I) = 2. - V(I) ZZ(I) = 1. Z1(I) = CB12(I) - CB13(I) Z2(I) = 0. Z3(I) = CB23(I)-CB33(I) Z4(I) = CB12(I) + CB13(I) Z5(I) = CB11(I)*ALPAX(I)+Z4(I)*ALPAR(I) Z6(I) = CB23(I) + CB33(I) ELSE SIGB(I) = ((CB23(I) - CB22(I))*ALPA0(I) + (CB33(I) - CB23(I))* $ ALPAR(I) + (CB13(I) - CB12(I))*ALPAX(I) + (CB36(I) - $ CB26(I))*ALPAX0(I)) * DT V(I) = DSQRT (CB22(I) / CB33(I)) V1(I) = 1. + V(I) V2(I) = 1. - V(I) V3(I) = 2. + V(I) V4(I) = 2. - V(I) ZZ(I) = CB33(I) - CB22(I) Z1(I) = (CB12(I) - CB13(I)) / ZZ(I) Z2(I) = (CB26(I) - 2*CB36(I)) / (3*CB33(I) + ZZ(I)) Z3(I) = SIGB(I) / ZZ(I) Z4(I) = CB12(I) + CB13(I) Z5(I) = CB26(I) + CB36(I) Z6(I) = CB23(I) + CB33(I) END IF 10 CONTINUE C C SET UP OF BC'S AT INNER RADIUS STRESS=0 IF(IFLAG(1).EQ.1) THEN BC(1,1) = CB13(1) BC(1,2) = 0. BC(1,3) = Z6(1) BC(1,4) = Z3(1)/ RI(1)**2 ET(1) = (CB13(1)*ALPAX(1)+CB23(1)*ALPA0(1)+CB33(1)*ALPAR(1)+ & CB36(1)*ALPAX0(1))*DT-(Z6(1)*DLOG(RI(1))+CB33(1))* & SIGB(1)/2./CB22(1)-PIN ELSE BC(1,1) = CB13(1) + Z6(1)*Z1(1) BC(1,2) = ((Z6(1) + CB33(1))*Z2(1) + CB36(1)) * RI(1) BC(1,3) = (CB23(1) + V(1)*CB33(1)) * RI(1)**(-V2(1)) BC(1,4) = (CB23(1) - V(1)*CB33(1)) * RI(1)**(-V1(1)) ET(1) = (CB13(1)*ALPAX(1) + CB23(1)*ALPA0(1) + CB33(1)*ALPAR(1) $ + CB36(1)*ALPAX0(1))*DT - Z6(1)*Z3(1)-PIN END IF DO 12 I1=5,2*J+2 BC(1,I1) = 0. 12 BC(2,(I1-2)) = 0. IF(IFLAG(J).EQ.1) THEN C SET UP OF BC'S FOR STRESS =0 AT OUTSIDE R BC(2,(2*J+1)) = Z6(J) BC(2,(2*J+2)) = Z3(J)/ RO(J)**2 ET(2)= (CB13(J)*ALPAX(J)+CB23(J)*ALPA0(J)+CB33(J)*ALPAR(J)+ & CB36(J)*ALPAX0(J))*DT-(Z6(J)*DLOG(RO(J))+CB33(J))* & SIGB(J)/2./CB22(J)-POUT BC(2,1)=C13(J) BC(2,2)=0. ELSE BC(2,1) = CB13(J) + Z6(J)*Z1(J) BC(2,2) = ((Z6(J) + CB33(J))*Z2(J) + CB36(J)) * RO(J) BC(2,(2*J+1)) = (CB23(J) + V(J)*CB33(J)) * RO(J)**(-V2(J)) BC(2,(2*J+2)) = (CB23(J) - V(J)*CB33(J)) * RO(J)**(-V1(J)) ET(2) = (CB13(J)*ALPAX(J) + CB23(J)*ALPA0(J) + CB33(J)*ALPAR(J) $ + CB36(J)*ALPAX0(J))*DT - Z6(J)*Z3(J)-POUT END IF BC(3,1) = 0. BC(3,2) = 0. BC(4,1) = 0. BC(4,2) = 0. ET(3) = P/(2.*PI) ET(4) = TW/(2.*PI) DO 30 I=1,J IF(IFLAG(I).EQ.1) THEN C SET UP OF BC'S FOR APPLIED LOADS NX AND TWIST C LOOP FOR ISOTROPIC LAYER BC(3,1) = BC(3,1) + CB11(I) * (RO(I)**2 - RI(I)**2)/2. BC(3,(1+2*I)) = Z4(I)*(RO(I)**2 - RI(I)**2 )/2. BC(3,(2+2*I)) = Z1(I)*(1./RO(I)- 1./RI(I)) ET(3) = ET(3) + (CB11(I)*ALPAX(I)+CB12(I)*ALPA0(I)+CB13(I)* & ALPAR(I)+CB16(I)*ALPAX0(I))*DT*(RO(I)**2-RI(I)**2)/2. & -(Z4(I)*(RO(I)**2*(DLOG(RO(I))-.5)/2.-RI(I)**2* & (DLOG(RI(I))-.5)/2.)+CB13(I)*(RO(I)**2-RI(I)**2)/2.)* & SIGB(I)/2./CB22(I) BC(4,2) = BC(4,2) + (CB66(I)*(RO(I)**4 -RI(I)**4)/4.) BC(4,(1+2*I)) = 0. BC(4,(2+2*I)) = 0. ET(4) = ET(4) + (CB16(I)*ALPAX(I)+CB26(I)*ALPA0(I)+ & CB36(I)*ALPAR(I)+CB66(I)*ALPAX0(I))*DT*(RO(I)**3- & RI(I)**3)/3.-((CB26(I)-CB36(I))*(RO(I)**3*(DLOG(RO(I))- & 1./3.)/3.-RI(I)**3*(DLOG(RI(I))-1./3.)/3.)+CB36(I)* & (RO(I)**3-RI(I)**3)/3.)*SIGB(I)/2./CB22(I) ELSE BC(3,1) = BC(3,1) + (CB11(I) + Z1(I)*Z4(I)) * (RO(I)**2 - $ RI(I)**2) / 2. BC(3,2) = BC(3,2) + ((CB13(I) + Z4(I)) * Z2(I) + CB16(I)) $ * (RO(I)**3 - RI(I)**3) / 3. BC(3,(1+2*I)) = (CB12(I) + V(I) * CB13(I)) * (RO(I)**V1(I) - $ RI(I)**V1(I)) / V1(I) BC(3,(2+2*I)) = (CB12(I) - V(I) * CB13(I)) * (RO(I)**V2(I) - $ RI(I)**V2(I)) / V2(I) ET(3) = ET(3) + (DT*(CB11(I)*ALPAX(I) + CB12(I)*ALPA0(I) + $ CB13(I)*ALPAR(I) + CB16(I)*ALPAX0(I)) - Z3(I)*Z4(I))* $ (RO(I)**2 - RI(I)**2) / 2. BC(4,1) = BC(4,1) + (CB16(I) + Z1(I)*Z5(I)) * (RO(I)**3 - $ RI(I)**3) / 3. BC(4,2) = BC(4,2) + (CB66(I) + Z2(I) * (Z5(I) +CB36(I)))* $ (RO(I)**4 - RI(I)**4) / 4. BC(4,(1+2*I)) = (CB26(I) + V(I) * CB36(I)) * (RO(I)**V3(I) - $ RI(I)**V3(I)) / V3(I) BC(4,(2+2*I)) = (CB26(I) - V(I) * CB36(I)) * (RO(I)**V4(I) - $ RI(I)**V4(I)) / V4(I) ET(4) = ET(4) + (DT*(CB16(I)*ALPAX(I) + CB26(I)*ALPA0(I) + $ CB36(I)*ALPAR(I) + CB66(I)*ALPAX0(I)) - Z3(I)*Z5(I))* $ (RO(I)**3 - RI(I)**3) / 3. END IF 30 CONTINUE C LOOPS TO SET BC'S AT INTERFACES BETWEEN LAYERS C FIRST MATCH DISPLACEMENT W THEN STRESS SIGR DO 34 KK=5,12 DO 34 LL=1,2*J+2 34 BC(KK,LL) = 0. DO 40 K=1,J-1 IF ((IFLAG(K).EQ.1).AND.(IFLAG(K+1).EQ.0)) THEN BC((4+K),1) = - Z1(K+1) * RO(K) BC((4+K),2) = - Z2(K+1) * RO(K)**2 BC((4+K),(1+2*K)) = RO(K) BC((4+K),(2+2*K)) = RO(K)**(-1) BC((4+K),(3+2*K)) = -RO(K)**V(K+1) BC((4+K),(4+2*K)) = -RO(K)**(-V(K+1)) ET(4+K) = Z3(K+1) * RO(K)-SIGB(K)*RO(K)*DLOG(RO(K))/2./CB22(K) BC((3+J+K),1) = CB13(K)-CB13(K+1)-Z6(K+1)*Z1(K+1) BC((3+J+K),2) = -((Z6(K+1)+CB33(K+1))*Z2(K+1)+CB36(K+1)) * RO(K) BC((3+J+K),(1+2*K)) = Z6(K) BC((3+J+K),(2+2*K)) = Z3(K)/RO(K)**2 BC((3+J+K),(3+2*K)) = -(CB23(K+1) + V(K+1)*CB33(K+1)) * $ RO(K)**(-V2(K+1)) BC((3+J+K),(4+2*K)) = -(CB23(K+1) - V(K+1)*CB33(K+1)) * $ RO(K)**(-V1(K+1)) ET(3+J+K) = DT *(CB13(K)*ALPAX(K) - CB13(K+1)*ALPAX(K+1)+ $ CB23(K)*ALPA0(K) - CB23(K+1)*ALPA0(K+1)+ $ CB33(K)*ALPAR(K) - CB33(K+1)*ALPAR(K+1)+ $ CB36(K)*ALPAX0(K) - CB36(K+1)*ALPAX0(K+1))+ & Z6(K+1)*Z3(K+1)-(Z6(K)*DLOG(RO(K))/2./CB22(K)+ & CB33(K)/2./CB22(K))*SIGB(K) END IF IF((IFLAG(K).EQ.0).AND.(IFLAG(K+1).EQ.1)) THEN BC((4+K),1) = Z1(K) * RO(K) BC((4+K),2) = Z2(K) * RO(K)**2 BC((4+K),(1+2*K)) = RO(K)**V(K) BC((4+K),(2+2*K)) = RO(K)**(-V(K)) BC((4+K),(3+2*K)) = - RO(K) BC((4+K),(4+2*K)) = -1./RO(K) ET(4+K) =SIGB(K+1)*RO(K)*DLOG(RO(K))/2./CB22(K+1) - Z3(K) * RO(K) BC((3+J+K),1) = CB13(K)+Z6(K)*Z1(K)-CB13(K+1) BC((3+J+K),2) = ((Z6(K) + CB33(K))*Z2(K) + CB36(K))*RO(K) BC((3+J+K),(1+2*K)) = (CB23(K) + V(K) *CB33(K) ) * $ RO(K)**(-V2(K) ) BC((3+J+K),(2+2*K)) = (CB23(K ) - V(K )*CB33(K )) * $ RO(K)**(-V1(K )) BC((3+J+K),(3+2*K)) =- Z6(K+1) BC((3+J+K),(4+2*K)) = -Z3(K+1)/RO(K)**2 ET(3+J+K) = DT *(CB13(K)*ALPAX(K) - CB13(K+1)*ALPAX(K+1)+ $ CB23(K)*ALPA0(K) - CB23(K+1)*ALPA0(K+1)+ $ CB33(K)*ALPAR(K) - CB33(K+1)*ALPAR(K+1)+ $ CB36(K)*ALPAX0(K) - CB36(K+1)*ALPAX0(K+1))+ & (Z6(K+1)*DLOG(RO(K))/2./CB22(K+1)+ & CB33(K+1)/2./CB22(K+1))*SIGB(K+1)-Z6(K)*Z3(K) END IF IF((IFLAG(K).EQ.1).AND.(IFLAG(K+1).EQ.1)) THEN BC((4+K),1)=0. BC((4+K),2)=0. BC((4+K),(1+2*K))=RO(K) BC((4+K),(2+2*K))=RO(K)**(-1) BC((4+K),(3+2*K))=-RO(K) BC((4+K),(4+2*K))=-RO(K)**(-1) ET(4+K)=(SIGB(K+1)/CB22(K+1)-SIGB(K)/CB22(K))*RO(K)*DLOG(RO(K))/2. ET(3+J+K) = DT *(CB13(K)*ALPAX(K) - CB13(K+1)*ALPAX(K+1)+ $ CB23(K)*ALPA0(K) - CB23(K+1)*ALPA0(K+1)+ $ CB33(K)*ALPAR(K) - CB33(K+1)*ALPAR(K+1)+ $ CB36(K)*ALPAX0(K) - CB36(K+1)*ALPAX0(K+1))+ & (Z6(K+1)*DLOG(RO(K))/2./CB22(K+1)+ & CB33(K+1)/2./CB22(K+1))*SIGB(K+1)-(Z6(K)* & DLOG(RO(K))/2./CB22(K)+CB33(K)/2./CB22(K))* & SIGB(K) BC((3+J+K),1)=CB13(K)-CB13(K+1) BC((3+J+K),2)=0. BC((3+J+K),(1+2*K))=Z6(K) BC((3+J+K),(2+2*K))=Z3(K)/RO(K)**2 BC((3+J+K),(3+2*K))=-Z6(K+1) BC((3+J+K),(4+2*K))=-Z3(K+1)/RO(K)**2 END IF IF((IFLAG(K).EQ.0).AND.(IFLAG(K+1).EQ.0)) THEN BC((4+K),1) = (Z1(K) - Z1(K+1)) * RO(K) BC((4+K),2) = (Z2(K) - Z2(K+1)) * RO(K)**2 BC((4+K),(1+2*K)) = RO(K)**V(K) BC((4+K),(2+2*K)) = RO(K)**(-V(K)) BC((4+K),(3+2*K)) = -RO(K)**V(K+1) BC((4+K),(4+2*K)) = -RO(K)**(-V(K+1)) ET(4+K) = (Z3(K+1) - Z3(K)) * RO(K) BC((3+J+K),1) = CB13(K)+Z6(K)*Z1(K)-CB13(K+1)-Z6(K+1)*Z1(K+1) BC((3+J+K),2) = ((Z6(K) + CB33(K))*Z2(K) + CB36(K) - (Z6(K+1) $ + CB33(K+1))*Z2(K+1) - CB36(K+1)) * RO(K) BC((3+J+K),(1+2*K)) = (CB23(K) + V(K)*CB33(K)) * RO(K)**(-V2(K)) BC((3+J+K),(2+2*K)) = (CB23(K) - V(K)*CB33(K)) * RO(K)**(-V1(K)) BC((3+J+K),(3+2*K)) = -(CB23(K+1) + V(K+1)*CB33(K+1)) * $ RO(K)**(-V2(K+1)) BC((3+J+K),(4+2*K)) = -(CB23(K+1) - V(K+1)*CB33(K+1)) * $ RO(K)**(-V1(K+1)) ET(3+J+K) = DT * (CB13(K)*ALPAX(K) - CB13(K+1)*ALPAX(K+1) + $ CB23(K)*ALPA0(K) - CB23(K+1)*ALPA0(K+1) + $ CB33(K)*ALPAR(K) - CB33(K+1)*ALPAR(K+1) + $ CB36(K)*ALPAX0(K) - CB36(K+1)*ALPAX0(K+1)) - $ Z6(K)*Z3(K) + Z6(K+1)*Z3(K+1) END IF 40 CONTINUE C WRITE(16,19) 19 FORMAT(///) DO 52 I=1,J WRITE(16,20) I, V(I) 20 FORMAT( 'LAMBDA(',I2,') =',D16.8/) 52 CONTINUE IF(IPRNT.EQ.1) THEN WRITE(16,22) 22 FORMAT(//'BC(I,J) EPS-O GAM-0 A(1 THRU 2J)'/) DO 35 K9=1,2*J+2 WRITE(16,45) (BC(K9,L9),L9=1,(2*J+2)) 45 FORMAT(6D16.8) 35 CONTINUE END IF WRITE(16,32) 32 FORMAT(//'ET(I)'//) WRITE(16,42) (ET(K8),K8=1,(2*J+2)) 42 FORMAT(D16.8) C C SOLVE BC*X=ET USING GASJON (REDDY) C CALL GASJON (2*J+2,BC,42,ET,42) C CALL LEQT1F(BC,1,2*J+2,42,ET,2,WKAREA,IER) C CALL PUTOUT(ITITLE,avgsigx,avgsigt,avgtauxt,Nx,Ny, $Nxy,sigtN) C IF (ITYPE .GT. 0 ) GO TO 100 CALL FAIL (J,INCR) C 100 CONTINUE 105 CONTINUE STOP END C C C SUBROUTINE INPUT (ANGLE) C IMPLICIT REAL*8 (A-H,M-Z) COMMON /IN/ THETA(20), T(20), CB11(20), CB12(20), CB13(20), $ CB16(20),CB22(20), CB23(20), CB26(20), CB33(20), CB36(20), $ CB66(20), RO(20), RI(20), ALPAX(20), ALPA0(20), ALPAR(20), $ ALPAX0(20), IFLAG(20),DT, P, J ,INCR,IPRNT COMMON /MATL/ E1(20), E2(20), E3(20), G12(20), NU23(20), NU32(20), $ NU13(20), NU31(20), NU12(20), NU21(20), ALPA1(20), ALPA2(20) $ , ALPA3(20), C11(20), C12(20), C13(20), C22(20), C23(20), $ C33(20), C66(20),TW,PIN,POUT DIMENSION ANGLE(20) c PI = 3.14159265 PI = dacos(-1.0d0) DO 20 I = 1, J THETA(I) = ANGLE(I) * PI / 180. IF (IPRNT.EQ.1) THEN WRITE(16,30) I,E1(I), E2(I), E3(I), G12(I), NU23(I), NU13(I), $ NU12(I), ALPA1(I), ALPA2(I), ALPA3(I), ANGLE(I), $ RI(I), RO(I), T(I) END IF 30 FORMAT(///'LAYER', I3///7X,'E1',14X,'E2',14X,'E3',13X,'G12',10X, $ 'NU23',6X,'NU13',6X,'NU12'/4D16.8,3F10.5//5X,'ALPHA 1', $ 10X,'ALPHA 2',10X,'ALPHA 3'/3D16.8//'PHI',10X,'RI',10X, $ 'RO',10X,'T'/F4.0,4X,3D12.5) C IF((E1(I).EQ.E2(I)).AND.(E1(I).EQ.E3(I))) THEN ANU=NU12(I) NU32(I)=ANU NU21(I)=ANU NU31(I)=ANU LAM=ANU*E1(I)/((1.+ANU)*(1.-2.*ANU)) IFLAG(I)=1 C11(I)=LAM+2.*G12(I) C12(I)=LAM C13(I)=LAM C22(I)=C11(I) C33(I)=C11(I) C23(I)=LAM C66(I)=G12(I) CB11(I)=C11(I) CB12(I)=C12(I) CB13(I)=C13(I) CB23(I)=C23(I) CB33(I)=C33(I) CB22(I)=C22(I) CB16(I)=0. CB26(I)=0. CB36(I)=0. CB66(I)=C66(I) ALPAX(I)=ALPA1(I) ALPA0(I)=ALPA1(I) ALPAR(I)=ALPA1(I) ALPAX0(I)=0. ELSE IFLAG(I)=0 NU32(I) = NU23(I) * E3(I) / E2(I) NU31(I) = NU13(I) * E3(I) / E1(I) NU21(I) = NU12(I) * E2(I) / E1(I) Z1 = (1. - NU23(I)*NU32(I) - NU13(I)*NU31(I) - NU12(I)*NU21(I) $ - 2*NU32(I)*NU13(I)*NU21(I)) / (E1(I)*E2(I)*E3(I)) C11(I) = (1. - NU23(I)*NU32(I)) / (E2(I)*E3(I)*Z1) C12(I) = (NU12(I) + NU32(I)*NU13(I)) / (E1(I)*E3(I)*Z1) C13(I) = (NU13(I) + NU12(I)*NU23(I)) / (E1(I)*E2(I)*Z1) C22(I) = (1. - NU13(I)*NU31(I)) / (E1(I)*E3(I)*Z1) C23(I) = (NU23(I) + NU21(I)*NU13(I)) / (E1(I)*E2(I)*Z1) C33(I) = (1. - NU12(I)*NU21(I)) / (E1(I)*E2(I)*Z1) C66(I) = G12(I) TAYTA = THETA(I) M = DCOS (TAYTA) N = DSIN (TAYTA) IF(ANGLE(I).EQ.0.) THEN M=1. N=0. END IF IF(ANGLE(I).EQ.90) THEN M=0. N=1. END IF Z2 = C12(I) + 2*C66(I) CB11(I) = C11(I)*M**4 + 2*M**2*N**2*Z2 + C22(I)*N**4 CB12(I) = (C11(I) + C22(I) - 4*C66(I))*M**2*N**2 + (N**4 + $ M**4)*C12(I) CB13(I) = M**2*C13(I) + N**2*C23(I) CB22(I) = C11(I)*N**4 + 2*M**2*N**2*Z2 + C22(I)*M**4 CB23(I) = C13(I)*N**2 + C23(I)*M**2 CB33(I) = C33(I) CB16(I) = ((C11(I) - Z2)*M**2 + (Z2 - C22(I))*N**2) * M * N CB26(I) = ((C11(I) - Z2)*N**2 + (Z2 - C22(I))*M**2) * M * N CB36(I) = (C13(I) - C23(I)) * M *N CB66(I) = M**2*N**2*(C11(I)-2*C12(I)+C22(I)) + C66(I)*(M*M-N*N)**2 ALPAX(I) = ALPA1(I)*M**2 + ALPA2(I)*N**2 ALPA0(I) = ALPA1(I)*N**2 + ALPA2(I)*M**2 ALPAR(I) = ALPA3(I) ALPAX0(I) = 2*N*M*(ALPA1(I) - ALPA2(I)) IF(ALPA1(I).EQ.ALPA2(I)) ALPAX0(I)=0.0 E2E3=E2(I)-E3(I) CB123=CB12(I)-CB13(I) IF ((DABS(E2E3).LT.1.D-06).AND.(DABS(CB123).LT.1.D-06))IFLAG(I)=1 END IF IF ( IPRNT . EQ . 1 ) THEN WRITE(16,50) C11(I), C12(I), C13(I), C22(I), C23(I), C33(I), $ C66(I), CB11(I), CB12(I),CB13(I), CB16(I), CB22(I), $ CB23(I), CB26(I), CB33(I), CB36(I), CB66(I), $ ALPAX(I), ALPA0(I), ALPAR(I), ALPAX0(I) END IF 50 FORMAT(////'C MATRIX'//3D16.8/16X,2D16.8/32X,D16.8/48X,D16.8/// $ 'CBAR MATRIX'//4D16.8/16X,3D16.8/32X,2D16.8/48X,D16.8/// $ 5X,'ALPHA-X',7X,'ALPHA-THETA',7X,'ALPHA-R' $ ,5X,'ALPHA-X,THETA'//4D16.8) 20 CONTINUE RETURN END C C C SUBROUTINE PUTOUT(ITITLE,avgsigx,avgsigt,avgtauxt, $Nx,Ny,Nxy,sigtN) IMPLICIT REAL*8 (A-H,M-Z) COMMON /IN/ THETA(20), T(20), CB11(20), CB12(20), CB13(20), $ CB16(20),CB22(20), CB23(20), CB26(20), CB33(20), CB36(20), $ CB66(20), RO(20), RI(20), ALPAX(20), ALPA0(20), ALPAR(20), $ ALPAX0(20), IFLAG(20),DT, P, J , INCR,IPRNT COMMON /MATL/ E1(20), E2(20), E3(20), G12(20), NU23(20), NU32(20), $ NU13(20), NU31(20), NU12(20), NU21(20), ALPA1(20), ALPA2(20) $ , ALPA3(20), C11(20), C12(20), C13(20), C22(20), C23(20), $ C33(20), C66(20),TW,PIN,POUT COMMON /BOUND/ BC(42,42), ET(42), V(20), V1(20),V2(20), $ Z1(20), Z2(20), Z3(20),SIGB(20) COMMON/OUT/ R(200),EPSX(11,20), $ EPS0(11,20),EPSR(11,20),GAMX0(11,20), $ SIGX(11,20),SIG0(11,20),SIGR(11,20),TAUX0(11,20), $ EPS1(11,20),EPS2(11,20),EPS3(11,20),GAM12(11,20), $ SIG1(11,20),SIG2(11,20),SIG3(11,20),TAU12(11,20), $ RR(11,20),TS(20),F1S1(20), $ F2S2(20),F2S3(20),F11S1(20), $ F22S2(20),F22S3(20),F66S6(20),F12S12(20),F12S13(20), $ F23S23(20) CHARACTER*4 ITITLE(18) DIMENSION W(11,20) PI = dacos(-1.0d0) Angle = Theta(1) * 180.0 / PI Area = PI * (RO(J)**2 - RI(1)**2) PolarI = PI * (RO(J)**4 - RI(1)**4)/2.0 IF (TW .EQ. 0) THEN Zeta = ET(2)* RI(1)/ET(1) ELSE Zeta = ET(1)/(ET(2)*RI(1)) END IF C Ex .. tube axial modulus C Gxt.. tube effective shear modulus Ex = 0.0 Gxt = 0.0 If (P .NE. 0.0 ) Ex = P / (Area * ET(1) ) If (TW .NE. 0.0 ) Gxt = TW / (PolarI * ET(2)) WRITE(16,10) ( ET(I1),I1=1,(2*J+2) ) 10 FORMAT(//10x,'LOAD CONSTANTS', $/5x,'EPSILON-0 = ',D12.6/5x,'GAMMA-0 = ',D12.6, $//7X,'A1',16X,'A2'/5(2D12.6,/)) DO 20 L=1,J TAYTA = THETA(L) M = DCOS(TAYTA) N = DSIN(TAYTA) T9 = T(L)/(INCR-1) DO 20 K=1,INCR R(K) = RI(L) + (K-1)*T9 EPSX(K,L) = ET(1) GAMX0(K,L) = ET(2)*R(K) IF( IFLAG(L).EQ.1 ) THEN ZK=CB26(L)-2.*CB36(L)/3./CB22(L)*ET(2) EPS0(K,L)=ET(1+2*L)+ET(2+2*L)/R(K)**2+ZK*R(K)+SIGB(L)/2./CB22(L) & *DLOG(R(K)) EPSR(K,L)=ET(1+2*L)-ET(2+2*L)/R(K)**2+2.*ZK*R(K)+SIGB(L)/2./ & CB22(L)*(DLOG(R(K))+1) END IF IF(IFLAG(L).EQ.0) THEN EPS0(K,L) = Z1(L)*ET(1) + Z2(L)*ET(2)*R(K) + Z3(L) + $ ET(1+2*L)*R(K)**(-V2(L)) + ET(2+2*L)* $ R(K)**(-V1(L)) EPSR(K,L) = Z1(L)*ET(1) + 2*Z2(L)*ET(2)*R(K) + Z3(L) + $ ET(1+2*L)*V(L)*R(K)**(-V2(L)) - ET(2+2*L)* $ V(L)*R(K)**(-V1(L)) END IF C W(K,L)=EPS0(K,L)*R(K) EA1 = ET(1) - ALPAX(L)*DT EA2 = EPS0(K,L) - ALPA0(L)*DT EA3 = EPSR(K,L) - ALPAR(L)*DT EA4 = GAMX0(K,L) - ALPAX0(L)*DT SIGX(K,L) = CB11(L)*EA1+CB12(L)*EA2+CB13(L)*EA3+CB16(L)*EA4 SIG0(K,L) = CB12(L)*EA1+CB22(L)*EA2+CB23(L)*EA3+CB26(L)*EA4 SIGR(K,L) = CB13(L)*EA1+CB23(L)*EA2+CB33(L)*EA3+CB36(L)*EA4 TAUX0(K,L) = CB16(L)*EA1+CB26(L)*EA2+CB36(L)*EA3+CB66(L)*EA4 EPS1(K,L) = EPSX(K,L)*M*M + EPS0(K,L)*N*N + GAMX0(K,L)*M*N EPS2(K,L) = EPSX(K,L)*N*N + EPS0(K,L)*M*M - GAMX0(K,L)*M*N EPS3(K,L) = EPSR(K,L) GAM12(K,L) = -EPSX(K,L)*2*M*N + EPS0(K,L)*2*M*N + GAMX0(K,L) $ *(M*M -N*N) SIG1(K,L) = SIGX(K,L)*M*M + SIG0(K,L)*N*N + TAUX0(K,L)*2*M*N SIG2(K,L) = SIGX(K,L)*N*N + SIG0(K,L)*M*M - TAUX0(K,L)*2*M*N SIG3(K,L) = SIGR(K,L) TAU12(K,L) = -SIGX(K,L)*M*N + SIG0(K,L)*M*N + TAUX0(K,L) $ *(M*M -N*N) 20 CONTINUE deltaRi = W(1,1)/RI(1) WRITE(17,900) Angle, ET(1), ET(2), Zeta, Ex, Gxt,deltaRi 900 FORMAT (x,F6.2,x,6(D12.5,x)) WRITE(16,30) 30 FORMAT (//'TOTAL STRAINS IN X-Y SYSTEM...LAYER...R/R...' $ ,'EPS-X...EPS-THETA...EPS-R...GAMMA-X,THETA'//) C C output global strains for through-thickness C DO 40 L=1,J DO 40 K=1,INCR R(K) = RI(L) + (K-1)*T(L)/(INCR-1) RR(K,L) = (R(K)-RI(1))/(RO(J)-RI(1)) WRITE(16,50)L,RR(K,L),EPSX(K,L),EPS0(K,L),EPSR(K,L),GAMX0(K,L) 50 FORMAT(I2,F8.5,4D15.7) 40 CONTINUE C C output global stresses for through-thickness C WRITE(16,60) 60 FORMAT (//'STRESSES IN X-Y SYSTEM...LAYER...R/R...', $ 'SIG-X...SIG-THETA...SIG-R...TAU-X,THETA'//) DO 70 L=1,J DO 70 K=1,INCR WRITE(16,50)L,RR(K,L),SIGX(K,L),SIG0(K,L),SIGR(K,L),TAUX0(K,L) 70 CONTINUE C C output tube properties WRITE(16,901) DO 80 L=1,J DO 80 K=1,INCR WRITE(16,50) L,RR(K,L),EPSX(K,L),GAMX0(K,L),W(K,L) 80 CONTINUE 901 FORMAT(//'DISPLACEMENTS ...LAYER....R/R....U(R,X)...V(R,X)... $ PER UNIT X ... W(R) '//) WRITE(16,90) 90 FORMAT (//'TOTAL STRAINS IN 1-2 SYSTEM...LAYER...R/R...', $ 'EPS-1...EPS-2...EPS-3...GAMMA-1,2'//) DO 100 L=1,J DO 100 K=1,INCR WRITE(16,50)L,RR(K,L),EPS1(K,L),EPS2(K,L),EPS3(K,L),GAM12(K,L) 100 CONTINUE WRITE(16,120) 120 FORMAT (//'STRESSES IN 1-2 SYSTEM...LAYER...R/R...', $ 'SIG-1...SIG-2...SIG-3...TAU-1,2'//) DO 130 L=1,J DO 130 K=1,INCR WRITE(16,50)L,RR(K,L),SIG1(K,L),SIG2(K,L),SIG3(K,L),TAU12(K,L) 130 CONTINUE C C output stresses, strains and displacements to files for plotting C DO 150 L = 1,J DO 150 K = 1,INCR WRITE(18,902)L,RR(K,L),SIGX(K,L),SIG0(K,L),SIGR(K,L), 1TAUX0(K,L),SIG1(K,L),SIG2(K,L),SIG3(K,L),TAU12(K,L) WRITE(19,902)L,RR(K,L),EPSX(K,L),EPS0(K,L),EPSR(K,L), 1GAMX0(K,L),EPS1(K,L),EPS2(K,L),EPS3(K,L),GAM12(K,L) 150 CONTINUE 902 FORMAT(I2,F8.5,8D15.7) RETURN END C C C SUBROUTINE FAIL(J,INCR) C C THIS SUBROUTINE ONLY WORKS FOR LAMINATES OF THE SAME MATERIAL C C READ FLAG 0 = NO FAILURE CRITERIA ; 1 = INPUT STRENGTHS ; C 2 = INPUT FIJ IMPLICIT REAL*8 (A-H,M-Z) COMMON/OUT/ R(200),EPSX(11,20), $ EPS0(11,20),EPSR(11,20),GAMX0(11,20), $ SIGX(11,20),SIG0(11,20),SIGR(11,20),TAUX0(11,20), $ EPS1(11,20),EPS2(11,20),EPS3(11,20),GAM12(11,20), $ SIG1(11,20),SIG2(11,20),SIG3(11,20),TAU12(11,20), $ RR(11,20),TS(20),F1S1(20), $ F2S2(20),F2S3(20),F11S1(20), $ F22S2(20),F22S3(20),F66S6(20),F12S12(20),F12S13(20), $ F23S23(20) C C THIS SUBROUTINE ONLY WORKS FOR LAMINATES OF THE SAME MATERIAL C C READ FLAG 0 = NO FAILURE CRITERIA ; 1 = TSAI-HILL:INPUT STRENGTHS ; C 2 = TSAI-WU:INPUT FI & FIJ C READ(15,*) FL IF (FL.EQ.1) GOTO 10 IF (FL.EQ.2) GOTO 20 GO TO 99 10 READ(15,*) XF,YF,SF F1 = 0. F2 = 0. F11 = 1./XF**2 F22 = 1./YF**2 F66 = 1./SF**2 F12 = -0.5*F11 F23 = -F22-F12 WRITE(16,70) 70 FORMAT(//'TSAI-HILL'//) GOTO 30 20 READ(15,*) F1,F2,F11,F22,F66 F12 = 0. F23 = 0. WRITE(16,45) 45 FORMAT(//'Tensor Polynomial'//) 30 DO 40 I=1,J F1S1(I) = F1 * SIG1(INCR,I) F2S2(I) = F2 * SIG2(INCR,I) F2S3(I) = F2 * SIG3(INCR,I) F11S1(I) = F11 * SIG1(INCR,I)**2 F22S2(I) = F22 * SIG2(INCR,I)**2 F22S3(I) = F22 * SIG3(INCR,I)**2 F66S6(I) = F66 * TAU12(INCR,I)**2 F12S12(I) = F12 * SIG1(INCR,I)*SIG2(INCR,I) F12S13(I) = F12 * SIG1(INCR,I)*SIG3(INCR,I) F23S23(I) = F23 * SIG2(INCR,I)*SIG3(INCR,I) TS(I) = F1S1(I) + F2S2(I) + F2S3(I) + F11S1(I) + F22S2(I) + $ F22S3(I) + F66S6(I) + F12S12(I) + F12S13(I) + F23S23(I) 40 TS(I) = DABS(TS(I)) WRITE(16,50) F1,F2,F11,F22,F66,F12,F23 50 FORMAT(7X,'F1',14X,'F2',13X,'F11',13X,'F22',13X,'F66',13X, $ 'F12',13X,'F23'/7D16.8) WRITE(16,65) 65 FORMAT(//'PLY',6X,'F1S1',12X,'F2S2',12X,'F2S3'/) DO 80 I=1,J 80 WRITE(16,60) I,F1S1(I),F2S2(I),F2S3(I) WRITE(16,110) 110 FORMAT(//'PLY',4X,'F11S1S1',9X,'F22S2S2',9X,'F22S3S3',9X, $ 'F66S6S6',9X,'F12S1S2',9X,'F12S1S3',9X,'F23S2S3'/) DO 90 I=1,J 90 WRITE(16,60) I,F11S1(I),F22S2(I),F22S3(I),F66S6(I),F12S12(I), $ F12S13(I),F23S23(I) WRITE(16,120) 120 FORMAT(//'PLY...FAILURE FRACTION,TS'/) DO 100 I=1,J 100 WRITE(16,60) I,TS(I) 60 FORMAT(I2,7D16.8) 99 RETURN END C C SUBROUTINE GASJON (N,A,NRMAX,B,NBMAX) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NRMAX,NRMAX),B(NBMAX),INDEX(42,2),IPVOT(42) DO 2 J=1,N 2 IPVOT(J)=0. DO 14 I=1,N T=0. DO 5 J=1,N IF (IPVOT(J).EQ.1) GO TO 5 DO 4 K=1,N IF (IPVOT(K)-1) 3,4,17 3 IF (DABS(T).GE.DABS(A(J,K))) GO TO 4 IROW=J ICOL=K T=A(J,K) 4 CONTINUE 5 CONTINUE IPVOT(ICOL)=IPVOT(ICOL)+1 IF (IROW.EQ.ICOL) GO TO 8 DO 6 L=1,N T=A(IROW,L) A(IROW,L)=A(ICOL,L) 6 A(ICOL,L)=T T=B(IROW) B(IROW)=B(ICOL) 7 B(ICOL)=T 8 INDEX(I,1)=IROW INDEX(I,2)=ICOL PIVOT=A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 9 L=1,N 9 A(ICOL,L)=A(ICOL,L)/PIVOT 10 B(ICOL)=B(ICOL)/PIVOT 11 DO 14 LI=1,N IF (LI.EQ.ICOL) GO TO 14 T=A(LI,ICOL) A(LI,ICOL)=0. DO 12 L=1,N 12 A(LI,L)=A(LI,L)-A(ICOL,L)*T 13 B(LI)=B(LI)-B(ICOL)*T 14 CONTINUE 17 RETURN END