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 [3], [4], [6], [17]–[25]. In a practical context the occurrence of those points in equilibrium paths of load–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 [3]. In direct methods the condition for occurrence of a critical point is included in the system to be solved [24], [25] . Arc–length methods, essentially a particular class of continuation methods, provide a useful tool for the tracking of the equilibrium path [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 [8], [14], [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–displacement stiffness relationship the critical load is finally obtained [12]–[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 [8], [14], [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 ${\textstyle {}^{0}V}$ is considered as a starting point (Figure 1). The body is submitted to different load levels, here denoted by the superscript ${\textstyle t}$, in a quasi-static analysis. Starting from a known actual configuration ${\textstyle {}^{t}V}$ the solution for an updated configuration at load level ${\textstyle t+\Delta t}$ is sought. Furthermore, the body is described in a generalized Lagrange formulation with respect to an arbitrary reference configuration at load level ${\textstyle {}^{r}V}$. This can be made coincident with the initial configuration ${\textstyle {}^{r}V={}^{0}V}$ (total Lagrange approach) or the current configuration ${\textstyle {}^{r}V={}^{t}V}$ (updated Lagrange approach)[1], [6], [13]. The configuration at the time ${\textstyle t+\Delta t}$ corresponds to a variation of the known configuration at the time ${\textstyle t}$.

Figure 1: Deformation process of a solid

The displacements at the time ${\textstyle t+\Delta t}$ are written as the sum of the known displacements ${\textstyle {{\boldsymbol {}}^{t}u}}$ and the sought displacement increments ${\textstyle \Delta {\boldsymbol {u}}}$, i.e.

${}^{t+\Delta t}{\boldsymbol {u}}\;=\;{}^{t}{\boldsymbol {u}}+\Delta {\boldsymbol {u}}$

(1)

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

${}_{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)$

(2)

${}_{r}^{t+\Delta t}\Delta \epsilon _{ij}\;=\;{}_{r}^{t+\Delta t}\epsilon _{ij}-{}_{r}^{t}\epsilon _{ij}={}_{r}e_{ij}+{}_{r}\eta _{ij}$

(3)

In above ${\textstyle {}_{r}e_{ij}}$ and ${\textstyle {}_{r}\eta _{ij}}$ contain terms with linear and quadratic dependence on the incremental displacements ${\textstyle \Delta {\boldsymbol {u}}}$, respectively, i.e.

${}_{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)$ (4) ${}_{r}\eta _{ij}\;=\;{\frac {1}{2}}{}_{r}\Delta u_{k,i}\;{}_{r}\Delta u_{k,j}$ (5)

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

$\delta {}^{t+\Delta t}u_{i}\;=\;\delta \Delta u_{i}$ (6) $\delta {}_{r}^{t+\Delta t}\Delta \epsilon _{ij}\;=\;\delta {}_{r}e_{ij}+\delta {}_{r}\eta _{ij}$ (7)

The first and second order virtual strain increments are

$\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)$ (8) $\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)$ (9)

The Piola-Kirchoff stress vector at ${\textstyle t+\Delta t}$ is obtained incrementally by

${}_{r}^{t+\Delta t}\sigma \;=\;{}_{r}^{t}\sigma +{}_{r}\Delta \sigma$

(10)

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

${}_{r}\Delta \sigma \;=\;{}_{r}{\boldsymbol {D}}({}_{r}{\boldsymbol {e}}+{}_{r}\eta )$

(11)

where ${\textstyle {}_{r}{\boldsymbol {D}}}$ is a fourth order constitutive tensor.

The expression of the PVW at time ${\textstyle t+\Delta t}$ can be written as

$\int _{{}^{r}V}\delta {}_{r}^{t+\Delta t}\epsilon _{ij}\;{}_{r}^{t+\Delta t}\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$

(12)

where ${\textstyle {}_{r}^{t+\Delta t}\sigma _{ij}}$ are the second Piola–Kirchhoff stress and ${\textstyle f_{i}^{B}}$, ${\textstyle f_{i}^{S}}$ denote the volume and surface forces, respectively. As usual ${\textstyle {}^{r}V}$ and ${\textstyle {}^{r}S}$ denote the volume and the boundary of the reference configuration.

Inserting equations (3) and (11) into (12) gives the full incremental form of the PVW as

$\int _{{}^{r}V}[\delta {}_{r}e_{ij}\;{}_{r}D_{ijkl}\;{}_{r}e_{kl}+\delta {}_{r}e_{ij}\;{}_{r}D_{ijkl}\;{}_{r}\eta _{kl}+\delta {}_{r}\eta _{ij}\;{}_{r}D_{ijkl}\;{}_{r}e_{kl}+\delta {}_{r}\eta _{ij}\;{}_{r}D_{ijkl}\;{}_{r}\eta _{kl}+$ $+\delta {}_{r}\eta _{ij}\;{}_{r}^{t}S_{ij}]dV\;=\;{}_{r}^{t+\Delta t}f-\int _{{}^{r}V}\delta {}_{r}e_{ij}\;{}_{r}^{t}\sigma _{ij}dV$

(13)

where

${}_{r}^{t+\Delta t}f\;=\;\int _{{}^{r}V}\delta {}^{t+\Delta t}u_{i}{}_{0}^{t+\Delta t}f_{i}^{B}dV+\int _{{}^{r}S}\delta {}^{t+\Delta t}u_{i}{}^{t+\Delta t}f_{i}^{S}dS$

(14)

Eq.(13) is the point of departure for the derivation of the secant stiffness matrix [8], [11]–[13].

As a next step the displacements and their increments are discretised in the standard finite element manner [26]

${}^{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}})$

