# Damage Detection of Bridges Using Vibration Data by Adjoint Variable Method

Authors: Akbar Mirzaee; Mohsenali Shayanfar; Reza Abbasnia
Journal: Shock and Vibration
 (2014)

---

## Abstract

This research entails a theoretical and numerical study on a new damage detection method for bridges, using response sensitivity in time domain. This method, referred to as “adjoint variable method,” is a finite element model updating sensitivity based method. Governing equation of the bridge-vehicle system is established based on finite element formulation. In the inverse analysis, the new approach is presented to identify elemental flexural rigidity of the structure from acceleration responses of several measurement points. The computational cost of sensitivity matrix is the main concern associated with damage detection by these methods. The main advantage of the proposed method is the inclusion of an analytical method to augment the accuracy and speed of the solution. The reliable performance of the method to precisely identify the location and intensity of all types of predetermined single, multiple, and random damages over the whole domain of moving vehicle speed is shown. A comparison study is also carried out to demonstrate the relative effectiveness and upgraded performance of the proposed method in comparison to the similar ordinary sensitivity analysis methods. Moreover, various sources of errors including the effects of noise and primary errors on the numerical stability of the proposed method are discussed.


---

## Body

## 1. Introduction

The main objective of developing the structural health monitoring (SHM) system for structures is to enhance structural safety. However, in bridges, SHM serves other economic benefits such as increased mission reliability, extended life of life-limited components, reduced tests, reduction in “down time,” increased equipment reliability, customization of maintenance actions, and greater awareness of operating personnel, resulting in fewer accidents. SHM also promises to help reduce maintenance costs [1].SHM algorithms are identified as static system identification (SI) and dynamic SI, according to the types of structural response used. Dynamics-based SI techniques assess the state of health of a structural component on the basis of the detection and analysis of its dynamic response. Such techniques can be classified on the basis of the type of response being considered for the investigations, on the frequency or time domain of interrogation and on the modality used to excite the component [2].The frequency domain SI and time-domain SI are more practical than static SI, as the static response, for instance, displacements of a structure are very difficult to measure in most cases.The developments in the field of SI using vibration data of civil engineering structures have been recently reviewed by several authors. Some recent studies are briefly described in the following.Doebling et al. [3, 4] have presented comprehensive review of literature mainly focusing on frequency-domain methods for damage detection in linear structures and declared that sufficient evidence exists to promote the use of measured vibration data for the detection of damage in structures, using both forced-response testing and long-term monitoring of ambient signals and there is a significant need in this field for research on the integration of theoretical algorithms with application-specific knowledge bases and practical experimental constraints. Another discussion on methods of damage detection and location using natural frequency changes has been presented by Salawu [5] and his study showed that damage detection using vibration frequencies is not very reliable.Alampalli and Fu [6] and Alampalli et al. [7] conducted laboratory and field studies on bridge structures to investigate the feasibility of measuring bridge vibration for inspection and evaluation. These studies focused on sensitivity of measured modal parameters to damage. Cross diagnosis using multiple signatures involving natural frequencies, mode shapes, modal assurance criteria, and coordinate modal assurance criteria was shown to be necessary to detect the damages. Casas and Aparicio studied concrete bridge structures and investigated dynamic response as an inspection tool to assess bearing conditions and girder cracking [8]. Their study showed the need to investigate more than one natural frequency and also to determine mode shapes in order that the damage could be successfully detected and located.The frequency-domain SI algorithms have been more widely developed and applied as the amount of measured data is reduced dramatically after the transform; thus, they can be handled easily. Unfortunately, the effects of local damages on the natural frequencies and mode shapes of higher modes are greater than lower ones, but they are usually difficult to measure from experiments. In addition, structural damping properties cannot be identified in frequency domain SI.The time-domain SI may be an attractive one to overcome the drawbacks of the frequency-domain SI. For time-domain SI, the forced vibration responses of the structure are needed in the identification. However, in some cases, it is either impractical or impossible to use artificial inputs to excite the civil engineering structures, so natural excitation must be measured along with the structural responses to assess the dynamic characteristics [9, 10]. In recent years, some researchers have investigated both the problem of load identification (moving load and impact load) and modal parameters identification under operational conditions [11, 12]. In addition, identification of the structural parameters applying a moving load has been considered in many papers. Law et al. [13] presented a novel moving force and prestress identification method based on the finite element and the wavelet-based methods for a bridge-vehicle system. Jiang et al. [14] identified the parameter of a vehicle moving on multispan continuous bridges. Zhu and Law [15] presented a method for damage detection of a simply supported concrete bridge structure in time domain using the interaction forces from the moving vehicles as excitation. Majumder and Manohar [16] proposed a time-domain approach for damage detection in beam structures using vibration data induced by a vehicle moving on a bridge deck.Sensitivity based methods allow a wide choice of physically meaningful parameters and this advantage has led to their widespread use in damage detection. The most important difficulty in sensitivity based SI methods is calculation of sensitivity matrix. Calculation of this massive matrix is repeated in each iteration and it is so time-consuming and has a significant effect on the efficiency of method. Despite the high importance of calculation method of sensitivity matrix and the optimization of its performance in SI procedure, there is no literature on this regard. In this paper, computational methods for sensitivity matrix are discussed and a novel sensitivity based damage detection method in time-domain referred to as “Adjoint variable method,” is developed. Computational algorithm of proposed method is presented and its performance is compared with the conventional methods and it is shown that the numerical cost is considerably reduced by using the concept of adjoint variable.The outline of the work is as follows: inverse problems along with model updating are briefly introduced in Section2. The basic theory of sensitivity analysis is addressed in Section 3 and the proposed algorithm will be presented in Section 4. Numerical simulations along with comparison studies are presented in Section 5 with studies on the effect of different factors which may affect the accuracy of the proposed analysis in practice. Conclusion will be drawn in the last section.


## 2. Finite Element Model Updating and Inverse Problem

Since many algorithms of damage detection are based on the difference between modified model before occurrence of damage and after that, problems such as parameter identification and damage detection are closely related to model updating. Discrepancy between two models is used for detection and quantification of damage.A key step in model-based damage identification is the updating of the finite element model of the structure in such a way that the measured responses can be reproduced by the FE model. A general flowchart of this operation is given in Figure1. The identification procedure presented in this paper is a sensitivity based model updating routine. Sensitivity coefficients are the derivatives of the system responses with respect to the physical parameters or input excitation force and are needed in the cost function of the flowchart of Figure 1.Figure 1
General flowchart of a FEM updating.


### 2.1. Finite Element Modeling of Bridge Vibration under Moving Loads

For a general finite element model of a linear elastic time-invariant structure, the equation of motion is given by(1)
[
M
]
{
z
,
t
t
}
+
[
C
]
{
z
,
t
}
+
[
K
]
{
z
}
=
[
B
]
{
F
}
,
where [
M
] and [
K
] are mass and stiffness matrices and [
C
] is damping matrix. Z
,
t
t, Z
,
t, and Z are the respective acceleration, velocity, and displacement vectors for the whole structure and {
F
} is a vector of applied forces with matrix [
B
] mapping these forces to the associated Dof’s of the structure. A proportional damping is assumed to show the effect of damping ratio on the dynamic magnification factor. Rayleigh damping, in which the damping matrix is proportional to the combination of the mass and stiffness matrices, is used. Consider
(2)
[
C
]
=
a
0
[
M
]
+
a
1
[
K
]
,
where a
0 and a
1 are constants to be determined from two modal damping ratios. If a more accurate estimation of the actual damping is required, a more general form of Rayleigh damping, the Caughey damping model, can be adopted.The dynamic responses of the structures can be obtained by direct numerical integration using Newmark method.


### 2.2. Objective Functions

The approach minimizes the difference between response quantities (usually acceleration response) of the measured data and model predictions. This problem may be expressed as the minimization ofJ, where
(3)
J
(
θ
)
=
∥
z
m
-
z
(
α
)
∥
2
=
ϵ
T
ϵ
,
ϵ
=
z
m
-
z
(
α
)
.
Here, z
m and z
(
α
) are the measured and computed response vectors, α is a vector of all unknown parameters, and ϵ is the response residual vector.


### 2.3. Penalty Function Methods

When the parameters of a model are unknown, they must be estimated using measured data. The measured response is a nonlinear function of the parameters. So, minimizing the error between the measured and predicted response will produce a nonlinear optimization problem.Penalty function method is generally used for modal sensitivity with a truncated Taylor series expansion in terms of the unknown parameters. In this paper, the truncated series of the dynamic responses in terms of the system parameterα are used to derive the sensitivity based formulation. The identification problem can be expressed as follows to find the vector {
α
} such that the calculated response best matches the measured response; that is,
(4)
[
Q
]
{
R
}
=
{
R
^
}
,
where the selection matrix [
Q
] is a matrix with elements of zeros or ones, matching the Dof’s corresponding to the measured response components. Vector {
R
} can be obtained from (4) for a given set of {
α
}.Let(5)
{
δ
z
}
=
{
R
^
}
-
[
Q
]
{
R
}
=
{
R
^
}
-
{
R
cal
}
,
where {
δ
z
} is the error vector in the measured output. In the penalty function method, we have
(6)
{
δ
z
}
=
[
S
]
{
δ
α
}
,
where {
δ
α
} is the perturbation in the parameters and [
S
] is the two-dimensional sensitivity matrix which is one of the matrices at time t in the three-dimensional sensitivity matrix shown in Figure 2 [17]. For a finite element model with N elements each with M system parameters, the number of unknown parameters is N
×
M, and N
×
M equations are needed to solve the parameters. Matrix [
S
] is on the parameter-t plane in Figure 2, and we can select any row of the three-dimensional sensitivity matrix, say, the ith row corresponding to the ith measurement for the purpose. When writing in full, (5) can be written as
(7)
{
δ
z
}
=
{
R
^
(
t
1
)
R
^
(
t
2
)
⋮
R
^
(
t
l
)
}
-
{
R
cal
(
t
1
)
R
cal
(
t
2
)
⋮
R
cal
(
t
l
)
}
with l
≥
N
×
M to make sure that the set of equation is overdetermined. Equation (6) can be solved by simple least-squares method as follows:
(8)
δ
α
=
[
S
T
S
]
-
1
S
T
δ
z
,
(9)
α
j
+
1
=
α
j
+
[
S
j
T
S
j
]
-
1
S
j
T
(
R
^
-
R
cal
)
.
The subscript j indicates the iteration number at which the sensitivity matrix is computed.Figure 2
Three-dimensional sensitivity matrix.One of the important difficulties in parameter estimation is ill-conditioning. In the worst case, this can mean that there is no unique solution to the estimation problem and many sets of parameters are able to fit the data. Many optimization procedures result in the solution of linear equations for the unknown parameters. The use of the singular value decomposition (SVD) [18] for these linear equations enables ill-conditioning to be identified and quantified. The options are then to increase the available data, which is often difficult and costly, or to provide extra conditions on the parameters. These can take the form of smoothness conditions (e.g., the truncated SVD), minimum norm parameter values (Tikhonov regularization), or minimum changes from the initial estimates of the parameters [19, 20].From experiences gained in model updating with simulated structures, Li and Law [21] found that Tikhonov regularization can give the optimal solution when there is no noise or very small noise in the measurement.


### 2.4. Tikhonov Regularization

Like many other inverse problems, (6) is an ill-conditioned problem. In order to provide bounds to the solution, the damped least-squares method (DLS) is used and singular-value decomposition is used in the pseudoinverse calculation. Equation (8) can be written in the following form:
(10)
δ
α
=
(
S
T
S
+
λ
R
I
)
-
1
S
T
δ
z
,
where λ is the nonnegative damping coefficient governing the participation of least-squares error in the solution. The solution of (10) is equivalent to minimizing the function:
(11)
J
(
{
δ
α
}
,
λ
R
)
=
∥
S
δ
α
-
δ
z
∥
2
+
λ
R
∥
δ
α
∥
2
with the second term in (11) that provides bounds to the solution. When the parameter λ
R approaches zero, the estimated vector {
δ
α
} approaches the solution obtained from the simple least-squares method. L-curve method is used in this paper to obtain the optimal regularization parameter λ
R.


### 2.5. Element Damage Index

In the inverse problem of damage identification, it is assumed that the stiffness matrix of the whole element decreases uniformly with damage, and the flexural rigidity,E
I
i of the ith finite element of the beam, becomes β
i
E
I
i, when there is damage [22]. The fractional change in stiffness of an element can be expressed as
(12)
Δ
K
b
i
=
(
K
b
i
-
K
~
b
i
)
=
(
1
-
β
i
)
K
b
i
,
where K
b
i and K
~
b
i are the ith element stiffness matrices of the undamaged and damaged beam, respectively. Δ
K
b
i is the stiffness reduction of the element. A positive value of β
i
∈
[
0,1
] will indicate a loss in the element stiffness. The ith element is undamaged when β
i
=
1 and the stiffness of the ith element is completely lost when β
i
=
0.The stiffness matrix of the damaged structure is the assemblage of the entire element stiffness matrixK
~
b
i:
(13)
K
b
=
∑
i
=
1
N
A
i
T
K
~
b
i
A
i
=
∑
i
=
1
N
β
i
A
i
T
K
b
i
A
i
,
where A
i is the extended matrix of element nodal displacement that facilitates assembling of global stiffness matrix from the constituent element stiffness matrix.


## 2.1. Finite Element Modeling of Bridge Vibration under Moving Loads

For a general finite element model of a linear elastic time-invariant structure, the equation of motion is given by(1)
[
M
]
{
z
,
t
t
}
+
[
C
]
{
z
,
t
}
+
[
K
]
{
z
}
=
[
B
]
{
F
}
,
where [
M
] and [
K
] are mass and stiffness matrices and [
C
] is damping matrix. Z
,
t
t, Z
,
t, and Z are the respective acceleration, velocity, and displacement vectors for the whole structure and {
F
} is a vector of applied forces with matrix [
B
] mapping these forces to the associated Dof’s of the structure. A proportional damping is assumed to show the effect of damping ratio on the dynamic magnification factor. Rayleigh damping, in which the damping matrix is proportional to the combination of the mass and stiffness matrices, is used. Consider
(2)
[
C
]
=
a
0
[
M
]
+
a
1
[
K
]
,
where a
0 and a
1 are constants to be determined from two modal damping ratios. If a more accurate estimation of the actual damping is required, a more general form of Rayleigh damping, the Caughey damping model, can be adopted.The dynamic responses of the structures can be obtained by direct numerical integration using Newmark method.


