Viscoplastic augmentation of the smooth cap model

:?;'All ELSEVIER

Nuclear Engineering and Design 150 (1994) 215-223

Nuclear Engineeri.ng and Design

Viscoplastic augmentation of the smooth cap model Leonard E. Schwer Schwer Engineering & Consulting Services, 455 South Frontage Road Suite 112, Burr Ridge, IL 60521-7104, USA

Abstract

The most common numerical viscoplastic implementations are formulations attributed to Perzyna. Although Perzyna-type algorithms are popular, they have several disadvantages relating to the lack of enforcement of the consistency condition in plasticity. The present work adapts a relatively unknown viscoplastic formulation attributed to Duvaut and Lions and generalized to multi-surface plasticity by Simo et al. The attraction of the Duvaut-Lions formulation is its ease of numerical implementation in existing elastoplastic algorithms. The present work provides a motivation for the Duvaut-Lions viscoplastic formulation, derivation of the algorithm and comparison with the Perzyna algorithm. A simple uniaxial strain numerical simulation is used to compare the results of the Duvaut-Lions algorithm, as adapted to the DYNA3D smooth cap model with results from a Perzyna algorithm adapted by Katona and Muleret to an implicit code.

1. Introduction

The dynamic behavior of geological materials such as rocks, concrete and soils is rate dependent. For many engineering applications, the loading rate in the geological media is slow enough to allow rate-dependent effects to be neglected. This significant class of applications, and a wealth of quasi-static laboratory data describing the behavior of geological materials, have led to the development of many sophisticated geological material constitutive models. For numerical simulations, perhaps the most popular of these geological constitutive models is the multi-surfaced cap model, whose numerical algorithm is usually based on the work of Sandier and Rubin (1979).

Compared with quasi-static laboratory data for geological materials, the quantity and quality of dynamic laboratory data are minimal. The dynamic data that do exist are usually limited to very high strain rate, uniaxial strain data obtained in gas gun experiments. While these data are useful, they do not provide a very interesting stress path for exercising sophisticated material models. Furthermore, these data are not typical of the stress paths of interest in geological material dynamic loading applications, such as underground explosions. In a companion work (Schwer, 1993), dynamic loading laboratory data obtained from a spherical wave test in dry limestone (Klopp, 1992) are compared with the numerical results of the proposed viscoplastic cap model implementation.

0029-5493/94/$07.00 © 1994 Elsevier Science S.A. All rights reserved SSDI 0029-5493(94)00734-G

216

L.E. Schwer / Nuclear Engineering anti Design 150 (1994) 215 223

The two most popular numerical viscoplastic formulations are the elastic-viscoplastic forms attributed to Perzyna (1963, 1966, 1971) and the endochronic forms usually attributed to Valanis (1971). Viscoplastic cap models based on the Perzyna formulation have been proposed by Katona and Muleret (1984) and by Simo et al. (1988a). The present work adapts a relatively unknown viscoplastic formulation attributed to Duvaut and Lions (1972) and generalized to multi-surface plasticity by Simo et al. (1988b). The attraction of the Duvaut-Lions viscoplastic formulation is in its ease of numerical implementation, especially in existing elasto-viscoplastic constitutive models. To convert an existing elasto-plastic constitutive algorithm to a viscoplastic constitutive algorithm requires a simple stress update loop through the integration points at the end of the existing routine. In the next section, motivation for the DuvautLions viscoplastic formulation is presented that follows loosely the detailed development presented by Simo and Hughes (in press). Section 3 presents a derivation of the Duvaut-Lions algorithm, while Section 4 provides a comparison of the Perzyna and Duvaut-Lions algorithms, listing the relative advantages and disadvantages of each algorithm. The numerical examples in Section 5 include a comparison of the Perzyna and Duvaut-Lions viscoplastic simulations of a hypothetical uniaxial strain loading of a McCormick Ranch sand sample. Details concerning the inviscid smooth cap model (Schwer and Murray, in preparation) used as the 'base' elastoplastic model for the DuvautALions formulation are not presented. The three-invariant inviscid cap model is based on the smooth cap model formulation of Pelessone (1989) with a modification suggested by Rubin (1991) for efficiently including the third stress invariant. The resulting viscoplastic three-invariant smooth cap model has been implemented into the public domain versions of the Lawrence Livermore National Laboratory explicit finite element codes DYNAZD (Whirley, 1992) and DYNA3D (Whirley, 1991). 2. Motivation

