Stability and error analysis for time integrators applied to strain-softening materials

__ _k!!ii

J?J&g

ELSEVIER

Computer methods in applied mechanics and englneerlng Comput. Methods Appl. Mech. Engrg. 124 (1995) 335-363

Stability and error analysis for time integrators strain-softening materials M. Kulkarni”,

T. Belytschkob’*,

applied to

A. Bayliss’

“Center for Computational Structures Technology, NASA Langley Research Center, Hampton, Virginia, USA bDepartment of Civil Engineering, Northwestern University, Evanston, Illinois, USA ‘Department of Engineering Science and Applied Mathematics, Northwestern University, Evanston, Illinois, USA

Received 4 October 1993

Abstract

The behavior of explicit time integrators is examined for the field equations of dynamic viscoplasticity in one dimension. The analysis proceeds by linearizing the field equations about the current state and freezing coefficients so that the solution to the system of ordinary differential equations can be computed explicitly. Amplification matrices for the central difference method (with rate tangent constitutive update) and fourth-order Runge-Kutta schemes are obtained by applying these schemes to the resulting linear, constant coefficient equations. These amplification matrices lead to local truncation error estimates for the temporal integrators as well as to estimates of the resulting growth rate. Analysis of the cental difference method with a forward Euler constitutive update indicates that the local truncation error involves as the product of the third power of the strain-rate and the square of the time step. In viscoplastic models, the viscoplastic strain-rate is an increasing (decreasing) function of the viscoplastic strain in softening (hardening). Also, the eigenvalues of the rate tangent method increase with strain-rate in softening. Thus, it is recommended that the time step be reduced in localization problems to keep the strain-increment within a prescribed tolerance. As a byproduct of the analysis, conditions for numerical stability for the central difference method with a rate tangent constitutive update are obtained for strain-hardening dynamic viscoplasticity. In the softening regime, the resulting linearized differential equations admit exponentially growing solution. Thus, conventional definitions of stability such as Neumann stability and absolute stability are inappropriate. In order to characterize stability of such systems, we introduce the concept of g-rel stability which combines the concept of relative and absolute stability. We show conditions under which the resulting system is g-rel stable and thus characterize the nature of the numerical stability that can be expected in such problems.

1. Introduction

Unstable material models such as strain-softening viscoplasticity have become increasingly important because they have been found to reproduce much of the localization behavior found in dynamic experiments such as those of Marchand and Duffy [l]. When the dynamic equations for solids with unstable material models are discretized in space, systems of first-order differential equations are obtained. When localization occurs, there are high gradients in both space and time which have to be resolved in order to successfully compute the behavior of these models. In this paper we consider the behavior of explicit time integrators when they are applied to such systems. We consider both the accuracy of such integrators and appropriate concepts of stability. We consider, in particular, two standard integrators, the central difference method with a rate tangent

* Corresponding

author.

0045-7825 /95/$09.50 0 1995 Elsevier Science S.A. All rights reserved SSDI 0045-7825(95)00786-S

336

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

I24 (1995) 335-363

constitutive update and the fourth-order Runge-Kutta method. By linearizing the field equations of dynamic viscoplasticity about the current state and then freezing the coefficients, we are able to obtain explicit representations for the local truncation error. Using these estimates, we obtain a criterion for adaptive choice of the time step in such problems. We illustrate by numerical examples the behavior of different integrators and the effectiveness of the adaptive time step procedure in following the solution through severe localization. We also consider various concepts of stability for such systems of equations. We show that the resulting linear systems of ordinary differential equations admit exponentially growing solutions in the softening regime. This behavior, while only for a linearized constant coefficient system of equations, does in fact model the rapid increase in viscoplastic strain and strain-rate which occurs during localization. As a result of this growth, standard concepts of stability, such as absolute and Neumann stability are inappropriate in characterizing the stability behavior of the temporal integrators. We analyze the behavior of temporal integrators for equations admitting exponentially growing solutions. Based on this analysis, we find that the stability properties of the time integrators can be characterized in terms of relative stability. We introduce the concept of g-rel stability and show that this notion of stability is satisfied by temporal integrators and that this concept characterizes the stability behavior of numerical simulations of severe localization. In Section 2 we provide a motivation for our results by analyzing the results of computations of severe localization in one and two dimensions. In Section 3 we introduce the concept of g-rel stability. In Section 4 we consider equations where all solutions exhibit exponential growth and show that the concept of g-rel stability is appropriate. In Section 5 we consider g-rel stability for systems with an exponential dichotomy, which are characteristic of viscoplastic equations in the softening regime. The local truncation error and adaptive time step procedures are considered in Section 6 and conclusions are drawn in Section 7. In the appendix, we present a stability analysis for the central difference in the hardening regime and obtain necessary conditions for Neumann stability.

2. Motivation In this section, numerical results are presented for the following problems with non-Drucker (unstable) materials: (1) Numerical simulation of the Hopkinson torsion experiment by a finite element model of a viscoplastic slab subjected to simple shear. (2) A two-dimensional transient analysis of a strain-softening cantilever beam. Deviations from energy balance are studied in the softening regime and it is found that conservation of energy is violated when deformations begin to localize. The objective of these calculations is to motivate the analysis which follows. 2.1. Strain-softening viscoplastic slab subjected to simple shear A viscoplastic slab subjected to simple shear was studied using one-dimensional linear finite elements as shown in Fig. 1. The response of this slab to a prescribed velocity was evaluated using the central difference and fourth-order Runge-Kutta integrators. The central difference method was used with three different constitutive updates. A power law model given in Eq. (6.4a) was used. The relevant dimensionless parameters are, E/v0 = 0.002, &L/u, = 6.25 x lo-‘, pci2L2/q, = 10m3, m = 0.05, where a, is the initial yield stress and p is the density of the material. The boundary conditions are as follows: ~(0, t> = 0, and u(L, t) = u&l(t). The initial conditions are stress and strain free with a linear initial velocity: ~(0, X) = u,x/ L. A dimensionless variable for time was used; t= &,/L and in what follows in this section the overbar is dropped and t implies the dimensionless time. A bilinear curve for g(&‘“) was used; g(EVP)= w,, + hleVP for E”~~0.002 and g(&‘“) = m0+ O.O02h, + h,(evP - 0.002) for E“’ 2 0.002 with h, = 25a, and h, = -2.5a,,. A small imperfection in the yield function was assumed with the minimum value of the yield taken to

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

Fig. 1. Viscoplastic

slab subjected

be 0.995~~ at the center of the slab. The imperfection yield strength is expressed as q,f(x) where

to simple

337

124 (1995) 335-363

shear.

is triangular in shape; the spatial variation of the

f(x) = 1.0

O.Ocx/L

GO.25

f(x)=l.O-6(4X-1)

0.25~x/LS0.50

(2.lb)

f(x)=l.O-6(3-4x)

0.5OSxlL~O.75

(2.lc)

f(x) = 1.0

0.75 CxlL

(2.ld)

(2.la)

s 1.0

with 6 = 0.005. 2.2. Response with strain-softening The response is characterized by a rapid growth in the maximum plastic strain and a sharp drop in stress at the center of the slab, beginning at about t = 0.057 (see Fig. 2). If the time step is constant throughout the simulation, the growth in the maximum viscoplastic strain and drop in stress is accompanied by a rapid rise in the energy error, which causes numerical overflow (Fig. 3). The energy error is defined as follows

-

stress -+I-- peak

strain

t

t Fig. 2. Plot of peak-strain at the same time. Fig. 3. Energy drop.

error

and stress at the center

of the slab. Note that the stress drop and rise if the peak strain in the slab occur

with time. Note that the energy

error

rises sharply

as soon as the stress in the center

of the slab begins

to

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

338

124 (1995) 335-363

(2.2) whereW,“t>We,,, WKE,are the internal work, work done by external forces and the kinetic energy, respectively; they are evaluated as described by Belytschko [2]. The energy error can be controlled by using a variable time step, determined by the following procedure: the time step At is reduced in proportion to the maximum strain rate in the slab, so if

AfW,,,)” > a

(2.3a)

then (2.3b)

