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 1 where appears the critical points ${\textstyle B}$, ${\textstyle C}$, ${\textstyle F}$, ${\textstyle H}$ and ${\textstyle K}$ 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 ${\textstyle B}$ and ${\textstyle F}$ of Figure 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 ${\textstyle G}$ and ${\textstyle H}$ of Figure 1).

 Figure 1: snap-back and snap through response [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)[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[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 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 [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[3]. In the non-linear finite element analysis, the system of non-linear algebraic equations to be solve take the form:

 ${\displaystyle {\boldsymbol {r}}({\boldsymbol {u}},\lambda )={\boldsymbol {f}}_{int}({\boldsymbol {u}})-\lambda {\boldsymbol {f}}_{ex}}$
(1)

where ${\textstyle {\boldsymbol {r}}={\boldsymbol {r}}({\boldsymbol {u}},\lambda )}$ is a function of the displacement vector ${\textstyle {\boldsymbol {u}}}$ and the load parameter ${\textstyle \lambda }$, and is the out-of-balance force vector which vanishes at equilibrium. The internal force vector ${\textstyle {\boldsymbol {f}}_{in}({\boldsymbol {u}})}$ is a function of the displacements and ${\textstyle {\boldsymbol {f}}_{ex}}$ is the prescribed external force vector and is scaled by the load intensity parameter ${\textstyle \lambda }$. Hence ${\textstyle \lambda {\boldsymbol {f}}_{ex}}$ represents the load applied at the current load increment. Basically, the arc-length method first consider the load factor ${\textstyle \lambda }$ in Eq.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 ${\textstyle \lambda }$ 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

 ${\displaystyle s=\int {\sqrt {d{\boldsymbol {u}}^{T}d{\boldsymbol {u}}+\psi ^{2}d\lambda ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}}}$
(2)

or increment form as

 ${\displaystyle \Psi ={\boldsymbol {\Delta u}}^{T}{\boldsymbol {\Delta u}}+\psi ^{2}\Delta \lambda ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}-\Delta l^{2}=0}$
(3)

where ${\textstyle {\boldsymbol {\Delta u}}}$ is vector of incremental displacement, ${\textstyle {\boldsymbol {f}}_{ex}}$ is specified external load vector, ${\textstyle \Delta \lambda }$ is incremental load-factor, ${\textstyle \Delta l}$ is fixed radius of desired intersection, and ${\textstyle \psi }$ 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 ${\textstyle \Delta l}$. Another constraint equations can also be found in the literature, e.g, Bathe and Dvorkin (1983).

Assuming that displacement increments ${\textstyle {\boldsymbol {\Delta u}}}$ satisfying constraint Eq 3, is possible to obtain a simplified linearized version of this constraint equation adopting the idea of Riks and Wempner(1971) for the ${\textstyle (i+1)th}$ iteration, i.e

 ${\displaystyle {\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}$
(4)

where ${\textstyle \delta {\boldsymbol {u}}_{i}}$ and ${\textstyle \delta \lambda _{i}}$ are the required corrections to ${\textstyle \Delta {\boldsymbol {u}}_{i}}$ and ${\textstyle \Delta \lambda _{i}}$ 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[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 [4].

Therefore, we can directly used a matrix form (Crisfield 1991), that with some simplifications in Eq 1 and 3, we compute the iterative change in displacement vector and load-factor in the ${\textstyle (i+1)th}$ iteration,

 ${\displaystyle \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}}_{i}^{T}&2\Delta \lambda _{i}\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}\end{array}}\right]_{i}^{-1}\left[{\begin{array}{c}{\boldsymbol {r}}\\\Psi \end{array}}\right]_{i}}$
(8)

where we define the effective general matrix of ${\textstyle n+1}$ dimension

 ${\displaystyle {\boldsymbol {\hat {K}}}_{T}=\left[{\begin{array}{cc}{\boldsymbol {K}}_{t}&-{\boldsymbol {f}}_{ex}\\2\Delta {\boldsymbol {p}}_{i}^{T}&2\Delta \lambda _{i}\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}\end{array}}\right]}$
(10)

where ${\textstyle {\boldsymbol {K}}_{T}}$ is the tangential stiffness matrix, and ${\textstyle {\boldsymbol {r}}_{i}}$ and ${\textstyle \Psi _{i}}$ are the previous values of out-of-balance load vector and arc-length. Note that ${\textstyle {\boldsymbol {\hat {K}}}_{T}}$ 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 ${\textstyle \delta {\boldsymbol {u}}_{i+1}}$ and ${\textstyle \delta \lambda _{i+1}}$ have been computed, the displacement vector and load-factor are updated, i.e

 ${\displaystyle \Delta {\boldsymbol {u}}_{i+1}=\Delta {\boldsymbol {u}}_{i}+\delta {\boldsymbol {u}}_{i+1}}$ (11) ${\displaystyle \Delta \lambda _{i+1}=\Delta \lambda _{i}+\delta \lambda _{i+1}}$ (12)

