Trusses Nonlinear Problems Solution with Numerical Methods of Cubic Convergence Order

ABSTRACT

A large part of the numerical procedures for obtaining the equilibrium path or load-displacement curve of structural problems with nonlinear behavior is based on the Newton-Raphson iterative scheme, to which is coupled the path-following methods. This paper presents new algorithms based on Potra-Pták, Chebyshev and super-Halley methods combined with the Linear Arc-Length path-following method. The main motivation for using these methods is the cubic order convergence. To elucidate the potential of our approach, we present an analysis of space and plane trusses problems with geometric nonlinearity found in the literature. In this direction, we will make use of the Positional Finite Element Method, which considers the nodal coordinates as variables of the nonlinear system instead of displacements. The numerical results of the simulations show the capacity of the computational algorithm developed to obtain the equilibrium path with force and displacement limits points. The implemented iterative methods exhibit better efficiency as the number of time steps and necessary accumulated iterations until convergence and processing time, in comparison with classic methods of Newton-Raphson and Modified Newton-Raphson.

Keywords:
Arc-Length; Positional Finite Element; Chebyshev; Potra-Pt´ak; Geometric Nonlinearity

RESUMO

Grande parte dos procedimentos numéricos para obter a trajetória de equilíbrio ou curva força-deslocamento de problemas estruturais com comportamento não-linear é baseado no esquema iterativo de Newton-Raphson ao qual são acoplados métodos de continuação. Este artigo apresenta novos algoritmos fundamentados nos métodos de Potra-Pták, Chebyshev e super-Halley, associados ao método de continuação Comprimento de Arco Linear. A principal motivação para utilizar tais métodos é a convergência de ordem cúbica. Para elucidar o potencial de nossa abordagem são apresentadas análises de problemas de treliças planas e espaciais com a não-linearidade geométrica encontrados na literatura. Nessa direção, o Método de Elementos Finitos Posicional é utilizado, cuja abordagem considera as coordenadas nodais como varíaveis do sistema não-linear ao invés dos deslocamentos. Os resultados numéricos das simulações mostram a capacidade do código computacional desenvolvido em obter o caminho de equilíbrio com pontos limites de força e deslocamento. Os métodos iterativos implementados exibem melhor eficiência quanto aos números de passos de força e iterações acumuladas necessárias até a convergência para a solução e o tempo de processamento, em comparação com os métodos clássicos de Newton-Raphson e Newton-Raphson Modificado.

Palavras-chave:
Comprimento de Arco; Elementos Finitos Posicional; Chebyshev; Potra- Pták; Não-Linearidade Geométrica

1 INTRODUCTION

In order to realize nonlinear analysis of structures with greater accuracy, it is extremely important that methods which adequately consider the effects of large rotations and displacements are employed. An efficient methodology for solving systems of nonlinear equations must be able to overcome the numerical problems associated with nonlinear behavior, tracing the entire equilibrium path (load versus displacement curve) of the structural system under analysis, identifying and passing through all critical points ^27),(21, In Figure 1, the nonlinear response of a structural system can be seen, in which a given displacement component may increase or decrease along the path. In this figure, load limit points (A, D), displacement limit points (B, C) and failure point (E) are identified ¹⁹.

Figure 1:
Equilibrium path of a structural system.

The Newton-Raphson method is one of the most used methods to solve non-linear problems in Structural Engineering ¹⁴. This method provides the solution of points in the equilibrium path by means of an incremental-iterative procedure. The solution near of a limit point in the path may diverge, due to the ill conditioning of the tangent stiffness matrix, or just because, for the established load level, there is no solution ²⁸.

In order to solve these problems of convergence, path-following techniques associated to this method have been developed, being possible to emphasize: Displacement Control ³, Generalized Displacement Control ^31),(33 and Arc-Length Control ^{25), (26), (24),(8}. In the Constant Load Control technique, the load parameter is kept invariable during the iterative cycle. The idea of path-following methods is to treat the load parameter as a variable, adding a constraint condition to the system of equations that describes the structural equilibrium for determination of oneself.

Modifications in the Newton-Raphson method can be made: solving the system of nonlinear equations in an inexact way, in other words, solve it by some iterative method imposing a precision, as inexact Newton method; approximating the Jacobian matrix using finite differences; and replacing the Jacobian with another matrix with some property, as in the Quasi-Newton methods. In Newton-Raphson method, a linear system is solved at each iteration, whose stiffness matrix is the Jacobian one, evaluated in the current iterated. One of the advantages of this method is the quadratic convergence rate (under suitable conditions). In addition, the optimum radius of convergence of this method is known in the literature; this means that, given a sequence generated by the Newton-Raphson method (whose initial point is outside the center ball in a solution and optimal radius), there is no guarantee that this sequence will converge to the respective solution. However, taking any initial point inside this ball, not only the convergence is guaranteed, but also the quadratic convergence rate.

Until the 1980s, iterative methods that had a higher order of convergence than Newton-Raphson method required the computation of higher order derivatives. Thus, the greater the computational cost making the practical use of these me-thods restricted to few cases. There are methods that have a cubic convergence rate, according to^5),(12, which are better than Newton-Raphson method in this aspect, such as methods belonging to the Chebyshev-Halley class and the Potra-Pták method¹.

In this paper we present algorithms for the incremental and iterative procedures based on Chebyshev, super-Halley and Potra-Pták methods associated with the Linear Arc-Length path-following technique (Riks-Wempner algorithm). These procedures are evaluated in a range of trusses problems with geometric nonlinearity, whose complex equilibrium paths have limit points. The Finite Element Positional formulation is used, which considers nodal positions as variables of the nonlinear system instead of displacements. The numerical results show that the greater computational effort of these methods (the cost of iteration of the Chebyshev, Potra-Pták and super-Halley methods is more expensive than the Newton-Raphson method) is compensated by a smaller number of time steps and iterations required by convergence to a given tolerance. The performance of the methods implemented with Matlab software is evaluated based on the following parameters: total number of time steps (TS), total number of iterations (k _t ), and processing time. The Finite Element Positional Formulation is used, which considers nodal positions as variables of the nonlinear system instead of displacements.

It is important to note that the classical literature deals with Newton-Raphson method combined with Riks Arc-Length techniques. However, by changing Newton-Raphson method by, for example, Chebyshev’s method, a new strategy following the Arc-Length philosophy must be adapted. In this sense, although the problems presented in this work are commonly found in the literature, the techniques of resolution are different from previous studies. Futhermore, the Chebyshev-Halley class depends on calculating tensor, which is usually expensive. In this direction, we use a simple approximation for the tensor based on ¹⁰ to reduce the computational costs. Hence we have the relevance of our exposure.

2 FORMULATION OF THE SPACE TRUSS FINITE ELEMENT

This section describes the truss element using the Positional Finite Element formulation ⁷. This element transmits only axial load and has a constant cross-sectional area A.

The coordinates (X ₁, Y ₁, Z ₁) and (X ₂, Y ₂, Z ₂) represent the initial configuration of the bar element (also known as reference coordinates). After a configuration change due to truss displacements, the bar will have new coordinates (x ₁, y ₁, z ₁) and (x ₂, y ₂, z ₂), according to the schematic drawing in Figure 2.

Figure 2:
Space truss finite element.

The initial (or reference) length L ₀ and the current length of the bar L are calculated, respectively, by:

$L_0=\sqrt{(X_2-X_1{)}^2+(Y_2-Y_1{)}^2+(Z_2-Z_1{)}^2},$ (2.1)

$L=\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2+(z_2-z_1{)}^2}.$ (2.2)

