Subscribe Now Subscribe Today
Research Article

Development of a Finite Element System for Vibration Based Damage Development of a Finite Element System for Vibration Based Damage Identification in Structures

Ashutosh Bagchi, Jag Humar and Ahmed Noman
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

Vibration based damage identification in a structure is based on the observation that damage causes changes in its dynamic characteristics, such as, frequencies and mode shapes. There are a number of Vibration Based Damage Identification (VBDI) algorithms which utilize change in vibration characteristics in a structure to determine the location and severity of possible damage. This article presents the development of a Finite Element (FE) program that implements the three dimensional beam and plate-shell elements and some of the vibration based damage identification algorithms. The FE system developed herein is also capable of performing model updating using the measured modal parameters such the natural frequencies and mode shape vectors. A three span continuous bridge in Canada has been presented as a case study. The three dimensional FE model of the bridge was correlated to the actual structure by updating the model with modal test data and the updated model is used for damage identification using simulated damage scenarios. The finite element system developed herein has been found to be robust and the VBDI algorithms, in spite of their sensitivity to incomplete modes and measurement noise, provide important clues about the location and severity of damage in a structure.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Ashutosh Bagchi, Jag Humar and Ahmed Noman , 2007. Development of a Finite Element System for Vibration Based Damage Development of a Finite Element System for Vibration Based Damage Identification in Structures. Journal of Applied Sciences, 7: 2404-2413.

DOI: 10.3923/jas.2007.2404.2413



Vibration based damage identification in a structure is based on the observation that damage causes changes in the dynamic characteristics, such as, natural frequencies and mode shapes of the structure. There are a number of Vibration Based Damage Identification (VBDI) algorithms which utilize the change in vibration characteristics of a structure to determine the location and severity of possible damage. Commercial Finite Element Analysis (FEA) packages do not provide enough flexibility to customize and incorporate the VBDI algorithms. This paper presents the development of a FEA system which implements some of the VBDI algorithms.

There are many challenges in the application of VBDI techniques as the measured values of the natural frequencies and mode shapes of a structure often contain random measurement errors. Finite element model of a civil engineering structure, such as a bridge or a dam, often contains a large number of degrees of freedom. However, only a limited number of sensors could be used for measuring the mode shapes which lead to an incomplete definition of a mode shape. It is sometimes possible to obtain the complete mode shapes from the incomplete ones obtained from field measurements, using some kind of modal expansion technique. This introduces additional error to the mode shapes. The incompleteness of mode shapes and random measurement errors degrade the performance of VBDI algorithms. Environmental effects on the structure’s behavior also affect the performance of damage detection algorithms.

The FEA system developed here is designed to construct finite element model of a structure using three dimensional beam and shell elements, correlate the model based on measured vibration quantities and simulate various damage scenarios. Based on a given distribution of the sensors in the structure, the performance of the damage detection algorithms implemented in the software can be studied by considering the incomplete modes of the damaged structure and simulated random measurement errors. Further work is necessary to model environmental effects.

The FEA system presented here has been developed in the MATLAB (2006) environment, which is a powerful tool for mathematical programming and quick prototyping. MATLAB is chosen here to take advantage of the availability of the functions for matrix operations, optimization and graphics. However, the programming languages provided with the MATLAB environment is interpretive in nature and it slows down the system when the problem size is very large (Bagchi, 2005).

The study also presents a study on a three span continuous steel free deck bridge located in Calgary, Canada. The bridge, opened to traffic in 1997, was constructed under the technical leadership of ISIS Canada Research Network and a number of sensors were installed in the bridge. Although the present study focuses on computer simulation, some data from the dynamic tests conducted on the bridge are also used.


A number of techniques for vibration based damage identification are available in the literature (Humar et al., 2006). They include (a) methods based on frequency changes, (b) methods based on mode shape changes, (c) mode shape curvature method, (d) methods based on change in flexibility matrix, (e) methods based on changes in uniform flexibility shape curvature, (f) strain energy based damage index method, (g) method based on modal residual vector, (h) matrix update methods and (i) semi analytical techniques, such as, neural network, genetic algorithm and statistical methods. Here, the strain energy based damage index method and the matrix update methods are used for their relative effectiveness as explained in Humar et al. (2003 and 2006).