Alternatively, instead of solving directly Eq 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 3. According to this technique, the iterative change of displacement for the new unknown load level ${\textstyle \Delta \lambda _{i+1}=\Delta \lambda _{i}+\delta \lambda _{i+1}}$ is written as [4]

 ${\displaystyle \delta {\boldsymbol {u}}_{i+1}={\boldsymbol {K}}_{T}^{-1}{\boldsymbol {r}}_{i}+\delta \lambda _{i}{\boldsymbol {K}}_{T}^{-1}{\boldsymbol {f}}_{ex}}$ (13) ${\displaystyle =\delta {\boldsymbol {u}}_{ri}+\delta \lambda _{i}\delta {\boldsymbol {u}}_{f}}$ (14)

where ${\textstyle \delta {\boldsymbol {u}}_{ri}}$ is the displacement vector resulting from eliminating the out-of-balance forces and ${\textstyle \delta {\boldsymbol {u}}_{f}}$ is the parallel displacement vector computed and stored after every update of the stiffness matrix ${\textstyle {\boldsymbol {K}}_{T}}$. Substituting values from 11 and 14 into the constraint equation yields in a quadratic equations in ${\textstyle \delta \lambda _{i}}$ so that

 ${\displaystyle a\delta \lambda _{i}^{2}+b\delta \lambda _{i}+c=0}$
(15)

where

 ${\displaystyle a=\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f}+\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (16.a) ${\displaystyle b=2\delta {\boldsymbol {u}}_{f}(\Delta {\boldsymbol {u}}_{i}+\delta {\boldsymbol {u}}_{ri})+2\Delta \lambda _{i}\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (16.b) ${\displaystyle 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}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (16.c)

This method is graphically shown in 2. Note that the parameter ${\textstyle \psi }$ convert everything to the same dimension and order of magnitude. Many analyst set this ${\textstyle \psi =0}$. This procedure is known as the indirect spherical AL method Crisfield(1891). A typical value for ${\textstyle \psi ^{2}}$ is taken as ${\textstyle {\boldsymbol {u}}_{0}^{T}{\boldsymbol {u}}_{0}^{T}/(\lambda _{0}^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}^{T})}$ [1]. Eq.15 is solved to get the value of ${\textstyle \delta \lambda }$, and so to define the iterative change completely. This equation leads to two results of ${\textstyle \delta \lambda }$. 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 ${\textstyle \Delta l}$ chosen.

 Figure 2: Arc-length procedure for specific iteration [8].

3 Modify arc-length procedure

Lam and Morley suggested decomposing the out-of-balance load ${\textstyle {\boldsymbol {r}}}$ into a parallel component ${\textstyle g_{i}{\boldsymbol {f}}_{ex}}$, which is in the specified applied load direction; and orthogonal component ${\textstyle {\boldsymbol {r}}^{\perp }=P^{\perp }(g_{i}{\boldsymbol {f}}_{ex})}$, with respect to fixed nodal load vector in case of divergence, where

 ${\displaystyle g_{i}={\frac {{\boldsymbol {r}}_{i}^{T}{\boldsymbol {f}}_{ex}}{{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}}}$ (17) ${\displaystyle {\boldsymbol {r}}^{\perp }={\boldsymbol {r}}_{i}-g_{i}{\boldsymbol {f}}_{ex}}$ (18)

In the equations above the, ${\textstyle P^{\perp }(\cdot )}$ 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 ${\textstyle \Delta {\boldsymbol {u}}_{i}}$ and the ${\textstyle \Delta \lambda _{i}}$, ${\textstyle \Delta {\boldsymbol {u}}_{i+1}}$ and ${\textstyle \delta \lambda _{i+1}}$ have the formation represented as

 ${\displaystyle \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 }}$ (19.a) ${\displaystyle =\Delta {\boldsymbol {u}}_{i}+\delta \lambda _{i}\delta {\boldsymbol {u}}_{f}+\eta _{i}\delta {\boldsymbol {u}}_{ri}}$ (19.b)

and

 ${\displaystyle \Delta \lambda _{i+1}=\Delta \lambda _{i}+g_{i}+\delta \lambda _{i}}$
(20)

where was introduced a scaling parameter ${\textstyle \eta }$ in Eqn 19.a which will be discussing later. Substituting Eqns 20 and 19.b into the constraint equations 3, gives a new quadratic equations in ${\textstyle \Delta \lambda _{i}}$, similar to Eqn 15 written as

 ${\displaystyle a(\delta \lambda _{i})^{2}+b(\delta \lambda _{i})+c=0}$
(21)

where, in this case

 ${\displaystyle a=\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f}+\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (22.a) ${\displaystyle 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}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (22.b) ${\displaystyle 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}}_{ex}^{T}{\boldsymbol {f}}_{ex}}$ (22.c)

This procedure becomes general arc-length procedure proposed by Crisfield when ${\textstyle \eta =1}$, ${\textstyle \psi ^{2}=0}$ and ${\textstyle g_{i}=0}$. 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 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 ${\textstyle \delta \lambda }$ without resolving any quadratic equations and avoiding the trace-back solution path when a wrong solution of ${\textstyle \delta \lambda }$ is found and choose.

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 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 14 for by a factor ${\textstyle \eta \neq {1}}$, 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 ${\textstyle \eta }$ in the usual constraint equation, that is, however, computed without introducing energy considerations. According to Eqn 21, in order to get the real root, the parameter ${\textstyle a}$, ${\textstyle b}$ and ${\textstyle c}$ in Eqn 22.a, 22.b and 22.c respectively have to satisfy that

 ${\displaystyle b^{2}-ac\geq 0}$
