Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2010 (2010), Article ID 142743, 22 pages
http://dx.doi.org/10.1155/2010/142743
Research Article

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

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

Academic Editor: Giuseppe Rega

Copyright © 2010 William Taylor Matias Silva and Luciano Mendes Bezerra. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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., [58]) 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 Φ=𝛼1||𝐹𝑥||𝐹𝑥𝑝𝛼2+𝛼3||𝐹𝑦||𝐹𝑦𝑝𝛼4+𝛼5||𝐹𝑧||𝐹𝑧𝑝𝛼6+𝛼7||𝑀𝑥||𝑀𝑥𝑝𝛼8+𝛼9||𝑀𝑦||𝑀𝑦𝑝𝛼10+𝛼11||𝑀𝑧||𝑀𝑧𝑝𝛼12+𝛼13||𝐹𝑥||𝐹𝑥𝑝𝛼14||𝑀𝑦||𝑀𝑦𝑝𝛼15+𝛼16||𝐹𝑥||𝐹𝑥𝑝𝛼17||𝑀𝑧||𝑀𝑧𝑝𝛼181=0,(2.1) where |()| is the absolute value of (), 𝐹𝑥 is an axial force, and 𝐹𝑦 and 𝐹𝑧 are shear forces. 𝑀𝑥 is the twisting moment, and 𝑀𝑦 and 𝑀𝑧 are bending moments. 𝐹𝑥𝑝 is the plastic axial force, 𝐹𝑦𝑝 and 𝐹𝑧𝑝 are plastic shear forces, respectively. 𝑀𝑥𝑝 is the plastic twisting moment, and 𝑀𝑦𝑝 and 𝑀𝑧𝑝 are plastic bending moments. The constants 𝛼𝑖 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 𝐹𝑗, can be written as follows:𝜕Φ𝜕𝐹𝑥=𝛼1𝛼2|||𝐹𝛼2𝑥1|||𝐹𝛼2𝑥𝑝+𝛼13𝛼14|||𝐹𝛼14𝑥1|||𝐹𝛼14𝑥𝑝||𝑀𝑦||𝑀𝑦𝑝𝛼15+𝛼16𝛼17|||𝐹𝛼17𝑥1|||𝐹𝛼17𝑥𝑝||𝑀𝑧||𝑀𝑧𝑝𝛼18𝐹sgn𝑥,𝜕Φ𝜕𝐹𝑦=𝛼3𝛼4|||𝐹𝛼4𝑦1|||𝐹𝛼4𝑦𝑝𝐹sgn𝑦;𝜕Φ𝜕𝐹𝑧=𝛼5𝛼6|||𝐹𝛼6𝑧1|||𝐹𝛼6𝑧𝑝𝐹sgn𝑧,𝜕Φ𝜕𝑀𝑥=𝛼7𝛼8|||𝑀𝛼8𝑥1|||𝑀𝛼8𝑥𝑝𝑀sgn𝑥,𝜕Φ𝜕𝑀𝑦=𝛼9𝛼10|||𝑀𝛼10𝑦1|||𝑀𝛼10𝑦𝑝+𝛼13𝛼15||𝐹𝑥||𝐹𝑥𝑝𝛼14|||𝑀𝛼15𝑦1|||𝑀𝛼15𝑦𝑝𝑀sgn𝑦,𝜕Φ𝜕𝑀𝑧=𝛼11𝛼12|||𝑀𝛼12𝑧1|||𝑀𝛼12𝑧𝑝+𝛼16𝛼18||𝐹𝑥||𝐹𝑥𝑝𝛼17|||𝑀𝛼18𝑧1|||𝑀𝛼18𝑧𝑝𝑀sgn𝑧,(2.2) where sgn() denotes the signal of the components of the vector 𝐹𝑗. 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 𝜕Φ𝜕𝐹𝑗𝑇1=𝜕Φ𝜕𝐹𝑥1𝜕Φ𝜕𝐹𝑦1𝜕Φ𝜕𝐹𝑧1𝜕Φ𝜕𝑀𝑥1𝜕Φ𝜕𝑀𝑦1𝜕Φ𝜕𝑀𝑧1,000000𝜕Φ𝜕𝐹𝑗𝑇2=000000𝜕Φ𝜕𝐹𝑥2𝜕Φ𝜕𝐹𝑦2𝜕Φ𝜕𝐹𝑧2𝜕Φ𝜕𝑀𝑥2𝜕Φ𝜕𝑀𝑦2𝜕Φ𝜕𝑀𝑧2.(2.3) It is noted that the sub indexes 𝑖,𝑗,𝑘,𝑙,, 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 𝜕Φ/𝜕𝐹𝑥 with respect to the generalized nodal forces are expressed as 𝜕2Φ𝜕𝐹2𝑥=𝛼1𝛼2𝛼2|||𝐹1𝛼2𝑥2|||𝐹𝛼2𝑥𝑝+𝛼13𝛼14𝛼14|||𝐹1𝛼14𝑥2|||𝐹𝛼14𝑥𝑝||𝑀𝑦||𝑀𝑦𝑝𝛼15+𝛼16𝛼17𝛼17|||𝐹1𝛼17𝑥2|||𝐹𝛼17𝑥𝑝||𝑀𝑧||𝑀𝑧𝑝𝛼18𝜕2Φ𝜕𝐹𝑥𝜕𝐹𝑦𝜕=0,2Φ𝜕𝐹𝑥𝜕𝐹𝑧𝜕=0,2Φ𝜕𝐹𝑥𝜕𝑀𝑥𝜕=0,2Φ𝜕𝐹𝑥𝜕𝑀𝑦=𝛼13𝛼14𝛼15|||𝐹𝛼14𝑥1|||𝐹𝛼14𝑥𝑝|||𝑀𝛼15𝑦1|||𝑀𝛼15𝑦𝑝𝑀sgn𝑦𝐹sgn𝑥,𝜕2Φ𝜕𝐹𝑥𝜕𝑀𝑧=𝛼16𝛼17𝛼18|||𝐹𝛼17𝑥1|||𝐹𝛼17𝑥𝑝|||𝑀𝛼18𝑧1|||𝑀𝛼18𝑧𝑝𝑀sgn𝑧𝐹sgn𝑥.(2.4) In an analogous procedure, the second derivatives of the other components of the vector {𝜕Φ/𝜕𝐹𝑗} 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𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘=1𝜕2Φ𝜕𝐹𝑥1𝜕𝐹𝑥1𝜕2Φ𝜕𝐹𝑥1𝜕𝐹𝑦1𝜕2Φ𝜕𝐹𝑥1𝜕𝐹𝑧1𝜕2Φ𝜕𝐹𝑥1𝜕𝑀𝑥1𝜕2Φ𝜕𝐹𝑥1𝜕𝑀𝑦1𝜕2Φ𝜕𝐹𝑥1𝜕𝑀𝑧1𝜕0000002Φ𝜕𝐹𝑦1𝜕𝐹𝑥1𝜕2Φ𝜕𝐹𝑦1𝜕𝐹𝑦1𝜕2Φ𝜕𝐹𝑦1𝜕𝐹𝑧1𝜕2Φ𝜕𝐹𝑦1𝜕𝑀𝑥1𝜕2Φ𝜕𝐹𝑦1𝜕𝑀𝑦1𝜕2Φ𝜕𝐹𝑦1𝜕𝑀𝑧1𝜕0000002Φ𝜕𝐹𝑧1𝜕𝐹𝑥1𝜕2Φ𝜕𝐹𝑧1𝜕𝐹𝑦1𝜕2Φ𝜕𝐹𝑧1𝜕𝐹𝑧1𝜕2Φ𝜕𝐹𝑧1𝜕𝑀𝑥1𝜕2Φ𝜕𝐹𝑧1𝜕𝑀𝑦1𝜕2Φ𝜕𝐹𝑧1𝜕𝑀𝑧1𝜕0000002Φ𝜕𝑀𝑥1𝜕𝐹𝑥1𝜕2Φ𝜕𝑀𝑥1𝜕𝐹𝑦1𝜕2Φ𝜕𝑀𝑥1𝜕𝐹𝑧1𝜕2Φ𝜕𝑀𝑥1𝜕𝑀𝑥1𝜕2Φ𝜕𝑀𝑥1𝜕𝑀𝑦1𝜕2Φ𝜕𝑀𝑥1𝜕𝑀𝑧1𝜕0000002Φ𝜕𝑀𝑦1𝜕𝐹𝑥1𝜕2Φ𝜕𝑀𝑦1𝜕𝐹𝑦1𝜕2Φ𝜕𝑀𝑦1𝜕𝐹𝑧1𝜕2Φ𝜕𝑀𝑦1𝜕𝑀𝑥1𝜕2Φ𝜕𝑀𝑦1𝜕𝑀𝑦1𝜕2Φ𝜕𝑀𝑦1𝜕𝑀𝑧1𝜕0000002Φ𝜕𝑀𝑧1𝜕𝐹𝑥1𝜕2Φ𝜕𝑀𝑧1𝜕𝐹𝑦1𝜕2Φ𝜕𝑀𝑧1𝜕𝐹𝑧1𝜕2Φ𝜕𝑀𝑧1𝜕𝑀𝑥1𝜕2Φ𝜕𝑀𝑧1𝜕𝑀𝑦1𝜕2Φ𝜕𝑀𝑧1𝜕𝑀𝑧100000000000000000000000000000000000000000000000000000000000000000000000000000012×12,𝜕(2.5)2Φ𝜕𝐹𝑗𝜕𝐹𝑘2=𝜕0000000000000000000000000000000000000000000000000000000000000000000000000000002Φ𝜕𝐹𝑥2𝜕𝐹𝑥2𝜕2Φ𝜕𝐹𝑥2𝜕𝐹𝑦2𝜕2Φ𝜕𝐹𝑥2𝜕𝐹𝑧2𝜕2Φ𝜕𝐹𝑥2𝜕𝑀𝑥2𝜕2Φ𝜕𝐹𝑥2𝜕𝑀𝑦2𝜕2Φ𝜕𝐹𝑥2𝜕𝑀𝑧2𝜕0000002Φ𝜕𝐹𝑦2𝜕𝐹𝑥12𝜕2Φ𝜕𝐹𝑦2𝜕𝐹𝑦2𝜕2Φ𝜕𝐹𝑦2𝜕𝐹𝑧2𝜕2Φ𝜕𝐹𝑦2𝜕𝑀𝑥2𝜕2Φ𝜕𝐹𝑦2𝜕𝑀𝑦2𝜕2Φ𝜕𝐹𝑦2𝜕𝑀𝑧2𝜕0000002Φ𝜕𝐹𝑧2𝜕𝐹𝑥2𝜕2Φ𝜕𝐹𝑧2𝜕𝐹𝑦2𝜕2Φ𝜕𝐹𝑧2𝜕𝐹𝑧2𝜕2Φ𝜕𝐹𝑧2𝜕𝑀𝑥2𝜕2Φ𝜕𝐹𝑧2𝜕𝑀𝑦2𝜕2Φ𝜕𝐹𝑧2𝜕𝑀𝑧2𝜕0000002Φ𝜕𝑀𝑥2𝜕𝐹𝑥2𝜕2Φ𝜕𝑀𝑥2𝜕𝐹𝑦2𝜕2Φ𝜕𝑀𝑥2𝜕𝐹𝑧2𝜕2Φ𝜕𝑀𝑥2𝜕𝑀𝑥2𝜕2Φ𝜕𝑀𝑥2𝜕𝑀𝑦2𝜕2Φ𝜕𝑀𝑥2𝜕𝑀𝑧2𝜕0000002Φ𝜕𝑀𝑦2𝜕𝐹𝑥2𝜕2Φ𝜕𝑀𝑦2𝜕𝐹𝑦2𝜕2Φ𝜕𝑀𝑦2𝜕𝐹𝑧2𝜕2Φ𝜕𝑀𝑦2𝜕𝑀𝑥2𝜕2Φ𝜕𝑀𝑦2𝜕𝑀𝑦2𝜕2Φ𝜕𝑀𝑦2𝜕𝑀𝑧2𝜕0000002Φ𝜕𝑀𝑧2𝜕𝐹𝑥2𝜕2Φ𝜕𝑀𝑧2𝜕𝐹𝑦2𝜕2Φ𝜕𝑀𝑧2𝜕𝐹𝑧2𝜕2Φ𝜕𝑀𝑧2𝜕𝑀𝑥2𝜕2Φ𝜕𝑀𝑧2𝜕𝑀𝑦2𝜕2Φ𝜕𝑀𝑧2𝜕𝑀𝑧212×12.(2.6)

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:𝐹𝑖=𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1,𝜆1>0𝐹trial𝑖=𝐹𝑖+𝐾𝑖𝑗𝑑𝑈𝑗,(3.1a)𝐾𝑖𝑗=𝐸𝐴𝐿00000𝐸𝐴𝐿0000012𝐸𝐼𝑧𝐿30006𝐸𝐼𝑧𝐿2012𝐸𝐼𝑧𝐿30006𝐸𝐼𝑧𝐿212𝐸𝐼𝑦𝐿306𝐸𝐼𝑦𝐿200012𝐸𝐼𝑦𝐿306𝐸𝐼𝑦𝐿20𝐺𝐽𝑥𝐿00000𝐺𝐽𝑥𝐿004𝐸𝐼𝑦𝐿0006𝐸𝐼𝑦𝐿202𝐸𝐼𝑦𝐿04𝐸𝐼𝑧𝐿06𝐸𝐼𝑧𝐿20002𝐸𝐼𝑧𝐿𝐸𝐴𝐿00000sym12𝐸𝐼𝑧𝐿30006𝐸𝐼𝑧𝐿212𝐸𝐼𝑦𝐿306𝐸𝐼𝑦𝐿30𝐺𝐽𝑥𝐿004𝐸𝐼𝑦𝐿04𝐸𝐼𝑧𝐿,(3.1b) where 𝐹𝑖 is the nodal force vector at the last load step that has achieved convergence.The term 𝐾𝑖𝑗𝑑𝑈𝑗 is the elastic predictor vector. 𝐾𝑖𝑗 is the stiffness matrix of the 3D beam element. 𝑑𝑈𝑗 is an increment in the vector 𝑈𝑗 which is defined as the nodal displacement vector. 𝐹trial𝑖 is the elastic trial nodal-force; {𝜕Φ/𝜕𝐹𝑗}1 is the plastic potential defined in the trial nodal-force; 𝜆1 is the plastic multiplier and the term 𝜆1𝐾𝑖𝑗{𝜕Φ/𝜕𝐹𝑗}1 is the plastic corrector. 𝐹𝑖 is the nodal-force after correction. In (3.1b), 𝐿 is the element length. 𝐴 is the cross-section area of the element. 𝐸 is the Young modulus. 𝐺 is the shear modulus. 𝐼𝑦 and 𝐼𝑧 are moments of inertia of the cross-section with respect to 𝑦 axis and 𝑧-axis, respectively. 𝐽𝑥 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).