The tangent stiffness matrix K _el and the elementary internal forces vector F _el are given by, respectively:

${\mathbf{K}}_{el}=\frac{EA}{L_0^3}\mathbf{B}+\frac{EA·{\epsilon }_G}{L_0}\mathbf{C},$ (2.3)

${\mathbf{F}}_{el}=\frac{EA·{\epsilon }_G}{L_0}\mathbf{d},$ (2.4)

where EA is the axial stiffness and ( _G is the Green Strain given by:

${\epsilon }_G=\frac{L^2-L_0^2}{2L_0^2}.$ (2.5)

In equations (2.3) and (2.4), the matrices B and C are defined, respectively, as follows:

$\mathbf{B}=\left[\begin{array}{cc}{\mathbf{I}}_3& -{\mathbf{I}}_3\\ -{\mathbf{I}}_3& {\mathbf{I}}_3\end{array}\right],$ (2.6)

were I ₃ is the identity matrix of order 3, C = dd ^T and d = [x ₁ - x ₂, y ₁ - y ₂, z ₁ - z ₂, x ₂ - x ₁, y ₂ - y ₁, z ₂ - z ₁]^T

3 NONLINEAR STRUCTURAL PROBLEM

The basic problem in nonlinear analysis is to find the equilibrium configuration of a structure that is under the action of applied forces. The system of nonlinear equations to be solved and that to governs the static equilibrium of a truss with geometrically nonlinear behavior is given by ²¹:

$\mathbf{g}(\mathbf{d},\lambda )={\mathbf{F}}_{int}\left(\mathbf{d}\right)-\lambda {\mathbf{F}}_r=\mathbf{0},$ (3.1)

being d the vector of coordinate in the nodal points of the structure, g the unba-lanced forces vector, F _int the internal forces vector (determined by vector d) and λ the load parameter responsible for the scaling of the reference vector F _r , that is a reference of arbitrary magnitude.

The equation (3.1) is a system of (n+1) unknowns, n nodal coordinate components (d) and one load parameter (λ), but only n equations. Therefore, an additional constraint equation must be added to the system, given by:

$c(\mathbf{d},\lambda )=0.$ (3.2)

The solution of the nonlinear system given by equations (3.1) and (3.2) is obtained by using an incremental and iterative scheme. Assuming the existence of solution (^t d, ^t λ) in time step t, we need to find (∆ d,∆λ ) such that:

$g\left({}^{t}d+\Delta d,{}^t\lambda +\Delta \lambda \right)=0,$ (3.3)

$c\left({}^{t}d+\Delta d,{}^t\lambda +\Delta \lambda \right)=0,$ (3.4)

Denoting

$g^{\left(k\right)}=g\left({}^{t}d+\Delta d^{\left(k\right)},{}^t\lambda +\Delta {\lambda }^{\left(k\right)}\right)=0,$ (3.5)

and

$c^{\left(k\right)}=g\left({}^{t}d+\Delta d^{\left(k\right)},{}^t\lambda +\Delta {\lambda }^{\left(k\right)}\right)=0,$ (3.6)

for all k =1,2,..., where ∆ d ^(k) is the incremental nodal coordinates vector and ∆λ ^(k) is the load increment. Note that the right superscript k is used herein to refer to the current iteration and the left superscript t refers to the previous time step. The application of the Newton-Raphson method for (3.3)-(3.4) results in the following iterative equations ²:

$\left[\begin{array}{cc}\frac{\partial {\mathbf{g}}^{(k-1)}}{\partial \mathbf{d}}& \frac{\partial {\mathbf{g}}^{(k-1)}}{\partial \lambda }\\ {\left(\frac{\partial c^{(k-1)}}{\partial \mathbf{d}}\right)}^T& \frac{\partial c^{(k-1)}}{\partial \lambda }\end{array}\right]\left[\begin{array}{c}\delta {\mathbf{d}}^{\left(k\right)}\\ \delta {\lambda }^{\left(k\right)}\end{array}\right]=-\left[\begin{array}{c}{\mathbf{g}}^{(k-1)}\\ c^{(k-1)}\end{array}\right],$ (3.7)

with

$\left[\begin{array}{c}\Delta {\mathbf{d}}^{\left(k\right)}\\ \Delta {\lambda }^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}\Delta {\mathbf{d}}^{(k-1)}\\ \Delta {\lambda }^{(k-1)}\end{array}\right]+\left[\begin{array}{c}\delta {\mathbf{d}}^{\left(k\right)}\\ \delta {\lambda }^{\left(k\right)}\end{array}\right].$ (3.8)

