Post-collapse analysis of plates and shells based on a rigid–plastic version of the TRIC element

Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775 www.elsevier.com/locate/cma

Post-collapse analysis of plates and shells based on a rigid–plastic version of the TRIC element q Leone Corradi a, Nicola Panzeri b

b,*

a Department of Nuclear Engineering, Politecnico, Via Ponzio 34/3, Milano 20133, Italy Department of Structural Engineering, Politecnico, Piazza Leonardo da Vinci 32, Milano 20133, Italy

Received 15 December 2002; received in revised form 22 May 2003; accepted 22 May 2003

Abstract The knowledge of the post-collapse response of structures is required in several situations. Typical examples are shells used as energy absorbers or bumpers, which must be able to undergo large plastic deformations by dissipating a sufficient amount of energy. The procedure known as sequential limit analysis can be employed to this purpose. The evolution of the structure after collapse is tracked by a sequence of rigid–plastic limit analyses, with the structural geometry progressively updated on the basis of the collapse mechanism computed in the previous step. A number of methods have been developed, mostly in recent years, for the finite element solution of the limit analysis problem. In this paper, a procedure proposed by one of the authors is employed, which a rather extensive computational experience proved reliable, stable and accurate. The procedure demands that the finite element be formulated on the basis of the natural approach, which suggests the use of a simple but well performing triangular shell element, named TRIC, recently developed within this framework by Argyris and co-workers. The formulation has to be modified to some extent to adapt to the rigid–plastic context, but the quality of the element performances is maintained. Some computations, referring to plates and shells, are presented, which assess that the ‘‘marriage’’ between the TRIC element and sequential limit analysis is successful.
2003 Elsevier B.V. All rights reserved. Keywords: Plates; Shells; Post-collapse behavior; Rigid–plastic analysis; Sequential limit analysis

1. Introduction The numerical solution of the limit analysis problem for perfectly plastic structures discretized via finite elements has experienced a growing interest in recent years, as witnessed by several papers on the subject

q Partial results were presented at the 6th International Conference on Computational Structures Technology (CST2002), Prague, 4–6 September 2002. * Corresponding author. Tel.: +39-02-2399-4351; fax: +39-02-2399-4220. E-mail address: [email protected] (N. Panzeri).

0045-7825/$ - see front matter
2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00373-6

3748

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

(see, e.g., [1–10]). Solution methods mainly exploit the kinematic (upper bound) theorem, reducing the problem to the search of the minimum of a convex, albeit non-smooth function. These methods can be employed also to obtain indications on the structural response subsequent to collapse, which must be considered in several instances. Examples are structures that must be able to develop large plastic deformations, such as energy absorbers or bumpers. The post-collapse evolution is followed by performing a sequence of rigid–plastic limit analyses with the structural geometry updated on the basis of the collapse mechanism obtained at the end of the previous step. Since elastic strains are neglected, only the rigid–plastic response is tracked, which however provides a meaningful piece of information on the energy that the structure can dissipate and on its deformation capabilities. Such a procedure is known as sequential limit analysis. To the authors knowledge, it was first proposed in the early 1960s to study the post-collapse response of plastic frames (see [11]). The method received new attention with the work of Yang on trusses [12] and recently Seitzberger and Rammerstorfer [13] used it successfully to simulate the large deformation crushing behavior of shells with polar symmetry. In a previous work [14] the authors performed a fairly extensive study by computing the post-collapse rigid–plastic curves for cylindrical shells, cones and hemispheres under different load conditions. Results assess the effectiveness of sequential limit analysis, which appears accurate, robust and reliable both for rising (stable) and decreasing (unstable) behaviors. In the paper above, sequential limit analysis was employed in conjunction with the limit analysis procedure established in [6]. In contrast to alternative approaches, mostly based on regularization techniques originated by a proposal of Huh and Yang [1], this procedure deals with the non-smooth nature of the function to be minimized by detecting and eliminating from the problem the finite elements that do not undergo plastic flow in the collapse mechanism, which are considered as rigid. To enforce this condition, it is essential that the finite element be formulated on the basis of the natural approach introduced by Argyris [15,16], separating deformation modes from rigid body motions: in fact, in a rigid element the first vanish, but the latter must survive. Only shell elements exhibiting polar symmetry were considered in [14] and, in this situation, the natural formulation does not entail difficulties. However, polar symmetry is a quite severe restriction, in that only particular shells can be handled and even axisymmetric shells might collapse according to diamond-shaped mechanisms or loose polar symmetry during their post-collapse evolution. To broaden the range of applicability of the method, more general elements are demanded. Well performing shell elements, however, are rare and usually involve high order shape functions, which makes the separation of deformation from rigid modes cumbersome. Recently, a simple triangular shell element, named TRIC, has been developed by Argyris and co-workers [17], and its formulation rests on the natural approach, making it a spontaneous candidate for the procedure. The element is simple but sophisticated (using authorsÕ terminology), non-conforming but converging [18], accounts for transverse shears without locking at the thin shell limit and several elastic computations assess its robustness and accuracy also in situations involving large displacements and rotations. Very recently the element was used for elastic–plastic analyses and results indicate that its performances are maintained also in the presence of material non-linearity [19]. In this study, a rigid–plastic version of the TRIC element is proposed, which appears suited for sequential limit analysis. To preserve the element performances, modifications with respect to the original formulation are reduced as much as possible, but the different context demands some change. In particular, the role of the elastic strain energy, a positive definite quadratic form, is played by the power of dissipation, a positively homogeneous of degree one function of strain rates, and this affects the formulation to some extent. On the other hand, other aspects are simpler, since limit analysis does not require the definition of the geometric stiffness and, when based on the kinematic theorem, of the stress state within the element. Moreover, the TRIC element has the capability of dealing with layered shells, which is not exploited in this paper, where only homogeneous shells are considered.

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3749

The outline of the paper is as follows. Some basic concepts of classical plasticity are summarized in Section 2. Section 3 describes the kinematics of the element, which remains as originally proposed in [17,18], except that deformations are purely plastic. The bulk of the paper consists of Section 4, computing the power of dissipation for the rigid–plastic TRIC element, Section 5, formulating the discrete limit analysis problem and outlining the solution strategy, and Section 6, where some computational aspects, specific of this study, are discussed. Finally, in Section 7 some examples are presented, which seem sufficient to establish that the ‘‘marriage’’ between TRIC and sequential limit analysis is a successful one.

2. Summary of basic results Limit analysis deals with rigid–perfectly plastic materials. Such a constitutive model assumes that stresses are confined within the convex domain uðrÞ ¼ f ðrÞ
r0 6 0;

ð2:1Þ

where u is the yield function (regular for simplicity) and r0 a yield limit. Deformations cannot occur as long as u < 0, while plastic flow may develop when equality holds. In this case, strain rates obey the normality rule ou k; k P 0: ð2:2Þ or In the equations above, r denotes the stress tensor, while e (as other geometric quantities) refers to rates. Eq. (2.2) might impose some restrictions on strain rates, by confining them within a (convex) domain De , the subspace spanned by the outward normals to the yield surface (typically, for von MisesÕ condition plastic flow is isochoric and De is the deviatoric subspace). The inclusion e 2 De defines the plastically admissible set for strain rates. b ðeÞ, known as power of dissipation (per unit volume) It is possible to associate to any e 2 De a function D and defined by HillÕs maximum principle as e¼

b ðeÞ ¼ rðeÞ : e ¼ max r : e subject to uðrÞ 6 0: D r

ð2:3Þ

The value rðeÞ of r solving the problem above is a stress state associated to e through the normality law, uniquely defined if the domain uðrÞ 6 0 is strictly convex. When the domain is convex but not strictly so, rðeÞ may not be unique, but the power of dissipation still is a uniquely defined function of strain rates. Such function is convex and positively homogeneous of degree one. Hence b b b ðeÞ ¼ rðeÞ : e ¼ o D : e ) rðeÞ ¼ o D : D oe oe

ð2:4a;bÞ

Eq. (2.3) implies DruckerÕs inequality ðrðeÞ
rÞ : e P 0 for all r such that uðrÞ 6 0. When imposed for all plastically admissible strain rates, the inequality enforces compliance with the admissible domain for stresses, i.e. ! b oD ðrðeÞ
rÞ : e ¼ ð2:5Þ
r : e P 0 8e 2 De ) uðrÞ 6 0: oe b ðeÞ completely describes the behavior of a rigid–perfectly plastic material: if its expression is Therefore, D known, the yield function needs not to be defined explicitly. The limit analysis problem considers a rigid–perfectly plastic solid subject to body forces kF on its volume X and surface tractions kf on the free portion oF X of its boundary. The constrained boundary oU X

3750

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

is fixed. Loads are defined as basic values F and f, affected by a load multiplier k, and the value s of k for which collapse is attained (collapse multiplier) is sought. The kinematic theorem of limit analysis states that s is the optimal value of the minimum problem Z b ðeÞ dX D ð2:6aÞ s ¼ min e;u

X

subject to e ¼ rS u in X;

u¼0

on oU X;

e 2 De

ð2:6b;cÞ ð2:6dÞ

in X; Z Z PðuÞ ¼ F  u dX þ

f  u dðoXÞ ¼ 1:

ð2:6eÞ

oF X

X

