On the arc length method: combining ideas and implementations aspects


==1 Introduction==

In non linear finite element analysis, traditional incremental or iterative method like the standard Newton-Raphson method (NR) is widely used to trace the right equilibrium path for nonlinear problems due to  large geometry change and/or material behaviour. Common examples appear in softening and buckling  instability problems involving snap-through or snap-back phenomena.  Typical equilibrium paths associated with snap-through and snap-back behavior are illustrated  in Figure [[#img-1|1]] where appears the critical points <math display="inline">B</math>, <math display="inline">C</math>, <math display="inline">F</math>, <math display="inline">H</math> and <math display="inline">K</math>   commonly categorized into limit points and into bifurcation points respectively. However, a further disadvantage of the ''NR method'' is that without some addition of special techniques  a falling load path cannot be handled and hence no convergence at specific locations along the trajectories of load-deflection is reached.  In other words, in such situations, the standard load-controlled ''NR method'' in which  incremental displacement vector is calculated by a prescribed incremental load factor lead to error  near the peak point ( see point <math display="inline">B</math> and <math display="inline">F</math> of Figure [[#img-1|1]]). To overcome these kinds of difficulties,  often, the ''NR method'' is accompanied with some ''load or displacement controlled strategy'' and a ''line searches'' methods. However, is well known that similar failures happens when equilibrium paths encounter the limit point of snap-back  (see points <math display="inline">G</math> and <math display="inline">H</math> of Figure [[#img-1|1]]).

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-fig_1a.png|600px|snap-back and snap through response <span id='citeF-7'></span>[[#cite-7|[7]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' snap-back and snap through response <span id='citeF-7'></span>[[#cite-7|[7]]].
|}

To overcome such problems, methods and techniques of switching between load and displacement  controls (Sabir and Lock, 1972), of using the artificial springs (Wright and Gaylord, 1968),  of abandoning the equilibrium iterations in the close vicinity of the limit point (Bergan and Soreide, 1978; Bergan et al., 1978),  or of allowing load reductions within load control (Cope and Rao, 1981; Bergan and Holand, 1979; Crisfield, 1982; Phillips and Zienkiewicz, 1976) have been successively proposed.  However, these methods and techniques usually require great care and cannot always give satisfactory results. Ramm (1981) reviewed all the available methods and suggested the most versatile is the arc-length method which  was originally developed by Riks (1972, 1979) and Wempner (1971) and later modified and widely used by Ramm (1981)  and Crisfield (1981,1982a, 1983,1984)<span id='citeF-1'></span>[[#cite-1|[1]]]. After that more and more attentions (Crisfield, 1986; Lam and Morley, 1992)  were paid to improve and apply this method.

The ''AL method'' is the most powerful path following method in the solution of equilibrium paths, bifurcation points and limit points by introducing a constraint condition which establishes the nonlinear equations in which the unknown load parameter is determined<span id='citeF-2'></span>[[#cite-2|[2]]].  The ''AL method'' has been accepted as the most effective numerical method for the post-buckling problems with ''NR method''  type algorithms. For such situations the ''AL method'' is useful to avoid bifurcation points and track unloading. The ''AL method'' method causes the Newton-Raphson equilibrium iterations  to converge along an arc, as can be seen in Figure [[#img-2|2]],  thereby often preventing divergence, even when the slope of the load vs. deflection curve becomes zero or negative. The constraint equation is forced to be satisfied at each iteration.

In this dissertation we described the basic aspect of the ''AL method'' and integrated ideas to improve the standard ''AL method'' . This combinations are based on the Crisfield arc-length method, but in this case using a spherical arc length and the strategies developed by <span id='citeF-1'></span>[[#cite-1|[1]]] to avoid complex roots. For analysis of softening material, an isotropic damage model is introduced with the objective energy mesh dissipation. Benchmarking examples are  tested for proving the implemented algorithm.        

==2 Basic formulation of Arc-Length solution procedure==

The Arc-Length method is described within the category of continuation methods and it is applied to obtain solution paths<span id='citeF-3'></span>[[#cite-3|[3]]]. In the non-linear finite element analysis, the system of non-linear algebraic equations to be solve take the form:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {r}(\boldsymbol {u},\lambda ) = \boldsymbol {f}_{int}(\boldsymbol {u}) - \lambda \boldsymbol {f}_{ex}      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\boldsymbol {r} = \boldsymbol {r}(\boldsymbol {u},\lambda )</math> is a function of the displacement vector <math display="inline">\boldsymbol {u}</math> and the load parameter <math display="inline">\lambda </math>, and is the out-of-balance force vector which vanishes at equilibrium.  The internal force vector <math display="inline">\boldsymbol {f}_{in}(\boldsymbol {u})</math> is a function of the displacements and <math display="inline">\boldsymbol {f}_{ex}</math> is the prescribed external force vector and is scaled by the load intensity parameter <math display="inline">\lambda </math>.  Hence <math display="inline">\lambda \boldsymbol {f}_{ex}</math> represents the load applied at the current load increment.  Basically, the arc-length method first consider the load factor <math display="inline">\lambda </math> in Eq.[[#eq-1|1]], as a variable in the residual equation. Then, an extra new constraint equation is added to the residual equilibrium equation for defining unequivocally  the next equilibrium point solution at an intersection between the solution path and the constrains equation. This representation enables introducing <math display="inline">\lambda </math> as an additional  "degree-of-freedom", thereby expanding the dimensions of the solution space by one, for tracing the snap-through and snap-back behaviour of the curves of load against displacements.  This constrains equation can be written in differential form

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s = \int \sqrt{d\boldsymbol {u}^T d\boldsymbol {u} + \psi ^2 d\lambda ^2\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

or increment form as

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Psi = \boldsymbol {\Delta u}^T\boldsymbol {\Delta u} + \psi ^2\Delta \lambda ^2\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}-\Delta l^2 = 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math display="inline">\boldsymbol {\Delta u}</math> is vector of incremental displacement, <math display="inline">\boldsymbol {f}_{ex}</math> is specified external load vector, <math display="inline">\Delta \lambda </math> is incremental load-factor, <math display="inline">\Delta l</math> is fixed radius of desired intersection,  and <math display="inline">\psi </math> is the scaling parameter. In load-displacement space, it is designed to search the solution  of the next iteration on a small ellipsoidal surface surrounding the last converged point,  of size determined by the prescribed arc length <math display="inline">\Delta l</math>. Another constraint equations can also be found in the literature, e.g, Bathe and Dvorkin (1983).

Assuming that displacement increments <math display="inline">\boldsymbol {\Delta u}</math>  satisfying constraint Eq [[#eq-3|3]], is possible to obtain a simplified linearized version of this constraint equation   adopting the idea of Riks and Wempner(1971) for the <math display="inline">(i+1)th</math> iteration, i.e

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\delta u_i}^T\boldsymbol {\Delta u_i} + \psi ^2\delta \lambda _i\Delta \lambda _i\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex} = 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">\delta \boldsymbol {u}_i</math> and <math display="inline">\delta \lambda _i</math> are the required corrections to <math display="inline">\Delta \boldsymbol {u}_i</math> and <math display="inline">\Delta \lambda _i</math> respectively. This results in this so called fixed normal plane constraint, with  the solution iterates all lying on the hyperplane normal to the tangent at last converged point in load-displacement space<span id='citeF-1'></span>[[#cite-1|[1]]]. However, the Risk and Wempner method was not suitable for standard finite  element analysis even with modified Newton-Raphson (mN-R) procedure, because equations proposed  by Riks destroy the banded nature of the stiffness matrix <span id='citeF-4'></span>[[#cite-4|[4]]].

Therefore, we can directly used a matrix form (Crisfield 1991), that with some simplifications in Eq [[#eq-1|1]]  and [[#eq-3|3]], we compute the iterative change in displacement  vector and load-factor in the <math display="inline">(i+1)th</math> iteration,

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\begin{array}{c}\delta \boldsymbol {u}  \\  \delta \lambda   \end{array} \right]_{i+1} = -\left[\begin{array}{cc}\boldsymbol {K}_t & -\boldsymbol {f}_{ex} \\ 2\Delta \boldsymbol {p}^T_i & 2\Delta \lambda _i\psi ^2\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex} \end{array} \right]^{-1}_i \left[\begin{array}{c}\boldsymbol {r} \\ \Psi  \end{array} \right]_i  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where we define the effective general matrix of <math display="inline">n+1</math> dimension

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\hat{K}}_T = \left[\begin{array}{cc}\boldsymbol {K}_t & -\boldsymbol {f}_{ex} \\ 2\Delta \boldsymbol {p}^T_i & 2\Delta \lambda _i\psi ^2\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex} \end{array} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

where <math display="inline">\boldsymbol {K}_T</math> is the tangential stiffness matrix, and <math display="inline">\boldsymbol {r}_i</math> and <math display="inline">\Psi _i</math> are the previous values of out-of-balance load vector and arc-length. Note that <math display="inline">\boldsymbol {\hat{K}}_T</math> is a non-symmetric matrix. However, as Crisfield (1981) stated, the direct full assembling and solving of this matrix  reduces the solution efficiency, since the symmetric  banded nature  of the system stiffness matrix can no longer be exploited. After the iterative change <math display="inline">\delta \boldsymbol {u}_{i+1}</math> and <math display="inline">\delta \lambda _{i+1}</math> have been computed, the displacement vector and load-factor  are updated, i.e

<span id="eq-11"></span>
<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \boldsymbol {u}_{i+1}  = \Delta \boldsymbol {u}_{i}  + \delta \boldsymbol {u}_{i+1}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|-
| style="text-align: center;" | <math> \Delta \lambda _{i+1} = \Delta \lambda _{i} + \delta \lambda _{i+1}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}
|}

Alternatively, instead of solving directly Eq [[#eq-8|8]],  constraint equation can be introduced by following the technique of Baltoz and Dhatt (1979) for displacement control at single point. Ramm (1981) and Crisfield (1981, 1982, 1983, 1984) applying this technique, gave the modification of the method and introduced  an indirect solution scheme for the constraint equation [[#eq-3|3]]. According to this technique,  the iterative change of displacement for the new unknown load level <math display="inline">\Delta \lambda _{i+1} =\Delta \lambda _{i} + \delta \lambda _{i+1}</math>   is written as <span id='citeF-4'></span>[[#cite-4|[4]]]

<span id="eq-13"></span>
<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \boldsymbol {u}_{i+1}  = \boldsymbol {K}_T^{-1}\boldsymbol {r}_i + \delta \lambda _{i}\boldsymbol {K}_T^{-1}\boldsymbol {f}_{ex}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|-
| style="text-align: center;" | <math>                     = \delta \boldsymbol {u}_{ri}     + \delta \lambda _{i}\delta \boldsymbol {u}_{f}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}
|}

where <math display="inline">\delta \boldsymbol {u}_{ri}</math> is the displacement vector resulting from eliminating the out-of-balance forces  and <math display="inline">\delta \boldsymbol {u}_{f}</math> is the parallel displacement vector computed and stored after every update of the stiffness matrix <math display="inline">\boldsymbol {K}_T</math>.  Substituting values from [[#eq-11|11]] and [[#eq-14|14]]  into the constraint equation yields in a quadratic equations in <math display="inline">\delta \lambda _i</math>  so that

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a\delta \lambda _{i}^2 + b\delta \lambda _{i} + c = 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a =  \delta \boldsymbol {u}_{f}^{T}\delta \boldsymbol {u}_{f}  + \psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16.a)
|-
| style="text-align: center;" | <math> b =  2\delta \boldsymbol {u}_{f}(\Delta \boldsymbol {u}_i + \delta \boldsymbol {u}_{ri}) + 2\Delta \lambda _i\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16.b)
|-
| style="text-align: center;" | <math> c =  (\Delta \boldsymbol {u}_i + \delta \boldsymbol {u}_{ri})^T(\Delta \boldsymbol {u}_i + \delta \boldsymbol {u}_{ri})-\Delta l^2 + \Delta \lambda ^2\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16.c)
|}
|}

This method is graphically shown in [[#img-2|2]].  Note that the parameter <math display="inline">\psi </math> convert everything to the same dimension and order of magnitude. Many analyst set this <math display="inline">\psi=0</math>. This procedure is known as the indirect  spherical AL method Crisfield(1891). A typical value for <math display="inline">\psi ^2</math> is taken as <math display="inline">\boldsymbol {u}_0^T\boldsymbol {u}_0^T/(\lambda _0^2\boldsymbol {f}^T_{ex}\boldsymbol {f}^T_{ex})</math> <span id='citeF-1'></span>[[#cite-1|[1]]].   Eq.[[#eq-15|15]] is solved to get the value of <math display="inline">\delta \lambda </math>, and so to define the iterative change completely.  This equation leads to two results of <math display="inline">\delta \lambda </math>. Despite the sophistication of the arc-length procedure  of Ramm and Crisfield certain difficulties remain in its applications.  The first concerns to the choice of proper  roots for the quadratic equations, which the method for selection of proper value will be discussed in subsequent sections; and secondly, complex roots of a quadratic equations something occurs and the reason are not clear yet, but is strongly related to the arc-length <math display="inline">\Delta l </math> chosen.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-fig_2.png|600px|Arc-length procedure for specific iteration <span id='citeF-8'></span>[[#cite-8|[8]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Arc-length procedure for specific iteration <span id='citeF-8'></span>[[#cite-8|[8]]].
|}

==3 Modify arc-length procedure==

Lam and Morley suggested decomposing the out-of-balance load <math display="inline">\boldsymbol {r}</math> into a  parallel component <math display="inline">g_i\boldsymbol {f}_{ex}</math>, which is in the specified applied load direction; and orthogonal component <math display="inline">\boldsymbol {r}^\perp = P^\perp (g_i\boldsymbol {f}_{ex})</math>, with respect to fixed nodal  load vector in case of divergence, where

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>g_i    = \frac{\boldsymbol {r}_i^T\boldsymbol {f}_{ex}}{\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|-
| style="text-align: center;" | <math> \boldsymbol {r}^\perp = \boldsymbol {r}_i-g_i\boldsymbol {f}_{ex}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}
|}

In the equations above the, <math display="inline">P^\perp (\cdot )</math> is the orthogonal operator introduce by convenience. In its work, they proposed a new strategy to deal with the complex roots arising from the solution of the quadratic equation, which will be be discussing later.  Given the <math display="inline">\Delta \boldsymbol {u}_i</math> and the <math display="inline">\Delta \lambda _i</math>, <math display="inline">\Delta \boldsymbol {u}_{i+1}</math> and <math display="inline">\delta \lambda _{i+1}</math> have the formation represented as

<span id="eq-19.a"></span>
<span id="eq-19.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \boldsymbol {u}_{i+1}  = \Delta \boldsymbol {u}_{i} + \delta \lambda _i\boldsymbol {K}^{-1}\boldsymbol {f}_{ex} + \eta _i\boldsymbol {K}^{-1}\boldsymbol {r}^\perp </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (19.a)
|-
| style="text-align: center;" | <math>                                                                                       = \Delta \boldsymbol {u}_{i} + \delta \lambda _i\delta \boldsymbol {u}_f + \eta _i\delta \boldsymbol {u}_{ri}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (19.b)
|}
|}

and

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \lambda _{i+1}  = \Delta \lambda _i   + g_i + \delta \lambda _i                                                                                                         </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

where was introduced a scaling  parameter <math display="inline">\eta </math> in Eqn [[#eq-19.a|19.a]]  which will be discussing later. Substituting   Eqns [[#eq-20|20]] and [[#eq-19.b|19.b]] into the constraint  equations [[#eq-3|3]], gives a new quadratic equations  in <math display="inline">\Delta \lambda _i</math>, similar to Eqn [[#eq-15|15]]   written as

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a(\delta \lambda _{i})^2 + b(\delta \lambda _{i}) + c = 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where, in this case

<span id="eq-22.a"></span>
<span id="eq-22.b"></span>
<span id="eq-22.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a =  \delta \boldsymbol {u}_{f}^{T}\delta \boldsymbol {u}_{f}  + \psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22.a)
|-
| style="text-align: center;" | <math> b =  2\delta \boldsymbol {u}_{f}(\Delta \boldsymbol {u}_i + \eta \delta \boldsymbol {u}_{ri}) + 2(\Delta \lambda _i-g_i)\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22.b)
|-
| style="text-align: center;" | <math> c =  (\Delta \boldsymbol {u}_i + \eta _i\delta \boldsymbol {u}_{ri})^T(\Delta \boldsymbol {u}_i + \eta _i\delta \boldsymbol {u}_{ri})-\Delta l^2 + (\Delta \lambda ^2-g_i)\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22.c)
|}
|}

This procedure becomes general arc-length procedure proposed by Crisfield when <math display="inline">\eta=1</math>, <math display="inline">\psi ^2=0</math> and <math display="inline">g_i=0</math>.  Thus far, the basic formulation of arc-length solution procedure is completed.  However, during calculation of practical problem, there still are certain difficulties remain in application  of arc-length method. The main problem is that quadratic equation [[#eq-21|21]]   sometimes has complex roots that cause failure of convergence. Methods to solve this problem are discussed in the next part. Even we can found in the literature that is possible to get the correct answer of <math display="inline">\delta \lambda </math> without resolving any quadratic equations and avoiding the trace-back solution path when a wrong solution of <math display="inline">\delta \lambda </math>   is found and choose.

==4 Adjustment for complex roots==

As we commented before, the main problem arises in this arc length procedure came from the complex solution obtained by the solutions of the quadratic equations [[#eq-21|21]].  Particularly, we get this situations in analysis reinforced concrete when a severe strain localizations occurs and unloading of adjacent elements on cracking, leading difficulties in solution convergence still persist, while the procedure repeatedly gives complex roots. Crisfield(1983,1984) for avoiding this problem, introduced a line searches, by multiplying the both terms of Eqn [[#eq-14|14]] for by a factor <math display="inline">\eta \neq{1}</math>, determining by minimization of potential energy, obtaining a considerable improvement.         However, exist in the literature another procedures to avoid this problem, that will be described in the following.

===4.1 Method 1. Lam and Morley procedure===

Lam and Morley (1992) proposed a method in a pseudo-line search based to avoid the complex roots by introducing a scaling parameter <math display="inline">\eta </math>  in the usual constraint equation, that is, however, computed without introducing energy considerations.  According to Eqn [[#eq-21|21]], in order to get the real root, the parameter <math display="inline">a</math>, <math display="inline">b</math> and <math display="inline">c</math> in Eqn [[#eq-22.a|22.a]], [[#eq-22.b|22.b]] and [[#eq-22.c|22.c]] respectively   have to satisfy that

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b^2 - ac \geq 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Which lead to

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vartheta (\eta ) = [\psi ^2(\Delta \lambda _i-g_i) + \delta \boldsymbol {u}_f^T(\Delta \boldsymbol {u}_i + \eta _i\delta \boldsymbol {u}_{ri})]^2 </math>
|-
| style="text-align: center;" | <math> - (\psi ^2 + \delta \boldsymbol {u}_f^T\delta \boldsymbol {u}_f)\times [\psi ^2(\Delta \lambda _i-g_i)^2 - \Delta l^2  </math>
|-
| style="text-align: center;" | <math> + (\Delta \boldsymbol {u}_i + \eta _i\delta \boldsymbol {u}_{ri})^T(\Delta \boldsymbol {u}_i + \eta _i\delta \boldsymbol {u}_{ri})\geq{0}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

After some algebra, we get the following quadratic equations in <math display="inline">\eta _i</math>

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vartheta (\eta _i) = \breve{a}\eta _i^2 + \breve{b}\eta _i + \breve{c} \leq 0                                                                                             </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

where

<span id="eq-26.a"></span>
<span id="eq-26.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\breve{a} =  (\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex} + \delta \boldsymbol {u}^T_{f}\delta \boldsymbol {u}_{f})\delta \boldsymbol {u}_{ri}^T\delta \boldsymbol {u}_{ri} - (\delta \boldsymbol {u}_{f}^T\delta \boldsymbol {u}_{ri})^2   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26.a)
|-
| style="text-align: center;" | <math> \breve{b} =  2(\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex} + \delta \boldsymbol {u}^T_{f}\delta \boldsymbol {u}_{f})(\Delta \boldsymbol {u}_i^T\delta \boldsymbol {u}_{ri}) - 2(\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}(\Delta \lambda _i-g_i) + \delta \boldsymbol {u}^T_{f}\Delta \boldsymbol {u}_i)\delta \boldsymbol {u}_{f}^T\delta \boldsymbol {u}_{ri}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26.b)
|-
| style="text-align: center;" | <math> \breve{c} =  (\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex} + \delta \boldsymbol {u}^T_{f}\delta \boldsymbol {u}_{f})(\Delta p^T_i\Delta p_i-\Delta l^2)  - (2\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}(\Delta \lambda _i-g_i) + \delta \boldsymbol {u}^T_{f}\Delta \boldsymbol {u}_i)\times \delta \boldsymbol {u}^T_{f}\Delta \boldsymbol {u}_i </math>
|-
| style="text-align: center;" | <math>               +\psi ^2\boldsymbol {f}^T_{ex}\boldsymbol {f}_{ex}(\Delta \lambda _i-g_i)^2\delta \boldsymbol {u}^T_{f}\delta \boldsymbol {u}_{f} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26.c)
|}
|}

