Geometric integrators for multiplicative viscoplasticity: Analysis of error accumulation

Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Geometric integrators for multiplicative viscoplasticity: Analysis of error accumulation A.V. Shutov *, R. Kreißig Institute of Mechanics and Thermodynamics, Chemnitz University of Technology, Str. d. Nationen 62, D-09111 Chemnitz, Germany

a r t i c l e

i n f o

Article history: Received 8 July 2009 Received in revised form 13 October 2009 Accepted 2 November 2009 Available online 10 November 2009 MSC: 74C20 65L20 Keywords: Viscoplasticity Finite strains Contractivity Exponential stability Inelastic incompressibility Integration algorithm Error accumulation

a b s t r a c t The inelastic incompressibility is a typical feature of metal plasticity/viscoplasticity. Over the last decade, there has been a great amount of research related to construction of numerical integration algorithms which exactly preserve this property. In this paper we examine, both numerically and mathematically, the excellent accuracy and convergence characteristics of such integrators. In order to simplify the considerations, we consider strain-driven processes without hardening effects. In terms of a classical model of multiplicative viscoplasticity, we illustrate the notion of exponential stability of the exact solution. We show that this property enables the construction of effective and stable numerical algorithms, if the incompressibility is exactly satisfied. On the other hand, if the incompressibility constraint is violated, spurious degrees of freedom are introduced. In general, this results in the loss of the exponential stability and a dramatic deterioration of convergence of numerical methods. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction The mechanical processing of materials may involve very large inelastic deformations. For instance, for equal channel angular extrusion of aluminium alloys, the introduced accumulated inelastic strain usually varies between 100 and 900% (depending on the number of extrusions [34]). Even larger deformations can be introduced by some incremental forming procedures like spin extrusion [20] (the accumulated inelastic strain ranges up to 1000%). Due to the highly nonlinear character of the underlying mechanical problem, a correct numerical simulation of such ‘‘long” processes is by no means a trivial task. It is desirable to have numerical algorithms which would be stable with respect to numerical errors, even if working with big time intervals and big time-steps. The assumption of exact inelastic incompressibility is widely implemented for construction of material models of metal plasticity and creep (see, for instance, [11]). Extensive studies were carried out concerning the construction of so-called geometrical

* Corresponding author. Tel.: +49 0 371 531 35024; fax: +49 0 371 531 23419. E-mail address: [email protected] (A.V. Shutov). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.11.002

numerical integrations which exactly preserve the incompressibility of the inelastic flow [4,10,12,19,22,26,27,30,31].1 In this paper, we assess those factors that result in more accurate computations, especially when integrating with big time-steps and for long times. To this end, we analyze the structural properties of the inelastic flow governed by a classical material model of multiplicative viscoplasticity. The material model is based on the multiplicative decomposition of the deformation gradient into inelastic and elastic parts. For simplicity, no hardening behavior is considered in this paper. However, the proposed methodology can be generalized to cover more complicated material behavior as well.2 Only strain-driven processes are considered in this paper.

1 The incompressibility condition is given by a linear invariant in the case of infinitesimal strains inelasticity. Since the linear invariants are exactly conserved by most of integration procedures (cf. [6]), the problem of the conservation of incompressibility only appears when working with finite strains. 2 By means of a series of numerical tests, it was shown for strain-controlled processes in [27] that the use of geometric integrators allows to eliminate the error accumulation even in the case of a complex material behaviour with nonlinear isotropic and kinematic hardening. In general, however, the construction of consistent integration procedures for the finite strain inelasticity is still an open problem (cf. [33]).

701

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Nomenclature inelastic right Cauchy–Green tensor (see (24)) Ci e T 2nd Piola–Kirchhoff tensor (see (26)) 1 second-rank identity tensor A  B ¼ AB product (composition) of two second-rank tensors A:B scalar product of two second-rank tensors AB tensor product of two second-rank tensors kAk l2 norm of a second-rank tensor (Frobenius norm) induced norm of a second-rank tensor (spectral norm) kAk (see (1)) deviatoric part of a tensor ðÞD transposition of a tensor ðÞT inverse of transposed tensor ðÞT trðÞ trace of a second-rank tensor ðÞ unimodular part of a tensor (see (2))

The analysis of a general initial-boundary value problem lies beyond the scope of the current work. We pay especial attention to the exponential stability of the inelastic flow, which is the key notion of the current study. We say that the solution to a Cauchy problem is exponentially stable, if for small perturbations of initial data an exponential decay estimate holds (see Section 2.1). From mechanical standpoint, the exponential stability implies fading memory behavior.3 Moreover, the exponential stability is deeply connected to contractivity (B-stability) of the system of equations, which can be used for stability analysis of numerical algorithms (see the monograph by Simo and Hughes [29]). The main conclusions of this paper regarding the problem of finite strain viscoplasticity are as follows.  The exact solution is exponentially stable for strain-driven processes with respect to small perturbations of initial data, if the incompressibility constraint is not violated.  In the case of exponential stability, the numerical error is uniformly bounded. In particular, there is no error accumulation even within large time periods.  If the incompressibility constraint is violated by some numerical algorithm, then, in general, the numerical error tends to accumulate over time. There exists a rich mathematical literature dealing with existence, uniqueness, regularity, and asymptotic behavior of solutions for certain plasticity/viscoplasticity problems in the context of infinitesimal strains (see [1,7,13] and references therein). A class of material models of monotone type which includes the class of generalized standard materials was defined and analyzed by Alber in [1]. In the context of finite strain viscoplasticity, however, only few theoretical works exist. Some preliminary investigations have been made by Neff in [23]. In this paper, we analyze the well-known material model of multiplicative viscoplasticity. The stability is proved analogously to the classical Lyapunov approach, based on the use of Lyapunov-candidate-functions. More precisely, the hyperelastic potential is used to construct a suitable Lyapunov-candidate (cf. [29]). The paper is organized as follows. In Section 2, we define the notion of exponential stability and prove the main theorem, which states that the numerical error is uniformly bounded if the exact solution is exponentially stable. A simple one-dimensional exam-

3

As Truesdell and Noll [32] have formulated, ‘‘Deformations that occurred in the distant past should have less influence in determining the present stress than those that occurred in the recent past”.

symðÞ hxi wel Distð; Þ K ki f F Sym M

qR k

l

symmetric part of a tensor MacCauley bracket (see (14)3) specific free energy density ‘‘distance” between two solutions (see (35)) yield stress proportionality factor (inelastic multiplier) (see (27)1) overstress (see (27)2) norm of the driving force (see (27)3) space of symmetric second-rank tensors invariant manifold (cf. (4), (28)) mass density in the reference configuration bulk modulus (see (30)) shear modulus (see (30))

ple is presented. In the next section, a classical material model of multiplicative viscoplasticity is formulated in the reference configuration. The change of the reference configuration is likewise discussed. Section 4 contains the definition and analysis of the distance between two solutions in terms of energy (Lyapunov-candidate). Next, the time-evolution of the distance is evaluated and the exponential stability of the exact solution is proved. Finally, the results of numerical tests are presented in Section 5, which illustrate the excellent accuracy and convergence characteristics of geometric integrators. We conclude this introduction with a few words regarding notation. Expression a :¼ b means a is defined to be another name for b. Throughout this article, bold-faced symbols denote first- and second-rank tensors in R3 . A coordinate-free tensor setting is used in this paper (cf. [14,28]). The scalar product of two second-rank tensors is defined by Ap: ffiffiffiffiffiffiffiffiffiffi B ¼ ffi trðABT Þ.4 This scalar product gives rise to the norm by kAk :¼ A : A. Moreover, we denote by k  k the induced norm of a tensor

kAk :¼ max kAxk2 ;

kxk2 :¼

kxk2 ¼1

pffiffiffiffiffiffiffiffiffi x  x:

ð1Þ

The overline ðÞ stands for the unimodular part of a tensor

A :¼ ðdet AÞ1=3 A:

ð2Þ D

1 trðAÞ1. 3

The deviatoric part of a tensor is defined as A :¼ A  The notation O stands for ‘‘Big-O” Landau symbol: f ðxÞ ¼ OðgðxÞÞ as x ! x0 if and only if there exists C < 1 such that kf ðxÞk 6 CkgðxÞk as x ! x0 . The inequality f ðxÞ 6 OðgðxÞÞ is understood as follows: there exists f ðxÞ ¼ OðgðxÞÞ such that f ðxÞ 6 f ðxÞ. 2. Differential equations on manifolds and exponential stability 2.1. General definitions Let us consider the Cauchy problem for a smooth function yðtÞ 2 Rn

_ yðtÞ ¼ f ðyðtÞ; dðtÞÞ;

yðt0 Þ ¼ y0 :

ð3Þ

Here, the initial value y0 and the function dðtÞ are supposed to be ~ðt; y0 ; t 0 Þ. In particular, given.5 Denote the exact solution to (3) by y we have