Eq. (2.6b,c) express compatibility, associating to a velocity field u, vanishing on the constrained boundary, the consequent strain rate distribution (rS u is the symmetric part of the velocity gradient). Eq. (2.6d) establishes the plastically admissible nature of strain rates and together with the compatibility conditions defines a mechanism. PðuÞ denotes the power of basic loads, which Eq. (2.6e) normalizes to unity. The minimum problem (2.6) is a classical and well-known result. Nevertheless, a proof of the statement is useful for future developments and is now outlined. Let Eqs. (2.6b) and (2.6e) be enforced in weak form by writing Z
 b ðeÞ
r : ðe
rS ðuÞÞ dX
jð1
PðuÞÞ ; Lðe; u; r; jÞ ¼ D ð2:7aÞ X

where rðxÞ and are Lagrangean multipliers. Then, the optimality condition reads [20] ð2:7bÞ

dL P 0 subject to e 2 De

in X;

u¼0

ð2:7c;dÞ

on oU X:

b ðeÞ is not differThe inequality reflects the non-stationary nature of the problem, due to the fact that D entiable when e ¼ 0. However, L is differentiable with respect to all other variables and Eq. (2.7b) must hold as an equality when variations with respect to them are taken. When referred to the Lagrangean multipliers, the condition enforces the constraints (2.6b), (2.6e). By imposing that L be stationary with respect to u one obtains Z  Z Z du L ¼ r : dðrS uÞ dX
j F  du dX þ f  dudðoXÞ ¼ 0 8du ¼ 0 on oU X: ð2:8Þ X

X

oF X

On the other hand, when e is varied dL is not necessarily zero and the condition reads Z
 b
r : de dX P 0 8de such that e þ de 2 De de L ¼ dD

ð2:9Þ

X

It is now necessary to distinguish between the regions Xp , where plastic flow actually develops, and Xr , b is differentiable and (2.9) holds as an equality. Since d D b ¼ ðo D b =oeÞ : de and which keeps rigid. In the first D because of Eq. (2.4b), one obtains r ¼ rðeÞ

in Xp :

ð2:10aÞ

b ¼ 0 and, hence, d D b ¼D b ðdeÞ ¼ rðdeÞ : de, rðdeÞ being the stress associated to variations de In Xr it is e ¼ 0, D through normality. Then, Eq. (2.9) establishes   ð2:10bÞ rðdeÞ
r : de P 0 8de 2 De in Xr :

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3751

By interpreting the Lagrangean multiplier field rðxÞ as a stress tensor, it is easily recognized that Eq. (2.8) enforces equilibrium under the basic loads multiplied by j and that Eqs. (2.10) imply uðrÞ 6 0 throughout (with equality holding in Xp ). Thus, the optimal value of j is a load multiplier which is statically admissible, or safe. In addition, it coincides with the optimal value of problem (2.6): in fact, by using the principle of virtual velocities and relations holding at solution (either a priori constraints or optimality conditions) one can write (a star marks values at solution) Z Z b ðe Þ dX ¼ s: j ¼ j Pðu Þ ¼ D r : e dX ¼ ð2:11Þ X

X

3. The TRIC element 3.1. Geometrical description TRIC is a triangular shell element with three nodes and six degrees of freedom per node. Its vertices are numbered from one to three and edges are labeled as a, b, and c (Fig. 1(a)). The length of each side is denoted by ‘k , k ¼ a; b; c. A local Cartesian reference system (x; y; z) is introduced, with origin in the element centroid and z ¼ 0 coinciding with the element mid-plane. The coordinates of the nodes are indicated with ðxi ; yi Þ, i ¼ 1; 2; 3. The position of a point on the mid-plane can also be expressed in terms of triangular coordinates Li , i ¼ 1; 2; 3 (L1 þ L2 þ L3 ¼ 1, Fig. 1(b)), related to x and y through the relations 1 1 1 ða1 þ b1 x þ c1 yÞ; L2 ¼ ða2 þ b2 x þ c2 yÞ; L3 ¼ ða3 þ b3 x þ c3 yÞ; ð3:1Þ L1 ¼ 2X 2X 2X where X is the element area and a1 ¼ x2 y3
x3 y2 ; a2 ¼ x3 y1
x1 y3 ; a3 ¼ x1 y2
x2 y1 ; ð3:2aÞ b1 ¼ y2
y3 ;

b2 ¼ y3
y1 ;

b3 ¼ y1
y2 ;

c1 ¼ x3
x2 ; c2 ¼ x1
x3 ; c3 ¼ x2
x1 : The following geometrical relations are easily established and will be used in the sequel: 2X ¼ ðb1 c2
b2 c1 Þ ¼ ðb2 c3
b3 c2 Þ ¼ ðb3 c1
b1 c3 Þ;

Fig. 1. The TRIC element.

ð3:2bÞ ð3:2cÞ ð3:3aÞ

3752

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

cos a ¼

x3
x2 c1 ¼ ; ‘a ‘a

sin a ¼

y3
y2 b1 ¼
; ‘a ‘a

cos b ¼

x1
x3 c2 ¼ ; ‘b ‘b

sin b ¼

cos c ¼

y1
y3 b2 ¼
; ‘b ‘b

x2
x1 c3 ¼ ; ‘c ‘c

sin c ¼

y2
y1 b3 ¼
; ‘c ‘c

ð3:3bÞ

ð3:3cÞ

where a, b, c are the angles indicated in Fig. 1(a). 3.2. Cartesian and natural strains The Cartesian strain components eij ¼ 12 ðsi;j þ sj;i Þ are grouped into two vectors em ¼ f e x

ey

exy gt ;

es ¼ f ezx

ezy gt

ð3:4a;bÞ

collecting membrane and transverse shear contributions. As usual in shell theory, ez is not considered as an independent variable: its value follows from the constitutive law and the condition rz ¼ 0. The element properties are established by replacing Cartesian strains with local quantities referred to the natural element system, defined by the three sides of the triangle. Natural membrane strains em ¼ fea eb ec gt measure the change in length of unit segments parallel to the sides, i.e. in the directions indicated as Yk (k ¼ a; b; c) in Fig. 1(b). The rules of tensor calculus express them as functions of Cartesian strains through the relation 8 9 2 38 9 < ea = cos2 a sin2 a 2 cos a sin a < ex = eb ¼ 4 cos2 b sin2 b 2 cos b sin b 5 ey ; ð3:5Þ : ; : ; 2 2 ec e cos c sin c 2 cos c sin c xy which can be inverted to give

em ¼ Tm em ;

2

b2 b3 ‘2a 1 6 c2 c3 ‘2a 6 Tm ¼ 24 4X b2 c3 þ c2 b3 2 ‘a 2

b3 b1 ‘2b c3 c1 ‘2b b3 c 1 þ c 3 b1 2 ‘b 2

3 b1 b2 ‘2c 7 c1 c2 ‘2c 7: b1 c2 þ c1 b2 2 5 ‘c 2

ð3:6Þ t

Three components of natural transverse shears are also introduced, collected in vector es ¼ fga gb gc g . Each component is the shearing strain between segments parallel to one edge and axis z. Simple geometrical considerations establish 8 9 2 3 cos a sin a   < ga = e gb ¼ 4 cos b sin b 5 zx : ð3:7Þ e : ; zy gc cos c sin c As the relation above shows, the three components of vector es are not independent. Nevertheless, when formulating the element properties, they are considered as such: as explained in detail in [18], this plays a key role in eliminating locking phenomena at the thin shell limit. Obviously, Eq. (3.7) cannot be inverted as it is. However, it is possible to solve three systems of two equations, by considering only two rows of the matrix. For instance, from the first two of them one obtains  
1        1
‘a b2 ‘b b1 ezx ga ga cos a sin a ¼ ¼ ezy 1 gb gb cos b sin b 2X
‘a c2 ‘b c1 or es1 ¼ Ts1 es

 1
‘a b2 Ts1 ¼ 2X
‘a c2

‘b b1 ‘ b c1

 0 : 0

ð3:8aÞ

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

By proceeding in the same way with reference to the other couples of rows, one also obtains   1 ‘a b3 0
‘c b1 es2 ¼ Ts2 es Ts2 ¼ ; 2X ‘a c3 0
‘c c1 es3 ¼ Ts3 es

Ts3 ¼

 1 0 2X 0


‘b b3
‘b c3

 ‘c b2 : ‘ c c2

3753

ð3:8bÞ

ð3:8cÞ

In deriving the expressions for matrices (3.6) and (3.8), use was made of Eqs. (3.3). 3.3. Model for natural strains According to the kinematics of the shell model, one can express natural strains as em ¼ g þ zv;

es ¼ t:

ð3:9a;bÞ

Eq. (3.9a) separates axial from bending contributions in the membrane strain components: g collects the membrane deformations of the mid-surface and v the natural curvatures. Eq. (3.9b) is just a change of name for transverse shear strains. Each contribution can be modeled in terms of 15 deformation modes, belonging to five groups of three (typical modes are illustrated in Fig. 2), which are now described. t Vector qm ¼ fqma qmb qmc g collects the axial modes, each corresponding to a change in length of a single side of the triangle, and governs the mid-plane stretching g. Curvatures are produced by six bending modes, collected in vectors qSb and qAb , producing transverse displacements, either symmetric or antisymmetric, on one edge of the triangle. Transverse shears are governed by the vector qs ¼ fqsa qsb qsc gt of shear modes, each predicting shearing of one side of the triangular prism, by preserving right angles in the two remaining t sides. Finally, qd ¼ fqd1 qd2 qd3 g is the vector of so-called drilling (or azimuth) modes, involving a rotation in the plane about a single vertex. In principle, they would contribute to g, but they are supposed not to affect local fields directly: their presence will be accounted for by means of a fictitious dissipation. Let b q t ¼ f qtm

