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

Authors: William Taylor Matias Silva; Luciano Mendes Bezerra
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 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
(2.1)Φ=α1(|Fx|Fxp)α2+α3(|Fy|Fyp)α4+α5(|Fz|Fzp)α6+α7(|Mx|Mxp)α8+α9(|My|Myp)α10+α11(|Mz|Mzp)α12+α13(|Fx|Fxp)α14(|My|Myp)α15+α16(|Fx|Fxp)α17(|Mz|Mzp)α18-1=0,
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:(2.2)∂Φ∂Fx=[α1α2|Fxα2-1|Fxpα2+α13α14|Fxα14-1|Fxpα14(|My|Myp)α15+α16α17|Fxα17-1|Fxpα17(|Mz|Mzp)α18]sgn(Fx),∂Φ∂Fy=α3α4|Fyα4-1|Fypα4sgn(Fy);∂Φ∂Fz=α5α6|Fzα6-1|Fzpα6sgn(Fz),∂Φ∂Mx=α7α8|Mxα8-1|Mxpα8sgn(Mx),∂Φ∂My=[α9α10|Myα10-1|Mypα10+α13α15(|Fx|Fxp)α14|Myα15-1|Mypα15]sgn(My),∂Φ∂Mz=[α11α12|Mzα12-1|Mzpα12+α16α18(|Fx|Fxp)α17|Mzα18-1|Mzpα18]sgn(Mz),
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
(2.3){∂Φ∂Fj}1T={∂Φ∂Fx1∂Φ∂Fy1∂Φ∂Fz1∂Φ∂Mx1∂Φ∂My1∂Φ∂Mz1000000},{∂Φ∂Fj}2T={000000∂Φ∂Fx2∂Φ∂Fy2∂Φ∂Fz2∂Φ∂Mx2∂Φ∂My2∂Φ∂Mz2}.
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
(2.4)∂2Φ∂Fx2=α1α2(α2-1)|Fxα2-2|Fxpα2+α13α14(α14-1)|Fxα14-2|Fxpα14(|My|Myp)α15+α16α17(α17-1)|Fxα17-2|Fxpα17(|Mz|Mzp)α18∂2Φ∂Fx∂Fy=0,∂2Φ∂Fx∂Fz=0,∂2Φ∂Fx∂Mx=0,∂2Φ∂Fx∂My=α13α14α15(|Fxα14-1|Fxpα14)(|Myα15-1|Mypα15)sgn(My)sgn(Fx),∂2Φ∂Fx∂Mz=α16α17α18(|Fxα17-1|Fxpα17)(|Mzα18-1|Mzpα18)sgn(Mz)sgn(Fx).
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(2.5)[∂2Φ∂Fj∂Fk]=1[∂2Φ∂Fx1∂Fx1∂2Φ∂Fx1∂Fy1∂2Φ∂Fx1∂Fz1∂2Φ∂Fx1∂Mx1∂2Φ∂Fx1∂My1∂2Φ∂Fx1∂Mz1000000∂2Φ∂Fy1∂Fx1∂2Φ∂Fy1∂Fy1∂2Φ∂Fy1∂Fz1∂2Φ∂Fy1∂Mx1∂2Φ∂Fy1∂My1∂2Φ∂Fy1∂Mz1000000∂2Φ∂Fz1∂Fx1∂2Φ∂Fz1∂Fy1∂2Φ∂Fz1∂Fz1∂2Φ∂Fz1∂Mx1∂2Φ∂Fz1∂My1∂2Φ∂Fz1∂Mz1000000∂2Φ∂Mx1∂Fx1∂2Φ∂Mx1∂Fy1∂2Φ∂Mx1∂Fz1∂2Φ∂Mx1∂Mx1∂2Φ∂Mx1∂My1∂2Φ∂Mx1∂Mz1000000∂2Φ∂My1∂Fx1∂2Φ∂My1∂Fy1∂2Φ∂My1∂Fz1∂2Φ∂My1∂Mx1∂2Φ∂My1∂My1∂2Φ∂My1∂Mz1000000∂2Φ∂Mz1∂Fx1∂2Φ∂Mz1∂Fy1∂2Φ∂Mz1∂Fz1∂2Φ∂Mz1∂Mx1∂2Φ∂Mz1∂My1∂2Φ∂Mz1∂Mz1000000000000000000000000000000000000000000000000000000000000000000000000000000]12×12,(2.6)[∂2Φ∂Fj∂Fk]2=[000000000000000000000000000000000000000000000000000000000000000000000000000000∂2Φ∂Fx2∂Fx2∂2Φ∂Fx2∂Fy2∂2Φ∂Fx2∂Fz2∂2Φ∂Fx2∂Mx2∂2Φ∂Fx2∂My2∂2Φ∂Fx2∂Mz2000000∂2Φ∂Fy2∂Fx12∂2Φ∂Fy2∂Fy2∂2Φ∂Fy2∂Fz2∂2Φ∂Fy2∂Mx2∂2Φ∂Fy2∂My2∂2Φ∂Fy2∂Mz2000000∂2Φ∂Fz2∂Fx2∂2Φ∂Fz2∂Fy2∂2Φ∂Fz2∂Fz2∂2Φ∂Fz2∂Mx2∂2Φ∂Fz2∂My2∂2Φ∂Fz2∂Mz2000000∂2Φ∂Mx2∂Fx2∂2Φ∂Mx2∂Fy2∂2Φ∂Mx2∂Fz2∂2Φ∂Mx2∂Mx2∂2Φ∂Mx2∂My2∂2Φ∂Mx2∂Mz2000000∂2Φ∂My2∂Fx2∂2Φ∂My2∂Fy2∂2Φ∂My2∂Fz2∂2Φ∂My2∂Mx2∂2Φ∂My2∂My2∂2Φ∂My2∂Mz2000000∂2Φ∂Mz2∂Fx2∂2Φ∂Mz2∂Fy2∂2Φ∂Mz2∂Fz2∂2Φ∂Mz2∂Mx2∂2Φ∂Mz2∂My2∂2Φ∂Mz2∂Mz2]12×12.


