A Radial Return Algorithm Application in Elastoplastic Frame Analysis Using Plastic Hinge Approach

Abstract

A new method is presented for first-order elastoplastic analysis of framed structures using a radial return predictor/corrector solution strategy. The proposed method assumes plastic hinge formation coupled with a yield surface. The yield surface is defined as a general function of axial force, shear forces, twisting, and biaxial bending moments on the cross-section of the frame. The material is regarded as linear and elastic-perfect plastic. The plastic deformations are governed by the normality criterion. Combining the Newton-Raphson method and the radial return algorithm, a consistent tangent modular matrix is proposed and fast and converging algorithms are presented. Examples demonstrate the accuracy and effectiveness of the proposed method.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Recommended by Giuseppe Rega

Department of Civil Engineering, University of Brasilia (UnB), 70910-900 Brasilia, DF, Brazil

Received 13 April 2010; Revised 15 July 2010; Accepted 16 August 2010

1. Introduction

Over the last few decades, in the context of computational plasticity, efficient algorithms have been developed for the integration of constitutive models for fragile and ductile materials. Excellent references on such integration schemes are the books of Simo and Hughes [1], Crisfield [2, 3], Doltsinis [4], among others. Such references describe in detail the implicit algorithms formulated in a continuum-based approach, and investigate the numerical performances of such algorithms in terms of efficiency, precision, robustness and convergence rate. Also in the literature there are research papers (see e.g., [5-8]) using implicit algorithms formulated in the stress resultant space for the collapse analysis of elastoplastic frames.

This paper presents a new method for a first-order elastoplastic analysis (i.e., small strains and small displacements) of framed structures under loading-unloading cycle based on the concepts of (a) radial return predictor/corrector algorithm and (b) limit analysis or the so-called "plastic hinge" approach. It is noted that the plastic hinge approach is an active area of research and able to deal with localization analysis [9], strong discontinuity, stress softening at failure, localized dissipative mechanism [10, 11], and collapse load of reinforced concrete frames [12]. For the plastic hinge analysis, this research presumes a generalized yield surface in the six-dimension space of stress resultants or generalized forces. The yield surface used in the proposed formulation is assumed as a continuous and convex function of the axial force, shear forces, twisting and biaxial bending moments (six generalized forces) acting on the structure cross-section. The proposed yield surface is a general expression that takes care of the different interactions of the generalized forces on the cross-sections. The material is considered to behave linearly elastic perfectly plastic and presents no strain hardening. The plastic deformations are governed by the normality principle and are confined to zero-length plastic zones at the element ends. The end sections can undergo an abrupt transition from a fully elastic to a fully plastic state. Combined stress resultants that initiate yielding on the cross-section are assumed to produce full plastification of the whole section. This paper describes in detail the proposed formulation presenting a new development for the Single-Vector Return Algorithm (hereafter named 1VRA) and an original Two-Vectors Return Algorithm (hereafter named 2VRA). The 1VRA is used for simulating one plastic hinge at one end of the beam element while the 2VRA simulates the appearance of two simultaneous plastic hinges at the ends of the beam element. The paper also brings the deduction of the consistent tangent modular matrix which together with the developed algorithms is fundamental for obtaining accuracy and convergence [1]. At the end of the paper, three examples are presented and discussed demonstrating the accuracy and effectiveness of the proposed method.

2. The Yield Surface Concept