Let us denote d ^(k) = ^t d + ∆ d ^(k) , for all k = 1,2,.... Consequently, we have that $\frac{\partial {\mathbf{g}}^{(k-1)}}{\partial \mathbf{d}}=\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)$ is the tangent stiffness matrix (determined by the vector d ^(k-1) = ^t d + ∆ d ^(k-1) ), $\frac{\partial {\mathbf{g}}^{(k-1)}}{\partial \lambda }=-{\mathbf{F}}_r$ is the reference vector, δλ ^(k) is the load sub-increment and δ d ^(k) is the sub-incremental nodal coordinates vector. The first equation of the system (3.7) provides the equation:

$\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\delta {\mathbf{d}}^{\left(k\right)}-\delta {\lambda }^{\left(k\right)}{\mathbf{F}}_r=-{\mathbf{g}}^{(k-1)}.$ (3.9)

Assuming that the inverse of the matrix K(d ^(k-1) ) exists, an expression can be obtained for δ d ^(k) from equation (3.9) ⁸:

$\delta {\mathbf{d}}^{\left(k\right)}=\delta {\mathbf{d}}_g^{\left(k\right)}+\delta {\lambda }^{\left(k\right)}\delta {\mathbf{d}}_r^{\left(k\right)},$ (3.10)

where δd _g ^(k) and δd _r ^(k) are obtained by, respectively:

$\delta {\mathbf{d}}_g^{\left(k\right)}=-{\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}{\mathbf{g}}^{(k-1)},$ (3.11)

$\delta {\mathbf{d}}_r^{\left(k\right)}={\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}{\mathbf{F}}_r,$ (3.12)

where g ^(k-1) = F _int (d ^(k-1) )-λ ^(k-1) F _r . The Newton-Raphson method starts with a first guess (∆d ⁽⁰⁾, ∆λ ⁽⁰⁾) for an approximated solution of the system. As such, for an sequence of the load increment ∆λ ^(k) , the respective nodal coordinate increment ∆d ^(k) is calculated. The process continues until a sufficiently accurate solution is reached. The total parameters of load (λ ^(k) ) and nodal coordinates (d ^(k) ) are update by, respectively:

${\lambda }^{\left(k\right)}={}^t\lambda +\Delta {\lambda }^{\left(k\right)},$ (3.13)

$d^{\left(k\right)}={}^{t}d+\Delta d^{\left(k\right)},$ (3.14)

For a truss consisting of ne finite elements of bar, the global stiffness matrix (K) and the global internal forces vector (F _int ) are obtained from the stiffness matrix (K _el ) and internal forces vector (F el) of each element, respectively, such that stated ²:

$\mathbf{K}=\sum _{i=1}^{ne}\mathbf{A}{\left[{\mathbf{K}}_{el}\right]}_i,$ (3.15)

${\mathbf{F}}_{int}=\sum _{i=1}^{ne}\mathbf{A}{\left[{\mathbf{F}}_{el}\right]}_i,$ (3.16)

where A is an assembly operator.

4 ARC-LENGTH METHOD

The methodology for the solution of nonlinear structural problems must be able to trace the complete equilibrium path, identifying and computing the limit points. For this, an incremental-iterative process is used and consists of two steps ¹⁸:

1. from the last equilibrium configuration of the structure, an increment of load is selected (defined as initial load increment - ∆λ ⁽⁰⁾), trying to satisfy some constraint equation imposed on the problem. After selecting this parameter, the initial increment of nodal coordinates (∆d ⁽⁰⁾) is determined; and

2. the second solution step seeks, by means of a path-following strategy, to correct the incremental solution, initially proposed in the previous step, in order to restore the structure equilibrium. If the iterations involve nodal coordinates (d) and the load parameter (λ ), then an additional constraint equation is required. The format of this equation is what distinguishes the various iteration strategies.