(23)

 ${\displaystyle \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}}$ ${\displaystyle -(\psi ^{2}+\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f})\times [\psi ^{2}(\Delta \lambda _{i}-g_{i})^{2}-\Delta l^{2}}$ ${\displaystyle +(\Delta {\boldsymbol {u}}_{i}+\eta _{i}\delta {\boldsymbol {u}}_{ri})^{T}(\Delta {\boldsymbol {u}}_{i}+\eta _{i}\delta {\boldsymbol {u}}_{ri})\geq {0}}$
(24)

After some algebra, we get the following quadratic equations in ${\textstyle \eta _{i}}$

 ${\displaystyle \vartheta (\eta _{i})={\breve {a}}\eta _{i}^{2}+{\breve {b}}\eta _{i}+{\breve {c}}\leq 0}$
(25)

where

 ${\displaystyle {\breve {a}}=(\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}+\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f})\delta {\boldsymbol {u}}_{ri}^{T}\delta {\boldsymbol {u}}_{ri}-(\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{ri})^{2}}$ (26.a) ${\displaystyle {\breve {b}}=2(\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}+\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f})(\Delta {\boldsymbol {u}}_{i}^{T}\delta {\boldsymbol {u}}_{ri})-2(\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}(\Delta \lambda _{i}-g_{i})+\delta {\boldsymbol {u}}_{f}^{T}\Delta {\boldsymbol {u}}_{i})\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{ri}}$ (26.b) ${\displaystyle {\breve {c}}=(\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}+\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f})(\Delta p_{i}^{T}\Delta p_{i}-\Delta l^{2})-(2\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}(\Delta \lambda _{i}-g_{i})+\delta {\boldsymbol {u}}_{f}^{T}\Delta {\boldsymbol {u}}_{i})\times \delta {\boldsymbol {u}}_{f}^{T}\Delta {\boldsymbol {u}}_{i}}$ ${\displaystyle +\psi ^{2}{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}(\Delta \lambda _{i}-g_{i})^{2}\delta {\boldsymbol {u}}_{f}^{T}\delta {\boldsymbol {u}}_{f}}$ (26.c)

Note that all this terms can be computed at the ${\textstyle i}$th iteration in ${\textstyle \delta \lambda _{i}}$. Newly computed the value of ${\textstyle \eta _{i}}$ 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 ${\textstyle \eta _{i}}$ should be chosen from ${\textstyle (0,1)}$. A suitable value of ${\textstyle \eta _{i}}$ proposed by Lam and Morley [1] is chosen as below:

 ${\displaystyle \eta _{i}=\eta _{2}-\xi {\hbox{ if }}\eta _{2}<1}$ (27.a) ${\displaystyle \eta _{i}=\eta _{2}+\xi {\hbox{ if }}-{\breve {b}}/{\breve {a}}<1<\eta _{2}}$ (27.b) ${\displaystyle \eta _{i}=\eta _{1}-\xi {\hbox{ if }}\eta _{1}<1<-{\breve {b}}/{\breve {a}}}$ (27.c) ${\displaystyle \eta _{i}=\eta _{1}+\xi {\hbox{ if }}1<\eta _{1}}$ (27.d)

where ${\textstyle \eta _{1}}$ and ${\textstyle \eta _{2}}$ are roots of ${\textstyle \vartheta (\eta _{i})}$; ${\textstyle \eta _{2}\geq \eta _{1}}$ and ${\textstyle \xi =\beta (\eta _{2}-\eta _{1})}$. Usual, a value of ${\textstyle \beta }$ is take as ${\textstyle 0.05}$. While we found in some certain situation in the practical calculations, ${\textstyle \eta _{i}}$ should be chosen carefully to avoid problems of numerical precision 1. 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 ${\textstyle \eta }$ gives real values, but in case it lead to complex values, the following strategy explain in Method ${\textstyle 2}$ most be adopted.

(1) 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.

4.2 Method 2. Elimination of r⊥

Generally the quadratic equation 25 solved for ${\textstyle \eta _{i}}$ 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 ${\textstyle {\boldsymbol {r}}^{\perp }}$ is too large. In case that Eqn 25 leads to complex values, the follow procedure proposed by Lam and Morley (1992) is adopted [4].

1. Compute the displacement due to the orthogonal component of out-of-balance load vector (orthogonal to applied load vector) and update the displacement
 ${\displaystyle \Delta {\boldsymbol {u}}^{cr}=\Delta {\boldsymbol {u}}_{i}+\delta {\boldsymbol {u}}_{ri}}$
(28)

where superscript ${\textstyle cr}$ is related to the complex roots.

2. Compute associated nodal forces ${\textstyle {\boldsymbol {f}}_{in}^{cr}}$.
3. Compute incremental loading parameter and the arc-length from orthogonal component of out-of-balance force vector.
 ${\displaystyle \lambda ^{cr}={\frac {({\boldsymbol {f}}_{in}^{cr})^{T}{\boldsymbol {f}}_{ex}}{{\boldsymbol {f}}_{ex}^{T}{\boldsymbol {f}}_{ex}}}}$ (29.a) ${\displaystyle \Delta \lambda ^{cr}=\lambda ^{cr}-\lambda _{0}}$ (29.b) ${\displaystyle (\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}}$ (29.c)