For practical problems there are approximate plastic interaction surfaces for each shape of beam cross-section, for instance, for rectangular and I-shaped sections of steel beams under bending and axial interaction; see [13, 14]. The derivation of the yield surfaces in terms of the generalized cross-sectional forces constitutes a difficult task. There are many ways of transforming the known yield conditions in terms of stress components to the space of the generalized forces [15]. All of them, however, are approximate in nature. There exists a vast literature on the subject (e.g., [16, 17]). In this paper, the yield surface is assumed to be a continuous and convex function of the generalized forces on a cross-section and may be mathematically represented as [figure omitted; refer to PDF] where |(·)| is the absolute value of (·) , _Fx is an axial force, and _Fy and _Fz are shear forces. _Mx is the twisting moment, and _My and _Mz are bending moments. _Fxp is the plastic axial force, _Fyp and _Fzp are plastic shear forces, respectively. _Mxp is the plastic twisting moment, and _Myp and _Mzp are plastic bending moments. The constants _αi represent positive real numbers which are functions of the geometric shape of the cross-section. As noted earlier, for practical purposes, the function Φ in (2.1) must be determined taking into account the cross-section shape. Several possibilities of Φ are considered later on in this paper, but of special interest, at this moment, is the basic observation that: (1 ) cross-sections for which the stress resultants rest inside the yield surface are considered elastic; (2 ) cross-sections with the stress resultants on the locus of the yield surface are considered fully plastic, and finally (3 ) the stress resultants are not allowed to be outside the yield surface as the material is of elastic-perfect behavior. To bring the outside stress resultants back to the yield surface a radial return predictor/corrector solution strategy is applied. Consequently, it is necessary to calculate the first and second derivatives of the yield function Φ . The derivatives are carried out with respect to the generalized forces on the cross-section. The following theory is presented for independent hinges located at ends of the beam. In case the yield surface exhibits corners (two yield surfaces are active), the 2VRA is used to bring the outside stress resultants back to the yield surface--for more details see chapter 14, in [3].

2.1. The First Derivative of the Yield Surface

To obtain the plastic potential function at each end of the beam element, the first derivative of Φ in (2.1) with respect to the generalized force vector, expressed here in indicial notation as _Fj , can be written as follows: [figure omitted; refer to PDF] where sgn (·) denotes the signal of the components of the vector _Fj . In order to obtain a compact format, these derivatives may be collected in a vector. This vector defines the plastic potential flow at the ends of the beam element. This vector may be written, respectively, for beam element ends 1 and 2 as [figure omitted; refer to PDF] It is noted that the sub indexes i, j, k, l, ..., and so forth are integers and they vary from 1 to 12; as 12 is the number of degrees of freedom for a 3D beam element. Also observe that in the last equations the left hand side is in indicial notation while the right-hand side is in expanded notation for a better understanding of the vector components.

2.2. The Second Derivative of the Yield Surface

The gradient of the plastic potential vector is obtained by the differentiation of each component of the vectors in (2.3) with respect to the generalized nodal forces. To illustrate, the second derivatives of the first component ∂Φ/∂_Fx with respect to the generalized nodal forces are expressed as [figure omitted; refer to PDF] In an analogous procedure, the second derivatives of the other components of the vector {∂Φ/∂_Fj } may be obtained. Collecting the second derivatives in a matrix, one can get the matrices that contain the gradient of the plastic potential flow at the ends of the beam element. Such matrices, for ends 1 and 2, are expressed as [figure omitted; refer to PDF] [figure omitted; refer to PDF]

3. A Backward Euler Algorithm

From now on we make use of the indicial notation for a better understanding of the algebraic operations. Suppose that there is a generalized force vector at node 1 outside the yield surface. The backward Euler algorithm (see chapter 6 of [2]) is based on the following equation: [figure omitted; refer to PDF] [figure omitted; refer to PDF] where _F...i is the nodal force vector at the last load step that has achieved convergence.The term _Kij d_Uj is the elastic predictor vector. _Kij is the stiffness matrix of the 3D beam element. d_Uj is an increment in the vector _Uj which is defined as the nodal displacement vector. ^Fitrial is the elastic trial nodal-force; _{{∂Φ/∂Fj }1} is the plastic potential defined in the trial nodal-force; _λ1 is the plastic multiplier and the term _{λ1Kij{∂Φ/∂Fj }1} is the plastic corrector. _F...i is the nodal-force after correction. In (3.1b), L is the element length. A is the cross-section area of the element. E is the Young modulus. G is the shear modulus. _Iy and _Iz are moments of inertia of the cross-section with respect to y axis and z -axis, respectively. _Jx is the polar moment of inertia of the cross-section. The element formulation is based on the Euler-Bernoulli theory. In (3.1b), the shear deformation is neglected. The element interpolation functions are linear for axial displacement and twisting. Hermite interpolation is used for bending. Figure 1 shows the geometric interpretation of the vectors in (3.1a).

Single-vector return algorithm--one plastic hinge.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Generally, the starting guess ^Fitrial does not satisfy the yield condition. Consequently, an iterative scheme is necessary to return the force vector outside the yield surface to the allowable stress resultants located on the yield surface.

3.1. Single-Vector Return Algorithm

