An accelerated incremental algorithm to trace the nonlinear equilibrium path of structures

This paper deals with the convergence acceleration of iterative nonlinear methods. An effective iterative algorithm, named the three–point method, is applied to nonlinear analysis of structures. In terms of computational cost, each iteration of the three–point method requires three evaluations of the function. In this study the effective functions have been proposed to accelerate the convergence process. The proposed method has a convergence order of eight, and it is important to note that its implementation does not require the computation of higher order derivatives compared to most other methods of the same order. To trace the equilibrium path beyond the limit point, a normal flow algorithm is implemented into a computer program. The three–point method is applied as an inner step in the normal flow algorithm. The procedure can be used for structures with complex behavior, including: unloading, snap–through, elastic post–buckling and inelastic post–buckling analyses. Several numerical examples are given to illustrate the efficiency and performance of the new method. Results show that the new method is comparable with the well–known existing methods and gives better results in convergence speed.


INTRODUCTION
A lot of research has focused on presenting formulations for the nonlinear analysis of structures.Various analysis methods have been proposed to estimate the real behavior of structures.Some of the more interesting structural forms include lattice domes or shallow reticulated caps that span long distances.These structures function as space frames, and are often used in place of continuous shell-type structures.The most common mode of failure for lattice domes or reticulated caps is instability, and the complex geometry used in their design usually prevents a closed form solution for the critical load.
Large deformation analysis requires the solution of a nonlinear system of equations.Nonlinear systems of equations are most commonly solved using iterative incremental techniques where small incremental changes in displacement are found by imposing small incremental changes in load on the structure.The resulting solutions are used to plot the equilibrium path for the structure.An excellent review of nonlinear analysis solution techniques and their implementation is given by Crisfield and Shi [10,11].The most common technique for solving nonlinear finite element equations is the Newton-Raphson (N-R) method.The Newton-Raphson method is famous for its rapid convergence but is known to fail at points (limit points) on the equilibrium path where the Jacobian (tangent stiffness) is singular or nearly singular.Bathe and Cimento [4] highlight some of the problems with the Newton-Raphson method and present various forms of the method that involve accelerations or line searches to maintain convergence during the solution process.
Arc-length methods have been used to overcome the problem of limit points.Arc-length methods are similar to the Newton-Raphson method except that the applied load increment becomes an additional unknown.A comparative study of arc-length methods is presented by Clarke and Hancock [9].The original idea behind the arc-length method was introduced by Riks [20] and Wempner [26].Eriksson [13] described a general setting for the evaluation of onedimensional generalized solution paths to augmented structural equilibrium problems.Saffari et al. [21] introduced a new algorithm for passing the equilibrium path, known as modified normal flow algorithm, for geometrically nonlinear analysis of space trusses.This algorithm can reduce both number of iterations and computing time.Also, other methods have been developed to surpass limit points: (a) the generalized displacement control method (GDC) which is one of the most reliable numerical procedures for nonlinear static analysis [6].Cardoso and Fonseca [6] shown that the GDC method can be seen as an orthogonal arc-length method with additional features.This method is originally proposed by Yang and Shieh [27]; (b) the work control method which was proposed by Bathe and Dvorkin [5] to enforce a constant value of work done in each iteration.This method avoids the limitations of the load control and displacement control methods in tracing the equilibrium path; (c) Chan [7] proposed the minimum residual displacement method (MRD) to remove the residual displacement or the unbalanced force in each iteration.A new approach for nonlinear analysis of structures, which accelerates the convergence rate, has been presented by Saffari and Mansouri [22].In [22], Saffari and Mansouri employed a mathematical method, namely two-point, to achieve the convergence state; however, in [22] the equilibrium path is traced until the limit point.A concept to accelerate the trend in nonlinear analysis and aimed to gain the ability for analysis of complex structures has been introduced by Saffari et al. [23].
Among numerical methods for solving nonlinear equations, the Newton-Raphson method dominates by reason of its quadratic convergence characteristics; however, in order to obtain such a convergence speed, a large amount of computational effort is needed.This effort is placed at each iteration in the construction of the tangent stiffness matrix and in the solving of the linear system with the new matrix.To reduce the global computational cost, a modified Newton-Raphson process with the same Jacobian at all iterations can be applied, but the quadratic convergence characteristic is then lost.After a few Newton-Raphson iterations, the Jacobian is recalculated only at a constant interval.These two approaches save the Jacobian calculation between the intervals.But at the same time the total number of Newton-Raphson iterations increases.In some cases, if less computational effort is required at the iterations which use the same Jacobian matrix than in the proceeding iterations, although the number of Newton-Raphson iterations will increase, the overall computational time can be reduced.
Higher order methods which require the second or higher order derivatives can be time consuming.One such method is the Newton method; however the Newton method can suffer from numerical instabilities.Thus, it is important to study the higher order variants of Newton's method, which require only function and first derivative calculations.In recent years, much attention has been given to develop several iterative methods for solving nonlinear equations.For instance, Two-point methods have been suggested by combining the well-known Newton's method with other one-point methods [18].Though the authors claimed the methods to be original, Petkovic and Petkovic [17] mentioned that these methods are found in Traub's book and are rediscovered in another way.Babajee and Dauhoo [2] analyzed the properties of these variants.Some methods were also derived and rediscovered from the Adomian Decomposition Method [8].
This paper presents an efficient method which accelerates the nonlinear analysis process.In this method, some path-following algorithms have been used to pass the limit points.A computer program is presented based on the concept developed herein which incorporates the three-point method.The number of function evaluations in this method is three, and new and effective functions have been introduced to evaluate the predictor phase.Numerical examples are presented and solved using the developed program.Results show that the proposed method is efficient in different kinds of path-following algorithms.