~ðt 0 ; y0 ; t0 Þ ¼ y0 : y 4 Although the notion of trace of a second-rank tensor is coordinate-free (cf. [28]), we may use a Cartesian coordinate system in order to define it as follows: trðAÞ :¼ A11 þ A22 þ A33 . 5 The system (3) is a system with input, and dðtÞ is interpreted as a forcing function.

702

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Suppose that all solutions lie on some manifold M  Rn

~ðt; y0 ; t 0 Þ 2 M; y

for all t P t 0 ; y0 2 M:

Next, using (5) we obtain for all k ¼ 1; 2; . . . ; n  1

ð4Þ

Then we say that (3) is a differential equation on the manifold M (cf. [5,6]). Next, we say that the solution yðtÞ to the problem (3) is locally exponentially stable on M, if there exist d > 0; c > 0; C 1 < 1, such that the following decay estimate holds

       ~  ð1Þ ð2Þ  cðtt0 Þ  ð1Þ ~ t; yð2Þ y0  y0 ; y t; y0 ; t0  y 0 ; t0  6 C 1 e ð1Þ ð2Þ 0; y0 ; y0

 n k k  y ~ð t; y; tÞ  y ~ðn t; k1 y; k1 tÞ  n k  ~ðk t; k1 y; k1 tÞ: ~ðk t; k y; k tÞ  y 6 C 1 ecð t tÞ y

ð9Þ

Substituting (9) in (8), we get

  ~ðn t; n1 y; n1 tÞ kn y  yðn tÞk 6 n y  y þ C1

ð5Þ

n1 X

ecð

n tk tÞ

 k k k  y ~ðk t; k1 y; k1 tÞ: ~ð t; y; tÞ  y

k¼1

ð10Þ

ð1Þ ky0

for all t0 P 2M such that  yðt 0 Þk 6 d; ð2Þ ky0  yðt0 Þk 6 d. We note that a somewhat different interpretation of the exponential stability can be met in the literature as well (cf., for example, Section 2.5 of [15]). Next, let us consider a numerical algorithm which solves (3) on the time interval ½0; T. Denote by n y the numerical solutions at time instances n t, where 0 ¼ 0 t < 1 t < 2 t <    < N t ¼ T, and 0 y ¼ y0 . Suppose that the error on the step is bounded by the second power of the step size. More precisely

 nþ1 n n  y ~ð t; y; tÞ  nþ1 y 6 C 2 ðnþ1 t  n tÞ2 ;

ð6Þ

where C 2 < 1 (cf. Fig. 1). For simplicity, we will consider constant time-steps only: Dt ¼ nþ1 t  n t ¼ const.

~ðk t; k y; k tÞ ¼ k y. Without loss of generality, we can asObviously, y sume that C 1 P 1. Next, substituting error estimation (6) into (10), we obtain 2

kn y  yðn tÞk 6 C 1 C 2 ðDtÞ

n X

ecð

n tk tÞ

ð11Þ

:

k¼1

But, n t  k t ¼ ðn  kÞDt. Thus, taking into account the well-known P1 i  expression for an infinite geometric series i¼0 r ¼ 1=ð1  rÞ for 1 ¼ 1 þ OðDtÞ, we get for small Dt jrj < 1 and that 1þOð DtÞ n X

n tk tÞ

ecð

6

1 X

eicDt ¼

i¼0

k¼1

¼ 2.2. Main theorem

1

1 1 ¼ 1  ecDt cDt þ OððDtÞ2 Þ

1

cDt 1 þ OðDtÞ

¼

1

cDt

ð1 þ OðDtÞÞ 6 2

1

cDt

:

ð12Þ

Finally, it follows from (11) and (12) that With definitions from the previous subsection we formulate the following theorem.

kn y  yðn tÞk 6

~ðt; y0 ; 0Þ be the exact solution. Suppose that Theorem 1. Let yðtÞ ¼ y conditions (5) and (6) hold. Moreover, suppose that the numerical solution of problem (3) lies exactly on M. Then there exists a constant C < 1 such that

Remark 1. The proof is essentially based on the assumption that k y 2 M. In general, if the numerical solution k y leaves the manifold M, the decay estimation (5) is not valid.

kn y  yðn tÞk 6 C Dt;

as Dt ! 0:

ð7Þ

Here, the constant C does not depend on the size of the time interval ½0; T. Proof. The proof is a modification of the standard error analysis (cf. [2]). In this paper we prove the theorem under the assumption that d ¼ 1. The proof can be easily generalized to cover arbitrary values of d > 0 by using mathematical induction and by assuming ~ðn t; 0 y; 0 tÞ ¼ yðn tÞ. The Dt 6 cd=ð2 maxðC 1 ; 1Þ C 2 Þ. First, note that y triangle inequality yields (cf. Fig. 1)

  ~ðn t; n1 y; n1 tÞ kn y  yðn tÞk 6 n y  y  n n1 n1  ~ð t; y; tÞ  y ~ðn t; n2 y; n2 tÞ þ    þ y  n 1 1  ~ð t; y; tÞ  y ~ðn t; 0 y; 0 tÞ: þ y

2C 1 C 2

c

Dt;



Remark 2. The theorem states that the error is uniformly bounded in the case of exponential stability. Thus, there is no error accumulation in the sense that the constant C in (7) does not depend on the overall time T. Moreover, let  > 0 be some small value. By choosing Dt 6 c=ð2C 1 C 2 Þ the numerical error kn y  yðn tÞk is guaranteed to be less than . Remark 3. If the exponential stability is replaced by the assumption that the right-hand side of (3) is a smooth function of y, a weaker error estimation is valid (cf. [2])

kn y  yðn tÞk 6 CeLT T Dt; ð8Þ

as Dt ! 0:

as Dt ! 0;

ð13Þ

where L ¼ sup kfy k. The effect of the growing multiplier on the righthand side of (13) is referred to as an effect of error accumulation. In that case, in order to guarantee a sufficient accuracy, the upper bound for Dt must depend on T. That makes the practical solution of some problems extremely expensive for large values of T. 2.3. One-dimensional example Let us consider a simple example which illustrates the notion of exponential stability. We examine the response of a one-dimensional viscoplastic device shown in Fig. 2a. The closed system of (constitutive) equations is as follows: The total strain is decomposed into the elastic part ee and the inelastic part ei

e ¼ ee þ ei : Fig. 1. Analysis of error accumulation.

The stress law

r on the elastic spring is governed by the elasticity

703

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

(a)

(b)

Fig. 2. Rheological model (a), and inelastic flow under monotonic loading (b).

r ¼ Eee ; E ¼ const > 0:

sponding inelastic flow of the Maxwell material is exponentially stable as well.

The time derivative of the inelastic strain is given by

1 r e_ i ¼ h f i ; f :¼ jrj  K; hxi :¼ maxðx; 0Þ; g jrj

ð14Þ

where the material constants K > 0 and g > 0 are referred to as yield stress and viscosity, respectively. In order to use the results of previous subsections, we rewrite the problem in the form

1

e_ i ¼ e_ i ðei ; eðtÞÞ ¼ hEjeðtÞ  ei j  K i signðeðtÞ  ei Þ: g ð1Þ

ð15Þ

ð2Þ

Let ei ðtÞ and ei ðtÞ be two solutions to (15). Following [29], we reqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ ð2Þ call that 12 Eðei  ei Þ2 defines an energy norm which is the natu-

3. Material model of multiplicative viscoplasticity In this section we discuss a classical material model of finite strain viscoplasticity. The main ingredients of this model are: multiplicative split of the deformation gradient, hyperelastic relations of Neo-Hooke type, classical Mises–Huber yield function in combination with the associative flow rule (J 2 flow theory), Perzyna-type viscoplasticity. The model is a special case among some relatively new models of viscoplasticity (cf. [10,18,27]). The rheological motivation for the constitutive equations is presented in Fig. 2a.

ral norm for the problem under consideration.6 Next, we consider a monotonic loading

3.1. Constitutive equations

eðtÞ ¼ e_ t; e_ ¼ const > 0:

Let us consider the multiplicative split of the deformation gradient F

Let us show that the exact solution satisfying the initial condition ei ¼ 0 is exponentially stable.7 Without loss of generality, we can asðkÞ sume that t 0 ¼ 0 in estimation (5). If jei ð0Þ  0j 6 d for k 2 f1; 2g, 0 0 then there exists a time instance t ¼ t ðdÞ such that the condition ð1Þ ð2Þ f P f0 > 0 holds for both solutions (ei ðtÞ and ei ðtÞ), if t P t 0 (see Fig. 2b). Then, under that assumption, we obtain

    @ e_ i ðei ; eðtÞÞ E ð1Þ ð2Þ ¼  ; e_ i ei ; eðtÞ  e_ i ei ; eðtÞ @ ei g  E  ð1Þ : ei  eð2Þ ¼ i

g

ð16Þ

g

ð17Þ

0