When one plastic hinge takes place at one end of the beam element, the Single-Vector Return Algorithm (1VRA) is applied. In the following equations it is assumes that the sub index "1" is related to node "1" of the beam element. The 1VRA brings the nodal-force vector back to the yield surface. To apply the 1VRA scheme, the residual force vector _ri is defined as [figure omitted; refer to PDF] This residual force vector represents the difference between the current force state _Fi and Backward Euler force_F...i . The trial force ^Fitrial in (3.2) is kept constant during the iteration process. A first-order Taylor's series expansion can be applied to (3.2) to obtain an expression for the new residual force vector ^rinew in terms of the old residual force vector ^riold , therefore [figure omitted; refer to PDF] where d_Fi is an infinitesimal increment of the generalized force vector _Fi ; d_λ1 is the variation of the plastic multiplier _λ1 and [^∂2 Φ/∂_Fj ∂_Fk ]d_Fk is the change in the potential plastic vector _{{∂Φ/∂Fj }1} . The goal is to achieve ^rinew =0 in (3.3). For that reason [figure omitted; refer to PDF] To determine an expression to d_Fi (the correcting force vector) in (3.4), define the square matrix _Qik as [figure omitted; refer to PDF] Now, considering (3.4) and (3.5) and after some algebraic manipulations to solve (3.4) for d_Fi , it follows that [figure omitted; refer to PDF] A first-order Taylor's series expansion of the yield function Φ around the final nodal-force vector _Fi is applied. This series expansion is necessary to get a linear approximation to the new value of the yield function ^Φnew , therefore [figure omitted; refer to PDF] Since the expression for d_Fi is available in (3.6) and for ^Φ1new =0 , the variation of the plastic multiplier d_λ1 is readily found [figure omitted; refer to PDF] This iterative procedure is continued until the yield criterion Φ=0 is satisfied at the final force state; that is, ^rnorm =|_ri |/|^Fitrial |<TOL and ^Φnorm =|Φ|<TOL , where ^rnorm is a norm for the residual force vector. ^Φnorm is defined as the residual yield norm and TOL is a tolerance for convergence. In this paper, it is assumed that TOL =1^0-10 .

3.2. Two-Vectors Return Algorithm

When two plastic hinges occur at the ends of the beam element, the Two-Vector Return Algorithm (2VRA) is applied. This algorithm is used so that the nodal-force vector at both ends of the beam element can be brought back to the yield surface. The geometric interpretation of this algorithm is illustrated in Figure 2. The Backward Euler force obtained with the 2VRA may be defined by the following expression: [figure omitted; refer to PDF] To derive a scheme for 2VRA, the residual force vector _ri is defined as [figure omitted; refer to PDF] This residual force vector represents the difference between the current force state _Fi and the Backward Euler force _F...i . Analogous to the 1VRA, the trial force ^Fitrial is kept constant during the iteration process. A first-order Taylor's series expansion can be applied to (3.10) to obtain an expression for the new residual force vector ^rinew in terms of the old residual force vector ^riold , therefore [figure omitted; refer to PDF] Making ^rinew =0 in (3.11), it follows that [figure omitted; refer to PDF] To determine an expression to d_Fi (the correcting force vector) in (3.12), define the square matrix _Qik as [figure omitted; refer to PDF] Now considering (3.13) and (3.12) and after some algebraic manipulations to solve (3.12) for d_Fi , it follows that [figure omitted; refer to PDF] First-order Taylor's series expansions of the yield function at node 1, _Φ1 , and at node 2, _Φ2 , around the final force vector _Fi are now developed. These series expansions generate linear approximations of the new values of the yield functions ^Φ1new and ^Φ2new as [figure omitted; refer to PDF] Since the expression for d_Fi is available in (3.14) and the goal is to achieve both ^Φ1new =0 and ^Φ2new =0 , then [figure omitted; refer to PDF] The two previous expressions form a system of equations in d_λ1 and d_λ2 . By introducing auxiliary variables [figure omitted; refer to PDF] and now using Cramer's rule for the previous system of equations, the variations of the plastic multiplier d_λ1 and d_λ2 are readily found by the expressions [figure omitted; refer to PDF] where the auxiliary variables are defined as [figure omitted; refer to PDF]

Two-Vectors return algorithm--two plastic hinges.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