(15)

where ${\textstyle {\boldsymbol {N}}}$ and ${\textstyle {\boldsymbol {a}}}$ are the shape function matrix and the nodal displacement vector, respectively.

The incremental secant stiffness expression is now obtained by substituting equations (4), (5), (8), (9), (10), (11) and (15) into (13) and collecting the terms depending on ${\textstyle {{\boldsymbol {\Delta }}a}}$ in the following form:

${}_{r}^{t}{\boldsymbol {K}}_{S}({}^{t}{\boldsymbol {a}},\Delta {\boldsymbol {a}})\;{{\boldsymbol {\Delta }}a}\;=\;{}_{r}^{t+\Delta t}{\boldsymbol {r}}$

(16)

where the incremental secant stiffness matrix is given by

${}_{r}^{t}{\boldsymbol {K}}_{S}(\Delta {\boldsymbol {a}})\;={}_{r}^{t}{\boldsymbol {K}}_{L}+{}_{r}^{t}{\boldsymbol {K}}_{M}(\Delta {\boldsymbol {a}})+{}_{r}^{t}{\boldsymbol {K}}_{N}(\Delta {\boldsymbol {a}}^{2})+{}_{r}^{t}{\boldsymbol {K}}_{\sigma }({\boldsymbol {\sigma }})$

(17)

Note that the secant matrix ${\textstyle {}_{r}^{t}{{\boldsymbol {K}}_{S}}}$ is split in terms according to their dependence on the incremental nodal displacements ${\textstyle \Delta {\boldsymbol {a}}}$. The form of the different matrices in equation (17) for truss and solid elements con be found in [8], [13] and [15] .

In eq.(16) ${\textstyle {}^{t+\Delta t}{\boldsymbol {r}}}$ is the residual force vector given by

${}_{r}^{t+\Delta t}{\boldsymbol {r}}\;=\;{}^{t+\Delta t}{\boldsymbol {f}}-\int _{{}^{r}V}{}_{r}^{t}{\boldsymbol {B}}_{L}^{T}{}_{r}^{t}\sigma _{ij}$

(18)

In eq.(16) ${\textstyle {}^{t+\Delta t}{\boldsymbol {f}}}$ is the equivalent nodal force vector deduced from eq.(14) and ${\textstyle {}_{r}^{t}{\boldsymbol {B}}_{L}}$ is the first order strain matrix relating first order strain increments ${\textstyle {}_{r}e_{ij}}$ and nodal displacements ${\textstyle \Delta u_{i}}$. The expression of this matrix can be found in [13] and [15].

