Combination of the critical displacement method with a damage model for structural instability analysis


==Abstract==

The paper describes the extension of the critical displacement method (CDM), presented by Oñate and Matias in 1996, to the instability analysis of structures with non-linear material behaviour using a simple damage model. The extended CDM is useful to detect instability points using a prediction of the critical displacement field and a secant load-displacement relationship accounting for material non-linearities. Examples of application of CDM to the instability analysis of structures using bar and solid finite elements are presented.

'''Keywords:''' Finite element method, Structural defects

==1 Introduction==

Mathematical phenomena such as bifurcation and limit points, commonly known as critical points, have received considerable  attention in the computational structural community <span id='citeF-3'></span>[[#cite-3|[3]]], <span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-6'></span>[[#cite-6|[6]]], <span id='citeF-17'></span>[[#cite-17|[17]]]&#8211;<span id='citeF-25'></span>[[#cite-25|[25]]]. In a practical context the occurrence of those points in equilibrium paths of load&#8211;deflection diagrams has serious consequences, as it means a loss of stability and a form change of the whole structure. Typical examples are the buckling of rods, plates and shell structures, diffuse necking bifurcation problems and the formation of shear bands in elastic-plastic solids. A detection and localization of critical points is therefore crucial as it permits a safe dimensioning of structures according to the expected loads.

Numerical methods for computation of critical points can be grouped into two categories, namely indirect and direct methods. With indirect methods the encounter of the critical point is judged with the help of a detection parameter, while the equilibrium path is being traced <span id='citeF-3'></span>[[#cite-3|[3]]]. In direct methods the condition for occurrence of a critical point is included in the system to be solved <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-25'></span>[[#cite-25|[25]]] . Arc&#8211;length methods, essentially a particular class of continuation methods, provide a useful tool for the tracking  of the equilibrium path <span id='citeF-4'></span>[[#cite-4|[4]]].

A blending of direct and indirect methods to compute critical points named as the ''critical displacement method ''(CDM) was developed recently for linear elastic bar and solid finite elements <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-14'></span>[[#cite-14|[14]]], <span id='citeF-15'></span>[[#cite-15|[15]]]. The basic idea of the CDM is to predict first the critical displacement pattern. This is done by using the tangent stiffness singularity condition at the critical point using a predicted perturbation of the last converged displacement field. The resulting eigenvalue problem is then linearized and solved. Via the secant load&#8211;displacement stiffness relationship the critical load is finally obtained <span id='citeF-12'></span>[[#cite-12|[12]]]&#8211;<span id='citeF-16'></span>[[#cite-16|[16]]].

In this paper the application of the CDM to the stability analysis of structures using a simple isotropic damage model will be presented. The damage model allows to estimate structural failure due to geometrical and material non linearities  under different stress and strain conditions. After a short introduction to the CDM and the damage model chosen the efficiency of the CDM is evaluated in the stability analysis of structures using truss and solid elements.

==2 Theoretical background of the CDM==

In this section a short introduction on the Critical Displacement Method (CDM) will be given. The approach follows closely the ideas presented in <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-14'></span>[[#cite-14|[14]]], <span id='citeF-15'></span>[[#cite-15|[15]]].

===2.1 Derivation of the secant stiffness matrix===

The secant stiffness equation relating total (or incremental) displacements and forces is derived via the principle of virtual work (PVW). In order to express the full incremental form of the PVW,  a body in its initial configuration with volume <math display="inline">{}^0V</math> is considered as a starting point (Figure [[#img-1|1]]). The body is submitted to different load levels,  here denoted by the superscript <math display="inline">t</math>, in a quasi-static analysis. Starting from a known actual configuration <math display="inline">{}^tV</math> the solution for an updated configuration at load level <math display="inline">t + \Delta t</math> is sought. Furthermore, the body is described in a generalized Lagrange formulation with respect to an arbitrary reference configuration at load level <math display="inline">{}^rV</math>. This can be made coincident with the initial configuration <math display="inline">{}^r V =  {}^0 V</math> (total Lagrange approach) or the current configuration <math display="inline">{}^r V  =  {}^t V</math> (updated Lagrange approach)<span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-6'></span>[[#cite-6|[6]]], <span id='citeF-13'></span>[[#cite-13|[13]]]. The configuration at the time <math display="inline">t + \Delta t</math> corresponds to a variation of the known configuration at the time <math display="inline">t</math>.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-fig1.png|500px|Deformation process of a solid]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Deformation process of a solid
|}

The displacements at the time <math display="inline">t + \Delta t</math> are written as the sum of the known displacements <math display="inline">{\boldsymbol {}^tu}</math> and the sought displacement increments <math display="inline">\Delta {\boldsymbol u}</math>, i.e.

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

Substituting the displacements in the Green strain tensor given by equation ([[#eq-2|2]]) and sorting the terms according to their dependence on the incremental displacements, leads to equation ([[#eq-3|3]]), where the strain increment vector is split into first and second order parts

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}_r^{t + \Delta t} \epsilon _{ij} \; = \; \frac{1}{2} \left({}^{t +\Delta t}    u_{i,j}+{}^{t + \Delta t} u_{j,i} + {}^{t + \Delta t}u_{k,i}    {}^{t + \Delta t} u_{k,j}\right)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}_r^{t + \Delta t} \Delta \epsilon _{ij} \; = \; {}_r^{t + \Delta t}     \epsilon _{ij} - {}_r^t \epsilon _{ij} = {}_r e_{ij} + {}_r \eta _{ij}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

In above <math display="inline">{}_r e_{ij}</math> and <math display="inline">{}_r \eta _{ij}</math> contain terms with linear and quadratic dependence on the incremental displacements <math display="inline">\Delta {\boldsymbol u}</math>, respectively, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}_r e_{ij} \; = \; \frac{1}{2} \left({}_r \Delta u_{i,j} +     {}_r \Delta u_{j,i} + {}_r^t u_{k,i} \; {}_r \Delta u_{k,j} +     {}_r \Delta u_{k,i} \; {}_r^t u_{k,j} \right)</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|-
| style="text-align: center;" | <math>     {}_r \eta _{ij} \; = \; \frac{1}{2} {}_r \Delta u_{k,i} \; {}_r \Delta u_{k,j} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}
|}

A comma in a lower right index denotes a derivative with respect to the Cartesian coordinate referred to by the index following (<math display="inline">u_{i,j}=\frac{\partial u_i}{\partial x_j}</math> with <math display="inline">i,j,k = 1,2,3</math>). A similar expression can be deduced for the virtual displacements. Taking into account the fact that <math display="inline">\delta {}^t u_i = 0</math> leads to the following expressions:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta {}^{t+\Delta t} u_i \; =\; \delta \Delta u_i  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|-
| style="text-align: center;" | <math>     \delta {}_r^{t + \Delta t} \Delta      \epsilon _{ij}  \; =  \; \delta {}_r e_{ij} + \delta {}_r \eta _{ij} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}
|}

The first and second order virtual strain increments are

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta {}_r e_{ij} \; = \; \frac{1}{2} \left(\delta {}_r \Delta u_{i,j} +     \delta {}_r \Delta u_{j,i} + {}_r^t u_{k,i} \; \delta {}_r \Delta u_{k,j} +     \delta {}_r \Delta u_{k,i} \; {}_r^t u_{k,j} \right)</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math>     \delta {}_r \eta _{ij} \; = \; \frac{1}{2} \left(\delta {}_r \Delta      u_{k,i} \; {}_r \Delta u_{k,j} + {}_r \Delta u_{k,i} \; \delta {}_r \Delta      u_{k,j} \right) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}
|}