4. 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.
 ${\displaystyle \mu ={\frac {\Delta l}{\Delta l^{cr}}}}$ (30.a) ${\displaystyle \Delta \lambda =\mu \Delta \lambda ^{cr}}$ (30.b) ${\displaystyle \Delta {\boldsymbol {u}}_{i+1}=\mu \Delta l^{cr}}$ (30.c)

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 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.

 ${\displaystyle \Delta {\boldsymbol {u}}=\Delta \lambda \Delta {\boldsymbol {u}}_{1}+\Delta {\boldsymbol {u}}_{2}}$
(31)

The first is the displacement ${\textstyle \Delta {\boldsymbol {u}}_{1}}$ due to a unit load factor multiplied by the incremental variation in the load level ${\textstyle \Delta \lambda }$ . The second ${\textstyle \Delta {\boldsymbol {u}}_{2}}$ 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:

 ${\displaystyle {\boldsymbol {t}}_{i}={\boldsymbol {u}}_{i}+\beta \lambda _{i}}$ (32) ${\displaystyle {\boldsymbol {n}}_{i}=\Delta {\boldsymbol {u}}+\beta \Delta \lambda }$ (33)

in which ${\textstyle {\boldsymbol {t}}_{i}}$ is the tangent of the current incremental load-displacement configuration, while ${\textstyle {\boldsymbol {n}}_{i}}$ is an arbitrary update direction chosen with reference to ${\textstyle {\boldsymbol {t}}_{i}}$, see Figure 3 .

The scalar product of these vectors yields a residual ${\textstyle {\boldsymbol {g}}_{i}}$. Form the scalar product using equations 31,32 and 33 a general expression can be derived for ${\textstyle \Delta \lambda }$

 ${\displaystyle \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}}}}$ (34) (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 ${\textstyle \Delta \lambda }$ as obtained by Crisfield (1981) but without solving quadratic equation and selection of proper root [4].

 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 ${\textstyle \Delta {\boldsymbol {u}}_{i-1}}$ and ${\textstyle \Delta {\boldsymbol {u}}_{i}}$ (Crisfield, 1991), i.e. ${\textstyle \delta \lambda _{i}}$ is the solution of 21 which gives the maximum product ${\textstyle \Delta {\boldsymbol {u}}_{i-1}^{T}\Delta {\boldsymbol {u}}_{i}}$ [5]:

 ${\displaystyle \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]}$
(36)

6 Solution formulation at beginning of the iteration

The arc-length increment ${\textstyle \Delta l_{i}}$ for the ${\textstyle i-th}$ loading step is determined by the following procedure:

 ${\displaystyle \Delta l_{i}=\Delta l_{i-1}{\frac {I_{d}}{I_{i-1}}}}$ (37) (38)

where ${\textstyle \Delta l_{i-1}}$ is the arc-length increment of the previous loading step, ${\textstyle I_{d}}$ is the desired number of iterations for the ${\textstyle i-th}$ loading step before convergence and ${\textstyle I_{i-1}}$ is the number of iterations required to converge in the ${\textstyle (i-1)th}$ loading step.

The first iteration from the previous converge point is a predictor solution in which the appropriate ${\textstyle \Delta \lambda }$ cannot be obtained by solving the quadratic equation showed above. Lam and Morley (1992) proposed a formation of and by introducing a scalar ${\textstyle \mu }$

 ${\displaystyle \mu ={\frac {\Delta l}{\sqrt {\Delta {\boldsymbol {u}}_{pr}^{T}\Delta {\boldsymbol {u}}_{pr}+\psi ^{2}\Delta \lambda _{pr}^{2}}}}}$ (39.a) ${\displaystyle \delta {\boldsymbol {u}}_{1}=\mu \Delta {\boldsymbol {u}}_{pr}}$ (39.b) ${\displaystyle \delta \lambda =\mu \Delta {\boldsymbol {\lambda }}_{pr}}$ (39.c) (39.d)

