Sensitivity analysis in plane stress elasto-plasticity and elasto-viscoplasticity

__ .-__ l!!J

2s

Computer methods in applied mechanics and englneerlng

4._

ELSEVIER

Comput.

Sensitivity

Methods

Appl.

Mech. Engrg.

137 (1996) 39549

analysis in plane stress elasto-plasticity elasto-viscoplasticity

and

Michal Kleiber *, Piotr Kowalczyk Institute of Fundamental

Technological Research, Polish Academy of Sciences, ul. $wi@okrzyska 21, 00-049 Warsaw, Poland Received

4 December

1995

Abstract Development of computational models for plane stress elasto-plasticity is known to make some difficulties non-existent in 3D and plane strain formulations. Tedious derivation of plane stress algorithmic tangent matrices which cannot be directly obtained as a special case of 3D matrices may serve here as an example. In this paper specific features of plane stress computations are discussed in the context of parameter sensitivity analysis. Elasto-plastic material model with isotropic/kinematic hardening as well as two viscoplastic models are considered. Explicit expressions for sensitivity gradients are derived. Instructive examples solved by FEM illustrate the paper.

1. Introduction

Sensitivity analysis for nonlinear structural mechanics problems has been emerging in the last decade as a fruitful area of engineering research. Theoretical formulations and computational algorithms developed for nonlinear system sensitivity, while still rarely effective enough to deal with the majority of practical problems, have now reached a level of maturity which allows one to solve numerically many illustrative examples that are very helpful in gaining further insight into the nature of nonlinear system sensitivity. Among the contributions which have been published so far on issues related to sensitivity of inelastic response one should mention [l-lo]. Basic aspects of the approach put forward in [4,5,9,10] can be summarized as follows: (i) Equilibrium (primal) analysis of the rate elastic-plastic problem is carried out in the typical time (or pseudo-time) interval [t, t + At] by using a specific integration scheme. Typically, the constitutive law is integrated using the radial return approach combined with appropriate linearization technique, resulting in the so-called consistent (or algorithmic) constitutive tangent matrix. (ii) Sensitivity evaluation follows as an add-on to such an analysis, i.e. sensitivity equations are derived on the basis of the time-discretized equilibrium problem. It turns out that the stiffness matrix to be used at the (non-iterative!) sensitivity stage is the one based on precisely the same algorithmic tangent consistent with the time integration algorithm used for integrating the rate equilibrium problem. Another consequence is that the same time step has to be used at both equilibrium and sensitivity stage of the analysis-not really a desired situation in view of typically more severe accuracy-implied step size limitations for sensitivity analysis as compared to equilibrium analysis. There is an alternative to the above approach. In it, sensitivity derivations are performed on the rate (rather than incremental) equilibrium problem and only then the resulting rate sensitivity formulation is subject to time integration by using any available algorithm (not necessarily the one used for integrating *Corresponding

author.

0045.7825/96/$15.00 @ 1996 Elsevier PII SOO45-7825(96)01072-9

Science

S.A. All rights reserved

396

M. Kleiber, P Kowalczyk/Cotnput.

Methods Appl. Mech. Engrg. 137 (1996) 395409

the rate equilibrium problem). The step size can also be selected independently for both the equilibrium and sensitivity analysis to account for possibly different accuracy indications. The above discussion is summarized in Fig. 1, in which q is the vector of nodal generalized coordinates while DDM stands for the direct differentiation method employed’ . As a matter of fact, the situation represented in Fig. 1, which takes reference to the semi-discretized formulation of elasto-plastic problem, is just a special case of a more general situation shown in Fig. 2. There, two alternate routes (A) and (B) are presented indicating how the original rate-type continuum based equilibrium problem can be transformed to yield time- and space-discretized sensitivity formulation. In regard to Fig. 1 we observe that in this article: (i) the semi-discretized formulation will be used as a starting point and (ii) the route (B) is followed in deriving time- and space-discretized sensitivity equations, which is in accord with our earlier papers [9,10]. Only to this approach will all our specific observations apply; nevertheless, many of the conclusions are believed to have a general significance. In spite of the progress so far, there are clearly still very many unsolved problems in the field of inelastic system sensitivity which need to be addressed. One of the issues not yet entirely clarified is