In the Arc-Length method proposed by Riks (presented in ^{25), (26}, the iteration path is always kept orthogonal to the initial tangent at each step. The expression for the initial load increment (predicted solution) is given by:

$\Delta {\lambda }^{\left(0\right)}=\frac{\Delta l}{{\Vert }^t\delta {\mathbf{d}}_r\Vert },$ (4.1)

where ∆l represents the arc-length increment and ^t δd _r is the part of δd ^(k) related to the vector F _r in the previous time step. As proposed by ⁸, the increase ∆l can be used as the control parameter in the current time step as follows:

$\Delta l={}^t\Delta l{\left(\frac{Nd}{{}^{t}k}\right)}^{1/2},$ (4.2)

where ^t ∆l represents the arc-length increment in the previous time step, Nd is the number of desired iterations for the convergence of the current iterative process, and ^t k is the number of iterations that were required to converge in the previous time step.

In the subsequent iterative process, the constraint equation c ^(k) used to calculate δλ ^(k) is obtained making the iterative solution (δd ^(k) , δλ ^(k) F _r ) results in an orthogonal vector to the predicted incremental solution (∆d ⁽⁰⁾, ∆λ ⁽⁰⁾ F _r ), in other words:

$c^{\left(k\right)}=\delta {{\mathbf{d}}^{\left(k\right)}}^T\Delta {\mathbf{d}}^{\left(0\right)}+\delta {\lambda }^{\left(k\right)}\Delta {\lambda }^{\left(0\right)}{{\mathbf{F}}_r}^T{\mathbf{F}}_r=0.$ (4.3)

Substituting equation (3.10) into equation (4.3) yields an expression for the determination of the load parameter correction (k > 1):

$\delta {\lambda }^{\left(k\right)}=-\frac{\Delta {{\mathbf{d}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}}^T\delta {\mathbf{d}}_g^{\left(k\right)}}{\Delta {{\mathbf{d}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}}^T\delta {\mathbf{d}}_r^{\left(k\right)}+\Delta {\lambda }^{\left(0\right)}{{\mathbf{F}}_r}^T{\mathbf{F}}_r}.$ (4.4)

Neglecting the second term in the denominator of equation (4.4), namely, ∆λ ⁽⁰⁾ F _r ^T F _r = 0, the parameter δλ ^(k) results:

$\delta {\lambda }^{\left(k\right)}=-\frac{\Delta {{\mathbf{d}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}}^T\delta {\mathbf{d}}_g^{\left(k\right)}}{\Delta {{\mathbf{d}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}}^T\delta {\mathbf{d}}_r^{\left(k\right)}}.$ (4.5)

The signal of the initial load increment ∆λ ⁽⁰⁾ can be positive or negative. The correct choice of signal is so important in definition of the sequences of solutions (d,λ ) that allow continuous advancement in load-displacement response. The procedure used consists of the analysis of the inner product between the incremental nodal coordinates obtained in the previous load step (^t ∆d) and the initial nodal coordinates sub-increment (δd _r ⁽⁰⁾), so if ^t ∆d ^T δd _r ⁽⁰⁾ > 0, the predictor ∆d ⁽⁰⁾ will have the same direction of δd _r ⁽⁰⁾, otherwise, the predictor will have opposite direction. The convergence criterion is expressed by the unbalanced forces vector norm and reference vector norm applied by:

$\Vert {\mathbf{g}}^{\left(k\right)}\Vert \le tol·\Vert {\mathbf{F}}_r\Vert ,$ (4.6)

where tol is a tolerance provided by the user. The displacements vector (u) is determined as follows:

$u^{\left(k\right)}=d^{\left(k\right)}-{}^{0}d,$ (4.7)

where ⁰ d is the nodal coordinate vector at time step 0 (undeformed structure).

5 SOLUTION METHODS

The basic problem of nonlinear analysis is to find the equilibrium configuration of a structure that is under the action of applied load. In this section, we present the formulations and algorithms for Newton-Raphson, Potra-Pták, Chebyshev and super-Halley methods adapted to the nonlinear structural problem described in Section 3.

The computational tests and algorithms developed for this problem were implemented using the Matlab software, version 8.6 R2015b ²⁰. The algorithms are implemented in a computer Core i7 with 8GB of memory.

5.1 The Newton-Raphson Method

The Newton-Raphson method (NR) only provides the solution of an one point in the of equilibrium path. For other points, the iterations of this method are combined with an incremental procedure. The iteration procedure of this method to determine the solution of the equations system given in (3.1) and (3.2) is described by ²:

$\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\delta {\mathbf{d}}^{\left(k\right)}=\delta {\lambda }^{\left(k\right)}{\mathbf{F}}_r-{\mathbf{g}}^{(k-1)}.$ (5.1)

Note that, isolating δd ^(k) in previous equation, we obtain the equation (3.10). The load sub-increment δλ ^(k) in equation (5.1) is determined as equation (4.5). The algorithm for the Standard Newton-Raphson method associated with Linear Arc-Length method is shown below ^{25), (26), (30}.

Algorithm 1:
Newton-Raphson Method associated with the Linear Arc-Length method (Riks-Wempner algorithm)

In the Modified Newton-Raphson method (MNR), the stiffness matrix K is evaluated at the beginning of the current time step (k = 0), remaining invariable along the iterative cycle (k = 1,...,i _max ).

5.2 Potra - Pták Method

The Potra-Pták method (PP), proposed in 1984, is a two-step method with cubic convergence according to ²³. The idea consists of two evaluations of the given system requiring only the calculation of first order derivatives. This method, adapted to the nonlinear structural problem, is described by the following equation:

$\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\delta {\mathbf{d}}^{\left(k\right)}=\delta {\lambda }_1^{\left(k\right)}{\mathbf{F}}_r-{\mathbf{g}}^{(k-1)}+\delta {\lambda }_2^{\left(k\right)}{\mathbf{F}}_r-\mathbf{g}\left({\mathbf{y}}^{(k-1)}\right).$ (5.2)

Assuming that the inverse of K(d ^(k-1) ) exists, the vector y ^(k-1) is calculated by:

${\mathbf{y}}^{(k-1)}={\mathbf{d}}^{(k-1)}+{\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}\left[\delta {\lambda }_1^{\left(k\right)}{\mathbf{F}}_r-{\mathbf{g}}^{(k-1)}\right].$ (5.3)

Isolating δd ^(k) in equation (5.2), we arrive at the expression:

$\delta {\mathbf{d}}^{\left(k\right)}=\delta {\mathbf{d}}_1^{\left(k\right)}+\delta {\mathbf{d}}_2^{\left(k\right)},$ (5.4)

where

$\delta {\mathbf{d}}_1^{\left(k\right)}={\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}\left[\delta {\lambda }_1^{\left(k\right)}{\mathbf{F}}_r-{\mathbf{g}}^{(k-1)}\right]=\delta {\lambda }_1^{\left(k\right)}\delta {\mathbf{d}}_r^{\left(k\right)}+\delta {\mathbf{d}}_g^{\left(k\right)},$ (5.5)

$\delta {\mathbf{d}}_2^{\left(k\right)}={\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}\left[\delta {\lambda }_2^{\left(k\right)}{\mathbf{F}}_r-\mathbf{g}\left({\mathbf{y}}^{(k-1)}\right)\right]=\delta {\lambda }_2^{\left(k\right)}\delta {\mathbf{d}}_r^{\left(k\right)}+\delta {\mathbf{y}}_g^{\left(k\right)}.$ (5.6)

The algorithm developed in this study for Potra-Pták method, adapted to the Linear Arc-Length method, is presented as follows.

Algorithm 2:
Potra-Pták method associated with the Linear Arc-Length method

5.2.1 Evaluation of δλ ₁ ^(k) and δλ ₂ ^(k)

A constraint equation, proposed by ³², must be respected in the two nonlinear solution steps - the predicted solution stage and the iteration cycle. The constraint equation imposed on the problem is of the form, adapted in this article:

${\mathbf{C}}_i^T\delta {\mathbf{d}}^{\left(k\right)}+k_i\delta {\lambda }_i^{\left(k\right)}=H_i,$ (5.7)

where C ₁ is a vector whose elements are constant, k _i is a scalar constant and H _i is an incremental parameter that can be displacement, arc-length or external work. Depending on the value adopted for these parameters, one can obtain different load and iteration strategies. The parameters δλ₁ ^(k) and δλ₂ ^(k) are calculated by, respectively:

$\delta {\lambda }_1^{\left(k\right)}=\frac{H_1-{\mathbf{C}}_1^T\delta {\mathbf{d}}_g^{\left(k\right)}}{{\mathbf{C}}_1^T\delta {\mathbf{d}}_r^{\left(k\right)}+k_1},$ (5.8)

$\delta {\lambda }_2^{\left(k\right)}=\frac{H_2-{\mathbf{C}}_2^T\delta {\mathbf{y}}_g^{\left(k\right)}}{{\mathbf{C}}_2^T\delta {\mathbf{d}}_r^{\left(k\right)}+k_2}.$ (5.9)

Considering the parameters k ₁ = k ₂ = 0, H ₁ = H ₂ = 0 and C ₁ = C ₂ = ∆u ⁽⁰⁾ in equations (5.8) and (5.9), which results in the Linear Arc-Length method.

5.3 Super-Halley and Chebyshev Methods

The iterative methods of Chebyshev-Halley family ²⁹, adapted to the nonlinear structural problem, are described by the equation:

${\mathbf{d}}^{\left(k\right)}={\mathbf{d}}^{(k-1)}+\delta {\mathbf{d}}^{\left(k\right)}+\frac{1}{2}{\mathbf{L}}^{(k-1)}{\left[\mathbf{I}-\gamma {\mathbf{L}}^{(k-1)}\right]}^{-1}\delta {\mathbf{d}}^{\left(k\right)},$ (5.10)

where γ is a real parameter, I is the identity matrix and L ^(k-1) is a matrix given by ^{12), (11), (9}:

${\mathbf{L}}^{(k-1)}={\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}{\left[{\mathbf{H}}^{(k-1)}\delta {\mathbf{d}}^{\left(k\right)}\right]}^T,$ (5.11)

where H ^(k-1) is a tensor. This family includes some methods, such as Chebyshev method (γ =0) and the super-Halley method (γ =1), which have cubic convergence. Making γ =0 in the equation (5.10), the Chebyshev method (Cbs) is given by the following equation ¹³:

${\mathbf{d}}^{\left(k\right)}={\mathbf{d}}^{(k-1)}+\delta {\mathbf{d}}^{\left(k\right)}+\frac{1}{2}{\mathbf{L}}^{(k-1)}\delta {\mathbf{d}}^{\left(k\right)}.$ (5.12)

Similarly, by making γ =1 in the equation (5.10), it is obtained the super-Halley method (sH) described by:

${\mathbf{d}}^{\left(k\right)}={\mathbf{d}}^{(k-1)}+\delta {\mathbf{d}}^{\left(k\right)}+\frac{1}{2}{\mathbf{L}}^{(k-1)}{\left(\mathbf{I}-{\mathbf{L}}^{(k-1)}\right)}^{-1}\delta {\mathbf{d}}^{\left(k\right)}.$ (5.13)

In this study, since we do not have access to the higher order derivatives that form the tensor H, we use an approximation for the matrix L ^(k-1) adapted from the research of ¹⁰:

${\mathbf{L}}^{(k-1)}\approx -{\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right]}^{-1}\frac{1}{p}\left[\mathbf{K}\left({\mathbf{d}}^{(k-1)}+p\delta {\mathbf{d}}^{\left(k\right)}\right)-\mathbf{K}\left({\mathbf{d}}^{(k-1)}\right)\right],$ (5.14)