NONLINEAR ANALYSIS OF STRUCTURES
In the following sections, large deflection inelastic analysis of structures including both geometric and material nonlinearities is briefly discussed.It is followed by a detailed description of the concept developed in this study.

Geometrically nonlinear analysis
Consider an arbitrary framed structure loaded at the nodes only, and let δ denote the corresponding deformed configuration.The equations of equilibrium of the system can be written as: in which {f (δ)} is the resultant of the nodal internal forces and citation{P } represents the external nodal forces.The member force deformation relationships denote that {f } is a highly nonlinear function of {δ}.
The load-deflection relations show that it is almost impossible to explicitly solve these equations.For computational purposes, it is useful to apply the differential form of the equation: {∆δ} in which {∆δ} stands for incremental displacements,{∆P } represents load increments, and ∂δ j ] constitutes the system tangent stiffness matrix which is explained in the following sections.

Truss element
The accuracy of the structure inelastic response depends on the accuracy of the member's load-displacement relationship used in the analysis.A number of models have been introduced in the literature to predict the nonlinear behavior of space trusses.In this study, a stressstrain relationship proposed by Hill et al. [15] is adopted to predict the inelastic post-buckling behavior of the trusses.
The curve can be expressed by the following relationship: -For tensile members: where F y denotes yield stress and u y = F y .L/ E.
-For compressive members: The figure of force-strain curve (Q-u/L) assumed applicable for steel materials both in tension and compression has been shown in [23].
Here Q cr = π 2 EI/L 2 (I =weak axis moment of inertia) and Q l is the asymptotic lower stress limit defined as Q l =r.Q cr .The corresponding critical buckling displacement is u cr = Q cr .L /(A.E ) while u' is defined as u' = u-u cr .Parameters X 1 and X 2 are constants depending on the slenderness ratio of the compressive members.
When a member is in compression and u ≥ u cr , the tangent modulus, E t , has to be used in place of E. The tangent modulus is obtained by: It can be seen that if loading path reaches point A, the member behavior follows relations in Eq. (3).

Frame element
A perfectly-plastic material assigned to plastic hinge zones is used in this study to consider material non-linear effects.In an elastic perfectly-plastic material, the effects of strain hardening are neglected.It further implies that once the yield moment M p is reached, the material yields and cannot withstand further stress.The material constitutive behavior is shown in Fig. 1.