discretization procedures (in space

and time)

continuum-based

design

design

differentiation

differentiation

Fig. 2.

‘In fact, as documented in [9], the adjoint system method that only the latter is subject to scrutiny in this article.

is unlikely

to perform

better

than DDM

for inelastic

problems

so

M. Kleiber, P Kowalczyk/Cotnput.

Methods Appl. Mech. Engrg. 137 (1996) 395409

397

related to somewhat ‘anomalous’ behaviour of design gradients at elastic-plastic transition points or, more precisely, to only directional properties of this behaviour at some points along the solution path. The problem was addressed in [8,10] and will be further taken up in this paper. Another unresolved problem (which may seem simple but in fact is not!) is the problem of sensitivity of plane stress elastoplasticity. The reason for this situation is the nature of return mapping algorithms now in standard use for solving elastic-plastic equilibrium problems. Such algorithms may be trivially modified for the plane strain problem but not at all for the plane stress problem in which satisfaction of the zero stress condition places a non-trivial constraint on the return algorithm. The problem was elegantly solved in [ll]; here we shall derive the appropriate consistent tangent using a slightly different approach, extend it to some viscoplastic models and then discuss the implications of the formulation for the study of sensitivity.

2. Finite element formulation Finite element aspects of the sensitivity analysis in the 3D case have been presented in [4,10]. Here, we only summarize the essential points of the approach, modifying it to apply to the plane stress problems. In the equilibrium finite element analysis the following equation is solved at the end of each time step, typically [t, t+At],

J

B;,q;dV=Qa,

c~=1,2 ,...,

N,

i=1,2

V

,...,

6,

(1)

where Bi, is the strain-nodal displacement matrix, a; denotes the stress vector related to nodal displacements qa by certain constitutive and geometric relations and Qu is the nodal load vector. N denotes the number of nodal parameters. For history-dependent problems the equation is usually formulated in an incremental form as

J

Bi, Aui dV = r+A’Qa- J

V

Bi, ‘vi dV 3 V

(2)

where Aai=AUi(‘uj,‘pl;,A&j;h),

1C=l,2,...

,K,

(3)

A&i = B;, Aqa denotes the strain vector, rpi is a set of internal parameters (like back stress or equivalent plastic strain) representing the state of the material, and h denotes the design parameter with respect to which the sensitivity is to be analyzed (for simplicity we restrict ourselves in this paper to considering a single parameter only; the subsequent considerations can easily be extended to a vector of such parameters h, though). The particular form of the incremental constitutive equation (3) must be derived from the rate-type constitutive equation for the material considered using the time integration scheme employed. For plane problems Eq. (2) can also be written as

J

J

b Bi, Aci da = ‘+” QabBi, ‘ai dR R R where b is the reference thickness of the element and fl its reference area (both assumed constant here). Obviously, index i in Eq. (4) runs over a reduced set of values corresponding to in-plane stresses only. In order to solve for Aq, this generally nonlinear set of equations an iteration scheme is needed. If the Newton-Raphson algorithm is applied, the subsequent corrections to the already obtained solutions are computed from equations with the following governing matrix

Km, =

JV

Bi, eij Bjp dV )

The matrix Cij is known as the algorithmic (or consistent with the time integration scheme) tangent matrix and will be shown below to have a fundamental importance in the sensitivity analysis. It may be shown that in order to compute design derivatives of any functional of our interest we have to know first the design derivative of the incremental nodal displacement vector, obtained as the solution

398

M. Kleiber, II Kowalczyk/Comput.

Methods Appl. Mech. Engrg. 137 (1996) 395-409

of the equilibrium problem (4). Its design differentiation