This iterative procedure is continued until the yield conditions _Φ1 =0 and _Φ2 =0 are satisfied at the final force vector state, that is, when ^rnorm =|_ri |/|^Fitrial |<TOL , ^Φ1norm =|_Φ1 |<TOL and ^Φ2norm =|_Φ2 |<TOL .

4. Consistent Tangent Stiffness Matrix

The final objective in deriving the backward Euler scheme for integration of elastoplastic constitutive equations is to use the previous algorithms (1VRA and 2VRA) in finite element computations. If the Newton-Raphson iterative scheme is used on a global equilibrium base, then the use of the so-called traditional elastoplastic stiffness matrix ^KijEP puts at risk the quadratic rate asymptotic convergence of the iterative process. In order to preserve the quadratic rate of convergence, a consistent stiffness matrix may be derived. In what follows, two derivations of consistent stiffness matrices are presented, respectively, for the 1VRA and 2VRA. Once the convergence criterion is satisfied and all the force points have returned to the yield surface, the consistent tangent stiffness matrices for any yielded element are updated before the start of the next loading cycle.

4.1. Single-Vector Return Algorithm

Starting from the Backward Euler equation (3.1a) and noting that vector _F...i is outside the yield surface and vector _Fi is on the yield surface, (3.1a) can be rewritten as [figure omitted; refer to PDF] The term d^Fitrial =_Kij d_Uj represents the elastic trial nodal-force increment and is associated with the nodal displacement increment d_Uj . The differential of (4.1) gives the following expression; [figure omitted; refer to PDF] After some algebraic treatment, [figure omitted; refer to PDF] Using matrix _Qik as defined by (3.5) and denoting the reduced stiffness matrix as _Rij =^{Qi[cursive l]-1}_{K[cursive l]j} , (4.3) changes into the following equation: [figure omitted; refer to PDF] It turns out that the form of (4.4) is similar to the nonconsistent form except for the change in _Kij to _Rij =^{Qi[cursive l]-1}_{K[cursive l]j} and for the fact that the normal to the yield surface is evaluated at the final force position. Assuming that the full consistency condition holds at the final force position, that is; Φ=0 , the differentiation of Φ with respect to the generalized force vector, considering (4.4), is [figure omitted; refer to PDF] and, therefore, the change in the plastic multiplier d_λ1 may be expressed as [figure omitted; refer to PDF] Finally, using (4.4) and (4.6), the elastoplastic consistent stiffness matrix ^KijAL may be determined by [figure omitted; refer to PDF]

4.2. Two-Vector Return Algorithm

Starting from the Backward Euler equation for the Two-Vector Return Algorithm, one can write the following expression analogous to (3.9): [figure omitted; refer to PDF] Differentiating (4.8) and performing algebraic manipulations, it follows that [figure omitted; refer to PDF] Using matrix _Qik as defined by (3.13) and denoting the reduced stiffness matrix as _Rij =^{Qi[cursive l]-1}_{K[cursive l]j} , we can rewrite the above equation as [figure omitted; refer to PDF] Again, it is assumed that the full consistency conditions should hold at the final force point, that is, _Φ1 =0 and _Φ2 =0 . Differentiating these two conditions and introducing the expression for d_Fi taken from (4.10), one can get [figure omitted; refer to PDF] Alternatively, by collecting terms algebraically, (4.11) are transformed in [figure omitted; refer to PDF] The previous two expressions represent a system of equations with d_λ1 and d_λ2 as unknowns. Introducing auxiliary variables _c1 , _c2 , _b11 , _b12 , _b21 and _b22 as [figure omitted; refer to PDF] With the help of Cramer's rule, the system in (4.12) is solved for d_λ1 and d_λ2 [figure omitted; refer to PDF] or alternatively [figure omitted; refer to PDF] If any plastic multiplier is negative, that is, d_λ1 <0 , or d_λ2 <0 , then this plastic multiplier is set to zero (i.e., d_λ1 =0 or d_λ2 =0 ) which corresponds to an initial plastic hinge. Again, introducing new auxiliary variables [figure omitted; refer to PDF] Equation (4.15) can be rewritten as [figure omitted; refer to PDF] Finally, using (4.10) and (4.17), the elastoplastic consistent stiffness matrix ^KijAL is obtained by the following expression: [figure omitted; refer to PDF]

5. Numerical Examples