where p ∈ (0, 1]. This approximation is relatively simple to calculate, and involves the stiffness matrices calculated in the points (d ^(k-1) + p δd ^(k) ) and d ^(k-1) .

The algorithm developed for the super-Halley and Chebyshev methods asso-ciated with the Linear Arc-Length method is presented as follows.

Algorithm 3:
Super-Halley and Chebyshev methods associated with the Linear

6 RESULTS AND DISCUSSION

In this section, we present the numerical results of space and plane trusses problems found in the literature, taking into account in the static analysis the geometric nonlinearity, in order to apply and compare the solution methods. The own weight of the structures is neglected in the analyzes. It is also worth noting that the generation of the mesh and the results are not counted in the processing time. In the simulations with the super-Halley and Chebyshev methods is considered p = 1 in equation (5.14) for the calculation of the matrix L.

6.1 Plane truss with two bars

Consider in Figure 3a the plane truss composed of two bars with axial stiffness EA = 1884.694 lbf, studied by¹⁹ and ³¹. The geometric imperfection is introduced by the horizontal (0.05P) and vertical (P) loads applied at node 2. The parameters considered for the path-following technique were: ∆l =1.0; Nd = 2; tol = 1.0×10^-6; i _max = 100; and ∆P = 10.0 lbf. Table 1 contains the numerical results (TS, k _t and processing time) of the simulations performed with the Newton-Raphson (NR), Modified Newton-Raphson (MNR), Chebyshev (Cbs), super-Halley (sH) and Potra-Pták (PP) methods. The vertical displacement at node 2 versus load P curves obtained with the Cbs, sH and PP methods, comparing them with that obtained by ¹⁹, are shown in the Figure 3b.