(17) over ½t ; t and taking the contractivity into account, we get

2 1  2 2E 0 1  ð1Þ ð2Þ ð1Þ ð2Þ E ei ðtÞ  ei ðtÞ 6 E ei ð0Þ  ei ð0Þ e g ðtt Þ : 2 2 Taking the square root of both sides we obtain the required expoEt0 nential decay estimation (5) with C 1 ¼ e g < 1 and c ¼ gE > 0. We note that the equations for the Maxwell material are covered by Eq. (14) as a special case (as K ¼ 0). Therefore, the correIt is known (see [29]) that

1 2E



2

ð2Þ eð1Þ ðtÞ  ei ðtÞ i

  b i :¼ 1 1  FT F1 : C i i 2

ð19Þ

b b T S :¼ b F 1 e SFe ;

e :¼ F1 SFT : T

ð20Þ

The hyperelastic potential relation for stresses and the classical Mises–Huber yield function (see, for example, [29]) are used in this paper8

b S ¼ 2qR

  be @wel C be @C

;

  rffiffiffi 2  b b D K: f :¼ ð C e SÞ   3

ð21Þ

The mass density qR > 0, the isotropic real-valued function wel , and the initial yield stress K are assumed to be known. Next, a covariant Oldroyd derivative is introduced as follows M

ðÞ :¼

d Li; ðÞ þ b L Ti ðÞ þ ðÞ b dt

b L i :¼ F_ i F1 i :

ð22Þ

is not increasing. This effect is

referred The contractivity can be easily proved by checking that   to as contractivity. 2 ð1Þ ð2Þ d 1 E e ðtÞ  e ðtÞ 6 0. dt 2 i i 7

b e :¼ F bT b C e Fe ;

S :¼ ðdet FÞT;

Due to the contractivity (for details see [29]),  2  2 ð1Þ 0 ð2Þ 0 ð1Þ ð2Þ 1 1 E e ðt Þ  e ðt Þ 6 E e ð0Þ  e ð0Þ . Moreover, integrating i i i i 2 2

6

ð18Þ

Here, b F e and Fi stand for elastic and inelastic parts, respectively (see [16,17]). The multiplicative split can be motivated by the idea of a local elastic unloading. A somewhat more consistent motivation can be derived from the concept of material isomorphism [3]. The elastic right Cauchy–Green tensor and inelastic Almansi tensor are defined, respectively by

Let T be the Cauchy stress tensor. The weighted Cauchy tensor (or Kirchhoff stress tensor) as well as the 2nd Piola–Kirchhoff tensors c and K f are defined respectively through operating on K

Therefore, we get from (16)

    d 1 ð1Þ ð2Þ ð1Þ ð2Þ _ ð2Þ e_ ð1Þ Eðei  ei Þ2 ¼ E ei  ei i  ei dt 2 2 E2  ð1Þ ei  eð2Þ ; for t P t 0 : ¼ i

F¼b F e Fi :

For the current example, the geometric property y 2 M is trivial: we put M ¼ R.

8 be b S The yield function f is formulated in terms of the Mandel stress tensor C c A more general operating on K.  definition  qffiffi of the Mises–Huber yield function can be  D 2 found in the literature as f :¼ S   3K, where S stands for the Kirchhoff stress tensor. Using (18) it can be shown that (21)2 follows from this general definition.

704

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Finally, the plastic flow is governed by the associative flow rule in combination with the Perzyna law (cf., for example, [10,11,26,27])9 M b b D b ¼ ki ð C e SÞ  ; C   i beb SÞD  ð C

ki ¼

1



1 f g k0

m ð23Þ

:

Here, g P 0 and m P 1 are known material parameters; k0 > 0 is used to get a dimensionless term in the bracket; the yield function f is interpreted as overstress within the context of Perzyna-type viscoplasticity [24,25]. In order to simplify the numerical treatment of the problem, we rewrite the constitutive equations in terms of tensor-valued varif Along with ables which operate on the reference configuration K. T the well-known right Cauchy-Green tensor C ¼ F F, we introduce an internal variable (inelastic right Cauchy–Green tensor) as

Ci :¼ FTi Fi :

ð24Þ

We rewrite the evolution Eq. (23)1 (cf., for example, [10,26,27]) and supplement it with the initial condition as follows

ki e D Ci ; C_ i ¼ 2 ðC TÞ F

Ci jt¼0 ¼ C0i ;

det C0i ¼ 1;

C0i 2 Sym:

ð25Þ

e the norm of the driving force Here, the 2nd Piola–Kirchhoff tensor T, F, and the inelastic multiplier ki are functions of ðC; Ci Þ, given by

@wel ðCC1 i Þ e ; T ¼ 2qR @C Ci ¼const rffiffiffi

m 1 1 2 K; ki ¼ f ; f ¼F 3 g k0

ð26Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e D 2 : F :¼ tr½ðC TÞ

ð27Þ

Next, we remark that the right-hand side in (25)1 is symmetric (cf. [28]). Moreover, taking into account the property trðABÞ ¼ ¼ detðAÞAT with trðBAÞ and combining the Jacobi formula @ detðAÞ @A the evolution Eq. (25)1, we get d @ detðCi Þ _ ki ki h e D i e D ¼ 0: : Ci ¼ 2 detðCi ÞC1 ðdetCi Þ ¼ i : ðC TÞ Ci ¼ 2 tr ðC TÞ F F dt @Ci Therefore, the exact solution of (25)–(27) has the following geometric property

M :¼ fB 2 Symj det B ¼ 1g:

ð28Þ

We note that the current material model is thermodynamically consistent. That means that the Clausius–Duhem inequality holds for arbitrary mechanical loadings

di :¼

1 e _ d   1  P 0: T :C w CCi 2qR dt el

ð29Þ

One mathematical interpretation of this inequality will be discussed in Section 4.1. To be definite, we use the following expression for the free energy density wel (generalized Neo–Hooke model [11])

9

k  pffiffiffiffiffiffiffiffiffiffiffiffi2 l ln det A þ ðtrA  3Þ; 2 2

ð30Þ

where k > 0; l > 0 are known material constants (bulk modulus and shear modulus, respectively).10 Substituting (30) in (26) we get the 2nd Piola–Kirchhoff stress tensor in the form

e ¼ k ln T

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi D detðCÞ C1 þ l C1 ðCC1 i Þ :

ð31Þ

In this paper we consider strain-driven processes. More precisely, we assume the deformation history CðtÞ to be given. In the following we analyze the exponential stability of the corresponding exact solution Ci ðtÞ.

We note that the constitutive equations of the previous subsection were formulated with respect to a certain fixed reference conf In order to simplify the analysis of the material figuration K. response, it may be useful to rewrite the constitutive equations with respect to some ‘‘new” local reference configuration f (see Fig. 3). We define ‘‘new” variables and corref new ¼ F0 K K sponding constitutive equations with respect to these variables in such a way that the same mechanical response (Cauchy stresses) is predicted. In what follows, we suppose that the mapping F0 is isochoric, i.e. detðF0 Þ ¼ 1. The ‘‘new” deformation gradient, inelastic and elastic parts of the deformation gradient, the right Cauchy tensor, and the inelastic right Cauchy tensor are given by

Fnew :¼ FF1 0 ; new

C

:¼

Moreover, the form of the evolution Eq. (23)1 can be uniquely determined by the principle of maximum plastic dissipation, if the yield function f is given by (21)2 (see [29]).

b F new e T F0 Ci F1 0 ;

Fnew :¼ Fi F1 i 0 ;

1 FT 0 CF0 ;

Cnew i

:¼

:¼ b Fe; 0 1 C0i new :¼ FT 0 Ci F0 :

ð32Þ

The situation is summarized in diagram Fig. 3. e the norm of the driving force The 2nd Piola–Kirchhoff tensor T, F, the inelastic multiplier ki , and the overstress f are transformed as follows

e new :¼ F0 TF e T; T 0

In particular, we get a reduced inequality for relaxation processes ðC ¼ constÞ

d   1  w CCi 6 0: dt el

qR wel ðAÞ :¼

3.2. Change of the reference configuration

Remark 4. The right Cauchy strain tensor C is symmetric. Since the function wel is isotropic, it makes no difference whether the derivative in (26) is understood as a general derivative or as a derivative with respect to a symmetric tensor (cf. [28]).

Ci 2 M;

Fig. 3. Commutative diagram: change of the reference configuration.

Fnew :¼ F;

knew :¼ ki i

f new :¼ f :

Since wel is isotropic, wel ðABÞ ¼ wel ðBAÞ. Using that property, it can be checked that wel ðCC1 i Þ is invariant under the change of the reference configuration

   1 : wel ðCnew Cnew Þ ¼ wel CC1 i i

ð33Þ

