- Email: [email protected]
A simple mixed finite element for static limit analysis
Computers and Structures 84 (2006) 1906–1917 www.elsevier.com/locate/compstruc
A simple mixed finite element for static limit analysis X. Yu a, F. Tin-Loi a
b,*
School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, NSW 2006, Australia b School of Civil and Environmental Engineering, University of New South Wales, NSW 2052, Australia Received 5 October 2005; accepted 7 August 2006
Abstract In this paper, we investigate the behavior of a simple mixed finite element for the limit analysis of plane structures. In particular, its ability to overcome incompressibility locking in plane strain situations is investigated. The element is constructed from a piecewise constant displacement field and a piecewise bilinear stress field, and is used within a mathematical programming based discrete representation of the classical static formulation. Several benchmark examples of both plane stress and plane strain situations are solved to illustrate the predictive accuracy and to assess the large-scale capability of the element. The results are compared with those obtained by a recent sophisticated enhanced strain mixed element formulation. 2006 Elsevier Ltd. All rights reserved. Keywords: Limit analysis; Mixed finite element; Plasticity; Collapse loads
1. Introduction Limit analysis is without doubt one of the most widely used instances of the so-called ‘‘direct’’ or ‘‘simplified’’ approaches for nonlinear analyses. As is well-known, it provides useful information regarding the collapse load of a perfectly-plastic structure, while avoiding the computationally expensive effort of an evolutive analysis that follows the whole time history of loading. In recent years, there has been renewed interest in discretized limit analysis, mostly motivated by three aims: application to specific, possibly non-conventional engineering situations; overcoming the volumetric locking behavior encountered in plane strain and 3D problems; and development of efficient computational strategies for solving practically motivated and large size problems. A comprehensive survey of such developments is not possible. We mention only a representative sample of recent work such as those by: Capsoni and Corradi [1] and Capsoni [2] who devel-
*
Corresponding author. Tel.: +61 2 9385 6223; fax: +61 2 9385 3169. E-mail address: [email protected] (F. Tin-Loi).
0045-7949/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.08.019
oped locking-free mixed 4-node plane strain elements with low distortion sensitivity and coarse mesh accuracy; Christiansen and Andersen [3] who devised a new computational approach that exploits the static–kinematic duality property of limit analysis for solving, without locking, large plane strain problems; Christiansen and Pedersen [4] who subsequently improved that method by incorporating a strategy for adaptive mesh refinement; Francescato et al. [5] who developed a combined static and kinematic approach to study the ductile failure of porous metallic materials; Krabbenhoft and Damkilde [6] who used an interior point method for the solution of general nonlinear programming problems arising from a static approach; Lyamin and Sloan [7,8] who also adopted an interior point method for both lower and upper limit analyses using nonlinear programming; and Tin-Loi and Ngo [9] who investigated the robustness of the p-version finite element method within the static framework. It is also notable that most approaches, at variance with the present work, use the kinematic approach to limit analysis. The present paper is primarily concerned with the assessment of a very simple mixed finite element, mentioned briefly in the classical paper of Anderheggen and
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Kno¨pfel [10], as to its capability to overcome incompressibility locking, its suitability for large-scale limit analysis, and also its ability to provide good indications of expected, not necessarily unique, collapse mechanisms. The element is a quadrilateral with a constant displacement field and a bilinear stress field. In spite of its simplicity and potential to overcome locking (as indicated by the favorable ratio of 3/2 of number of stress parameters to number of displacement parameters [10]), this particular element has surprisingly received almost no attention with regards to limit analysis. Its performance in terms of robustness and efficiency is compared with that of the sophisticated mixed enhanced strain element of Capsoni and Corradi [11], which was also implemented for limit analysis in this study, albeit within a static (rather than kinematic [11]) framework. The organization of this paper is as follows. In Section 2, we review the limit analysis model in a form suitable for its direct implementation via the simple element to be investigated. Section 3 deals with some key features of constructing the mixed quadrilateral finite element proposed. A summary of the discretized static problem is then presented in Section 4, before we describe (Section 5) various numerical examples dealing with classical benchmark problems, involving both plane stress and plane strain situations obeying von Mises or Tresca criterion to illustrate the performance of the simple mixed element. We conclude in Section 6. 2. The static limit analysis problem The classical static approach to limit analysis essentially requires that the assumed stress field must be in equilibrium with external loads (governed by a single monotonically increasing load parameter l) and must satisfy the yield condition everywhere within the continuum. The computed collapse limit is then a lower bound to the true limit. In the following, we briefly present a description of the limit analysis model in a form suitable for discretization by mixed elements, including in particular the one which forms the focus of this paper. The derivation of the equilibrium conditions follows the principle of virtual work, in much the same fashion as used by Anderheggen and Kno¨pfel [10] and Christiansen [12]. Typical indices i and j are used to refer to the Cartesian coordinates directions (x, y). Consider a bounded domain V occupied by a rigid perfectly-plastic material continuum. The boundary C of that body consists of two regions Ct and Cr. On Ct, the surface forces lti are prescribed while Cr refers to the region on which the (unknown) reactive surface tractions ri are developed. Volume (body) forces lgi act in V. The limit analysis problem is simply: given a fixed load distribution (gi, ti), it is required to calculate the maximum multiplier for l of this load (collapse load multiplier) that the body can sustain without collapsing. Collapse occurs through a mechanism formation in which the plastic flow is unrestricted. It is also
1907
of interest to find out the flow (mechanism) and stress fields at collapse. The external work rate W of forces (lgi, lti, ri), associated with a virtual plastic flow ui, is Z Z Z ui lgi dV þ ui lti dC þ ui ri dC: ð1Þ W ¼ V
Ct
Cr
The internal work or dissipation rate D for the pair of stresses rij and flow rates is given by Z 1 ðui;j þ uj;i Þrij dV D¼ ð2aÞ V 2 Z Z ð2bÞ ¼ ui rij;j dV þ ui rij mj dC V
C
where a subscripted comma indicates partial derivative, and mj are the components of a unit outward normal vector to the boundary. Note that (2b) has been obtained from (2a) by applying Green’s formula (assuming that the stress is continuous within the domain V). Differentiability of u is then no longer required to calculate D. We can now write the familiar equation of virtual work rate as W ¼ D:
ð3Þ
This equation also represents the equilibrium conditions for the body, as can easily be seen by equating (1) and (2b) to give Z Z ðrij;j þ lgi Þui dV þ ðlti rij mj Þui dC V Ct Z þ ðri rij mj Þui dC ¼ 0: ð4Þ Cr
We further assume that r is continuous across boundaries of domain V so that rij mj ¼ lti on C t ; rij mj ¼ ri on C r :
ð5aÞ ð5bÞ
Using these relations, we can now simplify (4) to give Z ðrij;j þ lgi Þui dV ¼ 0 in V ð6Þ V
which will be used to construct our mixed finite element, underpinned by, as assumed, independent (continuous) stress and (discontinuous) displacement fields. Finally, the stress field rij must satisfy the yield conditions for the assumed material, or in generic form f ðrij Þ 6 0
ð7Þ
which represents a convex (possibly unbounded) set of admissible stresses. The foregoing now allows the maximum limit multiplier ^ to be calculated from the following optimization l problem:
1908
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
^ ¼ max l
9 > > > =
l R
ðrij;j þ lgi Þui dV ¼ 0 in V ; V > rij vj ¼ lti on C t ; > > ; f ðrij Þ 6 0:
subject to :
ð8Þ
uym are arbitrary, we obtain the following series of equations that are decoupled across elements: Z ðrij;j þ lgi Þ dV ¼ 0; m ¼ 1 . . . M: ð11Þ Vm
It should be noted that we have eliminated (5b) since there are no constraints on ri.
By using Green’s theorem, we can rewrite (11) as Z Z rij vj dC þ l gi dV ¼ 0; m ¼ 1 . . . M: Cm
3. The simple mixed quadrilateral finite element Assume, as is conventional, that the two-dimensional domain under consideration is suitably subdivided into a mesh of M quadrilateral elements and N nodes. A generic element is shown in Fig. 1. The displacement field u is assumed to be constant within each element but it does not necessarily satisfy inter-element continuity. In particular, for the mth element, we have ux u¼ ; ðx; yÞ 2 V m ; m ¼ 1 . . . M ð9Þ uy m where Vm is the domain of the mth element. For the stress representation, we assume a bilinear stress field over each element domain. Thus, 2
rxx
3
6 7 4 ryy 5 ¼ rxy
4 X k¼1
2
rxx
3
6 7 N k 4 ryy 5 rxy k
ð10Þ
where the stress parameters [rxx, ryy, rxy]k are defined at each of the four nodes k = 1 . . . 4 of the element (see Fig. 1); and Nk is a local (element-level) weight function, identical to the shape function used in a classical 4-node isoparametric element. Stress continuity within and across elements is thus guaranteed. 4. The discrete limit analysis model We are now in a position to generate the discrete counterpart for the limit analysis problem (8). As a first step in constructing the discretized body force equilibrium relations for the model corresponding to (6), we substitute (9) into (6) and, after noting that uxm and 3
ð12Þ
Vm
We now substitute (10), for each element, into (12) to generate the assembled body force equilibrium relations for the whole finite element model. That is, 9 8 2 3 rxx > > = 4 < X t X ys yr Px 0 xr xs 6 7 ¼0 ; 4 ryy 5 þ l > >2 Py 0 xr xs y s y r ; m : k¼1 rxy k m ¼ 1 . . . M;
ð13Þ
where t is the element thickness (assumed to be constant for each element); x and y used with subscripts r, k, s are coordinates, with r, k, s representing nodal numbers pertaining to the mth element and defined in an anticlockwise convention (i.e. k = 1, 2, 3, 4 corresponds to r = 4, 1, 2, 3 and s = 2, 3, 4, 1); and lPx and lPy are, respectively, the x and y components of element volume forces. The surface force equilibrium equations (or stress boundary conditions) given by (5a) can be discretized in a number of ways, including, for example, a weighted integration method, or a boundary collocation method. When the load distribution is a priori known, the weighted integration method is written as (14a). It can be extended to the more general form of (14b) when the distribution of load is unknown but its total magnitude, however, is known, and thus can be related to the limit load, as in some rigid punch problems. In particular, 82 3 2 39 < rxx rxx = t yb ya 0 xa xb 4 ryy 5 þ 4 ryy 5 ; 0 xa xb y b y a : 2 rxy a rxy b Z tx ¼ l dC; ð14aÞ ty X Ca!b T ½ axx ayy axy i ½ axx ayy axy i ¼ l; ð14bÞ i
where [tx, ty] is the vector of distributed (per unit area) surface forces whose component directions coincide with the coordinates; a and b are two nodes at an element edge that is located at a stress boundary, with a ! b following the
4
uy
q
ux
[σ xx , σyy , σxy ]n Node n
2
y
q
q
b
a
a
a
1
x Fig. 1. Definition sketch of simple quadrilateral finite element.
Fig. 2. Stress boundary conditions: (a) linear continuous distribution along straight boundary; (b) with sharp change of boundary normal; (c) with sharp change of stress magnitude.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
anticlockwise direction for the element (Fig. 2a); Ca!b is the surface from a to b; and ½ axx ayy axy i is a weight vector, the explicit specification of which is determined by individual problems. The boundary collocation method only applies to the case where the load distribution is known beforehand. With this method, the constraints are imposed directly on the stress parameters of the boundary nodes, or rxx rxy vx tx ¼l ; n 2 ½1; N on C ð15Þ rxy ryy n vy n ty where [vx, vy]n is the directional vector of the surface normal at node n. Special care is needed to implement (15) when the normal to the stress boundary is discontinuous (Fig. 2b) or the load distribution q changes sharply (Fig. 2c). For these two cases, two sets of relations (15) need to be written to represent the stress boundary when the boundary node is approached from two opposite directions. In the case that these sets of equations conflict, which is always the case for Fig. 2c and sometimes also for Fig. 2b, one needs to make an arbitrary decision to drop one set from the conflicting pair. Obviously, this will introduce errors. To reduce such errors, we use the local mesh refinement scheme shown in Fig. 3. In a series of preliminary studies, we have tested both the integration and collocation methods and obtained very similar results both in accuracy and performance. Unless otherwise stated, the results reported in this paper are based on the boundary collocation method represented by (15). The yield constraints, typically either the well-known Tresca or von Mises criterion, are enforced at the nodes of the mesh. Since the stress distributions are piecewise bilinear, the yield conditions are automatically satisfied at every point in the discretized domain.
n div = 0
ndiv = 1
n div = 2
q
Fig. 3. Local mesh refinement: (a) before refinement; (b) after refinement with ndiv = 1; (c) after refinement with ndiv = 2; (d) used for sharp change of stress magnitude.
1909
It is useful to recall at this stage these yield criteria for plane situations; r0 is used to indicate the yield stress in simple tension. The material constant in these yield conditions can also be the pffiffiffi shear strength s0 of the material; for von Mises r0 ¼ s0 3 and for Tresca r0 = 2s0. Subscripts M and T denote Mises and Tresca, respectively. For plane strain, 4 2 fM ¼ ðrxx ryy Þ þ 4r2xy r20 ; 3 2
fT ¼ ðrxx ryy Þ þ 4r2xy r20 :
ð16aÞ ð16bÞ
For plane stress, 2
fM ¼ ðrxx ryy Þ þ rxx ryy þ 3r2xy r20 ;
ð16cÞ
fT ¼ maxðjr1 j; jr2 j; jr1 r2 jÞ r0 ;
ð16dÞ
where r1 and r2 are principal stresses. Finally, the discrete limit analysis problem can be set up as the following nonlinear programming problem: 9 ^ ¼ max l l > > > > = subject to: element equilibrium ð13Þ; stress boundary conditions ð14Þ or ð15Þ; > > > > ; yield conditions ð16Þ: ð17Þ It is useful to note that, with the above approach, the requirement of static admissibility is only satisfied in a ^ can only be considweak form. Therefore, the computed l ered to be an approximation, and not a strict bound, to the exact value. It may be argued that a strict bound solution, based on equilibrium elements [7] or compatible displacement elements [8], is more valuable in a limit analysis of a large-scale practical problem. Unfortunately, strict bounds can only be established in idealized situations. Even with the equilibrium element [7], static admissibility is not fully guaranteed. For example, (a) boundary equilibrium can be violated in the case of nonlinear distributions of boundary stresses; (b) inter-element equilibrium can be violated in the case of sharp changes of shear stresses at a boundary; and (c) element equilibrium can be violated in the case of non-constant body force. In addition, the strict bound property is also lost if the material strength has a nonlinear distribution across the domain, or if the material does not satisfy the associated flow rule. In view of these and also because our approach typically converges, as the model is refined, to the true value of collapse multiplier, we do not aim to obtain strict bounds. Obviously, whether the convergence is from above or below will in any case reveal the nature of the computed collapse limit. Finally, we note that solution of (17) will not only pro^, but also an idea of the (posvide the collapse multiplier l sibly non-unique) collapse mechanism obtained from the
1910
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Lagrange multipliers associated with the element equilibrium conditions (13) at optimality. 5. Numerical examples Our pilot numerical implementation consists of the following features: (a) formation of relevant quantities to set up the nonlinear programming problem in MATLAB; (b) production of an appropriate data (text) file for use by a generic GAMS (General Algebraic Modeling System) [13] model that implements the nonlinear programming formulation; and (c) solution by one of the state-of-the-art GAMS solvers (we typically use the well-known nonlinear programming solver GAMS/MINOS [14]). The use of GAMS provides a large number of advantages of which the most important for the work described herein are the availability of a number of industry standard mathematical programming solvers for linear and nonlinear programming, simplicity of modeling, large-scale facility, and automatic differentiation capability. In the following, we summarize some results related to four problems, namely a square plate with central circular hole under biaxial tension, a two-notched tensile specimen, a thick clamped beam, and the well-known Prandtl punch problem. For comparison, a sophisticated quadrilateral mixed element, initially developed for elastic–plastic analysis by Capsoni and Corradi [11], was also implemented within our static approach to limit analysis. Our implementation of this element involves stress parameters defined as unknown variables, and the yield conditions imposed at four vertices for each element; details are given in [15]. In terms of stress and velocity fields, this implementation is in effect the ‘‘dual’’ of that of Capsoni and Corradi [1] who, it must be stressed, used a kinematic approach. However, because of the differences in imposing yield conditions, the two implementations are not equivalent. For instance, Capsoni and Corradi [1] do not allow an element to become partly-plastic (and partly rigid) since no plasticrigid boundary was accommodated within their element. In addition, their iterative algorithm is to progressively add elements, deemed to be non-critical, to a ‘‘rigid’’ set. However, once included in that set, no element can become a candidate for future yielding. These limitations do not apply to our present approach. Therefore, given the same mesh and stopping criterion, the two implementations should yield similar, but not necessarily the same, predictions. In the following examples, unless otherwise noted, we compare the performance of our simple mixed element and with that of the mixed element of Capsoni and Corradi, as implemented by us within a static framework. Of course, we used identical testing conditions, e.g. same mesh, same solver, same feasibility and optimality tolerances (typically 106, the default in MINOS [14]), and same computer. All runs reported were carried out on a Pentium III 800 MHz PC.
5.1. Square plate with central circular hole This first example is a plane stress problem, involving the well-known reference case of a square plate with a hole diameter to plate width ratio of 0.2 (Fig. 4a). We analysed a quarter of the plate modeled with various meshes defined by parameters nh and nl, as in Fig. 4b. When p2 = 0, the problem becomes a uniaxial tension problem for which the exact collapse solution for the plane stress Tresca condition is p1/r0 = 0.8 [16] with the tensile stress at the plate edge arbitrarily distributed so that p1 refers to the averaged tensile stress at the plate edge. Three feasible failure mechanisms, namely in-plane symmetrical shearing, in-plane unsymmetrical shearing and out-ofplane shearing, were reported, any of which leads to the same solution. In our case, in view of the quarter model used, only a symmetrical failure mechanism is possible. We considered four cases as follows: • • • •
Case A: von Mises, p2 = p1, uniform stress. Case B: von Mises, p2 = 0.5p1, uniform stress. Case C: von Mises, p2 = 0, uniform stress. Case D: Tresca, p2 = 0, arbitrary distribution of tensile stress p1.
It is noted that the stress boundary for Case D was implemented based on (14b) because the load distribution is a priori unknown; for other cases, (15) was applied. Tables 1–4 summarize the results of these analyses for different numbers of nodes N and elements M. Some statistics related to the nonlinear programming problem (namely, number of variables nvar, and number of constraints ncon) and an idea of the computational cost (namely, run time in seconds) are also given. Table 5 also compares our results for the von Mises cases with some reported ones in the literature [9,16–21]. The results given in these Tables 1–5 show that: the calculated collapse pressures all converge with mesh refinement; our simple mixed element performed far better in terms of robustness and efficiency than our implementation
p2
p2
p1
p1
p1
nh
p2
nl
1.5 × n l
Fig. 4. Square plate: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 Table 1 Square plate: Case A (von Mises, p2 = p1) nh
nl
N
M
nvar
ncon
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
144 459 1611 6003
1911
Table 4 Square plate: Case D (Tresca, p2 = 0, arbitrary stress distribution) p1/r0
CPU (s)
nh
0.8486 0.8820 0.8906 0.8933
0.2 1.7 37.6 1491.1
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.9003 1.2 8 4 141 116 1973 2116 0.8964 25.0 16 8 513 464 7889 8408 0.8955 961.2 32 16 1953 1856 31,553 33,520 Faila –
nl
N
M
nvar
ncon
p1/r0
139 449 1591 5963
0.7762 0.7944 0.7984 0.7990
CPU (s) 0.2 1.9 23.7 618.1
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 500 537 0.8134 0.9 8 4 141 116 1984 2117 0.8042 9.2 16 8 513 464 7910 8409 0.8016 213.8 32 16 1953 1856 31,594 33,521 0.8006 6485.5
a
Run stopped after about 40 h of CPU time without convergence being achieved. Table 5 Square plate: comparison of limit multipliers Table 2 Square plate: Case B (von Mises, p2 = 0.5p1) nh
nl
N
M
nvar
ncon
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
144 459 1611 6003
p1/r0 0.8810 0.9007 0.9065 0.9086
CPU (s) 0.3 2.4 29.7 758.6
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.9254 1.1 8 4 141 116 1973 2116 0.9152 18.7 16 8 513 464 7889 8408 0.9121 370.2 32 16 1953 1856 31,553 33,520 0.9110 36077.2
Table 3 Square plate: Case C (von Mises, p2 = 0) nh
nl
N
M
nvar
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
Reference
Case A
Case B
Case C
Present work, fine mesh Analytical [16] Tin-Loi and Ngo [9] [p-version finite element] Belytschko [17] [triangular equilibrium element] Groß-Weege [18] [quadrilateral equilibrium element] Corradi and Zavelani [19] [triangular constant-stress element] Nguyen and Palgen [20] [quadrilateral equilibrium element] Liu et al. [21] [boundary element]
0.893 – 0.895
0.909 – 0.912
0.799 0.800 0.803
–
–
0.780
0.882
0.891
0.782
0.767
–
0.691
0.704
–
0.564
0.903
0.915
0.795
5.2. Two-notched tensile specimen
ncon 144 459 1611 6003
p1/r0 0.7640 0.7958 0.7984 0.7991
CPU (s) 0.2 2.4 34.9 969.0
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.8147 1.2 8 4 141 116 1973 2116 0.8057 20.4 16 8 513 464 7889 8408 0.8022 407.3 32 16 1953 1856 31,553 33,520 0.8008 16500.3
of the mixed element of Capsoni and Corradi [11], which is expected as the size of the nonlinear programming problem increases far more rapidly in the latter case; and our predictions for Cases C and D match particularly well the only analytical solution we are aware of. Finally, plots of the failure mechanisms we obtained for each of these four cases are shown in Fig. 5. The differences in the failure mechanisms between Case C and Case D are possibly due to the different yield conditions. For Case D, our failure mechanism agrees well with that of the analytical solution [16].
This type of problem, involving a plane strain doublenotched rectangular specimen under tension (Fig. 6a), is a popular benchmark used for testing volumetric locking in elastoplastic analyses [22]. We chose to carry out limit analyses on ‘‘infinitely’’ long specimens obeying von Mises criterion. Finite element models of a quarter of the plate (Fig. 6b) were constructed for the following three cases: • w/b = 3, h/b = 10. • w/b = 5, h/b = 10. • w/b = 8, h/b = 15. The indicated h/b dimensions were selected after some preliminary sensitivity analyses. Further, as indicated in Fig. 6b, we introduced some local mesh refinement, in the fashion of Fig. 3d, at the notch tip in an attempt to accommodate any stress discontinuities. In order to assess how extensive this refinement should be, we first carried out a series of analyses to investigate the influence of this local mesh refinement (measured through parameter ndiv) on the w/b = 3 case and using our proposed simple mixed element. Table 6 displays the results obtained; nb refers to
1912
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Case A
Case B
Case C
Case D
Fig. 5. Square plate: failure mechanisms for four load cases.
p
Table 6 Notch mesh refinement (w/b = 3, h/b = 10, nb = 10)
nb ( = 2)
b
2b
h
2h
ndiv
N
M
nvar
ncon
p/s0b
CPU (s)
0 1 2 3 4 5 6
3131 3136 3141 3146 3151 3156 3161
3000 3004 3008 3012 3016 3020 3024
9394 9409 9424 9439 9454 9469 9484
9542 9558 9574 9590 9606 9622 9638
3.5082 3.5208 3.5461 3.5506 3.5569 3.5582 3.5595
2416.0 2910.0 1925.6 1623.2 1590.8 1407.5 1439.4
2w
p
p
Fig. 6. Double-notched specimen: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
the number of elements, excluding those in the local refinement zone, underneath the half-width b of the unnotched segment (e.g. nb = 2 in the mesh shown in Fig. 6b). Interestingly, the computing times dropped dramatically with appropriate local refinement. The predicted collapse load factor varied less than about 0.1% when the local refinement was changed from ndiv = 4 to ndiv = 5. We adopted the latter discretization for our subsequent main runs.
Tables 7–9 summarize the limit analyses for the three cases. It is evident that our simple element performed reasonably well in providing collapse limits in all instances, whilst the models constructed using the mixed element of Capsoni and Corradi [11] could not be solved for even low values of nb. Moreover, our simple element did not exhibit any locking for this plane strain situation. We were unable to assess properly the accuracy of our computed limit loads since we could not find any corresponding values reported in the literature. However, the long double-notched tension specimen is equivalent to a rigid punch over a narrow block for which analytical lower and upper bounds are available [23]. Table 10 provides the appropriate comparisons. Plots of the failure mechanisms revealed by our simple finite element based model are given in Fig. 7 (with the
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
1913
Table 7 Double-notched specimen: w/b = 3, h/b = 10, ndiv = 5 nb
N
M
nvar
Current simple mixed element 2 172 140 517 4 558 500 1675 6 1184 1100 3553 8 2050 1940 6151 10 3156 3020 9469 12 4502 4340 13,507 Mixed 2 4 6
ncon
pw/s0b
CPU (s)
550 1738 3646 6274 9622 13,690
4.3495 3.8293 3.6754 3.6016 3.5582 3.5297
2.2 22.5 120.2 498.4 1407.5 3920.5
element of Capsoni and Corradi (present 172 140 2381 2555 558 500 8501 9065 1184 1100 18,701 19,895
implementation) 3.3758 21.9 3.3892 187.5 Fail –
Table 8 Double-notched specimen: w/b = 5, h/b = 10, ndiv = 5 nb
N
M
nvar
Current simple mixed element 2 256 220 769 4 886 820 2659 6 1916 1820 5749 8 3346 3220 10,039 10 5176 5020 15,529 12 7406 7220 22,219 Mixed 2 4 6
ncon
pw/s0b
CPU (s)
810 2738 5866 10,194 15,722 22,450
4.8701 4.5206 4.4109 4.3574 4.3256 4.3046
5.9 93.6 651.2 3392.9 14897.8 38769.1
element of Capsoni and Corradi (present 256 220 3741 4003 886 820 13,941 14,841 1916 1820 30,941 32,879
implementation) 4.1865 52.0 4.1957 7123.1 Faila –
a Run stopped after about 40 h of CPU time without convergence being achieved.
Fig. 7. Double-notched specimen: failure mechanisms for three geometries (notch is located at bottom left hand corner).
Table 9 Double-notched specimen: w/b = 8, h/b = 15, ndiv = 5 nb
N
Current simple 2 552 4 2038 6 4484 8 7890
M
nvar
mixed element 500 1657 1940 6115 4340 13,453 7700 23,671
ncon
pw/s0b
CPU (s)
5.3. Thick clamped beam
1720 6238 13,636 23,914
5.4248 5.2198 5.1531 5.1191
32.1 778.5 6271.0 34631.0
This problem, investigated recently by Capsoni and Corradi [1], involves a thick clamped beam in plane strain under a uniformly distributed pressure loading applied at its top central part, as shown in Fig. 8a. We modeled half of this beam (Fig. 8b) with mesh refinement (ndiv = 5) provided at the top pressure discontinuity point and the two corners of the fixed ends. We adopted a von Mises yield condition.
Mixed element of Capsoni and Corradi (present implementation) 2 552 500 8501 9065 4.9811 950.2 4 2038 1940 32,981 35,045 Faila – a Run stopped after about 40 h of CPU time without convergence being achieved.
notch oriented at the left hand bottom corner of the models).
p
p
2b h
Table 10 Double-notched specimen: comparison of limit multipliers
Present work Upper bound [23] Lower bound [23]
w/b = 3
w/b = 5
w/b = 8
3.530 4.000 3.000
4.305 4.899 3.333
5.119 6.000 3.556
2w Fig. 8. Thick clamped beam: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
1914
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Three geometries were analysed:
Table 13 Thick clamped beam: w/b = 5, h/b = 2, ndiv = 5
• w/b = 1, h/b = 2. • w/b = 2, h/b = 2. • w/b = 5, h/b = 2. The w/b = 1 case is particularly noteworthy since the failure mechanism associated with the exact solution requires shear sliding at the fixed end, and only models capable of displacement discontinuities will be able to simulate this failure mode. Tables 11–13 summarize some of the results obtained; np refers to the number of elements, excluding local refinement, in width b of material (i.e. under load p). Once again, it can be conclusively seen that our simple element is superior, both in terms of robustness and efficiency, as compared with the results provided by our static implementation of the more refined element of Capsoni and Corradi [1]. Contrary to expectation, however, the values of collapse load multiplier do not converge strictly monotonically for the w/b = 5 case (Table 13) as the mesh is refined.
Table 11 Thick clamped beam: w/b = 1, h/b = 2, ndiv = 5 np
N
nvar
ncon
p/s0
CPU (s)
Current 2 4 6 8 10 12
simple mixed element 45 28 136 75 52 226 121 92 364 183 148 550 261 220 784 355 308 1066
136 226 364 550 784 1066
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
0.3 0.6 1.7 3.3 7.6 13.1
M
N
np
M
nvar
ncon
Current simple mixed element 2 110 80 331 4 244 200 733 6 458 400 1375 8 752 680 2257 10 1126 1040 3379 12 1580 1480 4741 Mixed 2 4 6 8 10 12
element of Capsoni 110 80 244 200 458 400 752 680 1126 1040 1580 1480
p/s0
357 775 1433 2331 3469 4847
0.9461 1.0057 1.0142 1.0160 1.0162 1.0160
and Corradi (present 1361 1465 3401 3641 6801 7257 11,561 12,313 17,681 18,809 25,161 26,745
CPU (s) 1.5 8.1 35.5 98.6 242.3 587.4
implementation) 1.0961 9.7 1.0419 44.1 1.0307 173.2 1.0257 678.5 1.0226 3172.7 1.0207 10430.4
The failure mechanisms we obtained are depicted in Fig. 9. It is notable that the shear failure for the w/b = 1 case has been captured properly. The transition of shear failure to bending failure is also clear as w/b increases. Finally, Table 14 compares the present predictions with available [1] analytical bound solutions and computed results. As can be seen, our element predicts the exact solution for w/b = 1. It is noted that, in Table 14, the results attributed to the mixed element of Capsoni and Corradi were obtained by
Mixed element of Capsoni and Corradi (present implementation) 2 45 28 477 503 2.3992 1.5 4 75 52 885 935 Faila – a MINOS exited with ‘‘too many iterations’’ without achieving an optimal solution.
Table 12 Thick clamped beam: w/b = 2, h/b = 2, ndiv = 5
w/b = 1
np
N
nvar
ncon
p/s0
CPU (s)
Current 2 4 6 8 10 12
simple mixed element 80 56 241 136 104 409 224 184 673 344 296 1033 496 440 1489 680 616 2041
255 427 695 1059 1519 2075
1.5118 1.7264 1.7341 1.7395 1.7419 1.7421
0.4 2.8 7.3 21.4 40.8 97.8
M
Mixed element of Capsoni and Corradi (present 2 80 56 953 1021 4 136 104 1769 1889 6 224 184 3129 3333 a
implementation) 1.9604 5.4 1.8319 86.5 Faila –
MINOS exited with ‘‘too many iterations’’ without achieving an optimal solution.
w/b = 2
w/b = 5 Fig. 9. Thick clamped beam: failure mechanisms for three geometries.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 Table 14 Thick clamped beam: comparison of limit multipliers
Present work Upper bound [1] Lower bound [1] Capsoni and Corradi [1]
1915
Table 15 Prandtl’s punch: uniform distribution of vertical contact forces N
M
nvar
w/b = 1
w/b = 2
w/b = 5
np
2.000 2.000 2.000 2.286
1.742 2.000 1.333 1.828
1.016 2.000 0.500 1.071
Current simple mixed element 2 55 40 166 4 189 160 568 6 403 360 1210 8 697 640 2092 10 1071 1000 3214
them through their kinematic approach and are reported in [1], while the results attributed to the mixed element of Capsoni and Corradi in Tables 11–13 are from our present static implementation. Obviously, as is indeed the case, our results do not match those of Corradi and Capsoni since we used, as indicated earlier, a static approach (instead of a kinematic approach), different mesh configurations, a different method for ensuring yield conformity, and possibly also different convergence criteria. 5.4. Prandtl’s punch This classical plane strain problem consists of a semiinfinite rigid plastic medium under a punch load of width 2b (Fig. 10a). The exact collapse load [16], based on the plane strain Tresca criterion, is P = (2 + p)s0b = 5.1416s0b. We constructed a finite element model (Fig. 10b) involving half of a block of medium with dimensions 2w · h. The interaction between the isolated block and the rest of the medium was represented by fixed (no x or y displacement) boundary conditions. Initial trial simulations showed that w/b = 5 and h/b = 2 were large enough to simulate an ‘‘infinite’’ medium. Uniform meshes were used and local mesh refinement was not applied since the collapse load factor is independent of w/b. A plane strain Tresca criterion was used in all analyses. The interactions between the punch and the medium were represented by vertical force distributions and horizontal displacement restraints applied on the block. We considered two cases of vertical force distributions as follows:
ncon 159 555 1191 2067 3183
P/s0b
CPU (s)
5.1001 5.1170 5.1238 5.1277 5.1300
0.4 3.6 31.0 60.7 133.3
Mixed element of Capsoni and Corradi (present implementation) 2 55 40 681 714 5.3494 1.1 4 189 160 2721 2868 5.2485 17.0 6 403 360 6121 6462 5.2158 74.9 8 697 640 10,881 11,496 Faila – a
MINOS exited reporting ‘‘problem is unbounded or badly scaled’’.
Table 16 Prandtl’s punch: arbitrary distribution of vertical contact forces np
N
nvar
ncon
P/s0b
CPU (s)
Current 2 4 6 8 10
simple mixed element 55 40 166 189 160 568 403 360 1210 697 640 2092 1071 1000 3214
157 551 1185 2059 3173
5.1019 5.1194 5.1262 5.1298 5.1319
0.4 3.5 18.6 83.0 135.0
M
Mixed element of Capsoni and Corradi (present implementation) 2 55 40 684 715 5.8490 1.3 4 189 160 2726 2869 5.5069 17.2 6 403 360 6128 6463 Faila – a
MINOS exited reporting ‘‘problem is unbounded or badly scaled’’.
It is noted that the latter arbitrary distribution case is a more appropriate representation for the rigid punch problem but we also considered the uniform distribution as it is one that has been frequently quoted in the literature.
• Uniform distribution (using (15) for the stress boundary). • Arbitrary distribution (using (14b) for the stress boundary).
2P
2b
P
b
h
2w
Fig. 10. Prandtl’s punch: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
Fig. 11. Prandtl’s punch: (a) failure mechanism for uniform pressure distribution; (b) failure mechanism for arbitrary pressure distribution.
1916
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
The results of our limit analyses are summarized in Tables 15 and 16, representing, respectively, the uniform force distribution and arbitrary force distribution. Symbol np refers to the number of elements in contact with half of the punch. Corresponding collapse mechanisms are shown in Fig. 11. As can be seen, our simple mixed element performed well, and there does not appear to be any significant differences between the two load cases with both obtained mechanisms being similar. In line with our expectation, the collapse load predictions are both close to the theoretical value. 6. Concluding remarks Within the framework of a finite element based static approach to limit analysis, we present a simple 4-noded quadrilateral mixed finite element for reliably calculating collapse load multipliers. In particular, it is found that locking caused by volumetric incompressibility does not occur in plane strain problems. Furthermore, accurate solutions can be obtained in reasonable computing time, as demonstrated in the present paper through the results of four benchmark problems. Comparison of its performance with a sophisticated and recently developed [11] enhanced strain element favors our simple element. This is possibly not so surprising since, in our view, the element of Capsoni and Corradi is more suited for the solution of elastoplastic [11] rather than limit analysis [1] problems. Models using the enhanced strain element lead to very large increases in the sizes of the underlying nonlinear programming problems with mesh refinement. Capsoni and Corradi [1], it must be noted, used a kinematic approach and successfully reduced the problem size of their unconstrained but non-smooth optimization problems by identifying possible rigid regions which are then excluded in their optimization phase. Our prototype implementation consists essentially of a MATLAB code to provide appropriate equilibrium equations that a simple GAMS nonlinear programming model would use. As for the solution of the optimization problems, we used the state-of-the-art nonlinear solver GAMS/MINOS (which we found to be faster than another industrial-strength solver GAMS/CONOPT [24]). However, in some instances, we noted that the CPU time still grew very rapidly with increase in problem size. Adoption of the recently developed interior point nonlinear programming solvers is likely to provide an increase in computation speed, as they are based on quasi-Newton methods for which the iteration count is largely independent of the problem size, as has been reported [6–8]. Nevertheless, in the present study, we did not attempt to develop or implement special schemes to solve the optimization problems more efficiently as our primary concern was simplicity of finite element combined with a locking-free behavior. The use of mixed elements does not guarantee, as is well known, a strict bound to the collapse limit. Moreover, the estimates of limit multipliers typically converge mono-
tonically as the mesh is refined, thus providing a clear indication as to whether the actual collapse is being approached from above or below. The present simple mixed element appears, in addition to its acceptable accuracy, eminently suitable for providing good collapse mechanisms, even those involving displacement discontinuities. We did not notice any hourglass type displacement modes, as is often experienced with a 4-node isoparametric element that uses reduced (single centroid point) Gauss integration. Moreover, the stress field provided by the simple element can exhibit hourglass modes, as when the boundaries, subjected to rigid movement or fixed, have unknown stress distributions. However, the stress fields invariably become smooth in the yield region. Finally, worthwhile extensions of the present work would be develop a 3D version of the element, and to investigate adoption of piecewise linearized yield conditions, perhaps in concert with some adaptive scheme, so that a linear program needs to be solved rather than the typically more computationally expensive nonlinear optimization. Of course, the trade-off between increases in problem size with use of the more efficient linear programming codes needs to be carefully assessed. Acknowledgements This research was supported by a grant from the Australian Research Council (ARC). The work reported in this paper was carried out when Dr. Yu, the first named author, was employed, through ARC funds, as a Senior Research Associate at the University of New South Wales. Finally, the very constructive comments of two anonymous referees on an earlier version of this manuscript are gratefully acknowledged. References [1] Capsoni A, Corradi L. A finite element formulation of the rigidplastic limit analysis problem. Int J Numer Methods Eng 1997;40:2063–86. [2] Capsoni A. A mixed finite element model for plane strain limit analysis computations. Commun Numer Methods Eng 1999;15:101–12. [3] Christiansen E, Andersen KD. Computation of collapse states with von Mises type yield condition. Int J Numer Methods Eng 1999;46:1185–202. [4] Christiansen E, Pedersen OS. Automatic mesh refinement in limit analysis. Int J Numer Methods Eng 2001;50:1331–46. [5] Francescato P, Pastor J, Thai TH. Etude du crite`re de plasticite´ des mate´riaux poreux. C R Acad Sci, Se´r IIB 2001;329:753–60. [6] Krabbenhoft K, Damkilde L. A general non-linear optimization algorithm for lower bound limit analysis. Int J Numer Methods Eng 2003;56:165–84. [7] Lyamin AV, Sloan SW. Lower bound limit analysis using non-linear programming. Int J Numer Methods Eng 2002;55:573–611. [8] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Methods Geomech 2002;26:181–216. [9] Tin-Loi F, Ngo NS. Performance of the p-version finite element method for limit analysis. Int J Mech Sci 2003;45:1149–66.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 [10] Anderheggen E, Kno¨pfel H. Finite element limit analysis using linear programming. Int J Solids Struct 1972;8:1413–31. [11] Capsoni A, Corradi L. A mixed finite element model for plane strain elastic–plastic analysis: Part I – Formulation and assessment of the overall behavior; Part II – Application to the 4-node bilinear element. Comput Methods Appl Mech Eng 1997;141:67–93. [12] Christiansen E. Limit analysis of collapse states. In: Ciarlet PG, Lions JL, editors. Handbook of numerical analysis, vol. 4. Amsterdam: North-Holland; 1996. p. 193–312. [13] Brooke A, Kendrick D, Meeraus A, Raman R. GAMS: a user’s guide. Washington, DC 20007: Gams Development Corporation; 1998. [14] Gill PE, Murray W, Saunders MA, Wright MH. MINOS 5.5 user’s guide. Report 83-20R, Department of Operations Research, Stanford University, Stanford, CA, revised 1998. [15] Yu X, Tin-Loi F. Behavior of four finite elements for direct limit analysis. In: Loo YC, Chowdhury SH, Fragomeni S, editors. Advances in mechanics of structures and materials. Lisse: Balkema; 2002. p. 119–24. [16] Chen WF, Han DJ. Plasticity for structural engineers. New York: Springer-Verlag; 1988.
1917
[17] Belytschko T. Plane stress shakedown analysis by finite elements. Int J Mech Sci 1972;14:619–25. [18] Groß-Weege J. On the numerical assessment of the safety factor of elastic–plastic structures under variable loading. Int J Mech Sci 1997;39:417–33. [19] Corradi L, Zavelani A. A linear programming approach to shakedown analysis of structures. Comput Methods Appl Mech Eng 1974;3:37–53. [20] Nguyen DH, Palgen L. Shakedown analysis by displacement method and equilibrium finite elements. Trans Canad Soc Mech Eng 1980;61:34–40. [21] Liu Y, Zhang X, Cen Z. Numerical determination of limit loads for three-dimensional structures using boundary element method. Eur J Mech A/Solids 2004;23:127–38. [22] Nagtegaal JC, Parks DM, Rice JR. On numerically accurate finite element solutions in the fully plastic range. Comput Methods Appl Mech Eng 1974;4:153–77. [23] Calladine CR. Plasticity for engineers. Chichester: Ellis Horwood; 1985. [24] Drud A. CONOPT – a large-scale GRG code. ORSA J Comput 1994;6:207–16.
A simple mixed finite element for static limit analysis X. Yu a, F. Tin-Loi a
b,*
School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, NSW 2006, Australia b School of Civil and Environmental Engineering, University of New South Wales, NSW 2052, Australia Received 5 October 2005; accepted 7 August 2006
Abstract In this paper, we investigate the behavior of a simple mixed finite element for the limit analysis of plane structures. In particular, its ability to overcome incompressibility locking in plane strain situations is investigated. The element is constructed from a piecewise constant displacement field and a piecewise bilinear stress field, and is used within a mathematical programming based discrete representation of the classical static formulation. Several benchmark examples of both plane stress and plane strain situations are solved to illustrate the predictive accuracy and to assess the large-scale capability of the element. The results are compared with those obtained by a recent sophisticated enhanced strain mixed element formulation. 2006 Elsevier Ltd. All rights reserved. Keywords: Limit analysis; Mixed finite element; Plasticity; Collapse loads
1. Introduction Limit analysis is without doubt one of the most widely used instances of the so-called ‘‘direct’’ or ‘‘simplified’’ approaches for nonlinear analyses. As is well-known, it provides useful information regarding the collapse load of a perfectly-plastic structure, while avoiding the computationally expensive effort of an evolutive analysis that follows the whole time history of loading. In recent years, there has been renewed interest in discretized limit analysis, mostly motivated by three aims: application to specific, possibly non-conventional engineering situations; overcoming the volumetric locking behavior encountered in plane strain and 3D problems; and development of efficient computational strategies for solving practically motivated and large size problems. A comprehensive survey of such developments is not possible. We mention only a representative sample of recent work such as those by: Capsoni and Corradi [1] and Capsoni [2] who devel-
*
Corresponding author. Tel.: +61 2 9385 6223; fax: +61 2 9385 3169. E-mail address: [email protected] (F. Tin-Loi).
0045-7949/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.08.019
oped locking-free mixed 4-node plane strain elements with low distortion sensitivity and coarse mesh accuracy; Christiansen and Andersen [3] who devised a new computational approach that exploits the static–kinematic duality property of limit analysis for solving, without locking, large plane strain problems; Christiansen and Pedersen [4] who subsequently improved that method by incorporating a strategy for adaptive mesh refinement; Francescato et al. [5] who developed a combined static and kinematic approach to study the ductile failure of porous metallic materials; Krabbenhoft and Damkilde [6] who used an interior point method for the solution of general nonlinear programming problems arising from a static approach; Lyamin and Sloan [7,8] who also adopted an interior point method for both lower and upper limit analyses using nonlinear programming; and Tin-Loi and Ngo [9] who investigated the robustness of the p-version finite element method within the static framework. It is also notable that most approaches, at variance with the present work, use the kinematic approach to limit analysis. The present paper is primarily concerned with the assessment of a very simple mixed finite element, mentioned briefly in the classical paper of Anderheggen and
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Kno¨pfel [10], as to its capability to overcome incompressibility locking, its suitability for large-scale limit analysis, and also its ability to provide good indications of expected, not necessarily unique, collapse mechanisms. The element is a quadrilateral with a constant displacement field and a bilinear stress field. In spite of its simplicity and potential to overcome locking (as indicated by the favorable ratio of 3/2 of number of stress parameters to number of displacement parameters [10]), this particular element has surprisingly received almost no attention with regards to limit analysis. Its performance in terms of robustness and efficiency is compared with that of the sophisticated mixed enhanced strain element of Capsoni and Corradi [11], which was also implemented for limit analysis in this study, albeit within a static (rather than kinematic [11]) framework. The organization of this paper is as follows. In Section 2, we review the limit analysis model in a form suitable for its direct implementation via the simple element to be investigated. Section 3 deals with some key features of constructing the mixed quadrilateral finite element proposed. A summary of the discretized static problem is then presented in Section 4, before we describe (Section 5) various numerical examples dealing with classical benchmark problems, involving both plane stress and plane strain situations obeying von Mises or Tresca criterion to illustrate the performance of the simple mixed element. We conclude in Section 6. 2. The static limit analysis problem The classical static approach to limit analysis essentially requires that the assumed stress field must be in equilibrium with external loads (governed by a single monotonically increasing load parameter l) and must satisfy the yield condition everywhere within the continuum. The computed collapse limit is then a lower bound to the true limit. In the following, we briefly present a description of the limit analysis model in a form suitable for discretization by mixed elements, including in particular the one which forms the focus of this paper. The derivation of the equilibrium conditions follows the principle of virtual work, in much the same fashion as used by Anderheggen and Kno¨pfel [10] and Christiansen [12]. Typical indices i and j are used to refer to the Cartesian coordinates directions (x, y). Consider a bounded domain V occupied by a rigid perfectly-plastic material continuum. The boundary C of that body consists of two regions Ct and Cr. On Ct, the surface forces lti are prescribed while Cr refers to the region on which the (unknown) reactive surface tractions ri are developed. Volume (body) forces lgi act in V. The limit analysis problem is simply: given a fixed load distribution (gi, ti), it is required to calculate the maximum multiplier for l of this load (collapse load multiplier) that the body can sustain without collapsing. Collapse occurs through a mechanism formation in which the plastic flow is unrestricted. It is also
1907
of interest to find out the flow (mechanism) and stress fields at collapse. The external work rate W of forces (lgi, lti, ri), associated with a virtual plastic flow ui, is Z Z Z ui lgi dV þ ui lti dC þ ui ri dC: ð1Þ W ¼ V
Ct
Cr
The internal work or dissipation rate D for the pair of stresses rij and flow rates is given by Z 1 ðui;j þ uj;i Þrij dV D¼ ð2aÞ V 2 Z Z ð2bÞ ¼ ui rij;j dV þ ui rij mj dC V
C
where a subscripted comma indicates partial derivative, and mj are the components of a unit outward normal vector to the boundary. Note that (2b) has been obtained from (2a) by applying Green’s formula (assuming that the stress is continuous within the domain V). Differentiability of u is then no longer required to calculate D. We can now write the familiar equation of virtual work rate as W ¼ D:
ð3Þ
This equation also represents the equilibrium conditions for the body, as can easily be seen by equating (1) and (2b) to give Z Z ðrij;j þ lgi Þui dV þ ðlti rij mj Þui dC V Ct Z þ ðri rij mj Þui dC ¼ 0: ð4Þ Cr
We further assume that r is continuous across boundaries of domain V so that rij mj ¼ lti on C t ; rij mj ¼ ri on C r :
ð5aÞ ð5bÞ
Using these relations, we can now simplify (4) to give Z ðrij;j þ lgi Þui dV ¼ 0 in V ð6Þ V
which will be used to construct our mixed finite element, underpinned by, as assumed, independent (continuous) stress and (discontinuous) displacement fields. Finally, the stress field rij must satisfy the yield conditions for the assumed material, or in generic form f ðrij Þ 6 0
ð7Þ
which represents a convex (possibly unbounded) set of admissible stresses. The foregoing now allows the maximum limit multiplier ^ to be calculated from the following optimization l problem:
1908
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
^ ¼ max l
9 > > > =
l R
ðrij;j þ lgi Þui dV ¼ 0 in V ; V > rij vj ¼ lti on C t ; > > ; f ðrij Þ 6 0:
subject to :
ð8Þ
uym are arbitrary, we obtain the following series of equations that are decoupled across elements: Z ðrij;j þ lgi Þ dV ¼ 0; m ¼ 1 . . . M: ð11Þ Vm
It should be noted that we have eliminated (5b) since there are no constraints on ri.
By using Green’s theorem, we can rewrite (11) as Z Z rij vj dC þ l gi dV ¼ 0; m ¼ 1 . . . M: Cm
3. The simple mixed quadrilateral finite element Assume, as is conventional, that the two-dimensional domain under consideration is suitably subdivided into a mesh of M quadrilateral elements and N nodes. A generic element is shown in Fig. 1. The displacement field u is assumed to be constant within each element but it does not necessarily satisfy inter-element continuity. In particular, for the mth element, we have ux u¼ ; ðx; yÞ 2 V m ; m ¼ 1 . . . M ð9Þ uy m where Vm is the domain of the mth element. For the stress representation, we assume a bilinear stress field over each element domain. Thus, 2
rxx
3
6 7 4 ryy 5 ¼ rxy
4 X k¼1
2
rxx
3
6 7 N k 4 ryy 5 rxy k
ð10Þ
where the stress parameters [rxx, ryy, rxy]k are defined at each of the four nodes k = 1 . . . 4 of the element (see Fig. 1); and Nk is a local (element-level) weight function, identical to the shape function used in a classical 4-node isoparametric element. Stress continuity within and across elements is thus guaranteed. 4. The discrete limit analysis model We are now in a position to generate the discrete counterpart for the limit analysis problem (8). As a first step in constructing the discretized body force equilibrium relations for the model corresponding to (6), we substitute (9) into (6) and, after noting that uxm and 3
ð12Þ
Vm
We now substitute (10), for each element, into (12) to generate the assembled body force equilibrium relations for the whole finite element model. That is, 9 8 2 3 rxx > > = 4 < X t X ys yr Px 0 xr xs 6 7 ¼0 ; 4 ryy 5 þ l > >2 Py 0 xr xs y s y r ; m : k¼1 rxy k m ¼ 1 . . . M;
ð13Þ
where t is the element thickness (assumed to be constant for each element); x and y used with subscripts r, k, s are coordinates, with r, k, s representing nodal numbers pertaining to the mth element and defined in an anticlockwise convention (i.e. k = 1, 2, 3, 4 corresponds to r = 4, 1, 2, 3 and s = 2, 3, 4, 1); and lPx and lPy are, respectively, the x and y components of element volume forces. The surface force equilibrium equations (or stress boundary conditions) given by (5a) can be discretized in a number of ways, including, for example, a weighted integration method, or a boundary collocation method. When the load distribution is a priori known, the weighted integration method is written as (14a). It can be extended to the more general form of (14b) when the distribution of load is unknown but its total magnitude, however, is known, and thus can be related to the limit load, as in some rigid punch problems. In particular, 82 3 2 39 < rxx rxx = t yb ya 0 xa xb 4 ryy 5 þ 4 ryy 5 ; 0 xa xb y b y a : 2 rxy a rxy b Z tx ¼ l dC; ð14aÞ ty X Ca!b T ½ axx ayy axy i ½ axx ayy axy i ¼ l; ð14bÞ i
where [tx, ty] is the vector of distributed (per unit area) surface forces whose component directions coincide with the coordinates; a and b are two nodes at an element edge that is located at a stress boundary, with a ! b following the
4
uy
q
ux
[σ xx , σyy , σxy ]n Node n
2
y
q
q
b
a
a
a
1
x Fig. 1. Definition sketch of simple quadrilateral finite element.
Fig. 2. Stress boundary conditions: (a) linear continuous distribution along straight boundary; (b) with sharp change of boundary normal; (c) with sharp change of stress magnitude.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
anticlockwise direction for the element (Fig. 2a); Ca!b is the surface from a to b; and ½ axx ayy axy i is a weight vector, the explicit specification of which is determined by individual problems. The boundary collocation method only applies to the case where the load distribution is known beforehand. With this method, the constraints are imposed directly on the stress parameters of the boundary nodes, or rxx rxy vx tx ¼l ; n 2 ½1; N on C ð15Þ rxy ryy n vy n ty where [vx, vy]n is the directional vector of the surface normal at node n. Special care is needed to implement (15) when the normal to the stress boundary is discontinuous (Fig. 2b) or the load distribution q changes sharply (Fig. 2c). For these two cases, two sets of relations (15) need to be written to represent the stress boundary when the boundary node is approached from two opposite directions. In the case that these sets of equations conflict, which is always the case for Fig. 2c and sometimes also for Fig. 2b, one needs to make an arbitrary decision to drop one set from the conflicting pair. Obviously, this will introduce errors. To reduce such errors, we use the local mesh refinement scheme shown in Fig. 3. In a series of preliminary studies, we have tested both the integration and collocation methods and obtained very similar results both in accuracy and performance. Unless otherwise stated, the results reported in this paper are based on the boundary collocation method represented by (15). The yield constraints, typically either the well-known Tresca or von Mises criterion, are enforced at the nodes of the mesh. Since the stress distributions are piecewise bilinear, the yield conditions are automatically satisfied at every point in the discretized domain.
n div = 0
ndiv = 1
n div = 2
q
Fig. 3. Local mesh refinement: (a) before refinement; (b) after refinement with ndiv = 1; (c) after refinement with ndiv = 2; (d) used for sharp change of stress magnitude.
1909
It is useful to recall at this stage these yield criteria for plane situations; r0 is used to indicate the yield stress in simple tension. The material constant in these yield conditions can also be the pffiffiffi shear strength s0 of the material; for von Mises r0 ¼ s0 3 and for Tresca r0 = 2s0. Subscripts M and T denote Mises and Tresca, respectively. For plane strain, 4 2 fM ¼ ðrxx ryy Þ þ 4r2xy r20 ; 3 2
fT ¼ ðrxx ryy Þ þ 4r2xy r20 :
ð16aÞ ð16bÞ
For plane stress, 2
fM ¼ ðrxx ryy Þ þ rxx ryy þ 3r2xy r20 ;
ð16cÞ
fT ¼ maxðjr1 j; jr2 j; jr1 r2 jÞ r0 ;
ð16dÞ
where r1 and r2 are principal stresses. Finally, the discrete limit analysis problem can be set up as the following nonlinear programming problem: 9 ^ ¼ max l l > > > > = subject to: element equilibrium ð13Þ; stress boundary conditions ð14Þ or ð15Þ; > > > > ; yield conditions ð16Þ: ð17Þ It is useful to note that, with the above approach, the requirement of static admissibility is only satisfied in a ^ can only be considweak form. Therefore, the computed l ered to be an approximation, and not a strict bound, to the exact value. It may be argued that a strict bound solution, based on equilibrium elements [7] or compatible displacement elements [8], is more valuable in a limit analysis of a large-scale practical problem. Unfortunately, strict bounds can only be established in idealized situations. Even with the equilibrium element [7], static admissibility is not fully guaranteed. For example, (a) boundary equilibrium can be violated in the case of nonlinear distributions of boundary stresses; (b) inter-element equilibrium can be violated in the case of sharp changes of shear stresses at a boundary; and (c) element equilibrium can be violated in the case of non-constant body force. In addition, the strict bound property is also lost if the material strength has a nonlinear distribution across the domain, or if the material does not satisfy the associated flow rule. In view of these and also because our approach typically converges, as the model is refined, to the true value of collapse multiplier, we do not aim to obtain strict bounds. Obviously, whether the convergence is from above or below will in any case reveal the nature of the computed collapse limit. Finally, we note that solution of (17) will not only pro^, but also an idea of the (posvide the collapse multiplier l sibly non-unique) collapse mechanism obtained from the
1910
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Lagrange multipliers associated with the element equilibrium conditions (13) at optimality. 5. Numerical examples Our pilot numerical implementation consists of the following features: (a) formation of relevant quantities to set up the nonlinear programming problem in MATLAB; (b) production of an appropriate data (text) file for use by a generic GAMS (General Algebraic Modeling System) [13] model that implements the nonlinear programming formulation; and (c) solution by one of the state-of-the-art GAMS solvers (we typically use the well-known nonlinear programming solver GAMS/MINOS [14]). The use of GAMS provides a large number of advantages of which the most important for the work described herein are the availability of a number of industry standard mathematical programming solvers for linear and nonlinear programming, simplicity of modeling, large-scale facility, and automatic differentiation capability. In the following, we summarize some results related to four problems, namely a square plate with central circular hole under biaxial tension, a two-notched tensile specimen, a thick clamped beam, and the well-known Prandtl punch problem. For comparison, a sophisticated quadrilateral mixed element, initially developed for elastic–plastic analysis by Capsoni and Corradi [11], was also implemented within our static approach to limit analysis. Our implementation of this element involves stress parameters defined as unknown variables, and the yield conditions imposed at four vertices for each element; details are given in [15]. In terms of stress and velocity fields, this implementation is in effect the ‘‘dual’’ of that of Capsoni and Corradi [1] who, it must be stressed, used a kinematic approach. However, because of the differences in imposing yield conditions, the two implementations are not equivalent. For instance, Capsoni and Corradi [1] do not allow an element to become partly-plastic (and partly rigid) since no plasticrigid boundary was accommodated within their element. In addition, their iterative algorithm is to progressively add elements, deemed to be non-critical, to a ‘‘rigid’’ set. However, once included in that set, no element can become a candidate for future yielding. These limitations do not apply to our present approach. Therefore, given the same mesh and stopping criterion, the two implementations should yield similar, but not necessarily the same, predictions. In the following examples, unless otherwise noted, we compare the performance of our simple mixed element and with that of the mixed element of Capsoni and Corradi, as implemented by us within a static framework. Of course, we used identical testing conditions, e.g. same mesh, same solver, same feasibility and optimality tolerances (typically 106, the default in MINOS [14]), and same computer. All runs reported were carried out on a Pentium III 800 MHz PC.
5.1. Square plate with central circular hole This first example is a plane stress problem, involving the well-known reference case of a square plate with a hole diameter to plate width ratio of 0.2 (Fig. 4a). We analysed a quarter of the plate modeled with various meshes defined by parameters nh and nl, as in Fig. 4b. When p2 = 0, the problem becomes a uniaxial tension problem for which the exact collapse solution for the plane stress Tresca condition is p1/r0 = 0.8 [16] with the tensile stress at the plate edge arbitrarily distributed so that p1 refers to the averaged tensile stress at the plate edge. Three feasible failure mechanisms, namely in-plane symmetrical shearing, in-plane unsymmetrical shearing and out-ofplane shearing, were reported, any of which leads to the same solution. In our case, in view of the quarter model used, only a symmetrical failure mechanism is possible. We considered four cases as follows: • • • •
Case A: von Mises, p2 = p1, uniform stress. Case B: von Mises, p2 = 0.5p1, uniform stress. Case C: von Mises, p2 = 0, uniform stress. Case D: Tresca, p2 = 0, arbitrary distribution of tensile stress p1.
It is noted that the stress boundary for Case D was implemented based on (14b) because the load distribution is a priori unknown; for other cases, (15) was applied. Tables 1–4 summarize the results of these analyses for different numbers of nodes N and elements M. Some statistics related to the nonlinear programming problem (namely, number of variables nvar, and number of constraints ncon) and an idea of the computational cost (namely, run time in seconds) are also given. Table 5 also compares our results for the von Mises cases with some reported ones in the literature [9,16–21]. The results given in these Tables 1–5 show that: the calculated collapse pressures all converge with mesh refinement; our simple mixed element performed far better in terms of robustness and efficiency than our implementation
p2
p2
p1
p1
p1
nh
p2
nl
1.5 × n l
Fig. 4. Square plate: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 Table 1 Square plate: Case A (von Mises, p2 = p1) nh
nl
N
M
nvar
ncon
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
144 459 1611 6003
1911
Table 4 Square plate: Case D (Tresca, p2 = 0, arbitrary stress distribution) p1/r0
CPU (s)
nh
0.8486 0.8820 0.8906 0.8933
0.2 1.7 37.6 1491.1
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.9003 1.2 8 4 141 116 1973 2116 0.8964 25.0 16 8 513 464 7889 8408 0.8955 961.2 32 16 1953 1856 31,553 33,520 Faila –
nl
N
M
nvar
ncon
p1/r0
139 449 1591 5963
0.7762 0.7944 0.7984 0.7990
CPU (s) 0.2 1.9 23.7 618.1
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 500 537 0.8134 0.9 8 4 141 116 1984 2117 0.8042 9.2 16 8 513 464 7910 8409 0.8016 213.8 32 16 1953 1856 31,594 33,521 0.8006 6485.5
a
Run stopped after about 40 h of CPU time without convergence being achieved. Table 5 Square plate: comparison of limit multipliers Table 2 Square plate: Case B (von Mises, p2 = 0.5p1) nh
nl
N
M
nvar
ncon
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
144 459 1611 6003
p1/r0 0.8810 0.9007 0.9065 0.9086
CPU (s) 0.3 2.4 29.7 758.6
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.9254 1.1 8 4 141 116 1973 2116 0.9152 18.7 16 8 513 464 7889 8408 0.9121 370.2 32 16 1953 1856 31,553 33,520 0.9110 36077.2
Table 3 Square plate: Case C (von Mises, p2 = 0) nh
nl
N
M
nvar
Current simple mixed element 4 2 42 29 127 8 4 141 116 424 16 8 513 464 1540 32 16 1953 1856 5860
Reference
Case A
Case B
Case C
Present work, fine mesh Analytical [16] Tin-Loi and Ngo [9] [p-version finite element] Belytschko [17] [triangular equilibrium element] Groß-Weege [18] [quadrilateral equilibrium element] Corradi and Zavelani [19] [triangular constant-stress element] Nguyen and Palgen [20] [quadrilateral equilibrium element] Liu et al. [21] [boundary element]
0.893 – 0.895
0.909 – 0.912
0.799 0.800 0.803
–
–
0.780
0.882
0.891
0.782
0.767
–
0.691
0.704
–
0.564
0.903
0.915
0.795
5.2. Two-notched tensile specimen
ncon 144 459 1611 6003
p1/r0 0.7640 0.7958 0.7984 0.7991
CPU (s) 0.2 2.4 34.9 969.0
Mixed element of Capsoni and Corradi (present implementation) 4 2 42 29 494 536 0.8147 1.2 8 4 141 116 1973 2116 0.8057 20.4 16 8 513 464 7889 8408 0.8022 407.3 32 16 1953 1856 31,553 33,520 0.8008 16500.3
of the mixed element of Capsoni and Corradi [11], which is expected as the size of the nonlinear programming problem increases far more rapidly in the latter case; and our predictions for Cases C and D match particularly well the only analytical solution we are aware of. Finally, plots of the failure mechanisms we obtained for each of these four cases are shown in Fig. 5. The differences in the failure mechanisms between Case C and Case D are possibly due to the different yield conditions. For Case D, our failure mechanism agrees well with that of the analytical solution [16].
This type of problem, involving a plane strain doublenotched rectangular specimen under tension (Fig. 6a), is a popular benchmark used for testing volumetric locking in elastoplastic analyses [22]. We chose to carry out limit analyses on ‘‘infinitely’’ long specimens obeying von Mises criterion. Finite element models of a quarter of the plate (Fig. 6b) were constructed for the following three cases: • w/b = 3, h/b = 10. • w/b = 5, h/b = 10. • w/b = 8, h/b = 15. The indicated h/b dimensions were selected after some preliminary sensitivity analyses. Further, as indicated in Fig. 6b, we introduced some local mesh refinement, in the fashion of Fig. 3d, at the notch tip in an attempt to accommodate any stress discontinuities. In order to assess how extensive this refinement should be, we first carried out a series of analyses to investigate the influence of this local mesh refinement (measured through parameter ndiv) on the w/b = 3 case and using our proposed simple mixed element. Table 6 displays the results obtained; nb refers to
1912
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Case A
Case B
Case C
Case D
Fig. 5. Square plate: failure mechanisms for four load cases.
p
Table 6 Notch mesh refinement (w/b = 3, h/b = 10, nb = 10)
nb ( = 2)
b
2b
h
2h
ndiv
N
M
nvar
ncon
p/s0b
CPU (s)
0 1 2 3 4 5 6
3131 3136 3141 3146 3151 3156 3161
3000 3004 3008 3012 3016 3020 3024
9394 9409 9424 9439 9454 9469 9484
9542 9558 9574 9590 9606 9622 9638
3.5082 3.5208 3.5461 3.5506 3.5569 3.5582 3.5595
2416.0 2910.0 1925.6 1623.2 1590.8 1407.5 1439.4
2w
p
p
Fig. 6. Double-notched specimen: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
the number of elements, excluding those in the local refinement zone, underneath the half-width b of the unnotched segment (e.g. nb = 2 in the mesh shown in Fig. 6b). Interestingly, the computing times dropped dramatically with appropriate local refinement. The predicted collapse load factor varied less than about 0.1% when the local refinement was changed from ndiv = 4 to ndiv = 5. We adopted the latter discretization for our subsequent main runs.
Tables 7–9 summarize the limit analyses for the three cases. It is evident that our simple element performed reasonably well in providing collapse limits in all instances, whilst the models constructed using the mixed element of Capsoni and Corradi [11] could not be solved for even low values of nb. Moreover, our simple element did not exhibit any locking for this plane strain situation. We were unable to assess properly the accuracy of our computed limit loads since we could not find any corresponding values reported in the literature. However, the long double-notched tension specimen is equivalent to a rigid punch over a narrow block for which analytical lower and upper bounds are available [23]. Table 10 provides the appropriate comparisons. Plots of the failure mechanisms revealed by our simple finite element based model are given in Fig. 7 (with the
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
1913
Table 7 Double-notched specimen: w/b = 3, h/b = 10, ndiv = 5 nb
N
M
nvar
Current simple mixed element 2 172 140 517 4 558 500 1675 6 1184 1100 3553 8 2050 1940 6151 10 3156 3020 9469 12 4502 4340 13,507 Mixed 2 4 6
ncon
pw/s0b
CPU (s)
550 1738 3646 6274 9622 13,690
4.3495 3.8293 3.6754 3.6016 3.5582 3.5297
2.2 22.5 120.2 498.4 1407.5 3920.5
element of Capsoni and Corradi (present 172 140 2381 2555 558 500 8501 9065 1184 1100 18,701 19,895
implementation) 3.3758 21.9 3.3892 187.5 Fail –
Table 8 Double-notched specimen: w/b = 5, h/b = 10, ndiv = 5 nb
N
M
nvar
Current simple mixed element 2 256 220 769 4 886 820 2659 6 1916 1820 5749 8 3346 3220 10,039 10 5176 5020 15,529 12 7406 7220 22,219 Mixed 2 4 6
ncon
pw/s0b
CPU (s)
810 2738 5866 10,194 15,722 22,450
4.8701 4.5206 4.4109 4.3574 4.3256 4.3046
5.9 93.6 651.2 3392.9 14897.8 38769.1
element of Capsoni and Corradi (present 256 220 3741 4003 886 820 13,941 14,841 1916 1820 30,941 32,879
implementation) 4.1865 52.0 4.1957 7123.1 Faila –
a Run stopped after about 40 h of CPU time without convergence being achieved.
Fig. 7. Double-notched specimen: failure mechanisms for three geometries (notch is located at bottom left hand corner).
Table 9 Double-notched specimen: w/b = 8, h/b = 15, ndiv = 5 nb
N
Current simple 2 552 4 2038 6 4484 8 7890
M
nvar
mixed element 500 1657 1940 6115 4340 13,453 7700 23,671
ncon
pw/s0b
CPU (s)
5.3. Thick clamped beam
1720 6238 13,636 23,914
5.4248 5.2198 5.1531 5.1191
32.1 778.5 6271.0 34631.0
This problem, investigated recently by Capsoni and Corradi [1], involves a thick clamped beam in plane strain under a uniformly distributed pressure loading applied at its top central part, as shown in Fig. 8a. We modeled half of this beam (Fig. 8b) with mesh refinement (ndiv = 5) provided at the top pressure discontinuity point and the two corners of the fixed ends. We adopted a von Mises yield condition.
Mixed element of Capsoni and Corradi (present implementation) 2 552 500 8501 9065 4.9811 950.2 4 2038 1940 32,981 35,045 Faila – a Run stopped after about 40 h of CPU time without convergence being achieved.
notch oriented at the left hand bottom corner of the models).
p
p
2b h
Table 10 Double-notched specimen: comparison of limit multipliers
Present work Upper bound [23] Lower bound [23]
w/b = 3
w/b = 5
w/b = 8
3.530 4.000 3.000
4.305 4.899 3.333
5.119 6.000 3.556
2w Fig. 8. Thick clamped beam: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
1914
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
Three geometries were analysed:
Table 13 Thick clamped beam: w/b = 5, h/b = 2, ndiv = 5
• w/b = 1, h/b = 2. • w/b = 2, h/b = 2. • w/b = 5, h/b = 2. The w/b = 1 case is particularly noteworthy since the failure mechanism associated with the exact solution requires shear sliding at the fixed end, and only models capable of displacement discontinuities will be able to simulate this failure mode. Tables 11–13 summarize some of the results obtained; np refers to the number of elements, excluding local refinement, in width b of material (i.e. under load p). Once again, it can be conclusively seen that our simple element is superior, both in terms of robustness and efficiency, as compared with the results provided by our static implementation of the more refined element of Capsoni and Corradi [1]. Contrary to expectation, however, the values of collapse load multiplier do not converge strictly monotonically for the w/b = 5 case (Table 13) as the mesh is refined.
Table 11 Thick clamped beam: w/b = 1, h/b = 2, ndiv = 5 np
N
nvar
ncon
p/s0
CPU (s)
Current 2 4 6 8 10 12
simple mixed element 45 28 136 75 52 226 121 92 364 183 148 550 261 220 784 355 308 1066
136 226 364 550 784 1066
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000
0.3 0.6 1.7 3.3 7.6 13.1
M
N
np
M
nvar
ncon
Current simple mixed element 2 110 80 331 4 244 200 733 6 458 400 1375 8 752 680 2257 10 1126 1040 3379 12 1580 1480 4741 Mixed 2 4 6 8 10 12
element of Capsoni 110 80 244 200 458 400 752 680 1126 1040 1580 1480
p/s0
357 775 1433 2331 3469 4847
0.9461 1.0057 1.0142 1.0160 1.0162 1.0160
and Corradi (present 1361 1465 3401 3641 6801 7257 11,561 12,313 17,681 18,809 25,161 26,745
CPU (s) 1.5 8.1 35.5 98.6 242.3 587.4
implementation) 1.0961 9.7 1.0419 44.1 1.0307 173.2 1.0257 678.5 1.0226 3172.7 1.0207 10430.4
The failure mechanisms we obtained are depicted in Fig. 9. It is notable that the shear failure for the w/b = 1 case has been captured properly. The transition of shear failure to bending failure is also clear as w/b increases. Finally, Table 14 compares the present predictions with available [1] analytical bound solutions and computed results. As can be seen, our element predicts the exact solution for w/b = 1. It is noted that, in Table 14, the results attributed to the mixed element of Capsoni and Corradi were obtained by
Mixed element of Capsoni and Corradi (present implementation) 2 45 28 477 503 2.3992 1.5 4 75 52 885 935 Faila – a MINOS exited with ‘‘too many iterations’’ without achieving an optimal solution.
Table 12 Thick clamped beam: w/b = 2, h/b = 2, ndiv = 5
w/b = 1
np
N
nvar
ncon
p/s0
CPU (s)
Current 2 4 6 8 10 12
simple mixed element 80 56 241 136 104 409 224 184 673 344 296 1033 496 440 1489 680 616 2041
255 427 695 1059 1519 2075
1.5118 1.7264 1.7341 1.7395 1.7419 1.7421
0.4 2.8 7.3 21.4 40.8 97.8
M
Mixed element of Capsoni and Corradi (present 2 80 56 953 1021 4 136 104 1769 1889 6 224 184 3129 3333 a
implementation) 1.9604 5.4 1.8319 86.5 Faila –
MINOS exited with ‘‘too many iterations’’ without achieving an optimal solution.
w/b = 2
w/b = 5 Fig. 9. Thick clamped beam: failure mechanisms for three geometries.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 Table 14 Thick clamped beam: comparison of limit multipliers
Present work Upper bound [1] Lower bound [1] Capsoni and Corradi [1]
1915
Table 15 Prandtl’s punch: uniform distribution of vertical contact forces N
M
nvar
w/b = 1
w/b = 2
w/b = 5
np
2.000 2.000 2.000 2.286
1.742 2.000 1.333 1.828
1.016 2.000 0.500 1.071
Current simple mixed element 2 55 40 166 4 189 160 568 6 403 360 1210 8 697 640 2092 10 1071 1000 3214
them through their kinematic approach and are reported in [1], while the results attributed to the mixed element of Capsoni and Corradi in Tables 11–13 are from our present static implementation. Obviously, as is indeed the case, our results do not match those of Corradi and Capsoni since we used, as indicated earlier, a static approach (instead of a kinematic approach), different mesh configurations, a different method for ensuring yield conformity, and possibly also different convergence criteria. 5.4. Prandtl’s punch This classical plane strain problem consists of a semiinfinite rigid plastic medium under a punch load of width 2b (Fig. 10a). The exact collapse load [16], based on the plane strain Tresca criterion, is P = (2 + p)s0b = 5.1416s0b. We constructed a finite element model (Fig. 10b) involving half of a block of medium with dimensions 2w · h. The interaction between the isolated block and the rest of the medium was represented by fixed (no x or y displacement) boundary conditions. Initial trial simulations showed that w/b = 5 and h/b = 2 were large enough to simulate an ‘‘infinite’’ medium. Uniform meshes were used and local mesh refinement was not applied since the collapse load factor is independent of w/b. A plane strain Tresca criterion was used in all analyses. The interactions between the punch and the medium were represented by vertical force distributions and horizontal displacement restraints applied on the block. We considered two cases of vertical force distributions as follows:
ncon 159 555 1191 2067 3183
P/s0b
CPU (s)
5.1001 5.1170 5.1238 5.1277 5.1300
0.4 3.6 31.0 60.7 133.3
Mixed element of Capsoni and Corradi (present implementation) 2 55 40 681 714 5.3494 1.1 4 189 160 2721 2868 5.2485 17.0 6 403 360 6121 6462 5.2158 74.9 8 697 640 10,881 11,496 Faila – a
MINOS exited reporting ‘‘problem is unbounded or badly scaled’’.
Table 16 Prandtl’s punch: arbitrary distribution of vertical contact forces np
N
nvar
ncon
P/s0b
CPU (s)
Current 2 4 6 8 10
simple mixed element 55 40 166 189 160 568 403 360 1210 697 640 2092 1071 1000 3214
157 551 1185 2059 3173
5.1019 5.1194 5.1262 5.1298 5.1319
0.4 3.5 18.6 83.0 135.0
M
Mixed element of Capsoni and Corradi (present implementation) 2 55 40 684 715 5.8490 1.3 4 189 160 2726 2869 5.5069 17.2 6 403 360 6128 6463 Faila – a
MINOS exited reporting ‘‘problem is unbounded or badly scaled’’.
It is noted that the latter arbitrary distribution case is a more appropriate representation for the rigid punch problem but we also considered the uniform distribution as it is one that has been frequently quoted in the literature.
• Uniform distribution (using (15) for the stress boundary). • Arbitrary distribution (using (14b) for the stress boundary).
2P
2b
P
b
h
2w
Fig. 10. Prandtl’s punch: (a) geometry and loading; (b) finite element mesh of quadrilateral elements.
Fig. 11. Prandtl’s punch: (a) failure mechanism for uniform pressure distribution; (b) failure mechanism for arbitrary pressure distribution.
1916
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917
The results of our limit analyses are summarized in Tables 15 and 16, representing, respectively, the uniform force distribution and arbitrary force distribution. Symbol np refers to the number of elements in contact with half of the punch. Corresponding collapse mechanisms are shown in Fig. 11. As can be seen, our simple mixed element performed well, and there does not appear to be any significant differences between the two load cases with both obtained mechanisms being similar. In line with our expectation, the collapse load predictions are both close to the theoretical value. 6. Concluding remarks Within the framework of a finite element based static approach to limit analysis, we present a simple 4-noded quadrilateral mixed finite element for reliably calculating collapse load multipliers. In particular, it is found that locking caused by volumetric incompressibility does not occur in plane strain problems. Furthermore, accurate solutions can be obtained in reasonable computing time, as demonstrated in the present paper through the results of four benchmark problems. Comparison of its performance with a sophisticated and recently developed [11] enhanced strain element favors our simple element. This is possibly not so surprising since, in our view, the element of Capsoni and Corradi is more suited for the solution of elastoplastic [11] rather than limit analysis [1] problems. Models using the enhanced strain element lead to very large increases in the sizes of the underlying nonlinear programming problems with mesh refinement. Capsoni and Corradi [1], it must be noted, used a kinematic approach and successfully reduced the problem size of their unconstrained but non-smooth optimization problems by identifying possible rigid regions which are then excluded in their optimization phase. Our prototype implementation consists essentially of a MATLAB code to provide appropriate equilibrium equations that a simple GAMS nonlinear programming model would use. As for the solution of the optimization problems, we used the state-of-the-art nonlinear solver GAMS/MINOS (which we found to be faster than another industrial-strength solver GAMS/CONOPT [24]). However, in some instances, we noted that the CPU time still grew very rapidly with increase in problem size. Adoption of the recently developed interior point nonlinear programming solvers is likely to provide an increase in computation speed, as they are based on quasi-Newton methods for which the iteration count is largely independent of the problem size, as has been reported [6–8]. Nevertheless, in the present study, we did not attempt to develop or implement special schemes to solve the optimization problems more efficiently as our primary concern was simplicity of finite element combined with a locking-free behavior. The use of mixed elements does not guarantee, as is well known, a strict bound to the collapse limit. Moreover, the estimates of limit multipliers typically converge mono-
tonically as the mesh is refined, thus providing a clear indication as to whether the actual collapse is being approached from above or below. The present simple mixed element appears, in addition to its acceptable accuracy, eminently suitable for providing good collapse mechanisms, even those involving displacement discontinuities. We did not notice any hourglass type displacement modes, as is often experienced with a 4-node isoparametric element that uses reduced (single centroid point) Gauss integration. Moreover, the stress field provided by the simple element can exhibit hourglass modes, as when the boundaries, subjected to rigid movement or fixed, have unknown stress distributions. However, the stress fields invariably become smooth in the yield region. Finally, worthwhile extensions of the present work would be develop a 3D version of the element, and to investigate adoption of piecewise linearized yield conditions, perhaps in concert with some adaptive scheme, so that a linear program needs to be solved rather than the typically more computationally expensive nonlinear optimization. Of course, the trade-off between increases in problem size with use of the more efficient linear programming codes needs to be carefully assessed. Acknowledgements This research was supported by a grant from the Australian Research Council (ARC). The work reported in this paper was carried out when Dr. Yu, the first named author, was employed, through ARC funds, as a Senior Research Associate at the University of New South Wales. Finally, the very constructive comments of two anonymous referees on an earlier version of this manuscript are gratefully acknowledged. References [1] Capsoni A, Corradi L. A finite element formulation of the rigidplastic limit analysis problem. Int J Numer Methods Eng 1997;40:2063–86. [2] Capsoni A. A mixed finite element model for plane strain limit analysis computations. Commun Numer Methods Eng 1999;15:101–12. [3] Christiansen E, Andersen KD. Computation of collapse states with von Mises type yield condition. Int J Numer Methods Eng 1999;46:1185–202. [4] Christiansen E, Pedersen OS. Automatic mesh refinement in limit analysis. Int J Numer Methods Eng 2001;50:1331–46. [5] Francescato P, Pastor J, Thai TH. Etude du crite`re de plasticite´ des mate´riaux poreux. C R Acad Sci, Se´r IIB 2001;329:753–60. [6] Krabbenhoft K, Damkilde L. A general non-linear optimization algorithm for lower bound limit analysis. Int J Numer Methods Eng 2003;56:165–84. [7] Lyamin AV, Sloan SW. Lower bound limit analysis using non-linear programming. Int J Numer Methods Eng 2002;55:573–611. [8] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Methods Geomech 2002;26:181–216. [9] Tin-Loi F, Ngo NS. Performance of the p-version finite element method for limit analysis. Int J Mech Sci 2003;45:1149–66.
X. Yu, F. Tin-Loi / Computers and Structures 84 (2006) 1906–1917 [10] Anderheggen E, Kno¨pfel H. Finite element limit analysis using linear programming. Int J Solids Struct 1972;8:1413–31. [11] Capsoni A, Corradi L. A mixed finite element model for plane strain elastic–plastic analysis: Part I – Formulation and assessment of the overall behavior; Part II – Application to the 4-node bilinear element. Comput Methods Appl Mech Eng 1997;141:67–93. [12] Christiansen E. Limit analysis of collapse states. In: Ciarlet PG, Lions JL, editors. Handbook of numerical analysis, vol. 4. Amsterdam: North-Holland; 1996. p. 193–312. [13] Brooke A, Kendrick D, Meeraus A, Raman R. GAMS: a user’s guide. Washington, DC 20007: Gams Development Corporation; 1998. [14] Gill PE, Murray W, Saunders MA, Wright MH. MINOS 5.5 user’s guide. Report 83-20R, Department of Operations Research, Stanford University, Stanford, CA, revised 1998. [15] Yu X, Tin-Loi F. Behavior of four finite elements for direct limit analysis. In: Loo YC, Chowdhury SH, Fragomeni S, editors. Advances in mechanics of structures and materials. Lisse: Balkema; 2002. p. 119–24. [16] Chen WF, Han DJ. Plasticity for structural engineers. New York: Springer-Verlag; 1988.
1917
[17] Belytschko T. Plane stress shakedown analysis by finite elements. Int J Mech Sci 1972;14:619–25. [18] Groß-Weege J. On the numerical assessment of the safety factor of elastic–plastic structures under variable loading. Int J Mech Sci 1997;39:417–33. [19] Corradi L, Zavelani A. A linear programming approach to shakedown analysis of structures. Comput Methods Appl Mech Eng 1974;3:37–53. [20] Nguyen DH, Palgen L. Shakedown analysis by displacement method and equilibrium finite elements. Trans Canad Soc Mech Eng 1980;61:34–40. [21] Liu Y, Zhang X, Cen Z. Numerical determination of limit loads for three-dimensional structures using boundary element method. Eur J Mech A/Solids 2004;23:127–38. [22] Nagtegaal JC, Parks DM, Rice JR. On numerically accurate finite element solutions in the fully plastic range. Comput Methods Appl Mech Eng 1974;4:153–77. [23] Calladine CR. Plasticity for engineers. Chichester: Ellis Horwood; 1985. [24] Drud A. CONOPT – a large-scale GRG code. ORSA J Comput 1994;6:207–16.