Figure 3:
a) Plane truss with two bars and b) Equilibrium paths.

Table 1:
Numerical results of the plane truss with 2 bars.

Methods	TS	k t	Time (s)
NR	330	667	0.2246
MNR	1085	2191	0.4713
Cbs	144	288	0.1735
sH	129	257	0.2059
PP	111	221	0.1426

6.2 Plane truss with 22 bars

Consider the plane truss with 22 bars illustrated in Figure 4a, whose bars have elasticity modulus E = 3.0 × 10⁴ kip/in² and cross-sectional area equal to 20 and 40 in² for diagonals and the other members, respectively. This structure was stu-died by ²². The parameters used in the simulations were: ∆l =45.0; Nd = 2; tol = 1.0 × 10^-7; i _max = 100; and ∆P = 8.0 kip. The equilibrium paths (vertical displacement at node 13 versus load P curves) are shown in Figure 4b. The numerical results (TS, k _t and processing time) of the simulations with the implemented methods are shown in Table 2.

Figure 4:
a) Plane truss with 22 bars and b) Equilibrium paths.

Table 2:
Numerical results of the plane truss with 22 bars.

Methods	TS	k t	Time (s)
NR	1172	2355	2.4989
MNR	5336	10696	9.5401
Cbs	479	966	1.5305
sH	263	530	1.0473
PP	234	471	0.8647

6.3 Space truss with 12 bars

Consider the 12-bar space truss, shown in Figure 5a, whose three top nodes are free and the remaining six are fixed. The complete path was obtained in ¹⁵, ¹⁶ and ¹⁷. All bars have the same dimensionless axial stiffness EA = 1.0. At the top of the structure a load F is applied at the central node and loads of intensity 1.5F at the two adjacent nodes. For the Arc-Length strategy, the following parameters were used: ∆l = 2.0 × 10^-3; Nd = 2; i _max = 100; tol = 1.0 × 10^-6; and ∆F = 1.0. Figure 5b shows the vertical displacement at node 7 versus load F curves obtained with the Cbs, sH and PP methods, being close to the equilibrium path obtained by ¹⁵. The values of TS, k _t and processing time obtained from the simulations are presented in Table 3.

Figure 5:
a) Space truss with 12 bars and b) Equilibrium paths.

Table 3:
Numerical results of the space truss with 12 bars.

Methods	TS	k t	Time (s)
NR	247	485	0.5222
MNR	306	605	0.6147
Cbs	182	354	0.5146
sH	178	346	0.5742
PP	87	161	0.2794

6.4 Star Dome Space Truss with 24 bars

A star space dome truss with 24 bars and axial stiffness EA = 8.0 ×10⁷ N is shown in Figure 6a. In the top view of the dome, the outer circle has a radius of 50 cm and height zero, and the inner circle has a radius of 25 cm and a height of 6.216 cm. At the apex of the truss, with height 8.216 cm, a vertical load P is applied. Figure 6b shows the equilibrium paths (vertical displacement at the apex of the truss versus load P curves) obtained with the Cbs, sH and PP methods, being compatible with that obtained by ⁴. The parameters used in the simulations were: ∆l = 1.25; Nd = 2; i _max = 100; tol = 1.0 × 10^-6 ; and ∆F = 100.0 N. The values for the parameters TS, k _t and processing time obtained from the simulations are shown in Table 4.

Figure 6:
a) Star Dome Space Truss with 24 bars and b) Equilibrium paths.

Table 4:
Numerical results of the space truss with 24 bars.

Methods	TS	k t	Time (s)
NR	1131	2280	5.0963
MNR	No convergence
Cbs	357	726	1.5443
sH	162	331	0.9706
PP	179	366	0.8752

6.5 Star dome truss with 60 bars

The fifth structure analyzed corresponds to a trussed dome structure studied by ¹⁹ and ⁶, whose structural model is presented in Figure 7a. This structure has 25 nodes and 60 elements with dimensionless axial stiffness EA = 1.0 × 10⁴, being submitted to six vertical loads (P) of equal magnitude applied on nodes 13 to 18.

Figure 7:
a) Structural model of the star dome truss with 60 bars and b) Equilibrium paths.

