C C COMPLOT C C Fortran program for parametric studies of unidirectional C off-axis lamina and general laminates with variable C thickness 0, 90 and +/- theta layers. Output is stored C in files suitable for plotting. C The theory behind the program can be found in the book C Mechanics of Fibrous Composites by Carl T. Herakovich C (John Wiley and Sons, Inc, 1998, ISBN 0-471-10636-4) C C C COMPOSITE PLOTTING PROGRAM C General Laminates of type 0-T1, +-theta-T2, 90-T3 C Modified 8/28/93 for principal stresses in angle-ply C laminates under unit axial average stress (by CTH) C MODIFIED 12-19-92 FOR ALPAZ & NUZX FOR GENERAL LAMINATES C *** modified by CTH 12-9-92 to calculate & store for plotting C *** Ex, Ey, nuxy, Gxy for 0-T1, +-theta-T2, 90-T3 c *** laminates C *** modified 10/5/91 to consolidate programs *** C *** modified 11-11-96 to include output for ratio of C hoop to axial strain in angle-ply laminate tube under C internal pressure APTSR is the varialbe. Only printed C to unit 17 for plotting C C CALCULATES LAMINATE ALPHA Z AND NU XZ FOR ANGLE-PLY LAMINATES C CALCULATES Ex, Ey, NUxy, NUyx, AND Gxy FOR ANGLEPLY LAMINATES C C **** VALUES ARE STORED ON TAPE FOR PLOTTING C CALCULATES S, Q, SBAR, QBAR, E'S, G, NU'S, ETA'S, ALPHA'S C FOR OFF-AXIS LAMINAE C CALCULATES 21 VALUES FOR EACH ANGLE SPECIFIED C PROGRAM WRITTEN BY C. T. HERAKOVICH, APPLIED MECHANICS C UNIVERSITY OF VIRGINIA C CHARLOTTESVILLE, VIRGINIA 22903, PHONE (804) 24-3605 C email address: herak@virginia.edu IMPLICIT REAL*8(A-H,O-Z) CHARACTER*22 IFNAME,output,applt,Qijb,Sijb,EGNus,Etas, 1sig1s,prpplt,plotdat,D3Cbar,D3Sbar,zprpplt,etamis,off1sigs CHARACTER*8 MATL DOUBLE PRECISION DSIN, DCOS DIMENSION S(6,6),C(6,6),GALPX(181),GALPY(181),astar(4,181) DIMENSION Q(4),QB(1086),S2(4),SB(1086),PROP(1629) DIMENSION ITITLE(18),Y2(181),X2(181),Y4(181),X4(181),X(181),Y(181) 1,ANG(181),XX(181),YY(181),AN(181),A(3900),MQ(6),MP(9),APEX(181), 2APEY(181),APNXY(181),APNYX(181),APGXY(181),ALPX(181),ALPY(181), 3ALPXY(181),ALPXB(181),ALPYB(181),SIG1(181),SIG2(181),SIG12(181), 4EPS1(181),EPS2(181),EPS12(181),APTSR(181) DATA MQ/1,4,6,2,3,5/,MP/1,2,4,3,5,6,7,8,9/ C Prompt the user for the names of the input and output data files. C Read the names of the input and output data files from the screen. C Read name of material for use with standard plot output WRITE (*,31) 31 FORMAT(/'Enter your input data file name (default=input.data):') READ (*,51) IFNAME IF(IFNAME.EQ.' ') IFNAME='input.data' C WRITE (*,41) C 41 FORMAT(/'Enter your output data file name (default=output.data):') C READ (*,51) OFNAME 51 FORMAT(A22) C IF(OFNAME.EQ.' ') OFNAME='output.data' WRITE (*,61) 61 FORMAT(/'Enter material name 8 char. max (default=MATL):') READ (*,49) MATL 49 FORMAT (A8) do 6 i=8,2,-1 if (MATL(i:i).gt.' ') goto 7 6 continue 7 IF(MATL.EQ.' ') MATL='MATL' output = MATL(:i)//'.output' applt = MATL(:i)//'.applt' Qijb = MATL(:i)//'.Qijb' Sijb = MATL(:i)//'.Sijb' EGNus = MATL(:i)//'.EGNus' Etas = MATL(:i)//'.Etas' sig1s = MATL(:i)//'.sig1s' prpplt = MATL(:i)//'.prpplt' plotdat = MATL(:i)//'.plotdat' D3Cbar = MATL(:i)//'.D3Cbar' D3Sbar = MATL(:i)//'.D3Sbar' zprpplt = MATL(:i)//'.zprpplt' etamis = MATL(:i)//'.etamis' off1sigs = MATL(:i)//'.off1sigs' C C Open and rewind the input and output data files C OPEN (UNIT=15,FILE=IFNAME,STATUS='OLD',ACCESS='SEQUENTIAL', 1 FORM='FORMATTED') REWIND 15 OPEN (UNIT=16,FILE=output,STATUS='UNKNOWN') REWIND 16 OPEN (UNIT=17,FILE=applt,STATUS='UNKNOWN') REWIND 17 OPEN (UNIT=18,FILE=Qijb,STATUS='UNKNOWN') REWIND 18 OPEN (UNIT=11,FILE=Sijb,STATUS='UNKNOWN') REWIND 11 OPEN (UNIT=12,FILE=EGNus,STATUS='UNKNOWN') REWIND 12 OPEN (UNIT=13,FILE=Etas,STATUS='UNKNOWN') REWIND 13 OPEN (UNIT=14,FILE=sig1s,STATUS='UNKNOWN') REWIND 14 OPEN (UNIT=19,FILE=prpplt,STATUS='UNKNOWN') REWIND 19 OPEN (UNIT=20,FILE=plotdat,STATUS='UNKNOWN') REWIND 20 OPEN (UNIT=21,FILE=D3Cbar,STATUS='UNKNOWN') REWIND 21 OPEN (UNIT=22,FILE=D3Sbar,STATUS='UNKNOWN') REWIND 22 OPEN (UNIT=23,FILE=etamis,STATUS='UNKNOWN') REWIND 23 OPEN (UNIT=24,FILE=off1sigs,STATUS='UNKNOWN') REWIND 24 OPEN (UNIT=10,FILE=zprpplt,STATUS='UNKNOWN') REWIND 10 C C UNIT 15 - INPUT DATA FILE C UNIT 16 - STANDARD OUTPUT FILE C 17 LAMPLOT.APPLT - plot data for angle-plys engr. consts C 18 LAMPLOT.Qijb - plot data for off-axis Qijbar C 11 LAMPLOT.Sijb - plot data for off-axis Sijbar C 12 LAMPLOT.EGNus - pltdat: off-axis E,G,nu, etxy,x C 13 LAMPLOT.Etas - pltdat: off-axis etxy,y x,xy,y,xy C 14 ANGLPLY.SIG1s & EPS1s - sig1, sig2, sig12, eps1, C eps2, & eps12 for angle-ply laminates C 19 LAMPLOT.PRPPLT - off-axis plot data E, G, nu, etas, alphas C 10 LAMPLOT.ZPRPPLT - alpaz and NUxz for general laminates C 20 plotdat:EX, EY, NUXY, GXY,& ALPAs- GENERAL LAMINATE PLOTS C 21 genlam.D3Cbar, 13 off-axis 3-D stiffness coefficients C 22 genlam.D3Sbar, 13 off-axis 3-D compliance coefficients C WRITE (16,888) 888 FORMAT(////,10x,'***************************************',/, a10x,'*************************************************',/, 110X,'COMPLOT vJ98 - PARAMETRIC STUDIES FOR',/, 210X,'3D,OFF-AXIS & GENERAL LAMINATES',/, 210X,'CARL T. HERAKOVICH',/, 310X,'UNIVERSITY OF VIRGINIA',/, 410X,'email: HERAK@VIRGINIA.EDU',/ 510X,'For theory see MECHANICS OF FIBROUS COMPOSITES by ',/ 610X,'CARL T. HERAKOVICH, John Wiley & Sons, Inc., 1998 ',/ 510X,'******************************************************',//) READ (15,900) ITITLE 900 FORMAT (18A4) WRITE (16,901) ITITLE WRITE (10,901) ITITLE WRITE (11,901) ITITLE WRITE (12,901) ITITLE WRITE (13,901) ITITLE WRITE (14,901) ITITLE WRITE (17,901) ITITLE WRITE (18,901) ITITLE WRITE (19,901) ITITLE WRITE (20,901) ITITLE WRITE (21,901) ITITLE WRITE (22,901) ITITLE WRITE (23,901) ITITLE WRITE (24,901) 901 FORMAT (//,18A4,/) 902 FORMAT (6E9.2,/,4I5) 05 READ(15,*,END=1000) E1,E2,G12,GNU12,GNU23, 1AL1,AL2,TZERO,TINC,NOTINC,K1,K2,K3,K4,T1,T2,T3 C if K4 = 0, results for off-axis lamina C if K4 > 0, calculate properties for general laminate C where T1, T2, T3 are layer thicknesses C Laminate OF THE TYPE 0-T1/+-THETA-T2/90-T3 IF (E1.LT.0.0) go to 1000 C C LAMINA PROPERTIES AND PLOT PARAMETERS C K1=0 DON'T PRINT Q'S; K2=0 DON'T PRINT PROPS; C K3 - NOT IN USE AT THIS TIME C TZERO-SMALLEST ANGLE; TINC-INCREMENT; NOTINC- # OF INCREMENTS TMAX = TZERO + TINC * NOTINC WRITE (16,903) E1,E2,G12,GNU12,GNU23, 1AL1,AL2,TZERO,TMAX,NOTINC,TINC,K1,K2,K3,K4,T1,T2,T3 903 FORMAT (10X,'E1 =',E14.6,3X,'E2 = ',E14.6,3X,'G12 = ',E14.6,/10X,' 1NU12 =',E14.6,3X,'NU23 = ',E14.6, 2/10X,'ALPA1 = ',E14.6,3X,'ALPA2 = ',E14.6, 2/10X,'THETA RANGES FROM ',F7.3,' TO ',F10.3,/10X,'IN ',I5,3x,' 3INCREMENTS OF ',F10.3,' DEGREES',/,10X,'PRINT PARAMETERS K1 = 4',I2,3X,'K2 = ',I2,3X,'K3 = ', I2,3X,'K4 = ', I2,/,10X, 5'Layer Thicknesses T1 = ',F8.4,3X,'T2 = ',F8.4,3X,'T3 = ',F8.4) C CALL PROP3D (K,QB,E1,E2,GNU12,GNU23,G12,S,C) C C PROP3D CALCULATES 3-D CIJ AND SIJ FOR UNIDIRECTIONAL MATERIAL C NANGLS = NOTINC + 1 WRITE (24,904) 904 FORMAT (5x,'PRINCIPAL MATERIAL STRESSES - UNIT', 1' LOADING - OFF-AXIS LAMINA',/,3x,'Angle'10x,'SIG1',12x, 2'SGI2',12x,'TAU12',/,'*DATA') DO 17 I = 1,NANGLS ANG(I) = TZERO + TINC*(I-1) C C COMPUTE TRIG FUNCTIONS FOR ALL ANGLES C AN(I) = ANG(I) C AN ARE ANGLES IN DEGS., ANG ARE ANGLES IN RADIANS ANG(I) = ANG(I)* ACOS(-1.0)/180. X(I) = DSIN(ANG(I)) XX(I) = X(I)*X(I) Y(I) = DCOS(ANG(I)) YY(I) = Y(I)*Y(I) X2(I) = DSIN(2.*ANG(I)) X4(I) = DSIN(4.*ANG(I)) Y2(I) = DCOS(2.*ANG(I)) Y4(I) = DCOS(4.*ANG(I)) TWOMN = -X(I)*Y(I) WRITE (24,908) AN(I), YY(I), XX(I),TWOMN 908 FORMAT (3x,F6.1,2x,3(2x,E14.6)) C C CTE's ALPAX, ALPAY, & ALPAXY FOR EACH ANGLE C ALPX(I) = YY(I) * AL1 + XX(I) * AL2 ALPY(I) = XX(I) * AL1 + YY(I) * AL2 ALPXY(I) = 2.0 * X(I) * Y(I) * (AL1 - AL2) 17 CONTINUE C C 2-D STIFFNESS & COMPLIANCE MATRICES & INVARIANTS C GNU21 = E2*GNU12/E1 W = 1. - GNU12*GNU21 Q(1) = E1/W Q(2) = Q(1)*GNU21 Q(3) = E2/W Q(4) = G12 Q12Q = Q(1) + Q(3) Q16Q = 2.*Q(2) + 4.*Q(4) U1 = (3.*Q12Q + Q16Q)/8. U2 = (Q(1)-Q(3))/2. U3 = (Q12Q - Q16Q)/8. U4 = (Q12Q +6.*Q(2)-4.*Q(4))/8. U5 = (Q12Q - 2.*Q(2) + 4.*Q(4))/8. S2(1) = 1./E1 S2(2) = -GNU12/E1 S2(3) = 1./E2 S2(4) = 1./G12 IF (K1.EQ.0) GO TO 25 C C PRINT [Q] and [S] MATRICES C WRITE(16,915) 915 FORMAT (//10X,'Plane Str Q & S, PRINCIPAL MATERIAL COORD.') 916 FORMAT (2(10X,E14.6)/34X,E14.6,/48X,E14.6,/,10X,'S MATRIX'/) WRITE (16,916) (Q(J), J = 1,4) WRITE(16,917) (S2(J), J = 1,4) 917 FORMAT (2(10X,E14.6)/34X,E14.6,/48X,E14.6/) C C QBAR & SBAR FOR EACH PLY (ANGLE) C QB11=QB(L+1) QB(12)=QB(L+2) QB(16)=QB(L+3) C QB22=QB(L+4) QB26=QB(L+5) QB66=QB(L+6) C SAME FORMAT IS TRUE FOR THE SB MATRIX C 25 M = 0 L = 0 J = 0 write (16,930) 930 FORMAT (//,10x,'ANGLE-PLY ENGINEERING CONSTANTS',/, $3x,'ANG',10X,'Ex ',' Ey ',' NUxy ' $,' Gxy ',' ALPAx ',' ALPAy ',/) write (17,931) 931 FORMAT (//,10x,'ANGLE-PLY ENGINEERING CONSTANTS',/, 13x,'ANG',10X,'Ex ',' Ey ',' NUxy ', 2' Gxy ',' ALPAx ',' ALPAy ',/"*DATA",/) WRITE (14,961) 961 FORMAT (/,10x,'Angle-ply Material Principal stresses - 1unit stress',/,3x,'ANGLE',5x,'SIG1',10x,'SIG2',10x,'SIG12', 24x,'EPS1',4x,'EPS2',4x,'EPS12',/'*DATA') 30 DO 100 K = 1, NANGLS QB(L+1) = U1 + U2*Y2(K) + U3*Y4(K) QB(L+2) = U4 - U3*Y4(K) QB(L+3) = (U2*X2(K))/2. + U3*X4(K) IF (AN(K).EQ. 90.0 .OR. AN(K) .EQ. -90.0 ) QB(L+3)= 0.0 PA = S2(1) + S2(3) PB = S2(1) - S2(3) PC = 2.*S2(2) + S2(4) QB(L+4) = U1 - U2*Y2(K) + U3*Y4(K) QB(L+5) = (U2*X2(K))/2. - U3*X4(K) IF (AN(K).EQ. 90.0 .OR. AN(K) .EQ. -90.0 ) QB(L+5)= 0.0 QB(L+6) = U5 - U3*Y4(K) SB(L+1) = (3.*PA+PC)/8. + PB*Y2(K)/2. + (PA-PC)*Y4(K)/8. SB(L+2) = (6.+2. *Y4(K))*S2(2)/8. + (PA-S2(4))*(1-Y4(K))/8. SB(L+3) = PB/2. * X2(K) + X4(K)*(PA-PC)/4. IF (AN(K).EQ. 90.0 .OR. AN(K) .EQ. -90.0 ) SB(L+3)= 0.0 SB(L+4) = (3.*PA + PC)/8. - PB * Y2(K)/2. + Y4(K)*(PA-PC)/8. SB(L+5) = PB/2.*X2(K) + (PC-PA)*X4(K)/4. IF (AN(K).EQ. 90.0 .OR. AN(K) .EQ. -90.0 ) SB(L+5)= 0.0 SB(L+6) = PA/2. + S2(4)/2. - S2(2) + (PC-PA)*Y4(K)/2. C c CALCULATE ANGLE-PLY LAMINATE ELASTIC CONSTANTS C QDEL = QB(L+6)*(QB(L+1)*QB(L+4) - QB(L+2)*QB(L+2)) J = J + 1 APEX(J) = QDEL / (QB(L+4)*QB(L+6)) APEY(J) = QDEL / (QB(L+1)*QB(L+6)) APNXY(J) = QB(L+2) / QB(L+4) C WHY DOES APNXY = APNYX ??? APNYX(J) = QB(L+2) / QB(L+4) APGXY(J) = QB(L+6) C C CALCULATE LAMINATE CTE's ALPABX AND ALPABY C QD = QDEL/QB(L+6) ALPXB(J) = ALPX(K) 1 + (QB(L+4)*QB(L+3) - QB(L+2)*QB(L+5)) * ALPXY(K)/QD ALPYB(J) = ALPY(K) 2+ (-QB(L+2)*QB(L+3) + QB(L+1)*QB(L+5)) * ALPXY(K)/QD C Calculate and store material principal stresses for unit C axial stress on angle-ply laminate sig1(k) = (yy(k)*qb(l+1)+xx(k)*qb(l+2)+2.*Y(k)*X(k)*qb(l+3) 1 -APNXY(j)*(yy(k)*qb(l+2)+xx(k)*qb(l+4)+2.*Y(k)*X(k)*qb(l+5))) 2 / APEX(j) sig2(k) = (xx(k)*qb(l+1)+yy(k)*qb(l+2) - 2.*Y(k)*X(k)*qb(l+3) 1 -APNXY(j)*(xx(k)*qb(l+2)+yy(k)*qb(l+4)-2.*Y(k)*X(k)*qb(l+5))) 2 / APEX(j) sig12(k) =(y(k)*x(k)*(-qb(l+1)+qb(l+2)+APNXY(j)*(qb(l+2) 1 -qb(l+4)))+(yy(k)-xx(k))*(qb(l+3)-APNXY(j)*qb(l+5))) 2 / APEX(j) EPS1(k) = S2(1)*sig1(k) + s2(2)*sig2(k) EPS2(k) = S2(2)*sig1(k) + s2(3)*sig2(k) EPS12(k) = S2(4)*sig12(k) WRITE (14,960) AN(K),sig1(k),sig2(K),sig12(k), $EPS1(k),EPS2(k),EPS12(k) C C Print Angle-ply Engineering constants C IF (AN(K).GT.90.0) GO TO 100 IF (AN(K).LT.0.0) GO TO 100 WRITE (16,950)AN(K),APEX(J),APEY(J),APNXY(J),APGXY(J),ALPXB(J), 1ALPYB(J) C C Write Angle-ply Engineering constants to Unit 17 C APTSR(J) = (APEX(J)/APEY(J))*(2.0 - 1APNXY(J)*APEY(J)/APEX(J))/(1.0 - 2.0 * APNXY(J)) WRITE (17,950) AN(K),APEX(J),APEY(J),APNXY(J),APGXY(J),ALPXB(J), 1ALPYB(J),APTSR(J) 950 FORMAT (3X,F6.1,7(1X,E11.4)) 100 L = L + 6 960 FORMAT (F6.1,6(x,E10.3)) IF (K4 .LE. 0) WRITE (19,956) 956 FORMAT(4x,'OFF-AXIS LAMINA PROPERTIES',/, $4x'AN(I)',8x,'Ex',9x,'Ey',9x,'NUxy',7x,'Gxy',9x, 1,'ETAxy,x',4x,'ETAxy,y',5x,'ETAx,xy',5x,'ETAy,xy',6x, 2'ALPHAx',6x,'ALPHAy',6x,'ALPHAxy'/'*DATA') C C WRITE QB and SB if K1 not = 1 DO 150 I = 1, NANGLS LL = 9*(I-1) L = 6*(I-1) IF (K1.EQ.0) GO TO 50 WRITE(16,918) AN(I) 918 FORMAT (///,10X,'*********************************************', 1/,10X,'QBAR MATRIX -- ',3X,F7.3,2X,'DEGREE ANGLE'/) WRITE(16,919) (QB(L+J), J = 1,6) 919 FORMAT (10X,3(E14.6,5X)/29X,2(E14.6,5X)/48X,E14.6/) WRITE (16,920) AN(I) 920 FORMAT (10X,'SBAR MATRIX -- ',3X,F7.3,2X,'DEGREE ANGLE',/) WRITE(16,921) (SB(L+J), J = 1,6) 921 FORMAT (10X,3(E14.6,5X)/29X,2(E14.6,5X)/48X,E14.6,/) C C CALCULATE OFF-AXIS LAMINA ENGINEERING CONSTANTS C C EX = PROP(LL+1) EY = PROP(LL+2) NUXY = PROP(LL+3) C GXY = PROP(LL+4) NUYX = PROP(LL+5) ETAXY,X=PROP(LL+6) C ETAXY,Y = PROP(LL+7) ETAX,XY=PROP(LL+8) ETAY,XY=PROP(LL+9) C LL INCREASES BY NINE FOR EACH NEW ANGLE C 50 PROP(LL+1)=1./SB(L+1) PROP(LL+2)=1./SB(L+4) PROP(LL+3) = -SB(L+2)*PROP(LL+1) PROP(LL+4) = 1./SB(L+6) PROP(LL+5)=-SB(L+2)*PROP(LL+2) PROP(LL+6) = SB(L+3) * PROP(LL+1) PROP(LL+7) = SB(L+5) * PROP(LL+2) PROP(LL+8) = SB(L+3) * PROP(LL+4) PROP(LL+9) = SB(L+5) * PROP(LL+4) IF (K2.EQ.0) GO TO 150 C C store lamina properties for plotting C IF (K4 .LE. 0) WRITE (19,955) AN(I),PROP(LL+1),PROP(LL+2), $PROP(LL+3),PROP(LL+4),PROP(LL+6),PROP(LL+7),PROP(LL+8), 2PROP(LL+9),ALPX(I),ALPY(I),ALPXY(I) C C PRINT LAMINA ENGINEERING PROPERTIES C 955 FORMAT (3X,F7.3, 11(1X,E11.4)) WRITE(16,898) AN(I) 898 FORMAT (//,T15,'***LAMINA ENGINEERING CONSTANTS*** THETA =',F7.3) WRITE (16,928) PROP(LL+1),PROP(LL+2),PROP(LL+3),PROP(LL+5),PROP(LL+ 14),PROP(LL+6),PROP(LL+7),PROP(LL+8),PROP(LL+9),ALPX(I),ALPY(I), 2ALPXY(I) 928 FORMAT (/10X,'EXBAR = ',E14.6,7X,'EYBAR = ',E14.6,/,10X,'UXYBAR = 1 ',E14.6,6X,'UYXBAR = ',E14.6,/,10X,'GXYBAR = ',E14.6,6X, 2'ETAXY,X = ',E14.6,/,10X,'ETAXY,Y = ',E14.6,5X,'ETAX,XY = ',E14.6, 3/,10X,'ETAY,XY = ',E14.6,5X,'ALPX = ',E14.6,/,10X,'ALPY = ',E14.6 4,8X,'ALPXY = ',E14.6,/) 150 CONTINUE C C CALCULATE ENGINEERING PROPERTIES FOR SYMMETRIC LAMINATE C OF THE TYPE 0-T1/+-THETA-T2/90-T3 C if K4 > 0, T1, T2, T3 ARE LAYER THICKNESSES C if K4 < or = 0, SOLN IF FOR UNIDIRECTIONAL LAMINA C SAC3D CALCULATES 3-D TRANSFORMED C & S FOR ALL ANGLES C AND ALPHA Z AND NU XZ FOR ANGLE-PLY LAMINATES C if (K4 .GT. 0) call GENLAM (T1,T2,T3,QB,Q,TZERO,TINC,NOTINC 1,AN,NANGLS,ALPX,ALPY,ALPXY,AL1,AL2,GALPX,GALPY,astar) C CALL SAC3D(NANGLS,AN,ANG,C,S,ALPX,ALPY,ALPXY,K1, 1QB,GALPX,GALPY,astar,T1,T2,T3,Q,AL1,AL2,K4) C C STORE PLOT DATA IN THE [A] MATRIX C C IF (K3.EQ.0) GO TO 05 K = 1 C STORE QijBAR DO 250 I = 1, 6 DO 250 J = 1, NANGLS 01320 A(K)=QB(MQ(I) + (J-1)*6) 01330 250 K = K + 1 01340 C STORE SijBAR DO 300 I = 1, 6 01350 DO 300 J = 1, NANGLS 01360 A(K)=SB(MQ(I) + (J-1)*6) 01370 300 K = K + 1 01380 C STORE Ex, Ey, NUxy, Gxy, ETAxy,x , etc. DO 350 I = 1, 9 01390 DO 350 J = 1, NANGLS 01400 A(K)=PROP(MP(I) + (J-1)*9) 01410 350 K = K + 1 01420 C C Determine maximum value for each variable C C DETERMINE MAXIMUM VALUE (Z) FOR EACH VARIABLE C WRITE (16,908) C 908 FORMAT (/,10X,'MAX VALUES OF 21 OFF-AXIS PROPERTIES',/) DO 395 I = 1, 21 01440 LL = (I-1) * NANGLS 01450 Z = A(LL+1) 01460 C DETERMINE MAXIMUM VALUE (Z) DO 390 J = 1, NANGLS 01470 390 IF (ABS(A(LL+J)).GT.ABS(Z)) Z=(A(LL+J)) Z = ABS(Z) 01490 C WRITE VARIABLE NUMBER AND MAXIMUM VALUE C WRITE(16,907) I, Z C 907 FORMAT(5X,'VARIABLE NO. ',I5,5X,'MAX VALUE =',E14.6) C NORMALIZE w/r TO Z C DO 395 J = 1, NANGLS 01500 C 395 A(LL+J) = A(LL+J)/Z 01510 395 continue LL = NANGLS WRITE(18,929) 929 FORMAT(2x,'Angle',7X,'Q11',9X,'Q22',9X,'Q66',9x,'Q12',9x, 1'Q16',9X,'Q26',/'*DATA') WRITE (11,933) 933 FORMAT(2x,'Angle',7X,'S11',9X,'S22',9X,'S66',9x,'S12',9x, 1'S16',9X,'S26',/'*DATA') WRITE (12,934) 934 FORMAT(8x,'Angle',4X,'Ex',15X,'Ey',12X,'Gxy',12x,'NUxy',9X, 2'NUyx',9X,'ETAxy,x',/'*DATA') WRITE (13,932) 932 FORMAT (2X,'Angle',4X,'ETAxy,y',4X,'ETAx,xy',4X,'ETAy,xy',4X, 3/'*DATA') DO 410 J = 1, NANGLS C WRITE: ANGLE, Q11, Q22, Q66, Q12, Q16, Q26 WRITE(18,950) AN(J), A(J),A(LL+J),A(LL*2+J),A(LL*3+J),A(LL*4+J), 1A(LL*5+J) C WRITE: ANGLE, S11, S22, S66, S12, S16, S26 WRITE(11,950) AN(J), A(LL*6+J),A(LL*7+J),A(LL*8+J),A(LL*9+J), 1A(LL*10+J),A(LL*11+J) C WRITE: ANGLE, Ex, Ey, Gxy, NUxy, NUyx, ETAxy,x WRITE(12,950) AN(J), A(LL*12+J),A(LL*13+J),A(LL*14+J),A(LL*15+J), 1A(LL*16+J),A(LL*17+J) C WRITE: ANGLE, ETAxy,y , ETAx,xy , ETAy,xy WRITE(13,950) AN(J), A(LL*18+J),A(LL*19+J),A(LL*20+J) 410 CONTINUE 906 FORMAT (2X,F8.1,2(3X,E14.6)) 01530 IF (AN(1).NE.(-AN(NANGLS))) GO TO 05 WRITE(16,945) 945 FORMAT (//,20X,'ANGLE-PLY ETAXY,X & NUxy MISMATCH VALUES',//, 110X,'ANGLE',5X,'DELTA ETA',5X,'DELTA Nu',/) WRITE (23,936) 936 FORMAT (2x,'ANGLE-PLY LAMINATES ETA MISMATCH',/, $4X,'Angle',4x,'ETAxy,x mismatch',2x,'NUxy mismatch',/'*DATA') JJ = 21 * NANGLS MAX = (NANGLS-1)/2 + 1 DO 400 I = 1, MAX LL = 9*(I-1) KK = (NANGLS - I) * 9 MM = NANGLS + 1 - I A(JJ+I) = ABS(PROP(LL+6)-PROP(KK+6)) WRITE (16,906) AN(MM), A(JJ+I), 0.0 WRITE (23,906) AN(MM), A(JJ+I), 0.0 400 CONTINUE GO TO 05 1000 WRITE (16,905) 01550 905 FORMAT (///,10X,'PROBLEMS COMPLETE') 01560 STOP END C C **** SUBROUTINE GENLAM C SUBROUTINE GENLAM (T1,T2,T3,QB,Q,TZERO,TINC,NOTINC,AN, 1NANGLS,ALPX,ALPY,ALPXY,AL1,AL2,GALPX,GALPY,astar) C C CALCULATES Ex, Ey, nuxy, Gxy, Alphax and Alphay C for general laminates C IMPLICIT REAL*8(A-H,O-Z) DIMENSION QB(1086),Q(4),Ex(181),Ey(181),GLnuxy(181), 1Gxy(181),AN(181),ALPX(181),ALPY(181),GALPX(181),GALPY(181), 2ALPXY(181),astar(4,181) WRITE (20,900) L = 0 DO 100 I = 1, NANGLS A1 = Q(1)*T1 + T2*QB(L+1) + Q(3)*T3 A2 = Q(2)*(T1+T3) + T2*QB(L+2) A3 = Q(3)*T1 + T2*QB(L+4) + Q(1)*T3 A4 = Q(4)*(T1+T3) + T2*QB(L+6) THICK = T1 + T2 + T3 D1 = A1*A3 - A2*A2 D = THICK / D1 astar(1,I) = D * A3 astar(2,I) = - D * A2 astar(3,I) = D * A1 astar(4,I) = THICK / A4 Ex(I) = 1. / astar(1,I) Ey(I) = 1. / astar(3,I) GLnuxy(I) = - astar(2,I) / astar(1,I) Gxy(I) = 1. / astar(4,I) QA1 = T1*(Q(1)*AL1 + Q(2)*AL2) + T2*(QB(L+1)*ALPX(I) 1 + QB(L+2)*ALPY(I) + QB(L+3)*ALPXY(I)) 2+ T3*(Q(3)*AL2 + Q(2)*AL1) QA2 = T1*(Q(2)*AL1 + Q(3)*AL2) + T2*(QB(L+2)*ALPX(I) 1 + QB(L+4)*ALPY(I)+ QB(L+5)*ALPXY(I)) 2+ T3*(Q(2)*AL2 + Q(1)*AL1) GALPX(I) = (astar(1,I)*QA1 + astar(2,I)*QA2) / THICK GALPY(I) = (astar(2,I)*QA1 + astar(3,I)*QA2) / THICK WRITE (20,901) AN(I),Ex(I),Ey(I),GLnuxy(I), 1Gxy(I),GALPX(I),GALPY(I) 100 L = L + 6 RETURN 900 FORMAT (//,5X,'*** ENGINEERING PROPERTIES FOR GENERAL 0/+-th/90 1LAMINATES ***',/,7X,'ANG',10X,'Ex',14X,'Ey',14X, 2'nuxy',10X,'Gxy',14X,'ALPX',14X'ALPY',/,'*DATA') 901 FORMAT (5x,F6.1,2(3X,E14.6),3X,F7.3,3(3X,E14.6)) END C C C C ***** SUBROUTINE SAC3D ***** C SUBROUTINE SAC3D (NANGLS,AN,ANG,C,S,ALPX,ALPY, 1ALPXY,K1,QB,GAPLX,GAPLY,ASTAR,ZT1,ZT2,ZT3,Q,AL1,AL2,K4) C STIFFNESS & COMPLIANCE FOR TRANSVERSELY ISOTROPIC 3-D MATERIAL C PROGRAM NAME -- 3DSAC -- C 1,T1I,T2I IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DSIN,DCOS DIMENSION S(6,6),C(6,6),SB(6,6,181),CB(6,6,181),ANG(181), 1QB(1086),ALPX(181),ALPY(181),ALPXY(181), 2AN(181),T1(6,6),T2(6,6),A(6,6),B(6,6),T1I(6,6),T2I(6,6), 3astar(4,181),Q(4),GAPLX(181),GAPLY(181) C C BEGIN DO LOOP FOR EACH ANGLE C C WRITE HEADERS TO TAPE WRITE (21,908) 908 FORMAT (/,3X,'3-D SYMMETRIC STIFFNESS TERMS FOR PLOTTING', 1/,3X,'AN',6X,'C(1,1)',6X,'C(1,2)',6X,'C(1,6)',6X,'C(2,2)', 26X,'C(2,3)',6X,'C(2,6)',6X,'C(3,6)',6x,'C(4,4)',6x,'C(4,5)', 3'C(6,6)',6X,'C(1,3)',6X,'C(5,5)',6X,'C(3,3)'/,'*DATA') WRITE (22,909) 909 FORMAT (/,3X,'3-D SYMMETRIC COMPLIANCE TERMS FOR PLOTTING', 1/,3X,'AN S(1,1) S(1,2) S(1,6) S(2,2) S(2,3) 2S(2,6) S(3,6) S(4,4) S(4,5) S(6,6) S(1,3) 3S(5,5) S(3,3)'/,'*DATA') 20 DO 200 L = 1, NANGLS DO 120 I =1,6 DO 120 J = 1,6 A(I,J) = 0. B(I,J) = 0. CB(I,J,L) = 0. 120 SB(I,J,L) = 0. C C CONSTRUCT TRANSFORMATION MATRICES FOR ROTATION TO X-Y COORDINATES ANGR = ANG(L) XM = DCOS(ANGR) XX = XM * XM YN = DSIN(ANGR) YY = YN * YN XY = XM * YN T1(1,1) = XX T1(1,2) = YY T1(1,6) = 2. * XM * YN T1(2,1) = T1(1,2) T1(2,2) = T1(1,1) T1(2,6) = -T1(1,6) T1(3,3) = 1.0 T1(4,4) = XM T1(4,5) = -YN T1(5,4) = -T1(4,5) T1(5,5) = XM T1(6,1) = -XM * YN T1(6,2) = -T1(6,1) T1(6,6) = XX - YY C DO 50 I = 1,6 DO 50 J = 1,6 50 T2(I,J) = T1(I,J) C T2(1,6) = 0.5 * T2(1,6) T2(2,6) = 0.5 * T2(2,6) T2(6,1) = 2. * T2(6,1) T2(6,2) = 2. * T2(6,2) C DO 75 I = 1,6 DO 75 J = 1,6 T1I(I,J) = T1(I,J) 75 T2I(I,J) = T2(I,J) C T1I(1,6) = -T1I(1,6) T1I(2,6) = -T1I(2,6) T1I(4,5) = -T1I(4,5) T1I(5,4) = -T1I(5,4) T1I(6,1) = -T1I(6,1) T1I(6,2) = -T1I(6,2) T2I(1,6) = -T2I(1,6) T2I(2,6) = -T2I(2,6) T2I(4,5) = -T2I(4,5) T2I(5,4) = -T2I(5,4) T2I(6,1) = -T2I(6,1) T2I(6,2) = -T2I(6,2) C DO 100 I = 1,6 DO 100 J = 1,6 DO 100 K = 1,6 A(I,J) = A(I,J) + T1I(I,K) * C(K,J) 100 B(I,J) = B(I,J) + T2I(I,K) * S(K,J) C DO 125 I = 1,6 DO 125 J = 1,6 DO 125 K = 1,6 CB(I,J,L) = CB(I,J,L) + A(I,K) * T2(K,J) 125 SB(I,J,L) = SB(I,J,L) + B(I,K) * T1(K,J) C C C WRITE CB AND SB TO TAPE FOR PLOTTING C WRITE(21,910) AN(L),CB(1,1,L),CB(1,2,L),CB(1,6,L),CB(2,2,L), 1CB(2,3,L),CB(2,6,L),CB(3,6,L),CB(4,4,L),CB(4,5,L),CB(6,6,L), 2CB(1,3,L),CB(5,5,L),CB(3,3,L) WRITE(22,910) AN(L),SB(1,1,L),SB(1,2,L),SB(1,6,L),SB(2,2,L), 1SB(2,3,L),SB(2,6,L),SB(3,6,L),SB(4,4,L),SB(4,5,L),SB(6,6,L), 2SB(1,3,L),SB(5,5,L),SB(3,3,L) IF (K1.EQ.0) GO TO 200 WRITE(16,906) AN(L),((CB(I,J,L),J=1,6),I=1,6) WRITE(16,907) ((SB(I,J,L),J=1,6),I=1,6) 906 FORMAT (//6X,'3D-STIFFNESS & COMPLIANCE THETA =',F6.1,' DEG.', 1//20X,'**3D STIFFNESS CBA***',/,6(6(3X,E10.3)/)) 907 FORMAT (//20X,'**3D COMPLIANCE SBAR **',/,6(6(3X,E10.3)/)) 200 CONTINUE 910 FORMAT (2X,F6.1,13(X,E13.6)) CALL ALPAZ(ALPX,ALPY,ALPXY,QB,SB,NANGLS,ANG,AN, 1ZT1,ZT2,ZT3,S,Q,AL1,AL2,GAPLX,GAPLY,ASTAR,K4) RETURN END C C ************** SUBROUTINE ALPAZ *********************** C SUBROUTINE ALPAZ (ALPX,ALPY,ALPXY,QB,SB,NANGLS,ANG, 1AN,T1,T2,T3,S,Q,AL1,AL2,GALPX,GALPY,ASTAR,K4) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ALPZB(181),ALPX(181),ALPY(181),GNUYZ(181), 1QB(1086),SB(6,6,181),ALPXY(181),ANG(181),GNUXZ(181),AN(181), 2S(6,6),Q(4),GALPX(181),GALPY(181),F1(181),F2(181),ASTAR(4,181) C C CALCULATE LAMINATE ALPHA Z C WRITE (16,900) 900 FORMAT(/,5X,'LAMINA or LAMINATE THROUGH-THICKNESS PROPERTIES', 1/,7X,'THETA',6X,'ALPHAZ',11X,'NUXZ',9X,'NUYZ',/) WRITE (10,901) 901 FORMAT(/,5X,'LAMINA or LAMINATE THROUGH-THICKNESS PROPERTIES', 1/,7X,'THETA',6X,'ALPHAZ',11X,'NUXZ',9X,'NUYZ',/,'*DATA',/) THICK = T1 + T2 + T3 DO 200 I = 1, NANGLS IF ((AN(I) .LT. 0.0) .AND. (K4 .GT. 0) ) GO TO 200 L = 6 * (I-1) C K4 > 0, GENERAL LAMINATE, K4 < or = 0, Off-AXIS LAMINA IF (K4 .LE. 0) GO TO 100 ALPZB(I) = (((GALPX(I)-ALPX(I))*(SB(3,1,I)*QB(L+1) 1+SB(3,2,I)*QB(L+2) + SB(3,6,I)*QB(L+3)) + 2(GALPY(I)-ALPY(I))*(SB(3,1,I)*QB(L+2) 3+SB(3,2,I)*QB(L+4) + SB(3,6,I)*QB(L+5)) + 4(-ALPXY(I))*(SB(3,1,I)*QB(L+3)+SB(3,2,I)*QB(L+5) 5+ SB(3,6,I)*QB(L+6)) + AL2)*T2 + T1 *( S(3,1)* 6(Q(1)*(GALPX(I)-AL1) + Q(2)*(GALPY(I)-AL2)) + S(3,2)* 7(Q(2)*(GALPX(I)-AL1) + Q(3)*(GALPY(I)-AL2)) + AL2) + 8T3 * (S(3,2)* 9(Q(3)*(GALPX(I)-AL2) + Q(2)*(GALPY(I)-AL1)) + S(3,1)* A(Q(2)*(GALPX(I)-AL2) + Q(1)*(GALPY(I)-AL1)) + AL2))/THICK F1(I) = T1*(S(3,1)*Q(1) + S(3,2)*Q(2)) + 1T3 *(S(3,2)*Q(3) + S(3,1)*Q(2)) + 2T2 *(SB(1,3,I)*QB(L+1) + SB(2,3,I)*QB(L+2) + 3SB(3,6,I)*QB(L+3)) F2(I) = T1*(S(3,1)*Q(2) + S(3,2)*Q(3)) + 1T3*(S(3,2)*Q(2) + S(3,1)*Q(1)) + 2T2*(SB(1,3,I)*QB(L+2) + SB(2,3,I)*QB(L+4) + 3SB(3,6,I)*QB(L+5)) C GNUXZ(I) = -(ASTAR(1,I)*F1(I) + 1ASTAR(2,I)*F2(I))/(THICK * ASTAR(1,I)) C GNUYZ(I) = -(ASTAR(2,I)*F1(I) + 1ASTAR(3,I)*F2(I))/(THICK * ASTAR(3,I)) GO TO 300 100 GNUXZ(I) = -SB(3,1,I) / SB(1,1,I) GNUYZ(I) = -SB(2,3,I) / SB(2,2,I) ALPZB(I) = AL2 300 WRITE(10,905) AN(I),ALPZB(I),GNUXZ(I),GNUYZ(I) WRITE(16,905) AN(I),ALPZB(I),GNUXZ(I),GNUYZ(I) 905 FORMAT (5X,F6.1,3(2X,E14.6)) C WRITE ALPZB(I),GNUXZ(I),GNUYZ(I) 200 CONTINUE RETURN END C SUBROUTINE PROP3D (K,QB,E1,E2,GNU12,GNU23,G,S,C) C CALCULATES PRINCIPAL 3D STIFFNESS AND COMPLIANCE C for TRANSVERSELY ISOTROPIC MATERIAL C IMPLICIT REAL*8(A-H,O-Z) DIMENSION S(6,6),C(6,6) DO 10 I=1,6 DO 10 J = 1,6 S(I,J) = 0. C(I,J) = 0. 10 CONTINUE S(1,1) = 1./E1 S(1,2) = -GNU12/E1 S(1,3) = S(1,2) S(2,2) = 1./E2 S(2,3) = -GNU23/E2 S(3,3) = S(2,2) S(4,4) = 2.*(S(2,2) - S(2,3)) S(5,5) = 1./G S(6,6) = S(5,5) S(2,1) = S(1,2) S(3,2) = S(2,3) S(3,1) = S(1,3) SDET = S(1,1)*S(2,2)*S(3,3)-S(1,1)*S(2,3)*S(2,3)-S(2,2)*S(1,3)* 1S(1,3)-S(3,3)*S(1,2)*S(1,2)+2.*S(1,2)*S(2,3)*S(1,3) C(1,1)=(S(2,2)*S(3,3)-S(2,3)*S(2,3))/SDET C(1,2)=(S(1,3)*S(2,3)-S(1,2)*S(3,3))/SDET C(1,3)=(S(1,2)*S(2,3)-S(1,3)*S(2,2))/SDET C(2,2)=(S(3,3)*S(1,1)-S(1,3)*S(1,3))/SDET C(2,3)=(S(1,2)*S(1,3)-S(2,3)*S(1,1))/SDET C(3,3)=(S(1,1)*S(2,2)-S(1,2)*S(1,2))/SDET C(4,4)=1./S(4,4) C(5,5)=1./S(5,5) C(6,6)=1./S(6,6) C(2,1)=C(1,2) C(3,1)=C(1,3) C(3,2)=C(2,3) WRITE(16,904) ((S(I,J),J=1,6),I=1,6) WRITE (16,905) ((C(I,J),J=1,6),I=1,6) 904 FORMAT (//15X,'3D COMPLIANCE IN 1-2 COORDNATE SYSTEM',/ 1,6(6(3X,E10.3)/)) 905 FORMAT (//15X,'3D STIFFNESS IN 1-2 COORDINATE SYSTEM',/ 1,6(6(3X,E10.3)/)) 20 RETURN END