A simple extrapolated predictor for overcoming the starting and tracking issues in the arc-length method for nonlinear structural mechanics

05/20/2020
by   Chennakesava Kadapa, et al.
Swansea University
0

This paper presents a simplified implementation of the arc-length method for computing the equilibrium paths of nonlinear structural mechanics problems. In the proposed technique, the predictor is computed by extrapolating the solutions from two previously converged load steps. The extrapolation is a linear combination of the previous solutions; therefore, it is simple and inexpensive. Additionally, the proposed extrapolated predictor also serves as a means for identifying the forward movement along the equilibrium path without the need for any sophisticated techniques commonly employed for explicit tracking. The ability of the proposed technique to successfully compute complex equilibrium paths in static structural mechanics problems is demonstrated using five numerical examples involving truss and beam models.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

02/09/2020

An adaptive homotopy method for computing bifurcations of nonlinear parametric systems

In this paper, we present an adaptive step-size homotopy tracking method...
04/22/2019

A Simple Local Variational Iteration Method and Related Algorithm for Nonlinear Science and Engineering

A very simple and efficient local variational iteration method for solvi...
07/24/2020

The TR-BDF2 method for second order problems in structural mechanics

The application of the TR-BDF2 method to second order problems typical o...
09/11/2019

A Robust Numerical Path Tracking Algorithm for Polynomial Homotopy Continuation

We propose a new algorithm for numerical path tracking in polynomial hom...
11/10/2021

A Reverse Augmented Constraint preconditioner for Lagrange multiplier methods in contact mechanics

Frictional contact is one of the most challenging problems in computatio...
11/30/2021

Towards algorithm-free physical equilibrium model of computing

Our computers today, from sophisticated servers to small smartphones, op...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Nonlinear structures experience complex deformation behaviour beyond the limit points, for example, post-buckling, plastic yielding and damage. For computing the complex nonlinear response of structures, the arc-length method has become the de-facto incremental-iterative 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 arc-length 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 arc-length method, a constraint equation, called the arc-length 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 arc-length method is known as the spherical arc-length method WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 , cylindrical arc-length method CrisfieldCandS1981 ; RammArcLength1980 , and elliptical arc-length method ParkIJNME1982 . Linearised arc-length methods in which a linearised form of the arc-length equation is employed instead of the original arc-length equation, are also available WempnerIJSS1971 ; RiksJAM1972 ; CrisfieldIJNME1983 .

Some critical issues encountered in the computer implementation of the arc-length 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 arc-length 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 , Ritto-Correa and Camotim RittoCandS2008 , Al-Rasby AlRasbyCandS1991 , Krenk and Hededal KrenkCMAME1995 , and Kouhia KouhiaCMAME2008 , and references cited therein for the details regarding the commonly-encountered issues in the arc-length method, limitations of different flavours of the arc-length 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 arc-length 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 arc-length 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 arc-length implementations.

The rest of the paper is organised as follows. The governing equations for the arc-length 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 arc-length 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 post-buckling regime, can be quite complex with curves and loops in the load-deflection 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 Newton-Raphson 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 arc-length 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 arc-length 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 under-determined system of equations. To overcome this issue, the system of nonlinear equations in Eq. (4) is augmented with an additional equation, called as the arc-length equation, see WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 . The generic form of the arc-length equation is given as

(5)

where, is the arc-length parameter which parametrises the equilibrium path FengCS1996 , is the increment in the arc-length parameter and . Here, is a scalar parameter which helps to recover different arc-length schemes. For , the cylindrical arc-length method is recovered CrisfieldCandS1981 ; RammArcLength1980 ; for , the method becomes spherical arc length method WempnerIJSS1971 ; RiksIJSS1979 ; CrisfieldCandS1981 ; and for , the elliptical arc-length method is recovered ParkIJNME1982 . For the given , Eqs. (4) and (5) are solved together for increments and using the Newton-Raphson 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 Newton-Raphson 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 arc-length method is illustrated schematically in Fig. 1.

2.1 Issues associated with the arc-length method

