C 19-TH February, 1992 C C INPUT R* AND V INSTEAD MU AND BETA C TRY TO FIND THE CRITICAL GEOMETRIC RATIO FOR RUBBER C C EIGENVALUES SURFACE OF ANTISYMMETRIC DIFFUSE MODES C ************ FOR THE EI SUBREGIME *************** C C FOR INCOMPRESSIBLE ISOTROPIC PRESSURE-INSENSITIVE MATERIAL (RUBBER) C UNDER NO CONFINING STRESSES C C TRY TO COMPARE WITH BEATTY'S (1976) EXPERIMENTAL RESULTS C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C FOR PLOTTING ALL THOSE ELLIPTIC EIGENVALUES MODE C WITH CHANGING SZZ/(2.*GL) RATIO AND GAMMA (WAVELENGTH) C C DEVELOPED BY C K. T. CHAU C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C

C C HAVE TO INPUT FIVE MATERIAL PARAMETERS C AND ONE PARAMETER IS CHANGING WITH CURVES C C READ(5,10001)RG,RGK,PR,RR,SRE C READ(5,10002)REG C10001 FORMAT(5F10.0) C10002 FORMAT(F10.0) C C WE HAVE

C IN THIS PROGRAM C C DEFINE THE FOLLOWING RATIO : C C PR = POISSON RATIO C R = PRESSURE-SENSITIVE COEFFICIENT C RG = GT/GL C RGK = GL/K C SRE = SRR/E C SZE = SZZ/E (THIS THE UNKNOWN TO BE SOLVED) C RVV = R/2V (NON-NORMALITY PARAMETER) C WHERE C GL = LONGITUDINAL SHEAR MODULUS C GT = TRANSVERSE SHEAR MODULUS C E = INSTANTANEOUS TANGENT MODULUS C K = IN-PLANE BULK MODULUS C RR = PRESSURE-SENSITIVITY FACTOR C PR = EF C EFFCTIVE POISSON'S RATIO (V) C

C C KK1=1 FOR SYMMETRIC MODE C KK2=2 FOR ANTISYMMETRIC MODE C

C 20008 FORMAT(5X,' APPLIED COEFFICENT OF THE PROBLEM'/

C C

C C GAMMA = N(PI)A/L C WHERE C N = WAVE NUMBER C A = RADIUS OF THE CYLINDER C L = LENGTH OF THE CYLINDER C (GAMMA WILL BE THE X-AXIS OF THE PLOT) C

C C

2011 FORMAT(1X,' E/GL =',E20.6//) C

C C KK2 =1 FOR S<0 C KK2 =2 FOR S>0 C C

292 FORMAT(1X,' GAM SZZ/2GL',//)

C C ELLIPTIC COMPLEX Regime C

C C NB = NO. OF ROOTS WANTED C N= NO. DIVISION INPUT C

C C

12011 FORMAT(1X,'THERE IS NO ROOT (EC REGIME) AT GAM=',F20.5)

C XACC=1.E-6

C

776 CONTINUE 21111 FORMAT(2E20.6) C 221 CONTINUE C C C Elliptic-Imaginary C C C


C END OF THE EGL (VARYING) LOOP C


200 CONTINUE 100 CONTINUE 999 CONTINUE

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DEVELOPED BY C K. T. CHAU C REF: NUMERICAL RECIPES,PP.252 C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C Using Brent's method, find the root of a function FUNC knowm to c lie between X1 and X2. The root, returned as ZBREET, will be c refined until its accuracy is TOL. C

C

CC C

C C RENAME A, B AND C AND ADJUST BOUNDING INTERVAL D. C

CC C CONVERGENCE CHECK C

CC ATTEMPT INVERSE QUADRATIC INTERPOLATION C

C C CHECK WHETHER IN BOUNDS C

C C ACCEPT INTERPOLATION C

C C INTERPOLATION FAILED, USE BISECTION. C

C C BOUNDS DECREASES TOO SLOWLY, USE BISECTION. C

C C MOVE LAST BEST GUESS TO A C

C C EVALUATE NEW TRIAL ROOT C

C

C 11 CONTINUE

C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C SUBROUTINE TO BRACKET THE INTERVAL FROM X1 TO X2 C

C C GIVEN A FUNCTION FX DEFINED ON THE INTERVAL FROM X1 TO X2 C SUBDIVIDE THE INTERVAL INTO N EQUALLY SPACED SEGMENTS, AND C SEARCH FOR ZERO CROSSINGS OF THE FUNCTION. NB IS INPUT AS C THE MAXIMUM NUMBER OF ROOTS SOUGHT, AND IS RESET TO THE C NUMBER OF BRACKETING PAIRS XB1, XB2 THAT ARE FOUND. C C

C

C C LOOP OVER ALL INTERVALS C

C

C

C IF SIGN CHANGE OCCURS THEN RECORD VALUES FOR THE BOUNDS C

11 CONTINUE

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712 C

CC C NEW SUBROUTINE FOR ANTI-SYMMETRIC MODES C

C V3 FOR THE SECOND DIFFERENTIAL EQUATION

C C Elliptic-Complex Regime C C IF((VB*VB-4.*VA*VC).LT.0.0.AND.(VA*VC).GT.0.0) THEN C IF(VA.EQ.0.0) PAUSE 'A=0 IN FECA' C C CHECK POSITIVENESS OF P AND Q C C CALCULATION OF P AND Q

C

C

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

C

C

C C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

777 FEIA=U32*Y3*Y3*Y5+4.*Y2*Y3*Y6+U32*(Y7-3.*Y8)*(Y5*(4.*Y1+3.*Y3)

C ELSE 778 FEIA=WW1*(VQP*VRQ-VQQ*VRP)+WW2*(VPQ*VRP-VPP*VRQ)

999 RETURN

C


C EIGENVALUE EQUATION FOR SYMMETRIC GEOMETRIC DIFFUSE MODE C


C

C C Elliptic-Imaginary Subregime C

c PP=DSQRT((2.*DSQRT(VA*VC)-VB)/(4.*VA)) c QQ=DSQRT((2.*DSQRT(VA*VC)+VB)/(4.*VA))

C

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

C

C C

999 RETURN

C C SUBROUTINE FOR TIME CALCULATION FOR THE RUNNING TIME C **************************************************** C

$LARGE

300 CT = A

1000 FORMAT (41X, ' CURRENT TIME (SEC.)>>>>' ,F10.0,/)

100 PT = A

2000 FORMAT (41X, ' STARTING TIME (SEC.)>>>>',F10.0,/)

200 CT = A

3000 FORMAT (41X, ' ENDING TIME (SEC.)>>>>',F10.0,/)

C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

C


C---|


1


|


2


|


3


|


4


|


5


|


6


|


712

c C


C C FUNCTION RETURNS THE MODIFIED BESSEL FUNCTION OF ORDER N C FOR ANY REAL X AND N >= 2 C

11 CONTINUE