## 2.2. Objective Functions

The approach minimizes the difference between response quantities (usually acceleration response) of the measured data and model predictions. This problem may be expressed as the minimization ofJ, where
(3)
J
(
θ
)
=
∥
z
m
-
z
(
α
)
∥
2
=
ϵ
T
ϵ
,
ϵ
=
z
m
-
z
(
α
)
.
Here, z
m and z
(
α
) are the measured and computed response vectors, α is a vector of all unknown parameters, and ϵ is the response residual vector.


## 2.3. Penalty Function Methods

When the parameters of a model are unknown, they must be estimated using measured data. The measured response is a nonlinear function of the parameters. So, minimizing the error between the measured and predicted response will produce a nonlinear optimization problem.Penalty function method is generally used for modal sensitivity with a truncated Taylor series expansion in terms of the unknown parameters. In this paper, the truncated series of the dynamic responses in terms of the system parameterα are used to derive the sensitivity based formulation. The identification problem can be expressed as follows to find the vector {
α
} such that the calculated response best matches the measured response; that is,
(4)
[
Q
]
{
R
}
=
{
R
^
}
,
where the selection matrix [
Q
] is a matrix with elements of zeros or ones, matching the Dof’s corresponding to the measured response components. Vector {
R
} can be obtained from (4) for a given set of {
α
}.Let(5)
{
δ
z
}
=
{
R
^
}
-
[
Q
]
{
R
}
=
{
R
^
}
-
{
R
cal
}
,
where {
δ
z
} is the error vector in the measured output. In the penalty function method, we have
(6)
{
δ
z
}
=
[
S
]
{
δ
α
}
,
where {
δ
α
} is the perturbation in the parameters and [
S
] is the two-dimensional sensitivity matrix which is one of the matrices at time t in the three-dimensional sensitivity matrix shown in Figure 2 [17]. For a finite element model with N elements each with M system parameters, the number of unknown parameters is N
×
M, and N
×
M equations are needed to solve the parameters. Matrix [
S
] is on the parameter-t plane in Figure 2, and we can select any row of the three-dimensional sensitivity matrix, say, the ith row corresponding to the ith measurement for the purpose. When writing in full, (5) can be written as
(7)
{
δ
z
}
=
{
R
^
(
t
1
)
R
^
(
t
2
)
⋮
R
^
(
t
l
)
}
-
{
R
cal
(
t
1
)
R
cal
(
t
2
)
⋮
R
cal
(
t
l
)
}
with l
≥
N
×
M to make sure that the set of equation is overdetermined. Equation (6) can be solved by simple least-squares method as follows:
(8)
δ
α
=
[
S
T
S
]
-
1
S
T
δ
z
,
(9)
α
j
+
1
=
α
j
+
[
S
j
T
S
j
]
-
1
S
j
T
(
R
^
-
R
cal
)
.
The subscript j indicates the iteration number at which the sensitivity matrix is computed.Figure 2
Three-dimensional sensitivity matrix.One of the important difficulties in parameter estimation is ill-conditioning. In the worst case, this can mean that there is no unique solution to the estimation problem and many sets of parameters are able to fit the data. Many optimization procedures result in the solution of linear equations for the unknown parameters. The use of the singular value decomposition (SVD) [18] for these linear equations enables ill-conditioning to be identified and quantified. The options are then to increase the available data, which is often difficult and costly, or to provide extra conditions on the parameters. These can take the form of smoothness conditions (e.g., the truncated SVD), minimum norm parameter values (Tikhonov regularization), or minimum changes from the initial estimates of the parameters [19, 20].From experiences gained in model updating with simulated structures, Li and Law [21] found that Tikhonov regularization can give the optimal solution when there is no noise or very small noise in the measurement.


## 2.4. Tikhonov Regularization

Like many other inverse problems, (6) is an ill-conditioned problem. In order to provide bounds to the solution, the damped least-squares method (DLS) is used and singular-value decomposition is used in the pseudoinverse calculation. Equation (8) can be written in the following form:
(10)
δ
α
=
(
S
T
S
+
λ
R
I
)
-
1
S
T
δ
z
,
where λ is the nonnegative damping coefficient governing the participation of least-squares error in the solution. The solution of (10) is equivalent to minimizing the function:
(11)
J
(
{
δ
α
}
,
λ
R
)
=
∥
S
δ
α
-
δ
z
∥
2
+
λ
R
∥
δ
α
∥
2
with the second term in (11) that provides bounds to the solution. When the parameter λ
R approaches zero, the estimated vector {
δ
α
} approaches the solution obtained from the simple least-squares method. L-curve method is used in this paper to obtain the optimal regularization parameter λ
R.


## 2.5. Element Damage Index

In the inverse problem of damage identification, it is assumed that the stiffness matrix of the whole element decreases uniformly with damage, and the flexural rigidity,E
I
i of the ith finite element of the beam, becomes β
i
E
I
i, when there is damage [22]. The fractional change in stiffness of an element can be expressed as
(12)
Δ
K
b
i
=
(
K
b
i
-
K
~
b
i
)
=
(
1
-
β
i
)
K
b
i
,
where K
b
i and K
~
b
i are the ith element stiffness matrices of the undamaged and damaged beam, respectively. Δ
K
b
i is the stiffness reduction of the element. A positive value of β
i
∈
[
0,1
] will indicate a loss in the element stiffness. The ith element is undamaged when β
i
=
1 and the stiffness of the ith element is completely lost when β
i
=
0.The stiffness matrix of the damaged structure is the assemblage of the entire element stiffness matrixK
~
b
i:
(13)
K
b
=
∑
i
=
1
N
A
i
T
K
~
b
i
A
i
=
∑
i
=
1
N
β
i
A
i
T
K
b
i
A
i
,
where A
i is the extended matrix of element nodal displacement that facilitates assembling of global stiffness matrix from the constituent element stiffness matrix.


## 3. Sensitivity Analysis of Transient Dynamic Response

The objective of sensitivity analysis is to quantify the effects of parameter variations on calculated results. Terms such as influence, importance, ranking by importance, and dominance are all related to the sensitivity analysis.


### 3.1. Methods of Structural Sensitivity Analysis

When the parameter variations are small, the traditional way to assess their effects on calculated responses is the employment of perturbation theory, either directly or indirectly, via variational principles. The basic aim of perturbation theory is to predict the effects of small parameter variations without actually calculating the perturbed configuration but rather by using solely unperturbed quantities.Various methods employed in sensitivity analysis are listed in Figure3. Three approaches are used to obtain the sensitivity matrix: the approximation, discrete, and continuum approaches.Figure 3
Different approaches to sensitivity analysis.


### 3.2. Approximation Approach

In the approximation approach, sensitivity matrix is obtained by either the forward finite difference or by the central finite difference method.If the design is perturbed tou
+
Δ
u, where Δ
u represents a small change in the design, then the sensitivity of ψ
(
u
) can be approximated as
(14)
d
ψ
d
u
≈
ψ
(
u
+
Δ
u
)
-
ψ
(
u
)
Δ
u
.
Equation (14) is called the forward difference method since the design is perturbed in the direction of +
Δ
u. If -
Δ
u is substituted in (14) for Δ
u, then the equation is defined as the backward difference method. Additionally, if the design is perturbed in both directions, such that the design sensitivity is approximated by
(15)
d
ψ
d
u
≈
ψ
(
u
+
Δ
u
)
-
ψ
(
u
-
Δ
u
)
2
Δ
u
,
then the equation is defined as the central difference method.


### 3.3. Discrete Approach

In the discrete method, sensitivity matrix is obtained by design derivatives of the discrete governing equation. For this process, it is necessary to take the derivative of the stiffness matrix. If this derivative is obtained analytically using the explicit expression of the stiffness matrix with respect to the variable, it is an analytical method, since the analytical expressions of stiffness matrix are used. However, if the derivative is obtained using a finite difference method, the method is called a semianalytical method. The design represents a structural parameter that can affect the results of the analysis.The design sensitivity information of a general performance measure can be computed either with the direct differentiation method or with the adjoint variable method.


#### 3.3.1. Direct Differentiation Method