## 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:(2.2)∂Φ∂Fx=[α1α2|Fxα2-1|Fxpα2+α13α14|Fxα14-1|Fxpα14(|My|Myp)α15+α16α17|Fxα17-1|Fxpα17(|Mz|Mzp)α18]sgn(Fx),∂Φ∂Fy=α3α4|Fyα4-1|Fypα4sgn(Fy);∂Φ∂Fz=α5α6|Fzα6-1|Fzpα6sgn(Fz),∂Φ∂Mx=α7α8|Mxα8-1|Mxpα8sgn(Mx),∂Φ∂My=[α9α10|Myα10-1|Mypα10+α13α15(|Fx|Fxp)α14|Myα15-1|Mypα15]sgn(My),∂Φ∂Mz=[α11α12|Mzα12-1|Mzpα12+α16α18(|Fx|Fxp)α17|Mzα18-1|Mzpα18]sgn(Mz),
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
(2.3){∂Φ∂Fj}1T={∂Φ∂Fx1∂Φ∂Fy1∂Φ∂Fz1∂Φ∂Mx1∂Φ∂My1∂Φ∂Mz1000000},{∂Φ∂Fj}2T={000000∂Φ∂Fx2∂Φ∂Fy2∂Φ∂Fz2∂Φ∂Mx2∂Φ∂My2∂Φ∂Mz2}.
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
(2.4)∂2Φ∂Fx2=α1α2(α2-1)|Fxα2-2|Fxpα2+α13α14(α14-1)|Fxα14-2|Fxpα14(|My|Myp)α15+α16α17(α17-1)|Fxα17-2|Fxpα17(|Mz|Mzp)α18∂2Φ∂Fx∂Fy=0,∂2Φ∂Fx∂Fz=0,∂2Φ∂Fx∂Mx=0,∂2Φ∂Fx∂My=α13α14α15(|Fxα14-1|Fxpα14)(|Myα15-1|Mypα15)sgn(My)sgn(Fx),∂2Φ∂Fx∂Mz=α16α17α18(|Fxα17-1|Fxpα17)(|Mzα18-1|Mzpα18)sgn(Mz)sgn(Fx).
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(2.5)[∂2Φ∂Fj∂Fk]=1[∂2Φ∂Fx1∂Fx1∂2Φ∂Fx1∂Fy1∂2Φ∂Fx1∂Fz1∂2Φ∂Fx1∂Mx1∂2Φ∂Fx1∂My1∂2Φ∂Fx1∂Mz1000000∂2Φ∂Fy1∂Fx1∂2Φ∂Fy1∂Fy1∂2Φ∂Fy1∂Fz1∂2Φ∂Fy1∂Mx1∂2Φ∂Fy1∂My1∂2Φ∂Fy1∂Mz1000000∂2Φ∂Fz1∂Fx1∂2Φ∂Fz1∂Fy1∂2Φ∂Fz1∂Fz1∂2Φ∂Fz1∂Mx1∂2Φ∂Fz1∂My1∂2Φ∂Fz1∂Mz1000000∂2Φ∂Mx1∂Fx1∂2Φ∂Mx1∂Fy1∂2Φ∂Mx1∂Fz1∂2Φ∂Mx1∂Mx1∂2Φ∂Mx1∂My1∂2Φ∂Mx1∂Mz1000000∂2Φ∂My1∂Fx1∂2Φ∂My1∂Fy1∂2Φ∂My1∂Fz1∂2Φ∂My1∂Mx1∂2Φ∂My1∂My1∂2Φ∂My1∂Mz1000000∂2Φ∂Mz1∂Fx1∂2Φ∂Mz1∂Fy1∂2Φ∂Mz1∂Fz1∂2Φ∂Mz1∂Mx1∂2Φ∂Mz1∂My1∂2Φ∂Mz1∂Mz1000000000000000000000000000000000000000000000000000000000000000000000000000000]12×12,(2.6)[∂2Φ∂Fj∂Fk]2=[000000000000000000000000000000000000000000000000000000000000000000000000000000∂2Φ∂Fx2∂Fx2∂2Φ∂Fx2∂Fy2∂2Φ∂Fx2∂Fz2∂2Φ∂Fx2∂Mx2∂2Φ∂Fx2∂My2∂2Φ∂Fx2∂Mz2000000∂2Φ∂Fy2∂Fx12∂2Φ∂Fy2∂Fy2∂2Φ∂Fy2∂Fz2∂2Φ∂Fy2∂Mx2∂2Φ∂Fy2∂My2∂2Φ∂Fy2∂Mz2000000∂2Φ∂Fz2∂Fx2∂2Φ∂Fz2∂Fy2∂2Φ∂Fz2∂Fz2∂2Φ∂Fz2∂Mx2∂2Φ∂Fz2∂My2∂2Φ∂Fz2∂Mz2000000∂2Φ∂Mx2∂Fx2∂2Φ∂Mx2∂Fy2∂2Φ∂Mx2∂Fz2∂2Φ∂Mx2∂Mx2∂2Φ∂Mx2∂My2∂2Φ∂Mx2∂Mz2000000∂2Φ∂My2∂Fx2∂2Φ∂My2∂Fy2∂2Φ∂My2∂Fz2∂2Φ∂My2∂Mx2∂2Φ∂My2∂My2∂2Φ∂My2∂Mz2000000∂2Φ∂Mz2∂Fx2∂2Φ∂Mz2∂Fy2∂2Φ∂Mz2∂Fz2∂2Φ∂Mz2∂Mx2∂2Φ∂Mz2∂My2∂2Φ∂Mz2∂Mz2]12×12.


