- Email: [email protected]
Large displacement elastoplastic analysis of frames using an iterative LCP approach
Int. J. Mech. Sci. Vol. 33, No. 5, pp. 379 391, 1991
0020-7403/91 $3.00 + .00 ,~, 1991 Pergamon Press plc
Printed in Great Britain.
LARGE DISPLACEMENT ELASTOPLASTIC ANALYSIS OF FRAMES USING AN ITERATIVE LCP APPROACH R. R. WAKEFIELDand F. TIN-LOl School of Civil Engineering, Univerisity of New South Wales, P.O. Box 1, Kensington, N.S.W. 2033, Australia (Received 23 May 1990; and in revised form 15 October 1990)
Abstract--A method for the analysis of elastoplastic structures with large displacementsis presented in this paper. A parametric linear complementarity mathematical programming formulation, in combination with an iterative incremental approach, is used to solve the elastoplastic large displacement problem for assumed piecewise linearized yield conditions. Details of the computational procedure are described along with the steps required for implementation. Several well known examples are solved to demonstrate the accuracy of the method. The influence of different piecewise linearizations of the yield locus on the results is also discussed.
NOTATION C e F K~ KL Kr M M Mp N N Np q Q Qe Qp R S t u Z {~} {fl} {6} ~. t~ {a}
compatibility matrix elasticstrain vector applied force vector geometric stiffness matrix linear stiffnessmatrix = KC + KL = NtZN matrix bending moment plasticmoment capacity matrix of unit outward normals to piecewise linear yield surface axial force axial plastic capacity total generalized strain vector total generalized stress vector elasticstress vector plasticstress vector vectorof perpendicular distances of yield hyperplanes from the origin of stress space unassembledstructure stiffnessmatrix load multiplier nodal displacement vector matrix of influence coefficients subset of rows of LCP that have not been previously active subset of rows of LCP that are inactive but have been previously active subset of active rows of the LCP vectorof plastic multipliers vectorof plastic potentials subset of rows that reach the yield surface in a load step
INTRODUCTION Recent work by the a u t h o r s [1] has p r o d u c e d a m e t h o d for small displacement large scale n o n h o l o n o m i c (irreversible) elastoplastic analysis of structures using a linear complem e n t a r i t y formulation. The present paper is the first step in extending that work to the large scale elastoplastic analysis of structures to a c c o u n t for large displacements. The large displacement b e h a v i o u r of frames has received considerable a t t e n t i o n in the past. M a n y f o r m u l a t i o n s have been suggested for dealing with structures that are both materially n o n - l i n e a r a n d geometrically non-linear. Papers of interest include those of B a c k l u n d [2], W o n g a n d T i n - L o i [3], Kassimali [4] a n d Meek a n d L o g a n a t h a n [5] who solve the p r o b l e m from a traditional finite element approach; a n d those of C o r r a d i [6], C o n t r o et al. [7], F r a n c h i a n d C o h n [8], Maier a n d C o n t r o [9], De Freitas a n d Smith [10] 379
380
R.R. WAKEFIELDand F. T1N-Lo!
and Contro and Maier [ l l ] who use a mathematical programming approach. In their paper, De Freitas and Smith [10] give a unifying treatment of the formulations that have appeared in the literature. Of the mathematical programming formulations mentioned, only [11] has demonstrated the ability to solve truly large scale problems in the area of cable structures. This paper presents an alternative mathematical programming based incremental method which utilizes the robust and efficient large scale procedures detailed in Ref. [l]. It has therefore the potential to perform large elastoplastic analyses in the large displacement regime. A class of discretized skeletal structures made up of elastic-perfectly plastic components is considered. The non-linear yield functions are represented by approximate piecewise linear functions in generalized bending moment-axial force stress space. The materially non-linear behaviour of the structural components, under a small strain assumption, is handled by a parametric linear complementarity problem (PLCP) formulation, which has been used by several researchers to solve geometrically linear problems [12 15]. The large displacement behaviour of the frame is tackled by assuming that the response is geometrically linear within each incremental load step; a Newton-Raphson scheme is then used to ensure that joint equilibrium is preserved. The conventional geometric stiffness matrix is included in the stiffness evaluation. It is interesting to note that this matrix can be derived from the same minimum energy principle given for instance in Refs [-9, 1 I, 16]. Details of the geometrically linear analysis within each load step are presented, in particular, the equilibrium relations, the P L C P solution of the geometrically linear problem and the piecewise proportional loading scheme used are discussed. The computational procedure is then described along with the steps required for implementation. Operation of the method is illustrated in the first instance with an example taken from the literature, and the effect of the yield surface discretization adopted on the results produced is discussed in some detail. Several other well known examples are also analysed and the results compared with existing ones in order to assess the accuracy of the proposed method of analysis.
METHOD OF ANALYSIS
Preliminaries The basic assumptions made are: (a) The structure is plane and made up of straight prismatic elements. (b) Plastic deformations are lumped at the ends of each element while the element itself is considered to be elastic. (c) A generalized plastic hinge model is assumed valid. (d) Plastic deformations are consistent with the normality rule of the classical theory of plasticity. (e) The effect of strain-hardening on structure behaviour is not considered. (f) Deformations can be large but strains are assumed to be small. (g) All elements of the structure are assumed to have a form factor of unity. The beam model used in the analysis can be viewed as an elastic element with plastic deformations lumped at each end. The analysis is based on linear incremental techniques and within each load step the structure is assumed to behave in a geometrically linear manner. The material non-linearity (plasticity) is accounted for in each load step by using the P L C P formulation. Adjustments are made to the geometry at the end of each load step and the procedure is repeated iteratively until equilibrium is achieved. Each of the geometrically linear relations will now be described.
Equilibrium relations In the present work, "natural" generalized stresses [17] are used to describe the stress state of each element. The deformation pattern of the structure is then prescribed to depend upon a vector of q of generalized natural strains which are in turn defined by u the nodal displacement vector. The structural behaviour can then be described by the relation between q and the natural generalized stress vector Q. For the plane frames considered,
Large displacement elastoplastic analysis of frames
381
three natural stress resultants are required for each member, namely, the bending moments at each end and the axial force in the member. The equilibrium between the elements of the structure can be described in terms of the total, elastic and plastic strains represented by vectors q, e and p, respectively, the total stresses Q, the linear elastic response Qe to the external action F and the stresses QP corresponding to the plastic strains. At any instant, the above vectors are connected by the following relations: q = e + p
Q=Qe
+ Qp
CtQ e + CtQ p = F
(la)
(lb) (lc)
where C is the compatibility matrix for the structure and the superscript t represents transpose.
Mathematical programmin9 formulation In order to describe the mathematical programming formulation we need to elaborate on some of the relationships used. A piecewise linearized representation of the yield locus is adopted and the yield condition for the whole structure can be expressed in terms of hyperplanes [17] as follows: d0=N tQ-R~0
(2)
where dp is the vector of plastic potentials, N is the matrix of unit outward normals to the yield hyperplanes, R is the plastic capacity vector and 0 is the null vector. We assume that the material is elastic-perfectly plastic and that the plastic strain vector p can be represented as p = Nk (3) where k is a vector of non-negative plastic multipliers. Loading is assumed to be piecewise proportional. Interpreting the stresses as the sum of the linear-elastic response Q¢ to the loads and the self equilibrated stresses caused by the plastic strains [17], we have for each load multiplier t: Q = tQ ~ + Zp
(4)
where Z is a matrix of influence coefficients. The following relations then describe the nonholonomic (irreversible) elastoplastic problem: - ~ = R - tNtQ e - NtZNk
(5a)
- ~/> 0
(Sb)
~,/> 0
(5c)
dpt ~. ~> 0
(5d)
with the initial condition k = 0, t = 0; a superimposed dot represents rate. Problems described by these relations have been considered by m a n y researchers, with numerous references given in the reviews of Maier and Munro [18] and Maier and Smith [19]. The corresponding holonomic relations (non-linear elastic, reversible) are: - ~b = R - tNtQ ~ - N t Z N k
(6a)
- dp ~> 0
(6b)
k t> 0
(6c)
d~tk = 0
(6d)
Relations (6) differ from those in (5) in that Eqns (6c) and (6d) allow the value of 2 to decrease, an occurrence which cannot take place in the nonholonomic system. Also, at any stress point in the holonomic system, plastic deformation must become zero before elastic MS
33: 5 - E
382
R.R. WAKEFIELDand F. T1N-LOl
unloading can occur. We make use of both the nonholonomic and holonomic conditions in our solution procedure. The above equations fall into a special class of mathematical programming problems known as linear complementarity problems (LCP), which have received considerable attention in operations research; for example, see the paper by Cottle and Dantzig [20]. The system of equations represented by (6) can be solved parametrically for a given value of t to generate a solution ~.(t). Corresponding values of stresses and strains for the load multiplier t can then be determined using the following relationships: Q = tQ e + ZNk(t)
(Ta)
and q=S-IQ+p = t S - 1 Q e + ( S - 1 Z N + N)~,(t)
(7b)
where S is the unassembled structure stiffness matrix and an inverse is represented by the superscript - 1. The complete history of stresses and strains can thus be obtained. There are numerous methods for solving the LCP, e.g. Refs [13, 14, 21, 22, 23]. In this work, we use an elegant and direct method of solution proposed by Kaneko [14]. Before describing the solution procedure, some comments need to be made on the matrices and vectors involved. Namely, R ~> 0 and the matrix M = - N tZN is positive semi-definite. In the description of the procedure for solution we will adopt Kaneko's notation [ 14] and refer to a general iteration tableau:
0) where {3} is the subset of active rows of the system (i.e. those that have reached the yield surface and deformed plastically), {~} is the subset of rows that have not previously been active, and {/~} is the subset of inactive rows that have been active previously (i.e. those that have deformed plastically and then unloaded). The union {c~ w /~} then makes up the subset of inactive rows. Quantities k, NtQ e and 1VI represent the current values of the variables. The analysis proceeds as follows:
Step 1 Set up the initial tableau. Initialize the variables.
Step 2 If all - N tQe /> 0 then stop.
Step 3 Determine the minimum value of the load multiplier to cause any row of the tableau to be equal to zero, that is, for any element in the structure to reach yield. The minimum value is the next load step t* given by:
t* = m i n
i~a N]Q~
.. N - 4j Q je
> 0 .
(9)
Step 4 Work out {a}, the subset of rows that reach the yield surface as a result of the load step t*; {a} is the subset of {~ w/3} such that j ~ {a} if: Rj
+ t*NtQ *
+ l~/lji~,~ -~- O.
(10)
Large displacementelastoplastic analysis of frames
383
Form the following LCP sub-problem and solve: ~,~ = - NtQ,~ + 1Vl~d~ + 1('1~,~
(lla)
~b,, =
(11 b)
- N t Q,~ + l£'l,,,s dp,s + lgl,,,, ~,,,
~b >/0, k ~> 0, ~btk = 0.
(llc)
Use Graves' algorithm [24] to solve the LCP. It is able to solve the LCP if the M matrix is positive semi-definite, and identifies a solution ray if a solution exists [25]. Nonuniqueness of solution can also be detected here [26]. If the LCP has no solution then the structure has reached an unstable condition such as collapse or a limit point. This is indicated by the fact that the matrix M is no longer positive semi-definite. The analysis is then terminated. If a solution is found then use the pivots that gave a solution of the LCP sub-problem to perform pivots on the entire tableau. Then return to step 2. A detailed description of this algorithm can be found in Kaneko's original paper [14]. As can be seen from the above, similar with other schemes, the procedure automatically calculates t for which the next plastic event (activation of a yield plane or unloading) occurs. Of course the analysis can be carried out for a specified maximum value of t; this is useful in considering piecewise linear proportional loading stages. It must be mentioned that two computational problems arise in the implementation of Kaneko's method to solve large scale problems. In particular, the large size of the matrices involved and the resulting storage problems, and the propagation of roundoff errors and their effects on the accuracy of solutions generated. These problems have been addressed by the authors and details of their solution have been discussed in another paper [1].
Piecewise proportional loadiny and its application to geometrically non-linear analysis As mentioned in the previous section, the mathematical programming formulation can be used with piecewise proportional load paths. This is in fact utilized in the present method to keep account of the previous load history of the structure when the geometry is updated. Such a solution procedure, based on the theory of "stepwise holonomic" elastoplasticity with piecewise linearized yield conditions, was developed by Maier and his co-workers, eg. Refs. [12, 15, 27]. If we consider the structure at load level F i, the vector R in Eqn (6a) contains all the information about the previous structural behaviour, including whether any points on the structure have reached yield and how much further elastic deformation can take place before yielding occurs [8]. If we then apply a load increment AF to the structure, and perform the analysis, the pivoting within the P L C P formulation will occur at rows in the tableau corresponding to previously yielded nodes (i.e. those with an R value equal to zero). This has the effect of modifying the elastic problem that was set up at the start of the increment to include the effect of previous plastic and elastic deformations. Therefore, if we consider each new load increment as a new proportional load stage and use the value o f ~ at the end of the previous load stage for R in the next load stage, previous loading can be taken into account. This can be used to advantage in the geometrically non-linear analysis since the vector R can retain the state of the structure when geometric updating is undertaken.
Computational technique The non-linear response of the structure is generated by using an incremental load approach with a Newton-Raphson type iteration performed at each load level to satisfy joint equilibrium equations. If we consider a load level F i and assume that the corresponding deformed configuration of the system may be described symbolically as u i, then it is necessary for us to determine the configuration I! i + 1 corresponding to load level F i+l F ~ + AF in which AF represents a small increment in load. If we assume that the small load increment represents a new load step applied to the existing structure, we can proceed as follows: =
384
R.R. WAKEFIELDand F. TIN-LOl
Step 1 Calculate F i+1 = F i + AF, the new load level.
Step 2 Calculate the stiffness matrix K T of the structure as the sum of the linear stiffness matrix K L of the structure in its present configuration and the second order geometric stiffness matrix K~ [-28].
Step 3 Use the stiffness matrix and AF to solve for Au e the incremental elastic displacements and subsequently for AQ e the incremental elastic stresses.
Step 4 Use the P L C P method to solve for the increments in the elastic stresses and the plastic multipliers. Update the plastic multipliers as required and calculate the total displacements. Calculate the value of d~.
Step 5 Determine the new structure geometry, and calculate the new member lengths and direction cosines.
Step 6 Calculate the member end forces in local coordinates and, using the updated geometric information, work out the member end forces in global coordinates.
Step 7 Calculate the unbalanced joint forces by subtracting the member end forces from the applied nodal forces.
Step 8 Compare the unbalanced joint forces with a prescribed convergence criterion, for example, a force norm [29].
Step 9 Set the new value of the vector R equal to the value of - 4.
Step 10 If convergence has been achieved, return to Step 1 and apply the next load increment. If convergence has not been achieved, that is, the unbalanced force vector is larger than acceptable, use the unbalanced force vector as the new load increment and continue.
Step 11 Calculate the stiffness matrix of the structure. Solve for the elastic displacements and stresses and set up the PLCP.
Step 12 Perform a holonomic analysis. If all load multipliers are less than zero, reverse the out of balance forces (the force vector), solve the LCP, and reverse the results when updating 2 values, before doing subsequent stress calculations.
Step 13 Determine the new structure geometry and calculate updated member lengths and direction cosines.
Step 14 Return to Step 6. Before illustrating the operation of the algorithm on an example the following remark is appropriate. In the present scheme which combines a classical N e w t o n - R a p h s o n approach with a mathematical programming solution technique, the matrix M needs to be updated at every iteration. This appears to result in a computational burden of similar magnitude to that reported for the nonlinear programming approach given in [11]. In the latter direct energy minimization scheme, an iterative method is included in the nonlinear program
Large displacement elastoplastic analysis of frames
385
solution procedure. Direct comparison of the two approaches without detailed numerical testing on similar structures is difficult. ILLUSTRATIVE EXAMPLE
To illustrate the performance of the method, the analysis of a simply supported beam under uniform loading will be discussed. The beam is illustrated in Fig. l(a) and has previously been analysed by Backlund I-2] who used a finite element approach with hybrid Load Intensity w (MNIm)
1
~ pL
~ pL
i pL
~ pL
~ pL ~1
l =
0.5m
=
I-
-IIJ2
-IL
-I-
-I-
L
(a)
L
-IL
-I L/2
(b)
1.0
1.0
M/Mp
• , .
M/Mp
N/Np
1.0
N/Np
(c)
1.0
(d)
2
/,,"" °°" I
•,°
/
Z
7
I "~
~
"•~•
/.&
.o~
S ~
-
ElasticDeflection
.... ~---
LB Yield Locus
- - ~- -
UB Yield Locus
x
Backlund1976
0
0.00
0.01 Central Deflection (rn)
0.02
(e)
FIG. 1. Simply supported beam example: (a) beam layout; (b) half beam model; (c) lower bound yield locus; (d) upper bound yield locus; (e) load-deflection curves.
386
R.R. WAKEFIELDand F. TIN-LOl
beam elements. Contro et al. [7-] analysed the same structure both as being linear elastic and elastic-perfectly plastic. In the latter case, a non-linear programming approach was used. The beam is of uniform rectangular cross-section of 1 m width, 0.1 m depth, and made from a material with a Young's modulus of 220 x l09 N / m 2 and a yield stress of 300 x l06 N / m 2. The beam is subjected to a uniform load w. The non-linear yield locus, in terms of bending moment (M) and axial force (N), at any node is expressed by the equation +
- 1= 0
(12)
where Mp is the plastic capacity in pure bending and N o is the plastic capacity under pure axial load. Two piecewise linear approximations of this yield surface are used in our example to demonstrate the effect of linearization on the load deflection behaviour of the structure. The
W/8
15.984
,w/8
J w/8
1 w/8 .
I
¥~
T
I W/8
in~ L/16
I____1 I
I
L/8
L/8
I
L/8
I_
1_/16
L/2
__J
__1___1_
2-'--'--'--'--'-
2
I
L = 47.992 in
(a) 0.6" x
0.5"
0.4"
0.3"
~ 02"
01 /
x
Experimental (Maj=d, 1972)
00 0
1
Vertical Deflection of Centre Node (inches)
(b) FIG. 2. Pitched roof portal frame example: (a) dimensions and loading; (b) load-deflection curves.
Large displacement elastoplastic analysis of frames
387
approximations are shown in Figs l(c) and (d). Figure l(c) represents a lower bound approximation and Fig. l(d) represents an upper bound approximation. The structure is modelled using six beam elements as shown in Fig. l(b). The perfectly elastic structure was first analysed and the results, as shown in Fig. l(e), were found to be in good agreement with those of Backlund [2]. The results of the elastoplastic analyses with lower and upper bound polygons clearly show the influence of the yield surface discretization. The load-deflection response for both yield locus approximations give somewhat different results [Fig. l(e)]. The incidence of first yield is different in both cases. This is because the distance of the yield hyperplane from the origin, in the direction of the stress path taken during loading, is different. It is clear from Fig. l(d) that when bending moment is significant, the distance of the yield line from the origin is further in the upper bound case so that a greater load is required to cause first yield. As the load level is increased, the 4W W
4W
4W
==
r_~
L 4W
-I
L (a) 0.6
Ig 05"
04"
03"
_j
02Analysis
/ 0.1
0.0 000
•
Home,1963
x
De F r e i t a s
/ and Smith
i
i
i
0.02
0.04
0.06
1985
0.08
u/L
(b) FIG. 3.
Two storey frame example: (a) dimensions and loading; (b) load-deflection curves.
388
R . R . WAKEFIELD a n d F. TIN-Lol
linearization has a further effect through the normality rule of plasticity. The direction of the outward normal to the yield curve influences the plastic deformation taking place [7]. This indicates that in order to accurately model the deformation at any stress point, which is at yield, it is necessary to have the outward normal to the linearized yield locus at that point in a direction similar to the outward normal to the non-linear yield locus at that stress point. This is apparent in the cases illustrated [Fig. l(c)]. The upper bound yield locus models the deformation more accurately as the axial force in the beam becomes larger. This is because the upper bound yield line is tangent to the true yield locus at N = N p = 1 (i.e. the slopes of the normals to both curves are similar). When the lower bound yield locus is used, the deformation is more accurately modelled at lower load levels when both moment and axial force are significant. The problem of accurately representing the outward normal to the yield locus can be overcome by refining the piecewise linearization with particular consideration given to regions of the yield locus where there are large changes in the direction of the outward normal. FURTHER EXAMPLES The method presented has been used to analyse several frames that have appeared in the literature. These examples give an indication of the capability and accuracy of the method. The single bay pitched roof portal frame of Fig. 2(a) was experimentally tested by Majid [30] under the loading shown. The frame was modelled as 12 elements with the following
i2 D
L = 19.5 in
8 elements of equal length
(a)
(b)
+0
j
400
j
/
+.
'
-o m
A~
~
200
Elastic
o"
,c ,/.....
o
o 4p
.... t---
LB Yield Locus
- -,A- -
UB YieldLocus
o 0 0.0
ExperimenlalValues
i
i
i
i
0.2
0.4
0.6
08
1.0
Central Deflection (inches)
(c) FXG. 4. Fully fixed beam example: (a) dimensions and loading; (b) ball" beam model; (c) load deflection curves.
Large displacementelastoplasticanalysisof flames
389
properties: area 0.262in 2, moment of inertia 0.0057in 4, modulus of elasticity 28076 kips/in 2, yield stress 40.61 kips/in 2, plastic moment 1.36 kips-in and axial plastic capacity 10.64 kips. The flame has also been analysed by other researchers [3, 4]. In the present analyses, a lower bound yield locus [Fig. l(c)] was used. The obtained load-deflection diagram for the vertical deflection of the centre node is given in Fig. 2(b), along with experimental results reported by Majid [30]. The theoretical results are very similar to those reported in the literature [3, 4]. The second example is a two storey single bay flame first analysed by Home [31] and later analysed by Jennings and Majid [32], Corradi [6], De Freitas and Smith [10] and Wong and Tin-Loi [3]. The frame is shown in Fig. 3(a); member details were as follows: modulus of elasticity 30 x 1 0 6 lb/in 2, moment of inertia 51.2 in 4, area 3.2 in 2, length (L) 400 in, plastic moment 464 kip-in and axial plastic capacity 115 kips. The yield surface used in our analysis was the lower bound shown in Fig. l(c). Results of the analysis are presented in the form of a load-deflection curve, along with the results reported by Horne [31] and De Freitas and Smith [10]. Our results are in good agreement with reported ones. The final example is a fully fixed beam and subjected to lateral loading, as shown in Fig. 4(a). The beam, tested by Campbell and Charlton [33], was modelled using eight elements of equal lengths to represent half the beam [Fig. 4(b)]. The elements had the following properties: breadth 0.310 in, depth 0.304 in, span 19.5 in, axial stiffness 2.695 x 105 lb, flexural stiffness 20760 lb in 2, y!eld stress 16.1 tons/in 2, plastic moment 258 lb in and axial plastic capacity 3395 lb. Both the upper bound and lower bound yield surfaces shown in Figs. l(c) and (d) were used in the analysis. The load-deflection curves are given in Fig. 4(c) along with the experimental results of Campbell and Charlton [33]. The curves demonstrate similar behaviour to that found in the illustrative example. CONCLUDING REMARKS An incremental mathematical programming method for large displacement analysis, suitable for large structures, has been presented. The results of analyses of several well known examples from the literature demonstrate the accuracy of the method. The following remarks should be made: (a) In some cases the piecewise linearization of the yield locus can have a significant effect on the results obtained as was demonstrated in two examples described in this paper. An attempt should therefore be made to make the slope of the outward normal of the linearized yield surface close to the slope of the outward normal to the non-linear yield locus, in regions of the yield locus which are likely to be activated. (b) In the case of structures that exhibit a falling branch response after reaching a limit point, the method is only capable of tracing the load-deflection path up to that limit point. (c) While the solution of large scale problems has not been demonstrated explicitly in the examples given here, the method presented essentially only involves repeated analyses of the type detailed in Ref. [1]. One of the examples reported therein was that of the small displacement elastoplastic analysis under combined stresses of a 12 storey two-dimensional frame with 252 active degrees of freedom. The techniques for handling the large matrices involved and for controlling the propagation of round-off errors will still be effective for the present combined geometrically and materially nonlinear analysis. (d) The extension of the analysis to include higher dimensional yield surfaces, that is greater numbers of generalized stresses, is straightforward in terms of Eqn (1). (e) The number of incremental load steps needed in the analysis could be reduced by using a more accurate stiffness formulation, e.g. Refs [34, 35] than that used in the present method. Further work is at present being undertaken by the authors to widen the capabilities of the method to perform analysis of three-dimensional spatial structures. The basic formulation remains unaltered. However, the non-vectorial behaviour of large rotations in threedimensional space will need to be taken into account along with the effects of six generalized
390
R.R. WAKEFIELD and F. TIN-LOl
stress components. The more accurate stiffness matrices referred to in (e) above will also be implemented. Acknowledyements The authors would like to thank the reviewers for their constructive comments on an earlier version of this paper. REFERENCES I. R. R. WAKEFIELD and F. TIN-Lol, Large scale nonholonomic elastoplastic analysis using a linear complementarity formulation. Comput. Meth. Appl. Mech. Engn9 84, 229 (1990). 2. J. BACKLUND,Large deflection analysis of elasto-plastic beams and frames. Int. J. Mech. Sci. 18, 269 (1976). 3. M. B. WONG and F. TJN-LOl, Analysis of frames involving geometrical and material nonlinearitics. Comput. Struct. 34, 641 (1990). 4. A. KASSlMAL|, Large deformation analysis of elastic-plastic frames. J. Struct. Engng, ASCE 109, 1869 (1983). 5. J. L. MEEK and S. LOGANATHAN, Geometrically non-linear behaviour of space frame structures. Comput. Struct. 31, 35 (1989). 6. L. CORRADI, Large displacement elastic-plastic analysis, in Engineeriny Plasticity by Mathematical Pro.qrammin9 (edited by M. Z. Cohn and G. Maier). Pergamon Press, Oxford (1979). 7. R. CONTRO, L. CORRADI and R. NOVA, Large displacement analysis of elastoplastic structures a nonlinear programming approach. Solid Mech. Arch. 2, 433 (1977). 8. A. FRANCHIand M. Z. Cony, Strain-softening and large displacement analysis in structural plasticity. Comput. Struct. 11, 421 (1989). 9. G. MAIER and R. CONTRO, Energy approach to inelastic cable structure analysis. J. En.qn9 Mech. Div., ASCE 101,531 (19751. 10. J. A. T. DE FREITAS and D. L. SMITH, Elastoplastic analysis of planar structures for large displacements. J. Struct. Mech. 12, 419 (1984-85). 1I. R. CONTRO and G. MMER, Inelastic analysis of pretensioned cable structures by mathematical programming, in Steel Structures: Recent Research Advances and their Applications to Desiqn (edited by M. N. Pavlovic), pp. 85--105. Elsevier Applied Science, London (1986). 12. O. DE DONATO and G. MAIER, Historical deformation analysis of elastoplastic structures as a parametric linear complementarity problem. Meccanica 3, 166 (1976). 13. G. MAIER, S. GIACOMIN~ and F. PATERLIM, Combined elastoplastic and limit analysis via restricted basis linear programming. Comput. Meth. Appl. Mech. Engn.q 19, 21 (1979). 14. I. KANEKO. Complete solutions for a class of elastic-plastic structures. Comput. Meth. Appl. Mech. Enqno 21, 193 (1980). 15. O. DE DONATO and G. MAIER, Mathematical programming methods for the inelastic analysis of reinforced concrete frames allowing for limited rotation capacity. Int. J. Num. Meth. Engn9 4, 307 (1972). 16. L. CORRADI and G. MMER, Extremum theorems for large displacement analysis of discrete elastoplastic structures with piecewise linear yield surfaces. J. Optim. Theory Applic. 15, 50 (1975). 17. G. MINER, A matrix structural theory of piecewise-linear elastoplasticity with interacting yield planes. Meccanica 5, 54 (1970). 18. G. MA1ER and J. MUNRO, Mathematical programming applications to engineering plastic analysis. Appl. Mech. Rev. 35, 1631 (1982). 19. G. MAIER and D. L. SMITH, Update to "Mathematical programming applications to engineering plastic analysis" [AMR 35, 1631 (1982)]. Applied Mechanics Update, ASME (1986). 20. R. W. COTTLE and G. B. DANTZIG, Complementary pivot theory and mathematical programming. Linear Al,qebra Appl. 1, 103 (1968). 21. D. L. SMITH, The Wolfe Markowitz algorithm for non-holonomic elastoplastic analysis. Enonq Struct. 1, 9 (1978). 22. F. CAS('IATI, Elastic-plastic deformation analysis: a parametric linear programming method. J. Mec. Appl. 2, 269 (1978). 23. O. DE DONATO and G. MAIER, Local unloading in piecewise-linear plasticity. J. Engng. Mech. Div., ASCE 102, 383 (1976). 24. R. L. GRAVES,A principal pivoting simplex algorithm for linear and quadratic programming. Operations Res. 15. 482 (1967). 25. G. MAIER, Incremental plastic analysis in the presence of large displacements and physical instabilizing effects. Int. J. Solids Structs 7, 345 (1971). 26. R. R. WAKEFIELD and F. TIN-LOI, Mathematical programming and uniqueness in nonholonomic plasticity. Comput. Struct. 34, 477 (1990). 27. O. DF. DONATO and G. MAIER, Finite element elastoplastic analysis by quadratic programming: the multistage method, 2rid Int. Conf. on Struct. Mech. in Reactor Technology, 10-14 September, Paper M2/8, 1--12 (1973). 28. H. C. MARTIN, Finite elements and the analysis of geometrically non-linear problems. Proceedings, Japan-US Seminar on Matrix Methods of Structural Analysis and Design, 25-30 August, Tokyo, Japan (1969). 29. P. G. BERGAN and R. W. CLOUGH, Convergence criteria for iterative processes. AIAA J. 10, 1107 (1972). 30. K. I. MAJID, Non-Linear Structures. John Wiley, New York (1972). 31. M. R. HORNE, Elastic-plastic failure loads of plane frames. Proc. Roy. Sot'. London 274, 343 (1963).
Large displacement elastoplastic analysis of frames
391
32. A. JENNINGS and K. MAJID, An elastic-plastic analysis by computer for framed structures loaded up to collapse. The Structural Engineer 43, 407 (1965). 33. T. I. CAMPBELL and T. M. CHARLTON,Finite deformation of a fully fixed beam comprised of a non-linear material. Int. J. Mech. Sci. 15, 415 (1973). 34. C. ORAN,Tangent stiffness in space frames. J. Struct. Div., ASCE 99, 987 (1973). 35. J. L. MEEK and H. S. TAN, Geometrically non-linear analysis of space frames by an incremental iterative technique. Comput. Meth. Appl. Mech. Engng 47, 261 (1984).
0020-7403/91 $3.00 + .00 ,~, 1991 Pergamon Press plc
Printed in Great Britain.
LARGE DISPLACEMENT ELASTOPLASTIC ANALYSIS OF FRAMES USING AN ITERATIVE LCP APPROACH R. R. WAKEFIELDand F. TIN-LOl School of Civil Engineering, Univerisity of New South Wales, P.O. Box 1, Kensington, N.S.W. 2033, Australia (Received 23 May 1990; and in revised form 15 October 1990)
Abstract--A method for the analysis of elastoplastic structures with large displacementsis presented in this paper. A parametric linear complementarity mathematical programming formulation, in combination with an iterative incremental approach, is used to solve the elastoplastic large displacement problem for assumed piecewise linearized yield conditions. Details of the computational procedure are described along with the steps required for implementation. Several well known examples are solved to demonstrate the accuracy of the method. The influence of different piecewise linearizations of the yield locus on the results is also discussed.
NOTATION C e F K~ KL Kr M M Mp N N Np q Q Qe Qp R S t u Z {~} {fl} {6} ~. t~ {a}
compatibility matrix elasticstrain vector applied force vector geometric stiffness matrix linear stiffnessmatrix = KC + KL = NtZN matrix bending moment plasticmoment capacity matrix of unit outward normals to piecewise linear yield surface axial force axial plastic capacity total generalized strain vector total generalized stress vector elasticstress vector plasticstress vector vectorof perpendicular distances of yield hyperplanes from the origin of stress space unassembledstructure stiffnessmatrix load multiplier nodal displacement vector matrix of influence coefficients subset of rows of LCP that have not been previously active subset of rows of LCP that are inactive but have been previously active subset of active rows of the LCP vectorof plastic multipliers vectorof plastic potentials subset of rows that reach the yield surface in a load step
INTRODUCTION Recent work by the a u t h o r s [1] has p r o d u c e d a m e t h o d for small displacement large scale n o n h o l o n o m i c (irreversible) elastoplastic analysis of structures using a linear complem e n t a r i t y formulation. The present paper is the first step in extending that work to the large scale elastoplastic analysis of structures to a c c o u n t for large displacements. The large displacement b e h a v i o u r of frames has received considerable a t t e n t i o n in the past. M a n y f o r m u l a t i o n s have been suggested for dealing with structures that are both materially n o n - l i n e a r a n d geometrically non-linear. Papers of interest include those of B a c k l u n d [2], W o n g a n d T i n - L o i [3], Kassimali [4] a n d Meek a n d L o g a n a t h a n [5] who solve the p r o b l e m from a traditional finite element approach; a n d those of C o r r a d i [6], C o n t r o et al. [7], F r a n c h i a n d C o h n [8], Maier a n d C o n t r o [9], De Freitas a n d Smith [10] 379
380
R.R. WAKEFIELDand F. T1N-Lo!
and Contro and Maier [ l l ] who use a mathematical programming approach. In their paper, De Freitas and Smith [10] give a unifying treatment of the formulations that have appeared in the literature. Of the mathematical programming formulations mentioned, only [11] has demonstrated the ability to solve truly large scale problems in the area of cable structures. This paper presents an alternative mathematical programming based incremental method which utilizes the robust and efficient large scale procedures detailed in Ref. [l]. It has therefore the potential to perform large elastoplastic analyses in the large displacement regime. A class of discretized skeletal structures made up of elastic-perfectly plastic components is considered. The non-linear yield functions are represented by approximate piecewise linear functions in generalized bending moment-axial force stress space. The materially non-linear behaviour of the structural components, under a small strain assumption, is handled by a parametric linear complementarity problem (PLCP) formulation, which has been used by several researchers to solve geometrically linear problems [12 15]. The large displacement behaviour of the frame is tackled by assuming that the response is geometrically linear within each incremental load step; a Newton-Raphson scheme is then used to ensure that joint equilibrium is preserved. The conventional geometric stiffness matrix is included in the stiffness evaluation. It is interesting to note that this matrix can be derived from the same minimum energy principle given for instance in Refs [-9, 1 I, 16]. Details of the geometrically linear analysis within each load step are presented, in particular, the equilibrium relations, the P L C P solution of the geometrically linear problem and the piecewise proportional loading scheme used are discussed. The computational procedure is then described along with the steps required for implementation. Operation of the method is illustrated in the first instance with an example taken from the literature, and the effect of the yield surface discretization adopted on the results produced is discussed in some detail. Several other well known examples are also analysed and the results compared with existing ones in order to assess the accuracy of the proposed method of analysis.
METHOD OF ANALYSIS
Preliminaries The basic assumptions made are: (a) The structure is plane and made up of straight prismatic elements. (b) Plastic deformations are lumped at the ends of each element while the element itself is considered to be elastic. (c) A generalized plastic hinge model is assumed valid. (d) Plastic deformations are consistent with the normality rule of the classical theory of plasticity. (e) The effect of strain-hardening on structure behaviour is not considered. (f) Deformations can be large but strains are assumed to be small. (g) All elements of the structure are assumed to have a form factor of unity. The beam model used in the analysis can be viewed as an elastic element with plastic deformations lumped at each end. The analysis is based on linear incremental techniques and within each load step the structure is assumed to behave in a geometrically linear manner. The material non-linearity (plasticity) is accounted for in each load step by using the P L C P formulation. Adjustments are made to the geometry at the end of each load step and the procedure is repeated iteratively until equilibrium is achieved. Each of the geometrically linear relations will now be described.
Equilibrium relations In the present work, "natural" generalized stresses [17] are used to describe the stress state of each element. The deformation pattern of the structure is then prescribed to depend upon a vector of q of generalized natural strains which are in turn defined by u the nodal displacement vector. The structural behaviour can then be described by the relation between q and the natural generalized stress vector Q. For the plane frames considered,
Large displacement elastoplastic analysis of frames
381
three natural stress resultants are required for each member, namely, the bending moments at each end and the axial force in the member. The equilibrium between the elements of the structure can be described in terms of the total, elastic and plastic strains represented by vectors q, e and p, respectively, the total stresses Q, the linear elastic response Qe to the external action F and the stresses QP corresponding to the plastic strains. At any instant, the above vectors are connected by the following relations: q = e + p
Q=Qe
+ Qp
CtQ e + CtQ p = F
(la)
(lb) (lc)
where C is the compatibility matrix for the structure and the superscript t represents transpose.
Mathematical programmin9 formulation In order to describe the mathematical programming formulation we need to elaborate on some of the relationships used. A piecewise linearized representation of the yield locus is adopted and the yield condition for the whole structure can be expressed in terms of hyperplanes [17] as follows: d0=N tQ-R~0
(2)
where dp is the vector of plastic potentials, N is the matrix of unit outward normals to the yield hyperplanes, R is the plastic capacity vector and 0 is the null vector. We assume that the material is elastic-perfectly plastic and that the plastic strain vector p can be represented as p = Nk (3) where k is a vector of non-negative plastic multipliers. Loading is assumed to be piecewise proportional. Interpreting the stresses as the sum of the linear-elastic response Q¢ to the loads and the self equilibrated stresses caused by the plastic strains [17], we have for each load multiplier t: Q = tQ ~ + Zp
(4)
where Z is a matrix of influence coefficients. The following relations then describe the nonholonomic (irreversible) elastoplastic problem: - ~ = R - tNtQ e - NtZNk
(5a)
- ~/> 0
(Sb)
~,/> 0
(5c)
dpt ~. ~> 0
(5d)
with the initial condition k = 0, t = 0; a superimposed dot represents rate. Problems described by these relations have been considered by m a n y researchers, with numerous references given in the reviews of Maier and Munro [18] and Maier and Smith [19]. The corresponding holonomic relations (non-linear elastic, reversible) are: - ~b = R - tNtQ ~ - N t Z N k
(6a)
- dp ~> 0
(6b)
k t> 0
(6c)
d~tk = 0
(6d)
Relations (6) differ from those in (5) in that Eqns (6c) and (6d) allow the value of 2 to decrease, an occurrence which cannot take place in the nonholonomic system. Also, at any stress point in the holonomic system, plastic deformation must become zero before elastic MS
33: 5 - E
382
R.R. WAKEFIELDand F. T1N-LOl
unloading can occur. We make use of both the nonholonomic and holonomic conditions in our solution procedure. The above equations fall into a special class of mathematical programming problems known as linear complementarity problems (LCP), which have received considerable attention in operations research; for example, see the paper by Cottle and Dantzig [20]. The system of equations represented by (6) can be solved parametrically for a given value of t to generate a solution ~.(t). Corresponding values of stresses and strains for the load multiplier t can then be determined using the following relationships: Q = tQ e + ZNk(t)
(Ta)
and q=S-IQ+p = t S - 1 Q e + ( S - 1 Z N + N)~,(t)
(7b)
where S is the unassembled structure stiffness matrix and an inverse is represented by the superscript - 1. The complete history of stresses and strains can thus be obtained. There are numerous methods for solving the LCP, e.g. Refs [13, 14, 21, 22, 23]. In this work, we use an elegant and direct method of solution proposed by Kaneko [14]. Before describing the solution procedure, some comments need to be made on the matrices and vectors involved. Namely, R ~> 0 and the matrix M = - N tZN is positive semi-definite. In the description of the procedure for solution we will adopt Kaneko's notation [ 14] and refer to a general iteration tableau:
0) where {3} is the subset of active rows of the system (i.e. those that have reached the yield surface and deformed plastically), {~} is the subset of rows that have not previously been active, and {/~} is the subset of inactive rows that have been active previously (i.e. those that have deformed plastically and then unloaded). The union {c~ w /~} then makes up the subset of inactive rows. Quantities k, NtQ e and 1VI represent the current values of the variables. The analysis proceeds as follows:
Step 1 Set up the initial tableau. Initialize the variables.
Step 2 If all - N tQe /> 0 then stop.
Step 3 Determine the minimum value of the load multiplier to cause any row of the tableau to be equal to zero, that is, for any element in the structure to reach yield. The minimum value is the next load step t* given by:
t* = m i n
i~a N]Q~
.. N - 4j Q je
> 0 .
(9)
Step 4 Work out {a}, the subset of rows that reach the yield surface as a result of the load step t*; {a} is the subset of {~ w/3} such that j ~ {a} if: Rj
+ t*NtQ *
+ l~/lji~,~ -~- O.
(10)
Large displacementelastoplastic analysis of frames
383
Form the following LCP sub-problem and solve: ~,~ = - NtQ,~ + 1Vl~d~ + 1('1~,~
(lla)
~b,, =
(11 b)
- N t Q,~ + l£'l,,,s dp,s + lgl,,,, ~,,,
~b >/0, k ~> 0, ~btk = 0.
(llc)
Use Graves' algorithm [24] to solve the LCP. It is able to solve the LCP if the M matrix is positive semi-definite, and identifies a solution ray if a solution exists [25]. Nonuniqueness of solution can also be detected here [26]. If the LCP has no solution then the structure has reached an unstable condition such as collapse or a limit point. This is indicated by the fact that the matrix M is no longer positive semi-definite. The analysis is then terminated. If a solution is found then use the pivots that gave a solution of the LCP sub-problem to perform pivots on the entire tableau. Then return to step 2. A detailed description of this algorithm can be found in Kaneko's original paper [14]. As can be seen from the above, similar with other schemes, the procedure automatically calculates t for which the next plastic event (activation of a yield plane or unloading) occurs. Of course the analysis can be carried out for a specified maximum value of t; this is useful in considering piecewise linear proportional loading stages. It must be mentioned that two computational problems arise in the implementation of Kaneko's method to solve large scale problems. In particular, the large size of the matrices involved and the resulting storage problems, and the propagation of roundoff errors and their effects on the accuracy of solutions generated. These problems have been addressed by the authors and details of their solution have been discussed in another paper [1].
Piecewise proportional loadiny and its application to geometrically non-linear analysis As mentioned in the previous section, the mathematical programming formulation can be used with piecewise proportional load paths. This is in fact utilized in the present method to keep account of the previous load history of the structure when the geometry is updated. Such a solution procedure, based on the theory of "stepwise holonomic" elastoplasticity with piecewise linearized yield conditions, was developed by Maier and his co-workers, eg. Refs. [12, 15, 27]. If we consider the structure at load level F i, the vector R in Eqn (6a) contains all the information about the previous structural behaviour, including whether any points on the structure have reached yield and how much further elastic deformation can take place before yielding occurs [8]. If we then apply a load increment AF to the structure, and perform the analysis, the pivoting within the P L C P formulation will occur at rows in the tableau corresponding to previously yielded nodes (i.e. those with an R value equal to zero). This has the effect of modifying the elastic problem that was set up at the start of the increment to include the effect of previous plastic and elastic deformations. Therefore, if we consider each new load increment as a new proportional load stage and use the value o f ~ at the end of the previous load stage for R in the next load stage, previous loading can be taken into account. This can be used to advantage in the geometrically non-linear analysis since the vector R can retain the state of the structure when geometric updating is undertaken.
Computational technique The non-linear response of the structure is generated by using an incremental load approach with a Newton-Raphson type iteration performed at each load level to satisfy joint equilibrium equations. If we consider a load level F i and assume that the corresponding deformed configuration of the system may be described symbolically as u i, then it is necessary for us to determine the configuration I! i + 1 corresponding to load level F i+l F ~ + AF in which AF represents a small increment in load. If we assume that the small load increment represents a new load step applied to the existing structure, we can proceed as follows: =
384
R.R. WAKEFIELDand F. TIN-LOl
Step 1 Calculate F i+1 = F i + AF, the new load level.
Step 2 Calculate the stiffness matrix K T of the structure as the sum of the linear stiffness matrix K L of the structure in its present configuration and the second order geometric stiffness matrix K~ [-28].
Step 3 Use the stiffness matrix and AF to solve for Au e the incremental elastic displacements and subsequently for AQ e the incremental elastic stresses.
Step 4 Use the P L C P method to solve for the increments in the elastic stresses and the plastic multipliers. Update the plastic multipliers as required and calculate the total displacements. Calculate the value of d~.
Step 5 Determine the new structure geometry, and calculate the new member lengths and direction cosines.
Step 6 Calculate the member end forces in local coordinates and, using the updated geometric information, work out the member end forces in global coordinates.
Step 7 Calculate the unbalanced joint forces by subtracting the member end forces from the applied nodal forces.
Step 8 Compare the unbalanced joint forces with a prescribed convergence criterion, for example, a force norm [29].
Step 9 Set the new value of the vector R equal to the value of - 4.
Step 10 If convergence has been achieved, return to Step 1 and apply the next load increment. If convergence has not been achieved, that is, the unbalanced force vector is larger than acceptable, use the unbalanced force vector as the new load increment and continue.
Step 11 Calculate the stiffness matrix of the structure. Solve for the elastic displacements and stresses and set up the PLCP.
Step 12 Perform a holonomic analysis. If all load multipliers are less than zero, reverse the out of balance forces (the force vector), solve the LCP, and reverse the results when updating 2 values, before doing subsequent stress calculations.
Step 13 Determine the new structure geometry and calculate updated member lengths and direction cosines.
Step 14 Return to Step 6. Before illustrating the operation of the algorithm on an example the following remark is appropriate. In the present scheme which combines a classical N e w t o n - R a p h s o n approach with a mathematical programming solution technique, the matrix M needs to be updated at every iteration. This appears to result in a computational burden of similar magnitude to that reported for the nonlinear programming approach given in [11]. In the latter direct energy minimization scheme, an iterative method is included in the nonlinear program
Large displacement elastoplastic analysis of frames
385
solution procedure. Direct comparison of the two approaches without detailed numerical testing on similar structures is difficult. ILLUSTRATIVE EXAMPLE
To illustrate the performance of the method, the analysis of a simply supported beam under uniform loading will be discussed. The beam is illustrated in Fig. l(a) and has previously been analysed by Backlund I-2] who used a finite element approach with hybrid Load Intensity w (MNIm)
1
~ pL
~ pL
i pL
~ pL
~ pL ~1
l =
0.5m
=
I-
-IIJ2
-IL
-I-
-I-
L
(a)
L
-IL
-I L/2
(b)
1.0
1.0
M/Mp
• , .
M/Mp
N/Np
1.0
N/Np
(c)
1.0
(d)
2
/,,"" °°" I
•,°
/
Z
7
I "~
~
"•~•
/.&
.o~
S ~
-
ElasticDeflection
.... ~---
LB Yield Locus
- - ~- -
UB Yield Locus
x
Backlund1976
0
0.00
0.01 Central Deflection (rn)
0.02
(e)
FIG. 1. Simply supported beam example: (a) beam layout; (b) half beam model; (c) lower bound yield locus; (d) upper bound yield locus; (e) load-deflection curves.
386
R.R. WAKEFIELDand F. TIN-LOl
beam elements. Contro et al. [7-] analysed the same structure both as being linear elastic and elastic-perfectly plastic. In the latter case, a non-linear programming approach was used. The beam is of uniform rectangular cross-section of 1 m width, 0.1 m depth, and made from a material with a Young's modulus of 220 x l09 N / m 2 and a yield stress of 300 x l06 N / m 2. The beam is subjected to a uniform load w. The non-linear yield locus, in terms of bending moment (M) and axial force (N), at any node is expressed by the equation +
- 1= 0
(12)
where Mp is the plastic capacity in pure bending and N o is the plastic capacity under pure axial load. Two piecewise linear approximations of this yield surface are used in our example to demonstrate the effect of linearization on the load deflection behaviour of the structure. The
W/8
15.984
,w/8
J w/8
1 w/8 .
I
¥~
T
I W/8
in~ L/16
I____1 I
I
L/8
L/8
I
L/8
I_
1_/16
L/2
__J
__1___1_
2-'--'--'--'--'-
2
I
L = 47.992 in
(a) 0.6" x
0.5"
0.4"
0.3"
~ 02"
01 /
x
Experimental (Maj=d, 1972)
00 0
1
Vertical Deflection of Centre Node (inches)
(b) FIG. 2. Pitched roof portal frame example: (a) dimensions and loading; (b) load-deflection curves.
Large displacement elastoplastic analysis of frames
387
approximations are shown in Figs l(c) and (d). Figure l(c) represents a lower bound approximation and Fig. l(d) represents an upper bound approximation. The structure is modelled using six beam elements as shown in Fig. l(b). The perfectly elastic structure was first analysed and the results, as shown in Fig. l(e), were found to be in good agreement with those of Backlund [2]. The results of the elastoplastic analyses with lower and upper bound polygons clearly show the influence of the yield surface discretization. The load-deflection response for both yield locus approximations give somewhat different results [Fig. l(e)]. The incidence of first yield is different in both cases. This is because the distance of the yield hyperplane from the origin, in the direction of the stress path taken during loading, is different. It is clear from Fig. l(d) that when bending moment is significant, the distance of the yield line from the origin is further in the upper bound case so that a greater load is required to cause first yield. As the load level is increased, the 4W W
4W
4W
==
r_~
L 4W
-I
L (a) 0.6
Ig 05"
04"
03"
_j
02Analysis
/ 0.1
0.0 000
•
Home,1963
x
De F r e i t a s
/ and Smith
i
i
i
0.02
0.04
0.06
1985
0.08
u/L
(b) FIG. 3.
Two storey frame example: (a) dimensions and loading; (b) load-deflection curves.
388
R . R . WAKEFIELD a n d F. TIN-Lol
linearization has a further effect through the normality rule of plasticity. The direction of the outward normal to the yield curve influences the plastic deformation taking place [7]. This indicates that in order to accurately model the deformation at any stress point, which is at yield, it is necessary to have the outward normal to the linearized yield locus at that point in a direction similar to the outward normal to the non-linear yield locus at that stress point. This is apparent in the cases illustrated [Fig. l(c)]. The upper bound yield locus models the deformation more accurately as the axial force in the beam becomes larger. This is because the upper bound yield line is tangent to the true yield locus at N = N p = 1 (i.e. the slopes of the normals to both curves are similar). When the lower bound yield locus is used, the deformation is more accurately modelled at lower load levels when both moment and axial force are significant. The problem of accurately representing the outward normal to the yield locus can be overcome by refining the piecewise linearization with particular consideration given to regions of the yield locus where there are large changes in the direction of the outward normal. FURTHER EXAMPLES The method presented has been used to analyse several frames that have appeared in the literature. These examples give an indication of the capability and accuracy of the method. The single bay pitched roof portal frame of Fig. 2(a) was experimentally tested by Majid [30] under the loading shown. The frame was modelled as 12 elements with the following
i2 D
L = 19.5 in
8 elements of equal length
(a)
(b)
+0
j
400
j
/
+.
'
-o m
A~
~
200
Elastic
o"
,c ,/.....
o
o 4p
.... t---
LB Yield Locus
- -,A- -
UB YieldLocus
o 0 0.0
ExperimenlalValues
i
i
i
i
0.2
0.4
0.6
08
1.0
Central Deflection (inches)
(c) FXG. 4. Fully fixed beam example: (a) dimensions and loading; (b) ball" beam model; (c) load deflection curves.
Large displacementelastoplasticanalysisof flames
389
properties: area 0.262in 2, moment of inertia 0.0057in 4, modulus of elasticity 28076 kips/in 2, yield stress 40.61 kips/in 2, plastic moment 1.36 kips-in and axial plastic capacity 10.64 kips. The flame has also been analysed by other researchers [3, 4]. In the present analyses, a lower bound yield locus [Fig. l(c)] was used. The obtained load-deflection diagram for the vertical deflection of the centre node is given in Fig. 2(b), along with experimental results reported by Majid [30]. The theoretical results are very similar to those reported in the literature [3, 4]. The second example is a two storey single bay flame first analysed by Home [31] and later analysed by Jennings and Majid [32], Corradi [6], De Freitas and Smith [10] and Wong and Tin-Loi [3]. The frame is shown in Fig. 3(a); member details were as follows: modulus of elasticity 30 x 1 0 6 lb/in 2, moment of inertia 51.2 in 4, area 3.2 in 2, length (L) 400 in, plastic moment 464 kip-in and axial plastic capacity 115 kips. The yield surface used in our analysis was the lower bound shown in Fig. l(c). Results of the analysis are presented in the form of a load-deflection curve, along with the results reported by Horne [31] and De Freitas and Smith [10]. Our results are in good agreement with reported ones. The final example is a fully fixed beam and subjected to lateral loading, as shown in Fig. 4(a). The beam, tested by Campbell and Charlton [33], was modelled using eight elements of equal lengths to represent half the beam [Fig. 4(b)]. The elements had the following properties: breadth 0.310 in, depth 0.304 in, span 19.5 in, axial stiffness 2.695 x 105 lb, flexural stiffness 20760 lb in 2, y!eld stress 16.1 tons/in 2, plastic moment 258 lb in and axial plastic capacity 3395 lb. Both the upper bound and lower bound yield surfaces shown in Figs. l(c) and (d) were used in the analysis. The load-deflection curves are given in Fig. 4(c) along with the experimental results of Campbell and Charlton [33]. The curves demonstrate similar behaviour to that found in the illustrative example. CONCLUDING REMARKS An incremental mathematical programming method for large displacement analysis, suitable for large structures, has been presented. The results of analyses of several well known examples from the literature demonstrate the accuracy of the method. The following remarks should be made: (a) In some cases the piecewise linearization of the yield locus can have a significant effect on the results obtained as was demonstrated in two examples described in this paper. An attempt should therefore be made to make the slope of the outward normal of the linearized yield surface close to the slope of the outward normal to the non-linear yield locus, in regions of the yield locus which are likely to be activated. (b) In the case of structures that exhibit a falling branch response after reaching a limit point, the method is only capable of tracing the load-deflection path up to that limit point. (c) While the solution of large scale problems has not been demonstrated explicitly in the examples given here, the method presented essentially only involves repeated analyses of the type detailed in Ref. [1]. One of the examples reported therein was that of the small displacement elastoplastic analysis under combined stresses of a 12 storey two-dimensional frame with 252 active degrees of freedom. The techniques for handling the large matrices involved and for controlling the propagation of round-off errors will still be effective for the present combined geometrically and materially nonlinear analysis. (d) The extension of the analysis to include higher dimensional yield surfaces, that is greater numbers of generalized stresses, is straightforward in terms of Eqn (1). (e) The number of incremental load steps needed in the analysis could be reduced by using a more accurate stiffness formulation, e.g. Refs [34, 35] than that used in the present method. Further work is at present being undertaken by the authors to widen the capabilities of the method to perform analysis of three-dimensional spatial structures. The basic formulation remains unaltered. However, the non-vectorial behaviour of large rotations in threedimensional space will need to be taken into account along with the effects of six generalized
390
R.R. WAKEFIELD and F. TIN-LOl
stress components. The more accurate stiffness matrices referred to in (e) above will also be implemented. Acknowledyements The authors would like to thank the reviewers for their constructive comments on an earlier version of this paper. REFERENCES I. R. R. WAKEFIELD and F. TIN-Lol, Large scale nonholonomic elastoplastic analysis using a linear complementarity formulation. Comput. Meth. Appl. Mech. Engn9 84, 229 (1990). 2. J. BACKLUND,Large deflection analysis of elasto-plastic beams and frames. Int. J. Mech. Sci. 18, 269 (1976). 3. M. B. WONG and F. TJN-LOl, Analysis of frames involving geometrical and material nonlinearitics. Comput. Struct. 34, 641 (1990). 4. A. KASSlMAL|, Large deformation analysis of elastic-plastic frames. J. Struct. Engng, ASCE 109, 1869 (1983). 5. J. L. MEEK and S. LOGANATHAN, Geometrically non-linear behaviour of space frame structures. Comput. Struct. 31, 35 (1989). 6. L. CORRADI, Large displacement elastic-plastic analysis, in Engineeriny Plasticity by Mathematical Pro.qrammin9 (edited by M. Z. Cohn and G. Maier). Pergamon Press, Oxford (1979). 7. R. CONTRO, L. CORRADI and R. NOVA, Large displacement analysis of elastoplastic structures a nonlinear programming approach. Solid Mech. Arch. 2, 433 (1977). 8. A. FRANCHIand M. Z. Cony, Strain-softening and large displacement analysis in structural plasticity. Comput. Struct. 11, 421 (1989). 9. G. MAIER and R. CONTRO, Energy approach to inelastic cable structure analysis. J. En.qn9 Mech. Div., ASCE 101,531 (19751. 10. J. A. T. DE FREITAS and D. L. SMITH, Elastoplastic analysis of planar structures for large displacements. J. Struct. Mech. 12, 419 (1984-85). 1I. R. CONTRO and G. MMER, Inelastic analysis of pretensioned cable structures by mathematical programming, in Steel Structures: Recent Research Advances and their Applications to Desiqn (edited by M. N. Pavlovic), pp. 85--105. Elsevier Applied Science, London (1986). 12. O. DE DONATO and G. MAIER, Historical deformation analysis of elastoplastic structures as a parametric linear complementarity problem. Meccanica 3, 166 (1976). 13. G. MAIER, S. GIACOMIN~ and F. PATERLIM, Combined elastoplastic and limit analysis via restricted basis linear programming. Comput. Meth. Appl. Mech. Engn.q 19, 21 (1979). 14. I. KANEKO. Complete solutions for a class of elastic-plastic structures. Comput. Meth. Appl. Mech. Enqno 21, 193 (1980). 15. O. DE DONATO and G. MAIER, Mathematical programming methods for the inelastic analysis of reinforced concrete frames allowing for limited rotation capacity. Int. J. Num. Meth. Engn9 4, 307 (1972). 16. L. CORRADI and G. MMER, Extremum theorems for large displacement analysis of discrete elastoplastic structures with piecewise linear yield surfaces. J. Optim. Theory Applic. 15, 50 (1975). 17. G. MINER, A matrix structural theory of piecewise-linear elastoplasticity with interacting yield planes. Meccanica 5, 54 (1970). 18. G. MA1ER and J. MUNRO, Mathematical programming applications to engineering plastic analysis. Appl. Mech. Rev. 35, 1631 (1982). 19. G. MAIER and D. L. SMITH, Update to "Mathematical programming applications to engineering plastic analysis" [AMR 35, 1631 (1982)]. Applied Mechanics Update, ASME (1986). 20. R. W. COTTLE and G. B. DANTZIG, Complementary pivot theory and mathematical programming. Linear Al,qebra Appl. 1, 103 (1968). 21. D. L. SMITH, The Wolfe Markowitz algorithm for non-holonomic elastoplastic analysis. Enonq Struct. 1, 9 (1978). 22. F. CAS('IATI, Elastic-plastic deformation analysis: a parametric linear programming method. J. Mec. Appl. 2, 269 (1978). 23. O. DE DONATO and G. MAIER, Local unloading in piecewise-linear plasticity. J. Engng. Mech. Div., ASCE 102, 383 (1976). 24. R. L. GRAVES,A principal pivoting simplex algorithm for linear and quadratic programming. Operations Res. 15. 482 (1967). 25. G. MAIER, Incremental plastic analysis in the presence of large displacements and physical instabilizing effects. Int. J. Solids Structs 7, 345 (1971). 26. R. R. WAKEFIELD and F. TIN-LOI, Mathematical programming and uniqueness in nonholonomic plasticity. Comput. Struct. 34, 477 (1990). 27. O. DF. DONATO and G. MAIER, Finite element elastoplastic analysis by quadratic programming: the multistage method, 2rid Int. Conf. on Struct. Mech. in Reactor Technology, 10-14 September, Paper M2/8, 1--12 (1973). 28. H. C. MARTIN, Finite elements and the analysis of geometrically non-linear problems. Proceedings, Japan-US Seminar on Matrix Methods of Structural Analysis and Design, 25-30 August, Tokyo, Japan (1969). 29. P. G. BERGAN and R. W. CLOUGH, Convergence criteria for iterative processes. AIAA J. 10, 1107 (1972). 30. K. I. MAJID, Non-Linear Structures. John Wiley, New York (1972). 31. M. R. HORNE, Elastic-plastic failure loads of plane frames. Proc. Roy. Sot'. London 274, 343 (1963).
Large displacement elastoplastic analysis of frames
391
32. A. JENNINGS and K. MAJID, An elastic-plastic analysis by computer for framed structures loaded up to collapse. The Structural Engineer 43, 407 (1965). 33. T. I. CAMPBELL and T. M. CHARLTON,Finite deformation of a fully fixed beam comprised of a non-linear material. Int. J. Mech. Sci. 15, 415 (1973). 34. C. ORAN,Tangent stiffness in space frames. J. Struct. Div., ASCE 99, 987 (1973). 35. J. L. MEEK and H. S. TAN, Geometrically non-linear analysis of space frames by an incremental iterative technique. Comput. Meth. Appl. Mech. Engng 47, 261 (1984).