Figure 1 Material behavior
The yield moment is often defined through a yield criterion.A variety of yield criteria have been introduced in structural engineering.For the steel elements in this paper, the AISC-LRFD [1] criterion considering bending moment and axial force interaction is determined as follows: in which, M c = ϕ b ZF y , M pc is the reduced plastic moment capacity in the presence of axial force (Q c = ϕ c Q n );ϕ b andϕ c are the bending and axial resistance factors respectively; F y is the yield stress, and Z is the plastic section modulus.
In inelastic analysis, the tangent stiffness matrix introduced in section 2.2 has to be modified.A detailed discussion on the inelastic analysis can be found in [16].

THREE-POINT METHOD
For nonlinear problems, solving the nonlinear part of the system often requires the most computational effort.Since the Newton-Raphson method has quadratic convergence characteristics [22], it is often used for solving nonlinear equations; however, the Newton method can suffer from numerical instabilities.Because of this, it is important to study alternative methods, such as multi-point Newton-Like methods [3].
Multipoint methods without memory are methods that use new information at a number of points.It is an efficient way of generating higher order methods free from second and higher order derivatives.One classical problem in numerical analysis is the solution of nonlinear equations where f (x ) = 0. Efficiently finding zeros in a single variable nonlinear equation is an interesting and very old problem in numerical analysis with many applications to structural engineering.
The three-point iteration scheme for solving nonlinear equations is expressed as follows: It can be showed that the efficiency index of this approach is √ 2 = 1.4142 which is equal to that of Newton's method.In fact, Newton's method is applied three times.Recently, Dzunic et al. [12] have modified the above scheme in order to improve the computational efficiency of the iterative method.The derivatives f ′ (y)andf ′ (z) are approximated using following relations: in which It should be noted that the quantities s and t do not require new information since they are expressed by the already calculated quantities.Now the three-point scheme (7) has the form: where p and q are arbitrary real functions with Taylor's series of the form: It can be proved that the family of three-point methods (7) has the order eight, if p and q conditions (11,12) [12].Since the total number of function evaluations per iteration of the method (3) is four, the efficiency index of proposed method is 4  √ 8 = 1.6818 which is better than that of classic Newton's method.

Choice of functions p and q
The functions p and q have been searched in a general rational form: which satisfies conditions given in Eqs.(11,12) and its parameters (β, γ 1 , γ 2 , γ 3 ) have been selected using trial and error approach.
In this way many different functions can be generated.In this paper functions p 1 thru p 3 and q 1 thru q 4 have been chosen from reference [12].They are mentioned in the following: ) Another function for p and q is introduced in the proposed method, respectively: Among these choices, the results show that choice of p 4 and q 5 has more efficient than the other forms obtained by trial.

THREE-POINT SCHEME IN STRUCTURAL ANALYSIS
The three-point technique can be applied to analysis of structures with complex behaviors, including unloading, snap-through, elastic post-buckling and inelastic post-buckling analyses.In this paper a modified normal flow algorithm [21] is used to trace the equilibrium paths.The theory of this method is described in the following sections.

Modified normal flow algorithm
If i is the step number, j is the number of the modifying iteration, and the total load on the structure is {P } j i , or its equivalent, the product of a total ratio λ j toti and a given reference external load {P ref } is applied through a series of load increments.Mathematically, this is written as: A detailed discussion and the process of the normal flow algorithm have been represented in reference [21].
A nonlinear system of equilibrium equations can generally be formulated as follows: where J(λ j−1 i , δ j−1 i ) is the Jacobian matrix of order N × (N + 1) and S j i is the Newton-Raphson step size which is found through the following relationship: In the normal flow algorithm, the Newton-Raphson step size is the minimum solution of the system of Eq. ( 25) which is unique among an infinite number of solutions obtainable [21].A particular solution {V } should be obtained from: in which: where is the vector of the unbalanced forces, {δ I } j i is vector of tangential displacement in the converged point; and {∆δ R } j i is the vector of unbalanced displacements such that, in which {F } j i is the vector of resultant internal forces at the nodes.The vector of unbalanced force is computed through the following system of equations: Latin American Journal of Solids and Structures 9(2012) 425 -442 Using the following equation the minimum solution of the norm is calculated: in which {δ I } j i is the tangential displacement.A direct method of updating is adopted such that, the load increment is related to the number of iterations.Also sign of determinant of the tangential stiffness matrix of the previous step can be computed through the following relationship: where exponent γis a certain number, J D is the number of iterations assumed at the beginning of the calculations and J M is the number of iterations performed in the previous step.