for Duvaut-Lions

formulation

The simplest mechanical analog of a viscoplastic material is a spring, dashpot and Coulomb

I

1

I

E

I

I////1

~y 11

I

I

I

t

Viscoplastic

Elastic

Fig. 1. Mechanical analog of elastic viscoplastic material.

friction element. A schematic diagram of a one-dimensional Maxwell viscoplastic model is shown in Fig. 1. In the Maxwell model, the stress is the same in the elastic and viscoplastic components, and the strain rates are additive. The stress in the spring is given by the Hookean relationship O" ~-- ESspring = E e e

(1)

and the stress in the dashpot is given by the viscous damper relationship (2)

O" = q~dash = q~vp

The Coulomb friction device acts as a stress limiter, such that the stress never exceeds O-y. Furthermore, to make the model elastic-viscoplastic, as opposed to viscoelastic, the Coulomb device prevents straining in the viscoplastic element until the stress equals the limit stress, i.e. rigid viscoplastic. Thus, the stress-strain rate relationship for the viscoplastic component can be expressed as ~,vp =

t0

if O-< O-y if O- > O-y

O" -~- O ' y

(3)

The total strain rate across the viscoplastic Maxwell element is obtained by summing the elastic and viscoplastic rates of the components, i.e. = ~ + ~vp=~_~

O- - - O'y

~/

(4)