qtd

qtSb

qtAb

qts g

ð3:10Þ

be the vector collecting the modes above. Natural strains in the element are modeled as follows: q; bg b gðLi Þ ¼ qm ¼ b

b bg ¼ ½ I

0

b v ðLi Þb q; vðLi Þ ¼ a1 qSb þ a2 ðLi ÞqAb ¼ b tðLi Þ ¼ a3 ðLi Þqs ¼ b b t ðLi Þb q;

0 ;

0 0

ð3:11aÞ

b b v ðLi Þ ¼ ½ 0

b b t ðLi Þ ¼ ½ 0

0

0 0

0

a1

a3 ðLi Þ :

a2 ðLi Þ 0 ;

ð3:11bÞ ð3:11cÞ

Matrices an have been constructed in [17,18]. For the sake of completeness, their expressions are given in Appendix A. 3.4. Natural modes and relations to nodal displacements The TRIC element has 18 degrees of freedom, three displacements along the Cartesian axes and three rotations about them at each node (Fig. 3). According to the natural approach introduced by Argyris [15,16], it is possible to establish a one-to-one relation between nodal displacements U and the set of rigid body motions and natural deformation modes (or element generalized strains), collected in vectors q and q, respectively. Symbolically       q q Cq U ¼ ½ Aq A  U: ð3:12a;bÞ ; ¼ C q q

3754

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

Fig. 2. Typical deformation modes.

z

y

ψ

ϕ

w v u

x ϑ

Fig. 3. Nodal degrees of freedom.

Since six rigid body motions are present, the element has 12 natural modes. As it was discussed in [17,18], the relevant vector q reads qt ¼ f qtm

qtd

qtSb

qtA g;

ð3:13Þ

where qm , qSb and qd are as defined before and qA ¼ qAb þ qs ;

ð3:14Þ

i.e., antisymmetric bending and shear sum to produce the same natural modes. The relation between vector q and the vector (3.10) of expanded natural modes b q can be written as 2 3 I 0 0 0 0 60 I 0 0 07 6 7; q ¼ Lb q; L¼4 ð3:15Þ 0 0 I 0 05 0 0 0 I I where I and 0 denote the 3 · 3 identity and null matrices, respectively.

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3755

The second part q ¼ CU of Eq. (3.12b) establishes the compatibility relation for the finite element, by expressing its generalized strains as functions of nodal displacements (the expression for matrix C, which is as given in [17], is also reported in Appendix A). By combining this relation with Eq. (3.15), one obtains Lb q ¼ CU:

ð3:16Þ

3.5. Transformation to global coordinates and assemblage Eq. (3.16) relates the extended natural modes of an element to its nodal displacements, defined in the local Cartesian frame. These must be transformed into a global frame and assembled. The operation involves standard finite element procedures and is not discussed here. As a result, the compatibility condition for each element e can be written symbolically as Lb q e ¼ Ce U;

e ¼ 1; . . . ; N ;

ð3:17Þ

where L, the matrix defined in Eq. (3.15), is the same for all elements, b q e is the vector of extended natural modes for the element under consideration, Ce is the assembled compatibility operator for element e and U is the vector of free parameters, i.e., assemblage is meant to enforce not only interelement continuity but also displacement boundary conditions. Standard finite element procedures are also employed to construct the vector of nodal forces. Those equivalent to the basic loads are denoted by R and one can write PðuÞ ¼ Rt U:

ð3:18Þ

4. Power of dissipation for the TRIC element Consider a shell made of von MisesÕ material. Its power of dissipation per unit volume reads, in Cartesian coordinates sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ffi 
1=2 4 4 b ¼ r0 e2 þ ex ey þ e2y þ e2xy þ e2zx þ e2zy D etm hm em þ ets es ; ¼ r0 ð4:1aÞ 3 x 3 with 2

2 2 hm ¼ 4 1 3 0

1 2 0

3 0 0 5: 2

ð4:1bÞ

The condition that plastic flow be isochoric plays no role, since it merely establishes ez ¼
ex
ey and defines the single strain rate component that does not affect the formulation. In other words, any vectors em , es are plastically admissible. The power of dissipation is expressed in terms of the finite element extended modes by substituting, in sequence, the relations introduced in the preceding section. In terms of natural strains, one writes qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi b ¼ r0 ð4:2aÞ D etm Hm em þ ets Hs es : Matrix Hm is immediately computed by substituting Eq. (3.6) for em , to obtain Hm ¼ Ttm hm Tm :

ð4:2bÞ

3756

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

To define Hs , the contribution of transverse shear to the power of dissipation is computed as the average of the values provided by Eqs. (3.8). Namely  4   4 t 41  t es es ¼ es1 es1 þ ets2 es2 þ ets3 es3 ¼ ets Tts1 Ts1 þ Tts2 Ts2 þ Tts3 Ts3 es : 3 33 9 Hence  4 t T Ts1 þ Tts2 Ts2 þ Tts3 Ts3 : ð4:2cÞ 9 s1 The procedure is similar to that used in [17] to compute the shear strain contribution to the natural elastic stiffness of the element. Note that matrix Hs is positive definite for any non-degenerate triangle. By introducing Eqs. (3.9), the expression (4.2) of the power of dissipation becomes pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b ¼ r0 A þ Bz þ Cz2 D ð4:3Þ Hs ¼

with A ¼ gt Hm g þ tt Hs t;

B ¼ 2gt Hm v;

C ¼ vt Hm v:

ð4:4a–cÞ

The power of dissipation per unit shell surface is obtained by integrating the expression above over the element thickness, to obtain Z t=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dðg; v; tÞ ¼ r0 A þ Bz þ Cz2 dz: ð4:5Þ
t=2

D is a function of the point within the element, defined by its triangular coordinates Li ði ¼ 1; 2; 3Þ, and must be integrated over the element area, to obtain the power of dissipation for the finite element. The relation reads Z Z t=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z De ¼ DðLi Þ dX ¼ r0 A þ Bz þ Cz2 dz dX: ð4:6Þ Xe

Xe


t=2

The dependence on the point can be made explicit by introducing Eqs. (3.11), expressing the shell generalized strains in terms of the extended natural modes b q e of the element under consideration. Then
 qe; ð4:7aÞ A¼b q te b b tg Hm b bg þ b bt b b tt Hs b
 qe; B¼b q te b b tg Hm b bv þ b bg b b tv Hm b

ð4:7bÞ


 qe: C¼b q te b b tv Hm b bv b

ð4:7cÞ

By introducing the matrix
 Ke ¼ b b tg Hm b b tt Hs b b tg Hm b bg þ b bt þ z b bv þ b b g þ z2 b b tv Hm b bv; b tv Hm b q e and, hence one can write A þ Bz þ Cz2 ¼ b q te Ke b Z Z t=2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b q e dz dX ¼ De ðb De ¼ r0 q te Ke b qeÞ Xe

ð4:8Þ

ð4:9Þ


t=2

for all elements e ¼ 1; . . . ; N . Eq. (4.9) does not consider the effects of drilling modes, which are present in vector b q e but do not contribute to matrices b b k and, hence, to the value of De . Introducing corner rotations as additional

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3757

membrane degrees of freedom is possible, but not straightforward [21,22]. In any case, triangular shell elements with 15 nodal parameters are adequate in most cases and drilling degrees of freedom are required only to avoid singularities in the assembled matrix for plates or shells with flat portions [23]. Usually, in elastic formulations, as well as in [17,18] for the TRIC element, drilling modes are given a fictitious natural stiffness, governed by a parameter small enough to leave the solution unaffected but sufficient to rule out possible singularities. The same is done here for dissipation: the contribution of drilling modes is taken as proportional to that of symmetric bending by writing for each element qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dd ¼ r0 qtd Kd qd ; Kd ¼ dat1 Hm a1 : ð4:10Þ The constant d is taken as proportional to an estimate of the average dissipation of a finite element and computed as d ¼ kðs =N r0 Þ, where s is a reasonable a priori guess of the collapse multiplier and N is the number of finite elements. The computational experience gained so far indicates that results are unaffected by values of k ranging from 0.000001 to 0.1. The contribution above is included in the expression (4.8) of matrix Ke as a 3 · 3 diagonal block. Eq. (4.9) now defines, in terms of extended natural modes, the power of dissipation for the TRIC element. As such, it subsumes the constitutive behavior of the rigid–plastic element as a whole. In particular, the vector b eðqÞ ¼ oDe Q ob qe

ð4:11aÞ

defines extended stress modes at the finite element yield limit and the condition


 oDe b t b t b b q e Q eðqÞ
Q e ¼ b
Q e P 0 8b qe qe ob qe

ð4:11bÞ

enforces compliance with the element yield limit, in analogy with Eq. (2.5). Note that an explicit definition of the element yield limit is not required. It must be observed that the condition above merely establishes that stress modes are safe for the finite element, considered as a constituent of the discrete system. This fact does not imply by itself that stresses are locally safe in the original (continuous) shell.

5. The discrete limit analysis problem 5.1. Formulation By virtue of the relations above, the limit analysis problem for the finite element shell model reduces to the minimum problem X s ¼ min qeÞ De ð b ð5:1aÞ ^ qe ;U

e

subject to Lb q e ¼ Ce U; Rt U ¼ 1;

e ¼ 1; . . . ; N ;

ð5:1bÞ ð5:1cÞ