Let Ci ðtÞ be the solution of (25)–(27) for a given strain-controlled process fCðtÞgt2½0;T and given C0i . With the help of the relations (32), a one-to-one correspondence between the original and ‘‘new” quantities CðtÞ; Ci ðtÞ, and C0i is defined. It can be easily seen that the 10 As it will be clear in the following, the concrete ansatz for the volumetric part of the strain-energy function does not influence the considerations of the current paper. It was noted in [9] that the volumetric part of the strain-energy function (30) is not convex in det A.

705

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

closed system of equations which describes the evolution of Cnew ðtÞ i is obtained from (25)–(27) by formal replacing of original quantities by their ‘‘new” counterparts. We conclude that the mechanical response becomes invariant under the change of the reference configuration if the loading program CðtÞ and the initial value C0i are transformed according to (32).

Moreover, taking into account that kABk 6 kAk kBk , and that the norms k  k and k  k are equivalent, we get

 kABCk 6 CkAk kBk kCk;

ð42Þ

 < 1. Thus, with some constant C     1=2   ð2Þ 1=2 ð1Þ  ð2Þ 1=2   ð2Þ 1=2  ð1Þ  ð2Þ ð2Þ  C ¼ C  C C 1 C C C i i i i i  i   i    1=2   2  ð1Þ ð2Þ    Cð2Þ   6C  i  Ci Ci ;

ð42Þ

4. Analysis of exponential stability for multiplicative viscoplasticity

ð43Þ

4.1. Measuring the distance between solutions in terms of energy ð1Þ

ð2Þ

Suppose that Ci ðtÞ and Ci ðtÞ are two solutions to the problem (25)–(27) (with the same forcing function CðtÞ). Next, suppose that there exists a constant L < 1 such that

    ðkÞ 1=2   C ðtÞ  < L;  i

    ðkÞ 1=2   C ðtÞ  < L for all t > 0; k 2 f1; 2g:  i ð34Þ

We introduce the following measure of distance between two solutions in terms of energy11

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    ffi ð1Þ ð2Þ ð1Þ ð2Þ :¼ qR wel Ci ðCi Þ1 : Dist Ci ; Ci

ð35Þ

This measure has the following properties: (i) Invariance under the change of the reference configuration (cf. Section 3.2)

    ð1Þ ð2Þ ð1Þ ð2Þ Dist ðCi Þnew ; ðCi Þnew ¼ Dist Ci ; Ci : ð1Þ

ð36Þ

ð2Þ

(ii) For small Ci  Ci , there exist constants C 3 > 0 and C 4 < 1 such that

       ð1Þ  ð1Þ ð2Þ  ð1Þ ð2Þ ð2Þ  C 3 Ci  Ci  6 Dist Ci ; Ci 6 C 4 Ci  Ci :

ð37Þ

  ð1Þ ð2Þ ð1Þ ð2Þ P 0 and (iii) For all Ci ðtÞ; Ci ðtÞ 2 M we have Dist Ci ; Ci

  ð1Þ ð2Þ Dist Ci ; Ci ¼ 0;

ð1Þ

if and only if Ci

ð2Þ

¼ Ci :

ð38Þ

Proof (i) Identity (36) can be proved similarly to the invariance property (33). (ii) First, it follows from (30) that for small D we have (see Appendix A)

k l  2  qR wel ð1 þ DÞ ¼ ðtrDÞ2 þ tr DD þ OðD3 Þ; 8 4

  2  where D ¼ kDk. Note tr DD ¼ kDD k2 for D 2 Sym.  that   2 l D 2 k Thus, 8 ðtrDÞ þ 4 tr D is a norm on Sym. Since all norms on Sym are equivalent, there exist constants C 03 > 0 and C 04 < 1 such that for small D 2 Sym we have

C 03 kDk2 6 qR wel ð1 þ DÞ 6 C 04 kDk2 :

11



ð2Þ

Furthermore, substituting D ¼ Ci

1=2

ð44Þ  1=2 ð1Þ ð2Þ Ci Ci  1 in

(39), and combining it with (41), (43), and (44) we get

  1=2    2   ð1Þ ð2Þ  ð1Þ ð2Þ e 3  Cð2Þ   C  Ci  Ci  6 Dist Ci ; Ci  i   1=2  2   ð1Þ ð2Þ  e 4  Cð2Þ   6C  Ci  Ci :  i

ð45Þ

Finally, combining (45) with (34) we obtain (37). (iii) We note that wel ðAÞ P 0. Moreover, wel ðAÞ ¼ 0 if and only if A¼1 h In view of properties (i)–(iii), the function Dist is a natural measure of distance for the problem under consideration.12 Moreover, the dissipation inequality (29), which holds for all relaxation processes, can be interpreted as follows: During relaxation, the distance (measured in terms of energy) between any soluð2Þ ð1Þ tion Ci and a constant solution Ci C ¼ const is not increasing. 4.2. Sufficient condition for exponential stability Let us consider a loading program (strain-driven process) ð1Þ ð2Þ fCðtÞgt2½0;T on the time interval ½0; T. Let Ci ; Ci 2 M be two solutions. In order to prove the exponential stability, it is sufficient to prove that there exists t0 P 0 and c > 0 such that for all t P t 0 (cf. (17))

  2  2 d ð1Þ ð2Þ ð1Þ ð2Þ : Dist Ci ; Ci 6 c Dist Ci ; Ci dt

ð46Þ

Indeed, in that case, using the Gronwall’s inequality we get from (46) the following decay estimation

 2  2 0 ð1Þ ð2Þ ð1Þ ð2Þ Dist Ci ðtÞ; Ci ðtÞ 6 Dist Ci ðt 0 Þ; Ci ðt 0 Þ ecðtt Þ :

ð47Þ

Combining this result with (37), we get the required estimation of the type (5). Thus, the uniform error estimation of Theorem 1 follows immediately from (47). 4.3. Reduction of the stability analysis to a simplified problem with C¼1

ð40Þ

we have

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    1=2  1=2 ð1Þ ð2Þ ð2Þ ð1Þ ð2Þ Dist Ci ; Ci Ci Ci ¼ qR wel Ci :

   1=2  2  ð2Þ 1=2 ð1Þ  ð2Þ 1=2     Cð2Þ   C : 6C C C 1 i i   i   i

ð42Þ

ð39Þ

Next, due to the property

wel ðABÞ ¼ wel ðBAÞ;

     1=2    ð2Þ 1=2  ð2Þ 1=2 ð1Þ  ð2Þ 1=2  ð1Þ ð2Þ  ð2Þ  Ci Ci Ci 1 Ci Ci Ci  ¼    Ci

ð41Þ

The relation (35) can be seen as a generalization of the energy norm rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ð2Þ 2 1 E eð1Þ  ei (cf. Section 2.3). 2 i

Let t0 be an arbitrary time instance. In this section we discuss a procedure which helps to simplify the examination of the inequality (46) at the time t0 . 12 The function Dist is not symmetric: DistðA; BÞ – DistðB; AÞ. Symmetrized functions can p bffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e d e fi n e d b ffiy DistSym 1 ðA; BÞ :¼ 1=2ðDistðA; BÞ þ DistðB; AÞÞ; Sym Dist2 ðA; BÞ :¼ DistðA; BÞDistðB; AÞ. Nevertheless, none of these functions determine a metric on M, since the triangle inequality does not hold.

706

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

The first simplification of the problem is as follows. We  2 ð1Þ ð2Þ and note that the quantities Dist Ci ðt0 Þ; Ci ðt 0 Þ ,  ð1Þ 0 ð2Þ 0 2 ð1Þ 0 0 d DistðC ðt Þ; C ðt ÞÞ Þ depend solely on Cðt Þ; C ðt Þ, and i i i dt ð2Þ _ 0 Þ. Therefore, on the examination of (46) at Ci ðt 0 Þ but not on Cðt 0 t ¼ t we can replace the actual loading program fCgt2½0;T with a

constant loading (relaxation process): we take a constant Cðt 0 Þ instead of loading CðtÞ, where ðÞ stands for a unimodular part of a tensor. The second simplification is as follows. Let F0 be some ‘‘new” reference configuration and detðF0 Þ ¼ 1. There is a one-to-one corð1Þ ð2Þ respondence between the solutions Ci ðtÞ; Ci ðtÞ of the problem  new ð1Þ ðtÞ; with the forcing function CðtÞ and the solutions Ci  new ð2Þ ðtÞ with the forcing function Cnew ðtÞ (cf. Section 3.2) Ci 1 Cnew ðtÞ ¼ FT 0 CF0 ;



ðkÞ

Ci

new

ðkÞ

1 ðtÞ ¼ FT 0 Ci ðtÞF0 ;

k 2 f1; 2g:

It follows from (36) that

