PROGRAM edgesing C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C* PROGRAM : CALCULATES THE SINGULARITY OF THE STRESS FIELD AT THE * C* INTERFACE AND STRESS FREE EDGE OF A LAMINATE WHERE THE DISSIMILAR * C* MATERIALS HAVE ORTHORHOMBIC ELASTIC SYMMETRY. EACH LAYER IS A * C* FIBER-REINFORCED MATERIAL WHOSE FIBER ORIENTATION IS MEASURED BY * C* THE ANGLE THETA FROM THE LOAD AXIS. IN THIS PROGRAM WE USE THE * C* NOTATION OF T.C.T. TING WHERE ALL VARIABLES RELATED TO THE BOTTOM * C* LAYER ARE SUFFIXED BY THE LETTER P. * C* * C* EXAMPLE: STIFFNESSES OF TOP LAYER : CB(I,J) * C* STIFFNESS OF BOTTOM LAYER : CBP(I,J) * C* * C* THIS SOLUTION IS OUTLINED BY T.C.T.TING AND S.C. CHOU, FRACTURE * C* OF COMPOSITE MATERIALS, EDS. G.C.SIH AND V.P.TAMUZS, 1982, PROC. * C* 2ND USA-USSR SYMPOSIUM, LEIGH UNIV., BETHLEHEM, PA, MAR.9-12,1981 * C* * C* NIST LATEST UPDATE: **** 1/21/87 **** * C* (IMSL VERSION 1.1) * C* * C* PROGRAM WRITTEN BY R.D.KRIZ, FRACTURE AND DEFORMATION DIVISION, * C* INSTITUTE FOR MATERIALS SCIENCE AND ENGINEERING, * C* NATIONAL BUREAU OF STANDARDS, BOULDER, COLORADO 80303. * C* * C* Modified with new IMSL Version 2.0 routines **** 12/6/00 **** * C* Ron Kriz, Virginia Tech, Blacksburg, Virginia 24061 * C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INITIALLY ALL VARIABLES ARE DEFINED AS LOGICAL. THIS WILL GENERATE AN C ERROR IF A VARIABLE BELOW IS NOT REDEFINED AS REAL, INTEGER, OR COMPLEX. C HENCE ALL VARIABLES LISTED BELOW MUST BE REDEFINED APPROPRIATELY. C IMPLICIT LOGICAL (A-Z) C COMMON /L1/CB(6,6),P(3),V(3,3) COMMON /L2/CBP(6,6),PP(3),VP(3,3) COMMON /L3/T(3,3,3),TP(3,3,3) C INTEGER I,J,JJ,III,JJJ,L,IK,IR,maxfn C REAL A(6,6),AP(6,6),KT(100),KR(2) C REAL PI,E1,E2,E3,G12,G13,G23,NU21,NU31,NU32,NU12,NU13,NU23, *DEL,C11,C22,C33,C12,C13,C23,C44,C55,C66,THETA,THETAP,THETAR,M,N, *DEN,K,DET,F,CB,CBP,TEST real epirg,pindex,pindexp,errabs,errrel,aaa,bbb REAL S11,S12,S13,S22,S23,S33,S44,S55,S66 C COMPLEX W(6),WP(6),Z(6,6),ZP(6,6),B(3,3),BP(3,3) C COMPLEX FT,FTP,SUM,SUMP,V,VP,P,PP,T,TP,CMI COMPLEX TEST1,TEST2,TEST3,TEST1P,TEST2P,TEST3P c dimension title(48) C EXTERNAL FT,FTP,F C OPEN(5,FILE='edgesing.dat',STATUS='OLD',ERR=888) OPEN(6,FILE='edgesing.out',STATUS='UNKNOWN',ERR=888) OPEN(7,FILE='k.txt',STATUS='UNKNOWN',ERR=888) c c Read one line title c read(5,12)(title(i),i=1,12) 12 format(12a4) c c Write title line to output file c write(6,12)(title(i),i=1,12) C C DEFINE COMPLEX MINUS I "-SQRT(-1.)" CMI=(0.000000000000000,-1.000000000000000) C C GEOMETRIC PI(RADIANS)=180(DEGREES) PI=ACOS(-1.) C C CALCULATE STIFFNESSES FROM ENGINEERING ELASTIC PROPERTIES C (THE FIBER AXIS IS CHOOSEN AS X3) C C ELATIC PROPERTIES FROM REPORT VPI-E-77-16 C C E1=1.453E+06 C E2=E1 C E3=20.0E+06 C G12=0.4872E+06 C G13=0.8394E+06 C G23=G13 C NU21=0.4906 C NU31=0.3056 C NU32=NU31 C C ELASTIC PROPERTY APPROXIMATION AFTER PIPES-PAGANO C c E1=2.100E+06 c E2=E1 c E3=20.0E+06 c G12=0.850E+06 c G13=G12 c G23=G13 c NU21=0.2100 c NU31=NU21 c NU32=NU31 c c Read Engineering properties from data file c read(5,11)e1,e2,e3,g12,g13,g23,nu21,nu31,nu32 11 format(8(36x,e10.3,/),36x,e10.3) c c 1 2 3 4 c 1234567890123456789012345678901234567890 c Young's modulus in 1-direction, E1=+2.100E+06 c Young's modulus in 2-direction, E2=+2.100E+06 c Young's modulus in 3-direction, E3=+2.000E+07 c Shear modulus in 12-plane, G12=+0.850E+06 c Shear modulus in 13-plane, G13=+0.850E+06 c Shear modulus in 23-plane, G23=+0.850E+06 c Poission's ratio in 12-plane, NU21=+0.211E+00 c Poission's ratio in 13-plane, NU31=+0.211E+00 c Poission's ratio in 23-plane, NU32=+0.211E+00 C C CALCULATE OTHER POISSONS RATIOS C NU12=NU21*E1/E2 NU13=NU31*E1/E3 NU23=NU32*E2/E3 C C WRITE ENGINEERING ELASTIC CONSTANTS TO OUTPUT FILE C WRITE(6,251) WRITE(6,241) 241 FORMAT(' ENGINEERING ELASTIC CONSTANTS') WRITE(6,251) WRITE(6,678)E1,E2,E3,G12,G13,G23,NU12,NU13,NU23,NU21,NU31,NU32 678 FORMAT(' E1=',E11.4,' E2=',E11.4,' E3=',E11.4,/, * ' G12=',E11.4,' G13=',E11.4,' G23=',E11.4,/, * ' NU12=',E11.4,' NU13=',E11.4,' NU23=',E11.4,/, * ' NU21=',E11.4,' NU31=',E11.4,' NU32=',E11.4) C C CALCULATE STIFFNESSES CIJ FROM ENGINEERING PROPERTIES C DEL=1.-NU12*NU21-NU13*NU31-NU23*NU32-2.*NU21*NU32*NU13 C11=(1.-NU23*NU32)*E1/DEL C22=(1.-NU13*NU31)*E2/DEL C33=(1.-NU12*NU21)*E3/DEL C12=(NU21+NU31*NU23)*E1/DEL C13=(NU31+NU21*NU32)*E1/DEL C23=(NU32+NU12*NU31)*E2/DEL C44=G23 C55=G13 C66=G12 C C WRITE STIFFNESSES ON OUTPUT FILE C WRITE(6,251) WRITE(6,260) 260 FORMAT(' STIFFNESSES CALCULATED FROM ENGINERRING CONSTANTS') WRITE(6,251) WRITE(6,789)C11,C22,C33,C12,C13,C23,C44,C55,C66 C C CALCULATE COMPLIANCES SIJ FROM ENGINEERING PROPERTIES C S11=1./E1 S22=1./E2 S33=1./E3 S44=1./G23 S55=1./G13 S66=1./G12 S12=-NU21/E2 S13=-NU31/E3 S23=-NU32/E3 C C CALCULATE STIFFNESSES CIJ FROM COMPLIANCES SIJ C DEL=2.*S12*S23*S13+S11*S22*S33-S11*S23**2-S22*S13**2-S33*S12**2 C11=(S22*S33-S23**2)/DEL C22=(S33*S11-S13**2)/DEL C33=(S11*S22-S12**2)/DEL C12=(S13*S23-S12*S33)/DEL C13=(S12*S23-S13*S22)/DEL C23=(S12*S13-S23*S11)/DEL C44=1./S44 C55=1./S55 C66=1./S66 C C WRITE STIFFNESSES TO OUTPUT FILE C WRITE(6,251) WRITE(6,261) 261 FORMAT(' STIFFNESSES CALCULATED FROM COMPLIANCES WHICH WERE', *' CALC. FROM ENGR. CONSTANTS') WRITE(6,251) WRITE(6,789)C11,C22,C33,C12,C13,C23,C44,C55,C66 789 FORMAT(' C11=',E11.4,' C22=',E11.4,' C33=',E11.4,/, * ' C12=',E11.4,' C13=',E11.4,' C23=',E11.4,/, * ' C44=',E11.4,' C55=',E11.4,' C66=',E11.4) C C CLEAR ARRAYS A(I,J), CB(I,J), AP(I,J) CBP(I,J) C DO 2 J=1,6 DO 2 I=1,6 CB(I,J)=0.0 CBP(I,J)=0.0 A(I,J)=0.0 2 AP(I,J)=0.0 C C FIBER ANGLES(DEG.) (TOP LAYER: THETA) (BOTTOM LAYER: THETAP) c c Read fiber orientation in the top and bottom regions c read(5,222)theta,thetap 222 format(40x,f5.1,/,40x,f5.1) c 1 2 3 4 5 c 12345678901234567890123456789012345678901234567890 c Fiber orientation top region, theta=+75.0 c Fiber orientation bottom region, thetap=+00.0 c c THETA=+75.0 c THETAP=+00.0 c c THETA=+15.0 c THETAP=-15.0 C C WRITE OUT FIBER ANGLE FOR BOTH LAYERS C WRITE(6,251) WRITE(6,14)THETA,THETAP 14 FORMAT(' TOP LAYER, ANGLE (DEG.)=',E10.3,/, * ' BOTTOM LAYER, ANGLE (DEG.)=',E10.3,//, *' BELOW ARE STIFFNESSES FOR BOTH LAYERS') WRITE(6,251) C C CALCULATE STIFFNESSES IN THE ROTATED X3-X1 PLANE FOR THE TOP LAYER C THETAR=PI*THETA/180.0 M=COS(THETAR) N=SIN(THETAR) CB(1,1)=M**4*C11+2.*(M*N)**2*(C13+2.*C55)+C33*N**4 CB(1,3)=(N*M)**2*(C11+C33-4.*C55)+(N**4+M**4)*C13 CB(1,2)=M**2*C12+N**2*C23 CB(1,5)=M*N*((C13+2.*C55-C11)*M**2-(C13+2.*C55-C33)*N**2) CB(3,3)=N**4*C11+2.*(N*M)**2*(C13+2.*C55)+M**4*C33 CB(2,3)=N**2*C12+M**2*C23 CB(2,2)=C22 CB(2,5)=M*N*(C23-C12) CB(3,5)=M*N*((C13+2.*C55-C11)*N**2-(C13+2.*C55-C33)*M**2) CB(4,4)=C44*M**2+C66*N**2 CB(4,6)=M*N*(C44-C66) CB(6,6)=C66*M**2+C44*N**2 CB(5,5)=(M*N)**2*(C11+C33-2.*C13)+C55*(M**2-N**2)**2 CB(3,1)=CB(1,3) CB(2,1)=CB(1,2) CB(5,1)=CB(1,5) CB(3,2)=CB(2,3) CB(5,2)=CB(2,5) CB(5,3)=CB(3,5) CB(6,4)=CB(4,6) C C WRITE ROTATED STIFFNESSES FOR TOP LAYER C WRITE(6,251) DO 455 I=1,6 455 WRITE(6,456)(CB(I,J),J=1,6) 456 FORMAT(6(1X,E11.4)) C C CALCULATE STIFFNESSES IN THE ROTATED X3-X1 PLANE FOR BOTTOM LAYER C THETAR=PI*THETAP/180.0 M=COS(THETAR) N=SIN(THETAR) CBP(1,1)=M**4*C11+2.*(M*N)**2*(C13+2.*C55)+C33*N**4 CBP(1,3)=(N*M)**2*(C11+C33-4.*C55)+(N**4+M**4)*C13 CBP(1,2)=M**2*C12+N**2*C23 CBP(1,5)=M*N*((C13+2.*C55-C11)*M**2-(C13+2.*C55-C33)*N**2) CBP(3,3)=N**4*C11+2.*(N*M)**2*(C13+2.*C55)+M**4*C33 CBP(2,3)=N**2*C12+M**2*C23 CBP(2,2)=C22 CBP(2,5)=M*N*(C23-C12) CBP(3,5)=M*N*((C13+2.*C55-C11)*N**2-(C13+2.*C55-C33)*M**2) CBP(4,4)=C44*M**2+C66*N**2 CBP(4,6)=M*N*(C44-C66) CBP(6,6)=C66*M**2+C44*N**2 CBP(5,5)=(M*N)**2*(C11+C33-2.*C13)+C55*(M**2-N**2)**2 CBP(3,1)=CBP(1,3) CBP(5,1)=CBP(1,5) CBP(2,1)=CBP(1,2) CBP(3,2)=CBP(2,3) CBP(5,2)=CBP(2,5) CBP(5,3)=CBP(3,5) CBP(6,4)=CBP(4,6) C C WRITE ROTATED STIFFNESSES FOR TOP LAYER C WRITE(6,251) DO 457 I=1,6 457 WRITE(6,456)(CBP(I,J),J=1,6) WRITE(6,251) C C THE NEXT SECTION SOLVES FOR THE ROOTS OF EQN(10), BY T.C.T. TING(1982). C TING SOLVES EQN(10) WHERE EIGENVALUES EXIST ON THE OFF-DIAGONAL TERMS. C HERE THE PROBLEM IS REFORMULATED INTO A(6,6) WHERE THE THREE PAIRS OF C COMPLEX CONJUGATE ROOTS OF EQN(10) BECOME SIX UNIQUE EIGENVALUES W(6). C THE 1RST, 3RD, AND 5TH ELEMENTS OF W(6) ARE CHOOSEN WITH POSTIVE IMAGINARY C TERMS AS EIGENVALUES P(L=1,2,3). CORRESPONDING EIGENVECTORS V(I,L=1,2,3) C RECOVERED AS I=(4TH, 5TH, AND 6TH ELEMENTS) OF Z(I,J). A FINAL CHECK IS C IS OBTAINED BY SUBSTITUTING THE ABOVE RESULTS INTO EQN(10) OF TING(1982). C DETAILS OF THE ABOVE PROBLEM IS GIVEN IN THE AUTHORS PERSONAL NOTES. C C DEFINE A(I,J) NONZERO TERMS FOR THE TOP LAYER C DEN=CB(6,6)*CB(2,2)*CB(4,4)-CB(2,2)*CB(4,6)**2 A(1,2)=((CB(4,6)+CB(2,5))*CB(2,2)*CB(4,6)-(CB(1,2)+CB(6,6)) * *CB(2,2)*CB(4,4))/DEN A(2,1)=(CB(1,2)+CB(6,6))*(CB(4,6)**2-CB(4,4)*CB(6,6))/DEN A(2,3)=(CB(4,6)+CB(2,5))*(CB(4,6)**2-CB(4,4)*CB(6,6))/DEN A(3,2)=((CB(1,2)+CB(6,6))*CB(4,6)*CB(2,2)-(CB(4,6)+CB(2,5)) * *CB(2,2)*CB(6,6))/DEN A(1,4)=(CB(2,2)*CB(4,6)*CB(1,5)-CB(2,2)*CB(4,4)*CB(1,1))/DEN A(1,6)=(CB(2,2)*CB(4,6)*CB(5,5)-CB(2,2)*CB(4,4)*CB(1,5))/DEN A(2,5)=(CB(4,6)**2-CB(4,4)*CB(6,6))*CB(6,6)/DEN A(3,4)=(CB(4,6)*CB(2,2)*CB(1,1)-CB(2,2)*CB(6,6)*CB(1,5))/DEN A(3,6)=(CB(4,6)*CB(2,2)*CB(1,5)-CB(2,2)*CB(6,6)*CB(5,5))/DEN A(4,1)=1.0 A(5,2)=1.0 A(6,3)=1.0 C C DEFINE AP(I,J) NONZERO TERMS FOR THE BOTTOM LAYER C DEN=CBP(6,6)*CBP(2,2)*CBP(4,4)-CBP(4,6)*CBP(3,5)**2 AP(1,2)=((CBP(4,6)+CBP(2,5))*CBP(2,2)*CBP(4,6)-(CBP(1,2)+CBP(6,6)) * *CBP(2,2)*CBP(4,4))/DEN AP(2,1)=(CBP(1,2)+CBP(6,6))*(CBP(4,6)**2-CBP(4,4)*CBP(6,6))/DEN AP(2,3)=(CBP(4,6)+CBP(2,5))*(CBP(4,6)**2-CBP(4,4)*CBP(6,6))/DEN AP(3,2)=((CBP(1,2)+CBP(6,6))*CBP(4,6)*CBP(2,2)-(CBP(4,6)+CBP(2,5)) * *CBP(2,2)*CBP(6,6))/DEN AP(1,4)=(CBP(2,2)*CBP(4,6)*CBP(1,5)-CBP(2,2)*CBP(4,4)*CBP(1,1)) * /DEN AP(1,6)=(CBP(2,2)*CBP(4,6)*CBP(5,5)-CBP(2,2)*CBP(4,4)*CBP(1,5)) * /DEN AP(2,5)=(CBP(4,6)**2-CBP(4,4)*CBP(6,6))*CBP(6,6)/DEN AP(3,4)=(CBP(4,6)*CBP(2,2)*CBP(1,1)-CBP(2,2)*CBP(6,6)*CBP(1,5)) * /DEN AP(3,6)=(CBP(4,6)*CBP(2,2)*CBP(1,5)-CBP(2,2)*CBP(6,6)*CBP(5,5)) * /DEN AP(4,1)=1.0 AP(5,2)=1.0 AP(6,3)=1.0 C C WRITE A(I,J) AND AP(I,J) TO OUTPUT C WRITE(6,251) WRITE(6,281) 281 FORMAT(' MATRICES FOR EIGEN-VALUE PROBLEM FOR TOP/BOTTOM' * ,/,' LAYERS ARE WRITTEN BELOW [A]/[AP]') WRITE(6,251) DO 200 I=1,6 WRITE(6,250)(A(I,J),J=1,6) 200 CONTINUE WRITE(6,251) DO 201 I=1,6 WRITE(6,250)(AP(I,J),J=1,6) 201 CONTINUE WRITE(6,251) 250 FORMAT(6(1X,E11.4)) 251 FORMAT(1X,80(1H*)) C C SOLVE FOR EIGEN VALUES AND EIGEN VECTORS FOR BOTH A(6,6) AND AP(6,6) C c CALL EIGRF(A,6,6,2,W,Z,6,WK,IER) c CALL EIGRF(AP,6,6,2,WP,ZP,6,WKP,IERP) c call evcrg(6,a,6,w,z,6) pindex=epirg(6,6,a,6,w,z,6) c call evcrg(6,ap,6,wp,zp,6) pindexp=epirg(6,6,ap,6,wp,zp,6) C C WRITE PERFORMANCE PARAMETERS AND ARRAYS W,WP,Z, AND ZP C c WRITE(6,4)WK(1),IER c 4 FORMAT(' TOP LAYER, CONVERGED=',E10.3,' IER=',I6) WRITE(6,4)pindex 4 FORMAT(' TOP LAYER, CONVERGED=',E10.3) DO 41 J=1,6 WRITE(6,42)J,W(J),(I,Z(I,J),I=1,3) 41 WRITE(6,43)(I,Z(I,J),I=4,6) 42 FORMAT(' W(',I1,')',2E10.3,/,3(' Z(',I1,')',2E10.3)) 43 FORMAT(3(' Z(',I1,')',2E10.3)) WRITE(6,251) c WRITE(6,44)WKP(1),IERP c 44 FORMAT(' BOTTOM LAYER, CONVERGED=',E10.3,' IERP=',I6) WRITE(6,44)pindexp 44 FORMAT(' BOTTOM LAYER, CONVERGED=',E10.3) DO 441 J=1,6 WRITE(6,442)J,WP(J),(I,ZP(I,J),I=1,3) 441 WRITE(6,443)(I,ZP(I,J),I=4,6) 442 FORMAT(' WP(',I1,')',2E10.3,/,3(' ZP(',I1,')',2E10.3)) 443 FORMAT(3(' ZP(',I1,')',2E10.3)) C C THIS LOOP REASSIGNS INDICES FOR EIGENVALUES AND EIGENVECTORS, NORMALIZES C THE EIGENVECTORS, AND CHECKS RESUTLTS WITH TING'S ORIGINAL EQUATION. C JJ=0 C C CHOOSE COMPLEX ROOTS OR THEIR COMPLEX CONJUGATES C C LOOP THRU ARRAYS AND CHOOSE COMPLEX ROOTS(EIGENVALUES) WITH POSITIVE C IMAGINARY TERMS AND THEIR CORRESPONDING EIGENVECTORS DO 65 J=1,5,2 C C LOOP THRU ARRAYS AND CHOOSE COMPLEX CONJUGATE ROOTS(EIGENVALUES) WITH C NEGATIVE IMAGINARY TERMS AND THEIR CORRESPONDING EIGENVECTORS C DO 65 J=2,6,2 C C INCREMENT THE INDICE FOR THE ROOTS(EIGENVALUES) JJ=JJ+1 C C REASSIGN INDICES FROM Z(6,6),ZP(6,6) TO V(3,3),VP(3,3) C AND NORMALIZE V(3,3),VP(3,3) FOR TOP AND BOTTOM(P) LAYERS C SUM=(0.0,0.0) SUMP=(0.0,0.0) DO 50 I=4,6 SUM=SUM+Z(I,J)*CONJG(Z(I,J)) SUMP=SUMP+ZP(I,J)*CONJG(ZP(I,J)) 50 CONTINUE DO 51 I=4,6 V(I-3,JJ)=CMI*Z(I,J)/CSQRT(SUM) VP(I-3,JJ)=CMI*ZP(I,J)/CSQRT(SUMP) 51 CONTINUE C C REASSIGN INDICES FROM W(6),WP(6) TO P(3),PP(3) C FOR CRACKED AND UNCRACKED(P) LAYERS C P(JJ)=W(J) PP(JJ)=WP(J) C C CHECK BY SUBSTITUTING RESULTS INTO EQN(10) OF T.C.T TING(1982). C B(1,1)=CMPLX(CB(1,1))+P(JJ)*P(JJ)*CMPLX(CB(6,6)) B(1,2)=(CMPLX(CB(1,2))+CMPLX(CB(6,6)))*P(JJ) B(1,3)=CMPLX(CB(1,5))+P(JJ)*P(JJ)*CMPLX(CB(4,6)) B(2,2)=CMPLX(CB(6,6))+P(JJ)*P(JJ)*CMPLX(CB(2,2)) B(2,3)=(CMPLX(CB(4,6))+CMPLX(CB(2,5)))*P(JJ) B(3,3)=CMPLX(CB(5,5))+P(JJ)*P(JJ)*CMPLX(CB(4,4)) B(3,1)=B(1,3) B(3,2)=B(2,3) B(2,1)=B(1,2) BP(1,1)=CMPLX(CBP(1,1))+PP(JJ)*PP(JJ)*CMPLX(CBP(6,6)) BP(1,2)=(CMPLX(CBP(1,2))+CMPLX(CBP(6,6)))*PP(JJ) BP(1,3)=CMPLX(CBP(1,5))+PP(JJ)*PP(JJ)*CMPLX(CBP(4,6)) BP(2,2)=CMPLX(CBP(6,6))+PP(JJ)*PP(JJ)*CMPLX(CBP(2,2)) BP(2,3)=(CMPLX(CBP(4,6))+CMPLX(CBP(2,5)))*PP(JJ) BP(3,3)=CMPLX(CBP(5,5))+PP(JJ)*PP(JJ)*CMPLX(CBP(4,4)) BP(3,1)=BP(1,3) BP(3,2)=BP(2,3) BP(2,1)=BP(1,2) DO 62 III=1,3 SUM=(0.0,0.0) SUMP=(0.0,0.0) DO 60 JJJ=1,3 SUM=SUM+B(III,JJJ)*V(JJJ,JJ) SUMP=SUMP+BP(III,JJJ)*VP(JJJ,JJ) 60 CONTINUE WRITE(6,251) WRITE(6,61)III,SUM,III,SUMP 61 FORMAT(' SUM FOR EQUATION NO.',I1,' =',2E11.4,' TOP LAYER' * ,/,' SUM FOR EQUATION NO.',I1,' =',2E11.4,' BOTTOM LAYER') 62 CONTINUE WRITE(6,251) WRITE(6,63)JJ,P(JJ),(I,JJ,V(I,JJ),I=1,3) 63 FORMAT(' EIGEN-VALUE P(',I1,')=',2E11.4,' TOP LAYER',/, * ' NORMALIZED EIGEN-VECTOR V(',I1,',',I1,')=',2E11.4,/, * 24X,' V(',I1,',',I1,')=',2E11.4,/, * 24X,' V(',I1,',',I1,')=',2E11.4) WRITE(6,251) WRITE(6,64)JJ,PP(JJ),(I,JJ,VP(I,JJ),I=1,3) 64 FORMAT(' EIGEN-VALUE PP(',I1,')=',2E11.4,' BOTTOM LAYER',/, * ' NORMALIZED EIGEN-VECTOR VP(',I1,',',I1,')=',2E11.4,/, * 24X,' VP(',I1,',',I1,')=',2E11.4,/, * 24X,' VP(',I1,',',I1,')=',2E11.4) WRITE(6,251) 65 CONTINUE C C CALCULATE T(I,J,L) AND TP(I,J,L) MATRICES FOR TOP AND BOTTOM LAYERS C BY USING FUNCTION SUBROUTINES FT AND FTP. THESE MATRICES ARE DEFINED BY C EQN(9) OF T.C.T. TING(1982). C DO 73 I=1,3 DO 73 J=1,3 DO 73 L=1,3 T(I,J,L)=FT(I,J,L) TP(I,J,L)=FTP(I,J,L) 73 CONTINUE WRITE(6,251) C C WRITE T(I,J,L) AND TP(I,J,L) ARRAYS C DO 74 I=1,3 DO 76 J=1,3 76 WRITE(6,77)(I,J,L,T(I,J,L),L=1,3) 77 FORMAT(3(' T',I1,I1,I1,2E10.3)) 74 CONTINUE WRITE(6,251) DO 75 I=1,3 DO 78 J=1,3 78 WRITE(6,79)(I,J,L,TP(I,J,L),L=1,3) 79 FORMAT(3(' TP',I1,I1,I1,2E10.3)) 75 CONTINUE WRITE(6,251) C C TEST P(L) AND PP(L) ARRAYS WITH T(I,J,L) AND TP(I,J,L) ARRAYS C BY SUBSTITUTING ABOVE RESULTS INTO EQN(10) OF T.C.T TING(1982). C DO 80 L=1,3 TEST1=T(1,2,L)+P(L)*T(2,2,L) TEST2=T(1,3,L)+P(L)*T(2,3,L) TEST3=T(1,1,L)-P(L)*P(L)*T(2,2,L) TEST1P=TP(1,2,L)+PP(L)*TP(2,2,L) TEST2P=TP(1,3,L)+PP(L)*TP(2,3,L) TEST3P=TP(1,1,L)-PP(L)*PP(L)*TP(2,2,L) WRITE(6,81)L,TEST1,TEST2,TEST3 81 FORMAT(' L=',I1,' T1 =',2E10.3,' T2 =',2E10.3,' T3 =',2E10.3) WRITE(6,82)L,TEST1P,TEST2P,TEST3P 82 FORMAT(' L=',I1,' T1P=',2E10.3,' T2P=',2E10.3,' T3P=',2E10.3) WRITE(6,251) 80 CONTINUE WRITE(6,251) C C THIS NEXT SECTION DETERMINES THE THREE ROOTS OF AN ARRAY C(12,12). C THE ARRAY C(12,12) IS DEFINED BY T.C.T. TING(1982) BY PRESCRIBING C NECESSARY BOUNDARY CONDITIONS AT A STRESS FREE EDGE FOR A C BIMATERIAL INTERFACE WHERE THE DISSIMILAR MATERIALS ARE OTHOTROPIC. C C THIS LOOP INCREMENTS THE EXPONENT K FROM 0.0 TO 1.0 AND PROVIDES C AN APPROXIMATION FOR THE ROOT OF THE DETERMINATE OF C(12,12). C KR(1)=0.0 IR=1 K=0.00001 DO 750 IK=1,100 C C CALL THE FUNCTION "F" AND CALCULATE THE DETERMINATE AS A FUNCTION OF "K" DET=F(K) KT(IK)=DET IF(IK.EQ.1)GO TO 652 C C TEST FOR A SIGN CHANGE OF THE DETERMINANT TEST=KT(IK)*KT(IK-1) IF(TEST.GT.0.0)GO TO 652 C C TEST FOR MORE THAN ONE ROOT IN THE K INTERVAL 0.0 - 1.0 C IF MORE THAN ROOT EXISTS TERMINATE PROGRAM C aaa=kt(ik-1) bbb=kt(ik) IF(IR.LT.2)GO TO 655 WRITE(6,656)(I,KR(I),I=1,2) 656 FORMAT(' MORE THAN ONE ROOT IS FOUND BETWEEN 0.0 AND 1.0',/, *' THE ROOTS ARE',2(1X,'KR(',I1,')=',E12.5),/, *' PROGRAM TERMINATED') STOP C 655 KR(IR)=K-0.01 IR=IR+1 C C WRITE THE LOOP INCREMENT, EXPONENT, AND DETERMINANT C 652 WRITE(6,653)IK,K,DET 653 FORMAT(' IK=',I4,' K=',F12.9,' DET=',E16.9) C C INCREMENT ON THE EXPONENT AND END LOOP C K=K+0.01 750 CONTINUE C C TEST IF ROOT EXISTS OTHERWISE TERMINATE PROGRAM C IF(IR.EQ.1)WRITE(6,661) 661 FORMAT(' NO ROOT EXISTS FOR K BETWEEN 0.0 AND 1.0',/, *' PROGRAM TERMINATED') IF(IR.EQ.1)STOP C C WRITE ROOTS TO OUTPUT C WRITE(6,660)KR(1) 660 FORMAT(' APPROXIMATION FOR ROOT, K BETWEEN 0.0 AND 1.0' *,/,1X,' KR(1)=',E12.5) aaa=kr(1) bbb=kr(1)+0.01 write(6,666)aaa,bbb 666 format(' aaa=',e12.5,' bbb=',e12.5) C C IMSL SUBROUTINE CALCULATES IMPROVED ROOT KR(1) FROM THE INITIAL C APPROX.S FOR bbb GIVEN ABOVE. WE ASSUME THAT ONLY ONE ROOT EXIST. C c EPS=1.0E-5 c EPS2=1.0E-5 c ETA=1.0E-5 c NSIG=5 c ITMAX=100 c errabs=1.0e-5 errrel=1.0e-05 maxfn=100 C c CALL ZREAL1(F,EPS,EPS2,ETA,NSIG,1,KR,ITMAX,IER) c call zbren(f,errabs,errrel,aaa,bbb,maxfn) c c WRITE(6,989)KR(1) c 989 FORMAT(' EXACT ROOT TO WITHIN 5 SIGNIFICANT FIGURES',/, c *1X,' KR(1)=',E12.5) WRITE(6,989)bbb 989 FORMAT(' EXACT ROOT TO WITHIN 5 SIGNIFICANT FIGURES',/, *1X,' bbb=',E12.5) write(7,898)bbb 898 format(e12.5) C 888 STOP END C****************************************************************************** FUNCTION F(K) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * THIS FUNCTION CALCULATES THE NONZERO ELEMENTS OF THE C(12,12) ARRAY, * C * CALCULATES THE DETERMINANT OF THE ARRAY, AND RETURNS THE DETERMINANT * C * AS FUNTION K FOR A GIVEN VALUE OF THE EXPONENT K. DETAILS ARE GIVEN BY * C * T.C.T. TING IN 1982 PUBLICATION. A function subroutine is prescribed if * C * the IMSL routine used to find more accurate solution requires that the * C * determinate of the C matrix changes sign in the region: * C * * C * aaa=k(ik)=(+) -> bbb=k(ik)=(-) * C * where the initial quess is aaa * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON /L1/CB(6,6),P(3),V(3,3) COMMON /L2/CBP(6,6),PP(3),VP(3,3) COMMON /L3/T(3,3,3),TP(3,3,3) C COMPLEX P,PP,V,VP,T,TP,CN REAL F,K,CB,CBP,C(12,12),CC(12,12),DUM(12),D1,D2 real fac(12,12),ipvt(12) INTEGER I,J,L C C CALCULATE THE NONZERO COMPONENTS OF THE C(12,12) MATRIX C AS A FUNCTION OF THE EXPONENT K C C CLEAR THE C(I,J) MATRIX C DO 756 I=1,12 DUM(I)=0.0 DO 756 J=1,12 C(I,J)=0.0 756 CONTINUE C C CALCULATE THE NONZERO TERMS OF THE C(I,J) MATRIX C C BLOCKS 1 AND 2 DO 11 J=1,3 DO 11 L=1,3 CN=T(1,J,L)*(P(L))**(-K) C(J,L)=REAL(CN) C(J,L+3)=AIMAG(CN) 11 CONTINUE C C BLOCKS 3 AND 4 DO 13 J=1,3 DO 13 L=1,3 CN=TP(1,J,L)*(-PP(L))**(-K) C(J+3,L+6)=REAL(CN) C(J+3,L+9)=AIMAG(CN) 13 CONTINUE C C BLOCKS 5 AND 6 DO 15 J=1,3 DO 15 L=1,3 CN=T(2,J,L) C(J+6,L)=REAL(CN) C(J+6,L+3)=AIMAG(CN) 15 CONTINUE C C BLOCKS 7 AND 8 DO 17 J=1,3 DO 17 L=1,3 CN=TP(2,J,L) C(J+6,L+6)=-REAL(CN) C(J+6,L+9)=-AIMAG(CN) 17 CONTINUE C C BLOCKS 9 AND 10 DO 19 J=1,3 DO 19 L=1,3 CN=V(J,L) C(J+9,L)=REAL(CN) C(J+9,L+3)=AIMAG(CN) 19 CONTINUE C C BLOCKS 11 AND 12 DO 21 J=1,3 DO 21 L=1,3 CN=VP(J,L) C(J+9,L+6)=-REAL(CN) C(J+9,L+9)=-AIMAG(CN) 21 CONTINUE C C STORE ARRAY FOR LATER REFERENCE IF ERROR OCCURS IN SUBROUTINE LINV3F C DO 899 I=1,12 DO 899 J=1,12 CC(I,J)=C(I,J)/1.0E+03 C(I,J)=CC(I,J) 899 CONTINUE C C CALCULATE DETERMINANT OF THE C(12,12) MATRIX C c D1=1.0 c CALL LINV3F(C,DUM,4,12,12,D1,D2,WKS,IER) c IF(IER.EQ.0)GO TO 654 c WRITE(6,651)IER c 651 FORMAT('ERROR IN CALCULATING DETERMINATE, IER=',I4, c *' PROGRAM TERMINATED, ARRAY C(12,12) WRITTEN BELOW') c call lftrg(12,c,12,fac,12,ipvt) call lfdrg(12,fac,12,ipvt,d1,d2) c c GO TO 889 C C THE FUNCTION F IS RETURNED AS THE DETERMINATE OF ARRAY C(12,12) C c 654 F=D1*2.0**D2 c c F=D1*10.0**D2 F=D1 c RETURN C C OTHERWISE IF IER IS GREATER THAN ZERO; THE CC(12,12) ARRAY IS WRITTEN C TO OUTPUT AND THE PROGRAN IS TERMINATED. C c 889 DO 667 I=1,12 c 667 WRITE(6,668)(CC(I,J),J=1,6) c 668 FORMAT(6E10.3) c WRITE(6,251) c 251 FORMAT(1X,80(1H*)) c DO 669 I=1,12 c 669 WRITE(6,668)(CC(I,J),J=7,12) c STOP C END C****************************************************************************** FUNCTION FT(I,J,L) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * ELEMENTS OF THE T(I,J,L) ARRAY ARE CALCULATED FOR THE TOP LAYER * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON /L1/CB(6,6),P(3),V(3,3) COMPLEX P,V,T,FT REAL CB INTEGER IC,JC1,JC2,I,J,K,L IF((I.EQ.1).AND.(J.EQ.1))IC=1 IF((I.EQ.1).AND.(J.EQ.2))IC=6 IF((I.EQ.1).AND.(J.EQ.3))IC=5 IF((I.EQ.2).AND.(J.EQ.1))IC=6 IF((I.EQ.2).AND.(J.EQ.2))IC=2 IF((I.EQ.2).AND.(J.EQ.3))IC=4 IF((I.EQ.3).AND.(J.EQ.1))IC=5 IF((I.EQ.3).AND.(J.EQ.2))IC=4 IF((I.EQ.3).AND.(J.EQ.3))IC=3 T=(0.0,0.0) DO 10 K=1,3 IF(K.EQ.1)JC1=1 IF(K.EQ.1)JC2=6 IF(K.EQ.2)JC1=6 IF(K.EQ.2)JC2=2 IF(K.EQ.3)JC1=5 IF(K.EQ.3)JC2=4 T=T+(CMPLX(CB(IC,JC1))+P(L)*CMPLX(CB(IC,JC2)))*V(K,L) 10 CONTINUE FT=T RETURN END C****************************************************************************** FUNCTION FTP(I,J,L) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * ELEMENTS OF THE TP(I,J,L) ARRAY ARE CALCULATED FOR THE BOTTOM LAYER * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON /L2/CBP(6,6),PP(3),VP(3,3) COMPLEX PP,VP,TP,FTP REAL CBP INTEGER IC,JC1,JC2,I,J,K,L IF((I.EQ.1).AND.(J.EQ.1))IC=1 IF((I.EQ.1).AND.(J.EQ.2))IC=6 IF((I.EQ.1).AND.(J.EQ.3))IC=5 IF((I.EQ.2).AND.(J.EQ.1))IC=6 IF((I.EQ.2).AND.(J.EQ.2))IC=2 IF((I.EQ.2).AND.(J.EQ.3))IC=4 IF((I.EQ.3).AND.(J.EQ.1))IC=5 IF((I.EQ.3).AND.(J.EQ.2))IC=4 IF((I.EQ.3).AND.(J.EQ.3))IC=3 TP=(0.0,0.0) DO 10 K=1,3 IF(K.EQ.1)JC1=1 IF(K.EQ.1)JC2=6 IF(K.EQ.2)JC1=6 IF(K.EQ.2)JC2=2 IF(K.EQ.3)JC1=5 IF(K.EQ.3)JC2=4 TP=TP+(CMPLX(CBP(IC,JC1))+PP(L)*CMPLX(CBP(IC,JC2)))*VP(K,L) 10 CONTINUE FTP=TP RETURN END