where Eq. (5.1b) enforces compatibility (boundary conditions are included in the definition of U) and the counterpart of Eq. (2.6d) does not appear, since any b q e is plastically admissible. In Eq. (5.1c), R denotes the vector of nodal forces equivalent to basic loads. The problem (5.1) is convex, but its objective function is non-smooth and, hence, is not stationary at solution. Before writing the optimality conditions, it must be observed that De is a global value, referring to

3758

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

a finite element as a whole. In each element, De is either positive and differentiable (when b q e 6¼ 0) or equal to zero and not differentiable (if b q e ¼ 0). It follows that the partitioning of X in the two subdomains Xp and Xr now is replaced by the partition of the set E of the N finite elements into the two subsets Ep of the p 6 N elements that undergo plastic flow and Er of the N –p elements that keep rigid in the collapse mechanism. The Lagrangean function of problem (5.1) reads X X Lð b q e ; U; Qe ; jÞ ¼ qeÞ
De ð b Qte ðLb q e
Ce UÞ
jðRt U
1Þ; ð5:2Þ e

e

where vectors Qe ðe ¼ 1; . . . ; N Þ and the scalar j are Lagrangean multipliers. Since L is differentiable with respect to U, the optimality conditions read ! X t t dU L ¼ dU Ce Qe
jR ¼ 0 8dU; ð5:3aÞ e

q te Lt Qe P 0 dq L ¼ dDe
db

8db qe;

e ¼ 1; . . . ; N

ð5:3bÞ

in addition to the original constraints (5.1b and c). By following the same path of reasoning as in Section 2, it can be recognized that Eqs. (5.3) establish X Cte Qe ¼ jR ð5:4Þ e

and L t Qe ¼ db q te



oDe ob qe

if e 2 Ep ;

 oDe t
L Qe P 0 oðdb qeÞ

ð5:5aÞ

8db qe

if e 2 Er :

ð5:5bÞ

If the Lagrangean multiplier vectors Qe are interpreted as the natural (or generalized) stresses of the finite elements, Eq. (5.4) appears as the equilibrium equation for the discrete shell model. Moreover, the virtual btb velocity equivalence Qte qe ¼ Q e q e and Eq. (3.15) define the extended stress modes as b e ¼ L t Qe : Q

ð5:6Þ

Comparison with Eqs. (4.11) shows that Eqs. (5.5) imply compliance with the yield limit both in plastic and rigid elements. It can be concluded that the optimal value j of the Lagrangean multiplier j is a load multiplier which is statically admissible for the discrete system. Also, the discrete counterpart of Eq. (2.11) reads XZ XZ X t t b t b j ¼ j R U ¼ Qe qe dX ¼ De ð^qe Þ ¼ s: ð5:7Þ Q e q e dX ¼ e

Xe

e

Xe

e

Note that the optimal value of problem (5.1) defines the collapse multiplier for the discrete system. If a fully compatible finite element model were used, s would bound from above the collapse multiplier of the continuous shell. The TRIC element, however, is non-conforming, in that interelement continuity is not enforced completely, and the upper bound nature of the result cannot be a priori assessed. 5.2. Outline of the solution strategy The procedure employed for the numerical solution is summarized. Let us suppose, for illustration purposes, that all elements undergo plastic flow, so that L is differentiable everywhere and must be sta-

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3759

tionary at solution. Then, Eq. (5.3b) also hold as equalities and establish oDe =ob q e ¼ Lt Qe 8e. Because of the expression (4.9) of the power of dissipation for the finite element, one can write oDe b eb qe; ¼H ob qe

ð5:8aÞ

where b e ðb q e Þ ¼ r0 H

Z

Z

t=2
t=2

Xe

1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ke dz dX b qe q te Ke b

ð5:8bÞ

is a symmetric matrix depending on b q e , which is present in the square root term. On this basis, the optimality conditions reduce to the system of equations b eb q e ¼ Lt Qe ; e ¼ 1; . . . ; N ; H ð5:9aÞ Lb q e ¼ Ce U e ¼ 1; . . . ; N ; X Cte Qe ¼ jR;

ð5:9bÞ ð5:9cÞ

e

Rt U ¼ 1:

ð5:9dÞ

The problem is non-linear because Eq. (5.9a) are so and must be solved by subsequent iterations. At each iteration j, the conceptually simplest strategy goes through the following steps: b ej are evaluated on the basis of the vectors b 1. Matrices H q eðj
1Þ computed in the previous step (indications on the integration procedure are given in Section 6.1). 2. From Eq. (5.9a–c) one obtains, in sequence b
1 Lt Qe ; b q ¼H ð5:10aÞ e

ej

Qe ¼ Hej Ce Uðjþ1Þ

b
1 Lt  ; with Hej ¼ ½L H ej
1

" Hj Uðjþ1Þ ¼ jR

with Hj ¼

X

ð5:10bÞ

# Cte Hej Ce

:

ð5:10cÞ

e

3. Eq. (5.10c) is solved for Uðjþ1Þ under the condition (5.9d). This is accomplished by writing 1 jðjþ1Þ ¼ t ; Uðjþ1Þ ¼ jðjþ1Þ U : U ¼ H
1 j R; R U

ð5:11Þ

4. The result is substituted into Eqs. (5.10a and b) to obtain the vectors b
1 Lt Hej Ce Uðjþ1Þ : b q ej ¼ H ej

ð5:12Þ

On this basis, the subsequent iteration can be started. As Eq. (5.7) establishes, at solution the Lagrangean (load) multiplier j coincides with the value of the objective function. Therefore, the iteration process terminates when the equality jjþ1 ffi

N X

De ðb q ej Þ

e¼1

is satisfied to within a prescribed accuracy.

ð5:13Þ

3760

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

As described, the procedure does not consider that often the collapse mechanism entails plastic flow in some elements only. This occurrence is accounted for by means of the device established in previous work, which is now briefly outlined, referring to [6] for details. The procedure is started with a tentative vector U defined so as to induce plastic flow in all elements. At each iterative step the dissipation power is computed separately for each element and when it gets smaller than a prescribed tolerance, the relevant element is predicted to be rigid in the final mechanism. For that element, generalized strain rates must vanish and only rigid body motions survive in subsequent computations. The condition qe ¼ Ce U ¼ 0

ð5:14Þ

provides some constraints among the components of vector U, which can be replaced by a smaller size vector U1 by writing U ¼ G1 U1 . Then, the iteration process is continued with the rigid element ignored and the operation is repeated whenever the dissipation power of a new element gets sufficiently small. When the rth element is predicted as rigid, one writes Ur
1 ¼ Gr Ur and, hence with Gr ¼ Gr
1 Gr ¼ G1 G2    Gr :

U ¼ Gr Ur

ð5:15aÞ

Subsequent computations go through the same steps as before, except that U is replaced by the reduced vector Ur and the matrix in Eq. (5.10c) is replaced by " # X t t Hj ¼ Gr Ce Hej Ce Gr ; ð5:15bÞ e2Ep

where Ep is the current set of plastic elements. The procedure aims at identifying the finite elements which are not involved in the collapse mechanism, gradually transferring them from the set Ep of plastic elements (initially including all elements) to the (initially empty) set Er of rigid ones. A sequence of systems with a progressively decreasing number of elements and of free nodal parameters is thus considered. Each system consists of plastically deforming elements only, and its dissipation power is stationary at solution. An inherent limitation of the procedure is that an element cannot be removed from the rigid set Er once it has been introduced in it. The structure is progressively modified by the addition of constraints and, whenever this happens, the collapse multiplier of a different structure is searched. Strictly speaking, the solution obtained should be regarded merely as a kinematically admissible value, bounding from above the collapse multiplier of the original discrete system. However, the numerical experience gained so far seems to indicate that no wrong elements enter the rigid set. 5.3. Computation of matrix H Matrices H and He have the role played, in elastic analyses, by the overall stiffness matrix of the system and of the natural stiffness matrices of the finite elements. Once the latter are known, H is built by means of standard assemblage procedures (the introduction of rigid elements entails minor modifications only). However, constructing matrices He requires, for each element and for each iteration step, the inversion of b e and the subsequent inversion of the 12 · 12 matrix H
1 ¼ L H b
1 Lt . It is now demthe 15 · 15 matrix H e e onstrated how this operation can be reduced to a single inversion of a 3 · 3 matrix, with significant computational saving. be Note first that, in the strain model, shear modes are uncoupled from the remaining ones, i.e., matrix H has the form 2 3 H0 0 b e ¼ 4 1212 123 5 ð5:16Þ H 0 Hs 312

33

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3761

where Hs is the shear modes contribution. Matrix H0 , subsuming the effects of axial, drilling and both bending modes, is partitioned as follows: 2 3 h00 h0A ð5:17Þ H0 ¼ 4 99 93 5; hA0 hAA 33

39

where hAA is the contribution of antisymmetric bending (actually, matrix h00 is not full, since the contribution of drilling modes is uncoupled). From Eq. (5.16) one obtains 2 3 0 H
1 0 123 7 b
1 ¼ 6 H 4 1212 5: e 0 H
1 s 312

33

By virtue of the expression (3.15) of matrix L one writes 2 3 0 0 99 93
1 b
1 t 4 5: H
1 e ¼ L H e L ¼ H0 þ 0 H
1 s 39

ð5:18Þ

33

He is obtained by inverting the matrix above, which can be done by exploiting the formula for the inversion of a modified matrix [24]. This establishes that, given a symmetric matrix   0 0 M ¼ M0 þ ; ð5:19aÞ 0 m where M
1 0 ¼