Note that all this terms can be computed at the <math display="inline">i</math>th iteration in <math display="inline">\delta \lambda _i</math>.  Newly computed the value of <math display="inline">\eta _i</math> and chosen the more convenient value,  is then fed back into original constraint equation which is solved again. Suitable value of this scaling parameter are chosen close to 1. Zhou and Murray (1994) suggested that suitable value of <math display="inline">\eta _i</math> should be chosen from <math display="inline">(0,1)</math>. A suitable value of <math display="inline">\eta _i</math> proposed by Lam and Morley <span id='citeF-1'></span>[[#cite-1|[1]]]  is chosen as below:

<span id="eq-27.a"></span>
<span id="eq-27.b"></span>
<span id="eq-27.c"></span>
<span id="eq-27.d"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\eta _i = \eta _2 - \xi \hbox{ if }\eta _2<1    </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27.a)
|-
| style="text-align: center;" | <math> \eta _i = \eta _2 + \xi \hbox{ if }-\breve{b}/\breve{a}<1<\eta _2 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27.b)
|-
| style="text-align: center;" | <math> \eta _i = \eta _1 - \xi \hbox{ if }\eta _1<1<-\breve{b}/\breve{a} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27.c)
|-
| style="text-align: center;" | <math> \eta _i = \eta _1 + \xi \hbox{ if }1<\eta _1    </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27.d)
|}
|}

