PROGRAM Polar C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * PROGRAM GRAPH: CALCULATES POLAR DIAGRAMS FOR TRANSVERSELY * C * ISOTROPIC (HEXAGONAL) SYMMETRIC FIBER REINFORCED MATL'S. * C * USING EQUATIONS LISTED IN PAPER PRESENTED AT THE 2ND U.S. * C * - JAPAN SYMPOSIUM ON COMPOSITE MATERIALS, 1983, N.A.S.A. * C * LANGELY RES. CENTER, HAMPTON, VA.. DIAGRAMS WERE * C * ORIGINALLY DRAWN WITH D.I.S.S.P.L.A. SOFTWARE WRITTEN BY * C * I.S.S.C.O. CORP. SAN DIEGO, CA WHICH HAVE BEEN COMMENTED * C * OUT. FOR THE CRCD DATA IS STORED IN THREE FILES THAT WILL * C * BE PLOTTED USING PV-WAVE SOFTWARE. IMSL ROUTINES HAVE * C * ALSO BEEN CHANGED * C * FROM CALL LINV3F(S,NULL,1,6,6,-1.,D2,WK,IER) * C * TO call linds(6,c,6,s,6) * c * TO AVOID LINKING TO IMSL CALL INV(C,6,S) IS USED INSTEAD * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DIMENSION XNU13(361),YNU13(361),XNU32(361),YNU32(361) DIMENSION XG23(361),YG23(361),XG13(361),YG13(361) DIMENSION XE3(361),YE3(361),XT3(361),YT3(361),XK3(361),YK3(361) DIMENSION S(6,6),C(6,6) real negNU32(361,3), posNU32(361,3) REAL vfk(3),e3(361,3),g23(361,3),NU13,NU32(361,3),K3 c c Assign unit numbers to filenames c OPEN(5,FILE='polar.DAT',STATUS='OLD',ERR=888) OPEN(6,FILE='polar.out',STATUS='UNKNOWN',ERR=888) OPEN(7,FILE='e3.dat',STATUS='UNKNOWN',ERR=888) OPEN(8,FILE='g23.dat',STATUS='UNKNOWN',ERR=888) OPEN(9,FILE='nu32.dat',STATUS='UNKNOWN',ERR=888) C C CALC. RADIANS TO DEGREES CONVERSION FACTOR C PI=ACOS(-1.) CONV=PI/180. C C SET I.D. INFO. WITH ALL PLOT COMMANDS SENT TO A COMPRESSED C DIRECT ACCESS FILE NAMED " GPLT " C c CALL ID("R D KRIZ;DIV.562;X-3547",100) c CALL COMPRS(4HCONT,4) C C C READ MATERIAL TYPE AND PROPERTIES FROM INPUT FILE C CALL IN C C READ THREE FIBER VOLUME FRACTIONS, vfk(3) C read(5,13)(vfk(i),i=1,3) 13 format(23x,3f5.2) c c Write fiber volume fractions to polar.out c write(6,11)vfk(1),vfk(2),vfk(3) 11 format(' Three Fiber Volume Fractions: ',/,3(1x,f5.2)) c c Zero Poisson ratio arrays c do 50 k=1,3 do 50 j=1,361 posNU32(j,k)=0.0 negNU32(j,k)=0.0 50 continue c c Loop thru the three fiber volume fractions c do 100 k=1,3 VF=vfk(k) C C CALCULATE MATERIAL COMPLIANCES S(I,J) FROM Z. HASHIN EQN.S C GIVEN FIBER - EPOXY ELASTIC CONSTANTS AND FIBER VOLUME FRACTION,VF C ALSO CALCULATE MATERIAL DENSITY, RHO C CALL PROP(VF,S,C) C C CHOOSE TILT ANGLE, ALPHA (DEGREES) C ALPHA=00.0 AR=ALPHA*PI/180. CA=COS(AR) SA=SIN(AR) C C WRITE HEADER ON OUTPUT FILE WRITE(6,676)ALPHA,VF 676 FORMAT(' TILT ANGLE ALPHA=',E10.3,' VOLUME FRACTION=',E10.3,//, *5X,'THETA',9X,'E3',10X,'T3',9X,'G23',9X,'G13',10X,'K3',8X,'NU32', *9X,'NU13') C C CALCULATE POLAR DIAGRAMS FOR YOUNGS MODULI, POISSONS RATIOS, SHEAR C MODULI, AND LINEAR COMPRESSIBILITY ELASTIC CONSTANTS. C DO 100 J=1,361 C THETA=FLOAT(J)-1. THETAR=THETA*PI/180. CN=COS(THETAR) SN=SIN(THETAR) C C YOUNG'S MODULUS EQN. NO. (?E3) C E3(j,k)=1./(S(1,1)+(2.*S(1,3)+S(4,4)-2.*S(1,1))*CN*CN * +(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*CN**4) XE3(J)=E3(j,k)*SN YE3(J)=E3(j,k)*CN C C TORSIONAL MODULUS EQN. NO. (?T3) C T3=1./(S(4,4)+(S(1,1)-S(1,2)-S(4,4)/2.)*SN*SN * +2.*(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*CN*CN*SN*SN) XT3(J)=T3*SN YT3(J)=T3*CN C C SHEAR MODULUS IN THE X2-X3 PLANE EQN. NO. (?G23) C G23(j,k)=1./(S(4,4)+(S(1,1)-S(1,2)-S(4,4)/2.)*SN*SN *+2.*(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*CN*CN*SN*SN) XG23(J)=G23(j,k)*SN YG23(J)=G23(j,k)*CN C C SHEAR MODULUS IN THE X1-X3 PLANE EQN. NO. (?G13) C G13=1./(S(5,5)*CN*CN+2.*(S(1,1)-S(1,2))*SN*SN) XG13(J)=G13*SN YG13(J)=G13*CN C C LINEAR COMPRESSIBILITY EQN. NO. (?K3) C K3=S(1,1)+S(1,2)+S(1,3)-(S(1,1)-S(3,3)+S(1,2)-S(1,3))*CN*CN XK3(J)=K3*SN YK3(J)=K3*CN C C POISSONS RATIO NU32 EQN. NO. (?32) C NU32(j,k)=(-S(1,3)-(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*CN*CN*SN*SN) * /(S(1,1)+(2.*S(1,3)+S(4,4)-2.*S(1,1))*CN*CN * +(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*CN**4) XNU32(J)=NU32(j,k)*SN YNU32(J)=NU32(j,k)*CN C C POISSONS RATIO NU13 EQN. NO. (?13) C NU13=-1.*((S(1,1)+S(3,3)-S(4,4))*CA*CA*SA*SA*CN*CN * +S(1,2)*CA*CA*SN*SN * +S(1,3)*(CN*CN*(CA**4+SA**4)+SA*SA*SN*SN)) * /(S(1,1)+(2.*S(1,3)+S(4,4)-2.*S(1,1))*SA*SA * +(S(1,1)+S(3,3)-2.*S(1,3)-S(4,4))*SA**4) XNU13(J)=NU13*SN YNU13(J)=NU13*CN C WRITE(6,120)THETA,E3(j,k),T3,G23(j,k),G13,K3,NU32(j,k),NU13 120 FORMAT(8(2X,E10.3)) c c Sort out negative and positive Poisson's ratios c if(nu32(j,k).gt.0.0)posNU32(j,k)=NU32(j,k) if(nu32(j,k).lt.0.0)negNU32(j,k)=NU32(j,k) 100 continue c c Write results to files for polar plots c do 200 jj=1,361 THETA=FLOAT(jj)-1. THETAR=THETA*PI/180. WRITE(7,121)THETAR,(E3(jj,k),k=1,3) WRITE(8,121)THETAR,(G23(jj,k),k=1,3) WRITE(9,122)THETAR,(posNU32(jj,k),k=1,3),(negNU32(jj,k),k=1,3) 121 FORMAT(1X,E11.4,3(1X,E11.4)) 122 FORMAT(1X,E11.4,3(1X,E11.4),3(1X,E11.4)) 200 continue 888 stop end c c********1*********2*********3*********4*********5*********6*********7** SUBROUTINE IN C COMMON /HEAD1/MTYPE(50),ICE(25),ICF(50),IEM(30),INUM(30),IU(20) COMMON /HEAD2/IEFL(40),IEFT(40),IGFLT(40),IGFTT(40),IKFTT(40) COMMON /mf/EM,NUM,EFL,EFT,GFLT,GFTT,KFTT C REAL NUM,KFTT C C READ MATERIAL TYPE C READ(5,1)(MTYPE(I),I=1,50) 1 FORMAT(50A1) C C READ COMMENT AND ELASTIC EPOXY PROPERTIES C READ(5,2)(ICE(I),I=1,25),(IU(I),I=1,20),(IEM(I),I=1,30),EM, * (INUM(I),I=1,30),NUM 2 FORMAT(25A1,20A1,/,30A1,E12.5,/,30A1,E12.5) C C READ COMMENT AND ELASTIC FIBER PROPERTIES C READ(5,3)(ICF(I),I=1,50),(IEFL(I),I=1,40),EFL, * (IEFT(I),I=1,40),EFT, * (IGFLT(I),I=1,40),GFLT, * (IGFTT(I),I=1,40),GFTT, * (IKFTT(I),I=1,40),KFTT 3 FORMAT(50A1,/,4(40A1,E12.5,/),40A1,E12.5) C RETURN END c c********1*********2*********3*********4*********5*********6*********7** SUBROUTINE PROP(VF,S,C) COMMON /HEAD1/MTYPE(50),ICE(25),ICF(50),IEM(30),INUM(30),IU(20) COMMON /HEAD2/IEFL(40),IEFT(40),IGFLT(40),IGFTT(40),IKFTT(40) COMMON /mf/EM,NUM,EFL,EFT,GFLT,GFTT,KFTT DIMENSION C(6,6),S(6,6) REAL KFTT,NUFTT,NUFLT,NUM,LM,KM,NU12,NU13,NU23 REAL NU21,NU31,NU32,K23 C C CALCULATE OTHER EPOXY ELASTIC PROPERTIES C GM=EM/(2.*(1.+NUM)) LM=2.*GM*NUM/(1.-2.*NUM) KM=LM+GM C C CALCULATE OTHER FIBER ELASTIC PROPERTIES C NUFTT=EFT/(2.*GFTT)-1. CAL=EFL*(2.*KFTT-EFT-2.*KFTT*NUFTT)/(4.*KFTT*EFT) NUFLT=SQRT(CAL) C C WRITE EPOXY ELASTIC PROPERTIES ON OUTPUT FILE C WRITE(6,999) WRITE(6,4)(ICE(I),I=1,25),(IU(I),I=1,20),EM,GM,NUM 4 FORMAT(3X,25A1,20A1,/, * ' MATRIX YOUNGS MODULUS, EM=',E11.4,/, * ' MATRIX SHEAR MODULUS, GM=',E11.4,/, * ' MATRIX POISSON RATIO, NUM=',E11.4,/) C C WRITE FIBER ELASTIC PROPERTIES ON OUTPUT FILE C WRITE(6,5)(ICF(I),I=1,25),(IU(I),I=1,20),EFL,EFT,GFLT,GFTT,KFTT 5 FORMAT(3X,25A1,20A1,/, * ' FIBER LONGITUDINAL YOUNGS MODULUS, EFL=',E11.4,/, * ' FIBER TRANSVERSE YOUNGS MODULUS, EFT=',E11.4,/, * ' FIBER SHEAR MODULUS L-T PLANE, GFLT=',E11.4,/, * ' FIBER SHEAR MODULUS T PLANE, GFTT=',E11.4,/, * ' FIBER PLANE STRAIN BULK MODULUS, KFTT=',E11.4) C C WRITE MATERIAL TYPE C WRITE(6,7)(MTYPE(I),I=1,50),VF 7 FORMAT(/,2X,'MATERIAL: ',50A1,/,' FIBER VOL. FRACTION=',F5.2,//) C C EQUATIONS DEVELOPED BY ZVI HASHIN C REF: THEORY OF FIBER REINFORCED MATERIALS, NASA CR-1974, MARCH 74 C X1 FIBER AXIS ; X2-X3 PLANE TRANSVERSELY ISOTROPIC C C VM=1.-VF C E11=EFL*VF+EM*VM C GAMMA=GFTT/GM C BETA1=1./(3.-4.*NUM) C BETA2=KFTT/(KFTT+2.*GFTT) C ALPHA1=(BETA1-GAMMA*BETA2)/(1.+GAMMA*BETA2) C ALPHA2=(GAMMA+BETA1)/(GAMMA-1.) C G23=GM*((1.+ALPHA1*VF**3)*(ALPHA2+BETA1*VF)-3.*VF*VM**2*BETA1**2)/ C *((1.+ALPHA1*VF**3)*(ALPHA2-VF)-3.*VF*VM**2*BETA1**2) C G12=GM+VF/(1./(GFLT-GM)+VM/(2.*GM)) C G13=G12 C L1=2.*(1.-NUM) C L2=1.+EM*(1.-NUFLT-2.*NUFLT**2)/(EFL*(1.+NUM)) C NU12=(VF*NUFLT*L1+VM*NUM*L2)/(VF*L1+VM*L2) C NU13=NU12 C K23=KM+VF/(1./(KFTT-KM)+VM/(KM+GM)) C PSI=1.+4.*K23*NU12**2/E11 C NU23=(K23-PSI*G23)/(K23+PSI*G23) C E22=4.*K23*G23/(K23+PSI*G23) C E33=E22 C NU21=E22*NU12/E11 C NU31=E33*NU13/E11 C NU32=E33*NU23/E22 C C EQUATIONS DEVELOPED BY LEDBETTER-DATTA C VM=1.-VF E11=EFL*VF+EM*VM+(4.*VF*VM*(NUFLT-NUM)**2)/(1./GM+VF/KM+VM/KFTT) G23=GM+(GM*2.*VF*(KM+GM)*(GFTT-GM))/(2.*GM*(KM+GM)+ * VM*(KM+2.*GM)*(GFTT-GM)) G12=GM+(2.*VF*GM*(GFLT-GM))/((1.+VF)*GM+VM*GFLT) G13=G12 NU12=VM*NUM+VF*NUFLT+(VF*VM*(NUFLT-NUM)*(1./KM-1./KFTT))/ * (1./GM+VF/KM+VM/KFTT) NU13=NU12 K23=KM+(VF*(KM+GM)*(KFTT-KM))/((KFTT+GM)-VF*(KFTT-KM)) PSI=1.+4.*K23*NU12**2/E11 NU23=(K23-PSI*G23)/(K23+PSI*G23) E22=4.*K23*G23/(K23+PSI*G23) E33=E22 NU21=E22*NU12/E11 NU31=E33*NU13/E11 NU32=E33*NU23/E22 C C WRITE ENGINEERING ELASTIC CONSTANTS ON OUTPUT FILE C WRITE(6,8)(IU(I),I=1,20),E11,E22,(IU(I),I=1,20),G12,G23,NU12, *NU23 8 FORMAT(3X,'ENGINEERING ELASTIC PROPERTIES',/, * ' YOUNGS MODULI',20A1,/, * ' PARALLEL TO FIBER LONGITUDINAL AXIS, EL=',E10.3,/, * ' TRANSVERSE TO FIBER, ET=',E10.3,/, * ' SHEAR MODULI',20A1,/, * ' LONGITUDINAL-TRANSVERSE PLANE, GLT=',E10.3,/, * ' TRANSVERSE PLANE, GTT=',E10.3,/, * ' POISSONS RATIOS',/, * ' LONGITUDINAL-TRANSVERSE AXIS, NULT=',E10.3,/, * ' TRANSVERSE PLANE, NUTT=',E10.3,/) C C ZERO ALL STIFFNESS TERMS C DO 10 I=1,6 DO 10 J=1,6 C(I,J)=0.0 10 CONTINUE C C CALC. TRANSVERSELY ISOTROPIC STIFFNESSES, C(I,J) (HEXAGONAL SYMM.) C X3 FIBER AXIS ; X1-X2 PLANE TRANSVERSELY ISOTROPIC C DEL=1.-NU12*NU21-NU23*NU32-NU31*NU13-2.*NU21*NU32*NU13 C(3,3)=(1.-NU23*NU32)*E11/DEL C(1,1)=(1.-NU13*NU31)*E22/DEL C(2,2)=(1.-NU12*NU21)*E33/DEL C(1,3)=(NU21+NU31*NU23)*E11/DEL C(3,1)=C(1,3) C(2,3)=(NU31+NU21*NU32)*E11/DEL C(3,2)=C(2,3) C(1,2)=(NU32+NU12*NU31)*E22/DEL C(2,1)=C(1,2) C(6,6)=G23 C(4,4)=G13 C(5,5)=G12 C C WRITE STIFFNESSES ON OUTPUT FILE C WRITE(6,15)(IU(I),I=1,20) 15 FORMAT(18X,'STIFFNESSES',20A1,/) C WRITE(6,20) 20 FORMAT(15X,'C11 C12 C13 C14 C15 C16',/, * 15X,' C22 C23 C24 C25 C26',/, * 15X,' C33 C34 C35 C36',/, * 15X,' C44 C45 C46',/, * 15X,' SYMMETRIC C55 C56',/, * 15X,' C66',/) WRITE(6,21)(C(1,J),J=1,6) 21 FORMAT(1X,6(1X,E11.4)) WRITE(6,22)(C(2,J),J=2,6) 22 FORMAT(13X,5(1X,E11.4)) WRITE(6,23)(C(3,J),J=3,6) 23 FORMAT(25X,4(1X,E11.4)) WRITE(6,24)(C(4,J),J=4,6) 24 FORMAT(37X,3(1X,E11.4)) WRITE(6,25)C(5,5),C(5,6),C(6,6) 25 FORMAT(13X,' SYMMETRIC',24X,2(1X,E11.4),/,62X,E11.4) C C CALCULATE COMPLIANCES S(I,J) BY INVERTING STIFFNESSES C(I,J) C BY USING IMSL SUBROUTINE (LINV3F) old NIST routine C BY USING IMSL SUBROUTINE (LINDS) new VT routine C c DO 29 I=1,6 c DO 29 J=1,6 c S(I,J)=C(I,J) c 29 CONTINUE C c DO 32 I=1,6 c 32 WRITE(6,31)(S(I,J),J=1,6) c 31 FORMAT(1X,6(1X,E11.4)) c CALL LINV3F(S,NULL,1,6,6,-1.,D2,WK,IER) c call linds(6,c,6,s,6) call inv(c,6,s) C C WRITE COMPLIANCES S(I,J) ON OUTPUT FILE C WRITE(6,30) 30 FORMAT(//,25X,' COMPLIANCES',/, * 15X,'S11 S12 S13 S14 S15 S16',/, * 15X,' S22 S23 S24 S25 S26',/, * 15X,' S33 S34 S35 S36',/, * 15X,' S44 S45 S46',/, * 15X,' SYMMETRIC S55 S56',/, * 15X,' S66') WRITE(6,21)(S(1,J),J=1,6) WRITE(6,22)(S(2,J),J=2,6) WRITE(6,23)(S(3,J),J=3,6) WRITE(6,24)(S(4,J),J=4,6) WRITE(6,25)S(5,5),S(5,6),S(6,6) C WRITE(6,999) 999 FORMAT(1H1) C RETURN END c subroutine inv(a,n,b) dimension a(n,n),aa(n,n),b(n,n),indx(n) do 12 i=1,n do 11 j=1,n 11 aa(i,j)=a(i,j) 12 continue do 14 i=1,n do 13 j=i,n 13 b(i,j)=0.0 14 b(i,i)=1.0 call ludcmp(aa,n,n,indx,d) do 15 j=1,n 15 call lubksb(aa,n,n,indx,b(1,j)) return end c subroutine ludcmp(a,n,np,indx,d) parameter (nmax=100,tiny=1.0E-20) dimension a(np,np),indx(n),vv(nmax) d=1. do 12 i=1,n aamax=0. do 11 j=1,n if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if(aamax.eq.0.) pause 'singular matrix' vv(i)=1./aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0. do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if(dum.ge.aamax)then imax=i aamax=dum endif 16 continue if(j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.0)a(j,j)=tiny if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return end c subroutine lubksb(a,n,np,indx,b) dimension a(np,np),indx(n),b(n) ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0) then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return end