## 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:(3.1a)F̂i=Fitrial-λ1Kij{∂Φ∂Fj}1,λ1>0Fitrial=F¯i+KijdUj,(3.1b)Kij=[EAL00000-EAL0000012EIzL30006EIzL20-12EIzL30006EIzL212EIyL30-6EIyL2000-12EIyL30-6EIyL20GJxL00000-GJxL004EIyL0006EIyL202EIyL04EIzL0-6EIzL20002EIzLEAL00000sym12EIzL3000-6EIzL212EIyL306EIyL30GJxL004EIyL04EIzL], where F̅i is the nodal force vector at the last load step that has achieved convergence.The term KijdUj is the elastic predictor vector. Kij is the stiffness matrix of the 3D beam element. dUj 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)(b)Generally, the starting guessFitrial 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 vectorri is defined as
(3.2)ri=Fi-F̂i=Fi-(Fitrial-λ1Kij{∂Φ∂Fj}1).
This residual force vector represents the difference between the current force state Fi and Backward Euler forceF̂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
(3.3)rinew=riold+dFi+dλ1Kij{∂Φ∂Fj}1+λ1Kij[∂2Φ∂Fj∂Fk]1dFk,
where dFi is an infinitesimal increment of the generalized force vector Fi; dλ1 is the variation of the plastic multiplier λ1 and [∂2Φ/∂Fj∂Fk]dFk is the change in the potential plastic vector {∂Φ/∂Fj}1. The goal is to achieve rinew=0 in (3.3). For that reason
(3.4)0=riold+dλ1Kij{∂Φ∂Fj}1+(δik+λ1Kij[∂2Φ∂Fj∂Fk]1)dFk.
To determine an expression to dFi (the correcting force vector) in (3.4), define the square matrix Qik as
(3.5)Qik=δik+λ1Kij[∂2Φ∂Fj∂Fk]1.
Now, considering (3.4) and (3.5) and after some algebraic manipulations to solve (3.4) for dFi, it follows that(3.6)dFi=-Qiℓ-1(rℓold+dλ1Kℓj{∂Φ∂Fj}1).
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
(3.7)Φ1new=Φ1old+{∂Φ∂Fi}1dFi.
Since the expression for dFi is available in (3.6) and for Φ1new=0, the variation of the plastic multiplier dλ1 is readily found(3.8)dλ1=Φ1old-{∂Φ/∂Fi}1Qiℓ-1rℓold{∂Φ/∂Fi}1Qiℓ-1Kℓj{∂Φ/∂Fj}1.
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=10-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 Figure2. The Backward Euler force obtained with the 2VRA may be defined by the following expression:(3.9)F̂i=Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2,withλ1>0,λ2>0.
To derive a scheme for 2VRA, the residual force vector ri is defined as(3.10)ri=Fi-F̂i=Fi-(Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2).
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
(3.11)rinew=riold+dFi+dλ1Kij{∂Φ∂Fj}1+λ1Kij[∂2Φ∂Fj∂Fk]1dFk+dλ2Kij{∂Φ∂Fj}2+λ2Kij[∂2Φ∂Fj∂Fk]2dFk.
Making rinew=0 in  (3.11), it follows that(3.12)0=riold+dλ1Kij{∂Φ∂Fj}1+dλ2Kij{∂Φ∂Fj}2+(δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2)dFk.
To determine an expression to dFi (the correcting force vector) in (3.12), define the square matrix Qik as (3.13)Qik=δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2.
Now considering (3.13) and (3.12) and after some algebraic manipulations to solve (3.12) for dFi, it follows that
(3.14)dFi=-Qiℓ-1(rℓold+dλ1Kℓj{∂Φ∂Fj}1+dλ2Kℓj{∂Φ∂Fj}2).
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
(3.15)Φ1new=Φ1old+{∂Φ∂Fi}1dFi,Φ2new=Φ2old+{∂Φ∂Fi}2dFi.
Since the expression for dFi is available in (3.14) and the goal is to achieve both Φ1new=0 and Φ2new=0, then
(3.16)Φ1old-{∂Φ∂Fi}1Qiℓ-1rℓold=dλ1{∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}2,Φ2old-{∂Φ∂Fi}2Qiℓ-1rℓold=dλ1{∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}2.
The two previous expressions form a system of equations in dλ1 and dλ2. By introducing auxiliary variables(3.17)[a11a12a21a22]{dλ1dλ2}={b1b2}
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(3.18)dλ1=b1a22-b2a12a11a22-a12a21,dλ2=b2a11-b1a21a11a22-a12a21,
where the auxiliary variables are defined as
(3.19)b1=Φ1old-{∂Φ∂Fi}1Qiℓ-1rℓold,b2=Φ2old-{∂Φ∂Fi}2Qiℓ-1rℓold,a11={∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}1,a12={∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}2,a21={∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}1,a22={∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}2.Two-Vectors return algorithm—two plastic hinges.
(a)(b)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.