a bt

b c

 ð5:19bÞ

is known, M
1 is computed as   b
1 t
1
1 M ¼ M0
½p m b c

c

with p ¼ I þ mc:

ð5:19cÞ

Let M ¼ H
1 e , which implies, as comparison between Eqs. (5.18) and (5.19a) shows M0 ¼ H
1 0 ; Hence

M
1 0

m ¼ H
1 s :

ð5:20Þ

¼ H0 and (see Eqs. (5.19b) and (5.17)) b ¼ h0A ;

a ¼ h00 ;

c ¼ hAA :

It follows p ¼ I þ mc ¼ I þ also write

H
1 s hAA .

ð5:21Þ By exploiting the formula for the inverse of a matrix product, one can


1


1
1 p
1 m ¼ ½I þ H
1 s hAA  Hs ¼ ½Hs ðI þ Hs hAA Þ


1


1

¼ ½Hs þ hAA  :

ð5:22Þ

Hence He ¼ M  ¼


1



 h0A
1 ¼ H0
½Hs þ hAA  ½ hA0 hAA

h00
h0A ðHs þ hAA Þ
1 hA0 hA0
hAA ðHs þ hAA Þ
1 hA0

hAA 

h0A
h0A ðHs þ hAA Þ
1 hAA hAA
hAA ðHs þ hAA Þ
1 hAA



and only the inversion of the 3 · 3 matrix Hs þ hAA is required to construct He .

ð5:23Þ

3762

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

6. Computational aspects and comments 6.1. Details on the integration procedure To construct the problem to be solved in each iteration, the integrals (4.5) over the thickness and (4.6), (5.8b) over the element area must be computed. Numerical integration is spontaneous, but both integrals involve square root terms, with possible singularities in the integration domain; in this situation, Gauss integration exhibits poor accuracy even if a certain number of integration points is used. When integrating over the element area the problem has limited relevance. At collapse, most elements are completely plastic and no singularities show up in the integration domain. The small number of elements which are partially plastic are located on the boundary between plastic and rigid zones, dissipate little if at all and contribute marginally to the overall dissipation of the shell, so that possible numerical inaccuracies do not affect the result. On the other hand, the singularity in the square root in Eq. (4.5) is likely to occur within the shell thickness and standard numerical integration methods are not adequate. For homogeneous shells, closed form integration of (4.5) is possible. Note first that, matrix Ke being t positive definite, the equation A þ Bz þ Cz2 ¼ eðzÞ Ke eðzÞ ¼ 0 either has no real roots or two coincident roots, the latter case occurring if the three components of e simultaneously vanish for the same value of z (not necessarily within the plate). This implies 4AC
B2 P 0:

ð6:1Þ

In this situation one has (see, e.g., [25])  t=2 2Cz þ B pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi 4AC
B2 2Cz þ B pffiffiffiffi arcsinh pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D ¼ r0 A þ Bz þ Cz þ 4C 8C C 4AC
B2
t=2 with the second term vanishing when Eq. (6.1) holds as an equality and pffiffiffi D ¼ r0 t A if B ¼ 0; C ¼ 0:

ð6:2aÞ

ð6:2bÞ

At this point, the integral over the element area of the shell dissipation is written as (constant thickness is assumed, besides homogeneity) Z De ðge ; ve ; te Þ ¼ DðA; B; CÞ dX; ð6:3Þ Xe

where A, B and C are defined by Eq. (4.4) as functions of ge ; ve ; te and by Eq. (4.7) in terms of the extended modes b qe. The solution procedure outlined in Section 5 requires the computation, at each iteration step, of matrix b e in plastic elements, defined by Eq. (5.8a) so that H oDe b eb qe: ¼H ob qe By exploiting the chain rule for derivatives and by using Eqs. (3.11), one can write t  t  t   Z  Z  oDe oge oD ove oD ote oD t oD t oD t oD b b b ¼ þ þ bt bg þ dX ¼ þ bv dX: oge ove ote ob qe ob q e oge ob q e ove ob q e ote Xe Xe Because of Eqs. (4.4) and (3.11) one also has



 oD oD oA oD oB oD oD oD b oD b g þ v ¼ 2r0 Hm qe; ¼ þ ¼ 2r0 Hm bg þ bv b oge oA oge oB oge oA e oB e oA oB

ð6:4Þ

ð6:5Þ

ð6:6aÞ

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775



 oD oD oB oD oC oD oD oD b oD b ge þ ve ¼ 2r0 Hm qe; ¼ þ ¼ 2r0 Hm bg þ bv b ove oB ove oC ove oB oC oB oC oD oD oA oD oD b te ¼ 2r0 Hs qe: ¼ ¼ 2r0 Hs bt b ote oA ote oA oA b e reads Hence, matrix H  Z  oD bt oD bt oD bt t b b b b b b ð b Hm b g þ b t Hs b t Þ þ 2 H e ¼ 2r0 b Hm b v þ b Hm b v dX: oA g oB g oC v Xe

3763

ð6:6bÞ

ð6:6cÞ

ð6:7Þ

Note that the functions A, B and C (and, hence, the derivatives of D with respect to them) are configuration dependent. At each iteration, their values are computed for the configuration obtained in the previous step. The integration over the element area of both (6.3) and (6.7) is performed numerically. The use of three integration points located at Lg ¼ 2=3, Lh ¼ Lk ¼ 1=6 ðg ¼ 1; 2; 3; h 6¼ k 6¼ g 6¼ hÞ, with weights Wg ¼ ð1=3ÞXe , turns out to be adequate. 6.2. Antisymmetric bending and transverse shear The description given in Section 5 outlines the conceptual procedure, but computational convenience demands some modifications. The construction of matrices He requires the inversion of small (3 · 3) matrices and a few substitutions and matrix products. These operations are not computationally demanding in themselves, but must be performed for each plastic element at each iteration step. This suggests the search for approximations reducing the computational burden without jeopardizing the quality of results. The problem arises from the fact that, as Eq. (3.14) shows, antisymmetric bending and transverse shears contribute to the same natural modes. This condition may be enforced in an approximate way by writing qAb ¼ lqA ;

qs ¼ ð1
lÞqA

ð0 6 l 6 1Þ:

ð6:8Þ

Eq. (6.8) assumes that the ratio between bending and transverse shear contributions is the same for each antisymmetric mode within an element. l is a bending factor, in that l ¼ 1 implies the absence of transverse shears. If this assumption is introduced, the expression (4.9) of the element dissipation power becomes Z Z t=2 pffiffiffiffiffiffiffiffiffiffiffiffiffi qte Ce qe dz dX; ð6:9Þ De ¼ r0 Xe


t=2

where 2 Ce ¼ 4

lK0A

K00 99

93

lKA0

l2 KAA þ ð1
lÞ2 Ks

39

3 5

33

b e in (5.16) and (5.17), namely and matrix Ke is partitioned as H 2 3 K00 K0A 0 93 6 99 93 7 K K 0 7: 6 A0 AA Ke ¼ 4 39 33 33 5 0 0 Ks 39

ð6:10Þ

33

ð6:11Þ

33

The expression (6.9) of De now depends on the bending factor, which appears as an additional parameter, and the power of dissipation must be minimized also with respect to it. Let us assume that the bending factor le , constant within each element, has no correlation with the bending factors of adjacent elements. Then, minimization can be performed separately for each element by writing

3764

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

oDe ¼ 2r0 ole ¼

Z

Z

Xe

t=2


t=2

2r0 ðle qtA AqA

  1 pffiffiffiffiffiffiffiffiffiffiffiffiffi le qtA ðKAA þ Ks ÞqA
ðqtA Ks qA
qt0 K0A qA Þ dz dX qte Ce qe
qte Bqe Þ ¼ 0;

ð6:12Þ

where qe was partitioned as suggested by the partitioning Eq. (6.10) of Ce , namely qte ¼ f qt0

qtA g;

qt0 ¼ f qtm

qtd

qtSb g

ð6:13Þ

and the following matrices were introduced Z Z t=2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi ðKAA þ Ks Þ dz dX; A ¼ 33 qte Ce qe Xe
t=2 2 3 Z Z t=2 0
K0A 1 1 99 93 5 pffiffiffiffiffiffiffiffiffiffiffiffiffi 4 B ¼ dz dX: 1212 2 Xe
t=2 qte Ce qe
KA0 2Ks 39

ð6:14Þ

33

Note that the matrices above depend on le , affecting the square root terms. At each iteration, they are computed on the basis of the value obtained in the previous step, so that from Eq. (6.12) one obtains le ¼

qte Bqe : qtA AqA

ð6:15Þ

This permits the evaluation of matrices Ce and, hence, the expression of the element dissipation power Eq. (6.9) directly in terms of the element natural modes qe . The operations described in Section 5.2 for the solution of each iterative step now become: 1. Compute the matrices Z Z t=2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi Ce dz dX Hej ðqe Þ ¼ r0 t qe Ce qe Xe
t=2 on the basis of the vectors qeðj
1Þ obtained in the previous step and of the value (6.5) of le . 2. Assemble the matrix " # X t Hj ¼ Ce Hej Ce

ð6:16aÞ

ð6:16bÞ

e

of the linearized problem Hj Uðjþ1Þ ¼ jR for the jth step. 3. Solve Eq. (5.10c) for Uðjþ1Þ under the condition (5.9d), by using Eqs. (5.11). 4. Compute the vectors qej ¼ Ce Uðjþ1Þ