Equation (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–displacement pattern [13], [14]. An alternative procedure for deriving secant stiffness expressions can be found in [16].

Note that a “total” secant stiffness matrix can be obtained by making ${\textstyle \Delta {\boldsymbol {a}}={\boldsymbol {a}}^{t}}$. In a total Lagrangean description (i.e. ${\textstyle {}^{r}V={}^{0}V}$) the total secant matrix is given by [13]

${}_{0}^{t}{\boldsymbol {K}}_{S}={}_{0}^{t}{\boldsymbol {K}}_{L}+{}_{0}^{t}{\boldsymbol {K}}_{M}({}^{t}{\boldsymbol {a}})+{}_{0}^{t}{\boldsymbol {K}}_{N}({}^{t}{\boldsymbol {a}}^{2})$

(19)

Figure 2 shows schematically the difference between the incremental and total secant approaches. Further details can be found in [8], [13]–[15].

Figure 2: Total (left) and incremental (right) secant approaches

It is important to note that in the limiting case of ${\textstyle {{\boldsymbol {\Delta }}a}\rightarrow 0}$ the secant stiffness matrix reduces to the tangent stiffness matrix, i.e.

${}_{r}^{t}{\boldsymbol {K}}_{T}=\lim _{\Delta {\boldsymbol {a}}\rightarrow 0}{}_{r}^{t}{\boldsymbol {K}}_{S}={}_{r}^{t}{\boldsymbol {K}}_{L}+{}_{r}^{t}{\boldsymbol {K}}_{\sigma }$

(20)

Using ${\textstyle {}_{r}^{t}{\boldsymbol {K}}_{T}}$ the standard iteration techniques for the computation of stability points in structures can be obtained. For instance, the standard Newton–Raphson method is written from eq.(16) as

${}_{r}^{t}{\boldsymbol {K}}_{T}({}^{t}{\boldsymbol {a}})\Delta {\boldsymbol {a}}\;={}_{r}^{t+\Delta t}{\boldsymbol {r}}$

(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 [8], [13] and [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 ${\textstyle {}^{t}\lambda }$ is introduced. This parameter leaves the pattern of the external loads unchanged, only the size is varied (i.e. ${\textstyle {}^{t}{\boldsymbol {f}}={}^{t}\lambda {\boldsymbol {f}}_{0}}$).

The arc–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 ${\textstyle f({}^{t}\lambda ,{{\boldsymbol {}}^{t}u})=0}$, which determines the load parameter ${\textstyle {}^{t}\lambda }$. A variety of different control equations can be found in the literature [4], [6]. The ones used here are stated below in equation (22).

$f({}^{t}\lambda ,{{\boldsymbol {}}^{t}u})\;={\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}}$

(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 ${\textstyle {\boldsymbol {K}}_{T}}$ becomes singular at the point. Mathematically this means that equation (23) below has nontrivial solutions.

${\boldsymbol {K}}_{T}({\boldsymbol {a}}_{c})\Delta {\boldsymbol {a}}\;=\;0$

(23)

where ${\textstyle {\boldsymbol {a}}_{c}}$ 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 ${\textstyle {\boldsymbol {K}}_{T}}$ has some further consequences. The determinant of ${\textstyle {\boldsymbol {K}}_{T}}$ is equal to ${\textstyle 0}$ in the critical point, furthermore the smallest eigenvalue is ${\textstyle 0}$. This can be stated mathematically as

$\det {\boldsymbol {K}}_{T}({\boldsymbol {a}}_{c})\;=\;0$

(24a)

$\left[{\boldsymbol {K}}_{T}({\boldsymbol {a}}_{c})-\omega {\boldsymbol {I}}\right]\Phi \;=\;0\qquad {\mbox{with}}\quad \omega _{1}=0$

(24b)

Many indirect methods for the computation of critical points use condition (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 ${\textstyle det{\boldsymbol {K}}_{T}}$, the curve should become ${\textstyle 0}$ at the critical point according to eq.(24a).

Furthermore eq.(24b) indicates that the smallest eigenvalue turns zero at ${\textstyle {\boldsymbol {a}}_{c}}$. Keeping track of the number of negative diagonal elements of ${\textstyle {\boldsymbol {K}}_{T}}$ 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

$\left[{}_{r}^{t}{\boldsymbol {K}}_{L}+{\lambda }\;{}_{r}^{t}{\boldsymbol {K}}_{\sigma }\right]{\boldsymbol {a}}\;=\;0$

(25)

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

${\boldsymbol {f}}_{c}\;=\;{\lambda }\;{}^{t}{\boldsymbol {f}}$

(26)

The approach termed as the Critical Displacement Method (CDM) [8], [14]–[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 ${\textstyle {\boldsymbol {a}}_{c}}$ is expressed as

${}^{t+\Delta t_{c}}{\boldsymbol {a}}_{c}\;=\;{}^{t}{\boldsymbol {a}}+\Delta {\boldsymbol {a}}_{c}$

(27)

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

Inserting this approximation in the expression for the tangent stiffness matrix (20) together with the condition (24a), leads to

${det}\left[{\boldsymbol {K}}_{T}({}^{t}{\boldsymbol {a}})+\rho {\boldsymbol {K}}_{1}(\Phi )+\rho ^{2}{\boldsymbol {K}}_{2}(\Phi ^{2})\right]=0$

(28)

where the tangent matrix at the critical point has been split into terms according to their dependence on ${\textstyle \Phi }$. The detailed form of the matrices involved in eq.(28) can be seen in [8] and [15]. Note that the resultant equation (28) is a quadratic eigenvalue problem, that can be solved for ${\textstyle \rho }$.

Eq.(28) can be simplified to the following standard linear eigenvalue problem by neglecting second order terms in ${\textstyle \rho }$

${det}\left[{\boldsymbol {K}}_{T}({}^{t}{\boldsymbol {a}})+\rho {\boldsymbol {K}}_{1}(\Phi )\right]=0$

(29)

For the solution of eq.(29) the Inverse Iteration method has proved to be the best [1], [4]. The Rayleigh quotient ${\textstyle \rho }$ is used here as an approximation for the first eigenvalue. The convergence of this procedure can be enhanced by shifting the terms of matrix ${\textstyle {\boldsymbol {K}}_{T}}$. The Jacobi iteration method as an alternative is very time consuming and thus not very suitable for large problems.

Having obtained the smallest eigenvalue ${\textstyle \rho }$, the critical displacements are predicted from equation (27) with ${\textstyle {{\boldsymbol {\Delta }}a}_{c}=\rho {\boldsymbol {\Phi }}}$. For the determination of the corresponding critical loads ${\textstyle {\boldsymbol {f}}_{c}}$ these displacements are inserted into the incremental secant stiffness matrix (see eq.(16)), which gives the prediction for the critical load as

$\Delta {\boldsymbol {f}}_{c}\;=\;{}_{r}^{t}{\boldsymbol {K}}_{S}({}^{t}{\boldsymbol {a}},\rho \Phi )\;\rho \Phi$ (30) ${}^{t_{c}}{\boldsymbol {f}}\;=\;{}^{t}{\boldsymbol {f}}+\Delta {\boldsymbol {f}}_{c}$ (31)

The total secant stiffness matrix of eq.(19) can be used to directly calculate ${\textstyle {}^{t_{c}}{\boldsymbol {f}}}$ as shown in Figure 2. Further details can be found in [8], [13]–[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:

1. 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.
2. Starting from the displacements ${\textstyle {}^{t}{\boldsymbol {a}}}$, the critical displacements ${\textstyle {}^{t_{c}}{\boldsymbol {a}}}$ are found by solving equation (29) and using eq.(27) with ${\textstyle \Delta {\boldsymbol {a}}_{c}=\rho ^{t}{\boldsymbol {a}}}$.
3. The corresponding critical load is obtained from equations (20) and (31)
4. A further point on the equilibrium path is computed
5. Go to 2.

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–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 [10]. Indeed more sophisticated damage models can be used [9].

The damage parameter ${\textstyle d}$ is a measure of the loss of secant stiffness of the material. It ranges from ${\textstyle 0}$ (undamaged material) to ${\textstyle 1}$ (fully de-gradated material). The constitutive law is now written as

$\sigma \;=\;(1-d){\boldsymbol {D}}\;\epsilon ={\overline {\boldsymbol {D}}}\varepsilon$

(32)

where ${\textstyle {\overline {\boldsymbol {D}}}}$ 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 ${\textstyle {\overline {\sigma }}}$ and a damage threshold value ${\textstyle r(d)}$ as .

$F(\epsilon ,d)\;=\;{\overline {\sigma }}-r(d)\leq 0$

(33)

In our work we have taken

${\overline {\sigma }}(\epsilon )\;=\;{\sqrt {\epsilon :{\boldsymbol {D}}:\epsilon }}$ (34) $r(d)\;=\;{\frac {r_{0}}{1\;-\;(1+H)d}}$ (35)

In eq.(35) ${\textstyle r_{0}}$ is the initial damage tdshreshold value and ${\textstyle H}$ the hardening modulus. The stress norm ${\textstyle {\overline {\sigma }}}$ 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–strain curve is shown in Figure 3.

Figure 3: Constitutive plot for the single parameter damage model

If inequality (33) is violated, damage occurs. The damage parameter ${\textstyle d}$ is then updated according to the following evolution law

$d\;=\;{\frac {1}{1+H}}\left(1-{\frac {r_{0}}{\overline {\sigma }}}\right),$

(36)

which can be deduced from (35) and the damage criterion (33) when the equality sign is applied. In a damage case ${\textstyle d}$ is no longer constant, which affects the tangent stiffness matrix ${\textstyle {\boldsymbol {K}}_{T}}$. Using equations (32) and (35) the tangent constitutive matrix ${\textstyle {\boldsymbol {C}}}$ is obtained as

${\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 )&forF>0\\(1-d){\boldsymbol {D}}&forF\leq 0\end{cases}}.$

(37)

The damage model is described in an algorithmic form.

Compute: ${\overline {\sigma }}^{m+1}\;=\;{\sqrt {\epsilon ^{m+1}:{\boldsymbol {D}}:\epsilon ^{m+1}}}$ if ${\overline {\sigma }}\leq r^{m}$ ${\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}}$ else ${\overline {\sigma }}>r^{m}$ ${\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}}$

Note that initially the trial stresses ${\textstyle {\overline {\sigma }}^{m+1}}$ are computed and then compared with the actual threshold value ${\textstyle r^{m}}$. Following this comparison the values ${\textstyle r^{m+1}}$ and ${\textstyle d^{m+1}}$ are updated (damage, if–case) or not (undamaged, else–case). Index ${\textstyle m}$ 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 ${\textstyle {\boldsymbol {C}}^{m+1}}$ in order to solve equation (29). Index ${\textstyle m}$ 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 ${\textstyle d_{c}}$. Calculating first the critical Green stresses $\epsilon _{c}$, the stress norm ${\textstyle {\overline {\sigma }}}$ is obtained as

${\overline {\sigma }}(\epsilon _{c})\;=\;{\sqrt {\epsilon _{c}:{\boldsymbol {D}}:\epsilon _{c}}}$

(38)

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

$d_{c}={\begin{cases}d^{m}&for{\overline {\sigma }}(\epsilon _{c})\leq 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}}$

(39)

The option ${\textstyle d_{c}=d^{m}}$ 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.(38). The value of ${\textstyle d_{c}}$ 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–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 two-node truss elements and the last two ones eight-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 4. The product of Young's modulus with the cross area of the bar is ${\textstyle EA=5.0\times 10^{3}}$. A unit load ${\textstyle P=1}$ is applied at the top node. All the displacements plotted in Figures 5-8 are the vertical displacements of node four. The values are multiplied by ${\textstyle -1}$

Figure 4: Discretization of the arc structure into tow node bar element

Figure 5 shows the equilibrium path giving a limit load ${\textstyle P_{lim}=3.7}$ and a bifurcation load of ${\textstyle P_{bif}=3.47}$. 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 ${\textstyle P_{c}=3.0}$ (${\textstyle \sim 15\%}$ error) is predicted from the undeformed configuration.

Figures 6–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.

Figure 5: Arc structure. Equilibrium path and critical load prediction for nondamaged material

Figure 6: Arc with damage: $r_{0}=0.1$, $H=1$

Figure 7: Arc with damage: $r_{0}=0.1$, $H=0.5$

Figure 8: Arc with damage: $r_{0}=0.05$, $H=1$

The plots for different damage threshold values ${\textstyle r_{0}}$ and hardening modulus ${\textstyle H}$ 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 9) The arrow at node ten symbolizes again a unit load. A value of ${\textstyle EA=3.0\times 10^{6}}$ is used. In all figures the vertical displacements of node ten multiplied by ${\textstyle -1}$ versus the applied load are plotted.

Figure 9: Schematic picture of the bridge

Figure 10 shows that the CDM detects a critical point at a load level of ${\textstyle P_{crit}=88}$. This bifurcation point is very close to the maximum (limit) load point ${\textstyle P_{lim}=90}$.

Figure 10: Curved bridge. Equilibrium path and critical load prediction for nondamaged material

Introducing damage (see Figure 11) the critical point is shifted to a lower load value of ${\textstyle P_{crit}=83}$. ${\textstyle (\log \det({{\boldsymbol {K}}_{T}})}$ becomes negative at a load level of ${\textstyle P_{1}=82.8}$ and again positive at ${\textstyle P_{2}=83.4}$). The CDM prediction path converges to an intermediate value between these two points. The initial critical load prediction from the the undeformed configuration is ${\textstyle P_{c}=60}$ which is ${\textstyle \sim 38\%}$ 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.

Figure 11: Curved bridge with damage: $r_{0}=0.5$, $H=5$

5.3 Shallow circular arc

The shallow arc of Figure 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 Figures 13–16) show a plot of the negative vertical displacement of the upper node in the apex of the arc versus the applied load.

