PROGRAM lelpa_a REAL NU12,NU21,NULL,NUXY,NXYCX,NXYCY,NR,NB,MR,MB,LT,MEXP12,M *EXPXY,MXP12,MXPXY,MEXPB,NW,MW DIMENSION THETA(16),T(16),MTL(16),Z(16),ZAVG(16) DIMENSION E11(16),E22(16),G12(16),NU12(16),NU21(16) DIMENSION TEXP12(3,16),MEXP12(3,16) DIMENSION QK(3,3,16),SK(3,3,16),QB(3,3,16),SB(3,3,16) DIMENSION A(3,3),B(3,3),D(3,3),SM(6,6),V(6) DIMENSION A2(3,3),B2(3,3),C2(3,3),D2(3,3) DIMENSION AS(3,3),BS(3,3),CS(3,3),DS(3,3),RM(6,6),R(6) DIMENSION AR(3,3),BR(3,3),CR(3,3),DR(3,3),BM(6,6),W(6) DIMENSION AP(3,3),BP(3,3),CP(3,3),DP(3,3),CM(6,6),U(6) DIMENSION TXP12(3),MXP12(3),RT(3,3),TXPXY(3),MXPXY(3) DIMENSION TEXPXY(3,16),MEXPXY(3,16),TQ(3,3),AB(3),PP(3) DIMENSION BETA(3),GAMMA(3),PHI(3),PSI(3),TEXPB(3),MEXPB(3) DIMENSION NR(3),MR(3),NW(3),MW(3),NB(3),MB(3),LT(6),BLANK(6),VT(6) DIMENSION EP(3),SIG(3),NULL(3,3),TR(3,3),EP12(3),SIG12(3) DIMENSION EPF(3,16),EPF12(3,16),SIGF(3,16),SIGF12(3,16) DIMENSION SIGX(65),SIGY(65),SIGXY(65),SIGZ(65),X(65) CHARACTER ANGL(5,16)*1 CHARACTER HEAD(80)*1 WRITE(6,999) C OPEN(5,FILE='lelpa_a.dat',STATUS='OLD',ERR=888) OPEN(6,FILE='lelpa_a.out',STATUS='UNKNOWN',ERR=888) C C CONSTANTS PI=ACOS(-1.0) C C ZERO ARRAYS C DO 110 I=1,3 DO 110 J=1,3 A(I,J)=0.0 B(I,J)=0.0 D(I,J)=0.0 NULL(I,J)=0.0 GAMMA(I)=0.0 BETA(I)=0.0 PHI(I)=0.0 PSI(I)=0.0 DO 100 L=1,16 TEXP12(I,L)=0.0 MEXP12(I,L)=0.0 QK(I,J,L)=0.0 SK(I,J,L)=0.0 100 CONTINUE 110 CONTINUE C C READ HEAD ON FIRST DATA CARD C READ(5,10)(HEAD(I),I=1,80) 10 FORMAT(80A1) WRITE(6,12)HEAD 12 FORMAT(1X,80A1,/) C C READ AND WRITE OUT ORIENTATION AND THICKNESS C FOR EACH LAYER C READ(5,14)NL,NM 14 FORMAT(14X,I2,26X,I2) DO 120 ILAYER=1,NL READ(5,16)THETA(ILAYER),T(ILAYER),MTL(ILAYER),(ANGL(I,ILAYER),I=1, *5) 120 CONTINUE 16 FORMAT(18X,F5.1,11X,F9.6,14X,I2,12X,5A1) WRITE(6,18) DO 130 ILAYER =1,NL WRITE(6,20)ILAYER,THETA(ILAYER),T(ILAYER),MTL(ILAYER) 130 CONTINUE 18 FORMAT(31X,'LAMINAE STACKING',/,17X,'LAYER NO. ORIENTATION THICKNE *SS MTL. TYPE') 20 FORMAT(20X,I2,9X,F5.1,4X,F8.5,6X,I2) C C CALCULATE AND WRITE OUT TOTAL THICKNESS AND C DISTANCES FROM THE MIDPLANE TO EACH LAMINAE C TH=0.0 DO 140 K=1,NL TH=TH+T(K) 140 CONTINUE WRITE(6,22)TH 22 FORMAT(/,' TOTAL THICKNESS=',E10.3) Z(1)=-TH/2. NLP=NL+1 DO 150 K=2,NLP Z(K)=Z(K-1)+T(K-1) ZAVG(K-1)=(Z(K-1)+Z(K))/2. 150 CONTINUE C C Q*S ARE CALCULATED FROM MATERIAL PROPERTIES C C READ AND WRITE OUT MATERIAL PROPERTIES FOR EACH MATERIAL TYPE C DO 170 MAT=1,NM READ(5,24)K,E11(K),E22(K),G12(K),NU12(K),TEXP12(1,K),TEXP12(2,K),M *EXP12(1,K),MEXP12(2,K),CV WRITE(6,26)K,E11(K),E22(K),G12(K),NU12(K),TEXP12(1,K),TEXP12(2,K), *MEXP12(1,K),MEXP12(2,K),CV 24 FORMAT(9X,I2,4(5X,E10.3),/,6X,E10.3,3(7X,E10.3),/,30X,F4.1) 26 FORMAT(//,' MATERIAL TYPE',I2,/,26X,'ELASTIC CONSTANTS',/,6X,'EL = *',E10.3,' ET =',E10.3,' GLT=',E10.3,' NULT=',E10.3,/,11X,'THERMAL *COEFFICIENTS OF EXPANSION L-T COORD.',/,17X,'TEXPL=',E10.3,' TEXPT *=',E10.3,/,11X,'MOISTURE COEFFICIENTS OF EXPANSION L-T COORD.',/,1 *7X,'MEXPL=',E10.3,' MEXPT=',E10.3,/,11X,'MOISTURE VOLUME THRESHOLD *, CV=',E10.3,/) C C CALCULATE AND WRITE OUT STIFFNESSES AND COMPLIANCES FOR C EACH MATERIAL TYPE IN THE PRINCIPAL 1-2 DIRECTION C NU21(K)=NU12(K)*E22(K)/E11(K) QK(1,1,K)=E11(K)/(1.-NU12(K)*NU21(K)) QK(2,2,K)=E22(K)/(1.-NU12(K)*NU21(K)) QK(1,2,K)=NU21(K)*QK(1,1,K) QK(2,1,K)=QK(1,2,K) QK(3,3,K)=G12(K) SK(1,1,K)=1./E11(K) IF(E22(K).GT.0.0)SK(2,2,K)=1./E22(K) IF(E22(K).EQ.0.0)SK(2,2,K)=0.0 SK(1,2,K)=-NU12(K)/E11(K) SK(2,1,K)=SK(1,2,K) IF(G12(K).GT.0.0)SK(3,3,K)=1./G12(K) IF(G12(K).EQ.0.0)SK(3,3,K)=0.0 WRITE(6,28) DO 160 I=1,3 WRITE(6,30)(QK(I,J,K),J=1,3),(SK(I,J,K),J=1,3) 160 CONTINUE 28 FORMAT(10X,'Q STIFFNESS MATRIX',21X,'S COMPLIANCE MATRIX') 30 FORMAT(6(1X,E12.5)) 170 CONTINUE C C CALCULATE AND WRITE OUT STIFFNESSES AND COMPLIANCES FOR EACH C LAMINA IN THE X-Y DIRECTION C WRITE(6,999) DO 180 L=1,NL K=MTL(L) U1=(3.*QK(1,1,K)+3.*QK(2,2,K)+2.*QK(1,2,K)+4.*QK(3,3,K))/8. U2=(QK(1,1,K)-QK(2,2,K))/2. U3=(QK(1,1,K)+QK(2,2,K)-2.*QK(1,2,K)-4.*QK(3,3,K))/8. U4=(QK(1,1,K)+QK(2,2,K)+6.*QK(1,2,K)-4.*QK(3,3,K))/8. U5=(QK(1,1,K)+QK(2,2,K)-2.*QK(1,2,K)+4.*QK(3,3,K))/8. ANGLE=THETA(L)*PI/180. THETA4=4.*ANGLE THETA2=2.*ANGLE C1=COS(ANGLE) SN=SIN(ANGLE) CS2=C1**2 SN2=SN**2 CS3=C1**3 SN3=SN**3 CS4=C1**4 SN4=SN**4 QB(1,1,L)=U1+U2*COS(THETA2)+U3*COS(THETA4) QB(2,2,L)=U1-U2*COS(THETA2)+U3*COS(THETA4) QB(1,2,L)=U4-U3*COS(THETA4) QB(2,1,L)=QB(1,2,L) QB(3,3,L)=U5-U3*COS(THETA4) QB(1,3,L)=U2*SIN(THETA2)/2.+U3*SIN(THETA4) QB(3,1,L)=QB(1,3,L) QB(2,3,L)=U2*SIN(THETA2)/2.-U3*SIN(THETA4) QB(3,2,L)=QB(2,3,L) SB(1,1,L)=SK(1,1,K)*CS4+(2.*SK(1,2,K)+SK(3,3,K))*SN2*CS2+SK(2,2,K) **SN4 SB(1,2,L)=SK(1,2,K)*(SN4+CS4)+(SK(1,1,K)+SK(2,2,K)-SK(3,3,K))*SN2* *CS2 SB(2,2,L)=SK(1,1,K)*SN4+(2.*SK(1,2,K)+SK(3,3,K))*SN2*CS2+SK(2,2,K) **CS4 SB(1,3,L)=(2.*SK(1,1,K)-2.*SK(1,2,K)-SK(3,3,K))*SN*CS3-(2.*SK(2,2, *K)-2.*SK(1,2,K)-SK(3,3,K))*SN3*C1 SB(2,3,L)=(2.*SK(1,1,K)-2.*SK(1,2,K)-SK(3,3,K))*SN3*C1-(2.*SK(2,2, *K)-2.*SK(1,2,K)-SK(3,3,K))*SN*CS3 SB(3,3,L)=(4.*SK(1,1,K)+4.*SK(2,2,K)-8.*SK(1,2,K)-2.*SK(3,3,K))*SN *2*CS2+SK(3,3,K)*(SN4+CS4) SB(2,1,L)=SB(1,2,L) SB(3,1,L)=SB(1,3,L) SB(3,2,L)=SB(2,3,L) 180 CONTINUE C C CALCULATE AND WRITE OUT LAMINATE C A B D STIFFNESSES C DO 190 I=1,3 DO 190 J=1,3 DO 190 K=1,NL A(I,J)=A(I,J)+QB(I,J,K)*(Z(K+1)-Z(K)) B(I,J)=B(I,J)+(QB(I,J,K)*(Z(K+1)**2-Z(K)**2))/2. D(I,J)=D(I,J)+(QB(I,J,K)*(Z(K+1)**3-Z(K)**3))/3. 190 CONTINUE WRITE(6,12)HEAD WRITE(6,32) 32 FORMAT(25X,'LAMINATE STIFFNESS MATRIX',/,35X,'A I B',/,34X,'---+-- *-',/,35X,'B I D',/) DO 200 I=1,3 WRITE(6,36)A(I,1),A(I,2),A(I,3),B(I,1),B(I,2),B(I,3) 200 CONTINUE WRITE(6,34) 34 FORMAT(3X,34('-'),'+',34('-')) DO 210 I=1,3 WRITE(6,36)B(I,1),B(I,2),B(I,3),D(I,1),D(I,2),D(I,3) 210 CONTINUE 36 FORMAT(3X,1X,1PE10.3,2(1X,E10.3),' I ',3(1X,E10.3)) C C INVERT THE ABD MATRIX (LAMINATE COMPLIANCES) C CALL MINV(A,AS) CALL MULT(B2,AS,B) CALL SUBT(BS,NULL,B2) CALL MULT(CS,B,AS) CALL MULT(D2,B,B2) CALL SUBT(DS,D,D2) CALL MINV(DS,DP) CALL MULT(C2,DP,CS) CALL SUBT(CP,NULL,C2) CALL MULT(BP,BS,DP) CALL MULT(A2,BS,C2) CALL SUBT(AP,AS,A2) CALL MINV(D,DR) CALL MULT(B2,DR,B) CALL SUBT(CR,NULL,B2) CALL MULT(BR,B,DR) CALL MULT(D2,B,B2) CALL SUBT(AR,A,D2) C C CALCULATE LAMINATE RESULTANT STIFFNESSES, POISSON RATIO C AND COEFFICIENTS OF MUTUAL INFLUENCE C EXX=1./(AP(1,1)*TH) EYY=1./(AP(2,2)*TH) GXY=1./(AP(3,3)*TH) NUXY=-AP(1,2)/AP(1,1) NXYCX=AP(1,3)/AP(1,1) NXYCY=AP(2,3)/AP(2,2) WRITE(6,38)EXX,EYY,GXY,NUXY,NXYCX,NXYCY 38 FORMAT(//,23X,'LAMINATE RESULTANT STIFFNESSES',/,8X,'EXX=',1PE10.3 *,' EYY=',E10.3,' GXY=',E10.3,' NUXY=',E10.3,/,16X,'LAMINATE COEFFI *CIENTS OF MUTUAL INFLUENCE',/,20X,'NXY,X=',E10.3,' NXY,Y=',E10.3,/ *) C C CONVERT ABD & A*B*C*D* MATRICES C INTO SINGLE MATRICES SM & CM C DO 220 I=1,3 DO 220 J=1,3 SM(I,J)=A(I,J) SM(I,J+3)=B(I,J) SM(I+3,J)=B(I,J) SM(I+3,J+3)=D(I,J) CM(I,J)=AP(I,J) CM(I,J+3)=BP(I,J) CM(I+3,J)=CP(I,J) CM(I+3,J+3)=DP(I,J) 220 CONTINUE C C CONVERT AS BS CS DS MATRICES AND AR BR CR DR MATRICES C INTO SINGLE MATRICES RM AND BM DO 230 I=1,3 DO 230 J=1,3 RM(I,J)=AS(I,J) RM(I,J+3)=BS(I,J) RM(I+3,J)=CS(I,J) RM(I+3,J+3)=DS(I,J) BM(I,J)=AR(I,J) BM(I,J+3)=BR(I,J) BM(I+3,J)=CR(I,J) BM(I+3,J+3)=DR(I,J) 230 CONTINUE C C TRANSFORM 1-2 COEFF OF EXPANSION INTO X-Y COORD. FOR EACH LAYER C DO 260 L=1,NL K=MTL(L) DO 240 I=1,3 TXP12(I)=TEXP12(I,K) MXP12(I)=MEXP12(I,K) 240 CONTINUE C C CALCULATE TRANSFORMATION MATRIX C FROM 1-2 TO X-Y COORD. SYSTEM C ANGLE=THETA(L)*PI/180. RT(1,1)=COS(ANGLE)**2 RT(2,2)=RT(1,1) RT(1,2)=SIN(ANGLE)**2 RT(3,3)=RT(1,1)-RT(1,2) RT(1,3)=-2.*SIN(ANGLE)*COS(ANGLE) RT(2,3)=-RT(1,3) RT(2,1)=RT(1,2) RT(3,1)=-RT(1,3)/2. RT(3,2)=-RT(2,3)/2. C CALL TXY12(TXP12,TXPXY,RT) CALL TXY12(MXP12,MXPXY,RT) TXPXY(3)=TXPXY(3)*2. MXPXY(3)=MXPXY(3)*2. DO 250 I=1,3 TEXPXY(I,L)=TXPXY(I) MEXPXY(I,L)=MXPXY(I) 250 CONTINUE 260 CONTINUE C C CALCULATE C BETA & GAMMA TEMP COEFF OF EXPANSION C PHI & PSI MOISTURE COEFF OF EXPANSION C FOR THE LAMINATE IN THE X-Y COORD. C DO 290 K=1,NL DO 270 I=1,3 TXPXY(I)=TEXPXY(I,K) MXPXY(I)=MEXPXY(I,K) DO 270 J=1,3 TQ(I,J)=QB(I,J,K) 270 CONTINUE CALL SS(1,TXPXY,AB,TQ,NULL) CALL SS(1,MXPXY,PP,TQ,NULL) DO 280 I=1,3 BETA(I)=BETA(I)+AB(I)*(Z(K+1)-Z(K)) GAMMA(I)=GAMMA(I)+AB(I)*(Z(K+1)**2-Z(K)**2)/2. PHI(I)=PHI(I)+PP(I)*(Z(K+1)-Z(K)) PSI(I)=PSI(I)+PP(I)*(Z(K+1)**2-Z(K)**2)/2. 280 CONTINUE 290 CONTINUE WRITE(6,40)(BETA(I),I=1,3),(GAMMA(I),I=1,3),(PHI(I),I=1,3),(PSI(I) *,I=1,3) 40 FORMAT(18X,'THERMAL COEFFICIENTS FOR MIDPLANE LOADS',/,12X,' BETA- *1=',1PE10.3,' BETA-2=',E10.3,' BETA-3=',E10.3,/,12X,'GAMMA-1=',E *10.3,' GAMMA-2=',E10.3,' GAMMA-3=',E10.3,/,17X,'MOISTURE COEFFICIE *NTS FOR MIDPLANE LOADS',/,14X,'PHI-1=',E10.3,' PHI-2=',E10.3,' PHI *-3=',E10.3,/,14X,'PSI-1=',E10.3,' PSI-2=',E10.3,' PSI-3=',E10.3,// *) DO 310 I=1,3 VALUET=0.0 VALUEM=0.0 DO 300 J=1,3 VALUET=VALUET+AS(I,J)*BETA(J) VALUEM=VALUEM+AS(I,J)*PHI(J) 300 CONTINUE TEXPB(I)=VALUET MEXPB(I)=VALUEM 310 CONTINUE WRITE(6,42)TEXPB(1),MEXPB(1),TEXPB(2),MEXPB(2),TEXPB(3),MEXPB(3) 42 FORMAT(10X,'COEFFICIENTS OF EXPANSION FOR LAMINATE X-Y COORD.',/,1 *9X,'THERMAL COEFF. MOISTURE COEFF.',/,5X,'X-DIRECTION ',1PE1 *0.3,8X,E10.3,/,5X,'Y-DIRECTION ',E10.3,8X,E10.3,/,3X,'SHEAR XY- *PLANE ',E10.3,8X,E10.3) WRITE(6,999) C C READ IN NUMBER OF LOAD CASES C READ(5,44)NCASE 44 FORMAT(21X,I2) WRITE(6,12)HEAD C C DO LOOP THRU NUMBER OF LOAD CASES C DO 430 ICASE=1,NCASE READ(5,46)LDCASE 46 FORMAT(40X,I2) WRITE(6,48)LDCASE 48 FORMAT(3X,'LOAD CASE NUMBER=',I2) C C MIDPLANE RESPONSE WITHOUT THE INFLUENCE OF C TEMPERATURE OR MOISTURE C WRITE(6,50) 50 FORMAT(/,6X,'LAMINATE RESPONSE WITHOUT INFLUENCE OF TEMPERATURE OR * MOISTURE',//) C C READ IN AND WRITE OUT C V MIDPLANE STRAINS & CURVATURES MSW=-1 C R MIDPLANE STRESSES & CURVATURES MSW=0 C U MIDPLANE STRESSES & MOMENTS MSW=+1 C W MIDPLANE STRAINS & MOMENTS MSW=+2 C READ(5,54)MSW IF(MSW.EQ.0)READ(5,52)(R(I),I=1,6) IF(MSW.EQ.1)READ(5,52)(U(I),I=1,6) IF(MSW.EQ.-1)READ(5,52)(V(I),I=1,6) IF(MSW.EQ.2)READ(5,52)(W(I),I=1,6) 52 FORMAT(2(3X,E10.3),4X,E10.3,2(3X,E10.3),4X,E10.3) 54 FORMAT(I2) IF(MSW.EQ.0)WRITE(6,60)(R(I),I=1,6) IF(MSW.EQ.1)WRITE(6,56)(U(I),I=1,6) IF(MSW.EQ.-1)WRITE(6,58)(V(I),I=1,6) IF(MSW.EQ.2)WRITE(6,62)(W(I),I=1,6) 56 FORMAT(23X,'GIVEN LAMINATE MIDPLANE LOADS',/,16X,'NX=',1PE10.3,' N *Y=',E10.3,' NXY=',E10.3,/,16X,'MX=',E10.3,' MY=',E10.3,' MXY=',E10 *.3,//) 58 FORMAT(15X,'GIVEN LAMINATE MIDPLANE STRAINS AND CURVATURES',/,16X, *'EXO=',1PE10.3,' EYO=',E10.3,' EXYO=',E10.3,/,17X,'KX=',E10.3,' K *Y=',E10.3,' KXY=',E10.3,//) 60 FORMAT(14X,'GIVEN LAMINATE MIDPLANE STRESSES AND CURVATURES',/,16X *,'NX=',1PE10.3,' NY=',E10.3,' NXY=',E10.3,/,16X,'KX=',E10.3,' KY=' *,E10.3,' NXY=',E10.3,//) 62 FORMAT(16X,'GIVEN LAMINATE MIDPLANE STRAINS AND MOMENTS',/,16X,'EX *O=',1PE10.3,' EYO=',E10.3,' EXYO=',E10.3,/,17X,'MX=',E10.3,' MY=' *,E10.3,' MXY=',E10.3,//) C C AT THE LAMINATE MIDPLANE C CALCULATES U LOADS MSW=-1 C CALCULATES W STRAINS AND MOMENTS MSW=0 C CALCULATES V STRAINS AND CURVATURES MSW=+1 C CALCULATES R STRESSES AND CURVATURES MSW=+2 C CALL LMSS(MSW,R,W,U,V,RM,SM,CM,BM) C C WRITE OUT RESULTS C IF(MSW.EQ.1)WRITE(6,64)(V(I),I=1,6) 64 FORMAT(18X,'CALCULATED MIDPLANE STRAINS AND CURVATURE',/,16X,'EXO= *',1PE10.3,' EYO=',E10.3,' EXYO=',E10.3,/,17X,'KX=',E10.3,' KY=',E *10.3,' KXY=',E10.3,//) IF(MSW.EQ.-1)WRITE(6,66)(U(I),I=1,6) 66 FORMAT(21X,'CALCULATED MIDPLANE LOADS',/,16X,'NX=',1PE10.3,' NY=', *E10.3,' NXY=',E10.3,/,16X,'MX=',E10.3,' MY=',E10.3,' MXY=',E10.3,/ */) IF(MSW.EQ.0)WRITE(6,68)(W(I),I=1,6) 68 FORMAT(19X,'CALCULATED MIDPLANE STRAINS AND MOMENTS',/,16X,'EXO=', *1PE10.3,' EYO=',E10.3,' EXYO=',E10.3,/,17X,'MX=',E10.3,' MY=',E10 *.3,' MXY=',E10.3,//) IF(MSW.EQ.2)WRITE(6,70)(R(I),I=1,3) 70 FORMAT(17X,'CALCULATED MIDPLANE STRESSES AND CURVAUTRES',/,16X,'NX *=',1PE10.3,' NY=',E10.3,' NXY=',E10.3,/,16X,'KX=',E10.3,' KY=',E10 *.3,' KXY=',E10.3,//) C C GENERATE U AND V ARRAYS C IF ARRAYS W AND R WHERE USED C IF(MSW.EQ.-1.OR.MSW.EQ.1)GO TO 321 DO 320 I=1,3 U(I)=R(I) U(I+3)=W(I+3) V(I)=W(I) V(I+3)=R(I+3) 320 CONTINUE 321 CONTINUE C C READ IN C STRESS FREE TEMPERATURE AND OPERATING TEMPERATURE C AND C CHANGE IN PERCENT MOISTURE ABSORBED FROM C THE DRY STATE C C CALCULATE TEMPERATURE DIFFERENCE BETWEEN STRESS FREE TEMP C AND C OPERATING TEMPERATURE C READ(5,72)SFT,OPT,DW 72 FORMAT(17X,F6.1,16X,F6.1,28X,F6.3) DELT=OPT-SFT WRITE(6,74)SFT,OPT,DW 74 FORMAT(' STRESS FREE TEMP.=',F6.1,/,' OPERATING TEMP.=',F6.1,/,' P *ERCENT MOSITURE ABSORBED=',F6.3,//) IF(DW.LT.CV)WRITE(6,911) 911 FORMAT(12X,' PERCENT MOSITURE ABSORBED IS LESS THEN MOISTURE THRES *HOLD',/,12X,' THEREFORE PERCENT MOISTURE ABSORBED HAS BEEN DEFINED * AS ZERO',//) C C CALCULATE EFFECTIVE MIDPLANE LOADS DUE TO TEMP C CHANGE FROM STRESS FREE TEMP TO OPERATING TEMP (DELT3) C AND DUE TO PERCENT MOISTURE WEIGHT GAIN DW C DO 330 I=1,3 NR(I)=DELT*BETA(I) MR(I)=DELT*GAMMA(I) DWMCV=DW-CV IF(DW.LT.CV)DWMCV=0.0 NW(I)=(DWMCV)*PHI(I) MW(I)=(DWMCV)*PSI(I) 330 CONTINUE WRITE(6,76)DELT,(NR(I),I=1,3),(MR(I),I=1,3),DW,(NW(I),I=1,3),(MW(I *),I=1,3) 76 FORMAT(6X,'EFFECTIVE MIDPLANE LOADS DUE TO TEMP CHANGE FROM STRESS * FREE TEMP',/,21X,'TO OPERATING TEMP DELT=',F6.1,/,14X,'NX=',1PE1 *0.3,' NY=',E10.3,' NXY=',E10.3,/,14X,'MX=',E10.3,' MY=',E10.3,' MX *Y=',E10.3,/,9X,'EFFECTIVE MIDPLANE LOADS DUE TO A CHANGE IN MOISTU *RE DW=',F6.3,/,14X,'NX=',E10.3,' NY=',E10.3,' NXY=',E10.3,/,14X,'M *X=',E10.3,' MY=',E10.3,' MXY=',E10.3,//) C C CALCULATE MIDPLANE LOADS DUE TO CURING, MOISTURE, C AND APPLIED LOAD DO 340 I=1,3 NB(I)=NR(I)+NW(I)+U(I) MB(I)=MR(I)+MW(I)+U(I+3) 340 CONTINUE WRITE(6,78)(NB(I),I=1,3),(MB(I),I=1,3) 78 FORMAT(7X,'LAMINATE RESPONSE WITH THE INFLUENCE OF TEMPERATURE AND * MOISTURE',//,11X,' MIDPLANE LOADS DUE TO CURING, MOISTURE, AND AP *PLIED LOAD',/,14X,'NXBAR=',1PE10.3,' NYBAR=',E10.3,' NXYBAR=',E10. *3,/,14X,'MXBAR=',E10.3,' MYBAR=',E10.3,' MXYBAR=',E10.3,//) C C CALCULATE MIDPLANE STRAINS & CURVATURES FROM MIDPLANE C LOADS DUE TO CURING, MOISTURE, AND APPLIED LOAD C BUT FIRST COMBINE NB & MB INTO ONE MATRIX LT C DO 350 I=1,3 LT(I)=NB(I) LT(I+3)=MB(I) 350 CONTINUE CALL LMSS(1,BLANK,BLANK,LT,VT,RM,SM,CM,BM) WRITE(6,80)(VT(I),I=1,6) 80 FORMAT(4X,'MIDPLANE STRAINS & CURVATURES DUE TO CURING, MOISTURE, *AND APPLIED LOAD',/,17X,'EXO=',1PE10.3,' EYO=',E10.3,' EXYO=',E10. *3,/,18X,'KX=',E10.3,' KY=',E10.3,' KXY=',E10.3) WRITE(6,999) C C CALCULATE STRAINS & STRESSES IN EACH LAYER DUE TO C CURING, MOISTURE, AND APPLIED LOAD DO 380 K=1,NL DO 360 I=1,3 DWMCV=DW-CV IF(DW.LT.CV)DWMCV=0.0 EP(I)=VT(I)+VT(I+3)*ZAVG(K)-DELT*TEXPXY(I,K)-(DWMCV)*MEXPXY(I,K) DO 360 J=1,3 TQ(I,J)=QB(I,J,K) 360 CONTINUE CALL SS(1,EP,SIG,TQ,NULL) C C CALCULATE TRANSFORMATION MATRIX C FROM X-Y TO 1-2 COORD. SYSTEM C ANGLE=THETA(K)*PI/180. TR(1,1)=COS(ANGLE)**2 TR(2,2)=TR(1,1) TR(1,2)=SIN(ANGLE)**2 TR(3,3)=TR(1,1)-TR(1,2) TR(1,3)=2.*SIN(ANGLE)*COS(ANGLE) TR(2,3)=-TR(1,3) TR(2,1)=TR(1,2) TR(3,1)=-TR(1,3)/2. TR(3,2)=-TR(2,3)/2. C C TRANSFORM X-Y STRAINS INTO 1-2 COORD C EP(3)=EP(3)/2. CALL TXY12(EP,EP12,TR) EP(3)=EP(3)*2. EP12(3)=EP12(3)*2. C C TRANSFORM X-Y STRESSES INTO 1-2 COORD. C CALL TXY12(SIG,SIG12,TR) DO 370 I=1,3 EPF(I,K)=EP(I) EPF12(I,K)=EP12(I) SIGF(I,K)=SIG(I) SIGF12(I,K)=SIG12(I) 370 CONTINUE 380 CONTINUE C C GENERATE STRESS THRU THE THICKNESS C CALL SIGMAZ(TH,NL,T,SIGF,NP,X,SIGZ) C C WRITE OUT FINAL STRESSES C WRITE(6,12)HEAD WRITE(6,82)LDCASE,SFT,OPT,DW,(U(I),I=1,3),MSW,(U(I),I=4,6),(V(I),I *=1,6) 82 FORMAT(30X,'LOAD CASE NO.=',I2,//,' STRESS-FREE TEMP.=',1PE9.2,' *OPERATING TEMP.=',E9.2,' MOIST. CONT.=',E9.2,//,' MECHANICALIMID *PLANE LAMINATE LOADS',/,' LOADING I NX=',E9.2,' NY=',E9.2,' * NXY=',E9.2,/,' MSW=',I2,' I MX=',E9.2,' MY=',E9.2,' MXY *=',E9.2,/,11X,'IMIDPLANE LAMINATE STRAINS AND CURVATURES',/,11X,'I * EXO=',E9.2,' EYO=',E9.2,' EXYO=',E9.2,/,' PLY I KX=',E *9.2,' KY=',E9.2,' KXY=',E9.2,/,' ANGLE',/,' PLY SIG L * SIG T TAU LT SIG X SIG Y TAU XY SIGZ') DO 400 K=1,NL WRITE(6,84)SIGZ(4*K-3) 84 FORMAT(1X,68('-'),1X,1PE9.2) WRITE(6,86)K,SIGZ(4*K-2) 86 FORMAT(2X,I2,66X,1PE9.2) WRITE(6,88)THETA(K),(SIGF12(I,K),I=1,3),(SIGF(I,K),I=1,3),SIGZ(4*K *-1) 88 FORMAT(1X,F5.1,1X,3(1X,1PE9.2),3(1X,E9.2),3X,E9.2) WRITE(6,86)K,SIGZ(4*K) C C GENERATE STRESS ARRAYS TO BE PLOTTED C DO 390 I=1,4 SIGX(4*K-4+I)=SIGF(1,K) SIGY(4*K-4+I)=SIGF(2,K) SIGXY(4*K-4+I)=SIGF(3,K) 390 CONTINUE 400 CONTINUE SIGX(4*NL+1)=SIGF(1,NL) SIGY(4*NL+1)=SIGF(2,NL) SIGXY(4*NL+1)=SIGF(3,NL) WRITE(6,90) 90 FORMAT(1X,68('-')) C C PLOT STRESSES C CALL GRAPH1(NP,X,SIGX,SIGY,SIGXY,SIGZ,ANGL,HEAD,LDCASE) WRITE(6,999) C C WRITE OUT FINAL STRAINS C WRITE(6,12)HEAD WRITE(6,92)LDCASE,SFT,OPT,DW,(U(I),I=1,3),MSW,(U(I),I=4,6),(V(I),I *=1,6) 92 FORMAT(30X,'LOAD CASE NO.=',I2,//,' STRESS-FREE TEMP.=',1PE9.2,' *OPERATING TEMP.=',E9.2,' MOIST. CONT.=',E9.2,//,' MECHANICALIMIDP *LANE LAMINATE LOADS ',/,' LOADING I NX=',E9.2,' NY=',E9.2,' * NXY=',E9.2,/,' MSW=',I2,' I MX=',E9.2,' MY=',E9.2,' MXY *=',E9.2,/,11X,'IMIDPLANE LAMINATE STRAINS AND CURVATURES',/,11X,'I * EXO=',E9.2,' EYO=',E9.2,' EXYO=',E9.2,/,' PLY I KX=',E *9.2,' KY=',E9.2,' KXY=',E9.2,/,' ANGLE',/,' PLY STRN1 * STRN2 STRN12 STRNX STRNY STRNXY') DO 420 K=1,NL WRITE(6,94) 94 FORMAT(1X,68('-')) WRITE(6,96)K 96 FORMAT(2X,I2) WRITE(6,98)THETA(K),(EPF12(I,K),I=1,3),(EPF(I,K),I=1,3) 98 FORMAT(1X,F5.1,1X,3(1X,1PE9.2),3(1X,E9.2)) WRITE(6,96)K 420 CONTINUE 430 CONTINUE 999 FORMAT(1H1) 888 STOP END SUBROUTINE SS(IFLAG,STRN,STRS,Q,S) DIMENSION STRN(3),STRS(3),Q(3,3),S(3,3) C C IFLAG=1 - COMPUTES STRESSES GIVEN STRAINS C IFLAG=-1 - COMPUTES STRAINS GIVEN STRESSES C IF(IFLAG.LT.1)GO TO 1 DO 200 I=1,3 STRESS=0.0 DO 100 J=1,3 STRESS=STRESS+Q(I,J)*STRN(J) 100 CONTINUE STRS(I)=STRESS 200 CONTINUE GO TO 2 1 DO 400 I=1,3 STRAIN=0.0 DO 300 J=1,3 STRAIN=STRAIN+S(I,J)*STRS(J) 300 CONTINUE STRN(I)=STRAIN 400 CONTINUE 2 RETURN END SUBROUTINE LMSS(MSW,R,W,U,V,RM,SM,CM,BM) DIMENSION U(6),V(6),SM(6,6),CM(6,6),R(6),W(6),RM(6,6),BM(6,6) C C IN LAMINATE MIDPLANE C MSW=+1 CALCULATES STRAINS-CURVATURES GIVEN LOADS C MSW=0 CALCULATES STRAINS & MOMENTS GIVEN STRESSES & CURVATUTES C MSW=-1 CALCULATES LOADS GIVEN STRAINS-CURVATURES C MSW=2 CALCULATES STRESSES & CURVATURES GIVEN STRAINS & MOMENTS C IF(MSW.EQ.2)GO TO 4 IF(MSW.EQ.0)GO TO 3 IF(MSW.EQ.-1)GO TO 1 DO 200 I=1,6 VALUE=0.0 DO 100 J=1,6 VALUE=VALUE+CM(I,J)*U(J) 100 CONTINUE V(I)=VALUE 200 CONTINUE GO TO 2 1 DO 400 I=1,6 VALUE=0.0 DO 300 J=1,6 VALUE=VALUE+SM(I,J)*V(J) 300 CONTINUE U(I)=VALUE 400 CONTINUE GO TO 2 3 DO 600 I=1,6 VALUE=0.0 DO 500 J=1,6 VALUE=VALUE+RM(I,J)*R(J) 500 CONTINUE W(I)=VALUE 600 CONTINUE GO TO 2 4 DO 800 I=1,6 VALUE=0.0 DO 700 J=1,6 VALUE=VALUE+BM(I,J)*W(J) 700 CONTINUE R(I)=VALUE 800 CONTINUE 2 RETURN END SUBROUTINE MINV(X,XINV) DIMENSION X(3,3),XINV(3,3) DET=X(1,1)*(X(2,2)*X(3,3)-X(2,3)*X(3,2))+X(1,2)*(X(2,3)*X(3,1)-X(2 *,1)*X(3,3))+X(1,3)*(X(2,1)*X(3,2)-X(2,2)*X(3,1)) XINV(1,1)=(X(2,2)*X(3,3)-X(2,3)*X(3,2))/DET XINV(2,2)=(X(1,1)*X(3,3)-X(1,3)*X(3,1))/DET XINV(1,2)=(X(1,3)*X(3,2)-X(1,2)*X(3,3))/DET XINV(2,1)=(X(2,3)*X(3,1)-X(2,1)*X(3,3))/DET XINV(3,3)=(X(1,1)*X(2,2)-X(1,2)*X(2,1))/DET XINV(1,3)=(X(1,2)*X(2,3)-X(1,3)*X(2,2))/DET XINV(3,1)=(X(2,1)*X(3,2)-X(2,2)*X(3,1))/DET XINV(2,3)=(X(1,3)*X(2,1)-X(1,1)*X(2,3))/DET XINV(3,2)=(X(1,2)*X(3,1)-X(1,1)*X(3,2))/DET RETURN END SUBROUTINE SUBT(X3,X1,X2) DIMENSION X1(3,3),X2(3,3),X3(3,3) DO 100 I=1,3 DO 100 J=1,3 X3(I,J)=X1(I,J)-X2(I,J) 100 CONTINUE RETURN END SUBROUTINE MULT(X3,X1,X2) DIMENSION X1(3,3),X2(3,3),X3(3,3) DO 100 I=1,3 DO 100 J=1,3 X3(I,J)=0.0 100 CONTINUE DO 200 I=1,3 DO 200 J=1,3 DO 200 K=1,3 X3(I,J)=X1(I,K)*X2(K,J)+X3(I,J) 200 CONTINUE RETURN END SUBROUTINE TXY12(A,B,T) DIMENSION A(3),B(3),T(3,3) DO 200 I=1,3 VALUE=0.0 DO 100 J=1,3 VALUE=VALUE+T(I,J)*A(J) 100 CONTINUE B(I)=VALUE 200 CONTINUE RETURN END SUBROUTINE SIGMAZ(TH,NL,T,SIGY,NP,X,SIGZ) DIMENSION T(16),SIGY(3,16),X(65),SIGZ(65) REAL MZ X(1)=0.0 SIGZ(1)=0.0 NP=NL*4+1 DO 200 J=2,NP MZ=0.0 TL=0.0 DO 100 I=2,J K=(J-I+4)/4 TT=TL+T(K)/8.0 MZ=MZ+TT*SIGY(2,K)*T(K)/4.0 TL=TL+T(K)/4.0 100 CONTINUE X(J)=TL SIGZ(J)=45.0*MZ/14.0/(TH/2.0)**2 200 CONTINUE RETURN END SUBROUTINE GRAPH1(N,X,SIGX,SIGY,SIGXY,SIGZ,ANGL,HEAD,LDCASE) DIMENSION X(N),SIGX(N),SIGY(N),SIGXY(N),SIGZ(N),ZERO(65) CHARACTER ANGL(5,16)*1 CHARACTER HEAD(80)*1 CHARACTER ARRAY(71)*1 CHARACTER BLANK*1,H*1,VU*1,C*1,SX*1,SY*1,SXY*1,SZ*1,ZL*1 DATA BLANK/1H /,H/1H-/,VU/1H!/,C/1H+/,SX/1HX/,SY/1HY/,SXY/1HS/,SZ/ *1HZ/,ZL/1H0/ C C GENERATE ZERO ARRAY C DO 30 I=1,N ZERO(I)=0.0 30 CONTINUE C C MAXIMUM AND MINIMUM VALUES CHOOSEN C READ(5,42)ICHOSE 42 FORMAT(I2) IF(ICHOSE.EQ.1)READ(5,43)TAX,TIN 43 FORMAT(2(11X,E10.3,1X)) C C SCAN ARRAYS FOR MAXIMUM AND MINIMUM VALUES C IF(ICHOSE.EQ.1)GO TO 41 TAX=SIGX(1) TIN=SIGX(1) DO 40 ISTRS=1,4 DO 40 I=2,N IF(ISTRS.EQ.1)VAR=SIGX(I) IF(ISTRS.EQ.2)VAR=SIGY(I) IF(ISTRS.EQ.3)VAR=SIGXY(I) IF(ISTRS.EQ.4)VAR=SIGZ(I) IF(VAR.GT.TAX)TAX=VAR IF(VAR.LT.TIN)TIN=VAR 40 CONTINUE 41 CONTINUE C C PROVIDED FOR EDGE MARGINS BY APPROPRIATELY ROUNDING C OFF MAXIMUM AND MINIMUM VALUES C DO 120 IT=1,2 VAL=TIN IF(IT.EQ.1)VAL=TAX IF(VAL.EQ.0.0) GO TO 110 SIGN=1. IF(VAL.LT.0.0)SIGN=-1. ADD=+3. IF((VAL.GT.0.0).AND.(IT.EQ.2))ADD=-3. VAL=SIGN*VAL IF(VAL.GE.100.) GO TO 70 IF(VAL.GE.10.) GO TO 100 DO 50 I=1,20 K=VAL*10**I IF(K.GE.10.) GO TO 60 50 CONTINUE 60 VAL=SIGN*(K+ADD)/10**I GO TO 110 70 DO 80 I=1,20 K=VAL/10**I IF(K.LT.100.) GO TO 90 80 CONTINUE 90 VAL=SIGN*(K+ADD)*10**I GO TO 110 100 K=VAL+ADD VAL=SIGN*K 110 IF(IT.EQ.1)TAXG=VAL TING=VAL 120 CONTINUE DY=(TAXG-TING)/60. C C ASSUME X(N) ARRAY IS EVENLY INCREMENTED BY DX C DX=X(2)-X(1) C C PRINT GRAPH HEAD C WRITE(6,999) 999 FORMAT(1H1) WRITE(6,46)HEAD 46 FORMAT(80A1,/) WRITE(6,125) 125 FORMAT(15X,'SYMBOL SYMBOL CORRESPONDING',/,14X,'PRIORITY P *LOTTED',7X,'STRESS',/,11X,'LOW 5 -------- X ---- STRESS IN X DIR *ECTION',/,17X,'4 -------- Y ---- STRESS IN Y DIRECTION',/,17X,'3 - *------- S ---- SHEAR STRESS IN X-Y PLANE',/,17X,'2 -------- Z ---- * STRESS IN Z DIRECTION',/,11X,'HIGH 1 -------- 0 ---- ZERO LINE', */) WRITE(6,130)LDCASE,TIN,TAX,DY,DX,TING,TAXG 130 FORMAT(25X,'STRESS PLOT: LOAD CASE NO.=',I2,//,' MINIMUM STRESS=' *,E13.5,19X,' THICKNESS THICKNESS T(I)',/,' MAXIMUM STRESS=' *,E13.5,22X,'INCREMENT',8X,'LISTED DOWN',/,' STRESS INCREMENT=' *,E13.5,14X,E13.5,12X,'RIGHT',/,74X,'COLUMN',/,' I',9X,'MIN=',E13.5 *,2X,'<---------STRESS-------->',E13.5,'=MAX',4X,'I',/,5X,'PLY V' *,59X,'V',4X,'I',/,4X,' TYPE ',61('T'),4X,'I') C C FILL ARRAY TO BE PRINTED WITH APPROPRIATE GRAPH SYMBOLS C IMARK=4 C C LOOP THRU ARRAY Y(I) C DO 280 M=1,N DO 140 I=1,71 ARRAY(I)=BLANK 140 CONTINUE DO 150 I=1,8 IP=10*(I-1)+1 ARRAY(IP)=VU 150 CONTINUE IF(IMARK.LT.2.OR.IMARK.GT.2)GO TO 152 MK=(M+1)/4 DO 151 I=4,8 ARRAY(I)=ANGL(I-3,MK) 151 CONTINUE 152 CONTINUE IF(IMARK.LT.4) GO TO 180 DO 160 I=2,70 ARRAY(I)=H 160 CONTINUE DO 170 I=1,6 IC=10*I+1 ARRAY(IC)=C 170 CONTINUE IMARK=0 180 IMARK=IMARK+1 C C SCALE ARRAYS TO BE PLOTTED FROM 1 TO 70 C DO 111 ISTRS=1,5 IF(ISTRS.EQ.1)VAR=SIGX(M) IF(ISTRS.EQ.2)VAR=SIGY(M) IF(ISTRS.EQ.3)VAR=SIGXY(M) IF(ISTRS.EQ.4)VAR=SIGZ(M) IF(ISTRS.EQ.5)VAR=ZERO(M) IF(TING)190,230,240 190 P=-61.*TING/(TAXG-TING) IF(VAR) 200,210,220 200 S=P-P*VAR/TING GO TO 250 210 S=P GO TO 250 220 S=((61.-P)*VAR/TAXG)+P GO TO 250 230 S=(61./TAXG)*VAR GO TO 250 240 S=(61./(TAXG-TING))*(VAR-TING) 250 S=S+11 I=S Q=I I=I+1 IF(I.LE.11)GO TO 260 IF(I.GE.72)GO TO 260 IF(ISTRS.EQ.1)ARRAY(I)=SX IF(ISTRS.EQ.2)ARRAY(I)=SY IF(ISTRS.EQ.3)ARRAY(I)=SXY IF(ISTRS.EQ.4)ARRAY(I)=SZ IF(ISTRS.EQ.5)ARRAY(I)=ZL 260 CONTINUE 111 CONTINUE C C PLOT ARRAY AND PRINT VALUE FOR INDEPENDENT VARIABLE X(I) C WRITE(6,270) X(M),ARRAY 270 FORMAT(71X,1PE9.2,T2,71A1) 280 CONTINUE RETURN END