The above equation, valid when O-> O-y,is usually written in the alternative form (9"

E~ = ~ +

--

O'y

(5)

L.E. Schwer /Nuclear Engineering and Design 150 (1994) 215-223 I

strain rate formulations. Recasting the Coulomb friction element in terms of a plasticity yield function, i.e.

I

f(tr) = a

Yr, .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

......

217

-

try

Eq. (3) can be rewritten as .

.

.

.

.

~vp =

(tr)/n

iff(tr) < 0 (elastic) iff(tr) ~>0 (viscoplastic)

(8)

This can be written more compactly by introducing the ramp function notation Time Fig. 2. One-dimensional viscoplastic relaxation solution.

where z is called the relaxation time and is given by 1 - =

E -

(x + Ixl) 2

where 1" I is the absolute value. Now, we can write the above viscoplastic strain rate in the form typically associated with Perzyna-type viscoplastic formulations:

(6)

~vp _ (f(tr) > Of(tr) Otr

2.1. Stress relaxation response

(9)

2.3. D u v a u t - L i o n s - t y p e formulation

To examine the response of this simple viscoplastic model, let us consider the solution of Eq. (5) under the condition of a suddenly applied and then held constant strain, i.e. simple stress relaxation. Furthermore, we shall assume that the initial strain is large enough to cause an initial stress that exceeds the limit stress, i.e.

The Duvaut-Lions viscoplastic strain rate formulation follows directly from the ideas used in formulating the one-dimensional model problem. Let us recall the definition of the viscoplastic strain rate, i.e. ~vp _ (7 - - try

Eeo = a(0) > try Eq. (5) can be solved subject to this initial condition, using an integrating factor of exp(t/r) to yield the solution

Now, rewriting the above using the relaxation time, we obtain

tr( t) = ( Eeo - try)exp( t /z) + try

~vp __

E-l

(7)

Fig. 2 shows a schematic diagram of the above stress relaxation solution. The basic nature of the relaxation response is that the 'overstress' tr(0) decreases within increasing time, at an exponential rate controlled by the relaxation time parameter z. 2.2. Perzyna-type formulation

The simple one-dimensional viscoplastic model can be expressed in the form associated with the Perzyna (1963, 1966, 1971) type of viscoplastic

(tr - - try)

(lO)

T

which represents the Duvaut-Lions form of the viscoplastic strain rate. Although this forms seems trivial, it has great advantages over the Perzyna form in numerical implementations, as will be shown. To introduce the idea underlying the D u v a u t Lions form, and prepare the reader for the following sections, Eq. (lO) can be considered as stating that the viscoplastic strain rate is related to the difference between the current (viscoplastic) stress tr and the solution to the corresponding inviscid

L.E. Schwer / Nuclear Engineering and Design 150 (1994) 215 223

218

problem, represented by k; vp = E - l ( o

O'y,

such that

-

a.+~=a.exp

+

- - O'y)

~t

As time increases, the current stress relaxes to the inviscid solution and the viscoplastic strain rate decreases, at a rate controlled by the relaxation time. The algorithmic implications of the D u v a u t - L i o n s form of viscoplasticity are developed in the next section.

+ exp

z

+ C: ~ ds

(14)

By introducing the approximations

fl xp(

S),c ,as

3. Derivation of Duvaut-Lions algorithm The three-dimensional generalization of the D u v a u t - L i o n s viscoplastic strain rate formulation (Simo, in press) is given by

1 ~vp = _ C - l ( a _ 6) d =-

1

D-l(q

- q)

T

= 1 - exp( -- At/z) (C: Ann + ~)

(15)

_-{, expt

,16,

At/z

(11)

(12) the resulting viscoplastic stress update is obtained

where symbols in bold indicate matrices or equivalently tensors, e.g. a = a~/. The overbarred quantities 6 and q represent solutions of the inviscid elastoplastic problem. The quantities q represent internal hardening variables associated with motion and/or growth of the yield surface. The quantity ~ represents the evolution of the hardening quantities and is related through the diagonal matrix O by

as

a,+l=exp

+

-

-

a,+exp

1 - exp( - At/z) At/z

C:A~,+I

6,+

(17)

The corresponding hardening parameter update

q = D~ From the assumption that the strain rate can be separated into elastic and viscoplastic components, the stress rate can be written as = C : (~ - ~vp) = C:,~ - 1

(a - •)

(13)

T

where the tensor contraction is indicated by the colon (:) notation, e.g. G 0. = C : ~ = C , j , , ~ ,

As in the one-dimensional viscoplastic case, the stress rate equation (Eq. (13)) can be integrated in closed form over a time increment At, using an integrating factor, to yield

• ,,+]

= --D-lq.+

j

(19)

is obtained in a similar manner. Thus, the D u v a u t - L i o n s viscoplastic algorithm requires the solution of the corresponding inviscid elastoplastic problem to obtain the quantities 6, + 1 and q, + 1, followed by simple 'updates' using these quantities and the viscoplastic solution from the previous time step, i.e. a, and qn. 3.1. Alternative update f o r m

An even more remarkable viscoplastic update algorithm can be obtained by direct application of an implicit backward Euler algorithm to the D u v a u t - L i o n s viscoplastic strain rate (Eq. (11)),

L.E. Schwer / Nuclear Engineering and Design 150 (1994) 215-223

giving a.+l

219

3.3. More elaborate Duvaut-Lions models

--#n+l

=~

T

In Simo et al. (1988b), a remark is made, suggesting the construction of more elaborate Duvaut-Lions models. The example given is based on a family of monotonic C l functions g, used to express the viscoplastic strain rate as

C:A/~vP

= At C : (A~ -- At;e) ~7 [ trial -

At~a.+l - - a . + l )

~v~ = g ( l l ~ - ~ Ib) C _ l . ( a

where use has been made of the standard elastoplastic algorithm definitions an+l =an

The above update form may be rearranged and written as trial

a,+ l + (At/z)#n+ 1 1 + At/z

where

q trial

, + l + (At/z)~ +1 1 + At/z

is the energy norm. Mould and Levine (in press) proposed a practical modification of the Duvaut-Lions viscoplastic strain rate given by

(20)

Similarly, for the hardening parameters, we have q,+l =

~)

q- C : A / ~ e

trial ~'n+l = an + C:A~

a. + 1 =

_

(21)

The algorithmic implication of this alternative viscoplastic update is remarkable. The viscoplastic solution can be obtained from an inviscid model by temporarily saving the trail stress a trial and, after calculating the inviscid solution 6, applying the simple viscoplastic update given by Eq. (20).

C-l(a -a) '~vP = ~b(f)IIC-:~- =~11= where I1' I1~ is the Euclidean norm, e.g.

Ilall~ = (,~,~) l~ and ¢ ( f ) is a piecewise linear rate function that controls the magnitude of the viscoplastic strain rate as a function of the overstress. The function f is given by f=

II(a - ~)112

f'c 3.2. Elastic and inviscid solutions as limiting cases

Another advantage of the Duvaut-Lions viscoplastic formulation is that the elastic and inviscid solutions are preserved as limiting cases of the visoplastic solution. This is not necessarily the case when applying Perzyna-type formulations to multi-surface plasticity, such as traditional captype models (Simo, 1988b). It is easily shown from Eq. (20) that, in the limit as trial is obA t / z ~ O , the elastic case a.+l an+l tained; At/z ~ ~ , the inviscid case a.+l --*0.+ l is obtained.

where f'c is the unconfined compressive strength. This proposed piecewise linear viscoplastic strain rate model can be fit to experimental data covering the strain rate range of interest. Implementation of these more elaborate Duvaut-Lions models requires the solution of auxiliary non-linear equations to update the stresses. For example, the Mould and Levine model results in the non-linear scalar equation At IIC - l ( a t r i a l - - ~)I]2 ~b(ftrial) = 0

where the parameter e provides the stress update a . + l = (1 - ~)atnal + e a n + l

L.E. Schwer / Nuclear Engineering and Design 150 (1994) 215 223

220

4. Comparison of Perzyna and Duvaut-Lions algorithms

5. Numerical examples

In this section, a brief outline of the numerical algorithm for the Perzyna and D u v a u t - L i o n s viscoplastic algorithms is presented. The relative advantages and disadvantages of each algorithm are given.

4.1. Perzyna algorithm Let us assume that, at time t,, all the quantities are known and the applied strain increment A~, + l is given then an explicit Perzyna algorithm would proceed as follows: A ~ p = AtiVp Adr = C:(A/~n+ ] - Aenvp) drn + 1 = drn -]- Adr

(22)

L + J =f(dr,+ 1,/;vp)

~+ l

(49(f'+ l) 5 £~+ l q

-

where Eq. (22) is subject to stability criteria (Cormeau, 1975) for convergence of the explicit algorithm.

4.2. Duvaut- Lions algorithm

Katona and Mulert (1984) used their Perzynatype viscoplastic cap model to simulate the hypothetical uniaxial strain loading history shown in Fig. 3. The cap model material properties are those of McCormick Ranch sand, as given by Sandler and Rubin (1979). A summary of the properties used by Katona and Mulert is given in Table 2. The form of the Perzyna viscoplastic flow function used by Katona and Mulert is

where y is the fluidity paramter, f0 = 0.25 ksi and N=I. In their paper, Katona and Mulert presented viscoplastic axial stress histories calculated for three values of the fluidity parameter y, spanning two decades, i.e. y = 0.1-0.01 (from nearly inviscid to moderately viscous). Fig. 4 shows the viscoplastic axial stress history for y = 0.001 from Katona and Mulert, compared with the viscoplastic axial stress history for r = 0.25 calculated with the D u v a u t - L i o n s viscoplastic smooth cap model. Overall, the two viscoplastic models agree, with the possible exception of the times greater than 5.5. This agreement is quite gratifying considering

Let us assume that, at time t,, all the quantities are known and the applied strain increment A~n + ] is given. Then, an explicit D u v a u t - L i o n s algorithm would proceed as follows:

0.04

I

I

I

[

I

I

0.035

trial

dr,,+]

=

drn

+

C :A~n+ ]

#. + 1'-- from elasto-plastic algorithm trial dr.+] + ( At/~)#.+ ]

drn+ 1 --

1 + At/r

4.3. Comparison of Perzyna and Duvaut-Lions formulations

v ¢,-

0.025 0.02 0.015 0.01 0.005 0

The main advantages and disadvantages of the Perzyna and D u v a u t - L i o n s formulations are summarized in Table I.

\

0.03

,z-

I 1

I 2

I 3

I 4

I 5

I 6

Time

Fig. 3. Uniaxial strain loading history.

I

L.E. Schwer / Nuclear Engineering and Design 150 (1994) 2 1 5 - 2 2 3

Table 1 Advantages and disadvantages of Perzyna and D u v a u t - L i o n s formulations Advantages Perzyna formulation Has been in use longer and, hence, some model fits to c o m m o n metals and geological materials m a y be available

New functional relationships, i.e. $ ( f ) , can be easily added

Disadvantages

The resulting viscoplastic algorithm is not dependent on the material's corresponding elasto-plastic algorithm, because the consistancy condition is not enforced; in most cases, a new viscoplastic material model is introduced rather than augmenting the existing elasto-plastic model In multi-surface plasticity, such as traditional cap-type models, the Perzyna viscoplastic algorithm m a y not reduce to the rate-independent solution (Simo, 1988b) For certain material models and material properties additional stability criteria m a y need to be satisfied; these stability limits are important, since they are often independent of mesh size

Duvaut - L i o n s formulation Easily adapted to existing elasto-plastic material models

Converges to inviscid and elastic cases in the limit, including n o n - s m o o t h multisurface plasticity models

Introduction of more elaborate D u v a u t - L i o n s models requires modification of the algorithm to include at least the solution of a non-linear scalar equation There is no direct conversion of material model fits using the Perzyna formulation to D u v a u t - L i o n s formulations; this is a minor disadvantage that should be remedied as users adopt the more efficient D u v a u t Lions formulation and publish data for their material model fits

Table 2 McCormick Ranch sand cap model data

Elastic moduli K = 66.7 ksi G = 40.0 ksi

Cap surface X ( x ) = L(x) + RFe(x ) R =2.5 X = 0.189 ksi

Shear failure surface J'2 I/2 = ct - V exp( - flJ1 ) + OJl = 0.25 ksi ~, = 0.18 ksi fl = 0.67 ksi-i 0 = 0.0

Cap hardening e~(X) = IV[ 1 - exp{ -- D 1( X - Xo) -- D2( X -- Xo)2}] W = 0.066 D l = 0.67 ksi -1 32 = 0

Tension cutoff (not used)

221

222

L.E. Schwer / Nuclear Engineering anti Design 150 (1994) 215 223

1.5

I

l

I

I

[

^

.-~

I

I

f

I

.

.

.

.

t

3.5

Perzyna 7=0.001

=.5

3

[ I I - -

/...

2.5

- -

v

~

//''\

1.5

\

-

'<

-0.5

I

J

1

I

I

I

[

1

2

3

4

5

6

7

-

('~->~)

x = 1 , 0 0

x=0.500

.....

U~

<

-

.....

2

0.5

Elastic - -

z=0.250 -

-

.....

x = 0 . 1 2 5

Inviscid ('t->O)

0.5

-0.5 0

Time

I 1

I 2

I 3

I 4

J 5

/1

I 6

7

8

Time

Fig. 4. Comparison of Perzyna and Duvaut-Lions axial viscoplastic stress.

Fig. 5. Visoplastic axial stress response from Duvaut Lions smooth cap model.

the following: the two viscoplastic formulations are quite different, although the value of z = 0.25 was chosen to approximate the maximum stress magntiude calculated by Katona and Mulert; the cap models are slightly different; the integration schemes were different, Katona and Mulert using an implicit scheme, while an explicit scheme is used here. As mentioned previously in Section 3.2, the D u v a u t - L i o n s viscoplastic formulation preserves the elastic and inviscid solutions in the limit. Fig. 5 shows the elastic and inviscid axial stress histories for the uniaxial strain history (shown previously in Fig. 3), as calculated with the smooth cap model. These two stress histories were calculated using the viscoplastic augmented smooth cap model, but with the viscoplastic option turned off; the elastic case was calculated with the failure and cap surfaces shifted away from the stress trajectory to prevent intersection with the stress trajectory. Fig. 5 also shows the viscoplastic axial stress histories for four values of the viscosity parameter 3. The range of values of the viscosity parameter and the corresponding viscoplastic stress histories, compared with the elastic and inviscid stress histories, provide a guide to how the D u v a u t - L i o n s viscoplastic formulation approaches the limiting cases and the sensitivity of the viscosity parameter.

6.

Conclusions

The D u v a u t - L i o n s viscoplastic formulation provides the framework for an efficient numerical algorithm and an extremely easy implementation in existing elasto-plastic material models. The basic D u v a u t - L i o n s form implemented into the DYNA3D smooth cap model is shown to provide results very similar to those of the corresponding basic Perzyna-type form. Although not discussed in the present work, application of the D u v a u t - L i o n s viscoplastic smooth cap model to the more complex problem of simulating spherical wave propagation experimental results in limestone has also been very successful. Details of the experiment and simulation are provided in the companion work (Schwer, 1993).

Acknowledgments

The generosity of Professor Juan Simo in providing a draft copy of his book is gratefully acknowledged, as is the DYNA2D adaption of the DYNA3D viscoplastic three-invariant smooth cap model performed by Dr. Robert Whirley of LLNL. Support for this research was provided by a contract from the Defense Nuclear Agency under the technical supervision of Dr. Paul Senseny.

L.E. Schwer / Nuclear Engineering and Design 150 (1994) 215-223

References I. Cormeau, Numerical stability in quasi-static elasto/visco plasticity, Int. Numer. Methods Eng. 9 (1975) 109-127. G. Duvaut and J. L. Lions, Les Inequations en Mechanique et en Physique, Dunos, Paris, 1972. M.G. Katona and M.A. Mulert, A viscoplastic cap model for soils and rock, C.S. Desai and R.H. Gallagher (eds.), Mechanics of Engineering Materials, Wiley, New York, 1984, pp. 335-350. R.W. Klopp, A.L. Florence, J.K. Gran and J.W. Simons, Spherical wave tunnel (SWAT-I) test in a large Indiana limestone specimen, SRI Int. Tech. Rep. DNA-TR-92-000 (prepared for the Defense Nuclear Agency under Contract DNA-001-90-C-0032), August 1992. J.C. Mould Jr., and H.S. Levine, A rate-dependent three invariant softening model for concrete, in G.Z. Voyiadis, L.C. Bank and L.J. Jocobs (eds.), Mechanics of Materials and Structures, Elsevier, Amsterdam, in press. D. Pelessone, A modified formulation of the cap model, Gulf Atomics Rep. GA-C19579 (prepared for the Defense Nuclear Agency under Contract DNA-001086-C-0277), January 1989. P. Perzyna, The constitutive equations for rate-sensitive plastic materials, Quart. Appl. Math. 20 (1963) 321-332. P. Perzyna, Fundamental problems in viscoplasticity, Advances in Applied Mechanics, Vol. 9, Academic Press, New York, 1966, pp. 244-368. P. Perzyna, Thermodynamic theory of viscoplasticity, Advances in Applied Mecahnics, Vol. 11, Academic Press, New York, 1971. M.B. Rubin, Simple, convenient isotropic failure surface, Eng. Mech. Div. Am. Soc. Civil Eng. 117 (2) (1991) 348-369.

223

I.S. Sandier and D. Rubin, An algorithm and a modular subroutine for the cap model, Int. J. Numer. Anal. Methods Geomech. 3 (1979) 173-186. L.E. Schwer, Spherical wave propagation in limestone: simulation of the SRI SWAT-I experiment, Schwer Engineering & Consulting Services, Tech. Rep. SECS-TR-93-000, July 1993. L.E. Schwer and Y.D. Murray, A three-invariant smooth cap model with mixed hardening, in preparation. J.C. Simo, J.-W. Ju, K.S. Pister and R.L. Taylor, Assessment of cap model: consistent return algorithms and rate-dependent extensions, J. Eng. Mech. Div. Am Soc. Civil Eng. 114 (2) (1988a) 191 218. J.C. Simo, J.G. Kennedy and S. Govindjee, Non-smooth multisurfaee plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms, Int. J. Numer. Methods Eng. 26 (1988b) 2161-2185. J.C. Simo and T.J.R. Hughes, Plasticity Viscoplasticity and Viscoelasticity: Formulation and Numerical Aspects, Springer, Berlin in press. K.C. Valanis, A theory of viscoplasticity without a yield surface, Arch. Mech. 23 (1971) 517 555. R.G. Whirley, DYNA3D a nonlinear, explicit, three-dimensional finite element code for solid and structural mechanics users manual, University of California Lawrence Livermore National Laboratory, Tech. Rep. UCRL-MA107254, May 1991. R.G. Whirley and B.E. Engelmann, DYNA2D a nonlinear, explicit, two-dimensional finite element code for solid mechanics user manual, University of California Lawrence Livermore National Laboratory, Tech. Rep. UCRL-MA110630, April 1992.