where ${\textstyle ()_{pr}}$ corresponds to quantities of previous converged step. Another formation is also widely used in the cylindrical arc-length method (in which ${\textstyle \psi ^{2}}$ is set equal to zero), written as

 ${\displaystyle \mu ={\frac {\Delta l}{\sqrt {\Delta {\boldsymbol {u}}_{pr}^{T}\Delta {\boldsymbol {u}}_{pr}}}}}$ (40) (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:

1. Stiffness determinant. Follow the sign of the actual tangent stiffness determinant:
 ${\displaystyle sign(\delta \lambda )=sign(det({\boldsymbol {K}}_{T}({\boldsymbol {u}})))}$
(42)
2. Incremental work. Follow the sign of the predictor work increment:
 ${\displaystyle sign(\delta \lambda )=sign(\delta {\boldsymbol {u}}^{T}{\boldsymbol {f}}_{ex})}$
(43)
3. Secant path (Feng et al., 1995, 1996). The sign of ${\textstyle \delta \lambda }$ is determined as:
 ${\displaystyle sign(\delta \lambda )=sign(\Delta {\boldsymbol {u}}^{T}\delta {\boldsymbol {u}})}$
(44)

Criteria ${\textstyle (c)}$ 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 ${\textstyle {\boldsymbol {K}}_{T}}$ 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 ${\textstyle \delta \lambda }$ 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 ${\textstyle (c)}$ is the fact that ${\textstyle \Delta {\boldsymbol {u}}_{n}}$ carries with it information about the history of the current equilibrium path [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. 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 .

Arc length scheme.

1. In step ${\textstyle n+1}$, assuming that we converged in step ${\textstyle n}$, initialize iteration counter, ${\textstyle i:=0}$, and set initial guess for displacement and incremental load factor.
2.  ${\displaystyle \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}}$ ${\displaystyle \lambda _{n+1}^{0}=\lambda _{n}{\boldsymbol {r}}_{n+1}^{0}={\boldsymbol {f}}_{in}({\boldsymbol {u}}_{n})-\lambda _{n+1}^{0}{\boldsymbol {f}}_{ex}}$
3. Assemble the tangent stiffness matrix or a projection of it: ${\textstyle {\boldsymbol {K}}_{T}={\boldsymbol {K}}_{T}({\boldsymbol {u}}_{n+1})}$. It is important to calculate it because the quadratic convergence is ensured and also the snap-back point and bifurcation points can be localized.
4. Set ${\textstyle i=i+1}$. Solve the linear systems for perpendicular solution ${\textstyle {\boldsymbol {u}}_{ri}}$ and for the tangential solution ${\textstyle {\boldsymbol {u}}_{rf}}$:
 ${\displaystyle 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}}$ ${\displaystyle {\boldsymbol {u}}_{ri}={\boldsymbol {K}}_{T}^{-1}{\boldsymbol {r}}^{\perp }{\boldsymbol {u}}_{f}={\boldsymbol {K}}_{T}^{-1}{\boldsymbol {f}}_{ext}}$
5. Find iterative load factor ${\textstyle \delta \lambda ^{i}}$
1. if ${\textstyle i=1}$ (predictor solution) then compute:
 ${\displaystyle \delta \lambda _{n+1}^{1}=sign(\Delta {\boldsymbol {u}}_{n}^{T}\delta {\boldsymbol {u}}_{f})\mu \Delta \lambda _{n}}$
2. if ${\textstyle i\neq {1}}$ then solve the spherical arc-length constraint equation
 ${\displaystyle a(\delta \lambda _{i})^{2}+b(\delta \lambda _{i})+c=0}$

with coefficients defined by 22.a, 22.b and 22.c and choose root ${\textstyle \delta \lambda _{i}}$ according to 36 if real roots was found.

1. if complex root was found solve the quadratic equations for ${\textstyle \eta }$ in 25 and choose a suitable value for ${\textstyle \eta }$ according to 27.a, 27.b 27.c,27.d. If real roots was found, come back to b and try again to solve the quadratic equation for ${\textstyle \delta \lambda }$ in 21.
2. Is a complex value of ${\textstyle \eta }$ was found apply Method ${\textstyle 2}$, i.e, eliminate ${\textstyle {\boldsymbol {r}}^{\perp }}$ and update the new value of ${\textstyle \mu }$, ${\textstyle \Delta \lambda }$ and ${\textstyle {\boldsymbol {u}}}$ according to 30.a, 30.b and 30.c.
6. Apply correction to incremental load factor
 ${\displaystyle \Delta \lambda _{i+1}^{n+1}=\Delta \lambda _{i}-g_{i}+\delta \lambda _{i}}$
7. Compute iterative displacement
 ${\displaystyle \delta {\boldsymbol {u}}_{i+1}^{n+1}=\eta \delta {\boldsymbol {u}}_{ri}+\delta \lambda _{i}\delta {\boldsymbol {u}}_{f}}$
8. Correction to total and incremental displacements
 ${\displaystyle \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}}$
9. Update residual
 ${\displaystyle {\boldsymbol {r}}_{i+1}^{n+1}={\boldsymbol {f}}_{in}({\boldsymbol {u}}_{n+1}^{i+1})-\lambda _{i+1}^{n+1}{\boldsymbol {f}}_{ex}}$
10. Check for convergence with a appropriate norm. IF ${\textstyle \|{\boldsymbol {r}}\|/\|{\boldsymbol {f}}_{ex}\|<\epsilon _{tol}}$ THEN set ${\textstyle (\cdot )_{n+1}=(\cdot )_{i}^{n+1}}$ and EXIT. ELSE GO TO (ii).
11. If not convergence was not reached during the iteration, reduce the ${\textstyle \Delta l_{n+1}}$ with a suitable user factor ${\textstyle \beta }$ and try again.
 ${\displaystyle \Delta l_{n+1}=\beta \Delta l_{n+1}}$
12. If not convergence was not reached during the reduction of the ${\textstyle \Delta l_{n+1}}$ a bifurcation analysis must be done, for getting the new equilibrium path (Future work).

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 ${\textstyle {\boldsymbol {\sigma }}}$ is equivalent to the strain associated with its undamaged state under effective stress ${\textstyle {\overline {\boldsymbol {\sigma }}}}$. This damage state is monitored trough a single internal scalar variable, called damage or degradation , ${\textstyle d}$. This variable is a measure of the lost of secant stiffness of the material, and it ranges from ${\textstyle 0}$ for undamaged material to ${\textstyle 1}$ for fully degradated one ${\textstyle 0\leq d\leq 1}$.