Method based on modal strain energy: This method as originally developed by Stubbs et al. (1995) is applicable to beam type structures and is based on the comparison of modal strain energy before and after damage. The method can be extended to a structure of a general type. For the healthy structure the modal strain energy of member j in mode I is given by φTikjφi and the total strain energy of the structure deforming in its mode I is obtained from


where kj is the stiffness matrix of member j and K is the stiffness matrix of the structure. Similar expressions can be formed for the damaged structure. However, in deriving the strain energy of the damaged structure, Sdi, the global stiffness matrix of the damaged structure Kd is taken as being approximately equal to K. This is based on the reasoning that in practice only a few members are likely to suffer damage. The strain energy fractions of member j in mode I for the undamaged and damaged states, fij and fdij, respectively, are given by:


The ratio of the strain energy fractions as expressed in Eq. 3 is defined as the damage index for the jth element in mode I. When the information available from nm measured modes is used the damage index is given by


Members for which γj is relatively large are likely to have been damaged. This provides the location of damage. Numerical problems may arise in the evaluation of Eq. 3 when the denominator is very small in magnitude, which may be the case when the strain energy contributed by the jth member in the modes being considered is very small. In such a case a modified damage index, , for element j, as shown in Eq. 4 can be defined.


Matrix update methods: The matrix update methods constitute the largest class of methods developed for the identification of damage from measured vibration properties. The methods are based on the determination of perturbations in the property matrices caused by damage. The property matrices include the stiffness matrix K and the mass matrix M. The method as developed by Kabe (1985) has been briefly described here. For the undamaged structure the property matrices satisfy the following eigenvalue equation


where φi is the ith mass-orthonormal mode shape and λi is the associated eigenvalue (squared frequency). Damage in the structure may change both the stiffness and the mass matrices, altering the frequencies as well as the mode shapes. The eigenvalue equation for the damaged structure is given by:


Assuming that damage does not alter the mass matrix and using Eq. 6, we get:


where φdi = φi + δφi is the ith mode shape of the damaged structure. Pre-multiplying both sides of Eq.7 by φTdi and using the transpose of Eq. 7 we get:


Perturbation matrix δK can be expressed as


where βj is the damage factor for element j. Using Eq.8 and 9 and assuming that φdi contain mass ortho-normalized mode shape vectors, we get



where ne is the number of members and D is an m by ne matrix, β is the ne-vector of the unknown changes in element stiffness matrices and δλ is the m-vector of measured eigenvalue changes.

Solution of Eq. 11 provides both the location and severity of the damage. In general the number of unknown parameters in that equation is significantly greater than the number of measured frequencies and mode shapes; the problem is therefore under-determined and has an infinite number of probable solutions. A unique solution can be obtained through the minimization of an objective function subjected to appropriate constraints. In the present study, the objective function to be minimized is chosen to be the quadratic norm of the stiffness changes given by Eq. 12, subjected to the constraint specified in Eq. 13.



The problem now becomes a nonlinear optimization problem. In the present work, it is solved using the Sequential Quadratic Programming (SQP) method. Gill et al. (1981) provide an overview of the SQP methods.

M-FEM program structure
Description: M-FEM has been implemented in MATLAB for its simplicity and the availability of a number of basic functionalities for matrix operations, optimization and graphics. The programming environment in MATLAB is interpretive or scripting based, unlike other computer languages like FORTRAN or C/C++ which needs the program to be compiled before running. MATLAB code or program does not need to be compiled. The programming errors are detected at runtime. For fast prototyping of an algorithm, MATLAB is a useful tool. However, it provides flexibility in programming and debugging at the cost of computing efficiency.

The basic structure of the M-FEM system is given in Fig. 1. The program has a modular structure with following basic functionalities (a) Static Analysis (b) Modal Analysis (c) Vibration Based Damage Identification.