## 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 vectorri is defined as
(3.2)ri=Fi-F̂i=Fi-(Fitrial-λ1Kij{∂Φ∂Fj}1).
This residual force vector represents the difference between the current force state Fi and Backward Euler forceF̂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
(3.3)rinew=riold+dFi+dλ1Kij{∂Φ∂Fj}1+λ1Kij[∂2Φ∂Fj∂Fk]1dFk,
where dFi is an infinitesimal increment of the generalized force vector Fi; dλ1 is the variation of the plastic multiplier λ1 and [∂2Φ/∂Fj∂Fk]dFk is the change in the potential plastic vector {∂Φ/∂Fj}1. The goal is to achieve rinew=0 in (3.3). For that reason
(3.4)0=riold+dλ1Kij{∂Φ∂Fj}1+(δik+λ1Kij[∂2Φ∂Fj∂Fk]1)dFk.
To determine an expression to dFi (the correcting force vector) in (3.4), define the square matrix Qik as
(3.5)Qik=δik+λ1Kij[∂2Φ∂Fj∂Fk]1.
Now, considering (3.4) and (3.5) and after some algebraic manipulations to solve (3.4) for dFi, it follows that(3.6)dFi=-Qiℓ-1(rℓold+dλ1Kℓj{∂Φ∂Fj}1).
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
(3.7)Φ1new=Φ1old+{∂Φ∂Fi}1dFi.
Since the expression for dFi is available in (3.6) and for Φ1new=0, the variation of the plastic multiplier dλ1 is readily found(3.8)dλ1=Φ1old-{∂Φ/∂Fi}1Qiℓ-1rℓold{∂Φ/∂Fi}1Qiℓ-1Kℓj{∂Φ/∂Fj}1.
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=10-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 Figure2. The Backward Euler force obtained with the 2VRA may be defined by the following expression:(3.9)F̂i=Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2,withλ1>0,λ2>0.
To derive a scheme for 2VRA, the residual force vector ri is defined as(3.10)ri=Fi-F̂i=Fi-(Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2).
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
(3.11)rinew=riold+dFi+dλ1Kij{∂Φ∂Fj}1+λ1Kij[∂2Φ∂Fj∂Fk]1dFk+dλ2Kij{∂Φ∂Fj}2+λ2Kij[∂2Φ∂Fj∂Fk]2dFk.
Making rinew=0 in  (3.11), it follows that(3.12)0=riold+dλ1Kij{∂Φ∂Fj}1+dλ2Kij{∂Φ∂Fj}2+(δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2)dFk.
To determine an expression to dFi (the correcting force vector) in (3.12), define the square matrix Qik as (3.13)Qik=δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2.
Now considering (3.13) and (3.12) and after some algebraic manipulations to solve (3.12) for dFi, it follows that
(3.14)dFi=-Qiℓ-1(rℓold+dλ1Kℓj{∂Φ∂Fj}1+dλ2Kℓj{∂Φ∂Fj}2).
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
(3.15)Φ1new=Φ1old+{∂Φ∂Fi}1dFi,Φ2new=Φ2old+{∂Φ∂Fi}2dFi.
Since the expression for dFi is available in (3.14) and the goal is to achieve both Φ1new=0 and Φ2new=0, then
(3.16)Φ1old-{∂Φ∂Fi}1Qiℓ-1rℓold=dλ1{∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}2,Φ2old-{∂Φ∂Fi}2Qiℓ-1rℓold=dλ1{∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}2.
The two previous expressions form a system of equations in dλ1 and dλ2. By introducing auxiliary variables(3.17)[a11a12a21a22]{dλ1dλ2}={b1b2}
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(3.18)dλ1=b1a22-b2a12a11a22-a12a21,dλ2=b2a11-b1a21a11a22-a12a21,
where the auxiliary variables are defined as
(3.19)b1=Φ1old-{∂Φ∂Fi}1Qiℓ-1rℓold,b2=Φ2old-{∂Φ∂Fi}2Qiℓ-1rℓold,a11={∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}1,a12={∂Φ∂Fi}1Qiℓ-1Kℓj{∂Φ∂Fj}2,a21={∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}1,a22={∂Φ∂Fi}2Qiℓ-1Kℓj{∂Φ∂Fj}2.Two-Vectors return algorithm—two plastic hinges.
(a)(b)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 matrixKijEP 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
(4.1)Fi=Fitrial-λ1Kij{∂Φ∂Fj}1,withλ1>0.
The term dFitrial=KijdUj represents the elastic trial nodal-force increment and is associated with the nodal displacement increment dUj. The differential of (4.1) gives the following expression;
(4.2)dFi=KijdUj-dλ1Kij{∂Φ∂Fj}1-λ1Kij[∂2Φ∂Fj∂Fk]1dFk.
After some algebraic treatment,
(4.3)(δik+λ1Kij[∂2Φ∂Fj∂Fk]1)dFk=KijdUj-dλ1Kij{∂Φ∂Fj}1.
Using matrix Qik as defined by (3.5) and denoting the reduced stiffness matrix as Rij=Qiℓ-1Kℓj, (4.3) changes into the following equation:(4.4)dFi=Rij(dUj-dλ1{∂Φ∂Fj}1).
It turns out that the form of (4.4) is similar to the nonconsistent form except for the change in Kij to Rij=Qiℓ-1Kℓ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
(4.5){∂Φ∂Fi}1dFi=0⇒{∂Φ∂Fi}1Rij(dUj-dλ1{∂Φ∂Fj}1)=0,
and, therefore, the change in the plastic multiplier dλ1may be expressed as
(4.6)dλ1={∂Φ/∂Fi}1RijdUj{∂Φ/∂Fi}1Rij{∂Φ/∂Fj}1.
Finally, using (4.4) and (4.6), the elastoplastic consistent stiffness matrix KijALmay be determined by
(4.7)KijAL=Rij-Rim{∂Φ/(∂Fm)}1{∂Φ/(∂Fn)}1Rnj{∂Φ/(∂Fm)}1Rmn{∂Φ/(∂Fn)}1.