Three examples are presented in this section to demonstrate the accuracy and the effectiveness of the methods proposed and described in the previous sections of this paper. An incremental iterative method based on the Newton-Raphson method combined with constant arc length control method is employed for the solution of the nonlinear equilibrium equations. In the following examples, a tolerance for convergence criterion is set TOL =1^0-10 . The examples test the global performance of the Newton solution strategy and the local performance of the backward Euler integration schemes. The examples are concerned with a right-angle beam, a two-bay, two-storey frame, and a four-legged jacket frame, all made of steel. To verify the quadratic ratio of convergence for the radial return algorithm presented in this paper, the Euclidean norm (see Sections 3.1 and 3.2) of the residual forces vector at the nodes with plastic hinges is analyzed. This analysis is performed examining the error norm under a load step arbitrarily chosen as is usually done by many authors [1-3]. Moreover to verify the quadratic ratio of the convergences on global equilibrium (using the elastoplastic stiffness consistent matrix) the Euclidean norms of the residual forces and of the energy are adopted. As before, the monitoring of these two norms was carried out considering load steps arbitrarily chosen. The elements are analyzed following their number sequence. Since "n " elements may share a common node, the plastic hinge, if necessary, is attributed to the proper node of the element consistent with the element analysis succession. To avoid singularity at the global stiffness matrix, at most (n -1) plastic hinge may be assigned to the common node.

5.1. Right-... Angle Beam

In this example, a right-angle beam under a concentrated load P is analyzed. The geometry, material properties (EI, bending stiffness and GJ torsional stiffness), load location, and yield surface equation are shown in Figures 3(a) and 3(b). For this loading condition the right-angle beam is subjected to both bending and twisting moments. The same problem was analyzed by Ueda and Yao [18] for small deformations. In the present study, each of the members of the right-angle beam is modeled by linear grid elements. The first-order inelastic analysis of this beam is carried out for loading, unloading, and reversal loading. The progressive development of the plastic hinge formation is indicated by numbers in Figures 3(d). The load-displacement response for loading-unloading cycle is shown in Figures 3(e) while in [18, 19] only loading is reported. For the collapse load, the result is an excellent agreement with the value reported by Ueda and Yao [18] and Hodge [19].

A right angle bent of circular section. (a) Geometry and loading. (b) Material properties and yield function. (c) Finite element mesh. (d) Plastic hinge formation sequence. (e) Load-displacement response.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

The values of the residual force and residual yield norms for a typical load step are summarized in Table 1 for the Single-Vector Return Algorithm and in Table 2 for the Two-Vectors return algorithm. Note that, for Table 1, load steps 47 and 445 were arbitrarily chosen while for Table 2, load steps 106 and 450 were selected. The quadratic rate asymptotic convergence of the backward Euler integration scheme is exhibited in these results. Note that the Euclidian norm of the residual force and residual yield norm are very similar during the iterative process.

Table 1: Residual force and residual yield norms for typical load steps during the formation of the second plastic hinge.

Single-vector return algorithm
Residual force norm rnorm	Residual yield norm Φnorm	Residual force norm rnorm	Residual yield norm Φnorm
Step 47	Step 47	Step 445	Step 445
1.47597039E-03	1.3055981E-03	8.56514641E-05	9.37779227E-05
1.48013748E-06	4.26512984E-07	5.9731394E-09	2.21160446E-09
3.72157338E-13	4.59632332E-14	9.21947265E-17	2.22044605E-16

Table 2: Residual force and residual yield norms for typical load steps during the formation of first and third plastic hinges.

Two-vector return algorithm
Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm	Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm
Step 106	Step 106	Step 106	Step 450	Step 450	Step 450
1.65303075E-03	1.20357837E-05	1.47408979E-05	3.87694523E-03	6.19251038E-04	4.56960163E-03
1.89552293E-06	1.23009056E-07	5.42558567E-07	1.60549380E-05	7.05089222E-07	5.25607314E-06
5.81448484E-13	1.04360964E-14	7.37188088E-14	4.57042787E-11	2.61368704E-12	7.06168457E-12

The values of the energy and residual norms for typical load step are summarized in Table 1 for the 1VRA strategy and in Table 2 for the 2VRA approach. The quadratic rate of the asymptotic convergence of the backward Euler integration scheme is exhibited by the results presented here (see Table 2).