Element library: The element library in M-FEM includes a three dimensional beam element and a facet plate-shell element as shown in Fig. 2a and b, respectively. The plate-shell element is formed using the DKT plate bending element (Batoz et al., 1980) and a triangular plane stress element. The DKT plate element has three degrees of freedom (DOF) per node (Fig. 2e), two rotational DOFs (θx and θy) about the orthogonal horizontal directions and a translational DOFs (w) along the vertical direction. When combined with a triangular plane stress element, the plate bending element results in a facet shell element. The resulting shell element requires six DOFs per node as shown in Fig. 2b. Of these six DOFs, three are contributed by the DKT plate bending element and three must come from the plane stress element. Two versions of the plane stress element have been implemented in M-FEM:

Fig. 1: M-FEM system architecture

Fig. 2: M-FEM element library

The Constant Strain Triangle (CST) with two Degrees Of Freedom (DOF) per node corresponding to the in-plane displacements as shown in Fig. 2c (Zienkiewicz and Taylor, 1989),
The plane stress triangular element with three DOFs per node, two corresponding to the in-plane displacements and the other corresponding to the vertex rotation or drilling DOF as shown in Fig. 2d (Allman, 1984). This plane element is also referred in this study as the Allman’s triangle and linear stress triangle with drilling DOF (LSTD) element.

Since the CST element has only two DOFs (u and v in x and y directions, respectively) per node, the resulting shell element has zero stiffness associated with the remaining DOF corresponding to the rotation about the z axis (θz) or the drilling DOF and in the case of DKT-CST combination, a small stiffness of arbitrary magnitude may need to be added to the drilling DOF to avoid numerical difficulties when some of the adjacent elements are co-planar. Unlike the CST element, the drilling DOF stiffness arises naturally in the LSTD element and no adjustment is necessary for adjacent co-planer elements. Therefore, DKT-LSTD element has been used in the case studies presented here.

A quadrilateral plate-shell element is formed by combining two triangular elements as shown in Fig. 2f. In case of the CST element, a small stiffness corresponding to the drilling DOF needs to be added in order to construct the plate-shell element with six DOFs per node.

Operation of the software: The M-FEM software works partially in batch mode and partially in interactive mode. The data about the structure’s geometry, material properties, load and boundary conditions, sensor locations, the extent of measurement noise and the damage simulation scenario are provided in a data file in text format. The software reads the data and builds the mathematical model of the structure. After that the user needs to interact with the software based on a number of choices about the damage identification algorithms, plotting various results, viewing the structure’s geometry and mode shapes, using incomplete and noisy modes of the damaged structure, convergence criteria for model updating etc. The model once built by M-FEM resides in MATLAB workspace and the can be re-run with various damage conditions or measurement noise without re-building the base model. This saves time significantly when the analysis needs to be re-run with the same base model of the structure. The performance weakness of MATLAB is compensated in this case. M-FEM also utilizes a MATLAB feature to save the workspace for future use, in case a few more cases of simulations need to be done later using the base model of the structure already developed by M-FEM.

Validation: A number of structures involving 3D beam, plate, shell and a combination have been modeled with M-FEM and COSMOS (1995) to compare the static and modal analysis results. COSMOS implements the shell element as a combination of DKT plate bending element and Allman’s plane element. For validation of M-FEM, the global stiffness and mass matrices for a number of simple problems, not shown here to conserve space, have been obtained from both M-FEM and COSMOS for comparison. The coefficients of the global stiffness and mass matrices formed by M-FEM and COSMOS have also been compared and which have been found to have very little or no difference. After the correctness of the stiffness and mass matrices are established, the VBDI algorithms have been implemented on top of it.


