Elastic-plastic boundary element analysis as a linear complementarity problem

Elastic-plastic boundary element analysis as a linear complementarity problem G. Maier and G. Novati Department Italy (Received

of Structural Engineering,

August

Technical University (Politecnico),

Milan,

1982)

Elastic-plastic constitutive laws allowing for nonassociation, ‘corners’ and hardening (or softening) are formulated both in the usual terms of rates and in terms of piecewise-linear approximation (local or global). The direct boundary element method combined with such descriptions of nonlinear material behaviour, is shown to reduce elasto-plastic structural analysis, both by infinitesimal and finite steps, to a linear complementary problem. On this basis, a general extremum characterization of the solution and conditions for solution uniqueness are pointed out. Various solution procedures are envisaged in the new unitary framework. Key words: boundary

Introduction Boundary element methods have only recently developed from their mathematical origin to being an important tool of engineering mechanics. Their increasing popularity is primarily due to some peculiarities, which are, in certain circumstances, computationally attractive compared to other well-established methods, primarily to the finite element methods. Among these features those frequently advocated are: simpler data preparation, smaller equation systems, selective and accurate evaluation of stresses and displacements in points internal to the domain, fitness for problems over infinite domains.1’2 Elastoplastic analysis by boundary integral equations and boundary elements has been the object of pioneering investigations by Swedlow, Cruse, Riccardella, Mendelson et al. in the ~OS.~,~In the last few years, some conceptual and mathematical difficulties and shortcomings have been overcome and corrected; the computational methodology of boundary elements in inelasticity have been notably developed and successfully tested in a number of papers.‘-” The contributions presented in this paper concern the following aspects of plastic analysis by the direct boundary element method. A broad class of elastoplastic constitutive laws (covering nonnormality, nonregular yield surfaces and hardening or softening behaviour) is first formulated in terms of rates; then it is given an approximate description in view of its use in evolutive (‘historical’) elastoplastic This paper was originally presented at the 4th International Conference on Boundary element methods in engineering, Southampton,

UK, 21-23 September 1982

elements,

elastoplasticity,

complementarity

analysis by finite steps. Both descriptions, developed from that proposed elsewhere, r2-15 have the form of a linear complementarity problem. Such mathematical structure is shown to be preserved when the boundary value problem in the finite step is subsequently formulated by the direct boundary element method. A similar formulation had been studied by Maier on the basis of other methods of discrete structural modelling. ‘6-18 In the present context, the complementarity formulation appears fruitful as a unifying framework for establishing extremum theorems and conditions for uniqueness of solution and overall stability. Moreover, a number of alternative mumerical techniques are related to the recognition of the step elastoplastic analysis by boundary elements as a linear complementarity problem. This paper is intended to briefly point out some direct consequences of this circumstance, which appear to have been overlooked so far.

Constitutive

laws

Elastoplastic incremental

laws at corners

A fairly general description of the elastic-plastic behaviour of materials in terms of rates (indicated by dots and to be interpreted as derivatives with respect to some increasing function t of the physical time) can be given the following traditional expression:r7*r9

(1) C%=~

(Oij,n) a.

a, Z 0

(~ = 1, . . . ,v>

(2)

4

0307-904X/83/02074-09/$03.00

74

Appl. Math. Modelling,

1983,

Vol.

7, April

0 1983 Butterworth &Co. (Publishers) Ltd

Elastic-plastic BE analysis: G. Maier and G. Novati for given stress rates tiij at a known stress state Uij after a known yielding history, can be expressed by the following relation set: if&
&,,a=0

@,fl=

I,...,y’)

(lob)

We will term the ‘inverse’ flow rule that which the incremental plastic response of the material above sense for given strain rates e, (it? in the imposed strain rates). The relevant relation set formulated by substituting into equation (lOa) sion for ciij provided by (1) and (2):

governs in the presence of is readily the expres-

(5) In equation (2) the ‘plastic potentials’ $, are functions of the current stress states and the past plastic histories (through some parameters denoted by 71,e.g. plastic work per unit volume), but are independent from incremental quantities (or rates); a potential is defined for each of the JJ ‘yield functions’ @,(cij, tr) or ‘yield modes’, which define in the stress space the current elastic domain (whose boundary is the current ‘yield surface’ or yield locus) by the conditions: $a(Uij,7T)GO (Ci=l,...f_Y)

(6)

The ‘plastic multiplier rate’ L is a nonnegative scalar variable which specifies the entity of the strain rate contribution provided by the oath yield mode, as equation (2) stipulates that this contribution is proportional to the stress gradient of function 11/,(in other terms is directed as the outward normal to the surface J/, = constant in the current stress point oij). In equation (3) the ‘hardening coefficients’ hap govern the way in which the incremental yielding process (defined by &) affects the elastic domain and the yield locus, i.e. express the instantaneous ‘hardening rule’. Their expression is obtained from the total differential of & as soon as the meaning of n is specified. E.g., if 71means plastic work density SO that ?i = Uije$, then: (7) Equation (4) expresses Prager’s ‘consistency rule’, according to which the yield mode (Ycontributes to the plastic strain rates (a, > 0) only if the relevant yield surface & = 0 contains the stress point Uij during the incremental process involving Goij = tiij 6t and 6eij = tij 6t. It is readily seen that the requirements stated in (4) can be equivalently expressed as follows: cp&=o &a,=0

(a=

l,...,y)

(11) kijrs being the tensor of elastic moduli; thus: kijrs ids-

i, = ~

~

kijrs ~

aUJj

ij

+ h,p

A+0

rs Wa)