The values of the energy and residual norms, for each iteration, are summarized in Table 3. The quadratic rate of the asymptotic convergence of the Newton iterative scheme is exhibited by the results here presented (see Table 3). In addition, it is observed that the Euclidian norm of the residuals lags behind the energy norm in the iterative process - see Table 3. This fact is explained because energy is a scalar product between residual forces and displacement increments. Those quantities decreases along the convergence process and accordingly their product (energy) decrease in a faster ratio.

Table 3: Error norms for Newton iterative scheme on global equilibrium.

Single-vector return algorithm
Residual norm	Energy norm	Residual norm	Energy norm
Step 16	Step 16	Step 250	Step 250
0.11634E-04	0.10000E+01	0.21960E-03	0.10000E+01
0.18054E-08	0.41898E-12	0.18174E-07	0.19603E-10
0.74031E-16	0.93424E-24	0.18630E-14	0.65129E-23

5.2. Two-Bay, Two-Storey Fame

A two-bay, two-storey rigid frame is now analyzed. Its geometric and material properties, the external loading conditions and the yield function are shown in Figure 4. The same problem was analyzed by Argyris et al. [17] for small deformation and large deformation using the computer program LARSTRAN; and by Halder and Ming [20] for large deformation. In the present study, each member is modeled using linear 2D beam elements. The first-order inelastic analysis is undertaken for loading and unloading conditions. The progressive development of the plastic hinges and the load-displacement response are shown in Figure 4(d). The results are in good agreement with the results presented by Argyris et al. [17].

A two-bay, two-storey rigid frame. (a) Geometry and loading. (b) Material properties and yield function. (c) Finite element mesh. (d) Plastic hinge formation sequence. (e) Load-displacement response.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

The values of the residual force and residual yield norms for load steps arbitrarily chosen are summarized in Table 4 for the Single-Vector Return Algorithm. Again the quadratic rate of the asymptotic convergence of the backward Euler integration scheme is exhibited by these results.

Table 4: Residual force and residual yield norms for typical load steps during the formation of the first plastic hinge.

Single-vector return algorithm
Residual force norm rnorm	Residual yield norm Φnorm	Residual force norm rnorm	Residual yield norm Φnorm
Step 11	Step 11	Step 60	Step 60
1.20935377E-04	2.41865226E-07	8.93804269E-04	6.25601990E-06
3.45076307E-11	6.60402844E-11	6.22404505E-09	9.11217879E-09
8.02626977E-17	0.0	5.12166402E-17	0.0

The values of the energy and residual norms, for each iteration, are summarized in Table 5. Once more, the quadratic rate of the asymptotic convergence of the Newton iterative scheme is exhibited by theses results. Note that the Euclidian norm of the residual lags behind the energy norm in the iterative process for the same reason as explained in the former example.

Table 5: Error norms for Newton iterative scheme on global equilibrium.

Single-vector return algorithm
Residual norm	Energy norm	Residual norm	Energy norm
Step 16	Step 16	Step 250	Step 250
0.80063E-02	0.10000E+01	0.23954E-01	0.10000E+01
0.18017E-03	0.24336E-06	0.14595E-01	0.82672E-03
0.15248E-08	0.50571E-13	0.15491E-15	0.11855E-08
0.49717E-13	0.39726E-23	0.87519E-13	0.28146E-19

5.3. Four-Legged Jacket Type

This example is concerned with a four-legged jacket. It is a type of platform structure often used for off-shore industries. The geometry and dimensions of the structure, loading condition, and material properties are specified in Figure 5. Two different structural systems denoted here as bracing type 1 and bracing type 2 are considered. The dimensions of the four-legged jackets, type 1 and type 2, are also present in Figure 5. In this work, each member of the four-legged jacket is modeled by linear 3D beam elements with first-order inelastic analysis performed considering loading and unloading conditions. This structure was dicretized with 8 nodes and 18 elements. Each structural member corresponds to one finite element. The development of plastic hinges and load-displacement response for both systems are also shown in Figure 5. The same example was studied by Shi and Atluri [21]. They performed a second-order analysis of this problem; therefore, their curves (load-displacement) and the collapse load represent lower bound limits. Our results are in good qualitative agreement with them. As expected, the curves in Figure 5 are slightly superior to the curves reported in [21].

A Four-legged jacket. (a) Geometry and loading. (b) Material properties, and yield function. (c) Plastic hinge formation sequence. (d) Load-displacement response.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