### 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):
(4.8)Fi=Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2,withλ1>0,λ2>0.
Differentiating (4.8) and performing algebraic manipulations, it follows that(4.9)(δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2)dFk=KijdUj-dλ1Kij{∂Φ∂Fj}1-dλ2Kij{∂Φ∂Fj}2.
Using matrix Qik as defined by (3.13) and denoting the reduced stiffness matrix as Rij=Qiℓ-1Kℓj, we can rewrite the above equation as(4.10)dFi=Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2).
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 dFi taken from (4.10), one can get(4.11){∂Φ∂Fi}1dFi=0⇒{∂Φ∂Fi}1Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2)=0,{∂Φ∂Fi}2dFi=0⇒{∂Φ∂Fi}2Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2)=0.
Alternatively, by collecting terms algebraically, (4.11) are transformed in
(4.12){∂Φ∂Fi}1RijdUj=dλ1{∂Φ∂Fi}1Rij{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}1Rij{∂Φ∂Fj}2,{∂Φ∂Fi}2RijdUj=dλ1{∂Φ∂Fi}2Rij{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}2Rij{∂Φ∂Fj}2.
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(4.13)c1={∂Φ∂Fi}1RijdUj,c2={∂Φ∂Fi}2RijdUj,b11={∂Φ∂Fi}1Rij{∂Φ∂Fj}1,b12={∂Φ∂Fi}1Rij{∂Φ∂Fj}2,b21={∂Φ∂Fi}2Rij{∂Φ∂Fj}1,b22={∂Φ∂Fi}2Rij{∂Φ∂Fj}2.
With the help of Cramer’s rule, the system in  (4.12) is solved for dλ1 and dλ2(4.14)dλ1=c1b22-c2b12b11b22-b12b21,dλ2=c2b11-c1b21b11b22-b12b21
or alternatively
(4.15)dλ1=(b22{∂Φ/∂Fm}1-b12{∂Φ/∂Fm}2b11b22-b12b21)RmndUn,dλ2=(b11{∂Φ/∂Fm}2-b21{∂Φ/∂Fm}1b11b22-b12b21)RmndUn.
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(4.16)β1=b22b11b22-b12b21,β2=b12b11b22-b12b21,β3=b11b11b22-b12b21,β4=b21b11b22-b12b21,
Equation (4.15) can be rewritten as
(4.17)dλ1=(β1{∂Φ∂Fm}1Rmn-β2{∂Φ∂Fm}2Rmn)dUn,dλ2=(β3{∂Φ∂Fm}2Rmn-β4{∂Φ∂Fm}1Rmn)dUn.
Finally, using  (4.10) and (4.17), the elastoplastic consistent stiffness matrix KijALis obtained by the following expression:(4.18)KijAL=Rij-(β1Rim{∂Φ∂Fm}1{∂Φ∂Fn}1Rnj-β2Rim{∂Φ∂Fm}1{∂Φ∂Fn}2Rnj)-(β3Rim{∂Φ∂Fm}2{∂Φ∂Fn}2Rnj-β4Rim{∂Φ∂Fm}2{∂Φ∂Fn}1Rnj).