&a,=0

$20

(q/3=

l,...,y’)

(12b)

For y’> 1 and 4, # IJ, the preceding flow rules concern, in the jargon of plasticity theory, a ‘corner’ of the current yield locus (i.e. in oij more than one yield mode is activatable) and a ‘nonassociative’ constitution (i.e. the ‘normality law’ is not fulfilled). The latter feature is exhibited by most nonmetallic materials: soils, rocks, concrete and in general materials with ‘internal friction’. The coverage of the above analytical descriptions can be further clarified by considering some special cases: (a) If y’ = 1: Oij is a ‘smooth’, regular point of the current yield surface, i.e. the tangent plane and the outward normal to it are uniquely defined. (b) When @a = $, (a! = 1 , . . . , y’): the law is associated, i.e. exhibits normality (in a generalized sense at corners) of the plastic strain rate tensor to the yield surface. (c) If hM = 0 (q fl = 1, . . . , y’): the instantaneous behaviour is perfectly plastic, nonhardening, i.e. corresponds to the material model of ideal plasticity (unlimited plastic strain, yield locus not modified by yielding). (d) If hog = 0 for OL# 0: noninteractive yield modes, according to Koiter’s hardening rule. (e) If h, < 0: strain softening of mode (Y,whose yield surface moves inward (locally shrinks) at yielding, so that negative second-order plastic work is performed. Other specializations and detailed discussion of the above notions can be found elsewhere.r3-15~‘9

(8) (9)

Equation (8) selects among the y modes, they ‘< y which can be activated in the incremental process (i.e. the relevant L can be >O), since for them @a = 0; equation (9) can be interpreted equivalently either as a sum with respect to CYor componentwise for each.&, the equivalence being ensured by the sign constraints I& < 0 required by the yield condition (6) for ally’ ‘activable’ modes. Therefore, the ‘direct’ plastic flow rule which governs the incremental plastic response (i.e. the h, of the y’ activable yield modes)

Piecewise-linearized increments

elas toplas tic laws in finite

Elastic-plastic structural analysis is generally a timestepping, ‘historical’ procedure along the time-history of the external actions. Each step is necessarily finite and implies to find a structural response in terms of finite increments such as Aeij and AU,. A straightforward method of dealing with step At is to assume that the static situation at to is completely known (together with the previous yielding history), to solve the nonlinear rate prob-

Appl.

Math.

Modelling,

1983,

Vol. 7, April

75

Elastic-plastic BE analysis: G. Maier and G. Novati lem at r0 and to integrate e.g.:

along the interval At, so that

Uii(to+ At) = Uii(to) + 6ij(to) At

(13)