where <math display="inline">\eta _1</math> and <math display="inline">\eta _2</math> are roots of <math display="inline">\vartheta (\eta _i)</math>; <math display="inline">\eta _2\geq \eta _1</math> and <math display="inline">\xi = \beta (\eta _2-\eta _1)</math>. Usual, a value of <math display="inline">\beta </math> is take as <math display="inline">0.05</math>. While we found in some certain situation in the practical  calculations, <math display="inline">\eta _i</math> should be chosen carefully to avoid problems of numerical precision <span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]]. The value of parameter is initially set to 1; with the consequent divergence of solution process due to complex. This way the problem of complex roots can be avoided successfully. The authors further report that generally the quadratic equation solved for <math display="inline">\eta </math> gives real values, but in case it lead to complex values, the following strategy explain in Method <math display="inline">2</math> most be adopted.

<span id="fn-1"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) By author's experiences, it necessary to take care about computer precision, due to sometimes the value of the quantities in quadratic equations are really small and can be lost significant digits.</span>

===4.2 Method 2. Elimination of <math>\boldsymbol {r}^{\perp }</math>===

Generally the quadratic equation [[#eq-25|25]] solved for <math display="inline">\eta _i</math> gives real roots. Since it is not possible to find a real roots to give the chosen arc length,  Lam and Morley (1992) suggested that the main problem is that the orthogonal  error <math display="inline">\boldsymbol {r}^{\perp }</math> is too large. In case that Eqn [[#eq-25|25]] leads to complex values,  the follow procedure proposed by Lam and Morley (1992) is adopted <span id='citeF-4'></span>[[#cite-4|[4]]].

<ol>

<li>Compute the displacement due to the orthogonal component of out-of-balance load vector (orthogonal to applied load vector) and update the displacement

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\Delta \boldsymbol {u}^{cr} = \Delta \boldsymbol {u}_i + \delta \boldsymbol {u}_{ri}                                                                                             </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

where superscript <math display="inline">cr</math> is related to the complex roots.
<li>Compute associated nodal forces <math display="inline">\boldsymbol {f}_{in}^{cr}</math>. </li>

<li>Compute incremental loading parameter and the arc-length from orthogonal component of out-of-balance force vector.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\lambda ^{cr}         = \frac{(\boldsymbol {f}_{in}^{cr})^T\boldsymbol {f}_{ex}}{\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.a)
|-
| style="text-align: center;" | <math> \Delta \lambda ^{cr}   = \lambda ^{cr}-\lambda _0  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.b)
|-
| style="text-align: center;" | <math> (\Delta l^{cr})^2  = (\Delta \boldsymbol {u}^{cr})^T\Delta \boldsymbol {u}^{cr} + \psi ^2(\Delta \lambda ^{cr})^2\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}                                                                                              </li>

</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.c)
|}
|}
<li>Reduce the loading parameter by multiplying the loading parameter by the ratio of original arc-length to arc-length associated with complex roots, and iterate the solution process.

<span id="eq-30.a"></span>
<span id="eq-30.b"></span>
<span id="eq-30.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mu                = \frac{\Delta l}{\Delta l^{cr}}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (30.a)
|-
| style="text-align: center;" | <math> \Delta \lambda      = \mu \Delta \lambda ^{cr}          </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (30.b)
|-
| style="text-align: center;" | <math> \Delta \boldsymbol {u}_{i+1} = \mu \Delta l^{cr}               </li>

</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (30.c)
|}
|}

</ol>

Compared with reducing arc-length directly, this method is more targeted to eliminate residuals and hence should be more effective.

===4.3 Method 3. No solution of quadratic equations===

Forde and Stiemer (1987) proposed another arc-length method in which the direct solving of quadratic equation  [[#eq-21|21]] is not required, and thus the problem of dealing with complex roots is avoided. These kind  of arc-length methods are called arc-length orthogonality methods. In the standard incremental iterative solution of equilibrium equation for proportional loading,  the incremental displacement can be written in two components.

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \boldsymbol {u} = \Delta \lambda \Delta \boldsymbol {u}_1 + \Delta \boldsymbol {u}_2                                                                                                                </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

The first is the displacement <math display="inline">\Delta \boldsymbol {u}_1</math> due to a unit load factor multiplied by the incremental variation  in the load level <math display="inline">\Delta \lambda </math> . The second <math display="inline">\Delta \boldsymbol {u}_2</math> is the displacement update for a conventional  load controlled procedure due to the unbalanced loads. Therefore, A general arc length procedure can be derived from orthogonality  principles:

<span id="eq-32"></span>
<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {t}_i = \boldsymbol {u}_i + \beta \lambda _i   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|-
| style="text-align: center;" | <math> \boldsymbol {n}_i = \Delta \boldsymbol {u} + \beta \Delta \lambda  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}
|}