## 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
(4.1)Fi=Fitrial-λ1Kij{∂Φ∂Fj}1,withλ1>0.
The term dFitrial=KijdUj represents the elastic trial nodal-force increment and is associated with the nodal displacement increment dUj. The differential of (4.1) gives the following expression;
(4.2)dFi=KijdUj-dλ1Kij{∂Φ∂Fj}1-λ1Kij[∂2Φ∂Fj∂Fk]1dFk.
After some algebraic treatment,
(4.3)(δik+λ1Kij[∂2Φ∂Fj∂Fk]1)dFk=KijdUj-dλ1Kij{∂Φ∂Fj}1.
Using matrix Qik as defined by (3.5) and denoting the reduced stiffness matrix as Rij=Qiℓ-1Kℓj, (4.3) changes into the following equation:(4.4)dFi=Rij(dUj-dλ1{∂Φ∂Fj}1).
It turns out that the form of (4.4) is similar to the nonconsistent form except for the change in Kij to Rij=Qiℓ-1Kℓ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
(4.5){∂Φ∂Fi}1dFi=0⇒{∂Φ∂Fi}1Rij(dUj-dλ1{∂Φ∂Fj}1)=0,
and, therefore, the change in the plastic multiplier dλ1may be expressed as
(4.6)dλ1={∂Φ/∂Fi}1RijdUj{∂Φ/∂Fi}1Rij{∂Φ/∂Fj}1.
Finally, using (4.4) and (4.6), the elastoplastic consistent stiffness matrix KijALmay be determined by
(4.7)KijAL=Rij-Rim{∂Φ/(∂Fm)}1{∂Φ/(∂Fn)}1Rnj{∂Φ/(∂Fm)}1Rmn{∂Φ/(∂Fn)}1.


## 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):
(4.8)Fi=Fitrial-λ1Kij{∂Φ∂Fj}1-λ2Kij{∂Φ∂Fj}2,withλ1>0,λ2>0.
Differentiating (4.8) and performing algebraic manipulations, it follows that(4.9)(δik+λ1Kij[∂2Φ∂Fj∂Fk]1+λ2Kij[∂2Φ∂Fj∂Fk]2)dFk=KijdUj-dλ1Kij{∂Φ∂Fj}1-dλ2Kij{∂Φ∂Fj}2.
Using matrix Qik as defined by (3.13) and denoting the reduced stiffness matrix as Rij=Qiℓ-1Kℓj, we can rewrite the above equation as(4.10)dFi=Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2).
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 dFi taken from (4.10), one can get(4.11){∂Φ∂Fi}1dFi=0⇒{∂Φ∂Fi}1Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2)=0,{∂Φ∂Fi}2dFi=0⇒{∂Φ∂Fi}2Rij(dUj-dλ1{∂Φ∂Fj}1-dλ2{∂Φ∂Fj}2)=0.
Alternatively, by collecting terms algebraically, (4.11) are transformed in
(4.12){∂Φ∂Fi}1RijdUj=dλ1{∂Φ∂Fi}1Rij{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}1Rij{∂Φ∂Fj}2,{∂Φ∂Fi}2RijdUj=dλ1{∂Φ∂Fi}2Rij{∂Φ∂Fj}1+dλ2{∂Φ∂Fi}2Rij{∂Φ∂Fj}2.
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(4.13)c1={∂Φ∂Fi}1RijdUj,c2={∂Φ∂Fi}2RijdUj,b11={∂Φ∂Fi}1Rij{∂Φ∂Fj}1,b12={∂Φ∂Fi}1Rij{∂Φ∂Fj}2,b21={∂Φ∂Fi}2Rij{∂Φ∂Fj}1,b22={∂Φ∂Fi}2Rij{∂Φ∂Fj}2.
With the help of Cramer’s rule, the system in  (4.12) is solved for dλ1 and dλ2(4.14)dλ1=c1b22-c2b12b11b22-b12b21,dλ2=c2b11-c1b21b11b22-b12b21
or alternatively
(4.15)dλ1=(b22{∂Φ/∂Fm}1-b12{∂Φ/∂Fm}2b11b22-b12b21)RmndUn,dλ2=(b11{∂Φ/∂Fm}2-b21{∂Φ/∂Fm}1b11b22-b12b21)RmndUn.
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(4.16)β1=b22b11b22-b12b21,β2=b12b11b22-b12b21,β3=b11b11b22-b12b21,β4=b21b11b22-b12b21,
Equation (4.15) can be rewritten as
(4.17)dλ1=(β1{∂Φ∂Fm}1Rmn-β2{∂Φ∂Fm}2Rmn)dUn,dλ2=(β3{∂Φ∂Fm}2Rmn-β4{∂Φ∂Fm}1Rmn)dUn.
Finally, using  (4.10) and (4.17), the elastoplastic consistent stiffness matrix KijALis obtained by the following expression:(4.18)KijAL=Rij-(β1Rim{∂Φ∂Fm}1{∂Φ∂Fn}1Rnj-β2Rim{∂Φ∂Fm}1{∂Φ∂Fn}2Rnj)-(β3Rim{∂Φ∂Fm}2{∂Φ∂Fn}2Rnj-β4Rim{∂Φ∂Fm}2{∂Φ∂Fn}1Rnj).


