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

Abstract

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

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+𝛼₃||𝐹^𝑦||𝐹_𝑦𝑝^𝛼4+𝛼₅||𝐹^𝑧||𝐹_𝑧𝑝^𝛼6+𝛼₇||𝑀^𝑥||𝑀_𝑥𝑝^𝛼8+𝛼₉||𝑀^𝑦||𝑀_𝑦𝑝^𝛼10+𝛼₁₁||𝑀^𝑧||𝑀_𝑧𝑝^𝛼12+𝛼₁₃||𝐹^𝑥||𝐹_𝑥𝑝^𝛼14||𝑀^𝑦||𝑀_𝑦𝑝^𝛼15+𝛼₁₆||𝐹^𝑥||𝐹_𝑥𝑝^𝛼17||𝑀^𝑧||𝑀_𝑧𝑝^𝛼18−1=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:𝜕Φ𝜕𝐹_𝑥=⎡⎢⎢⎢⎣𝛼₁𝛼₂|||𝐹_{𝛼2𝑥−1}|||𝐹^𝛼2_𝑥𝑝+𝛼₁₃𝛼₁₄|||𝐹_{𝛼14𝑥−1}|||𝐹^𝛼14_𝑥𝑝||𝑀^𝑦||𝑀_𝑦𝑝^𝛼15+𝛼₁₆𝛼₁₇|||𝐹_{𝛼17𝑥−1}|||𝐹^𝛼17_𝑥𝑝||𝑀^𝑧||𝑀_𝑧𝑝^𝛼18⎤⎥⎥⎥⎦sgn𝐹_𝑥𝜕Φ,𝜕𝐹_𝑦=𝛼₃𝛼₄|||𝐹_{𝛼4𝑦−1}|||𝐹^𝛼4_𝑦𝑝sgn𝐹_𝑦;𝜕Φ𝜕𝐹_𝑧=𝛼₅𝛼₆|||𝐹_{𝛼6𝑧−1}|||𝐹^𝛼6_𝑧𝑝𝜕Φsgn𝐹_𝑧,𝜕𝑀_𝑥=𝛼₇𝛼₈|||𝑀_{𝛼8𝑥−1}|||𝑀^𝛼8_𝑥𝑝𝜕Φsgn𝑀_𝑥,𝜕𝑀_𝑦=⎡⎢⎢⎢⎣𝛼₉𝛼₁₀|||𝑀_{𝛼10𝑦−1}|||𝑀^𝛼10_𝑦𝑝+𝛼₁₃𝛼₁₅||𝐹^𝑥||𝐹_𝑥𝑝^𝛼14|||𝑀_{𝛼15𝑦−1}|||𝑀^𝛼15_𝑦𝑝⎤⎥⎥⎥⎦sgn𝑀_𝑦𝜕Φ,𝜕𝑀_𝑧=⎡⎢⎢⎢⎣𝛼₁₁𝛼₁₂|||𝑀_{𝛼12𝑧−1}|||𝑀^𝛼12_𝑧𝑝+𝛼₁₆𝛼₁₈||𝐹^𝑥||𝐹_𝑥𝑝^𝛼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𝑥=𝛼₁𝛼₂𝛼₂−1|||𝐹^{𝛼2𝑥−2}|||𝐹^𝛼2_𝑥𝑝+𝛼₁₃𝛼₁₄𝛼₁₄−1|||𝐹^{𝛼14𝑥−2}|||𝐹^𝛼14_𝑥𝑝||𝑀^𝑦||𝑀_𝑦𝑝^𝛼15+𝛼₁₆𝛼₁₇𝛼₁₇−1|||𝐹^{𝛼17𝑥−2}|||𝐹^𝛼17_𝑥𝑝||𝑀^𝑧||𝑀_𝑧𝑝^𝛼18𝜕²Φ𝜕𝐹_𝑥𝜕𝐹_𝑦=0,𝜕²Φ𝜕𝐹_𝑥𝜕𝐹_𝑧=0,𝜕²Φ𝜕𝐹_𝑥𝜕𝑀_𝑥𝜕=0,₂Φ𝜕𝐹_𝑥𝜕𝑀_𝑦=𝛼₁₃𝛼₁₄𝛼₁₅⎛⎜⎜⎜⎝|||𝐹^{𝛼14𝑥−1}|||𝐹^𝛼14_𝑥𝑝⎞⎟⎟⎟⎠⎛⎜⎜⎜⎝|||𝑀^{𝛼15𝑦−1}|||𝑀^𝛼15_𝑦𝑝⎞⎟⎟⎟⎠sgn𝑀_𝑦𝜕₂sgn𝐹_𝑥,Φ𝜕𝐹_𝑥𝜕𝑀_𝑧=𝛼₁₆𝛼₁₇𝛼₁₈⎛⎜⎜⎜⎝|||𝐹^{𝛼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𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘=₁⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜕²Φ𝜕𝐹_𝑥1𝜕𝐹_𝑥1𝜕²Φ𝜕𝐹_𝑥1𝜕𝐹_𝑦1𝜕²Φ𝜕𝐹_𝑥1𝜕𝐹_𝑧1𝜕²Φ𝜕𝐹_𝑥1𝜕𝑀_𝑥1𝜕²Φ𝜕𝐹_𝑥1𝜕𝑀_𝑦1𝜕²Φ𝜕𝐹_𝑥1𝜕𝑀_𝑧1𝜕000000₂Φ𝜕𝐹_𝑦1𝜕𝐹_𝑥1𝜕²Φ𝜕𝐹_𝑦1𝜕𝐹_𝑦1𝜕²Φ𝜕𝐹_𝑦1𝜕𝐹_𝑧1𝜕²Φ𝜕𝐹_𝑦1𝜕𝑀_𝑥1𝜕²Φ𝜕𝐹_𝑦1𝜕𝑀_𝑦1𝜕²Φ𝜕𝐹_𝑦1𝜕𝑀_𝑧1𝜕000000₂Φ𝜕𝐹_𝑧1𝜕𝐹_𝑥1𝜕²Φ𝜕𝐹_𝑧1𝜕𝐹_𝑦1𝜕²Φ𝜕𝐹_𝑧1𝜕𝐹_𝑧1𝜕²Φ𝜕𝐹_𝑧1𝜕𝑀_𝑥1𝜕²Φ𝜕𝐹_𝑧1𝜕𝑀_𝑦1𝜕²Φ𝜕𝐹_𝑧1𝜕𝑀_𝑧1𝜕000000₂Φ𝜕𝑀_𝑥1𝜕𝐹_𝑥1𝜕²Φ𝜕𝑀_𝑥1𝜕𝐹_𝑦1𝜕²Φ𝜕𝑀_𝑥1𝜕𝐹_𝑧1𝜕²Φ𝜕𝑀_𝑥1𝜕𝑀_𝑥1𝜕²Φ𝜕𝑀_𝑥1𝜕𝑀_𝑦1𝜕²Φ𝜕𝑀_𝑥1𝜕𝑀_𝑧1𝜕000000₂Φ𝜕𝑀_𝑦1𝜕𝐹_𝑥1𝜕²Φ𝜕𝑀_𝑦1𝜕𝐹_𝑦1𝜕²Φ𝜕𝑀_𝑦1𝜕𝐹_𝑧1𝜕²Φ𝜕𝑀_𝑦1𝜕𝑀_𝑥1𝜕²Φ𝜕𝑀_𝑦1𝜕𝑀_𝑦1𝜕²Φ𝜕𝑀_𝑦1𝜕𝑀_𝑧1𝜕000000₂Φ𝜕𝑀_𝑧1𝜕𝐹_𝑥1𝜕²Φ𝜕𝑀_𝑧1𝜕𝐹_𝑦1𝜕²Φ𝜕𝑀_𝑧1𝜕𝐹_𝑧1𝜕²Φ𝜕𝑀_𝑧1𝜕𝑀_𝑥1𝜕²Φ𝜕𝑀_𝑧1𝜕𝑀_𝑦1𝜕²Φ𝜕𝑀_𝑧1𝜕𝑀_𝑧1000000000000000000000000000000000000000000000000000000000000000000000000000000⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦_12×12𝜕,(2.5)₂Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₂⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣00000000000000000000000000000000000=0000000𝜕₂000000000000000000000000000000000000Φ𝜕𝐹_𝑥2𝜕𝐹_𝑥2𝜕²Φ𝜕𝐹_𝑥2𝜕𝐹_𝑦2𝜕²Φ𝜕𝐹_𝑥2𝜕𝐹_𝑧2𝜕²Φ𝜕𝐹_𝑥2𝜕𝑀_𝑥2𝜕²Φ𝜕𝐹_𝑥2𝜕𝑀_𝑦2𝜕²Φ𝜕𝐹_𝑥2𝜕𝑀_𝑧2000000𝜕²Φ𝜕𝐹_𝑦2𝜕𝐹_𝑥12𝜕²Φ𝜕𝐹_𝑦2𝜕𝐹_𝑦2𝜕²Φ𝜕𝐹_𝑦2𝜕𝐹_𝑧2𝜕²Φ𝜕𝐹_𝑦2𝜕𝑀_𝑥2𝜕²Φ𝜕𝐹_𝑦2𝜕𝑀_𝑦2𝜕²Φ𝜕𝐹_𝑦2𝜕𝑀_𝑧2000000𝜕²Φ𝜕𝐹_𝑧2𝜕𝐹_𝑥2𝜕²Φ𝜕𝐹_𝑧2𝜕𝐹_𝑦2𝜕²Φ𝜕𝐹_𝑧2𝜕𝐹_𝑧2𝜕²Φ𝜕𝐹_𝑧2𝜕𝑀_𝑥2𝜕²Φ𝜕𝐹_𝑧2𝜕𝑀_𝑦2𝜕²Φ𝜕𝐹_𝑧2𝜕𝑀_𝑧2000000𝜕²Φ𝜕𝑀_𝑥2𝜕𝐹_𝑥2𝜕²Φ𝜕𝑀_𝑥2𝜕𝐹_𝑦2𝜕²Φ𝜕𝑀_𝑥2𝜕𝐹_𝑧2𝜕²Φ𝜕𝑀_𝑥2𝜕𝑀_𝑥2𝜕²Φ𝜕𝑀_𝑥2𝜕𝑀_𝑦2𝜕²Φ𝜕𝑀_𝑥2𝜕𝑀_𝑧2000000𝜕²Φ𝜕𝑀_𝑦2𝜕𝐹_𝑥2𝜕²Φ𝜕𝑀_𝑦2𝜕𝐹_𝑦2𝜕²Φ𝜕𝑀_𝑦2𝜕𝐹_𝑧2𝜕²Φ𝜕𝑀_𝑦2𝜕𝑀_𝑥2𝜕²Φ𝜕𝑀_𝑦2𝜕𝑀_𝑦2𝜕²Φ𝜕𝑀_𝑦2𝜕𝑀_𝑧2000000𝜕²Φ𝜕𝑀_𝑧2𝜕𝐹_𝑥2𝜕²Φ𝜕𝑀_𝑧2𝜕𝐹_𝑦2𝜕²Φ𝜕𝑀_𝑧2𝜕𝐹_𝑧2𝜕²Φ𝜕𝑀_𝑧2𝜕𝑀_𝑥2𝜕²Φ𝜕𝑀_𝑧2𝜕𝑀_𝑦2𝜕²Φ𝜕𝑀_𝑧2𝜕𝑀_𝑧2⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦_12×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_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁,𝜆₁>0𝐹^trial_𝑖=𝐹_𝑖+𝐾_𝑖𝑗𝑑𝑈_𝑗,(3.1a)𝐾_𝑖𝑗=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐸𝐴𝐿12𝐸0𝐼_𝑧0000−𝐸𝐿𝐴00000𝐿³0006𝐸𝐼^𝑧𝐿²0−12𝐸𝐼^𝑧𝐿³0006𝐸𝐼^𝑧𝐿²12𝐸𝐼_𝑦𝐿³0−6𝐸𝐼^𝑦𝐿²000−12𝐸𝐼^𝑦𝐿³0−6𝐸𝐼^𝑦𝐿²𝐺𝐽_𝑥0𝐿4𝐸0𝐼_𝑦0000−𝐺𝐿𝐽^𝑥00𝐿0006𝐿𝐸²𝐼^𝑦02𝐸𝐼^𝑦4𝐸𝐼_𝑧𝐿0𝐿0−6𝐿𝐸²𝐼^𝑧0002𝐸𝐼^𝑧𝐸𝐴𝐿sym𝐿12𝐸0𝐼_𝑧0000𝐿³000−6𝐸𝐼^𝑧𝐿²12𝐸𝐼_𝑦𝐿³06𝐸𝐼^𝑦𝐿³𝐺𝐽_𝑥0𝐿4𝐸0𝐼_𝑦0𝐿4𝐸0𝐼_𝑧𝐿⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(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; {𝜕Φ/𝜕𝐹_𝑗}₁ is the plastic potential defined in the trial nodal-force; 𝜆₁ is the plastic multiplier and the term 𝜆₁𝐾_𝑖𝑗{𝜕Φ/𝜕𝐹_𝑗}₁ 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).

(a)

(b)

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_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁.(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_𝑖+𝑑𝐹_𝑖+𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁𝑑𝐹_𝑘,(3.3) where 𝑑𝐹_𝑖 is an infinitesimal increment of the generalized force vector 𝐹_𝑖; 𝑑𝜆₁ is the variation of the plastic multiplier 𝜆₁ and [𝜕²Φ/𝜕𝐹_𝑗𝜕𝐹_𝑘]𝑑𝐹_𝑘 is the change in the potential plastic vector {𝜕Φ/𝜕𝐹_𝑗}₁. The goal is to achieve 𝑟^new_𝑖=0 in (3.3). For that reason 0=𝑟^old_𝑖+𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝛿_𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁𝑑𝐹^𝑘.(3.4) To determine an expression to 𝑑𝐹_𝑖 (the correcting force vector) in (3.4), define the square matrix 𝑄_𝑖𝑘 as 𝑄_𝑖𝑘=𝛿_𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁.(3.5) Now, considering (3.4) and (3.5) and after some algebraic manipulations to solve (3.4) for 𝑑𝐹_𝑖, it follows that𝑑𝐹_𝑖=−𝑄^−1𝑖ℓ𝑟^oldℓ+𝑑𝜆₁𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁.(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 Φ^new₁=Φ^old₁+𝜕Φ𝜕𝐹_𝑖₁𝑑𝐹_𝑖.(3.7) Since the expression for 𝑑𝐹_𝑖 is available in (3.6) and for Φ^new₁=0, the variation of the plastic multiplier 𝑑𝜆₁ is readily found𝑑𝜆₁=Φ^old1−𝜕Φ/𝜕𝐹_𝑖₁𝑄^−1𝑖ℓ𝑟^old_ℓ𝜕Φ/𝜕𝐹_𝑖₁𝑄^−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=10⁻¹⁰.

3.2. Two-Vectors Return Algorithm

When two plastic hinges occur at the ends of the beam element, the Two-Vector Return Algorithm (2VRA) is applied. This algorithm is used so that the nodal-force vector at both ends of the beam element can be brought back to the yield surface. The geometric interpretation of this algorithm is illustrated in Figure 2. The Backward Euler force obtained with the 2VRA may be defined by the following expression:𝐹_𝑖=𝐹^trial_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁−𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂,with𝜆₁>0,𝜆₂>0.(3.9) To derive a scheme for 2VRA, the residual force vector 𝑟_𝑖 is defined as𝑟_𝑖=𝐹_𝑖−𝐹_𝑖=𝐹_𝑖−𝐹^trial_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁−𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂.(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_𝑖+𝑑𝐹_𝑖+𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁𝑑𝐹_𝑘+𝑑𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂+𝜆₂𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₂𝑑𝐹(_𝑘3..11) Making 𝑟^new_𝑖=0 in (3.11), it follows that0=𝑟^old_𝑖+𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂+𝛿_𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁+𝜆₂𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₂𝑑𝐹^𝑘.(3.12) To determine an expression to 𝑑𝐹_𝑖 (the correcting force vector) in (3.12), define the square matrix 𝑄_𝑖𝑘 as 𝑄_𝑖𝑘=𝛿_𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁+𝜆₂𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₂.(3.13) Now considering (3.13) and (3.12) and after some algebraic manipulations to solve (3.12) for 𝑑𝐹_𝑖, it follows that 𝑑𝐹_𝑖=−𝑄^−1𝑖ℓ𝑟^oldℓ+𝑑𝜆₁𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₂.(3.14) First-order Taylor’s series expansions of the yield function at node 1, Φ₁, and at node 2, Φ₂, around the final force vector 𝐹_𝑖 are now developed. These series expansions generate linear approximations of the new values of the yield functions Φ^new₁ and Φ^new₂ as Φ^new₁=Φ^old₁+𝜕Φ𝜕𝐹_𝑖₁𝑑𝐹_𝑖,Φ^new₂=Φ^old₂+𝜕Φ𝜕𝐹_𝑖₂𝑑𝐹_𝑖.(3.15) Since the expression for 𝑑𝐹_𝑖 is available in (3.14) and the goal is to achieve both Φ^new₁=0 and Φ^new₂=0, then Φ^old₁−𝜕Φ𝜕𝐹_𝑖₁𝑄^−1𝑖ℓ𝑟^old_ℓ=𝑑𝜆₁𝜕𝜕𝐹Φ_𝑖₁𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝜕𝜕𝐹Φ_𝑖₁𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₂Φ_old2,−𝜕Φ𝜕𝐹_𝑖₂𝑄^−1𝑖ℓ𝑟^old_ℓ=𝑑𝜆₁𝜕𝜕𝐹Φ_𝑖₂𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝜕𝜕𝐹Φ_𝑖₂𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₂.(3.16) The two previous expressions form a system of equations in 𝑑𝜆₁ and 𝑑𝜆₂. By introducing auxiliary variables𝑎¹¹𝑎₁₂𝑎₂₁𝑎₂₂𝑑𝑑𝜆𝜆¹₂=𝑏𝑏¹₂(3.17) and now using Cramer's rule for the previous system of equations, the variations of the plastic multiplier 𝑑𝜆₁ and 𝑑𝜆₂ are readily found by the expressions𝑑𝜆₁=𝑏¹𝑎₂₂−𝑏₂𝑎₁₂𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁,𝑑𝜆₂=𝑏²𝑎₁₁−𝑏₁𝑎₂₁𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁,(3.18) where the auxiliary variables are defined as 𝑏₁=Φ^old₁−𝜕Φ𝜕𝐹_𝑖₁𝑏₂𝑄^−1𝑖ℓ𝑟^old_ℓ,=Φ^old₂−𝜕Φ𝜕𝐹_𝑖₂𝑎₁₁𝑄^−1𝑖ℓ𝑟^old_ℓ,=𝜕Φ𝜕𝐹_𝑖₁𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁𝑎₁₂,=𝜕Φ𝜕𝐹_𝑖₁𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₂𝑎₂₁,=𝜕Φ𝜕𝐹_𝑖₂𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₁𝑎₂₂,=𝜕Φ𝜕𝐹_𝑖₂𝑄^−1𝑖ℓ𝐾_ℓ𝑗𝜕𝜕𝐹Φ_𝑗₂.(3.19)

(a)

(b)

Figure 2 

Two-Vectors return algorithm—two plastic hinges.

This iterative procedure is continued until the yield conditions Φ₁=0 and Φ₂=0 are satisfied at the final force vector state, that is, when 𝑟^norm=|𝑟_𝑖|/|𝐹^trial_𝑖|<TOL, Φ^1norm=|Φ₁|<TOL and Φ^2norm=|Φ₂|<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_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁,with𝜆₁>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; 𝑑𝐹_𝑖=𝐾_𝑖𝑗𝑑𝑈_𝑗−𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁−𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁𝑑𝐹_𝑘.(4.2) After some algebraic treatment, 𝛿^𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁𝑑𝐹^𝑘=𝐾_𝑖𝑗𝑑𝑈_𝑗−𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁.(4.3) Using matrix 𝑄_𝑖𝑘 as defined by (3.5) and denoting the reduced stiffness matrix as 𝑅_𝑖𝑗=𝑄^−1𝑖ℓ𝐾_ℓ𝑗, (4.3) changes into the following equation:𝑑𝐹_𝑖=𝑅_𝑖𝑗𝑑𝑈^𝑗−𝑑𝜆₁𝜕𝜕𝐹Φ_𝑗₁.(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 𝜕𝜕𝐹Φ_𝑖₁𝑑𝐹_𝑖=0⟹𝜕Φ𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝑑𝑈^𝑗−𝑑𝜆₁𝜕𝜕𝐹Φ_𝑗₁=0,(4.5) and, therefore, the change in the plastic multiplier 𝑑𝜆₁may be expressed as 𝑑𝜆₁=𝜕Φ/𝜕𝐹^𝑖₁𝑅_𝑖𝑗𝑑𝑈_𝑗𝜕Φ/𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝜕Φ/𝜕𝐹_𝑗₁.(4.6) Finally, using (4.4) and (4.6), the elastoplastic consistent stiffness matrix 𝐾^{𝐴𝐿𝑖𝑗}may be determined by 𝐾^{𝐴𝐿𝑖𝑗}=𝑅_𝑖𝑗−𝑅^𝑖𝑚𝜕Φ/(𝜕𝐹_𝑚)₁𝜕Φ/(𝜕𝐹_𝑛)₁𝑅_𝑛𝑗𝜕Φ/(𝜕𝐹_𝑚)₁𝑅_𝑚𝑛𝜕Φ/(𝜕𝐹_𝑛)₁.(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_𝑖−𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁−𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂,with𝜆₁>0,𝜆₂>0.(4.8) Differentiating (4.8) and performing algebraic manipulations, it follows that𝛿^𝑖𝑘+𝜆₁𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₁+𝜆₂𝐾_𝑖𝑗𝜕²Φ𝜕𝐹_𝑗𝜕𝐹_𝑘₂𝑑𝐹^𝑘=𝐾_𝑖𝑗𝑑𝑈_𝑗−𝑑𝜆₁𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁−𝑑𝜆₂𝐾_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂.(4.9) Using matrix 𝑄_𝑖𝑘 as defined by (3.13) and denoting the reduced stiffness matrix as 𝑅_𝑖𝑗=𝑄^−1𝑖ℓ𝐾_ℓ𝑗, we can rewrite the above equation as𝑑𝐹_𝑖=𝑅_𝑖𝑗𝑑𝑈^𝑗−𝑑𝜆₁𝜕𝜕𝐹Φ_𝑗₁−𝑑𝜆₂𝜕𝜕𝐹Φ_𝑗₂.(4.10) Again, it is assumed that the full consistency conditions should hold at the final force point, that is, Φ₁=0 and Φ₂=0. Differentiating these two conditions and introducing the expression for 𝑑𝐹_𝑖 taken from (4.10), one can get𝜕𝜕𝐹Φ_𝑖₁𝑑𝐹_𝑖=0⟹𝜕Φ𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝑑𝑈^𝑗−𝑑𝜆₁𝜕𝜕𝐹Φ_𝑗₁−𝑑𝜆₂𝜕Φ𝜕𝜕𝐹Φ_𝑗₂=0,𝜕𝐹_𝑖₂𝑑𝐹_𝑖=0⟹𝜕Φ𝜕𝐹_𝑖₂𝑅_𝑖𝑗𝑑𝑈^𝑗−𝑑𝜆₁𝜕𝜕𝐹Φ_𝑗₁−𝑑𝜆₂𝜕𝜕𝐹Φ_𝑗₂=0.(4.11) Alternatively, by collecting terms algebraically, (4.11) are transformed in 𝜕𝜕𝐹Φ_𝑖₁𝑅_𝑖𝑗𝑑𝑈_𝑗=𝑑𝜆₁𝜕𝜕𝐹Φ_𝑖₁𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝜕𝜕𝐹Φ_𝑖₁𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂𝜕Φ,𝜕𝐹_𝑖₂𝑅_𝑖𝑗𝑑𝑈_𝑗=𝑑𝜆₁𝜕𝜕𝐹Φ_𝑖₂𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁+𝑑𝜆₂𝜕𝜕𝐹Φ_𝑖₂𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂.(4.12) The previous two expressions represent a system of equations with 𝑑𝜆₁ and 𝑑𝜆₂ as unknowns. Introducing auxiliary variables 𝑐₁, 𝑐₂, 𝑏₁₁, 𝑏₁₂, 𝑏₂₁ and 𝑏₂₂ as𝑐₁=𝜕Φ𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝑑𝑈_𝑗,𝑐₂=𝜕Φ𝜕𝐹_𝑖₂𝑅_𝑖𝑗𝑑𝑈_𝑗,𝑏₁₁=𝜕Φ𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁𝑏₁₂,=𝜕Φ𝜕𝐹_𝑖₁𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂,𝑏₂₁=𝜕Φ𝜕𝐹_𝑖₂𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₁,𝑏₂₂=𝜕Φ𝜕𝐹_𝑖₂𝑅_𝑖𝑗𝜕𝜕𝐹Φ_𝑗₂.(4.13) With the help of Cramer’s rule, the system in (4.12) is solved for 𝑑𝜆₁ and 𝑑𝜆₂𝑑𝜆₁=𝑐¹𝑏₂₂−𝑐₂𝑏₁₂𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁,𝑑𝜆₂=𝑐²𝑏₁₁−𝑐₁𝑏₂₁𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁(4.14) or alternatively 𝑑𝜆₁=𝑏²²𝜕Φ/𝜕𝐹_𝑚₁−𝑏₁₂𝜕Φ/𝜕𝐹_𝑚₂𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁𝑅^𝑚𝑛𝑑𝜆₂𝑑𝑈_𝑛,=𝑏¹¹𝜕Φ/𝜕𝐹_𝑚₂−𝑏₂₁𝜕Φ/𝜕𝐹_𝑚₁𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁𝑅^𝑚𝑛𝑑𝑈_𝑛.(4.15) If any plastic multiplier is negative, that is, 𝑑𝜆₁<0, or 𝑑𝜆₂<0, then this plastic multiplier is set to zero (i.e., 𝑑𝜆₁=0 or 𝑑𝜆₂=0) which corresponds to an initial plastic hinge. Again, introducing new auxiliary variables𝛽₁=𝑏²²𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁,𝛽₂=𝑏¹²𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁,𝛽₃=𝑏¹¹𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁,𝛽₄=𝑏²¹𝑏₁₁𝑏₂₂−𝑏₁₂𝑏₂₁,(4.16) Equation (4.15) can be rewritten as 𝑑𝜆₁=𝛽₁𝜕𝜕𝐹Φ_𝑚₁𝑅_𝑚𝑛−𝛽₂𝜕𝜕𝐹Φ_𝑚₂𝑅_𝑚𝑛𝑑𝑈^𝑛𝑑𝜆₂,=𝛽₃𝜕𝜕𝐹Φ_𝑚₂𝑅_𝑚𝑛−𝛽₄𝜕𝜕𝐹Φ_𝑚₁𝑅_𝑚𝑛𝑑𝑈^𝑛.(4.17) Finally, using (4.10) and (4.17), the elastoplastic consistent stiffness matrix 𝐾^{𝐴𝐿𝑖𝑗}is obtained by the following expression:𝐾^{𝐴𝐿𝑖𝑗}=𝑅_𝑖𝑗−𝛽₁𝑅_𝑖𝑚𝜕𝜕𝐹Φ_𝑚₁𝜕𝜕𝐹Φ_𝑛₁𝑅_𝑛𝑗−𝛽₂𝑅_𝑖𝑚𝜕𝜕𝐹Φ_𝑚₁𝜕𝜕𝐹Φ_𝑛₂𝑅_𝑛𝑗−𝛽₃𝑅_𝑖𝑚𝜕𝜕𝐹Φ_𝑚₂𝜕𝜕𝐹Φ_𝑛₂𝑅_𝑛𝑗−𝛽₄𝑅_𝑖𝑚𝜕𝜕𝐹Φ_𝑚₂𝜕𝜕𝐹Φ_𝑛₁𝑅_𝑛𝑗.(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=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 “𝑛” 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].

(a)

(b)

(c)

(d)

(e)

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.

Table 1 

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

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.

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

(a)

(b)

(c)

(d)

(e)

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.

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.

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

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

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.

Table 8 

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

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.