in which <math display="inline">\boldsymbol {t}_i</math> is the tangent of the current incremental load-displacement configuration, while <math display="inline">\boldsymbol {n}_i</math> is an arbitrary update direction chosen with reference to <math display="inline">\boldsymbol {t}_i</math>, see Figure [[#img-3|3]] .

The scalar product of these vectors yields a residual <math display="inline">\boldsymbol {g}_i</math>.  Form the scalar product using equations [[#eq-31|31]],[[#eq-32|32]] and [[#eq-33|33]]  a general expression can be derived for <math display="inline">\Delta \lambda </math>

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \lambda = \frac{\boldsymbol {g}_i-\boldsymbol {u}_i^T\Delta \boldsymbol {u}_2}{\beta ^2\lambda _i + \boldsymbol {u}_i^T\Delta \boldsymbol {u}_1}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|-
| style="text-align: center;" | 
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}
|}

The expression can be simplified for particular cases of orthogonality. The use of the method simplifies the solution process. The method reveals exactly similar results for <math display="inline">\Delta \lambda </math> as obtained by Crisfield (1981) but without solving quadratic equation and selection of proper root <span id='citeF-4'></span>[[#cite-4|[4]]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-fig_3.png|600px|arc-length orthogonality method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' arc-length orthogonality method.
|}

==5 The choice of the appropriate root==

The iterative load factor is normally chosen as the solution to the quadratic equation that yields the minimum angle between <math display="inline">\Delta \boldsymbol {u}_{i-1}</math> and <math display="inline">\Delta \boldsymbol {u}_{i}</math> (Crisfield, 1991),  i.e. <math display="inline">\delta \lambda _i</math> is the solution of [[#eq-21|21]] which gives the maximum product <math display="inline">\Delta \boldsymbol {u}_{i-1}^T\Delta \boldsymbol {u}_{i}</math> <span id='citeF-5'></span>[[#cite-5|[5]]]:

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \lambda _i = arg\left[\max _{\delta \lambda |(\delta \lambda _{i})^2 + b(\delta \lambda _{i}) + c = 0}(\Delta \boldsymbol {u}_{i-1} + \eta \delta \boldsymbol {u}_{ri} + \delta \lambda \boldsymbol {u}_f)^T \Delta \boldsymbol {u}_{k-1} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

==6 Solution formulation at beginning of the iteration ==

The arc-length increment <math display="inline">\Delta l_i</math> for the <math display="inline">i-th</math> loading step is determined by the following procedure:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta l_i = \Delta l_{i-1} \frac{I_d}{I_{i-1}}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|-
| style="text-align: center;" | 
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}
|}

where <math display="inline">\Delta l_{i-1}</math> is the arc-length increment of the previous loading step, <math display="inline">I_d</math> is the desired number of iterations for the <math display="inline">i-th</math> loading step before convergence and <math display="inline">I_{i-1}</math>   is the number of iterations required to converge in the <math display="inline">(i-1)th</math> loading step.

The first iteration from the previous converge point is a predictor solution in which  the appropriate <math display="inline">\Delta \lambda </math> cannot be obtained by solving the quadratic equation showed above.  Lam and Morley (1992) proposed a formation of  and  by introducing a scalar <math display="inline">\mu </math>

<span id="eq-39.a"></span>
<span id="eq-39.b"></span>
<span id="eq-39.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu                  = \frac{\Delta l}{\sqrt{\Delta \boldsymbol {u}_{pr}^T\Delta \boldsymbol {u}_{pr} + \psi ^2\Delta \lambda _{pr}^2}}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.a)
|-
| style="text-align: center;" | <math> \delta \boldsymbol {u}_1       = \mu \Delta \boldsymbol {u}_{pr}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.b)
|-
| style="text-align: center;" | <math> \delta \lambda        = \mu \Delta \boldsymbol {\lambda }_{pr}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.c)
|-
| style="text-align: center;" | 
| style="width: 5px;text-align: right;white-space: nowrap;" | (39.d)
|}
|}

where <math display="inline">()_{pr}</math> corresponds to quantities of previous converged step. Another formation is also widely used in the cylindrical arc-length method  (in which <math display="inline">\psi ^2</math> is set equal to zero),   written as

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu                  = \frac{\Delta l}{\sqrt{\Delta \boldsymbol {u}_{pr}^T\Delta \boldsymbol {u}_{pr}}}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|-
| style="text-align: center;" | 
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}
|}

While in both formation, the sign of the  should be choose carefully to avoid track back  on the current path  (Crisfield, 1991). Many procedures are developed to predict the continuation direction, i.e. to choose the sign of   that carries on tracing the current solution path. Some criteria are listed below:

<ol>

<li>''Stiffness determinant''. Follow the sign of the actual tangent stiffness determinant:

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

sign(\delta \lambda ) = sign(det(\boldsymbol {K}_T(\boldsymbol {u})))   </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}
<li>''Incremental work''. Follow the sign of the predictor work increment:

<span id="eq-43"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

sign(\delta \lambda ) = sign(\delta \boldsymbol {u}^T\boldsymbol {f}_{ex})   </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}
<li>''Secant path'' (Feng et al., 1995, 1996). The sign of <math display="inline">\delta \lambda </math> is determined as:

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

sign(\delta \lambda ) = sign(\Delta \boldsymbol {u}^T\delta \boldsymbol {u})   </li>

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

</ol>

Criteria <math display="inline">(c)</math> are implemented in program NOLM. Procedure (a) is widely used in commercial finite element codes and works well in the absence of bifurcations. In the presence of bifurcations, however, it is known not to be appropriate and fails in most cases. As pointed out by Crisfield (1991), its ill-conditioned behaviour stems from the fact that the sign of <math display="inline">\boldsymbol {K}_T</math> changes either when a limit point or when a bifurcation point is passed. In this case, the predictor cannot distinguish between these two quite different situations, unless further (usually computationally expensive) analyses are undertaken. In the presence of a bifurcation, instead of following the current equilibrium path, the solution will oscillate about the bifurcation point. This property is mathematically proved by Feng et al. (1997). Procedure (b), on the other hand, is blind  to bifurcations and can continue to trace an equilibrium path after passing a bifurcation point. However, this criterion proves ineffective in the descending branch of the load-deflection curve in snap-back  problems, where the predicted positive slope  will provoke a back tracking  load increase.  One important feature shared by the criteria (a) and (b) is the fact that they rely exclusively  on information relative to the current equilibrium point (at the beginning of the increment).  The decision on the sign of <math display="inline">\delta \lambda </math> is made without considering the history of the currently  traced equilibrium path. In situations such as the ones pointed out above, this may result in wrong direction prediction. In contrast, a key point concerning criterion <math display="inline">(c)</math> is the fact that <math display="inline">\Delta \boldsymbol {u}_n</math> carries with  it information about the history of the current equilibrium path <span id='citeF-5'></span>[[#cite-5|[5]]].

===6.1 Summary Arc Length Algorithm Procedure===

A basic arc length algorithm implemented in NOLM program is summarized in the following box depicted in [[#6.1 Summary Arc Length Algorithm Procedure|6.1]]. This algorithm correspond to the spherical arc-length method, however, with small change in the    code, a cylindrical arc-length procedure can also be used .

boxed   BOXHbox BOXBox

<span style="text-align: center; font-size: 75%;"> Arc length scheme.

<ol>

<li> In step <math display="inline">n+1</math>, assuming that we converged in step <math display="inline">n</math>, initialize iteration counter, <math display="inline">i:=0</math>, and set initial  guess for displacement and incremental load factor. </li>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mu = \frac{\Delta l_{n+1}}{\sqrt{\Delta \boldsymbol {u}_{n}^T\Delta \boldsymbol {u}_{n} + \psi ^2\Delta \lambda _{n}^2}}  \Delta \boldsymbol {u}_{n+1}^0 = \Delta \boldsymbol {u}_{n} </math>
|-
| style="text-align: center;" | <math> \lambda _{n+1}^0 = \lambda _{n}  \boldsymbol {r}^0_{n+1} = \boldsymbol {f}_{in}(\boldsymbol {u}_n)-\lambda _{n+1}^0\boldsymbol {f}_{ex}  

</math>
|}
|}
<li> Assemble the tangent stiffness matrix or a projection of it: <math display="inline">\boldsymbol {K}_T = \boldsymbol {K}_T(\boldsymbol {u}_{n+1})</math>. It is important to calculate it because the quadratic convergence  is ensured and also the snap-back point and bifurcation points can be localized.  </li>
<li> Set <math display="inline">i= i + 1</math>. Solve the linear systems for perpendicular solution <math display="inline">\boldsymbol {u}_{ri}</math> and for the tangential solution <math display="inline">\boldsymbol {u}_{rf}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

g_i    = \frac{\boldsymbol {r}_i^T\boldsymbol {f}_{ex}}{\boldsymbol {f}_{ex}^T\boldsymbol {f}_{ex}}  \boldsymbol {r}^\perp = \boldsymbol {r}_i-g_i\boldsymbol {f}_{ex} </math>
|-
| style="text-align: center;" | <math> \boldsymbol {u}_{ri} = \boldsymbol {K}_T^{-1}\boldsymbol {r}^{\perp }  \boldsymbol {u}_{f} = \boldsymbol {K}_T^{-1}\boldsymbol {f}_{ext}    </li>

</math>
|}
|}
<li> Find iterative load factor <math display="inline">\delta \lambda ^i</math>