## 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 setTOL=10-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 loadP 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)(b)(c)(d)(e)The values of the residual force and residual yield norms for a typical load step are summarized in Table1 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 algorithmResidual force norm rnormResidual yieldnorm ΦnormResidual force norm rnormResidual yieldnorm ΦnormStep 47Step 47Step 445Step 4451.47597039E-031.3055981E-038.56514641E-059.37779227E-051.48013748E-064.26512984E-075.9731394E-092.21160446E-093.72157338E-134.59632332E-149.21947265E-172.22044605E-16Table 2
Residual force and residual yield norms for typical load steps during the formation of first and third plastic hinges.
Two-vector return algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 106Step 106Step 106Step 450Step 450Step 4501.65303075E-031.20357837E-051.47408979E-053.87694523E-036.19251038E-044.56960163E-031.89552293E-061.23009056E-075.42558567E-071.60549380E-057.05089222E-075.25607314E-065.81448484E-131.04360964E-147.37188088E-144.57042787E-112.61368704E-127.06168457E-12The values of the energy and residual norms for typical load step are summarized in Table1 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 Table3. 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 algorithmResidual normEnergy normResidual normEnergy normStep 16Step 16Step 250Step 2500.11634E-040.10000E+010.21960E-030.10000E+010.18054E-080.41898E-120.18174E-070.19603E-100.74031E-160.93424E-240.18630E-140.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 Figure4. 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)(b)(c)(d)(e)The values of the residual force and residual yield norms for load steps arbitrarily chosen are summarized in Table4 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 algorithmResidual forcenorm rnormResidual yieldnorm ΦnormResidual forcenorm rnormResidual yieldnorm ΦnormStep 11Step 11Step 60Step 601.20935377E-042.41865226E-078.93804269E-046.25601990E-063.45076307E-116.60402844E-116.22404505E-099.11217879E-098.02626977E-170.05.12166402E-170.0The values of the energy and residual norms, for each iteration, are summarized in Table5. 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 algorithmResidual normEnergy normResidual normEnergy normStep 16Step 16Step 250Step 2500.80063E-020.10000E+010.23954E-010.10000E+010.18017E-030.24336E-060.14595E-010.82672E-030.15248E-080.50571E-130.15491E-150.11855E-080.49717E-130.39726E-230.87519E-130.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 Figure5. 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.The values of the residual force and residual yield norms for arbitrary load steps (12 and 41) are summarized in Tables6 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=10-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 algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 12Step 12Step 12Step 41Step 41Step 418.03978006E-039.81043834E-039.79441075E-031.36989462E-021.79089005E-021.78810582E-026.24153357E-052.38396940E-052.37928885E-051.96679153E-047.88306625E-057.86893201E-058.85460742E-101.48586699E-101.44924517E-109.31316656E-091.60199565E-091.58390967E-098.20520831E-170.03.33066907E-162.40162583E-170.00.0Table 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 algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 12Step 12Step 12Step 41Step 41Step 418.03359176E-039.80846569E-039.78495657E-031.36815406E-021.78952610E-021.78530190E-026.23470796E-052.38312272E-052.37622611E-051.96249346E-047.87644610E-057.85560337E-058.84047353E-101.42486467E-101.41957335E-109.29161239E-091.76712267E-091.64043112E-093.12120416E-181.11022302E-161.11022302E-168.12816025E-174.44089210E-160.0Table 8
Four-legged jacket type-1. Error norms for Newton iterative scheme on global equilibrium.
Single-vector return algorithmResidual normEnergy normResidual normEnergy normStep 12Step 12Step 52Step 520.11997E-030.10000E+010.20903E-030.10000E+010.16535E-090.27673E-120.41175E-070.43179E-100.18777E-150.41131E-250.36921E-150.24472E-22Table 9
Four-legged jacket type-2. Error norms for Newton iterative scheme on global equilibrium.
Single-vector return algorithmResidual normEnergy normResidual normEnergy normStep 25Step 25Step 45Step 450.34837E-030.10000E+010.75999E-030.10000E+010.10245E-050.44342E-080.81462E-060.89862E-080.95123E-110.12589E-150.13896E-110.88152E-16