new  new     ð1Þ ð2Þ ð1Þ ð2Þ Dist Ci ðt 0 Þ; Ci ðt 0 Þ ¼ Dist Ci ðt 0 Þ; Ci ðt0 Þ ;

  2  new  new 2  d d ð1Þ ð2Þ ð1Þ ð2Þ ; Dist Ci Dist Ci ðtÞ; Ci ðtÞ ðtÞ; Ci ðtÞ ¼ 0 dt dt t¼t0 t¼t

   ð2Þ ð1Þ  ð1Þ ð2Þ We abbreviate D :¼ Ci  Ci . Since Ci ; Ci 2 M, the following estimations hold (cf. Appendix B)

   ð2Þ1 ð1Þ1 ð1Þ tr Ci  Ci Ci ¼ OðD2 Þ;    ð1Þ ð1Þ1 ð2Þ1 tr Ci Ci ¼ OðD2 Þ;  Ci

D  D   ð2Þ1 ð1Þ ð2Þ1 ð1Þ1 ð1Þ Ci Ci ¼ Ci  Ci Ci þ 1   ð2Þ1 ð1Þ1 ð1Þ Ci þ OðD2 Þ: ¼ Ci  Ci   ð1Þ ð2Þ1 T   T D @wel Ci Ci l ð2Þ1 ð1Þ ð2Þ1 qR  ð1Þ ð2Þ1  ¼ Cð1Þ C C C i i 2 i i @ Ci Ci l ð1Þ1 ð2Þ  ð2Þ1 ð1Þ1  ð1Þ ¼ Ci Ci C i  Ci Ci þ OðD2 Þ 2 l  ð2Þ1 ð1Þ1  ð1Þ ¼  Ci Ci þ OðD2 Þ: C 2 i 

Without loss of generality we assume detðCðt 0 ÞÞ ¼ 1. By choosing F0 ¼ ðCðt 0 ÞÞ1=2 the problem can be reduced to the simplified problem with Cðt 0 Þ ¼ 1.13   2 ð1Þ ð2Þ 4.4. Evaluation of dtd Dist Ci ; Ci   2 ð1Þ ð2Þ at some fixed In this section we evaluate dtd Dist Ci ; Ci time instance t 0 . Without loss of generality (cf. the previous section) it can be assumed that Cðt 0 Þ ¼ 1. In that reduced case, the evolution Eq. (25) takes the form

  D  D       þ OðD2 Þ að1Þ  að2Þ ¼ a  Cið1Þ1    Cð2Þ1 i 

*

pffiffiffiffiffiffiffiffi +m

lx  2=3K

glx

k0

:

   ð1Þ1 D  q2ffiffi  K is poIt can be assumed that the overstress f ¼ l C   i 3 qffiffi q ffiffi 2 2 K. Thus, we suppose K=l < sitive and bounded by 3 3  qffiffi  ð1Þ1 D   6 2 2K=l.  C   i 3 Using the property A : ðBCDÞ ¼ ðBT ADT Þ : C it follows from (49) and (53) that

  2 ð35Þ d    d ð1Þ ð2Þ ð2Þ1 Dist Ci ; Ci ¼ qR wel Cð1Þ i Ci dt dt   ð1Þ ð2Þ1 @wel Ci Ci d  ð1Þ ð2Þ1   : ¼ qR  Ci C i ð1Þ ð2Þ1 dt @ C i Ci l  ð1Þ  ð1Þ1 ð2Þ1  ð49Þ and ð53Þ ¼  Ci  Ci C 2 i   D  D ð1Þ1 ð2Þ1 : að1Þ Ci  að2Þ Ci

ð48Þ

l  ð1Þ  ð1Þ1 ð2Þ1  ð51Þ C þ OðD3 Þ ¼ 2  Ci  Ci 2 i   ð1Þ1 ð2Þ1  að2Þ Ci : að1Þ Ci l  ð1Þ  ð1Þ1 ð2Þ1   Ci þ OðD3 Þ ¼  að1Þ Ci Ci 2   l ð1Þ1 ð2Þ1  ðað1Þ  að2Þ Þ  Ci : Ci 2    ð1Þ ð1Þ1 ð2Þ1 ð2Þ1 : Ci  Ci

Ci Ci

Next, using the product rule we get from (48)1

  d  ð1Þ ð2Þ1  _ ð1Þ ð2Þ1 ð1Þ ð2Þ1 ¼ C i Ci þ Ci Ci Ci C i dt   ð1Þ ð1Þ1 _ ð1Þ ð2Þ1 ð1Þ ð2Þ1 ð2Þ ð2Þ1 ¼ C i Ci þ Ci Ci Ci Ci Ci Ci

 D  D  ð48Þ1 ð1Þ ð1Þ1 ð2Þ1 ð2Þ1 ¼ Ci að1Þ Ci  að2Þ Ci ; Ci where aðkÞ

ð49Þ

   ðkÞ1 D   for k 2 f1; 2g. :¼ a  C   i

ð54Þ

Furthermore, we compute the derivative of wel ðAÞ using a coordinate-free tensor setting (see, for example, [14,28])

 T D @w ðAÞ k pffiffiffiffiffiffiffiffiffiffiffiffi l qR el ¼ ln det A AT þ AT A : @A 2 2

ð50Þ

þ OðD3 Þ ¼ F I þ F II þ OðD3 Þ;

Alternatively, the problem can be reduced to the case  1=2 ð1Þ F0 ¼ Ci ðt0 Þ .

ð1Þ Ci ðt 0 Þ

¼ 1 by choosing

ð55Þ

where F I and F II are given by F I :¼ 

l

F II :¼  13

ð53Þ

  ð1Þ1 D  .  the derivative of aðxÞ at x ¼  Next, denote by a   Ci Therefore,

 a   ð1Þ1 D  ð1Þ1 ð2Þ1 þ OðD2 Þ: ð54Þ : Ci  Ci ¼  ð1Þ1 D  Ci   C   i

 new  new 2  d ð1Þ ð2Þ ðtÞ; Ci ðtÞ Dist Ci 0 dt t¼t new  new  2 ð1Þ ð2Þ 6 cDist Ci ðt 0 Þ; Ci ðt0 Þ :

aðxÞ :¼

ð52Þ

Combining (50) with (52) we obtain

1

ð51Þ

Thus, using (51)1 we get

Therefore, estimation (46) is equivalent to

   1 D   1 D  C C Ci ; C_ i ¼ a  i   i

as D ! 0:

2







ð2Þ1 ð1Þ ð1Þ1 ð2Þ1 Ci C i að1Þ tr Cð1Þ1  Ci  Ci i

l

a

 D   2  Cð1Þ1    i



ð1Þ1 Ci



; D     ð1Þ1 ð2Þ1 ð1Þ1 ð2Þ1 1 : ðCi : Ci  Ci  Ci Þ : ð56Þ

707

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Now, for any pair of real positive numbers ðh; DÞ let us define a subset of M M by

Sðh; DÞ :¼



ð1Þ

ð2Þ

Ci ; Ci



      ð1Þ1 D  ð1Þ ð2Þ   6 h;  2 M Mj Ci  Ci  6 D :   Ci ð57Þ

By definition, put   2að1Þ F II ð2Þ :¼  U Cð1Þ i ; Ci a F I

! D      ðCð1Þ1 Þ  ð1Þ1 ð2Þ1 ð1Þ1 ð2Þ1 i : C 2   C  Ci 1 : Ci i i ðCð1Þ1 ÞD  i     : ¼ ð1Þ1 ð2Þ1 ð1Þ ð1Þ1 ð2Þ1  Ci  Ci Ci Ci tr Ci ð58Þ

There exists a function qðhÞ > 0 such that

    ð1Þ ð2Þ ð1Þ ð2Þ þ OðDÞ for all Ci ; Ci 2 Sðh; DÞ: qðhÞ P U Ci ; Ci

ð59Þ

The numerical evaluation of the function qðhÞ is discussed in the Appendix C. The function can be approximated by qðhÞ  h3 =12 for 0 6 h 6 0:2 (see Appendix C). Moreover, suppose that (cf. Section 4.5) ð1Þ

a

! rffiffiffi 2 : K=l a Pq 2 3

ð60Þ

This condition will be discussed in the next subsection. Multiplying  F I OðDÞ ¼ OðD3 Þ we both sides of (60) by F I a1ð1Þ < 0 and noting that aað1Þ    qffiffi  ð1Þ ð2Þ get for all Ci ; Ci 2 S 2 23K=l; D

! rffiffiffi  a  2 a ð59Þ   ð1Þ ð2Þ  K=l FI 6 q 2 F I 6 U Ci ; Ci þ OðDÞ ð1Þ F I ð1Þ a a 3

Finally, combining (63) with (64) we get the required estimation qffiffi (46) if the assumptions hold: 0 < f0 6 f 6 23K, qffiffi following  að1Þ P q 2 23K=l a . 4.5. Analysis of the sufficient stability condition In this subsection we analyze the condition (60) which was used previously to prove the inequality (46). For brevity, we write ð1Þ subsection. The condition (60) can be Ci instead of C i  in this  qffiffi   1 D   P q 2 2K=l daðxÞ   C rewritten as a  D  (see the   i dx 3 x¼ðC1 i Þ   sketch in Fig. 4).  1 D  q2ffiffi > K=l to ensure that the overNext, we suppose    Ci 3 stress is larger than zero. Using (48)2 it can be easily shown that (60) is equivalent to

   1 D   P xcr ;  C   i where the critical value xcr is given by

