PROGRAM plysing C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C* PROGRAM : CALCULATES THE SINGULARITIES OF THE DISPLACEMENT FIELD * C* FOR A CRACK TIP PERPENDICULAR TO A BIMATERIAL INTERFACE WHERE THE * C* DISSIMILAR MATERIALS HAVE ORTHORHOMBIC ELASTIC SYMMETRY. EACH * C* LAYER IS A FIBER-REINFORCED MATERIAL WHOSE FIBER ORIENTATION IS * C* MEASURED BY THE ANGLE THETA FROM THE LOAD AXIS. IN THIS PROGRAM * C* ALL VARIABLES RELATED TO THE LAYER ADJACENT TO THE CRACKED LAYER * C* ARE SUFFIXED BY THE LETTER P. * C* * C* EXAMPLE: STIFFNESSES OF LAYER WITH CRACK : CB(I,J) * C* STIFFNESS OF ADJACENT LAYER WITH NO CRACK: CBP(I,J) * C* * C* THIS SOLUTION IS OUTLINED BY T.C.T.TING AND R.H.HOANG,INT.J.SOLID * C* AND STRUCTURES,VOL.20,NO.5,PP 439-454,1984. IN THIS PROGRAM * C* CORRECTED ELASTIC PROPERTIES ARE COMPARED WITH PREVIOUS APPROX. * C* * C* LATEST UPDATE: **** 2/16/87 **** * 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* 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(3) 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='plysing.dat',STATUS='OLD',ERR=888) OPEN(6,FILE='plysing.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 NUMBER MINUS "I" 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 Poisson's ratio in 12-plane, NU21=+0.211E+00 c Poisson's ratio in 13-plane, NU31=+0.211E+00 c Poisson'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.) (CRACKED LAYER: THETA) (UNCRACKED 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=+90.0 c THETAP=+00.0 C C WRITE OUT FIBER ANGLE FOR LAYERS WITH AND WITHOUT CRACK C WRITE(6,251) WRITE(6,14)THETA,THETAP 14 FORMAT(' LAYER WITH CRACK, ANGLE (DEG.)=',E10.3,/, * ' LAYER WITHOUT CRACK, ANGLE (DEG.)=',E10.3,//, *' BELOW ARE STIFFNESSES FOR CRACKED/UNCRACKED LAYERS') WRITE(6,251) C C CALCULATE STIFFNESSES IN THE ROTATED X3-X1 PLANE FOR THE CRACKED 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 CRACKED 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 UNCRACKED 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 UNCRACKED 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(11), BY T.C.T. TING(1984). C TING SOLVES EQN(11) 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(11) 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(8) OF TING(1984). C DETAILS OF THE ABOVE PROBLEM IS GIVEN IN THE AUTHORS PERSONAL NOTES. C C DEFINE A(I,J) NONZERO TERMS FOR THE CRACKED LAYER C DEN=CB(3,3)*CB(4,4)*CB(5,5)-CB(4,4)*CB(3,5)**2 A(1,2)=((CB(2,3)+CB(4,4))*CB(4,4)*CB(3,5)-(CB(4,6)+CB(2,5)) * *CB(4,4)*CB(3,3))/DEN A(2,1)=(CB(4,6)+CB(2,5))*(CB(3,5)**2-CB(3,3)*CB(5,5))/DEN A(2,3)=(CB(2,3)+CB(4,4))*(CB(3,5)**2-CB(3,3)*CB(5,5))/DEN A(3,2)=((CB(4,6)+CB(2,5))*CB(3,5)*CB(4,4)-(CB(2,3)+CB(4,4)) * *CB(4,4)*CB(5,5))/DEN A(1,4)=(CB(4,4)*CB(3,5)*CB(4,6)-CB(4,4)*CB(3,3)*CB(6,6))/DEN A(1,6)=(CB(3,5)*CB(4,4)**2-CB(4,4)*CB(3,3)*CB(4,6))/DEN A(2,5)=(CB(3,5)**2-CB(3,3)*CB(5,5))*CB(2,2)/DEN A(3,4)=(CB(3,5)*CB(4,4)*CB(6,6)-CB(4,4)*CB(5,5)*CB(4,6))/DEN A(3,6)=(CB(3,5)*CB(4,4)*CB(4,6)-CB(5,5)*CB(4,4)**2)/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 UNCRACKED LAYER C DEN=CBP(3,3)*CBP(4,4)*CBP(5,5)-CBP(4,4)*CBP(3,5)**2 AP(1,2)=((CBP(2,3)+CBP(4,4))*CBP(4,4)*CBP(3,5)-(CBP(4,6)+CBP(2,5)) * *CBP(4,4)*CBP(3,3))/DEN AP(2,1)=(CBP(4,6)+CBP(2,5))*(CBP(3,5)**2-CBP(3,3)*CBP(5,5))/DEN AP(2,3)=(CBP(2,3)+CBP(4,4))*(CBP(3,5)**2-CBP(3,3)*CBP(5,5))/DEN AP(3,2)=((CBP(4,6)+CBP(2,5))*CBP(3,5)*CBP(4,4)-(CBP(2,3)+CBP(4,4)) * *CBP(4,4)*CBP(5,5))/DEN AP(1,4)=(CBP(4,4)*CBP(3,5)*CBP(4,6)-CBP(4,4)*CBP(3,3)*CBP(6,6)) * /DEN AP(1,6)=(CBP(3,5)*CBP(4,4)**2-CBP(4,4)*CBP(3,3)*CBP(4,6))/DEN AP(2,5)=(CBP(3,5)**2-CBP(3,3)*CBP(5,5))*CBP(2,2)/DEN AP(3,4)=(CBP(3,5)*CBP(4,4)*CBP(6,6)-CBP(4,4)*CBP(5,5)*CBP(4,6)) * /DEN AP(3,6)=(CBP(3,5)*CBP(4,4)*CBP(4,6)-CBP(5,5)*CBP(4,4)**2)/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 CRACKED/UNCRACKED' * ,/,' 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 IER AND IERP AND ARRAYS W,WP,Z, AND ZP C c WRITE(6,4)WK(1),IER c 4 FORMAT(' LAYER WITH CRACK, CONVERGED=',E10.3,' IER=',I6) 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(' LAYER WITHOUT CRACK CONVERGED=',E10.3,' IERP=',I6) 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 CRACKED AND UNCRACKED(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(8) OF T.C.T TING(1984). C B(1,1)=CMPLX(CB(6,6))+P(JJ)*P(JJ)*CMPLX(CB(5,5)) B(1,2)=(CMPLX(CB(4,6))+CMPLX(CB(2,5)))*P(JJ) B(1,3)=CMPLX(CB(4,6))+P(JJ)*P(JJ)*CMPLX(CB(3,5)) B(2,2)=CMPLX(CB(2,2))+P(JJ)*P(JJ)*CMPLX(CB(4,4)) B(2,3)=(CMPLX(CB(2,3))+CMPLX(CB(4,4)))*P(JJ) B(3,3)=CMPLX(CB(4,4))+P(JJ)*P(JJ)*CMPLX(CB(3,3)) B(3,1)=B(1,3) B(3,2)=B(2,3) B(2,1)=B(1,2) BP(1,1)=CMPLX(CBP(6,6))+PP(JJ)*PP(JJ)*CMPLX(CBP(5,5)) BP(1,2)=(CMPLX(CBP(4,6))+CMPLX(CBP(2,5)))*PP(JJ) BP(1,3)=CMPLX(CBP(4,6))+PP(JJ)*PP(JJ)*CMPLX(CBP(3,5)) BP(2,2)=CMPLX(CBP(2,2))+PP(JJ)*PP(JJ)*CMPLX(CBP(4,4)) BP(2,3)=(CMPLX(CBP(2,3))+CMPLX(CBP(4,4)))*PP(JJ) BP(3,3)=CMPLX(CBP(4,4))+PP(JJ)*PP(JJ)*CMPLX(CBP(3,3)) 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,' CRACKED LAYER' * ,/,' SUM FOR EQUATION NO.',I1,' =',2E11.4,' UNCRACKED 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,' CRACKED 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,' LAYER NO CRACK',/, * ' 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 CRACKED AND UNCRACKED LAYERS C BY USING FUNCTION SUBROUTINES FT AND FTP. THESE MATRICES ARE DEFINED BY C EQN(9) OF T.C.T. TING(1984). 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 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(12) OF T.C.T TING(1984). C DO 80 L=1,3 TEST1=T(1,2,L)+P(L)*T(1,3,L) TEST2=T(2,2,L)+P(L)*T(2,3,L) TEST3=T(2,2,L)-P(L)*P(L)*T(3,3,L) TEST1P=TP(1,2,L)+PP(L)*TP(1,3,L) TEST2P=TP(2,2,L)+PP(L)*TP(2,3,L) TEST3P=TP(2,2,L)-PP(L)*PP(L)*TP(3,3,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(18,18). C THE ARRAY C(18,18) IS DEFINED IN SECTION 4 OF T.C.T. TING(1984) BY C PRESCRIBING NECESSARY BOUNDARY CONDITIONS FOR A CRACK NORMAL TO 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 ROOTS OF THE DETERMINATE OF C(18,18). C KR(1)=0.0 KR(2)=0.0 KR(3)=0.0 IR=1 K=0.0000001 DO 750 IK=1,100 C 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 THREE ROOTS IN THE INTERVAL K(0.0,1.0) c if more than three roots exist terminate program c aaa=kt(ik-1) bbb=kt(ik) IF(IR.LT.4)GO TO 655 WRITE(6,656)(I,KR(I),I=1,3) 656 FORMAT(' MORE THAN THREE ROOTS EXIST FOR K BETWEEN 0.0 AND 1.0' *,/,' BELOW WE WRITE ONLY THE FIRST THREE ROOTS' *,/,1X,3(' KR(',I1,')=',F12.9)) 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 FOR LESS THAN THREE ROOTS AND WRITE ROOTS TO OUTPUT C IF(IR.LT.3)WRITE(6,661)IR 661 FORMAT(' BETWEEN 0.0 AND 0.1 ONLY ',I1,' ROOT/S WERE FOUND') WRITE(6,660)(I,KR(I),I=1,3) 660 FORMAT(' APPROXIMATION FOR THREE ROOTS FOR K BETWEEN 0.0 AND 1.0' *,/,1X,3(' KR(',I1,')=',F12.9)) IF(IR.LT.3)STOP C C WRITE ROOTS TO OUTPUT C WRITE(6,662)KR(1) 662 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 ZREAL1 CALCULATES IMPROVED ROOTS KR(3) FROM THE INITIAL C APPROX.S FOR KR(3) GIVEN ABOVE. WE ASSUME THAT ONLY THREE ROOTS EXIST. C c EPS=1.0E-5 c EPS2=1.0E-5 c ETA=1.0E-2 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,3,KR,ITMAX,IER) c call zbren(f,errabs,errrel,aaa,bbb,maxfn) c c WRITE(6,989)(I,KR(I),I=1,3) c 989 FORMAT(' EXACT ROOTS TO WITHIN 5 SIGNIFICANT FIGURES',/, c *1X,3(' KR(',I1,')=',F12.9)) 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(18,18) ARRAY, * C * CALCULATES THE DETERMINANT OF THE ARRAY, AND RETURNS THE DETERMINANT * C * AS FUNTION OF THE EXPONENT K. DETAILS ARE GIVEN IN SECTION 4 OF 84 PUB. * C * A function subroutine is prescribed if the IMSL routine used to find * C * more accurate solution requires that the determinate of the C matrix * C * changes sign in the region: * C * * C * aaa=k(ik)=(+) -> bbb=k(ik)=(-) * C * where the initial quess is aaa * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 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 COMPLEX P,PP,V,VP,T,TP,CM1,CN REAL F,K,CB,CBP,C(18,18),CC(18,18),DUM(18),D1,D2 real fac(12,12),ipvt(12) INTEGER I,J,L C C DEFINE COMPLEX MINUS ONE C CM1=(-1.000000000000000,0.000000000000000) C C CALCULATE THE NONZERO COMPONENTS OF THE C(18,18) MATRIX C AS A FUNCTION OF THE EXPONENT K C C CLEAR THE C(I,J) MATRIX C DO 756 I=1,18 DUM(I)=0.0 DO 756 J=1,18 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(3,J,L)*CM1**(-K) C(J,L)=+REAL(CN) C(J,L+3)=+AIMAG(CN) 11 CONTINUE C BLOCKS 5 AND 6 DO 15 J=1,3 DO 15 L=1,3 CN=T(2,J,L)*(P(L))**(-K) C(J+3,L)=+REAL(CN) C(J+3,L+3)=+AIMAG(CN) 15 CONTINUE C BLOCKS 7 AND 8 DO 17 J=1,3 DO 17 L=1,3 CN=TP(2,J,L)*(PP(L))**(-K) C(J+3,L+6)=-REAL(CN) C(J+3,L+9)=-AIMAG(CN) 17 CONTINUE C BLOCKS 9 AND 10 DO 19 J=1,3 DO 19 L=1,3 CN=V(J,L)*(P(L))**(1.-K) C(J+6,L)=+REAL(CN) C(J+6,L+3)=+AIMAG(CN) 19 CONTINUE C BLOCKS 11 AND 12 DO 21 J=1,3 DO 21 L=1,3 CN=VP(J,L)*(PP(L))**(1.-K) C(J+6,L+6)=-REAL(CN) C(J+6,L+9)=-AIMAG(CN) 21 CONTINUE C BLOCKS 13 AND 14 DO 23 J=1,3 DO 23 L=1,3 CN=TP(2,J,L)*(-PP(L))**(-K) C(J+9,L+6)=-REAL(CN) C(J+9,L+9)=-AIMAG(CN) 23 CONTINUE C BLOCKS 15 AND 16 DO 25 J=1,3 DO 25 L=1,3 CN=T(2,J,L)*(-P(L))**(-K) C(J+9,L+12)=+REAL(CN) C(J+9,L+15)=+AIMAG(CN) 25 CONTINUE C BLOCKS 17 AND 18 DO 27 J=1,3 DO 27 L=1,3 CN=VP(J,L)*(-PP(L))**(1.-K) C(J+12,L+6)=+REAL(CN) C(J+12,L+9)=+AIMAG(CN) 27 CONTINUE C BLOCKS 19 AND 20 DO 29 J=1,3 DO 29 L=1,3 CN=V(J,L)*(-P(L))**(1.-K) C(J+12,L+12)=-REAL(CN) C(J+12,L+15)=-AIMAG(CN) 29 CONTINUE C BLOCKS 3 AND 4 DO 13 J=1,3 DO 13 L=1,3 CN=T(3,J,L)*CM1**(-K) C(J+15,L+12)=+REAL(CN) C(J+15,L+15)=+AIMAG(CN) 13 CONTINUE C C STORE ARRAY FOR LATER REFERENCE IF ERROR OCCURS IN SUBROUTINE LINV3F C DO 899 I=1,18 DO 899 J=1,18 CC(I,J)=C(I,J)/1.0E+03 C(I,J)=CC(I,J) 899 CONTINUE C C CALCULATE DETERMINANT OF THE C(18,18) MATRIX C BY C IMSL SUBROUTINE LINV3F c D1=1.0 c CALL LINV3F(C,DUM,4,18,18,D1,D2,WKS,IER) c IF(IER.NE.0)GO TO 654 c call lftrg(12,c,12,fac,12,ipvt) call lfdrg(12,fac,12,ipvt,d1,d2) C C THE FUNCTION F IS RETURNED AS THE DETERMINATE OF ARRAY C(18,18) c F=D1*2.0**D2 c c F=D1*10.0**D2 F=D1 c RETURN C C IF IER IS GREATER THAN ZERO; THE CC(18,18) ARRAY IS WRITTEN C TO OUTPUT AND THE PROGRAN IS TERMINATED. C c 654 WRITE(6,651)IER c 651 FORMAT('ERROR IN CALCULATING DETERMINATE, IER=',I4, c *' PROGRAM TERMINATED',/,' ARRAY C(18,18) WRITTEN BELOW') c DO 667 I=1,18 c 667 WRITE(6,668)(CC(I,J),J=1,9) c 668 FORMAT(9E10.3) c WRITE(6,251) c 251 FORMAT(1X,80(1H*)) c DO 669 I=1,18 c 669 WRITE(6,668)(CC(I,J),J=10,18) 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 CRACKED LAYER * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT LOGICAL (A-Z) 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=6 IF(K.EQ.1)JC2=5 IF(K.EQ.2)JC1=2 IF(K.EQ.2)JC2=4 IF(K.EQ.3)JC1=4 IF(K.EQ.3)JC2=3 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 UNCRACKED LAYER * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT LOGICAL (A-Z) 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=6 IF(K.EQ.1)JC2=5 IF(K.EQ.2)JC1=2 IF(K.EQ.2)JC2=4 IF(K.EQ.3)JC1=4 IF(K.EQ.3)JC2=3 TP=TP+(CMPLX(CBP(IC,JC1))+PP(L)*CMPLX(CBP(IC,JC2)))*VP(K,L) 10 CONTINUE FTP=TP RETURN END