ð6:16cÞ

and go to the subsequent iteration. 6.3. Comparison with the original version of the TRIC element The formulation presented simply amounts at rephrasing for rigid–plastic behavior that given in [17–19]. To preserve its effectiveness, modifications were reduced as much as possible. However, the different context demanded some changes, which are now discussed. To establish comparison, it must be noted that matrices He play the same role as the elastic natural stiffness of individual elements, while H corresponds to the stiffness matrix of the assembled discrete system.

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3765

With this in mind, it is immediately recognized that the way of including transverse shear contributions and of accounting for drilling modes are essentially the same. When constructing matrix H some artifice is required to handle the strain modes qA , containing both antisymmetric bending and transverse shear, which contribute with a priori unknown fractions. The procedure described in Section 5.3 is different from the original proposal. However, it does not entail any approximation and the computational burden looks comparable. Actually, the procedure could be applied for elastic analyses as well, even if advantages are not evident. In any case, it is worth observing that it does not require that membrane strains due to antisymmetric bending be artificially uncoupled from those induced by axial and symmetric bending modes. The procedure is feasible computationally and was applied in some cases. However, the iterative nature of the solution process requires that matrix H be constructed and the problem solved a number of times. Moreover matrices, being configuration dependent, must be updated at each step of the iteration process, motivating the introduction of additional simplifying assumptions. To this end, the exact procedure presented in Section 5.3 is replaced with that described in Section 6.2. This entails some drastic approximations, of more or less empirical nature, and their consequences are hard to assess in general terms (in particular, the assumption of values of le independent in each element entails an additional relaxation of interelement continuity). A few problems solved with both methods produced negligible differences in the results. However, it must be observed that, in contrast to the original formulation, only homogeneous shells are considered here. In this situation, the effects of transverse shears are minor and the bending factors turn out to be close to unity throughout, drastically limiting interelement discontinuities induced by strongly different values in adjacent elements. At this stage, it seems too daring state that the approximation is adequate for laminated shells as well.

7. Numerical examples 7.1. Sequential limit analysis The procedure, used to trace the post-collapse curves of plates and shells, is as follows. The limit analysis problem for the rigid–plastic structure is first solved. Besides the value of the collapse multiplier s, the solution provides the vector U of the nodal velocities corresponding to the collapse mechanism. They are transformed into displacements by multiplying them by a fictitious time Dt, defined by the condition
; maxkUDtk ¼ U ð7:1Þ nodes


is an assigned value. The displacements so obtained are added to the previous nodal coordinates, where U to produce the mesh for the subsequent limit analysis problem.
governs the step amplitude. Clearly, the smaller is U
, the more accurate is the solution. The value of U However, each step solves a standard limit analysis problem, which does not suffer of possible instabilities such as those experienced by incremental analyses when the response is unstable (decreasing post-collapse behavior). In other words, the selection of a suitable value is not conditioned by numerical stability requirements and comparatively large steps are adequate in most cases. To follow the post-critical response, a few problems are to be solved. Moreover, since the mechanism changes during the evolution, each problem must be started from scratch, with all elements potentially deformable. This raises the question of computational efficiency, in particular of the convenience of eliminating rigid elements. Whenever a new element is found to be rigid, the rigidity constraints must be enforced, which entails additional computations. However, after few iterations the number of plastic elements is decreased noticeably and the solution speeds up significantly. The procedure assumes that the elimination of rigid elements is legitimate within a single step. Again, each step involves a standard limit

3766

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

analysis solution and, as it was discussed at the end of Section 5.2, such an expedient, even if hard to justify on a theoretical ground, can be considered correct on the basis of computational experience. Several problems, referring to plates and shells of various shape and subject to different loading, were solved and complete results are discussed in detail in [26]. Only two of them are presented here to demonstrate the effectiveness of the proposed model and solution procedure. 7.2. Simply supported plate A first set of examples refers to square (side a, thickness t), homogeneous plates, simply supported on their boundary and subject to uniform transverse pressure, with basic value p¼4

M0 ; a2

ð7:2Þ

where M0 ¼ r0 t2 =4, r0 being the uniaxial yield limit. Thickness variations are accounted for by means of the slenderness parameter a ð7:3Þ b¼ : t The performances of the rigid–plastic version of the TRIC element were tested first by computing the collapse multiplier of a Kirchhoff plate. The exact solution is not known, but available results bracket it rather stringently within the interval ð7:4Þ

6:216 6 s 6 6:256:

The first result was obtained numerically by Hodge and Belytschko [27]: even if more than 30 years old, it was still quoted as the best lower bound available in a comparatively recent collection compiled by Save [28]. The upper bound corresponds to a finite element solution based on a conforming Kirchhoff element obtained in [8]. Computations were performed with different meshes involving from 8 to 148 elements for a quarter of plate (24–445 degrees of freedom), on purpose arranged in a rather irregular pattern. As slenderness ratio, b ¼ 100 was assumed, corresponding to a plate thin enough to make transverse shear effects negligible. Results converged toward the value s ¼ 6:220, contained in the interval (7.4), which was obtained with the finest mesh. It is also of interest to compare the results with those obtained in [8], where the same problem was solved by using square elements, one based on Kirchhoff model (see [23]), the second resting on Mindlin–Reissner formulation, with 12 degrees of freedom. Some figures are indicated in Table 1. Obviously, Kirchhoff model is more efficient and a smaller number of degrees of freedom is required to reach convergence. Results obtained with the present approach compare favorably with those produced by the Mindlin plate element, since equivalent meshes give equivalent differences with respect to the converged Kirchhoff value and the TRIC element appears more flexible. In addition, TRIC applies to any shells, while the Mindlin element Table 1 Numerical results for the Kirchhoff plate Degrees of freedom

Kirchhoff [8]

66 102 209 217 445 457

6.258 6.256

Mindlin [8]

Present

6.286 6.178 6.220 6.281

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3767

Fig. 4. Collapse multiplier as function of the slenderness ratio.

was specifically devised for the problem at hand and hardly can be used in a different context. It must be observed, however, that the TRIC element is non-conforming and convergence is attained from below. Therefore, in spite of the kinematic approach used, results cannot be considered as upper bounds. Values of the slenderness parameter b ranging from 1 to 100 are next considered to assess the element performances with different transverse shear influence. The computed collapse multipliers are compared in Fig. 4 with those obtained in [8] with the square Mindlin element and with the following upper bound estimate 8 < 6:256 if b P 2:709; ^s ¼ p4ffiffiffi ð7:5Þ b if b 6 2:709: : 3 The first value refers to a Kirchhoff plate, the second assumes that the plate undergoes a rigid vertical translation with shear dissipation along its edges only. For thick plates, the collapse multipliers provided by the TRIC element seem somewhat too low, but results get very good already for b ¼ 4. To this purpose, it must be observed that b ¼ 1 corresponds to a cube and very small b to solids that hardly can be considered as plates, even if thick. The post-collapse behavior of a thin plate is next examined by sequential limit analysis. A quarter of the plate is discretized with 98 elements, involving 298 degrees of freedom. For comparison, the equivalent elastic–plastic problem was solved incrementally with the code ABAQUS [29], by employing two of the available shell elements. The first, labeled as S3R, is triangular and, if the same mesh as for sequential limit analysis is adopted, gives rise to a problem of comparable size; the second, named S9R5, is definitely more refined, in that a subdivision into 64 elements produces 1734 degrees of freedom. Table 2 lists the significant figures for the meshes used in computations. The following geometrical and mechanical properties were adopted a ¼ 20 mm;

h ¼ 0:2 mm;

r0 ¼ 200 MPa;

E ¼ 200 GPa;

m ¼ 0:3

ð7:6Þ

(elastic properties are required by ABAQUS analyses only). The post-collapse curves obtained are depicted in Fig. 5(a), where the load multiplier is plotted, up to a value more than 30 time the collapse level, versus the (dimensionless) central displacement. The solid line refers to sequential limit analysis, while elastic plastic solutions are dashed. Curve 1, referring to the S3R triangle, appears exceedingly stiff. A good

3768

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

Table 2 Models and solution data for the post-collapse analysis of the square plate

SLA ABAQUS ABAQUS

Element type

Number of elements

Nodes

Degrees of freedom

Number of steps (1)

CPU time (1)

TRIC S3R S9R5

98 98 64

64 64 289

298 384 1734

15 47 58

23 00 29 00 0 1 24

00

(1) for a maximum displacement of 0.2a.

(b) 120

(a) 120

80

load multiplier

load multiplier

100

1 2

40

in plane constraints

80 60 simply supported

40 20

β = 100 0 0.00

0.05

0.10

0.15

0.20

central displacement / side length

β = 100 0.25

0 0.00

0.05

0.10

0.15

0.20

0.25

central displacement / side length

Fig. 5. Thin plate: post-collapse evolution.