<ol>

</li>
<li> if <math display="inline">i=1</math> (predictor solution) then compute:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\delta \lambda _{n+1}^1 = sign(\Delta \boldsymbol {u}_n^T\delta \boldsymbol {u}_f)\mu \Delta \lambda _{n}  </li>

</math>
|}
|}
<li> if <math display="inline">i\neq{1}</math> then solve the spherical arc-length constraint equation

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

a(\delta \lambda _{i})^2 + b(\delta \lambda _{i}) + c = 0  </li>

</math>
|}
|}

with coefficients defined by [[#eq-22.a|22.a]], [[#eq-22.b|22.b]] and [[#eq-22.c|22.c]]  and choose root <math display="inline">\delta \lambda _i</math> according to  [[#eq-36|36]] if real roots was found.

<ol>

<li> if complex root was found solve the quadratic equations for <math display="inline">\eta </math> in [[#eq-25|25]] and choose a suitable value for <math display="inline">\eta </math> according to [[#eq-27.a|27.a]], [[#eq-27.b|27.b]] [[#eq-27.c|27.c]],[[#eq-27.d|27.d]]. If real roots was found, come back to b and try again to solve the quadratic equation for <math display="inline">\delta \lambda </math>  in [[#eq-21|21]]. </li>
<li> Is a complex value of <math display="inline">\eta </math> was found apply Method <math display="inline">2</math>, i.e, eliminate <math display="inline">\boldsymbol {r}^{\perp }</math> and update the new value of <math display="inline">\mu </math>, <math display="inline">\Delta \lambda </math> and <math display="inline">\boldsymbol {u}</math> according to [[#eq-30.a|30.a]], [[#eq-30.b|30.b]] and [[#eq-30.c|30.c]].     </li>

</ol>

</ol>

<li> Apply correction to incremental load factor

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\Delta \lambda ^{n+1}_{i+1}  = \Delta \lambda _i - g_i + \delta \lambda _i  </li>

</math>
|}
|}
<li> Compute iterative displacement

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\delta \boldsymbol {u}^{n+1}_{i+1}  = \eta \delta \boldsymbol {u}_{ri} + \delta \lambda _{i}\delta \boldsymbol {u}_{f}  </li>

</math>
|}
|}
<li> Correction to total and incremental displacements

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\Delta \boldsymbol {u}_{i+1}^{n+1}  = \Delta \boldsymbol {u}_{i} + \delta \boldsymbol {u}_{i+1}  \boldsymbol {u}_{i+1}^{n+1}  = \boldsymbol {u}_{i}^{n+1} + \delta \boldsymbol {u}_{i+1}    </li>

</math>
|}
|}
<li> Update residual

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\boldsymbol {r}_{i+1}^{n+1} =  \boldsymbol {f}_{in}(\boldsymbol {u}_{n+1}^{i+1})-\lambda ^{n+1}_{i+1}\boldsymbol {f}_{ex}  </li>

</math>
|}
|}
<li> Check for convergence with a appropriate norm.  

IF <math display="inline">\|  \boldsymbol {r}  \| /\|  \boldsymbol {f}_{ex}  \| <\epsilon _{tol}</math>  THEN set <math display="inline">(\cdot )_{n+1}=(\cdot )^{n+1}_i</math>   and EXIT.

ELSE GO TO (ii). </li>

<li> If not convergence was not reached during the iteration, reduce the <math display="inline">\Delta l_{n+1}</math> with a suitable user factor <math display="inline">\beta </math> and try again.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\Delta l_{n+1} = \beta \Delta l_{n+1}  </li>

</math>
|}
|}

<li> If not convergence was not reached during  the reduction of the <math display="inline">\Delta l_{n+1}</math> a bifurcation analysis must be done, for getting the new equilibrium path (Future work). </li>

</ol>

</span>

==7 Smeared Isotropic Rankine Damage Models==

===7.1 Constitutive Equation===

In the last years, the so-called ''continuum damage model'' have been widely accepted  as an alternative to deal with complex constitutive behavior. This approach has proved simple and versatile, as well a rigorously based in fundamental constitutive theory. The ''continuum damage model'' explains the process of gradual material deterioration induced  by nucleation and grown of multitude of micro-defects. It is based on the definition of effective stress,  which is introduced in connection with the hypothesis strain equivalence. In others words, the strain associated with a damage state under the applied  stress <math display="inline">\boldsymbol {\sigma }</math>  is equivalent to the strain associated  with its undamaged state under effective  stress <math display="inline">\overline{\boldsymbol {\sigma }}</math>. This damage state is monitored trough a single internal scalar variable,  called ''damage'' or ''degradation'' , <math display="inline">d</math>. This variable is a measure  of the lost of secant stiffness of the material, and it ranges from <math display="inline">0</math> for undamaged material to <math display="inline">1</math> for fully degradated one <math display="inline">0 \leq d \leq 1</math>.

In the present work, the effective stresses <math display="inline">\overline{\boldsymbol {\sigma }}</math> are computed in terms of the total strain tensor <math display="inline">\boldsymbol {\varepsilon }</math> as

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\boldsymbol {\sigma }}  = \boldsymbol {C} : \boldsymbol {\varepsilon }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where <math display="inline">\boldsymbol {C}</math> is the secant fourth order constitutive tensor and <math display="inline">(:)</math> denotes double contraction. The constitutive equations for the scalar isotropic damage model used in this work is:

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\sigma }  = (1-d)\overline{\boldsymbol {\sigma }} =(1-d) \boldsymbol {C}_0 : \boldsymbol {\varepsilon }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

where <math display="inline">\boldsymbol {C}_0</math> is the linear elastic constitutive tensor. This equations reveals that the material stiffness  is only affected by a scalar factor and therefore isotropic is preserve and the explicit integration of the constitutive equation is achieved, as only current values of strains and damage.

===7.2 Isotropic damage model ===