rffiffiffi 2 K þ ðm  1ÞlqðhÞ 3 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi1 , !2 u rffiffiffi rffiffiffi u 2 2 C K þ ðm  1ÞlqðhÞ þ 4lqðhÞ K A ð2lÞ; þt 3 3

xcr :¼

qffiffi  qffiffi  with h ¼ 2 23K=l. For small values of q 2 23K=l , a simplified estimation for xcr is valid

0 ! !! 2 1 rffiffiffi rffiffiffi rffiffiffi 2 2 2 A: K=l þ mq 2 K=l þ O@ q 2 K=l xcr ¼ 3 3 3

ð60Þ

ð58Þ

¼ 2F II þ OðD3 Þ:

This result is visualized in Fig. 4 for m ¼ 1.

ð61Þ

Multiplying both sides of (61) by 1/2 and adding 1=2F I þ F II , we get

  2 d ð1Þ ð2Þ Dist Ci ; Ci 6 1=2F I þ OðD3 Þ: dt

ð62Þ

Next, if f P f0 for some f0 > 0, then there exists C 5 > 0 such that

  l  ð1Þ1 ð2Þ1  ð1Þ 1=2 2 ð56Þ 2 F I ¼  að1Þ  C  C C i i  i  6 C 5 D : 2

2

 

ð2Þ1 að1Þ  Cð1Þ1  Ci i

 1=2  2 ð1Þ  Ci 

ð63Þ

 1=2   1=2  2  ð1Þ 1=2  ð1Þ1 2 ð2Þ1 ð1Þ   C  Ci Ci  Ci i    2     2  ð1Þ 1=2 2  ð1Þ 1=2 ð2Þ1  ð1Þ 1=2 l   C ¼  að1Þ C 6  Ci Ci Ci  1 i     2 2      1=2 ð39Þ 1=2  ð1Þ 1=2  l  ð2Þ1 ð1Þ  q wel Cð1Þ 6  að1Þ C 6 =C 04  Ci Ci Ci R i   2       ð1Þ 1=2 2 l ð40Þ  q wel Cð1Þ Cð2Þ1 : ¼  að1Þ C 6 =C 04  ð64Þ C R i i  i  2 6

l

 

ð66Þ

The situation is summarized in Fig. 5. Let us consider a concrete in order to estimate the or qffiffi example 

l . For instance, we put for a typical

2 K= 3

l  25; 000 MPa. Thus, aluminium alloy K  300 MPa, qffiffi ð84Þ 2 23K=l  0:014. Next, qð0:014Þ  ð0:014Þ3 =12  2:3  107 (see

Similarly to the proof of (37) we obtain with some C 6 > 0

l

0 ! !!2 1 rffiffiffi rffiffiffi 2 2 A: K=l þ O@ q 2 K=l fcr ¼ mlq 2 3 3

der of magnitude of q 2

Therefore, for small D, inequality (62) yields

  2 d ð1Þ ð2Þ Dist Ci ; Ci 6 1=4F I : dt

ð65Þ

where the critical overstress fcr is estimated by

Combining this result with (55) we obtain

ð56Þ

Therefore, in terms of the overstress, the condition (60) is equivalent to

f P fcr ;

F I þ F II 6 1=2F I þ OðD3 Þ:

FI ¼ 

   1 D  q2ffiffi  K. We recall that the overstress is given by f ¼ l C   i 3

Appendix C). Therefore, the critical overstress is given by ð66Þ

fcr  m 0:0057 MPa. For physically reasonable values of m ðm 6 100Þ this critical qffiffi value is negligible compared to the size of the elastic domain 23K  245 MPa.

að1Þ C 6  Cð1Þ i

Fig. 4. Sketch of condition (60) for m ¼ 1.

708

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711 nþ1

Fig. 5. Domain of exponential stability.

Remark 5. Since the overstress f is isolated from zero due to the sufficient stability condition (65), the current theory can not be applied to exactly quasistatic processes with f 6 0. On the other hand, the theory is directly applicable to nearly quasistatic processes with the overstress f larger than fcr but still negligible qffiffi compared to the size of the elastic domain 23K.

The numerical implementation of the material model (25)–(27) within a displacement-based Finite Element Method (FEM) with implicit time stepping is based on the implicit integration of the evolution Eq. (25) (see, for example, [29]). This procedure should provide the stresses as a function of the strain history. More precisely, suppose that the right Cauchy–Green tensor nþ1 C at the time tnþ1 ¼ tn þ Dt is known and assume that the internal variable Ci at the time tn is given by n Ci . We need to compute the internal variable Ci at the time t nþ1 in order to evaluate the e ¼ Tð e nþ1 C; nþ1 Ci Þ. stress tensor nþ1 T Note that the norm of the driving force F and the overstress f can be represented as functions of nþ1 C and nþ1 Ci :

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i2 nþ1 nþ1 e ðnþ1 C; nþ1 Ci ÞÞD ; Fð C; Ci Þ ¼ tr ðnþ1 C T

rffiffiffi 2 nþ1 nþ1 nþ1 nþ1 K: f ð C; Ci Þ ¼ Fð C; Ci Þ  3

ð67Þ

For what follows it is useful to introduce the incremental inelastic parameter

n :¼ Dt nþ1 ki : Thus, according to the Perzyna rule, we get the following equation with respect to nþ1 C; nþ1 Ci and n

n¼

m Dt f ðnþ1 C; nþ1 Ci Þ : k0 g

ð68Þ

The remaining equation for finding unknown nþ1 Ci and n is obtained through the time discretization of (25), which will be discussed in the next subsection. 5.1. Euler-Backward Method and geometric implicit integrators We introduce a nonlinear operator Bðnþ1 C; nþ1 Ci ; nÞ as (cf. the right-hand side of (25)1)

Bðnþ1 C; nþ1 Ci ; nÞ :¼ 2

 D n nþ1 e nþ1 C Tð C; nþ1 Ci Þ ; nþ1 nþ1 Fð C; Ci Þ

e nþ1 C; nþ1 Ci Þ and Fðnþ1 C; nþ1 Ci Þ are given by (31) and (67), where Tð respectively. Let us consider the classical Euler-Backward Method (EBM) (see, for example, [4,8,29]) being applied to the evolution problem (25) nþ1

Ci ¼ ½1  Bð

nþ1

C;

nþ1

1 n

Ci ; nÞ

Ci :

ð69Þ

Since the symmetry of the internal variable nþ1 Ci is exactly preserved by the EBM14, this equation is equivalent to 14

Moreover, it was shown in [27] that the symmetry is exactly preserved by the Euler-Backward Method and Exponential Method even in a more general case of a nonlinear kinematic hardening.

ð70Þ

The main disadvantage of the EBM is that the incompressibility condition is violated. In other words, in general, detðnþ1 Ci Þ – 1. In order to avoid this problem, a modification of EBM was proposed by Helm in [12]. The Modified Euler-Backward Method (MEBM) (see [12,27]) uses the following equation which is obtained by applying the projection operator ðÞ to the right-hand side of (69) nþ1

Ci ¼ ½1  Bðnþ1 C; nþ1 Ci ; nÞ

1 n

Ci :

ð71Þ

Next, since the symmetry property is exactly preserved by MEBM ðnþ1 Ci 2 SymÞ, this method can be equivalently rewritten as nþ1

5. Accuracy testing of implicit integrators

Ci ¼ symð½1  Bðnþ1 C; nþ1 Ci ; nÞ1 n Ci Þ:

Ci ¼ symð½1  Bðnþ1 C; nþ1 Ci ; nÞ1 n Ci Þ:

ð72Þ

Finally, the Exponential Method (EM) (see, for instance, [4,21,22,35]) is based on the use of the tensor exponential expðÞ nþ1

Ci ¼ exp½Bðnþ1 C; nþ1 Ci ; nÞn Ci ;

ð73Þ

where

exp½B :¼ 1 þ B þ

1 2 1 3 B þ B þ : 2! 3!

ð74Þ

Since both the symmetry and unimodularity of the numerical solution are exactly preserved ðnþ1 Ci 2 MÞ, the Exponential Method can be rewritten as follows (cf. [27]): nþ1

Ci ¼ symðexp½Bðnþ1 C; nþ1 Ci ; nÞn Ci Þ:

ð75Þ