Figure 12: Geometry of the arc structure

The equilibrium path in Figure 13 shows a limit point at a load level of ${\textstyle P_{crit}=161}$. After the two first predictions the CDM approximates this value almost exactly.

A characteristic phenomena of all the diagrams with damage (Figures 14 - 16) 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.

Figure 13: Shallow circular arc. Equilibrium path and critical load prediction for non damaged material

Figure 14: Shallow circular arc with damage: $r_{0}=0.3$, $H=0.5$

Figure 15: Shallow circular arc with damage: $r_{0}=0.2$, $H=1$

Figure 16: Shallow circular arc with damage: $r_{0}=0.2$, $H=0.5$

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 ${\textstyle R=100}$ and a span of ${\textstyle S=160}$, analyzed using 10 eight node solid elements. The exact geometrical and material data can be seen in Figure 17. A vertical load is placed in the apex of the arc.

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 18 shows results for the non damaged model. A bifurcation point is found at a load level of ${\textstyle P_{bif}=1'070}$ and a limit load point at ${\textstyle P_{limit}=1'263}$. The first one is predicted quite accurately by the CDM even for the initial undeformed configuration (${\textstyle \sim 25\%}$ error).

Figure 18: Hollow arc. Equilibrium path and critical load predictions for non damaged material

Figure 19: Hollow arc. Equilibrium path and critical load predictions for non damaged material