The model defined in equation ([[#eq-45|45]]) is fully determined if the value of <math display="inline">d</math> can be evaluated at every time of the deformation process. To that extends one must define a equivalent effective stress  or a suitable norm, <math display="inline">\tau </math>. It is used to compare different states of the stress, so that is possible to define such concepts as  loading , unloading  and reloading , well known as Kuhn-Tucker conditions. The norm is a scalar and positive function and have a zero value for the undeformed state.   In this work, the ''equivalent stress'' will assume the following form:

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau  = \left(\theta + \frac{1-\theta }{n}\sqrt{\boldsymbol {\sigma }_0:\boldsymbol {C}_0^{-1}: \boldsymbol {\sigma }_0}  \right)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

where <math display="inline">n</math> is defined as ratio compressive strength/tensile strength and <math display="inline">\theta </math> is a weighting factor depending on the state of stress <math display="inline">\sigma _0</math>. A satisfactory definition of <math display="inline">\theta </math> is

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\theta  = \frac{\sum _{i=1}^3\langle  \sigma _0^i \rangle }{\sum _{i=1}^3 |   \sigma _0^i  |}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

where <math display="inline">\langle  \cdot  \rangle </math> are the Macaulay brackets <math display="inline">(\langle  x \rangle  = x \hbox{, if } x\geq{0} \hbox{ and } \langle  x \rangle =0\hbox{, if } x<0)</math>. Note that the values of <math display="inline">\theta </math> range from <math display="inline">0</math> for triaxial compression to one for triaxial tension and for intermediate stress a value betweens <math display="inline">0<\theta{<1}</math>.   With the above definitions for the equivalent effective stress, the damage criterion, <math display="inline">\mathcal{F}(\tau , r)\leq{0}</math>, formulated in the undamaged stress spaces, is introduced as:

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{F}(\tau ^t, r^t) = \tau ^t - r^t \leq{0} \hbox{ } \forall t \geq 0     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

where the <math display="inline">\tau ^t</math> is the norm described before and <math display="inline">r^t</math> is the current damage threshold at the time <math display="inline">t</math>. This last value is an internal stress-like variable whose controls the size of damage surface. If <math display="inline">r_0=\sigma _0</math> is the initial threshold value, in this case, the initial uniaxial damage stress(tensile strength),  it must be <math display="inline">r^t\geq r_0=\sigma _0</math>. Damage occurs when the norm <math display="inline">\tau </math> exceeds the current threshold value. Fig [[#img-4|4]] shows the elastic domain of the function defined in Eqn [[#eq-49|49]].

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-fig_4.png|600px|Elastic domain selected for the isotropic damage model.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Elastic domain selected for the isotropic damage model.
|}

The evolution of the damage bonding surface for loading, unloading and reloading conditions is controlled by the Kuhn-Tucker relations and the damage consistency conditions, which are

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{r}\geq 0  \mathcal{F}(\tau , r) \leq{0}  \dot{r}\mathcal{F}(\tau , r) = 0  </math>
|-
| style="text-align: center;" | <math> \hbox{if } \mathcal{F}(\tau , r) = 0 \hbox{ then } \dot{r}\dot{\mathcal{F}}(\tau , r) = 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

leading, in view of equation ([[#eq-49|49]]) , to the loading condition

<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\tau } = \dot{r}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

This, in turn, the evolution of the internal variable <math display="inline">r</math> may be explicitly integrated to render:

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r = \hbox{max}\{ r_0, \hbox{max}(\tau (t))\}  0\leq t< \infty   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

Note that ([[#eq-51|51]]) allows to compute the current values for <math display="inline">r</math> in terms of the current value of <math display="inline">\tau </math>, which depends explicitly on the current total strains. Stress softening is controlled by the  ''softening function'' <math display="inline">q = q(r)</math>. In this work, the following exponential softening law is used:

<span id="eq-53"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q(r) = r_0 e^{(-2H_S(\frac{r-r_0}{r_0}))} \hbox{ , } r\geq r_0    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

where <math display="inline">H_S</math> is the softening parameter which defines the softening response, dependent on the element size and  also sets a maximum size for the elements that can be used in the analysis.  The figure shows a schematic representation of this function.  Finally, the damage index <math display="inline">d = d(r)</math> is explicitly defined in terms of the corresponding current value of damage threshold,

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d(r) = 1-\frac{q(r)}{r} \hbox{ , } r\geq r_0    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

so that it is monotonically increasing function such that <math display="inline">0\leq d \leq 1</math>.

===7.3 Mechanical Dissipation===

The mechanical free energy <math display="inline">\Psi </math> for an isotropic damage model is defined as follows:

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Psi = (1-d)\Psi _0(\boldsymbol {\varepsilon }) = (1-d)[\frac{1}{2}\boldsymbol {\varepsilon }:\boldsymbol {C}_0:\boldsymbol {\varepsilon }] \geq 0     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

Thus, the rate of energy dissipation may be written as

<span id="eq-56"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\mathcal{D}} = \Psi _0 \dot{d} =   \geq 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

Consider now an ''uniaxial'' experiments in which the tensile strain <math display="inline">\varepsilon </math> increases monotonically and quasi-statically from an initial unstressed state to another in which full degradation takes place. In this case, <math display="inline">\Psi _0 = \frac{\sigma _0^2}{2E} = \mathcal{U}_0</math>, where <math display="inline">E</math> is the Young's modulus and <math display="inline">\sigma _0 = E\varepsilon </math>. Therefore, for any process the total specific dissipated energy (per unit volume) is

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{D} = \int _0^\infty \dot{\mathcal{D}} dt =  \int _0^\infty \Psi _0 \dot{d}dt </math>
|-
| style="text-align: center;" | <math> \mathcal{D} = \frac{1}{2E}(\int _{r_0}^\infty q(r)dr-\int _{q_0}^0 rdq)      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

Combining equation ([[#eq-57|57]]) with  ([[#eq-54|54]]) results:

<span id="eq-58"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{D} = (1+\frac{1}{H_S})\frac{\sigma _0^2}{2E} = (1+\frac{1}{H_S})\mathcal{U}_0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

Finally, as it is usual in smeared crack models, it is possible to relate the specific dissipated energy <math display="inline">\mathcal{D}</math>

<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{D}\mathcal{L}_e = \mathcal{G}_f   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

where <math display="inline">\mathcal{G}_f</math> is the mode I fracture energy per unit area (assumed to be a material property)  and <math display="inline">\mathcal{L}_e</math> is the computational width of the fracture zone well knows as ''element characteristic length'', that is computed depending on the geometric dimension of the element. Therefore, the <math display="inline">H_S</math> parameter can be expressed

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>H_S = \frac{\overline{H}_S\mathcal{L}_e}{1-\overline{H}_S\mathcal{L}_e} \geq 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

where the ''specific softening parameter'' <math display="inline">\overline{H}_S = \mathcal{U}_0 \mathcal{G}_f</math> measures the brittleness of the material and <math display="inline">H_S</math>, the ''elemental softening parameter'' measures the  brittleness of the element. For linear simplex element, the characteristic length can be taken as the representative size of the element, <math display="inline">\mathcal{L}_e = l_e</math>. In this work, and assuming that the elements are almost equilateral, the size of the element can be computed <math display="inline">l_e^2 = \sqrt{A_e}</math> for quadrilateral elements  being <math display="inline">A_e</math> the area of the finite element. It is interesting to comment that for a linear irreducible formulation, the discrete localization band <math display="inline">\mathcal{w}</math> is only one element across. The box [[#7.3 Mechanical Dissipation|7.3]] shows the very simple algorithm has to be used to evaluate stresses using this smeared crack models. Note the updating of the threshold value of <math display="inline">\tau ^t</math> as the nonlinear behavior proceeds:

boxed   BOXBox

<span style="text-align: center; font-size: 75%;"> Smeared Isotropic damage model.

<ol>

<li> Initial Data for time <math display="inline">t+1</math>

<ol>

</li>
<li>Material Properties <math display="inline">\sigma _t</math>,<math display="inline">\sigma _c</math>, <math display="inline">n</math>, <math display="inline">E_0</math>, <math display="inline">\nu </math>, <math display="inline">\mathcal{G}_f</math> </li>
<li>Current values: <math display="inline">\boldsymbol {\varepsilon }^{t+1}</math>,<math display="inline">d^t</math> and <math display="inline">r^t</math> </li>

</ol>

<li>Determine <math display="inline">H_S</math>. Eqn [[#eq-60|60]]. </li>
<li>If <math display="inline">t=0</math> the initialize <math display="inline">\tau _0=\sigma _0</math>. </li>
<li>Evaluate undamaged stresses.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\boldsymbol {\sigma }_0^{t+1} = \boldsymbol {C}_0:\boldsymbol {\varepsilon }^{t+1}   </li>

</math>
|}
|}
<li>Evaluate <math display="inline">\tau ^{n+1}</math>. Eqn [[#eq-47|47]].

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\tau  = \left(\theta + \frac{1-\theta }{n}\sqrt{\boldsymbol {\sigma }_0:\boldsymbol {C}_0^{-1}: \boldsymbol {\sigma }_0}  \right) </li>

</math>
|}
|}
<li>Update internal variables

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

r^{n+1} = max(r^t,\tau ^{n+1}) </math>
|-
| style="text-align: center;" | <math> d^{n+1} = d(r^{t+1})            </li>

</math>
|}
|}
<li>Update stresses

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\boldsymbol {\sigma ^{t+1}} = (1-d^{n+1}) \boldsymbol {\sigma _0^{t+1}}  </li>

</math>
|}
|}

</ol>

</span>

==8 Validation Examples==

The application of the arc-length method is presented in this work to the problem of softening material and catch the equilibrium path   of the structure. The arc length method was integrated in NOLM program, together to a smeared isotropic damage model and no linear truss element. We begin to solve a very simple problem using this non linear  truss element as a benchmarking. Results reveals the strong  snap through and snap back in the equilibrium path of this examples.

Other examples concerns to the continuum models discretized in quadrilateral elements using the smeared crack model coded in NOLM. Another constitutive laws, implemented in NOLM is used. However, no mechanical dissipation and objectivity exist, so a very difficult convergence is found is some situations. We proposed an objective mechanical dissipation for future work to take account the dissipated energy, related strongly with the plastic work of the model. For quadratic convergence for the described  isotropic smeared crack model is necessary to apply the tangent stiffness matrix.  However, it is really difficult to compute and other alternatives need to be adopted. The most common way to get a projection of the tangent stiffness matrix is applying the perturbation method. But, in the present work we used the secant stiffness matrix defined as <math display="inline">\boldsymbol {C}_{sec} = (1-d)\boldsymbol {C}_0</math>. This can be effectively used, but not quadratic convergence is guaranteed. In any case is more faster than using the non-damage stiffness matrix. However, in presence of snap back points, sometimes is encountered some difficult to converge.  

It was coded in NOLM subroutines allow connect the results for future reading in Gid, the Pre and post-processing  developed at CIMNE <span id='citeF-6'></span>[[#cite-6|[6]]].

===8.1 Two 2D nonlinear truss element===

The first numerical example consist of two-member nonlinear truss element shown in the figure [[#img-5|5]] and has been widely used as a benchmark problem for comparison of numerical solution algorithm. The following material are assumed: Axial rigidity <math display="inline">EA=1884.694lb</math>, <math display="inline">h=25.874in</math> and <math display="inline">\alpha=63.4</math> degrees.  The same example was tested by Yang and Shieh (1990). Three loading cases will be studied herein. In the first case (Case I), only a vertical load is imposed; secondly (Case II),  the horizontal load is considered as an imperfection, i.e <math display="inline">P_u=0.05P_v</math> and the third one (Case III), the vertical load is treated as an imperfection, i.e <math display="inline">P_v = 0.05P_u</math>. The results obtained by the present method for the three considered cases are shown in  Figs [[#img-6|6]], fig:Sec4-3 and [[#img-9|9]] respectively,  which are according  with solution presented by Pecknold et ''al'' (1985).

Note that, in the first case an snap-through is obtained. As can be seen in the figure fig:Sec4-3, there are four limits points and two snap back points.  The present analysis has demonstrated the self-adaptive capability of this implemented arc-length method. At the third case, the observation for the first load can also be applied here.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_32a.png|600px|Two non-linear truss element. Geometry and boundary conditions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Two non-linear truss element. Geometry and boundary conditions.
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_8a.png|600px|Analysis result for case I.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Analysis result for case I.
|}

<div id='img-7a'></div>
<div id='img-7b'></div>
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_3a.png|600px|Vertical Load P<sub>v</sub> vs horizontal displacements. ]]
|[[Image:draft_Samper_836349612-Fig_4a.png|600px|Vertical Load P<sub>v</sub> vs vertical displacements. ]]
|- style="text-align: center; font-size: 75%;"
| (a) Vertical Load <math>P_v</math> vs horizontal displacements. 
| (b) Vertical Load <math>P_v</math> vs vertical displacements. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' Analysis results for Case II.
|}

<div id='img-8a'></div>
<div id='img-8b'></div>
<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_5a.png|600px|Vertical Load P<sub>u</sub> vs horizontal displacements. ]]
|[[Image:draft_Samper_836349612-Fig_6a.png|600px|Vertical Load P<sub>u</sub> vs vertical displacements. ]]
|- style="text-align: center; font-size: 75%;"
| (a) Vertical Load <math>P_u</math> vs horizontal displacements. 
| (b) Vertical Load <math>P_u</math> vs vertical displacements. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' Analysis results for Case III.
|}

===8.2 A simple <math>2D</math> benchmarking example===

In this example, we test our code of arc-length method with a problem of square non-linear material under tension.  A simple <math display="inline">2D</math> <math display="inline">1\times{1} m\times m</math>  block discretized in <math display="inline">1</math>, <math display="inline">4</math> and <math display="inline">16</math> quadrilateral bilinear element and drown in figure [[#img-9|9]], [[#img-10a|10a]] and [[#img-10b|10b]] respectively, are analyzed both Von-Mises linear/exponential softening model and Drucker-Prager linear/exponential model. A plain strain analysis is performed. Material are assumed and is depicted in table [[#table-1|1]]. The softening modulus <math display="inline">H</math> for linear softening was taken <math display="inline">1MPa</math>.  For exponential softening, we take the following formula for evolution of shear stress

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _s({\kappa }) = \tau _{s0}e^{-h\kappa } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

where <math display="inline">\kappa </math> the equivalent plastic strain and <math display="inline">h</math> is a constant that controls the softening curve, here assumed to be <math display="inline">h=3000</math>. Note that this models, is not objective due to no dissipation of energy is controlled. For define a DP-Model, the parameter <math display="inline">\alpha </math>, that controls the pressure is assumed <math display="inline">\alpha=0.15</math>. Take into account that the objective of this examples is testing the ''AL method'' implemented in NOLM.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 ''Material data of 2D benchmarking example.''
|- style="border-top: 2px solid;"
|  Material Data 
| Value(<math display="inline">N-m</math>) 
|- style="border-top: 2px solid;"
|  Young's Modulus <math display="inline">E</math>              
| 10 <math display="inline">MPa</math>   
|-
| Poisson's ratio <math display="inline">\upsilon </math>       
| 0.33         
|- style="border-bottom: 2px solid;"
| Shear strength <math display="inline">\tau _{s0}</math>      
| 1 <math display="inline">KPa</math>    

|}

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_35a.png|600px|Mesh one.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Mesh <math>one</math>.
|}

<div id='img-10a'></div>
<div id='img-10b'></div>
<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_33a.png|600px|Mesh two. ]]
|[[Image:draft_Samper_836349612-Fig_34a.png|600px|Mesh three. ]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh <math>two</math>. 
| (b) Mesh <math>three</math>. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 10:''' Analysis meshes cases.
|}

The results of the analysis is plotted in figures [[#img-11|11]] and [[#img-12|12]]. Observe that for both material model, identical responses is obtained for the three different meshes when is used the same softening low. However, although not was used a energy dissipated criteria, all element have the same geometric properties making the softening responses objective.

<div id='img-11a'></div>
<div id='img-11b'></div>
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_25a.png|600px|Vertical Load vs vertical displacement using DP-model with exponential softening. ]]
|[[Image:draft_Samper_836349612-Fig_26a.png|600px|Vertical Load vs vertical displacement using DP-model with linear softening. ]]
|- style="text-align: center; font-size: 75%;"
| (a) Vertical Load vs vertical displacement using DP-model with exponential softening. 
| (b) Vertical Load vs vertical displacement using DP-model with linear softening. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 11:''' Analysis cases for DP-Model.
|}

<div id='img-12a'></div>
<div id='img-12b'></div>
<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_35b.png|600px|Vertical Load vs vertical displacement using Von-Mises model with exponential softening. ]]
|[[Image:draft_Samper_836349612-Fig_36b.png|600px|Vertical Load vs vertical displacement using Von-Mises model with linear softening. ]]
|- style="text-align: center; font-size: 75%;"
| (a) Vertical Load vs vertical displacement using Von-Mises model with exponential softening. 
| (b) Vertical Load vs vertical displacement using Von-Mises model with linear softening. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 12:''' Analysis cases for Von Mises -Model.
|}

===8.3 A <math>2D</math> pull test benchmarking example===

In this example, analysis of a cracked flat under surface tension is taken to test our code of arc-length method.  We selected three different sizes of meshes of the same model depicted in figure [[#img-13|13]]  and calculated  with perfect plastic Drucker-Prager model, perfect plastic Von-Mises model and exponential softening Von-Mises model.  We also tested the influence of different softening  parameters with exponential softening Drucker-Prager model.  The same material used in the previous example are assumed again.  A plain strain analysis is performed. The analysis cases are

<ol>

<li> For Drucker-Prager model with exponential softening taken <math display="inline">\alpha=0.15</math> and <math display="inline">h_1=100</math> <math display="inline">h_2=300</math> and <math display="inline">h_3=180</math>. Case I. </li>
<li>For Von-Mises model with perfect plastic. Case II. </li>
<li>For Von-mises model with with exponential softening <math display="inline">h=1000</math>. Case III. </li>

</ol>

The Case I is represented is the figure [[#img-14|14]].  Observe how the ''AL method'' obtain the softening behavior and thre different responses and post-pick load   are obtained when is used different values of <math display="inline">h</math>. Case II represent a perfect plastic  behavior and the analysis results is plotted in fig [[#img-15|15]]. Note how the ''AL method'' implemented obtained the expected perfect plastic behavior.  The  Case III analysis is depicted in figure [[#img-16|16]]. Note how the results is not objective due to no dissipation is controlled. As a future work, is necessary to implement a mesh objective dissipation in the standard model material implemented in NOLM.   

<div id='img-13a'></div>
<div id='img-13b'></div>
<div id='img-13c'></div>
<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_12a.png|600px|Mesh one. 50 nodes and 34 elements.]]
|[[Image:draft_Samper_836349612-Fig_10a.png|600px|Mesh two. 605 nodes and 544 elements.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh one. 50 nodes and 34 elements.
| (b) Mesh two. 605 nodes and 544 elements.
|-
| colspan="2"|[[Image:draft_Samper_836349612-Fig_11a.png|600px|Mesh three. 167 nodes and 136 elements. ]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Mesh three. 167 nodes and 136 elements. 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 13:''' Finite element meshes for analysis.
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_37b.png|600px|Responses of pull test using three different parameter in exponential softening.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Responses of pull test using three different parameter in exponential softening.
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_29a.png|600px|Perfect plastic response for three meshes using Von-Mises model.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Perfect plastic response for three meshes using Von-Mises model.
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_38b.png|600px|Non objective response in Von-Mises material with exponential softening.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Non objective response in Von-Mises material with exponential softening.
|}

===8.4 A <math>2D</math> notched beam===

In this example, we  have selected for analysis a notch beam tested by Petterson(1981) . The major adantages of this experiments are that it has been repeated several  times and that necessary material parameter, such as <math display="inline">\mathcal{G}_f</math>, has been carefully specified. Finite element meshes, dimension and properties are shown in figure [[#img-17|17]]. <math display="inline">702</math> nodes and <math display="inline">612</math> quadrilateral bi-linear element was used in the analysis. We are assuming a plane stress state for this model. The beam is subjected to a point load in the top center, producing a direct traction in the notch. An exponential softening curve was used according to the average element size. The analysis were performed by full arc-length method implemented and the smeared  isotropic crack model described. The figure [[#img-18|18]] compare experimental results and numerical load deflection curves. Note  how the softening behavior can be satisfactorily simulated  together with the maximum load.

<div id='img-17a'></div>
<div id='img-17b'></div>
<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_14a.png|600px|Geometry and material properties.]]
|[[Image:draft_Samper_836349612-Fig_40b.png|600px|Mesh for analysis.]]
|- style="text-align: center; font-size: 75%;"
| (a) Geometry and material properties.
| (b) Mesh for analysis.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 17:''' Finite element meshes for analysis and basic properties.
|}

<div id='img-18a'></div>
<div id='img-18b'></div>
<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_15a.png|600px|Experimental and numerical results reported by Rots et ''al'' (1985).]]
|[[Image:draft_Samper_836349612-Fig_1b.png|600px|Numerical result using smeared crack model with objective dissipation.]]
|- style="text-align: center; font-size: 75%;"
| (a) Experimental and numerical results reported by Rots et ''al'' (1985).
| (b) Numerical result using smeared crack model with objective dissipation.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 18:''' Experimental and numerical result reported and numerical result using smeared crack model with objective dissipation.
|}