Algorithm of nonlinear analysis based on three-point method
Here, a general flow chart for the iterative process is presented as follows: Step 1) Initialization of variables and parameters.
Step 2) Input first incremental load, λ 1 1 , J D , J M , boundary conditions, connectivity, material properties and structural geometry.
Step 3) For the first iteration at each increment step i : 3.1 Form the system tangent stiffness matrix [τ ] Step 4) Compute the element force vector.Update the tangent modulus accounting for inelastic effects.Update global stiffness matrix Step 5) Solve for tangential displacement:[τ ] Solve for {∆δ R } j i from Eq. (30).
Step 6) Calculate {d 1 } Step 7) Update global nodal displacements and determine the unbalanced load: Step 8) Compute s: and construct p(s).
Step 9) Determine the new node locations: Step 10) Update global nodal displacements and determine the unbalanced load: Latin American Journal of Solids and Structures 9(2012) 425 -442 Step 11) Compute t: and construct q(s,t).
Step 12) Determine the new node locations: Step 13) Determine the load increment parameter using Eq. ( 28) and update it.
Step 14) Check the convergence criteria by the following relation: where ε is a predefined error in the course of calculations.

NUMERICAL RESULTS
A computer program is developed in this study which incorporates the three-point algorithm.CPU times taken by calculation process can be measured conveniently when the corresponding command prompt on the screen is set.All examples are solved using a 32-bit Pentium 2.00 GHz processor (Dual core).To examine the performance of the method proposed here and to evaluate the results obtained, some examples are presented in this section.Newton-Raphson method and other sub-stepping algorithms including three-point method are applied to three cases of elastic, elastic post-buckling (EPB), and inelastic post-buckling (IPB) analyses of typical and well-known structures.A predefined tolerance of ε=10 −5 is also assigned to all cases defined in the analysis.Parameters r, X 1 , X 2 are assigned fixed values as r = 0.4, X 1 = 50 and X 2 = 100 [25] and remain unchanged throughout the analysis.Newton-Raphson method or three point procedure are used in inner steps of the path-following algorithms.
Comparison between these two approached are listed in the following Tables.

Example 1
The geometric dimensions of the geodesic dome truss shown in Fig. 2 are adopted from Ramesh and Krishnamoorthy [19].This truss has a total of 61 nodes and 156 members with identical cross-section (A = 6.5 cm 2 and I = 1 cm 4 ) located at the outer supports are restrained as pin supports.A concentrated vertical load of P = 8 kN is applied to the center of dome as shown in Fig. 2. Parameters of Eq. ( 32) are selected such that∆λ 1 1 = 0.01, λ max = 0.5, γ = 0.1, J D = 10, J max = 100.The elevation of the truss is defined by the following equation: The material properties are E = 6895 kN/cm 2 and F y = 400 kN/cm 2 .The loaddisplacement curves obtained from four cases of analysis are shown in Fig. 3.
Table 1 presents a summary of the CPU time taken for application of classic Newton-Raphson method and the three-point method in three cases of analysis for various pathfollowing methods.Number of iterations is shown in Table 2. Results show that for all cases analyzed, threepoint approach has better performance as compared to the Newton-Raphson method.

