C CLAP C C Fortran program for analysis of unidirectional lamina, C laminates and 3-D properties of anisotropic materials C thermal-mechanical loading, interlaminar forces and C moments, and failure analysis. C The theory behind the program can be found in Chapters C 3,4,5, 8 and 9 of 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 Modified Feb. 15, 1997 to print strain distribution to a C file for plotting C Modified Dec. 13, 1995 for Interlaminar Force & Moment C when loading is moment Mx. C Modified Oct. 20, 1994 to calculate G23 C Modified Nov. '93 for improvement to thru-thickness C Force and Moment Calculations and plotting. C Modified Dec. 4, 1995 for plot file naming C Modified Dec. 6, 1992 to include stress & strain C calculations at top & bottom of each ply C clapv2 - modified for incremental TTT interlaminar C forces and moments including plot file Jan. 30, '91. C C **** CLAPALPAZ IS MODIFIED VERSION TO ADD CTE ALPAZ 5-4-88 C AND THROUGH-THE-THICKNESS NU,XZ 2-14-84 C MODIFIED 3-5-88 FOR UNSYMMETRIC LAMINATES C COMPOSITE LAMINATE ANALYSIS PROGRAM (CLAP) C THIS VERSION HAS CORRECT TEMPERATURE DEPENDENT PROPERTY SOLN C AS OF 3-12-83 C PROGRAM WRITTEN BY C. T. HERAKOVICH, APPLIED MECHANICS PROGRAM C UNIVERSITY OF VIRGINIA C CHARLOTTESVILLE, VIRGINIA 22903, PHONE (804) 924-3605 C EMAIL: HERAK@VIRGINIA.EDU C IMPLICIT REAL*8(A-H,O-Z) C DOUBLE PRECISION X,Y,Y2,X2,Y4,X4,DSIN,DCOS,AN,G,A,AI,AD,ANG,AP,AW, C 1B,BS,BP,BW,C,CS,CP,CW,D,DI,DSI,DS,XX,YY,ALPABA,BETABA CHARACTER*4 ITITLE(18) DIMENSION A(9),ANG(10),B(9),D(9),Q(20),QB(120),S(20),NA(20),KEY(25 1),H(21),MATL(20),EL(5),ET(5),ULT(5),UTL(5),G(5),U1(5),U2(5),U3(5), 2U4(5),X(10),Y(10),XX(10),YY(10),ALPA(10),ALPAX(60),CURV(3), 3EPS(60),EPSX(60),SIG(60),SIGX(60),SIGRX(60),ENAUT(3), 4FORC(3),BEND(3),EPSXT(60),SIGXT(60),TFORC(3),TBEND(3), HFORC(3),HB 5END(3),TOBEND(3),TENAUT(3),TCURV(3),TOENUT(3),TOCURV(3),BETA(10),B 6ETAX(60),U23(5), HENAUT(3),HCURV(3),SIGMX(60),EPSMX(60),XT(5),XC( 75),YT(5),YC(5),SS(5),F12(5),ELM(5),ETM(5),GM(5),XTM(5),XCM(5),YTM( 85),YCM(5),SSM(5),F12M(5),ALPAM(10),BETAM(10),ETT(5),ULTT(5), 9GT(5),XTT(5),XCT(5),YTT(5),YCT(5),SST(5),F12T(5),EPSRX(60),ULTM( A5),TOFORC(3),RTFORC(3),RTBEND(3),ALPAT(10),BETAT(10),RTCURV(3),RTE BNAT(3),ELT(5),SB(120),FTSTRN(60),FMSTRN(60),HRFORC(3),HRBEND C(3),HRCURV(3),HRENAT(3),ALPAMX(60),BETAMX(60),ALPABA(3),BETABA(3) DIMENSION AI(9),AN(10),AP(9),AW(9),BS(9),BP(9),BW(9),CS(9),CP(9) 1,CW(9),DI(9),DSI(9),DS(9),DP(9),DUMMY(60),U23M(5),G23(5),G23M(5) DIMENSION Y2(10),X2(10),Y4(10),X4(10),F(25) C quantitites with leading B correspond to Bottom of ply DIMENSION BEPSX(60),BSIGMX(60),BSIGXT(60),BEPSXT(60), 2BEPSMX(60),BEPSRX(60),BSIGX(60),BSIGRX(60),BEPSRX(60),HB(21) CHARACTER*22 IFNAME,output,intFMplt,strplt,stainplt CHARACTER*8 PROBLM C C Prompt the user for the names of the input files and C standard problem name lead-in for output files C Read the names from the screen. C WRITE (*,3) 3 FORMAT(/'Enter your input data file name (default=input.data):') READ (*,8) IFNAME IF(IFNAME.EQ.' ') IFNAME='input.data' 8 FORMAT (A22) WRITE (*,4) 4 FORMAT(/'Enter problem name 8 character max (default=PROBLM):') READ (*,5) PROBLM 5 FORMAT(A8) do 6 i=8,2,-1 if (PROBLM(i:i).gt.' ') goto 7 6 continue 7 IF(PROBLM.EQ.' ') PROBLM='PROBLM' output = PROBLM(:i)//'.output' strplt = PROBLM(:i)//'.strplt' intFMplt = PROBLM(:i)//'.intFMplt' stainplt = PROBLM(:i)//'.stainplt' C WRITE (*,6) C 6 FORMAT(/'Enter interlaminar plot output file(deflt=plot.output):') C READ (*,5) PFNAME C IF(PFNAME.EQ.' ') PFNAME='plot.output' C WRITE (*,7) C 7 FORMAT(/'Enter stress plot output file (deflt=stress.plot):') C READ (*,5) SFNAME C IF(SFNAME.EQ.' ') SFNAME='stress.plot' 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=intFMplt,STATUS='UNKNOWN') REWIND 17 OPEN (UNIT=18,FILE=strplt,STATUS='UNKNOWN') REWIND 18 OPEN (UNIT=19,FILE=stainplt,STATUS='UNKNOWN') REWIND 19 1111 READ (15,900) ITITLE KWRITE = 0 DO 01 I = 1, 60 01 DUMMY(I) = 0.0 DO 02 I = 1,3 TBEND(I) = 0.0 HFORC(I) = 0.0 HBEND(I) = 0.0 02 TFORC(I) = 0.0 900 FORMAT (18A4) WRITE (16,888) 888 FORMAT(////,10x,'****************************************',/, a10x,'*************************************************',/, 110X,'CLAP vJ98 - COMPOSITE LAMINATE ANALYSIS PROGRAM',/, 210X,'CARL T. HERAKOVICH',/, 310X,'UNIVERSITY OF VIRGINIA (804) 924-3605',/, 410X,'EMAIL: HERAK@VIRGINIA.EDU',/ 510X,'For theory see MECHANICS OF FIBROUS COMPOSITES by ',/ 610X,'CARL T. HERAKOVICH, John Wiley & Sons, Inc., 1998 ',/ 710X,'******************************************************',//) WRITE (16,901) ITITLE 901 FORMAT (//,18A4,///) 902 FORMAT (5I5,F10.0) READ (15,*) KEY WRITE(16,955) (KEY(I), I=1,25) 955 FORMAT (//T10,'KEYS',/T10,13I3,/T10,12I3,/) IF (KEY(2).LT.0) STOP WRITE (17,900) ITITLE WRITE (18,900) ITITLE WRITE (19,900) ITITLE READ(15,*,END=1000) NPLYS,NMATLS,NANGLS,ISYN,KDH,HO 899 FORMAT (25I2) C C LAMINATE CONFIGURATION C WRITE (16,903) NPLYS,NMATLS,NANGLS,ISYN,KDH,HO 903 FORMAT (10X,'NPLYS = ',I5,5X,'NMATLS = ',I5,5X,'NANGLS = ',I5,/ 1,10X,'ISYM = ',I3,3X,'KDH = ',I3,3X,'HALF THICKNESS = ',E14.6) 15 READ(15,*) (NA(I),H(I),I = 1,NPLYS) C C H(I) IS Z COORDINATE TO TOP OF I TH PLY C 904 FORMAT (6(I2,F10.0)) READ(15,*) (ANG(I),MATL(I),I = 1,NANGLS) 905 FORMAT (8(F8.0,I2)) WRITE (16,906) 906 FORMAT(//'0',T30,'***LAMINATE GEOMETRY***',//, 111X,'PLY',5X,'MATL',5X,'ANG',5X,'ANG NO.',10X,'TH', 210X,'BH',/) 907 FORMAT (11X,I2,7X,I2,5X,F6.1,5X,I2,4X,E14.6,4X,E14.6/) H(NPLYS+1) = HO IF (ISYN.EQ.0) GO TO 16 DO 17 I = 1,NPLYS NA(NPLYS + I) = NA(NPLYS+1-I) 17 H(2*(NPLYS+1)-I) = -H(I) H(NPLYS+1) = 0.0 NPLYS = NPLYS * 2 C WRITE (16,897) C 897 FORMAT(/T30,'*** SYMMETRIC LAMINATE ***',//) C C HO IS TH HALF THICKNESS OF THE LAMINATE C 16 DO 19 I = 1, NPLYS 19 HB(I) = H(I + 1) C HB(I) IS Z TO BOTTOM OF PLY I WRITE(16,907)(I,MATL(NA(I)),ANG(NA(I)),NA(I),H(I), 1HB(I),I =1,NPLYS) DO 20 I =1,NANGLS AN(I) = ANG(I) ANG(I) = ANG(I)*3.14159265/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)) 20 CONTINUE L = 0 C C READ STIFFNESS, THERMAL ,MOISTURE & STRENGTH PARAMETERS C K = 0 DO 30 I = 1, NMATLS L1 = L + 1 L2 = L + 2 READ(15,*) EL(I),ET(I),ULT(I),G(I),ALPA(L1),ALPA(L2),BETA(L1), 1BETA(L2),XT(I),XC(I),YT(I),YC(I),SS(I),F12(I),U23(I) READ(15,*) ELM(I),ETM(I),ULTM(I),GM(I),ALPAM(L1),ALPAM(L2), 1BETAM(L1),BETAM(L2),XTM(I),XCM(I),YTM(I),YCM(I),SSM(I), 2F12M(I),U23M(I) 975 FORMAT (8E10.0) G23(I) = ET(I) / (2.0 * (1.0 + U23(I))) G23M(I) = ETM(I) / (2.0 * (1.0 + U23M(I))) WRITE(16,970) ELM(I),EL(I),ETM(I),ET(I),GM(I),G(I),ULTM(I),ULT(I), 1U23M(I),U23(I),G23M(I),G23(I),ALPAM(L1),ALPA(L1),ALPAM(L2), 2ALPA(L2),BETAM(L1),BETA(L1),BETAM(L2),BETA(L2),XTM(I),XT(I), 3XCM(I),XC(I),YTM(I),YT(I),YCM(I),YC(I),SSM(I),SS(I), 4F12M(I),F12(I) F(K+1) = 1./XT(I) + 1./XC(I) F(K+2) = 1./YT(I) + 1./YC(I) F(K+3) = -1./(XT(I)*XC(I)) F(K+4) = -1./(YT(I)*YC(I)) F(K+5) = 1./(SS(I)*SS(I)) K = K + 5 30 L = L + 2 908 FORMAT (8E10.0) 970 FORMAT (//T10,'**** LINEAR PROPERTY DEPENDENCE FROM RM TEMP ***',/ O T10,'E11 = ',E14.6,' T + ',E14.6,/, 1 T10,'E22 = ',E14.6,' T + ',E14.6,/, 2 T10,'G12 = ',E14.6,' T + ',E14.6,/, 3 T10,'U12 = ',E14.6,' T + ',E14.6,/, E T10,'U23 = ',E14.6,' T + ',E14.6,/, E T10,'G23 = ',E14.6,' T + ',E14.6,/, 4 T10,'ALPA1 = ',E14.6,' T + ',E14.6,/, 5 T10,'ALPA2 = ',E14.6,' T + ',E14.6,/, 6 T10,'BETA1 = ',E14.6,' T + ',E14.6,/, 7 T10,'BETA2 = ',E14.6,' T + ',E14.6,/, 8 T10,'XT = ',E14.6,' T + ',E14.6,/, 9 T10,'XC = ',E14.6,' T + ',E14.6,/, A T10,'YT = ',E14.6,' T + ',E14.6,/, B T10,'YC = ',E14.6,' T + ',E14.6,/, C T10,'SS = ',E14.6,' T + ',E14.6,/, D T10,'F12 = ',E14.6,' T + ',E14.6,/) C IF (KEY(18).EQ.0) WRITE (16,1221) 1221 FORMAT (/10X,'KEY(18) = 0,CONSTANT ROOM TEMPERATURE PROPERTIES',/) IF (KEY(18).GT.0) WRITE(16,1222) 1222 FORMAT (/10X,'KEY(18) > 0, PROPERTIES TEMPERATURE DEPENDENT',/) TTT = 0.0 CALL PRPRTY(X2,X4,Y2,Y4,TTT,HO,AI,AN,AP,AW,BS,BP,BW,CS,CP,CW, 1DI,DSI,DS,DP,A,B,D,Q,QB,S,NA,KEY,H,MATL,EL,ET,ULT,G,U1,U2,U3,U4, 2X,Y,XX,YY,ALPA,ALPAX,ALPABA,BETA,BETAX,BETABA,XT,XC,YT,YC,SS, 3F12,ELM,ETM,GM,XTM,XCM,YTM,YCM,SSM,F12M,ALPAM,BETAM,ETT,ULTT,GT, 4XTT,XCT,YTT,YCT,SST,F12T,ULTM,ALPAT,BETAT,ELT,SB,NANGLS,NMATLS, 5ALPAM,BETAMX,NPLYS,UTL) IF (KEY(18).EQ.0.AND.KWRITE.GT.0) GO TO 31 KWRITE = KWRITE + 1 CALL RITPRP (NMATLS,EL,ET,ULT,G,ALPA,BETA,XT,XC,YT,YC,SS,F12) IF (KEY(2).EQ.0) GO TO 31 C C 3-D PROBLEM IF KEY(2) .GT. 0 C CALL SAC3D (KEY,H,AI,QB,NANGLS,EL,ET,ULT,G,U23,NPLYS,ANG,NA,HO, 1ALPABA,ALPAX,ALPA,MATL) C 31 IF (KEY(4).EQ.0) GO TO 2000 IF (KEY(4).GT.0) GO TO 5000 GO TO 1111 C KEY(4) > 0 IMPLIES LAMINA SOLUTION ONLY C KEY(4) < 0 IMPLIES LAMINATE CONSTANTS ONLY C KEY(5) = NUMBER OF LAMINATE LOAD CASES 5000 CALL LAMINA(QB,SB,ANG,X,Y,XX,YY,KEY,EPSX,SIGX,EPS,SIG,NPLYS,NA) GO TO 1111 2000 K2 = KEY(5) LL = 3 DO 4000 K1 = 1,K2 LL = LL + 3 DO 400 I = 1,3 ENAUT(I) = 0. CURV(I) = 0. FORC(I) = 0. 400 BEND(I) = 0. J = KEY(LL) WRITE(16,951) K1 951 FORMAT ('0',T20,'******************************************* ',//, 1T35,'LOADCASE NO',I4,//) GO TO (2010,2020,2030,2040), J C C GIVEN MIDPLANE STRAINS AND CURVATURES C 2010 READ(15,*) ENAUT,CURV,CURET,RMT,OPRT,HMC WRITE (16,895) CURET,RMT,OPRT,HMC CALL MULT2N(FORC,BEND,A,B,B,D,ENAUT,CURV) WRITE(16,971) (ENAUT(I),CURV(I),FORC(I),BEND(I),I=1,3) 971 FORMAT ('0',T20,'***GIVEN MIDPLANE STRAINS AND CURVATURES',/, 1T12,'ENAUT',12X,'KAPPA',15X,'N',12X,'MOMENTS',/,4(3X,E14.6)) 895 FORMAT (//T10,'***ROOM TEMPERATURE RESULT FOR MECHANICAL LOADING O 1NLY ***',/,T10,'CURE TEMPERATURE =',F7.1,5X,'ROOM TEMPERATURE = 2',F7.1,/,T10,'OPERATING TEMPERATURE =',F7.1,5X,'MOISTURE CONTENT = 3',F7.3,/) IF (ENAUT(1).EQ.0.0.AND.ENAUT(2).EQ.0.0.AND.ENAUT(3).EQ.0.0.AND.CU 1RV(1).EQ.0.0.AND.CURV(2).EQ.0.0.AND.CURV(3).EQ.0.0) GO TO 3000 CALL RESULT (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV,F,F1 12,MATL,DUMMY,DUMMY,BSIGX,BEPSX,HB) GO TO 3000 C C GIVEN FORCES AND CURVATURES C 2020 READ (15,*) FORC, CURV,CURET,RMT,OPRT,HMC WRITE (16,895) CURET,RMT,OPRT,HMC 941 FORMAT ('0',T20,'***GIVEN FORCES N & CURVATURES KAPPA ',/ 1,T10,'N',15X,'KAPPA',10X,'STRAINS E0',8X,'MOMENTS',/, 24(3X,E14.6)) CALL MULT2N (ENAUT,BEND,AI,BS,CS,DS,FORC,CURV) WRITE(16,941) (FORC(I),CURV(I),ENAUT(I),BEND(I),I = 1,3) IF (FORC(1).EQ.0.0.AND.FORC(2).EQ.0.0.AND.FORC(3).EQ.0.0.AND.CURV( 11).EQ.0.0.AND.CURV(2).EQ.0.0.AND.CURV(3).EQ.0.0) GO TO 3000 CALL RESULT (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV,F,F1 12,MATL,DUMMY,DUMMY,BSIGX,BEPSX,HB) GO TO 3000 C C GIVEN FORCES AND MOMENTS C C 2030 READ (15,*) FORC, BEND,CURET,RMT,OPRT,HMC WRITE (16,895) CURET,RMT,OPRT,HMC 940 FORMAT (10F8.0) CALL MULT2N (ENAUT,CURV,AP,BP,CP,DSI,FORC,BEND) WRITE (16,949) 949 FORMAT ('0',T20,'*** GIVEN FORCES AND MOMENTS ***',//) WRITE(16,942) (FORC(I),BEND(I),ENAUT(I),CURV(I),I = 1,3) 942 FORMAT(T10,'N',15X,'MOMENT',10X,'STRAINS',10X,'KAPPA',/,4(3X,E14.6 1)) DO 2038 I = 1,3 IF (FORC(I).NE.0.0 .OR. BEND(I).NE.0.0) GO TO 2039 2038 CONTINUE GO TO 3000 2039 CALL RESULT (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV,F,F1 12,MATL,DUMMY,DUMMY,BSIGX,BEPSX,HB) GO TO 3000 C C GIVEN MIDPLANE STRAINS AND MOMENTS C 2040 READ(15,*) ENAUT,BEND,CURET,RMT,OPRT,HMC WRITE (16,895) CURET,RMT,OPRT,HMC CALL MULT2N (FORC,BEND,AW,BW,CW,DI,ENAUT,BEND) WRITE(16,943) (ENAUT(I),CURV(I),FORC(I),BEND(I),I = 1,3) 943 FORMAT ('0',T20,'***GIVEN MID-PLANE STRAINS AND MOMENTS ***',/, 1T10,'ENAUT',12X,'KAPPA',12X,'FORCES',12X,'MOMENTS',/, 24(3X,E14.6)) IF (ENAUT(1).EQ.0.0.AND.ENAUT(2).EQ.0.0.AND.ENAUT(3).EQ.0.0.AND.BE 1ND(1).EQ.0.0.AND.BEND(2).EQ.0.0.AND.BEND(3).EQ.0.0) GO TO 3000 CALL RESULT (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV,F,F1 12,MATL,DUMMY,DUMMY,BSIGX,BEPSX,HB) C C THERMAL MOISTURE LOADING C 3000 TEMP =-(CURET - OPRT) FABTD = -(CURET-RMT) WRITE (16,917) TEMP,FABTD 917 FORMAT(//////T15,'*********** THE THERMAL PROBLEM *********',/ 1,T15,'CURE TO OPERATING TEMP CHANGE =',F7.1,/,T15,'CURE TEMP TO RO 1OM TEMPERATURE CHANGE =',F7.1,/) C C FABTD & TEMP = 0 IF NO THERMAL EFFECTS C FABTD TEMP IS THE TEMPERATURE DROP DURING CURE AND TEMP IS THE C TEMPERATURE CHANGE T0 TO T FROM CURE TO OPERATING TEMPERATURE C IF (DABS(FABTD).LT.1.0) GO TO 3010 C C RESIDUAL STRESS PROBLEM C TTT = -FABTD T4 = 0.0 C WRITE (16,960) NPLYS C 960 FORMAT (/,5x,'NPLYS = ',I5,/) CALL LTDSTN(TTT,T4,HMC,NMATLS,MATL,NPLYS,ALPAX,ALPAMX,BETAX,BETAM 1X,FTSTRN,FMSTRN,X,Y,XX,YY,NA,KEY(18)) C WRITE (16,960) NPLYS CALL EQFAM(NPLYS,H,QB,NA,RTFORC,RTBEND,HRFORC,HRBEND,FTSTRN,FMSTRN 1) C WRITE (16,960) NPLYS CALL MULT2N(RTENAT,RTCURV,AP,BP,CP,DSI,RTFORC,RTBEND) CALL MULT2N (HRENAT,HRCURV,AP,BP,CP,DSI,HRFORC,HRBEND) C C THE RESIDUAL STRESS PROBLEM AT ROOM TEMPERATURE C C C RETURNS NT, MT, NH, MH FOR TEMPERATURE CHANGE DURING FABRICATION C WRITE (16,950) C 890 FORMAT(////,T20,'***RESIDUAL STRESSES AND STRAINS***',/ 1T20,'***AT ROOM TEMPERATURE ****',/) C THE RESIDUAL STRESS PROBLEM AT ROOM TEMPERATURE C 950 FORMAT (///,T10,'**** EQUIVALENT ROOM TEMP. THERMAL LOADS **'/) WRITE(16,942) (RTFORC(I),RTBEND(I),RTENAT(I),RTCURV(I), I =1,3) WRITE(16,890) C WRITE (16,960) NPLYS CALL RESULT (QB,H,NPLYS,SIGRX,EPSRX,NA,X,Y,XX,YY,KEY,RTENAT,RTCURV 1,F,F12,MATL,FTSTRN,DUMMY,BSIGRX,BEPSRX,HB) 891 FORMAT(////,T20,'***RESIDUAL STRESSES AND STRAINS***',/ 1T15,'***AT OPRT TEMPERATURE =',F6.0,'****',/) C C MOISTURE LOADS & STRESSES AT ROOM TEMPERATURE C WRITE (16,1946) WRITE(16,942) (HRFORC(I),HRBEND(I),HRENAT(I),HRCURV(I),I = 1,3) IF (HMC.EQ.0.0) GO TO 3045 WRITE(16,1890) HMC C WRITE (16,960) NPLYS CALL RESULT (QB,H,NPLYS,SIGMX,EPSMX,NA,X,Y,XX,YY,KEY,HRENAT,HRCURV 1,F,F12,MATL,DUMMY,FMSTRN,BSIGMX,BEPSMX,HB) 1946 FORMAT (///,T20,'*** EQUIVALENT MOISTURE LOADS ***',/) 3045 DO 3050 I = 1,3 TFORC(I) = RTFORC(I) HFORC(I) = HRFORC(I) TBEND(I) = RTBEND(I) 3050 HBEND(I) = HRBEND(I) C C OPERATING TEMPERATURE PROBLEM C 3010 IF (DABS(TEMP).LT.1.0) GO TO 3020 IF (RMT.EQ.OPRT) GO TO 3020 C TTT = OPRT - RMT C WRITE (16,960) NPLYS CALL PRPRTY(X2,X4,Y2,Y4,TTT,HO,AI,AN,AP,AW,BS,BP,BW,CS,CP,CW, 1DI,DSI,DS,DP,A,B,D,Q,QB,S,NA,KEY,H,MATL,EL,ET,ULT,G,U1,U2,U3,U4, 2X,Y,XX,YY,ALPA,ALPAX,ALPABA,BETA,BETAX,BETABA,XT,XC,YT,YC,SS, 3F12,ELM,ETM,GM,XTM,XCM,YTM,YCM,SSM,F12M,ALPAM,BETAM,ETT,ULTT,GT, 4XTT,XCT,YTT,YCT,SST,F12T,ULTM,ALPAT,BETAT,ELT,SB,NANGLS,NMATLS, 5ALPAMX,BETAMX,NPLYS,UTL) TTC = CURET - RMT TTO = OPRT - RMT C WRITE (16,960) NPLYS CALL LTDSTN(TTC,TTO,HMC,NMATLS,MATL,NPLYS,ALPAX,ALPAMX,BETAX,BETAM 1X,FTSTRN,FMSTRN,X,Y,XX,YY,NA,KEY(18)) C WRITE (16,960) NPLYS CALL EQFAM (NPLYS,H,QB,NA,TFORC,TBEND,HFORC,HBEND,FTSTRN,FMSTRN) CALL MULT2N(TENAUT,TCURV,AP,BP,CP,DSI,TFORC,TBEND) CALL MULT2N(HENAUT,HCURV,AP,BP,CP,DSI,HFORC,HBEND) WRITE(16,891) OPRT WRITE (16,952) 952 FORMAT (///,T10,'**** EQUIVALENT OPRT TEMP. THERMAL LOADS ***') WRITE(16,942) (TFORC(I),TBEND(I),TENAUT(I),TCURV(I), I =1,3) C WRITE (16,960) NPLYS CALL RESULT (QB,H,NPLYS,SIGRX,EPSRX,NA,X,Y,XX,YY,KEY,TENAUT,TCURV, 1F,F12,MATL,FTSTRN,DUMMY,BSIGRX,BEPSRX,HB) C C ***** MOISTURE PROBLEM AT OPERATING TEMPERATURE C WRITE (16,1946) WRITE(16,942) (HFORC(I),HBEND(I),HENAUT(I),HCURV(I),I = 1,3) IF (HMC.EQ.0.0) GO TO 3020 WRITE(16,1890) HMC C WRITE (16,960) NPLYS CALL RESULT (QB,H,NPLYS,SIGMX,EPSMX,NA,X,Y,XX,YY,KEY,HENAUT,HCURV, 1F,F12,MATL,DUMMY,FMSTRN,BSIGMX,BEPSMX,HB) C RETURN NT,MT,NH,MH FOR OPERATING TEMPERATURE C HMC = 0 IF ZERO MOISTURE EFFECT C C CALCULATE N & M DUE TO MECHANICAL, THERMAL & MOISTURE LOADS C 3020 DO 2031 I = 1,3 TOFORC(I) = FORC(I) + TFORC(I)+ HFORC(I) 2031 TOBEND(I) = BEND(I) + TBEND(I) + HBEND(I) C C CALCULATE MIDPLANE STRAINS & CURVATURES C C C THERMAL RESULTS AT OPERATING TEMPERATURE C C C EQUIVALENT MOISTURE LOADS AT OPERATING TEMP(I THINK) C 947 FORMAT(///,T20,'*** MECHANICAL, THERMAL & MOISTURE LOADING ***',/ 1) CALL MULT2N (TOENUT,TOCURV,AP,BP,CP,DSI,TOFORC,TOBEND) C C THE COMBINED LOADS PROBLEM C WRITE (16,947) WRITE(16,942) (TOFORC(I),TOBEND(I),TOENUT(I),TOCURV(I),I = 1,3) 1890 FORMAT(////,T20,'***MOISTURE STRESSES AND STRAINS***',/ 1T20,'*** MOISTURE CONTENT = ',F7.3,'% ***',/) 500 IF (KEY(LL+2).EQ.0) GO TO 4000 WRITE(16,889) OPRT 889 FORMAT(////,T10,'*** TOTAL STRESSES & STRAINS DUE TO COMBINED ***' 1,/10X,'*** MECHANICAL AND THERMAL LOADING***',/,T10,'*** OPERATING 2TEMP = ',F5.1,' ***',/) C WRITE (16,960) NPLYS CALL RESULT (QB,H,NPLYS,SIGXT,EPSXT,NA,X,Y,XX,YY,KEY,TOENUT,TOCURV 1,F,F12,MATL,FTSTRN,FMSTRN,BSIGXT,BEPSXT,HB) 4000 CONTINUE GO TO 1111 1000 STOP END C C ***** SUBROUTINE INVERT ***** C SUBROUTINE INVERT(A,B) C INVERTS A 3X3 MATRIX `A INTO A 3X3 MATRIX `B = `A INVERSE C DOUBLE PRECISION A,B,D IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(9),B(9) D = A(1)*(A(5)*A(9) - A(8)*A(6)) - A(4)*(A(2)*A(9) 1 - A(3)*A(8)) + A(7)*(A(2)*A(6) - A(3)*A(5)) B(1) = (A(5)*A(9) - A(8)*A(6))/D B(2) = -(A(4)*A(9) - A(6)*A(7))/D B(3) = (A(4)*A(8) - A(5)*A(7))/D B(4) = -(A(2)*A(9) - A(3)*A(8))/D B(5) = (A(1)*A(9) - A(3)*A(7))/D B(6) = -(A(1)*A(8) - A(2)*A(7))/D B(7) = (A(2)*A(6) - A(3)*A(5))/D B(8) = -(A(1)*A(6) - A(3)*A(4))/D B(9) = (A(1)*A(5) - A(2)*A(4))/D RETURN END C C ***** SUBROUTINE EQFAM ***** C SUBROUTINE EQFAM (NPLYS,H,QB,NA,TFORC,TBEND,HFORC,HBEND,FTSTRN,FMS 1TRN) C C CALCULATES EQUIVALENT THERMAL & HYGROSCOPIC FORCE & MOMENT C FOR UNIFORM DELTA T IMPLICIT REAL*8(A-H,O-Z) DIMENSION H(1),NA(1),QB(1),TFORC(1),TBEND(1),HFORC(1),FTSTRN(1),FM 1STRN(1),HBEND(1),DELTA(3),HBETA(3) DO 20 I = 1,3 TFORC(I) = 0. TBEND(I) = 0. HFORC(I) = 0. 20 HBEND(I) = 0. DO 50 I =1, NPLYS DH = H(I+1) - H(I) DH2 = (H(I+1) ** 2 - H(I)**2)/2. L = 6 * NA(I) - 6 J = 3 * (I - 1) DELTA(1) = (QB(L+1)*FTSTRN(J+1)+QB(L+2)*FTSTRN(J+2) 1 + QB(L+3) * FTSTRN(J+3)) DELTA(2) = (QB(L+2)*FTSTRN(J+1) + QB(L+4)*FTSTRN(J+2) 1 + QB(L+5)*FTSTRN(J+3)) DELTA(3) = (QB(L+3)*FTSTRN(J+1) + QB(L+5)*FTSTRN(J+2) 1 + QB(L+6)*FTSTRN(J+3)) HBETA(1) = (QB(L+1)*FMSTRN(J+1)+QB(L+2)*FMSTRN(J+2) 1 + QB(L+3) * FMSTRN(J+3)) HBETA(2) = (QB(L+2)*FMSTRN(J+1) + QB(L+4)*FMSTRN(J+2) 1 + QB(L+5)*FMSTRN(J+3)) HBETA(3) = (QB(L+3)*FMSTRN(J+1) + QB(L+5)*FMSTRN(J+2) 1 + QB(L+6)*FMSTRN(J+3)) DO 50 K = 1,3 TFORC(K) = TFORC(K) + DELTA(K) * DH TBEND(K) = TBEND(K) + DELTA(K) * DH2 HFORC(K) = HFORC(K) + HBETA(K) * DH HBEND(K) = HBEND(K) + HBETA(K) * DH2 50 CONTINUE RETURN END C C ***** SUBROUTINE SANDS ***** C SUBROUTINE SANDS (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV 1,FTSTRN,FMSTRN,BSIGX,BEPSX,HB) C DOUBLE PRECISION X,Y,XX,YY IMPLICIT REAL*8(A-H,O-Z) C CALCULATES STRESS & STRAIN GIVEN MID-PLANE STRAINS C AND CURVATURES, TEMPERATURE AND MOISTURE C DIMENSION QB(120),H(21),SIGX(60),EPSX(60),ENAUT(3),CURV(3), 1BSIG1(60),BEPS1(60) DIMENSION NA(20),SIG1(60),EPS1(60),X(10),Y(10),XX(10),YY(10), 1KEY(25),FTSTRN(60),FMSTRN(60),BSIGX(60),BEPSX(60),HB(21) DO 100 I = 1, NPLYS C DH2 = (H(I+1)-H(I))/2. + H(I) L = 3 * I - 3 DO 50 J = 1,3 BEPSX(L+J) = ENAUT(J) + HB(I)*CURV(J) - FTSTRN(L+J) - FMSTRN(L+J) 50 EPSX(L+J) = ENAUT(J) + H(I)*CURV(J) - FTSTRN(L+J) - FMSTRN(L+J) K = NA(I)*6 - 6 SIGX(L+1) = QB(K+1)*EPSX(L+1)+QB(K+2)*EPSX(L+2)+QB(K+3)*EPSX(L+3) SIGX(L+2) =QB(K+2)*EPSX(L+1)+QB(K+4)*EPSX(L+2)+QB(K+5)*EPSX(L+3) SIGX(L+3)=QB(K+3)*EPSX(L+1)+QB(K+5)*EPSX(L+2)+QB(K+6)*EPSX(L+3) BSIGX(L+1) = QB(K+1)*BEPSX(L+1)+QB(K+2)*BEPSX(L+2) 1+QB(K+3)*BEPSX(L+3) BSIGX(L+2) =QB(K+2)*BEPSX(L+1)+QB(K+4)*BEPSX(L+2) 1+QB(K+5)*BEPSX(L+3) BSIGX(L+3)=QB(K+3)*BEPSX(L+1)+QB(K+5)*BEPSX(L+2) 1+QB(K+6)*BEPSX(L+3) DO 75 J = 1,3 BEPSX(L+J) = ENAUT(J) + HB(I)*CURV(J) 75 EPSX(L+J) = ENAUT(J) + H(I)*CURV(J) 100 CONTINUE CALL RITESS (NPLYS,SIGX,EPSX,BSIGX,BEPSX,H,HB,1) C the last variable in the CALL RITESS statement C is used for column labeling only C1 = 0.5 C2 = 2.0 C3 = 1.0 CALL ROTATE (EPS1,EPSX,X,Y,XX,YY,NPLYS,NA,C1,C2) CALL ROTATE (SIG1,SIGX,X,Y,XX,YY,NPLYS,NA,C3,C3) CALL ROTATE (BEPS1,BEPSX,X,Y,XX,YY,NPLYS,NA,C1,C2) CALL ROTATE (BSIG1,BSIGX,X,Y,XX,YY,NPLYS,NA,C3,C3) WRITE (19,912) 125 WRITE (18,911) DO 200 I = 1, NPLYS L = 3 * I - 3 WRITE (18,903)I,H(I),SIGX(L+1),SIGX(L+2),SIGX(L+3),SIG1(L+1), 1SIG1(L+2),SIG1(L+3) WRITE (18,903) I, HB(I),BSIGX(L+1),BSIGX(L+2),BSIGX(L+3), 1BSIG1(L+1),BSIG1(L+2),BSIG1(L+3) WRITE(19,903)I,H(I),EPSX(L+1),EPSX(L+2),EPSX(L+3),EPS1(L+1), 1EPS1(L+2),EPS1(L+3) WRITE (19,903) I, HB(I),BEPSX(L+1),BEPSX(L+2),BEPSX(L+3), 1BEPS1(L+1),BEPS1(L+2),BEPS1(L+3) 903 FORMAT (2X,I5,3X,E14.6,6(2X,E14.6)) 910 FORMAT (//,10X,'STRESS & STRAIN IN NATURAL COORDINATES',/) 200 CONTINUE WRITE (16,910) CALL RITESS (NPLYS,SIG1,EPS1,BSIG1,BEPS1,H,HB,4) 911 FORMAT (//,10X,'THROUGH THICKNESS STRESS DISTRIBUTION',/, 15X,'PLY',8X,'H',14X,'SIGX',13X,'SIGY',10X,'TAUXY',12X, 2'SIG1',12X,'SIG2',12X,'TAU12',/'*DATA') 912 FORMAT (//,10X,'THROUGH THICKNESS STRAIN DISTRIBUTION',/, 15X,'PLY',8X,'H',14X,'EPSX',13X,'EPSY',10X,'EPSXY',12X, 2'EPS1',12X,'EPS2',12X,'GAM12',/'*DATA') RETURN END C C ***** SUBROUTINE RITESS ***** SUBROUTINE RITESS (NPLYS,SIGX,EPSX,BSIGX,BEPSX,H,HB,K) C WRITES STRESS AND STRAIN VALUES FOR ALL PLIES C AT TOP & BOTTOM OF PLY C the last variable K in the SUBROUTINE RITESS statement C is used for column labeling only IMPLICIT REAL*8(A-H,O-Z) DIMENSION SIGX(60),EPSX(60),BSIGX(60),BEPSX(60), 1N(6),H(21),HB(21) DATA N(1),N(2),N(3),N(4),N(5),N(6)/'X','Y','XY','1','2','12'/ K1 = K + 2 WRITE (16,902) (N(I), I = K,K1) 902 FORMAT (//,30X,'***PLY TOP & BOTTOM STRESSES***',//,5X,'PLY', 18X,'H',14X,'SIG',A1,13X,'SIG',A1,10X,'TAU',A2,/) L = 0 903 FORMAT (2X,I5,3X,E14.6,3(2X,E14.6)) 905 FORMAT (2X,I5,3X,E14.6,3(2X,E14.6),/) DO 150 I = 1, NPLYS WRITE(16,903) I,H(I), SIGX(L+1),SIGX(L+2),SIGX(L+3) WRITE(16,905) I,HB(I),BSIGX(L+1),BSIGX(L+2),BSIGX(L+3) 150 L = L + 3 WRITE (16,904) (N(I), I = K,K1) 904 FORMAT (//30X,'***PLY TOP & BOTTOM STRAINS***', 1//,5X,'PLY',8X,'H',14X,'EPS',A1, 113X,'EPS',A1,10X,'GAMMA',A2,/) L = 0 DO 200 I = 1,NPLYS WRITE(16,903) I,H(I),EPSX(L+1),EPSX(L+2),EPSX(L+3) WRITE(16,905) I,HB(I),BEPSX(L+1),BEPSX(L+2),BEPSX(L+3) 200 L = L + 3 RETURN END C C C ***** SUBROUTINE ABD ***** C C CALCULATES ALL FORMS OF A, B, D MATRICES FOR USE IN C CONSTITUTIVE EQUATIONS C SUBROUTINE ABD(A,B,D,AI,AP,AW,BS,BP,BW,CS,CP,CW,DS,DSI,DI) C DOUBLE PRECISION A,B,D,AI,AP,AW,BS,BP,BW,CS,CP,CW,DS,DSI,DI C DOUBLE PRECISION G IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(1),B(1),D(1),AI(1),AP(1),AW(1),BS(1),BP(1),BW(1), 1CS(1),CP(1),CW(1),DS(1),DSI(1),DI(1),G(9) DO 50 I = 1,9 AI(I) = -AI(I) 50 CONTINUE CALL MULT4N (AI,B,BS) DO 75 I = 1,9 AI(I) = - AI(I) 75 CONTINUE CALL MULT4N (B,AI,CS) CALL MULT4N (CS,B,G) DO 100 I = 1,9 100 DS(I) = D(I) - G(I) CALL INVERT (DS,DSI) DO 150 I = 1,9 150 DSI(I) = - DSI(I) CALL MULT4N(DSI,CS,CP) DO 175 I = 1,9 175 DSI(I) = - DSI(I) CALL MULT4N(BS,DSI,BP) CALL MULT4N(BS,CP,G) DO 200 I = 1,9 200 AP(I) = AI(I) + G(I) CALL INVERT(D,DI) CALL MULT4N(B,DI,BW) DO 250 I = 1,9 250 DI(I) = - DI(I) CALL MULT4N(DI,B,CW) CALL MULT4N(B,CW,G) DO 275 I = 1,9 275 DI(I) = - DI(I) DO 300 I = 1,9 300 AW(I) = A(I) - G(I) RETURN END C C ***** SUBROUTINE MULT2N ***** C SUBROUTINE MULT2N(A,B,C,D,E,F,X,Y) C MULTIPLIES A 9X9 MATRIX INTO A 9X1 COLUMN MATRIX C DOUBLE PRECISION C,D,E,F IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3),B(3),C(9),D(9),E(9),F(9),X(3),Y(3) L = 0 DO 100 J = 1,3 A(J) = 0. 100 B(J) = 0. DO 200 J = 1,3 A(J) = A(J) + C(L+1)*X(1)+C(L+4)*X(2)+C(L+7)*X(3) 1 + D(L+1)*Y(1) + D(L+4)*Y(2) + D(L+7)*Y(3) B(J) = B(J) + E(L+1)*X(1)+E(L+4)*X(2)+E(L+7)*X(3) 1 + F(L+1)*Y(1) + F(L+4)*Y(2) + F(L+7)*Y(3) 200 L = L + 1 RETURN END C C ***** SUBROUTINE MULT4N ***** C SUBROUTINE MULT4N(B,C,A) C MULTIPLIES A 3X3 MATRIX INTO A 3X3 MATRIX C DOUBLE PRECISION A,B,C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(9),B(9),C(9) DO 50 I = 1,9 50 A(I) = 0. L = 1 DO 100 I = 1,9 IF (I.LT.4) J = 1 IF (I.GT.3.AND.I.LT.7) J = 4 IF (I.GT.6) J = 7 A(I) = A(I) + B(L)*C(J)+B(L+3)*C(J+1)+B(L+6)*C(J+2) L = L + 1 IF (L.GT.3) L = 1 100 CONTINUE RETURN END C C ***** SUBROUTINE STRESS ***** C SUBROUTINE STRESS(E,S,Q,NPLYS,NA) C CALCULATES STRESSES ASSOCIATED WITH STRAINS E FOR EACH PLY C OF STIFFNESS Q IMPLICIT REAL*8(A-H,O-Z) DIMENSION E(60),S(60),Q(1),NA(1) DO 100 I = 1,NPLYS K = 3*I - 3 L = 6 * NA(I)-6 S(K+1) = Q(L+1)*E(K+1) + Q(L+2)*E(K+2) + Q(L+3)*E(K+3) S(K+2) = Q(L+2)*E(K+1)+Q(L+4)*E(K+2) + Q(L+5)*E(K+3) S(K+3) = Q(L+3)*E(K+1) + Q(L+5)*E(K+2) + Q(L+6)*E(K+3) 100 CONTINUE RETURN END C C ***** SUBROUTINE ROTATE ***** C SUBROUTINE ROTATE(A,B,X,Y,XX,YY,NPLYS,NA,C1,C2) C DOES TRANSFORMATION A = T B FOR EACH PLY OF LAMNATE C DOUBLE PRECISION X,Y,XX,YY IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(60),B(60),X(1),Y(1),XX(1),YY(1),NA(1) DO 100 I = 1,NPLYS K = 3*I - 3 J = NA(I) XY = X(J)*Y(J) A(K+1) = YY(J)*B(K+1) + XX(J)*B(K+2) + 2.0*C1*XY*B(K+3) A(K+2) = XX(J)*B(K+1) + YY(J)*B(K+2) -2.0*C1* XY*B(K+3) A(K+3) = -C2*XY*B(K+1) + C2*XY*B(K+2) + (YY(J)-XX(J))*B(K+3) 100 CONTINUE RETURN END C C ***** SUBROUTINE LAMINA ***** C SUBROUTINE LAMINA (QB,SB,ANG,X,Y,XX,YY,KEY,EPSX,SIGX,EPS1,SIG1,NPL 1YS,NA) C CALCULATES STRESS OR STRAIN FOR OFF-AXIS LAMINA C DOUBLE PRECISION X,Y,XX,YY,ANG IMPLICIT REAL*8(A-H,O-Z) DIMENSION QB(1),SB(1),ANG(1),X(1),Y(1),XX(1),YY(1),KEY(1),EPSX(1), 1EPS1(1),SIGX(1),SIG1(1),NA(1) C KEY(5) = NO OF LOAD CASES C KEY(6) THRU KEY(6)+KEY(5) INDICATE LOAD CASE TYPES C IF KEY(6) = 1 GIVEN STRAINS C = 2 GIVEN STRESSES C = 3 DETERMINE AXIAL FAILURE STRESS AS FUNCTION OF C THETA FOR TSAI-WU CRITERIA C = 4 DETERMINE SHEAR FAILURE STRESS AS FUNCTION OF C THETA FOR TSAI-WU CRITERIA C = 5 PLOT TSAI-WU FAILURE SURFACE IN 1,2 STRESS SPACE K1 = KEY(5) DO 100 I = 1,K1 J = KEY(5+I) GO TO (1000,2000,3000,4000,5000), J 1000 READ(15,*) (EPSX(K),K=1,3) WRITE(16,903) (EPSX(K),K=1,3) 903 FORMAT (4X,'***',3(E14.6,3X)) CALL STRESS(EPSX,SIGX,QB,NPLYS,NA) WRITE(16,900) 900 FORMAT ('0',T20,'*** GIVEN LAMINA STRAINS ***',/) CALL RITESS (NPLYS,SIGX,EPSX,1) IF (KEY(3).EQ.0) GO TO 100 CALL ROTATE(EPS1,EPSX,X,Y,XX,YY,NPLYS,NA,.5,2.) CALL ROTATE(SIG1,SIGX,X,Y,XX,YY,NPLYS,NA,1.,1.) CALL RITESS(NPLYS,SIG1,EPS1,4) GO TO 100 2000 READ(15,*) (SIGX(K),K=1,3) WRITE(16,903) (SIGX(K),K=1,3) CALL STRESS(SIGX,EPSX,SB,NPLYS,NA) WRITE(16,901) CALL RITESS (NPLYS,SIGX,EPSX,1) 901 FORMAT ('0',T20,'*** GIVEN LAMINA STRESSES ***',/) 902 FORMAT (3F9.0) IF (KEY(3).EQ.0) GO TO 100 CALL ROTATE(EPS1,EPSX,X,Y,XX,YY,NPLYS,NA,.5,2.) CALL ROTATE (SIG1,SIGX,X,Y,XX,YY,NPLYS,NA,1.,1.) CALL RITESS (NPLYS,SIG1,EPS1,4) 3000 GO TO 100 4000 GO TO 100 5000 GO TO 100 100 CONTINUE RETURN END C C ******* SUBROUTINE INTWCURV C C SUBROUTINE CALCULATES INTERFACE MOMENT MZ AND INTERFACE FORCES C FZX AND FYZ, AT SPECIFIED LOCATIONS, IN THE PRESENCE OF C CURVATURE. C SUBROUTINE INTWCURV (NPLYS,CURV,QB,H,KEY,NA) IMPLICIT REAL*8(A-H,O-Z) DIMENSION H(1), KEY(21), CURV(1),QB(1),NA(1) WRITE (16,901) 901 FORMAT (//T20,'***INTERLAMINAR RESULTS - WITH CURVATURE***',/ 19X,'Z', 7X,'SIGMA Z MOM',5X,'TAUZY FORCE',5X,'TAUZX FORCE',/) WRITE (17,904) FZX = 0.0 FYZ = 0.0 ZMZ = 0.0 WRITE (16,902) H(NPLYS+1), ZMZ, FYZ, FZX WRITE (17,902) H(NPLYS+1), ZMZ, FYZ, FZX NINC = key(25) C KEY(25) IS THE NUMBER OF INCREMENTS W/I LAYER DO 30 I = NPLYS, 1, -1 DZ = (H(I+1) - H(I))/NINC J = 6 * NA(I) -6 QD2 = -(QB(J+2)*CURV(1) + QB(J+4)*CURV(2) + QB(J+5)* 1CURV(3)) QD6 = -(QB(J+3)*CURV(1) + QB(J+5)*CURV(2) + QB(J+6)* 1CURV(3)) DO 40 LL = 1, NINC zj = h(I+1) - DZ * (LL-1) zj1 = zj - DZ zstar = zj1 zsq = (zj*zj - zj1*zj1)/2.0 FZX = FZX + QD6 * zsq FYZ = FYZ + QD2 * zsq zcube = ((h(I+1)**3 - zj1**3)/3.0 - zj1*(h(I+1)**2 - 1zj1*zj1)/2.00) ZMZ = qd2 * zcube IF (I.EQ.NPLYS) go to 60 do 50 JJ = I+1, NPLYS K = 6 * NA(JJ) -6 QD2M = -(QB(K+2)*CURV(1) + QB(K+4)*CURV(2) + QB(K+5)* 1CURV(3)) zjj = h(jj+1) zjj1 = H(jj) zcube = ((zjj**3 - zjj1**3)/3.0 - zstar*(zjj*zjj - 1zjj1*zjj1)/2.00) ZMZ = ZMZ + qd2m * zcube 50 continue 60 WRITE (16,902) zstar, ZMZ, FYZ, FZX WRITE (17,902) zstar, ZMZ, FYZ, FZX 40 continue 30 CONTINUE 902 FORMAT (5X,F7.5,2X,E14.6,2X,E14.6,2X,E14.6) 904 FORMAT (5x,'INTERLAMINAR MOMENTS & FORCES',/, a8x,'CURVATURE LOADING',/,8X,'Z',13x,'Mz', a10x,'Fyz',15x,'Fzx',/,'*DATA') RETURN END C C ***** SUBROUTINE MOMENT ***** C SUBROUTINE MOMENT (NPLYS,SIGX,BSIGX,H,KEY) C SUBROUTINE CALCULATES INTERFACE MOMENT MZ AND INTERFACE FORCES C XF AND YF AT EACH INTERFACE. IMPLICIT REAL*8(A-H,O-Z) DIMENSION SIGX(1), H(1),YM(20),XF(20),YF(20),KEY(25),BSIGX(60) WRITE (16,901) 901 FORMAT (//T20,'***INTERLAMINAR STRESS RESULT***',/5X,'INTERFACE', 14X,'SIGMA Z MOM',8X,'TAUZY FORCE',9X,'TAUZX FORCE',/) WRITE (17,904) DO 10 I = 1,NPLYS+1 YF(I) = 0. XF(I) = 0. 10 YM(I) = 0. WRITE (16,902) 1, YM(1),YF(1),XF(1) DO 30 J = 1, NPLYS L = 3*(J-1) DH = H(J+1)-H(J) XF(J+1) = XF(J) - (SIGX(L+3)+BSIGX(L+3))*DH/2 YF(J+1) = YF(J) - (SIGX(L+2)+BSIGX(L+2))*DH/2 DO 50 I = 1, J L = 3*(I-1) DH = H(I+1)-H(I) YM(J+1) =YM(J+1)-(SIGX(L+2)+BSIGX(L+2))*(DH/2) 1*((H(I+1)+H(I))/2.-H(J+1)) 50 continue 29 WRITE (16,902) J+1, YM(J+1),YF(J+1),XF(J+1) 30 CONTINUE IF (KEY(24) .GT. O ) THEN write (16,905) 905 Format (/,5x,'Through-Thickness Interlaminar Forces', 1' & Moment'/,8X,'Z',10x,'Mz',12x,'Fyz',10x,'Fzx',/) z = -H(1) WRITE (17,903) z, 0.0, 0.0, 0.0 WRITE (16,903) z, 0.0, 0.0, 0.0 NINC = KEY(25) DO 100 I = 1, NPLYS DZ = (H(I+1) - H(I))/NINC DO 200 M = 1, NINC ZMZ = 0.0 Z = H(I) + DZ * M LL = 3 * (I-1) C IF (I-1 .EQ. 0) THEN C LLL = 0 C ELSE C LLL = 1 C END IF FZX = XF(I) - SIGX(LL+3)*(Z - H(I)) FYZ = YF(I) - SIGX(LL+2)*(Z - H(I)) IF (I .EQ. 1) GO TO 61 JJ = I - 1 DO 60 J = 1, JJ L = 3*(J-1) DH = H(J+1)-H(J) 60 ZMZ = ZMZ - SIGX(L+2)* DH *((H(J+1)+H(J))/2 - Z) 61 ZMZ = ZMZ - SIGX(LL+2)*(Z - H(I))*((Z + H(I))/2 - Z) zt = -z WRITE (17, 903) ZT, ZMZ, FYZ, FZX WRITE (16, 903) ZT, ZMZ, FYZ, FZX 200 CONTINUE 100 CONTINUE ELSE END IF 902 FORMAT (5X,I5,6X,E14.6,6X,E14.6,6X,E14.6) 903 FORMAT (5X,F7.5,3(2X,E14.6)) 904 FORMAT (5x,'INTERLAMINAR MOMENTS & FORCES',/,8x,'Z',13x,'Mz', a10x,'Fyz',15x,'Fzx',/,'*DATA') RETURN END C C ***** SUBROUTINE SAC3D ***** C SUBROUTINE SAC3D (KEY,H,AI,QB,NANGLS,EL,ET,ULT,G,U23,NPLYS,ANG,NA, 1HO,ALPABA,ALPAX,ALPA,MATL) C STIFFNESS & COMPLIANCE FOR TRANSVERSELY ISOTROPIC 3-D MATERIAL C PROGRAM NAME -- 3DSAC -- C DOUBLE PRECISION ANGR,T1,T2,XM,YN,XX,YY,XY,S,C,SB,CB,A,B,DSIN,DCOS C 1,T1I,T2I IMPLICIT REAL*8(A-H,O-Z) DIMENSION S(6,6),C(6,6),SB(6,6,10),CB(6,6,10),ANG(10),H(21),KEY(25 1),AI(9),QB(120),EL(5),ET(5),ULT(5),G(5),U23(5),NA(20),ALPABA(1), 2ALPAX(1),ALPA(1),MATL(1), 1T1(6,6),T2(6,6),A(6,6),B(6,6),T1I(6,6),T2I(6,6),P(6,6) DO 10 I=1,6 DO 10 J = 1,6 S(I,J) = 0. C(I,J) = 0. T1(I,J) = 0. T2(I,J) = 0. P(I,J) = 0. 10 CONTINUE S(1,1) = 1./EL(1) S(1,2) = -ULT(1)/EL(1) S(1,3) = S(1,2) S(2,2) = 1./ET(1) S(2,3) = -U23(1)/ET(1) S(3,3) = S(2,2) S(4,4) = 2.*(S(2,2) - S(2,3)) S(5,5) = 1./G(1) 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) IF (KEY(1).EQ.0) GO TO 20 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,'COMPLIANCE MATRIX S IN 1-2 COORDNATE SYSTEM',/ 1,6(6(3X,E10.3)/)) 905 FORMAT (//15X,'STIFFNESS MATRIX C IN 1-2 COORDINATE SYSTEM',/ 1,6(6(3X,E10.3)/)) C C BEGIN DO LOOP FOR EACH ANGLE C 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 IF (KEY(1).EQ.0) GO TO 200 BETA = ANG(L)*180/ACOS(-1.0) WRITE(16,906) BETA,((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,'STIFFNESS & COMPLIANCE FOR THETA = ',F6.2,'DEGREES', 1//20X,'****STIFFNESS CBAR****',/,6(6(3X,E10.3)/)) 907 FORMAT (//20X,'****COMPLIANCE SBAR****',/,6(6(3X,E10.3)/)) 200 CONTINUE CALL NUXZ (AI,NPLYS,NA,QB,SB,H,HO) CALL ALPAZ (NPLYS,ALPABA,ALPAX,QB,NA,SB,H,ALPA,HO,MATL) RETURN END C C ************** SUBROUTINE ALPAZ *********************** C SUBROUTINE ALPAZ (NPLYS,ALPABA,ALPAX,QB,NA,SB,H,ALPA,HO,MATL) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ALPABA(1),ALPAX(1),NA(1),QB(120),H(21),SB(6,6,10), 1ALPA(1),MATL(1) ALPAZB = 0.0 DO 200 I = 1, NPLYS L = 6 * (NA(I)-1) DH = H(I+1) - H(I) K = NA(I) L2 = 2 * (MATL(K) - 1) LL = 3 * (NA(I) - 1 ) DALX = ALPABA(1) - ALPAX(LL+1) DALY = ALPABA(2) - ALPAX(LL+2) DALXY = ALPABA(3) - ALPAX(LL+3) ALPAZB = ALPAZB +(SB(3,1,K)*(QB(L+1)*DALX+QB(L+2)*DALY+QB(L+3)* 1DALXY) + SB(3,2,K)*(QB(L+2)*DALX+QB(L+4)*DALY+QB(L+5)*DALXY) 2+ SB(3,6,K)*(QB(L+3)*DALX + QB(L+5)*DALY + QB(L+6)*DALXY) 3+ ALPA(L2+2))*DH C WRITE(16,905) ALPA(L2+2), ALPAZB C 905 FORMAT (/,10X,2(E14.6,5X)) 200 CONTINUE ALPAZB = ALPAZB / (2. * HO ) WRITE(16,900) ALPAZB 900 FORMAT (//,10X,'LAMINATE CTE ALPAZBAR = ',E14.6) RETURN END C C ************** SUBROUTINE NUXZ *********************** C SUBROUTINE NUXZ (AI,NPLYS,NA,QB,SB,H,HO) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AI(9),NA(1),QB(120),H(21),SB(6,6,10),SQ(3) DO 100 I = 1, 3 100 SQ(I) = 0.0 DO 200 I = 1, NPLYS L = 6 * (NA(I)-1) DH = H(I+1) - H(I) K = NA(I) SQ(1) = SQ(1) +(SB(1,3,K)*QB(L+1)+SB(2,3,K)*QB(L+2)+SB(3,6,K)*QB(L 1+3)) * DH SQ(2) = SQ(2) +(SB(1,3,K)*QB(L+2)+SB(2,3,K)*QB(L+4)+SB(3,6,K)*QB(L 1+5)) * DH SQ(3) = SQ(3) +(SB(1,3,K)*QB(L+3)+SB(2,3,K)*QB(L+5)+SB(3,6,K)*QB(L 1+6)) * DH 200 CONTINUE ZNUXZ = -(SQ(1) + SQ(2)*AI(2)/AI(1) + SQ(3)*AI(3)/AI(1))/(2*HO) WRITE(16,900) ZNUXZ 900 FORMAT (//,10X,'LAMINATE POISONS RATIO NUXZ = ',E14.6) ZNUYZ = -(SQ(1)*AI(4)/AI(5) + SQ(2) + SQ(3)*AI(6)/AI(5))/(2*HO) WRITE(16,901) ZNUYZ 901 FORMAT (/,10X,'LAMINATE POISONS RATIO NUYZ = ',E14.6,//) RETURN END C C ***** SUBROUTINE TSAIWU ***** C SUBROUTINE TSAIWU (F,F12,SIG1,BSIG1,NPLYS,NA,MATL,X,Y,XX,YY) C C *** CALCULATES TSAI-WU FUNCTION FOR EACH PLY *** C C DOUBLE PRECISION X,Y,XX,YY IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(1),F12(1),SIG1(1),TW(20),NA(1),MATL(1),BTW(20) DIMENSION X(1),Y(1),XX(1),YY(1),SIGP(60),BSIG1(60),BSIGP(60) CALL ROTATE (SIGP,SIG1,X,Y,XX,YY,NPLYS,NA,0.1D01,0.1D01) CALL ROTATE (BSIGP,BSIG1,X,Y,XX,YY,NPLYS,NA,0.1D01,0.1D01) WRITE (16,900) K = 0 DO 100 I = 1, NPLYS J = NA(I) L = (MATL(J) - 1 ) * 5 TW(I) = F(L+1)*SIGP(K+1) + F(L+2)*SIGP(K+2) + F(L+3)*SIGP(K+1)* 1 SIGP(K+1) + F(L+4)*SIGP(K+2)*SIGP(K+2) + F(L+5)*SIGP(K+3) 2 * SIGP(K+3) + 2. * F12(L+1)*SIGP(K+1)*SIGP(K+2) BTW(I) =F(L+1)*BSIGP(K+1)+F(L+2)*BSIGP(K+2)+F(L+3)*BSIGP(K+1)* 1BSIGP(K+1)+F(L+4)*BSIGP(K+2)*BSIGP(K+2)+F(L+5)*BSIGP(K+3) 2 * BSIGP(K+3)+2.*F12(L+1)*BSIGP(K+1)*BSIGP(K+2) K = K + 3 100 WRITE (16,901) I,TW(I),BTW(I) 900 FORMAT (///T20,'*** TSAI-WU FUNCTION FOR EACH PLY ***',// 1T20,'PLY NO',7X,'TOP PLY T-W ',9X,'BOTTOM PLY T-W',/) 901 FORMAT (T23,I2,7X,E14.6,7X,E14.6) RETURN END C C ***** SUBROUTINE FAILUR ***** C SUBROUTINE FAILUR (NPLYS,NA,QB,YY,XX,X,Y,A,B,Z,P,R) C *** CALCULATES COEFFICIENTS FOR INSERTION INTO FAILURE C *** CRITERIA FOR PARAMETRIC STUDY IN IN-PLANE OR MOMENT C *** LOAD SPACE. P IS THE VECTOR OF COEFFICIENTS AND R C *** IS THE VECTOR OF LOAD RATIOS C C DOUBLE PRECISION YY,XX,X,Y,A,B IMPLICIT REAL*8(A-H,O-Z) DIMENSION NA(1),QB(1),YY(1),XX(1),X(1),Y(1),A(1),B(1),Z(1), 1Q(9),T(9),E(9),F(9),G(9),P(9),R(3) DO 100 I = 1, NPLYS K = 6 * NA(I) - 6 Q(1) = QB(K+1) Q(2) = QB(K+2) Q(3) = QB(K+3) Q(4) = Q(2) Q(5) = QB(K+4) Q(6) = QB(K+5) Q(7) = Q(3) Q(8) = Q(6) Q(9) = QB(K+6) M = NA(I) T(1) = YY(M) T(2) = XX(M) T(3) = 2. * X(M) * Y(M) T(4) = T(2) T(5) = T(1) T(6) = -T(3) T(7) = T(6)/2. T(8) = -T(7) T(9) = YY(M)- XX(M) CALL MULT4N(T,Q,E) DO 200 J = 1,9 F(J) = A(J) + Z(I) * B(J) 200 CONTINUE CALL MULT4N (E,F,G) N = (I-1) * 3 DO 300 L = 1,3 P(N+L) = G(J) * R(1) + G(J+1)*R(2) + G(J+2)*R(3) 300 J = J + 3 100 CONTINUE RETURN END C C ***** SUBROUTINE TRNSFR ***** C SUBROUTINE TRNSFR (NANGLS,MATL,XX,YY,X,Y,ALPA,ALPAX) C C TRANSFORMS ALPA 1-2 TO ALPA X-Y C C DOUBLE PRECISION XX,YY,X,Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION MATL(1),XX(1),YY(1),X(1),Y(1),ALPA(1),ALPAX(1) 175 M = 0 DO 200 I = 1, NANGLS N = 2 * MATL(I) - 1 ALPAX(M+1) = YY(I)*ALPA(N) + XX(I)*ALPA(N+1) ALPAX(M+2) = XX(I)*ALPA(N) + YY(I)*ALPA(N+1) ALPAX(M+3) = -2. * X(I)*Y(I) * (ALPA(N+1)-ALPA(N)) M = M + 3 200 CONTINUE RETURN END C C ***** SUBROUTINE RESULT ***** C SUBROUTINE RESULT (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CUR 1V,F,F12,MATL,FTSTRN,FMSTRN,BSIGX,BEPSX,HB) C DOUBLE PRECISION X,Y,XX,YY IMPLICIT REAL*8(A-H,O-Z) C C CALL SUBROUTINESS TO CALCULATE STRESSES, Y MOMENT & TSAI-WU C DIMENSION QB(1),H(1),SIGX(1),EPSX(1),NA(1),X(1),Y(1),XX(1),YY(1), 1KEY(1),ENAUT(1),CURV(1),F(1),F12(1),MATL(1),FTSTRN(1),FMSTRN(1), 2BSIGX(1),BEPSX(1),HB(1) CALL SANDS (QB,H,NPLYS,SIGX,EPSX,NA,X,Y,XX,YY,KEY,ENAUT,CURV, 1FTSTRN,FMSTRN,BSIGX,BEPSX,HB) IF (KEY(24).EQ.1) CALL MOMENT (NPLYS,SIGX,BSIGX,H,KEY) IF (KEY(24).EQ.2) CALL INTWCURV (NPLYS,CURV,QB,H,KEY,NA) CALL TSAIWU (F,F12,SIGX,BSIGX,NPLYS,NA,MATL,X,Y,XX,YY) RETURN END C C ***** SUBROUTINE LTDSTN ***** C SUBROUTINE LTDSTN (T0,T,HMC,NMATLS,MATL,NPLYS,ALPA,ALPAM,BETA,BETA 1M,FTSTRN,FMSTRN,X,Y,XX,YY,NA,MF) C ROUTINE TO CALCULATE FREE THERMAL & MOISTURE STRAINS FOR C LINEAR TEMPERATURE DEPENDENT PROPERTIES C C DOUBLE PRECISION X,Y,XX,YY IMPLICIT REAL*8(A-H,O-Z) DIMENSION ALPA(1),ALPAM(1),BETA(1),BETAM(1),FTSTRN(1),FMSTRN(1),M 1ATL(1),X(1),Y(1),XX(1),YY(1),NA(1) DT1 = (T + T0)/2. DT2 = (T - T0) DO 25 I = 1,NPLYS L = 3*I - 3 DO 25 J = 1,3 FTSTRN(L+J) = 0. 25 FMSTRN(L+J) = 0. DO 50 I =1, NPLYS KK = NA(I)*3 - 3 L = 3*I - 3 DO 50 J = 1,3 C MF = 0 FOR T INDEPENDENT PROPERTIES C MF = 1 FOR T DEPENDENT PROPERTIES C FTSTRN(L+J) = (ALPAM(KK+J)*DT1*MF + ALPA(KK+J)) * DT2 FMSTRN(L+J) = (BETAM(KK+J)*HMC*MF + BETA(KK+J)) * HMC 50 CONTINUE 150 WRITE(16,900) L = 0 DO 200 I = 1, NPLYS WRITE(16,901) I,FTSTRN(L+1),FTSTRN(L+2),FTSTRN(L+3) 200 L = L + 3 900 FORMAT (/T15,'*** FREE THERMAL STRAINS ***',/,8X,'PLY',10X,'EPSX', 115X,'EPSY',13X,'GAMMAXY') 901 FORMAT (5X,I5,3(5X,E14.6)) 902 FORMAT (/T15,'*** FREE MOISTURE STRAINS ***',/,8X,'PLY',10X,'EPSX' 1,15X,'EPSY',13X,'GAMMAXY') WRITE(16,902) L = 0 DO 300 I = 1, NPLYS WRITE(16,901) I,FMSTRN(L+1),FMSTRN(L+2),FMSTRN(L+3) 300 L = L + 3 100 RETURN END C C ***** SUBROUTINE PRPRTY ***** C SUBROUTINE PRPRTY(X2,X4,Y2,Y4,TTT,HO,AI,AN,AP,AW,BS,BP,BW,CS,CP,CW 1,DI,DSI,DS,DP,A,B,D,Q,QB,S,NA,KEY,H,MATL,EL,ET,ULT,G,U1,U2,U3,U4, 2X,Y,XX,YY,ALPA,ALPAX,ALPABA,BETA,BETAX,BETABA,XT,XC,YT,YC,SS, 3F12,ELM,ETM,GM,XTM,XCM,YTM,YCM,SSM,F12M,ALPAM,BETAM,ETT,ULTT,GT, 4XTT,XCT,YTT,YCT,SST,F12T,ULTM,ALPAT,BETAT,ELT,SB,NANGLS,NMATLS, 5ALPAMX,BETAMX,NPLYS,UTL) C C CALCULATES BASIC LAMINATE PROPERTIES GIVEN E'S, G'S, ETC C DOUBLE PRECISION X,Y,Y2,X2,Y4,X4,DSIN,DCOS,AN,G,A,AI,AD,ANG,AP,AW, C 1B,BS,BP,BW,C,CS,CP,CW,D,DI,DSI,DS,XX,YY,QA,QBB,ALPABA,BETABA IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(9),B(9),D(9),Q(20),QB(120),S(20),NA(20),KEY(25),H(21), 1MATL(20),EL(5),ET(5),ULT(5),G(5),U1(5),U2(5),U3(5),U4(5),UTL(5), 2X(10),Y(10),XX(10),YY(10),ALPA(10),ALPAX(60),ALPABA(3), 6BETA(10),BETAX(60),BETABA(3), 7XT(5),XC(5),YT(5),YC(5),SS(5),F12(5),ELM(5),ETM(5),GM(5),XTM(5),XC 8M(5),YTM(5),YCM(5),SSM(5),F12M(5),ALPAM(10),BETAM(10),ETT(5),ULTT( 95),GT(5),XTT(5),XCT(5),YTT(5),YCT(5),SST(5),F12T(5), AULTM(5),ALPAT(10),BETAT(10),ELT(5),U5(5),SB(120),ALPAMX(60), CBETAMX(60) DIMENSION X2(1),X4(1),Y2(1),Y4(1) DIMENSION AI(9),AN(10),AP(9),AW(9),BS(9),BP(9),BW(9),CS(9),CP(9) 1,CW(9),DI(9),DSI(9),DS(9),DP(9),QA(3),QBB(3) C C CALCULATE PROPERTIES AT GIVEN TEMPERATURE - LINEAR DEPENDENCE C WRITE (16,850) TTT 850 FORMAT(//T15,'PROPERTIES AT ',F5.1,'DEGREES F ABOVE ROOM TEMP.',/) DO 10 I = 1,9 A(I) = 0. B(I) = 0. 10 D(I) = 0. DO 11 I = 1,3 QBB(I) = 0.0 11 QA(I) = 0.0 DO 50 I = 1,NMATLS ELT(I) = ELM(I)*TTT + EL(I) ETT(I) = ETM(I) *TTT + ET(I) ULTT(I) = ULTM(I)*TTT + ULT(I) GT(I) = GM(I)*TTT + G(I) XTT(I) = XTM(I)*TTT + XT(I) XCT(I) = XCM(I)*TTT + XC(I) YTT(I) = YTM(I)*TTT + YT(I) YCT(I) = YCM(I)*TTT + YC(I) SST(I) = SSM(I)*TTT + SS(I) J = 2 * (I-1) BETAT(J+1) = BETAM(J+1) *TTT + BETA(J+1) ALPAT(J+1) = ALPAM(J+1) * TTT + ALPA(J+1) BETAT(J+2) = BETAM(J+2) *TTT + BETA(J+2) ALPAT(J+2) = ALPAM(J+2) *TTT + ALPA(J+2) 50 F12T(I) = F12M(I)*TTT+ F12(I) IF (KEY(21).EQ.0) GO TO 31 CALL RITPRP(NMATLS,ELT,ETT,ULTT,GT,ALPAT,BETAT,XTT,XCT,YTT,YCT,SST 1,F12T) 31 L = 0 C C STIFFNESS & COMPLIANCE MATRICES & INVARIANTS C DO 51 I = 1, NMATLS UTL(I) = ETT(I)*ULTT(I)/ELT(I) W = 1. - ULTT(I)*UTL(I) Q(L+1) = ELT(I)/W Q(L+2) = Q(L+1)*UTL(I) Q(L+3) = ETT(I)/W Q(L+4) = GT(I) Q12Q = Q(L+1) + Q(L+3) Q16Q = 2.*Q(L+2) + 4.*Q(L+4) U1(I) = (3.*Q12Q + Q16Q)/8. U2(I) = (Q(L+1)-Q(L+3))/2. U3(I) = (Q12Q - Q16Q)/8. U4(I) = (Q12Q +6.*Q(L+2)-4.*Q(L+4))/8. U5(I) = (Q12Q - 2.*Q(L+2) + 4.*Q(L+4))/8. S(L+1) = 1./ELT(I) S(L+2) = -ULTT(I)/ELT(I) S(L+3) = 1./ETT(I) S(L+4) = 1./GT(I) IF (KEY(21).EQ.0) GO TO 51 WRITE(16,915) I 915 FORMAT (//10X,'Q AND S MATRIX FOR MATERIAL NO.',I5// 110X,'Q MATRIX') K =I*4 - 4 916 FORMAT (2(10X,E14.6)/34X,E14.6,/48X,E14.6,//10X,'S MATRIX'/) WRITE (16,916) (Q(K+J), J = 1,4) WRITE(16,917) (S(K+J), J = 1,4) 917 FORMAT (2(10X,E14.6)/34X,E14.6,/48X,E14.6/) C C QBAR & SBAR FOR EACH PLY C 51 L = L + 4 M = 0 L = 0 DO 100 K = 1, NANGLS J = MATL(K) JJ = 4*MATL(K) - 4 QB(L+1) = U1(J) + U2(J)*Y2(K) + U3(J)*Y4(K) QB(L+2) = U4(J) - U3(J)*Y4(K) QB(L+3) = (U2(J)*X2(K))/2. + U3(J)*X4(K) PA = S(JJ+1) + S(JJ+3) PB = S(JJ+1) - S(JJ+3) PC = 2.*S(JJ+2) + S(JJ+4) QB(L+4) = U1(J) - U2(J)*Y2(K) + U3(J)*Y4(K) QB(L+5) = (U2(J)*X2(K))/2. - U3(J)*X4(K) QB(L+6) = U5(J) - U3(J)*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))*S(JJ+2)/8. + (PA-S(JJ+4))*(1-Y4(K))/8. SB(L+3) = PB/2. * X2(K) + X4(K)*(PA-PC)/4. 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. SB(L+6) = PA/2. + S(JJ+4)/2. - S(JJ+2) + (PC-PA)*Y4(K)/2. 100 L = L + 6 IF (KEY(1).EQ.0) GO TO 175 DO 150 I = 1, NANGLS L = I*6 - 6 WRITE(16,918) AN(I) 918 FORMAT (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,/) 150 CONTINUE C C CALCULATE THERMAL COEFFICIENTS FOR EACH PLY C 175 CALL TRNSFR (NANGLS,MATL,XX,YY,X,Y,ALPAT,ALPAX) CALL TRNSFR (NANGLS,MATL,XX,YY,X,Y,BETAT,BETAX) CALL TRNSFR (NANGLS,MATL,XX,YY,X,Y,BETAM,BETAMX) CALL TRNSFR (NANGLS,MATL,XX,YY,X,Y,ALPAM,ALPAMX) IF (KEY(21).EQ.0) GO TO 255 WRITE (16,922) 922 FORMAT (///25X,'*** THERMAL COEFFICIENTS ***' 1 ,//10X,'ANGLE',10X,'ALPHAX',13X,'ALPHAY',10X,'ALPHAXY'/) L = 0 DO 250 I = 1, NANGLS WRITE (16,923) AN(I),ALPAX(L+1),ALPAX(L+2),ALPAX(L+3) 923 FORMAT (7X,F7.2,3(5X,E14.6)) 250 L = L + 3 WRITE (16,1922) 1922 FORMAT (///25X,'*** MOISTURE COEFFICIENTS ***' 1 ,//10X,'ANGLE',10X,'BETAX',14X,'BETAY',11X,'BETAXY'/) L = 0 DO 251 I = 1, NANGLS WRITE (16,923) AN(I),BETAX(L+1),BETAX(L+2),BETAX(L+3) 251 L = L + 3 C C CALCULATE [A],[B],[C],ALPABAR & BETABAR FOR LAMINATE C 255 DO 300 I = 1, NPLYS DH = H(I+1) - H(I) DH2 =(H(I+1)**2 - H(I)**2)/2. DH3 =(H(I+1)**3 - H(I)**3)/3. J = 6 * NA(I) - 6 A(1) = A(1) + QB(J+1)*DH B(1) = B(1) + QB(J+1)*DH2 D(1) = D(1) + QB(J+1)*DH3 A(2) = A(2) + QB(J+2) * DH B(2) = B(2) + QB(J+2) * DH2 D(2) = D(2) + QB(J+2) * DH3 A(3) = A(3) + QB(J+3)*DH B(3) = B(3) + QB(J+3) * DH2 D(3) = D(3) + QB(J+3) * DH3 A(5) = A(5) + QB(J+4) * DH B(5) = B(5) + QB(J+4) * DH2 D(5) = D(5) + QB(J+4) * DH3 A(6) = A(6) + QB(J+5) * DH B(6) = B(6) + QB(J+5) * DH2 D(6) = D (6) + QB(J+5) * DH3 A(9) = A(9) + QB(J+6) * DH B(9) = B(9) + QB(J+6) * DH2 D(9) = D(9) + QB(J+6) * DH3 K = 3 * (NA(I) - 1) QA(1)=QA(1)+(QB(J+1)*ALPAX(1+K)+QB(J+2)*ALPAX(2+K)+QB(J+3)* 1ALPAX(3+K))*DH QA(2)=QA(2)+(QB(J+2)*ALPAX(1+K)+QB(J+4)*ALPAX(2+K)+QB(J+5)* 1ALPAX(3+K))*DH QA(3)=QA(3)+(QB(J+3)*ALPAX(1+K)+QB(J+5)*ALPAX(2+K)+QB(J+6)* 1ALPAX(3+K))*DH QBB(1)=QBB(1)+(QB(J+1)*BETAX(1+K)+QB(J+2)*BETAX(2+K)+QB(J+3)* 1BETAX(3+K))*DH QBB(2)=QBB(2)+(QB(J+2)*BETAX(1+K)+QB(J+4)*BETAX(2+K)+QB(J+5)* 1BETAX(3+K))*DH QBB(3)=QBB(3)+(QB(J+3)*BETAX(1+K)+QB(J+5)*BETAX(2+K)+QB(J+6)* 1BETAX(3+K))*DH 300 CONTINUE WRITE (16,930) QA(1),QA(2),QA(3) 930 FORMAT (//,10X,'EQUIVALENT UNIT THERMAL FORCES',/,10X,3(E14.6,3X)) A(4) = A(2) B(4) = B(2) D(4) = D(2) A(7) = A(3) B(7) = B(3) D(7) = D(3) A(8) = A(6) B(8) = B(6) D(8) = D(6) WRITE(16,924) 924 FORMAT (////10X,'A MATRIX'/) 925 FORMAT (10X,3(5X,E14.6)/29X,2(5X,E14.6)/53X,E14.6) WRITE(16,925) A(1),A(4),A(7),A(5),A(8),A(9) WRITE(16,926) 926 FORMAT (//10X,'B MATRIX'/) WRITE(16,925) B(1),B(4),B(7),B(5),B(8),B(9) WRITE (16,927) 927 FORMAT (//10X,'D MATRIX'/) WRITE(16,925) D(1),D(4),D(7),D(5),D(8),D(9) 305 CALL INVERT (A,AI) C C CALCULATE LAMINATE ENGINEERING CONSTANTS C TT = 1./(2.*HO) EXBAR = 1./AI(1)*TT EYBAR = 1./AI(5)*TT UXYBAR = -AI(4)/AI(1) GXYBAR = 1./AI(9)*TT UYXBAR = -AI(2)/AI(5) ETXYCX = AI(3)/AI(1) ETXYCY = AI(6)/AI(5) ETXCXY = AI(3)/AI(9) ETYCXY = AI(6)/AI(9) 2000 CALL ABD(A,B,D,AI,AP,AW,BS,BP,BW,CS,CP,CW,DS,DSI,DI) WRITE (16,931) 931 FORMAT (//10X,'INVERTED A MATRIX'/) WRITE (16,925) AI(1),AI(4),AI(7),AI(5),AI(8),AI(9) WRITE (16, 940) 940 FORMAT (///10X,' A PRIME MATRIX'/) WRITE (16,941) (AP(I), I = 1,9) WRITE (16,942) 942 FORMAT (///10X,' B PRIME MATRIX'/) WRITE (16, 941) (BP(I), I = 1,9) 941 FORMAT ((10X, 3(5X,E14.6))/) WRITE (16, 943) 943 FORMAT (///10X,' C PRIME MATRIX'/) WRITE (16,941) (CP(I), I = 1,9) WRITE (16,944) 944 FORMAT (///' D PRIME MATRIX'/) WRITE (16, 941) (DSI(I), I = 1,9) WRITE(16,898) 898 FORMAT ('0',T20,'***LAMINATE ENGINEERING CONSTANTS***') WRITE (16,928) EXBAR,EYBAR,UXYBAR,UYXBAR,GXYBAR,ETXYCX,ETXYCY, 1ETXCXY,ETYCXY 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,/) C CALCULATE LAMINATE ALPHABAR AND BETABAR CALL MULT3N (AI,QA,ALPABA) CALL MULT3N (AI,QBB,BETABA) WRITE (16,929) ALPABA(1),BETABA(1),ALPABA(2),BETABA(2), 1ALPABA(3),BETABA(3) 929 FORMAT (/10X,'ALPAX =',E14.6,7X,'BETAX =',E14.6,/,10X,'ALPAY =', 1E14.6,7X,'BETAY =',E14.6,/,10X,'ALAPXY =',E14.6,7X,'BETAXY =', 2E14.6,/) 311 RETURN END C C ********* SUBROUTINE MULT3N *********** C MULTIPLIES A (3X3) MATRIX [A] INTO A (3X1) COLUMN VECTOR B C THE PRODUCT IS THE (3X1) COLUMN VECTOR C C SUBROUTINE MULT3N (A,B,C) C DOUBLE PRECISION A,B,C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(9),B(3),C(3) C(1) = A(1)*B(1)+A(2)*B(2)+A(3)*B(3) C(2) = A(4)*B(1)+A(5)*B(2)+A(6)*B(3) C(3) = A(7)*B(1)+A(8)*B(2)+A(9)*B(3) RETURN END C C ***** SUBROUTINE RITPRP ***** C SUBROUTINE RITPRP (NMATLS,EL,ET,ULT,G,ALPA,BETA,XT,XC,YT,YC,SS,F12 1) IMPLICIT REAL*8(A-H,O-Z) DIMENSION EL(1),ET(1),ULT(1),G(1),ALPA(1),BETA(1),XT(1),YT(1),SS(1 1),F12(1),XC(1),YC(1) WRITE(16,910) (I,EL(I),ET(I),ULT(I),G(I),I =1,NMATLS) 910 FORMAT ('0',T30,'***MATERIAL PROPERTIES***',//, 110X,'MATL NO ',I3,/10X,'E11 = ',E14.6,5X,'ET(1)2 = ',E14.6,/, 210X,'U12 = ',E14.6,5X,'G12 = ',E14.6,/) L = 0 DO 40 J = 1, NMATLS WRITE(16,912) J, ALPA(L+1), ALPA(L+2),BETA(L+1),BETA(L+2) 40 L = L + 2 912 FORMAT(/10X,'MATL NO. ',I3,/10X,'ALPHA1 =',E14.6,4X,'ALPHA2 =', 1E14.6,/,10X,'BETA1 =',E14.6,4X,'BETA2 =',E14.6,/) DO 45 I =1, NMATLS WRITE (16,913) I 913 FORMAT (//10X,'STRENGTH PARAMETERS FOR MATERIAL NO.',I3,/) WRITE(16,914) XT(I),XC(I),YT(I),YC(I),SS(I),F12(I) 914 FORMAT (10X,'XT = ',E14.6,3X,'XC = ',E14.6,/ 110X,'YT = ',E14.6,3X,'YC = ',E14.6,/ 210X,' S = ',E14.6,3X,'F12 = ',E14.6,//) 45 CONTINUE RETURN END