===8.5 Double cantilever beam benchmark===

Fig shows  [[#img-19|19]] a  Double cantilever beam  (DBC) tested by Sock, Baron and Fracois(1979). In this work,   we used the smeared crack models using two different fracture energy to guarantee objective dissipation. A total number of <math display="inline"> 608</math> nodes and <math display="inline">548</math> linear quadrilateral elements was used for discretization.    This example is selected  because represent a pure mode I of fracture. However, some simplifications  is done in the model with respect to the original model, due to this example is aimed to test the arc length method.  Material properties is drown in table [[#table-2|2]]. Note that for define the smeared crack model is only required <math display="inline">3</math> parameter. In fig,  the load is plotted against the crack mouth opening displacement (CMOD).  The computation was performed for exponential softening and assuming a plane stress state.Development of the fracture Zone is represented in fig. 4.8  which represent the evolution of damage (fracture process ) For three different process in the calculation.

<div id='img-19a'></div>
<div id='img-19b'></div>
<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_16a.png|600px|Geometry properties.]]
|[[Image:draft_Samper_836349612-Fig_13a.png|600px|Mesh for analysis.]]
|- style="text-align: center; font-size: 75%;"
| (a) Geometry properties.
| (b) Mesh for analysis.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 19:''' Finite element meshes for analysis and basic properties.
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 ''Material data of 2D benchmarking example.''
|- style="border-top: 2px solid;"
|  Material Data 
| Value(<math display="inline">N-m</math>) 
|- style="border-top: 2px solid;"
|  Young's Modulus <math display="inline">E</math>                   
| 40 <math display="inline">GPa</math>   
|-
| Poisson's ratio <math display="inline">\upsilon </math>            
| 0.2           
|-
| Tensile strength <math display="inline">\sigma _0</math>            
| 4 <math display="inline">MPa</math>    
|-
| Fracture Energy <math display="inline">\mathcal{G}_{f1}</math>    
| 350 <math display="inline">N/m</math>    
|- style="border-bottom: 2px solid;"
| Fracture Energy <math display="inline">\mathcal{G}_{f2}</math>    
| 250 <math display="inline">N/m</math>    

|}

Due to the boundaries condition apply, the forces tends to open the beam in front of the notch and all damage process begin in this place. Damage progresses gradually as the loading process increases  and finally the beam is separated in two part. The evolution of damage to the end  of simulation is potted in [[#img-21b|21b]] . Note that when the damage reaches the value of <math display="inline">1</math>, means that this zone is totally ''smeared'' cracked.

<div id='img-20a'></div>
<div id='img-20b'></div>
<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_17a.png|600px|Experimental and numerical results reported by Rots et ''al'' (1985.]]
|[[Image:draft_Samper_836349612-Fig_9a.png|600px|Numerical result using smeared crack model with objective dissipation.]]
|- style="text-align: center; font-size: 75%;"
| (a) Experimental and numerical results reported by Rots et ''al'' (1985.
| (b) Numerical result using smeared crack model with objective dissipation.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 20:''' Experimental and numerical result reported for DCB and numerical result using smeared crack model with objective dissipation.
|}