Combining (68) with one of the discretization methods (Eqs. (70), (72) or (75)) a closed system of nonlinear algebraic equations for finding unknown nþ1 Ci and n is obtained. One possible solution strategy for the resulting system of equations was discussed in [27], and the application of a coordinate-free tensor formalism to the numerical solution was analyzed in [28]. We do not discuss the concrete solution strategies in this paper, since all robust procedures yield the same numerical result. We recall that the geometric property of the exact flow ðCi 2 MÞ is exactly satisfied by MEBM and EM. Therefore, we refer to these two methods as to geometric integrators. For all the three methods, the error on the step is bounded by the second power of the step size (cf. estimation (6)), if the righthand side is a smooth function. On the other hand, strong local nonlinearities due to the distinction into elastic and inelastic material behavior or due to the non-smoothness of the loading function CðtÞ may increase the error on the corresponding time step. 5.2. Testing results The theoretical results obtained in this study are validated via a series of numerical tests. Let us simulate the material behavior under strain-controlled, nonproportional and non-monotonic loading in the time interval t 2 ½0; 300. Suppose that the deformation gradient is defined by

FðtÞ ¼ F0 ðtÞ; where F0 ðtÞ is a piecewise linear function of time t such that F0 ð0Þ ¼ F1 ; F0 ð100Þ ¼ F2 ; F0 ð200Þ ¼ F3 , and F0 ð300Þ ¼ F4 with

Table 1 Material parameters. k (MPa)

l (MPa)

K (MPa)

m ()

g ðs1 Þ

k0 (Mpa)

73,500

28,200

270

3.6

2 106

1

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

709

Fig. 6. Accuracy analysis concerning Euler-Backward Method (EBM), Modified Euler-Backward Method (MEBM), and Exponential Method (EM).

Fig. 7. Summary of the main results.

0

2

B F2 :¼ @ 0 0 0 1 1 1 0 B C F3 :¼ @ 0 1 0 A; 0 0 1

F1 :¼ 1;

0

1

0C A; 0 p1ffiffi2 0 1 pffiffi 2 B F4 :¼ @ 0 2 0 0

p1ffiffi 2

0

1

C 0 A: p1ffiffi 2

Thus, we put

8 > < ð1  t=100ÞF1 þ ðt=100ÞF2 if t 2 ½0; 100; F0 ðtÞ :¼ ð2  t=100ÞF2 þ ðt=100  1ÞF3 if t 2 ð100; 200; > : ð3  t=100ÞF3 þ ðt=100  2ÞF4 if t 2 ð200; 300: The material parameters used in simulations are summarized in Table 1. Next, we suppose that the reference configuration is stress free at t ¼ 0. Therefore, we put

Ci jt¼0 ¼ 1: The numerical solution obtained with an extremely small time step . ðDt ¼ 0:01sÞ will be named the exact solution and denoted by Cexact i Next, the numerical solutions with Dt ¼ 1s and Dt ¼ 0:5s are de  is plotted in Fig. 6. . The error Cnumer  Cexact noted by Cnumer i i i For all three methods the error is proportional to Dt. Moreover, in accordance with Theorem 1 (cf. Section 2.2), the error is uniformly bounded for geometric integrators (MEBM and EM). More precisely, the error is bounded by C Dt where the constant C does not depend on the size of the entire time interval. Next, since the incompressibility condition is violated by EBM, the geometric property (28) is lost and some spurious degrees of freedom are introduced. that  numer In exact  case, only a weaker error estimation is valid: e ðTÞDt, where C e ðTÞ depends on the size T of C  Ci  6 C i the entire time interval. 6. Discussion and conclusion

evolution equations of finite strain plasticity/viscoplasticity, which exactly preserve the inelastic incompressibility condition. The excellent accuracy and convergence properties of such algorithms were tested by numerical computations. Particularly, the long term accuracy of geometric integrators was analyzed in the paper [27], and the absence of error accumulation was numerically verified for strain-driven processes. In the current study, a rigorous mathematical formulation of this phenomenon is proposed. The main result of the current paper is as follows: the numerical error is uniformly bounded by C Dt if the incompressibility condition is satisfied. In terms of a classical model of multiplicative viscoplasticity we prove that all first-order accurate geometric integrators are equivalent in that sense. This theoretical result corresponds with the numerical tests. Indeed, MEBM and EM are equivalent concerning the accuracy and convergence (cf. Fig. 6). The main results are summarized diagrammatically in Fig. 7. The property of the exponential stability of the exact plastic flow was mathematically analyzed in this paper. Obviously, this property must be utilized during the development of new material models and corresponding algorithms in order to improve the accuracy and convergence of numerical computations. Acknowledgements This research was supported by the German Science Foundation (DFG) within the collaborative research center SFB 692 ”Highstrength aluminium based light weight materials for reliable components”.

Appendix A Suppose D ¼ kDk ! 0. Let us show that

k l  2  qR wel ð1 þ DÞ ¼ ðtrDÞ2 þ tr DD þ OðD3 Þ: 8

4

ð76Þ

First, recall the Taylor expansion of detð1 þ DÞ up to second order In the last decade, intensive research has been carried out concerning the development of so-called geometric integrators for the

detð1 þ DÞ ¼ 1 þ trðDÞ þ 1=2ðtrðDÞÞ2  1=2trðD2 Þ þ OðD3 Þ:

ð77Þ

710

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Therefore,

Substituting (82) in (58) and taking (83) into account, we obtain

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi detð1 þ DÞ ¼ 1 þ 1=2trðDÞ þ OðD2 Þ;  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 k k ln detð1 þ DÞ ¼ ðtrDÞ2 þ OðD3 Þ: 2 8 Next, note that for small

ð1 þ eÞ

1=3

ð78Þ

e we have 2

!

ð1Þ1 D



ðCi



ð2Þ ¼ U Cð1Þ i ; Ci

Þ

 2  ðCð1Þ1 ÞD  : X ð1 : XÞ i ð1Þ

trðXCi XÞ

þ OðDÞ:

Thus, we define qðhÞ as 3

¼ 1  1=3e þ 2=9e þ Oðe Þ:

ð1Þ

^ðCi Þ; qðhÞ :¼  max  q ðCð1Þ1 ÞD 6h i

Combining this with (77), we get

l

 l  1=3 trð1 þ DÞ  3 trð1 þ DÞ  3 ¼ ðdet ð1 þ DÞÞ 2 2 l  ¼ 1  1=3trD þ 1=18ðtrDÞ2 2   þ 1=6trðD2 Þ þ OðD3 Þ ð3 þ trDÞ  3  l   2 tr D  1=3ðtrDÞ2 þ OðD3 Þ ¼ 4 l  2  þ OðD3 Þ: ¼ tr DD 4

!

ð1Þ1 D

^ðCið1Þ Þ :¼ q

ð79Þ

Finally, (76) follows from (30), using (78) and (79).

ðCi Þ  2  ðCð1Þ1 ÞD  : X ð1 : XÞ i   : ð1Þ tr XCi X

max

X2T

C

ð1Þ1 M i

  ^ Cð1Þ The function q can be evaluated as follows. First, for each X i  1=2 ð1Þ . Next, define a vector space introduce Y ¼ X Ci    1=2 ð1Þ T :¼ Y 2 Symj Ci : Y ¼ 0 . Thus,

X 2 T Cð1Þ1 M () Y 2 T; i

  ð1Þ tr XCi X ¼ kYk;

Appendix B Let A; B 2 M and kA  Bk ! 0. Let us prove, for instance, that

B1 : ðA  BÞ ¼ OðkA  Bk2 Þ:

ð80Þ

Indeed, since detðÞ is a smooth function, we have

ðCi ÞD  2     ð1Þ1 D  : X ¼ B1 : Y;   C   i ð1Þ1

where

 D  1=2 Cð1Þ1 i ð1Þ   B1 :¼ 2 Ci  ð1Þ1 D  ;  C   i 

@ detðBÞ detðAÞ ¼ detðBÞ þ : ðA  BÞ þ OðkA  Bk2 Þ: @B Next, using the Jacobi formula, we get

detðAÞ ¼ detðBÞ þ detðBÞBT : ðA  BÞ þ OðkA  Bk2 Þ:

1 : X ¼ B2 : Y;

 1=2 ð1Þ B2 :¼ Ci :

Therefore, T

Finally, taking into account that detðAÞ ¼ detðBÞ ¼ 1 and B we obtain (80).

1

¼B ,

  ^ Cð1Þ q ¼ max ½ðB1 : YÞðB2 : YÞ: i Y2T;kYk¼1

Remark 6. Note that for the tangential space T B M to the manifold M in Sym we have

T B M ¼ fX 2 SymjB1 : X ¼ 0g: Thus, relation (80) implies that limA!B ððA  BÞ=kA  BkÞ 2 T B M (if the limit exists).

Next, we compute the orthogonal projections of B1 and B2 on T:

  1=2  1=2 1 ð1Þ ð1Þ B0k :¼ Bk  Bk : Ci Ci  ;   ð1Þ 1=2 2   C   i

k 2 f1; 2g:

Thus, Appendix C

^ðCið1Þ Þ ¼ max q

Y2T;kYk¼1