agreement is obtained with curve 2, referring to the S9R5 element, with minor discrepancies arising from elastic deformations, which sequential limit analysis does not include. However, this result was obtained with nearly six times as many degrees of freedom, stressing the superiority of the TRIC element. Some data referring to the number of load increments and to the CPU time required to attain a displacement of 4 mm ¼ 0.2a are also included in Table 2 (all computations were performed on the same machine, a PCAthlon 1200 MHz processor). It is worth noticing the small number of load steps involved by sequential limit analysis, which can be taken as an indication of the numerical stability of the method. Finally, the effect of the in-plane sliding of supports is analyzed, by comparing previous results with those obtained when in-plane velocities are constrained to be zero on the boundary (Fig. 5(b)). The collapse multiplier does not change but, as well expected, the post-collapse response is much stiffer initially. The slope of the curve does not keep increasing because of a sort of compensation between stiffening geometric effects and the presence of severe membrane forces reducing the bending capacity of the plate. 7.3. Cylindrical shell A cylindrical shell with the geometric and material properties listed in Table 3 is next examined. Both a vertical load and an external pressure are considered. For the first case experimental results are available and are reported in [13], where the problem was studied by sequential limit analysis, as also done in [14]. Both papers employed finite elements exhibiting polar symmetry. The values of the squash load P0 and of

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3769

Table 3 Cylindrical shell properties Radius (mm)

Thickness (mm)

Height (mm)

r0 (MPa)

E (GPa)

P0 (kN)

PE (kN)

k


19.5

1.0

25.0

200

200

24.5

760

0.18

Table 4 Cylindrical shell under axial load: features of the finite element models and computation figures (the number of steps and the CPU time refer to a shortening of 8 mm) Model

Element

A B C D E F G

Axisymmetric Axisymmetric TRIC S3R S9R5 TRIC S9R5

Notes

Method

Number of elements

Degrees of freedom

P0 (kN)

Number of steps

CPU time

Quarter Quarter Quarter Whole Whole

SLA ABAQUS SLA ABAQUS ABAQUS SLA ABAQUS

50 50 272 272 120 400 200

148 148 751 960 3156 1081 5046

24.67 24.50 24.75 24.50 24.45 24.85 24.43

29 1091 18 1529 1475 18 1477

15 00 40 02 00 54 00 230 12 00 250 16 00 0 1 31 00 400 09

00

the elastic buckling load Pp the E for ffiffiffiffiffiffiffiffiffiffiffi ffi axially compressed cylinder are also indicated in Table 3. The dimensionless parameter k
¼ P0 =PE gives an indication on the slenderness of the shell and the value of 0.18 states that it is very stocky. For the axial load case, the different models used in computations are listed in Table 4. Models A and B, referring to axisymmetric shell elements, were employed in [14] and are reported for the sake of completeness. Models C–E consider one quarter of the cylinder only, while the two remaining consider the whole shell, meshed to a slightly lower refinement. The post-collapse curves are depicted in Fig. 6(a) and (b). The first compares results obtained by sequential limit analysis with those presented in [13] and with the experimental curve. The latter refers to a considerably longer specimen (height of 155 mm) and, because of this fact and of the boundary constraints induced by the test equipment, comparatively large

Fig. 6. Cylindrical shell under axial load: post-collapse evolution.

3770

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

Fig. 7. Cylindrical shell under axial load: mechanism evolution up to contact.

displacements are present prior to collapse (the complete curve can be found in [13]). However, the process zone is very nearly the same and, as far as the post-collapse response is concerned, comparison is meaningful. In Fig. 6(a) the post-critical experimental plot was shifted toward left, so as to depart from the collapse load. Fig. 6(b) compares the curve produced by model C (the most refined mesh with TRIC element) with those obtained with incremental analyses by ABAQUS. As for the preceding example, to match the accuracy of TRIC the rather refined S9R5 shell element is required, involving a considerably larger number of degrees of freedom. In the presence of an unstable post-collapse response, the superior stability of sequential limit analysis becomes dramatically evident, as assessed by the significantly reduced number of steps required and the consequent reduction in computation time (see Table 4). This substantiates the fact that, in contrast to incremental computations, sequential limit analysis is insensitive to the rising or decreasing nature of the response. Fig. 7(a) illustrates the evolution of the mechanism from collapse up to contact. The deformed mesh of model C is depicted in Fig. 7(b) for a configuration close to contact. Although no alternative results were available for comparison, the investigation was extended to the case of external pressure, which exhibits a curious aspect: the collapse mechanism and the initial, slightly increasing, post-collapse evolution conform to the polar symmetry of the problem, but at a certain point deformations concentrate in few narrow zones of the shell, detectable in Fig. 8, and the load drops suddenly. In the first phase the most significant displacement is the radial component at mid-height; however, when plastic strain concentrate this quantity becomes meaningless and, for this reason, the post-collapse curves in Fig. 9 refer to the top edge displacement, measuring the cylinder elongation (vertical displacements are constrained to be zero at the bottom and to have the same value at the top). The cylinder actually elongates, because of two effects more than compensating the shortening induced by inflection: the strong Poisson effect associated to incompressible plastic strains and, for significant inflections, the tension induced by the follower nature of the external pressure. Three curves are plotted in Fig. 9. Curve 1 is obtained on the basis of axisymmetric elements and is unable to represent strain concentrations: its behavior is only influenced by the relative weights of Poisson effect, inflection induced shortening and tensile components of follower pressure. Results provided by an

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3771

Fig. 8. Cylindrical shell under external pressure: strain concentration (magnification factor 3).

16

Pressure (MPa)

12

2

1 3

8

4

1 Axisymmetric 2 S.L.A. (TRIC) 3 ABAQUS (S9R5)

0 0.0

0.2

0.4

0.6

0.8

Top edge displacement (mm) Fig. 9. Cylindrical shell under external pressure: post-collapse evolution.

incremental, elastic–plastic analysis, performed with model G in Table 4, follow the same path up to an elongation of about 0.3 mm, when polar symmetry is abandoned and a sudden local collapse occurs (curve 3, dashed). The sequential limit analysis plot (obtained with model F) has the same qualitative trend, even if it is shifted upward by about 5% (curve 2). This result is unexpected and, in fact, anomalous, representing the single significant discrepancy with respect to incremental analyses experienced so far (the discrepancy does not disappear with mesh refinement). In any case, the departure from polar symmetric behavior is predicted at the correct elongation level, stressing the need for general shell elements even for axisymmetric structures. An additional example, referring to an axially compressed square tube with non-symmetric collapse mechanism and post-collapse evolution, confirms the superior efficiency of the proposed procedure when dealing with unstable responses. Its description is omitted because of space limitations, but details can be found in [26,30].

3772

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

8. Conclusions The procedure proposed in this paper to compute the post-collapse curves of plates and shells appears robust and reliable. In particular, comparison between sequential limit analysis used in conjunction with the TRIC element and standard elastic-plastic incremental analyses performed with ABAQUS establishes the following points: 1. to match the performances of the TRIC element, ABAQUS requires models involving 4–6 times as many degrees of freedom; 2. the significant reduction in the number of load steps required by the computations, up to a factor of about 50 when the response is unstable, clearly stresses the superior stability of sequential limit analysis. Point 1 simply assesses that the modifications introduced to adapt the original TRIC formulation to the rigid–plastic context do not jeopardize the element performances. The second point is worth a comment. In a sense, comparison is unfair in that ABAQUS automatically selects the arc length (and, hence, the step amplitude) and the selection criterion might be too severe. In other words, it is likely that comparable results could be obtained with a reduced number of steps if their amplitudes were assigned by the user, as done with the present approach. In any case, the very fact that ABAQUS algorithm spontaneously selects an extremely large number of small steps indicates that incremental computations are, at least potentially, subject to numerical instability problems, a drawback which sequential limit analysis does not suffer. The procedure is flexible enough for several extensions to be envisaged. The simplest one refers to hardening behaviors, which can be included by updating the yield limit on the basis of the amount of plastic strains developed in the previous steps. Dynamic effects, which play a key role in crash analyses, can also be considered by including inertia forces within the framework of rigid–plastic dynamics. These and other developments are presently under study. Acknowledgements This study is a part of the project Molecular level instruments for biomaterial interface design, within the framework of the Large-Scale Computing program of the Politecnico di Milano. The financial support of the Institution is gratefully acknowledged. The finite element code ABAQUS was run at the Department of Structural Engineering of the Politecnico di Milano under academic license. Appendix A The expressions for the natural strain models, Eqs. (3.11), as well as that of the compatibility matrices Ce , are as in the original formulation of the TRIC element and their derivation is clearly explained in [18]. To make the presentation self-contained, they are now cast in the present symbolism and listed, referring to the paper above for details. Eq. (3.11a) is obvious, being evident that membrane strains due to axial modes are constant. Natural curvatures depend on bending modes, governing the transverse (in the z-direction) displacement, which is expresses as 1 wðLi Þ ¼ ð‘a L2 L3 qbSa þ ‘b L3 L1 qSbb þ ‘c L1 L2 qSbc Þ 2 1 þ ð‘a L2 L3 ðL2
L3 ÞqbAa þ ‘b L3 L1 ðL3
L1 ÞqAbb þ ‘c L1 L2 ðL1
L2 ÞqAbc Þ; ðA:1Þ 2

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

3773

where the two addends correspond to symmetric and antisymmetric bending, respectively. By considering local coordinates such as Ya , running along side a of the element (Fig. 1(b)), natural curvatures have the expression va ¼

o2 w oYa2

ðA:2Þ

3 ow X ow oLi ¼ and the relations oYa oLi oYa i¼1 oL2 1 oL3 1 ¼ 0; ¼
; ¼ ; ‘a oYa oYa ‘a 1 oL2 oL3 1 ¼ ; ¼ 0; ¼
; ‘b ‘b oYb oYb 1 oL2 1 oL3 ¼
; ¼ ; ¼0 ‘c oYc ‘c oYc

(and analogous). Then, the chain rule oL1 oYa oL1 oYb oL1 oYc