The Crowchild Bridge located in Calgary, Alberta, Canada, is a two-lane traffic overpass with three continuous spans. The details of the bridge superstructure are shown in Fig. 3. The superstructure is said to be the first continuous-span steel-free bridge deck in the world (Bakht and Mufti, 1998). It is composed of five longitudinal steel girders (900 mm deep), a polypropylene fiber reinforced concrete slab deck and prefabricated glass fiber reinforced concrete barriers. The five longitudinal girders are spaced at 2 m. Four evenly spaced cross-frames in each span and steel girder diaphragms at the supports hold the main girders in place. The main girders are also connected by evenly spaced steel straps placed across the top of the girders to provide lateral restraint to them. The girders and straps are connected to the deck slab by stainless steel stubs. The deck is 9030 mm wide and does not contain any internal steel reinforcement. The slab thickness is 275 mm along the girders and 185 mm elsewhere. A monitoring program for the bridge has been developed by ISIS Canada. Static and ambient vibrations tests have been conducted on the bridge by a team from the University of British Columbia in 1997 (Ventura et al., 2000) and a team from the University of Alberta in 1998-99 (Cheng and Afhami, 1999).

An analytical model of the Crowchild Bridge is constructed using three dimensional beam elements for the piers, girders, diaphragms and the cross frames including the steel straps and shell elements for the deck and side barriers. The deck elements are connected to the girder elements by rigid beam elements. The piers are assumed to be fixed at their base, while roller and pin supports are assumed to exist at the north and south abutments, respectively. The FE model contains 351 elements, 247 nodes and 1399 active degrees of freedom. The density of steel and concrete is assumed to be 76 and 24 kN m‾3, respectively. The concrete compressive strength is taken as 35 MPa. The modulus of elasticity for concrete is assumed to be 30 GPa for the deck and 27 GPa for the barrier and pier; for steel it is assumed to be 200 GPa. The FE model of Crowchild Bridge is shown in Fig. 4a and the damage area assumed for computer simulation is shown in Fig. 4b.

Initially, the concrete is assumed to be un-cracked. The model that is based on this assumption is updated and correlated with the data obtained from vibration test conducted in 1997 by the University of British Columbia team (Ventura et al., 2000) and is used as the base model for the undamaged structure as of 1997. This model will be referred to as CC-97 model. The first four natural frequencies were reported to be 2.78, 3.41, 3.44 and 3.89 Hz. The model update and correlation were carried out by using the optimization process defined by Eq. 10 and 11. During the process of model updating, the stiffness coefficients of individual elements are modified to fine tune the resulting modal frequencies of the system.

Fig. 3: Details of the Crowchild Bridge (a) west elevation, (b) Cross section

The mode shapes of the undamaged structure derived from model CC-97 are shown in Fig. 5.

Field test conducted by the University of Alberta team in 1998-99 (Cheng and Afhami, 1999) revealed that the frequencies of the bridge reduced slightly due to some cracking of concrete. To consider the effect of cracking the stiffness of the concrete elements is reduced and the model is correlated with the 1999 test results. The stiffness of deck elements is assumed to be 90% of un-cracked stiffness in the positive moment region and 70% in the negative moment region. The stiffness of the barrier is assumed to be 80% in the positive moment region and 70% in the negative moment region. The pier stiffness is reduced by 10%. The correlated base model will be referred to as CC-99. The natural frequencies at this stage were reported to be 2.60, 2.90, 3.63 and 3.85 Hz. The mode shapes are similar to those of CC-97 (Fig. 5).

The vibration based damage identification techniques implemented in M-FEM has been studied by applying the techniques to detect simulated damage in the bridge models. The damage is assumed to be concentrated in the area identified in Fig. 4b.

Fig. 4: (a) The Finite Element model of Crowchild Bridge, (b) Damaged area for computer simulation

Fig. 5: Mode shapes: (a) mode 1, (b) mode 2, (c) mode 3 and (d) mode 4

The following scenario is considered: the longitudinal girders (corresponding to element numbers 174, 177, 180, 183 and 186 in the finite element model) suffer 70% damage. The damage magnitude of 70% is rather high, but it is chosen only to demonstrate the program. The damage simulated herein by reducing the stiffness of the corresponding elements in the model. The mode shape and frequencies obtained from the damaged model by modal analysis. In practice however, such frequencies and mode shapes would be obtained from modal tests.