The Piola-Kirchoff stress vector at <math display="inline">t + \Delta t</math> is obtained incrementally by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}_r^{t + \Delta t}\sigma  \;=\;    {}_r^t \sigma  + {}_r \Delta \sigma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

The constitutive equation chosen assumes a linear dependence between the second Piola-Kirchhoff stress increments and the Green-Lagrange strain increments as

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}_r \Delta \sigma  \; = \; {}_r {\boldsymbol D} ({}_r {\boldsymbol e} +    {}_r \eta )  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">{}_r {\boldsymbol D}</math> is a fourth order constitutive tensor.

The expression of the PVW at time <math display="inline">t + \Delta t</math> can be written as

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{{}^{r}V} \delta {}^{t + \Delta t}_r    \epsilon _{ij}\; {}^{t + \Delta t}_r \sigma _{ij} dV =    \int _{{}^{r}V} \delta u_i \;{}^{t+    \Delta t}f_i^B dV +\int _{{}^{r}S} \delta u_i \; {}^{t+    \Delta t}f_i^S dS  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math display="inline">{}^{t + \Delta t}_r \sigma _{ij}</math> are the second Piola&#8211;Kirchhoff stress and <math display="inline">f_i^B</math>, <math display="inline">f_i^S</math> denote the volume and surface forces, respectively. As usual <math display="inline">{}^rV</math> and <math display="inline">{}^rS</math> denote the volume and the boundary of the reference configuration.