fig1
Figure 1: Single-vector return algorithm—one plastic hinge.

Generally, the starting guess 𝐹trial𝑖 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 𝑟𝑖 is defined as 𝑟𝑖=𝐹𝑖𝐹𝑖=𝐹𝑖𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1.(3.2) This residual force vector represents the difference between the current force state 𝐹𝑖 and Backward Euler force𝐹𝑖. The trial force 𝐹trial𝑖 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 𝑟new𝑖 in terms of the old residual force vector 𝑟old𝑖, therefore 𝑟new𝑖=𝑟old𝑖+𝑑𝐹𝑖+𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1𝑑𝐹𝑘,(3.3) where 𝑑𝐹𝑖 is an infinitesimal increment of the generalized force vector 𝐹𝑖; 𝑑𝜆1 is the variation of the plastic multiplier 𝜆1 and [𝜕2Φ/𝜕𝐹𝑗𝜕𝐹𝑘]𝑑𝐹𝑘 is the change in the potential plastic vector {𝜕Φ/𝜕𝐹𝑗}1. The goal is to achieve 𝑟new𝑖=0 in (3.3). For that reason 0=𝑟old𝑖+𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1𝑑𝐹𝑘.(3.4) To determine an expression to 𝑑𝐹𝑖 (the correcting force vector) in (3.4), define the square matrix 𝑄𝑖𝑘 as 𝑄𝑖𝑘=𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1.(3.5) Now, considering (3.4) and (3.5) and after some algebraic manipulations to solve (3.4) for 𝑑𝐹𝑖, it follows that𝑑𝐹𝑖=𝑄1𝑖𝑟old+𝑑𝜆1𝐾𝑗𝜕Φ𝜕𝐹𝑗1.(3.6) A first-order Taylor’s series expansion of the yield function Φ around the final nodal-force vector 𝐹𝑖 is applied. This series expansion is necessary to get a linear approximation to the new value of the yield function Φnew, therefore Φnew1=Φold1+𝜕Φ𝜕𝐹𝑖1𝑑𝐹𝑖.(3.7) Since the expression for 𝑑𝐹𝑖 is available in (3.6) and for Φnew1=0, the variation of the plastic multiplier 𝑑𝜆1 is readily found𝑑𝜆1=Φold1𝜕Φ/𝜕𝐹𝑖1𝑄1𝑖𝑟old𝜕Φ/𝜕𝐹𝑖1𝑄1𝑖𝐾𝑗𝜕Φ/𝜕𝐹𝑗1.(3.8) This iterative procedure is continued until the yield criterion Φ=0 is satisfied at the final force state; that is, 𝑟norm=|𝑟𝑖|/|𝐹trial𝑖|<TOL and Φnorm=|Φ|<TOL, where 𝑟norm is a norm for the residual force vector. Φnorm is defined as the residual yield norm and 𝑇𝑂𝐿 is a tolerance for convergence. In this paper, it is assumed that TOL=1010.

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:𝐹𝑖=𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2,with𝜆1>0,𝜆2>0.(3.9) To derive a scheme for 2VRA, the residual force vector 𝑟𝑖 is defined as𝑟𝑖=𝐹𝑖𝐹𝑖=𝐹𝑖𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2.(3.10) This residual force vector represents the difference between the current force state 𝐹𝑖 and the Backward Euler force 𝐹𝑖. Analogous to the 1VRA, the trial force 𝐹trial𝑖 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 𝑟new𝑖 in terms of the old residual force vector 𝑟old𝑖, therefore 𝑟new𝑖=𝑟old𝑖+𝑑𝐹𝑖+𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1𝑑𝐹𝑘+𝑑𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2+𝜆2𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘2𝑑𝐹𝑘.(3.11) Making 𝑟new𝑖=0 in (3.11), it follows that0=𝑟old𝑖+𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2+𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1+𝜆2𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘2𝑑𝐹𝑘.(3.12) To determine an expression to 𝑑𝐹𝑖 (the correcting force vector) in (3.12), define the square matrix 𝑄𝑖𝑘 as 𝑄𝑖𝑘=𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1+𝜆2𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘2.(3.13) Now considering (3.13) and (3.12) and after some algebraic manipulations to solve (3.12) for 𝑑𝐹𝑖, it follows that 𝑑𝐹𝑖=𝑄1𝑖𝑟old+𝑑𝜆1𝐾𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝐾𝑗𝜕Φ𝜕𝐹𝑗2.(3.14) First-order Taylor’s series expansions of the yield function at node 1, Φ1, and at node 2, Φ2, around the final force vector 𝐹𝑖 are now developed. These series expansions generate linear approximations of the new values of the yield functions Φnew1 and Φnew2 as Φnew1=Φold1+𝜕Φ𝜕𝐹𝑖1𝑑𝐹𝑖,Φnew2=Φold2+𝜕Φ𝜕𝐹𝑖2𝑑𝐹𝑖.(3.15) Since the expression for 𝑑𝐹𝑖 is available in (3.14) and the goal is to achieve both Φnew1=0 and Φnew2=0, then Φold1𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝑟old=𝑑𝜆1𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗2,Φold2𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝑟old=𝑑𝜆1𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗2.(3.16) The two previous expressions form a system of equations in 𝑑𝜆1 and 𝑑𝜆2. By introducing auxiliary variables𝑎11𝑎12𝑎21𝑎22𝑑𝜆1𝑑𝜆2=𝑏1𝑏2(3.17) and now using Cramer's rule for the previous system of equations, the variations of the plastic multiplier 𝑑𝜆1 and 𝑑𝜆2 are readily found by the expressions𝑑𝜆1=𝑏1𝑎22𝑏2𝑎12𝑎11𝑎22𝑎12𝑎21,𝑑𝜆2=𝑏2𝑎11𝑏1𝑎21𝑎11𝑎22𝑎12𝑎21,(3.18) where the auxiliary variables are defined as 𝑏1=Φold1𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝑟old,𝑏2=Φold2𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝑟old,𝑎11=𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗1,𝑎12=𝜕Φ𝜕𝐹𝑖1𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗2,𝑎21=𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗1,𝑎22=𝜕Φ𝜕𝐹𝑖2𝑄1𝑖𝐾𝑗𝜕Φ𝜕𝐹𝑗2.(3.19)