## 5.1. Right-ِAngle Beam

In this example, a right-angle beam under a concentrated loadP 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)(b)(c)(d)(e)The values of the residual force and residual yield norms for a typical load step are summarized in Table1 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 algorithmResidual force norm rnormResidual yieldnorm ΦnormResidual force norm rnormResidual yieldnorm ΦnormStep 47Step 47Step 445Step 4451.47597039E-031.3055981E-038.56514641E-059.37779227E-051.48013748E-064.26512984E-075.9731394E-092.21160446E-093.72157338E-134.59632332E-149.21947265E-172.22044605E-16Table 2
Residual force and residual yield norms for typical load steps during the formation of first and third plastic hinges.
Two-vector return algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 106Step 106Step 106Step 450Step 450Step 4501.65303075E-031.20357837E-051.47408979E-053.87694523E-036.19251038E-044.56960163E-031.89552293E-061.23009056E-075.42558567E-071.60549380E-057.05089222E-075.25607314E-065.81448484E-131.04360964E-147.37188088E-144.57042787E-112.61368704E-127.06168457E-12The values of the energy and residual norms for typical load step are summarized in Table1 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 Table3. 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 algorithmResidual normEnergy normResidual normEnergy normStep 16Step 16Step 250Step 2500.11634E-040.10000E+010.21960E-030.10000E+010.18054E-080.41898E-120.18174E-070.19603E-100.74031E-160.93424E-240.18630E-140.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 Figure4. 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)(b)(c)(d)(e)The values of the residual force and residual yield norms for load steps arbitrarily chosen are summarized in Table4 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 algorithmResidual forcenorm rnormResidual yieldnorm ΦnormResidual forcenorm rnormResidual yieldnorm ΦnormStep 11Step 11Step 60Step 601.20935377E-042.41865226E-078.93804269E-046.25601990E-063.45076307E-116.60402844E-116.22404505E-099.11217879E-098.02626977E-170.05.12166402E-170.0The values of the energy and residual norms, for each iteration, are summarized in Table5. 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 algorithmResidual normEnergy normResidual normEnergy normStep 16Step 16Step 250Step 2500.80063E-020.10000E+010.23954E-010.10000E+010.18017E-030.24336E-060.14595E-010.82672E-030.15248E-080.50571E-130.15491E-150.11855E-080.49717E-130.39726E-230.87519E-130.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 Figure5. 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.The values of the residual force and residual yield norms for arbitrary load steps (12 and 41) are summarized in Tables6 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=10-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 algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 12Step 12Step 12Step 41Step 41Step 418.03978006E-039.81043834E-039.79441075E-031.36989462E-021.79089005E-021.78810582E-026.24153357E-052.38396940E-052.37928885E-051.96679153E-047.88306625E-057.86893201E-058.85460742E-101.48586699E-101.44924517E-109.31316656E-091.60199565E-091.58390967E-098.20520831E-170.03.33066907E-162.40162583E-170.00.0Table 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 algorithmResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normResidual force normrnormResidual yield normΦ1normResidual yield normΦ2normStep 12Step 12Step 12Step 41Step 41Step 418.03359176E-039.80846569E-039.78495657E-031.36815406E-021.78952610E-021.78530190E-026.23470796E-052.38312272E-052.37622611E-051.96249346E-047.87644610E-057.85560337E-058.84047353E-101.42486467E-101.41957335E-109.29161239E-091.76712267E-091.64043112E-093.12120416E-181.11022302E-161.11022302E-168.12816025E-174.44089210E-160.0Table 8
Four-legged jacket type-1. Error norms for Newton iterative scheme on global equilibrium.
Single-vector return algorithmResidual normEnergy normResidual normEnergy normStep 12Step 12Step 52Step 520.11997E-030.10000E+010.20903E-030.10000E+010.16535E-090.27673E-120.41175E-070.43179E-100.18777E-150.41131E-250.36921E-150.24472E-22Table 9
Four-legged jacket type-2. Error norms for Newton iterative scheme on global equilibrium.
Single-vector return algorithmResidual normEnergy normResidual normEnergy normStep 25Step 25Step 45Step 450.34837E-030.10000E+010.75999E-030.10000E+010.10245E-050.44342E-080.81462E-060.89862E-080.95123E-110.12589E-150.13896E-110.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.



---

*Source: 142743-2010-11-07.xml*