The direct differentiation method (DDM) is a general, accurate, and efficient method to compute finite element response sensitivities to the model parameters. This method directly solves for the design dependency of a state variable and then computes performance sensitivity using the chain rule of differentiation. This method clearly shows the implicit dependence on the design, and a very simple sensitivity expression can be obtained.Consider a structure in which the generalized stiffness and mass matrices have been reduced by accounting for boundary conditions. Let the damping force be represented in the form ofC
(
b
)
z
,
t where z
,
t
=
d
z
/
d
t denotes the velocity vector. Under these conditions, Lagrange’s equation of motion becomes the second-order differential equation, as [23]
(16)
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
=
F
(
t
,
b
)
with the initial conditions
(17)
z
(
0
)
=
z
0
,
z
,
t
(
0
)
=
z
,
t
0
.
If design parameters are just related to stiffness matrix, we have
(18)
[
M
]
{
∂
z
,
t
t
∂
b
i
}
+
[
C
]
{
∂
z
,
t
∂
b
i
}
+
[
K
]
{
∂
z
∂
b
i
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
in which {
∂
z
/
∂
b
i
}, {
∂
z
,
t
/
∂
b
i
}, and {
∂
z
,
t
t
/
∂
b
i
} are sensitivity vectors of displacement, velocity, and acceleration with respect to design parameter b
i, respectively. Assume that(19a)
Y
,
t
t
=
∂
z
,
t
t
∂
b
i
,
(19b)
Y
,
t
=
∂
z
,
t
∂
b
i
,
(19c)
Y
=
∂
z
∂
b
i
.So, by replacing (19a), (19b), and (19c) to (18), we have
(20)
[
M
]
{
Y
,
t
t
}
+
[
C
]
{
Y
,
t
}
+
[
K
]
{
Y
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
.
The right side of (20) can be considered as an equivalent force, so (20) is similar to (16) and sensitivity vectors can be obtained by Newmark method.


#### 3.3.2. Adjoint Variable Method

Sensitivity analysis can be performed very efficiently by using deterministic methods based on adjoint functions. The use of adjoint functions for analyzing the effects of small perturbations in a linear system was introduced by Wigner [24].This method constructs an adjoint problem that solves the adjoint variable, which contains all implicit dependent terms.For the dynamic response of structure, the following form of a general performance measure will be considered:(21)
ψ
=
g
(
z
(
T
)
,
b
)
+
∫
0
T
G
(
z
,
b
)
d
t
,
where the final time T is determined by a condition in the form:
(22)
Ω
(
z
(
T
)
,
z
,
t
(
T
)
,
b
)
=
0
.
It is presumed that (22) uniquely determines T, at least locally. This requires that the time derivative of Ω is nonzero at T as follows:
(23)
Ω
,
t
=
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
z
,
t
t
(
T
)
≠
0
.
When final time T is prescribed before the response analysis, the relation in (22) needs not be considered.To obtain the design sensitivity ofΨ, define a design variation in the form:
(24)
b
τ
=
b
+
τ
δ
b
.
Design b is perturbed in the direction of δ
b with the parameter τ. Substituting b
τ into (21), the derivative of (21) can be evaluated with respect to τ at τ
=
0. Leibnitz’s rule of differentiation of an integral may be used to obtain the following expression:
(25)
ψ
′
=
∂
g
∂
b
δ
b
+
∂
g
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
G
(
z
(
T
)
,
b
)
T
′
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
,
where
(26)
z
′
=
z
′
(
b
,
δ
b
)
≡
d
d
τ
z
(
t
,
b
+
τ
δ
b
)
|
τ
=
0
=
d
d
b
[
z
(
t
,
b
)
]
δ
b
,
T
′
=
T
′
(
b
,
δ
b
)
≡
d
d
τ
T
(
b
+
τ
δ
b
)
|
τ
=
0
=
d
T
d
b
δ
b
.
Note that since the expression in (21), that determines T, depends on the design, T will also depend on the design. Thus, terms arise in (25), that involve the derivative of T with respect to the design. In order to eliminate these terms, differentiate (22) with respect to τ and evaluate it at τ
=
0 in order to obtain
(27)
∂
Ω
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
∂
Ω
∂
z
,
t
[
z
,
t
′
(
T
)
+
z
,
t
t
(
T
)
T
′
]
+
∂
Ω
∂
b
δ
b
=
0
.
This equation may also be written as
(28)
Ω
,
t
T
′
=
[
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
t
(
T
)
]
T
′
=
-
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Since it is presumed by (23) that Ω
,
t
≠
0, then
(29)
T
′
=
-
1
Ω
,
t
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Substituting the result of (29) into (25), the following is obtained:
(30)
ψ
′
=
[
∂
g
∂
z
-
(
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
)
1
Ω
,
t
∂
Ω
∂
z
]
z
′
(
T
)
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
+
∂
g
∂
b
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
.
Note that ψ
′ depends on z
′ and z
,
t
′ at T, as well as on z
′ within the integration.In order to writeΨ
′ in (29) explicitly in terms of a design variation, the adjoint variable technique can be used. In the case of a dynamic system, all terms in (16) can be multiplied by λ
T
(
t
) and integrated over the interval [
0
,
T
] to obtain the following identity in λ:
(31)
∫
0
T
λ
T
[
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
-
F
(
t
,
b
)
]
d
t
=
0
.
Since this equation must hold for arbitrary λ, which is now taken to be independent of the design, substitute b
τ into (31) and differentiate it with respect to τ in order to obtain the following relationship:
(32)
∫
0
T
[
λ
T
M
(
b
)
z
,
t
t
′
+
λ
T
C
(
b
)
z
,
t
′
+
λ
T
K
(
b
)
z
′
-
∂
R
∂
b
δ
b
]
d
t
=
0
,
where
(33)
R
=
λ
~
T
F
(
t
,
b
)
-
λ
~
T
M
(
b
)
z
~
,
t
t
-
λ
~
T
C
(
b
)
z
~
,
t
-
λ
~
T
K
(
b
)
z
~
with the superposed tilde (
~
) denoting variables that are held constant during the differentiation with respect to the design in (32).Since (32), contains the time derivatives of z
′, integrate the first two integrands by parts in order to move the time derivatives to λ, as follows:
(34)
λ
T
M
(
b
)
z
,
t
′
(
T
)
-
λ
,
t
T
(
T
)
M
(
b
)
z
′
(
T
)
+
λ
T
C
(
b
)
z
′
(
T
)
+
∫
0
T
{
[
λ
,
t
t
T
M
(
b
)
-
λ
,
t
T
C
(
b
)
+
λ
T
K
(
b
)
]
z
′
-
∂
R
∂
b
δ
b
}
d
t
=
0
.
The adjoint variable method expresses the unknown terms in (30), in terms of the adjoint variable (
λ
). Since (34) must hold for arbitrary functions λ
(
t
), λ may be chosen so that the coefficients of terms involving z
′
(
T
), z
,
t
′
(
T
), and z
′ in (30) and (34) are equal. If such a function λ
(
t
) can be found, then the unwanted terms in (30), involving z
′
(
T
), z
,
t
′
(
T
) and z
′, can be replaced by terms that explicitly depend on δ
b in (34), and to be more specific, choose a λ
(
t
) that satisfies the following:
(35)
M
(
b
)
λ
(
T
)
=
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
t
,
(36)
M
(
b
)
λ
,
t
(
T
)
=
C
T
(
b
)
λ
(
T
)
-
∂
g
T
∂
z
+
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
(37)
M
(
b
)
λ
,
t
t
-
C
T
(
b
)
λ
,
t
+
K
(
b
)
λ
=
∂
G
T
∂
z
,
0
≤
t
≤
T
.
Note that once the dynamic equations of (16) and (17) is solved and (22) is used to determine T, then z
(
T
), z
,
t
(
T
), ∂
Ω
/
∂
z, ∂
Ω
/
∂
z
,
t, and Ω
,
t may be evaluated. Equation (23) can then be solved for λ
(
T
) since the mass matrix M
(
b
) is nonsingular. Having determined λ
(
T
), all terms on the right of (36) can be evaluated, and the equation can be solved for λ
,
t
(
T
). Thus, a set of terminal conditions on λ has been determined. Since M
(
b
) is nonsingular, (37) may then be integrated from T to 0, yielding the unique solution λ
(
t
). Taken as a whole, (35), through (37), may be thought of as a terminal value problem.Since the terms involving a variation in the state variable in (30) and (34) are identical, substitute (34) into (30) to obtain
(38)
ψ
′
=
∂
g
∂
b
δ
b
+
∫
0
T
[
∂
G
∂
b
+
∂
R
∂
b
]
d
t
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
≡
∂
ψ
∂
b
δ
b
.
Every term in this equation can now be calculated. The terms ∂
g
/
∂
b, ∂
G
/
∂
b, and ∂
Ω
/
∂
b represent explicit partial derivatives with respect to the design. The term ∂
R
/
∂
b, however, must be evaluated from (33), thus requiring λ
(
t
). Note also that since design variation δ
b does not depend on time, it is taken outside the integral in (38).Since (38) must hold, for all δ
b, the design derivative vector of ψ is
(39)
d
ψ
d
b
=
∂
g
∂
b
(
z
(
T
)
,
b
)
+
∫
0
T
[
∂
G
∂
b
(
z
,
b
)
+
∂
R
∂
b
(
λ
(
t
)
,
z
(
t
)
,
z
,
t
(
t
)
,
z
,
t
t
(
t
)
,
b
)
]
d
t
-
1
Ω
,
t
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
∂
Ω
∂
b
.


### 3.4. Continuum Approach

In the continuum approach, the design derivative of the variational equation is taken before it is discretized. If the structural problem and sensitivity equations are solved as a continuum problem, then it is called the continuum-continuum method. The continuum sensitivity equation is solved by discretization in the same way that structural problems are solved. Since differentiation is taken at the continuum domain and is then followed by discretization, this method is called the continuum-discrete method.


### 3.5. Sensitivity Method Selection

The advantage of the finite difference method is obvious. If structural analysis can be performed and the performance measure can be obtained as a result of structural analysis, then the expressions in (14) and (15) are virtually independent of the problem types considered.Major disadvantage of the finite difference method is the accuracy of its sensitivity results. Depending on perturbation size, sensitivity results are quite different. For a mildly nonlinear performance measure, relatively large perturbation provides a reasonable estimation of sensitivity results. However, for highly nonlinear performances, a large perturbation yields completely inaccurate results. Thus, the determination of perturbation size greatly affects the sensitivity result. And even though it may be necessary to choose a very small perturbation, numerical noise becomes dominant for a too-small perturbation size. That is, with a too-small perturbation, no reliable difference can be found in the analysis results.The continuum-continuum approach is so limited and is not applicable in complex engineering structures because very simple, classical problems can be solved analytically.The discrete and continuum-discrete methods are equivalent under the conditions given below, using a beam as the structural component. It has also been argued that the discrete and continuum-discrete methods are equivalent under the conditions given below [23].First, the same discretization (shape function) used in the FEA method must be used for continuum design sensitivity analysis. Second, an exact integration (instead of a numerical integration) must be used in the generation of the stiffness matrix and in the evaluation of continuum-based design sensitivity expressions. Third, the exact solution (and not a numerical solution) of the finite element matrix equation and the adjoint equation should be used to compare these two methods. Fourth, the movement of discrete grid points must be consistent with the design parameterization method used in the continuum method.In this paper, two different analytical discrete methods, including direct differential method (DDM) and adjoint variable method (ADM) are presented and efficiency of proposed method is investigated when compared with DDM method.


## 3.1. Methods of Structural Sensitivity Analysis

When the parameter variations are small, the traditional way to assess their effects on calculated responses is the employment of perturbation theory, either directly or indirectly, via variational principles. The basic aim of perturbation theory is to predict the effects of small parameter variations without actually calculating the perturbed configuration but rather by using solely unperturbed quantities.Various methods employed in sensitivity analysis are listed in Figure3. Three approaches are used to obtain the sensitivity matrix: the approximation, discrete, and continuum approaches.Figure 3
Different approaches to sensitivity analysis.


## 3.2. Approximation Approach

In the approximation approach, sensitivity matrix is obtained by either the forward finite difference or by the central finite difference method.If the design is perturbed tou
+
Δ
u, where Δ
u represents a small change in the design, then the sensitivity of ψ
(
u
) can be approximated as
(14)
d
ψ
d
u
≈
ψ
(
u
+
Δ
u
)
-
ψ
(
u
)
Δ
u
.
Equation (14) is called the forward difference method since the design is perturbed in the direction of +
Δ
u. If -
Δ
u is substituted in (14) for Δ
u, then the equation is defined as the backward difference method. Additionally, if the design is perturbed in both directions, such that the design sensitivity is approximated by
(15)
d
ψ
d
u
≈
ψ
(
u
+
Δ
u
)
-
ψ
(
u
-
Δ
u
)
2
Δ
u
,
then the equation is defined as the central difference method.


## 3.3. Discrete Approach

In the discrete method, sensitivity matrix is obtained by design derivatives of the discrete governing equation. For this process, it is necessary to take the derivative of the stiffness matrix. If this derivative is obtained analytically using the explicit expression of the stiffness matrix with respect to the variable, it is an analytical method, since the analytical expressions of stiffness matrix are used. However, if the derivative is obtained using a finite difference method, the method is called a semianalytical method. The design represents a structural parameter that can affect the results of the analysis.The design sensitivity information of a general performance measure can be computed either with the direct differentiation method or with the adjoint variable method.


### 3.3.1. Direct Differentiation Method

The direct differentiation method (DDM) is a general, accurate, and efficient method to compute finite element response sensitivities to the model parameters. This method directly solves for the design dependency of a state variable and then computes performance sensitivity using the chain rule of differentiation. This method clearly shows the implicit dependence on the design, and a very simple sensitivity expression can be obtained.Consider a structure in which the generalized stiffness and mass matrices have been reduced by accounting for boundary conditions. Let the damping force be represented in the form ofC
(
b
)
z
,
t where z
,
t
=
d
z
/
d
t denotes the velocity vector. Under these conditions, Lagrange’s equation of motion becomes the second-order differential equation, as [23]
(16)
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
=
F
(
t
,
b
)
with the initial conditions
(17)
z
(
0
)
=
z
0
,
z
,
t
(
0
)
=
z
,
t
0
.
If design parameters are just related to stiffness matrix, we have
(18)
[
M
]
{
∂
z
,
t
t
∂
b
i
}
+
[
C
]
{
∂
z
,
t
∂
b
i
}
+
[
K
]
{
∂
z
∂
b
i
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
in which {
∂
z
/
∂
b
i
}, {
∂
z
,
t
/
∂
b
i
}, and {
∂
z
,
t
t
/
∂
b
i
} are sensitivity vectors of displacement, velocity, and acceleration with respect to design parameter b
i, respectively. Assume that(19a)
Y
,
t
t
=
∂
z
,
t
t
∂
b
i
,
(19b)
Y
,
t
=
∂
z
,
t
∂
b
i
,
(19c)
Y
=
∂
z
∂
b
i
.So, by replacing (19a), (19b), and (19c) to (18), we have
(20)
[
M
]
{
Y
,
t
t
}
+
[
C
]
{
Y
,
t
}
+
[
K
]
{
Y
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
.
The right side of (20) can be considered as an equivalent force, so (20) is similar to (16) and sensitivity vectors can be obtained by Newmark method.


### 3.3.2. Adjoint Variable Method

Sensitivity analysis can be performed very efficiently by using deterministic methods based on adjoint functions. The use of adjoint functions for analyzing the effects of small perturbations in a linear system was introduced by Wigner [24].This method constructs an adjoint problem that solves the adjoint variable, which contains all implicit dependent terms.For the dynamic response of structure, the following form of a general performance measure will be considered:(21)
ψ
=
g
(
z
(
T
)
,
b
)
+
∫
0
T
G
(
z
,
b
)
d
t
,
where the final time T is determined by a condition in the form:
(22)
Ω
(
z
(
T
)
,
z
,
t
(
T
)
,
b
)
=
0
.
It is presumed that (22) uniquely determines T, at least locally. This requires that the time derivative of Ω is nonzero at T as follows:
(23)
Ω
,
t
=
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
z
,
t
t
(
T
)
≠
0
.
When final time T is prescribed before the response analysis, the relation in (22) needs not be considered.To obtain the design sensitivity ofΨ, define a design variation in the form:
(24)
b
τ
=
b
+
τ
δ
b
.
Design b is perturbed in the direction of δ
b with the parameter τ. Substituting b
τ into (21), the derivative of (21) can be evaluated with respect to τ at τ
=
0. Leibnitz’s rule of differentiation of an integral may be used to obtain the following expression:
(25)
ψ
′
=
∂
g
∂
b
δ
b
+
∂
g
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
G
(
z
(
T
)
,
b
)
T
′
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
,
where
(26)
z
′
=
z
′
(
b
,
δ
b
)
≡
d
d
τ
z
(
t
,
b
+
τ
δ
b
)
|
τ
=
0
=
d
d
b
[
z
(
t
,
b
)
]
δ
b
,
T
′
=
T
′
(
b
,
δ
b
)
≡
d
d
τ
T
(
b
+
τ
δ
b
)
|
τ
=
0
=
d
T
d
b
δ
b
.
Note that since the expression in (21), that determines T, depends on the design, T will also depend on the design. Thus, terms arise in (25), that involve the derivative of T with respect to the design. In order to eliminate these terms, differentiate (22) with respect to τ and evaluate it at τ
=
0 in order to obtain
(27)
∂
Ω
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
∂
Ω
∂
z
,
t
[
z
,
t
′
(
T
)
+
z
,
t
t
(
T
)
T
′
]
+
∂
Ω
∂
b
δ
b
=
0
.
This equation may also be written as
(28)
Ω
,
t
T
′
=
[
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
t
(
T
)
]
T
′
=
-
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Since it is presumed by (23) that Ω
,
t
≠
0, then
(29)
T
′
=
-
1
Ω
,
t
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Substituting the result of (29) into (25), the following is obtained:
(30)
ψ
′
=
[
∂
g
∂
z
-
(
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
)
1
Ω
,
t
∂
Ω
∂
z
]
z
′
(
T
)
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
+
∂
g
∂
b
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
.
Note that ψ
′ depends on z
′ and z
,
t
′ at T, as well as on z
′ within the integration.In order to writeΨ
′ in (29) explicitly in terms of a design variation, the adjoint variable technique can be used. In the case of a dynamic system, all terms in (16) can be multiplied by λ
T
(
t
) and integrated over the interval [
0
,
T
] to obtain the following identity in λ:
(31)
∫
0
T
λ
T
[
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
-
F
(
t
,
b
)
]
d
t
=
0
.
Since this equation must hold for arbitrary λ, which is now taken to be independent of the design, substitute b
τ into (31) and differentiate it with respect to τ in order to obtain the following relationship:
(32)
∫
0
T
[
λ
T
M
(
b
)
z
,
t
t
′
+
λ
T
C
(
b
)
z
,
t
′
+
λ
T
K
(
b
)
z
′
-
∂
R
∂
b
δ
b
]
d
t
=
0
,
where
(33)
R
=
λ
~
T
F
(
t
,
b
)
-
λ
~
T
M
(
b
)
z
~
,
t
t
-
λ
~
T
C
(
b
)
z
~
,
t
-
λ
~
T
K
(
b
)
z
~
with the superposed tilde (
~
) denoting variables that are held constant during the differentiation with respect to the design in (32).Since (32), contains the time derivatives of z
′, integrate the first two integrands by parts in order to move the time derivatives to λ, as follows:
(34)
λ
T
M
(
b
)
z
,
t
′
(
T
)
-
λ
,
t
T
(
T
)
M
(
b
)
z
′
(
T
)
+
λ
T
C
(
b
)
z
′
(
T
)
+
∫
0
T
{
[
λ
,
t
t
T
M
(
b
)
-
λ
,
t
T
C
(
b
)
+
λ
T
K
(
b
)
]
z
′
-
∂
R
∂
b
δ
b
}
d
t
=
0
.
The adjoint variable method expresses the unknown terms in (30), in terms of the adjoint variable (
λ
). Since (34) must hold for arbitrary functions λ
(
t
), λ may be chosen so that the coefficients of terms involving z
′
(
T
), z
,
t
′
(
T
), and z
′ in (30) and (34) are equal. If such a function λ
(
t
) can be found, then the unwanted terms in (30), involving z
′
(
T
), z
,
t
′
(
T
) and z
′, can be replaced by terms that explicitly depend on δ
b in (34), and to be more specific, choose a λ
(
t
) that satisfies the following:
(35)
M
(
b
)
λ
(
T
)
=
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
t
,
(36)
M
(
b
)
λ
,
t
(
T
)
=
C
T
(
b
)
λ
(
T
)
-
∂
g
T
∂
z
+
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
(37)
M
(
b
)
λ
,
t
t
-
C
T
(
b
)
λ
,
t
+
K
(
b
)
λ
=
∂
G
T
∂
z
,
0
≤
t
≤
T
.
Note that once the dynamic equations of (16) and (17) is solved and (22) is used to determine T, then z
(
T
), z
,
t
(
T
), ∂
Ω
/
∂
z, ∂
Ω
/
∂
z
,
t, and Ω
,
t may be evaluated. Equation (23) can then be solved for λ
(
T
) since the mass matrix M
(
b
) is nonsingular. Having determined λ
(
T
), all terms on the right of (36) can be evaluated, and the equation can be solved for λ
,
t
(
T
). Thus, a set of terminal conditions on λ has been determined. Since M
(
b
) is nonsingular, (37) may then be integrated from T to 0, yielding the unique solution λ
(
t
). Taken as a whole, (35), through (37), may be thought of as a terminal value problem.Since the terms involving a variation in the state variable in (30) and (34) are identical, substitute (34) into (30) to obtain
(38)
ψ
′
=
∂
g
∂
b
δ
b
+
∫
0
T
[
∂
G
∂
b
+
∂
R
∂
b
]
d
t
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
≡
∂
ψ
∂
b
δ
b
.
Every term in this equation can now be calculated. The terms ∂
g
/
∂
b, ∂
G
/
∂
b, and ∂
Ω
/
∂
b represent explicit partial derivatives with respect to the design. The term ∂
R
/
∂
b, however, must be evaluated from (33), thus requiring λ
(
t
). Note also that since design variation δ
b does not depend on time, it is taken outside the integral in (38).Since (38) must hold, for all δ
b, the design derivative vector of ψ is
(39)
d
ψ
d
b
=
∂
g
∂
b
(
z
(
T
)
,
b
)
+
∫
0
T
[
∂
G
∂
b
(
z
,
b
)
+
∂
R
∂
b
(
λ
(
t
)
,
z
(
t
)
,
z
,
t
(
t
)
,
z
,
t
t
(
t
)
,
b
)
]
d
t
-
1
Ω
,
t
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
∂
Ω
∂
b
.


## 3.3.1. Direct Differentiation Method

The direct differentiation method (DDM) is a general, accurate, and efficient method to compute finite element response sensitivities to the model parameters. This method directly solves for the design dependency of a state variable and then computes performance sensitivity using the chain rule of differentiation. This method clearly shows the implicit dependence on the design, and a very simple sensitivity expression can be obtained.Consider a structure in which the generalized stiffness and mass matrices have been reduced by accounting for boundary conditions. Let the damping force be represented in the form ofC
(
b
)
z
,
t where z
,
t
=
d
z
/
d
t denotes the velocity vector. Under these conditions, Lagrange’s equation of motion becomes the second-order differential equation, as [23]
(16)
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
=
F
(
t
,
b
)
with the initial conditions
(17)
z
(
0
)
=
z
0
,
z
,
t
(
0
)
=
z
,
t
0
.
If design parameters are just related to stiffness matrix, we have
(18)
[
M
]
{
∂
z
,
t
t
∂
b
i
}
+
[
C
]
{
∂
z
,
t
∂
b
i
}
+
[
K
]
{
∂
z
∂
b
i
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
in which {
∂
z
/
∂
b
i
}, {
∂
z
,
t
/
∂
b
i
}, and {
∂
z
,
t
t
/
∂
b
i
} are sensitivity vectors of displacement, velocity, and acceleration with respect to design parameter b
i, respectively. Assume that(19a)
Y
,
t
t
=
∂
z
,
t
t
∂
b
i
,
(19b)
Y
,
t
=
∂
z
,
t
∂
b
i
,
(19c)
Y
=
∂
z
∂
b
i
.So, by replacing (19a), (19b), and (19c) to (18), we have
(20)
[
M
]
{
Y
,
t
t
}
+
[
C
]
{
Y
,
t
}
+
[
K
]
{
Y
}
=
-
∂
[
K
]
∂
b
i
{
z
}
-
α
2
∂
[
K
]
∂
b
i
{
z
,
t
}
.
The right side of (20) can be considered as an equivalent force, so (20) is similar to (16) and sensitivity vectors can be obtained by Newmark method.


## 3.3.2. Adjoint Variable Method

Sensitivity analysis can be performed very efficiently by using deterministic methods based on adjoint functions. The use of adjoint functions for analyzing the effects of small perturbations in a linear system was introduced by Wigner [24].This method constructs an adjoint problem that solves the adjoint variable, which contains all implicit dependent terms.For the dynamic response of structure, the following form of a general performance measure will be considered:(21)
ψ
=
g
(
z
(
T
)
,
b
)
+
∫
0
T
G
(
z
,
b
)
d
t
,
where the final time T is determined by a condition in the form:
(22)
Ω
(
z
(
T
)
,
z
,
t
(
T
)
,
b
)
=
0
.
It is presumed that (22) uniquely determines T, at least locally. This requires that the time derivative of Ω is nonzero at T as follows:
(23)
Ω
,
t
=
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
z
,
t
t
(
T
)
≠
0
.
When final time T is prescribed before the response analysis, the relation in (22) needs not be considered.To obtain the design sensitivity ofΨ, define a design variation in the form:
(24)
b
τ
=
b
+
τ
δ
b
.
Design b is perturbed in the direction of δ
b with the parameter τ. Substituting b
τ into (21), the derivative of (21) can be evaluated with respect to τ at τ
=
0. Leibnitz’s rule of differentiation of an integral may be used to obtain the following expression:
(25)
ψ
′
=
∂
g
∂
b
δ
b
+
∂
g
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
G
(
z
(
T
)
,
b
)
T
′
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
,
where
(26)
z
′
=
z
′
(
b
,
δ
b
)
≡
d
d
τ
z
(
t
,
b
+
τ
δ
b
)
|
τ
=
0
=
d
d
b
[
z
(
t
,
b
)
]
δ
b
,
T
′
=
T
′
(
b
,
δ
b
)
≡
d
d
τ
T
(
b
+
τ
δ
b
)
|
τ
=
0
=
d
T
d
b
δ
b
.
Note that since the expression in (21), that determines T, depends on the design, T will also depend on the design. Thus, terms arise in (25), that involve the derivative of T with respect to the design. In order to eliminate these terms, differentiate (22) with respect to τ and evaluate it at τ
=
0 in order to obtain
(27)
∂
Ω
∂
z
[
z
′
(
T
)
+
z
,
t
(
T
)
T
′
]
+
∂
Ω
∂
z
,
t
[
z
,
t
′
(
T
)
+
z
,
t
t
(
T
)
T
′
]
+
∂
Ω
∂
b
δ
b
=
0
.
This equation may also be written as
(28)
Ω
,
t
T
′
=
[
∂
Ω
∂
z
z
,
t
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
t
(
T
)
]
T
′
=
-
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Since it is presumed by (23) that Ω
,
t
≠
0, then
(29)
T
′
=
-
1
Ω
,
t
(
∂
Ω
∂
z
z
′
(
T
)
+
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∂
Ω
∂
b
δ
b
)
.
Substituting the result of (29) into (25), the following is obtained:
(30)
ψ
′
=
[
∂
g
∂
z
-
(
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
)
1
Ω
,
t
∂
Ω
∂
z
]
z
′
(
T
)
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
z
,
t
z
,
t
′
(
T
)
+
∫
0
T
[
∂
G
∂
z
z
′
+
∂
G
∂
b
δ
b
]
d
t
+
∂
g
∂
b
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
.
Note that ψ
′ depends on z
′ and z
,
t
′ at T, as well as on z
′ within the integration.In order to writeΨ
′ in (29) explicitly in terms of a design variation, the adjoint variable technique can be used. In the case of a dynamic system, all terms in (16) can be multiplied by λ
T
(
t
) and integrated over the interval [
0
,
T
] to obtain the following identity in λ:
(31)
∫
0
T
λ
T
[
M
(
b
)
z
,
t
t
+
C
(
b
)
z
,
t
+
K
(
b
)
z
-
F
(
t
,
b
)
]
d
t
=
0
.
Since this equation must hold for arbitrary λ, which is now taken to be independent of the design, substitute b
τ into (31) and differentiate it with respect to τ in order to obtain the following relationship:
(32)
∫
0
T
[
λ
T
M
(
b
)
z
,
t
t
′
+
λ
T
C
(
b
)
z
,
t
′
+
λ
T
K
(
b
)
z
′
-
∂
R
∂
b
δ
b
]
d
t
=
0
,
where
(33)
R
=
λ
~
T
F
(
t
,
b
)
-
λ
~
T
M
(
b
)
z
~
,
t
t
-
λ
~
T
C
(
b
)
z
~
,
t
-
λ
~
T
K
(
b
)
z
~
with the superposed tilde (
~
) denoting variables that are held constant during the differentiation with respect to the design in (32).Since (32), contains the time derivatives of z
′, integrate the first two integrands by parts in order to move the time derivatives to λ, as follows:
(34)
λ
T
M
(
b
)
z
,
t
′
(
T
)
-
λ
,
t
T
(
T
)
M
(
b
)
z
′
(
T
)
+
λ
T
C
(
b
)
z
′
(
T
)
+
∫
0
T
{
[
λ
,
t
t
T
M
(
b
)
-
λ
,
t
T
C
(
b
)
+
λ
T
K
(
b
)
]
z
′
-
∂
R
∂
b
δ
b
}
d
t
=
0
.
The adjoint variable method expresses the unknown terms in (30), in terms of the adjoint variable (
λ
). Since (34) must hold for arbitrary functions λ
(
t
), λ may be chosen so that the coefficients of terms involving z
′
(
T
), z
,
t
′
(
T
), and z
′ in (30) and (34) are equal. If such a function λ
(
t
) can be found, then the unwanted terms in (30), involving z
′
(
T
), z
,
t
′
(
T
) and z
′, can be replaced by terms that explicitly depend on δ
b in (34), and to be more specific, choose a λ
(
t
) that satisfies the following:
(35)
M
(
b
)
λ
(
T
)
=
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
t
,
(36)
M
(
b
)
λ
,
t
(
T
)
=
C
T
(
b
)
λ
(
T
)
-
∂
g
T
∂
z
+
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
T
∂
z
,
(37)
M
(
b
)
λ
,
t
t
-
C
T
(
b
)
λ
,
t
+
K
(
b
)
λ
=
∂
G
T
∂
z
,
0
≤
t
≤
T
.
Note that once the dynamic equations of (16) and (17) is solved and (22) is used to determine T, then z
(
T
), z
,
t
(
T
), ∂
Ω
/
∂
z, ∂
Ω
/
∂
z
,
t, and Ω
,
t may be evaluated. Equation (23) can then be solved for λ
(
T
) since the mass matrix M
(
b
) is nonsingular. Having determined λ
(
T
), all terms on the right of (36) can be evaluated, and the equation can be solved for λ
,
t
(
T
). Thus, a set of terminal conditions on λ has been determined. Since M
(
b
) is nonsingular, (37) may then be integrated from T to 0, yielding the unique solution λ
(
t
). Taken as a whole, (35), through (37), may be thought of as a terminal value problem.Since the terms involving a variation in the state variable in (30) and (34) are identical, substitute (34) into (30) to obtain
(38)
ψ
′
=
∂
g
∂
b
δ
b
+
∫
0
T
[
∂
G
∂
b
+
∂
R
∂
b
]
d
t
δ
b
-
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
1
Ω
,
t
∂
Ω
∂
b
δ
b
≡
∂
ψ
∂
b
δ
b
.
Every term in this equation can now be calculated. The terms ∂
g
/
∂
b, ∂
G
/
∂
b, and ∂
Ω
/
∂
b represent explicit partial derivatives with respect to the design. The term ∂
R
/
∂
b, however, must be evaluated from (33), thus requiring λ
(
t
). Note also that since design variation δ
b does not depend on time, it is taken outside the integral in (38).Since (38) must hold, for all δ
b, the design derivative vector of ψ is
(39)
d
ψ
d
b
=
∂
g
∂
b
(
z
(
T
)
,
b
)
+
∫
0
T
[
∂
G
∂
b
(
z
,
b
)
+
∂
R
∂
b
(
λ
(
t
)
,
z
(
t
)
,
z
,
t
(
t
)
,
z
,
t
t
(
t
)
,
b
)
]
d
t
-
1
Ω
,
t
[
∂
g
∂
z
z
,
t
(
T
)
+
G
(
z
(
T
)
,
b
)
]
∂
Ω
∂
b
.


## 3.4. Continuum Approach

In the continuum approach, the design derivative of the variational equation is taken before it is discretized. If the structural problem and sensitivity equations are solved as a continuum problem, then it is called the continuum-continuum method. The continuum sensitivity equation is solved by discretization in the same way that structural problems are solved. Since differentiation is taken at the continuum domain and is then followed by discretization, this method is called the continuum-discrete method.


## 3.5. Sensitivity Method Selection

The advantage of the finite difference method is obvious. If structural analysis can be performed and the performance measure can be obtained as a result of structural analysis, then the expressions in (14) and (15) are virtually independent of the problem types considered.Major disadvantage of the finite difference method is the accuracy of its sensitivity results. Depending on perturbation size, sensitivity results are quite different. For a mildly nonlinear performance measure, relatively large perturbation provides a reasonable estimation of sensitivity results. However, for highly nonlinear performances, a large perturbation yields completely inaccurate results. Thus, the determination of perturbation size greatly affects the sensitivity result. And even though it may be necessary to choose a very small perturbation, numerical noise becomes dominant for a too-small perturbation size. That is, with a too-small perturbation, no reliable difference can be found in the analysis results.The continuum-continuum approach is so limited and is not applicable in complex engineering structures because very simple, classical problems can be solved analytically.The discrete and continuum-discrete methods are equivalent under the conditions given below, using a beam as the structural component. It has also been argued that the discrete and continuum-discrete methods are equivalent under the conditions given below [23].First, the same discretization (shape function) used in the FEA method must be used for continuum design sensitivity analysis. Second, an exact integration (instead of a numerical integration) must be used in the generation of the stiffness matrix and in the evaluation of continuum-based design sensitivity expressions. Third, the exact solution (and not a numerical solution) of the finite element matrix equation and the adjoint equation should be used to compare these two methods. Fourth, the movement of discrete grid points must be consistent with the design parameterization method used in the continuum method.In this paper, two different analytical discrete methods, including direct differential method (DDM) and adjoint variable method (ADM) are presented and efficiency of proposed method is investigated when compared with DDM method.


## 4. Proposed Method

While structural vibration responses are used for damage detection, assumingG
=
0, (37), is a free vibration of beam with terminal conditions. Solving (37), for a single degree of freedom system is as follows:
(40)
m
λ
,
t
t
-
c
λ
,
t
+
k
λ
=
0
with terminal conditions:
λ
(
T
)
,
λ
˙
(
T
)
,
λ
T
(
t
)
=
e
ξ
ω
(
t
-
T
)
(
A
1
sin
(
ω
D
t
)
+
B
1
cos
⁡
(
ω
D
t
)
)
,
A
1
=
(
λ
,
t
(
T
)
ω
D
-
ξ
1
-
ξ
2
λ
(
T
)
)
cos
⁡
(
ω
D
T
)
+
λ
(
T
)
sin
⁡
(
ω
D
T
)
,
B
1
=
λ
(
T
)
cos
⁡
⁡
(
ω
D
t
)
-
A
1
tan
⁡
(
ω
D
T
)
,
in which
(41)
ξ
=
c
2
m
ω
=
c
c
cr
<
1
,
ω
D
=
ω
1
-
ξ
2
.
When time T is known, the coefficients of the characteristic equation of T
′ and thereupon Ω will be zero, so the terminal conditions are as follows:
(42)
λ
(
T
)
=
0
,
(43)
λ
,
t
(
T
)
=
M
-
1
(
b
)
×
(
-
∂
g
T
∂
z
)
.
Substitute (42) into (43) to obtain
(44)
A
1
=
λ
,
t
(
T
)
ω
D
cos
⁡
⁡
(
ω
D
T
)
,
B
1
=
-
λ
,
t
(
T
)
ω
D
sin
⁡
(
ω
D
T
)
.
Note that ∂
g
/
∂
z like A
1 and B
1 is dependent on time T, so terminal values for different amounts of T are not similar and adjoint equation should be calculated for all amounts of T separately. So,
(45)
λ
T
(
t
)
=
e
ξ
ω
(
t
-
T
)
(
λ
,
t
(
T
)
ω
D
cos
⁡
(
ω
D
T
)
sin
⁡
(
ω
D
t
)
m
m
m
m
m
m
-
λ
,
t
(
T
)
ω
D
sin
⁡
(
ω
D
T
)
cos
⁡
(
ω
D
t
)
)
=
P
T
f
(
t
)
+
Q
T
g
(
t
)
in which
(46)
P
T
=
e
-
ξ
ω
T
λ
,
t
(
T
)
ω
D
cos
⁡
(
ω
D
T
)
,
f
(
t
)
=
e
ξ
ω
t
sin
⁡
(
ω
D
t
)
,
Q
T
=
-
e
-
ξ
ω
T
λ
,
t
(
T
)
ω
D
sin
(
ω
D
T
)
,
g
(
t
)
=
e
ξ
ω
t
cos
⁡
⁡
(
ω
D
t
)
.


### 4.1. Sensitivity Matrix for Physical Parameter

Using (39) and assuming T is known and G
=
0 because of using structural vibration data, (47) can be obtained
(47)
d
ψ
d
b
=
∫
0
T
∂
R
∂
b
d
t
.In this equation,(48)
R
=
λ
~
T
F
(
t
)
-
λ
~
T
M
z
~
,
t
t
-
λ
~
T
C
(
b
)
z
~
,
t
-
λ
~
T
K
(
b
)
z
~
,
C
=
a
0
K
(
b
)
+
a
1
M
,
is Rayleigh damping matrix, so,
(49)
∂
R
∂
b
=
-
λ
T
~
a
0
∂
K
∂
b
z
,
t
~
-
λ
T
~
∂
K
∂
b
z
~
.
And finally component of sensitivity matrix in time T is
(50)
d
ψ
d
b
(
T
)
=
∫
0
T
(
-
λ
T
~
a
0
∂
K
∂
b
z
,
t
~
-
λ
T
~
∂
K
∂
b
z
~
)
d
t
.
In a multidegree of freedom problem, solving the above equations directly is not possible, and, for this purpose, change the variables as follows:
(51)
{
λ
}
=
[
ϕ
]
{
Y
}
.
In this equation, matrix [
ϕ
] forms vibration modes (modal matrix) and terminal conditions of above equations are
(52)
{
Y
(
T
)
}
=
M
-
1
[
ϕ
]
T
[
m
]
{
λ
(
T
)
}
,
(53)
{
Y
,
t
(
T
)
}
=
M
-
1
[
ϕ
]
T
[
m
]
{
λ
,
t
(
T
)
}
.
By inserting (51) in (37) and multiplying [
ϕ
]
T in both sides, the new equation in modal space is
(54)
[
M
]
{
Y
,
t
t
}
-
[
C
]
{
Y
,
t
}
+
[
K
]
{
Y
}
=
{
0
}
.
Each of [
M
], [
C
], and [
K
] matrices are diagonal, so,
(55)
M
i
{
Y
,
t
t
i
}
-
C
i
{
Y
,
t
i
}
+
K
i
{
Y
i
}
=
{
0
}
,
(56)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
Y
〉
×
[
ϕ
]
T
×
a
0
[
∂
k
∂
b
]
×
{
z
,
t
}
m
m
m
+
〈
Y
〉
×
[
ϕ
]
T
×
[
∂
k
∂
b
]
×
{
z
}
d
t
.
Consider
(57)
[
ϕ
]
T
×
a
0
[
∂
k
∂
b
]
×
{
z
,
t
}
=
{
z
z
,
t
}
,
[
ϕ
]
T
×
[
∂
k
∂
b
]
×
{
z
}
=
{
z
z
}
.
Equation (56) can be reduced to the following equation:
(58)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
Y
〉
×
{
z
z
,
t
}
+
〈
Y
〉
×
{
z
z
}
d
t
.
From (45), variable Y in modal space can be written as
(59)
{
Y
}
=
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
.
Replacing (59) in (58), a new expression is derived to calculate the sensitivity as follows:
(60)
d
ψ
d
b
(
T
)
=
-
∫
0
T
(
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
)
T
×
{
z
z
,
t
}
m
m
mm
+
(
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
)
T
×
{
z
z
}
d
t
.
Equation (60) can be rewritten as follows:
(61)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
P
(
T
)
〉
×
(
{
f
(
t
)
}
·
{
z
z
,
t
}
+
{
f
(
t
)
}
·
{
z
z
,
t
}
)
m
m
m
n
+
〈
Q
(
T
)
〉
×
(
{
g
(
t
)
}
·
{
z
z
,
t
}
+
{
g
(
t
)
}
·
{
z
z
,
t
}
)
d
t
.
Consider following parameters:
(62)
A
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
,
B
=
∫
0
T
{
g
(
t
)
}
·
{
z
z
,
t
}
d
t
,
C
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
}
d
t
,
D
=
∫
0
T
{
g
(
t
)
}
·
{
z
z
}
d
t
.
So, (61) is presented as
(63)
d
ψ
d
b
(
T
)
=
-
〈
P
(
T
)
〉
×
(
{
A
}
+
{
C
}
)
-
〈
Q
(
T
)
〉
×
(
{
B
}
+
{
C
}
)
.
The solution of (63) is directly too time-consuming, because in each time step all terms in (63) should be recalculated. Therefore, an incremental solution is developed as follows:
(64)
{
A
T
+
Δ
T
}
=
∫
0
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
+
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
,
(65)
{
A
T
+
Δ
T
}
=
{
A
T
}
+
{
δ
A
}
,
{
δ
A
}
=
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
≅
{
f
(
T
+
Δ
T
2
)
}
·
{
z
z
,
t
(
T
+
Δ
T
2
)
}
.
Similar to (65) for other parameters, we have
(66)
{
δ
B
}
=
∫
T
T
+
Δ
T
{
g
(
t
)
}
·
{
z
z
,
t
}
d
t
≅
{
g
(
T
+
Δ
T
2
)
}
·
{
z
z
,
t
(
T
+
Δ
T
2
)
}
,
{
δ
C
}
=
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
}
d
t
≅
{
f
(
T
+
Δ
T
2
)
}
·
{
z
z
(
T
+
Δ
T
2
)
}
,
{
δ
D
}
=
∫
T
T
+
Δ
T
{
g
(
t
)
}
·
{
z
z
}
d
t
≅
{
g
(
T
+
Δ
T
2
)
}
·
{
z
z
(
T
+
Δ
T
2
)
}
.
And finally the sensitivity expression in time T
+
Δ
T is as follows:
(67)
d
ψ
d
b
(
T
+
Δ
T
)
=
-
〈
P
(
T
+
Δ
T
)
〉
×
(
{
A
T
+
Δ
T
}
+
{
C
T
+
Δ
T
}
)
-
〈
Q
(
T
+
Δ
T
)
〉
×
(
{
B
T
+
Δ
T
}
+
{
D
T
+
Δ
T
}
)
.


### 4.2. Computational Algorithm

The computational algorithm that leads to the determination of sensitivity matrix is as follows.Step 1.
Calculateλ
,
t
(
T
) from (43).Step 2.
Calculateω, ω
D, and ϕ from and consider i
=
1.Step 3.
For theith element, calculate ∂
K
/
∂
b, z
z
,
t, and z
z and consider j
=
1.Step 4.
For thejth sensor and the corresponding Dof, calculate λ
,
t
(
T
) from Step 1 and Y
,
t
(
T
) from (53) and T
n
=
Δ
t and T
o
=
0.Step 5.
ConsiderA
=
B
=
C
=
D
=
0.Step 6.
CalculateT
m
=
T
0
+
(
Δ
t
/
2
) and calculate P
(
T
n
)
-
Q
(
T
n
)
-
f
(
T
m
)
-
g
(
T
m
) from  (45).Step 7.
Calculate{
δ
A
}, {
δ
B
}, {
δ
C
}, and {
δ
D
} from ((65)~(66)).Step 8.
Calculated
ψ
/
d
b
(
T
n
) from (67).Step 9.
IfT
n
<
T
final, consider T
0
=
T
n and T
n
=
T
n
+
Δ
t and go to Step 5, otherwise go to next step.Step 10.
Ifj
< number of sensors, consider j
=
j
+
1 and go to Step 4, otherwise go to next step.Step 11.
Ifi
< number of elements, consider i
=
i
+
1 and go to Step 3, otherwise finish.


### 4.3. Procedure of Iteration for Damage Detection

The initial analytical model of a structure deviates from the true model and measurement from the initial intact structure is used to update the analytical model. The improved model is then treated as a reference model, and measurement from the damaged structure will be used to update the reference model.When response measurement from the intact state of the structure is obtained, the sensitivities are computed from the proposed algorithm or direct differentiate method (20) based on the analytical model of the structure and the well knowing input force and velocity. The vector of parameter increments is then obtained from (8) or (10) using the computed and experimentally obtained responses. The analytical model is then updated and the corresponding response and its sensitivity are again computed for the next iteration. When measurement from the damaged state is obtained, the updated analytical model is used in the iteration in the same way as that using measurement from the intact state. Convergence is considered to be achieved when the following criteria are met as follows:
(68)
∥
E
i
+
1
-
E
i
∥
∥
E
i
∥
×
100
%
≤
Tol
1
,
∥
Response
i
+
1
-
Response
i
∥
∥
Response
i
∥
×
100
%
≤
Tol
2
.
The final vector of identified parameter increments corresponds to the changes occurring in between the two states of the structure. The tolerance is set equal to 1
×
10
-
6 in this study except otherwise specified.Equation (6) has been popularly used in the form of the first-order approximation of the increment on the left-hand side of the equation. The higher-order term of the Taylor expansion has been omitted in the computation. The iterative computation described above on the updating of the sensitivity and the system aims at reducing error due to such an omission, particularly with large local damages.


## 4.1. Sensitivity Matrix for Physical Parameter

Using (39) and assuming T is known and G
=
0 because of using structural vibration data, (47) can be obtained
(47)
d
ψ
d
b
=
∫
0
T
∂
R
∂
b
d
t
.In this equation,(48)
R
=
λ
~
T
F
(
t
)
-
λ
~
T
M
z
~
,
t
t
-
λ
~
T
C
(
b
)
z
~
,
t
-
λ
~
T
K
(
b
)
z
~
,
C
=
a
0
K
(
b
)
+
a
1
M
,
is Rayleigh damping matrix, so,
(49)
∂
R
∂
b
=
-
λ
T
~
a
0
∂
K
∂
b
z
,
t
~
-
λ
T
~
∂
K
∂
b
z
~
.
And finally component of sensitivity matrix in time T is
(50)
d
ψ
d
b
(
T
)
=
∫
0
T
(
-
λ
T
~
a
0
∂
K
∂
b
z
,
t
~
-
λ
T
~
∂
K
∂
b
z
~
)
d
t
.
In a multidegree of freedom problem, solving the above equations directly is not possible, and, for this purpose, change the variables as follows:
(51)
{
λ
}
=
[
ϕ
]
{
Y
}
.
In this equation, matrix [
ϕ
] forms vibration modes (modal matrix) and terminal conditions of above equations are
(52)
{
Y
(
T
)
}
=
M
-
1
[
ϕ
]
T
[
m
]
{
λ
(
T
)
}
,
(53)
{
Y
,
t
(
T
)
}
=
M
-
1
[
ϕ
]
T
[
m
]
{
λ
,
t
(
T
)
}
.
By inserting (51) in (37) and multiplying [
ϕ
]
T in both sides, the new equation in modal space is
(54)
[
M
]
{
Y
,
t
t
}
-
[
C
]
{
Y
,
t
}
+
[
K
]
{
Y
}
=
{
0
}
.
Each of [
M
], [
C
], and [
K
] matrices are diagonal, so,
(55)
M
i
{
Y
,
t
t
i
}
-
C
i
{
Y
,
t
i
}
+
K
i
{
Y
i
}
=
{
0
}
,
(56)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
Y
〉
×
[
ϕ
]
T
×
a
0
[
∂
k
∂
b
]
×
{
z
,
t
}
m
m
m
+
〈
Y
〉
×
[
ϕ
]
T
×
[
∂
k
∂
b
]
×
{
z
}
d
t
.
Consider
(57)
[
ϕ
]
T
×
a
0
[
∂
k
∂
b
]
×
{
z
,
t
}
=
{
z
z
,
t
}
,
[
ϕ
]
T
×
[
∂
k
∂
b
]
×
{
z
}
=
{
z
z
}
.
Equation (56) can be reduced to the following equation:
(58)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
Y
〉
×
{
z
z
,
t
}
+
〈
Y
〉
×
{
z
z
}
d
t
.
From (45), variable Y in modal space can be written as
(59)
{
Y
}
=
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
.
Replacing (59) in (58), a new expression is derived to calculate the sensitivity as follows:
(60)
d
ψ
d
b
(
T
)
=
-
∫
0
T
(
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
)
T
×
{
z
z
,
t
}
m
m
mm
+
(
{
P
(
T
)
}
·
{
f
(
t
)
}
+
{
Q
(
T
)
}
·
{
g
(
t
)
}
)
T
×
{
z
z
}
d
t
.
Equation (60) can be rewritten as follows:
(61)
d
ψ
d
b
(
T
)
=
-
∫
0
T
〈
P
(
T
)
〉
×
(
{
f
(
t
)
}
·
{
z
z
,
t
}
+
{
f
(
t
)
}
·
{
z
z
,
t
}
)
m
m
m
n
+
〈
Q
(
T
)
〉
×
(
{
g
(
t
)
}
·
{
z
z
,
t
}
+
{
g
(
t
)
}
·
{
z
z
,
t
}
)
d
t
.
Consider following parameters:
(62)
A
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
,
B
=
∫
0
T
{
g
(
t
)
}
·
{
z
z
,
t
}
d
t
,
C
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
}
d
t
,
D
=
∫
0
T
{
g
(
t
)
}
·
{
z
z
}
d
t
.
So, (61) is presented as
(63)
d
ψ
d
b
(
T
)
=
-
〈
P
(
T
)
〉
×
(
{
A
}
+
{
C
}
)
-
〈
Q
(
T
)
〉
×
(
{
B
}
+
{
C
}
)
.
The solution of (63) is directly too time-consuming, because in each time step all terms in (63) should be recalculated. Therefore, an incremental solution is developed as follows:
(64)
{
A
T
+
Δ
T
}
=
∫
0
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
=
∫
0
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
+
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
,
(65)
{
A
T
+
Δ
T
}
=
{
A
T
}
+
{
δ
A
}
,
{
δ
A
}
=
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
,
t
}
d
t
≅
{
f
(
T
+
Δ
T
2
)
}
·
{
z
z
,
t
(
T
+
Δ
T
2
)
}
.
Similar to (65) for other parameters, we have
(66)
{
δ
B
}
=
∫
T
T
+
Δ
T
{
g
(
t
)
}
·
{
z
z
,
t
}
d
t
≅
{
g
(
T
+
Δ
T
2
)
}
·
{
z
z
,
t
(
T
+
Δ
T
2
)
}
,
{
δ
C
}
=
∫
T
T
+
Δ
T
{
f
(
t
)
}
·
{
z
z
}
d
t
≅
{
f
(
T
+
Δ
T
2
)
}
·
{
z
z
(
T
+
Δ
T
2
)
}
,
{
δ
D
}
=
∫
T
T
+
Δ
T
{
g
(
t
)
}
·
{
z
z
}
d
t
≅
{
g
(
T
+
Δ
T
2
)
}
·
{
z
z
(
T
+
Δ
T
2
)
}
.
And finally the sensitivity expression in time T
+
Δ
T is as follows:
(67)
d
ψ
d
b
(
T
+
Δ
T
)
=
-
〈
P
(
T
+
Δ
T
)
〉
×
(
{
A
T
+
Δ
T
}
+
{
C
T
+
Δ
T
}
)
-
〈
Q
(
T
+
Δ
T
)
〉
×
(
{
B
T
+
Δ
T
}
+
{
D
T
+
Δ
T
}
)
.


## 4.2. Computational Algorithm

The computational algorithm that leads to the determination of sensitivity matrix is as follows.Step 1.
Calculateλ
,
t
(
T
) from (43).Step 2.
Calculateω, ω
D, and ϕ from and consider i
=
1.Step 3.
For theith element, calculate ∂
K
/
∂
b, z
z
,
t, and z
z and consider j
=
1.Step 4.
For thejth sensor and the corresponding Dof, calculate λ
,
t
(
T
) from Step 1 and Y
,
t
(
T
) from (53) and T
n
=
Δ
t and T
o
=
0.Step 5.
ConsiderA
=
B
=
C
=
D
=
0.Step 6.
CalculateT
m
=
T
0
+
(
Δ
t
/
2
) and calculate P
(
T
n
)
-
Q
(
T
n
)
-
f
(
T
m
)
-
g
(
T
m
) from  (45).Step 7.
Calculate{
δ
A
}, {
δ
B
}, {
δ
C
}, and {
δ
D
} from ((65)~(66)).Step 8.
Calculated
ψ
/
d
b
(
T
n
) from (67).Step 9.
IfT
n
<
T
final, consider T
0
=
T
n and T
n
=
T
n
+
Δ
t and go to Step 5, otherwise go to next step.Step 10.
Ifj
< number of sensors, consider j
=
j
+
1 and go to Step 4, otherwise go to next step.Step 11.
Ifi
< number of elements, consider i
=
i
+
1 and go to Step 3, otherwise finish.


## 4.3. Procedure of Iteration for Damage Detection

The initial analytical model of a structure deviates from the true model and measurement from the initial intact structure is used to update the analytical model. The improved model is then treated as a reference model, and measurement from the damaged structure will be used to update the reference model.When response measurement from the intact state of the structure is obtained, the sensitivities are computed from the proposed algorithm or direct differentiate method (20) based on the analytical model of the structure and the well knowing input force and velocity. The vector of parameter increments is then obtained from (8) or (10) using the computed and experimentally obtained responses. The analytical model is then updated and the corresponding response and its sensitivity are again computed for the next iteration. When measurement from the damaged state is obtained, the updated analytical model is used in the iteration in the same way as that using measurement from the intact state. Convergence is considered to be achieved when the following criteria are met as follows:
(68)
∥
E
i
+
1
-
E
i
∥
∥
E
i
∥
×
100
%
≤
Tol
1
,
∥
Response
i
+
1
-
Response
i
∥
∥
Response
i
∥
×
100
%
≤
Tol
2
.
The final vector of identified parameter increments corresponds to the changes occurring in between the two states of the structure. The tolerance is set equal to 1
×
10
-
6 in this study except otherwise specified.Equation (6) has been popularly used in the form of the first-order approximation of the increment on the left-hand side of the equation. The higher-order term of the Taylor expansion has been omitted in the computation. The iterative computation described above on the updating of the sensitivity and the system aims at reducing error due to such an omission, particularly with large local damages.


## 5. Numerical Results

To illustrate the formulations presented in the previous sections, we consider the system shown in Figures4 and 8, and capabilities of the proposed method are investigated.Figure 4
Multispan bridge model used in detection procedure.The relative percentage error (RPE) in the identified results is calculated from (69), where ∥
·
∥ is the norm of matrix and E
Identified and E
True are the identified and the true elastic modulus, respectively. Consider
(69)
RPE
=
∥
E
Identified
-
E
True
∥
∥
E
True
∥
×
100
%
.
Since the true value of elastic modulus is unknown, RPE can just be used for investigating the efficiency of method.


### 5.1. Multispan Model

A three-span bridge as shown in Figure4 is studied to illustrate the proposed method. It consists of 30 Euler-Bernoulli beam elements with 31 nodes each one with two Dof’s. The mass density of material is 7.8
×
10
3 kg/m3 and the elastic modulus of material is 2.1
×
10
7 N/cm2. The total length of bridge is 30 m and height and width of the frame section are, respectively, 200 and 200 mm. The first five undamped natural frequencies of the intact bridge are 37.73, 55.17, 66.97, 134.2, and 196.485 Hz. Rayleigh damping model is adopted with the damping ratios of the first two modes taken to be equal to 0.05. The equivalent Rayleigh coefficients, a
0 and a
1, are, respectively, 0.1 and 4.804
×
10
-
5.The transverse point loadP has a constant velocity, V
=
L
/
T, where T is the traveling time across the bridge and L is the total length of the bridge.For the forced vibration analysis, an implicit time integration method, called “the Newmark integration method” is used with the integration parametersβ
=
1
/
4 and γ
=
1
/
2, which leads to the constant-average acceleration approximation.Speed parameter is defined as(70)
α
v
=
V
V
cr
,
in which V
cr is critical speed (
V
cr
=
(
π
/
l
)
EI
/
ρ
), V is moving load speed, and ρ is mass per unit length of beam.


#### 5.1.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement noise are studied and they are shown in Table1.Table 1
Damage scenarios for multispan bridge.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M1-1
Single
23
5%
Nil
M1-2
Multi
8, 13, and 29
11%, 4%, and 7%
Nil
M1-3
Multi
3, 7, 19, 25, and 28
12%, 6%, 5%, 2%, and 18%
Nil
M1-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M1-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M1-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilLocal damage is simulated with a reduction in the elastic modulus of material of an element. The sampling rate is 10000 Hz and 450 data of the acceleration response (degree of indeterminacy is 15) collected along the z-direction at nodes 5, 15, and 25 are used in the identification.Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.088 in DDM method and 0.0354 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as in Scenario 6. This scenario consists of no simulated damage in the structure but with the initial elastic modulus of material of all the elements underestimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all scenarios (Figure5) and the RPE parameter is shown in Tables 2 and 3.Table 2
RPE of DDM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0465
0.0461
0.0457
0.0454
0.045
0.0743
0.0416
0.0471
0.088
M1-2
0.3135
0.317
0.3165
0.3157
0.315
0.2937
0.291
0.2967
0.0038
M1-3
0.0273
0.0268
0.0265
0.0262
0.0259
0.0281
0.0007
0.0007
0.0007
M1-4
0.052
0.0525
0.0516
0.0522
0.0531
0.0382
0.0576
0.0346
0.0155
M1-5
0.0411
0.0395
0.0408
0.0367
0.0403
0.06
0.0542
0.0207
0.0091
M1-6
0.0502
0.0546
0.0485
0.0471
0.0431
0.046
0.0422
0.041
0.0007Table 3
RPE of ADM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0354
0.0346
0.0342
0.0338
0.0337
0.0003
0.0214
0.0107
0.0079
M1-2
0.0496
0.0338
0.0493
0.0585
0.0575
0.0294
0.024
0.0396
0.0214
M1-3
0.0008
0.0005
0.0005
0.0005
0.0005
0.0007
0.0007
0.0007
0.0006
M1-4
0.0271
0.0247
0.0222
0.0077
0.0071
0.0161
0.0006
0.0287
0.0007
M1-5
0.0051
0.0047
0.0028
0.0035
0.0031
0.1971
0.0171
0.001
0.0134
M1-6
0.0526
0.0237
0.0156
0.0009
0.0008
0.065
0.001
0.0008
0.0007Figure 5
Detection of damage location and amount in elements 3, 7, 19, 25, and 28 and distribution of error in different elements with ADM scheme.Further studies on Scenario 6 shows that both methods are sensitive to the initial model error and for the maximum 20% initial error can be converged and a relatively good finite element model is therefore needed for the damage detection procedure.


#### 5.1.2. Effect of Noise

Noise is the random fluctuation in the value of measured or input that causes random fluctuation in the output value. Noise at the sensor output is due to either internal noise sources, such as resistors at finite temperatures, or externally generated mechanical and electromagnetic fluctuations [6].To evaluate the sensitivity of results to such measurement noise, noise-polluted measurements are simulated by adding to the noise-free acceleration vector a corresponding noise vector whose root-mean-square (r.m.s.) value is equal to a certain percentage of the r.m.s. value of the noise-free data vector. The components of all the noise vectors are of Gaussian distribution, uncorrelated, and with a zero mean and unit standard deviation. Then, on the basis of the noise-free accelerationZ
,
t
t
n
f, the noise-polluted acceleration Z
,
t
t
n
p of the bridge at location x can be simulated by
(71)
Z
,
t
t
n
p
=
Z
,
t
t
n
f
+
RMS
(
Z
,
t
t
n
f
)
×
N
level
×
N
unit
,
where RMS
(
Z
,
t
t
n
f
) is the r.m.s value of the noise-free acceleration vector, Z
,
t
t
n
f
×
N
level is the noise level, and N
unit is a randomly generated noise vector with zero mean and unit standard deviation [14].In order to study the effect of noise on stability of sensitivity methods, Scenario 2 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure have been studied.Results are illustrated in Figure6 for DDM and ADM methods, respectively.RPE contours with respect to noise level and loops.
(a)
(b)These contours show that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 1.3%, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than 1.3%, a denoising tool alongside sensitivity methods should be used.


#### 5.1.3. Efficiency of Proposed Method

In order to compare and quantify the performance of different methods and evaluate the proposed method, relative efficiency parameter (REP) is defined as follows:(72)
REP
=
ST
DDM
ST
ADM
in which ST is the solution time of SI method. In fact, this parameter represents the computation cost of method.Figure7 shows REP changes with respect to the speed parameter in different scenarios.Figure 7
REP changes in different scenarios with respect to speed parameter.Figure 8
Plane grid bridge model used in detection procedure.Table4 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 2.1599 and 12.3739 and its average is 4.6580; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 21.5% of other sensitivity based finite element model updating methods.Table 4
REP ranges in different scenarios.
Damage scenario
Max REP
Min REP
Average
M1-1
12.3739
4.9093
7.6744
M1-2
3.5953
2.2271
2.7166
M1-3
5.4912
4.5801
4.9990
M1-4
6.0214
2.287
4.6553
M1-5
3.8383
2.1599
3.1221
M1-6
7.6027
3.2449
4.7804
Total
12.3739
2.1599
4.6580


### 5.2. Plane Grid Model

A plane grid model of bridge is studied as another numerical example to illustrate the effectiveness of the proposed method. The finite element model of the structure is shown in Figure8. The structure is modeled by 46 frame elements and 32 nodes with three Dof at each node for the translation and rotational deformations. The mass density of material is 7.8
×
10
3 kg/m3 and the elastic modulus of material is 2.1
×
10
7 N/cm2. The first five undamped natural frequencies of the intact bridge are 45.59, 92.77, 181.74, 259.73, and 399.07 Hz. Rayleigh damping model is adopted with the damping ratios of the first two modes taken to be equal to 0.05. The equivalent Rayleigh coefficients a
0 and a
1 are, respectively, 0.1 and 2.364
×
10
-
5.


#### 5.2.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement of the noise are studied and they are shown in Table5.Table 5
Damage scenarios for grid model.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M2-1
Single
41
7%
Nil
M2-2
Multi
3, 26, 35, and 40
9%, 14%, 3%, and 8%
Nil
M2-3
Multi
5, 7, 12, 15, 24, and 37
4%, 11%, 6%, 2%, 10%, and 16%
Nil
M2-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M2-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M2-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilThe sampling rate is 14000 Hz and 460 data of the acceleration response (degree of indeterminacy is 10) collected along thez-direction at nodes 4, 11, 21, and 27 are used.Similar to the previous model, Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.0006 in DDM method and 0.0011 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as Scenario 6. This scenario consists of no simulated damage in the structure, but with the initial elastic modulus of material of all the elements under-estimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all the scenarios (Figure9) and the RPE parameter is shown in Tables 6 and 7.Table 6
RPE of DDM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0003
0.0003
0.0003
0.0003
0.0004
0.0004
0.0003
0.0006
0.0006
M2-2
0.0005
0.0006
0.0005
0.0004
0.0003
0.0004
0.0004
0.0005
0.0006
M2-3
0.0004
0.0004
0.0006
0.0003
0.0006
0.0005
0.0003
0.0005
0.0003
M2-4
0.0006
0.0006
0.0004
0.0005
0.0005
0.0005
0.0004
0.0002
0.0004
M2-5
0.0005
0.0006
0.0006
0.0004
0.0004
0.0003
0.0005
0.0004
0.0003
M2-6
0.0004
0.0004
0.0003
0.0005
0.0004
0.0004
0.0006
0.0003
0.0004Table 7
RPE of ADM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0002
0.0002
0.0001
0.0009
0.0011
0.0034
0.0014
0.0007
0.0007
M2-2
0.0006
0.0008
0.0004
0.0011
0.001
0.0014
0.0015
0.0012
0.0007
M2-3
0.0005
0.0007
0.0097
0.001
0.0011
0.0013
0.0018
0.001
0.0007
M2-4
0.0003
0.0003
0.0007
0.0013
0.0007
0.001
0.0012
0.0008
0.0011
M2-5
0.001
0.001
0.0008
0.0009
0.001
0.001
0.0014
0.0007
0.0012
M2-6
0.0007
0.0007
0.0007
0.0009
0.0011
0.0011
0.0011
0.0011
0.0011Figure 9
Detection of damage location and amount in elements 5, 7, 12, 15, 24, and 37 and distribution of error in different elements with ADM scheme.


#### 5.2.2. Effect of Noise

In order to study effect of noise on stability of sensitivity methods, scenario 3 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure has been studied.Figure10 shows that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 2% and 1.7% for ADM method and DDM method, respectively, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than mentioned values, a denoising tool such as wavelet transform alongside sensitivity methods should be used. The wavelet transform is mainly attractive because of its ability to compress and encode information to reduce noise or to detect any local singular behavior of a signal [25].RPE contours with respect to noise level and loops.
(a)
(b)


#### 5.2.3. Efficiency of Proposed Method

Figure11 shows REP changes with respect to the speed parameter in different scenarios. Table 8 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 1.4998 and 3.1370 and its average is 2.1173; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 47.2% of other sensitivity based finite element model updating methods.Table 8
REP ranges in different scenarios for model 2.
Damage scenario
Max REP
Min REP
Average
M2-1
2.423
1.4998
1.9089
M2-2
3.0713
1.8519
2.442633
M2-3
3.137
1.7166
2.443778
M2-4
2.8168
1.748
2.153533
M2-5
2.5382
1.6374
2.0865
M2-6
2.2976
1.4389
1.859456
Total
3.137
1.4389
2.117258Figure 11
REP changes in different scenarios with respect to speed parameter for model 2.


## 5.1. Multispan Model

A three-span bridge as shown in Figure4 is studied to illustrate the proposed method. It consists of 30 Euler-Bernoulli beam elements with 31 nodes each one with two Dof’s. The mass density of material is 7.8
×
10
3 kg/m3 and the elastic modulus of material is 2.1
×
10
7 N/cm2. The total length of bridge is 30 m and height and width of the frame section are, respectively, 200 and 200 mm. The first five undamped natural frequencies of the intact bridge are 37.73, 55.17, 66.97, 134.2, and 196.485 Hz. Rayleigh damping model is adopted with the damping ratios of the first two modes taken to be equal to 0.05. The equivalent Rayleigh coefficients, a
0 and a
1, are, respectively, 0.1 and 4.804
×
10
-
5.The transverse point loadP has a constant velocity, V
=
L
/
T, where T is the traveling time across the bridge and L is the total length of the bridge.For the forced vibration analysis, an implicit time integration method, called “the Newmark integration method” is used with the integration parametersβ
=
1
/
4 and γ
=
1
/
2, which leads to the constant-average acceleration approximation.Speed parameter is defined as(70)
α
v
=
V
V
cr
,
in which V
cr is critical speed (
V
cr
=
(
π
/
l
)
EI
/
ρ
), V is moving load speed, and ρ is mass per unit length of beam.


### 5.1.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement noise are studied and they are shown in Table1.Table 1
Damage scenarios for multispan bridge.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M1-1
Single
23
5%
Nil
M1-2
Multi
8, 13, and 29
11%, 4%, and 7%
Nil
M1-3
Multi
3, 7, 19, 25, and 28
12%, 6%, 5%, 2%, and 18%
Nil
M1-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M1-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M1-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilLocal damage is simulated with a reduction in the elastic modulus of material of an element. The sampling rate is 10000 Hz and 450 data of the acceleration response (degree of indeterminacy is 15) collected along the z-direction at nodes 5, 15, and 25 are used in the identification.Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.088 in DDM method and 0.0354 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as in Scenario 6. This scenario consists of no simulated damage in the structure but with the initial elastic modulus of material of all the elements underestimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all scenarios (Figure5) and the RPE parameter is shown in Tables 2 and 3.Table 2
RPE of DDM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0465
0.0461
0.0457
0.0454
0.045
0.0743
0.0416
0.0471
0.088
M1-2
0.3135
0.317
0.3165
0.3157
0.315
0.2937
0.291
0.2967
0.0038
M1-3
0.0273
0.0268
0.0265
0.0262
0.0259
0.0281
0.0007
0.0007
0.0007
M1-4
0.052
0.0525
0.0516
0.0522
0.0531
0.0382
0.0576
0.0346
0.0155
M1-5
0.0411
0.0395
0.0408
0.0367
0.0403
0.06
0.0542
0.0207
0.0091
M1-6
0.0502
0.0546
0.0485
0.0471
0.0431
0.046
0.0422
0.041
0.0007Table 3
RPE of ADM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0354
0.0346
0.0342
0.0338
0.0337
0.0003
0.0214
0.0107
0.0079
M1-2
0.0496
0.0338
0.0493
0.0585
0.0575
0.0294
0.024
0.0396
0.0214
M1-3
0.0008
0.0005
0.0005
0.0005
0.0005
0.0007
0.0007
0.0007
0.0006
M1-4
0.0271
0.0247
0.0222
0.0077
0.0071
0.0161
0.0006
0.0287
0.0007
M1-5
0.0051
0.0047
0.0028
0.0035
0.0031
0.1971
0.0171
0.001
0.0134
M1-6
0.0526
0.0237
0.0156
0.0009
0.0008
0.065
0.001
0.0008
0.0007Figure 5
Detection of damage location and amount in elements 3, 7, 19, 25, and 28 and distribution of error in different elements with ADM scheme.Further studies on Scenario 6 shows that both methods are sensitive to the initial model error and for the maximum 20% initial error can be converged and a relatively good finite element model is therefore needed for the damage detection procedure.


### 5.1.2. Effect of Noise

Noise is the random fluctuation in the value of measured or input that causes random fluctuation in the output value. Noise at the sensor output is due to either internal noise sources, such as resistors at finite temperatures, or externally generated mechanical and electromagnetic fluctuations [6].To evaluate the sensitivity of results to such measurement noise, noise-polluted measurements are simulated by adding to the noise-free acceleration vector a corresponding noise vector whose root-mean-square (r.m.s.) value is equal to a certain percentage of the r.m.s. value of the noise-free data vector. The components of all the noise vectors are of Gaussian distribution, uncorrelated, and with a zero mean and unit standard deviation. Then, on the basis of the noise-free accelerationZ
,
t
t
n
f, the noise-polluted acceleration Z
,
t
t
n
p of the bridge at location x can be simulated by
(71)
Z
,
t
t
n
p
=
Z
,
t
t
n
f
+
RMS
(
Z
,
t
t
n
f
)
×
N
level
×
N
unit
,
where RMS
(
Z
,
t
t
n
f
) is the r.m.s value of the noise-free acceleration vector, Z
,
t
t
n
f
×
N
level is the noise level, and N
unit is a randomly generated noise vector with zero mean and unit standard deviation [14].In order to study the effect of noise on stability of sensitivity methods, Scenario 2 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure have been studied.Results are illustrated in Figure6 for DDM and ADM methods, respectively.RPE contours with respect to noise level and loops.
(a)
(b)These contours show that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 1.3%, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than 1.3%, a denoising tool alongside sensitivity methods should be used.


### 5.1.3. Efficiency of Proposed Method

In order to compare and quantify the performance of different methods and evaluate the proposed method, relative efficiency parameter (REP) is defined as follows:(72)
REP
=
ST
DDM
ST
ADM
in which ST is the solution time of SI method. In fact, this parameter represents the computation cost of method.Figure7 shows REP changes with respect to the speed parameter in different scenarios.Figure 7
REP changes in different scenarios with respect to speed parameter.Figure 8
Plane grid bridge model used in detection procedure.Table4 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 2.1599 and 12.3739 and its average is 4.6580; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 21.5% of other sensitivity based finite element model updating methods.Table 4
REP ranges in different scenarios.
Damage scenario
Max REP
Min REP
Average
M1-1
12.3739
4.9093
7.6744
M1-2
3.5953
2.2271
2.7166
M1-3
5.4912
4.5801
4.9990
M1-4
6.0214
2.287
4.6553
M1-5
3.8383
2.1599
3.1221
M1-6
7.6027
3.2449
4.7804
Total
12.3739
2.1599
4.6580


## 5.1.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement noise are studied and they are shown in Table1.Table 1
Damage scenarios for multispan bridge.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M1-1
Single
23
5%
Nil
M1-2
Multi
8, 13, and 29
11%, 4%, and 7%
Nil
M1-3
Multi
3, 7, 19, 25, and 28
12%, 6%, 5%, 2%, and 18%
Nil
M1-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M1-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M1-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilLocal damage is simulated with a reduction in the elastic modulus of material of an element. The sampling rate is 10000 Hz and 450 data of the acceleration response (degree of indeterminacy is 15) collected along the z-direction at nodes 5, 15, and 25 are used in the identification.Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.088 in DDM method and 0.0354 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as in Scenario 6. This scenario consists of no simulated damage in the structure but with the initial elastic modulus of material of all the elements underestimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all scenarios (Figure5) and the RPE parameter is shown in Tables 2 and 3.Table 2
RPE of DDM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0465
0.0461
0.0457
0.0454
0.045
0.0743
0.0416
0.0471
0.088
M1-2
0.3135
0.317
0.3165
0.3157
0.315
0.2937
0.291
0.2967
0.0038
M1-3
0.0273
0.0268
0.0265
0.0262
0.0259
0.0281
0.0007
0.0007
0.0007
M1-4
0.052
0.0525
0.0516
0.0522
0.0531
0.0382
0.0576
0.0346
0.0155
M1-5
0.0411
0.0395
0.0408
0.0367
0.0403
0.06
0.0542
0.0207
0.0091
M1-6
0.0502
0.0546
0.0485
0.0471
0.0431
0.046
0.0422
0.041
0.0007Table 3
RPE of ADM method for model 1.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M1-1
0.0354
0.0346
0.0342
0.0338
0.0337
0.0003
0.0214
0.0107
0.0079
M1-2
0.0496
0.0338
0.0493
0.0585
0.0575
0.0294
0.024
0.0396
0.0214
M1-3
0.0008
0.0005
0.0005
0.0005
0.0005
0.0007
0.0007
0.0007
0.0006
M1-4
0.0271
0.0247
0.0222
0.0077
0.0071
0.0161
0.0006
0.0287
0.0007
M1-5
0.0051
0.0047
0.0028
0.0035
0.0031
0.1971
0.0171
0.001
0.0134
M1-6
0.0526
0.0237
0.0156
0.0009
0.0008
0.065
0.001
0.0008
0.0007Figure 5
Detection of damage location and amount in elements 3, 7, 19, 25, and 28 and distribution of error in different elements with ADM scheme.Further studies on Scenario 6 shows that both methods are sensitive to the initial model error and for the maximum 20% initial error can be converged and a relatively good finite element model is therefore needed for the damage detection procedure.


## 5.1.2. Effect of Noise

Noise is the random fluctuation in the value of measured or input that causes random fluctuation in the output value. Noise at the sensor output is due to either internal noise sources, such as resistors at finite temperatures, or externally generated mechanical and electromagnetic fluctuations [6].To evaluate the sensitivity of results to such measurement noise, noise-polluted measurements are simulated by adding to the noise-free acceleration vector a corresponding noise vector whose root-mean-square (r.m.s.) value is equal to a certain percentage of the r.m.s. value of the noise-free data vector. The components of all the noise vectors are of Gaussian distribution, uncorrelated, and with a zero mean and unit standard deviation. Then, on the basis of the noise-free accelerationZ
,
t
t
n
f, the noise-polluted acceleration Z
,
t
t
n
p of the bridge at location x can be simulated by
(71)
Z
,
t
t
n
p
=
Z
,
t
t
n
f
+
RMS
(
Z
,
t
t
n
f
)
×
N
level
×
N
unit
,
where RMS
(
Z
,
t
t
n
f
) is the r.m.s value of the noise-free acceleration vector, Z
,
t
t
n
f
×
N
level is the noise level, and N
unit is a randomly generated noise vector with zero mean and unit standard deviation [14].In order to study the effect of noise on stability of sensitivity methods, Scenario 2 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure have been studied.Results are illustrated in Figure6 for DDM and ADM methods, respectively.RPE contours with respect to noise level and loops.
(a)
(b)These contours show that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 1.3%, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than 1.3%, a denoising tool alongside sensitivity methods should be used.


## 5.1.3. Efficiency of Proposed Method

In order to compare and quantify the performance of different methods and evaluate the proposed method, relative efficiency parameter (REP) is defined as follows:(72)
REP
=
ST
DDM
ST
ADM
in which ST is the solution time of SI method. In fact, this parameter represents the computation cost of method.Figure7 shows REP changes with respect to the speed parameter in different scenarios.Figure 7
REP changes in different scenarios with respect to speed parameter.Figure 8
Plane grid bridge model used in detection procedure.Table4 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 2.1599 and 12.3739 and its average is 4.6580; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 21.5% of other sensitivity based finite element model updating methods.Table 4
REP ranges in different scenarios.
Damage scenario
Max REP
Min REP
Average
M1-1
12.3739
4.9093
7.6744
M1-2
3.5953
2.2271
2.7166
M1-3
5.4912
4.5801
4.9990
M1-4
6.0214
2.287
4.6553
M1-5
3.8383
2.1599
3.1221
M1-6
7.6027
3.2449
4.7804
Total
12.3739
2.1599
4.6580


## 5.2. Plane Grid Model

A plane grid model of bridge is studied as another numerical example to illustrate the effectiveness of the proposed method. The finite element model of the structure is shown in Figure8. The structure is modeled by 46 frame elements and 32 nodes with three Dof at each node for the translation and rotational deformations. The mass density of material is 7.8
×
10
3 kg/m3 and the elastic modulus of material is 2.1
×
10
7 N/cm2. The first five undamped natural frequencies of the intact bridge are 45.59, 92.77, 181.74, 259.73, and 399.07 Hz. Rayleigh damping model is adopted with the damping ratios of the first two modes taken to be equal to 0.05. The equivalent Rayleigh coefficients a
0 and a
1 are, respectively, 0.1 and 2.364
×
10
-
5.


### 5.2.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement of the noise are studied and they are shown in Table5.Table 5
Damage scenarios for grid model.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M2-1
Single
41
7%
Nil
M2-2
Multi
3, 26, 35, and 40
9%, 14%, 3%, and 8%
Nil
M2-3
Multi
5, 7, 12, 15, 24, and 37
4%, 11%, 6%, 2%, 10%, and 16%
Nil
M2-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M2-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M2-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilThe sampling rate is 14000 Hz and 460 data of the acceleration response (degree of indeterminacy is 10) collected along thez-direction at nodes 4, 11, 21, and 27 are used.Similar to the previous model, Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.0006 in DDM method and 0.0011 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as Scenario 6. This scenario consists of no simulated damage in the structure, but with the initial elastic modulus of material of all the elements under-estimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all the scenarios (Figure9) and the RPE parameter is shown in Tables 6 and 7.Table 6
RPE of DDM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0003
0.0003
0.0003
0.0003
0.0004
0.0004
0.0003
0.0006
0.0006
M2-2
0.0005
0.0006
0.0005
0.0004
0.0003
0.0004
0.0004
0.0005
0.0006
M2-3
0.0004
0.0004
0.0006
0.0003
0.0006
0.0005
0.0003
0.0005
0.0003
M2-4
0.0006
0.0006
0.0004
0.0005
0.0005
0.0005
0.0004
0.0002
0.0004
M2-5
0.0005
0.0006
0.0006
0.0004
0.0004
0.0003
0.0005
0.0004
0.0003
M2-6
0.0004
0.0004
0.0003
0.0005
0.0004
0.0004
0.0006
0.0003
0.0004Table 7
RPE of ADM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0002
0.0002
0.0001
0.0009
0.0011
0.0034
0.0014
0.0007
0.0007
M2-2
0.0006
0.0008
0.0004
0.0011
0.001
0.0014
0.0015
0.0012
0.0007
M2-3
0.0005
0.0007
0.0097
0.001
0.0011
0.0013
0.0018
0.001
0.0007
M2-4
0.0003
0.0003
0.0007
0.0013
0.0007
0.001
0.0012
0.0008
0.0011
M2-5
0.001
0.001
0.0008
0.0009
0.001
0.001
0.0014
0.0007
0.0012
M2-6
0.0007
0.0007
0.0007
0.0009
0.0011
0.0011
0.0011
0.0011
0.0011Figure 9
Detection of damage location and amount in elements 5, 7, 12, 15, 24, and 37 and distribution of error in different elements with ADM scheme.


### 5.2.2. Effect of Noise

In order to study effect of noise on stability of sensitivity methods, scenario 3 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure has been studied.Figure10 shows that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 2% and 1.7% for ADM method and DDM method, respectively, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than mentioned values, a denoising tool such as wavelet transform alongside sensitivity methods should be used. The wavelet transform is mainly attractive because of its ability to compress and encode information to reduce noise or to detect any local singular behavior of a signal [25].RPE contours with respect to noise level and loops.
(a)
(b)


### 5.2.3. Efficiency of Proposed Method

Figure11 shows REP changes with respect to the speed parameter in different scenarios. Table 8 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 1.4998 and 3.1370 and its average is 2.1173; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 47.2% of other sensitivity based finite element model updating methods.Table 8
REP ranges in different scenarios for model 2.
Damage scenario
Max REP
Min REP
Average
M2-1
2.423
1.4998
1.9089
M2-2
3.0713
1.8519
2.442633
M2-3
3.137
1.7166
2.443778
M2-4
2.8168
1.748
2.153533
M2-5
2.5382
1.6374
2.0865
M2-6
2.2976
1.4389
1.859456
Total
3.137
1.4389
2.117258Figure 11
REP changes in different scenarios with respect to speed parameter for model 2.


## 5.2.1. Damage Scenarios

Five damage scenarios of single, multiple, and random damages in the bridge without measurement of the noise are studied and they are shown in Table5.Table 5
Damage scenarios for grid model.
Damage scenario
Damage type
Damage location
Reduction in elastic modulus
Noise
M2-1
Single
41
7%
Nil
M2-2
Multi
3, 26, 35, and 40
9%, 14%, 3%, and 8%
Nil
M2-3
Multi
5, 7, 12, 15, 24, and 37
4%, 11%, 6%, 2%, 10%, and 16%
Nil
M2-4
Random
All elements
Random damage in all elements with an average of 5%
Nil
M2-5
Random
All elements
Random damage in all elements with an average of 15%
Nil
M2-6
Estimation of undamaged state
All elements
5% reduction in all elements
NilThe sampling rate is 14000 Hz and 460 data of the acceleration response (degree of indeterminacy is 10) collected along thez-direction at nodes 4, 11, 21, and 27 are used.Similar to the previous model, Scenario 1 studies the single damage scenario. The iterative solution converges in all speed parameter ranges with a maximum RPE of 0.0006 in DDM method and 0.0011 in ADM method.Scenarios 2 and 3 are on multiple damages with different amount of measured responses for the identification and Scenarios 4-5 are on random damages with different average for the identification. These scenarios also converge in all speed parameter ranges. One more scenario with model error is also included as Scenario 6. This scenario consists of no simulated damage in the structure, but with the initial elastic modulus of material of all the elements under-estimated by 5% in the inverse identification.Using both described methods, including DDM and proposed method, the damage locations and amount are identified correctly in all the scenarios (Figure9) and the RPE parameter is shown in Tables 6 and 7.Table 6
RPE of DDM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0003
0.0003
0.0003
0.0003
0.0004
0.0004
0.0003
0.0006
0.0006
M2-2
0.0005
0.0006
0.0005
0.0004
0.0003
0.0004
0.0004
0.0005
0.0006
M2-3
0.0004
0.0004
0.0006
0.0003
0.0006
0.0005
0.0003
0.0005
0.0003
M2-4
0.0006
0.0006
0.0004
0.0005
0.0005
0.0005
0.0004
0.0002
0.0004
M2-5
0.0005
0.0006
0.0006
0.0004
0.0004
0.0003
0.0005
0.0004
0.0003
M2-6
0.0004
0.0004
0.0003
0.0005
0.0004
0.0004
0.0006
0.0003
0.0004Table 7
RPE of ADM method for model 2.
Damage scenario
Speed parameter
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
M2-1
0.0002
0.0002
0.0001
0.0009
0.0011
0.0034
0.0014
0.0007
0.0007
M2-2
0.0006
0.0008
0.0004
0.0011
0.001
0.0014
0.0015
0.0012
0.0007
M2-3
0.0005
0.0007
0.0097
0.001
0.0011
0.0013
0.0018
0.001
0.0007
M2-4
0.0003
0.0003
0.0007
0.0013
0.0007
0.001
0.0012
0.0008
0.0011
M2-5
0.001
0.001
0.0008
0.0009
0.001
0.001
0.0014
0.0007
0.0012
M2-6
0.0007
0.0007
0.0007
0.0009
0.0011
0.0011
0.0011
0.0011
0.0011Figure 9
Detection of damage location and amount in elements 5, 7, 12, 15, 24, and 37 and distribution of error in different elements with ADM scheme.


## 5.2.2. Effect of Noise

In order to study effect of noise on stability of sensitivity methods, scenario 3 (speed ratio of moving load is considered to be constant and equal to 0.5) is considered and different levels of noise pollution are investigated, and RPE changes with increasing number of loops for the iterative procedure has been studied.Figure10 shows that both ADM and DDM methods are sensitive to the noise and if the noise level becomes greater than 2% and 1.7% for ADM method and DDM method, respectively, these methods lose their effectiveness and are not able to detect damage. So, in cases with noise level greater than mentioned values, a denoising tool such as wavelet transform alongside sensitivity methods should be used. The wavelet transform is mainly attractive because of its ability to compress and encode information to reduce noise or to detect any local singular behavior of a signal [25].RPE contours with respect to noise level and loops.
(a)
(b)


## 5.2.3. Efficiency of Proposed Method

Figure11 shows REP changes with respect to the speed parameter in different scenarios. Table 8 shows that in different scenarios and for different speed parameters, the efficiency parameter is between 1.4998 and 3.1370 and its average is 2.1173; therefore, the adjoint variable method is extremely successful and computational cost for this method is about 47.2% of other sensitivity based finite element model updating methods.Table 8
REP ranges in different scenarios for model 2.
Damage scenario
Max REP
Min REP
Average
M2-1
2.423
1.4998
1.9089
M2-2
3.0713
1.8519
2.442633
M2-3
3.137
1.7166
2.443778
M2-4
2.8168
1.748
2.153533
M2-5
2.5382
1.6374
2.0865
M2-6
2.2976
1.4389
1.859456
Total
3.137
1.4389
2.117258Figure 11
REP changes in different scenarios with respect to speed parameter for model 2.


## 6. Conclusion

A new damage detection method based on finite element model updating and sensitivity technique using acceleration time history data of a bridge deck affected by a moving vehicle with specified load, named “ADM” method, is presented. The updating procedure can be regarded as a parameter identification technique which aims to fit the unknown parameters of an analytical model such that the model behaviour corresponds as closely as possible to the measured behaviour.Newmark method is used to calculate the structural dynamic response and its dynamic response sensitivity matrix is calculated by adjoint variable method. In order to solve ill-posed inverse problem Tikhonov regularization method is used and L-curve method is implemented to find optimum value of the regularization parameter.In proposed method, an incremental solution for adjoint variable equation developed that calculates each element of sensitivity matrix separately. The main advantage is inclusion of an analytical method to augment the accuracy and speed of the solution.Numerical simulations demonstrate the efficiency and accuracy of the method to identify location and intensity of single, multiple, and random damages in different bridge models.Comparison studies confirmed that computational cost for this method is much lower than other traditional sensitivity methods. For modern, practical engineering applications, the cost of damage detection analysis is expensive. So, this method is feasible for large-scale problems.Similar to other sensitivity methods, the drawback of proposed method is its low stability against input measurement noise, which can be easily improved by using low-pass denoising tools such as wavelets.



---

*Source: 698658-2014-06-04.xml*