The equilibrium paths (vertical displacement at the top of the dome versus load P) obtained with the Cbs, sH and PP methods are shown in Figure 7b, showing good agreement with the response obtained by ¹⁹. The parameters used for the path-following strategy were: ∆l = 3.0; Nd = 2; i _max = 100; ∆P = 0.1; and tol = 1.0 × 10^-5. Table 5 shows the numerical results of the simulations.

Table 5:
Numerical results of the space truss with 60 bars.

Methods	TS	k t	Time (s)
NR	832	1674	5.2363
MNR	No convergence
Cbs	287	578	3.2959
sH	193	388	2.5548
PP	212	427	2.5765

6.6 Discussion of results

The analyzed structures are characterized by a strongly nonlinear behavior, in which the equilibrium paths exhibit force limit points (maximum and minimum points) and displacement limit points (the tangent at these points is vertical). When the structure reaches the limit points, it can become unstable; therefore, their identification is of great importance for an engineering design.

The procedure used to change the signal of the initial load increment (∆λ ⁽⁰⁾) was efficient, since it was able to identify and overcome the limit points in the equilibrium paths, defining its correct direction. In addition, the device for verification of signal change is easy to implement computationally. There are other strategies for checking the signal changes, such as the variation in the signal of the stiffness matrix determinant, or the change in the sign of some eigenvalue of this matrix.

The super-Halley, Chebyshev and Potra-Pták methods proved to be more computationally efficient than the Newton-Raphson method, since the number of time steps and cumulative iterations required for convergence to the approximate solution of the problem were smaller, although the computational cost of the iteration of these methods is greater than the iteration of the Newton-Raphson method. With this, there was a reduction in the amount in which the systems of linear equations (generated from the Finite Element discretization) are solved in the incremental process, and also in the updates of the stiffness matrix (K) and the internal force vector (F _int ) during the iterative cycle.

It is noted that in the iteration of the super-Halley, Chebyshev and Potra-Pták methods the same stiffness matrix is used for the solution of the systems of li-near equations; thus, these systems can be solved via decomposition (for example, LU decomposition), since a single factorization at the beginning of the iteration is necessary.

The Chebyshev and super-Halley methods are very attractive for solving nonlinear problems, because they have the cubic convergence rate. However, they are computationally expensive for two reasons: (i) the need to obtain the matrix L at each iteration; and (ii) the resolution of two linear systems. The analyzes with the Chebyshev and super-Halley methods with the approximation for the matrix L adapted from ¹⁰ were promising, since they were able to reach the response (equilibrium path) with less number of time steps and iterations, if compared to the Newton-Raphson and Modified Newton-Raphson methods in the numerical examples studied.

In the analysis performed with the Modified Newton-Raphson method, there was a numerical instability occurring the non-convergence for the solution of the last two problems. In fact, problems of convergence with this method can occur during the analysis, since the stiffness matrix K is computed only at the beginning of each time step, remaining invariable along the iterative cycle. However, the convergence can be achieved by increasing the tolerance value (tol), decreasing the value of the arc-length increment (∆l) or decreasing the value of load increment (∆P or ∆F).

The matrix of the structural system is characterized by a high sparsity index. It is possible to obtain a better numerical efficiency of the presented models, by means of algorithms that store the non-null coefficients present in the matrix and perform operations between matrices and vectors with these coefficients, thus avoiding redundant calculations involving null elements.

7 CONCLUSION

In nonlinear simulations of structural problems in an incremental and iterative process, solving the system of linear equations generated at each iteration is, in general, the most expensive step, and which demands greater computational time and effort during processing. Even with the impact of the microelectronics industry on the development of computational components, with the emphasis on more compact memory systems and faster processors, these powerful machines alone cannot always adequately handle various structural models, either for lack of memory or for excessive response time.

From the numerical examples studied, we highlight the good agreement between the results obtained and those of the literature on the achievement of equilibrium paths, validating the developed computational codes. The increasing simulation of complex structural models - through the Finite Element Method - has required the manipulation of a large amount of data, which is intrinsic to the method, as well as the search for a decrease in the response time to solve the structural problem.

It is observed a better computational performance in all the analyzes with the super-Halley, Chebyshev and Potra-Pták methods in comparison with the classic methods of Newton-Raphson and Newton-Raphson Modified, associated with the Linear Arc-Length path-following technique proposed by ^{25), (26}. Moreover, in the simulations with these methods, no numerical instability occurred, converging to the solution in all problems of nonlinear trusses studied. The Chebyshev method, whose iteration is computationally cheaper than super-Halley method, is a good alternative in the solution of nonlinear problems.

As future research, it is suggested to implement algorithms that consider the physical nonlinearity and to adapt the codes implemented for studies in dynamic analysis. In addition, numerical studies with other approximations for the matrix L in the Chebyshev and super-Halley methods can be carried out.

ACKNOWLEDGMENTS

The authors wish to thank the Graduate Program in Mathematics of the State University of Maringá and the Federal Technological University of Paraná for their support in the development of this research. The authors also wish to thank the two anonymous referees whose comments helped shape the final version of this paper.

REFERENCES

• ¹
  S. Amat, S. Busquier & J.M. Gutiérrez. Geometric constructions of iterative functions to solve nonlinear equations. Journal of Computational and Applied Mathematics, 157 (2003), 197-205.
• ²
  K. Bathe. Finite Element Procedures. Prentice Hall (2006). URL https://books.google.com. br/books?id=rWvefGICfO8C
  » URL https://books.google.com. br/books?id=rWvefGICfO8C
• ³
  J.L. Batoz & G. Dhat. Incremental displacement algorithms for nonlinear problems. International Journal for Numerical Methods in Engineering, 14 (1979), 1262-1267.