(.I

Bi, C?iiBip dV

V

yields

-d&p > dh

(‘5) As we can see here, we have to do with a linear set of equations with the same governing matrix analysis. In order to form the r.h.s. vector some additional derivatives of incremental stresses (3) have to be computed and design increments of stress and internal parameters have to be accumulated at each integration point. K,p as in the last iteration of the equilibrium

3. Plane stress problem formulation The idea of equilibrium computations in the case of plane stress is clearly the same as in the general 3D case. There, the computations are based on the return mapping algorithm and iterative scheme involving the consistent elasto-plastic tangent moduli obtained by linearization of the algorithm. However, due to a somewhat numerically inconvenient constraint of zero stress in the prescribed direction (which results in the fact that the strain component in this direction cannot be explicitly computed from the displacement field), the procedure that leads to the solution looks quite different. The general formulation of the return mapping algorithm follows here the one presented by Simo and Taylor [ll]. Assume that the tensors of stress 0 and strain E in the plane stress problem are represented as vectors {a,,, a,,, &gX,,} and {cXx, +,, X&E,,}, respectively. The factor fi is used here for convenience of representing some scalar product operations. The component a,, equals zero while gZzzhas to be computed ‘additionally’. Elastic isotropic constitutive equation in this notation reads

u=DE,

(7)

and the transversal elastic strain is expressed as &

zz=-L(& 1-V

+Eyy) = -~(c& +uyy).

In the elasto-plastic range strain is the sum of the elastic and plastic components, E =

Ee + EP

and only the elastic one fulfills Eqs. (7)-(S). The plastic incompressibility plastic strain, EZz= -(El, + $),

(9) condition is imposed on the (W

i.e. the plastic strain is a deviatoric tensor. A rate-type constitutive equation is postulated for Ed, its form being dependent on the particular constitutive model considered. Here, we present the formulae for the following three material types: l Rate-independent elasto-plasticity:

where

M. Kleiber, I! Kowalczyk/Comput.

z=

I - flol)(u-a).

@=

@-K=O,

$

k = K’iP,

Methods Appl. Mech. Engrg. 137 (1996) 395-409

3

II4

399

(12) (13)

>

2 z $+p,

(14)

with ~‘(8’) and q’(8) being given hardening functions, c” the equivalent plastic strain rate and (Ythe back stress. Note that the evolution equation for (Ydefines actually only its deviator rate*, while the volumetric part can be computed from the plane stress condition cr,, = 0. 0 Overstress viscoplasticity [13]:

where

(16) and the positive monotonic overstress function 4 may be assumed as, for instance, (17)

l

Here, p is the so-called viscosity parameter while m is another material constant. Power-law strain and strain-rate hardening viscoplasticity [14]:

f(ti,e”)

=

io

’ (-E-- )lJrn dep)

g(P) = E&O

l/n

I+”

( go>* (!I )

(18)

1+

4

Here, .s~ and go are scalar reference values of plastic strain and its rate, respectively, while m, II are additional dimensionless material constants. E denotes, as before, the elastic Young modulus. The deviator of the plane stress tensor can be expressed as

(19)

Note that the vector s represented by {sXx,s,,r, &s,,} contains only in-plane components of the stress deviator. The non-zero transversal component sZZcan be computed easily. It is straightforward to check that the deviator norm in the 3D space (i.e. sjX + s;~ + s:~ + 2s$) can be expressed as X/X. The following formulae will be helpful in performing further transformations. orthogonal operator

*In most papers (and, consequently, its volumetric part deal with the back

Let Q be the following

devoted to 3D elasto-plastic computational algorithms (see [12.10], for instance) the notion of the ‘back stress’ the symbol a) is actually related to the deviator which appears in Eq. (14). There is then no need to compute at all. Here, the volumetric part of the back stress is introduced only for the notation simplicity (enabling us to stress as a ‘plane stress’ object). The notion appears somewhat ‘artificial’ in the context of its physical meaning.

400

M Kleiber, R KowalczyklComput.

i-

Djl =

i

--

Methods Appl. Mech. Engrg. 137 (1996) 395-409

h ho 1 Jz A O. --

0

1

1

0

We can then write

J

(20)

P=QpQT,

(21) [Bij] = diag

D=QbQT,

?.-

-.&

E

1-Y’

l+Y’

l+v I

.

(22)

3.1. Return mapping algorithm To solve the problem of elasto-plastic scheme:

deformation

in time we introduce

the following integration

u=‘u+Au=‘u+D(A~-Acp)

(23)

AE~= hP&,

(24)

where &=U--(Y,

(25) A=

AZQ

3655~ =T’

@=

(26)

J---STPS The scheme above is valid for both rate-independent elasto-plasticity and elasto-viscoplasticity. In the latter case no kinematic hardening is here assumed, i.e. (Y= 0, 77’= 0 and g = u. Combining Eqs. (23)(26) yields

where the ‘trial’ elastic stress (explicitly known for a given in-plane strain) has been defined as 5” = ‘u - ‘~+DAE. Eq. (27) may be transformed z&J = 5” )

(28)

to the form

Z=yJ+hDP,

y, =1+&j

(29)

which, after inverting Z, allows to find the value of 6 at known A. The value of A (or At?Q,as both are related by Eq. (26)) must be computed from the equation of the type (30)

F(A, 5) = 0 which for rate-independent

elasto-plasticity

reads

F=c?-K(p)=0

(31)

while for the viscoplastic models considered in this paper it becomes F=Atf((T,t?)-At?=O.

(32)

Since the equation is nonlinear in A, the Newton-Raphson scheme is typically used to find the solution. Subsequent corrections aA to an approximate solution A(‘) (the latter assumed to have been already

M. Kleiber, II Kowalczyk/Comput.

obtained--i

Methods Appl. Mech. Engrg. 137 (1996) 395-409

401

is the iteration number) are computed from the equation A(‘+‘) = A(‘) + 6A(‘)

(33)

with dF(A(‘)) pzzz dh where for the rate-independent ,=l-;Ai’

>

(34) elasto-plasticity

ii’ = K’

(35)

while for the viscoplastic models considered in this paper af

2 hi? ,

K’ = 1 - At af (36) 8ep ’ respectively. The value of d&/dA is needed in both the cases. We substitute it from Eq. (29) differentiated with respect to A to obtain y2

=

At

ab

-

3

(37) which can be solved for the unknown dg/dA. (Dependence of q’ on 8’ and thus on A has been neglected for simplicity but it may also be considered without affecting the general idea of the algorithm.) The result, when substituted in Eq. (34), allows to solve for the correction SA(‘) and Eq. (37) can be formed again based on the new values of A and 6 (the latter computed from Eq. (29)). After reaching convergence the final stress value may be computed as u = g + LY. In order to avoid frequent inverting of the generally non-symmetric matrix 2 one may use transformations (20)-(22). Let us denote (38) It can easily be shown that (39) By substituting Eqs. (20)-(22) into Eqs. (29) and (37) and premultiplying following sets of equations

them by QT we obtain the

Z&5”,

(40)

(41) in which i?=y,I+Ab~

(42)

is a diagonal matrix and $’ dAfQTg”. It can be shown that after some algebraic manipulations the solution 8 can be expressed as (43) while dt/dA

as

dg

-&I?

dh=

1+AK’

-5;G l+hc’

l+AG

(5xX+ Syy)K (-6Xx + S&G 2&G l+AZ? ’ l+AG z-z

(44)

M. Kleiber, II KowalczyklComput.

402

Methods Appl. Mech. Engrg. 137 (1996) 395-409

where K=

E

277’

E ;3

&

giQq+3

l+v

(45)

3 .

The value of 5 is given as t=Q& The consistent tangent matrix for the above time integration subsection.

scheme will be derived in the next

3.2. Design variations In order to obtain the expression for dAo/dh as a linear function of the variation dAs/dh (which is necessary to perform the sensitivity analysis of the equilibrium solution) we have to differentiate with respect to h all the equations developed in the previous subsection. By using the notion of total (8) and partial (a) design variations, which formally replace the total and explicit design derivatives according to

s(.) =

4.)

ah,

dh

a(-) = $

Sh,

(46)

we can write, see Eqs. (23)-(26), ~A~==~~-~~~+SA~,

(47)

$Acu= ; (6A~‘J+hdq’~+h&&,

(48)

1.

from which we obtain St=;

[

sy+sA~-~(Shl7’+laq’)~

Further, by differentiating zsg+;

(49)

Eq. (29) we obtain (; $I+DP

hdn15+SA

) ~-~D(AE-AP&)-~‘~-D&AE=O _

(50)

A&P

where 80 is represented

as

(51) Combining it with Eq. (49),and premultiplying

(I-l+;

by D-’ yields

P) $A

Pt -D-l

The unknown variation &Acan be derived by differentiation model considered. 3.2.1. Rate-independent

~D(AE -AE~) - $AE = 0.

(52)

of Eq. (30) for any specific constitutive

elasto-plasticity

By design differentiating

Eq. (30) with Eqs. (31) and (26) we have

~f~Ptr~_i3K__~~r~_~OK18h=0

(53)

with y2 defined by Eq. (35). Substituting in it Eq. (49) and solving with respect to $A yields

6h =

3’)‘~ 2(K’+?7’)

gTP(fi’J.+&Aa) -0-V

1 K’+7f

which by utilizing Eq. (26) and observing that

y2h

87’

+ $$

(8K

+

K’s

(54)

‘2’) I

K’+?f

=

3/(2h) (7, - y2) becomes

M. Kleiber, I? Kowalczyk/Comput.

Methods Appl. Mech. Engrg. 1.37 (1996) 395-409

403

(55)

After substituting it into Eq. (52) and some rearrangements

A

-

(y,

aK -y2)a(K’6 ‘7

+

+ A?

a$)P&

-D-l

dD(A&

we get

- A&)

-

6AE= 0,

where the notation (PC) ~3(PJ) stands for the tensor product, i.e. Pt(PoT. solved for 6Aa yielding the following result

Eq. (56) may easily be

(57)

where A dAu = (n _ y2)a (aK + AP a@Pg

+ i?D-’

dD(A& - A&) ,

(58) (59) (60)

dAu -& aA&

(61)

)

and

(62) is the consistent (algorithmic) tangent matrix for the integration scheme considered. REMARK. The above form of the consistent matrix i? looks quite different (and a bit more simply) than that obtained in [ll]. They are fully equivalent, though, which may be shown by multiplying &-’ in Eq. (62) by its counterpart & derived in [ll] to get the result which is equal to I. Unfortunately, no use of the transformations (20)-(22) can be made at this stage of the analysis-the matrix product (Pt)@(Pg) will never have a diagonal form. Having derived Eq. (57) we are able to form and solve the global sensitivity problem for SAq and determine ~A.E at the integration points. Then, to perform the calculations at the next time step, design variations of certain parameters must be updated. By knowing 6Au we can calculate sg and 6A making use of Eqs. (49) and (55), respectively. Then, by using Eq. (48), we can update the design variation of the back stress (Ywhile $A? may be computed from (63) see Eq. (26). 3.2.2. Viscoplasticity In the viscoplastic models considered Eqs. (49) and (52) become $5 = &a = $Au + &,

no kinematic hardening is assumed. Thus, (Y= 0, 7’ = 0 and (64)

404

M. Kleiber, l? KowalczyklComput.

Methods Appl. Mech. Engrg. 137 (1996) 395-W

(D-’ + AP) 6Aa + 6h Pa + AP 6la - D-’ aD(Ae -A&‘) - RAE = 0. By differentiating

(65)

Eqs. (32) and (26) and applying notation (36) we get

yz6ir+(l-r~-‘)G’d-~K’r3-Sh+Ati)f=O,

(66)

which can be solved for 6h and substituted into (65). By observing that + 6Ac-r)

$a = &‘P(8-u we can write

‘u+sAa)+(l-i?)k?P+At6’f and +hP + 4ii’a2 9y2 (Pa) 8 (Pu)] $A@ + [,p + -& +

3(1 - ii’) s,_pPu

(67)

1 (Pa) @ (Pu)] 8’0

+ 3At G’f PU - D-’ dD(Ac -a~~) - RAE = 0.

23a

2K'(+

(68)

(69)

The latter equation can be solved for &Au yielding (70)

where dAu = -z

i?Pu + eD_’ ~D(AE-AE~)

aau -=eD-’

(72)

-1

815 aAu

-

a k3

=-

3(1 - K’) .. 2._,e

CPU

aAu -e dAE

(71)

)

(73) (74)

and

is the consistent (algorithmic) tangent matrix for the constitutive model and the integration scheme considered. As in the rate-independent case, Eq. (70) allows to build and solve the global sensitivity problem for 6Aq and determine $AE at the integration points. Then, to perform the calculations at the next time step, design variations of u and S’ must be updated. This can be done by utilizing Eqs. (70), (68) and (63). 4. Computational examples A plane stress finite element model of a uniaxially stretched elasto-plastic rectangular strip with a circular hole was considered. The finite element mesh and nodal boundary conditions for displacements as well as applied loads (only a quarter of the strip is considered due to symmetry conditions) are

M. Kleiber, I? Kowalczyk/Comput.

Methods Appl. Mech. Engrg. 137 (1996) 395-409

405

presented in Fig. 3. The values of the load forces were increased linearly in time t according to the formula FCn,= Fcn,t, with the values of F(n) listed in Table 1. The three material models discussed in the previous sections have been examined. For all the three cases the values of elastic constants E and v were assumed as 2.01 x10” N/m2 and 0.25, respectively. For the rate-independent and overstress models a linear hardening function K(i?)

=

Uy +

K'p ,

K’

=

const.,

(76)

was adopted. The remaining material constants assumed for the models were: l rate-independent elasto-plasticity: uY = 2 x lo9 N/m2, K’ = 6x lo9 N/m2, v’ E 0, l overstress viscoplasticity: aY = 2 x lo9 N/m2, K’ = 6 x 10” N/m2, p = 0.5 s, m = 2, l power-law strain and strain-rate hardening viscoplasticity: co = 0.01, E. = 0.01 SC’, m = 0.25, n = 2.5. The strip thickness b = 1 mm was assumed in all the cases.

Table 1 Reference

nodal

load values

(n)

FC,I [N/s1

1 2 3 4 5 6 7

1000.0 2060.2 2245.6 2516.8 2832.1 3103.0 1632.6

17mm il/_L

Fig. 3. Finite element

c 34mm

model of the stretched

strip with a circular

hole.

406

M. Kleiber, I? Kowalczyk/Comput.

Methods Appl. Mech. Engrg. 137 (1996) 395-409

For each material model two loading tests were performed: one with a coarse and another with a fine time discretization (At = 0.05 s and 0.001 s, respectively). The tests were completed when the load exceeded by a few per cent the level at which plastic flow occurred in the total strip cross-section (for the models with the yield limit). The deformation of the specimen boundary and the development of the plastic zones for the case of rate-independent elasto-plasticity is shown in Fig. 4. In the subsequent figures presenting results, individual symbols correspond to the solution obtained with the coarse time discretization while solid lines represent the solution obtained with the fine disretization. Horizontal displacement of the point A (see Fig. 3) was monitored and is presented in Fig. 5. Sensitivity gradients of this displacement with respect to variations of thickness and various material parameters were calculated and are shown in Figs. 6-8. The results for different design parameters are scaled by multiplying them by the nominal values of the respective parameters-the so-scaled values become then the measure of sensitivity with respect to relative changes of the parameter values. This makes it easy to compare the influence of changes in different design parameters on the specimen response. Similarly to the 3D case (see [lo]), it can be seen that the sensitivity graphs for rate-independent material behaviour are discontinuous. Small but numerous jumps that refer to crossing the yield limit by stresses at subsequent integration points are clearly visible for the fine time discretization test (they are hardly noticeable when larger At is chosen), especially for h = b and h = cry, The problem of the discontinuities was discussed recently in the papers [C&8];its more detailed discussion in the context of the return mapping algorithm can be found in [lo]. Speaking briefly, a discontinuity in the sensitivity gradient may occur at the very yield limit point if its value is different for different design variations, depending on whether the structure perturbed by such a variation is in the state already after or still before crossing the yield limit at a certain time integration point. Then, the left-hand limit value at the discontinuity point refers to the perturbed states that are still elastic while the right-hand limit value describes the sensitivity with respect to the design variations at which the structure is already plastic. Such discontinuities cause no substantial difficulties in the approach presented in this paper. The equations derived in the previous sections properly account for the jumps caused by possible discontinuities occurring within the time step. A problem may only appear if a discontinuity occurs at the very end of

h1.70

t=1.30

t=1.10

‘\

:;, I&

k2.10

:::: :.: L.:

:. J t-2.50

,’

,/.:

.:

.:

.:;::.:.. ..:::,

:::.

‘.

:::

:.

‘_

,::,:...::;

.;.

L--‘.,:i;:

::,.

e, :::::

Fig. 4. Deformation of the strip and plastic zone evolution (dots denote Gauss independent elasto-plasticity, At = 0.05, displacements scaled by the coefficient

::..

G;,j;’ ::::, -.:: ..:.::

integration 5.0.

points

where

plastic

flow occurs);

rate

h4. Kleiber, l? Kowalczyk/Comput.

Methods Appl. Mech. Engrg. 137 (1996) 395409

407

a time step-then either one of the two possible sensitivity limits should be taken as the instantaneous sensitivity, in correspondence with the assumed state of the material (elastic or plastic) at this time instant (the value of the sensitivity gradient output as a result will then be valid only for a certain class of design variations-such that the structure perturbed by these variations is at the considered time instant in the elastic or plastic state, respectively). Consequently, depending on the assumption, the jump caused

0.0006

ate-independenl

0.50

(Huber-Mises) model

1.oo

1.50

2.00

2.50

time I [s] Fig. 5. Horizontal

E -c 3

displacement

of node A (absolute

value).

0.003

5 s

.s 2

.z

0.002

i E E 8 m

0.001

P e

0.000 0.50

1 .oo

1.50

2.00

2.50

pseudo-time f [s] Fig. 6. Sensitivities

of the horizontal

displacement

of node A, rate-independent

(Huber-Mises)

elasto-plasticity.

408

M. Kleiber, I! Kowalczyk/Comput.

Methods

Appl. Mech. Engrg. 137 (1996) 395409

0.0010 r-

0

h=b

cl

h-E

-E: s^

A V

h = oy

s s

X

h=p

0.0006

0 .z ._ i! w %

0.0004

E

0.0008

h=m

A..

j % 5

0.0002

_...

L

-7-z-F i-

0.0000

0.00

0.50

/ , , 1.00

1.50

2.00

2.50

time 1 [s] Fig. 7. Sensitivities of the horizontal displacement of node A, overstress viscoplasticity.

i 0.00

0.50

1.00

1.50

2.00

2.50

time f [sl Fig. 8. Sensitivities of the horizontal displacement of node A, power-law strain and strain-rate hardening viscoplasticity.

M. Kleiber, P Kowalczyk/Comput.

409

Methods Appl. Mech. Engrg. 137 (1996) 395-409

by a discontinuity occuring at this point is accounted for either at the current or at the next time step. The uncertainty about the result at such a point has thus no influence on the results at subsequent time instants corresponding to local loading. As a matter of fact, it is highly unlikely that such a situation will ever appear in the computational practice. Moreover, it may be shown that possible errors caused by choosing the wrong limit value of the sensitivity gradient at a discontinuity point for an arbitrary design variation are not larger than those encountered when exact (continuous) results obtained in a close neighbourhood of such a discontinuity are checked against finite design perturbations. Note that in computational applications of sensitivity results we always have to do with finite perturbations not necessarily complying with the perturbation ‘smallness’ assumption adopted when developing sensitivity equations. If the equilibrium solution of the perturbed system happens to be on the opposite side of this discontinuity (for instance, if the primary structure is elastic and the perturbed one is plastic), the sensitivity results are erroneous to the same extent as those at the discontinuity point for infinitesimal variations.

5. Conclusions The idea of sensitivity evaluation in the plane stress elasto-plastic as well as elasto-viscoplastic problems appears to be essentially the same as in the 3D case (see [lo]). Obviously, due to a different type of the stress constraint condition the formulae derived in this paper differ from those for 3D case and the computation of the consistent tangent matrix e seems more complicated and tedious. The algorithm for sensitivity computations can easily be implemented in a finite element code. The computational examples show similar characteristic features of the sensitivity graphs to those typically obtained in the 3D case (discontinuities at yield limit points, for instance). References [l] J.J. Tsay and J.S. Arora, Nonlinear structural design sensitivity analysis for path-dependent problems. Part 1: General theory, Comput. Methods Appl. Mech. Engrg. 81 (1990) 183-208. [2] J.J. Tsay, J.E.B. Cardoso and J.S. Arora, Nonlinear structural design sensitivity analysis for path-dependent problems. Part 2: Analytical examples, Comput. Methods Appl. Mech. Engrg. 81 (1990) 209-228. [3] Q. Zhang, S. Mukherjee and A. Chandra, Design sensitivity coefficients for elasto-viscoplastic problems by boundary element method, Int. J. Numer. Methods Engrg. 34 (1992) 947-966. elasto-plasticity, Comput. Methods Appl. Mech. [41 C.A. Vidal and R.B. Haber, Design sensitivity analysis for rate-independent Engrg. 107 (1993) 393-431. sensitivity analysis for problems with any material and kinematic nonlinearity, Comput. PI M. Kleiber, Shape and non-shape Methods Appl. Mech. Engrg. 108 (1993) 73-97. D.A. Tortorelli and C.A. Vidal, Tangent operators and design sensitivity formulations for transient nonlinear PI P Michaleris, coupled problems with applications to elastoplasticity, Int. J. Numer. Methods Engrg. 37 (1994) 2471-2499. structures, Int. J. Numer. Methods Engrg. 37 (1994) I71 M. Ohsaki and J.S. Arora, Design sensitivity analysis of elasto-plastic 737-762. method for design sensitivity analysis of elasto-plastic 181TH. Lee and J.S. Arora, A computational Appl. Mech. Engrg. 122 (1995) 27-50. and P. Kowalczyk, Parameter sensitivity of elasto-plastic I91 M. Kleiber, T.D. Hien, H. Antunez

structures, response,

Comput. Engrg.

Methods

Comput.

I2

(1995) 263-280. Constitutive parameter sensitivity in elasto-plasticity, Comput. Mech. 16 (1995) 297-306. [lOI M. Kleiber and P. Kowalczyk, for plane stress elastoplasticity, Int. J. Numer. Methods Engrg. 22 illI J.C. Simo and R.L. Taylor, A return mapping algorithm (1986) 649670. for rate-independent elastoplasticity, Comput. Methods Appl. Mech. PI J.C. Simo and R.L. Taylor, Consistent tangent operators Engrg.

48 (1985) 101-118. The constitutive equation for rate-sensitive plastic C.E Shih and A. Needleman, A tangent modulus

I131 P Perzyna. P41 D. Peirce, 875-887.

materials, Quart. Appl. Math. 20 (1963) 321-332. method for rate-dependent solids, Comput. Struct.

18 (1984)