fig2
Figure 2: Two-Vectors return algorithm—two plastic hinges.

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 𝑟norm=|𝑟𝑖|/|𝐹trial𝑖|<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 𝐾𝐸𝑃𝑖𝑗 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 𝐹𝑖 is outside the yield surface and vector 𝐹𝑖 is on the yield surface, (3.1a) can be rewritten as 𝐹𝑖=𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1,with𝜆1>0.(4.1) The term 𝑑𝐹trial𝑖=𝐾𝑖𝑗𝑑𝑈𝑗 represents the elastic trial nodal-force increment and is associated with the nodal displacement increment 𝑑𝑈𝑗. The differential of (4.1) gives the following expression; 𝑑𝐹𝑖=𝐾𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1𝑑𝐹𝑘.(4.2) After some algebraic treatment, 𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1𝑑𝐹𝑘=𝐾𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1.(4.3) Using matrix 𝑄𝑖𝑘 as defined by (3.5) and denoting the reduced stiffness matrix as 𝑅𝑖𝑗=𝑄1𝑖𝐾𝑗, (4.3) changes into the following equation:𝑑𝐹𝑖=𝑅𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝜕Φ𝜕𝐹𝑗1.(4.4) It turns out that the form of (4.4) is similar to the nonconsistent form except for the change in 𝐾𝑖𝑗 to 𝑅𝑖𝑗=𝑄1𝑖𝐾𝑗 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 𝜕Φ𝜕𝐹𝑖1𝑑𝐹𝑖=0𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝜕Φ𝜕𝐹𝑗1=0,(4.5) and, therefore, the change in the plastic multiplier 𝑑𝜆1may be expressed as 𝑑𝜆1=𝜕Φ/𝜕𝐹𝑖1𝑅𝑖𝑗𝑑𝑈𝑗𝜕Φ/𝜕𝐹𝑖1𝑅𝑖𝑗𝜕Φ/𝜕𝐹𝑗1.(4.6) Finally, using (4.4) and (4.6), the elastoplastic consistent stiffness matrix 𝐾𝐴𝐿𝑖𝑗may be determined by 𝐾𝐴𝐿𝑖𝑗=𝑅𝑖𝑗𝑅𝑖𝑚𝜕Φ/(𝜕𝐹𝑚)1𝜕Φ/(𝜕𝐹𝑛)1𝑅𝑛𝑗𝜕Φ/(𝜕𝐹𝑚)1𝑅𝑚𝑛𝜕Φ/(𝜕𝐹𝑛)1.(4.7)

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): 𝐹𝑖=𝐹trial𝑖𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2,with𝜆1>0,𝜆2>0.(4.8) Differentiating (4.8) and performing algebraic manipulations, it follows that𝛿𝑖𝑘+𝜆1𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘1+𝜆2𝐾𝑖𝑗𝜕2Φ𝜕𝐹𝑗𝜕𝐹𝑘2𝑑𝐹𝑘=𝐾𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗1𝑑𝜆2𝐾𝑖𝑗𝜕Φ𝜕𝐹𝑗2.(4.9) Using matrix 𝑄𝑖𝑘 as defined by (3.13) and denoting the reduced stiffness matrix as 𝑅𝑖𝑗=𝑄1𝑖𝐾𝑗, we can rewrite the above equation as𝑑𝐹𝑖=𝑅𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝜕Φ𝜕𝐹𝑗1𝑑𝜆2𝜕Φ𝜕𝐹𝑗2.(4.10) 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 𝑑𝐹𝑖 taken from (4.10), one can get𝜕Φ𝜕𝐹𝑖1𝑑𝐹𝑖=0𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝜕Φ𝜕𝐹𝑗1𝑑𝜆2𝜕Φ𝜕𝐹𝑗2=0,𝜕Φ𝜕𝐹𝑖2𝑑𝐹𝑖=0𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝑑𝑈𝑗𝑑𝜆1𝜕Φ𝜕𝐹𝑗1𝑑𝜆2𝜕Φ𝜕𝐹𝑗2=0.(4.11) Alternatively, by collecting terms algebraically, (4.11) are transformed in 𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝑑𝑈𝑗=𝑑𝜆1𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗2,𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝑑𝑈𝑗=𝑑𝜆1𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗1+𝑑𝜆2𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗2.(4.12) The previous two expressions represent a system of equations with 𝑑𝜆1 and 𝑑𝜆2 as unknowns. Introducing auxiliary variables 𝑐1, 𝑐2, 𝑏11, 𝑏12, 𝑏21 and 𝑏22 as𝑐1=𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝑑𝑈𝑗,𝑐2=𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝑑𝑈𝑗,𝑏11=𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗1,𝑏12=𝜕Φ𝜕𝐹𝑖1𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗2,𝑏21=𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗1,𝑏22=𝜕Φ𝜕𝐹𝑖2𝑅𝑖𝑗𝜕Φ𝜕𝐹𝑗2.(4.13) With the help of Cramer’s rule, the system in (4.12) is solved for 𝑑𝜆1 and 𝑑𝜆2𝑑𝜆1=𝑐1𝑏22𝑐2𝑏12𝑏11𝑏22𝑏12𝑏21,𝑑𝜆2=𝑐2𝑏11𝑐1𝑏21𝑏11𝑏22𝑏12𝑏21(4.14) or alternatively 𝑑𝜆1=𝑏22𝜕Φ/𝜕𝐹𝑚1𝑏12𝜕Φ/𝜕𝐹𝑚2𝑏11𝑏22𝑏12𝑏21𝑅𝑚𝑛𝑑𝑈𝑛,𝑑𝜆2=𝑏11𝜕Φ/𝜕𝐹𝑚2𝑏21𝜕Φ/𝜕𝐹𝑚1𝑏11𝑏22𝑏12𝑏21𝑅𝑚𝑛𝑑𝑈𝑛.(4.15) If any plastic multiplier is negative, that is, 𝑑𝜆1<0, or 𝑑𝜆2<0, then this plastic multiplier is set to zero (i.e., 𝑑𝜆1=0 or 𝑑𝜆2=0) which corresponds to an initial plastic hinge. Again, introducing new auxiliary variables𝛽1=𝑏22𝑏11𝑏22𝑏12𝑏21,𝛽2=𝑏12𝑏11𝑏22𝑏12𝑏21,𝛽3=𝑏11𝑏11𝑏22𝑏12𝑏21,𝛽4=𝑏21𝑏11𝑏22𝑏12𝑏21,(4.16) Equation (4.15) can be rewritten as 𝑑𝜆1=𝛽1𝜕Φ𝜕𝐹𝑚1𝑅𝑚𝑛𝛽2𝜕Φ𝜕𝐹𝑚2𝑅𝑚𝑛𝑑𝑈𝑛,𝑑𝜆2=𝛽3𝜕Φ𝜕𝐹𝑚2𝑅𝑚𝑛𝛽4𝜕Φ𝜕𝐹𝑚1𝑅𝑚𝑛𝑑𝑈𝑛.(4.17) Finally, using (4.10) and (4.17), the elastoplastic consistent stiffness matrix 𝐾𝐴𝐿𝑖𝑗is obtained by the following expression:𝐾𝐴𝐿𝑖𝑗=𝑅𝑖𝑗𝛽1𝑅𝑖𝑚𝜕Φ𝜕𝐹𝑚1𝜕Φ𝜕𝐹𝑛1𝑅𝑛𝑗𝛽2𝑅𝑖𝑚𝜕Φ𝜕𝐹𝑚1𝜕Φ𝜕𝐹𝑛2𝑅𝑛𝑗𝛽3𝑅𝑖𝑚𝜕Φ𝜕𝐹𝑚2𝜕Φ𝜕𝐹𝑛2𝑅𝑛𝑗𝛽4𝑅𝑖𝑚𝜕Φ𝜕𝐹𝑚2𝜕Φ𝜕𝐹𝑛1𝑅𝑛𝑗.(4.18)

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=1010. 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 [13]. 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 “𝑛” 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 (𝑛−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 𝑃 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].

