1 Introduction
Nonlinear structures experience complex deformation behaviour beyond the limit points, for example, postbuckling, plastic yielding and damage. For computing the complex nonlinear response of structures, the arclength method has become the defacto incrementaliterative numerical technique in computational structural mechanics using the finite element method (FEM). Based on the techniques originally pioneered by Wempner WempnerIJSS1971 , Riks RiksJAM1972 ; RiksIJSS1979 and Crisfield CrisfieldCandS1981 , numerous flavours of the arclength method have been proposed for computing complex equilibrium paths in nonlinear static structural mechanics problems, see Bergan et al. BerganIJNME1978 , Batoz and Dhatt BatozIJNME1979 , Ramm RammArcLength1980 , Powell and Simons PowellIJNME1981 , Fried FriedCMAME1984 , Gierlinski and Smith GierlinskiCandS1985 , Schweizerhof and Wriggers SchweizerhofCMAME1986 , and Krenk KrenkIJNME1995 .
In the arclength method, a constraint equation, called the arclength equation, is added to the original (discrete) nonlinear equations of the problem. The augmented system of equations thus obtained is solved for the incremental load factor along with the incremental displacements. Depending upon the type of the constraint equation, the arclength method is known as the spherical arclength method WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 , cylindrical arclength method CrisfieldCandS1981 ; RammArcLength1980 , and elliptical arclength method ParkIJNME1982 . Linearised arclength methods in which a linearised form of the arclength equation is employed instead of the original arclength equation, are also available WempnerIJSS1971 ; RiksJAM1972 ; CrisfieldIJNME1983 .
Some critical issues encountered in the computer implementation of the arclength methods are those associated with computing the solution of a matrix system that is increased in size; finding the solution of a quadratic equation for the load parameter increment, especially the case of complex roots of the quadratic equation; computation of the arclength radius at the beginning of each load step; and the determination of the sign of the load parameter increment at the first iteration of each load step to ensure forward movement along the equilibrium path.
Over the years, several researchers have proposed different techniques for overcoming these issues towards computing complex equilibrium paths in static nonlinear structures. These techniques consist of one or a combination of (a) imposing the orthogonality conditions, (b) computing the angle between residuals and displacement increments, (c) imposing certain constraints, and (d) computing the sign of the determinant of the stiffness matrix. The discussion of all of these techniques is beyond the scope and the interest of this article. The reader is referred to Crisfield CrisfieldCandS1981 , Bergan et al. BerganIJNME1978 , Krenk KrenkIJNME1995 , Bellini and Chulya BelliniCandS1987 , Lam and Morley LamJSE1992 , Carrera CarreraCandS1994 , Feng et al. FengCS1996 , RittoCorrea and Camotim RittoCandS2008 , AlRasby AlRasbyCandS1991 , Krenk and Hededal KrenkCMAME1995 , and Kouhia KouhiaCMAME2008 , and references cited therein for the details regarding the commonlyencountered issues in the arclength method, limitations of different flavours of the arclength method and various techniques proposed for overcoming them.
Despite their tremendous success in computing the complex equilibrium paths in the response of nonlinear structures, the existing techniques for overcoming the starting and tracking in the arclength are cumbersome and expensive. In this work, a simple extrapolation operator is proposed for overcoming the critical issues related to the starting and tracking in the arclength method. In the proposed technique, the solution at the predictor step is computed as a linear combination of solutions at the two previously converged load steps. Therefore, the proposed technique is simple, inexpensive and easier to implement. Furthermore, the proposed technique also helps in tracking the solution in the forward direction along the equilibrium path, without the need for any sophisticated techniques commonly employed in the classical arclength implementations.
The rest of the paper is organised as follows. The governing equations for the arclength method and the details of the proposed technique are discussed in Section 2. The suitability of the proposed technique in successfully computing complex equilibrium paths is illustrated in Section 3 using five benchmark examples involving nonlinear truss and beam models. The paper is concluded with Section 4 with a summary of observations made and conclusions drawn from the present work.
2 The arclength method
By adapting the finite element method (FEM) for computing numerical solutions, the governing discrete system of equations for the nonlinear elasticity problem can be written as,
(1) 
where,
is the nodal displacement vector,
is the load factor, is the internal force vector, is the external force vector, and is the residual vector.Often, it is impossible to obtain the numerical solutions of Eq. (1) when the external load is applied all at once. Moreover, the response of the structures, especially in the postbuckling regime, can be quite complex with curves and loops in the loaddeflection curves. Therefore, numerical solutions of Eq. (1) are obtained by using an incremental approach, in which the solutions are computed using an iterative technique, for example, the NewtonRaphson scheme, at each increment.
In the incremental approach, the displacement and load factor at the current load step, and , respectively, are computed as increments, and , from their respective values at the previously converged load step, and , as
(2)  
(3) 
where, the subscripts and , respectively, denote the current and previously converged load steps. Using (2) and (3), the residual vector for the current load step can be written as,
(4) 
If the region of interest of the response of the structure under consideration does not include any limit points, then either load control method (LCM) or the displacement control method (DCM) can be used. However, if one is interested in tracking the response of the structure beyond the limit points, then the arclength method (ALM) must be employed, see Leon et al. LeonAMR2011 and references cited therein for the comprehensive details on various iterative techniques for numerical solutions of nonlinear structures.
In the arclength method, the system of nonlinear equations in Eq. (4) is solved for both the displacement, , as well as the loading parameter,
. This approach increases the number of degrees of freedom (DOFs) by one, making Eq. (
4) an underdetermined system of equations. To overcome this issue, the system of nonlinear equations in Eq. (4) is augmented with an additional equation, called as the arclength equation, see WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 . The generic form of the arclength equation is given as(5) 
where, is the arclength parameter which parametrises the equilibrium path FengCS1996 , is the increment in the arclength parameter and . Here, is a scalar parameter which helps to recover different arclength schemes. For , the cylindrical arclength method is recovered CrisfieldCandS1981 ; RammArcLength1980 ; for , the method becomes spherical arc length method WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 ; and for , the elliptical arclength method is recovered ParkIJNME1982 . For the given , Eqs. (4) and (5) are solved together for increments and using the NewtonRaphson scheme.
Starting with an initial guess , the displacement increment, , and the load increment, , are computed by iterative updates, and , as
(6) 
where, is the iteration counter and is the maximum number of iterations. Using the expressions in Eq. (6), and at the current iteration can be written as,
(7) 
By applying the NewtonRaphson scheme to solve Eqs. (4) and (5), we get the following matrix system,
(8) 
where,
(9)  
(10)  
(11)  
(12) 
The matrix system in Eq. (8) is solved for until a predefined convergence criterion is satisfied. Typical evolution of solution increments in the arclength method is illustrated schematically in Fig. 1.
2.1 Issues associated with the arclength method
When solving the coupled matrix system (8), one quickly runs into four main difficulties:

solving a matrix system of equations that is increased in size,

starting the arclength algorithm at the first load step,

predicting the solution increments at the first iteration for the subsequent load steps, and

identifying the correct (forward) direction of evolution along the equilibrium path.
These difficulties, along with the techniques for overcoming them, are discussed below.
2.1.1 Issue 1: solving a matrix system of equations that is increased in size
The size of the matrix system given by Eq. (8) is one higher than the stiffness matrix system obtained with the FEM with the loadcontrol method. While this might not be a big issue in inhouse codes, it does require modifications to the code to account for the increased matrix size. Moreover, the increased bandwidth due to the offdiagonal terms and deteriorates the performance of the matrix solver. To overcome this issue, different approaches were proposed in the literature, c.f. Crisfield CrisfieldCandS1981 , Batoz and Dhatt BatozIJNME1979 , and Schweizerhof SchweizerhofCMAME1986 , all of which solve Eq. (8) by splitting it into a suitable form such that the matrix system of the original size can be used.
In the present work, the splitting approach based on the Schur complement of the stiffness matrix () BatozIJNME1979 , which is also valid for the cylindrical arclength method (), is adapted. Accordingly, the coupled matrix system given in Eq. (8) is solved for and as
(13)  
(14) 
where,
(15) 
This approach is similar to the one proposed by Crisfield CrisfieldCandS1981 in which a quadratic equation is solved for . Note that, in this approach, the case of complex roots as well as the ambiguity associated with choosing the correct solution of the quadratic equation, are completely avoided.
Remark I: The disadvantage of the splitting scheme is that two matrix solves are required for computing and at every iteration. The associated cost can be minimised by factorising the matrix once and then using the factorisation for computing the solution of multiple righthand sides. In fact, many matrix solvers, for example, SuperLU superlu , PARDISO pardiso6.0a and MUMPS mumps1 support the solution of matrix systems for multiple righthand sides at once. An alternative is to compute at the first iteration only; this, however, deteriorates the convergence of iterations. In this work, both and are solved for at every iteration.
2.1.2 Issue 2: starting the arclength algorithm at the first load step
The lack of information regarding the arclength increment makes it difficult to start the algorithm at the first load step. However, since the limit points are not encountered in the first few load steps, the difficulty of starting the arclength method can be easily overcome by computing displacement at the first load step using the load control method in which the load increment is specified as the user input. Based on the author’s experience, a value of that produces a noticeable deformation of the structure from the original configuration, is sufficient enough.
For the first load step, the values of , and in the matrix system (8) are taken as,
(16) 
so that
(17) 
Once is obtained at the first load step using the load control method, the arclength increment can be computed using Eq. (5) as
(18) 
The value of thus computed is used in the subsequent load steps, either with or without a multiplication factor ().
2.1.3 Issue 3: predicting the solution increments at the first iteration for the subsequent load steps
This issue, together with issue 4, are the most crucial and difficult issues in the successful implementation of the arclength method. In the loadcontrol method, the solution at the previously converged load step is often employed as the predictor at the current load step. Using this choice, we get,
(19a)  
(19b) 
Although the predicted displacement increment in Eq. (19a) works for the load control method, the choice results in all kinds of problems in the arclength method since it yields and , which makes it impossible to compute from Eq. (8). A variety of techniques have been proposed for computing suitable nonzero values as the predictors at the first iteration for the arclength methods. These techniques vary in the complexity of implementation, computational cost incurred and their success in computing the complex equilibrium paths.
In the proposed work, an extrapolation operator based on the solutions at the two previously converged load steps is employed for predicting the solution increments at the first iteration in each load step (except the first load step) of the arclength method. Accordingly,
(20) 
where, and are the increments of arclength parameter, respectively, for the current and previously converged load steps. The parameter accounts for the adaptive load stepping. For uniform load increments, .
Using the predictor given by Eq. (20), the solution increments at the predictor step become
(21a)  
(21b) 
where, and , respectively, are the displacement increment and load factor increment at the previously converged load step.
The above equations indicate that the initial guess for the solution increment at the current load step is a constant multiple of the solution increments at the previously converged load step. Therefore, the predictor given by Eq. (21
) is, in fact, a better estimate than zero, and it brings the initial guess
closer to the vicinity of the actual solution , as illustrated in Fig. 2.It is worth pointing that the proposed predictor step is simple and also inexpensive when compared with the techniques employed in the classical arclength implementations in the literature RiksIJSS1979 ; CrisfieldCandS1981 ; RammArcLength1980 ; SchweizerhofCMAME1986 ; KrenkIJNME1995 .
2.1.4 Issue 4: identifying the correct (forward) direction of evolution along the equilibrium path
Identification of the correct direction of movement along the equilibrium path is another crucial issue in the successful implementation of the arclength method. Popular techniques proposed for overcoming this issue consist of sophisticated methods which require comparison of the sign of vector products, evaluation of the sign of the determinant the stiffness matrix, enforcement of orthogonal conditions, or a combination of these. Although these techniques have been proven to be successful in computing the complex equilibrium paths in the response of static nonlinear structures, their computer implementation is quite cumbersome.
In the present work, no special technique is employed for determining and/or identifying the correct direction of evolution along the equilibrium path. This is taken care by the extrapolation operator used for predicting the solution increments at the first iteration given by Eqs. (21a) and (21b). Towards understanding the reason(s) behind the ability of the proposed technique to successfully compute complex equilibrium paths, the predicted solutions for the 1DOF problem are illustrated schematically in Fig. 3 for four different scenarios along the equilibrium path for the case with uniform increments of the arclength parameter. A similar illustration is presented in Fig. 4 for the case with an adaptive cutting. As illustrated in Figs. 3 and 4, the solution predicted using the proposed technique is always located on the same side of the actual solution . Another interpretation is that the direction of is always opposite to that of , thereby forcing the solution increments at the predictor step in the correct direction along the equilibrium path.
Therefore, the proposed seemingly simple and lowcost extrapolation operator given by Eq. (20) not only predicts the solution at the first iteration but also serves to identify the correct direction along the equilibrium path, without the need for any sophisticated techniques. This ability of the proposed technique to successfully compute complex equilibrium paths in nonlinear structural mechanics problems is illustrated with numerical examples in Section 3.
3 Numerical examples
The ability of the proposed arclength implementation in capturing complex equilibrium paths is illustrated using five benchmark examples consisting of nonlinear truss and beam models. The pseudocode for the arclength method using the proposed predictor is presented in Algorithm. 1. The nonlinear space truss finite element models are discussed briefly in A. For the examples modelled with beams, the geometrically exact beam element is used; the reader is referred to Chapter 17 in Zienkiewicz and Taylor bookfemZienkiewiczVol2 for the details.
The spherical arclength method () is used in all the simulations reported in this work. The convergence tolerance () is assumed to be , and the maximum number of iterations, , is set to 10. When the convergence is not achieved according to these two criteria, the arclength increment () is changed adaptively, see lines 4248 in Algorithm. 1. The specified value of the point load is one, unless stated otherwise explicitly.
3.1 Plane truss with three members
The first numerical example consists of a planar truss with 3 bars which are arranged in the configuration shown in Fig. (a)a. For this problem, the truss model based on the engineering strain () is used, see A. This truss structure experiences a highly nonlinear deformation behaviour, as illustrated with loaddisplacement curves for node 2 in Fig. (b)b. For , the bar 12 acts like a rigid structure and it undergoes significantly less deformation. For this case, the loadcontrol method results in a snapthrough behaviour but accurate solutions can be obtained using the displacementcontrol method. However, for , the displacementcontrol method also fails to track the correct equilibrium path as it results in a snapback response. To accurately compute the response of the structure over a wide range of parameters, the arclength method is essential.
Numerical solutions are computed using the proposed arclength implementation for four different values of and the loaddisplacement curves are presented in Fig. 6 together with the analytical solution. As shown, the proposed technique captures equilibrium paths for all four cases quite accurately. It is also worth mentioning that the proposed implementation does not require prohibitively smaller increments for successful computation of numerical solutions.
3.2 Space truss with 12 members
This example consists of a 12bar space truss structure whose geometry and boundary conditions are shown in Fig. 7. This example has been previously studied by Yang and Leu YangCMAME1991 , Yang et al. YangMAMS2007 , Krenk and Hededal KrenkCMAME1995 , Leon et al. LeonAMR2011 , and Habibi and Bidmeshki HabibiIJSSD2019 . The parameters are: and . Each bar is discretised with one nonlinear truss element based on the GreenLagrange strain measure, see A. The analysis is started using a load increment of , which corresponds to an arclength of . The response of the structure is presented in terms of loaddisplacement curves in Figs. 8 and 9 and displacementdisplacement curve in Fig. 10. These graphs illustrate that the results obtained using the proposed technique match well with the solution obtained by Krenk and Hededal KrenkCMAME1995 . The results obtained for the present example show that the proposed arclength implementation captures complex nonlinear response of the structure quite well.
3.3 Lee frame
This example consists of a planar frame whose geometry and boundary conditions are as shown in Fig. (a)a. The geometric and material parameters are: cm, cm, kN/cm, and . The point load is taken as 1 kN. In this example, the snapback behaviour occurs once the frame undergoes a significantly large deformation. The analysis is performed using 20 nonlinear beam elements and with an initial load increment , which corresponds to an arclength increment of . The loaddisplacement curve for the node at which the point load is applied is shown in Fig. 12 along with a reference solution from Schweizerhof and Wriggers SchweizerhofCMAME1986 . Deformed shapes of the frame at four different points a, b, c and d in Fig. 12 are shown in (b)b. The loaddisplacement graphs presented for the example illustrate that the proposed technique successfully captures the load limit point as well as the displacement limit point.
3.4 Hingedclamped 215degree arch
This is another widelyused benchmark example for demonstrating the instabilities in structural mechanics problems. This example consists of a 215 degree circular arch of radius cm which is hinged on its one end and clamped on the other end, as depicted in Fig. (a)a. A point load is applied at the crown. The parameters are : , , , , and , in consistent units. The problem is discretised with 60 nonlinear beam elements. The initial load increment used for this problem is . The loaddeflection curve in terms of normalised load and normalised displacement is presented in Fig. (b)b. The loaddisplacement curve for this problem consists of an almost vertical jump in , and this very difficult equilibrium path is captured quite successfully using the proposed technique. The result obtained with the proposed arclength implementation is in good agreement with the reference solutions taken from Han et al. HanIJNLM2008 and Kreja and Schmidt. KrejaIJNLM2006 . The buckling load obtained in the present work has less than 1% error relative to the analytical buckling value of . Deformed shapes of the arch at eight different load instants are shown in Fig. 14.
3.5 Semicircular hinged arch
The last example consists of a semicircular arch of radius cm that is hinged on its two ends, as shown in Fig. 15. The parameters are : cm, cm, N/cm, , and . The problem is discretised with 50 nonlinear beam elements. Similar to the studies conducted in Yang and Shieh YangAIAA1990 , the analysis is performed for two different loading conditions, as shown in Fig. 15: (i) a point load at the crown of the arch which yields symmetrical deformation of the arch and (ii) a point load at an offset angle of which yields an asymmetrical deformation of the arch.
The equilibrium paths for the symmetric and asymmetric loading are presented, respectively, in Figs. (a)a and (b)b. The deformed shapes at 15 different instants are shown in Figs. 17 and 18, respectively, for the symmetric and asymmetric loading. The loaddisplacement curves in Fig. 16 show that the proposed technique captures the looping paths of the loaddisplacement curve quite well, and that the solution obtained with the proposed technique matches well with those presented in Yang and Shieh YangAIAA1990 . The ability of the proposed technique in computing such complex equilibrium paths is quite remarkable considering especially that it is completely devoid of any sophisticated techniques used for explicitly tracking the forward movement along the equilibrium path in the classical implementations of the arclength methods.
4 Summary and conclusions
This contribution presents a simple extrapolation operator for overcoming the wellknown issues associated with the predictor step as well as with identifying the correct direction along the equilibrium path using the arclength method. The proposed technique computes the predictor as a linear combination of two solutions at the previously converged load steps, which makes it simple and inexpensive when compared with the other techniques proposed in the literature on the arclength method. Moreover, the proposed approach is also applicable to the adaptive load stepping strategy. Another attractive feature of the present scheme is that it is free from adhoc parameters. The simplistic nature of the proposed technique renders it suitable for adaption in the existing computer implementations of the arclength method with some minor modifications.
The ability of the proposed scheme in successfully computing complex equilibrium paths consisting of load and displacement points as well as complex loops is demonstrated using five benchmark examples in nonlinear structural mechanics. The presented numerical results show excellent agreement between the results obtained with the proposed approach and the analytical/reference solution. It is also worth pointing out that the present approach does not require prohibitively smaller increments for its success. The capability of the proposed seemingly simple and economical predictor to capture the complex response of the structure is quite remarkable especially considering that it does not involve any sophisticated computations and comparisons commonly employed for tracking the forward movement along equilibrium path in the classical implementation of arclength methods.
In the present work, the standard NewtonRaphson method is employed. As the future work, the proposed technique can be explored with the modified NewtonRaphson method either for only or for the whole iteration. Other possibilities include combining the proposed technique with the existing strategies for improving the convergence in the vicinity of limit points towards enhancing the computational efficiency further. The ongoing work focuses on the application of the proposed technique for finite element analysis using continuum finite elements.
Supplementary material
The computer implementation of the arclength method using the proposed technique is available as GNU Octave scripts at the GitHub repository https://github.com/chennachaos/ArcLengthMethod. It is possible to use these scripts in MATLAB with some minor modifications.
Acknowledgements
The author acknowledges the support of the Supercomputing Wales project, which is partfunded by the European Regional Development Fund (ERDF) via the Welsh Government.
Appendix A Nonlinear truss element
For the twonoded space truss element, the length of the element in the original and the deformed configurations, and , respectively, are given by
(A.1)  
(A.2) 
where, and , respectively, are the coordinates of nodes 1 and 2 in the original configuration, and and are the coordinates in the current configuration. Nodal coordinates in the current configuration are related to their respective values in original configuration via the relations
(A.3)  
(A.4) 
where, and are the displacement of nodes 1 and 2, respectively.
For the truss model based on the engineering strain (), which is defined as
(A.5) 
the internal force vector () and the stiffness matrix () for an element are given by,
(A.6) 
where, is the area of the element and is the Young’s modulus,
(A.7) 
(A.8) 
For the truss model based on the GreenLagrange strain (), which is defined as
(A.9) 
the internal force vector and the stiffness matrix for an element are given by,
(A.10) 
References
References
 (1) G. A. Wempner. Discrete approximates related to nonlinear theories of solids. International Journal of Solids and Structures, 7(11):1581–1599, 1971.
 (2) E. Riks. The application of Newton’s method to the problem of elastic stability. Journal of Applied Mechanics, 39:1060–1066, 1972.
 (3) E. Riks. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures, 15:529–551, 1979.
 (4) M. A. Crisfield. A fast incremental/iterative solution procedure that handles snapthrough. Computers and Structures, 13:55–62, 1981.
 (5) P. G. Bergan, G. Horrigmoe, B. Krakeland, and T. H. Soreide. Solution technioues for nonlinear finite element problems. International Journal for Numerical Methods in Engineering, 12:1677–1696, 1978.
 (6) J. L. Batoz and G. Dhatt. Incremental displacement algorithms for nonlinear problems. International Journal for Numerical Methods in Engineering, 14(8):1262–1267, 1979.
 (7) E. Ramm. Strategies for tracing the nonlinear response near limit points. In EuroUS Workshop on Nonlinear Finite Element Analysis in Structural Mechanics, editor, Bathe, K. J. and Stein, E. and Wunderlich, W., pages 63–89, RuhrUniverstät Bochum, Germany, 1980. Springer.
 (8) G. Powell and J. Simons. Improved iterative strategy for nonlinear structures. International Journal for Numerical Methods in Engineering, 17:1455–1467, 1981.
 (9) I. Fried. Orthogonal trajectory accession to the nonlinear equilibrium curve. Computer Methods in Applied Mechanics and Engineering, 47:283–297, 1984.
 (10) J. T. Gierlinski and T. R. G. Smith. A variable load iteration procedure for thinwalled structures. Computers and Structures, 21:1085–1094, 1985.
 (11) K. H. Schweizerhof and P. Wriggers. Consistent linearization for path following methods in nonlinear FE analysis. Computer Methods in Applied Mechanics and Engineering, 59(3):261–279, 1986.
 (12) S. Krenk. An orthogonal residual procedure for nonlinear finite element equations. International Journal for Numerical Methods in Engineering, 38(5):823–839, 1995.
 (13) K. C. Park. A family of solution algorithms for nonlinear structural analysis based on relaxation equations. International Journal for Numerical Methods in Engineering, 18(9):1337–1347, 1982.
 (14) M. A. Crisfield. An arclength method including line searches and accelerations. International Journal for Numerical Methods in Engineering, 19(9):1269–1289, 1983.
 (15) P. X. Bellini and A. Chuyla. An improved automatic incremental algorithm for the efficient solution of nonlinear finite element equations. Computers and Structures, 26(12):99–110, 1987.
 (16) W. F. Lam and C. T. Morley. Arclength method for passing limit points in structural calculation. Journal of Structural Engineering, 118(1):169–185, 1992.
 (17) E. Carrera. A study on arclengthtype methods and their operation failures illustrated by a simple model. Computers and Structures, 50(2):217–229, 1994.
 (18) Y. T. Feng, D. Perić, and D. R. J. Owen. A new criterion for determination of initial loading parameter in arclength methods. Computers and Structures, 58(3):479–485, 1996.
 (19) M. RittoCorrea and D. Camotim. On the arclength and other quadratic control methods: established, less known and new implementation procedures. Computers and Structures, 86(1112):1353–1368, 2008.
 (20) S. N. AlRasby. Solution techniques in nonlinear structural analysis. Computers and Structures, 40(4):985–993, 1991.
 (21) S. Krenk and O. Hededal. A dual orthogonality procedure for nonlinear finite element equations. Computer Methods in Applied Mechanics and Engineering, 123(14):95–107, 1995.
 (22) R. Kouhia. Stabilized forms of orthogonal residual and constant incremental work control path following methods. Computer Methods in Applied Mechanics and Engineering, 197(1316):1389–1396, 2008.
 (23) S. E. Leon, G. H. Paulino, A. Pereira, I. F. M. Menezes, and E. N. Lages. A unified library of nonlinear solution schemes. Applied Mechanics Reviews, 64:040803, 2011.
 (24) X. S. Li, J. W. Demmel, J. R. Gilbert, L. Grigori, M. Shao, and I. Yamazaki. SuperLU Users’ Guide, August 2011.
 (25) A. De Coninck, B. De Baets, D. Kourounis, F. Verbosio, O. Schenk, S. Maenhout, and J. Fostier. Needles: toward largescale genomic prediction with markerbyenvironment interaction. Genetics, 203(1):543–555, 2016.
 (26) P. R. Amestoy, I. S. Duff, J. Koster, and J. Y. L’Excellent. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications, 23(1):15–41, 2001.
 (27) O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method for Solid and Structural Mechanics. Elsevier Butterworth and Heinemann, Oxford, England, Sixth edition, 2005.
 (28) Y. B. Yang and L. J. Leu. Constitutive laws and force recovery procedures in nonlinear analysis of trusses. Computer Methods in Applied Mechanics and Engineering, 92(1):121–131, 1991.
 (29) Y. B. Yang, L. J. Leu, and J. P. Yang. Key considerations in tracing the postbuckling response of structures with multi winding loops. Mechanics of Advanced Materials and Structures, 14(3):175–189, 2007.
 (30) A. Habibi and S. Bidmeshki. An optimized approach for tracing pre and postbuckling equilibrium paths of space trusses. International Journal of Structural Stability and Dynamics, 19(3):1950040, 2019.
 (31) S. C. Han, H. D. Ham, and W. KanokNukulchai. Geometrically nonlinear analysis of arbitrary elastic supported plates and shells using an elementbased Lagrangian shell element. International Journal of NonLinear Mechanics, 43:53–64, 2008.
 (32) I. Kreja and R. Schmidt. Large rotations in firstorder shear deformation FE analysis of laminated shells. International Journal of NonLinear Mechanics, 41:101–123, 2006.
 (33) Y. B. Yang and M. S. Shieh. Solution method for nonlinear problems with multiple critical points. AIAA Journal, 28(12):2110–2115, 1990.
Comments
There are no comments yet.