The values of the residual force and residual yield norms for arbitrary load steps (12 and 41) are summarized in Tables 6 and 7 using the 2VRA for the two different structural systems considered in this example. Again, the quadratic rate asymptotic convergence of the Backward Euler integration scheme is exhibited by the reported result. Note that the Euclidian norm of the residual force and residual yield norm are very similar during the iterative process. The values of the energy and residual norms along the iterations are summarized in Table 8 and Table 9 (for load steps 12 and 52) for the two different structural systems studied. Again for both systems the quadratic rate of the asymptotic convergence of the Newton iterative scheme is exhibited. In Tables 6 and 7, note the fast convergence of the radial return algorithms in just 4 iterations. For obtaining the structure equilibrium, see Tables 8 and 9, only 3 iterations were necessary even using a tight tolerance (TOL =1^0-10 ) parameter. The same observations can be verified in the examples studied before.

Table 6: Four-legged jacket type-1. Residual force and residual yield norms for typical load steps during the formation of the first and second plastic hinges.

Two-vector return algorithm
Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm	Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm
Step 12	Step 12	Step 12	Step 41	Step 41	Step 41
8.03978006E-03	9.81043834E-03	9.79441075E-03	1.36989462E-02	1.79089005E-02	1.78810582E-02
6.24153357E-05	2.38396940E-05	2.37928885E-05	1.96679153E-04	7.88306625E-05	7.86893201E-05
8.85460742E-10	1.48586699E-10	1.44924517E-10	9.31316656E-09	1.60199565E-09	1.58390967E-09
8.20520831E-17	0.0	3.33066907E-16	2.40162583E-17	0.0	0.0

Table 7: Four-legged jacket type-1. Residual force and residual yield norms for typical load steps during the formation of the first and second plastic hinges.

Two-vector return algorithm
Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm	Residual force norm rnorm	Residual yield norm Φ1norm	Residual yield norm Φ2norm
Step 12	Step 12	Step 12	Step 41	Step 41	Step 41
8.03359176E-03	9.80846569E-03	9.78495657E-03	1.36815406E-02	1.78952610E-02	1.78530190E-02
6.23470796E-05	2.38312272E-05	2.37622611E-05	1.96249346E-04	7.87644610E-05	7.85560337E-05
8.84047353E-10	1.42486467E-10	1.41957335E-10	9.29161239E-09	1.76712267E-09	1.64043112E-09
3.12120416E-18	1.11022302E-16	1.11022302E-16	8.12816025E-17	4.44089210E-16	0.0

Table 8: Four-legged jacket type-1. Error norms for Newton iterative scheme on global equilibrium.

Single-vector return algorithm
Residual norm	Energy norm	Residual norm	Energy norm
Step 12	Step 12	Step 52	Step 52
0.11997E-03	0.10000E+01	0.20903E-03	0.10000E+01
0.16535E-09	0.27673E-12	0.41175E-07	0.43179E-10
0.18777E-15	0.41131E-25	0.36921E-15	0.24472E-22

Table 9: Four-legged jacket type-2. Error norms for Newton iterative scheme on global equilibrium.

Single-vector return algorithm
Residual norm	Energy norm	Residual norm	Energy norm
Step 25	Step 25	Step 45	Step 45
0.34837E-03	0.10000E+01	0.75999E-03	0.10000E+01
0.10245E-05	0.44342E-08	0.81462E-06	0.89862E-08
0.95123E-11	0.12589E-15	0.13896E-11	0.88152E-16

6. Conclusions

In this paper, the radial return algorithms were tested with different yield surfaces given good results. It is noted that the algorithm formulated provides a way to keep the generalized force vector always on the yield surface or inside it. The proposed method deals simultaneously with two plastic hinges located at the ends of 3D beam finite elements under load-unloading cycles. Combining the Newton-Raphson method and the radial return method provides a "consistent" tangent modular matrix, and robust, accurate, and fast converging algorithms. The use of the "consistent" tangent modular matrix is fundamental for achieving the quadratic rate of the asymptotic convergence with the Newton's method. The Single- and the Two-Vector Return Algorithms have been proposed for simulating, respectively, one and two plastic hinges at the ends of the beam element.

References

[1] J. C. Simo, T. J. R. Hughes Computational Inelasticity , vol. 7, of Interdisciplinary Applied Mathematics, pp. xiv+392, Springer, New York, NY, USA, 1998.

