module15a/ 0040755 0000431 0000012 00000000000 07217532150 0013201 5 ustar 00rkriz staff 0000040 0000141 module15a/figs/ 0040755 0000431 0000012 00000000000 07175276644 0014151 5 ustar 00rkriz staff 0000040 0000141 module15a/ARCHIVE/ 0040755 0000431 0000012 00000000000 07217532754 0014274 5 ustar 00rkriz staff 0000040 0000141 module15a/ARCHIVE/MOD_F77/ 0040755 0000431 0000012 00000000000 07275042660 0015332 5 ustar 00rkriz staff 0000040 0000141 module15a/ARCHIVE/MOD_F77/README 0100644 0000431 0000012 00000004346 07275043452 0016216 0 ustar 00rkriz staff 0000040 0000141 --- This is the README file for ../module15a/ARCHIVE/MOD_F77 --- The original Fortran77 "legacy" code created at NIST is modified to create files necessary to calculate the singularity exponent on the radius vector, r, at an interface seperating two dissimilar anisotropic materials where this interface terminates normal to the stress-free edge of a laminated composite. This is a working archive directory where all results are confirmed prior to making this into an interactive NPIB interactive module by creating results in this directory. Instructions for creating results are outlined below. All programs in this directory work so that the viewer can download, compile and run this programs on their computers. This section is divided into two sections: 1. Compile and execute FEM programs % source /usr/local/vni/CTT3.0/ctt/bin/cttsetup.csh % $FC $FFLAGS edgesing.f -o edgesing.x $LINK_FNL % edgesing.x 2. Cleanup files in this directory % remove_excess_files.sh (Removes excess files for archiving) --------------------------------------------------------------------- CREATION OF INTERACTIVE NPIB INTERFACE FOR THIS MODULE Once verified, these working programs are copied to the NPIB directory /www/jwave/applets98/Receiver/commands/strohedge_.d except for *.out, *.dat such that the only remaining files in the strohedge_.d directory are *.f, *.x, and strohedge.html DO NOT COPY README file from the MOD_F77 directory to strohedge_.d DO NOT DELETE *.sh file FROM the MOD_F77 directory. The remove_excess_files.sh shell script file will remove the appropriate files from MOD_F77 for archiving. A shell script file named "strohedge_" recreates results in a unique directory located in http://www.jwave.vt.edu/output/"app_date-time"/.. /www/jwave/applets98/Receiver/commands/strohedge_: ------------------------------------------------------------------- #!/bin/csh #-------------------- strohedge_ ------------------------ cp /www/jwave/applets98/Receiver/commands/strohedge_.d/* $1 cd $1 #source /usr/local/vni/ipt/bin/.solaris.csh #source /usr/local/vni/wave/bin/wvsetup ./edgesing.x rm *.x exit 0 ---------------------------------------------------------------- Ronald D. Kriz, Virginia Tech, 10-07-00 module15a/ARCHIVE/MOD_F77/remove_excess_files.sh 0100755 0000431 0000012 00000000236 07314213105 0021704 0 ustar 00rkriz staff 0000040 0000141 #!/bin/csh # Go to the ../module15a/ARCHIVE/MOD_F77 directory cd ~rkriz/crcd/modules/module15a/ARCHIVE/MOD_F77 # Remove excess files for Archiving rm *.x module15a/ARCHIVE/MOD_F77/edgesing.dat 0100644 0000431 0000012 00000001064 07217451302 0017600 0 ustar 00rkriz staff 0000040 0000141 ******* Pipes-Pagano Approximation Gr/Ep ******* Young's modulus in 1-direction, E1=+2.100E+06 Young's modulus in 2-direction, E2=+2.100E+06 Young's modulus in 3-direction, E3=+2.000E+07 Shear modulus in 12-plane, G12=+0.850E+06 Shear modulus in 13-plane, G13=+0.850E+06 Shear modulus in 23-plane, G23=+0.850E+06 Poisson's ratio in 12-plane, NU21=+0.210E+00 Poisson's ratio in 13-plane, NU31=+0.210E+00 Poisson's ratio in 23-plane, NU32=+0.210E+00 Fiber orientation top region, theta =+75.0 Fiber orientation bottom region, thetap=+00.0 module15a/ARCHIVE/MOD_F77/edgesing.f 0100644 0000431 0000012 00000063525 07275042535 0017277 0 ustar 00rkriz staff 0000040 0000141 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 ,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 module15a/ARCHIVE/MOD_F77/k.txt 0100644 0000431 0000012 00000000015 07275042563 0016320 0 ustar 00rkriz staff 0000040 0000141 0.29388E-01 module15a/ARCHIVE/MOD_F77/edgesing.out 0100644 0000431 0000012 00000041520 07275042563 0017651 0 ustar 00rkriz staff 0000040 0000141 ******* Pipes-Pagano Approximation Gr/Ep ******* ******************************************************************************** ENGINEERING ELASTIC CONSTANTS ******************************************************************************** E1= 0.2100E+07 E2= 0.2100E+07 E3= 0.2000E+08 G12= 0.8500E+06 G13= 0.8500E+06 G23= 0.8500E+06 NU12= 0.2100E+00 NU13= 0.2205E-01 NU23= 0.2205E-01 NU21= 0.2100E+00 NU31= 0.2100E+00 NU32= 0.2100E+00 ******************************************************************************** STIFFNESSES CALCULATED FROM ENGINERRING CONSTANTS ******************************************************************************** C11= 0.2213E+07 C22= 0.2213E+07 C33= 0.2024E+08 C12= 0.4771E+06 C13= 0.5648E+06 C23= 0.5648E+06 C44= 0.8500E+06 C55= 0.8500E+06 C66= 0.8500E+06 ******************************************************************************** STIFFNESSES CALCULATED FROM COMPLIANCES WHICH WERE CALC. FROM ENGR. CONSTANTS ******************************************************************************** C11= 0.2213E+07 C22= 0.2213E+07 C33= 0.2024E+08 C12= 0.4771E+06 C13= 0.5648E+06 C23= 0.5648E+06 C44= 0.8500E+06 C55= 0.8500E+06 C66= 0.8500E+06 ******************************************************************************** TOP LAYER, ANGLE (DEG.)= 0.750E+02 BOTTOM LAYER, ANGLE (DEG.)= 0.000E+00 BELOW ARE STIFFNESSES FOR BOTH LAYERS ******************************************************************************** ******************************************************************************** 0.1791E+08 0.5590E+06 0.1685E+07 0.0000E+00 0.4193E+07 0.0000E+00 0.5590E+06 0.2213E+07 0.4830E+06 0.0000E+00 0.2193E+05 0.0000E+00 0.1685E+07 0.4830E+06 0.2300E+07 0.0000E+00 0.3132E+06 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8500E+06 0.0000E+00 0.0000E+00 0.4193E+07 0.2193E+05 0.3132E+06 0.0000E+00 0.1970E+07 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8500E+06 ******************************************************************************** 0.2213E+07 0.4771E+06 0.5648E+06 0.0000E+00 0.0000E+00 0.0000E+00 0.4771E+06 0.2213E+07 0.5648E+06 0.0000E+00 0.0000E+00 0.0000E+00 0.5648E+06 0.5648E+06 0.2024E+08 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8500E+06 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8500E+06 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.8500E+06 ******************************************************************************** ******************************************************************************** MATRICES FOR EIGEN-VALUE PROBLEM FOR TOP/BOTTOM LAYERS ARE WRITTEN BELOW [A]/[AP] ******************************************************************************** 0.0000E+00 -0.1658E+01 0.0000E+00 -0.2107E+02 0.0000E+00 -0.4933E+01 -0.6368E+00 0.0000E+00 -0.9913E-02 0.0000E+00 -0.3842E+00 0.0000E+00 0.0000E+00 -0.2581E-01 0.0000E+00 -0.4933E+01 0.0000E+00 -0.2318E+01 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 ******************************************************************************** 0.0000E+00 -0.1561E+01 0.0000E+00 -0.2603E+01 0.0000E+00 0.0000E+00 -0.5998E+00 0.0000E+00 0.0000E+00 0.0000E+00 -0.3842E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 -0.1000E+01 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.1000E+01 0.0000E+00 0.0000E+00 0.0000E+00 ******************************************************************************** TOP LAYER, CONVERGED= 0.852E-02 W(1)-0.477E-06 0.461E+01 Z(1) 0.937E+00 0.000E+00 Z(2)-0.121E-07 0.132E+00 Z(3) 0.243E+00 0.159E-07 Z(4) 0.195E-07-0.203E+00 Z(5) 0.287E-01 0.139E-08 Z(6) 0.470E-08-0.527E-01 W(2)-0.477E-06-0.461E+01 Z(1) 0.937E+00 0.000E+00 Z(2)-0.121E-07-0.132E+00 Z(3) 0.243E+00-0.159E-07 Z(4) 0.195E-07 0.203E+00 Z(5) 0.287E-01-0.139E-08 Z(6) 0.470E-08 0.527E-01 W(3)-0.298E-07 0.100E+01 Z(1)-0.177E+00 0.246E-06 Z(2)-0.253E-07-0.172E+00 Z(3) 0.664E+00 0.000E+00 Z(4) 0.256E-07 0.177E+00 Z(5)-0.171E+00 0.139E-07 Z(6) 0.232E-08-0.662E+00 W(4)-0.298E-07-0.100E+01 Z(1)-0.177E+00-0.246E-06 Z(2)-0.253E-07 0.172E+00 Z(3) 0.664E+00 0.000E+00 Z(4) 0.256E-07-0.177E+00 Z(5)-0.171E+00-0.139E-07 Z(6) 0.232E-08 0.662E+00 W(5) 0.939E-41 0.663E+00 Z(1) 0.711E-01-0.147E-06 Z(2) 0.168E-07 0.517E+00 Z(3)-0.182E+00-0.384E-07 Z(4)-0.239E-07-0.107E+00 Z(5) 0.780E+00 0.000E+00 Z(6) 0.491E-07 0.274E+00 W(6) 0.939E-41-0.663E+00 Z(1) 0.711E-01 0.147E-06 Z(2) 0.168E-07-0.517E+00 Z(3)-0.182E+00 0.384E-07 Z(4)-0.239E-07 0.107E+00 Z(5) 0.780E+00 0.000E+00 Z(6) 0.491E-07-0.274E+00 ******************************************************************************** BOTTOM LAYER, CONVERGED= 0.186E-01 WP(1) 0.000E+00 0.112E+01 ZP(1) 0.590E+00 0.000E+00 ZP(2) 0.000E+00 0.456E+00 ZP(3) 0.000E+00 0.000E+00 ZP(4) 0.000E+00-0.527E+00 ZP(5) 0.408E+00 0.000E+00 ZP(6) 0.000E+00 0.000E+00 WP(2) 0.000E+00-0.112E+01 ZP(1) 0.590E+00 0.000E+00 ZP(2) 0.000E+00-0.456E+00 ZP(3) 0.000E+00 0.000E+00 ZP(4) 0.000E+00 0.527E+00 ZP(5) 0.408E+00 0.000E+00 ZP(6) 0.000E+00 0.000E+00 WP(3) 0.000E+00 0.100E+01 ZP(1) 0.000E+00 0.000E+00 ZP(2) 0.000E+00 0.000E+00 ZP(3) 0.707E+00 0.000E+00 ZP(4) 0.000E+00 0.000E+00 ZP(5) 0.000E+00 0.000E+00 ZP(6) 0.000E+00-0.707E+00 WP(4) 0.000E+00-0.100E+01 ZP(1) 0.000E+00 0.000E+00 ZP(2) 0.000E+00 0.000E+00 ZP(3) 0.707E+00 0.000E+00 ZP(4) 0.000E+00 0.000E+00 ZP(5) 0.000E+00 0.000E+00 ZP(6) 0.000E+00 0.707E+00 WP(5) 0.000E+00 0.894E+00 ZP(1) 0.408E+00 0.000E+00 ZP(2) 0.000E+00 0.527E+00 ZP(3) 0.000E+00 0.000E+00 ZP(4) 0.000E+00-0.456E+00 ZP(5) 0.590E+00 0.000E+00 ZP(6) 0.000E+00 0.000E+00 WP(6) 0.000E+00-0.894E+00 ZP(1) 0.408E+00 0.000E+00 ZP(2) 0.000E+00-0.527E+00 ZP(3) 0.000E+00 0.000E+00 ZP(4) 0.000E+00 0.456E+00 ZP(5) 0.590E+00 0.000E+00 ZP(6) 0.000E+00 0.000E+00 ******************************************************************************** SUM FOR EQUATION NO.1 =-0.4438E+01 0.3642E+01 TOP LAYER SUM FOR EQUATION NO.1 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.2 =-0.3716E+00 0.2631E+02 TOP LAYER SUM FOR EQUATION NO.2 = 0.0000E+00-0.6250E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.3 =-0.5000E+00 0.9027E+00 TOP LAYER SUM FOR EQUATION NO.3 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** EIGEN-VALUE P(1)=-0.4768E-06 0.4612E+01 TOP LAYER NORMALIZED EIGEN-VECTOR V(1,1)=-0.9591E+00-0.9223E-07 V(2,1)= 0.6545E-08-0.1354E+00 V(3,1)=-0.2488E+00-0.2218E-07 ******************************************************************************** EIGEN-VALUE PP(1)= 0.0000E+00 0.1119E+01 BOTTOM LAYER NORMALIZED EIGEN-VECTOR VP(1,1)=-0.7911E+00 0.0000E+00 VP(2,1)= 0.0000E+00-0.6117E+00 VP(3,1)= 0.0000E+00 0.0000E+00 ******************************************************************************** ******************************************************************************** SUM FOR EQUATION NO.1 = 0.1250E+01-0.6260E+00 TOP LAYER SUM FOR EQUATION NO.1 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.2 = 0.4621E-01-0.9785E+00 TOP LAYER SUM FOR EQUATION NO.2 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.3 =-0.2500E+00-0.1075E+00 TOP LAYER SUM FOR EQUATION NO.3 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** EIGEN-VALUE P(2)=-0.2980E-07 0.1003E+01 TOP LAYER NORMALIZED EIGEN-VECTOR V(1,2)= 0.2505E+00-0.3619E-07 V(2,2)= 0.1974E-07 0.2423E+00 V(3,2)=-0.9373E+00-0.3279E-08 ******************************************************************************** EIGEN-VALUE PP(2)= 0.0000E+00 0.1000E+01 BOTTOM LAYER NORMALIZED EIGEN-VECTOR VP(1,2)= 0.0000E+00 0.0000E+00 VP(2,2)= 0.0000E+00 0.0000E+00 VP(3,2)=-0.1000E+01 0.0000E+00 ******************************************************************************** ******************************************************************************** SUM FOR EQUATION NO.1 =-0.2500E+00 0.2547E+00 TOP LAYER SUM FOR EQUATION NO.1 =-0.1250E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.2 =-0.2589E-01 0.4595E+00 TOP LAYER SUM FOR EQUATION NO.2 = 0.0000E+00-0.5000E+00 BOTTOM LAYER ******************************************************************************** SUM FOR EQUATION NO.3 =-0.6250E-01 0.2592E-01 TOP LAYER SUM FOR EQUATION NO.3 = 0.0000E+00 0.0000E+00 BOTTOM LAYER ******************************************************************************** EIGEN-VALUE P(3)= 0.9387E-41 0.6632E+00 TOP LAYER NORMALIZED EIGEN-VECTOR V(1,3)=-0.1286E+00 0.2863E-07 V(2,3)= 0.0000E+00-0.9355E+00 V(3,3)= 0.3292E+00-0.5896E-07 ******************************************************************************** EIGEN-VALUE PP(3)= 0.0000E+00 0.8936E+00 BOTTOM LAYER NORMALIZED EIGEN-VECTOR VP(1,3)=-0.6117E+00 0.0000E+00 VP(2,3)= 0.0000E+00-0.7911E+00 VP(3,3)= 0.0000E+00 0.0000E+00 ******************************************************************************** ******************************************************************************** T111-0.179E+08-0.169E+01 T112 0.421E+06-0.655E+00 T113-0.575E+06 0.265E+00 T121 0.756E+00-0.387E+07 T122 0.413E-01 0.420E+06 T123-0.161E-01-0.868E+06 T131-0.450E+07-0.428E+00 T132-0.801E+06-0.158E+00 T133 0.123E+06 0.387E-02 T211 0.756E+00-0.387E+07 T212 0.413E-01 0.420E+06 T213-0.161E-01-0.868E+06 T221 0.840E+06 0.158E+00 T222-0.418E+06 0.752E-02 T223 0.131E+07 0.147E-01 T231 0.188E+00-0.975E+06 T232 0.265E-01-0.799E+06 T233 0.332E-01 0.186E+06 T311-0.450E+07-0.428E+00 T312-0.801E+06-0.158E+00 T313 0.123E+06 0.387E-02 T321 0.188E+00-0.975E+06 T322 0.265E-01-0.799E+06 T323 0.332E-01 0.186E+06 T331-0.139E+07-0.117E+00 T332 0.112E+05-0.559E-01 T333 0.186E+06 0.298E-01 ******************************************************************************** TP111-0.142E+07 0.000E+00 TP112 0.000E+00 0.000E+00 TP113-0.102E+07 0.000E+00 TP121 0.000E+00-0.127E+07 TP122 0.000E+00 0.000E+00 TP123 0.000E+00-0.114E+07 TP131 0.000E+00 0.000E+00 TP132-0.850E+06 0.000E+00 TP133 0.000E+00 0.000E+00 TP211 0.000E+00-0.127E+07 TP212 0.000E+00 0.000E+00 TP213 0.000E+00-0.114E+07 TP221 0.114E+07 0.000E+00 TP222 0.000E+00 0.000E+00 TP223 0.127E+07 0.000E+00 TP231 0.000E+00 0.000E+00 TP232 0.000E+00-0.850E+06 TP233 0.000E+00 0.000E+00 TP311 0.000E+00 0.000E+00 TP312-0.850E+06 0.000E+00 TP313 0.000E+00 0.000E+00 TP321 0.000E+00 0.000E+00 TP322 0.000E+00-0.850E+06 TP323 0.000E+00 0.000E+00 TP331-0.602E+05 0.000E+00 TP332 0.000E+00 0.000E+00 TP333 0.538E+05 0.000E+00 ******************************************************************************** L=1 T1 =-0.372E+00 0.270E+02 T2 = 0.000E+00 0.903E+00 T3 = 0.118E+03 0.536E+01 L=1 T1P= 0.000E+00-0.625E+00 T2P= 0.000E+00 0.000E+00 T3P=-0.625E+00 0.000E+00 ******************************************************************************** L=2 T1 = 0.462E-01-0.100E+01 T2 =-0.188E+00-0.108E+00 T3 =-0.188E+00-0.672E+00 L=2 T1P= 0.000E+00 0.000E+00 T2P= 0.000E+00 0.000E+00 T3P= 0.000E+00 0.000E+00 ******************************************************************************** L=3 T1 =-0.259E-01 0.375E+00 T2 =-0.391E-01 0.259E-01 T3 = 0.000E+00 0.272E+00 L=3 T1P= 0.000E+00-0.625E+00 T2P= 0.000E+00 0.000E+00 T3P=-0.688E+00 0.000E+00 ******************************************************************************** ******************************************************************************** IK= 1 K= 0.000010000 DET= 0.171677327E+01 IK= 2 K= 0.010009999 DET= 0.111024725E+01 IK= 3 K= 0.020009998 DET= 0.211749864E+01 IK= 4 K= 0.030009998 DET=-0.311088443E+01 IK= 5 K= 0.040009998 DET=-0.929969311E+01 IK= 6 K= 0.050009996 DET=-0.277408648E+01 IK= 7 K= 0.060009994 DET=-0.582746029E+01 IK= 8 K= 0.070009992 DET=-0.103270364E+01 IK= 9 K= 0.080009989 DET=-0.164812732E+01 IK= 10 K= 0.090009987 DET=-0.244692111E+01 IK= 11 K= 0.100009985 DET=-0.344394541E+01 IK= 12 K= 0.110009983 DET=-0.465089035E+01 IK= 13 K= 0.120009981 DET=-0.607640314E+01 IK= 14 K= 0.130009979 DET=-0.772587872E+01 IK= 15 K= 0.140009984 DET=-0.960176849E+01 IK= 16 K= 0.150009990 DET=-0.117033947E+01 IK= 17 K= 0.160009995 DET=-0.140270019E+01 IK= 18 K= 0.170010000 DET=-0.165662038E+01 IK= 19 K= 0.180010006 DET=-0.193117595E+01 IK= 20 K= 0.190010011 DET=-0.222518635E+01 IK= 21 K= 0.200010017 DET=-0.253723025E+01 IK= 22 K= 0.210010022 DET=-0.286567736E+01 IK= 23 K= 0.220010027 DET=-0.320868683E+01 IK= 24 K= 0.230010033 DET=-0.356424475E+01 IK= 25 K= 0.240010038 DET=-0.393017888E+01 IK= 26 K= 0.250010043 DET=-0.430419064E+01 IK= 27 K= 0.260010034 DET=-0.468386698E+01 IK= 28 K= 0.270010024 DET=-0.506674194E+01 IK= 29 K= 0.280010015 DET=-0.545025253E+01 IK= 30 K= 0.290010005 DET=-0.583183289E+01 IK= 31 K= 0.300009996 DET=-0.620890141E+01 IK= 32 K= 0.310009986 DET=-0.657890797E+01 IK= 33 K= 0.320009977 DET=-0.693933105E+01 IK= 34 K= 0.330009967 DET=-0.728771687E+01 IK= 35 K= 0.340009958 DET=-0.762169361E+01 IK= 36 K= 0.350009948 DET=-0.793902588E+01 IK= 37 K= 0.360009938 DET=-0.823755360E+01 IK= 38 K= 0.370009929 DET=-0.851532078E+01 IK= 39 K= 0.380009919 DET=-0.877049541E+01 IK= 40 K= 0.390009910 DET=-0.900140572E+01 IK= 41 K= 0.400009900 DET=-0.920660210E+01 IK= 42 K= 0.410009891 DET=-0.938481617E+01 IK= 43 K= 0.420009881 DET=-0.953497982E+01 IK= 44 K= 0.430009872 DET=-0.965624332E+01 IK= 45 K= 0.440009862 DET=-0.974796009E+01 IK= 46 K= 0.450009853 DET=-0.980972481E+01 IK= 47 K= 0.460009843 DET=-0.984134007E+01 IK= 48 K= 0.470009834 DET=-0.984280968E+01 IK= 49 K= 0.480009824 DET=-0.981437588E+01 IK= 50 K= 0.490009815 DET=-0.975646782E+01 IK= 51 K= 0.500009835 DET=-0.966974449E+01 IK= 52 K= 0.510009825 DET=-0.955502701E+01 IK= 53 K= 0.520009816 DET=-0.941332340E+01 IK= 54 K= 0.530009806 DET=-0.924582672E+01 IK= 55 K= 0.540009797 DET=-0.905389309E+01 IK= 56 K= 0.550009787 DET=-0.883899975E+01 IK= 57 K= 0.560009778 DET=-0.860278797E+01 IK= 58 K= 0.570009768 DET=-0.834698105E+01 IK= 59 K= 0.580009758 DET=-0.807337475E+01 IK= 60 K= 0.590009749 DET=-0.778394413E+01 IK= 61 K= 0.600009739 DET=-0.748061848E+01 IK= 62 K= 0.610009730 DET=-0.716545582E+01 IK= 63 K= 0.620009720 DET=-0.684048843E+01 IK= 64 K= 0.630009711 DET=-0.650777817E+01 IK= 65 K= 0.640009701 DET=-0.616940308E+01 IK= 66 K= 0.650009692 DET=-0.582738781E+01 IK= 67 K= 0.660009682 DET=-0.548372889E+01 IK= 68 K= 0.670009673 DET=-0.514038563E+01 IK= 69 K= 0.680009663 DET=-0.479923773E+01 IK= 70 K= 0.690009654 DET=-0.446207523E+01 IK= 71 K= 0.700009644 DET=-0.413061285E+01 IK= 72 K= 0.710009634 DET=-0.380645823E+01 IK= 73 K= 0.720009625 DET=-0.349108815E+01 IK= 74 K= 0.730009615 DET=-0.318586063E+01 IK= 75 K= 0.740009606 DET=-0.289202237E+01 IK= 76 K= 0.750009596 DET=-0.261065626E+01 IK= 77 K= 0.760009587 DET=-0.234271193E+01 IK= 78 K= 0.770009577 DET=-0.208899641E+01 IK= 79 K= 0.780009568 DET=-0.185015452E+01 IK= 80 K= 0.790009558 DET=-0.162669408E+01 IK= 81 K= 0.800009549 DET=-0.141897297E+01 IK= 82 K= 0.810009539 DET=-0.122718787E+01 IK= 83 K= 0.820009530 DET=-0.105140173E+01 IK= 84 K= 0.830009520 DET=-0.891532707E+01 IK= 85 K= 0.840009511 DET=-0.747362232E+01 IK= 86 K= 0.850009501 DET=-0.618543196E+01 IK= 87 K= 0.860009491 DET=-0.504606724E+01 IK= 88 K= 0.870009482 DET=-0.404968643E+01 IK= 89 K= 0.880009472 DET=-0.318940544E+01 IK= 90 K= 0.890009463 DET=-0.245738339E+01 IK= 91 K= 0.900009453 DET=-0.184492910E+01 IK= 92 K= 0.910009444 DET=-0.134260845E+01 IK= 93 K= 0.920009434 DET=-0.940354538E+01 IK= 94 K= 0.930009425 DET=-0.627586460E+01 IK= 95 K= 0.940009415 DET=-0.393321848E+01 IK= 96 K= 0.950009406 DET=-0.226293993E+01 IK= 97 K= 0.960009396 DET=-0.115068686E+01 IK= 98 K= 0.970009387 DET=-0.481595325E+01 IK= 99 K= 0.980009377 DET=-0.141388774E+01 IK= 100 K= 0.990009367 DET=-0.174782121E+01 APPROXIMATION FOR ROOT, K BETWEEN 0.0 AND 1.0 KR(1)= 0.20010E-01 aaa= 0.20010E-01 bbb= 0.30010E-01 EXACT ROOT TO WITHIN 5 SIGNIFICANT FIGURES bbb= 0.29388E-01 VECTOR V(1,1)=-0.9591E+00-0.9223E-07 V(2,1)= 0.6545E-08-0.1354E+00 V(3,1)=-0.2488E+00-0.2218E-07 ****************************module15a/ARCHIVE/ORIG_NIST/ 0040755 0000431 0000012 00000000000 07226440623 0015662 5 ustar 00rkriz staff 0000040 0000141 module15a/ARCHIVE/ORIG_NIST/EDGE 0100644 0000431 0000012 00000057501 07145326737 0016327 0 ustar 00rkriz staff 0000040 0000141 PROGRAM EDGE 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* LATEST UPDATE: **** 1/21/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* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 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 IER,IERP,I,J,JJ,III,JJJ,L,IK,IR,NSIG,ITMAX C REAL A(6,6),AP(6,6),WK(50),WKP(50),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,EPS,EPS2,ETA 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 EXTERNAL FT,FTP,F C OPEN(6,FILE='TAPE6',STATUS='UNKNOWN',ERR=888) 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 E1=2.100E+06 E2=E1 E3=20.0E+06 G12=0.850E+06 G13=G12 G23=G13 NU21=0.2100 NU31=NU21 NU32=NU31 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 THETA=+75.0 THETAP=+00.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 CALL EIGRF(A,6,6,2,W,Z,6,WK,IER) CALL EIGRF(AP,6,6,2,WP,ZP,6,WKP,IERP) C C WRITE PERFORMANCE PARAMETERS IER AND IERP AND ARRAYS W,WP,Z, AND ZP C WRITE(6,4)WK(1),IER 4 FORMAT(' TOP LAYER, 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) WRITE(6,44)WKP(1),IERP 44 FORMAT(' BOTTOM LAYER, 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 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" 654 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 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) C C IMSL SUBROUTINE ZREAL1 CALCULATES IMPROVED ROOT KR(1) FROM THE INITIAL C APPROX.S FOR KR(1) GIVEN ABOVE. WE ASSUME THAT ONLY ONE ROOT EXIST. C EPS=1.0E-5 EPS2=1.0E-5 ETA=1.0E-5 NSIG=5 ITMAX=100 C CALL ZREAL1(F,EPS,EPS2,ETA,NSIG,1,KR,ITMAX,IER) WRITE(6,989)KR(1) 989 FORMAT(' EXACT ROOT TO WITHIN 5 SIGNIFICANT FIGURES',/, *1X,' KR(1)=',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. * 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,WKS(20) INTEGER I,J,L,IER 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 BY C IMSL SUBROUTINE LINV3F D1=1.0 CALL LINV3F(C,DUM,4,12,12,D1,D2,WKS,IER) IF(IER.EQ.0)GO TO 654 WRITE(6,651)IER 651 FORMAT('ERROR IN CALCULATING DETERMINATE, IER=',I4, *' PROGRAM TERMINATED, ARRAY C(12,12) WRITTEN BELOW') GO TO 889 C C THE FUNCTION F IS RETURNED AS THE DETERMINATE OF ARRAY C(12,12) 654 F=D1*2.0**D2 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 889 DO 667 I=1,12 667 WRITE(6,668)(CC(I,J),J=1,6) 668 FORMAT(6E10.3) WRITE(6,251) 251 FORMAT(1X,80(1H*)) DO 669 I=1,12 669 WRITE(6,668)(CC(I,J),J=7,12) 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 LX(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)module15a/ARCHIVE/ORIG_NIST/Other_Programs/ 0040755 0000431 0000012 00000000000 07217703237 0020620 5 ustar 00rkriz staff 0000040 0000141 module15a/ARCHIVE/ORIG_NIST/Other_Programs/MONE.FOR 0100644 0000431 0000012 00000001400 07145326740 0021716 0 ustar 00rkriz staff 0000040 0000141 PROGRAM MONE INTEGER N0,N1,N2 REAL K COMPLEX I,V1,CV,CV0,CV1,CV2 C OPEN(6,FILE='MONE.OUT',STATUS='UNKNOWN',ERR=888) C PI=ACOS(-1.) I=(0.0,1.0) N0=0 N1=1 N2=2 V1=(-1.0,0.0) WRITE(6,4)PI,V1,I 4 FORMAT(' PI=',E11.4,' V1=',2E11.4,' I=',2E11.4) K=0.0 DO 5 J=1,101 CV=V1**K CV0=(CEXP(-I*PI-2*N0*PI))**K CV1=(CEXP(-I*PI-2*N1*PI))**K CV2=(CEXP(-I*PI-2*N2*PI))**K K=K+0.01 C C WRITE RESULTS C 5 WRITE(6,10)K,CV,CV0,CV1,CV2 10 FORMAT(' K=',F5.2,' CV=',2E13.6,' CV0=',2E13.6,/, * 39X,' CV1=',2E13.6,/, * 40X,' CV2=',2E13.6) 888 STOP END module15a/ARCHIVE/ORIG_NIST/Other_Programs/MONE.OUT 0100644 0000431 0000012 00000052116 07145326741 0021752 0 ustar 00rkriz staff 0000040 0000141 PI= 0.3142E+01 V1=-0.1000E+01 0.0000E+00 I= 0.0000E+00 0.1000E+01 K= 0.01 CV= 0.100000E+01 0.000000E+00 CV0= 0.100000E+01 0.000000E+00 CV1= 0.100000E+01 0.000000E+00 CV2= 0.100000E+01 0.000000E+00 K= 0.02 CV= 0.999507E+00 0.314108E-01 CV0= 0.999507E+00 0.314108E-01 CV1= 0.938638E+00 0.294979E-01 CV2= 0.881476E+00 0.277015E-01 K= 0.03 CV= 0.998027E+00 0.627905E-01 CV0= 0.998027E+00 0.627905E-01 CV1= 0.880171E+00 0.553757E-01 CV2= 0.776233E+00 0.488364E-01 K= 0.04 CV= 0.995562E+00 0.941083E-01 CV0= 0.995562E+00 0.941083E-01 CV1= 0.824529E+00 0.779409E-01 CV2= 0.682878E+00 0.645510E-01 K= 0.05 CV= 0.992115E+00 0.125333E+00 CV0= 0.992115E+00 0.125333E+00 CV1= 0.771635E+00 0.974801E-01 CV2= 0.600153E+00 0.758169E-01 K= 0.06 CV= 0.987688E+00 0.156434E+00 CV0= 0.987688E+00 0.156434E+00 CV1= 0.721410E+00 0.114260E+00 CV2= 0.526920E+00 0.834559E-01 K= 0.07 CV= 0.982287E+00 0.187381E+00 CV0= 0.982287E+00 0.187381E+00 CV1= 0.673773E+00 0.128529E+00 CV2= 0.462156E+00 0.881609E-01 K= 0.08 CV= 0.975917E+00 0.218143E+00 CV0= 0.975917E+00 0.218143E+00 CV1= 0.628637E+00 0.140517E+00 CV2= 0.404937E+00 0.905141E-01 K= 0.09 CV= 0.968583E+00 0.248690E+00 CV0= 0.968583E+00 0.248690E+00 CV1= 0.585918E+00 0.150438E+00 CV2= 0.354435E+00 0.910034E-01 K= 0.10 CV= 0.960294E+00 0.278991E+00 CV0= 0.960294E+00 0.278991E+00 CV1= 0.545527E+00 0.158490E+00 CV2= 0.309905E+00 0.900357E-01 K= 0.11 CV= 0.951057E+00 0.309017E+00 CV0= 0.951057E+00 0.309017E+00 CV1= 0.507377E+00 0.164857E+00 CV2= 0.270680E+00 0.879492E-01 K= 0.12 CV= 0.940881E+00 0.338738E+00 CV0= 0.940881E+00 0.338738E+00 CV1= 0.471381E+00 0.169707E+00 CV2= 0.236161E+00 0.850234E-01 K= 0.13 CV= 0.929776E+00 0.368125E+00 CV0= 0.929776E+00 0.368124E+00 CV1= 0.437450E+00 0.173199E+00 CV2= 0.205815E+00 0.814881E-01 K= 0.14 CV= 0.917755E+00 0.397148E+00 CV0= 0.917755E+00 0.397148E+00 CV1= 0.405498E+00 0.175475E+00 CV2= 0.179164E+00 0.775312E-01 K= 0.15 CV= 0.904827E+00 0.425779E+00 CV0= 0.904827E+00 0.425779E+00 CV1= 0.375440E+00 0.176669E+00 CV2= 0.155781E+00 0.733050E-01 K= 0.16 CV= 0.891007E+00 0.453990E+00 CV0= 0.891007E+00 0.453990E+00 CV1= 0.347191E+00 0.176902E+00 CV2= 0.135287E+00 0.689320E-01 K= 0.17 CV= 0.876307E+00 0.481754E+00 CV0= 0.876307E+00 0.481754E+00 CV1= 0.320668E+00 0.176289E+00 CV2= 0.117342E+00 0.645096E-01 K= 0.18 CV= 0.860742E+00 0.509041E+00 CV0= 0.860742E+00 0.509041E+00 CV1= 0.295791E+00 0.174930E+00 CV2= 0.101648E+00 0.601142E-01 K= 0.19 CV= 0.844328E+00 0.535827E+00 CV0= 0.844328E+00 0.535827E+00 CV1= 0.272481E+00 0.172921E+00 CV2= 0.879347E-01 0.558050E-01 K= 0.20 CV= 0.827081E+00 0.562083E+00 CV0= 0.827081E+00 0.562083E+00 CV1= 0.250660E+00 0.170348E+00 CV2= 0.759664E-01 0.516267E-01 K= 0.21 CV= 0.809017E+00 0.587785E+00 CV0= 0.809017E+00 0.587785E+00 CV1= 0.230254E+00 0.167289E+00 CV2= 0.655325E-01 0.476121E-01 K= 0.22 CV= 0.790155E+00 0.612907E+00 CV0= 0.790155E+00 0.612907E+00 CV1= 0.211190E+00 0.163816E+00 CV2= 0.564464E-01 0.437843E-01 K= 0.23 CV= 0.770513E+00 0.637424E+00 CV0= 0.770513E+00 0.637424E+00 CV1= 0.193399E+00 0.159994E+00 CV2= 0.485432E-01 0.401585E-01 K= 0.24 CV= 0.750111E+00 0.661312E+00 CV0= 0.750111E+00 0.661312E+00 CV1= 0.176812E+00 0.155881E+00 CV2= 0.416773E-01 0.367434E-01 K= 0.25 CV= 0.728969E+00 0.684547E+00 CV0= 0.728969E+00 0.684547E+00 CV1= 0.161365E+00 0.151531E+00 CV2= 0.357197E-01 0.335430E-01 K= 0.26 CV= 0.707107E+00 0.707107E+00 CV0= 0.707107E+00 0.707107E+00 CV1= 0.146993E+00 0.146993E+00 CV2= 0.305568E-01 0.305568E-01 K= 0.27 CV= 0.684547E+00 0.728969E+00 CV0= 0.684547E+00 0.728969E+00 CV1= 0.133637E+00 0.142309E+00 CV2= 0.260887E-01 0.277816E-01 K= 0.28 CV= 0.661312E+00 0.750111E+00 CV0= 0.661312E+00 0.750111E+00 CV1= 0.121239E+00 0.137519E+00 CV2= 0.222269E-01 0.252115E-01 K= 0.29 CV= 0.637424E+00 0.770513E+00 CV0= 0.637424E+00 0.770513E+00 CV1= 0.109743E+00 0.132657E+00 CV2= 0.188941E-01 0.228391E-01 K= 0.30 CV= 0.612907E+00 0.790155E+00 CV0= 0.612907E+00 0.790155E+00 CV1= 0.990961E-01 0.127754E+00 CV2= 0.160221E-01 0.206555E-01 K= 0.31 CV= 0.587785E+00 0.809017E+00 CV0= 0.587785E+00 0.809017E+00 CV1= 0.892469E-01 0.122838E+00 CV2= 0.135509E-01 0.186512E-01 K= 0.32 CV= 0.562083E+00 0.827081E+00 CV0= 0.562083E+00 0.827080E+00 CV1= 0.801471E-01 0.117933E+00 CV2= 0.114281E-01 0.168159E-01 K= 0.33 CV= 0.535827E+00 0.844328E+00 CV0= 0.535827E+00 0.844328E+00 CV1= 0.717503E-01 0.113060E+00 CV2= 0.960778E-02 0.151394E-01 K= 0.34 CV= 0.509042E+00 0.860742E+00 CV0= 0.509042E+00 0.860742E+00 CV1= 0.640125E-01 0.108239E+00 CV2= 0.804964E-02 0.136112E-01 K= 0.35 CV= 0.481754E+00 0.876307E+00 CV0= 0.481754E+00 0.876307E+00 CV1= 0.568918E-01 0.103486E+00 CV2= 0.671852E-02 0.122209E-01 K= 0.36 CV= 0.453991E+00 0.891006E+00 CV0= 0.453991E+00 0.891006E+00 CV1= 0.503482E-01 0.988138E-01 CV2= 0.558368E-02 0.109586E-01 K= 0.37 CV= 0.425779E+00 0.904827E+00 CV0= 0.425780E+00 0.904827E+00 CV1= 0.443439E-01 0.942355E-01 CV2= 0.461831E-02 0.981440E-02 K= 0.38 CV= 0.397148E+00 0.917755E+00 CV0= 0.397148E+00 0.917755E+00 CV1= 0.388431E-01 0.897611E-01 CV2= 0.379906E-02 0.877910E-02 K= 0.39 CV= 0.368125E+00 0.929776E+00 CV0= 0.368125E+00 0.929776E+00 CV1= 0.338119E-01 0.853990E-01 CV2= 0.310559E-02 0.784381E-02 K= 0.40 CV= 0.338738E+00 0.940881E+00 CV0= 0.338738E+00 0.940881E+00 CV1= 0.292180E-01 0.811561E-01 CV2= 0.252022E-02 0.700016E-02 K= 0.41 CV= 0.309017E+00 0.951056E+00 CV0= 0.309017E+00 0.951056E+00 CV1= 0.250312E-01 0.770381E-01 CV2= 0.202760E-02 0.624029E-02 K= 0.42 CV= 0.278991E+00 0.960294E+00 CV0= 0.278992E+00 0.960294E+00 CV1= 0.212228E-01 0.730492E-01 CV2= 0.161441E-02 0.555683E-02 K= 0.43 CV= 0.248690E+00 0.968583E+00 CV0= 0.248690E+00 0.968583E+00 CV1= 0.177657E-01 0.691928E-01 CV2= 0.126913E-02 0.494294E-02 K= 0.44 CV= 0.218144E+00 0.975917E+00 CV0= 0.218144E+00 0.975917E+00 CV1= 0.146346E-01 0.654711E-01 CV2= 0.981784E-03 0.439224E-02 K= 0.45 CV= 0.187382E+00 0.982287E+00 CV0= 0.187382E+00 0.982287E+00 CV1= 0.118053E-01 0.618853E-01 CV2= 0.743748E-03 0.389885E-02 K= 0.46 CV= 0.156435E+00 0.987688E+00 CV0= 0.156435E+00 0.987688E+00 CV1= 0.925541E-02 0.584361E-01 CV2= 0.547592E-03 0.345735E-02 K= 0.47 CV= 0.125334E+00 0.992115E+00 CV0= 0.125334E+00 0.992115E+00 CV1= 0.696374E-02 0.551234E-01 CV2= 0.386916E-03 0.306274E-02 K= 0.48 CV= 0.941089E-01 0.995562E+00 CV0= 0.941089E-01 0.995562E+00 CV1= 0.491041E-02 0.519463E-01 CV2= 0.256215E-03 0.271045E-02 K= 0.49 CV= 0.627911E-01 0.998027E+00 CV0= 0.627912E-01 0.998027E+00 CV1= 0.307679E-02 0.489037E-01 CV2= 0.150764E-03 0.239630E-02 K= 0.50 CV= 0.314114E-01 0.999507E+00 CV0= 0.314114E-01 0.999507E+00 CV1= 0.144544E-02 0.459936E-01 CV2= 0.665137E-04 0.211645E-02 K= 0.51 CV= 0.655387E-06 0.100000E+01 CV0= 0.699099E-06 0.100000E+01 CV1= 0.302108E-07 0.432140E-01 CV2= 0.130553E-08 0.186745E-02 K= 0.52 CV=-0.314102E-01 0.999507E+00 CV0=-0.314101E-01 0.999507E+00 CV1=-0.127469E-02 0.405623E-01 CV2=-0.517300E-04 0.164611E-02 K= 0.53 CV=-0.627899E-01 0.998027E+00 CV0=-0.627899E-01 0.998027E+00 CV1=-0.239298E-02 0.380357E-01 CV2=-0.911985E-04 0.144957E-02 K= 0.54 CV=-0.941077E-01 0.995562E+00 CV0=-0.941076E-01 0.995562E+00 CV1=-0.336811E-02 0.356312E-01 CV2=-0.120545E-03 0.127524E-02 K= 0.55 CV=-0.125333E+00 0.992115E+00 CV0=-0.125333E+00 0.992115E+00 CV1=-0.421248E-02 0.333454E-01 CV2=-0.141583E-03 0.112075E-02 K= 0.56 CV=-0.156434E+00 0.987688E+00 CV0=-0.156434E+00 0.987688E+00 CV1=-0.493761E-02 0.311750E-01 CV2=-0.155849E-03 0.983996E-03 K= 0.57 CV=-0.187381E+00 0.982287E+00 CV0=-0.187381E+00 0.982287E+00 CV1=-0.555423E-02 0.291164E-01 CV2=-0.164635E-03 0.863052E-03 K= 0.58 CV=-0.218142E+00 0.975917E+00 CV0=-0.218142E+00 0.975917E+00 CV1=-0.607228E-02 0.271659E-01 CV2=-0.169030E-03 0.756199E-03 K= 0.59 CV=-0.248689E+00 0.968583E+00 CV0=-0.248689E+00 0.968583E+00 CV1=-0.650101E-02 0.253198E-01 CV2=-0.169944E-03 0.661889E-03 K= 0.60 CV=-0.278990E+00 0.960294E+00 CV0=-0.278990E+00 0.960294E+00 CV1=-0.684898E-02 0.235744E-01 CV2=-0.168137E-03 0.578732E-03 K= 0.61 CV=-0.309016E+00 0.951057E+00 CV0=-0.309016E+00 0.951057E+00 CV1=-0.712410E-02 0.219258E-01 CV2=-0.164240E-03 0.505481E-03 K= 0.62 CV=-0.338737E+00 0.940881E+00 CV0=-0.338737E+00 0.940881E+00 CV1=-0.733372E-02 0.203702E-01 CV2=-0.158776E-03 0.441020E-03 K= 0.63 CV=-0.368124E+00 0.929777E+00 CV0=-0.368124E+00 0.929777E+00 CV1=-0.748459E-02 0.189040E-01 CV2=-0.152174E-03 0.384350E-03 K= 0.64 CV=-0.397147E+00 0.917755E+00 CV0=-0.397147E+00 0.917755E+00 CV1=-0.758294E-02 0.175232E-01 CV2=-0.144785E-03 0.334580E-03 K= 0.65 CV=-0.425778E+00 0.904827E+00 CV0=-0.425778E+00 0.904827E+00 CV1=-0.763454E-02 0.162243E-01 CV2=-0.136893E-03 0.290914E-03 K= 0.66 CV=-0.453990E+00 0.891007E+00 CV0=-0.453990E+00 0.891007E+00 CV1=-0.764465E-02 0.150035E-01 CV2=-0.128727E-03 0.252641E-03 K= 0.67 CV=-0.481753E+00 0.876307E+00 CV0=-0.481753E+00 0.876307E+00 CV1=-0.761813E-02 0.138574E-01 CV2=-0.120468E-03 0.219131E-03 K= 0.68 CV=-0.509040E+00 0.860743E+00 CV0=-0.509040E+00 0.860743E+00 CV1=-0.755943E-02 0.127823E-01 CV2=-0.112260E-03 0.189822E-03 K= 0.69 CV=-0.535826E+00 0.844329E+00 CV0=-0.535826E+00 0.844329E+00 CV1=-0.747262E-02 0.117750E-01 CV2=-0.104213E-03 0.164214E-03 K= 0.70 CV=-0.562082E+00 0.827081E+00 CV0=-0.562082E+00 0.827081E+00 CV1=-0.736142E-02 0.108320E-01 CV2=-0.964102E-04 0.141864E-03 K= 0.71 CV=-0.587784E+00 0.809018E+00 CV0=-0.587784E+00 0.809018E+00 CV1=-0.722923E-02 0.995021E-02 CV2=-0.889132E-04 0.122379E-03 K= 0.72 CV=-0.612906E+00 0.790156E+00 CV0=-0.612906E+00 0.790156E+00 CV1=-0.707914E-02 0.912639E-02 CV2=-0.817649E-04 0.105411E-03 K= 0.73 CV=-0.637423E+00 0.770514E+00 CV0=-0.637423E+00 0.770514E+00 CV1=-0.691396E-02 0.835756E-02 CV2=-0.749939E-04 0.906523E-04 K= 0.74 CV=-0.661311E+00 0.750112E+00 CV0=-0.661311E+00 0.750112E+00 CV1=-0.673623E-02 0.764078E-02 CV2=-0.686165E-04 0.778304E-04 K= 0.75 CV=-0.684546E+00 0.728970E+00 CV0=-0.684546E+00 0.728970E+00 CV1=-0.654827E-02 0.697322E-02 CV2=-0.626399E-04 0.667049E-04 K= 0.76 CV=-0.707106E+00 0.707108E+00 CV0=-0.707106E+00 0.707108E+00 CV1=-0.635215E-02 0.635217E-02 CV2=-0.570634E-04 0.570635E-04 K= 0.77 CV=-0.728968E+00 0.684548E+00 CV0=-0.728968E+00 0.684548E+00 CV1=-0.614975E-02 0.577501E-02 CV2=-0.518808E-04 0.487194E-04 K= 0.78 CV=-0.750110E+00 0.661313E+00 CV0=-0.750110E+00 0.661313E+00 CV1=-0.594274E-02 0.523924E-02 CV2=-0.470813E-04 0.415078E-04 K= 0.79 CV=-0.770512E+00 0.637425E+00 CV0=-0.770512E+00 0.637425E+00 CV1=-0.573263E-02 0.474246E-02 CV2=-0.426508E-04 0.352840E-04 K= 0.80 CV=-0.790154E+00 0.612908E+00 CV0=-0.790154E+00 0.612908E+00 CV1=-0.552075E-02 0.428235E-02 CV2=-0.385731E-04 0.299205E-04 K= 0.81 CV=-0.809016E+00 0.587786E+00 CV0=-0.809016E+00 0.587786E+00 CV1=-0.530831E-02 0.385672E-02 CV2=-0.348301E-04 0.253057E-04 K= 0.82 CV=-0.827080E+00 0.562085E+00 CV0=-0.827080E+00 0.562085E+00 CV1=-0.509635E-02 0.346348E-02 CV2=-0.314029E-04 0.213415E-04 K= 0.83 CV=-0.844327E+00 0.535828E+00 CV0=-0.844327E+00 0.535828E+00 CV1=-0.488579E-02 0.310063E-02 CV2=-0.282721E-04 0.179421E-04 K= 0.84 CV=-0.860741E+00 0.509043E+00 CV0=-0.860741E+00 0.509043E+00 CV1=-0.467745E-02 0.276625E-02 CV2=-0.254183E-04 0.150324E-04 K= 0.85 CV=-0.876306E+00 0.481755E+00 CV0=-0.876306E+00 0.481755E+00 CV1=-0.447203E-02 0.245853E-02 CV2=-0.228220E-04 0.125466E-04 K= 0.86 CV=-0.891006E+00 0.453992E+00 CV0=-0.891006E+00 0.453992E+00 CV1=-0.427014E-02 0.217575E-02 CV2=-0.204646E-04 0.104273E-04 K= 0.87 CV=-0.904826E+00 0.425781E+00 CV0=-0.904826E+00 0.425781E+00 CV1=-0.407230E-02 0.191629E-02 CV2=-0.183279E-04 0.862451E-05 K= 0.88 CV=-0.917754E+00 0.397149E+00 CV0=-0.917754E+00 0.397150E+00 CV1=-0.387894E-02 0.167857E-02 CV2=-0.163945E-04 0.709459E-05 K= 0.89 CV=-0.929776E+00 0.368126E+00 CV0=-0.929776E+00 0.368126E+00 CV1=-0.369043E-02 0.146115E-02 CV2=-0.146479E-04 0.579956E-05 K= 0.90 CV=-0.940880E+00 0.338740E+00 CV0=-0.940880E+00 0.338740E+00 CV1=-0.350708E-02 0.126263E-02 CV2=-0.130725E-04 0.470640E-05 K= 0.91 CV=-0.951056E+00 0.309019E+00 CV0=-0.951056E+00 0.309019E+00 CV1=-0.332912E-02 0.108170E-02 CV2=-0.116534E-04 0.378646E-05 K= 0.92 CV=-0.960293E+00 0.278993E+00 CV0=-0.960293E+00 0.278993E+00 CV1=-0.315675E-02 0.917127E-03 CV2=-0.103771E-04 0.301485E-05 K= 0.93 CV=-0.968583E+00 0.248692E+00 CV0=-0.968583E+00 0.248692E+00 CV1=-0.299010E-02 0.767733E-03 CV2=-0.923070E-05 0.237006E-05 K= 0.94 CV=-0.975916E+00 0.218145E+00 CV0=-0.975916E+00 0.218145E+00 CV1=-0.282927E-02 0.632422E-03 CV2=-0.820230E-05 0.183345E-05 K= 0.95 CV=-0.982287E+00 0.187383E+00 CV0=-0.982287E+00 0.187383E+00 CV1=-0.267431E-02 0.510158E-03 CV2=-0.728092E-05 0.138892E-05 K= 0.96 CV=-0.987688E+00 0.156436E+00 CV0=-0.987688E+00 0.156436E+00 CV1=-0.252526E-02 0.399967E-03 CV2=-0.645644E-05 0.102261E-05 K= 0.97 CV=-0.992114E+00 0.125335E+00 CV0=-0.992114E+00 0.125335E+00 CV1=-0.238210E-02 0.300935E-03 CV2=-0.571952E-05 0.722556E-06 K= 0.98 CV=-0.995562E+00 0.941103E-01 CV0=-0.995562E+00 0.941104E-01 CV1=-0.224481E-02 0.212202E-03 CV2=-0.506164E-05 0.478476E-06 K= 0.99 CV=-0.998027E+00 0.627925E-01 CV0=-0.998027E+00 0.627926E-01 CV1=-0.211332E-02 0.132964E-03 CV2=-0.447497E-05 0.281551E-06 K= 1.00 CV=-0.999506E+00 0.314128E-01 CV0=-0.999506E+00 0.314129E-01 CV1=-0.198757E-02 0.624661E-04 CV2=-0.395238E-05 0.124217E-06 K= 1.01 CV=-0.100000E+01 0.205979E-05 CV0=-0.100000E+01 0.214721E-05 CV1=-0.186745E-02 0.400981E-08 CV2=-0.348737E-05 0.748812E-11 4 0.144957E-02 K= 0.54 CV=-0.941077E-01 0.995562E+00 CV0=-0.941076E-01 0.995562E+00 CV1=-0.336811E-02 0.356312E-01 CV2=-0.120545E-03 0.127524E-02 K= 0.55 CV=-0.125333E+00 0.992115E+00 CV0=-0.125333E+00 0.992115E+00 CV1=-0.421248E-02 0.333454E-01 CV2=-0.141583E-03 0.1120module15a/ARCHIVE/ORIG_NIST/Other_Programs/README 0100644 0000431 0000012 00000000505 07175311133 0021466 0 ustar 00rkriz staff 0000040 0000141 -- This is the README file for ../module??/ARCHIVE/ORIG_NIST -- These are the original DATA & Fortran77 "legacy" code created at NIST in the 1980s by Ron Kriz. These files are modfied here for web-based CRCD access. ---------------------------------------------------------------- Ronald D. Kriz, Virginia Tech, 10-07-00 module15a/ARCHIVE/ORIG_NIST/Other_Programs/ROOTS 0100644 0000431 0000012 00000005264 07145326742 0021457 0 ustar 00rkriz staff 0000040 0000141 PROGRAM ROOTS INTEGER NDEG,IER REAL COEF(7),M,N REAL NU12,NU13,NU23,NU21,NU31,NU32 COMPLEX Z(6) C OPEN(6,FILE='TAPE6',STATUS='UNKNOWN',ERR=888) C NDEG=6 PI=ACOS(-1.) C C CALCULATE STIFFNESSES FROM ENGINEERING ELASTIC PROPERTIES C REPORTED IN VPI-E-77 WHERE THE FIBER AXIS IS CHOOSEN AS X3 C (MARCH,77) E1=1.453E+06 E2=E1 E3=20.0E+06 G12=0.4872E+06 G13=0.8394E+06 G23=G13 NU21=0.4906 NU31=0.3056 NU32=NU31 NU12=NU21*E1/E2 NU13=NU31*E1/E3 NU23=NU32*E2/E3 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 LOOP THRU FIBER ANGLES(DEG.) 0 TO 90 IN STEPS OF 10 DEG. C DTHETA=10.0 THETA=0.0 DO 1000 K=1,10 C C CALCULATE STIFFNESSES IN THE ROTATED X3-X1 PLANE C THETAR=PI*THETA/180.0 M=COS(THETAR) N=SIN(THETAR) CB11=M**4*C11+2.*(M*N)**2*(C13+2.*C55)+C33*N**4 CB13=(N*M)**2*(C11+C33-4.*C55)+(N**4+M**4)*C13 CB12=M**2*C12+N**2*C23 CB15=-M*N*(C11*M**2-C33*N**2-(C13+2.*C55)*(M**2-N**2)) CB33=N**4*C11+2.*(N*M)**2*(C13+2.*C55)+M**4*C33 CB23=N**2*C12+M**2*C23 CB22=C22 CB25=M*N*(C23-C12) CB35=-M*N*(N**2*C11-M**2*C33-(C13+2.*C55)*(N**2-M**2)) CB44=C44*M**2+C66*N**2 CB46=M*N*(C44-C66) CB66=C66*M**2+C44*N**2 CB55=(M*N)**2*(C11+C33-2.*C13)+C55*(M**2-N**2)**2 C C CLEAR ARRAY COEF(I) C DO 2 I=1,NDEG 2 COEF(I)=0.0 C C DEFINE COEF(I) NONZERO TERMS C A=CB66 B=CB55 C=CB46+CB25 D=CB46 E=CB35 F=CB22 G=CB44 H=CB23+CB44 M=CB44 N=CB33 COEF(1)=B*G*N-G*E**2 COEF(3)=A*G*N+B*G*M+B*F*N+2.*C*H*E-B*H**2-N*C**2-2.*E*D*G-F*E**2 COEF(5)=A*G*M+A*F*N+F*M*B+2.*C*H*D-A*H**2-M*C**2-G*D**2-2.*E*D*F COEF(7)=A*F*M-F*D**2 C C WRITE COEF(I,J) TO TAPE6 C DO 200 I=1,NDEG WRITE(6,250)I,COEF(I) 200 CONTINUE 250 FORMAT(' COEF(',I1,')=',E11.4) C C SOLVE FOR ROOTS C CALL ZRPOLY(COEF,NDEG,Z,IER) C C WRITE OUT PERFORMANCE PARAMETER C WRITE(6,4)IER 4 FORMAT(' PERFORMANCE IER=',I6) C C WRITE RESULTS C WRITE(6,14)THETA 14 FORMAT(' ANGLE (DEG.)=',E10.3) DO 25 J=1,NDEG WRITE(6,15)J,Z(J) 15 FORMAT(' ROOT Z(',I1,')=',2E10.3) 25 CONTINUE THETA=THETA+DTHETA 1000 CONTINUE 888 STOP END module15a/ARCHIVE/ORIG_NIST/Other_Programs/RUN.BAT 0100644 0000431 0000012 00000000400 07145326742 0021605 0 ustar 00rkriz staff 0000040 0000141 REM - filename is RUN.BAT PATH C:\FORTRAN SET PROFORT.ERR=\FORTRAN\PROFORT.ERR PROFORT %1 %2 > %1.LST LINK %1,,NUL,\FORTRAN\LIBRARY\PROFORT.LIB %1 DEL %1.OBJ PATH C:\;C:\DOS;C:\KRM ORTRAN\PROFORT.ERR PROFORT %1 %2 > %1.LST LINK %1,,NUL,\FORTRAN\ module15a/ARCHIVE/ORIG_NIST/Other_Programs/SING 0100644 0000431 0000012 00000061304 07145326743 0021307 0 ustar 00rkriz staff 0000040 0000141 PROGRAM SING 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* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 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 IER,IERP,I,J,JJ,III,JJJ,L,IK,IR,NSIG,ITMAX C REAL A(6,6),AP(6,6),WK(50),WKP(50),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,EPS,EPS2,ETA 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 EXTERNAL FT,FTP,F C OPEN(6,FILE='TAPE6',STATUS='UNKNOWN',ERR=888) 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 E1=2.100E+06 E2=E1 E3=20.0E+06 G12=0.850E+06 G13=G12 G23=G13 NU21=0.2100 NU31=NU21 NU32=NU31 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 THETA=+90.0 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 CALL EIGRF(A,6,6,2,W,Z,6,WK,IER) CALL EIGRF(AP,6,6,2,WP,ZP,6,WKP,IERP) C C WRITE PERFORMANCE PARAMETERS IER AND IERP AND ARRAYS W,WP,Z, AND ZP C WRITE(6,4)WK(1),IER 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) WRITE(6,44)WKP(1),IERP 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) 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 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 EPS=1.0E-5 EPS2=1.0E-5 ETA=1.0E-2 NSIG=5 ITMAX=100 C CALL ZREAL1(F,EPS,EPS2,ETA,NSIG,3,KR,ITMAX,IER) WRITE(6,989)(I,KR(I),I=1,3) 989 FORMAT(' EXACT ROOTS TO WITHIN 5 SIGNIFICANT FIGURES',/, *1X,3(' KR(',I1,')=',F12.9)) 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 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 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,WKS(20) INTEGER I,J,L,IER 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 D1=1.0 CALL LINV3F(C,DUM,4,18,18,D1,D2,WKS,IER) IF(IER.NE.0)GO TO 654 C C THE FUNCTION F IS RETURNED AS THE DETERMINATE OF ARRAY C(18,18) F=D1*2.0**D2 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 654 WRITE(6,651)IER 651 FORMAT('ERROR IN CALCULATING DETERMINATE, IER=',I4, *' PROGRAM TERMINATED',/,' ARRAY C(18,18) WRITTEN BELOW') DO 667 I=1,18 667 WRITE(6,668)(CC(I,J),J=1,9) 668 FORMAT(9E10.3) WRITE(6,251) 251 FORMAT(1X,80(1H*)) DO 669 I=1,18 669 WRITE(6,668)(CC(I,J),J=10,18) 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 ',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,/, * module15a/ARCHIVE/ORIG_NIST/Other_Programs/SINGT 0100644 0000431 0000012 00000002621 07145326743 0021430 0 ustar 00rkriz staff 0000040 0000141 PROGRAM SING INTEGER N,IA,IJOB,IZ,IER REAL A(4,4),WK(24) COMPLEX W(4),Z(4,4),ZN C OPEN(6,FILE='TAPE6',STATUS='UNKNOWN',ERR=888) C IA=4 IZ=4 N=4 IJOB=2 C C CLEAR ARRAY A(I,J) C DO 2 I=1,N DO 2,J=1,N 2 A(I,J)=0.0 C C DEFINE A(I,J) NONZERO TERMS C A(1,1)=4.0 A(2,2)=4.0 A(3,3)=4.0 A(4,4)=4.0 A(1,4)=3.0 A(4,1)=3.0 A(1,2)=-5.0 A(3,1)=5.0 A(4,3)=5.0 A(2,4)=-5.0 A(2,3)=-3.0 A(3,2)=-3.0 C C SOLVE FOR EIGEN VALUES AND EIGEN VECTORS C CALL EIGRF(A,N,IA,IJOB,W,Z,IZ,WK,IER) C C WRITE OUT PERFORMANCE PARAMETER C WRITE(6,4)WK(1) 4 FORMAT(' PERFORMANCE=',E10.3) C C NORMALIZE EIGEN-VECTORS C C DO 5 J=1,N C ZN=Z(1,J) C DO 5 I=1,N C Z(I,J)=Z(I,J)/ZN C 5 CONTINUE C C WRITE RESULTS C DO 25 J=1,N WRITE(6,15)J,W(J) 15 FORMAT(' EIGEN-VALUES W(',I1,')=',2E10.3) WRITE(6,20)(I,J,Z(I,J),I=1,N) C 20 FORMAT(' NORMALIZED EIGEN-VECTORS Z(',I1,',',I1,')=',2E10.3,/, 20 FORMAT(' EIGEN-VECTORS Z(',I1,',',I1,')=',2E10.3,/, * 25X,' Z(',I1,',',I1,')=',2E10.3,/, * 25X,' Z(',I1,',',I1,')=',2E10.3,/, * 25X,' Z(',I1,',',I1,')=',2E10.3) 25 CONTINUE 888 STOP END module15a/ARCHIVE/ORIG_NIST/Other_Programs/STIFF1.FOR 0100644 0000431 0000012 00000002400 07145326743 0022120 0 ustar 00rkriz staff 0000040 0000141 PROGRAM STIFF1 REAL NU12,NU13,NU23,NU21,NU31,NU32 OPEN(6,FILE='STIFF1.OUT',STATUS='UNKNOWN',ERR=888) E1=20.0E+06 E2=1.453E+06 E3=E2 G12=0.8394E+06 G13=G12 G23=0.4872E+06 NU12=0.3056 NU13=NU12 NU23=0.4906 NU21=NU12*E2/E1 NU31=NU13*E3/E1 NU32=NU23*E3/E2 DEL=1.-NU12*NU21-NU13*NU31-NU23*NU32-2.*NU21*NU32*NU13 C33=(1.-NU23*NU32)*E1/DEL C11=(1.-NU13*NU31)*E2/DEL C22=(1.-NU12*NU21)*E3/DEL C13=(NU21+NU31*NU23)*E1/DEL C23=(NU31+NU21*NU32)*E1/DEL C12=(NU32+NU12*NU31)*E2/DEL C66=G23 C44=G13 C55=G12 WRITE(6,101)C11,C22,C33,C12,C13,C23,C44,C55,C66,DEL 101 FORMAT(' STIFFNESSES C11=',E11.4,/, * 12X,' C22=',E11.4,/, * 12X,' C33=',E11.4,/, * 12X,' C12=',E11.4,/, * 12X,' C13=',E11.4,/, * 12X,' C23=',E11.4,/, * 12X,' C44=',E11.4,/, * 12X,' C55=',E11.4,/, * 12X,' C66=',E11.4,/, * 12X,' DEL=',E11.4) 888 STOP END C55=',E11.4,/, * 12X,' C66=',E11.4,/, * 12X,' D module15a/ARCHIVE/ORIG_NIST/Other_Programs/STIFF1.OUT 0100644 0000431 0000012 00000000456 07145326744 0022153 0 ustar 00rkriz staff 0000040 0000141 STIFFNESSES C11= 0.1953E+07 C22= 0.1953E+07 C33= 0.2055E+08 C12= 0.9778E+06 C13= 0.8955E+06 C23= 0.8955E+06 C44= 0.8394E+06 C55= 0.8394E+06 C66= 0.4872E+06 DEL= 0.7391E+00 module15a/ARCHIVE/ORIG_NIST/Other_Programs/STIFF2.FOR 0100644 0000431 0000012 00000002400 07145326744 0022122 0 ustar 00rkriz staff 0000040 0000141 PROGRAM STIFF2 REAL NU12,NU13,NU23,NU21,NU31,NU32 OPEN(6,FILE='STIFF2.OUT',STATUS='UNKNOWN',ERR=888) E1=1.453E+06 E2=E1 E3=20.0E+06 G23=0.8394E+06 G13=G23 G12=0.4872E+06 NU21=0.4906 NU31=0.3056 NU32=NU31 NU12=NU21*E1/E2 NU13=NU31*E1/E3 NU23=NU32*E2/E3 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 WRITE(6,101)C11,C22,C33,C12,C13,C23,C44,C55,C66,DEL 101 FORMAT(' STIFFNESSES C11=',E11.4,/, * 12X,' C22=',E11.4,/, * 12X,' C33=',E11.4,/, * 12X,' C12=',E11.4,/, * 12X,' C13=',E11.4,/, * 12X,' C23=',E11.4,/, * 12X,' C44=',E11.4,/, * 12X,' C55=',E11.4,/, * 12X,' C66=',E11.4,/, * 12X,' DEL=',E11.4) 888 STOP END C55=',E11.4,/, * 12X,' C66=',E11.4,/, * 12X,' D module15a/ARCHIVE/ORIG_NIST/Other_Programs/STIFF2.OUT 0100644 0000431 0000012 00000000456 07145326746 0022156 0 ustar 00rkriz staff 0000040 0000141 STIFFNESSES C11= 0.1953E+07 C22= 0.1953E+07 C33= 0.2055E+08 C12= 0.9778E+06 C13= 0.8955E+06 C23= 0.8955E+06 C44= 0.8394E+06 C55= 0.8394E+06 C66= 0.4872E+06 DEL= 0.7391E+00 module15a/ARCHIVE/ORIG_NIST/Other_Programs/TEST 0100644 0000431 0000012 00000031054 07145326746 0021330 0 ustar 00rkriz staff 0000040 0000141 PROGRAM SING 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 * PROGRAM * 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* PROGRAM WRITTEN BY R.D.KRIZ, FRACTURE AND DEFORMATION DIVISION, * C* NATIONAL BUREAU OF STANDARDS, BOULDER, COLORADO 80303. 5/14/86 * C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * COMMON /L1/CB(6,6),PL(3),V(3,3) COMMON /L2/CBP(6,6),PLP(3),VP(3,3) INTEGER NN,IA,IJOB,IZ,IER REAL A(6,6),AP(6,6),WK(50),WKP(50),WKS(20),KK REAL NU12,NU13,NU23,NU21,NU31,NU32,N,M,C(19,19),CB,CBP COMPLEX W(6),U(3),Z(6,6),SUM,V,B(3,3) COMPLEX WP(6),UP(3),ZP(6,6),SUMP,VP,BP(3,3) COMPLEX FTT,FTP,CI,CV,PL,PLP,TT(3,3,3),TP(3,3,3),CN EXTERNAL FTP,FTT C OPEN(6,FILE='TAPE6',STATUS='UNKNOWN',ERR=888) C C DEFINE INITIAL PARAMETERS C C PARAMETERS USED IN CALCULATING EIGEN-VALUES IA=6 IZ=6 NN=6 IJOB=2 IFAIL=1 C C GEOMETRIC PI(RADIANS)=180(DEGREES) PI=ACOS(-1.) C C IMAGINARY NUMBER "I" CI=(0.0,1.0) C C INTEGER MULTIPLIER FOR CALCULATING HIGHER ORDER C IMAGINARY NUMBERS CV, PL, AND PLP NR=0.0 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 E1=1.453E+06 E2=E1 E3=20.0E+06 G12=0.4872E+06 G13=0.8394E+06 G23=G13 NU21=0.4906 NU31=0.3056 NU32=NU31 C C CALCULATE OTHER POISSONS RATIOS C NU12=NU21*E1/E2 NU13=NU31*E1/E3 NU23=NU32*E2/E3 C C ELASTIC PROPERTY APPROXIMATION AFTER PIPES 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 CALCULATE OTHER POISSONS RATIOS C NU12=NU21*E1/E2 C NU13=NU31*E1/E3 C NU23=NU32*E2/E3 C C SAME APPROXIMATION FOR OTHER POISSONS RATIOS C NU13=NU31*E1/E3 C NU12=NU13 C NU23=NU12 C C WRITE ENGINEERING ELASTIC CONSTANTS 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 OTHER ELASTIC 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 UNROTATED STIFFNESSES C WRITE(6,251) WRITE(6,261) 261 FORMAT(' STIFFNESSES FOR PRINCIPAL MATERIAL AXES') 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,NN DO 2 I=1,NN 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 THETA=30.0 THETAP=60.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*(C11*M**2-C33*N**2-(C13+2.*C55)*(M**2-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*(N**2*C11-M**2*C33-(C13+2.*C55)*(N**2-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*(C11*M**2-C33*N**2-(C13+2.*C55)*(M**2-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*(N**2*C11-M**2*C33-(C13+2.*C55)*(N**2-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 DEFINE A(I,J) NONZERO TERMS 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 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 TAPE6 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,NN WRITE(6,250)(A(I,J),J=1,NN) 200 CONTINUE WRITE(6,251) DO 201 I=1,NN WRITE(6,250)(AP(I,J),J=1,NN) 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 C CALL EIGRF(A,NN,IA,IJOB,W,Z,IZ,WK,IER) CALL EIGRF(AP,NN,IA,IJOB,WP,ZP,IZ,WKP,IERP) C C WRITE OUT PERFORMANCE PARAMETER C WRITE(6,4)WK(1),IER 4 FORMAT('LAYER WITH CRACK, CONVERGED=',E10.3,' IER=',I6) DO 41 J=1,NN 41 WRITE(6,42)J,W(J),(I,J,Z(I,J),I=4,NN) 42 FORMAT(' W(',I1,')=',2E10.3,3(' Z(',I1,',',I1,')=',2E10.3)) WRITE(6,44)WKP(1),IERP 44 FORMAT('LAYER WITHOUT CRACK CONVERGED=',E10.3,' IERP=',I6) DO 441 J=1,NN 441 WRITE(6,442)J,WP(J),(I,J,ZP(I,J),I=4,NN) 442 FORMAT(' WP(',I1,')=',2E10.3,3(' ZP(',I1,',',I1,')=',2E10.3)) C C NORMALIZE EIGEN-VECTORS AND CHECK RESULTS C JJ=0 DO 65 J=1,5,2 JJ=JJ+1 C C NORMALIZE EIGEN-VECTORS AND REASSIGN INDEX NUMBERS (1=1,2=3,3=5) C RECOVER ONLY THE COMPLEX CONJUGATES WITH POSITIVE IMAGINARY TERMS C SUM=0.0 SUMP=0.0 DO 50 I=4,NN SUM=SUM+Z(I,J)*CONJG(Z(I,J)) SUMP=SUMP+ZP(I,J)*CONJG(ZP(I,J)) 50 CONTINUE DO 51 I=4,NN V(I-3,JJ)=Z(I,J)/CSQRT(SUM) VP(I-3,JJ)=ZP(I,J)/CSQRT(SUMP) 51 CONTINUE C C REASSIGN INDEX NUMBERS FOR EIGEN-VALUES (1=1,2=3,3=5) C RECOVER ONLY THE COMPLEX CONJUGATES WITH POSITIVE IMAGINARY TERMS C U(JJ)=W(J) UP(JJ)=WP(J) C C CHECK RESULTS C B(1,1)=CMPLX(CB(6,6))+U(JJ)*U(JJ)*CMPLX(CB(5,5)) B(1,2)=(CMPLX(CB(4,6))+CMPLX(CB(2,5)))*U(JJ) B(1,3)=CMPLX(CB(4,6))+U(JJ)*U(JJ)*CMPLX(CB(3,5)) B(2,2)=CMPLX(CB(2,2))+U(JJ)*U(JJ)*CMPLX(CB(4,4)) B(2,3)=(CMPLX(CB(2,3))+CMPLX(CB(4,4)))*U(JJ) B(3,3)=CMPLX(CB(4,4))+U(JJ)*U(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))+UP(JJ)*UP(JJ)*CMPLX(CBP(5,5)) BP(1,2)=(CMPLX(CBP(4,6))+CMPLX(CBP(2,5)))*UP(JJ) BP(1,3)=CMPLX(CBP(4,6))+UP(JJ)*UP(JJ)*CMPLX(CBP(3,5)) BP(2,2)=CMPLX(CBP(2,2))+UP(JJ)*UP(JJ)*CMPLX(CBP(4,4)) BP(2,3)=(CMPLX(CBP(2,3))+CMPLX(CBP(4,4)))*UP(JJ) BP(3,3)=CMPLX(CBP(4,4))+UP(JJ)*UP(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,U(JJ),(I,JJ,V(I,JJ),I=1,3) 63 FORMAT(' EIGEN-VALUE U(',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,UP(JJ),(I,JJ,VP(I,JJ),I=1,3) 64 FORMAT(' EIGEN-VALUE UP(',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 888 STOP END C C CALCULATE OTHER ELASTIC 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 UNROTATED STIFFNESSES C WRITE(6,251) WRITE(6,261) 261 FORMAT(' STIFFNESSES FOR PRImodule15a/ARCHIVE/ORIG_NIST/README 0100644 0000431 0000012 00000000506 07226440661 0016542 0 ustar 00rkriz staff 0000040 0000141 -- This is the README file for ../module15a/ARCHIVE/ORIG_NIST -- These are the original DATA & Fortran77 "legacy" code created at NIST in the 1980s by Ron Kriz. These files are modfied here for web-based CRCD access. ---------------------------------------------------------------- Ronald D. Kriz, Virginia Tech, 10-07-00 module15a/ARCHIVE/README 0100644 0000431 0000012 00000001106 07314213022 0015125 0 ustar 00rkriz staff 0000040 0000141 ------ This is the README file for ../module15a/ARCHIVE ------ This is a "working" ARCHIVE that demonstrates how existing "legacy" code is modified to run interactively on the CRCD web-site as an NPIB applications. FIRST modification: modify the Fortran legacy code to run on the CRCD web server as a Fortran program. ~rkriz/crcd/modules/module15a/ARCHIVE/MOD_F77 FINAL modification: Copy files from the MOD_F77 directory to the ../Receiver/commands/app_.d directory -------------------------------------------------------------- Ronald D. Kriz, Virginia Tech, 10-07-00 module15a/start.html 0100644 0000431 0000012 00000000443 07313002241 0015210 0 ustar 00rkriz staff 0000040 0000141
Instructions for using this form:
Below is a brief description of a method to calculate singuariltes that exists at an interface seperating two dissimilar anisotropic materials that terminates at a laminate stress free edge. A more complete description is give in Section 7. "Cracks Near Interfaces Between Dissimilar Anisotropic Materials"
It is important to note that singularities coincide with the coordinate system origin. Hence each interface singularity has a unique coordinate system. Since the objective here is to reproduce published results the same coordinate system established by Ting et al. will be used. For this module we have expanded on the solution introduced by Ting et al. and developed the numerical approach to solve for these singularities. The solution is organized into three parts: 1) transformation of fourth order stiffness tensor, Cijkl, for arbitrary rotation in the 1-3 plane, see eqn 35 Ref.[25], 2) calcuate eigenvalues, pk, and eigen vectors, ni and tij, from equations (9) and (10) Ref.[25], and 3) calcuate exponent, k (singularity) on radius, r, from the determinate of equations (30 and (31) of Ref.[25] using appropriate boundary conditions. From these boundary conditions 12 homogeneous equations are generated. Equation coefficients are functions of elastic properties and the exponent, k on the radius vector whose origin is located at the stress free edge where the singularity exists. With the elastic properties prescribed in the NPIB form the determinate of the coefficient matrix is only a function of the exponent, k. For comparison results from Ting et al. Ref.[25] are provided below in Table 1 using the Pipes-Pagano elastic property approximation which are the default elastic properties in the NPIB form.
Exponent of r-k |
q'=0 | q'=90 | q'=-q |
---|---|---|---|
q=0 | - | 3.3388 x 10-2 | - |
q=15 | 1.3528 x 10-4 | 3.2814 x 10-2 | 6.4322 x 10-4 |
q=30 | 2.6286 x 10-3 | 2.8682 x 10-2 | 1.1658 x 10-2 |
q=45 | 9.6461 x 10-3 | 2.0575 x 10-2 | 2.5575 x 10-2 |
q=60 | 1.9866 x 10-2 | 1.0519 x 10-2 | 2.3346 x 10-2 |
q=75 | 2.9388 x 10-2 | 2.6785 x 10-3 | 8.9444 x 10-3 |
q=90 | 3.3388 x 10-2 | - | - |