In the present work, the effective stresses ${\textstyle {\overline {\boldsymbol {\sigma }}}}$ are computed in terms of the total strain tensor ${\textstyle {\boldsymbol {\varepsilon }}}$ as

 ${\displaystyle {\overline {\boldsymbol {\sigma }}}={\boldsymbol {C}}:{\boldsymbol {\varepsilon }}}$
(45)

where ${\textstyle {\boldsymbol {C}}}$ is the secant fourth order constitutive tensor and ${\textstyle (:)}$ denotes double contraction. The constitutive equations for the scalar isotropic damage model used in this work is:

 ${\displaystyle {\boldsymbol {\sigma }}=(1-d){\overline {\boldsymbol {\sigma }}}=(1-d){\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}}$
(46)

where ${\textstyle {\boldsymbol {C}}_{0}}$ 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 (45) is fully determined if the value of ${\textstyle d}$ can be evaluated at every time of the deformation process. To that extends one must define a equivalent effective stress or a suitable norm, ${\textstyle \tau }$. 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:

 ${\displaystyle \tau =\left(\theta +{\frac {1-\theta }{n}}{\sqrt {{\boldsymbol {\sigma }}_{0}:{\boldsymbol {C}}_{0}^{-1}:{\boldsymbol {\sigma }}_{0}}}\right)}$
(47)

where ${\textstyle n}$ is defined as ratio compressive strength/tensile strength and ${\textstyle \theta }$ is a weighting factor depending on the state of stress ${\textstyle \sigma _{0}}$. A satisfactory definition of ${\textstyle \theta }$ is

 ${\displaystyle \theta ={\frac {\sum _{i=1}^{3}\langle \sigma _{0}^{i}\rangle }{\sum _{i=1}^{3}|\sigma _{0}^{i}|}}}$
(48)

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

 ${\displaystyle {\mathcal {F}}(\tau ^{t},r^{t})=\tau ^{t}-r^{t}\leq {0}{\hbox{ }}\forall t\geq 0}$
(49)

where the ${\textstyle \tau ^{t}}$ is the norm described before and ${\textstyle r^{t}}$ is the current damage threshold at the time ${\textstyle t}$. This last value is an internal stress-like variable whose controls the size of damage surface. If ${\textstyle r_{0}=\sigma _{0}}$ is the initial threshold value, in this case, the initial uniaxial damage stress(tensile strength), it must be ${\textstyle r^{t}\geq r_{0}=\sigma _{0}}$. Damage occurs when the norm ${\textstyle \tau }$ exceeds the current threshold value. Fig 4 shows the elastic domain of the function defined in Eqn 49.

 Figure 4: Elastic domain selected for the isotropic damage model.

 ${\displaystyle {\dot {r}}\geq 0{\mathcal {F}}(\tau ,r)\leq {0}{\dot {r}}{\mathcal {F}}(\tau ,r)=0}$ ${\displaystyle {\hbox{if }}{\mathcal {F}}(\tau ,r)=0{\hbox{ then }}{\dot {r}}{\dot {\mathcal {F}}}(\tau ,r)=0}$
(50)

 ${\displaystyle {\dot {\tau }}={\dot {r}}}$
(51)

This, in turn, the evolution of the internal variable ${\textstyle r}$ may be explicitly integrated to render:

 ${\displaystyle r={\hbox{max}}\{r_{0},{\hbox{max}}(\tau (t))\}0\leq t<\infty }$
(52)

Note that (51) allows to compute the current values for ${\textstyle r}$ in terms of the current value of ${\textstyle \tau }$, which depends explicitly on the current total strains. Stress softening is controlled by the softening function ${\textstyle q=q(r)}$. In this work, the following exponential softening law is used:

 ${\displaystyle q(r)=r_{0}e^{(-2H_{S}({\frac {r-r_{0}}{r_{0}}}))}{\hbox{ , }}r\geq r_{0}}$
(53)

where ${\textstyle H_{S}}$ 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 ${\textstyle d=d(r)}$ is explicitly defined in terms of the corresponding current value of damage threshold,

 ${\displaystyle d(r)=1-{\frac {q(r)}{r}}{\hbox{ , }}r\geq r_{0}}$
(54)

so that it is monotonically increasing function such that ${\textstyle 0\leq d\leq 1}$.

7.3 Mechanical Dissipation

The mechanical free energy ${\textstyle \Psi }$ for an isotropic damage model is defined as follows:

 ${\displaystyle \Psi =(1-d)\Psi _{0}({\boldsymbol {\varepsilon }})=(1-d)[{\frac {1}{2}}{\boldsymbol {\varepsilon }}:{\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}]\geq 0}$
(55)

Thus, the rate of energy dissipation may be written as

 ${\displaystyle {\dot {\mathcal {D}}}=\Psi _{0}{\dot {d}}=\geq 0}$
(56)

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

 ${\displaystyle {\mathcal {D}}=\int _{0}^{\infty }{\dot {\mathcal {D}}}dt=\int _{0}^{\infty }\Psi _{0}{\dot {d}}dt}$ ${\displaystyle {\mathcal {D}}={\frac {1}{2E}}(\int _{r_{0}}^{\infty }q(r)dr-\int _{q_{0}}^{0}rdq)}$
