- Email: [email protected]
Limit analysis for a general class of yield conditions
Eur. J. Mech. A/Solids 19 (2000) 401–421 2000 Éditions scientifiques et médicales Elsevier SAS. All rights reserved S0997-7538(00)00170-4/FLA
Limit analysis for a general class of yield conditions Alan R.S. Ponter a,1 , Paolo Fuschi b , Markus Engelhardt a a Department of Engineering, University of Leicester, Leicester, LE1 7RH, UK b DASTEC - Departimento Arte Scienza Tecnica del Costruire, University of Reggio Calabria, Via Melissari Feo di Vito, 89124
Reggio Calabria, Italy Abstract – The paper describes a generalisation of the programming method described by Ponter and Carter (1997) for the evaluation of optimal upper bounds on the limit load of a body composed of a rigid/perfectly plastic material. The method is based upon similar principles to the ‘Elastic Compensation’ method which has been used for design calculations for some years but re-interpreted as a non-linear programming method. A sufficient condition for convergence is derived which relates properties of the yield surface to those of the linear solutions solved at each iteration. The method is demonstrated through an application to a Drucker–Prager yield condition in terms of the Von Mises effective stress and the hydrostatic pressure. Implementation is shown to be possible using the user routines in a commercial finite element code, ABAQUS. The examples chosen indicate that stable convergent solutions may be obtained. There are, however, limits to the application of the method if isotropic linear solutions are used for an isotropic yield surface. In an accompanying paper (Ponter and Engelhardt, 2000) the method is extended to shakedown and related problems. 2000 Éditions scientifiques et médicales Elsevier SAS plasticity / limit / programming / finite elements / Drucker–Prager / life assessment / design
1. Introduction There has been a growth of interest in recent times in methods of directly computing performance indicators for structures subjected to complex loading as an aid to design and for the implementaion of life assessment methods (Goodall et al., 1991). The simplest of these indicators is the limit load for a rigid/perfectly plastic body as the first step towards shakedown limits and related solutions. In previous papers (Ponter and Carter, 1997a, b) methods were described for the computation of optimal upper bounds on the limit state and shakedown limit of a perfectly plastic body, with a Von Mises yield condition, subjected to variable load and temperature. The purpose of this paper is to provide a generalisation of this method to a general yield condition. The general characteristics of the method are the same as the ‘Elastic Compensation’ method which has been described by a number of authors (see Marriott, 1985; Mackenzie and Boyle, 1993; Nadarajah et al., 1993; Shi et al., 1993; Mackenzie et al., 1993; Seshadri and Fernando, 1991). A sequence of linear problems with spatial varying moduli are generated where, at each iteration, the stress is brought within yield by changes of the moduli, assuming a fixed strain rate distribution. Interpreting the method as an upper bound programming technique, Ponter and Carter (1997a and b) demonstrated that the method provides a monotonically reducing sequence of upper bounds which converges to the least upper bound for a class of displacement fields. When implemented in a finite element code, the method proves to be very stable and converged solutions can usually be obtained after about 10 iterations. The lower bound derived from the resulting stress distributions, which formed the primary interest of the ‘Elastic Compensation’ method, were shown to converge, not necessarily monotonically, to the least upper bound. Other programming methods, both linear and non-linear, have been described in the literature of which Zhang et al. (1993) and Liu et al. (1995) provide recent examples. 1 E-mail: [email protected]
402
A.R.S. Ponter et al.
In this paper we discuss the generalisation of the method to an arbitrary convex yield surface. The arguments are confined to limit analysis and the extension to shakedown analysis and related cyclic problems will be discussed elsewhere (Ponter and Engelhardt, 2000). The generalisation of the method may be understood in two stages. The first involves a re-interpretation of the basis of the method as a general programming technique. At each iteration the moduli of the linear material are adjusted so that, for the current strain rate distribution, the linear material is matched to the yield condition. However, this condition alone is insufficient to uniquely define the moduli or for convergence to occur and additional conditions need to be applied. In the following sections we derive a sufficient condition for convergence and then interpret this condition in terms of an algorithm for the evaluation of the moduli of the linear material. Essentially, the surface of complementary energy density for the linear problem which touches the yield surface at the position of the current strain rate must otherwise lie outside the yield surface in stress space. The method is illustrated with reference to two classes of isotropic yield conditions which depend upon the Von Mises effective stress σ¯ and the hydrostatic pressure p. The first is for yield conditions which describe an ellipse in (σ¯ , p) space where the yield conditions and contours of constant complementary energy density may coincide and where the Von Mises yield condition is a special case. We find, as expected, that the method converges when the sufficient condition applies but also that convergence does not occur for all those cases studied where the condition does not apply. The method is then applied to two cases of a Drucker–Prager model for granular material where the application of the method in its general case can be demonstrated. We find that where the yield surface has regions of low curvature an isotropic yield condition may require the definition of an anisotropic linear material. All the numerical examples given in the paper have been generated using the user routines of a general purpose commercial finite element code, ABAQUS. These methods may, therefore, be viewed as general purpose methods which are capable of having significant use in design and where very accurate values of significant structural performance indicators are required.
2. Convexity conditions 2.1. The linear material We consider a class of linear materials defined by the strain rate potential U (˙εij )given by,
U (˙εij ) = 12 Dij kl ε˙ ij − r˙ij ε˙ kl − r˙kl ,
(1)
where Dij kl is the compliance tensor of coefficients which satisfy the usual symmetry properties and r˙ij is a predefined initial strain rate. The stress corresponding to a strain rate is given by, σijL =
∂U = Dij kl ε˙ kl − r˙kl . ∂ ε˙ ij
(2)
The complementary potential U (σij ) is given by, U (σij ) = 12 Cij kl (σij − qij )(σkl − qkl ),
(3)
Limit analysis for a general class of yield conditions
403
where Cij kl is the inverse of Dij kl and qij is an initial stress given by, r˙ij = −Cij kl qkl ,
qij = −Dij kl r˙kl .
or equivalently
(4)
and hence, ε˙ ij =
∂U L L = Cij kl σkl − qkl . ∂σij
(5)
The two potential functions are related by,
U σijL + U ε˙ ij = σijL ε˙ ij − qij r˙ij .
(6)
Above and in the following we denote stresses derived from such a linear material with a superscript L, to distinguish them from stresses associated with the yield condition which will be denoted by a superscript p. At this stage we note that, for uniqueness of the solution of continuum problems associated with the linear material, we require that both U and U are convex, i.e. U
ε˙ ij1
−U
ε˙ ij2
∂U − ∂ ε˙ ij
2 εij
ε˙ ij1 − ε˙ ij2 > 0
(7a)
with equality only when ε˙ ij1 = ε˙ ij2 , and
U σij1 − U σij2 −
∂U ∂σij
2 εij
σij1 − σij2 > 0
(7b)
with equality only when σij1 = σij2 . 2.2. The yield condition Consider a rigid plastic material with a yield condition, f (σij ) 6 0
(8)
and associated flow rule, ∂f , f (σij ) = 0. ∂σij We assume that the yield function is convex and that the maximum work principle holds, ε˙ ij = λ˙ p
σij − σij∗ ε˙ ij > 0, p
p
(9)
(10)
where σij is the stress associated with ε˙ ij at yield and σij∗ is any state of stress which satisfies the yield condition (8). p
p
2.3. The relationship between the linear material and the associated flow rule In (Ponter and Carter, 1997a), the convergence proof was constructed by assuming that the appropriate class of linear materials for a Von Mises yield condition was incompressible and isotropic. In the general case we
404
A.R.S. Ponter et al.
Figure 1. The relative convexity condition, equation (11), for a Von Mises yield condition requires that the shaded areas should have a positive area as shown.
require a method of constructing linear solutions so that convergence is assured. The key to this choice can be understood from a simple example given by Ponter and Carter (1997a), where we show that the iterative process can be understood as a programming method where the slope of the potential energy of the linear material (in the rate form) is matched to the slope of the upper bound to the limit load at each strain rate distribution. The solution of the linear problem produces the minimum of the potential energy and, at the same time reduces the upper bound to the limit load. In the following we describe a generalisation of this argument, in two stages. Initially we define a relationship between the convexity conditions (7) and the maximum work principle (10) which places constraints on the class of linear materials. We then show that this relationship is sufficient for convergence to take place in an iterative process for the limit load. 2.4. The relative convexity condition Consider a Von Mises yield condition f = σ − σy = 0 where σ is the Von Mises effective stress and σy is the uniaxial yield stress. In this case the linear material is an incompressible isotropic material described by a single material constant E and ε˙ = σ /E, where ε˙ denotes the effective strain rate. The matching process can be described graphically, as shown in figure 1. Suppose, at some stage in the iterative process, the current strain rate distribution is given by ε˙ i . The linear material may be matched to the yield condition by choosing E = σy /˙εi so that both the linear material and the associated flow rule yield the same state of effective stress for effective strain rate ε˙ i . The solution of the subsequent linear problem for this choice of E will produce a final strain rate ε˙ f which will generally differ from ε˙ i and may be greater or less. From figure 1 it can be seen that the following inequality always holds for the potential U corresponding to the linear material matched in this way;
U ε˙ f − U ε˙ i > σy ε˙ f − ε˙ i .
(11)
Limit analysis for a general class of yield conditions
405
i
For inequality (11) to remain true, the hatched areas, for ε˙ f > ε˙ i and ε˙ f 6 ε˙ , must be positive and are evidently i so. Equality occurs if, and only if, ε˙ f = ε˙ . In the general case for a general class of yield condition and the class of linear problems described by (2) or (5), the process requires that for an initial strain rate ε˙ iji the linear material is matched to the associated flow rule,
pi
σij = σijLi = Dij kl ε˙ iji − r˙ij = pi
∂U ∂ ε˙ ij
,
(12)
i ε˙ ij
p
where σij is the stress, at yield, associated with ε˙ ij = ε˙ iji . Generally this equation will not be sufficient to define Dij kl and r˙ij uniquely. Equation (12) merely states that a contour of constant U is tangential to the yield pi surface at σij (for a yield surface without corners) in stress space. A wide class of linear materials satisfy this condition. We now assume that, for some members of this class the following generalisation of inequality (11) holds, f
pf f
pi
U ε˙ ij − U ε˙ iji > σij ε˙ ij − σij ε˙ iji , pf
p
(13)
f
where σij is the stress at yield associated with ε˙ ij = ε˙ ij . The methods of constructing solutions to equation (12) for which (13) are satisfied are discussed in section 4. Inequality (13) may be rewritten on noting (12) and (2) as the extended comparative convexity condition; f
U ε˙ ij − U ε˙ iji −
∂U ∂ ε˙ ij
f
i ε˙ ij
pf
pi f
ε˙ ij − ε˙ iji > σij − σij ε˙ ij > 0,
(14)
where the second inequality follows from the maximum work principle (10). Inequality (14) relates the degree of convexity of the linear material and the rigid plastic material. Note that any linear material which satisfies both the matching condition and inequality (13) possesses a convex U and will give rise to unique solutions to the continuum problem associated with the linear material. In the following we assume equality in (14) occurs f if, and only if, ε˙ ij = ε˙ iji , an assumption which may be tested in individual cases. 3. The convergence proof Consider a class of problems where a body of a rigid plastic material occupies a volume V with surface S. Over part of S, namely Su , displacement rates u˙ i = 0 whereas on the remainder of S, namely ST , either surface tractions are zero or loads Ppi (xi ) are applied where P is a positive scalar loading parameter. For an arbitrary distribution of compatible strain rates ε˙ ijc , where the corresponding displacement rates u˙ ci satisfy the boundary conditions on Su , the following upper bound PUc B on the limit load PL exists, Z
PUc B
ST
pi u˙ ci
dS =
Z V
pc
σij ε˙ ijc dV ;
PUc B > PL ,
(15)
pc
where σij is the stress at yield associated with strain rate ε˙ ijc . We now prove the following result; For an initial strain rate field ε˙ ijc = ε˙ iji , a matching linear material is defined by equation (12) which satisfies inequality (13). The continuum problem for this linear material is then solved with load parameter P = PUi B , f the upper bound corresponding to ε˙ iji . The solution to this linear problem produces a (unique) solution ε˙ ij . If
406
A.R.S. Ponter et al.
f
f
PU B is the upper bound corresponding to ε˙ ijc = ε˙ ij then, f
PU B 6 PUi B ,
(16)
f
where equality holds if, and only if, ε˙ iji = ε˙ ij everywhere within V . The proof of this result follows from the potential energy theorem for the linear problem in its rate form. The f solution ε˙ ij to the linear problem minimises the total potential energy of the body (in the rate form) and hence, Z
f
V
U ε˙ ij dV − PUi B
Z
f
ST
pi u˙ i dS 6
Z V
U ε˙ iji dV − PUi B
Z ST
pi u˙ ii dS
(17)
f
with equality only when ε˙ iji = ε˙ ij . Hence from inequality (13) and the definition of the upper bound (15); Z
pf f
pi
f
σij ε˙ ij − σij ε˙ iji dV = PU B
V
6
Z
Z
U V
f
ST
pi u˙ i dS − PUi B
f ε˙ ij dV
−
Z
Z
ST
pi u˙ ii dS
V
U ε˙ iji dV .
(18)
Combining inequalities (17) and (18) gives f
PU B f
Z ST
f
pi u˙ i dS − PUi B R
Z ST
pi u˙ ii dS 6 PUi B
Z ST
f
pi u˙ i dS − PUi B
Z ST
pi u˙ ii dS
(19)
f
and hence PU B 6 PUi B as ST pi u˙ i dS > 0. The convergence of the iterative process follows from the repeated application of this result. The process may be initiated by the solution of a linear problem for arbitrary material constants and load parameter P producing an initial strain rate field ε˙ ijc = ε˙ ij0 . A sequence of solutions ε˙ ijc = ε˙ ijk are then produced where, for the kth solution the material constants are generated by the matching process from the previous solution ε˙ ijk−1 and the load parameter is given by P = PUk−1 ˙ ijk−1 . Then the general result gives B , the upper bound generated by ε that, PUk B 6 PUk−1 B
(20)
and the process generates a sequence of monotonically reducing upper bounds where equality occurs if and only if ε˙ ijk ≡ ε˙ ijk+1 . If the linear solutions are evaluated exactly within a certain function space then the sequence converges to the limit load value defined within that function space. Similarly, if the elastic solutions are derived as the minimum of the potential energy within some restricted class of functional descriptions of the displacement rate field then inequality (17) remains true and the sequence will converge to the minimum upper bound corresponding to that class of displacement rate fields. In the case of finite element solution, for a fixed mesh and p value, the process will converge to the lowest upper bound within the class of displacement rate fields so defined. 4. Conditions under which inequalities (13) and (14) are satisfied by the matching linear material The matching of the linear material to the associated flow rule requires a choice of material constants so that f inequality (13) or, equivalently inequality (14) is satisfied for an arbitrary ε˙ ij . In addition, there is an implicit
Limit analysis for a general class of yield conditions
407
Figure 2. Schematic representation of the quantities in the proof in section 4. f
consistancy assumption that the class of strain rate fields ε˙ ij possesses a location on the yield surface where pf there is an associated stress σij . This is not a trivial assumption as the yield condition and the linear material are independently defined. An example of a situation where this condition would not be satisfied would be a linear material which allows a rate of volume change, whereas the yield condition requires incompressibility. We may prove the following result for the class of matching linear materials for which U and U are convex; f
pi
Inequality (13) is satisfied by an arbitrary ε˙ ij if the contour of constant U which passes through σij entirely encloses the yield surface. f
The proof of this result is assisted by the schematic diagram in figure 2. Corresponding to ε˙ ij there is a strain rate, f
ε˙ ijm = α ε˙ ij for some value of α > 0 so that,
(21)
U σijLm = U σijLi .
(22)
It can easily be shown, from equations (1)–(5) that
α=
U (σ Li ) f f 1 D ε˙ ε˙ 2 ij kl ij kl
1/2
.
408
A.R.S. Ponter et al.
By definition the stress state σijLi for the matching linear material lies on the yield surface where the associated strain rates is ε˙ iji , pi
σijLi = σij .
(23)
If, momentarily, we regard the surface U (σij ) = U (σijLm) as a yield surface which encloses the yield surface f = 0, then; pf
σijLm − σij ε˙ ijm > 0 or, on noting equation (21);
pf f
σijLm − σij ε˙ ij > 0.
(24)
On noting the relationship between U and U , equation (6), and (22),
U ε˙ ijm − σijLm ε˙ ijm = U ε˙ iji − σijLi ε˙ iji .
(25)
Hence on noting (23) and (25) we obtain, f
U ε˙ ij − U ε˙ ijm −
∂U ∂ ε˙ ij
f
m ε˙ ij
f
f
pi
ε˙ ij − ε˙ ijm = U ε˙ ij − U ε˙ iji − σijLmε˙ ij + σij ε˙ iji .
(26)
As we assume U is convex then f
U ε˙ ij − U ε˙ ijm −
∂U ∂ ε˙ ij
f
m ε˙ ij
ε˙ ij − ε˙ ijm > 0.
(27)
Hence, on noting the inequality (24) we obtain from (26) and (27), f
pf f
pi
U ε˙ ij − U ε˙ iji − σij ε˙ ij − σij ε˙ iji > 0
(28)
and inequality (13) is proven. This condition includes, of course, the situation where a contour of constant U coincides exactly with the yield condition. In the case of a Von Mises yield condition this requires the linear material to satisfy the incompressibility condition, i.e. Poissons ratio ν = 0.5. A lower value of ν could result in the violation of inequality (28) in some circumstances, as a contour of U = constant which touches the yield surface would otherwise lie inside; a higher value would violate the condition that U is convex. At the same time the f consistancy assumption, that the strain rate field ε˙ ij given by the linear material may always be associated with the yield condition, is also violated. In summary, in this section we have shown that, in principle, it is possible to construct a sequence of monotonically reducing upper bounds for a general yield condition. In the following section we apply this theory to an important class of isotropic yield conditions which depend upon the Von Mises effective stress and the hydrostatic pressure. 5. Isotropic yield conditions which depend upon the hydrostatic pressure Consider a class of isotropic materials where the yield condition f (σ , p) = 0 depends upon the Von q Mises effective stress σ = 32 σij0 σij0 , where σij0 = σij − pδij denotes the deviatoric stress and p = − 13 σkk
Limit analysis for a general class of yield conditions
409
the (compressive) hydrostatic pressure. The associated flow rule becomes, ε˙ = λ˙
∂f ∂σ
and
ε˙ v = λ˙
∂f , ∂p
(29)
q
where ε˙ = 23 ε˙ ij ε˙ ij and ε˙ v = −˙εkk , the effective strain rate and the rate of volume decrease respectively. The upper bound equation (15) may be simplified to; Z
PUc B
ST
pi u˙ ci dS
=
Z V
σ cy ε˙ c + pyc ε˙ vc dV ,
(30)
where (σ cy , pyc ) denotes the effective stress components at yield associated with (˙ε c , ε˙ vc ). In the following we consider the limit state for a range of yield conditions of this general form, chosen to demonstrate the convergence of the method. In all cases we choose a class of isotropic linear materials defined by the complementary energy function U (σij ) =
1 2
σ 2 /E + (p − pL )2 /K ,
(31)
where E and K are a tensile and a bulk moduli, respectively, and K = Eg(νL ),
where g(νL ) =
1 , 3(1 − 2νL )
(32)
where νL is a Poisson’s ratio. From equation (31) follows the relationships between the stress and strain invariants; ε˙ = σ /E,
ε˙ v = (p − pL )/Eg(νL )
(33)
Hence the matching condition, at the kth iteration, is given by any set of moduli (Ek+1 , νL(k+1), pL(k+1)) which satisfy the equations, k
ε˙ = σ yk /Ek+1
and
ε˙ vk = pyk − pL(k+1) /Ek+1 g νL(k+1) ,
(34)
where (σ yk , pyk ) are the stress values at yield associated with (˙εk , ε˙ vk ). As equations (34) are insufficient to define all three moduli we need to apply, in addition, the condition that a contour of constant U corresponding to (Ek+1 , νL(k+1), pL(k+1)) which touches the yield surface at (σ yk , pyk ) coincides or lies outside the yield surface elsewhere. In the following we first investigate the method for a class of elliptic yield conditions where the possibility exists that the yield surface and a contour of constant complementary energy coincide. We then discuss more general types of yield conditions and show that the class of linear isotropic linear materials given by equations (31) are not able to provide convergent solutions for all isotropic yield conditions of the general form of f (σ , p) = 0. 5.1. A class of elliptic yield conditions Consider a class of isotropic yield conditions of the general form; f = σy
2
+ (py − p0 )2 /g(νp ) − σ02 = 0,
(35)
410
A.R.S. Ponter et al.
Figure 3. The convergence criterion requires, for an elliptical yield surface, that the Poissons ratio for the linear material should be greater than that of the associated flow rule.
where g(νp ) = 1/3(1 − 2νp ) and νp is a plastic Poisson’s ratio. In (p, σ ) space the yield condition describes an ellipse with its centre at p0 and the major and minor axes coincident with space axes. The associated flow rule (9) for the yield condition (35) is given by; ˙ y ε˙ = 2λσ
and
˙ y − p0 )/g(νp ), ε˙ v = 2λ(p
(36)
where (σ y , py ) denotes the point on the yield surface associated with (˙ε , ε˙ v ). The inverse of equations (36) is given by; σ0 ε˙ σy = q 2 ε˙ + ε˙ v2 g(νp )
and
σ0 ε˙ v g(νp ) py = q + p0 . 2 ε˙ + ε˙ v2 g(νp )
(37)
As there are three independant moduli (Ek+1 , νL(k+1) and pL(k+1) ) and two determining equations (34), one modulus may be predefined. If we pre-assign the Poissons ratio νL , we obtain from (34) the following equations for the remaining moduli, Ek+1 = σ yk /˙ε yk
and
pL(k+1) = pyk − σ yk g(νL ) ε˙ νk /˙εk .
(38)
When νL = νp , a contour of constant U coincides with the yield condition and, from the general result in section 4, convergence will occur. For 0.5 > νL > νp > 0 any contour of constant U which touches the yield surface will otherwise lie outside and convergence will also take place. However for 0.5 > νp > νL > 0 any contour of constant U which touches the yield surface will otherwise lie within the yield surface and convergence cannot be assured. These three cases are illustrated in figure 3. 5.2. Numerical examples for elliptic yield surfaces The following examples of the technique using finite element methods have been implemented in the commercial finite element code ABAQUS. The advantages of using a widely available code are evident;
Limit analysis for a general class of yield conditions
411
Figure 4. Uniform block subjected to a uniform compressive stress and plane strain. Eight noded quadrilateral elements.
the method is transferable and the extensive libraries of element types, pre-processor and post-processor systems and efficiency of such codes become available. There is, however, the disadvantage that it is not possible to change the standard procedures within such codes. The implementation involves the solution of a sequence of linear problems where the material moduli are changed at each iteration according to the algorithm equations (38). The upper bound, equation (30), will then converge to the least upper bound associated with the class of displacement fields described by the element structure. This result presumes that all integrations in the evalution of the stiffness matrix and the upper bound are carried out exactly. In the implementation, in ABAQUS integration, routines use Gaussian integration of an order which corresponds to the exact evaluation of the stiffness matrix integrand assuming the material constants are constant within an element. As the algorithm (36) is applied, independently, at each Gauss point there are higher order terms which introduce an error in the stiffness matrix. From extensive experience with the application of the method we are confident that such errors are not significant. In the examples in this and the following section we study two plane strain problems, that of a uniform block of material under uniform loading, shown in figure 4 and the indentation of a half space indentation, shown in figure 5. The elements are eight noded quadrilateral elements. Figure 6 shows the convergence of the upper bound for the indentation problem for the case when νp = νL = 0.5, i.e. for a Von Mises yield condition. The linear solutions are incompressible and hybrid elements were used with a value of νL = 0.499999999. Two solutions are shown in figure 6, corresponding the mesh shown in figure 5 and a more dense mesh where each element has been subdivided into 16 elements of the same type, i.e. 576 elements in total. The finer mesh provides a converged solution which differs from the well known Prandtl solution (Hill, 1950) by 2.5 %. This illustrates the point that the method is capable of producing extremely accurate solutions for sufficiently refined meshes. All the following examples use the less refined mesh of figure 5 as we are primarily concerned with understanding the relationship between the behaviour of the numerical solutions and the convergence criteria. Figures 7 and 8 show the behaviour of the upper bound for the uniform block. Although the problem involves a uniform state of stress and strain, the method needs to seek out the appropriate point on the yield surface corresponding to the correct solution. For νL 6 νp = 0.32818, when convergence is not proven, the solutions
412
A.R.S. Ponter et al.
Figure 5. Symmetric half-space of indentor problems. Eight noded quadriclateral elements.
Figure 6. Convergence of upper bound load parameter P for the indentor problem for a Von Mises yield condition (νL = 0.5). The upper curve A corresponds to the mesh shown in figure 5, whereas the lower curve B corresponds to a mesh where each element has been subdivided into 16 equal elements. The analytical Prandtl solution is shown for comparison.
are non-convergent and stabilise at values which depend upon the choice of νL as shown in figure 7. For νL > νp = 0.32818 when the convergence proof applies the two solutions for νL = 0.32818 and 0.4 converge
Limit analysis for a general class of yield conditions
413
Figure 7. Variation of the load parameter with iteration for the plane strain uniform block of figure 4 for the elliptic yield condition with νp = 0.32818, for νL = 0.1 (curve A) and νL = 0.01 (curve B) when the sufficient condition for convergence is violated. In both cases the apparent converged values are considerably higher that the minimum upper bound (see figure 8).
Figure 8. Variation of the load parameter with iteration for the plane strain uniform block of figure 4 for the elliptical yield condition with νp = 0.32818, for νL = νp = 0.32818 (curve A), νL = 0.4 (curve B) and νL = 0.499 (curve C). In all cases the sufficient condition for convergence is satisfied but the rate of convergence deteriorates as νL increases.
414
A.R.S. Ponter et al.
Figure 9. Variation of the upper bound load parameter P for the indentor problem for elliptic yield surface for νp = 0.32818, νL = 0.1 (curve A) and νL = 0.01 (curve B). The sufficient condition for convergence is not satisfied and convergence does not occur.
to a common, correct, value of the limit load as shown in figure 8. For νL = 0.499 where the aspect ratio of the elliptical contours of constant U is large compared to that of the yield surface, convergence is very slow and has not been approached after 30 iterations. The most rapid convergence occurs when νL = νp . For the indentation for νL < νp = 0.32818, figure 10, convergence to a common value, independent of νL occurs, but convergence becomes progressively slower as νL increases. For νL > νp = 0.32818, figure 9, convergence does not occur and the upper bound values for large numbers of iterations depend upon the value of νL . These solutions demonstrate that the conditions of the convergence theorem are generally required to ensure convergence, but that the rate of convergence is influenced by how closely contours of U match the yield surface. In this case it is possible to obtain a perfect match with νL = νp , the convergence rate is optimised and the number of iterations required for a convergent solution is similar to that for a Von Mises yield condition. For the yield surfaces discussed in the next section such a match is not possible. 5.3. Prager–Drucker yield conditions In this section we consider two yield conditions of the Prager–Drucker type which are used to model the failure behaviour of soils and granular material. The yield surface consists of three district regions, shown in figure 11, a cone section AB which represents shearing failure with dilatation, a cap section BC where the deformation is predominantly volumetric compression and, a corner, B, where the strain rate is indeterminate in direction but lies between the normal vectors to the adjacent sections.
Limit analysis for a general class of yield conditions
415
Figure 10. Variation of the load parameter with iteration for the plane strain indentor problem of figure 10 for the elliptical yield condition with νp = 0.32818, for νL = νp = 0.32818 (curve A), νL = 0.4 (curve B) and νL = 0.499 (curve C). In all cases the sufficient condition for convergence is satisfied but the rate of convergence deteriorates as νL increases.
Figure 11. The two versions of modified Drucker–Prager yield condition used in the examples, where the cone section AB is described by a parabola (upper curve) or a hyperbola (lower curve). The equations for the yield surface are given in the appendix.
For reasons which will become evident when we discuss numerical solutions, we choose two alternative descriptions of the cone section, a parabolic and a hyperbolic form as shown in figure 11. The yield functions and the stress resultants at yield are given in the appendix.
416
A.R.S. Ponter et al.
Figure 12. Convergence for uniform block of figure 4 for Drucker–Prager yield condition with parabolic cone. The sufficient condition for convergence is satisfied.
Figure 13. Complementary strain energy contours for iterations 1–6 for the problem of figure 12.
The matching process has the same characteristics as in the case of an elliptic yield surface, equation (34) repeated here; ε˙ = σ y /E
and
ε˙ v = (py − pL )/Eg(νL ).
(34)
Limit analysis for a general class of yield conditions
417
Figure 14. Variation of load parameter for uniform block for Drucker–Prager yield condition with hyperbolic cone. Non-convergence occurs after the 4th iteration.
Equations (34) are sufficient to define two of the three material parameters (E, νL , pa ). The remaining condition is the requirement that the contour of U = constant which touches the yield surface should otherwise lie outside the yield surface. This condition can be satisfied for both yield surfaces for the corner and the cap by ensuring that the U = constant contour passes through the apex p = −pt and σ = 0. For the cone section U = constant approaches the yield surface for sufficiently large values of pL in the parabolic case, but there exists no values of the moduli so that the contour always lies outside the yield surface in the hyperbolic case. Hence, an isotropic linear material may not always exist which allows convergent solutions for an isotropic yield condition and there are limits to the range of yield conditions which may be used. This is not an inherent restriction on the application of the method as it is possible to overcome this problem by using an anisotropic linear material, but that possibility will not be pursued here. This example demonstrates that the main restriction on the application of the method derives from the condition that the U = constant contour lies outside the yield surface and this condition can be difficult to apply in cases where the curvature of the yield surface is locally small. Hence the method is most appropriate where the yield condition is similar in form to the appropriate U = constant surface. Figures 12 and 14 show the solutions for the uniform block for each of the yield conditions and figures 13 and 15 show U = constant surfaces for the first few iterations in each case. The initial solution is generated by choosing the moduli so that a U = constant contour passes through the two yield points on the p axis. In accordance with the convergence results in section 4, the solution for the parabolic cone converges but the solution for the hyperbolic cone fails to converge and this is borne out by the numerical results. Figure 16 shows the convergence of the plane strain indentor problem of figure 5 for the parabolic cone with the mesh as shown (curve A) and a refined mesh with each element subdivided into 16 elements (curve B). In both cases the changes in the load parameter are confined to the 5th significant figure after 20 iterations
418
A.R.S. Ponter et al.
Figure 15. Complementary strain energy contours for iterations 1–6 for the problem of figure 14. At the 4th iteration the sufficient condition for convergence is violated and the solution fails to converge subsequently.
Figure 16. Convergence of the load parameter for the plain strain indentation problem of figure 5 for the Drucker–Prager yield condition with a parabolic cone section. Curve A corresponds to the mesh shown in figure 5 and curve B to the same mesh where each element has been subdivided into 16 equal elements.
and convergence occurs. Figure 17 shows the variation of the upper bound with iterations for the hyperbolic cone. The sufficient condition for convergence is not satisfied and the solutions fail to converge with increasing iterations.
Limit analysis for a general class of yield conditions
419
Figure 17. Non-convergence of the load parameter for the indentation problem of figure 5 for the Drucker–Prager yield condition with a parabolic cone section.
6. Conclusions The paper discusses the convergence criteria for the method previously discussed by Carter and Ponter (1997a), for a Von Mises yield condition, but here extended to a general yield condition. The method is re-interpreted as a programming method where a linear material is matched to the current strain field and the yield condition in an iterative process and a sufficient condition for convergence is derived. The set of numerical examples provide the following indications. (1) The convergence of the method is most rapid when contours of constant complementary energy for the linear material are close or identical to the yield surface. (2) Although the convergence criterion is a sufficient condition, there is clear evidence that violation of the condition is likely to lead to non-convergence or apparent convergence to a non-optimal value of a load parameter. (3) The method is capable of providing efficient converged solutions for complex yield surfaces although there are cases where an isotropic yield condition may require an anisotropic linear material. (4) The method requires the evaluation of a sequence of linear finite element problems and implimentaion is possible within commercial finite element codes, in this case ABAQUS, by the use of user routines. As a result the method is more assessable than others which require specialist software. The method is applied here to limit analysis problems. The extension to shakedown limits and related solutions requires essentially the same algorithms (Ponter and Engelhardt, 2000), and hence the work of this paper may be regarded as the first step to the generation of such solutions for a general yield surface. Appendix. Yield conditions and stress resultants (σ y , py ) at yield A.1. Cap section BC Yield function:
fc = (p − pa )2 + σ 2 (pb − pa )2 /σ02
1/2
− (pb − pa ) = 0.
(A1)
420
A.R.S. Ponter et al.
Stress resultants at yield:
σ y = σ0 ε˙ / ε˙ 2 +
pb − p a σ0
2
1/2
ε˙ v2
py = pa + (pb − pa ) ε˙ v / ε˙ v2 +
,
σ0 pb − p a
2
1/2
ε˙ 2
.
(A2)
A.2. Parabolic section AB Yield function:
fp = (pa + pt ) Stress resultants at yield:
ε˙ v − ε˙
For
> tan α =
σ σ0
2
− (p + pt ) = 0.
σ0 ; 2(pa + pt )
(A3)
(A4)
ε˙ ε˙ σ02 σ02 , py = −pt + 2(pa + pt ) ε˙ v 4(pa + pt ) ε˙ v ε˙ v σ0 and for 06− 6 tan α = ; 2(pa + pt ) ε˙ σ y = σ0 , p y = pa .
2
σy = −
(A5) (A6) (A7)
A.3. Hyperbolic section BC Yield function: fh = l02 + σ 2 l0 =
where
1/2
− (p tan β − d),
σ02 − (pa + pt )2 tan2 β 2 (pa + pt ) tan β
and
(A8)
d = pt tan β − l0 .
(A9)
Stress resultants at yield
For
ε˙ v − ε˙
σ02 − l02 > tan β, σ02
l0 tan β ε˙ σy = 2 2 , (˙ε tan β − ε˙ v2 )1/2
and for σ y = σ0 ,
1 l0 ε˙ v py = d− 2 2 tan β (˙ε tan β − ε˙ v2 )1/2
(A10)
ε˙ v σ02 − l02 06− 6 tan β σ02 ε˙ p y = pa .
(A11)
Acknowledgements The authors gratefully acknowledge the support during the course of this work of Nuclear Electric Ltd. and the University of Leicester (for ME) and of a European Commission Grant HCM Network Contract ERBCHRXCT940629 (for PF).
Limit analysis for a general class of yield conditions
421
References Collins I.F., Boulbibane, Application of Kinematic Shakedown Theorems to Pavement Design, Proc. Euromech 385, Inelastic Analysis of Structures under Variable Loading: Theory and Applications - Technical University of Aachen, September 1998. Goodall I.W., Goodman A.M., Chell G.C., Ainsworth R.A., Williams J.A., R5: An Assessment Procedure for the High Temperature Response of Structures, Nuclear Electric Ltd., Barnwood, Gloucester, Report, 1991. Hill R., The Mathematical Theory of Plasticity, Oxford University Press, London, 1950. Liu Y.H., Cen Z.Z., Xu B.Y., A numerical method for plastic limit analysis of 3-D structures, Int. J. Sol. Struct. 32 (1995) 1645–1658. Marriott D.L., Evaluation of deformation on load control of stresses under inelastic conditions using elastic finite element stress analysis, ASME Pressure Vessels and Piping Conference, 1985, Pittsburgh, ‘Codes and Standards and Applications for Design and Analysis of Pressure Vessel and Piping Components’, eds: R. Seshadri et al., ASME PVP, 1985, 136, p. 3. Mackenzie D., Boyle T., A method of estimating limit loads by iterative elastic analysis. I-simple examples, Int. J. Pres. Ves. & Piping 53 (1993) 77. Mackenzie D., Nadarajah C., Shi J., Boyle J.T., Simple bounds on limit loads by elastic finite element analysis, ASME J. Press. Ves. Tech. 115 (1993) 27. Nadarajah C., Mackenzie D., Boyle J.T., A method of estimating limit loads by iterative elastic analysis II-nozzle sphere intersections with internal pressure and radial load, Int. J. Pres. Ves & Piping 53 (1993) 97. Ponter A.R.S., Carter K.F., Limit state solutions, based uponlinear solutions with a spatially varying elastic modulus, Meth, Appl. Mech. Eng. 140 (1997a) 237–258. Ponter A.R.S., Carter K.F., Shakedown state simulation techniques based upon linear elastic solutions, Meth. Appl. Mech. Eng. 140 (1997b) 259– 279. Ponter A.R.S., Engelhardt M., Shakedown limits for a general yield condition: implementation and application for a Von Mises yield condition, Eur. J. Mech. A/Solids 19 (2000) 423–445 (this issue). Seshadri R., Fernando C.P.D., Limit loads of mechanical components and structures based on linear elastic solutions using the GLOSS R-Node method, Trans ASME PVP, Vol 210-2, San Diego, 1991, pp. 125–134. Shi J., Mackenzie D., Boyle J.T., A method of estimating limit loads by iterative elastic analysis. I-simple examples, Int. J. Pres. Ves. & Piping 53 (1993) 77. Zhang Y.G., Zhang P., Lu M.W., Computational limit analysis of rigid plastic bodies in plane strain, Acta Mech. Sol. Sinica 6 (1993) 341–348.
Limit analysis for a general class of yield conditions Alan R.S. Ponter a,1 , Paolo Fuschi b , Markus Engelhardt a a Department of Engineering, University of Leicester, Leicester, LE1 7RH, UK b DASTEC - Departimento Arte Scienza Tecnica del Costruire, University of Reggio Calabria, Via Melissari Feo di Vito, 89124
Reggio Calabria, Italy Abstract – The paper describes a generalisation of the programming method described by Ponter and Carter (1997) for the evaluation of optimal upper bounds on the limit load of a body composed of a rigid/perfectly plastic material. The method is based upon similar principles to the ‘Elastic Compensation’ method which has been used for design calculations for some years but re-interpreted as a non-linear programming method. A sufficient condition for convergence is derived which relates properties of the yield surface to those of the linear solutions solved at each iteration. The method is demonstrated through an application to a Drucker–Prager yield condition in terms of the Von Mises effective stress and the hydrostatic pressure. Implementation is shown to be possible using the user routines in a commercial finite element code, ABAQUS. The examples chosen indicate that stable convergent solutions may be obtained. There are, however, limits to the application of the method if isotropic linear solutions are used for an isotropic yield surface. In an accompanying paper (Ponter and Engelhardt, 2000) the method is extended to shakedown and related problems. 2000 Éditions scientifiques et médicales Elsevier SAS plasticity / limit / programming / finite elements / Drucker–Prager / life assessment / design
1. Introduction There has been a growth of interest in recent times in methods of directly computing performance indicators for structures subjected to complex loading as an aid to design and for the implementaion of life assessment methods (Goodall et al., 1991). The simplest of these indicators is the limit load for a rigid/perfectly plastic body as the first step towards shakedown limits and related solutions. In previous papers (Ponter and Carter, 1997a, b) methods were described for the computation of optimal upper bounds on the limit state and shakedown limit of a perfectly plastic body, with a Von Mises yield condition, subjected to variable load and temperature. The purpose of this paper is to provide a generalisation of this method to a general yield condition. The general characteristics of the method are the same as the ‘Elastic Compensation’ method which has been described by a number of authors (see Marriott, 1985; Mackenzie and Boyle, 1993; Nadarajah et al., 1993; Shi et al., 1993; Mackenzie et al., 1993; Seshadri and Fernando, 1991). A sequence of linear problems with spatial varying moduli are generated where, at each iteration, the stress is brought within yield by changes of the moduli, assuming a fixed strain rate distribution. Interpreting the method as an upper bound programming technique, Ponter and Carter (1997a and b) demonstrated that the method provides a monotonically reducing sequence of upper bounds which converges to the least upper bound for a class of displacement fields. When implemented in a finite element code, the method proves to be very stable and converged solutions can usually be obtained after about 10 iterations. The lower bound derived from the resulting stress distributions, which formed the primary interest of the ‘Elastic Compensation’ method, were shown to converge, not necessarily monotonically, to the least upper bound. Other programming methods, both linear and non-linear, have been described in the literature of which Zhang et al. (1993) and Liu et al. (1995) provide recent examples. 1 E-mail: [email protected]
402
A.R.S. Ponter et al.
In this paper we discuss the generalisation of the method to an arbitrary convex yield surface. The arguments are confined to limit analysis and the extension to shakedown analysis and related cyclic problems will be discussed elsewhere (Ponter and Engelhardt, 2000). The generalisation of the method may be understood in two stages. The first involves a re-interpretation of the basis of the method as a general programming technique. At each iteration the moduli of the linear material are adjusted so that, for the current strain rate distribution, the linear material is matched to the yield condition. However, this condition alone is insufficient to uniquely define the moduli or for convergence to occur and additional conditions need to be applied. In the following sections we derive a sufficient condition for convergence and then interpret this condition in terms of an algorithm for the evaluation of the moduli of the linear material. Essentially, the surface of complementary energy density for the linear problem which touches the yield surface at the position of the current strain rate must otherwise lie outside the yield surface in stress space. The method is illustrated with reference to two classes of isotropic yield conditions which depend upon the Von Mises effective stress σ¯ and the hydrostatic pressure p. The first is for yield conditions which describe an ellipse in (σ¯ , p) space where the yield conditions and contours of constant complementary energy density may coincide and where the Von Mises yield condition is a special case. We find, as expected, that the method converges when the sufficient condition applies but also that convergence does not occur for all those cases studied where the condition does not apply. The method is then applied to two cases of a Drucker–Prager model for granular material where the application of the method in its general case can be demonstrated. We find that where the yield surface has regions of low curvature an isotropic yield condition may require the definition of an anisotropic linear material. All the numerical examples given in the paper have been generated using the user routines of a general purpose commercial finite element code, ABAQUS. These methods may, therefore, be viewed as general purpose methods which are capable of having significant use in design and where very accurate values of significant structural performance indicators are required.
2. Convexity conditions 2.1. The linear material We consider a class of linear materials defined by the strain rate potential U (˙εij )given by,
U (˙εij ) = 12 Dij kl ε˙ ij − r˙ij ε˙ kl − r˙kl ,
(1)
where Dij kl is the compliance tensor of coefficients which satisfy the usual symmetry properties and r˙ij is a predefined initial strain rate. The stress corresponding to a strain rate is given by, σijL =
∂U = Dij kl ε˙ kl − r˙kl . ∂ ε˙ ij
(2)
The complementary potential U (σij ) is given by, U (σij ) = 12 Cij kl (σij − qij )(σkl − qkl ),
(3)
Limit analysis for a general class of yield conditions
403
where Cij kl is the inverse of Dij kl and qij is an initial stress given by, r˙ij = −Cij kl qkl ,
qij = −Dij kl r˙kl .
or equivalently
(4)
and hence, ε˙ ij =
∂U L L = Cij kl σkl − qkl . ∂σij
(5)
The two potential functions are related by,
U σijL + U ε˙ ij = σijL ε˙ ij − qij r˙ij .
(6)
Above and in the following we denote stresses derived from such a linear material with a superscript L, to distinguish them from stresses associated with the yield condition which will be denoted by a superscript p. At this stage we note that, for uniqueness of the solution of continuum problems associated with the linear material, we require that both U and U are convex, i.e. U
ε˙ ij1
−U
ε˙ ij2
∂U − ∂ ε˙ ij
2 εij
ε˙ ij1 − ε˙ ij2 > 0
(7a)
with equality only when ε˙ ij1 = ε˙ ij2 , and
U σij1 − U σij2 −
∂U ∂σij
2 εij
σij1 − σij2 > 0
(7b)
with equality only when σij1 = σij2 . 2.2. The yield condition Consider a rigid plastic material with a yield condition, f (σij ) 6 0
(8)
and associated flow rule, ∂f , f (σij ) = 0. ∂σij We assume that the yield function is convex and that the maximum work principle holds, ε˙ ij = λ˙ p
σij − σij∗ ε˙ ij > 0, p
p
(9)
(10)
where σij is the stress associated with ε˙ ij at yield and σij∗ is any state of stress which satisfies the yield condition (8). p
p
2.3. The relationship between the linear material and the associated flow rule In (Ponter and Carter, 1997a), the convergence proof was constructed by assuming that the appropriate class of linear materials for a Von Mises yield condition was incompressible and isotropic. In the general case we
404
A.R.S. Ponter et al.
Figure 1. The relative convexity condition, equation (11), for a Von Mises yield condition requires that the shaded areas should have a positive area as shown.
require a method of constructing linear solutions so that convergence is assured. The key to this choice can be understood from a simple example given by Ponter and Carter (1997a), where we show that the iterative process can be understood as a programming method where the slope of the potential energy of the linear material (in the rate form) is matched to the slope of the upper bound to the limit load at each strain rate distribution. The solution of the linear problem produces the minimum of the potential energy and, at the same time reduces the upper bound to the limit load. In the following we describe a generalisation of this argument, in two stages. Initially we define a relationship between the convexity conditions (7) and the maximum work principle (10) which places constraints on the class of linear materials. We then show that this relationship is sufficient for convergence to take place in an iterative process for the limit load. 2.4. The relative convexity condition Consider a Von Mises yield condition f = σ − σy = 0 where σ is the Von Mises effective stress and σy is the uniaxial yield stress. In this case the linear material is an incompressible isotropic material described by a single material constant E and ε˙ = σ /E, where ε˙ denotes the effective strain rate. The matching process can be described graphically, as shown in figure 1. Suppose, at some stage in the iterative process, the current strain rate distribution is given by ε˙ i . The linear material may be matched to the yield condition by choosing E = σy /˙εi so that both the linear material and the associated flow rule yield the same state of effective stress for effective strain rate ε˙ i . The solution of the subsequent linear problem for this choice of E will produce a final strain rate ε˙ f which will generally differ from ε˙ i and may be greater or less. From figure 1 it can be seen that the following inequality always holds for the potential U corresponding to the linear material matched in this way;
U ε˙ f − U ε˙ i > σy ε˙ f − ε˙ i .
(11)
Limit analysis for a general class of yield conditions
405
i
For inequality (11) to remain true, the hatched areas, for ε˙ f > ε˙ i and ε˙ f 6 ε˙ , must be positive and are evidently i so. Equality occurs if, and only if, ε˙ f = ε˙ . In the general case for a general class of yield condition and the class of linear problems described by (2) or (5), the process requires that for an initial strain rate ε˙ iji the linear material is matched to the associated flow rule,
pi
σij = σijLi = Dij kl ε˙ iji − r˙ij = pi
∂U ∂ ε˙ ij
,
(12)
i ε˙ ij
p
where σij is the stress, at yield, associated with ε˙ ij = ε˙ iji . Generally this equation will not be sufficient to define Dij kl and r˙ij uniquely. Equation (12) merely states that a contour of constant U is tangential to the yield pi surface at σij (for a yield surface without corners) in stress space. A wide class of linear materials satisfy this condition. We now assume that, for some members of this class the following generalisation of inequality (11) holds, f
pf f
pi
U ε˙ ij − U ε˙ iji > σij ε˙ ij − σij ε˙ iji , pf
p
(13)
f
where σij is the stress at yield associated with ε˙ ij = ε˙ ij . The methods of constructing solutions to equation (12) for which (13) are satisfied are discussed in section 4. Inequality (13) may be rewritten on noting (12) and (2) as the extended comparative convexity condition; f
U ε˙ ij − U ε˙ iji −
∂U ∂ ε˙ ij
f
i ε˙ ij
pf
pi f
ε˙ ij − ε˙ iji > σij − σij ε˙ ij > 0,
(14)
where the second inequality follows from the maximum work principle (10). Inequality (14) relates the degree of convexity of the linear material and the rigid plastic material. Note that any linear material which satisfies both the matching condition and inequality (13) possesses a convex U and will give rise to unique solutions to the continuum problem associated with the linear material. In the following we assume equality in (14) occurs f if, and only if, ε˙ ij = ε˙ iji , an assumption which may be tested in individual cases. 3. The convergence proof Consider a class of problems where a body of a rigid plastic material occupies a volume V with surface S. Over part of S, namely Su , displacement rates u˙ i = 0 whereas on the remainder of S, namely ST , either surface tractions are zero or loads Ppi (xi ) are applied where P is a positive scalar loading parameter. For an arbitrary distribution of compatible strain rates ε˙ ijc , where the corresponding displacement rates u˙ ci satisfy the boundary conditions on Su , the following upper bound PUc B on the limit load PL exists, Z
PUc B
ST
pi u˙ ci
dS =
Z V
pc
σij ε˙ ijc dV ;
PUc B > PL ,
(15)
pc
where σij is the stress at yield associated with strain rate ε˙ ijc . We now prove the following result; For an initial strain rate field ε˙ ijc = ε˙ iji , a matching linear material is defined by equation (12) which satisfies inequality (13). The continuum problem for this linear material is then solved with load parameter P = PUi B , f the upper bound corresponding to ε˙ iji . The solution to this linear problem produces a (unique) solution ε˙ ij . If
406
A.R.S. Ponter et al.
f
f
PU B is the upper bound corresponding to ε˙ ijc = ε˙ ij then, f
PU B 6 PUi B ,
(16)
f
where equality holds if, and only if, ε˙ iji = ε˙ ij everywhere within V . The proof of this result follows from the potential energy theorem for the linear problem in its rate form. The f solution ε˙ ij to the linear problem minimises the total potential energy of the body (in the rate form) and hence, Z
f
V
U ε˙ ij dV − PUi B
Z
f
ST
pi u˙ i dS 6
Z V
U ε˙ iji dV − PUi B
Z ST
pi u˙ ii dS
(17)
f
with equality only when ε˙ iji = ε˙ ij . Hence from inequality (13) and the definition of the upper bound (15); Z
pf f
pi
f
σij ε˙ ij − σij ε˙ iji dV = PU B
V
6
Z
Z
U V
f
ST
pi u˙ i dS − PUi B
f ε˙ ij dV
−
Z
Z
ST
pi u˙ ii dS
V
U ε˙ iji dV .
(18)
Combining inequalities (17) and (18) gives f
PU B f
Z ST
f
pi u˙ i dS − PUi B R
Z ST
pi u˙ ii dS 6 PUi B
Z ST
f
pi u˙ i dS − PUi B
Z ST
pi u˙ ii dS
(19)
f
and hence PU B 6 PUi B as ST pi u˙ i dS > 0. The convergence of the iterative process follows from the repeated application of this result. The process may be initiated by the solution of a linear problem for arbitrary material constants and load parameter P producing an initial strain rate field ε˙ ijc = ε˙ ij0 . A sequence of solutions ε˙ ijc = ε˙ ijk are then produced where, for the kth solution the material constants are generated by the matching process from the previous solution ε˙ ijk−1 and the load parameter is given by P = PUk−1 ˙ ijk−1 . Then the general result gives B , the upper bound generated by ε that, PUk B 6 PUk−1 B
(20)
and the process generates a sequence of monotonically reducing upper bounds where equality occurs if and only if ε˙ ijk ≡ ε˙ ijk+1 . If the linear solutions are evaluated exactly within a certain function space then the sequence converges to the limit load value defined within that function space. Similarly, if the elastic solutions are derived as the minimum of the potential energy within some restricted class of functional descriptions of the displacement rate field then inequality (17) remains true and the sequence will converge to the minimum upper bound corresponding to that class of displacement rate fields. In the case of finite element solution, for a fixed mesh and p value, the process will converge to the lowest upper bound within the class of displacement rate fields so defined. 4. Conditions under which inequalities (13) and (14) are satisfied by the matching linear material The matching of the linear material to the associated flow rule requires a choice of material constants so that f inequality (13) or, equivalently inequality (14) is satisfied for an arbitrary ε˙ ij . In addition, there is an implicit
Limit analysis for a general class of yield conditions
407
Figure 2. Schematic representation of the quantities in the proof in section 4. f
consistancy assumption that the class of strain rate fields ε˙ ij possesses a location on the yield surface where pf there is an associated stress σij . This is not a trivial assumption as the yield condition and the linear material are independently defined. An example of a situation where this condition would not be satisfied would be a linear material which allows a rate of volume change, whereas the yield condition requires incompressibility. We may prove the following result for the class of matching linear materials for which U and U are convex; f
pi
Inequality (13) is satisfied by an arbitrary ε˙ ij if the contour of constant U which passes through σij entirely encloses the yield surface. f
The proof of this result is assisted by the schematic diagram in figure 2. Corresponding to ε˙ ij there is a strain rate, f
ε˙ ijm = α ε˙ ij for some value of α > 0 so that,
(21)
U σijLm = U σijLi .
(22)
It can easily be shown, from equations (1)–(5) that
α=
U (σ Li ) f f 1 D ε˙ ε˙ 2 ij kl ij kl
1/2
.
408
A.R.S. Ponter et al.
By definition the stress state σijLi for the matching linear material lies on the yield surface where the associated strain rates is ε˙ iji , pi
σijLi = σij .
(23)
If, momentarily, we regard the surface U (σij ) = U (σijLm) as a yield surface which encloses the yield surface f = 0, then; pf
σijLm − σij ε˙ ijm > 0 or, on noting equation (21);
pf f
σijLm − σij ε˙ ij > 0.
(24)
On noting the relationship between U and U , equation (6), and (22),
U ε˙ ijm − σijLm ε˙ ijm = U ε˙ iji − σijLi ε˙ iji .
(25)
Hence on noting (23) and (25) we obtain, f
U ε˙ ij − U ε˙ ijm −
∂U ∂ ε˙ ij
f
m ε˙ ij
f
f
pi
ε˙ ij − ε˙ ijm = U ε˙ ij − U ε˙ iji − σijLmε˙ ij + σij ε˙ iji .
(26)
As we assume U is convex then f
U ε˙ ij − U ε˙ ijm −
∂U ∂ ε˙ ij
f
m ε˙ ij
ε˙ ij − ε˙ ijm > 0.
(27)
Hence, on noting the inequality (24) we obtain from (26) and (27), f
pf f
pi
U ε˙ ij − U ε˙ iji − σij ε˙ ij − σij ε˙ iji > 0
(28)
and inequality (13) is proven. This condition includes, of course, the situation where a contour of constant U coincides exactly with the yield condition. In the case of a Von Mises yield condition this requires the linear material to satisfy the incompressibility condition, i.e. Poissons ratio ν = 0.5. A lower value of ν could result in the violation of inequality (28) in some circumstances, as a contour of U = constant which touches the yield surface would otherwise lie inside; a higher value would violate the condition that U is convex. At the same time the f consistancy assumption, that the strain rate field ε˙ ij given by the linear material may always be associated with the yield condition, is also violated. In summary, in this section we have shown that, in principle, it is possible to construct a sequence of monotonically reducing upper bounds for a general yield condition. In the following section we apply this theory to an important class of isotropic yield conditions which depend upon the Von Mises effective stress and the hydrostatic pressure. 5. Isotropic yield conditions which depend upon the hydrostatic pressure Consider a class of isotropic materials where the yield condition f (σ , p) = 0 depends upon the Von q Mises effective stress σ = 32 σij0 σij0 , where σij0 = σij − pδij denotes the deviatoric stress and p = − 13 σkk
Limit analysis for a general class of yield conditions
409
the (compressive) hydrostatic pressure. The associated flow rule becomes, ε˙ = λ˙
∂f ∂σ
and
ε˙ v = λ˙
∂f , ∂p
(29)
q
where ε˙ = 23 ε˙ ij ε˙ ij and ε˙ v = −˙εkk , the effective strain rate and the rate of volume decrease respectively. The upper bound equation (15) may be simplified to; Z
PUc B
ST
pi u˙ ci dS
=
Z V
σ cy ε˙ c + pyc ε˙ vc dV ,
(30)
where (σ cy , pyc ) denotes the effective stress components at yield associated with (˙ε c , ε˙ vc ). In the following we consider the limit state for a range of yield conditions of this general form, chosen to demonstrate the convergence of the method. In all cases we choose a class of isotropic linear materials defined by the complementary energy function U (σij ) =
1 2
σ 2 /E + (p − pL )2 /K ,
(31)
where E and K are a tensile and a bulk moduli, respectively, and K = Eg(νL ),
where g(νL ) =
1 , 3(1 − 2νL )
(32)
where νL is a Poisson’s ratio. From equation (31) follows the relationships between the stress and strain invariants; ε˙ = σ /E,
ε˙ v = (p − pL )/Eg(νL )
(33)
Hence the matching condition, at the kth iteration, is given by any set of moduli (Ek+1 , νL(k+1), pL(k+1)) which satisfy the equations, k
ε˙ = σ yk /Ek+1
and
ε˙ vk = pyk − pL(k+1) /Ek+1 g νL(k+1) ,
(34)
where (σ yk , pyk ) are the stress values at yield associated with (˙εk , ε˙ vk ). As equations (34) are insufficient to define all three moduli we need to apply, in addition, the condition that a contour of constant U corresponding to (Ek+1 , νL(k+1), pL(k+1)) which touches the yield surface at (σ yk , pyk ) coincides or lies outside the yield surface elsewhere. In the following we first investigate the method for a class of elliptic yield conditions where the possibility exists that the yield surface and a contour of constant complementary energy coincide. We then discuss more general types of yield conditions and show that the class of linear isotropic linear materials given by equations (31) are not able to provide convergent solutions for all isotropic yield conditions of the general form of f (σ , p) = 0. 5.1. A class of elliptic yield conditions Consider a class of isotropic yield conditions of the general form; f = σy
2
+ (py − p0 )2 /g(νp ) − σ02 = 0,
(35)
410
A.R.S. Ponter et al.
Figure 3. The convergence criterion requires, for an elliptical yield surface, that the Poissons ratio for the linear material should be greater than that of the associated flow rule.
where g(νp ) = 1/3(1 − 2νp ) and νp is a plastic Poisson’s ratio. In (p, σ ) space the yield condition describes an ellipse with its centre at p0 and the major and minor axes coincident with space axes. The associated flow rule (9) for the yield condition (35) is given by; ˙ y ε˙ = 2λσ
and
˙ y − p0 )/g(νp ), ε˙ v = 2λ(p
(36)
where (σ y , py ) denotes the point on the yield surface associated with (˙ε , ε˙ v ). The inverse of equations (36) is given by; σ0 ε˙ σy = q 2 ε˙ + ε˙ v2 g(νp )
and
σ0 ε˙ v g(νp ) py = q + p0 . 2 ε˙ + ε˙ v2 g(νp )
(37)
As there are three independant moduli (Ek+1 , νL(k+1) and pL(k+1) ) and two determining equations (34), one modulus may be predefined. If we pre-assign the Poissons ratio νL , we obtain from (34) the following equations for the remaining moduli, Ek+1 = σ yk /˙ε yk
and
pL(k+1) = pyk − σ yk g(νL ) ε˙ νk /˙εk .
(38)
When νL = νp , a contour of constant U coincides with the yield condition and, from the general result in section 4, convergence will occur. For 0.5 > νL > νp > 0 any contour of constant U which touches the yield surface will otherwise lie outside and convergence will also take place. However for 0.5 > νp > νL > 0 any contour of constant U which touches the yield surface will otherwise lie within the yield surface and convergence cannot be assured. These three cases are illustrated in figure 3. 5.2. Numerical examples for elliptic yield surfaces The following examples of the technique using finite element methods have been implemented in the commercial finite element code ABAQUS. The advantages of using a widely available code are evident;
Limit analysis for a general class of yield conditions
411
Figure 4. Uniform block subjected to a uniform compressive stress and plane strain. Eight noded quadrilateral elements.
the method is transferable and the extensive libraries of element types, pre-processor and post-processor systems and efficiency of such codes become available. There is, however, the disadvantage that it is not possible to change the standard procedures within such codes. The implementation involves the solution of a sequence of linear problems where the material moduli are changed at each iteration according to the algorithm equations (38). The upper bound, equation (30), will then converge to the least upper bound associated with the class of displacement fields described by the element structure. This result presumes that all integrations in the evalution of the stiffness matrix and the upper bound are carried out exactly. In the implementation, in ABAQUS integration, routines use Gaussian integration of an order which corresponds to the exact evaluation of the stiffness matrix integrand assuming the material constants are constant within an element. As the algorithm (36) is applied, independently, at each Gauss point there are higher order terms which introduce an error in the stiffness matrix. From extensive experience with the application of the method we are confident that such errors are not significant. In the examples in this and the following section we study two plane strain problems, that of a uniform block of material under uniform loading, shown in figure 4 and the indentation of a half space indentation, shown in figure 5. The elements are eight noded quadrilateral elements. Figure 6 shows the convergence of the upper bound for the indentation problem for the case when νp = νL = 0.5, i.e. for a Von Mises yield condition. The linear solutions are incompressible and hybrid elements were used with a value of νL = 0.499999999. Two solutions are shown in figure 6, corresponding the mesh shown in figure 5 and a more dense mesh where each element has been subdivided into 16 elements of the same type, i.e. 576 elements in total. The finer mesh provides a converged solution which differs from the well known Prandtl solution (Hill, 1950) by 2.5 %. This illustrates the point that the method is capable of producing extremely accurate solutions for sufficiently refined meshes. All the following examples use the less refined mesh of figure 5 as we are primarily concerned with understanding the relationship between the behaviour of the numerical solutions and the convergence criteria. Figures 7 and 8 show the behaviour of the upper bound for the uniform block. Although the problem involves a uniform state of stress and strain, the method needs to seek out the appropriate point on the yield surface corresponding to the correct solution. For νL 6 νp = 0.32818, when convergence is not proven, the solutions
412
A.R.S. Ponter et al.
Figure 5. Symmetric half-space of indentor problems. Eight noded quadriclateral elements.
Figure 6. Convergence of upper bound load parameter P for the indentor problem for a Von Mises yield condition (νL = 0.5). The upper curve A corresponds to the mesh shown in figure 5, whereas the lower curve B corresponds to a mesh where each element has been subdivided into 16 equal elements. The analytical Prandtl solution is shown for comparison.
are non-convergent and stabilise at values which depend upon the choice of νL as shown in figure 7. For νL > νp = 0.32818 when the convergence proof applies the two solutions for νL = 0.32818 and 0.4 converge
Limit analysis for a general class of yield conditions
413
Figure 7. Variation of the load parameter with iteration for the plane strain uniform block of figure 4 for the elliptic yield condition with νp = 0.32818, for νL = 0.1 (curve A) and νL = 0.01 (curve B) when the sufficient condition for convergence is violated. In both cases the apparent converged values are considerably higher that the minimum upper bound (see figure 8).
Figure 8. Variation of the load parameter with iteration for the plane strain uniform block of figure 4 for the elliptical yield condition with νp = 0.32818, for νL = νp = 0.32818 (curve A), νL = 0.4 (curve B) and νL = 0.499 (curve C). In all cases the sufficient condition for convergence is satisfied but the rate of convergence deteriorates as νL increases.
414
A.R.S. Ponter et al.
Figure 9. Variation of the upper bound load parameter P for the indentor problem for elliptic yield surface for νp = 0.32818, νL = 0.1 (curve A) and νL = 0.01 (curve B). The sufficient condition for convergence is not satisfied and convergence does not occur.
to a common, correct, value of the limit load as shown in figure 8. For νL = 0.499 where the aspect ratio of the elliptical contours of constant U is large compared to that of the yield surface, convergence is very slow and has not been approached after 30 iterations. The most rapid convergence occurs when νL = νp . For the indentation for νL < νp = 0.32818, figure 10, convergence to a common value, independent of νL occurs, but convergence becomes progressively slower as νL increases. For νL > νp = 0.32818, figure 9, convergence does not occur and the upper bound values for large numbers of iterations depend upon the value of νL . These solutions demonstrate that the conditions of the convergence theorem are generally required to ensure convergence, but that the rate of convergence is influenced by how closely contours of U match the yield surface. In this case it is possible to obtain a perfect match with νL = νp , the convergence rate is optimised and the number of iterations required for a convergent solution is similar to that for a Von Mises yield condition. For the yield surfaces discussed in the next section such a match is not possible. 5.3. Prager–Drucker yield conditions In this section we consider two yield conditions of the Prager–Drucker type which are used to model the failure behaviour of soils and granular material. The yield surface consists of three district regions, shown in figure 11, a cone section AB which represents shearing failure with dilatation, a cap section BC where the deformation is predominantly volumetric compression and, a corner, B, where the strain rate is indeterminate in direction but lies between the normal vectors to the adjacent sections.
Limit analysis for a general class of yield conditions
415
Figure 10. Variation of the load parameter with iteration for the plane strain indentor problem of figure 10 for the elliptical yield condition with νp = 0.32818, for νL = νp = 0.32818 (curve A), νL = 0.4 (curve B) and νL = 0.499 (curve C). In all cases the sufficient condition for convergence is satisfied but the rate of convergence deteriorates as νL increases.
Figure 11. The two versions of modified Drucker–Prager yield condition used in the examples, where the cone section AB is described by a parabola (upper curve) or a hyperbola (lower curve). The equations for the yield surface are given in the appendix.
For reasons which will become evident when we discuss numerical solutions, we choose two alternative descriptions of the cone section, a parabolic and a hyperbolic form as shown in figure 11. The yield functions and the stress resultants at yield are given in the appendix.
416
A.R.S. Ponter et al.
Figure 12. Convergence for uniform block of figure 4 for Drucker–Prager yield condition with parabolic cone. The sufficient condition for convergence is satisfied.
Figure 13. Complementary strain energy contours for iterations 1–6 for the problem of figure 12.
The matching process has the same characteristics as in the case of an elliptic yield surface, equation (34) repeated here; ε˙ = σ y /E
and
ε˙ v = (py − pL )/Eg(νL ).
(34)
Limit analysis for a general class of yield conditions
417
Figure 14. Variation of load parameter for uniform block for Drucker–Prager yield condition with hyperbolic cone. Non-convergence occurs after the 4th iteration.
Equations (34) are sufficient to define two of the three material parameters (E, νL , pa ). The remaining condition is the requirement that the contour of U = constant which touches the yield surface should otherwise lie outside the yield surface. This condition can be satisfied for both yield surfaces for the corner and the cap by ensuring that the U = constant contour passes through the apex p = −pt and σ = 0. For the cone section U = constant approaches the yield surface for sufficiently large values of pL in the parabolic case, but there exists no values of the moduli so that the contour always lies outside the yield surface in the hyperbolic case. Hence, an isotropic linear material may not always exist which allows convergent solutions for an isotropic yield condition and there are limits to the range of yield conditions which may be used. This is not an inherent restriction on the application of the method as it is possible to overcome this problem by using an anisotropic linear material, but that possibility will not be pursued here. This example demonstrates that the main restriction on the application of the method derives from the condition that the U = constant contour lies outside the yield surface and this condition can be difficult to apply in cases where the curvature of the yield surface is locally small. Hence the method is most appropriate where the yield condition is similar in form to the appropriate U = constant surface. Figures 12 and 14 show the solutions for the uniform block for each of the yield conditions and figures 13 and 15 show U = constant surfaces for the first few iterations in each case. The initial solution is generated by choosing the moduli so that a U = constant contour passes through the two yield points on the p axis. In accordance with the convergence results in section 4, the solution for the parabolic cone converges but the solution for the hyperbolic cone fails to converge and this is borne out by the numerical results. Figure 16 shows the convergence of the plane strain indentor problem of figure 5 for the parabolic cone with the mesh as shown (curve A) and a refined mesh with each element subdivided into 16 elements (curve B). In both cases the changes in the load parameter are confined to the 5th significant figure after 20 iterations
418
A.R.S. Ponter et al.
Figure 15. Complementary strain energy contours for iterations 1–6 for the problem of figure 14. At the 4th iteration the sufficient condition for convergence is violated and the solution fails to converge subsequently.
Figure 16. Convergence of the load parameter for the plain strain indentation problem of figure 5 for the Drucker–Prager yield condition with a parabolic cone section. Curve A corresponds to the mesh shown in figure 5 and curve B to the same mesh where each element has been subdivided into 16 equal elements.
and convergence occurs. Figure 17 shows the variation of the upper bound with iterations for the hyperbolic cone. The sufficient condition for convergence is not satisfied and the solutions fail to converge with increasing iterations.
Limit analysis for a general class of yield conditions
419
Figure 17. Non-convergence of the load parameter for the indentation problem of figure 5 for the Drucker–Prager yield condition with a parabolic cone section.
6. Conclusions The paper discusses the convergence criteria for the method previously discussed by Carter and Ponter (1997a), for a Von Mises yield condition, but here extended to a general yield condition. The method is re-interpreted as a programming method where a linear material is matched to the current strain field and the yield condition in an iterative process and a sufficient condition for convergence is derived. The set of numerical examples provide the following indications. (1) The convergence of the method is most rapid when contours of constant complementary energy for the linear material are close or identical to the yield surface. (2) Although the convergence criterion is a sufficient condition, there is clear evidence that violation of the condition is likely to lead to non-convergence or apparent convergence to a non-optimal value of a load parameter. (3) The method is capable of providing efficient converged solutions for complex yield surfaces although there are cases where an isotropic yield condition may require an anisotropic linear material. (4) The method requires the evaluation of a sequence of linear finite element problems and implimentaion is possible within commercial finite element codes, in this case ABAQUS, by the use of user routines. As a result the method is more assessable than others which require specialist software. The method is applied here to limit analysis problems. The extension to shakedown limits and related solutions requires essentially the same algorithms (Ponter and Engelhardt, 2000), and hence the work of this paper may be regarded as the first step to the generation of such solutions for a general yield surface. Appendix. Yield conditions and stress resultants (σ y , py ) at yield A.1. Cap section BC Yield function:
fc = (p − pa )2 + σ 2 (pb − pa )2 /σ02
1/2
− (pb − pa ) = 0.
(A1)
420
A.R.S. Ponter et al.
Stress resultants at yield:
σ y = σ0 ε˙ / ε˙ 2 +
pb − p a σ0
2
1/2
ε˙ v2
py = pa + (pb − pa ) ε˙ v / ε˙ v2 +
,
σ0 pb − p a
2
1/2
ε˙ 2
.
(A2)
A.2. Parabolic section AB Yield function:
fp = (pa + pt ) Stress resultants at yield:
ε˙ v − ε˙
For
> tan α =
σ σ0
2
− (p + pt ) = 0.
σ0 ; 2(pa + pt )
(A3)
(A4)
ε˙ ε˙ σ02 σ02 , py = −pt + 2(pa + pt ) ε˙ v 4(pa + pt ) ε˙ v ε˙ v σ0 and for 06− 6 tan α = ; 2(pa + pt ) ε˙ σ y = σ0 , p y = pa .
2
σy = −
(A5) (A6) (A7)
A.3. Hyperbolic section BC Yield function: fh = l02 + σ 2 l0 =
where
1/2
− (p tan β − d),
σ02 − (pa + pt )2 tan2 β 2 (pa + pt ) tan β
and
(A8)
d = pt tan β − l0 .
(A9)
Stress resultants at yield
For
ε˙ v − ε˙
σ02 − l02 > tan β, σ02
l0 tan β ε˙ σy = 2 2 , (˙ε tan β − ε˙ v2 )1/2
and for σ y = σ0 ,
1 l0 ε˙ v py = d− 2 2 tan β (˙ε tan β − ε˙ v2 )1/2
(A10)
ε˙ v σ02 − l02 06− 6 tan β σ02 ε˙ p y = pa .
(A11)
Acknowledgements The authors gratefully acknowledge the support during the course of this work of Nuclear Electric Ltd. and the University of Leicester (for ME) and of a European Commission Grant HCM Network Contract ERBCHRXCT940629 (for PF).
Limit analysis for a general class of yield conditions
421
References Collins I.F., Boulbibane, Application of Kinematic Shakedown Theorems to Pavement Design, Proc. Euromech 385, Inelastic Analysis of Structures under Variable Loading: Theory and Applications - Technical University of Aachen, September 1998. Goodall I.W., Goodman A.M., Chell G.C., Ainsworth R.A., Williams J.A., R5: An Assessment Procedure for the High Temperature Response of Structures, Nuclear Electric Ltd., Barnwood, Gloucester, Report, 1991. Hill R., The Mathematical Theory of Plasticity, Oxford University Press, London, 1950. Liu Y.H., Cen Z.Z., Xu B.Y., A numerical method for plastic limit analysis of 3-D structures, Int. J. Sol. Struct. 32 (1995) 1645–1658. Marriott D.L., Evaluation of deformation on load control of stresses under inelastic conditions using elastic finite element stress analysis, ASME Pressure Vessels and Piping Conference, 1985, Pittsburgh, ‘Codes and Standards and Applications for Design and Analysis of Pressure Vessel and Piping Components’, eds: R. Seshadri et al., ASME PVP, 1985, 136, p. 3. Mackenzie D., Boyle T., A method of estimating limit loads by iterative elastic analysis. I-simple examples, Int. J. Pres. Ves. & Piping 53 (1993) 77. Mackenzie D., Nadarajah C., Shi J., Boyle J.T., Simple bounds on limit loads by elastic finite element analysis, ASME J. Press. Ves. Tech. 115 (1993) 27. Nadarajah C., Mackenzie D., Boyle J.T., A method of estimating limit loads by iterative elastic analysis II-nozzle sphere intersections with internal pressure and radial load, Int. J. Pres. Ves & Piping 53 (1993) 97. Ponter A.R.S., Carter K.F., Limit state solutions, based uponlinear solutions with a spatially varying elastic modulus, Meth, Appl. Mech. Eng. 140 (1997a) 237–258. Ponter A.R.S., Carter K.F., Shakedown state simulation techniques based upon linear elastic solutions, Meth. Appl. Mech. Eng. 140 (1997b) 259– 279. Ponter A.R.S., Engelhardt M., Shakedown limits for a general yield condition: implementation and application for a Von Mises yield condition, Eur. J. Mech. A/Solids 19 (2000) 423–445 (this issue). Seshadri R., Fernando C.P.D., Limit loads of mechanical components and structures based on linear elastic solutions using the GLOSS R-Node method, Trans ASME PVP, Vol 210-2, San Diego, 1991, pp. 125–134. Shi J., Mackenzie D., Boyle J.T., A method of estimating limit loads by iterative elastic analysis. I-simple examples, Int. J. Pres. Ves. & Piping 53 (1993) 77. Zhang Y.G., Zhang P., Lu M.W., Computational limit analysis of rigid plastic bodies in plane strain, Acta Mech. Sol. Sinica 6 (1993) 341–348.