INTRODUCTION
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.
VIBRATION BASED DAMAGE IDENTIFICATION
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 φ^{T}_{i}k_{j}φ_{i} and the total strain energy of the structure deforming in its mode I is obtained from
where k_{j} 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, S^{d}_{i},
the global stiffness matrix of the damaged structure K^{d} 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, f_{ij} and
f^{d}_{ij}, respectively, are given by:
The ratio of the strain energy fractions as expressed in Eq.
3 is defined as the damage index for the j^{th} 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 massorthonormal 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. Premultiplying both sides
of Eq.7 by φ^{T}_{di} 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 orthonormalized
mode shape vectors, we get
where ne is the number of members and D is an m by ne matrix, β is the
nevector of the unknown changes in element stiffness matrices and δλ
is the mvector 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 underdetermined 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.
MFEM program structure
Description: MFEM 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 MFEM 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 MFEM includes a three dimensional
beam element and a facet plateshell element as shown in Fig.
2a and b, respectively. The plateshell 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 MFEM:

Fig. 1: 
MFEM system architecture 

Fig. 2: 
MFEM element library 
• 
The Constant Strain Triangle (CST) with two Degrees Of Freedom
(DOF) per node corresponding to the inplane displacements as shown in Fig.
2c (Zienkiewicz and Taylor, 1989), 
• 
The plane stress triangular element with three DOFs per node, two corresponding
to the inplane 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 DKTCST 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 coplanar. Unlike the CST
element, the drilling DOF stiffness arises naturally in the LSTD element and
no adjustment is necessary for adjacent coplaner elements. Therefore, DKTLSTD
element has been used in the case studies presented here.
A quadrilateral plateshell 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
plateshell element with six DOFs per node.
Operation of the software: The MFEM 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 MFEM resides in MATLAB workspace and the can be rerun with various damage conditions or measurement noise without rebuilding the base model. This saves time significantly when the analysis needs to be rerun with the same base model of the structure. The performance weakness of MATLAB is compensated in this case. MFEM 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 MFEM.
Validation: A number of structures involving 3D beam, plate, shell and a combination have been modeled with MFEM 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 MFEM, the global stiffness and mass matrices for a number of simple problems, not shown here to conserve space, have been obtained from both MFEM and COSMOS for comparison. The coefficients of the global stiffness and mass matrices formed by MFEM 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.
CASE STUDY: THE CROWCHILD BRIDGE
The Crowchild Bridge located in Calgary, Alberta, Canada, is a twolane 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 continuousspan steelfree 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 crossframes 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 199899 (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 uncracked. 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 CC97 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 CC97 are shown in Fig.
5.
Field test conducted by the University of Alberta team in 199899 (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 uncracked
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 CC99. 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 CC97 (Fig. 5).
The vibration based damage identification techniques implemented in MFEM 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 CC97 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 CC97 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 CC97
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 pseudoinverse 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 CC97 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 CC99 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.
DISCUSSION
A finite element analysis system, MFEM for vibration based damage identification in structures has been developed in the MATLAB environment. MFEM 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 compilerbased language like FORTRAN or C/C++, the efficient use of MATLAB workspace can be advantageous when it is necessary to reuse the base model built using MFEM.
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.
ACKNOWLEDGMENTS
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.