Figure 20: Hollow arc with damage: $r_{0}=8.$, $H=1$

Figure 21: Hollow arc with damage: $r_{0}=8.$, $H=0.5$

Results accounting for the effect of damage (Figures 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 ${\textstyle 40\%}$ 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

[1] Bathe, K. J.: Finite element procedures, Prentice Hall, Englewood Cliffs, New Jersey, 1996

[2] Bertsekas, D. P.: Constrained optimisation and Lagrange multiplier methods, Academic Press, 1982

[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

[4] Crisfield, M. A.: Non–linear finite element analysis of solids and structures, Vol1, John Wiley & Sons, 1991

[5] Dawe, D.J.: “Finite deflection analysis of shallow arches by discrete element method”, International Journal of Numerical Methods in Engineering, Vol. 3 (1971)

[6] Kleiber, M.: Incremental Finite Element modeling in non-linear mechanics of solids and structures, Ellis Horwood Limited, Chichester

[7] Lueneberger, D. G.: Linear and nonlinear programming, Addison–Wesley, Massachusetts, 1984

[8] Matias, W.T.: The critical displacement method for structural instability analysis (in Spanish), PhD Thesis, Universitat Politecnica de Catalunya, 1996

[9] Oñate, E.: “Reliability analysis of concrete structures and numerical experimental studies”, Research Report no. 107, CIMNE, Barcelona, 1997

[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

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

[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

[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

[14] Oñate, E., Matias, W. T.: “Enhanced prediction of structural instability points using a critical displacement method”, Advances in finite element technology, N.–E. Wiberg (Eds.), CIMNE, Barcelona, 1995

[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

[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–236, 1998

[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

[18] Riks, E.: “Some computational aspects of the stability analysis of nonlinear structures”, in Computer Methods in Applied Mechanics and EngineeringVol. 47 (1984), p. 219-259

[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

[20] Seydel, R.: “Numerical computation of branch points in nonlinear structural analysis”, Num. Math., Vol. 33 (1979), p. 339-352

[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

[22] Waszczysztn, Z.: “Numerical problems of nonlinear stability analysis of elastic structures”, Computers and Structures, Vol. 17 (1983), p. 13-24

[23] Wagner, W.: Zur Behandlung von Stabilitätsproblemen der Elastostatik mit der Methode der Finiten Elemente, Habilitation Thesis, IBNM Universität Hannover, 1990

[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

[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

[26] Zienkiewicz, O. C., Taylor: The Finite Element Method,Vol. I (1989) and Vol. II (1991), McGraw–Hill