We need to construct and evaluate a function qðhÞ such that for small D

    ð1Þ ð2Þ ð1Þ ð2Þ þ OðDÞ for all Ci ; Ci 2 Sðh; DÞ: qðhÞ P U Ci ; Ci

  ð1Þ ð2Þ 2 Sðh; DÞ (cf. (57)). It follows from (80), that Let Ci ; Ci

  ð1Þ1 ð2Þ1 ð1Þ Ci  Ci : Ci ¼ OðD2 Þ:

ð81Þ ð1Þ1

By X we denote the orthogonal projection of Ci tangential space T Cð1Þ1 M. Using (81), we get for X

ð2Þ1

 Ci

on the

2 3   1 6 ð1Þ1 7 ð1Þ ð1Þ1 ð2Þ1 ð2Þ1 ð1Þ : Ci  X ¼ Ci  Ci  4 Ci  Ci  5C  ð1Þ 2 i Ci 

h

  i    Y : sym B01  B02 : Y ¼ kmax sym B01  B02 ;

   where kmax sym B01  B02 is the maximal eigenvalue of the sym  0 metric operator sym B1  B02 : Sym ! Sym. Obviously, the same n o maximal eigenvalue has its restriction on T 0 ¼ Span B01 ; B02 . It can be easily seen that

       sym B01  B02 B01 ¼ 1=2 B01 : B02 B01 þ 1=2 B01 : B01 B02 ;       sym B01  B02 ðB02 Þ ¼ 1=2 B02 : B02 B01 þ 1=2 B01 : B02 B02 :

i

¼

ð1Þ1 Ci



ð2Þ1 Ci

2

þ OðD Þ:

      ð1Þ 1=2 2 ð1Þ  , we have Moreover, since tr XCi X ¼  X C i  

Therefore, the n o matrix of the restricted operator with respect to the basis B01 ; B02 has the following form:

ð82Þ

1 2

B01 : B02

B02 : B02

B01 : B01

B01 : B02

!

:

Both eigenvalues of A are real, since A represents a symmetric operator. Finally,

3

OðD Þ   ¼ OðDÞ: ð1Þ tr XCi X

A :¼

ð83Þ

     ^ Cð1Þ ¼ kmax sym B01  B02 ¼ kmax ðAÞ: q i

A.V. Shutov, R. Kreißig / Comput. Methods Appl. Mech. Engrg. 199 (2010) 700–711

Fig. 8. Function qðhÞ and its approximation by h3 =12.

  ð1Þ ^ Cð1Þ is a continuous function of Ci . Therefore, the Note that q i ð1Þ  ^ðC Þ is well defined. We commaximum qðhÞ ¼ max ðCð1Þ1 ÞD 6h q i

pute it by the brutal force method. Moreover, the following parametrization can be used to simplify the computations. For any tensor ð1Þ Ci there exist a Cartesian coordinate system and real numbers ð1Þ 1 k ; k2 > 0 such that the matrix of Ci takes the diagonal form 1 2 1 1 diagðk ; k ; 1=ðk k ÞÞ. The function qðhÞ is plotted in Fig. 8 for 0 6 h 6 0:2. The following simplified estimation approximates the function qðhÞ with a high precision for 0 6 h 6 0:2 (cf. Fig. 8)

qðhÞ  h3 =12:

ð84Þ

References [1] H.D. Alber, Materials With Memory. Initial-Boundary Value Problems for Constitutive Equations with Internal Variables, Springer, 1997. [2] N.S. Bahvalov, N.P. Jidkov, G.M. Kobelkov, Numerical methods, Lab of base knowledge, 2001 (in Russian). [3] A. Bertram, Elasticity and Plasticity of Large Deformations, Springer, 2005. [4] W. Dettmer, S. Reese, On the theoretical and numerical modelling of Armstrong–Frederick kinematic hardening in the finite strain regime, Comput. Meth. Appl. Mech. Engrg. 193 (2004) 87–116. [5] E. Hairer, Geometric integration of ordinary differential equations on manifolds, BIT Numer. Math. 41 (5) (2001) 996–1007. [6] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure Preserving Algorithms for Ordinary Differential Equations, Springer, 2006. [7] W. Han, B.D. Reddy, Plasticity. Mathematical Theory and Numerical Analysis, Springer, 1999. [8] S. Hartmann, G. Lührs, P. Haupt, An efficient stress algorithm with applications in viscoplasticity and plasticity, Int. J. Numer. Meth. Engrg. 40 (1997) 991– 1013. [9] S. Hartmann, P. Neff, Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility, Int. J. Solid Struct. 40 (2003) 2767–2791. [10] S. Hartmann, K.J. Quint, M. Arnold, On plastic incompressibility within timeadaptive finite elements combined with projection techniques, Comput. Meth. Appl. Mech. Engrg. 198 (2008) 173–178.

711

[11] P. Haupt, Continuum Mechanics and Theory of Materials, second ed., Springer, 2002. [12] D. Helm, Stress computation in finite thermoviscoplasticity, Int. J. Plasticity 22 (2006) 1699–1721. [13] I.R. Ionescu, M. Sofonea, Functional and Numerical Methods in Viscoplasticity, Oxford Mathematical Monographs, first ed., Oxford University Press, 1993. [14] M. Itskov, Tensor Algebra and Tensor Analysis for Engineers: With Applications to Continuum Mechanics, Springer, 2007. [15] J. Kato, A.A. Martynyuk, A.A. Shestakov, Stability of motion of nonautonomous systems (method of limiting equations), Stab. Control: Theory Meth. Appl., vol. 3, Gordon and Breach Publishers, 1996. [16] E. Kröner, Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen, Arch. Rat. Mech. Anal. 4 (1959) 273–334. [17] E.H. Lee, Elastic–plastic deformation at finite strains, J. Appl. Mech. 36 (1969) 1–6. [18] A. Lion, Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological elements, Int. J. Plasticity 16 (2000) 469–494. [19] G. Lührs, S. Hartmann, P. Haupt, On the numerical treatment of finite deformations in elastoviscoplasticity, Comput. Meth. Appl. Mech. Engrg. 144 (1997) 1–21. [20] R. Michel, R. Kreißig, H. Ansorge, Thermomechanical finite element analysis (FEA) of spin extrusion, Forschung im Ingenieurwesen 68 (2003) 19–24. [21] C. Miehe, E. Stein, A canonical model of multiplicative elasto-plasticity: formulation and aspects of the numerical implementation, Eur. J. Mech. A – Solid 11 (1992) 25–43. [22] C. Miehe, Exponential map algorithm for stress updates in anisotropic multiplicative elastoplasticity for single crystals, Int. J. Numer. Meth. Engrg. 39 (1996) 3367–3390. [23] P. Neff, Mathematische Analyse multiplikativer Viskoplastizität, Ph.D. Thesis TU Darmstadt, Shaker Verlag, 2000. [24] P. Perzyna, The constitutive equations for rate sensitive plastic materials, Quart. Appl. Math. 20 (1963) 321–331. [25] P. Perzyna, Fundamental problems in visco-plasticity, in: G. Kuerti (Ed.), Advances in Applied Mechanics, vol. 9, Academic Press, New York, 1966, pp. 243–377. [26] S. Reese, D. Christ, Finite deformation pseudo-elasticity of shape memory alloys – constitutive modelling and finite element implementation, Int. J. Plasticity 24 (2008) 455–482. [27] A.V. Shutov, R. Kreißig, Finite strain viscoplasticity with nonlinear kinematic hardening: phenomenological modeling and time integration, Comput. Meth. Appl. Mech. Engrg. 197 (2008) 2015–2029. [28] A.V. Shutov, R. Kreißig, Application of a coordinate-free tensor formalism to the numerical implementation of a material model, ZAMM 88 (11) (2008) 888–909. [29] J. Simo, T. Hughes, Computational Inelasticity, Springer, 1998. [30] J.C. Simo, C. Miehe, Associative coupled thermoplasticity at finite strains: formulation, numerical analysis and implementation, Comput. Meth. Appl. Mech. Engrg. 98 (1992) 41–104. [31] B. Svendsen, A logarithmic-exponential backward-Euler-based split of the flow rule for anisotropic inelastic behaviour at small elastic strain, Int. J. Numer. Meth. Engrg. 70 (2007) 496–504. [32] C. Truesdell, W. Noll, The Non-Linear Field Theories of Mechanics, second ed., Springer, 1992. [33] G. Vadillo, R. Zaera, J. Fernandez-Saez, Consistent integration of the constitutive equations of Gurson materials under adiabatic conditions, Comput. Meth. Appl. Mech. Engrg. 197 (2008) 1280–1295. [34] R.Z. Valiev, R.K. Islamgaliev, I.V. Alexandrov, Bulk nanostructured materials from severe plastic deformation, Prog. Mater. Sci. 45 (2000) 103–189. [35] G. Weber, L. Annand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoelastic solids, Comput. Meth. Appl. Mech. Engrg. 79 (1990) 173–202.