where4,,, is the maximum strain rate in the slab and (Y is a user defined parameter. For the above parameter values, the Runge-Kutta and central difference method give an energy conserving response for (Y= 1.0 X 10-4. The relative performance of the Runge-Kutta and central difference method for various constitutive updates is summarized in Table 1. Note that the forward Euler, Crank-Nicolson and backward Euler updates correspond to the Rate Tangent Modulus update [3] for 8 = 0, l/2 and 1, respectively. The Crank-Nicolson and backward Euler methods can compute the response in the stress drop for larger values of (Y and consequently need fewer time steps than the forward Euler update and the Runge-Kutta method. Here, N’“‘“’ and Ndrop are the total number of time steps and the number of steps required in the stress drop, respectively. The parameter (Ydepends on the exponent m and g(&‘“). The energy error for the various constitutive updates coupled with the central difference method is almost identical. The energy error at t = 0.059 for the various updates for (Y= 1.0 x 10m4 is summarized in Table 2. The Runge-Kutta method exhibits smaller energy errors than the central difference method, except in the stress drop. Fig. 4 shows a comparison between the energy error with the Runge-Kutta method and the central difference method coupled with the forward Euler update. Note that the rapid rise in the peak strain, the sharp drop in the stress and the rise in the energy error occur at the same time. The onset of the rapid rise in the peak strain is independent of shape of the imperfection, but depends on the minimum value of the yield stress in the imperfection. The spatial profiles of viscoplastic strain at various times depend on the spatial character of the imperfection [4]. Spatial profiles of the viscoplastic strain for the imperfection of Eq. (2.1) at various times are shown in Fig. 5. In contrast to this unstable behavior in the softening regime, the response in hardening is wellbehaved. Fig. 6 shows the energy error for both methods in the hardening regime. Note that the energy error decreases in the hardening regime. Table 1 Comparison of integrators for the strain-softening Runge-Kutta

Atmma, ;,,,..I N drcp

3.0 x 1om5 1.0 x 1o-4 4857 3137

slab

Central difference Forward Euler

Crank-Nicolson

Backward Euler

2.0 x 1o-5 6021 1.0 x 1o-4

2.0 x 1om5 3425x 1o-4 5.0

2.0 x 1om5 5.0 x 1o-4 3425

3391

1165

1165

Table 2 Energy error at t = 0.059 for the various constitutive updates for cx= 1.0 x lo-”

6 (t = 0.059)

Forward Euler

Crank-Nicolson

Backward Euler

7.51 x 1o-3

7.59 x 1o-3

7.27 x lo-’

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

o’ 0

339

124 (1995) 335-363

I 0.2

0.4

0.6

1

0.8

XIL

Fig. 4. Comparison of energy error for the central difference and Runge-Kutta

time integrators.

Fig. 5. Viscoplastic strain profiles in the slab at various times.

2.0

2.0

10”

0

0.2

0.6

0.4

0.8

P J

I 1 4

I

25.0

Fig. 6. Energy error in the hardening regime for the central difference and Runge-Kutta Fig. 7. Tip loaded strain-softening

2.3.

P

I c

integrators.

cantilever beam.

Strain-softening cantilever beam

As another example of a problem with material instabilities, results for a tip loaded cantilever beam are presented (see Fig. 7). An 8 x 10 mesh of four-node quadrilateral elements with one-point quadrature and hourglass control [5] was used. The relevant parameters are as follows: Young’s modulus = 5.0 x 104; elastic wave speed = 7.071 x 103; initial yield = 300.0; Poisson’s ratio = 0.45; load P = 200.0 H(t). The uniaxial stress-strain curve is bilinear and is described in Table 3. The material is governed by an isotropic von Mises 52 flow model in plane strain. A large deformation analysis was performed and the Jaumann rate was used as a frame indifferent stress rate. Energy error as defined in Eq. (2.2) was studied for both the hardening and softening regimes. Figs. 8-10 show the displacement, plastic strain at the root of the cantilever, and energy error versus time, respectively. For comparison, the energy error for hardening is shown in Fig. 11. Table 3 Uniaxial stress-strain

curve for the strain-softening

Effective plastic strain

Effective stress

0.0 0.20 20.0

300.0 95.0 10.0

cantilever beam

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

340

124 (1995) 335-363

0.5

1.5

I

time

Fig. 8. Tip displacement

2

2.5

3

time

with time for the tip loaded

Fig. 9. Viscoplastic

strain

strain-softening

at the root of the beam

cantilever

beam

shown

in Fig. 7.

of Fig. 7 with time.

10’

100

lQ

10-l

lo"

10-3

0.5

u

1

1.5

2

time Fig. 10. Energy Fig. 11. Comparison

of energy

2.5

U

3

4

e

10

time error

with time for the cantilever

error

for strain

softening

beam

and hardening

of Fig. 7. for the beam

of Fig. 7.

It is seen that in the softening regime, energy error grows rapidly when the effective plastic strain at the root of the cantilever reaches a certain point in the softening curve; within few time steps the energy error reaches unacceptable levels. In the hardening regime, the energy error diminishes with time.

3. Stability and exponentially growing solutions Many notions of stability have been presented in the literature, but most of these definitions relate to bounded solutions. The purpose of this section is to present some definitions in the context of both bounded and unbounded solutions; in the following sections these are examined in detail. Consider a solution y(t) to a differential equation jl +.@,

Y) = 0

(3.la) (3.lb)

Y(0) = Yo We will use the following definitions of stability [6]: DEFINITION.

Ill

y(t) is stable if for every E sufficiently small, there exists a 6 so that if

- Y(O)ll
(3.2a)

M. Kulkarni et al. I Compur. Methods Appl. Mech. Engrg.

341

123 (1995) 335-363

then

IIJW-

Y@)II < E >

(3.2b)

Vt

where j(t) is the solution to (3.1) with initial data j(0). DEFINITION.

Ill

y(t) is relatively

-Y(O)ll

stable if for every E sufficiently small, there exists a S so that if (3.3a)


then

lb-w-Y(M

< EIlY(Qll I

(3.3b)

Vt

A stable solution does not necessarily satisfy the conditions of relative stability as defined in Eq. (3.3a)-(3.3b). T o circumvent this, the following definition of generalized relative stability is introduced: this is termed g-rel stability. DEFINITION.

y(t) is g-rel stable if for every E sufficiently small, there exists a 6 so that if

Il~m- Y(O)ll

(3.4a)

< 6

then

IIW -Y(Qll

<&(llY(9ll

+ IlJT~)II + I)?

Vt

(3.4b)

This definition also symmetrizes the concept of relative stability. In applications it is not always clear which solution represents the true solution and which solution represents the perturbed solution. With this definition, the difference of the two solutions is required to be small relative to both solutions. Also, there may be cases, e.g. strain-softening viscoplasticity, where two solutions j(l) and y(t) can start growing at different times. In this case the behavior of the relative difference can differ dramatically depending on which solution is used in the denominator. A nice feature of g-rel stability is that a stable solution is always g-rel stable, whereas it is not always relatively stable. While the above definitions of stability are given in terms of initial conditions, stability can also be defined in terms of perturbations of coefficients of the solutions. The stability of numerical techniques for the solution of ordinary differential equations is usually studied in terms of (3.2), but for the purpose of analysis, the equations are linearized. As a consequence, the perturbed solution is governed by the same linearized equation. To develop conditions for stability, a characteristic equation is developed with roots ri. The solutions to the ordinary differential equations at time step n are ry. Obviously, ri must lie within the unit circle for the solutions to be bounded. This leads to the definition of absolute stability. A method is said to be absolutely stable for a continued given time step if all roots of the stability equation lie within the unit circle, i.e. Ir,I < 1 for all i.

In studying the stability of partial differential equations, the notion of von Neumann stability is often used. This is closely related to absolute stability. For example, in finite difference methods, in a von Neumann stability analysis, a Fourier transform is applied in space to the partial differential equations, which leads to a set of ordinary differential equations which are then checked for stability (i.e. lril < 1). Similarly, in finite element methods, the ordinary differential equations resulting from the semidiscretization are required to have characteristic solutions Ir,l s 1. Thus, in these methods of analysis, von Neumann stability is closely related to absolute stability. The properties of numerical approximations of unbounded solutions are described next. Clearly, many of the definitions of numerical stability are not applicable when the true solution (and the numerical approximation) are unbounded; when a differential equation has an unbounded solution, then the numerical solution must also be unbounded, so lril > 1 for some i and the solution will not be absolutely stable. Therefore a broader concept of stability (or acceptability of numerical solutions) is needed for unbounded solutions.

342

M. Kulkarni et al. I Comput. Methods Appl. Mech.

Engrg.

124 (1995)

335-363

Consider a discretization, w,, corresponding to a grid size At. Proceeding by analogy with the concept of relative stability, the relative error may be defined as

WAI

K = IIY, n

(3.5a)

IlYA

and the asymptotic relative error as follows K

Y

IlYn - %II

=lim

12-=

(3.5b)

llYnII

where y, = y(n At) and w, is the numerical solution. The discretization there exists a k such that

is relatively stable if for At d S,

IlYn - %I1 ~~llY”ll~>O and the discretization

(3.6)

w, is g-rel stable if for At G 6, there exists a k such that

IIY, - w,II WIYnII

+ l),n>O

Consider relative stability for systems of differential equations wheref(t, of the form

(3.7) y) is linear in y, i.e. systems

j+Ay=O

(3.8a)

Y(0) = Y0

(3.8b)

where A is an m x m matrix, assumed to be real symmetric and y is an m vector. A superposed dot implies a derivative with respect to the independent variable (this notation will be adhered to throughout this paper). When the matrix A has both positive and negative eigenvalues, the system (3.7) is said to possess an exponential dichotomy [7]. Let Ai be the eigenvalues of A, with A, < A, < . . . < A,, and A, CO. The effect of initial conditions on the relative stability of solutions to Eqs. (3.7) is now examined. Let m Y0 = YPX, (3.9)

c

i=l

where xi are the orthogonal eigenvectors of the matrix A and yp are the components basis. The solution to the system of Eqs. (3.7) is then

y(t) =

L$, YP e-*k

of y” with x, as a

(3.10)

Consider two initial conditions y” and j” such that yp=Yp,

i=2,.

(3.11a)

. . ,m

y; =o

(3.11b)

YY=6

(3.11c)

so that IIY” -BOIL

(3.12)

= 8

The solutions corresponding

to these initial conditions are

Y(t) = *$ Yp e-*“Xi

(3.13a)

,zlyP e-*“xi

(3.13b)

f(t) =

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

343

Thus, for sufficiently large t, it follows that IIy(r) -Ml, Ilu(W

= 6 e-*I’

= l~ile-~*’

IIYW -mllm IIYW

8 pt =Ir’:ie

(3.14a) (3.14b) (3.14c)

where 11.I(_ is the L, norm and /3 = A, - A,. Note that since A, < A,, p > 0 and e” grows exponentially. Thus, for this choice of initial conditions, the solution y(t) cannot satisfy Eqs. (3.2) or (3.3) and is consequently relatively unstable. Note that when yr # 0 IluWlL = 1~3 emAl’ IIYW - WL IIYW x

=- 6 Iv3

and the solution y(t) is relatively stable. Thus, solutions with non-zero components growth direction are relatively stable.

(3.15a) (3.15b) along the fastest

4. Analysis for first-order systems In this section, various types of stability are examined for forward Euler and Runge-Kutta integration of first-order systems. Although much of what follows is available in the literature, it is repeated here to give a comprehensive picture. Consider systems of the form i+Ay=O

(4.la)

Y(0) = Yl,

(4.lb)

where A is real symmetric. The above system of equations can be uncoupled by expanding y in terms of the eigenvectors of A to yield scalar equations of the form Y+Ay=O

(4.2a)

Y(0) = Yn

(4.2b)

where A is an eigenvalue of A. In one-dimensional strain-softening viscoplasticity, at least one eigenvalue is negative (see Section 6). It will be shown that (1) If the eigenvalue A in Eq. (2a) is negative, then Euler integration is unstable but relatively stable regardless of the size of the time step; however, the relative error grows unacceptably large (tends to unity). (2) If A > 0, then Euler integration is stable and g-rel stable if and only if the time step is smaller than a critical time step, i.e. if At < 2/A. Thus, when A < 0, Euler integration is g-rel stable regardless of the size of the time step and when A > 0, Euler integration is g-rel stable when At < 2/A. These results will then be generalized to the Runge-Kutta method and will be verified by numerical results. 4.1. Stability of the first-order equation j + Ay = 0 The closed form solution to Eqs. (4.2) is given by Y(t) = Y,,e-*’

(4.3)

The notion of stability defined in Eqs. (3.2a)-(3.2b) can be applied to the solution (4.3) as follows; consider two initial conditions y, and Y, and a small parameter 6 such that

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

344

124 (1995) 335-363

(4.4)

90= Yo + 6 and let y(t) and Y(L)be the solutions corresponding exists an R such that

to y, and y0. Then the solution y(t) is stable if there

Iv(t) - J?l(t>I=G& Vt

(4.5)

The stability of the solution (4.3) depends on the sign of the exponent A. The difference and j-(t) as a function of time is given by

between y(t)

p(t)=(y,+S)e-“’ y(t) -y(t)

(4.6a)

= S ep*’

(4.6b)

Note that when A 0, the difference in the two solutions decays exponentially and Eq. (4.5) is satisfied for F = 6. Thus, the solution y(t) defined in Eq. (4.3) is unstable in the sense of Eqs. (3.2a)-(3.2b) for A < 0 and stable for A > 0. The solution y(t) is relatively stable irrespective of the sign of A as

IN>-

Y@)l =

$ IY(4

4.2. Relative error and stability for the forward Euler method The forward Euler approximation

to Eq. (4.2a) is given by (4.7a)

Y,,+,=y,+j,At=y,(l-AAt) The amplification

factor for the forward Euler method (i.e. the ratio y,,

I /y,)

is (4.7b)

CY FE=(l-AAf) The solution to the difference

is given by recursively applying Eq. (4.7a)

y, = y,,( 1 - A At)”

(4.8)

For absolute stability, 1~FEI< 1, which yields the conditions A > 0 and At < 2/A. The closed form amplification factor may be defined analogously as the ratio of the closed form solution at time t + At to the closed form solution at time t. The closed form amplification factor is given bY a

CL

-h At

(4.9)

=e

For the forward Euler integration Yo K,

=

e-AaAt-

of Eq. (4.2), the relative error defined in Eq. (3.5a) is

yo( 1 - A At)” -An

Ar

=l-(‘$,3’)”

Yo e When A < 0, (4.11a)

l-AAt
Fig. 12 shows the behavior of e-* ” and (1 - A At) with A At for -1 s A At s 1. It can be seen that < 1 for all At > 0

(4.11b)

It follows that, as n+ ~0, ((1 - A At)/(e-” “‘))” + 0 and K, = 1. The growth of relative error with time for the case A < 0 for the forward Euler and Runge-Kutta integrators is shown in Fig. 13 for A = +O.2, At=O.l. When A > 0, the relative error may be written as

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

345

124 (1995) 335-363

IO0

-0.5

-1

10.101

-0.6

-0.2

0.2

0.6

0

I

7

20

40

60

80

I00

t

)iat Fig. 12. Behavior Fig. 13. Relative

error

of em””

and l-AAI

for the forward

Euler

with AAr in the interval and Runge-Kutta

-l
time integrators

for Eq. (4.la).

(4.12)

~,=1-(1-AAt)“e^“~’

Thus, when At > Atcrit, where Atcrit = 2/h and A > 0, (1 - A At) < - 1 and consequently the relative error grows unboundedly. This corresponds to the situation where the absolute stability condition is violated. The relative error is bounded when 11- A AtI 0, there is a transition between bounded and unbounded relative error depending on the ratio I(1 - A At)1/eA” “, which defines the critical time step Atcrit (often called the stable time step). For A < 0, however, there is no sharp transition between boundedness and exponential growth and the asymptotic relative error equals unity irrespective of the time step used. For larger values of the time step, this asymptote is reached faster. This is an important feature of the forward Euler method for A < 0; the relative error is at most unity regardless of the time step. Consequently, according to the definition of Eq. (3.5a) the method is relatively stable for any time step for A < 0. For A > 0, the relative error error is unbounded when A At > 1.278, so in the regime 1.278 < A At < 2, the forward Euler method is absolutely stable, but not relatively stable. This occurs because the error of the numerical solution becomes large compared to the magnitude of the decaying exact solution. However, according to the definition of g-rel stability, Euler integration for A > 0 is g-rel stable in the same regime that it possesses absolute stability. The conclusions for the forward Euler method for A < 0 are equally valid for other one-step integration schemes. The relative error for one-step schemes is obtained from Eq. (4.12) with (1 - A At) replaced by the amplification factor of the integration method aNUM a K,=l-

-

( e

NUM -A

Ar

n

(4.13)

)

For an explicit forth-order

Runge-Kutta

scheme, the amplification factor is

ff RK=l-AAt+;AAtZ-;AAt3+&AAt4 i.e. aRK agrees with a Taylor series expansion for e-*” up to fourth order in At. When A < 0, aRK < e-h ‘I so, using the same arguments as for the forward Euler, the asymptotic relative error will be unity. The amplification factor of a higher-order method agrees with a Taylor series expansion for ePh ” to a greater number of terms than the first-order accurate forward Euler method. Consequently ((a N”“) /(e-” Ar))n+ 0 s1ower and the relative error is improved (Fig. 13). It can be seen that for A < 0, neither the forward Euler method nor the Runge-Kutta method exhibits the unbounded growth relative to the exact solution which characterizes the numerical instability which occurs when At > 2/A and A > 0. Instead, the major difficulty encountered with the forward Euler method is that the relative error becomes unacceptably large. With higher-order integrators, such as the

346

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

Runge-Kutta method, the relative error is much smaller, so these integrators improve the results for this class of problems. Note that absolute stability requires that A > 0, so an analysis for absolute stability would indicate that they cannot be used for non-Drucker materials. However, the relative error remains bounded for non-Drucker materials and the numerical methods are in fact acceptable in this regime; the major difficulty is the rapid growth of error in low order schemes such as forward Euler. The boundedness of the relative error for the forward Euler for A < 0 has been previously noted by Quinney [8, p. 761.

5. Analysis of second-order systems with exponential dichotomy In this section, second-order

systems of the form

Mj+Ky=O

(5.la)

are considered with the restriction that M and K are real symmetric, and M is positive definite. Exploiting the symmetry properties of M and K (and, consequently the orthogonality of the eigenvectors of these matrices) the system of Eqs. (5.la) can be uncoupled to yield scalar equations of the form j; +/A,y=o

(5.lb)

where hi satisfy the following eigenvalue problem Kx = AMx

(5.2)

Attention is restricted to a system with an exponential dichotomy where at least one eigenvalue is negative and at least one is positive. The independent variable (time) is normalized as follows, t= t/(-A;)l’z where Ai is the smallest negative eigenvalue. Dropping the overbar for simplicity, Eq. (5.lb) may be written as jj=y

(5.3a)

Y(0) = Y,

(5.3b)

Y(0) = P,

(5.3c)

where superposed dots indicate derivatives with respect to & the bars will be omitted henceforth for simplicity. This second-order equation may be recast as a first-order system, see Eq. (5.19). It will be subsequently shown that this system has eigenvalues of +l and -1, i.e. it possesses an exponential dichotomy. Eqs. (5.3) also describe a one-dimensional oscillator with a negative stiffness. If the coefficients are frozen at the onset of strain-softening, this equation provides a useful model for strain-softening. This problem was previously studied by Hildebrand [9]. Conditions under which the relative error is bounded are now evaluated for the central difference and Runge-Kutta methods. The closed form solution for this equation is given by

Consider initial conditions such that Y, + Yo = P

(5.5)

The solution (5.4) for initial conditions with /3 = 0 is a decaying solution. For small values of /3, the solution exhibits a delayed growth. This is a consequence of the fact that there is a small quantity multiplying an exponentially growing term and the exponential term does not dominate the solution from the beginning. This feature mimics the behavior of the viscoplastic strain-rate in strain-softening, see Section 6. A slight perturbation of either Y. or y, will result in the appearance of the exponentially growing term in (5.4) causing large changes in the solution y(t). Thus, the solution (5.4) with /3 = 0 is an

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

(4

124 (1995) 335-363

0.35

347

1

p =0.5x IO”

0.3

,,.----:

0.25

0.1 0.05 0‘

.

0

4

2

6

8

10

12

14

t @‘I

2.0 Id

0.0 IO0

-2.0 Id :*fi: -4.0 ld

: : : : : :

-6.0 ld

-8.0 Id * 0

5

I5

10

20

25

t

Fig. 14. (a) Behavior of the g-rel stability parameter for Eq. (5.3); (b) Behavior of the central difference and Runge-Kutta schemesfor Eq. (5.3). Note that the central difference approximation diverges as soon as the exact solution begins to grow.

unstable solution in the sense of Eqs. (3.2). The solution with p = 0 is also relatively unstable according to condition (3.3b). For /3 #O, the solution (5.4) is unstable but relatively stable. Fig. 14(a) shows the behavior of the relative difference tREL = ([Y(l) - y(l)l)/( ly(t)l + Iy(l(t>l) of the solution y(t) corresponding to /3 = lo-’ and solutions Y(c) corresponding to p = 5.0 x lo-‘, 6.0 x lo-* and 7.0 X lo-‘, respectively. It can be seen that after a sharp rise, the relative difference asymptotically reaches the relative difference in /3. This sharp rise in the relative difference arises because the solution y(t) starts to grow earlier than the solutions jr(~). Thus, for small /3, the solution of Eq. (5.4) is relatively stable. In contrast to the behavior of SREL, the difference between the two solutions, 5 = &((t) - y(t)1 grows exponentially. Fig. 14(b) shows the behavior of the Runge-Kutta and central difference integrators for initial conditions p = lo-‘, i.e. a slight perturbation of Eq. (5.5) with p = 0. Observe that as soon as the closed form solution begins to grow, the central difference scheme diverges from the closed form solution, whereas the fourth-order Runge-Kutta gives acceptable results. An analysis of the relative error defined in Eq. (3.5a) sheds some light on the behavior of these integrators. It is shown that the relative error for a central difference approximation to Eq. (5.3a) for initial conditions for B = 0 or small fi tends to be unacceptably large, even for very small values of the time increment. For the Runge-Kutta method, relative errors are within acceptable limits. 5.1. Central difference method When the second derivative in Eq. (5.3a) is replaced by a centered finite difference the difference approximation to Eq. (5.3a) is given by Yn+, -(At*+2)y,,+~,z-,

=O

approximation,

(5.6)

348

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

The solution to this difference

124 (1995) 335-363

equation is given by

y, = c, A; + C,AY

(5.7)

where A, and A, are roots of the equation A2 - (At’ + 2)h + 1 = 0

(5.8)

The roots of the above quadratic are given by A r~,=l+$kAt

For At-O,

(5.9a)

the square root in the above equation may be approximated

by

(5.9b) The two roots of the quadratic equation (5.8) may then be written as a series in At, which agree with the expansion for eA’ and emA’, respectively, to second order in At. Thus * At3 A,=l+At+~+~=e”‘+~(bt’)

(5,lOa)

2 At3 h,=l-At+~-s~e-“+o(A~‘)

(5.10b)

Substituting

for A, and A, in Eq. (5.7) the solution to the difference

y,* = C,

1-f At + $

equation may be written as

+ O(n At’)

(5.11a)

This may also be written as y, = C,(eA’)” + C2(e-Ar)n + O(n At’)

(5.11b)

The constants C, and C, are determined from the consistency of the initial conditions. A supplementary starting value yr is assumed and consistency is enforced by requiring that the error in the initial value of the derivative go to zero as the time increment goes to zero, i.e. let Yl -Yo ~ = Y. + At

(r

(At)

(5.12)

where a(Af) is the error in approximating the initial derivative, conditions (5.3b)-(5.3c) and Eq. (5.11a) give

which is first order in At. The initial

y, = c, + c 2

(5.13) 2

yr=C,

l+A,+-$

+

O(Alc’)

b

Eliminating the supplementary starting value y1 and expressing C, and C, in terms of yr and am, substituting in Eq. (5.11b), the difference solution may be written as Y, =i((y,

+Yo) -$

(5.14)

(5.15)

and

yo+a(At))e”A’+~((yo-jo)+~yo-a(At))e~”A’+O(nAt3) (5.16)

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

349

124 (1995) 335-363

Note that the coefficients multiplying the exponential terms in Eq. (5.16) differ from the closed form solution by the term -y, At/2 + a(At) which tends to zero as At-to. Expanding y, in a Taylor series expansion about t = 0, it can be shown that -y, At/2 + a(At) = i,, At*/6. Thus, the coefficients of the central difference method differ from those of the closed form solution by a second-order term in At. The relative error as defined by Eq. (3Sa) for this problem is given by n AI

e

K,

=

4At) > /? e” ‘*

~(1-

+

The behavior of two cases:

-e

-n

Al +

(5.17)

O(n At’)

+ (yO - %,) e-’ ”

as a function of time depends on the initial conditions and will be subdivided into

K,

Case 1: p=O $ Kn

=

y, (y,

_

a(At) j.)

(e*”

AI -

1)

The relative error grows exponentially unstable.

+

(5.18a)

O(n At”)

in this case, so the numerical solution is unstable and relatively

Case 2: /3 # 0 Since edfl ” << e” A’ and

(5.18b) Since the numerator and denominator unacceptably large for small p. 5.2. Runge-Kutta

are constants,

the relative

error

is bounded,

but may be

method

In analyzing error growth for the Runge-Kutta first-order system

method,

Eq. (5.3a) is treated

as an equivalent

(5.19a) subject to the initial conditions (5.19b) The closed form solution for the system of Eqs. (5.19) is given by {;;;;}=;[;y; The amplification

(5.20)

$;$]{;:J matrix for the Runge-Kutta

method is given by Lambert

[lo]

G RK=I+AtA+;A,2.2++At3A3+&At4A4

(5.21a)

where (5.21b) Substituting Eq. (5.21b) in Eq. (5.21a), the amplification matrix for the Runge-Kutta written as

method may be

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

350

G Rx =

1 1 1 +-i_At2+zAt4

I

124 (1995) 335-363

At+;Al’

At +$it3

1 +~At2+zAt4

The recursive relation for the Runge-Kutta 1+~At2+~At4 1 1

(521c)

1

1

I

may then be written as

At+;At”

Af+;At3

1 +TAt’ 1

(5.22)

+zAt4 1

The eigenvalues of A are + 1 and - 1. The corresponding eigenvalues of the amplification Eq. (5.21a) and the spectral mapping theorem [ll, p. 3501, are

matrix from

(5.23a) A2=1-At+;At2-;At3++4At4 and the eigenvectors

(5.23b)

are given by

+T= [l,

+11

(5.24a)

+T=[l,

-11

(5.24b)

The recursive relation (5.22) can therefore =P diag[h;

be written as

hi]P-’

(5.25a)

where p=

1 [1

1 -1

1

(5.25b)

is the matrix of eigenvectors of A. Carrying out the multiplication recursive relation (5.22) may be written as A;+h; A;-A;

Al-h; A?+A;

From Eq. (5.26) the Runge-Kutta

in Eq. (5.25a) the solution of the

I{I y, Y,

approximation

(5.26) to y(n At) is given by (5.27)

A;+;(y,-&)A;

Note that the eigenvalues A, and A, are series expansions for ear and ePA’ accurate to fourth order in At, i.e. A, = eAr + O(At’)

(5.28a)

A, = e-*’ + O(At”)

(5.28b)

The Runge-Kutta y,=+p

approximation

en At +$(y,-YO)e-“*’

to the closed form solution is + O(n At”)

(5.29)

Unlike in the central difference method, see Eq. (5.16), the coefficients of the exponential terms are exact and consequently, errors arise only from truncation of the expansions for the exponential terms. Numerical results for three values of p for the relative error of the central difference and

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

351

124 (1995) 335-363

10”

10’2 10’0 IO’ &06 IO’ IO’ 100 lo”

lo”’

lOA

t Fig. 15. Behavior growth of relative

Table 4 Relative error Kutta method

integrator

integrator

applied

applied

to Eq. (5.3) for various

to Eq. (5.3) for various

for the central differencemethod and Rungeapplied

P

Central

1o-x 10mh 10-’ 1om4 1om3 lo-’ 10-l 1.0

66.613 0.6626 0.06267 0.00267 0.00332 0.00392 0.00398 0.00399

to (5.4)-(5.5) difference

5

IO

IS

20

i

t

of central difference error for /3 = 0.0.

Fig. 16. Behavior of Runge-Kutta is accurate for all values of p.

1

0

with At = 0.004

values

Table 5 Relative error for the central (5.4)-(4.5) for /3 = lo-’ At

K,

1.87 1.89 1.77 1.66 1.55 1.55 1.55 1.55

0.001 0.002 0.003 0.004 0.005

0.4149 1.664 3.747 6.6625 10.414

lo-* lo-“’ 1o-‘O 1o-‘O 10-l” lo-lo lo-lo 1o-‘O

is an exponential

values of p. Note that the Runge-Kutta

Runge-Kutta x X x x x x x x

of @ Note that there

difference

method

integrator

applied

to

Runge-Kutta integration are shown in Figs. 15 and 16, respectively. The parameters used were At = 0.004, and j. = -1.0 and y, = 1.0 + p. When /3 is small, the relative error of the central difference method tends to a relatively large constant. This constant increases as p decreases. When /? = 10e7, the asymptotic relative error K, = 6.0, when /3 = lo-*, K, = 66.0. When p = 2.0, the relative error for the central difference method was constant at lo-*. The relative error for the Runge-Kutta is much smaller than that for the central difference method. The variation of K, for the central difference method for different values of /3 with At = 0.004 is summarized in Table 4. Thus, for growing solutions, the Runge-Kutta method is characterized by smaller relative errors just like in the absolutely stable regimes. The central difference method tends to exhibit a larger asymptotic error for larger values of the time step for fixed values of /3. The effect of the size of the time step on the relative error in the central difference method is summarized in Table 5 for /3 = 10e7.

6. Analysis of viscoplasticity for central difference and Runge-Kutta time integrators In this section, the central difference method and Runge-Kutta method are analyzed in the softening regime for a one-dimensional semi-discretization. The term central difference method in this section refers to central difference integration of equations of motion with forward Euler integration of the constitutive equations. The field equations of one-dimensional viscoplasticity are recast as a system of first-order ODE’s in time. A closed form solution of this system of equations is then obtained in terms of matrix

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

352

124 (1995) 335-363

exponentials of the coefficient matrix. The amplification matrices for the Runge-Kutta and central difference methods are obtained and compared with the closed form amplification matrix. The method of frozen coefficients is used, so the amplification matrices are assumed to be constant during a time increment. Unlike the case of linear first- and second-order systems, the amplification matrices are not real symmetric. The softening regime is shown to admit growing solutions, and the hardening regime is shown to admit only decaying solutions. Both methods are shown to be absolutely unstable in the softening regime. Using the method of frozen coefficients, the behavior of the solutions in the softening regime may be modelled by the first- and second-order equations dealt with in Section 4 and Section 5. It has been shown that numerical integration of these equations is relatively stable. It will be shown that the amplification matrix of the Runge-Kutta method agrees with that of the closed form amplification matrix up to fourth order in At, whereas the central difference method has an erroneous second-order term for both hardening and softening. Consequently, the central difference method is only first-order accurate and its eigenvectors differ from those of the closed form amplification matrix. The Runge-Kutta method has the same eigenvectors as the coefficient matrix of the original first-order system and errors arise from truncation alone. 6.1.

Field equations of one-dimensional

The field equations for one-dimensional momentum

dynamic viscoplasticity

viscoplasticity

are as follows

equation

O-,=pti

(6.1)

viscoplastic constitutive

equation

b = m,,- +(a,7))

(6.2)

I bbldt

(6.3a)

4 = $Jsgn (T

(6.3b)

Y=

and

Here, m denotes the stress, 4 the viscoplastic strain rate, u is the velocity, E the elastic modulus and A superposed dot represents a material time derivative and a comma represents a derivative with respect to the variable which follows. These equations describe a one-dimensional bar in uniaxial stress or simple shearing of a viscoplastic slab. The strain rate 4 depends on the sign of the hardening modulus. This plays an important role in error propagation for the central difference method. Consider a typical functional form for the strain rate [12] p the density of the material.

H

l/RI

( >

4(a,r> =d go

sign(m)

where b and m are material parameters.

The derivatives of +(u, y) are

(6.4~)

[r$ sign(a)],, where

= - ’ zgT$)h

(6.4d)

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

h=F

124 (1995) 335-363

ds

353

(6.4e)

is the plastic modulus. Note that in the softening regime, h < 0 so [C#J sign(a)] ,y > 0, whereas in the hardening regime h > 0, so [C#Jsign(o)] ,y < 0. The viscoplastic strain rate c,t~increases in magnitude with increasing viscoplastic strain in the softening regime. 6.2.

Semidiscrete

equations

In what follows, attention is restricted to a program of loading where u > 0, and 4 > 0. The differential equations (6.1)-(6.3~) are now discretized in space. The spatial derivatives in the governing equations (6.1)-(6.3) are replaced by centered finite difference approximations, as follows ,j

=

uj+l/2

-

.x

uj-l12

Ax

vj.x =

v,+1/2

-

‘j-112

(6.6)

Ax

Assume the following form for the field variables Uj = d(t)/_2

(6.7)

vi = v)(t)/_bj

(6.8)

4’ = @(t)jJ

(6.9) (6.10)

(6.11) where p = e”” ‘I, i = a, m is the wave number and (Y= 2 sin@ AX/~). Substituting these equations in the governing equations gives (6.12) (6.13) (6.14) variables are used: V = v’/c, ~7= a/E, i= tc/Ax, 4 = 4’ Ax/c, f, = 4,,, Ax/c, f, = E4,, Ax/c; c = (E/P)“~ is the elastic wave speed. Substituting these in the semidiscrete equations and dropping the overbars, the linearized equations may be written as The following dimensionless

Y,, + AY = 0

(6.15a)

YT = [v, (+741

(6.15b)

where

and A=

0 -ia ( -iafU

-iff 0 0

0 (6.15~) -(&f,,

)

A time interval t c t* 6 t + At is considered

and the error per step is evaluated

at time t + AC. The

354

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg. 124 (1995) 335-363

matrix A is frozen at time t and is assumed constant in the interval t 6 t* d t + At. The solution for the above equation within the time step is [13] y(t*> = e-Nr*V y(t)

(6.16a)

tdt*st+At

where eWA’is the matrix with the same eigenvectors as A, and eigenvalues given by eeAi’ where A, are the eigenvalues of A. It is assumed that A has distinct eigenvalues and linearly independent eigenvectors; the matrix eeA’ is then given by e -A’ = p diag[e-“I’

e-*2’

e-“3’]~-’

(6.16b)

where P is the matrix whose columns are the eigenvectors of A. Thus, the values of the field variables at t + At can be written in terms of the values at t as follows y(t f At) = eeAA’ y(t)

The amplification G

(6.17)

matrix CCL = emA*’ may be expanded in an infinite series as follows [13]

cL=Z-AfA+~Az-~A3+~A4+O(Afi)

(6.18)

The solutions to the system of the ordinary differential equations (6.15) decay if all the eigenvalues of of A satisfy the following characteristic equation

A have positive real parts (see [14]). The eigenvalues

ef + (f, -f,)e’

(6.19a)

+ Ly2ei+ u*f, = 0

Let e; = -e;

(6.19b)

The t?; satisfy (6.19~)

e3-(fy-f,)ef+ry2e,-a%=0

The Descartes rule of signs is used to obtain the signs of the roots of the polynomial (6.19~); ‘The number of positive real roots of an equation with real coefficients is either equal to the number of variations in sign or is less than that number by a positive even integer. A root of multiplicity m is counted as m roots’ [15]. Note that f, > 0, f, > 0 for softening (6.4a)-(6.4e). If f, >f,, then Eq. (6.19~) has three variations of sign, thus there are either three real positive roots or one real positive root. If f,
matrix for the Runge-Kutta 3

2

method

algorithm can be written in terms of the matrix A as

4

G RK=Z-AtA+$A2_$A3+$A4

(6.20)

where Z is the identity matrix. Note that the amplification matrix for the Runge-Kutta method is a matrix polynomial of the matrix A and consequently has the same eigenvectors as A [13, 141. Therefore, the eigenvalues of GRK, hi are related to the eigenvalues of A, ei, by the following 2

Aj=l-Atei+~e~-~ej+~e~

3

4

(6.21)

Since at least one eigenvalue of A is real and negative in the softening regime, the corresponding eigenvalue for G RK will lie outside the unit circle. Thus, from the viewpoint of absolute stability theory, the Runge-Kutta method is unconditionally unstable in softening.

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

6.4.

Amplification

matrix for the central difference

124 (1995) 335-363

355

method

Details of the derivation of the amplification matrix are given in Appendix A, Eqs. (Al)-(A25). Normalizing and ordering the solution vector as in Eq. (6.15), and setting the integration parameter 8 = 0, the central difference amplification matrix derived in Appendix A, may be written as

G

CD

=

l--‘At* ia At ( iafWAt

icrht 1 0

The above amplification GcD=(;

;

icu At2 -At 1+ Ar(f, -f,) 1

(6.22)

matrix may be written as a series in Ai

8)+At(iP

f

(!$,)+$(2;

i

-;&)

(6.23)

A Neumann analysis for the central difference method is performed in Appendix A and indicates that at least one eigenvalue lies outside the unit circle for the softening regime. Necessary conditions for numerical stability in the hardening regime are obtained as a natural consequence of the Neumann analysis. The series expansion for the closed form amplification matrix derived in Eq. (6.18) is compared with that of the Runge-Kutta method equation (6.20) and the central difference method equation (6.23). The Runge-Kutta method truncates the series after four terms and the amplification matrix has the same eigenvectors as the closed form amplification matrix (GRK can be expressed as a matrix polynomial of A). Thus, errors in the Runge-Kutta method arise from truncation alone and can be made sufficiently small provided sufficiently small time increments are used. The amplification matrix for the central difference method agrees with the closed form amplification matrix up to first order in At, but has an erroneous second-order term. Conse uently, this amplification %D matrix has eigenvectors that differ from those of the coefficient matrix A. (G cannot be expressed a matrix polynomial of the matrix A.) To obtain an estimate of the error in the incremental values of the solution vector, the difference between the closed form incremental solution and the incremental solution for the central difference method is examined. The solution vector at time t + At, y(t + At) may be written in terms of y(t) as (6.24)

y(t + At) = G CDy(t)

Using the series expansion for the closed form amplification matrix and the central difference amplification matrix and retaining terms up to second order in At, the error term may be written as yerr(t + At) = $

(A* - B)y(t) + O(At’)

(6.25)

where (6.26) The error vector may be written as a

2

y”“(t + At) = +

2

-i&f, iaf,(f,

-L)

For the power law model equation

0 -a

-a’f,

2

(6.27) (f, -f,)’

(6.4a) the term (f, -f,)

is given by (6.28)

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

356

124 (1995) 335-363

In real materials E/a S- h/g; moreover in the softening regime (+ decreases and E is constant, so the ratio E/u increases and 4 increases with viscoplastic strain according to Eqs. (6.4a)-(6.4e). In the hardening regime, the ratio E/c+ decreases for constant and 4 decreases with viscoplastic strain. Substituting Eq. (6.28) in to the last row of Eq. (6.27) it can be seen that the error in 4 is of order 4” At*. Consequently, the error increases in the softening regime and decreases in the hardening regime. Although the above conclusions are based on a linearized analysis, they are borne out by deterioration in the performance of the numerical method when the strain rates become very large and decreasing energy error in the hardening regime. It is nevertheless puzzling as to why the time step needs to be reduced so much in the latter stages of the simulations, for although the linearized analysis indicate a growth of error with increasing strains, it indicates nothing pathological. A family of semi-implicit constitutive updates is given by Peirce et al. [3] (6.29) where 0 < 8 G 1. The forward Euler, Crank-Nicolson and backward Euler updates correspond to 8 = 0, l/2 and 1, respectively. From Eqs. (6.4b)-(6.4c) the strain rate 4 and the derivatives 4,, and 4,,, grow very rapidly at the onset of the localization, which is accompanied by a sharp drop in stress. As a consequence, the time step must be reduced progressively so that the stress drop may be tracked till the stress vanishes in the slab. In order to examine the various definitions of stability in the context of numerical solutions of the field equations of viscoplasticity, the behavior of the solutions to the slab problem were considered. These solutions were effected using the central difference method coupled with the backward Euler constitutive update (0 = 1) with various values of 6 in Eq. (2.1). Stability was examined by evaluating L, norms of the viscoplastic strain rate 4 and examining their behavior with time. For this purpose, the difference norm and relative difference norm defined in Section 4 are generalized as follows

5(t)= Ilf, -

(6.30a)

fAl2

where fi is the vector of values of $+ in the elements. The relative difference

IlfI SREL(t)

=

IMl2

norm is defined as by

f2112

(6.30b)

+ llf2112

The solutions f, and f2 are stable if s(t) satisfies Eqs. (3.3) and are g-rel stable if they satisfy Eqs. (3.4), the constant (=l) that appears in Eqs. (3.3) plays no role in this particular problem and hence is dropped from the definitions (6.30). In all cases, f, corresponds to 6 = 5.1 X 10e3. t(t) and SREL(t) were

t

I

Fig. 17.

e(t) for S = 5.1 x 10m3 for the field equations Fig. 18. (,,,

with time for various

of one-dimensional values

of 6.

viscoplasticity.

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

357

evaluated with f2 corresponding to 6 = 5.075 x lo-‘, 5.05 x 10-3, 5.025 X 10e3 and 5.0 X 10m3, respectively. In all cases, the difference norm grows unboundedly in the last stages of the computation. Fig. 17 shows the behavior of l(t) with time for f, corresponding to S = 5.1 X 10-3. Fig. 18 shows the behavior of the relative difference norm eRri_ with time. tREL(t) shows an initial sharp rise and then decreases. Furthermore, the maximum value of tREL(t) decreases linearly with difference in 6. We can conclude that the solutions are g-rel stable but not stable. Note that this behavior is similar to that of the second-order system with an exponential dichotomy, Fig. 14(a).

7. Conclusions The evaluation of the effectiveness of various numerical integrators and the selection of a time step poses a challenging problem because the usual notion of von Neumann stability is not applicable. This is a consequence of the fact that the ordinary differential equations that result from semi-discretization admit unbounded solutions. The numerical solutions are therefore absolutely unstable. A von Neumann analysis was made of the field equations of one-dimensional viscoplasticity using the method of frozen coefficients. It was shown that this analysis predicts an unconditional instability in the softening regime. Usually an instability predicted by a von Neumann analysis is a characteristic of the numerical method and manifests itself in a spurious growth of the numerical solution which is large compared to the solution of the differential equation. In strain-softening viscoplasticity however, the von Neumann instability is a manifestation of the instability of the underlying equations, the growth of the numerical solution can be said to be inherited from the partial differential equations. The numerical solution in fact grows at a rate similar to the growth of the exact solutions. Another way to look at this dilemma is that the von Neumann analysis employs the method of frozen coefficients, which in general is not a good model for long-term behavior of the solution. In order to understand the behavior of numerical approximations to unbounded solutions, a new definition of g-rel stability was proposed. These notions of relative stability were applied to linear first-order systems with unbounded solutions and it was shown that: (1) For negative eigenvalues, there is no critical time step which marks a transition to unbounded growth (2) The forward Euler method and Runge-Kutta method are relatively stable for systems with negative eigenvalues and the major practical concern is the growth in error. For second-order equations with an exponential dichotomy, the truncation error for the central difference method is of second order, whereas that of the Runge-Kutta is of fourth order. The central difference method also has a second-order error in the coefficients which causes a deterioration in performance for certain initial conditions. However, when the field equations of viscoplasticity are integrated using the central difference method coupled with a semi-implicit update of the constitutive equations, the central difference performs better than the Runge-Kutta method. This superiority is apparent in the later stages of the localization problems when the strain rate grows rapidly as the stress begins to drop in the slab. This behavior cannot be explained by a truncation error analysis of the linearized system and indicates strongly non-linear behavior in the late stages of localization. Numerical studies of a one-dimensional slab problem with softening and localization show that the solutions are Neumann unstable but satisfy the conditions of generalized relative stability [16]. However, these studies have shown that the error grows quite large due to round-off errors. Therefore the time step must be controlled so that the strain increment in each step is below a specified value.

Acknowledgment

The support of the U.S. Air Force Office of Scientific Research and the National Science Foundation is gratefully acknowledged.

358

M. Kulkarni et al.

I Comput. Methods Appl. Mech. Engrg. 124 (1995) 335-363

Appendix A: Absolute stability of the central difference method for dynamic viscoplasticity The governing equations (6.1)-(6.3) are replaced by central difference approximations with a uniform spatial grid interval AX and time interval At. A superscript represents the time step and a subscript represents a grid point. The difference forms of the governing equations are given below. These difference equations are then linearized and the amplification matrix, which relates the field variables at time step y1+ 1 to the values at time step n is derived A. 1. Difference

equations

Eqs. (6.1)-(6.6) momentum

yield the following difference equations

equation

n

?I+112

n

uj+l/2

&

‘Tj-li2

-

aj-1/2

”

@j+lI2

=P n

=

’j

-

v2

(Ala)

At

pr(vl+“’

- f,~y-l’~)

@lb)

where Ax (AIc)

‘=ht viscoplastic constitutive n

n-1/2 ‘j+l

n-l

uj+1/2

-

(+j+112

At n-1 @j+ll2

n

uj+li2

equation

-

=E

Y”-“2 I

Ax

- 4$z2

=~(U;,‘12-y;-1’2)-EAt~:I,11~~

(A24 @2b)

To relate the value of the strain rate at time step (n + 1) to the value at time step n, the strain rate is expanded in a Taylor series about the previous step and terms of first order are retained. +;;‘:i:’ = +rL,‘:,’ + 4.c ‘5+1/z + 4.y AYi+1/2

(A3a)

where +,V and $.., are the derivatives of $(o; y) with respect to u and y, respectively. In the above expression, A5 + 1,2 = a,:;;;;

- a;,-,::2

Wb)

= a;f1,2 - a,=-,,

Similarly, the viscoplastic strain increment

(A3c) may be written as

AY/+l/z = ~j’Zl’;z’- ~yLl’/‘2’

(A3d)

A generalized trapezoidal rule is used to obtain an estimate for the effective viscoplastic strain AY,, l ,2 =

‘Yj+l/*

At[@;;‘:/,

+ (1 - +#I;;;;;]

(A3e)

where 0 represents an integration parameter; 0 G 0 < 1. This represents a family of constitutive updates, and 8 = 0 represents the forward Euler constitutive update. Let F, =

+,,/(l -

5 = &,/(l Substituting

‘WJ,~At)

(A4a)

-WJ,~ At)

(A4b)

Eq. (ARC), Eq. (A3e) in Eq. (A3a) gives

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

= F,(u;+,

+;:;i;

,2

-

a,=-;,,)+ (1 + F,

359

124 (1995) 335-363

(A3

At,+:,-,‘;;

Combining the above with (A2a) yields =F,

4;:::;

1

~(U~~:‘2--Ur-1’2)-EAt~;;:,:l’

I

+ (1 + F, A+&;;

(Aha)

n+1/2 FCTE 4 I+112 =--($;;‘2-~;-“2)+(l+FyAt-EF,At)~;,;;;

(A6b)

From Eq. (A2a) n gj+l12

” oj~l/2

n ‘+j+ll2

=rT:1;11,2+~(~~~:i2-U;1i2)-EAt~:’;11:L2

(A7a)

=~:‘,‘,2+~(ul-1’2-~~~~12)--EAt~;’_;’::2

(A7b)

-

n ‘+j-l/2

=

n-1 u~+112

-

n-1 aj-li2

+

E

r

(v,=_:”

-

2~;~”

+

v;:,“~)

-

E At(r&,‘;

- #I;:~‘,?)

(A7c)

Substituting Eq. (A7c) in Eq. (Ala) yields an expression for the updated value uY+~‘~of the velocity in terms of the historical variables. $+I/2

_

ql12

=~(~,:2-“~~~,2)+~(,;;:12-2u”-l/2+u;_11;2) pr - $$

I

<+;+::2’ - c#q;;;)

(A@

Expressions for the updated value of the stress (T;+,,~ and the updated viscoplastic strain rate $TTI1:,’ in terms of the historical variables are as follows n

n-l -

uj+l/2

A.2.

=~(U:‘,“2-U:1~1’2)-EAt9:+~~~

uj+ll2

Amplification

(A9)

matrix

Eqs. (A8)-(AlO) relate the velocities at step II to the velocities at step n - 1; the values of the stresses and strain rates at step (n + l/2) to the values at step (n - l/2). To derive the amplification matrix: (1) all the quantities are shifted forward by one step; (2) quantities evaluated at the midpoint of the time step, then, are shifted back by half a step. The relevant update equations at a generic time step n + 1 are summarized below @;;,2 = q ?I+1 uj+1!2

n+l

ui

(u;+, - u;) + (1 + F, At - EF, At)+,=,,,

=u;,,_,+;(u;+,

(All)

-u;)-EAt+;+,,2

1 E = u; + pr (u;+,,2 - UjnPli2)+ 2 (UT+,- 2u; + r&J

Following a Neumann procedure, form

G412) - y

(4;+,l, -

4jn_,/2)

(A13)

pr

assume that the difference equations admit solutions of the following

360

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

(T,? = &An

(AW

v,: = &”

(AW

Here, p = elkI Ax, where i = fl and k is the wave number. The update equations (All)-(A13) are recast in terms of the variables 5, A, v to obtain the amplification matrix which relates these variables at the it + lth step to their values at the nth step. Using the substitutions (A14)-(A16), the updated independent variables may now be written as follows 5 “+‘=(1+F,At-EF,At)~“+~(~1’2---i’2)~~ V

A

l

II+1

n+l

vfl + -.@ (p

zz

-E

Recognizing l/2 F

(A17)

5” +;

At

(p112 - p-“‘)vn

+ A”

- /p2)A”

W-9 (A19)

that -112 =

--CL

ia

(A20)

~+cL-l-2=-(y2

(A21)

where a =2sin-

khw 2

(A221

Eqs. (A17)-(A19)

may be written in the following form

d n+l = Gd”

(~23)

dT = (5, A, v>

(A24)

where

and G is the amplification matrix defined as follows

G425)

We seek stability conditions for the viscoplastic equations (6.1)-(6.3). Mathematically, we seek conditions under which all the eigenvalues of the amplification matrix lie within the unit circle. Since the closed form solution is growing, we expect, a priori, that in the softening regime at least one eigenvalue will lie outside the unit circle. The eigenvalues ei of the amplification matrix, which is derived in Appendix A.2 satisfy the following equation

(A26)

The characteristic

polynomial

is then a cubic in the ei and is given as follows

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

361

e; - ((3 - R) + At(F, - EF,))ef - ((R - 3) + At{2EF, + (R - 2)})e, - (( 1 + At(F, - EF,)) = 0 (A27) where R = (w/r)*,

c = (E/P)“~,

the elastic wave speed.

REMARK.

Recall from Eqs. (Ale) and (A22) that r = Ax/At, and CY= 2 sin@ AX/~). The maximum value of (Y is +2 and consequently at this value of (Y, R = 4 represents the Courant Frederich Levy condition for stability of the difference equations of the one-dimensional wave equation.

A .3. Z-transform

To find the appropriate to the left-hand plane.

stability condition, the following transformation

1+z e-I =1-z The characteristic

is used to map the unit circle

6428) equation may be written in terms of z as

c,z3 + c,z* + c,z2 + c, = 0

(~29)

C,) = (8 - 2R) + At((4 - R)F, - 4EF,}

(A30)

C, = At{(R - 4)Fy + 4EF,}

(A31)

C, = 2R i- RF, At

~432)

C, = -RF_, At

(A33)

where

C,C, - C&I, = 8ERF,

At

6434)

The necessary and sufficient conditions that the roots of the characteristic equation lie in the left-hand z-plane, which are necessary for absolute stability is given by the Hurwitz matrix criterion, which may be stated as

A.4.

C,, > 0

(A35)

c, >o

(A36)

c,>o

(A37)

c,>o

(A38)

CC*-c,c,>o

(A39)

Practical stability criteria for the hardening regime

Eqs. (6.4a)-(6.4e) indicate that for power law rate dependent models, +,,, > 0, for softening, +..y < 0 for hardening and +,, > 0 irrespective of the sign of the hardening modulus. From Eqs. (A4a)-(A4b), it follows that in the hardening regime, F, > 0 and F, < 0. Thus, stability conditions (A33) and (A34) are satisfied in he hardening regime. As a first step in evaluating the stable time step, we show that the stable time step must be less than the Courant number. It is shown that R = 4 represents the Courant condition and C,=(8-2R)-C,.

WO)

Furthermore, for stability C,, > 0 and C, > 0, thus if R > 4, i.e. if the time step is above the Courant number, C,, < 0. Thus a necessary condition for stability is R < 4. We next consider the stability condition (A37), i.e. C, >O. This leads to

M. Kulkami et al.

362

I Comput. Methods Appl. Mech. Engrg. 124 (1995) 335-363

6441) The stability condition (A35) may be rearranged

in the following manner

C, = 8 + 4 At@, - EF,) - C,

(~42)

Since C, 2 0, C,, 2 0 for stability, a necessary condition which follows from Eq. (A42) is that 8 + 4At(F,, - EF,) z=0 Since F,, < 0 and F, > 0 in hardening,

(A43a) the above condition may be written as (A43b)

Substituting for F, and F, from Eq. (A4a) and Eq. (A4b), this condition may be written in terms of the field variables as

Ars if,l(l-ie)+%

(A43c) Condition (A43c) is more restrictive than Eq. (A41) and thus provides a necessary condition for stability in the hardening regime. Sufficient conditions for stability can only be obtained by solving Eq. (A30). AS.

Softening regime

For this case, necessary conditions for the eigenvalues to lie within the unit circle are that F, > 0, and F, < 0. As shown in Eqs. (6.4a)-(6,4d), +,., >O, for softening, +,., O irrespective of the sign of the hardening modulus. In the softening regime the conditions F, > 0, and F, < 0 lead to two mutually exclusive conditions; 1 - 0$,, At > 0 and 1 - #+,, At < 0. These conditions cannot be satisfied for a positive At. This leads to the expected result that in the softening regime, at least one eigenvalue lies outside the unit circle. This implies that the central difference method cannot be absolutely stable for any value of the time step.

References [II A. Marchand

and J. Duffy, An experimental study of the formation process of adiabatic shear bands in a structural steel, J. Mech. Phys. Solids 54 (1988) 806-812. and time integration procedures, in: T. Belytschko and T.J.R. Hughes, 121T. Belytschko, An overview of semidiscretization eds., Computational Methods for Transient Analysis (North Holland, Amsterdam, 1987) l-65. A tangent modulus method for rate dependent solids, Comput. Struct. 18(5) (1984) 131 D. Peirce, C.F. Shih and A. Needleman, 875-887. B. Moran and M. Kulkarni, On the crucial role of imperfections in quasi-static solutions, J. App. Mech., [41 T. Belytschko, ASME 58 (1991) 658-665. and quadrilateral with orthogonal hourglass control, JJNME PI D.P. Flanagan and T. Belytschko, A uniform strain hexahedron 17 (1981) 679-706. of Differential Equations (Academic Press, San Diego, 1989). 161 D. Zwillinger, Handbook 171 A. Bayliss, A double shooting scheme for certain unstable and singular boundary value problems, Math. Comput. 32(141) (1978) 61-71. to the Numerical Solution of Differential Equations (Research Studies Press Ltd., Letchworth, PI D. Quinney, An Introduction Hertfordshire, England, 1985). Finite Difference Equations and Simulations (Prentice Hall, Englewood Cliffs, NJ, 1968). [91 F-B. Hildebrand, Methods in Ordinary Differential Equations (Wiley, London and New York, 1973). [lfll J.D. Lambert, Computational Mathematics (Wiley Eastern Limited, New Delhi, 1983). 1111 E. Kreyszig, Advanced Engineering Material rate dependence and mesh sensitivity in localization problems, CMAME 67 (1988) 69-85. [121 A. Needleman, 1131 M.C. Pease, Methods of Matrix Algebra (Academic Press, New York and London, 1965).

M. Kulkarni et al. I Comput. Methods Appl. Mech. Engrg.

124 (1995) 335-363

[14] J. Ortega, Matrix Theory; A Second Course (Plenum Press, New York and London, 1987). (151 L.E. Dickson, First Course in the Theory of Equations (John Wiley, London, 1922). [16] H.L. Schreyer and Z. Chen, One-dimensional softening with localization, J. App. Mech., ASME 53 (1986) 791-797.

363