The experimental frequencies and mode shapes obtained from the modal tests are often affected by random measurement noise. Another factor to be considered is that the finite element model usually contains a large number of degrees of freedom, but in the field only a limited number of sensors could be deployed. In addition, measurement of rotational degrees of freedom is not practical. Thus, only a few translation degrees of freedom are measured and the experimental mode shapes are far from complete. This implies that either the incomplete mode shapes need to be expanded to the full size of the finite element model, or the finite element model is condensed to the measured degrees of freedom.

In the present study, incomplete mode shapes are constructed by selecting those elements of the analytical damaged mode shapes that correspond to a set of degrees of freedom along which measurements would be probably be made during a field test. The incomplete mode shapes are corrupted by adding a small amount of random error to simulate the noise in measurement; they are then expanded using a method proposed by O’Callahan et al. (1989) to obtain a complete set of mode shapes for the damaged structure. In this study, the incomplete mode shapes have been generated assuming that all three translation degrees of freedom at every node have been measured, while in practice, perhaps fewer degrees of freedom would be measured. To simulate the random measurement noise, the eigenvalues and mode shapes of the damaged structure are corrupted by a random noise of 0.5 and 1%.

Figure 6 and 7 show the damage identification results for CC-97 model. It is observed from Fig. 6a that the damage index identifies the damage locations correctly when complete and noise free mode shapes of the damaged structure are used. Figure 6b and c show the severity of damage as expressed by damage factors (β vector in Eq. 11) evaluated using the matrix update methods. Estimates of β obtained from Eq. 11 by taking the pseudo inverse of D provide rough estimates of the severity of damage and are shown in Fig. 6b. It should be noted that in this case β can take negative values.

Fig. 6:
Results for the CC-97 model with complete and noise free modes of the damaged structure-(a) Damage index, (b) Damage factors by pseudo inverse solution of Eq. 11 and (c) Damage factors by optimized solution of Eq. 11

The solution of Eq. 11 obtained from the optimization of the objective function in Eq. 12 with the constraints defined by Eq. 13 give a very good measure of the damage extent (Fig.6c).

Figure 7 shows the damage indices and damage factors CC-97 model with noisy (E1) and incomplete mode shapes (M1). Determination of the location and the extent of damage are affected both by the measurement errors and incomplete modes. The general area of damage is identified correctly, while the severity of damage is not as accurate. The pseudo-inverse solution produces the damage factors similar to that in Fig. 6, while the optimized solution of Eq. 11 produces spurious damage in some elements, which, of course, can be eliminated when the damage indices in Fig. 7a are considered.

Fig. 7:
Results for the CC-97 model with incomplete and noisy modes of the damaged structure (a) Damage index, (b) Damage factors by pseudo inverse solution of Eq. 11 and (c) Damage factors by optimized solution of Eq. 11

From the results of damage identification for CC-99 model with incomplete and noisy modes, not shown here to conserve space, it is observed that the damage index method identifies the damage location quite effectively, but the estimates of the damage extent by matrix update method are not very good. However, when the results of the damage index method and the matrix update method are considered together, it is possible to get important clues about the severity and location of the damage.


A finite element analysis system, M-FEM for vibration based damage identification in structures has been developed in the MATLAB environment. M-FEM utilizes the flexibility of the MATLAB framework and the availability of many useful functions necessary for a FEM system. Although software built using MATLAB performs poorly as compared to that developed in a compiler-based language like FORTRAN or C/C++, the efficient use of MATLAB workspace can be advantageous when it is necessary to re-use the base model built using M-FEM.

Vibration based damage identification is a useful tool for structural health monitoring. In the past its success has been demonstrated by applications in the damage detection of simple structures. However, its application to real life and complex civil engineering structures is still beset with considerable challenge. The first challenge manifests itself in creating a realistic mathematical model of the structure. This is because structural properties, such as the effective moment of inertia and modulus of elasticity, are quite uncertain and variable, especially in structures built from reinforced concrete where concrete cracking and quality of concrete may have a marked influence on the physical properties. Most modeling attempts, including those used in the present study, use arbitrary estimates of the effective moments of inertia, based on engineering judgment. Correlation of a model which must be based on such choices often requires unrealistic adjustments in the element stiffness coefficients, even when the differences between the analytical and experimental frequencies are small. In the present case, such adjustments could be kept within reasonable bounds only after the selection of the physical properties was fine tuned by a process of trial and error.