Inserting equations ([[#eq-3|3]]) and ([[#eq-11|11]]) into ([[#eq-12|12]]) gives the full incremental form of the PVW as

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{{}^rV} [ \delta {}_r e_{ij} \; {}_rD_{ijkl} \;{}_r e_{kl} +     \delta {}_r e_{ij}\; {}_r D_{ijkl} \; {}_r \eta _{kl} + \delta {}_r     \eta _{ij} \; {}_rD_{ijkl} \; {}_r e_{kl} + \delta {}_r \eta _{ij} \;{}_r     D_{ijkl} \; {}_r \eta _{kl}+  </math>
|-
| style="text-align: center;" | <math>     + \delta {}_r \eta _{ij} \; {}_r^t S_{ij}     ] dV \; = \; {}_r^{t + \Delta t}  f -     \int _{{}^rV} \delta {}_r e_{ij} \; {}_r^t \sigma _{ij} dV  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t + \Delta t}_r f \; = \; \int _{{}^rV} \delta {}^{t + \Delta t}     u_i {}_0^{t + \Delta t} f^B_i dV +     \int _{{}^r S} \delta {}^{t + \Delta t} u_i {}^{t +        \Delta t} f^S_i dS </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Eq.([[#eq-13|13]]) is the point of departure for the derivation of the secant stiffness matrix <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-11'></span>[[#cite-11|[11]]]&#8211;<span id='citeF-13'></span>[[#cite-13|[13]]].

As a next step the displacements and their increments are discretised in the standard finite element manner <span id='citeF-26'></span>[[#cite-26|[26]]]

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t {\boldsymbol u} \; = \; {\boldsymbol N} {}^t {\boldsymbol a} \quad , \quad          \Delta {\boldsymbol u} \;  =  \; {\boldsymbol N} \Delta {\boldsymbol a} \quad , \quad          \delta (\Delta {\boldsymbol u}) \; = \; {\boldsymbol N} \delta (\Delta {\boldsymbol a})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">{\boldsymbol N}</math> and <math display="inline">{\boldsymbol a}</math> are the shape function matrix and the nodal displacement vector, respectively.

''The incremental secant stiffness'' expression is now obtained by substituting equations ([[#eq-4|4]]), ([[#eq-5|5]]), ([[#eq-8|8]]), ([[#eq-9|9]]), ([[#eq-10|10]]), ([[#eq-11|11]]) and ([[#eq-15|15]]) into ([[#eq-13|13]]) and collecting the terms depending on <math display="inline">{\boldsymbol \Delta a}</math> in the following form:

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t_r{\boldsymbol K}_S({}^t {\boldsymbol a},\Delta {\boldsymbol a}) \; {\boldsymbol \Delta a} \;=\;       {}^{t + \Delta t}_r {\boldsymbol r}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

where the incremental secant stiffness matrix is given by

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t_r{\boldsymbol K}_S(\Delta {\boldsymbol a}) \; =     {}^t_r{\boldsymbol K}_L + {}^t_r{\boldsymbol K}_M(\Delta {\boldsymbol a}) +     {}^t_r{\boldsymbol K}_N(\Delta {\boldsymbol a}^2) +     {}^t_r {\boldsymbol K}_\sigma (\boldsymbol \sigma )  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

Note that the secant matrix <math display="inline">{}_r^t{\boldsymbol K_S}</math> is split in terms according to their dependence on the incremental nodal displacements <math display="inline">\Delta {\boldsymbol a}</math>. The form of the different matrices in equation ([[#eq-17|17]]) for truss and solid elements con be found in <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-13'></span>[[#cite-13|[13]]] and <span id='citeF-15'></span>[[#cite-15|[15]]] .

In eq.(16) <math display="inline">{}^{t + \Delta t} {\boldsymbol r}</math> is the residual force vector given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t + \Delta t}_r {\boldsymbol r} \; = \; {}^{t + \Delta t} {\boldsymbol f}      - \int _{{}^{r}V} {}^t_r{\boldsymbol B}^T_L      {}^t_r \sigma  _{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

In eq.([[#eq-16|16]]) <math display="inline">{}^{t + \Delta t}{\boldsymbol f}</math> is the equivalent nodal force vector deduced from eq.([[#eq-14|14]]) and <math display="inline">{}^t_r{\boldsymbol B}_L</math> is the first order strain matrix relating first order strain increments <math display="inline">{}_r e_{ij}</math> and nodal displacements <math display="inline">\Delta u_i</math>. The expression of this matrix can be found in <span id='citeF-13'></span>[[#cite-13|[13]]] and <span id='citeF-15'></span>[[#cite-15|[15]]].

Equation ([[#eq-16|16]]) is a non linear system of algebraic equations relating the increment of forces with the displacement increments. Indeed this secant equation can be used as a basis for deriving different non linear solution schemes to track the load&#8211;displacement pattern <span id='citeF-13'></span>[[#cite-13|[13]]], <span id='citeF-14'></span>[[#cite-14|[14]]]. An alternative procedure for deriving secant stiffness expressions can be found in <span id='citeF-16'></span>[[#cite-16|[16]]].

Note that a “total” secant stiffness matrix can be obtained by making <math display="inline">\Delta {\boldsymbol a}={\boldsymbol a}^t</math>. In a total Lagrangean description (i.e. <math display="inline">{}^r V = {}^0 V</math>) the total secant matrix is given by <span id='citeF-13'></span>[[#cite-13|[13]]]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t_0{\boldsymbol K}_S={}^t_0 {\boldsymbol K}_L + {}^t_0{\boldsymbol K}_M({}^t{\boldsymbol a}) +     {}^t_0{\boldsymbol K}_N({}^t{\boldsymbol a}^2) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Figure [[#img-2|2]] shows schematically the difference between the incremental and total secant approaches. Further details can be found in <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-13'></span>[[#cite-13|[13]]]&#8211;<span id='citeF-15'></span>[[#cite-15|[15]]].

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-fig2.png|500px|Total (left) and incremental (right) secant approaches]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Total (left) and incremental (right) secant approaches
|}

It is important to note that in the limiting case of <math display="inline">{\boldsymbol \Delta a} \rightarrow 0</math> the secant stiffness matrix reduces to the tangent stiffness matrix, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t_r {\boldsymbol K}_T = \lim _{\Delta {\boldsymbol a} \rightarrow 0} {}^t_r     {\boldsymbol K}_S = {}^t_r {\boldsymbol K}_L + {}^t_r {\boldsymbol K}_\sigma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

Using <math display="inline">{}^t_r {\boldsymbol K}_T</math> the standard iteration techniques for the computation of stability points in structures can be obtained. For instance, the standard Newton&#8211;Raphson method is written from eq.([[#eq-16|16]]) as

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^t_r {\boldsymbol K}_T({}^t{\boldsymbol a}) \Delta {\boldsymbol a} \;     =  {}^{t+ \Delta t}_r{\boldsymbol r}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

In our work this theory was implemented for two-node truss elements and 8 node solid elements. The detailed form of all the relevant matrices can be seen in <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-13'></span>[[#cite-13|[13]]] and <span id='citeF-15'></span>[[#cite-15|[15]]].

===2.2 Computing the equilibrium path===

In order to trace the equilibrium path the load terms have to be scalable. For this purpose, the scalar load parameter <math display="inline">{}^t \lambda </math> is introduced. This parameter leaves the pattern of the external loads unchanged, only the size is varied (i.e. <math display="inline">{}^t{\boldsymbol f}={}^t\lambda {\boldsymbol f}_0</math>).

The arc&#8211;length method is commonly used to calculate equilibrium solutions for different loads and displacements. Therefore the equation defining the load increment is formally extended by one equation <math display="inline">f({}^t\lambda ,{\boldsymbol {}^tu})=0</math>, which determines the load parameter <math display="inline">{}^t \lambda </math>. A variety of different control equations can be found in the literature <span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-6'></span>[[#cite-6|[6]]]. The ones used here are stated below in equation ([[#eq-22|22]]).

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f({}^t\lambda ,{\boldsymbol {}^tu}) \; = \begin{cases}   \qquad {}^{t+\Delta t} u_i - c & \mbox{ displacement control}\\   ({\boldsymbol u_m} - {}^t {\boldsymbol u})^T ({}^{t+\Delta t}{\boldsymbol u}-{\boldsymbol u_m})+(   \lambda _m - {}^t \lambda )({}^{t+\Delta t}\lambda - \lambda _m) & \mbox{ normal plane} \\ \\   \sqrt{\|{}^{t+\Delta t}{\boldsymbol u} - {}^t {\boldsymbol u}\| ^2 +({}^{t+\Delta t}\lambda - {}^t \lambda )^2}   - ds & \mbox{ cylindrical}\end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

For the control equation the normal plane iteration was found to be rather  stable and reliable in this work. Note however that the type of control equation can be problem dependent.

===2.3 Prediction of instability points===

The main criterion for existance of a critical point is that the tangent stiffness matrix <math display="inline">{\boldsymbol K}_T</math> becomes singular at the point. Mathematically this means that equation ([[#eq-23|23]]) below has nontrivial solutions.

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}_T({\boldsymbol a}_c) \Delta{\boldsymbol a} \; = \; 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

where <math display="inline">{\boldsymbol a}_c</math> is the critical displacement vector.

This situation corresponds to either bifurcation points, where the equilibrium path splits into a second branch, or limit points, where a snap through occurs.

The singularity of the tangent matrix <math display="inline">{\boldsymbol K}_T</math> has some further consequences. The determinant of <math display="inline">{\boldsymbol K}_T</math> is equal to <math display="inline">0</math> in the critical point, furthermore the smallest eigenvalue is <math display="inline">0</math>. This can be stated mathematically as

<span id="eq-24.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\det {\boldsymbol K}_T({\boldsymbol a}_c) \; =  \; 0   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24a)
|}

<span id="eq-24.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[{\boldsymbol K}_T({\boldsymbol a}_c) - \omega {\boldsymbol I}\right]\Phi \;  = \; 0    \qquad \mbox{with} \quad \omega _1 = 0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24b)
|}

Many indirect methods for the computation of critical points use condition ([[#eq-24.a|24a]]) to determine their location on the equilibrium path. Nevertheless, in order to verify if a point indeed is a critical point these criteria can serve as test functions. Plotting <math display="inline">det {\boldsymbol K}_T</math>, the curve should become <math display="inline">0</math> at the critical point according to eq.([[#eq-24.a|24a]]).

Furthermore eq.([[#eq-24.b|24b]]) indicates that the smallest eigenvalue turns zero at <math display="inline">{\boldsymbol a}_c</math>. Keeping track of the number of negative diagonal elements of <math display="inline">{\boldsymbol K}_T</math> during the computation of equilibrium points on the path permits as well to detect critical points.

As a rather common method for critical point computation, the standard limit load analysis shall be mentioned briefly. This it is carried out by solving the following eigenvalue problem

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[{}^{t}_r {\boldsymbol K}_{L} + {\lambda }\; {}^{t}_r    {\boldsymbol K}_\sigma \right]{\boldsymbol a} \;=\; 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

The estimated critical load is then obtained as a multiple of the known load at the equilibrium point.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol f}_{c}\;=\; {\lambda } \; {}^{t}{\boldsymbol f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

The approach termed as the Critical Displacement Method (CDM) <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-14'></span>[[#cite-14|[14]]]&#8211;<span id='citeF-16'></span>[[#cite-16|[16]]] is based on approximating the tangent matrix singularity condition at the critical point in terms of known values at the last converged solution. The unknown critical displacement field <math display="inline">{\boldsymbol a}_c</math> is expressed as

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{t + \Delta t_c}{\boldsymbol a}_c\; = \; {}^t{\boldsymbol a} + \Delta {\boldsymbol a}_c   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">{}^t{\boldsymbol a}</math> are the nodal displacements at the actual equilibrium configuration from which the prediction has to be made and <math display="inline">\Delta {\boldsymbol a}_c</math> is an estimate of the critical displacement increment pattern yielding structural instability at <math display="inline">t_c = t+\Delta t_c</math>. It is now assumed <math display="inline">\Delta {\boldsymbol a}_c</math> to be of the form <math display="inline">\Delta {\boldsymbol a}_c = \rho \Phi </math> with <math display="inline">\Phi </math> being an approximation for the critical displacement pattern and <math display="inline">\rho </math> a multiplier. As a convenient starting value of <math display="inline">\Phi </math> the last known value of the displacement pattern <math display="inline">{}^t{\boldsymbol a}</math> can be chosen (i.e. <math display="inline">\Phi ={}^t{\boldsymbol a}</math>).

Inserting this approximation in the expression for the tangent stiffness matrix (20) together with the condition ([[#eq-24.a|24a]]), leads to

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{det} \left[{\boldsymbol K}_T({}^t{\boldsymbol a}) + \rho {\boldsymbol K}_1(\Phi ) +      \rho ^2 {\boldsymbol K}_2(\Phi ^2)\right]= 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

where the tangent matrix at the critical point has been split into terms according to their dependence on <math display="inline">\Phi </math>. The detailed form of the matrices involved in eq.([[#eq-28|28]]) can be seen in <span id='citeF-8'></span>[[#cite-8|[8]]] and <span id='citeF-15'></span>[[#cite-15|[15]]]. Note that the resultant equation ([[#eq-28|28]]) is a quadratic eigenvalue problem, that can be solved for <math display="inline">\rho </math>.

Eq.([[#eq-28|28]]) can be simplified to the following standard linear eigenvalue problem by neglecting second order terms in <math display="inline">\rho </math>

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{det}   \left[{\boldsymbol K}_T({}^t{\boldsymbol a}) + \rho {\boldsymbol K}_1(\Phi )      \right]= 0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

For the solution of eq.([[#eq-29|29]]) the Inverse Iteration method has proved to be the best <span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-4'></span>[[#cite-4|[4]]]. The Rayleigh quotient <math display="inline">\rho </math> is used here as an approximation for the first eigenvalue. The convergence of this procedure can be enhanced by shifting the terms of matrix <math display="inline">{\boldsymbol K}_T</math>. The Jacobi iteration method as an alternative is very time consuming and thus not very suitable for large problems.

Having obtained the smallest eigenvalue <math display="inline">\rho </math>, the critical displacements are predicted from equation ([[#eq-27|27]]) with <math display="inline">{\boldsymbol \Delta a}_c = \rho {\boldsymbol \Phi }</math>. For the determination of the corresponding critical loads <math display="inline">{\boldsymbol f}_c</math> these displacements are inserted into the incremental secant stiffness matrix (see eq.([[#eq-16|16]])), which  gives the prediction for the critical load as

<span id="eq-30"></span>
<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta {\boldsymbol f}_c \;  =  \; {}^t_r{\boldsymbol K}_S({}^t{\boldsymbol a},\rho      \Phi ) \; \rho \Phi  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|-
| style="text-align: center;" | <math>     {}^{t_c}{\boldsymbol f} \;  =  \; {}^t{\boldsymbol f} + \Delta {\boldsymbol f}_c  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}
|}

The total secant stiffness matrix of eq.([[#eq-19|19]]) can be used to directly calculate <math display="inline">{}^{t_c}{\boldsymbol f}</math> as shown in Figure [[#img-2|2]]. Further details can be found in <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-13'></span>[[#cite-13|[13]]]&#8211;<span id='citeF-15'></span>[[#cite-15|[15]]].

The approach proposed above can be applied in different ways so as to obtain the different approximations for the critical load value. For the prediction of critical points the following strategy has been used:

<ol>

<li>A point on the equilibrium path is              computed using any standard iterative procedure. As above mentioned, arc-length method based on a normal plane iteration has been chosen in this work.         </li>
<li>Starting from the displacements <math display="inline">{}^t {\boldsymbol a}</math>,               the critical displacements <math display="inline">{}^{t_c}{\boldsymbol a}</math> are               found by solving equation ([[#eq-29|29]]) and using eq.([[#eq-27|27]]) with <math display="inline">\Delta {\boldsymbol a}_c=\rho ^t {\boldsymbol a}</math>.         </li>
<span id="point2"></span>
<li>The corresponding critical load is obtained from               equations ([[#eq-20|20]]) and ([[#eq-31|31]])         </li>
<li>A further point on the equilibrium path is computed </li>

<li>Go to [[#point2|2]]. </li>

</ol>

Compared to the standard limit load analysis where the critical load is obtained directly, the CDM predicts the critical displacements first. Then in an additional step the critical load is determined.  Note however that the cost of the CDM is similar to that of the standard limit load computations as both involve the solution of an eigenvalue problem.

When assembling the secant stiffness matrices special care has to be taken if  a damage model is included. As it will be shown next the critical displacements might belong to a damaged region in the load&#8211;deflection diagram, so that a different constitutive relation has to be applied.

==3 Combination of the CDM with a damage model==

In order to study the behavior of the CDM when a more complex constitutive model is used, a  damage model has been implemented. Such a model permits to reproduce material failure under different stress or strain conditions. For simplicity reasons the single parameter isotropic damage model is chosen <span id='citeF-10'></span>[[#cite-10|[10]]]. Indeed more sophisticated damage models can be used <span id='citeF-9'></span>[[#cite-9|[9]]].

The damage parameter <math display="inline">d</math> is a measure of the loss of secant stiffness  of the material. It ranges from <math display="inline">0</math> (undamaged material) to <math display="inline">1</math> (fully de-gradated material). The constitutive law is now written as

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

where <math display="inline">\overline {\boldsymbol D}</math> is the damage constitutive matrix.

Similar to plasticity analysis, a damage initiation criterion and damage evolution laws have to be defined.

A damage criterion is formulated in terms of a stress norm <math display="inline">\overline{\sigma }</math> and a damage threshold value <math display="inline">r(d)</math> as .

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>F(\epsilon ,d) \; = \; \overline{\sigma }   - r(d) \leq 0 </math>
|-
| style="text-align: center;" | 
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

In our work we have taken

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\sigma }(\epsilon ) \; = \;   \sqrt{\epsilon : {\boldsymbol D}: \epsilon } </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|-
| style="text-align: center;" | <math>   r(d) \; = \; \frac{r_0}{1\; - \;(1+H)d}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}
|}

In eq.([[#eq-35|35]]) <math display="inline">r_0</math> is the initial damage tdshreshold value and <math display="inline">H</math> the hardening  modulus. The stress norm <math display="inline">\overline{\sigma }</math> is the square root of the product of the undamaged stress tensor with the strain tensor. For the simplest one dimensional form of this damage model, the stress&#8211;strain curve is shown in Figure [[#img-3|3]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-pic08.png|300px|Constitutive plot for the single parameter damage model]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Constitutive plot for the single parameter damage model
|}

If inequality ([[#eq-33|33]]) is violated, damage occurs. The damage parameter <math display="inline">d</math> is then updated according to the following evolution law

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

which can be deduced from ([[#eq-35|35]]) and the damage criterion ([[#eq-33|33]]) when the equality sign is applied. In a damage case <math display="inline">d</math> is no longer constant, which affects the tangent stiffness matrix <math display="inline">{\boldsymbol K}_T</math>. Using equations ([[#eq-32|32]]) and ([[#eq-35|35]]) the tangent constitutive matrix <math display="inline">{\boldsymbol C}</math> is obtained as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol C} \; = \begin{cases} \; (1-d) {\boldsymbol D} - \frac{r_0}{(1+H)} \;   \frac{1}{\overline{\sigma } \;^3} ({\boldsymbol D}   \epsilon ) \otimes    ({\boldsymbol D}\epsilon ) & for F > 0 \\   (1-d) {\boldsymbol D} & for F \le 0\end{cases}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

The damage model is described in an algorithmic form.
{| class="floating_imageSCP" style="text-align: left; border: 1px solid #BBB; margin: 1em auto; width: 50%;max-width: 50%;font-size: 100%;"
|-
|style="text-align: center;"|Compute: <math>\overline{\sigma }^{m+1} \; = \; \sqrt{\epsilon ^{m+1}:{\boldsymbol D}:\epsilon ^{m+1} }</math>
|-
|style="text-align: center;"|if <math>\overline{\sigma } \le r^m</math>
|-
|-
| style="text-align: center;" | <math>\begin{array}{l}      r^{m+1} \; & = & \; r^m \\      d^{m+1} \; & = & \; d^m \\      \sigma ^{m+1} \; & = & \; (1-d^{m+1})      {\boldsymbol D} \epsilon ^{m+1} \\      {\boldsymbol C}^{m+1}  \; & = & \; (1-d^{m+1}) {\boldsymbol D} \\ \end{array}</math>
|-
|-
|sstyle="text-align: center;"|else <math>\overline{\sigma } >  r^m</math>
|-
|-
| style="text-align: center;" | <math>\begin{array}{l}      r^{m+1} \; & = & \; \overline{\sigma }^{m+1} \\      d^{m+1} \; & = & \; d(r^{m+1}) \\      \sigma ^{m+1} \; & = & \; (1-d^{m+1})      {\boldsymbol D} \epsilon ^{m+1} \\      {\boldsymbol C}^{m+1}  \;& = &\; (1-d^{m+1}) {\boldsymbol D} - \frac{r_0}{(1+H)} \;      \frac{1}{\overline{\sigma } \;^3} ({\boldsymbol D}      \epsilon ) \otimes       ({\boldsymbol D}\epsilon ) \end{array}</math>
|}

Note that initially the trial stresses <math display="inline">\overline{\sigma }^{m+1}</math> are computed and then compared with the actual threshold value <math display="inline">r^m</math>. Following this comparison the values <math display="inline">r^{m+1}</math> and <math display="inline">d^{m+1}</math> are updated (damage, ''if''&#8211;case) or not (undamaged, ''else''&#8211;case). Index <math display="inline">m</math> indicates here the different load levels.

==4 Prediction of instability points with the CDM and a damage model==

Having computed an equilibrium point with this algorithm a prediction of the critical load is made with the CDM as described in a previous section. When the approximate tangent stiffness matrix is  assembled the damaged constitutive tangent matrix <math display="inline">{\boldsymbol C}^{m+1}</math> in order to solve equation (29). Index <math display="inline">m</math> is used here to emphasize that the damage parameter of the current time load step corresponding to the last convergent solution in the equilibrium path is chosen. Solution of eq.(29) gives the critical displacemente pattern as previously explained.

As a next step the critical load is obtained from eqs.(30) and (31) For this purpose the incremental secant stiffness matrix has to be constructed. It is desirable to perform this computation taking into account that at the critical point elements might be in a further damaged state. The solution to this problem lies in the use of a critical damage parameter <math display="inline">d_c</math>. Calculating first the critical Green stresses <math>\epsilon _c</math>, the stress norm <math display="inline">\overline{\sigma }</math> is obtained as

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

The critical damage parameter is now chosen according to the following scheme:

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d_c = \begin{cases} d^m  & for \overline{\sigma }(\epsilon _c) \le r^m \\   \frac{1}{1+H} \left(1 - \frac{r_0}{\overline{\sigma }(\epsilon _c)} \right)& for   \overline{\sigma }(\epsilon _c) > r^m\end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

The option <math display="inline">d_c=d^m</math> means that the damage level of the last equilibrium position is applied. Alternatively a new damage parameter is computed using the stress norm from eq.([[#eq-38|38]]). The value of <math display="inline">d_c</math> is then used to estimate the secant stiffness matrix at the critical point. As the critical point prediction is independent from the computation of equilibrium points on the load&#8211;deflection path, by no means an update of the damage parameter or the threshold value should be made at this prediction step.

==5 Examples==

In this section examples of critical point predictions using the CDM and the damage model previously described will be presented. The first two examples use 2 node truss elements and the last two ones 8 node solid elements to illustrate the potential of the method.

===5.1 Simple arc structure with bar elements===

An outline of the first example is shown in Figure [[#img-4|4]]. The product of Young's modulus with the cross area of the bar is <math display="inline">E A = 5.0 \times 10^3</math>. A unit load <math display="inline">P = 1</math> is applied at the top node. All the displacements plotted in Figures [[#img-5|5]]-[[#img-8|8]] are the vertical displacements of node 4. The values are multiplied by <math display="inline">-1</math>

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-pic05.png|600px|Discretization of the arc structure into tow node bar element]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Discretization of the arc structure into tow node bar element
|}

Figure [[#img-5|5]] shows the equilibrium path giving a limit load <math display="inline">P_{lim} = 3.7</math> and a bifurcation load of <math display="inline">P_{bif} = 3.47</math>. This path has been obtained by a standard incremental load process using a normal plane iteration control procedure. The detailed results correspond to values of the critical displacement and the critical load predicted from different points of the equilibrium path using the CDM with a standard (undamaged) linear elastic constitutive model. Note that a value of <math display="inline">P_c = 3.0</math> (<math display="inline">\sim 15 % </math> error) is predicted from the undeformed configuration.

Figures [[#img-6|6]]&#8211;[[#img-8|8]] show that introducing the damage model of previous section diminishes the value of the critical load of the structure while the values of the critical displacement are not so different. For a better comparison of results, the undamaged equilibrium path is drawn as a dotted line.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-9b3149.png|600px|Arc structure. Equilibrium path and critical load prediction for nondamaged material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Arc structure. Equilibrium path and critical load prediction for nondamaged material
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-888c36.png|600px|Arc with damage: r₀= 0.1, H = 1]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Arc with damage: <math>r_0 = 0.1</math>, <math>H = 1</math>
|}

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-30306b.png|600px|Arc with damage: r₀= 0.1, H = 0.5]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Arc with damage: <math>r_0 = 0.1</math>, <math>H = 0.5</math>
|}

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-c8d086.png|600px|Arc with damage: r₀= 0.05, H = 1]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Arc with damage: <math>r_0 = 0.05</math>, <math>H = 1</math>
|}

The plots for different damage threshold values <math display="inline">r_0</math> and hardening modulus <math display="inline">H</math> demonstrate that the prediction follows well the new equilibrium path. The critical point predicted by the CDM is now the maximum load point. It is quite remarkable that the values of the critical load are predicted very accurately from the undeformed configurations in all cases. This is very useful for practical purposes as ideally a single computation from the original configuration should provide a good estimate of the instability load.

===5.2 Curved bridge===

The next example is a curved bridge formed by 35 bar elements connecting 19 nodes. (Figure [[#img-9|9]]) The arrow at node 10 symbolizes again a unit load. A value of <math display="inline">E A = 3.0 \times 10^6</math> is used. In all figures the vertical displacements of node 10 multiplied by <math display="inline">-1</math> versus the applied load are plotted.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-pic07.png|600px|Schematic picture of the bridge]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Schematic picture of the bridge
|}

Figure [[#img-10|10]] shows that the CDM detects a critical point at a load level of <math display="inline">P_{crit} = 88</math>. This bifurcation point is very close to the maximum (limit) load point <math display="inline">P_{lim} = 90</math>.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-01ffa2.png|600px|Curved bridge. Equilibrium path and critical load prediction for nondamaged material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Curved bridge. Equilibrium path and critical load prediction for nondamaged material
|}

Introducing damage (see Figure [[#img-11|11]]) the critical point is shifted to a lower load value of <math display="inline">P_{crit} = 83</math>. (<math display="inline">\log \det ({\boldsymbol K_T})</math> becomes negative at a load level of <math display="inline">P_1 = 82.8</math> and again positive at <math display="inline">P_2 = 83.4</math>). The CDM prediction path converges to an intermediate value between these  two points. The initial critical load prediction from the the undeformed configuration is <math display="inline">P_c = 60</math> which is <math display="inline">\sim 38 % </math> deviated from the actual critical load value. Progressing on the equilibrium path the predictions of the critical load and the critical displacement become more exact.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-fd728b.png|600px|Curved bridge with damage: r₀= 0.5, H = 5]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Curved bridge with damage: <math>r_0 = 0.5</math>, <math>H = 5</math>
|}

===5.3 Shallow circular arc===

The shallow arc of Figure [[#img-12|12]] has been studied with 2D solid elements. For symmetry reasons only one half of the arc was discretised using 10 eigth node solid elements and 53 nodes.  A point load acts vertically in the middle of the arc. Both ends of the arc are clamped. The diagrams of Figures 13&#8211;16 show a plot of the negative vertical displacement of the upper node in the apex of the arc versus the applied load.

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-pic15.png|600px|Geometry of the arc structure]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Geometry of the arc structure
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-fcad2f.png|600px|Shallow circular arc. Equilibrium path and critical load prediction for non damaged material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Shallow circular arc. Equilibrium path and critical load prediction for non damaged material
|}

The equilibrium path in Figure [[#img-13|13]] shows a limit point at a load level of <math display="inline">P_{crit}= 161</math>. After the two first predictions the CDM approximates this value almost exactly.

A characteristic phenomena of all the diagrams with damage (Figures [[#img-16|16]] - [[#img-14|14]]) is that the critical point is situated on the declining path after the maximum load has been reached. The value of the current stiffness parameter in these points indicates a critical point corresponding to a singularity of the tangent stiffness matrix.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-5c4599.png|600px|Shallow circular arc with damage: r₀= 0.3, H = 0.5]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Shallow circular arc with damage: <math>r_0 = 0.3</math>, <math>H = 0.5</math>
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-bdc55f.png|600px|Shallow circular arc with damage: r₀= 0.2, H = 1]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Shallow circular arc with damage: <math>r_0 = 0.2</math>, <math>H = 1</math>
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-ddcdf9.png|600px|Shallow circular arc with damage: r₀= 0.2, H = 0.5]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Shallow circular arc with damage: <math>r_0 = 0.2</math>, <math>H = 0.5</math>
|}

The CDM results confirm the critical points identified by standard limit point procedures. A comparison with the undamaged predictions demonstrate the influence  of the damage model in the CDM in the early states of the equilibrium path.

===5.4 Hollow arc===

The last example is a hollow arc with a radius of <math display="inline">R=100</math> and a span of <math display="inline">S=160</math>, analyzed using 10 eight node solid elements. The exact geometrical and material data can be seen in Figure [[#img-17|17]]. A vertical load is placed in the apex of the arc.

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-pic16.png|600px|Geometry of the hollow arc]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Geometry of the hollow arc
|}

In the diagrams the vertical displacements of the lower apex node versus the applied load is plotted. Figure [[#img-18|18]] shows results for the non damaged model. A bifurcation point is found at a load level of <math display="inline">P_{bif}=1070</math> and a limit load point at <math display="inline">P_{limit}=1263</math>. The first one is predicted quite accurately by the CDM even for the initial undeformed configuration (<math display="inline">\sim 25 %</math> error).

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-diag121.png|600px|Hollow arc. Equilibrium path and critical load predictions for non damaged material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Hollow arc. Equilibrium path and critical load predictions for non damaged material
|}
<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-80689a.png|600px|Hollow arc. Equilibrium path and critical load predictions for non damaged material]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Hollow arc. Equilibrium path and critical load predictions for non damaged material
|}

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-cffac3.png|600px|Hollow arc with damage: r₀= 8., H = 1]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' Hollow arc with damage: <math>r_0 = 8.</math>, <math>H = 1</math>
|}

<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_863111864-test-picture-1a8d00.png|600px|Hollow arc with damage: r₀= 8., H = 0.5]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Hollow arc with damage: <math>r_0 = 8.</math>, <math>H = 0.5</math>
|}

Results accounting for the effect of damage (Figures&nbsp;19 and 20) show that the critical point again lies on the declining path after the maximum load. The behavior of the structure has changed and the limit and bifurcation points are closer together now. The critical load path shows a deviation when damage occurs. Note again that the CDM yielded the correct critical load value in all cases studied.

==6 Conclusions==

In this paper the extension of the Critical Displacement Method with a damage constitutive model for structural instability analysis using truss and solid elements was presented.

The results using the CDM were quite good for the undamaged and damaged material models and the critical point was always encountered with high precision. The values of the critical load predicted from the undeformed configuration at the cost of a single eignvalue analysis were quite accurate in most problems, with <math display="inline">40 %</math> deviation in the worst case. This value can be a useful first estimate of the instability load for practical purposes.

==7 Acknowledgment==

The authors acknowledge the help of Dr. W.T. Matias in providing advice and the initial code for implementation of the damage model in the CDM.

==REFERENCES==

<div id="cite-1"></div>
[[#citeF-1|[1]]]  Bathe, K. J.: ''Finite element procedures'', Prentice Hall, Englewood Cliffs, New Jersey, 1996 <div id="cite-2"></div>
[2]  Bertsekas, D. P.: ''Constrained optimisation and Lagrange multiplier methods'', Academic Press, 1982 <div id="cite-3"></div>
[[#citeF-3|[3]]]  Choong, K. K., Hangai, Y.: “Review on methods of bifurcation analysis for geometrically nonlinear structures”, ''Bulletin of IASS'', Vol. 34, no. 112, p. 133-149 <div id="cite-4"></div>
[[#citeF-4|[4]]]  Crisfield, M. A.: ''Non&#8211;linear finite element analysis of solids and structures, Vol1'', John Wiley & Sons, 1991 <div id="cite-5"></div>
[5]  Dawe, D.J.: “Finite deflection analysis of shallow arches by discrete element method”, ''International Journal of Numerical Methods in Engineering'', Vol. 3 (1971) <div id="cite-6"></div>
[[#citeF-6|[6]]]  Kleiber, M.: ''Incremental Finite Element modeling in non-linear mechanics of solids and structures'', Ellis Horwood Limited, Chichester <div id="cite-7"></div>
[7]  Lueneberger, D. G.: ''Linear and nonlinear programming'', Addison&#8211;Wesley, Massachusetts, 1984 <div id="cite-8"></div>
[[#citeF-8|[8]]]  Matias, W.T.: ''The critical displacement method for structural instability analysis'' (in Spanish), PhD Thesis, Universitat Politecnica de Catalunya, 1996 <div id="cite-9"></div>
[[#citeF-9|[9]]]  Oñate, E.:  “Reliability analysis of concrete structures and numerical experimental studies”, ''Research Report'' no. 107, CIMNE, Barcelona, 1997 <div id="cite-10"></div>
[[#citeF-10|[10]]]  Oliver, J., Cervera, M., Oller, S., Lubliner, J.: “Isotropic damage models and smeared crack analysis of concrete”, ''Second International Conference on Computer Aided Analysis and Design of Concrete Structures'', Zell am See (Austria), 1990, Vol. 2, p.945-958 <div id="cite-11"></div>
[[#citeF-11|[11]]]  Oñate, E., Oliver, J., Miquel, J., Suárez, B.: “A finite element formulation for geometrically nonlinear problems using a secant matrix”, ''Proceedings of International Conference on Computational Mechanics'', Tokyo May 25-29,  Springer Verlag (1986) <div id="cite-12"></div>
[[#citeF-12|[12]]]  Oñate, E.:  “Possibilities of the secant stiffness matrix of non linear finite element analysis”, ''Non Linear Engineering Computations'', N. Bicanic ''et al.'' (Ed.), Pineridge Press, 1991 <div id="cite-13"></div>
[[#citeF-13|[13]]]  Oñate, E.:  “On the derivation and possibilities of the secant stiffness matrix for non linear finite element analysis”, ''Computational Mechanics'', Vol.15, p. 572-593, 1995 <div id="cite-14"></div>
[[#citeF-14|[14]]]  Oñate, E., Matias, W. T.: “Enhanced prediction of structural instability  points using a critical displacement method”, ''Advances in finite element technology'', N.&#8211;E. Wiberg (Eds.), CIMNE, Barcelona, 1995 <div id="cite-15"></div>
[[#citeF-15|[15]]]  Oñate, E., Matias, W. T.: “A critical displacement approach for predicting structural instability”,   ''Computer Methods in Applied Mechanics and Engineering'', 134 (1996), p. 135-161 <div id="cite-16"></div>
[[#citeF-16|[16]]]  Moran, A., Oñate, E. and Miquel, J., “A general procedure for deriving symmetric expressions for the secant and tangent stiffness matrices in finite element analysis”, ''International Journal for Numerical Methods in Engineering'', Vol. 42 , pp. 219&#8211;236, 1998 <div id="cite-17"></div>
[[#citeF-17|[17]]]  Riks, E.: “An incremental approach to the solution of snapping and buckling problems”, in ''International Journal of Solids and Structures'', Vol. 15 (1979), p. 529-551 <div id="cite-18"></div>
[18]  Riks, E.: “Some computational aspects of the stability analysis of nonlinear structures”, in ''Computer Methods in Applied Mechanics and Engineering''Vol. 47 (1984), p. 219-259 <div id="cite-19"></div>
[19]  Riks, E.: “Bifurcation and stability, a numerical approach”, in ''Proceedings on Innovaive Methods for Nonlinear Problems'', Eds. Liu, Belytschko  and Park, Pineridge Presss, p. 313-344, 1984 <div id="cite-20"></div>
[20]  Seydel, R.: “Numerical computation of branch points in nonlinear structural analysis”, ''Num. Math.'', Vol. 33 (1979), p. 339-352 <div id="cite-21"></div>
[21]  Skeie, G., Felippa, C.A.: “Detecting and traversing bifurcation points in nonlinear structural analysis”, ''International Journal of Space Structures'', Vol. 6 (1991), p. 77-98 <div id="cite-22"></div>
[22]  Waszczysztn, Z.: “Numerical problems of nonlinear stability analysis of elastic structures”, ''Computers and Structures'', Vol. 17 (1983), p. 13-24 <div id="cite-23"></div>
[23]  Wagner, W.:  ''Zur Behandlung von Stabilitätsproblemen der Elastostatik mit der Methode der Finiten Elemente'', Habilitation Thesis, IBNM Universität Hannover, 1990 <div id="cite-24"></div>
[[#citeF-24|[24]]]  Wriggers, P., Wagner, W., Miehe, C.:  “A quadratically convergent procedure for the calculation of stability points in finite element analysis”, ''Computer Methods in Applied Mechanics and Engineering'', 70 (1988), p. 329-347 <div id="cite-25"></div>
[[#citeF-25|[25]]]  Wriggers, P., Simo, J.C.: “A general procedure for the direct computation of turning and bifurcation points”, ''International Journal for Numerical Methods in Engineering'', Vol. 30 (1990), p. 155-176 <div id="cite-26"></div>
[[#citeF-26|[26]]]  Zienkiewicz, O. C., Taylor:  ''The Finite Element Method'',Vol. I (1989) and Vol. II (1991), McGraw&#8211;Hill