[2] M. A. Crisfield Non-Linear Finite Element Analysis of Solids and Structures , vol. 1, of Essentials, John Wiley & Sons, Chichester, UK, 1991.

[3] M. A. Crisfield Non-Linear Finite Element Analysis of Solids and Structures , vol. 2, of Advanced Topics, John Wiley & Sons, Chichester, UK, 1997.

[4] I. Doltsinis Elements of Plasticity: Theory and Computation , vol. 1, of High Performance Structures and Materials, pp. x+311, WIT, Southampton, UK, 2000.

[5] A. S. Gendy, A. F. Saleeb, "Generalized yield surface representations in the elasto-plastic three-dimensional analysis of frames," Computers and Structures , vol. 49, no. 2, pp. 351-362, 1993.

[6] F. G. A. Al-Bermani, K. Zhu, "Nonlinear elastoplastic analysis of spatial structures under dynamic loading using kinematic hardening models," Engineering Structures , vol. 18, no. 8, pp. 568-576, 1996.

[7] B. Skallerud, B. Haugen, "Collapse of thin shell structures--stress resultant plasticity modeling within a co-rotated ANDES finite element formulation," International Journal for Numerical Methods in Biomedical Engineering , vol. 46, pp. 1961-1986, 1999.

[8] S. Krenk, C. Vissing-Jørgensen, L. Thesbjerg, "Efficient collapse analysis techniques for framed structures," Computers and Structures , vol. 72, no. 4, pp. 481-496, 1999.

[9] M. E. Marante, R. Picón, J. Flórez-López, "Analysis of localization in frame members with plastic hinges," International Journal of Solids and Structures , vol. 41, no. 14, pp. 3961-3975, 2004., [email protected]

[10] F. Armero, D. Ehrlich, "Numerical modeling of softening hinges in thin Euler-Bernoulli beams," Computers & Structures , vol. 84, no. 10-11, pp. 641-656, 2006.

[11] D. Ehrlich, F. Armero, "Finite element methods for the analysis of softening plastic hinges in beams and frames," Computational Mechanics , vol. 35, no. 4, pp. 237-264, 2005., [email protected]

[12] W. A. Thanoon, D. K. Paul, M. S. Jaafar, D. N. Trikha, "Influence of torsion on the inelastic response of three-dimensional r.c. frames," Finite Elements in Analysis and Design , vol. 40, no. 5-6, pp. 611-628, 2004., [email protected]

[13] M. R. Horn Plastic Theory of Structures , MIT Press, Cambridge, UK, 1972.

[14] B.G. Neal The Plasctic Methods of Structural Analysis , Chapman and Hall, New York, NY, USA, 1977.

[15] M. S. Park, B. C. Lee, "Geometrically non-linear and elastoplastic three-dimensional shear flexible beam element of von-mises-type hardening material," International Journal for Numerical Methods in Engineering , vol. 39, no. 3, pp. 383-408, 1996.

[16] J. G. Orbison, W. McGuire, J. F. Abel, "Yield surface applications in nonlinear steel frame analysis," Computer Methods in Applied Mechanics and Engineering , vol. 33, no. 1-3, pp. 557-573, 1982.

[17] J. H. Argyris, B. Boni, U. Hindenlang, M. Kleiber, "Finite element analysis of two- and three-dimensional elasto-plastic frames-the natural approach," Computer Methods in Applied Mechanics and Engineering , vol. 35, no. 2, pp. 221-248, 1982.

[18] Y. Ueda , T. Yao, "The plastic node method: a new method of plastic analysis," Computer Methods in Applied Mechanics and Engineering , vol. 34, pp. 1089-1104, 1982.

[19] P. G. Hodge, Jr. Plastic Analysis of Structures , of McGraw-Hill Series in Engineering Sciences, pp. xiv+364, McGraw-Hill, New York, NY, USA, 1959.

[20] A. Haldar, N. K. Ming, "Elasto-plastic large deformation analysis of PR steel frames for LRFD," Computers and Structures , vol. 31, no. 5, pp. 811-823, 1989.

[21] G. Shi, S. N. Atluri, "Elasto-plastic large deformation analysis of space-frames: a plastic-hinge and stress-based explicit derivation of tangent stiffnesses," International Journal for Numerical Methods in Engineering , vol. 26, no. 3, pp. 589-615, 1988.