Other challenges include the practical limitation on the degrees of freedom along which modal measurements can be made, the number of modes that can be measured and the measurement errors that are inevitably present in field tests. As shown here all of these have profound impact on the success of damage identification. In the present study, in order for the methods to succeed, measurements were required along all of the translational degrees of freedom.

Of the many VBDI techniques available only two are used here because of their relative merit in comparison to others. The results obtained form these methods show that in spite of all of the problems, VBDI methods may be successful in identifying at least the location of damage in a structure, such as the Crowchild Bridge. It is generally not possible to detect damages of small magnitudes through VBDI techniques since the changes in modal properties due to small amounts of damage are not appreciable. Also, if the modal strain energy contribution of an element is very small, even a large degree of damage in it could not be detected by the VBDI techniques. Finally, it should be noted that computer simulation studies can not quite predict the effect of environmental factors, such as temperature changes, changes in boundary conditions and nonlinearities caused by the opening and closing of cracks and slip in bolted joints. These are potential areas of further research.


The financial support provided by ISIS Canada Research Network to the second author and Concordia University (through the SRT and FRDP grants) to the first author are gratefully acknowledged.

1:  Allman, D.J., 1984. A compatible triangular element including vertex rotations for plane elasticity analysis. Comput. Struct., 9: 1-8.

2:  Bakht, B. and A. Mufti, 1998. Five steel-free bridge deck slabs in Canada. SEI Struct. Eng. Int., 8: 196-200.
Direct Link  |  

3:  Batoz, J.L., K.J. Bathe and L.W. Ho, 1980. A study of three-node triangular plate bending element. Int. J. Numeric. Methods Eng., 15: 1771-1812.

4:  Bagchi, A., 2005. Development of a finite element system for vibration based damage identification in structures. Proceedings of the 2nd International Conference on Structural Health Monitoring of Intelligent Infrastructure (SHMII'2), Nov. 15-18, Shenzen, China.

5:  Cheng, J.J.R. and S. Afhami, 1999. Field instrumentation and monitoring of crowchild bridge in calgary, alberta. University of Alberta Report for ISIS Canada.

6:  COSMOS/M, 1995. Finite element analysis system user guide. Struct. Res. Anal. Corp., Santa Monica, CA.

7:  Gill, P.E., W. Murray and M.H. Wright, 1981. Practical Optimization. Academic Press, London.

8:  Humar, J.L., A. Bagchi and H. Xu, 2003. Challenges in vibration based structural health monitoring. Proceedings of the 1st International Structural Health Monitoring II-1'2003 Conference, Nov. 13-15, Toyko, Japan.

9:  Humar, J.L., A. Bagchi and H. Xu, 2006. Performance of vibration-based techniques for the identification of structural damage. Struct. Health Monitoring-An Int. J., 5: 215-241.
Direct Link  |  

10:  Kabe, A.M., 1985. Stiffness matrix adjustment using mode data. AIAA J., 23: 1431-1436.

11:  MATLAB, 2006. User Guide. Mathworks Inc., Natick, MA.

12:  O`Callahan, J., P. Avitabile and R. Riemer, 1989. System equivalent reduction expansion process (SEREP). Proceedings of the 7th International Modal Analysis Conference, January 30-February 2, 1989, Las Vegas, Nevada, pp: 29-37.

13:  Stubbs, N., J.T. Kim and C.R. Farrar, 1995. Field verification of a non-destructive damage localization and severity estimation algorithm. Proceedings of the 13th International Modal Analysis Conference, February 13-16, 1995, USA., pp: 210-218.

14:  Ventura, C.E., T. Onur and R.C. Tsai, 2000. Dynamic characteristics of the Crowchild Trail Bridge. Can. J. Civil Eng., 27: 1046-1056.
Direct Link  |  

15:  Zienkiewicz, O.C. and R.L. Taylor, 1989. The Finite Element Method. 4th Edn., Vol. 1-2, McGraw Hill, UK.

©  2021 Science Alert. All Rights Reserved