(57)

Combining equation (57) with (54) results:

 ${\displaystyle {\mathcal {D}}=(1+{\frac {1}{H_{S}}}){\frac {\sigma _{0}^{2}}{2E}}=(1+{\frac {1}{H_{S}}}){\mathcal {U}}_{0}}$
(58)

Finally, as it is usual in smeared crack models, it is possible to relate the specific dissipated energy ${\textstyle {\mathcal {D}}}$

 ${\displaystyle {\mathcal {D}}{\mathcal {L}}_{e}={\mathcal {G}}_{f}}$
(59)

where ${\textstyle {\mathcal {G}}_{f}}$ is the mode I fracture energy per unit area (assumed to be a material property) and ${\textstyle {\mathcal {L}}_{e}}$ 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 ${\textstyle H_{S}}$ parameter can be expressed

 ${\displaystyle H_{S}={\frac {{\overline {H}}_{S}{\mathcal {L}}_{e}}{1-{\overline {H}}_{S}{\mathcal {L}}_{e}}}\geq 0}$
(60)

where the specific softening parameter ${\textstyle {\overline {H}}_{S}={\mathcal {U}}_{0}{\mathcal {G}}_{f}}$ measures the brittleness of the material and ${\textstyle H_{S}}$, 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, ${\textstyle {\mathcal {L}}_{e}=l_{e}}$. In this work, and assuming that the elements are almost equilateral, the size of the element can be computed ${\textstyle l_{e}^{2}={\sqrt {A_{e}}}}$ for quadrilateral elements being ${\textstyle A_{e}}$ the area of the finite element. It is interesting to comment that for a linear irreducible formulation, the discrete localization band ${\textstyle {\mathcal {w}}}$ is only one element across. The box 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 ${\textstyle \tau ^{t}}$ as the nonlinear behavior proceeds:

boxed BOXBox

Smeared Isotropic damage model.

1. Initial Data for time ${\textstyle t+1}$
1. Material Properties ${\textstyle \sigma _{t}}$,${\textstyle \sigma _{c}}$, ${\textstyle n}$, ${\textstyle E_{0}}$, ${\textstyle \nu }$, ${\textstyle {\mathcal {G}}_{f}}$
2. Current values: ${\textstyle {\boldsymbol {\varepsilon }}^{t+1}}$,${\textstyle d^{t}}$ and ${\textstyle r^{t}}$
2. Determine ${\textstyle H_{S}}$. Eqn 60.
3. If ${\textstyle t=0}$ the initialize ${\textstyle \tau _{0}=\sigma _{0}}$.
4. Evaluate undamaged stresses.
 ${\displaystyle {\boldsymbol {\sigma }}_{0}^{t+1}={\boldsymbol {C}}_{0}:{\boldsymbol {\varepsilon }}^{t+1}}$
5. Evaluate ${\textstyle \tau ^{n+1}}$. Eqn 47.
 ${\displaystyle \tau =\left(\theta +{\frac {1-\theta }{n}}{\sqrt {{\boldsymbol {\sigma }}_{0}:{\boldsymbol {C}}_{0}^{-1}:{\boldsymbol {\sigma }}_{0}}}\right)}$
6. Update internal variables
 ${\displaystyle r^{n+1}=max(r^{t},\tau ^{n+1})}$ ${\displaystyle d^{n+1}=d(r^{t+1})}$
7. Update stresses
 ${\displaystyle {\boldsymbol {\sigma ^{t+1}}}=(1-d^{n+1}){\boldsymbol {\sigma _{0}^{t+1}}}}$

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 ${\textstyle {\boldsymbol {C}}_{sec}=(1-d){\boldsymbol {C}}_{0}}$. 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 [6].

8.1 Two 2D nonlinear truss element

The first numerical example consist of two-member nonlinear truss element shown in the figure 5 and has been widely used as a benchmark problem for comparison of numerical solution algorithm. The following material are assumed: Axial rigidity ${\textstyle EA=1884.694lb}$, ${\textstyle h=25.874in}$ and ${\textstyle \alpha =63.4}$ 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 ${\textstyle P_{u}=0.05P_{v}}$ and the third one (Case III), the vertical load is treated as an imperfection, i.e ${\textstyle P_{v}=0.05P_{u}}$. The results obtained by the present method for the three considered cases are shown in Figs 6, fig:Sec4-3 and 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.

 Figure 5: Two non-linear truss element. Geometry and boundary conditions.
 Figure 6: Analysis result for case I.
 (a) Vertical Load ${\displaystyle P_{v}}$ vs horizontal displacements. (b) Vertical Load ${\displaystyle P_{v}}$ vs vertical displacements. Figure 7: Analysis results for Case II.
 (a) Vertical Load ${\displaystyle P_{u}}$ vs horizontal displacements. (b) Vertical Load ${\displaystyle P_{u}}$ vs vertical displacements. Figure 8: Analysis results for Case III.