The possible significant errors and error accumulation arising in this method have been discussed.” Remedies were proposed, e.g. resting on the notion of ‘subincremerits’.... A substantially different approach, to be adopted here, rests on the concept of ‘piecewise-linearization’ and can be formulated in two alternative ways: (i) A priori piecewise-linearization. One chooses y yield functions linear in the stresses, so that the yield conditions (6) define a polyhedron which approximates the original elastic domain. By a similar choice for the plastic potentials, each yielding mode, i.e. each face of the polyhedrical yield surface, when activated provides a plastic strain contribution with fixed direction in the space of the strain components (with axes superposed to those for stress components). These assumptions, usually combined with association (@, = $3, are rather familiar in limit analysis of perfectly plastic systems.‘5*‘9 For hardening systems parallel hypotheses may be adopted for the hardening rule, by assuming that the yield functions & and the plastic potentials 4, (if distinct from &J are linear functions also of the plastic multipliers X, (which become measures of the total yielding according to the relevant modes).12 Clearly, the above approximation of the constitutive law can be limited to the region of the stress space where the stress point may be reasonably conjectured to stay along the mechanical process under consideration. In elastoplastic analysis only the ‘potentially active’ yield planes will be considered in each step, i.e. one will ignore those which, on the basis of the elastic stress response, appear unlikely to be attained by the stress point in the step. (ii) A posteriori piecewise linearization. In elastoplastic analysis the construction of a piecewise-linear approximation to the constitutive law according to the very same criteria outlined in (i) can be conveniently carried out at each step ‘locally’, i.e. in the vicinity of the stress point at the starting time to. A single, suitably chosen yield plane at each loading step is often sufficient. Automatic procedures, easy to implement in computer codes, have been proposed.2o>22 The underlying philosophy and motivation of piecewiselinearization by either approach is that, after this approximation of the constitutive law has been done, the finite-step boundary value problem can be solved without further errors besides those due to rounding-off and discretization. Approximating constitutive laws without violation of their basic features and with computational advantages, appears reasonable if kept within the margin of uncertainty exhibited by the experimental data on nonlinear material behaviour.20 Let y’ denote the linearized yielding modes regarded as potentially active in a loading step At starting at to. The cuth of they’ yield planes will be represented by the equality in the yield condition: $0, c Q$[Uii(to) + AU,]

-R,

-hqAb

< 0

(14)

where the stress gradient @z (outward normal to the plane), the ‘yield limit’ R, (the distance at to of the plane from the origin of the stress axes, if @$ is normalized) and the hardening coefficients hap are now constants. The inequalities (14), for ol,/3 = 1, . . . , y’, define the current

76

Appl.

Math.

Modelling,

1983,

Vol.

7, April

(for the step) polyherical approximation of the elastic domain (globally y’ = y or locally y’
-R,

6

0)

(15)

the finite incremental plastic process is readily seen to be governed by the relations: (16a) Ax,>0

& Ahx, = 0

and by the finite-step

counterparts

((YJ = 1, . . . ,y’) to equations

(16b)

(1) and

(2): Aeii = Cijrs Aurs + \k$ Ax, + Ae$

(17)

where \k denotes the (constant) stress-gradient of the &h linearized plastic potential. The ‘inverse’ piecewise-linearized laws for finite increments, counterparts to the rate relations (11) and (12), read: (18)

Remarks (a) For At = 6t -+ 0, the above piecewise-linear laws for finite increments formally reduce (up to the factor At) to those in terms of rates previously examined. In fact, in the limit 4: = 0 and, hence, Go,= r&6 t, for (Y= 1, . . . , y’. In this sense the finite relations can be regarded as a generalization of the rate relations. (b) For uii(to) = 0 (i.e. when the starting situation at to is the stressless, ‘original’ state of the material), 4: = -R,, the step symbols A can be removed and the finite laws become a piecewise-linear, holonomic (path independent, reversible, nonlinear-elastic) model for a single-step description of material nonlinear behaviour.13T’6 (c) Within the step At, the piecewise-linear laws are obviously path-independent (holonomic), but the irreversible nature of plastic deformations is accounded for in passing from a step to the subsequent one: e.g. by updating R, in equation (15) on the basis of the increments Alb, in the preceding step. (d) Since the hardening moduli h,p are constant, at least in a step, the yield planes translate at yielding (unless hap = 0). The special structures of the matrix hag which correspond to kinematic, isotropic and other hardening behaviours have been established.13 (e) The rate and finite direct and inverse relationships (lo), (12), (16) and (19) are all analogous: in the sense that they allconsist of a linear relation, in 5 sign-constrained variations, and a nonlinear complementarity relation requiring that in each pair of corresponding variables at least one must vanish. This common peculiar mathematical structure is known as the ‘linear complementarity problem’. In view of its important role in mathematical programming and in other areas of applied mathematics, it has been the object of many investigations. Some of the relevant results will be used later, with reference to the finite step formulations. These formulations alone will be referred to here, because of their generality and their computational usefulness; but it is worth stressing that the rate formulations are covered as a limit case (for At + 0) in what follows.

Elastic-plastic BE analysis: G. Maier and G. Novati

Boundary element formulation Equations at the boundary A typical boundary value problem over the domain a for a finite step At in the course of an evolutive (historical) elasto-plastic analysis, consists of the equilibrium equations (bi being body forces): Aoii,i + Abi = 0 the compatibility

in R

equations

(20)

(Ui being ‘small’displacements):

Aeii = f (aUi, j + AU,, i)

in s2

(21)

and the constitutive relations which will be assumed in the piecewise-linearized form, equations (16), (17) or (1 S), (19). The kinematic boundary conditions on the constrained part r, of the boundary r and the static conditions on the free part r,, read, respectively: aUi= Aiii

on

Aoiini = A,Di

r, on

rf

(229 Wb)

where iii and pi denote given displacements and assigned tractions, ni direction cosines of the outward normal vector to r. For the sake of brevity and simplicity, the boundary r will be assumed here as ‘smooth’, i.e. such that ni are everywhere uniquely defined. Similarly, the body considered will be homogeneous over the domain a as for the elastic properties, but the extension of what follows to piecewise-constant elastic tensor would be straightforward through the familiar notion of ‘assembling subregions’.‘!2 The direct integral approach to the solution of the above problem is centred on the Somigliana identity generalized to imposed, ‘anelastic’ strains eyj, namely (see references 1,2,6,23): * Cij(E)Au,(~2) f P$Q, X) Aui(x) dr J

constrained boundary I’,; the fictitious unit forces in E and the tractions p$ on I’ pertaining to the fundamental solution (see Ref. 1). Alternatively, equation (23) can be generated through a weighted residual approach and intcgration by parts.‘pl’ The discretization of Somigliana identity (23) written for a boundary point or ‘node’ &, is usually performed by the following basic operations, with clear similarity to the finite element method:25 the boundary r is subdivided into, say L, elements r,; on each r,, suitable points or ‘nodes’ are chosen in view of modelling fields over it and its geometry; over each r, the vector fields G(x) and p(x) are ‘modelled’ so that they depend on the nodal values collected in vectors u, and pe through an interpolation matrix N,(x) of assigned shape functions; the domain a (or only part of it) is subdivided into, say Lc, cells !&; like for boundary elements, on each fi2, points called henceforth ‘cell nodes’ or ‘strain points’ are chosen and the fields b(x) and e(x) are modelled through assigned functions in N,(x) interpolating the nodal values contained in b,, eC. To fix ideas, one may refer to a two-dimensional structure in plane stress, with L straight boundary elements and Lc cells, all modelled by linear interpolation functions; at a vertex between elements the outward normal ni is (uniquely) defined with respect to the original boundary l?. The discretized version of the integral equation (23) with the above symbols and in matrix notation+ reads:

L

-dSP*(tn,X)~e(x) I

dr be

1e _

l’e

r

+f

u*G,,xPc(x) lC [s

=

da

I

bc

ac

(23)

In this equation t and x are the vectors of the coordinates of ‘load point’ and ‘field point’, respectively; all integrations are meant with respect to x. If [ is internal to the domain s2, cij = 1 for i = j, cij = 0 for i #j; if t belongs to the (‘smooth’) boundary, cij = i for i =i, cij = 0 for i #j. The starred symbols concern known functions obtained from a ‘fundamental’ solution for a fictitious homogeneous linear-elastic unlimited domain a* encompassing !& subjected to unit forces Fi concentrated in [: pz represents the jth component of traction at x on the surface element dI’ of l?; U$ the ith displacement component in x; u& the component kZof the stress tensor in x. Kelvin fundamental solutions (a* = elastic space or plane) are referred to here;1*2 however, the applicability and usefulness of Mindlin’s (half-space) or Melan’s (half-plane) fundamental solutions have been demonstrated.‘Oy’l Equation (23) for internal load-points can be established or interpreted as an application of Betti’s reciprocity theorem of linear elasticity,% to the body fi acted upon by the two following systems of external actions: the actual ones (including Ae$) and the unknown reactions &i on the

The integrations of known functions in square brackets have to be carried out mostly numerically, taking into account the singularities of the fundamental solutions at x = ~n*rJ,9Jo After writing equation (24) for all boundary nodes n, let us gather in vectors Apr and Aur all independent nodal tractions and displacements, respectively, along the boundary F (performing suitable identifications in coinciding nodes), and in vectors Ab and Aa all independent nodal body forces and anelastic strains. The resulting assembled matrix equation can be written as follows: HAur=GApr+BAb+D(Aet+AeP)

(25)

Strain increments in the plastic zones If the anelastic strains on the r.h.s. of (25) were known (e.g. as in thermoelastic analysis), equation (25) would be t Bold face posed tilde ponents are corresponding summation

symbols denote matrices (and column-vectors); a superindicates tranpose; 0 is a matrix or vector whose comall zero; equalities and inequalities apply to each pair of entries. For indices concerning matrix symbols the rule does not apply

Appl.

Math.

Modelling,

1983,

Vol. 7, April

77

Elastic-plastic BE analysis: G. Maier and G. Novati sufficient for solving the problem in the sense of the boundary element method, i.e. for computing the subvector Au” (contained in Au,) of unknown displacements on the free part of r and the subvector Ap” (in Apr) of unknown tractions on the constrained part of r. Thereafter, other quantities of possible interest within L? could be calculated making direct or indirect use of equation (23) for internal points.1*2 The presence in equation (25) of the unknown vector A@ requires to generate other equations coupled with (25). Precisely, the independent? strain increments Ae in the same cell nodes (or strain points) C;, already concerned by AeP will be computed as functions of Aur, Ap and AeP. To this purpose the following methods have been considered in the recent literature: (a) The compatibility equation (21) is applied in the finite difference sense, after having calculated displacements in suitably chosen points by equation (23). (b) As for finite elements, the displacement field is modelled over the cells considered; the strain field model, which follows through compatibility (21), is then used to compute the desired strains from nodal displacements calculated via (23). (c) Equation (23) written for internal points 5:is considered in its original form (prior to discretization) and its derivatives with respect to Ei according to (21) are constructed analytically. Because of the difficulties related to the singularity of c;“,, for x = [, the correct expressions for the main fundamental solutions were achieved only recently.‘-’ This method, though laborious, turns out to be preferable for computational efficiency and accuracy. Particular provisions are required for calculating strains (or stresses) in points on the boundary;’ however, this requirement is not so severe if all strain points are internal to the cells. The discrete version of the integral expression for Aeij(&) can be constructed by the same line of thoughts followed for generating equation (24). Semianalytical integration techniques were proposed and shown fruitful for the integral concerning the &derivatives of c&i1 After integration and assemblage, the desired equation can be cast into the form: A~=~~Au,-GAP~+~A~+(~+C)(A&+A~~) (26) where C denotes the diagonal matrix which contains the ‘free terms’ related to the singular integrals of the derivatives of ui*,,. The use for the capped matrices in (26) of the same symbols as in (25) is merely intended to emphasize the correspondence of roles. The vector of stresses in the strain points can be derived through the block-diagonal matrix k which contains as repeated diagonal block the matrix of elastic moduli: Au = k(Ae - Aer - AeP)

(27)

Substituting (26) into (27) one obtains, with obvious meaning of the new symbols: Au = H’Aq

- G’ Apr + B’Ab + (D’+ C’)(A&+

The final ingredient of the formulation by the plastic constitutive relations.

is provided, clearly,

t This implies 6 (3 in the plane) ‘engineering strain’ components 2C,,2e,,2e,,)=z in each point: {E11E22E33

78

A@) /? 0,

Appl. Math. Modelling, 1983, Vol. 7, April

The plastic strain increments in all strain points will be expressed as linear function of the plastic multiplier increments, according to equation (17). In matrix notation: AeP=\kAL

(29)

where Ah is the vector of all Aka (cu = 1, . . . , y’) in all strain points considered in a given order; matrix \k repressents a block-diagonal matrix whose blocks in diagonal positions correspond to strain points and contain as columns the plastic potential gradients ‘k$ (expressed as vectors in the space of the independent strain components) of the relevant strain point. Analogously we will form: the vectors $ and #’ of all yield functions involved in the step calculation; the blockdiagonal matrix @ of all yield function gradients; the matrix h E diag[h,p] of all hardening coefficients. By means of these new symbols, the ‘direct’ relationship (16) for all strain points can be formulated compactly as follows: @=&Ao-hAx+$‘O,

(3Ba)

SAA=0

(30b)

The parallel compact representation (19) reads:

of all ‘inverse’ relations

(31a) (31b) Thus, the direct boundary element method for the finitestep elasto-plastic analysis has been led to two alternative well-posed formulations. The former (I) associates: the usual discrete version (25) of the boundary integral equation; the expression (28) for the stress increments in all strain points; the flow rule (29); the direct laws (30) which govern the dependence of plastic strain increments on stress increments. The latter (II) combines with equation (25): the expression (26) for the strain increments; the flow rule (29); the inverse laws (31) which define the dependence of plastic strain increments on (elasto-plastic) strain increments Aeep.

Step analysis as a complementarity problem: some theoretical results Formulation of the complementarity problem A typical move in the direct boundary element method is to rearrange equation (25J so that the variable unknown quantities {A?, Ap”} = AX on the boundary be separate from the boundary data (superscript d). Thus (25) becomes: A AX = Afr + B Ab + D(Ae’+ having set after partitioning 4-[;1::],

AeP)

(32)

of H and G: _

Afr= tz5]

(3

(33)

Through same easy manipulations and defining -4’ and Aft! by means of submatrices of H’ and G’ precisely like in (33) by those of H and G, equation (28) becomes: Aa = A’AX + Af; + B' Ab + (C’+D’)(Ae’+

AeP) (34)

By inverting the nonsingular, nonsymmetric matrix A, (32) can be solved with respect to AX and sub-

equation

Elastic-plastic BE analysis: G. Maier and G. Novati stituted in (34). This leads one to express the stress increments (in the same strain points considered in forming A&‘) in terms of external action and plastic strain increments alone :

A general extremum property Consider the optimization problem: minimize : x(Ah) = A”h(h - &If*)

Ao=A#'+Ad'

(35)

where :

(36)

Au+A'A-'Afr + Af;

(37)

A& = (A’A-‘B + B’) Ab

(38)

A$ =MAe’

(39)

M f A’A-‘D + (C’ + D’)

(40)

The mechanical inter retation of the new symbols is straightforward. A oEpcontains the stress increments in the fictitious linear elastic response of the structure to external actions (fictitious, as there was no plastic yielding in the step). This response is expressed as superposition of the effects of the loads and displacements assigned on the boundary (37) of the body forces (38) and of the thermal strains assigned on the domain (39). The square matrixM (40) operates the transformation of strains imposed in strain points (and, hence, through the cell interpolation matrices, over the domain a) into the consequent stresses there. Clearly, this transformation rests not on the constitutive laws in the strain points, but on the elastic properties, geometry and constraints of the whole structure. Equations (29) and (30) instead, ignore the overall system and reflect only the constitutive laws separately for each strain point and relates the plastic strain increments there to the stress increments Au. Substituting (29) into (35) and (35) into (30) one obtains:

I$=~A~+@~-(~-&M'P)AKO

(414

AX>0

(41b)

iAh=

In mathematical terms, this is a ‘linear complementarity problem’ in the variable vectors AX and 9. Mechanically, this relation set governs the distribution of the plastic yielding (vector Ax) throughout the system, when its fictitious elastic response to external action increments is known. The formal analogy with the relation sets (10) and (12) concerning the plastic constitutive laws is worth noting. As expected, an analogous mathematical structure can be arrived at in the finite element method (and in other approaches to discrete elasto-plastic analysis); in that context it was previously pointed out in Refs 13, 16, and 18. The difference of meaning of the variables (true stresses and strains in strain points rather than ‘generalized’ variables) is not substantial. However, the generation of vector AoE and matrix M is substantially different and the procedure used here is bound to benefit from the peculiar features of the boundary element method, primarily from the smaller size of the matrix to invert A. In view of equation (27), it can be easily realized that the very same relations set (41) would be constructed, if the same path of reasoning were applied to the latter (II) formulation at the end of the preceding section, i.e. to equations (2.5) (29) and (26) (31) instead of (28) (30).

(42a)

subject to: (h-&M’P)Ax>&Ati+#’

AhA&+A~+A&

Ax - A%b A8

AL>0

(42b)

If problem (41) has a solution, this solution also solves problem (42); if the minimum of problem (42) is zero, then any solution of this problem also solves the complementarity problem (41); otherwise, problem (42) has no solution. The equivalence, in the above sense, between the two problems can be demonstrated straightforwardly. The extremum formulation (42) is a problem in (generally nonconvex) quadratic programming and hardly has computational value. However, when matrix (h -&M*) is symmetric positive definite, two much stronger dual extremum characterizations of the solution can be established (by means of Kuhn-Tucker theorem and duality theory of constrained optimization).26 These are convex quadratic programs, each with only one of the two constraint sets (42b). These developments, parallel to those established by Maier in a finite element context,13j18 will be pursued elsewhere. Existence and uniqueness of solution The following statements can be directly derived from known results in complementarity theory (see, Refs. 27, 28 and references given there): (a) a unrque solution to problem (41) exists for any vector of data @ At? + Go, if and only if the matrix (h - @MQ) has all principal minors positive; (b) the number of solutions to (41) is finite for all vectors of data, if and only if all the principal minors of the above marix are nonzero, (c) if matrix (/r -@Mq) is positive semidefinite and the linear constraints (42b) are consistent, then the problem is solvable; if it is positive definite, the constraints (42b) are always consistent (i.e. they define in the Ai space a nonempty domain). It is worth noting that in statement (a) ‘for any vector of data’ implies ‘for any set of external action increments’, but the converse is not necessarily true. Uniqueness of the Ax solution to problem (41) means also uniqueness of the incremental response in terms of stresses, strains and displacements, since all these quantities are linear elastic consequences of Ah and the given external action increments.

Some solution procedures The procedures outlined below for solving the finite-step problem in elasto-plastic analysis with piecewise-linearized incremental constitutive laws, flow more or less directly from the linear complementarity formulation (41). The first two methods (so called ‘initial stress’ and ‘initial strain’) appear already in use, with marginal variants, in the boundary element literature,6,10,” but acquire here a unitary framework seemingly lacking so far. Other more elaborate techniques, which are out of the scope of this paper, might be transferred from the domain of mathematical programming, once the linear complementarity

Appt. Math. Modelling, 1983, Vol. 7, April

79

Elastic-plastic BE analysis: G. Maier and G. Novati

nature of the mechanical problem in point is recognized. Convergence conditions, comparative assessment of computational merits and illustrative examples cannot be dealt with here.

(2) Solve (e.g. by familiar Gauss elimination or some variant of it) the linear system consequent to the above changes: (h-&Mq)Ah=&A#

(46)

‘Initial strain ’ method

(1) Choose an initial vector Ah or A&’ (e.g. = 0, or, for proportional loading with equal steps, = the one obtained in the preceding step); calculate the selfstresses Aa’ = MAeQ = M\kAx. (2) Setting 4 = 0 in the direct plastic law (30a) or in (41a) compute Ax by solving: h AA = i(AoE

+ Au”) + $J’

(43)

(3) For each component AX, of Ai choose the maximum Axi between zero and its value obtained from (43) thus satisfying the constitutive sign constraint. (4) Perform a convergence test on AA’ or on Aep = *Ah’; if not satisfied go to (1) and replace Ax by the new vector Ax’, or AeP by Aep . Phases (2) and (3) can be replaced by a single statement: solve in each strain point the linear complementarity problem of the ‘direct’ law (30a). ‘Initial stress’ method (1) Choose an initial Ah or a AeP; calculate Ad = M AeP = M\k Ax and Aeep = AeP + k-’ (AoE + A#). (2) Setting $ = 0 in the inverse plastic low (31a), compute

Ah by solving (&k\k + h) AA = @k A&P + @JO

(44)

Note that this expression can be obtained from (31a) for $ = 0 if @k A$’ is added and subtracted in (3 la). (3) Evaluate Ahi = max(0, Ahi from (44)) and, if needed, Aen’ = \k Ax’. (4) Perform a convergence test; if not satisfied go to (1) and replace AA or Aep by the new (primed) vector. Phases (2) and (3) can be replaced by: solve in each strainpoint the linear complementarity problem of the ‘inverse’ law (31a). Over-relaxation

(1) Choose an initial vector Ax (e.g. Ax = 0). (2) Set $O = 0 in equation (41a); then, with reference to the linear equation: (h-~MQ)Ax=&AaaE+$ro

(45)

Calculate Ax?) at the nth iteration step according to the familiar Gauss-Seidel recurrent formula for linear equations, then compute: Ai?)

= max{O, A@-‘)

+ w(Av)-

A@-l))}

(46)

where o is the over-relaxation factor to be chosen in the interval 0 + 2. (3) After a convergence test go to (2) or stop. Tangent modulus (1) Inspect vector #JOQ 0 (yield functions

at the beginning of the step interval): set go= 0 the components whose absolute values are below an assigned tolerance; consider as not activatable in the step the yield modes corresponding to the other components i.e. ignore them, in forming they’ yield conditions in each strain point so that 4’= 0.

80

Appl.

Math.

Modelling,

1983,

Vol.

7, April

(3) After a convergence test go to (1) considering not activatable also the yield modes i corresponding to negative AX, in the solution of (46). Quadratic programming

Any traditional algorithm for solving convex quadratic programming problems (such as those of Beale, Wolfe, etc.,26 as well as most nonlinear programming algorithms, are known to generate at least a Kuhn-Tucker point, if not the (global optimum) solution. Hence, any commercial code inplementing such algorithm when applied to the optimization problem (42) will generate: the solution of (42) if matrix (h - @Mk) is positive definite; the/a solution, if it is positive semidefinite; the/a solution whatever this matrix may be, provided the minimum be zero. The equivalence between problems (42) and (41) in the sense of the preceding section permits one to derive the/a solution to the mechanical problem in point. For a quadratic programming in finite element elastoplastic analysis (see references 12 and 30). Remarks

(a) Both matrices h in (43) and (&k\k + h) in (44) are easy to invert, because of their block-diagonal nature. But the latter can be inverted also for h = 0, i.e. for perfectlyplastic behaviour, while the former cannot. This fact is in favour of the ‘initial stress’ method, which uses the constitution to obtain stresses from strains instead of in the converse direction (see references 15, 29-3 1 for the FE context). (b) The over-relaxation method represents Cryer’s generalization to linear complementarity of the accelerated version of the Gauss-Seidel iteration scheme. Convergence is guaranteed with 0 < w < 2 for a positive definite matrix. The optimal over-relaxation factor w proves to be almost constant for a certain class of problems; therefore this factor can be estimated for that class by numerical tests. However, the lack of sparsity in the matrix is likely to make the ‘direct’ (elimination of pivoting) methods preferable (see Refs 26, 27 and relevant references). (c) The above procedure called tangent modulus bears a close resemblance to that proposed years ago in the finite element context.32 Under proportional loading the loop 3 + 1 might be given up, by accepting the possible error of some negative Axi. The mechanical interpretation of (46) and the name rest on the remark that in each strain point an elastic-plastic compliance is attributed to the material, but otherwise this is dealt with as linear. In fact setting r#~= 0 in the direct plastic law (16) (17) one obtains: Ae = (k-’ + +h-‘a)

Au

(47)

where the block-diagonal matrix in parenthesis contains in diagonal position the ‘tangent’ local compliance matrices. (d) The linear complementarity formulation (41) in the variables Ai and 4 alone (having substituted all other unknowns), have been established and used in ‘he finite element and other contexts. But it is worth hi ssing that the

Elastic-plastic BE analysis: G. Maier and G. Novati

boundary element method is particularly suitable and advantageous in this respect. In fact the order of the inverse matrix A-’ needed in constructing the ‘strain to selfstress’ transformation matrix, is equal to the number of the unknowns on the boundary alone; in the finite element approach the counterpart would be the system elasticstiffness matrix, of the order equal to the number of degrees of freedom. The direct use of the ‘strain-to-selfstress’ influence functions for the continuum in the boundary integral equation approach has been proposed by Polizzotto.33 (e) Clearly, the size of the complementarity problem is given by the sum of the numbers y ’ of yield conditions (or activatable modes in the load step) in all strain points. Actually, one will set y’ # 0 (mostly y’ = 1 or a small number) in all strain points within the ‘plastic zones’, i.e. those parts of the domain fi where plastic yielding in the current step is reasonably expected. This conjecture and selection can be automatically performed on the basis of the linear elastic stress response A#, and corrected if Au resulting at the end of the calculations for the step, leads to the violation of some of the neglected yield conditions. It is worth noting that, the dimensionalities of the vectors A& (plastic strains), Ae’ (thermal strains) and Ab (body forces) are quite different: the last two vectors, generating data only, might concern cells over the whole domain a; hence, they may be much larger than A# which, on the contrary, will concern plastic zones alone. For the sake of simplicity of notation this computationally important circumstance has not been explicitly expressed in the foregoing developments .

give rise to some alternative

Acknowledgement Support from the CNR, PAdIS, is gratefully acknowledged.

References

4

5

6

7

8 9 10

Conclusions The findings of this paper can be summarized

as follows:

(a) In view of the numerical step-by-step analysis of elastoplastic systems, an approximation of the constitutive laws in terms of finite increments has been proposed with the following feature: mathematically, it is a linear complementarity problem with the dimensions equal to the number of ‘potentially active’ current yield plans; mechanically, it covers lack of normality, corners and hardening (or softening) and formally reduces to the traditional rate relations as the step amplitude tends to zero; computationally, it can help reducing error accumulation due to the finite nature of the load step. (b) The direct boundary element method has been used in association with the above approximate description of nonlinear behaviour, and the elastoplastic analysis by finite (as well as infinitesimal) steps has been cast into the form of a linear complementarity problem. This involves the plastic multiplier increments in all strain points over the potentially yielding zones of the structure. The fitness of the boundary element method to generate such compact formulation has been demonstrated. (c) As straightforward consequences of the complementarity formulation pointed out, a very general extremum characterization and conditions for solution uniqueness have been established; stronger statements for less general classes of material behaviours appear possible (and will bC considered separately). (d) In the unitary framework developed for the problem in point, various procedures for the numerical solution have been envisaged in a version matching the new complementarity formulation. This formulation is found to

solution techniques, the of which may be worth exploring elsewhere.

potentialities

11

12

13

14

15

16 17

18

19 20 21

22 23

Banerjee, P. K. and Butterfield, R. ‘Boundary elements methods in engineering science’, McGraw-Hill, London, 1981 Brebbia, C. A. and Walker, S. ‘The boundary element techniques in engineering’, Newnes-Butterworths, London, 1980 Mendelson, A. and Albers, L. U. ‘Application of boundary integral equations to elastoplastic problems’, in ‘Boundary Integral Equation Method: Computational Applications in Applied Mechanics’ (Cruse and Rizzo, Eds), ASME, New York, 1975,47-84 Swedlow, J. L. and Cruse, T. A. ‘Formulation of boundary integral equations for three dimensional elastoplastic flow’, Inr. J. Solids Struct. 1971,7,1673 Banerjee, P. K. and Mustoe, G. G. ‘The boundary element method for two dimensional problems of elastoplasticity’, in ‘Recent Advances in Boundary Element Methods’ (Brebbia, ed.), University of Southampton, 1978,283-300 Banerjee, P. K: and Cathie, D. N. ‘A direct formulation and numerical implementation of the boundary element method for two dimensional problems of elastoplasticity’, Int .I. Mech. Sci. 1980, 22,233 Bui, H. J. L. ‘Some remarks about the formulation of three dimensional thermoelastoplastic problems by integral equations’, Int. J. Solids Struct. 1978, 14, 935 Mukherjee, S. ‘Corrected boundary integral equations in planar thermoplasticity’, Int. J. Solids S&t. i977, i3, 331 Telles. J. C. F. and Brebbia. C. A. ‘On the aunlication of the boun&ry element method to plasticity’, Ai&. Math. Modelling 1979,3,466 Telles, J. C. F. and Brebbia, C. A. ‘Elastoplastic boundary element analysis’, Proc. Europe-U.S. workshop on nonlinear finite element analysis in structural mechanics (Wunderlich et al., Eds), Ruhr-University Bochum, Springer Verlag, 1980, 403-434 Telles, J. C. F. and Brebbia, C. A. ‘New developments in elastoplastic analysis’, Proc. 3rd Znt. Seminar on Recent Advances in Boundary Element Methods (Brebbia, ed.), Irvine, 1981,350-370 De Donato, 0. and Maier, G. ‘Finite element elastoplastic analysis by quadratic programming: the multistage-method’, 2nd Int. Conj: on Struct. Mech. Reactor in Technology, SMiRT, Berlin, Vol. V, Part M, 1973 Maier, G. ‘A matrix structural theory of piecewise-linear plasticity with interacting yield planes’, Meccanica 1970, 5 (2), 1 Maier, G. ‘Sui legami associati tra sforzi e deformazioni incrementali in elastoplasticiti’, Rend. Ist. Lomb. Scienze e Lettere 1966,100,809 Maier, G. ‘Mathematical programming methods in structural analysis’, Variational Methods in Engineering, Proc. Int. Symp. (Eds H. Tottenham and C. Brebbia), Southampton Univ. Press, 1973,1-32 Maier, G. ‘Quadratic programming and theory of elasticperfectly plastic structures’, Meccanica 1968, 3 (4), 1968, 265 Maier, G. ‘A minimum principle for incremental elastoplasticity with non associated flow laws’, J. Mech. Phys. Solids 1970, 18, 319 Maier, G. ‘Incremental plastic analysis in the presence of large displacements and physical instabilizing effects’, Int. J. Sol. Struct. 197 1,7,345 Martin, J. B. ‘Plasticity’, MIT Press, 1975. Hodge, P. G., Jr. ‘Automatic piecewise linearization in ideal plasticity’, Comp. Meth. Appl. Mech. Eng. 1977, 10, 249 Schreyer, H. L., Kulak, R. F. and Kramer, J. M. ‘Accurate numerical solutions for elastic-plastic models’, Trans. ASME, J. Pressure Vessel Technol. 1979, 101, 226 Cannarozzi, A. A. ‘A nontraditional linearizing procedure in limit analysis’, J. Struct. Mech. 1980, 8 (4). 449 Somigliana, C. ‘Sopra l’equilibrio di dn c&o elastic0 isotropo’, fl nuovo cimento 1886,17-19

Appl.

Math.

Modelling,

1983,

Vol. 7, April

81

Elastic-plastic BE analysis: G. Maier and G. Novati 24 25 26 21 28

Betti, E. ‘Teoria dell’elasticita’,n nno~o cimento, 1872,7-10 Zienkiewicz, 0. C. ‘The finite element method’ (3rd ed.), McGraw-Hill, London, 1977 Kunzi, H. P. and Krelle, W. ‘Non-linear programming’, Blaisdell, 1966 Cottle, R. W. ‘On a problem in linear inequalities’, J. London Math. Sot. 1968,43,378 Murthy, K. G. ‘On the number of solutions to the complementarity quadratic programming problem’, Linear Algebra

30

and its Applic.

33

1975

29

Argyris, G. H. and Sharpf, D. W. ‘Methods of elasto-plastic analysis’, Froc. Symp. ‘Finite Element Techniques’, Stuttgart, I.S.D., 1969

82

Appl.

Math.

Modelling,

1983,

Vol.

7, April

31

32

Franchi, A. and Genna, F. ‘Plastic collapse in plane problems: results of a numerical approach’, Proc. VI’ Congress0 Naz. AIMETA, Genova, 7-9 October 1982 Zienkiewicz, 0. C., Valliappan, S. and King, I. P. ‘Elastoplastic solutions of engineering problems, initial stress finite element approach’, ht. .I. Num. Meth. Eng. 1969, 1,1969,75 Marcal, P. V. and King, I. P. ‘Elastic-plastic analysis of two dimensional stress systems by the finite element method’, Int. J. Mech. Sci. 1967,9,143

Polizzotto, C. ‘An alternative formulation of the boundary element method’, Appl. Math. ModelZing 1982,6 (21, 73