Example 2
A circular dome truss taken from reference [25] is shown in Fig. 4.This structure has 168 elements and 73 nodes for a total of 147 degrees-of-freedom.The structures is subjected to a vertical load P = 500 kN at the apex.The out-of-plane motion has been constrained with pin supports appended to each end of the truss.
The cross-sectional properties for area (A) and moment of inertia (I ) are 50.431cm 2 and 52.94 cm 4 respectively for all members.The elastic modulus (E ) for all members is 2.04×10 4 kN/cm 2 and F y = 25 kN/cm 2 .
For Eq. (32), parameters∆λ 1 1 = 0.01,λ max = 2,J D = 5, J max = 100,γ = 0.1 are selected.Fig. 5 illustrates the numerical responses obtained from the proposed formulation for four cases of analyses in a manner similar to Example 1.The comparisons between the results of applying the two methods (Newton-Raphson and three-point) are listed in Table 3 and Table 4. Normal flow GDC Arc-length MRD algorithm [21] method [6] method [20] method [ Fig. 5 shows the significant shortcoming of Newton-Raphson method in passing limit point.
From Table 3, the three-point algorithm requires the least computational time for convergence compared to Newton-Raphson method.

Example 3
Schewdeler's truss [14], shown in Fig. 6 having 264 elements and 97 nodes with pin supports at the outer nodes provides a good opportunity to evaluate the efficiency of the method discussed herein.The axial stiffness (EA) for all members is taken as 640 ×10 3 kN, F y = 25 kN/cm 2 , and I = 30.04cm 4 .The external loading was a typical equipment loading of the magnitude P = 50 kN at the crown node.For Eq. (32) parameters∆λ 1 1 = 0.01, λ max = 1, γ = 0.1,J D = 2, J max = 100 are used.Fig. 7 shows variation of vertical displacement at central node with the load P obtained in three cases of analysis.When each of these four cases is analyzed using the four different iteration algorithms, a total of 16 sets of results are obtained.For the purpose of this study, it is interesting to investigate the summary of the results and compare the computational time taken by all of methods.Table 5 provides a summary of all results for convenience.
It is important to note that all computations were carried out with the smallest step possible and with high precision.From Table 6, the minimum number of iterations is associated with the three-point scheme.

Example 4
The last example shown in Fig. 8 has been analyzed by Spiliopoulos and Patsios [24].For all members E = 220000 kN/m 2 .Section properties for columns are I = 0.000313 m 4 , A = 0.022 m 2 .Beam sections (including the inclined members) given as: I = 0.0000249 m 4 , A = 0.00543 m 2 .The load increment, ∆P, is 0.795 kN and P = 31.80kN.The frame was analyzed using the three-point method developed in this paper.Fig. 9 shows the load-displacement curve obtained by applying the method developed in the present study and Newton-Raphson method, respectively.
The results show that the number of iterations and the computing time used in deploying the classic Newton-Raphson approach is more than that used by applying the developed method.The results are presented in Table 7. Table 8 shows that three-point scheme also can reduce the computational effort in nonlinear analysis of frames.

CONCLUSIONS
This paper provides a simple and accurate high order approach which can be used in complex structural analysis.A mathematical method which was named three-point method can reduce the global computational time of nonlinear analysis and number of iterations by near 40-50% when compared with other more commonly used methods.In the three-point method using two functions, the acceleration of convergence was increased.Several functions and a new function were proposed and results were compared with those obtained by deploying the classic Newton-Raphson method.A computer program was developed to numerically solve a system of nonlinear equations in an incremental form.The numerical examples represented a good compromise between reducing the computing time, reducing the number of the iterations and obtaining results sufficient accuracy.

LatinFigure 2
Figure 2 Geodesic dome truss, dimensions are given in cm

Figure 3
Figure 3 Load-displacement curves of geodesic dome truss at apex

Figure 4 Figure 5
Figure 4 Circular dome truss solved in example 2 (dimensions are given in cm)

Figure 6
Figure 6 Schewdeler's dome truss, dimensions are given in cm

Figure 7
Figure 7 Central node vertical load-displacement

Figure 8 Figure 9 N
Figure 8 Five-storey frame

Table 1
Comparison of CPU time (sec) for example 1

Table 2
Comparison of number of iterations for example 1

Table 3
Comparison of CPU time (sec) for example 2

Table 4
Comparison of number of iterations for example 2

Table 5
Comparison of CPU time (sec) for example 3

Table 7
Comparison of CPU time (sec) for example 4