8.2 A simple ${\displaystyle 2D}$ benchmarking example8.2 A simple 2 D {\displaystyle 2D} 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 ${\textstyle 2D}$ ${\textstyle 1\times {1}m\times m}$ block discretized in ${\textstyle 1}$, ${\textstyle 4}$ and ${\textstyle 16}$ quadrilateral bilinear element and drown in figure 9, 10a and 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 1. The softening modulus ${\textstyle H}$ for linear softening was taken ${\textstyle 1MPa}$. For exponential softening, we take the following formula for evolution of shear stress

 ${\displaystyle \tau _{s}({\kappa })=\tau _{s0}e^{-h\kappa }}$
(61)

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

 Material Data Value(${\textstyle N-m}$) Young's Modulus ${\textstyle E}$ 10 ${\textstyle MPa}$ Poisson's ratio ${\textstyle \upsilon }$ 0.33 Shear strength ${\textstyle \tau _{s0}}$ 1 ${\textstyle KPa}$
 Figure 9: Mesh ${\displaystyle one}$.
 (a) Mesh ${\displaystyle two}$. (b) Mesh ${\displaystyle three}$. Figure 10: Analysis meshes cases.

The results of the analysis is plotted in figures 11 and 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.

 (a) Vertical Load vs vertical displacement using DP-model with exponential softening. (b) Vertical Load vs vertical displacement using DP-model with linear softening. Figure 11: Analysis cases for DP-Model.
 (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. Figure 12: Analysis cases for Von Mises -Model.

8.3 A ${\displaystyle 2D}$ pull test benchmarking example8.3 A 2 D {\displaystyle 2D} 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 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

1. For Drucker-Prager model with exponential softening taken ${\textstyle \alpha =0.15}$ and ${\textstyle h_{1}=100}$ ${\textstyle h_{2}=300}$ and ${\textstyle h_{3}=180}$. Case I.
2. For Von-Mises model with perfect plastic. Case II.
3. For Von-mises model with with exponential softening ${\textstyle h=1000}$. Case III.

The Case I is represented is the figure 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 ${\textstyle h}$. Case II represent a perfect plastic behavior and the analysis results is plotted in fig 15. Note how the AL method implemented obtained the expected perfect plastic behavior. The Case III analysis is depicted in figure 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.

 (a) Mesh one. 50 nodes and 34 elements. (b) Mesh two. 605 nodes and 544 elements. (c) Mesh three. 167 nodes and 136 elements. Figure 13: Finite element meshes for analysis.
 Figure 14: Responses of pull test using three different parameter in exponential softening.
 Figure 15: Perfect plastic response for three meshes using Von-Mises model.
 Figure 16: Non objective response in Von-Mises material with exponential softening.

8.4 A ${\displaystyle 2D}$ notched beam8.4 A 2 D {\displaystyle 2D} 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 ${\textstyle {\mathcal {G}}_{f}}$, has been carefully specified. Finite element meshes, dimension and properties are shown in figure 17. ${\textstyle 702}$ nodes and ${\textstyle 612}$ 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 18 compare experimental results and numerical load deflection curves. Note how the softening behavior can be satisfactorily simulated together with the maximum load.

 (a) Geometry and material properties. (b) Mesh for analysis. Figure 17: Finite element meshes for analysis and basic properties.
 (a) Experimental and numerical results reported by Rots et al (1985). (b) Numerical result using smeared crack model with objective dissipation. 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 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 ${\textstyle 608}$ nodes and ${\textstyle 548}$ 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 2. Note that for define the smeared crack model is only required ${\textstyle 3}$ 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.

 (a) Geometry properties. (b) Mesh for analysis. Figure 19: Finite element meshes for analysis and basic properties.

 Material Data Value(${\textstyle N-m}$) Young's Modulus ${\textstyle E}$ 40 ${\textstyle GPa}$ Poisson's ratio ${\textstyle \upsilon }$ 0.2 Tensile strength ${\textstyle \sigma _{0}}$ 4 ${\textstyle MPa}$ Fracture Energy ${\textstyle {\mathcal {G}}_{f1}}$ 350 ${\textstyle N/m}$ Fracture Energy ${\textstyle {\mathcal {G}}_{f2}}$ 250 ${\textstyle N/m}$

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 21b . Note that when the damage reaches the value of ${\textstyle 1}$, means that this zone is totally smeared cracked.

 (a) Experimental and numerical results reported by Rots et al (1985. (b) Numerical result using smeared crack model with objective dissipation. Figure 20: Experimental and numerical result reported for DCB and numerical result using smeared crack model with objective dissipation.
 (a) Contour Sxx Cauchy Stress field at the end of simulation. (Gid pre-post processor[6]). (b) Evolution of damage for ${\displaystyle {\mathcal {G}}_{f1}}$ at the end of simulation. (Gid pre-post processor[6]). (c) Direccion of principal stress at the end of simulation. (Gid pre-post processor[6]). (d) Deformed mesh ${\displaystyle (\times {500)}}$. (Gid pre-post processor[6]). 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 ${\textstyle (0.5-0.75)}$ 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[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

[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

[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

[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–2 133–149

[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

[5] de Souza Neto, E. A. and Peri, D. and Owen, D. R. J. (2008) "Computational Methods for Plasticity". John Wiley & Sons, Ltd 1–15

[6] GiD. (2009) "The Personal Pre and Post Processor"

[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

[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–2 216–233

Document information

Published on 05/03/18
Submitted on 05/03/18