Latest revision as of 13:55, 20 March 2018

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:

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

or increment form as

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

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,

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

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

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]

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

where

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

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

and

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

where, in this case

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.

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

Which lead to

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

where

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:

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

(¹) 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}</li>}$

   (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}</li>}$ (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}</li>}$ (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.

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:

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 }$

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]:

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:

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 }$

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

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}})))</li>}$

   (42)

2. Incremental work. Follow the sign of the predictor work increment:

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

   (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}})</li>}$

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

${\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}}$

2. 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.
3. 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}</li>}$

4. 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}</li>}$

   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</li>}$

      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.
5. Apply correction to incremental load factor

   ${\displaystyle \Delta \lambda _{i+1}^{n+1}=\Delta \lambda _{i}-g_{i}+\delta \lambda _{i}</li>}$

6. Compute iterative displacement

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

7. 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}</li>}$

8. 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}</li>}$

9. 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).
10. 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}</li>}$

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

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:

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:

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

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:

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.

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

leading, in view of equation (49) , to the loading condition

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

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:

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,

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:

Thus, the rate of energy dissipation may be written as

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

Combining equation (57) with (54) results:

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

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

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}</li>}$

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)</li>}$

6. Update internal variables

   ${\displaystyle r^{n+1}=max(r^{t},\tau ^{n+1})}$ ${\displaystyle d^{n+1}=d(r^{t+1})</li>}$

7. Update stresses

   ${\displaystyle {\boldsymbol {\sigma ^{t+1}}}=(1-d^{n+1}){\boldsymbol {\sigma _{0}^{t+1}}}</li>}$