define the following 2 1 6 ‘a 0 6 6 1 a1 ¼ 6 6 0 ‘b 6 4 0 0

non-zero matrices in Eq. (3.11b) 3 3 2 3 ‘b ‘c
ðL3
L2 Þ
2 L1 L 1 07 7 6 ‘a ‘a ‘2a 7 7 6 ‘a 3 ‘c 7 7 6 7: 6 L
ðL
L Þ
L 07 2 1 3 2 ; a ¼ 2 2 2 7 7 6 ‘b ‘b ‘b 7 7 6 5 5 4 1 ‘a ‘b 3
2 L3 L
ðL
L Þ 3 2 1 2 ‘c ‘c ‘c ‘c

ðA:3Þ

ðA:4a;bÞ

As for shearing modes, it can be observed from the relevant representation in Fig. 2 that each of them implies a linear variation of transverse shear strain rates over the element, with maximum value along a side and vanishing at the opposite vertex. The expression of matrix a3 in Eq. (3.11c) is 2 3 1
L1 0 0 1 a3 ¼ 4 0 ðA:4cÞ 1
L2 0 5: 4 0 0 1
L3 The coefficient is justified by the fact that vector t collects tensorial strain rate components, one half of the engineering shear strains, which in turn are one half of the shear mode amplitude (Fig. 2). The relation (3.12b) between rigid-deformation modes and nodal displacements is considered now. At each vertex i ði ¼ 1; 2; 3Þ of the element, the six degrees of freedom are labeled as indicated in Fig. 3. Rigid body motions consist of three translations and three rotations. The first are q1 ¼

u1 þ u2 þ u3 ; 3

q2 ¼

v1 þ v2 þ v3 ; 3

q3 ¼

w1 þ w2 þ w3 ; 3

ðA:5aÞ

while rigid rotations about the three axes of the local reference system read c1 w1 þ c2 w2 þ c3 w3 b1 w1 þ b2 w2 þ b3 w3 ; q5 ¼
; 2X 2X b1 v 1 þ b2 v 2 þ b3 v 3
c 1 u1
c 2 u2
c 3 u3 q6 ¼ ; 4X q4 ¼

ðA:5bÞ

where ci ; bi are the projections of the element sides on the cartesian axes defined by Eqs. (3.2). Eqs. (A.5) permit the expression of matrix Cq in Eq. (3.12b). The compatibility matrix C is obtained on the basis of the following relations:

3774

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775

Axial modes

2

0




0

8 9 6 < qma = 6 6 cos b qmb ¼ 6 : ; 6 6 ‘b qmc 4 cos c
‘c

sin b ‘b sin c
‘c

cos a ‘a




sin a ‘a

0

0

cos c ‘c

sin c ‘c

cos a ‘a cos b
‘b 0

Symmetric bending modes 8 9 2 0 < qSba = qSbb ¼ 4
sin b : ; qSbc sin c

0 cos b
cos c

sin a
cos a
sin a 0 0 sin b
sin c cos c 0

Antisymmetric bending plus shear modes 8 9 2 0 < qAa = qAb ¼ 4 sin b : ; qAc sin c

0
cos b
cos c

sin a 0 sin c


cos a 0
cos c

sin a sin b 0

38 9 sin a > u1 > > > > > > v1 > ‘a 7 > > > 7> cos b 7< u2 = 7
: ‘b 7 > > v2 > 7> > > > 5> > > > u3 > 0 ; : v3

8 9 > > > #1 > > > 3> > > u1 > cos a > = < > # 2 :
cos b 5 u2 > > > > 0 > > > # > > > > ; : 3> u3

9 8 #1
q4 > > > > > > 3> u1
q5 > > > > >
cos a < = #
q 2 4 5 :
cos b u2
q5 > > > > 0 > > > > #3
q > > 4 > > ; : u3
q5

Drilling modes 8 9 8 9 < qd1 = < w1
q6 = ¼ w
q6 : q : d2 ; : 2 ; qd3 w3
q6 With the aid of Eqs. (A.5), the relations above permit the construction of the 12 · 18 matrix C, expressing deformation modes in terms of the element nodal displacements in the local reference system. Matrices Ce are obtained on this basis, by following the procedure outlined in Section 3.5. References [1] H. Huh, W.H. Yang, A general algorithm for limit solutions of plane stress problems, Int. J. Solids Struct. 28 (1991) 727–738. [2] G.L. Jiang, Nonlinear finite element formulation of kinematic limit analysis, Int. J. Numer. Methods Engrg. 38 (1995) 2775–2807. [3] Y.H. Liu, Z.Z. Cen, B.Y. Xu, A numerical method for plastic limit analysis of 3-D structures, Int. J. Solids Struct. 32 (1995) 1645– 1658. [4] S.W. Sloan, P.W. Kleeman, Upper bound limit analysis using discontinuous velocity fields, Comput. Methods Appl. Mech. Engrg. 127 (1995) 293–314. [5] R.S. Ponter, K.F. Carter, Limit state solutions based upon linear elastic solutions with a spatially varying elastic modulus, Comput. Methods Appl. Mech. Engrg. 140 (1997) 237–258. [6] A. Capsoni, L. Corradi, A finite element formulation for the rigid–plastic limit analysis problem, Int. J. Numer. Methods Engrg. 40 (1997) 2063–2086. [7] E. Christiansen, K.D. Andersen, Computation of collapse states with von Mises type yield condition, Int. J. Numer. Methods Engrg. 46 (1999) 1185–1202. [8] A. Capsoni, L. Corradi, Limit analysis of plates: a finite element formulation, Struct. Engrg. Mech. 8 (1999) 325–341. [9] R.S. Ponter, P. Fuschi, M. Engelhart, Limit analysis for a general class of yield conditions, Eur. J. Mech. A/Solids 19 (2000) 401– 421.

L. Corradi, N. Panzeri / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3747–3775 [10] [11] [12] [13]

[14] [15] [16] [17]

[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

3775

E. Christiansen, O.S. Pedersen, Automatic mesh refinement in limit analysis, Int. J. Numer. Methods Engrg. 50 (2001) 1331–1346. M.R. Horne, W. Merchant, The Stability of Frames, Maxwell, London, 1965. W.H. Yang, Large deformations of structures by sequential limit analysis, Int. J. Solids Struct. 30 (1993) 1001–1013. M. Seitzberger, F.G. Rammerstorfer, On the application of the exact Ilyushin yield surface for plastic collapse analysis of shell structures, in: IASS-IACM 2000, Fourth International Colloquium on Computation of Shells and Spatial Structures, Chania, Greece, June 2000. L. Corradi, N. Panzeri, C. Poggi, Post-critical behavior of moderately thick axisymmetric shells: a sequential limit analysis approach, Int. J. Struct. Stabil. Dynam. 1 (2001) 293–311. J.H. Argyris, Continua and discontinua, in: Proceedings of The First Conference on Matrix Methods in Structural Mechanics, AFFDL, TR 66–88, Dayton, OH, 1966, pp. 11–92. J.H. Argyris, H. Balmer, J.St. Doltsinis, P.C. Dunne, M. Haase, M. Muller, D.W. Scharpf, Finite element method: The natural approach, Comput. Methods Appl. Mech. Engrg. 17/18 (1979) 1–106. J.H. Argyris, L. Tenek, L. Olofsson, TRIC: a simple but sophisticated 3-node triangular element based on six rigid-body and 12 straining modes for fast computational simulations of arbitrary isotropic and laminated composite shells, Comput. Methods Appl. Mech. Engrg. 145 (1997) 11–85. J.H. Argyris, M. Papadrakakis, C. Apostolopoulou, S. Koutsourelakis, The TRIC shell element: theoretical and numerical investigation, Comput. Methods Appl. Mech. Engrg. 182 (2000) 217–245. J.H. Argyris, M. Papadrakakis, L. Karapitta, Elasto-plastic analysis of shells with triangular element TRIC, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3613–3636. P.D. Panagiotopoulus, Inequality Problems in Mechanics and Applications: Convex and Nonconvex Energy Functions, Birkhouser, Boston, 1985. D.J. Allman, A compatible triangular element including vertex rotations for plane elasticity analysis, Comput. Struct. 19 (1984) 1– 8. P.G. Bergan, C.A. Felippa, A triangular membrane element with rotational degrees of freedom, Comput. Methods Appl. Mech. Engrg. 50 (1985) 25–69. O.C. Zienckiewicz, R.L. Taylor, The Finite Element Method, McGraw-Hill, London, 1989. J.S. Prezminiecki, Theory of matrix structural analysis, McGraw-Hill, New York, 1968. M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972. N. Panzeri, Post-critical behavior of shells: a sequential limit analysis approach, Ph.D. thesis, Department of Structural Engineering, Politecnico di Milano, 2001. Ph.G. Hodge Jr., T. Belytschko, Numerical methods for the limit analysis of plates, Trans. ASME, J. Appl. Mech. 35 (1968) 796– 802. M. Save, Atlas of Limit Loads of Metal Plates, Shells and Disks, Elsevier, Amsterdam, 1995. ABAQUS/Standard UserÕs Manuals, Release 6.1, Hibbit, Karlsson, Sorensen Inc., Pawtucket, RI, USA, 2000. L. Corradi, N. Panzeri, A triangular finite element for the study of the post-critical behaviour of shells, in: The Sixth International Conference on Computational Structures Technology, Prague, Czech Republic, 4–6 September 2002, paper CST79 (CD Rom).