When solving the coupled matrix system (8), one quickly runs into four main difficulties:

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

  2. starting the arc-length algorithm at the first load step,

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

  4. 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 load-control method. While this might not be a big issue in in-house codes, it does require modifications to the code to account for the increased matrix size. Moreover, the increased bandwidth due to the off-diagonal 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 arc-length 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 right-hand sides. In fact, many matrix solvers, for example, SuperLU superlu , PARDISO pardiso-6.0a and MUMPS mumps1 support the solution of matrix systems for multiple right-hand 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.

Figure 1: A typical evolution of solution increments in the arc-length method.

2.1.2 Issue 2: starting the arc-length algorithm at the first load step

The lack of information regarding the arc-length 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 arc-length 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 arc-length 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 arc-length method. In the load-control 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 arc-length method since it yields and , which makes it impossible to compute from Eq. (8). A variety of techniques have been proposed for computing suitable non-zero values as the predictors at the first iteration for the arc-length 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 arc-length method. Accordingly,

(20)

where, and are the increments of arc-length 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 arc-length implementations in the literature RiksIJSS1979 ; CrisfieldCandS1981 ; RammArcLength1980 ; SchweizerhofCMAME1986 ; KrenkIJNME1995 .

Figure 2: Evolution of solution increments in the arc-length method with .

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 arc-length 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 1-DOF problem are illustrated schematically in Fig. 3 for four different scenarios along the equilibrium path for the case with uniform increments of the arc-length 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 low-cost 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.

(a)
(b)
(c)
(d)
Figure 3: Illustration of predictors for different scenarios along the equilibrium path. and •  represent the predicted and actual solutions, respectively. The direction of the predicted solution is denoted with a thick red dashed line.
(a) Uniform increment
(b) Non-uniform increment (adaptive cutting)
Figure 4: Illustration of predictors with uniform and non-uniform increments in the arc-length parameter.
1:Set , , and . converged=False. . Initialise variables.
2:Compute
3:for  to  do
4:     #1 Predictor step:
5:     if  then
6:         
7:         
8:         
9:     end if
10:     
11:     
12:     convergedPrev = converged
13:     converged = False
14:     #2 Corrector step:
15:     for  to  do
16:         Compute: , , , and in Eq. (8)
17:         if  then
18:              converged = True
19:              Exit iteration loop
20:         end if
21:         Solve: and from Eqs. (13) and (14)
22:         
23:         
24:         
25:         
26:     end for
27:     #3 Solution update:
28:     if converged then
29:         if n == 1 then
30:              
31:              
32:              
33:         end if
34:         
35:         
36:         
37:         if convergedPrev then
38:              
39:         end if
40:         
41:         
42:     else
43:         if convergedPrev then
44:              
45:         else
46:              
47:         end if
48:     end if
49:end for
Algorithm 1 Algorithm for the arc-length method

3 Numerical examples

The ability of the proposed arc-length implementation in capturing complex equilibrium paths is illustrated using five benchmark examples consisting of nonlinear truss and beam models. The pseudocode for the arc-length 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 book-fem-ZienkiewiczVol2 for the details.

The spherical arc-length 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 arc-length increment () is changed adaptively, see lines 42-48 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 load-displacement curves for node 2 in Fig. (b)b. For , the bar 1-2 acts like a rigid structure and it undergoes significantly less deformation. For this case, the load-control method results in a snap-through behaviour but accurate solutions can be obtained using the displacement-control method. However, for , the displacement-control method also fails to track the correct equilibrium path as it results in a snap-back response. To accurately compute the response of the structure over a wide range of parameters, the arc-length method is essential.

Numerical solutions are computed using the proposed arc-length implementation for four different values of and the load-displacement 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.

(a)
(b)
Figure 5: 3-member planar truss: (a) geometry and boundary conditions and (b) analytical load-displacement curves for node 2 for various values of with .
(a)
(b)
(c)
(d)
Figure 6: 3-member planar truss: load-displacement curves for the free DOFs with and different values of . The markers represent the converged load steps.

3.2 Space truss with 12 members

This example consists of a 12-bar 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 Green-Lagrange strain measure, see A. The analysis is started using a load increment of , which corresponds to an arc-length of . The response of the structure is presented in terms of load-displacement curves in Figs. 8 and 9 and displacement-displacement 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 arc-length implementation captures complex nonlinear response of the structure quite well.

Figure 7: 12-member space truss: geometry and boundary conditions.
Figure 8: 12-member space truss: load-displacement curve.
Figure 9: 12-member space truss: load-displacement curve.
Figure 10: 12-member space truss: graph of displacement versus displacement .

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 snap-back 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 arc-length increment of . The load-displacement 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 load-displacement graphs presented for the example illustrate that the proposed technique successfully captures the load limit point as well as the displacement limit point.

(a)
(b)
Figure 11: Lee frame: (a) geometry and boundary conditions and (b) deformed configurations of the frame at points a, b, c and d marked in Fig. 12.
Figure 12: Lee frame: load-displacement curve. The markers for the present work represent the converged load steps.

3.4 Hinged-clamped 215-degree arch

This is another widely-used 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 load-deflection curve in terms of normalised load and normalised displacement is presented in Fig. (b)b. The load-displacement 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 arc-length 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.

(a)
(b)
Figure 13: Hinged-clamped 215 arch: (a) geometry and boundary conditions and (b) the load-displacement curve. The markers for the present work represent the converged load steps.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 14: Hinged-clamped 215 arch: deformed shapes at eight different load steps.   : original configuration and    : deformed configuration.

3.5 Semi-circular hinged arch

The last example consists of a semi-circular 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 off-set 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 load-displacement curves in Fig. 16 show that the proposed technique captures the looping paths of the load-displacement 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 arc-length methods.

Figure 15: Semi-circular arch: geometry and boundary conditions. Length units are centimeters.
(a) Symmetric loading
(b) Asymmetric loading
Figure 16: Semi-circular arch: load-displacement curves for symmetric and asymmetric loading.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
Figure 17: Semi-circular arch: deformed shapes at different load steps for symmetric loading.   : original configuration and    : deformed configuration.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
Figure 18: Semi-circular arch: deformed shapes at different load steps for asymmetric loading.   : original configuration and    : deformed configuration.

4 Summary and conclusions

This contribution presents a simple extrapolation operator for overcoming the well-known issues associated with the predictor step as well as with identifying the correct direction along the equilibrium path using the arc-length 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 arc-length 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 ad-hoc parameters. The simplistic nature of the proposed technique renders it suitable for adaption in the existing computer implementations of the arc-length 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 arc-length methods.

In the present work, the standard Newton-Raphson method is employed. As the future work, the proposed technique can be explored with the modified Newton-Raphson 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 arc-length 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 part-funded by the European Regional Development Fund (ERDF) via the Welsh Government.

Appendix A Nonlinear truss element

For the two-noded 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 Green-Lagrange 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 snap-through. Computers and Structures, 13:55–62, 1981.
  • (5) P. G. Bergan, G. Horrigmoe, B. Krakeland, and T. H. Soreide. Solution technioues for non-linear 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 Euro-US Workshop on Nonlinear Finite Element Analysis in Structural Mechanics, editor, Bathe, K. J. and Stein, E. and Wunderlich, W., pages 63–89, Ruhr-Universtä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 non-linear 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 thin-walled 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 non-linear 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 arc-length 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(1-2):99–110, 1987.
  • (16) W. F. Lam and C. T. Morley. Arc-length method for passing limit points in structural calculation. Journal of Structural Engineering, 118(1):169–185, 1992.
  • (17) E. Carrera. A study on arc-length-type 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 arc-length methods. Computers and Structures, 58(3):479–485, 1996.
  • (19) M. Ritto-Correa and D. Camotim. On the arc-length and other quadratic control methods: established, less known and new implementation procedures. Computers and Structures, 86(11-12):1353–1368, 2008.
  • (20) S. N. Al-Rasby. 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 non-linear finite element equations. Computer Methods in Applied Mechanics and Engineering, 123(1-4):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(13-16):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 large-scale genomic prediction with marker-by-environment 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 post-buckling 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. Kanok-Nukulchai. Geometrically non-linear analysis of arbitrary elastic supported plates and shells using an element-based Lagrangian shell element. International Journal of Non-Linear Mechanics, 43:53–64, 2008.
  • (32) I. Kreja and R. Schmidt. Large rotations in first-order shear deformation FE analysis of laminated shells. International Journal of Non-Linear 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.