<div id='img-21a'></div>
<div id='img-21b'></div>
<div id='img-21c'></div>
<div id='img-21d'></div>
<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_836349612-Fig_18a.png|600px|Contour Sxx Cauchy Stress field at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).]]
|[[Image:draft_Samper_836349612-Fig_19a.png|600px|Evolution of damage for \mathcalG<sub>f1</sub> at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).]]
|- style="text-align: center; font-size: 75%;"
| (a) Contour Sxx Cauchy Stress field at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).
| (b) Evolution of damage for <math>\mathcal{G}_{f1}</math> at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).
|-
|[[Image:draft_Samper_836349612-Fig_41b.png|600px|Direccion of principal stress at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).]]
|[[Image:draft_Samper_836349612-Fig_42b.png|600px|Deformed mesh (×500). (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).]]
|- style="text-align: center; font-size: 75%;"
| (c) Direccion of principal stress at the end of simulation. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).
| (d) Deformed mesh <math>(\times{500)}</math>. (Gid pre-post processor<span id='citeF-6'></span>[[#cite-6|[6]]]).
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 21:''' Double cantilever beam analysis.
|}

==9 Conclusion, recommendations, actual and future works==

===9.1 Conlcusion===

In this work we introduced  and implemented the ''AL method'' in NOLM program. Basically  the actual implementation in NOLM is based in the work proposed by  Lam and Morley (1992) and was coded with take into account the sign for no ''tracking back'' the solution. An efficient method was included for avoid complex roots and if the problem persist, two method is available to get the convergence. In case of no convergence is reached, we introduce  a user factor, for reducing the arc-length, that should be taken betweens <math display="inline">(0.5-0.75)</math> per cent of the actual arc length.

A linear truss-element, together to the isotropic smeared crack model was introduced for testing and benchmarking the ''AL method'' . During the work, we found several difficulties to converge in the standard constitutive law implemented  in NOLM when a softening parameter was introduced. It is caused by the no mechanical  dissipation in the model. However, for a perfect plastic an hardening they get the corresponding results.

All meshing process and post-processing was made with Gid<span id='citeF-6'></span>[[#cite-6|[6]]], proving that it is possible to connect the developed owner software and this Pre-and post processing program.

===9.2 Recommendations===

By author experience, it is recommended, as a future work, the redesign of the basis of  the NOLM code. Difficulties are found in debugging process when several static variables, jumping ( GO TO sentences) and writing process was used. Some bugs was fixed, one of them related to the solver and solution flag processes.

Its really necessary the checking and testing of the constitutive models implemented in NOLM with its corresponding benchmarking. Also recommended the change of the data base of NOLM  (Common Block A), due to lost of memory and a lot of them are not using in a real run time calculations and the difficulties that could be encountered  for parallelization.

The storage of the global stiffness matrix can be done in a parallel way.  However, the needed to know the exactly position of its terms make difficult a future parallelization. In other words, not only be considered the parallelization of the solver using MKL Pardiso, but also all looping in elements, conditions and nodes. A suitable way to reaches this objectives is a modular or object programing inside NOLM, now available in the new INTEL COMPILER versions.

===9.3 Actual and future works===

The actual and future works are aimed to checking standard constitutive laws integrated in NOLM and coded and objective dissipation in the finite element.  This can be done with a ''small'' relative change in the softening parameter curve adding the characterized length of the elements.

Integrating an analysis of bifurcations points when convergence can not reaches using  the ''AL method'' for stopping  or tracing another possible equilibrium path for the next iteration step. It is done by checking the determinant of the global stiffness matrix using the LU decomposition.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Lam, W. and Morley, C. (1992) "Arc‐Length Method for Passing Limit Points in Structural Calculation", Volume 118. Journal of Structural Engineering 1 169-185

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Cinthia A. G. Sousa and Paulo M. Pimenta. (2010) "A new parameter to Arc-Length method in non-linear structural analysis", Volume XXIX. Mecánica Computacional 1841-1848

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' J.V. Ferreira and A.L. Serpa. (2005) "Application of the arc-length method in nonlinear frequency response", Volume 284. Journal of Sound and Vibration 1&#8211;2 133&#8211;149

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' MEMON Bashir-Ahmed, SU Xiao-zu. (2004) "Arc-length technique for nonlinear finite element analysis.", Volume 5. J Zhejiang Univ Sci 5 618-28

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' de Souza Neto, E. A. and Peri, D. and Owen, D. R. J. (2008) "Computational Methods for Plasticity". John Wiley & Sons, Ltd 1&#8211;15

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' GiD. (2009) "The Personal Pre and Post Processor"

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' M. A. Crisfield. (1986) "Snap-through and snap-back response in concrete structures and the dangers of under-integration", Volume 22. International Journal for Numerical Methods in Engineering 3 751-767

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Kyoung Soo Lee and Sang Eul Han and Taehyo Park. (2011) "A simple explicit arc-length method using the dynamic relaxation method with kinetic damping", Volume 89. Computers & Structures  1&#8211;2 216&#8211;233