• ⁴
  J. Bonet, A.J. Gil&R.D.Wood. Worked examples in nonlinear continuum mechanics for finite element analysis. Cambridge University Press (2012).
• ⁵
  V. Candela & A. Marquina. Recurrence relations for rational cubic methods II: the Chebyshev method. Computing, 45(4) (1990), 355-367.
• ⁶
  K.K. Choong & Y. Hangai. Review on methods of bifurcation analysis for geometrically nonlinear structures. Bulletin of the International Association for Shell and Spatial Structures, 34(112) (1993), 133-149.
• ⁷
  H.B. Coda & M. Greco. A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering, 193 (2004), 3541-3557.
• ⁸
  M.A. Crisfield. Non-Linear Finite Element Analysis of Solids and Structures: Essentials. John Wiley & Sons, Inc., New York, NY, USA (1991).
• ⁹
  J. Ezquerro & M. Hernández. A modification of the convergence conditions for Picards iteration. Computational & Applied Mathematics, 23(1) (2004), 55-65.
• ¹⁰
  J. Ezquerro & M. Hernández. An optimization of Chebyshev’s method. Journal of Complexity, 25(4) (2009), 343-361.
• ¹¹
  M. Frontini & E. Sormani. Third-order methods from quadrature formulae for solving systems of nonlinear equations. Applied Mathematics and Computation, 149(3) (2004), 771-782.
• ¹²
  J.M. Gutierrez & M.A. Hernández. An acceleration of Newtons method: Super-Halley method. Applied Mathematics and Computation, 117(2) (2001), 223-239.
• ¹³
  J.M. Gutíerrez & M.A. Hernández. A family of Chebyshev-Halley type methods in Banach spaces. Bulletin of the Australian Mathematical Society, 55(1) (1997), 113-130.
• ¹⁴
  W.E. Haisler, J.A. Stricklin & F.J. Stebbins. Development and evaluation of solution procedures for geometrically nonlinear structural analysis. AIAA Journal, 10(3) (1972), 264-272.
• ¹⁵
  S. Krenk. Non-linear modeling analysis of solids and structures. Cambridge, Cambridge, UK (2009).
• ¹⁶
  S. Krenk & O. Hededal. A dual orthogonality procedure for non-linear finite element equations. Computer Methods in Applied Mechanics and Engineering, 123 (1995), 95-107.
• ¹⁷
  E.G.M. Lacerda, D.N. Maciel & A.C. Scudelari. Geometrically static analysis of trusses using the Arc-Length Method and the positional formulation of Finite Element Method. In XXXV Iberian Latin-American Congress on Computational Methods in Engineering (2014).
• ¹⁸
  S.E. Leon, G.H. Paulino, A. Pereira, I.F.M. Menezes & E.N. Lages. A Unified Library of Nonlinear Solution Schemes. Applied Mechanics Reviews, 64 (2011), 1-26.
• ¹⁹
  W.T. Matias. El control variable de los desplazamientos en el análisis no lineal elástico de estructuras de barras. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 18(4) (2002), 549-572.
• ²⁰
  MATLAB. version 8.6.0 (R2015b). The MathWorks Inc., Natick, Massachusetts (2015).
• ²¹
  D.P. Maximiano, A.R.D. Silva & R.A.M. Silveira. Iterative strategies associated with the normal flow technique on the nonlinear analysis of structural arches. Revista Escola de Minas(Impresso), 67 (2014), 143-150.
• ²²
  M.R. Pajand & M.T. Hakkak. Nonlinear analysis of truss structures using dynamic relaxation. IJE Transactions B: Applications, 19(1) (2006), 11-22.
• ²³
  F.A. Potra & V. Pták. Nondiscrete induction and iterative processes, volume 103. Pitman Advanced Publishing Program (1984).
• ²⁴
  E. Ramm. Strategies for tracing the non-linear response near limit-points, nonlinear finite element analysis in structural mechanics. Springer-Verlag, Berlin, GE (1981).
• ²⁵
  E. Riks. The application of Newtons method to the problem of elastic stability. Journal of Applied Mechanics, 39(4) (1972), 1060-1065.
• ²⁶
  E. Riks. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, 15 (1979), 529-551.
• ²⁷
  P.F.N. Rodrigues, W.D. Varela & R.A. Souza. Análise de Estratégias de Solução do Problema Nãolinear. Revista Ciência e Tecnologia, 8(2) (2008), 36-49.
• ²⁸
  R.A.M. Silveira, G. Rocha & P.B. Gonçalvez. Estratégias numéricas para análises geometricamente não lineares. In XV Congresso Brasileiro de Engenharia Mecânica, Águas de Lindóia , pp. 117-120 (1999).
• ²⁹
  T. Steihaug & S. Suleiman. Rate of convergence of higher order methods. Applied Numerical Mathematics, 67 (2013), 230-242.
• ³⁰
  G.A.Wempner. Discrete approximations related to nonlinear theories of solids. International Journal of Solids and Structures, 7(11) (1971), 1581-1599.
• ³¹
  Y.B. Yang & S.R. Kuo. Theory & Analysis of Nonlinear Framed Structures. Prentice-Hall, Singapore (1994).
• ³²
  Y.B. Yang & S.R. Kuo. Theory and analysis of nonlinear framed structures. Prentice Hall PTR (1994).
• ³³
  Y.B. Yang & M.S. Shieh. Solution method for nonlinear problems with multiple critical points. AIAA Journal, 28(12) (1990), 2110-2116.

Publication Dates

• Publication in this collection
  Jan-Apr 2018

History

• Received
  28 Mar 2017
• Accepted
  02 Mar 2018