fig3
Figure 3: 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.

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.

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

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.

tab3
Table 3: Error norms for Newton iterative scheme on global equilibrium.
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].

fig4
Figure 4: 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.

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.

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

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.

tab5
Table 5: Error norms for Newton iterative scheme on global equilibrium.
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].

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

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=1010) parameter. The same observations can be verified in the examples studied before.

tab6
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.
tab7
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.
tab8
Table 8: Four-legged jacket type-1. Error norms for Newton iterative scheme on global equilibrium.
tab9
Table 9: Four-legged jacket type-2. Error norms for Newton iterative scheme on global equilibrium.

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 and T. J. R. Hughes, Computational Inelasticity, vol. 7 of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 1998. View at Zentralblatt MATH
  2. M. A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, vol. 1 of Essentials, John Wiley & Sons, Chichester, UK, 1991. View at Zentralblatt MATH
  3. M. A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, vol. 2 of Advanced Topics, John Wiley & Sons, Chichester, UK, 1997. View at Zentralblatt MATH
  4. I. Doltsinis, Elements of Plasticity: Theory and Computation, vol. 1 of High Performance Structures and Materials, WIT, Southampton, UK, 2000. View at Zentralblatt MATH
  5. A. S. Gendy and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  6. F. G. A. Al-Bermani and 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. View at Publisher · View at Google Scholar · View at Scopus
  7. B. Skallerud and 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. View at Google Scholar
  8. S. Krenk, C. Vissing-Jørgensen, and L. Thesbjerg, “Efficient collapse analysis techniques for framed structures,” Computers and Structures, vol. 72, no. 4, pp. 481–496, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  9. M. E. Marante, R. Picón, and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  10. F. Armero and D. Ehrlich, “Numerical modeling of softening hinges in thin Euler-Bernoulli beams,” Computers & Structures, vol. 84, no. 10-11, pp. 641–656, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  11. D. Ehrlich and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  12. W. A. Thanoon, D. K. Paul, M. S. Jaafar, and 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. View at Publisher · View at Google Scholar · View at Scopus
  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 and 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. View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  16. J. G. Orbison, W. McGuire, and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  17. J. H. Argyris, B. Boni, U. Hindenlang, and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  18. Y. Ueda and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  19. P. G. Hodge, Jr., Plastic Analysis of Structures, McGraw-Hill Series in Engineering Sciences, McGraw-Hill, New York, NY, USA, 1959.
  20. A. Haldar and 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. View at Publisher · View at Google Scholar · View at Scopus
  21. G. Shi and 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus