NASA Technical Memorandum 102338 

Refinements in a Viscoplastic Model 

A.D. Freed 

Lewis Research Center 

Cleveland, Ohio 

and 

K. P. Walker 

Engineering Science Sofhyare 

SmithJeld, Rhode Island 

Prepared for the 

1989 Winter Annual Meeting 

of the American Society of Mechanical Engineers 

San Francisco, California, December 10-15, 1989 

(NASA-T&-102338) REPINEHENTS IB A rS 89-28036 VISCOPLASTIC HODEL (NASA. Levis Research 

enter) 18 p CSCL 20K 

Unclas 

G3/39 0225965 
REFINEMENTS IN A VISCOPLASTIC MODEL 

A. D. Freed 

National Aeronautics and Space Adminstration 

Lewis Research Center 

Cleveland, OH 44135 

K. P. Walker 

Engineering Science Software, Inc. 

Smithfield, RI 02917 

ABSTRACT 

A thermodynamically admissible theory of viscoplasticity with two internal variables (a 

back stress and a drag strength) is presented. Six material functions characterize a specific 

viscoplastic model. In the pursuit of compromise between accuracy and simplicity, a model 

is developed that is a hybrid of two existing viscoplastic models. A limited number of applications of the model to Al, Cu, and Ni are presented. A novel implicit integration method is 

also discussed, and applications are made to obtain solutions using this viscoplastic model. 

1. INTRODUCTION 

The design of structures exposed to high-temperature environments requires mathematical models capable 

of accurately predicting short-term plastic strains, long-term creep strains, and the interactions between them. 

Multiaxial, cyclic, and nonisothermal loading histories are normal service conditions, not exceptional ones. In 

the pursuit to achieve this objective, a hybrid model is developed using WALKER’S [l] theoretical structure for 

the evolution of internal state and MILLER’S [l] methodology for characterizing material functions. The main reason for doing inelastic analysis is not so much to determine the deformation response of the structure per se, but 

rather to assess its useful service life. Results from inelastic analysis are input to life assessing schemes. These 
schemes typically require such quantities as: stress range, mean stress, inelastic strain range, and the percentages 

of time-dependent creep and time-independent plasticity that make up the inelastic strain range [3]. It is therefore important that the mathematical model used in high-temperature structural analysis be capable of predicting 

these parameters with reasonable accuracy. 

A wide range of materials are used in structures that are designed to operate in high-temperature environments. The designer's choice of material is usually dictated by cost and by such factors as its ability to: carry 

load, conduct heat, resist crack initiation and growth, and resist oxidation and corrosion. A typical metal may 

have some desirable properties, and some undesirable ones as well. With the technology of metal-matrix fabrication coming of age, the engineer can now begin to tailor materials that provide optimal characteristics for a 

given application. At NASA, composite matrices of aluminum, copper, nickel, and their alloys are being considered for a variety of applications. Along with ceramics, tungsten and tungsten alloys are being considered as 

fiber materials. An objective of this paper is to develop a model for characterizing their viscoplastic behavior. 

These material representations will be used in the near future in a mesomechanics model [4,5] to predict the 

viscoplastic response of metal matrix composites. 

The need for reliable and expedient algorithms to numerically integrate viscoplastic models is of 

paramount importance. Such models are composed of systems of stiff, first-order, linear, differential equations 

whose coefficients are coupled and strongly nonlinear; hence, they are inherently difficult to integrate. The 

implicit integration algorithm of WALKER [6,7], which is presented in appendix A and extended there to include 

the second-order term in its difference formulation, is applied to our viscoplastic model. 

2. VISCOPLASTIC THEORY 

The structure of our viscoplastic theory is a special case of a much more general structure that FREED & 

CHABOCHE [8] have shown to be physically motivated and thermodynamically admissible. The theory admits two 

fundamental internal variables: they are: the (positive scalar-valued) drag suength D > 0, and the (deviatoric 

tensor-valued) back stress Bi,. The drag strength accounts for isotropic hardening effects. The back stress 

accounts for kinematic (or flow-induced anisotropic) hardening effects, and is considered to be the sum of two 

separately evolving constituents. One constituent Bij represents short-range back stress effects, while the other 

constituent Bi: represents long-range back stress effects [9,10]. These internal variables are considered to 

evolve phenomenologically through competitive processes associated with strain hardening, strain-induced 

dynamic recovery, and time-induced thermal recovery. 

The strain &ij is taken to be the sum of an elastic (or thermodynamically reversible) strain E$ and an 

inelastic (or thermodynamically irreversible) strain E;, i.e., 

(1) 

with no inelastic strain occurring in the stress-free virgin state. Small strains, displacements, and rotations are 

considered to make up the deformation of a material element. 

The constitutive equations of the theory are given by: 

E.. = &e. + &P IJ IJ lJ 

D =Do+h5 (2) 

where E is the elastic modulus, v is the Poisson ratio, a is the coefficient of thermal expansion, T-To is 

the change in temperature, qj is the stress, and 6, is the Kronecker delta. Here p$ and 6 2 0 are the thermodynamic displacements conjugate to the thermodynamic forces Bi! and D-Do 2 0, with H, > 0 and 

h > 0 denoting their hardening coefficients, and Do > 0 defining the annealed value of drag strength. Repeated 

Latin indices are summed over from 1 to 3 in the usual manner. 

The evolution equations of the theory are given by: 

2 
which, when combined with the constitutive equations (2b and 2c), give rise to the governing viscoplastic equations of our theory; they are: 

2 

3 Bij = -(H1+H2) &$ - 

d=h[[l-f]Wv-8.] +-I; D-Do 

D-Do h 

where 

with directions (unit or projection) defined by: 

and with norms (or magnitudes) defined by: 

I&Pl=~~ , Vl=dW 

where Jij€ (Sij,Bij,Sij-Bij,Xij,Sij-Xij), Sij = oij - ~~6~~13 is the deviatoric stress, and Sij-Bij is the viscous 

stress that drives inelastic flow. This choice for the directions and norms scales the theory for tension. The 

notation Xij=Bij. is used for easier visual distinction between Xij and Bii. The derivatives of the hardening 

coefficients h , H1 and H2 represent their change with temperature, if any. 

The material functions of the theory are: 9(T) > 0 is the thermal diffusivity, Z(IS -B 110) = l/e 2 0 is 

the Zener [ll] parameter, L > 0 and 1 > 0 are the limiting states associated with dynamic recovery, and 

R 2 0 and r 1 0 are the thermal recovery functions. At this point in the development, the limit and thermalrecovery functions are general functions of state. Specific forms for the functions 0,Z, I, L , r , and R will 

vary from model to model. There is one material constant in these evolution equations; it is the nonproportionality parameter a E [O,l] found in eqn. (5d). This parameter provides added flexibility in modeling nonproportional kinematic behavior; it has no effect on proportional behavior. What it does (more or less) is proportion 

the nonlinear strain-induced characteristics of the back stress between the mechanisms of hardening and recovery 

[81. The value of this parameter has a profound impact on predicted response under nonproportional ratchetting 

conditions [12]. 

Notice that the drag strength D evolves with W’l(D-Do). It Cannot evolve with either IEPpl or Sij&$ 

(the choices in existing models), and satisfy the thermodynamic constraint of intrinsic dissipation throughout the 

entire domain of state space [8] defined by the set (T, oij, D , Bij Therefore, in the absence of recovery, the 

was introduced phenomenologically by MILLER & SHERBY [131 into their model; here it is introduced as a direct 

consequence of the thermodynamics. Because D-Do is in the denominator of this forcing function, it appears 

to be singular in the virgin state, D =Do. But actually the forcing function is zero valued there, because the 

numerator, ~’, is zero valued in the virgin state, too. 

drag strength parabolically hardens with viscous work, Le., D Q. & W’. Parabolic hardening of the drag strength 

3 
The basic form of the evolution equations (3b and 3c) for the kinematic variables and bz was proposed by WALKER [l], and is a compromise between accuracy and simplicity. By assuming that the long-range 

back stress Ei; does not thermally recover, which is physically incorrect (a compromise), it can be integrated 

exactly and combined with the short-range back stress Ei) = Xi,. The result is a single governing equation (4b) 

for the back stress (instead of two) which has both short- and long-range features providing improved accuracy. 

The direction for dynamic recovery di, of the back stress, as established in eqn. (5d), presents a refinement to 

the original theory. Dynamic recovery is taken to be proportioned by a factor a between directions opposing 

the short-range back stress Xi, and the stress difference Si, -Xi,. The scaling factor uiiXi,/w I in eqn. (Sd) is 

introduced to make the model insensitive to the direction of dynamic recovery under proportional loading conditions @e., the predicted response is the Same for all values of the parameter a E [OJ]). The model only becomes 

sensitive to the direction of dynamic recovery (Le., the value of a) under nonproportional loading conditions. 

This scaling factor is also required for the theory to be thermodynamically admissible [8]. 

3. VISCOPLASTIC MODEL 

Our model is predicated upon the assumption that the magnitude of inelastic strain rate (in particular, Z,,) 

can be represented as a sum of two power functions in stress (cf- appendix B). One to govern at the highertemperatures/lower-stresses, while the other to govern at the lower-temperaweshigher-stresses. An analogous 

expression for Z has recently been proposed by CHABOCHE [I41 to model the transient strain-rate sensitivity of 

a stainless steel. This parameter is a generalization of WALKER'S [l] power-law model, and a simplification 

(especially from a numerical integration viewpoint) of MILLER'S [2] hyperbolic sine model, i.e., a compromise 

between accuracy and simplicity. We motivate its choice with the steady-state creep data of fig. 1, where 

2, = I&"VB. By applying Miller's hypothesis (cJ appendix B) to our relationship for the stress dependence of 

inelastic flow at steady state, we are able to obtain explicit relationships for the transient Zener parameter Z, 

and both of the thermal recovery functions r and R. Details on the derivation of the material functions, and a 

procedure for determining values for their associated constants, are presented in appendix B. 

The material functions of OUT viscoplastic model are: 

L= CD (W 

Hl+H2 R= --- [ H, VTL] [ [=In + [fq] 

where eqn. (6e) implies no dynamic recovery of the drag strength. (A particular limit function 1 is under consideration to account for monotonic/cyclic interactions, which are not incorporated in the current model, i.e., 

1 = -.) The material constants for this model are comprised of: the elastic constants: E, a, and v; the thermal 

activation energy and transition temperature: Q and T,; the steady-state constants: A, A', n , and n'; the 

constant-structure constant: m; the threshold constants: b and c; the transient constants: h, HI, and H2; the 

nonproportional constant: a; and the initial values: Do and D(t=O). The parameter k is the Boltzmann constant, and H(D-Do) is the Heaviside function. Values for these constants are given in table 1 for aluminum, 

copper, nickel, and tungsten. If one chooses, a reference temperature could be introduced (thereby setting 

B(Tref)=l) which would, in effect, scale the drag strength. 

4 
This model has several deficiencies that the authors are aware of, and are working on. It does not account 

for monotonic/cyclic interaction effects. It overpredicts normal ratchetting behavior. And it does not account for 

inverse strain rate sensitivity effects exhibited by many alloys. These deficiencies are common amongst most 

viscoplastic models found in the literature. 

4. NUMERIC INTEGRATION 

and obtain the following system of stiff, linear, first-order, difference equations: 

One can combine the constitutive equations (2) and the evolution equations (3) of our viscoplastic theory 

AK + KAc = Az (7c) 

with 

where the coefficients are nonlinear and coupled, and are defined by: 

and where the nonhomogeneous terms are also nonlinear and coupled, and are defined by: 

Bij 

CL Axij = 2p Aeij - ABij + -Ap 

AZ = h (S-8 I 111.~~1 (W 

Here Ap, Ah, and AH1 represent the incremental change in these moduli due to an incremental change AT in 

temperature. This system of difference equations characterizes the deviatoric response of the material, while the 

constitutive equation a, = ~K[E~~ - a(T-To)Sii] characterizes the hydrostatic response of the material. Here 

p = E/2(l+v) is the elastic shear modulus, and K = E/3(1-2v) is the elastic bulk modulus. 

The system of difference equations (7) can be solved by classical methods, such as the explicit methods of 

Euler and Runge-Kutta. Such methods require numerous time steps in order to maintain numeric stability and 

accuracy because the coefficients and nonhomogeneous terms make this system of equations mathematically 

stiff. This is not a detriment when numemus data pairs are needed to plot smooth accurate response curves. 

However a new method, the uniformly-valid implicit integration scheme of WALKER [6], may be preferred for 

finite element analyses where computational times are large, and the ability to create smooth response curves at 

the nodal points is not the typical required output [7]. An advantage of this method is that it is not constrained 

by time step size in order to maintain numeric stability and accuracy since the scheme is uniformly valid. 

Walker's method transforms the system of difference equations (7) into recursive integral equations, which are 

then solved using an asymptotic expansion technique (cf appendix A). The resulting integration scheme for the 

5 
general viscoplastic theory presented in eqn. (4) is given by the following system of stiff, coupled, implicit, 

recursive relations: 

Ayi;(t+At) [ 1 - e-Ab(t+At) + 1 Xi, (t +At) = Xi, (t ) e”b(t+At) + Ab (t +At ) 

where t 2 0 and where Af > 0 may be arbitrarily large or small. These equations must be solved simultaneously using a technique such as Newton-Raphson iteration, whe;e Aa, Ab, and Ac are the iterated variables 

[6,7]. For viscoplastic models that have a back stress and a drag strength (or yield strength) which evolve, 

Walker’s implicit integration scheme has only three variables to iterate over; whereas, the classical implicit 

integration schemes (e.g. backward Euler difference) have thirteen variables to iterate over (Le., six for inelastic 

strain, six for back stress, and one for drag strength). 

5. APPLICATIONS 

A few applications for aluminum, copper, and nickel (data are insufficient to include tungsten) are 

presented to demonstrate some of the capabilities and limitations of the model, and to show feasibility of the 

Walker integration method. 

The model’s capability in predicting the monotonic behavior of aluminum as a function of temperature is 

demonstrated in fig. 2. The model produces acceptable predictions above approximately .2 Th (or about - 

100°C, where T,, is the material’s homologous temperature). This is an improvement over a single power-law 

for the Zener parameter with an exponent of n-5, which produces acceptable predictions above approximately 

.5 Th (or about 200OC for aluminum). In applications requiring modeling at both elevated temperatures and 

temperatures below .2 Th , the hyperbolic Zener parameter of MILLER [2] is recommended. Numeric integration 

using the technique of Walker described in the previous section and the second-order method of Runge-Kutta 

with automatic time-step control was done here with good agreement between the two methods (less than 1 MPa 

difference, with a maximum error of 20% which occurred at the highest temperatures where the stress is on the 

order of 5 MPa in fig. 2). The method of Walker was found to be superior to that of Runge-Kutta in both accuracy and computation time. This was especially true when the model became mathematically stiff (Le., in the 

regions where hardening and recovery mechanisms are of comparable magnitude). Here the difference in cpu 

times exceeded an order in magnitude, and the Runge-Kutta solutions oscillated. 

A thermornechanical prediction of the model (in this case, for nickel) is presented in fig. 3a. The model 

produces the correct shape of the hysteresis loop with two exceptions. Of course, it cannot model the serrated 

yielding (or strain avalanching, [28]) observed in the experimental response. And second, the model presently 

has no provision to account for monotonic/c yclic interactions, and therefore it continues to isotropically harden 

beyond that which is experimentally observed. This deficiency is being worked on. 

Figure 3b presents a comparison between predictions obtained using the methods of Walker (with several 

different time steps) and second-order Runge-Kutta (with automatic time stepping). The Walker analysis with 

the largest time step size (50°C / step) ran about 2.5 times faster than the Runge-Kutta analysis, the one with the 

intermediate time step (approximately 15°C / step) ran about 1.5 times faster, while the one with the smallest 

time step (1.5OC / step) ran about 3 times slower. The greatest variation in the three Walker solutions occurs at 

the knees of the hysteresis curve. A benefit of the Walker (implicit integral) method is that this error is not 

6 
propagated, as all three solutions agree at the turnaround points. This is not true of explicit differential methods, 

such as Runge-Kutta, where an error initiated is also propagated. The effect of a slowly accumulated propagating error is found in the Runge-Kutta solution, and can be observed by looking at successive peaks in fig. 3b. 

The model’s capability in predicting nonproportional material behavior (for copper) is demonstrated in fig. 

4. Here the influence of the nonproportionality parameter a on predicted response is shown. The limiting 

case of a = 0 represents the nonlinear kinematic hardening rule of Armstrong and Frederick [9,14], while the 

other extreme of a = 1 denotes a nonlinear Prager-like hardening rule [8,12]. The prediction for a = 1 produces an unwanted ratchetting down the axial stress axis; a result also observed by LAMBA & SIDEBOTrOM [29] 

in plasticity analyses using the Prager and Ziegler hardening rules. There is little difference in the predicted 

response over a wide range in the parameter a (a = 0 to a = 0.9 in this experiment). Therefore a more 

severe nonproportional loading history (such as advocated in ref. 12) is necessary to characterize the value of the 

material constant a. This figure illustrates another limitation of the model. Our model, like most, is formulated with a von Mises type norm, whereas copper behaves like a Tresca material [29]. That is why the 

predicted axial stress response is too soft, as the value of the drag strength was correlated with the shear 

response in this figure. 

6. CONCLUSIONS 

Within the framework of a thermodynamically admissible theory of viscoplasticity with an evolving back 

stress and drag strength, a model is presented that is a hybrid of the WALKER [ll and MILLER [2] models. The 

model is for pure metals, and has been applied to aluminum, copper, and nickel. Tungsten has been partially 

characterized, using rules of thumb instead of experimental data (they do not exist) to establish values for 

roughly half of the constants. Both strengths and weaknesses of the model are discussed. The implicit integration scheme of WALKER [6,7] has been further developed and applied to our model, and was found to be superior 

to the second-order method of Runge-Kutta. A detailed numerical study of the Walker integration technique is 

currently underway. 

REFERENCES 

Dl 

r21 

r31 

r41 

PI 

[61 

r71 

181 

r91 

r 101 

Walker, K. P, 1981, “Research and Development Program for Nonlinear Structural Modeling With 

Advanced Time-Temperature Dependent Constitutive Relationships,” NASA CR-165533. 

Miller, A., 1976, “An Inelastic Constitutive Model for Monotonic, Cyclic, and Creep Deformation: Part I 

- Equations Development and Analytical Procedures,” J. Eng. Mater. Technol., Vol. 96, p. 97. 

Saltsman, J. F., and Halford, G. R., 1988, “Life Prediction of Thermomechanical Fatigue Using Total 

Strain Version of Strainrange Phoning (SRP): A Proposal,” NASA TP-2779. 

Walker, K. P., Jordan, E. H., and Freed, A. D., 1989, “Nonlinear Mesomechanics of Composites With 

Periodic Microstructure: First Report,” NASA TM-102051. 

Walker, K. P., Jordan, E. H.. and Freed, A. D., 1989. “Equivalence of Green’s Function and the Fourier 

Series Representation of Composites With Periodic Microstructure,” Micromechanics and Inhomogeneity - 

The Toshio Mura Anniversary Volume, Springer-Verlag, p. 535. 

Walker, K. P., 1987, “A Uniformly Valid Asymptotic Integration Algorithm for Unified Viscoplastic Constitutive Models,” Advances in Inelastic Analysis, s. Nakazawa et al., eds., AMD - Vol. 88 I PED - Vol. 

28, ASME, p. 13. 

Chulya, A., and Walker, K. P., 1989, “A New Uniformly Valid Asymptotic Integration Algorithm for 

Elasto-Plastic-Creep and Unified Viscoplastic Theories Including Continuum Damage,” to appear: NASA 

TM. 

Freed, A. D., and Chaboche, J. L., “Viscoplasticity: A Thermodynamic Formulation,” submitted to: J. 

Appl. Mech. 

Chaboche, J. L., and Rousselier, G., 1983, “On the Plastic and Viscoplastic Constitutive Equations - Part 

I: Rules Developed With Internal Variable Concept,” J. Pressure Vessel Technol., Vol. 105, p. 153. 

Lowe, T. C., and Miller, A. K., 1986, “Modeling Internal Stresses in the Nonelastic Deformation of 

Metals,” J. Eng. Mater. Technol., Vol. 108, p. 365. 

7 
Zener, C., and Hollomon. J. H., 1944, “Effect of Strain Rate Upon Plastic Flow of Steel,” J. Appl. Phys., 

Vol. 15, p. 22. 

Burlet, H., and Cailletaud, G., 1987, “Modeling of Cyclic Plasticity in Finite Element Codes,” Constitutive Laws for Engineering Materials: Theory and Applications, C. S. Desai et al., eds., Elsevier Applied 

Science, New York, p. 1157. 

Miller, A. K., and Sherby, 0. D., 1978, “A Simplified Phenomenological Model for Non-Elastic Deformation: Predictions of Pure Aluminum Behavior and Incorporation of Solute Strengthening Effects,” Acta 

Metall. Vol. 26, p. 289. 

Chaboche, J. L., 1989, “Constitutive Equations for Cyclic Plasticity and Cyclic Viscoplasticity,” to 

appear: Int. J. Plast.. 

Weertman, J., 1956, “Creep of Polycrystalline Aluminum as Determined from Strain Rate Tests,” J. 

Mech. Phys. Solids, Vol. 4, p. 230. 

Servi, I. S., and Grant, N. J., 1951, “Creep and Smss Rupture Behavior of Aluminum as a Function of 

Purity,” J. Metals, Vol. 3, p. 909. 

Luthy, H., Miller, A. K., and Sherby, 0. D., 1980, “The Stress and Temperature Dependence of SteadyState Flow at Intermediate Temperatures for Pure Polycrystalline Aluminum,” Acta Metall., Vol. 28, p. 

169. 

Jenkins, W. D., and Digges, T. G., 1951, “Creep of Annealed and Cold-Drawn High-Purity Copper,” J. 

Res. NBS, Vol. 47, p. 272. 

Barren, C. R., and Sherby, 0. D., 1964, “Steady-State Creep Characteristics of Polycrystalline Copper in 

the Temperature Range 400” to 950”C,” Trans. AIME, Vol. 230, p. 1322. 

Pahutovl, M., cadek, J., and Ryi, P., 1971, “High Temperature Creep in Copper,” Philos. Mag., Vol. 23, 

p. 509. 

Jenkins, W. D., Digges, T. G., and Johnson, C. R., 1954, “Creep of High-Purity Nickel,” J. Res. NBS, 

Vol. 53, p. 329. 

Weertman, J., and Shahinian, P., 1956, “Creep of Polycrystalline Nickel,” J. Metals, Vol. 8, p. 1223. 

King, G. W., and Sell, H. G., 1965, “The Effect of Thoria on the Elevated-Temperature Tensile Properties 

of Recrystallized High-Purity Tungsten,” Trans. AIME, Vol. 233, p. 1104. 

Klopp, W. D., Witzke, W. R., and Raffo, P. L., 1965, “Effect of Grain Size on Tensile and Creep Properties of Arc-Melted and Electron-Beam-Melted Tungsten at 2250” to 4140°F;’ Trans. AIME, Vol. 233, p. 

1860. 

Flagella, P. N., 1967, “High Temperature Creep-Rupture Behavior of Unalloyed Tungsten,” Proceedings 

of the Third International Symposium on High Temperature Technology, Butterworths, p. 279. 

Henshall, G. A., 1987, “Solute-Enhanced Back Stresses and Their Role in a Simplified Phenomenological 

Constitutive Model for the Nonelastic Deformation of Metals and Alloys,” PhD. Thesis, Stanford University. 

Verrilli, M. J., 1988, unpublished experimental data, NASA. 

Yan, B.-D., and Laird, C., 1985, “The Phenomenon of Strain Avalanches in Cyclic Deformation,” Acta 

Metall., Vol. 33, p. 2023. 

Lamba, H. S., and Sidebottom, 0. M., 1978, “Cyclic Plasticity for Nonproportional Paths: Part 2 - Comparison With Predictions of Three Incremental Plasticity Models,” J. Eng. Mater. Technol., Vol. 100, p. 

104. 

Il’yushin, A. A., 1954, “On the Relations Between Stresses and Small Deformations in Mechanics of Continuous Media,” PPM J. Appl. Math. Mech., Vol. 18, p. 641. 

Sherby, 0. D., and Miller, A. K., 1979, “Combining Phenomenology and Physics in Describing the High 

Temperature Mechanical Behavior of Crystalline Solids,” J. Eng. Mater. Technol., Vol. 101, p. 387. 

Walker, K. P., and Jordan, E. H., 1985, “Biaxial Constitutive Modelling and Testing of a Single Crystal 

Superalloy at Elevated Temperatures,” presented at: Biaxial and Multiaxial Fatigue Conference, 

8 
Sheffield, England. Published in: Biaxial and Multiaxial Fatigue, EGF3, M. W. Brown and K. J. Miller, 

eds., Mechanical Engineering Publications, 1989, p. 145. 

Marquis, D., 1979, "Modklisation et Identification de L'Ccrouissage Anisotrope des MCtaux," PhD. 

Thesis, Universitk Pierre et Marie Curie, Paris. 

APPENDIX A 

WALKER IMPLICIT INTEGRATION METHOD 

A recursive integration relation for a linear. first-order, differential equation is given. The uniformly-valid 

asymptotic expansion of WALKER [61 is applied to this integral equation. In the past, this method has been 

applied to specific viscoplastic models retaining only the first term in the expansion [6,7]. Here we apply it to a 

general differential equation so that its structure is more fully revealed, and it is extended to include the second 

term in the expansion, 

Consider the linear, first-order, differential equation 

X+UX=3 (All 

where li and + are the parameters, which may be nonlinear functions of X, and X is the variable to be 

solved for. One can integrate this differential equation over the time interval [O,t +At 3 and obtain 

where Xo = X(0). The integral in this equation is the nonhomogeneous contribution to the solution of the 

differential equation (Al). This integral equation can be cast into a recursive relation by dividing the interval of 

integration into two parts, viz., 

Substituting the identity 1 = e-'(l)eu(l) into the first integral of this equation results in 

which simplifies to the desired recursive integral equation 

where t 2 0 and u (0) = X (0) = 0. The derivation of this relation made use of eqn. (A2) when integrated over 

the time interval [O,t]. 

The recursive relation given by eqn. (AS) for the differential relation in eqn. (Al) becomes practical only 

when a solution strategy exists for evaluating the integral 

When u(t) is a monotonically increasing function of f, this integral possesses an evanescent memory of the 

forcing function av(t)/ag, causing most of the contribution to the integral to come from the region near the 

9 
upper limit 4 = t+At , where the exponential term is approximately equal to its maximum value of 1. Using the 

method of Laplace, WALKER [6] derived a uniformly-valid asymptotic expansion for this integral equation by 

first expanding the arguments of the exponential and derivative into a Taylor series, and then integrating. In particular, by writing u(4) = u[t+At - (t+At-t)] and expanding it by Taylor's theorem to give 

du (t+At) + . . . , the integral equation (A6) becomes 

d(t+At) ~(4) = u(t+At) - (t+&<) 

which can be rewritten in equivalent form as 

tj introducing the change in variable z = t+At-€,. If one now expands the argument of the derivative in this 

integrand using Taylor's theorem to give v (t +At -z ) = v (t +At ) - z dv(t+At) + . . . , then this derivative can 

be represented by the series = - i + zV - . . . , and the integral in eqn. (A8) 

d(t+&) aV(t+At-z) - - av(t+At-z) - aZ act +& 

can be expressed as 

AI 

I(A~) = j e-(&- '.. )(i - ZV + . . . )dz 

Zd 

Integrating this relationship gives a uniformly valid asymptotic expansion for the integral in eqn. (A@; viz., 

where the third and higher-order terms in the expansion have been neglected. Here we recall that the derivatives of u and v are evaluated at time t+At. 

A uniformly-valid, implicit, numeric, integration method providing a solution strategy to the recursive 

integral equation (A5) for the linear, first-order, differential equation (Al) is therefore given by the difference 

relation 

where bx(t+At)=x(t+At)-x(r) =i(t+At)At, Ai(t+At) =x(t+At)-i(t) =i(t+At)At, and hence 

hK(t+At) - &(t) = [x(t+&) - x(t)]Af = i(t+At)(At)2 for a constant time step size At. This implicit 

integration scheme must be solved iteratively, since Au(t+At) and Av(t+At) may be nonlinear functions of 

X (t +At ). 

The implicit integration method of eqn. (All) is valid for an arbitrary time step size At > 0, be it large 

or small, even when only the first term in the expansion is retained, i.e., it is uniformly valid. This follows 

because whenever Au>l, or equivalently Au XBAX, eqn. (All) reduces to 

Av (t +At ) X (t +At) = 

Au(t+At) 

and whenever Au 4, or equivalently Au X aAX , it reduces to 

X(t+At) = X(t) + Av(t+At) = X(t) + AX(t+&) 

when keeping only the first term in the asymptotic expansion, or to 

X(t+At) =X(r) +%[Av(t+Af) + Av(t)] 

= X(t) + I/z[AX(t+At) + AX(t)I 

10 
when keeping the fust two terms in the asymptotic expansion, as in eqn. (All). Equation (A12) is the asymptotic solution to eqn. (Al) whenever Au Xs-AX, and eqns. (A13 and A14) are valid whenever Au XaAX. In 

fact, eqn, (A13) is the Euler backward difference method, while eqn. (A14) is the trapezoidal difference method. 

Because keeping the second term of the asymptotic expansion comes with negligible computational expense per 

iteration (one must only retain Av(f) from the previous time step), and because numerical experiments show 

that this reduces the number of iterative steps required to attain convergence, it is therefore recommended that 

the second term in the asymptotic expansion be retained in numeric algorithms. 

Viscoplastic models are composed of systems of stiff, linear, first-order, differential equations with coupled 

nonlinear parameters. An algorithm for the numeric solution of these models using the implicit integration 

scheme of eqn. (All) (where only the first term in the asymptotic expansion was kept) and the iterative technique of Newton-Raphson is given in ref. 6. An earlier algorithm of WALKER’S [ll (where the derivatives of u 

and v are evaluated at time t instead of time t+At) is not as efficient. 

Even though this integration scheme is stable and accurate for time steps of arbitrary size from the numerins viewpoint, it will still be necessary to moderate the time step size from a physics viewpoint in order to properly trace the loading history, especially when it is nonproportional such that the principal axes of tensor-valued 

variables rotate in a material reference frame. The evolutionary equations in unified viscoplastic formulations 

introduce an evanescent strain memory effect (evanescent along the inelastic strain path). This effect is consistent with L’WSHI”s [30] delay-trace hypothesis, that a material’s response depends not on the whole of the 

previous strain path, but only on the most recent part of it. This physical fact is assimilated in the integration 

algorithm by expanding the exponential integrand in eqn. (A6) about the current state. 

APPENDIX B 

MODEL DEVELOPMENT AND CHARACTERIZATION 

Here we derive the material functions given in eqn. (6) that are unique to this model, and present a methodology for determining values for the material constants of the model, 

MILLER’S HYPOTHESIS 

The Zener parameter Z in the flow equation (3a) is a temperature normalized function describing the 

magnitude of inelastic strain rate, i.e., Z = IiP.Pv8. At steady state, one must be able to express Z as only a 

function of stress, since the internal variables do not evolve at steady state. Therefore, one can propose two distinct Zener parameters, each differing in argument and function; one for steady states, and the other for transient 

states. If one defines these arguments as 

z,=z,[!$ , .=.[?I IS-B I 

then the hypothesis of Miller [21 given by 

establishes the functional dependence of the transient Zener parameter Z 

(B2) 

(the harder of the two to characterize) 

in terms of the steady-state Zener parameter Z,. Here A > 0 is a material constant introduced to normalize 

the units, and rn 2 1 is the constant-structure exponent. 

The advantages of this hypothesis go beyond just providing the functional form for the transient Zener 

parameter. One can now equate arguments between the Z of eqn. (B2) and the Z, of eqn. (Bl) under steady 

state conditions. Assuming that the steady-state value for back stress is proportional to the value of applied 

stress, i.e., 

PI, = b IS1 033) 

then it immediately follows that 

- m-1 

D, = A(1-b)[ F] 034) 

11 
independent of the choice of functional form for Zss. Actually, Miller proposed eqn. (B4) and then derived eqn. 

(B3), where the proposition of eqn. (B4) was based on data from warm-worked material. These relationships 

play an important role in characterizing the thermal recovery functions. Note that m = 1 implies 

D, = A (1-b) = constant, but in the model of eqn. (6) m > 1 because of the recovery function r . 

MATERIAL FUNCTIONS 

A steady-state Zener parameter that is the sum of two power-law functions, i.e., 

n’ 

z,, = [ y] + [ -4 ”[ 9 

is proposed, and we motivate this choice with the steady-state creep data of fig. 1. Applying Miller’s 

hypothesis, eqn. (B2), to this form for Z, results in the transient Zener parameter Z given in eqn. (6b). The 

function for thermal diffusivity, eqn. (6a), was derived by MILLER [21. 

The limit function L for the short-range back stress controls the degree of curvature exhibited at the 

knee of the stress-strain curve. Equation (6c) assumes that L is proportional to the value of the drag strength. 

This relationship is in accordance with the physical observation that a material’s capacity to sustain an internal 

stress (of which L is the measure) is directly related to the density and distribution of its dislocation substructures (of which D is the measure). 

By definition Bij = 0 at steady state, and therefore from the constitutive equation (2b) one obtains 

from which one derives the thermal recovery function R for the back stress. Here it is considered that f = 0 

at steady state, which implies that HI = H2 = 0 at steady state. Combining this relationship with the evolution 

equations (3) for and p;, and using the fact that Xi,, Bi;, and Si, are coaxial with each other at steady 

state (a direct consequence of the chosen evolution equations), leads to 

which when combined with eqns. (B3 and B5) results in the expression for the thermal recovery function R 

given in eqn. (a). 

At steady state, D = 0 by definition, and therefore the evolution equation (3d) for drag strength reduces 

to 

where eqns. (B3 and B4) are used to reduce the coefficient IS-B I/D in the middle expression. Combining 

eqns. (B4. B5, and B8) results in the thermal recovery function r for drag strength presented in eqn. (Sf). 

MATERIAL CONSTANTS 

For a complete material description using the model presented herein, there are: three elastic constants: 

E, a, and v; two thermal constants: Q and TI; four steady state constants: A , A ’, n , and n ’; the constantstructure constant m ; two threshold constants: b and c ; three transient constants: h , H and H2; the nonproportionality constant: a; and the initial value: DO. Their values are given in table 1 for aluminum, copper, 

nickel, and tungsten. There are seventeen independent material constants in all, which to a good approximation 

are temperature insensitive. There are a number of rules of thumb (or estimates) given below that one can use to 

assign values to the various constants if data are unavailable. 

The activation energy for self diffusion Q > 0 and the melting temperature T, > 0 are typical handbook data. A useful estimate is that the transition temperature Tt > 0 associated with the activation energy for 

inelastic flow is given by TI = O.6Tm with actual values ranging between OST, and O.8Tm [31]. 

12 
With the thermal diffusivity now characterized, one can determine the steady-state Zener parameter 

Z, = IEp:PISs/€l from creep data and plot it against its associated flow stress, as done in fig. 1. The lines in this 

figure represent the fit of eqn. (B5) to the data. Values for the constants A > 0, A’ > 0, n > 0, and n‘ > 0 are 

readily obtained by eye, or by the method of least squares. Since creep data are probably the most abundant of 

all high-temperature experimental data, one should have no difficulty in characterizing the Norton creep constants A and n (n = 5, typically). High-stress creep data, however, are sometimes scarce or unavailable. If 

this is the case, one can use the estimates given by: n‘ = 3.51 and Zco = (A’/A )nn”(n’-n) = 10’os-l, with 

values ranging between 3n to 4n for n’ and lo9 to 10” s-l for the Zener parameter Z,, at cross over (cf. 

The remaining steady state constant, the Bauschinger coefficient bE(0,l) defined in eqn. (B3), can be 

obtained from a tensile test where a specimen is loaded to its saturation (or steady state) strength qat, and then 

unloaded sufficient enough to induce some reversed plasticity. One then obtains a value for the Bauschinger 

coefficient from the relationship b = (os,+ o,)/2a,,, where ov is the reversed yield strength. A word of 

caution. A small offset strain should be used to establish 6, to ensure a realistic value for 6. Values for b 

may range from 0.2 to 0.7, with a typical value being 0.4. The larger the value for 6, the greater the predicted 

Bauschinger effect. 

There is one more material constant that one can ascertain from creep testing. It is the constant structure 

exponent m > 1 which can be determined from stress-drop tests. Fortunately, this parameter is nearly insensitive to material selection. Typical values for the constant structure exponent are: m = 1.8 for pure metals, and 

m = 1.5 for alloys [13]. 

The next step is to obtain those constants one can from stable stress-strain hysteresis loops, where the isotropic variable has saturated. Ratios of plastic strain range to elastic strain range of 1:l to 3:l generally work 

well. The constant H2 > 0 represents the tangent modulus in the region of linear strain hardening after going 

around the knee of the curve, with H2=Hl/100 serving as a rule of thumb. The constant H1 > 0 defines the 

tangent modulus at the onset of inelastic deformation. Our experience from optimization calculations is that 

H1 = E15 provides a good estimate for its value. The larger the value for the material constant c , the more 

rounded the knee in the hysteresis loop will be. The method of nonlinear least squares [32] is a useful tool, but 

not a necessary one, for determining c and HI (and D,, since the drag strength is an unknown constant at this 

time) since there are three constants to be determined simultaneously at this step of material characterization. 

The constants Do > 0 and h > 0 are characterized next. The constant h establishes the rate of isotropic hardening in monotonic tension to large strains. The initial value Do for drag strength is obtained by comparing predicted and experimental values for the yield strength. Some iteration between these two constants 

may be necessary to obtain a pmper fit. Again, the method of nonlinear least squares may be helpful, but it is 

certainly not necessary. 

The final material constant to quantify is the nonproportionality parameter a~[O,l]. This requires 

nonproportional experimental data. If they do not exist, which is likely to be the case, setting a = 0 implies that 

the theory resembles a two-surface plasticity model with Mroz hardening C331. Further experimental study is 

required before typical values for a can be given. 

fig. 1). 

13 
TABLE 1 

Al I cu I NI 

~ 

Elastic Constants$ 

w 

EO 

El 

E2 

a0 

a1 

V 

GPa 71 

MPUC -36 

.34 

MPapc2 -.011 

1pC 25x10" 

1pc2 

127 

41 

-.027 

.34 

16~10~ 

5x104 

$ 

210 

-2.5 

-.063 

.3 1 

13x10" 

E =Eo+EIT+E2T2 

a = + a1T 

T is in degrees centigrade. 

395 

-2.3 

-.027 

.28 

4.4~ 1 0" 

.16x104 

I I Material 

* 

Constant 

A 

A' 

b 

DO 

Hl 

H2 

mt 

n 

n' 

Q 

Tt 

c 

h 

constant a, nor is there adequate experience to 

assign values via a rule of thumb. 

Determined by engineering judgment because of a lack 

of data. 

Units 

MPa 

MPa 

MPa 

MPa 

GPa 

MPa 

kJ/mole 

OK 

Al 

.076 

4.4 

.6 

23 

.1 

.OS 

10 

100 

1.8 

4 

14 

140 

560 

cu 

1.2 

18 

.35 

7 

1.2 

3 

20 

200 

1.8 

5.5 

19 

200 

820 

I Data are not available to characte 

Ni 

.55 

19 

.5 

15 

.5 

2 

70 

700 

1.8 

5 

15 

290 

looot 

w 

.29 

26 

.4t 

lot 

.3t 

5t 

loot 

1Wt 

1.8 

4.5 

19 

590 

2600 

ze the material 

14 
,- MATERIAL REFERENCE , 

100 

80 

m 60 

P 

v) 

v) 

w z

40 

20 

FIGURE 1. - STRESS DEPENDENCE OF ZENER PARAMETER 

FOR ALUMINUM, COPPER, NICKEL. AND TUNGSTEN. 

STRAIN. 

E 

MODEL PREDICTIONS 

USED TO DETERMINE 

MATERIAL CONSTANT h 

0 

-200 -100 0 100 200 300 400 500 

TEMPERATURE. OC 

FIGURE 2. - MONOTONIC STRESS-STRAIN RESPONSE OF 

ALUMINUM AS A FUNCTION OF TEMPERATURE. STRAIN 

RATE. h I 0.000128/s. 

15 
MODEL PREDICTION 

EXPERIMNTAL DATA 

(REF. 27) 

100 - 

-50 - 

-100 T 

n" 

z 

k 

W oz 

100 - 

50 - 

(a) FIRST TWO CYCLES AND SATURATED CYCLE. 

WALKER METHOD FOR 

VARIOUS TIME STEP SIZES, 

DEG C/STEP 

1.5 

0 15 

0 50 

_____----- RUNGE-KUTTA ETHOD 

T 

-.003 -.002 -.001 0 ,001 .002 .003 

PECHAN I CAL STRAIN 

(b) THEORETICAL PREDICTIONS USING INTEGRATION 

METHODS OF WALKER AND RUNGE-KUTTA. 

FIGURE 3. - INPHASE THERmmECHANICAL RESPONSE 

OF NICKEL. CONTROL IS VIA LINEAR RAMPS IN 

TEMPERATURE AND MECHANICAL STRAIN WITH A 

600s CYCLE TIE. TEMPERATURE RANGE, 400 TO 

700 'C. ECHANICAL STRAIN RANGE. fo.0026. 

ao2 r 

-0.2 

-.008 ,004 0 * 004 .008 

AXIAL STRAIN. 

(a) STRAIN CONTROL PROGRAM. PATH SEGPENTS SEQUENCE, 0 - 1 - 2 - ... 4 8. EACH PATH 

SEGMENT IS APPLIED IN 60s. 

NONPROPORT I ONAL 

MTERIAL CONSTANT. 

a 

.75 

.9 

. . . . . . . . . 0 

1 

EXPERIMENT (REF. 29) 

200 - T 

71 

L 

- 

-200 + I I I 

-200 -100 0 100 200 

AXIAL STRESS, MPa 

(b) EXPERIMENTAL AND PREDICTED STRESS RESPONSE, 

D = 6.5 MPa IN THESE CALCULATIONS. 

FIGURE 4. - CYCLIC, NONPROPORTIONAL, AXIALTORSIONAL TEST OF COPPER AT ROOM TEMPERATURE. THE TEST SPECIMEN WAS CYCLED TO 

ISOTROPIC SATURATION PRIOR TO TESTING. 

16 
1. Report No. 

NASA TM-102338 

4. Title and Subtitle 

Refinements in a Viscoplastic Model 

2. Government Accession No. 

7. Author@) 

A.D. Freed and K.P. Walker 

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No of pages 

Unclassified Unclassified 18 

9. Performing Organization Name and Address 

National Aeronautics and Space Administration 

Lewis Research Center 

Cleveland, Ohio 44135-3191 

22. Price' 

A03 

12. Sponsoring Agency Nfme and Address 

National Aeronautics and Space Administration 

Washington, D.C. 20546-0001 

15. Supplementary Notes 

3. Recipient's Catalog No. 

5. Report Date 

6. Performing Organization Code 

8. Performing Organization Report No. 

E-4988 

10. Work Unit No. 

505-63-1B 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 

Prepared for the 1989 Winter Annual Meeting of the American Society of Mechanical Engineers, San Francisco, 

California, December 10-15, 1989. 

16. Abstract 

A thermodynamically admissible theory of viscoplasticity with two internal variables (a back stress and a drag 

strength) is presented. Six material functions characterize a specific viscoplastic model. In the pursuit of 

compromise between accuracy and simplicity, a model is developed that is a hybrid of two existing viscoplastic 

models. A limited number of applications of the model to Al, Cu, and Ni are presented. A novel implicit 

integration method is also discussed, and applications are made to obtain solutions using this viscoplastic model. 

17. Key Words (Suggested by Author@)) 

Viscoplasticity; Copper; Nickel; 

Aluminum; Tungsten 

18. Distribution Statement 

Unclassified - Unlimited 

Subject Category 39 

*For sale by the National Technical Information Service, Springfield, Virginia 221 61 NASA FORM 1626 OCT 86