C Updated 19.50 17 December 1984 PMS1 C C********************************************************************| C | C This subroutine was originally written by Harold Frost, | C but has been completely rewritten, changed and updated | C by Philip Sargent. | C | C This version: #SIC4 11 Dec 1984 | C New version from #SIC3 in Fortran 77 | C | C********************************************************************| C | SUBROUTINE EQNS(G,GAMMA,JFIELD,T,TAU,ATVOL,BURGV,TMELT, 1 SHMODO,TCNORM,DOV,QVOL,WDOB,QBDRY,ACDOC,QCORE,QCREEP, 2 DGRAIN,A2,EN,GAMZRO,TOBST,DELFOB,TPNORM,GAMPLS,DELPLS, 3 TSTOP,BETA,PHONON,ELECTR,PHSCHG,PLSCUT,CHKNUM) C DIMENSION G(9) C REAL LARGE,SMALL C SAVE TLAST2 C COMMON /C1/ GZ(9),THETA C | C Eqns to match versions 186 and 166, 786, 986 and 1086. C | C BETA Power-law breakdown normalised stress. PHONON | C Phonon drag. ELECTR Electron drag (pure guess). PHSCHG | C alpha-Ti and alpha-Zr transformation temps. | C CHKNUM is a check number to ensure that the CONST values | C are correct for the particular version of subroutine | C EQNS in use. | C | C | C SIC1 PLSCUT is the normalised stress below which Peierls | C SIC1 stress darg control is disabled. | C | IF(CHKNUM.EQ.-900.0) GO TO 10 STOP 999 C 10 CONTINUE C---- IF (T.GT.200. .AND. T.LT.800.)WRITE(6,11) T,TAU 11 FORMAT(' EQNS:',F5.0,E10.3) C | C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| C calculate constants | C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| C | BOLTZ = 1.38054E-23 GASC = 8.31434E-03 IF(T.EQ.TLAST2) GO TO 20 C | C Calculate the shear modulus as a function of temperature. | C | SHMOD = SHMODO*(1. - (T - 300.)*(TCNORM/TMELT)) SHMD0K= SHMODO*(1 + 300*TCNORM/TMELT) C | C Special trap for alpha Ti which has TCNORM > 1.0 and | C therefore would have a negative modulus at its melting | C point. Actually, it transforms to beta at 1155K anyway. | C | C PHSCHG is the temperature (degrees Absolute) above which | C the map should not be used. This fudge is to enable us | C to run data for materials with extremely strong | C temperature sensitivit of their shear moduli, such that | C were they not to transform to another phase, they would | C have negative moduli at temperatures before they melted. | C | IF ( T.GT.PHSCHG) SHMOD = SHMODO* 1 (1. - (PHSCHG - 300.)*(TCNORM/TMELT)) C | C Calculate activation energies for obstacle-controlled | C and Peierls' stress controlled glide. | C | C | C Versions 1 to 22 of EQNS use: | C | C DELF = DELFOB*SHMOD*BURGV**3 C DEL = DELPLS * SHMODO * BURGV**3 C | C We use the following lines in Version 23 because | C although we should really use temperature dependent | C moduli, we should also really be using temperature | C dependent burgers vectors. Since we have not got the | C data for the expansion of the lattice with temperature, | C we partially compensate by omitting the temperature | C dependence of the modulus. FUDGE!!! | C | DELF = DELFOB*SHMD0K*BURGV**3 DEL = DELPLS*SHMD0K*BURGV**3 C C Calculate dvol, dbdry, and dcore. | C | RLGLIM = 139. SMALL = EXP(-RLGLIM) LARGE = 1.0E+30 QVOLKT = QVOL/(GASC*T) IF(QVOLKT.LT.RLGLIM) DVOL = DOV*EXP(-QVOLKT) IF(QVOLKT.GE.RLGLIM) DVOL = SMALL WDBE = LOG(WDOB) - QBDRY/(GASC*T) IF (WDBE.LT.RLGLIM) THEN WDBDRY = EXP(WDBE) ELSE WDBDRY = SMALL ENDIF ADE = LOG(ACDOC) - QCORE/(GASC*T) IF (ADE.LT.RLGLIM) THEN ADCORE = EXP(ADE) ELSE ADCORE = SMALL ENDIF QCRPKT = QCREEP/(GASC*T) IF(QCRPKT.LT.RLGLIM) DCREEP = EXP(-QCRPKT) IF(QCRPKT.GE.RLGLIM) DCREEP = SMALL C | C--------------------------------------------------------------------| C | TLAST2 = T 20 CONTINUE C | C--------------------------------------------------------------------| C | C TAUG is the (normalised) stress below which it is | C assumed that no glide can occur. This can lead to a | C discontinuity in strain-rate at this (normalised) | C stress. | C This is because TAUG is the 'athermal' strength below | C which no glide can occur at any temperature. | C | TAUG = TOBST * 0.5 C | C---- IF (T.GT.200. .AND. T.LT.800.)WRITE(6,21) DVOL,WDBDRY,ADCORE,DCREEP C---- 21 FORMAT('+',28X,4E10.3) C---- 22 FORMAT('+',82X,4E10.3) IF(T.EQ.TLAST) GO TO 50 C | C--------------------------------------------------------------------| C | C This next line calculates the Nabarro-Herring creep term | C without the stress term. | C | C IF(DVOL.GT.1.0E-20) GO TO 43 C X = SMALL C GO TO 44 43 X = 42.*DVOL*SHMOD*ATVOL/(BOLTZ*T*DGRAIN**2) C | C This next line calculates the Coble creep term, again | C without the stress term. It includes the Nabarro-Herring | C contribution and from no on these two mechanisms will be | C taken together. | C | 44 CONTINUE C1 = X*(1.+3.14159*WDBDRY/(DGRAIN*DVOL)) C | C This next line calculates the Power-Law creep term | C (again without the stress term) for the "lattice | C diffusion" contribution. Actually it uses DCREEP and not | C DVOL to allow for those cases where the activatione | C energy is not that for lattice diffusion (e.g. alpha | C Ti). | C This has the effect of replacing QVOL with QCREEP and | C DOV with 1.0, therefore the 'A' used here is different | C from the 'A' used in the Frost and Ashby book (1982). | C You have to multiply their 'A' by DOV to get the correct | C map output. | C | C2 = (DCREEP*SHMOD*BURGV/(BOLTZ*T))*A2 C | C This next line calculates the Power-Law creep term | C (again without the stress term) for the dislocation-core | C diffusion contribution. | C | C Remember, because we have a new A (and A2) we are also | C using a different ADCORE from before (ADCORE := ADCORE / | C DOV) from when we used volume diffusion as the rate | C determining step for HT power law creep. Now we use | C DCREEP as the rate determiner we use this new ADCORE if | C we want to get the same maps out. | C | C3 = (10./BURGV**2)*(ADCORE/DCREEP) C | C---- IF (T.GT.200. .AND. T.LT.800.)WRITE(6,22) X,C1,C2,C3 TLAST = T 50 CONTINUE C | C--------------------------------------------------------------------| C DIFFUSION CREEP | C--------------------------------------------------------------------| C | C A threshold stress is now introduced below which no | C diffusion controlled creep will be permitted. This means | C that a couple of safety features have to be included to | C prevent negative strain rates. The effective stress for | C diffusion creep will be the applied stress minus the | C threshold stress (Arzt and Ashby, 1982). BB is the | C grain-boundary dislocation burger's vector and is | C assumed always to be about Burger's-vector/3. C4 is the | C energy factor minus one. The energy factor is the | C proportional increase in the energy of the grain | C boundary dislocations as they climb from their minimum | C energy positions to their maximum energy positions. If | C only the increase in length is considered, then the | C energy factor must be between 1.08 (equiaxed grains, min | C to max perimeter values) and 1.7 (equiaxed grains, min | C perimeter to diameter value, certainly an overestimate). | C Here we have taken a factor of 1.3 and hence C4=0.3 . | C The grain boundary dislocations are the sources and | C sinks for the vacancies in the grain boundaries and must | C climb if these sources and sinks operate. | C | BBURGV = BURGV/3.0 C4 = 0.3 TR = C4 * BBURGV / DGRAIN TAU3 = TAU -TR IF(TAU3.LE.0.0) TAU3 = 1.0 E-30 C C Tuesday 11 December 1984 | C TSTOP is QVOL/(139 * GASC) and cuts off diffusion | C coefficient calculations below this temperature. Past | C versions jump from here to Obstacle control creep bit. | C | C C Thursday 6 December 1984 | C This line is likely to undeflow, so inserted print | C statement. | C | GDIFF = C1 * TAU3 C GZ(1) = GDIFF C C--------------------------------------------------------------------| C POWER LAW CREEP | C--------------------------------------------------------------------| C | C For shear stresses greater than about one | C -eight-hundredth of the shear modulus the power-law | C description of creep breaks down and a different | C empirical description must be used. A hyperbolic-sine | C description is identical to the power-law at low stress | C levels and thus can be used over both the conventional | C and the breakdown regions if a suitable change-over | C stress is chosen. High stresses in the breakdown regime | C are encountered in hot-working conditions and good data | C is obnly available for copper and aluminium. | C Nevertheless, the breakdown point (normalised stress TAU | C = BETA) is the best guess we have for the behaviour of | C any metal and the hyperbolic-sine description will be | C used for these high stresses wherever | C obstacle-controlled glide or Peierls'-stress controlled | C glide is not dominant. | C | PLBD = 1.0/BETA C | C We will take a threshold stress for power law creep to | C be the same as that for diffusion creep. If the gb | C dislocations cannot move, then the lattice dislocations | C cannot annihilate at grain boundaries and power law | C creep cannot occur. | C | ARG = PLBD * TAU3 C | C The next line protects against exponential overflow. It | C effectively limits the maximum normalised stress to | C TAU=139./(PLBD*EN) i.e. usually about 3.98E-02, it does | C this without warning. This contour and any at higher | C stresses will all appear coincident on the maps. | C | IF(ARG.LT.139./EN) GO TO 60 ARG = 139./EN 60 CONTINUE C IF(TAU.GT.TR) GO TO 70 GPOWLL = 1.0 E-20 GPOWLC = 1.0 E-20 GO TO 80 C 70 CONTINUE C | GPOWLL = C2*((SINH(ARG)/PLBD)**EN) C | GPOWLC = GPOWLL * C3 * TAU3**2 C 80 CONTINUE C GZ(2) = GPOWLL GZ(3) = GPOWLC C | C--------------------------------------------------------------------| C OBSTACLE CONTROL | C--------------------------------------------------------------------| C | C These next two lines check that the stress level is | C greater than half the theoretical strength of the | C material. | C If it is not, then there is no possibility of glide and | C the program does not bother to calculate GPLS or GOBST. | C The straining is wholly due to creep (power-law and | C diffusion). | C That, at least, is the logic behind this fudge. The | C actual reason is that the constant strain-rate contours | C (on the stress/temperature map) for power-law creep and | C for obstacle glide do not actually intersect at any | C temperature for a number of low strain-rate contours. | C This is because the deviation from the power-law at high | C stresses is actually much greater than is predicted by | C the formulae here, even the formula which pretends to | C describe power-law breakdown. The contours really must | C intersect as power-law breakdown is only an intermediate | C stage between power-law creep and simple | C obstacle-controlled glide. Since the power-law breakdown | C equation does actually predict quite high strain-rates | C we must check to see that these are below the limiting | C Peierls' stress limiting strain rate. | C Apart from all that, TAUG really does represent a | C 'completely athermal' resistance to glide in the | C material (well, it scales as the shear modulus anyway) | C which is revealed only at very low strain rates. | C | 90 IF(TAU3.GT.TAUG) GO TO 100 GOBST = 1.0E-50 GZ(4) = GOBST GO TO 110 C | C This next bit calculates the strain-rate caused by | C obstacle-glide. DELF is the activation energy and X is | C the term which is exponentiated, a check is made that | C this exponentiation will not go over the integer limit | C of the machine. GOBST is the strain-rate due to | C obstacle-glide. | C | 100 X2 = DELF*(1.-TAU3/TOBST)/(BOLTZ*T) IF(X2.GT.139.) X2 = 139. IF(X2.LT.-139.) X2 = -139. GOBST = GAMZRO * EXP(-X2) GZ(4) = GOBST C | C--------------------------------------------------------------------| C PEIERLS' STRESS DRAG CONTROL | C--------------------------------------------------------------------| C | C A check is now made on whether to continue to use the | C obstacle-glide strain rate or whether to calculate what | C the exponential term for the Peierls' mechansism would | C be, in any case, the Peierls' strain-rate is calculated. | C | 110 X3 = TAU3*SHMOD/(TPNORM*SHMODO) IF(X3.GE.1.0) GO TO 130 C | C This next bit calculates the strain-rate due only to the | C Peierls' stress resistance to glide. DEL is the | C activation energy. Y is the term to be exponentiated, a | C check is made to see whether it is too large for the | C machine to cope with. | C GPLS is the strain-rate due to this Peierls' mechanism. | C | Y = -DEL*(1. - X3**0.75)**1.33333/(BOLTZ*T) C Thursday 6 December 1984 | C This line may possibly undeflow. | C | IF(Y.GT.-140) GPLS=EXP(Y)*GAMPLS*TAU3*TAU3 IF(Y.LE.-140) GPLS=SMALL GO TO 140 C | 130 GPLS = GAMPLS * TAU3*TAU3 140 CONTINUE IF(TAU3.LT.PLSCUT) GPLS = SMALL GZ(5) = GPLS C | C--------------------------------------------------------------------| C PHONON DRAG CONTROL | C--------------------------------------------------------------------| C | C This next bit calculates the strain rate which is | C limited by non-relativistic phonon drag. We know that | C B/rho (B divided by the mobile dislocation density) is | C between 1E-16 and 1E-15 Newton seconds for various forms | C of pure aluminium (Kumar et al. 1968). We will estimate | C the mobile dislocation density as 1.E+11 m(-2). | C | C Actually, we will take a modulus-corrected 'B-type' | C quantity incorporating the shear modulus at 300 K, the | C burger's vector (squared), B and the mobile dislocation | C density. This we will call a normalised Viscous Mobility | C Coefficient (VMC). We will take it to have a value of | C 5.2 E +06 per second (derived from aluminium data). This | C constant needs only a single lg(strain-rate)/stress data | C point in the viscous regime to be defined. We do not | C actually need to know B and rho separately. VMC varies | C quite widely over at least an order of magnitude for the | C data we have for just aluminium and copper alone. | C VMC = rho * shmod(300K) * b**2 / B | C | C It is observed in Titanium that linear-viscous | C dependence of flow-stress on strain-rate is observed at | C strain rates as low as 1E-3 at room temperature. This | C contradicts the maps that we draw up using phonon-drag | C data for Aluminium and Copper. Apart from fiddling with | C TOBST in an attempt to resolve this, we have introduced | C a temperature dependence of the phonon drag and also a | C guess at a constant electronic contribution to prevent | C there being no drag at cryogenic temperatures. | C | VMC = 1.0/(PHONON*T + ELECTR) GNRPD = VMC * TAU3 GZ(6) = GNRPD C | C We also need a relativistic-type relationship to ensure | C that the dislocation velocities only approach the speed | C of sound asymtotically. Since we do not know anything | C about dislocation densities at such high stresses, we | C will pick a limiting strain rate, rather than a limiting | C velocity. 1.0 E +6 per second is a good limit. | C | GLIMIT = 1.0 E+06 GLHALF = GLIMIT/2.0 C GRPD = GNRPD/SQRT(1.+ (GNRPD/GLIMIT)**2) GZ(7) = GRPD GRPDG = VMC*TAUG/SQRT(1.+ (VMC*TAUG/GLIMIT)**2) C | C--------------------------------------------------------------------| C FINAL STRAIN RATE DECISIONS | C--------------------------------------------------------------------| C | C This next line combines the strain-rates due to the | C lattice and core diffusion contributions to the | C Power-Law creep strain-rate. | C | GPOWL = GPOWLL + GPOWLC C | C This next bit sorts out how to organise the calculation | C of the overall strain-rate. | C | C #22 There are only a fixed number of dislocations and they | C #22 can either glide past obstacles by thermal activation | C #22 (obstacle-glide) or they can climb past them (power-law | C #22 creep). They will follow whichever process they can do | C #22 fastest. However, there is an overall limitation that | C #22 neither of these processes can proceed at a rate faster | C #22 than that controlled by the Peierls' stress resistance | C #22 to dislocation movement. Therefore, if the Peierls' | C #22 stress controlled strain rate (GPLS) is lower than | C #22 either the obstacle-glide rate (GOBST) or the power-law | C #22 creep rate (GPOWL) then GPLS will control the overall | C #22 rate. | C | GO TO 210 C | C #22 Here ther is no glide (either TAU > TAUG or power law | C #22 creep is faster than obstacle controlled glide or both). | C | 160 CONTINUE GAMMA = GDIFF + GPOWL IF(GDIFF.LT.GPOWL) GO TO 180 170 JFIELD = 1 RETURN C | C Here the diffusional flow is faster than power law | C creep. | C | 180 CONTINUE IF(GPOWLC.GT.GPOWLL) GO TO 190 JFIELD = 2 RETURN C | C This next bit covers the case where obstacle-controlled | C glide is faster than Power-law creep. | C | 190 CONTINUE JFIELD = 3 RETURN C | C SIC1 Here we take the HIGHEST strain rate out of the | C SIC1 drag-control and the diffusion/power law creep control. | C SIC1 This is because we assume that as soon as glide occurs | C SIC1 it becomes Peierls-limited. | C | 210 CONTINUE GAMMA = GDIFF + GPOWL GB = AMIN1(GRPD,GPLS) IF(GB.GT.GAMMA) GO TO 215 C | C This next bit is only in SIC3 and it applies phonon drag | C restrictions to the power law creep field. | C | IF(GRPD.LE.GAMMA) GO TO 230 GO TO 160 C 215 IF(GRPD.LT.GPLS) GO TO 230 GAMMA = GPLS 220 JFIELD = 5 RETURN C 230 CONTINUE GAMMA = GRPD 240 IF(GAMMA.GT.GLHALF) GO TO 250 JFIELD = 6 RETURN C 250 CONTINUE JFIELD = 7 RETURN C | C JFIELD is set to the number flag for the dominant | C mechanism: (1) Nabarro-Herring/Coble, (2) | C Power-law (lattice diffusion), (3) Power-law (core | C diffusion) (4) Obstacle-controlled glide, (5) | C Peierls' stress controlled glide, (6) Non-Relativistic | C Phonon Drag or (7) Relativistic Phonon Drag. | C | END