Published in Engineering Computations Vol. 18 (3-4), pp. 642-662, 2001
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
Mathematical phenomena such as bifurcation and limit points, commonly known as critical points, have received considerable attention in the computational structural community , , , –. 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 . In direct methods the condition for occurrence of a critical point is included in the system to be solved ,  . Arc–length methods, essentially a particular class of continuation methods, provide a useful tool for the tracking of the equilibrium path .
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 , , . 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 –.
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.
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 is considered as a starting point (Figure 1). The body is submitted to different load levels, here denoted by the superscript , in a quasi-static analysis. Starting from a known actual configuration the solution for an updated configuration at load level is sought. Furthermore, the body is described in a generalized Lagrange formulation with respect to an arbitrary reference configuration at load level . This can be made coincident with the initial configuration (total Lagrange approach) or the current configuration (updated Lagrange approach), , . The configuration at the time corresponds to a variation of the known configuration at the time .
|Figure 1: Deformation process of a solid|
The displacements at the time are written as the sum of the known displacements and the sought displacement increments , i.e.
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
In above and contain terms with linear and quadratic dependence on the incremental displacements , respectively, i.e.
A comma in a lower right index denotes a derivative with respect to the Cartesian coordinate referred to by the index following ( with ). A similar expression can be deduced for the virtual displacements. Taking into account the fact that leads to the following expressions:
The first and second order virtual strain increments are
The Piola-Kirchoff stress vector at is obtained incrementally by
The constitutive equation chosen assumes a linear dependence between the second Piola-Kirchhoff stress increments and the Green-Lagrange strain increments as
where is a fourth order constitutive tensor.
The expression of the PVW at time can be written as
where are the second Piola–Kirchhoff stress and , denote the volume and surface forces, respectively. As usual and denote the volume and the boundary of the reference configuration.
As a next step the displacements and their increments are discretised in the standard finite element manner 
where and are the shape function matrix and the nodal displacement vector, respectively.
where the incremental secant stiffness matrix is given by
Note that the secant matrix is split in terms according to their dependence on the incremental nodal displacements . The form of the different matrices in equation (17) for truss and solid elements con be found in ,  and  .
In eq.(16) is the residual force vector given by
In eq.(16) is the equivalent nodal force vector deduced from eq.(14) and is the first order strain matrix relating first order strain increments and nodal displacements . The expression of this matrix can be found in  and .
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 , . An alternative procedure for deriving secant stiffness expressions can be found in .
Note that a “total” secant stiffness matrix can be obtained by making . In a total Lagrangean description (i.e. ) the total secant matrix is given by 
|Figure 2: Total (left) and incremental (right) secant approaches|
It is important to note that in the limiting case of the secant stiffness matrix reduces to the tangent stiffness matrix, i.e.
Using 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
In order to trace the equilibrium path the load terms have to be scalable. For this purpose, the scalar load parameter is introduced. This parameter leaves the pattern of the external loads unchanged, only the size is varied (i.e. ).
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 , which determines the load parameter . A variety of different control equations can be found in the literature , . The ones used here are stated below in equation (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.
The main criterion for existance of a critical point is that the tangent stiffness matrix becomes singular at the point. Mathematically this means that equation (23) below has nontrivial solutions.
where 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 has some further consequences. The determinant of is equal to in the critical point, furthermore the smallest eigenvalue is . This can be stated mathematically as
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 , the curve should become at the critical point according to eq.(24a).
Furthermore eq.(24b) indicates that the smallest eigenvalue turns zero at . Keeping track of the number of negative diagonal elements of 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
The estimated critical load is then obtained as a multiple of the known load at the equilibrium point.
The approach termed as the Critical Displacement Method (CDM) , – 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 is expressed as
where are the nodal displacements at the actual equilibrium configuration from which the prediction has to be made and is an estimate of the critical displacement increment pattern yielding structural instability at . It is now assumed to be of the form with being an approximation for the critical displacement pattern and a multiplier. As a convenient starting value of the last known value of the displacement pattern can be chosen (i.e. ).
Inserting this approximation in the expression for the tangent stiffness matrix (20) together with the condition (24a), leads to
where the tangent matrix at the critical point has been split into terms according to their dependence on . The detailed form of the matrices involved in eq.(28) can be seen in  and . Note that the resultant equation (28) is a quadratic eigenvalue problem, that can be solved for .
Eq.(28) can be simplified to the following standard linear eigenvalue problem by neglecting second order terms in
For the solution of eq.(29) the Inverse Iteration method has proved to be the best , . The Rayleigh quotient is used here as an approximation for the first eigenvalue. The convergence of this procedure can be enhanced by shifting the terms of matrix . The Jacobi iteration method as an alternative is very time consuming and thus not very suitable for large problems.
Having obtained the smallest eigenvalue , the critical displacements are predicted from equation (27) with . For the determination of the corresponding critical loads these displacements are inserted into the incremental secant stiffness matrix (see eq.(16)), which gives the prediction for the critical load as
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:
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.
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 . Indeed more sophisticated damage models can be used .
The damage parameter is a measure of the loss of secant stiffness of the material. It ranges from (undamaged material) to (fully de-gradated material). The constitutive law is now written as
where 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 and a damage threshold value as .
In our work we have taken
In eq.(35) is the initial damage tdshreshold value and the hardening modulus. The stress norm 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 is then updated according to the following evolution law
which can be deduced from (35) and the damage criterion (33) when the equality sign is applied. In a damage case is no longer constant, which affects the tangent stiffness matrix . Using equations (32) and (35) the tangent constitutive matrix is obtained as
The damage model is described in an algorithmic form.
Note that initially the trial stresses are computed and then compared with the actual threshold value . Following this comparison the values and are updated (damage, if–case) or not (undamaged, else–case). Index indicates here the different load levels.
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 in order to solve equation (29). Index 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 . Calculating first the critical Green stresses , the stress norm is obtained as
The critical damage parameter is now chosen according to the following scheme:
The option 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 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.
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.
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 . A unit load 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
|Figure 4: Discretization of the arc structure into tow node bar element|
Figure 5 shows the equilibrium path giving a limit load and a bifurcation load of . 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 ( 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: ,|
|Figure 7: Arc with damage: ,|
|Figure 8: Arc with damage: ,|
The plots for different damage threshold values and hardening modulus 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.
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 is used. In all figures the vertical displacements of node ten multiplied by 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 . This bifurcation point is very close to the maximum (limit) load point .
|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 . becomes negative at a load level of and again positive at ). The CDM prediction path converges to an intermediate value between these two points. The initial critical load prediction from the the undeformed configuration is which is 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: ,|
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 . 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: ,|
|Figure 15: Shallow circular arc with damage: ,|
|Figure 16: Shallow circular arc with damage: ,|
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.
The last example is a hollow arc with a radius of and a span of , 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 and a limit load point at . The first one is predicted quite accurately by the CDM even for the initial undeformed configuration ( 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: ,|
|Figure 21: Hollow arc with damage: ,|
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.
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 deviation in the worst case. This value can be a useful first estimate of the instability load for practical purposes.
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.
 Zienkiewicz, O. C., Taylor: The Finite Element Method,Vol. I (1989) and Vol. II (1991), McGraw–Hill