Finding buckling points for nonlinear structures by dynamic relaxation scheme

Abstract

Dynamic Relaxation Method (DRM) is an explicit approach for solving the simultaneous systems of equations. In this tactic, the fictitious mass and damping are added to the static governing equations, and an artificial dynamic system is constructed. By using DRM for nonlinear analysis, the structural static equilibrium path is obtained. This outcome is extremely valuable, since it leads to the behavior of structures. Among the finding related to the structural static path are the critical and buckling points for nonlinear structures. In this paper, a new way for calculating the load factor is proposed by setting the external work zero. Mixing the dynamic relaxation scheme with external work technique has not been formulated so far. In all incremental-iterative methods, the load factor increment sign should be determinated by extra calculations. This sign leads to increase or decrease of the load increment. It is worth emphasizing that sign of the load factor increment changes at the load limit points. Therefore, the sign determinations are required in the common work control methods. These disadvantages are eliminated in the proposed algorithm. In fact, the suggested load factor depends only on the Dynamic Relaxation (DR) fictitious parameters. Besides, all DR calculations are performed via vector operation. Moreover, the load factor is calculated only by one formula, and it has the same relation in the all solution processes. In contrast to the arc length techniques, which requires the sign determined, the proposed scheme does not need any sign finding. It is shown that author’s technique is quicker than the other dynamic relaxation strategies. To prove the capability and efficiency of the presented scheme, several numerical tests are performed. The results indicate that the suggested approach can trace the complex structural static paths, even in the snap-back and snap-through parts.

Keywords

load factor / external work / dynamic relaxation / static equilibrium path / large displacement

Introduction

To analyze structures under the static loads, by using stiffness tactics, the coming system of equations should be solved:

(1)\mathit{K}\mathit{X}=\mathit{F}=\mathit{P}\mathrm{.}

In this equality, the stiffness matrix, the displacement vector, the internal force vector, and the external load vector are shown by K, X, F, and P, respectively. It should be reminded that the internal forces are related to nodal displacements in the structural geometric nonlinear analysis. To solve Eq. (1), numerical methods can be deployed. In general, direct schemes, such as Gaussian elimination and Cholesky factorization are applicable in solving a linear system of equations. Nevertheless, these techniques cannot be merely utilized for solving a set of nonlinear equations. In this case, these strategies in the companion with iterative tactics are employed.

Iterative algorithms are divided into two categories: explicit approaches and implicit approaches. In the former group, the internal forces are applied to achieve the responses. As a result, vector calculations are used in these methods. The simplicity and high efficiency are the main characteristics of these tactics. On the other hand, the implicit schemes are formulated based on the stiffness matrix. Due to utilizing matrix operation in these algorithms, the implicit techniques are complex and time consuming. To clarify this issue, the load and displacement limit points, in which the stiffness matrix is zero or undefined, lead to difficulties in solution procedure. Note that the convergence rate of the implicit ways is more, in comparison with the explicit ones [¹].

Newton-Raphson method is not able to trace the structural static path, having the limit load points [²]. To overcome this drawback, the control displacement scheme was introduced [³]. However, this approach is weak in passing the snap-back points. For the first time, Wempner [⁴] and Riks [⁵] proposed the arc-length procedure. The arc-length methods have different categories. One of them is the normal plane tactic, in which; the locus of the obtained iterative solutions is perpendicular to the tangent of the first equilibrium point in the current step [⁶]. The other type is updated normal plane technique, which was presented by Ramm [⁷]. Later, Crisfield [⁸] formulated the cylindrical arc-length scheme. The residual load minimization method was suggested by Bergan [⁹].

In recent years, various strategies are introduced for tracing structural static path. Quasi-newton methods and orthogonal condition were used by Krenk and Hededal for nonlinear analysis of structures [¹⁰]. Rezaiee-Pajand and Boroshaki [¹¹] utilized variable arc length technique to trace the structural static path. Kim and Kim [¹²] deployed neural networks and Newton-Raphson tactic to analyze nonlinear structures. Toklu [¹³] analyzed nonlinear trusses by minimizing the total strain energy. Ligarò and Valvo [¹⁴] minimized the total strain energy of structures to analyze the nonlinear regular pyramidal trusses. Safari and Mansouri proposed a two-point method for solving the nonlinear system of equations [¹⁵].

One of the well-known schemes for nonlinear analysis of structures is Dynamic Relaxation Methock (DRM). This is an explicit procedure. At first, Otter and Day [¹⁶–¹⁸] used these approaches in 1960s. DRM employed the Second-Order Richardson tactic. This way was developed by Frankel [¹⁹]. Initially, Rushton [²⁰] applied DRM in the structural nonlinear analysis. He used this method to geometric nonlinear behavior of the bending plate. DRM can also be utilized for tracing static equilibrium path [²¹–²⁵]. Moreover, this scheme has been widely used for form-finding of structures [²⁶–³⁴].

Recently, researchers have investigated the sensitivity and uncertainty analysis for various structures. Vu-Bac et al. [³⁵] provide a sensitivity analysis toolbox consisting of a set of Matlab functions that offer utilities for quantifying the influence of uncertain input parameters on uncertain model outputs. Furthermore, Hamdia et al. [³⁶] present a methodology for stochastic modeling of the fracture in polymer/particle nanocomposites. In 2018, Hamdia et al. [³⁷] performed a sensitivity analysis to identify the key input parameters influencing the energy conversion factor of flexoelectric materials by a NURBS-based formulation.

It is worth emphasizing that researchers have not been so far successful in presenting an approach capable of tracing all kinds of the static equilibrium paths. In other words, the most robust tactics are not able to trace the equilibrium path of all problems. This paper aims to present a formulation to improve DRM for perfectly tracing the static equilibrium path. A new formula will be proposed for the load factor calculation. In this way, the work increment will be established based on the Dynamic Relaxation (DR) fictitious parameters. Based on the external work approach, the work increment is equal to zero in each iteration. By setting work increment zero, a new relationship will be obtained for the load factor. In fact, the work increment depends on the load factor and the displacement increment. Moreover, in the DRM, the displacement increment is related to the velocity. Also, the velocity ​​depends on the fictitious mass and damping. As a result, the proposed load factor depends only on the DR fictitious parameter. By using of this load factor, the structural static paths can be traced.

Following the introduction section, the basis of DRM will be presented. Afterwards, previously proposed techniques for tracing structural static path will be briefly reviewed. Then, a new method for calculating the load factor will be formulated. Finally, several numerical samples will be solved to assess the efficiency and robustness of the suggested strategy. These include a variety of the geometrical nonlinear problems, such as 2D and 3D trusses, 2D frames, plane arcs and shells.

Dynamic Relaxation Method (DRM)

DRM is an explicit tactic for solving the simultaneous system of equations. In this approach, an artificial mass and damping are added to the structural static equilibrium equations. Afterwards, the fictitious dynamic system will be solved.

(2)M^n{\stackrel{\mathrm{..}}{X}}^n+C^n{\stackrel{.}{X}}^n+K^n{X}^n=F^n=P^n\mathrm{.}

To set up the former equations, M and C are utilized, which denote the fictitious mass matrix and the damping matrix, respectively. Note that these matrices are diagonal. The displacement, velocity and acceleration are demonstrated by X,\stackrel{.}{X}, and \stackrel{\mathrm{..}}{X}, respectively. Moreover, the number of iterations is denoted by n. To solve Eq. ‎(2), the inertia and damping forces should become zero. In other words, the steady-state response of Eq. ‎(2) is the solution of Eq. ‎(1). It is worth emphasizing; the residual force vector leads to artificial oscillation of the structure. This force vector can be computed as below:

(3)R=M\stackrel{\mathrm{..}}{X}+C\stackrel{.}{X}=P-F\mathrm{.}

Obviously, the response of the structure can be achieved under the condition that the residual force vector is equal to zero. In DRM, it is presumed that the velocity changes linearly in each time step. Hence, the acceleration is constant. Accordingly, the iterative relationships utilized in the strategy are obtained with the help of central finite difference tactic as follows:

(4){\stackrel{.}{X}}_i^{n+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{2}{m}_{ii}^n-C_{ii}^n{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}{\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{2}{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}(p_i^n-f_i^n)\begin{array}{l},\end{array}i=\mathrm{1},\mathrm{2},\mathrm{...},ndof,

(5)X_i^{n+\mathrm{1}}=X_i^n+\delta X_i^n=X_i^n+t^{n+\mathrm{1}}{\stackrel{.}{X}}_i^{n+\frac{\mathrm{1}}{\mathrm{2}}}\begin{array}{l},\end{array}i=\mathrm{1},\mathrm{2},\mathrm{...},ndof\mathrm{.}

The fictitious time step, the ith entries of internal force vector, the ith diagonal entries of artificial mass and damping are shown by t^n,f_i^n,m_{ii}^n andC_{ii}^n, respectively. The external load of the static system is denoted byp_i. Furthermore, ndof is the number of degrees of freedom. The displacement increment is demonstrated by\delta X_i^n. Equations ‎(4) and ‎(5) are repeated to reach the stable solution. By assuming that the fictitious mass and damping matrices are diagonal leads to explicit formulations in DRM. It is clear from presenting formulas that the vector operators are merely applied in this technique.

For the first time, Welsh and Cassell proposed the idea of employing the fictitious mass. In addition, they improved the DR iterations by using numerical methods [³⁸,³⁹]. Wood [⁴⁰] specified the artificial mass by deploying the upper bound of the spectral radius of the coefficient matrix. He proved that the convergence rate of DRM is more than the degenerate Chebyshev in the structural analysis. Brew and Brotton showed that DRM is capable of solving a system of linear equations, resulting from the planar frame. Moreover, they concluded that this procedure can consider nonlinear effects. Additionally, these investigators suggested that the mass of each degree of freedom should be proportional to the corresponding diagonal elements of the stiffness matrix [⁴¹]. In another paper, a new approach for estimating the critical damping was proposed by Bunce [⁴²]. Note that Rayleigh's principle was applied in this algorithm.

Initial stability conditions for this solution technique were presented by Cassell and Hobbs. They approximated the fictitious mass with the help of Gerschgörin circle theory. By this way, these researchers assessed the convergence and stability of the nonlinear analysis, based on DRM [⁴³]. Felippa [⁴⁴] formulated the implicit DR method. In this scheme, the mass and damping parameters are computed by varying the structural stiffness in fictitious time. In 1983, Underwood [⁴⁵] suggested a well-known formulation for iterations of the DRM. He assumed that the time step is constant during the iterations process. Qiang [⁴⁶] determined the artificial time and damping by utilizing Rayleigh’s principle. The modified adaptive DR method was proposed by Zhang et al. [⁴⁷]. They suggested the nodal damping template and initial displacement vector. Munjiza et al. [⁴⁸,⁴⁹] presumed that the damping is proportional to the exponent of mass and stiffness matrix. These investigators proved that all modes will be damped critically when the damping matrix is expressed in the form of \mathrm{2}M{(M^{-\mathrm{1}}S)}^{\mathrm{0}.\mathrm{5}}.

To reduce the number of iterations, increase the error rank and improve the convergence rate of DRM, Rezaiee-Pajand and Taghavian Hakkak [⁵⁰] applied Taylor series expansion. With the help of the first three terms of the Taylor series, they obtained the displacement. These researchers assumed that the fictitious time step is constant in all iterations. In 2008, Kadkhodayan et al. [⁵¹] formulated the fictitious time step for DRM by minimizing the residual forces of iterations. Rezaiee-Pajand and Alamatian [⁵²] deployed the DRM for nonlinear dynamic analysis of structures. They were used DR with modified time step. The integration error was reduced in their scheme.

The relation between the optimum time step and the critical damping was determined by Rezaiee-Pajand and Sarafrazi [⁵³]. They proved that the constant time step does not have any effect on the convergence rate. Moreover, Rezaiee-Pajand and Alamatian [⁵⁴] proposed new formulas for the fictitious mass and damping in 2010. Based on their work, it was deduced that the convergence rate can be increased without any additional calculations. Furthermore, Rezaiee-Pajand et al. [⁵⁵] suggested a new technique for computing the fictitious damping. This goal was reached by error minimization between two successive steps and also using Stodola iterative process.

Mass and damping matrix along with the artificial time step play the important role in DRM. By setting the damping zero and specifying the time step ratio, Rezaiee-Pajand and Sarafrazi [⁵⁶] presented a new formulation. To demonstrate this approach's application, they performed the nonlinear dynamic structural analysis. In 2011, based on minimization of the residual energy, Rezaiee-Pajand et al. [⁵⁷] proposed another new method. Recently, Rezaiee-Pajand et al. [⁵⁸] investigated extensively the capability of 12 well-known DRM in the finite element analysis of frames and trusses. Rezaiee-Pajand and Estiri [⁵⁹] proposed a procedure for finding the load factor by imposing the work increment of the external forces to zero. They suggested a new formula for calculating the load factor by minimizing external work [⁶⁰]. Moreover, these researchers calculated another variable load factor by minimizing the unbalanced displacement [⁶¹]. To evaluate the DRM other skills, Rezaiee-Pajand and Estiri [⁶²–⁶⁴] compared the abilities of DR approaches in analyzing bending plates and frames. In another event, Labbafi et al. [⁶⁵] Compared viscous and kinetic DRMs in form-finding of membrane structures.

In DR formulations, the fictitious diagonal mass and damping matrices, the initial displacement vector and the fictitious time step are unknown. Usually, it is presumed that the entries of the initial displacement vector equal zero or one. Furthermore, the fictitious time step is considered to be one. Based on Rayleigh’s principle, other researchers determined the time step [⁴⁶]. The most appropriate approach for calculating the artificial mass matrix is Gerschgorin's circle theory [⁴³,⁴⁵]. Rezaiee-Pajand and Alamatian [⁵⁴] suggested another efficient strategy in which the fictitious mass is calculated automatically. The most common way for identifying the damping parameter is based on Rayleigh’s principle [⁴²,⁴⁷]. Two famous algorithms which utilize the aforesaid principle were proposed by Zhang and Underwood [⁴⁵,⁶⁶]. In another paper, Rezaiee-Pajand et al. [⁵³–⁵⁵] presented a technique with high convergence rate. In the coming lines; the well-known formulas applied to determine the unknown parameters related to the mass and damping will be discussed. In 1972, one of the famous relations used in computing the mass matrix is proposed by Bunce [⁴²]. This relationship has the next shape:

(6)m_{ii}^n={\displaystyle \frac{{(t^n)}^{\mathrm{2}}}{\mathrm{4}}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}\left|k_{ij}\right|}\mathrm{.}

The entries of the stiffness matrix are denoted by k_{ij}. It is obvious that nodal fictitious mass calculated in any iteration from the current values of the coefficients in the stiffness matrix. By minimizing the displacement error between two successive iterations and utilizing Gerschgorin's circle theory, the entries of the fictitious mass matrix can be obtained as follows [⁵⁴]:

(7)m_{ii}^n=\max \left[{\displaystyle \frac{{(t^n)}^{\mathrm{2}}}{\mathrm{2}}}{k}_{ii},{\displaystyle \frac{{(t^n)}^{\mathrm{2}}}{\mathrm{4}}}{\displaystyle \sum _{j=\mathrm{1}}^{ndof}\left|k_{ij}\right|}\right]\mathrm{.}

In general, the minimum frequency of the artificial dynamic system is achieved with the help of Rayleigh’s principle. The minimum frequency can be calculated by using the succeeding formula:

(8){\omega }_{\mathrm{m}\mathrm{i}\mathrm{n}}=\sqrt{{\displaystyle \frac{X^{T}F}{X^{T}MX}}}\mathrm{.}

In this equality, {\omega }_{\min } is the minimum frequency. Zhang assumed that c=\mathrm{2}{\omega }_{\min } [⁶⁶]. Based on structural dynamic theories, Rezaiee-Pajand and Alamatian [⁵⁴] presented the coming relations for the determination of viscous damping:

(9)C=\sqrt{{\omega }_{min}^{\mathrm{2}}[\mathrm{4}-t^{\mathrm{2}}{\omega }_{min}^{\mathrm{2}}]}M\mathrm{.}

Usually, it is presumed that the time step is equal to one. Besides this, several algorithms were proposed to estimate it. One of the famous procedures suggested for detecting the time step is based on minimization of the residual forces, which leads to the coming result [⁵¹]:

(10)t={\displaystyle \frac{R^T\stackrel{.}{F}}{{\stackrel{.}{F}}^T\stackrel{.}{F}}}

In this relation, the internal force increment of the ith degree of freedom is denoted by \stackrel{.}{F}. Based on minimization of the residual energy, other researchers presented a formula to calculate the time step [⁵⁷].

Tracing structural static path

Nonlinear behavior of structures includes various properties. Among them, buckling loads, post-buckling regions, load and displacement limit point in the stable and unstable path can be mentioned. In addition, these points have a key role in post-buckling behavior of structures. In the stable path of structures, the load and displacement increase simultaneously. On the other hand, unstable path of the structure may occur when the force decreases and displacement increases. Existence of the limit point leads to the difficulties in tracing the equilibrium path. These issues indicate that choosing an appropriate method for evaluating the nonlinear behavior of structures is considerably significant.

Note that load and displacement limit points cause snap-through and snap-back in the equilibrium path. For tracing the equilibrium path of structures, DRM can be deployed. In most of the DR formulations, it is presumed that the external loads are constant in each load increment. As a result, the common DR procedures cannot trace the aforesaid regions. Around these limit points, the structural static path mostly jumps to another location.

In Fig. 1, snap-through and snap-back parts of the equilibrium path are demonstrated. The displacement and load limit points are shown by A and C, respectively. By jumping, the structural static path has been transferred from A to B with slight increment in displacement. Furthermore, increasing slightly the load leads to the movement of point C to D. In other words, the common DRM can only trace the increasing branch of the equilibrium path. To remove this limitation, automatic techniques were proposed by researchers. In these strategies, variable external forces are considered to trace snap-through and snap-back parts.

Rushton employed DRM for post-buckling of the plates [⁶⁷]. The post-buckling behavior of laminated composite plates with large deformation was assessed by Turvey and Wittrick [⁶⁸]. Hook and Rushton [⁶⁹] investigated the buckling of beams and plates with intermediate support. Kadkhodayan et al. [⁷⁰] analyzed the buckling and post-buckling of plates with the help of DR algorithm and dynamic criterion of the stability. On the other hand, by applying the displacement increment in all iterations, Ramesh and Krishnamoorthy [⁷¹] traced the structural static path, which includes snap-through regions. This process was not efficient, when the snap-back occurs. To remedy this difficulty, these researchers mixed the variable arc length method with DRM [⁷²]. It is worth emphasizing; the arc length scheme is not automatic. Additionally, the efficiency of the arc length approach is greatly depended on the selected reference degree of freedom and the quantity of arc length. Lee et al. [²²,²³] mixed the kinetic DRM with the arc length strategy. In tracing the equilibrium path, which has snap-through and snap-back regions, Rezaiee-Pajand and Alamatian [²⁴] deployed the minimization of the unbalanced force and displacement for calculating the load factor in DR iteration. By this way, the analysis converges to the near static points instead of the far ones. In the minimum residual force technique, the load factor is obtained with the help of coming equation:

(11)\lambda ={\displaystyle \frac{P_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{T}}F}{P_{\mathrm{r}\mathrm{e}\mathrm{f}}^T{P}_{\mathrm{r}\mathrm{e}\mathrm{f}}}}\mathrm{.}

Herein, l and P_ref are the load factor and the reference external load, respectively. Another study proposed a formula by minimizing the out-of-balance energy [²⁵]. This relation has the succeeding appearance:

(12){\lambda }^n={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\left[\mathrm{4}{t}^n{f}_i^n-(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{4}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}\frac{{(p_{i,\mathrm{r}\mathrm{e}\mathrm{f}})}^{\mathrm{2}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}}}\mathrm{.}

By utilizing displacement techniques, Alamatian also presented two relations for computing the load factor [²¹]. The first one is based on the minimization of the unbalanced displacement. This equation has the following shape:

(13){\lambda }^n={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{{(\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n)}^{\mathrm{2}}}}\left[\mathrm{2}{t}^n{f}_i^n-(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\left({\displaystyle \frac{p_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\right)}^{\mathrm{2}}}}}\mathrm{.}

The second criterion is established by minimizing the kinetic energy in each step:

(14){\lambda }^n={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{p}_{i,\mathrm{r}\mathrm{e}\mathrm{f}}{\left({\displaystyle \frac{m_{ii}^n}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\right)}^{\mathrm{2}}\left[\mathrm{2}{t}^n{f}_i^n-(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\left({\displaystyle \frac{m_{ii}^n{p}_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\right)}^{\mathrm{2}}}}}\mathrm{.}

Rezaiee-Pajand and Estiri proposed the following formulas for finding the load factor [⁵⁹–⁶¹]:

(15){\lambda }^n={\lambda }^{n-\mathrm{1}}+{\displaystyle \frac{-{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\left[\mathrm{2}{t}^n{r}_i^{n-\mathrm{1}}+(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}\frac{{(p_i)}^{\mathrm{2}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}}},

(16){\lambda }^n={\lambda }^{n-\mathrm{1}}+{\displaystyle \frac{-{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_{i,\mathrm{r}\mathrm{e}\mathrm{f}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\left[\mathrm{2}{t}^n{r}_i^{n-\mathrm{1}}+(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{4}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}\frac{{(p_i)}^{\mathrm{2}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}}},

(17){\lambda }^n={\lambda }^{n-\mathrm{1}}+{\displaystyle \frac{-{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_i}{{(\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n)}^{\mathrm{2}}}}.\left((\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{2}{t}^n{r}_i^{n-\mathrm{1}}\right)}}{\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\left({\displaystyle \frac{p_i}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\right)}^{\mathrm{2}}}}}\mathrm{.}

It is worth mentioning, for tracing the structural static path; the residual force is obtained from Eq. ‎(18), instead of Eq. ‎(3):

(18)R=\lambda P-F\mathrm{.}

In this paper, a new method for tracing the static equilibrium path of structures is suggested. This approach is based on external work.

New load factor

To find better nonlinear equation solver, researchers have combined DR with the other methods, and published them in the engineering journals, so far. For example, combination of the methods, such as, DR and residual load minimization [²⁵], DR and arc-length [²²,²³], DR and residual displacement minimization [²¹] have been done up to now. In the same category but in the different way, the authors will propose another new formulation. Based on this technique and in each iteration, the external work increment is equal to zero. The work is increment calculated by using the following formula:

(19)\delta w^n={\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}\delta {\lambda }^n{P}_i\delta X_i}\mathrm{.}

Substituting displacement increment obtained from Eq. ‎(5) and deploying Eq. ‎(18) into Eq. ‎(19) lead to the next result:

(20)\begin{array}{l}\delta w^n={\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}\delta {\lambda }^n{p}_i{t}^{n+\mathrm{1}}{\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}}=\delta {\lambda }^n{t}^{n+\mathrm{1}}{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{p}_i\left[{\displaystyle \frac{\mathrm{2}{m}_{ii}^n-C_{ii}^n{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}{\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{2}{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}{r}_i^n\right]}\\ \begin{array}{l}\end{array}\begin{array}{l}\end{array}\begin{array}{l}\end{array}=\delta {\lambda }^n{t}^{n+\mathrm{1}}{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{p}_i\left[{\displaystyle \frac{\mathrm{2}{m}_{ii}^n-C_{ii}^n{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}{\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{2}{t}^n}{\mathrm{2}{m}_{ii}^n+C_{ii}^n{t}^n}}({\lambda }^n{p}_i-f_i^n)\right]}\mathrm{.}\end{array}

As previously stated, the work increment is zero in this approach. Hence, by setting Eq. ‎(20) zero, two values are achieved for load factor. One of them is\delta {\lambda }^n=\mathrm{0}. The second one is shown in Eq. ‎(21). Since the equilibrium path is dependent on {\lambda }^n, the zero value is not acceptable.

(21)A^n{\lambda }^n+B^n=\mathrm{0}\begin{array}{ll}& \end{array}\Rightarrow \begin{array}{ll}& \end{array}{\lambda }^n=-{\displaystyle \frac{B^n}{A^n}},

whereA^n and B^n are calculated from the coming relationships:

(22)A^n=\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}\frac{{(p_i)}^{\mathrm{2}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\begin{array}{l}\end{array}\begin{array}{l}\end{array},\begin{array}{l}\end{array}\begin{array}{l}\end{array}{B}^n={\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{-p_i}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\left[\mathrm{2}{t}^n{f}_i^n-(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}\mathrm{.}

Inserting the aforementioned values into Eq. ‎(21) leads to the new load factor:

(23){\lambda }^n={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^{ndof}{\displaystyle \frac{p_i}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}\left[\mathrm{2}{t}^n{f}_i^n-(\mathrm{2}{m}_{ii}^n-t^n{C}_{ii}^n){\stackrel{.}{X}}_i^{n-\frac{\mathrm{1}}{\mathrm{2}}}\right]}}{\mathrm{2}{t}^n{\displaystyle \sum _{i=\mathrm{1}}^{ndof}\frac{{(p_i)}^{\mathrm{2}}}{\mathrm{2}{m}_{ii}^n+t^n{C}_{ii}^n}}}}\mathrm{.}

This load factor is based on the external work technique, which will be utilized in the geometric nonlinear analysis. It is worth emphasizing; mixing dynamic relaxation with external work technique was not so far used for tracing static equilibrium path.

Including the formula, the DRM steps for tracing the equilibrium path are as follows:

Step 1: Select initial values for velocity (zero vector), displacement (zero vector or convergence displacement of the previous increment), fictitious time step (t=\mathrm{1}), and convergence criterion for the residual force (e_R={\mathrm{10}}^{-\mathrm{4}}) and the kinetic energy (e_K={\mathrm{10}}^{-\mathrm{4}}).

Step 2: Calculate the elemental tangent stiffness matrix and internal force vector.

Step 3: Construct fictitious mass and damping matrices.

Step 4: Find the load factor from Eqs. ‎(11)–(14) and Eqs. ‎(23). Note that these relations are used for minimizing residual force, out-of-balance energy, unbalanced displacement, kinetic energy and external work, respectively.

Step 5: Calculate the residual force vector by using Eq. ‎(18).

Step 6: If \sqrt{{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{(r_i^n)}^{\mathrm{2}}}/{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{(r_i^{\mathrm{1}})}^{\mathrm{2}}}}\le e_R, go to Step 9; Step otherwise, update the fictitious velocity vector using Eq. ‎(4).

Step 7: If \sqrt{{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{\left({\stackrel{.}{X}}_i^{n+\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}}/{\displaystyle \underset{i=\mathrm{1}}{\overset{ndof}{\mathrm{\Sigma }}}{\left({\stackrel{.}{X}}_i^{\mathrm{1}}\right)}^{\mathrm{2}}}}\le e_K, go to Step 9; otherwise, continue.

Step 8: By employing Eq. ‎(5), update the displacement vector and go to Step 2.

Step 9: Print the displacements and load factor of the current increment.

Step 10: If \lambda >{\lambda }_{\max }, stop; otherwise, \lambda =\lambda +\delta \lambda and go to Step 2.

Where, \delta \lambda is the load factor increment, and it is selected by the user. On the other hand, the maximum load factor is{\lambda }_{\max }. In all analyses, {\lambda }_{\max }is assumed to be 10. Moreover, r¹ and Ẋ¹ are the residual force and velocity at the first iteration each step loading, respectively.

To obtain structural load, the load factor is multiplied to the reference load. By utilizing this load, the structure will be analyzed, and the responses are found. The next increment will be started with these results. Both predictor and corrector steps are required in the common work control formulations. Hence, the load increments are different for the first iteration and corrector ones. In the general incremental-iterative methods, load factor increment sign should be specified. In other words, it is important that the increment will increase or decrease the load. It is worth emphasizing; the sign of load factor increment changes at the load limit points [⁷³]. Therefore, the sign determination will be needed for the common work control method. These disadvantages are eliminated in the proposed algorithm. In fact, the suggested load factor depends only on the DR fictitious parameters. Accordingly, all calculations are performed via vector operation. Moreover, the load factor is calculated only by one formula. Thus, it has the same relation in the all process. In contrast to the arc length techniques, which required the sign determined, the proposed scheme does not need any sign finding. This is another advantage of the new formulation.

The linear equilibrium equations for a flat triangular shell element in its local coordinate system are first perturbed to yield the in-plane geometric stiffness matrix. Then out-of-plane considerations that involve the effect of rigid body rotations on member forces yield an out-of-plane geometric stiffness matrix. The shell element that was chosen for that purpose combined the constant stress triangle (CST) flat triangular membrane element and of the discrete Kirchhoff theory (DKT) flat triangular plate element. Finally, a computer program, featuring incremental analysis and DR method, geometric effects, pure deformation isolation, internal stress retrieval and updating of nodal forces and coordinates is presented and used to solve a number of problems. Authors modified and used the reference [⁷⁴] for analyzing the shell structure. Moreover, this code can be adjusted by iso-geometric analysis (IGA), which is an appropriate method for modeling the Kirchhoff-Love shell model [⁷⁵].

Numerical samples

To utilize the presented method, a FORTRAN program is prepared by the authors. By using of this program, the geometric analyses of various structures are performed. Among them are 2D and 3D trusses, planer arcs, frames and shells. All data unites are consistence. The proposed scheme is independent of element types. Hence, this technique can be applied to any other finite elements too. For this purpose, in addition to the truss and frame elements, a four-node quadrilateral element is used for geometric nonlinear analysis of the bending plates. In all of problems, the reference degree of freedom is shown by D, and the related structural static path is drawn. In the previous sections, the formulas required for calculating the load factor of DR tactic’s iterations were presented. Table 1 demonstrates all symbols for the used schemes.

To analyze the structures, it is presumed that the initial displacements are zero. For other increments, the achieved static point in the previous step is utilized for the initial guess. Moreover, the fictitious time step is equal to one in all iterations. It is worth emphasizing; the load factor increment affects the number of points of the equilibrium path. In other words, if the increment decreases, the load-displacement curve is traced with more points. With the help of this property, the structural behavior can be assessed more accurately. For comparing the efficiency of the methods, three criteria are proposed. In the following line, these criteria or indicators are introduced:

(24)S\mathrm{1}={\displaystyle \frac{NumberofConvergedPoints}{TotalIterations}}\begin{array}{l}\end{array},\begin{array}{l}\end{array}S\mathrm{2}={\displaystyle \frac{TotalIterations}{TimeIterations}}\begin{array}{l}\end{array},\begin{array}{l}\end{array}S\mathrm{3}={\displaystyle \frac{NumberofConvergedPoints}{TotalIterations}}\mathrm{.}

The greater value for S1 indicates that the method has obtained more convergence points in lesser iteration. S2 and S3 are directly related with the consuming time of the algorithms. On the other hand, the accuracy of the obtained static equilibrium path can be improved by increasing the number of convergence points. Additionally, the consuming time of the procedure is boosted by increasing the convergence points. In addition to improve the accuracy and efficiency of DR tactic, this paper aims to reduce the consuming time of the aforementioned technique. The above-cited criteria are applied to compare the convergence rate of the presented approaches.

Two-member truss 1

A planar truss which has two members is demonstrated in Fig. 2. This structure is under a concentrated downward load. Previously, Greco et al. [⁷³] analyzed this truss. He used two different relations; positional and corotational formulations. Their results are used for verification of authors’ method and the related program. The elasticity modulus, the reference force, the cross-sectional areas of member and the increment of the load factor are 7.17 × 10⁴ N/mm², 100 kN, 60 mm², and 0.1, respectively. The analytical solution for this problem is available. It was obtained by imposing the equilibrium at the deformed position. Equations ‎(25) and ‎(26) are nonlinear equilibrium equation based on the applied vertical force and deformed position [⁷⁶]. The subscript numbers are denoted bars' number. Other parameters are shown in Fig. 2.

(25)-{\displaystyle \frac{\cos {\gamma }_i}{\cos {\beta }_i}}{\displaystyle \frac{P}{\cos {\gamma }_i.\tan {\beta }_i+\sin {\gamma }_i}}=E_{\mathrm{1}}{A}_{\mathrm{1}}\left[\mathrm{2}{\displaystyle \frac{L-x}{L}}{\displaystyle \frac{\cos {\beta }_{\mathrm{0}}}{\cos {\beta }_i}}-\mathrm{1}\right]\mathrm{.}

(26)-{\displaystyle \frac{P}{\cos {\gamma }_i.\tan {\beta }_i+\sin {\gamma }_i}}=E_{\mathrm{2}}{A}_{\mathrm{2}}\left[\mathrm{2}{\displaystyle \frac{x}{L}}{\displaystyle \frac{\cos {\beta }_{\mathrm{0}}}{\cos {\gamma }_i}}-\mathrm{1}\right]\mathrm{.}

Moreover, Young’s modulus and cross section are shown by E and A, correspondingly. This truss has only one degree of freedom, if geometry and loading are symmetric. In other words, the tip horizontal displacement equal to zero. Fellippa [⁷⁷] derived the following tangent stiffness matrix for this structure. Equation ‎(27) shows this matrix. Moreover, he calculated the vertical force for this truss by using Eq. ‎(28), too. In these formulas, vertical motion is positive upwards.

(27)K={\displaystyle \frac{\mathrm{8}EA}{\sqrt{{(\mathrm{4}{H}^{\mathrm{2}}+L^{\mathrm{2}})}^{\mathrm{3}}}}}[\mathrm{3}{D}_y^{\mathrm{2}}+\mathrm{6}H{D}_y+\mathrm{2}{H}^{\mathrm{2}}],

(28)\lambda P={\displaystyle \frac{\mathrm{4}EA}{\sqrt{{(\mathrm{4}{H}^{\mathrm{2}}+L^{\mathrm{2}})}^{\mathrm{3}}}}}[\mathrm{2}(H+D_y)(\mathrm{2}H{D}_y+D_y^{\mathrm{2}})]\mathrm{.}

where, D_y, H and L are the tip deflection, the height and span length of truss, respectively. By setting stiffness matrix determinant zero, two roots are reached. The following roots are load limit points, which are always real value.

(29)\mathrm{3}{D}_y^{\mathrm{2}}+\mathrm{6}H{D}_y+\mathrm{2}{H}^{\mathrm{2}}=\mathrm{0}\begin{array}{ll}& \end{array}\to \begin{array}{ll}& \end{array}{D}_y=-\mathrm{1}\pm {\displaystyle \frac{H}{\sqrt{\mathrm{3}}}}\mathrm{.}

By substituting these values into Eqs. ‎(28) and ‎(29), the exact displacements and loads of the limit points are achieved. The displacements are 0.4226 m and 1.5774 m. Moreover; the limit point loads are±148.102 kN, correspondingly. Now, this truss will be analyzed numerically, too. The proposed method is utilized for analyzing two-member truss. Figurer 3 demonstrates the structural response for various load levels. In this figure, the equilibrium path of the center node is plotted. The obtained results are compatible with the outcomes of the other researchers [⁷³,⁷⁷]. Note that all the schemes lead to the same load-displacement curve. At the first limit point, the load factor is 1.478. This value is multiplied by the reference load of 100 kN to be equal to 147.8 kN, which is the maximum force that the structure can support before snap-through buckling. Afterward, the load decreases until another limit point is attained. This value is -1.469. Following this, the load and displacement increases simultaneously. It is worth emphasizing; the limit points are similar to exact solution.

There are two zero-load configurations for this truss. At the first state, the top node deflection is 1 m. At this case, the structure is placed in a horizontal position. At the second, it is 2 m. In other words, the configuration is symmetric to the initial situation with respect to the horizontal plane. The initial situation is when the truss is unloaded. Moreover, Fig. 3 shows that the accuracy of MRF, MDI, MKE, and WCM are the same. In other words, these procedures can trace the structural static path with the same converged points and number of total iterations.

Table 2 shows the values of analysis indicators. Based on S1, MRF technique converges more rapidly than the other strategies. WCM, MKE, and MDI have the next ranking. The slowest scheme is MRE. Based on S2 and S3, it is obvious that WCM is the most appropriate one. In other words, in this tactic, the consumed time for iterations is less than the other approaches. On the other hand, for all the algorithms, the locations of the convergence points in the snap-through part (between 0.5 and 1.5 m), are approximately similar. By investigating these indicators, it can be deduced that WCM is more efficient, in comparison with other techniques.

Two-member truss 2

Here, the truss shown in Fig. 4 is analyzed. The axial rigidity and the reference load are 8366.888 N and 270 N, respectively. This structure is analyzed under two load pattern [⁷⁸]. In first case, horizontal load is treated as an imperfection, i.e., P_u = 0.05P_v. In second case, vertical load considered as an imperfection, i.e., P_v = 0.05P_u.

The equilibrium paths related to the vertical and horizontal displacements of the tip node under the first loading case are illustrated in Fig. 5. Based on these figures, four load limit points emerge. At the first point, the relevant load, the vertical and horizontal displacements are 1863.356 N, 183.41 and 145.67 mm, respectively. This is the maximum load of structure before it becomes unstable. Other load limit points are occurred in±2264.904 and -1863.356 N. For the horizontal displacement, the snap-back points are placed at displacement±463.55 mm. This truss has two zero-load situations. These positions are placed in the snap-back points. In these cases, the vertical deflection of the tip node is 656.59 mm.

To compare these techniques, number of increments and iterations of each method under the first load pattern are presented in Table 3. Based on Table 3, MRF and WCM are the most appropriate methods.

Figure 6 shows the load-displacement curves of the tip node under the second loading pattern. At the first stage, the obtained results are compared with the responses of other references [⁷⁸]. Both the achieved results and the responses of the reference are similar. It is worthwhile to remark that the force-displacement curves of all the methods are the same. However, their number of convergence points and analysis duration are not similar.

Based on Fig. 6, it is obvious that the post-buckling behavior begins after the load value reaches 2410.487 N. Hence, this force is the load limit point. The corresponding displacements of this load are 419.84 and 592.4 mm. The first snap-back point is occurred when the tip deflection is 661.92 mm. The next limit displacement points are placed in 651 and 1312.93 mm. In the last point, the load is zero. The second snap-through appears at the load of -2410.487 N. The vertical and horizontal displacements are 893.06 and -592.46 mm, respectively. It is worth emphasizing; although the WCM strategy has fewer converged points in comparison with MDI and MKE, its accuracy is acceptable. Note that the converged increments are very close in MDI and MKE procedures. In authors' approach, the number of convergence increments increase near the limit points. On the stable branch of the equilibrium path, the number of these points is less than that of MDI and MKE schemes.

In Table 4, the number of iterations and the score of the methods are shown. The numbers presented in a parenthesis denote the grade of each scheme. Since the number of convergence points of MDI and MKE methods are more than the proposed strategy, their analysis duration is more than the other schemes. Based on Table 4, the most suitable procedures for the analysis of this truss are MRF and MRE processes. Furthermore, the suggested algorithm is ranked third.

Yang and Shieh [⁷⁸] achieved an improved method for nonlinear problems with multiple critical points. It is worth emphasizing; the number of iterations and convergence points of the arc length algorithms are dependent on the value of the arc length. In other words, the analyzer should set the initial arc length. Therefore, a slight change in the magnitude of the arc length may lead to divergence. As a consequence, the wrong equilibrium path may be achieved. On the other hand, analyzer's role in the dynamic relation technique is considerably less than the aforesaid method. In the DR strategies, as the solution procedure approaches the limit points, the iteration number increases usually. As a result, this algorithm can easily trace the structural equilibrium path.

Sixteen-member shallow truss

The next sample shows the geometric nonlinear behavior of a space truss. This structure has 16 members and 9 nodes. Four corner nodes of this truss are simply supported. The geometric characteristics of this structure are demonstrated in Fig. 7. Hrinda [⁷⁹] analyzed this truss by using an arc length approach. The height of the middle nodes and the top node are 20 and 40 mm, respectively. Member’s cross-sectional areas are 100 mm². Young’s modulus of the members is 10⁷ N/mm². Additionally, the reference load is equal to 25 kN. To analyze this truss; it is assumed that the load factor increment of each load step equals 1.

Hrinda reached the limit points at the load factor of±3.84. Moreover, the displacements of these points were 0.97 and 3.09 cm, respectively. The equilibrium path of the reference degree of freedom is illustrated in Fig. 8. Based on this curve, it is obvious that the convergence points of the MDI, MKE, and WCM techniques are the same. In other words, they have the same accuracy. These three methods have traced more appropriately the structural static path between two load limit points. In the first limit point, the displacement and the load are 9.506 mm and 94.91 kN, respectively. This is the maximum load that the structure can support while still remaining stable. Then, the load decreases and deflection increases until the second limit point occurs. The next limit point is placed in the snap-through part of the equilibrium path. The load and the displacement of second point are -95.38 kN and 30.64 mm, respectively. Note that, these points are analogous to the ones mentioned in Ref [⁸⁰]. Results show that the proposed method has a suitable accuracy. It is worth emphasizing; this structure is unstable as long as the displacement is equal to 4.45 cm. Then, the structure becomes stable. The numbers of converged points of the MRF and MRE procedures between two limit points are 14 and 19, correspondingly. It should be added, MDI, MKE and WCM include 24 convergence points. In other words, these methods can trace more accurately the structural static path. Although the accuracies of these schemes are the same, but the consumed time for WCM is lesser than the others.

According to Fig. 8, if the vertical displacement of tip node is 20 or 40 cm, then the corresponding loads are zero. It is noted that in this state, nodes' deflections located in the middle plane are equal to zero. It should be reminded that the top node is located 20 mm above the middle horizontal plane. Hence, at the first zero-load configuration (when tip deflection is 20 mm), top node is coplanar with middle nodes. Because it is located 20 mm below middle nodes' plane in the second zero-load position, this state is symmetric to the initial unloaded situation with respect to this plane.

In Table 5, the results of this structure are inserted. Based on indicators S1, S2, and S3, the efficiency ranks of three methods are as follows, MRF, WCM and MRF. The indicator S1 of the MDI, MKE and WCM, are the same. Based on S2, authors’ scheme is the most suitable strategy. By utilizing the other criteria, the suggested algorithm is ranked third. Nevertheless, this scheme traces the static equilibrium path more appropriately, in comparison with MRF and MRE. Therefore, based on the criterion S3, WCM has the first rank.

Cross-shaped truss

In this part, the truss shown in Fig. 9 is analyzed [⁵⁷]. This three dimensional structure has three degrees of freedom. The modulus of elasticity, the cross-sectional areas and the reference load are 10⁴ N/mm², 10 mm², and 10 kN, respectively. To analyze this truss, it is presumed that the load factor increment is equal to 0.1. This structure is analyzed by the large deformation theory.

The force and tangent stiffness relationships were obtained as follows [⁵⁷]:

(30)K={\displaystyle \frac{\mathrm{2}EA}{L^{\mathrm{3}}}}[\mathrm{2}{(D_y+H)}^{\mathrm{2}}+(D_y^{\mathrm{2}}+\mathrm{2}H{D}_y)],

(31)\lambda P={\displaystyle \frac{\mathrm{2}EA}{L^{\mathrm{3}}}}[(D_y+H)(D_y^{\mathrm{2}}+\mathrm{2}H{D}_y)],

where, D_y, H and Lare the tip deflection, the height of truss and initial length of each member, respectively. By setting the stiffness matrix zero, the critical displacement can be calculated from the coming relationship:

(32)\mathrm{3}{D}_y^{\mathrm{2}}+\mathrm{6}H{D}_y+\mathrm{2}{H}^{\mathrm{2}}=\mathrm{0}\begin{array}{ll}& \end{array}\to \begin{array}{ll}& \end{array}{D}_y=-\mathrm{1}\pm {\displaystyle \frac{H}{\sqrt{\mathrm{3}}}}

Based on Eqs.(27) and (30), the displacements of limit points are the same for these trusses. However, the loads are not similar. By substituting the values into Eqs. ‎(31) and ‎(32), the exact solutions are obtained. Hence, the displacements are 65.13 and 243.09 mm. Moreover, the limit point loads are±45.442 kN, correspondingly. Following this sample, the structure is analyzed numerically, and the results are achieved.

All four presented methods are used to determine the load factor-displacement curve. The obtained structural static paths are shown in Fig. 10. It should be reminded that MRE technique diverges at the 24th increment. As a result, this method is not able to trace the static equilibrium path. Based on the finding curve, it is obvious that MRF, MDI, MKE, and WCM converged to the same point. This analogy is rooted in the fact that the load factor increment is assumed to be a small value. However, the analysis duration is different. The proposed scheme can trace the structural equilibrium path faster related to the other strategies. Two snap-through points occur in the equilibrium path of this truss. Limit points' displacements are 58.95 and 234.4 mm, respectively. Moreover, the limit point loads are 45.11 and 44.82 kN, correspondingly. These results are similar to the exact solution. Similar to example 1, there are two zero-load configurations for this problem. At the first, the tip deflection is 154.15 mm. This value is the height of truss. Therefore, the top node is located in a horizontal location. At the second state, it top point deflection is 308 mm. Hence, the configuration is symmetric to the initial unloaded position with respect to the horizontal plane.

Table 6 includes the indicator values and the ranking. The WCM consumes less time than the other methods. Consequently, it is ranked first, based on the criterion S2 and S3. On the other hand, MKE and MDI techniques lead to the similar results. In other words, their consuming time, number of iterations and number of convergence points are the same. By comparing these schemes, it can be deduced that WCM is the best strategy for analyzing this truss.

Dome truss

Yang et al. [⁸⁰] previously analyzed this truss. They considered the material of structure elastic or elasto-plastic. Moreover, other researchers solved it with geometric nonlinear behavior, too [²⁵]. This structure includes 73 nodes, 168 members, and 147 degrees of freedom. The dome truss is shown in Fig. 11. The height of the truss is 1790.22 mm. Based on the pervious solution, this truss behaves extremely nonlinearly. The elasticity modulus and members’ cross-sectional areas are 10⁵ N/mm² and 1 mm², respectively. Furthermore, the reference load is 100 N. it should be added that this load is applied to the truss tip. Additionally, it is assumed that the load factor increment is 1. The lowest level of this structure is simple supported.

The answers belong to the presented techniques are inserted in Table 7. It should be added that the number of the convergence points of the MRF method is less than the other tactics. As a result, this approach is the least efficient scheme. Moreover, fewer numbers of the convergence points in MRF and MRE algorithms lead to reduce the analysis time. Based on the criterion S2, the WCM is the most appropriate tactic. The MRF is not practical in snap-through parts of the static equilibrium path. The local jumps have been occurred in the MRF snap-through part, and therefore, the WCM is ranked first based on indicator S3. In other words, WCM is the most efficient approach for tracing the static equilibrium path of this truss. In addition to its high accuracy, its convergence rate is appropriate.

In Fig. 12, the static path of the truss is plotted. In this curve, the vertical displacement of the top node is shown. First, the load and displacement increase simultaneously. The process is continued until the deflection is equal to 74.53 mm. The corresponding load is 849.6 N. This force is the first snap-through point. Afterward, the structure is entered to the unstable region, in which the load is backward and displacement increases. The latter limit point occurs in the load and deflection of -664.7 N and 324.9 mm, respectively. Again, the load and displacement increase simultaneously. It should be reminded; the structure is unstable between 75.53 and 502 mm. The second point is where the load factor is equal to the first limit point.

Note that all methods are able to trace the equilibrium path of this truss. However, the numbers of convergence points and their locations are different. Herein, MDI, MKE, and WCM perform appropriately and trace the complete structural static path. On the other hand, MRF and MRE algorithms are ranked second and third, respectively. The numbers of the convergence points of MRF and MRE scheme between two limit points are 40 and 98, correspondingly. This value is 155 for all three other methods. The dome truss static path in Fig. 12 has complete agreement with the previously published results [²⁵].

Star truss under symmetric loading

Figure 13 shows the star truss. This structure has been previously analyzed by other researchers [²⁵]. Some investigators solved this truss with various loading [⁸¹,⁸²]. At first; this shallow shell was studied by Hangai and Kawamata [⁸³]. This three-dimensional truss has 13 nodes and 21 degrees of freedoms. The modulus of elasticity and the cross-sectional areas are 303000 N/mm² and 3.17 mm², respectively. Furthermore, the load factor increment is assumed to be 1. The reference load equals 50 N for node 1. Moreover, the reference loads are 25 N for nodes 2, 3, 4, 5, 6, and 7.

The structural static path related to the vertical displacement of the nodes 1 and 2 are shown in Fig. 14. In the path associated with node 1, two limit load points exist. The load factor at the first point is 7.967. Additionally, the relevant displacement is 7.841 mm. At this limit point, the structure becomes unstable. Afterwards, the load increment decreases until another limit point is reached. The corresponding load factor and displacement are -4.508 and 29.06 mm, correspondingly. Note that the truss is still unstable. At this state, the load reverses and begins to increase. The structure continues to deform until a new stable equilibrium point is obtained, where the load factor is equal to the first limit load point, and the structure becomes stable. Hence, it can sustain an additional load. The star truss has two zero-load situations. At the first case, the vertical deflections of nodes 1 and 2 are 18.8 mm and -1.15 mm, respectively. Therefore, the node 1 was placed in the middle horizontal plane. At the latter state, the displacements are 40 mm and zero mm, correspondingly. At the second configuration, it was 20 mm below the internal plane. Hence, the second situation is symmetric to the initial unloaded position with respect to this plane.

Note that MRF and MKE are not able to trace the equilibrium path of this truss. It is worth emphasizing; MRE, MDI, and WCM are the most suitable techniques. These approaches are capable of completely tracing the static path of the star truss. However, their numbers of convergence points are not similar. MDI, and WCM schemes require more increments to trace the static equilibrium path, in comparison to the MRE. In MRE, MDI and WCM algorithms, the numbers of convergence points between the load limit points are 27, 53 and 38, respectively. It should be added that these three methods can trace the equilibrium path. The best three methods are MDI, WCM and MRE.

To compare these techniques, the number of the increments, iterations and the analysis duration of the aforementioned approaches are inserted in Table 8. In Table 9, the ranking of these strategies are presented. It is obvious that MRE scheme is the most appropriate method. However, in tracing the path between two limit points, this approach is less accurate than WCM and MDI algorithms.

Star truss under asymmetric load

This truss is demonstrated in Fig. 13. The aforesaid structure is subjected to the asymmetric loading. This truss was previously analyzed by other researchers [²¹]. Young's modulus, the cross-sectional areas and the load factor increment are 5000 N/mm², 2 mm², and 1, respectively. The reference loads related to nodes 1, 2, and 3 are 0.4 N, correspondingly.

The vertical and horizontal displacements of node 1 are plotted versus load factor in Fig. 15. In the equilibrium path relevant to vertical displacement of node 1, two load limit points are existed. The corresponding loads are 3.87 N and -2.467 N, respectively. Additionally, the displacements of these points are 7.777 and 29.49 mm, respectively. It should be reminded that the structure is unstable between 7.777 and 47.04 mm. In other words, the first point is the initial snap-through point. Moreover, the second point is where the load factor is equal to the first limit point. The present load-displacement curves are in good agreement with the previously results [²¹]. At loads 3.827 N and 2.45 N, the snap-back points emerge in the horizontal displacement of the node 1 static path. Note that the related displacements are 0.2767 mm and -0.2446 mm, correspondingly. Both MRF and MKE are not able to trace the equilibrium path of the aforementioned truss. Conversely, MRE, MDI and WCM are the most appropriate tactics. These algorithms are able to trace the equilibrium path completely. Nevertheless, their numbers of convergence points are different. MDI, and WCM require more increments to trace the structural static path. It should be added that, in MRE, MDE, and WCM methods, the numbers of convergence points, between load limit points, are 24, 55, and 35, respectively. For the aforesaid tactics, the numbers of convergence points, between displacement limit points, are 23, 56, and 40, correspondingly.

In Table 10, the number of increments, iterations and consumed time of each algorithm are inserted. Based on the previously introduced criteria, the schemes are ranked in Table 11. According to the obtained results, it is concluded that MRE strategy is the most suitable technique. However, this approach is the least accurate scheme. Note that authors’ algorithm is ranked first, based on the criterion S1. It is worth emphasizing that the given conclusions are compatible with the results of the previous sample.

Schwedler dome

Schwedler's dome is a 240-member truss which has 219 degrees of freedom. This structure is shown in Fig. 16. Krishnamoorthy et al. [⁸²] analyzed this truss. Young's modulus and members’ cross-sectional areas are 2.1 × 10⁵ N/mm² and 450 mm², respectively. The reference load is 600 N. It is applied to truss' tip. All based nodes are hinged.

Figure 17 plots the tip load against the vertical and horizontal displacements. The vertical displacement curve of the tip node (node 1) has snap-through points. Furthermore, the horizontal displacement (node 2) has a snap-back point. Clearly, all methods are able to trace the structural static path completely. However, their numbers of convergence points are different. When MRF is utilized, local jumping occurs in the equilibrium path. The equilibrium path obtained by MRF, MRE, MDI, MKE, and WCM technique includes 110, 429, 746, 746, and 746 convergence points between the limit points. The snap-back point related to the horizontal displacement occurs when the displacement reaches 4.336 mm. The results achieved from these schemes are inserted in Table 12.

William Toggle frame

A two-member frame is shown in Fig. 18. This structure was previously analyzed by other researchers [²³,²⁵,⁸⁴]. William [⁸⁵] solved this frame analytically by considering the effects of the axial force and the flexural shortening. He proposed the force-deflection relation as follows;

(33)\lambda P\approx \mathrm{2}(F_{\mathrm{1}}+F_{\mathrm{2}}\sin {\beta }_{\mathrm{0}}),

(34)F_{\mathrm{1}}={\displaystyle \frac{EA}{L}}\left[D_y\sin {\beta }_{\mathrm{0}}-\mathrm{0.6}{\displaystyle \frac{D_y^{\mathrm{2}}}{L}}\right],

where, F₁ and F₂ are depended on the axial and flexural rigidity, correspondingly. Moreover, the initial length, the member angle relative to the horizontal and tip deflection are showen by L, {\beta }_{\mathrm{0}} and D_y, respectively. To solve the Toggle frame, Wood and Zienkiewicz [⁸⁴] used a total lagrangian finite element formulation.

To solve this frame, the finite element method is applied. A nonlinear solution by using large deformations theory is pursued. Each frame member modeled with five elements. In this sample, it is assumed that reference load, elasticity modulus, cross-sectional and moment of inertia are 44.28 N, 7.1 × 10⁴ N/mm², 118.06 mm², and 374.61 mm⁴, respectively. Moreover, the load factor increment is equal to 1.

The tip node displacement of the toggle frame is plotted versus the load factor in Fig. 19. The results are in good agreement with Williams’ curve, which includes flexural shortening. All the presented tactics are able to completely trace the equilibrium path of this structure. However, the numbers of convergence points belong to these tactics are different. In comparison to other schemes, MDI, MKE and WCM methods require more increments. The load and displacement of the first load limit point are 155.4 N and 6.1 mm, respectively. In this diagram, the snap-through part of the static equilibrium path is gradually formed. In other words, after the first limit point, the structure softens slightly with increasing deformations as the load decreases and then the hardening state occurs. For the second point, the related load and displacement are 143.16 N and 10.28 mm, correspondingly. It should be noted that the toggle frame has an unstable region between 6.1 mm and 12.44 mm.

The findings are inserted in Table 13. Note that the numbers of convergence points are similar in MDI, MKE and WCM approaches. Based on this table, obviously, WCM algorithm is the fastest technique. In other words, although the number of increment is similar in three methods, but the convergence rate of the WCM is more than other schemes. Consequently, the suggested technique is ranked first, based on the indicators S2 and S3.

Shallow arc

A two-dimensional arc is shown in Fig. 20. The geometry of this structure is defined by the coming function:

(35)y=\mathrm{5}\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{\mathrm{\pi }x}{\mathrm{1}\mathrm{0}\mathrm{0}}}\right)\begin{array}{l}\end{array}\mathrm{c}\mathrm{m}\mathrm{.}

This arc is subjected to the asymmetric loading. Note that eccentric loads are applied. The eccentricity of the both loads is assumed to be 10 cm relative to the center point. To analyze this structure, 20 frame elements are used. Young’s modulus is equal to 10⁵ N/mm². Moreover, the cross-sectional area and moment of inertia are 32 mm² and 10⁴ mm⁴, respectively. In this problem, it is presumed that the load factor increment is 1. Additionally, loads P₁, P₂ and P₃ are 600, 400, and 400 N, correspondingly.

It should be mentioned; this shallow arc has not been so far analyzed. In Fig. 21, the equilibrium path of this structure is plotted. All the presented techniques can completely trace the equilibrium path of this structure. In other words, the accuracies of all the methods are approximately the same. Nevertheless, the numbers of converged points of these methods are different. For Fig. 21(a), the load is increased until the load factor, and the central deflection become 3.87 and 3.04 cm, respectively. Afterward, the load decreases until another limit point is achieved. This load factor is 2.282. The corresponding deflection is 7.07 cm. Following this state, the load and displacement increases simultaneously. After the first limit point, the arc is unstable until the load factor is equal to 3.87 again. The rankings of these algorithms are presented in Table 14. Based on these outcomes, it is deduced that WCM is the best approach for the analysis of the aforesaid arc.

Hinged circular arch

In this part, a hinged circular arch is analyzed. The geometry of this structure is shown in Fig. 22. Previously, Yang and Shieh [⁷⁸] investigated the behavior of this arch. This arch is divided by 100 straight elements with equal length. The load is applied in a small distance to the tip node. The elasticity modulus, cross-sectional areas and moment of inertia are 1378 kPa, 64.52 cm², and 41.62 cm⁴, respectively. Note that the lowest nodes of this structure are simply supported.

Figure 23 shows the load-displacement curve. At the first stage, the obtained results are compared with the responses of other references [⁷⁸]. Both the achieved results and the responses of the reference are similar. It is worthwhile to remark that the force-displacement curves of all the methods are the same. However, their number of convergence points and analysis duration are not similar. Based on Fig. 23, it is obvious that the post-buckling behavior begins after the load value reaches 25.874016 N. Hence, this force is the load limit point. The other limit points are occurred in -38.083776, 71.168, -99.243776, 170.429568, -231.260416, and 278.231296 N.

In Table 15, the scores of the methods are shown. Based on this table, all indicators S1 for the MDI, MKE and WCM are the same. Based on S2 and S3, authors’ scheme is the most suitable strategy. Based on Table 15, the most suitable procedures for the analysis of this arch are WCM and MKE processes.

Cylindrical roof

The three-node triangular element is used for geometric nonlinear analysis of a shell [⁷⁴]. The structure is a cylindrical roof shown in Fig. 24 [⁷⁵,⁸⁶]. This problem is the part of a cylinder with the angle of 0.1 rad. The modulus of elasticity and Poisson’s ratio of the shell are 3102.75 N/mm² and 0.3, respectively. It is assumed that the structural thickness is 12.7 mm. Furthermore, the reference load is 300 N. To solve this structure; the load factor increment of 1 is utilized. Due to symmetry, a 10×10 mesh is used to discretize a quarter of the roof.

The load-displacement curve for the center node is plotted in Fig. 25. This node has maximum deflection. Moreover, the static path for the point B is drawn. According to the findings, the accuracy of MRF tactic is low. The other methods are approximately the same. The number of the convergence points of MRF scheme between two limit points is 4. This is the lowest value among all techniques. It should be added that the obtained results are compatible with the previously reported outcomes [⁸⁶]. Therefore, the proposed formulation can be used for the shell analysis, as well.

The number of the convergence points, iterations and the analysis duration of the approaches are inserted in Table 16. Based on these results, the MDI & MKE are approximately the same. Moreover, the WCM analysis duration is about %10 less than the MDI & MKE schemes. The grades belong to the presented techniques are inserted in Table 17. It should be added that the number of the converged points for the MRF method is fewer than the other tactics. Furthermore, fewer numbers of the converged points in the MRF algorithm lead to reduction of the analysis time. Based on the criterion S3, the MRF is the most appropriate tactic. Nevertheless, it has the minimum accuracy. On the other hand, the proposed procedure has the first ranking based on S2. The WCM is the efficient approach for finding the load-deflection curve of this shell. In addition to its high accuracy, its convergence rate and time are appropriate. The deformed geometry under the maximum load is portrayed in Fig. 26.

Cantilever plate with end bending moment

In this example, the cantilever plate will be studied, which has a moment at the free end [⁸⁷]. The structural mesh is shown in Fig. 27. This plate has Young’s modulus 1.2 × 10⁶, Poisson’s ratio 0, and the reference load 5.236. It is assumed that the load is applied by a pair of concentrated moments at the end nodes.

Figure 28 depicts the horizontal (U) and vertical displacement (W) curves against the end moment. Based on this figure, the snap-back occurs in the structural static path. This point related to the vertical deflection when the load reaches 8.701 and 19.37. On the other hand, the horizontal displacement curve has the limit point where displacement reaches 14.55. Furthermore, the corresponding load is 37.3. Based on the analysis results, the accuracy of the MRF method is low. On the other hand, the precisions of other approaches are approximately the same. The numbers of the convergence points of MRF, MRE, MDI, MKE, and WCM tactics are 18, 127, 110, 221, and 270, respectively. Figure 29 shows the initial position and two deformed configurations.

In Table 18, the ranks of all schemes are inserted. The Table 18 shows that MRE strategy is the most suitable technique for the cantilever plate. The WCM has the next ranking. Based on the criterion S1 & S3, authors’ algorithm is ranked second. In addition to its high accuracy, its convergence rate is appropriate.

Multi spans cylindrical shell

Now, a two bay cylindrical shell with a 60° central angle will be analyzed. It should be mentioned; this structure has not been so far analyzed. The shell static path has snap-through and snap-back points. The young's modulus, Poisson's ratio and shell thickness are 20 kN/mm², 0.3, and 30 mm, respectively. Additionally, the reference load P is equal to 2.5 kN. The geometry of this structure is shown in Fig. 30. The radius of cylindrical is 1 m. Owing to symmetry; one-quarter of the shell is modeled with 15×10 elements. Note that, four corner nodes (C, D, E and F) are simple support, and these points can’t move at X, Y and Z directions.

The vertical deflection static paths for the points A & B are plotted in Fig. 31(a). In these equilibrium paths, two load limit points are existed. The corresponding loads are 17725 and 11341 N, respectively. Additionally, the deflections are 11.55 and 22.66 cm for point A, correspondingly. These values are 4.954 cm and 8.802 cm for point B, respectively. Figure 31(b) shows the load-displacement curve for the horizontal displacement of point B. This static path has one snap-back point, in addition to two snap-through points. The related load and horizontal displacement are 14451 N and 3.98 mm, respectively. Between two limit load points, the numbers of converged points of the MRF, MRE, MDI, MKE, and WCM procedures are 11, 79, 73, 74, and 48, correspondingly. Nevertheless, all the presented techniques can completely trace the equilibrium path of this structure. In fact, the accuracies of all methods are approximately the same. It should be added, the numbers of converged points of the MRF, MRE, MDI, MKE and WCM approaches, between the first limit load and snap-back points, are 8, 56, 52, 50, and 33, respectively. Figure 32 demonstrates the initial position and the final structural situation. The deformed geometry is portrayed under the maximum load (25 kN).

The number of iterations, convergence points and the analysis time are inserted in Table 19. Based on this table, the MDI & MKE are approximately same. The WCM analysis duration is about 55% of the MRE, MDI and MKE methods. The results are inserted in Table 20. Based on the criteria of S1, S2 & S3, the WCM is the most appropriate tactic. In other words, in addition to its high accuracy, its convergence rate and time are appropriate.

Open hemispherical shell

Figure 33 shows a hemispherical shell with an 18° opening at the top [⁸⁷,⁸⁸]. The shell is loaded by two orthogonal pairs of diametrically opposite equal forces. One pair is pinched along one direction whereas the other pair is pulled along the perpendicular direction. Owing to symmetry, only one-quarter of the shell is loaded and modeled using 380 elements. The modulus of elasticity, thickness, radius’ and maximum load are 6.825 × 10⁶, 0.04, 10, and 400, respectively.

The load-displacement curves are plotted in Figs. 34(a) and 34(b) for the horizontal and vertical displacements, correspondingly. Based on Fig. 34(b), a snap-back point exists in the vertical deflection curve of point A. In this point; the load and displacement are 113.9256 and 0.939, respectively. These results are compatible with the previously reported one [⁸⁹]. For this curve, The MDI, MKE and WCM methods are more accurate. Based on Ref. [⁸⁹], the vertical deflection at the maximum load is 0.0927. Figure 35 shows the deformed structural configurations. All methods are able to trace the structural static path completely. However, their numbers of convergence points are different. Table 21 lists the number of convergence points, iterations and time whereas Table 22 lists the scores.

Semi-cylindrical shell

This example is a well-known benchmark problem. Figure 36 shows the geometry of this structure [⁸⁷]. The maximum load is 1600. The load is applied to the middle of the free end. The other circumferential periphery is fully clamped. The two straight edges parallel to Y are restrained. In other words, the nodal Z translations and rotation about Y are zero, i.e., W= q_Y = 0. Young's modulus, Poisson's ratio and shell thickness are 2.0685 × 10⁷, 0.3, and 0.03, respectively. Owing to symmetry, one-half of the semi-cylindrical shell is modeled by 20 × 12 shell elements.

Figure 37 shows the equilibrium static path for A & B points. The obtained results are compatible with the outcomes of the other researchers [⁸⁷]. Based on the aforementioned figures, the snap-back point exists in the curves related to the point B. For the vertical deflection, it is occurred when the load is reached to 875.27. The corresponding deflection is 3.33. For the horizontal displacement, the related load and displacement are 946.9 and 24.1, respectively. The local jumps have been occurred in the MRF, but the other methods can appropriately trace the static path. The deformed configuration at the maximum load is portrayed in the Fig. 38. All the schemes are ranked in Table 23. According to the obtained results, it is concluded that WCM strategy is the most suitable technique. Moreover, the MRF method is the worst tactic for this structure.

Imperfection sample

As it is shown in Fig. 39, the final example is a planar, static, conservative system with two degrees of freedom [⁹⁰]. Since buckling is very sensitive to the imperfections, this structure is studied to illustrate the special situation of the zero-stiffness post buckling. In this study, the load is applied to the node D at vertical direction. The maximum load is 8896.8 kg, and Young’s modulus is 20685.1 kg/cm². It is assumed that the cross-sectional areas are 25.81, 25.81, 19.35, and 48.39 cm² for AB, BC, CD, and DE rods, respectively.

Figure 40 shows the equilibrium static path for C point at horizontal direction. All presented tactics are able to completely trace the structural equilibrium path. However, the numbers of convergence points belong to these schemes are different. For accurate validation, this truss is solved by Crisfield’s arc-length, as well. The arc-length response is the same as DR curves. Based on Fig. 40, the snap-back points exist in the load-displacement curve. The first point is occurred when the load is reached the value of 197531.2167 kg, with the corresponding displacement of 109.8 cm. Besides, the second limit displacement point is occurred in 2925144.241 kg and with displacement 56.52 cm. For better investigation, the details of the converged points near the buckling region are portrayed in Fig. 41. This structure is deformed in a low load factor about 8896.8 kg. The horizontal displacement is 52.7 cm in this force. Under application of more loads, the truss resists versus deformation. In other words, the deformation is produced by more applied load. After the first limit point occurs, the displacement returns and the force increases. All the used schemes are ranked in Table 24. According to the obtained results, it is concluded that MKE and WCM strategies are the most suitable techniques for solving this structure. Moreover, the MDI method is the worst approach for this sample.

Ranking evaluation

Based on the aforementioned indicators, a number is assigned to each technique. This number is ranged from 1 to 5. Number one is ascribed to the best scheme. Number 5 represents the worst tactic. To quantify the efficiency of the presented algorithms, the ranking of the i^th method is defined by Q_ij. For instance, based on indicator S2, the WCM strategy is ranked first in 13 solved problems. Hence, Q_i1 is equal to 13 for this approach. The number of structures which cannot be analyzed by using of ith tactic is shown by Q_0i. Based on Q_ij values, the score of the ith scheme is computed as below:

(36)G_{ij}=\mathrm{1}\mathrm{0}\mathrm{0}\times {\displaystyle \underset{j=\mathrm{1}}{\overset{\mathrm{5}}{\mathrm{\Sigma }}}\frac{Q_{ij}(\mathrm{6}-j)}{\mathrm{85}}}\mathrm{.}

When a method is ranked first in all the samples, its related score is 85. Consequently, G_ij of this technique equals 100. The scores of the presented strategies are illustrated in Tables 25–27. When a procedure is not able to trace the static equilibrium path, it is not considered in Eq. ‎(36).

According to the obtained results, the rankings of all tactics are presented in Table 28. This table shows that the worst method is MKE tactic. Based on the two criteria, authors’ algorithm is ranked first.

Conclusion

In early part of the paper, DRM and its related formulas were introduced. Afterwards, the existing approaches for determining the load factor were discussed briefly. Finally, a new formulation for load factor was suggested, which took advantage of the external work. This procedure assumes that the work increment is zero for all iterations. By using the suggested scheme, it is possible to trace the equilibrium path of structure with extremely nonlinearity, even in the snap-back and snap-through parts.

The proposed tactic is entirely automatic and only uses DR fictitious parameters. The efficiency of the suggested formulation was evaluated with the analysis of 2D and 3D trusses, arcs, 2D frames and shells. All of these structures had geometrical nonlinear behavior. If the number of convergence points is low, then the time of analysis reduces. On the other hand, for more accurate structural static path, analysis time increases. Therefore, analysis consuming time is an important factor. In this study, the accuracy, analysis time and no jumping or deviations of the structural static path are considered. These three criteria were used to compare the suggested technique with other strategies. The numerical examples have proven that the proposed formulation is more accurate than MRE and MRF. Moreover, the accuracy of the suggested algorithm is similar with MDI and MKE. Based on the numbers of convergence points, total numbers of iterations, and the analysis consuming time, it is concluded that authors’ technique is more efficient, in comparison to the other presented methods.

References

[1] Rezaiee H. Nonlinear Structural Analysis Using Dynamic Relaxation Method. Thesis for the Master's Degree. Mashhad: Ferdowsi University, 2012

[2] Chen W F, Lui E M. Stability design of steel frames. Boca Raton: CRC press, 1991

[3] Zienkiewicz O C. Incremental displacement in non-linear analysis. International Journal for Numerical Methods in Engineering, 1971, 3(4): 587–588

[4] Wempner G A. Discrete approximations related to nonlinear theories of solids. International Journal of Solids and Structures, 1971, 7(11): 1581–1599

[5] Riks E. The application of Newton’s method to the problem of elastic stability. Journal of Applied Mechanics, 1972, 39(4): 1060–1065

[6] Riks E. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, 1979, 15(7): 529–551

[7] Ramm E. Strategies for tracing the nonlinear response near limit points. In: Wunderlich W, Stein E, Bathe K J, eds. Nonlinear Finite Element Analysis in Structural Mechanics. Heidelberg: Springer Berlin Heidelberg, 1981, 63–89

[8] Crisfield M A. A fast incremental/iterative solution procedure that handles “snap-through”. Computers & Structures, 1981, 13(1–3): 55–62

[9] Bergan P G. Solution algorithms for nonlinear structural problems. Computers & Structures, 1980, 12(4): 497–509

[10] Krenk S, Hededal O. A dual orthogonality procedure for non-linear finite element equations. Computer Methods in Applied Mechanics and Engineering, 1995, 123(1–4): 95–107

[11] Rezaiee-Pajand M, Boroshaki F. A variable arc-length method. Asian Journal of Structural Engineering, 1999, 3: 21–44

[12] Kim J H, Kim Y H. A predictor–corrector method for structural nonlinear analysis. Computer Methods in Applied Mechanics and Engineering, 2001, 191(8–10): 959–974

[13] Toklu Y C. Nonlinear analysis of trusses through energy minimization. Computers & Structures, 2004, 82(20–21): 1581–1589

[14] Ligarò S S, Valvo P S. Large displacement analysis of elastic pyramidal trusses. International Journal of Solids and Structures, 2006, 43(16): 4867–4887

[15] Saffari H, Mansouri I. Non-linear analysis of structures using two-point method. International Journal of Non-linear Mechanics, 2011, 46(6): 834–840

[16] Day A S. An introduction to dynamic relaxation. Engineer, 1965, 219: 218–221

[17] Otter J R H. Computations for prestressed concrete reactor pressure vessels using dynamic relaxation. Nuclear Structural Engineering, 1965, 1(1): 61–75

[18] Otter J R H, Day A S. Tidal computations. Engineer, 1960, 289: 177–182

[19] Frankel S P. Convergence rates of iterative treatments of partial differential equations. Mathematical Tables and Other Aids to Computation, 1950, 4(30): 65–75

[20] Rushton K R. Large deflection of variable-thickness plates. International Journal of Mechanical Sciences, 1968, 10(9): 723–735

[21] Alamatian J. Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method. Computers & Structures, 2013, 114–115: 84–97

[22] Lee K S, Han S E, Park T. A simple explicit arc-length method using the dynamic relaxation method with kinetic damping. Computers & Structures, 2011, 89(1–2): 216–233

[23] Lee K S, Han S E, Hong J W. Post-buckling analysis of space frames using concept of hybrid arc-length methods. International Journal of Non-linear Mechanics, 2014, 58: 76–88

[24] Rezaiee-Pajand M, Alamatian J. Dynamic relaxation method for tracing the statical path of truss structures. Journal of Modeling in Engineering, 2009, 3(17): 27–39

[25] Rezaiee-Pajand M, Alamatian J. Automatic DR structural analysis of snap-through and snap-back using optimized load increments. Journal of Structural Engineering, 2011, 137(1): 109–116

[26] Bagrianski S, Halpern A B. Form-finding of compressive structures using Prescriptive Dynamic Relaxation. Computers & Structures, 2014, 132: 65–74

[27] Barnes M R. Form-finding and analysis of tension space structures by dynamic relaxation. Dissertation for the Doctoral Degree. Ann Arbor: The City University, 1977

[28] Barnes M R. Form-finding and analysis of prestressed nets and membranes. Computers & Structures, 1988, 30(3): 685–695

[29] Barnes M R. Form and stress engineering of tension structures. Structural Engineering Review, 1994, 6(3): 175–202

[30] Barnes M R. Form finding and analysis of tension structures by dynamic relaxation. International Journal of Space Structures, 1999, 14(2): 89–104

[31] Han S E, Lee K S. A study of the stabilizing process of unstable structures by dynamic relaxation method. Computers & Structures, 2003, 81(17): 1677–1688

[32] Hegyi D, Sajtos I, Geiszter G, Hincz K. Eight-node quadrilateral double-curved surface element for membrane analysis. Computers & Structures, 2006, 84(31–32): 2151–2158

[33] Lewis W J, Lewis T S. Application of formian and dynamic relaxation to the form-finding of minimal surfaces. Journal of the International Association for Shell and Spatial Structures, 1996, 37(3): 165–186

[34] Wood R D. A simple technique for controlling element distortion in dynamic relaxation form-finding of tension membranes. Computers & Structures, 2002, 80(27–30): 2115–2120

[35] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[36] Hamdia K M, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. International Journal of Fracture, 2017, 206(2): 215–227

[37] Hamdia K M, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 95–109

[38] Cassell A, Kinsey P, Sefton D. Cylindrical shell analysis by dynamic relaxation. In: Proceedings of the Institution of Civil Engineers. Ice Virtual Library, 1968,75–84

[39] Otter J, Pippard A, Lane R, Welch A, King I, Wood W, Cubitt N, Hayes R, Hobbs R, Zienkiewicz O. Discussion: dynamic relaxation. In: Proceedings of the Institution of Civil Engineers. Ice Virtual Library, 1967, 723–750

[40] Wood W L. Note on dynamic relaxation. International Journal for Numerical Methods in Engineering, 1971, 3(1): 145–147

[41] Brew J S, Brotton D M. Nonlinear structural analysis by dynamic relaxation. International Journal for Numerical Methods in Engineering, 1971, 3(4): 463–483

[42] Bunce J W. A note on the estimation of critical damping in dynamic relaxation. International Journal for Numerical Methods in Engineering, 1972, 4(2): 301–303

[43] Cassell A C, Hobbs R E. Numerical stability of dynamic relaxation analysis of non-linear structures. International Journal for Numerical Methods in Engineering, 1976, 10(6): 1407–1410

[44] Felippa C A. Dynamic relaxation under general increment control. Mathematical Programming, 1982, 24: 103–133

[45] Underwood P. Dynamic Relaxation (in Structural Transient Analysis). Computational Methods for Transient Analysis. Amsterdam: North-Holland, 1983: 245–265

[46] Qiang S. An adaptive dynamic relaxation method for nonlinear problems. Computers & Structures, 1988, 30(4): 855–859

[47] Zhang L C, Kadkhodayan M, Mai Y W. Development of the maDR method. Computers & Structures, 1994, 52(1): 1–8

[48] Munjiza A A. K^m proportional damping for dynamic relaxation. International Journal for Engineering Modelling, 1996, 9(1–4): 1–9

[49] Munjiza A, Owen D R J, Crook A J L. An M( M^-1K)^m proportional damping in explicit integration of dynamic structural systems. International Journal for Numerical Methods in Engineering, 1998, 41(7): 1277–1296

[50] Rezaiee-Pajand M, Taghavian Hakkak M. Nonlinear analysis of truss structures using dynamic relaxation. International Journal of Engineering, 2006, 19(1): 11–22

[51] Kadkhodayan M, Alamatian J, Turvey G J. A new fictitious time for the dynamic relaxation (DXDR) method. International Journal for Numerical Methods in Engineering, 2008, 74(6): 996–1018

[52] Rezaiee-Pajand M, Alamatian J. Nonlinear dynamic analysis by dynamic relaxation method. Structural Engineering and Mechanics, 2008, 28(5): 549–570

[53] Rezaiee-Pajand M, Sarafrazi S R. Nonlinear structural analysis using dynamic relaxation method with improved convergence rate. International Journal of Computational Methods, 2010, 7(4): 627–654

[54] Rezaiee-Pajand M, Alamatian J. The dynamic relaxation method using new formulation for fictitious mass and damping. Structural Engineering and Mechanics, 2010, 34(1): 109–133

[55] Rezaiee-Pajand M, Kadkhodayan M, Alamatian J, Zhang L C. A new method of fictitious viscous damping determination for the dynamic relaxation method. Computers & Structures, 2011, 89(9–10): 783–794

[56] Rezaiee-Pajand M, Sarafrazi S R. Nonlinear dynamic structural analysis using dynamic relaxation with zero damping. Computers & Structures, 2011, 89(13–14): 1274–1285

[57] Rezaiee-Pajand M, Kadkhodayan M, Alamatian J. Timestep selection for dynamic relaxation method. Mechanics Based Design of Structures and Machines, 2012, 40(1): 42–72

[58] Rezaiee-Pajand M, Sarafrazi S R, Rezaiee H. Efficiency of dynamic relaxation methods in nonlinear analysis of truss and frame structures. Computers & Structures, 2012, 112–113(0): 295–310

[59] Rezaiee-Pajand M, Estiri H. Computing the structural buckling limit load by using dynamic relaxation method. International Journal of Non-linear Mechanics, 2016, 81: 245–260

[60] Rezaiee-Pajand M, Estiri H. Finding equilibrium paths by minimizing external work in dynamic relaxation method. Applied Mathematical Modelling, 2016, 40(23–24): 10300–10322

[61] Rezaiee-Pajand M, Estiri H. Mixing dynamic relaxation method with load factor and displacement increments. Computers & Structures, 2016, 168: 78–91

[62] Rezaiee-Pajand M, Estiri H. A comparison of large deflection analysis of bending plates by dynamic relaxation. Periodica Polytechnica. Civil Engineering, 2016, 60(4): 619–645

[63] Rezaiee-Pajand M, Estiri H. Comparative analysis of three-dimensional frames by dynamic relaxation methods. Mechanics of Advanced Materials and Structures, 2017, 25(6): 451–466

[64] Rezaiee-Pajand M, Estiri H. Geometrically nonlinear analysis of shells by various dynamic relaxation methods. World Journal of Engineering, 2017, 14(5): 381–405

[65] Labbafi S F, Sarafrazi S R, Kang T H K. Comparison of viscous and kinetic dynamic relaxation methods in form-finding of membrane structures. Advances In Computational Design, 2017, 2(1): 71–87

[66] Zhang L G, Yu T X. Modified adaptive dynamic relaxation method and its application to elastic-plastic bending and wrinkling of circular plates. Computers & Structures, 1989, 33(2): 609–614

[67] Rushton K R. Post-buckling of tapered plates. International Journal of Mechanical Sciences, 1969, 11(5): 461–480

[68] Turvey G, Wittrick W. The large deflection and post-buckling behaviour of some laminated plates. Aeronautical Quarterly, 1973, 24(2): 77–86

[69] Hook P M, Rushton K R. Buckling of beams and plates onto an intermediate support studied by the dynamic relaxation method. Journal of Strain Analysis for Engineering Design, 1975, 10(3): 153–158

[70] Kadkhodayan M, Zhang L C, Sowerby R. Analyses of wrinkling and buckling of elastic plates by DXDR method. Computers & Structures, 1997, 65(4): 561–574

[71] Ramesh G, Krishnamoorthy C S. Post-buckling analysis of structures by dynamic relaxation. International Journal for Numerical Methods in Engineering, 1993, 36(8): 1339–1364

[72] Ramesh G, Krishnamoorthy C S. Inelastic post-buckling analysis of truss structures by dynamic relaxation method. International Journal for Numerical Methods in Engineering, 1994, 37(21): 3633–3657

[73] Greco M, Menin R, Ferreira I, Barros F. Comparison between two geometrical nonlinear methods for truss analyses. Structural Engineering and Mechanics, 2012, 41(6): 735–750

[74] Levy R, Spillers W R. Analysis of geometrically nonlinear structures. Chapman & Hall, 1995

[75] Vu-Bac N, Duong T X, Lahmer T, Zhuang X, Sauer R A, Park H S, Rabczuk T. A NURBS-based inverse analysis for reconstruction of nonlinear deformations of thin shell structures. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 427–455

[76] Greco M, Vicente C E R. Analytical solutions for geometrically nonlinear trusses. REM. Revista Escola de Minas, 2009, 62(2): 205–214

[77] Felippa C A. Nonlinear Finite Element Methods (ASEN 5017). Colorado: University of Colorado, 2001

[78] Yang Y B, Shieh M S. Solution method for nonlinear problems with multiple critical points. AIAA Journal, 1990, 28(12): 2110–2116

[79] Hrinda G. Snap-through instability patterns in truss structures. In: The 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference: American Institute of Aeronautics and Astronautics, 2010, 1313–1324

[80] Yang Y B, Yang C T, Chang T P, Chang P K. Effects of member buckling and yielding on ultimate strengths of space trusses. Engineering Structures, 1997, 19(2): 179–191

[81] Jagannathan D, Epstein H, Christiano P. Snap-through buckling of reticulated shells. Ice Proceedings, 1975, 59(4): 727–742

[82] Krishnamoorthy C S, Ramesh G, Dinesh K U. Post-buckling analysis of structures by three-parameter constrained solution techniques. Finite Elements in Analysis and Design, 1996, 22(2): 109–142

[83] Hangai Y, Kawamata S. Nonlinear analysis of space frames and snap-through buckling of reticulated shell structures. Proceedings of IASS Pacific Symposium on Tension Structures and Space Frames, 1972, 803–816

[84] Wood R D, Zienkiewicz O. Geometrically nonlinear finite element analysis of beams, frames, arches and axisymmetric shells. Computers & Structures, 1977, 7(6): 725–735

[85] Williams F. An approach to the non-linear behaviour of the members of a rigid jointed plane framework with finite deflections. Quarterly Journal of Mechanics and Applied Mathematics, 1964, 17(4): 451–469

[86] Sze K Y, Chan W K, Pian T H H. An eight-node hybrid-stress solid-shell element for geometric non-linear analysis of elastic shells. International Journal for Numerical Methods in Engineering, 2002, 55(7): 853–878

[87] Jeon H M, Lee Y, Lee P S, Bathe K J. The MITC3+ shell element in geometric nonlinear analysis. Computers & Structures, 2015, 146: 91–104

[88] Sze K Y, Liu X H, Lo S H. Popular benchmark problems for geometric nonlinear analysis of shells. Finite Elements in Analysis and Design, 2004, 40(11): 1551–1569

[89] Mohan P.Development and applications of a flat triangular element for thin laminated shells. Aiaa Journal, 1997, 36(2), 273–281

[90] Jia X, Hoefinger G, Mang H A. Imperfection sensitivity or insensitivity of zero-stiffness postbuckling … that is the question. Proceedings in Applied Mathematics and Mechanics, 2009, 9(1): 241–242

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature