A dynamic stiffness-based framework for harmonic input estimation and response reconstruction considering damage

Abstract

In structural health monitoring (SHM), the measurement is point-wise but structures are continuous. Thus, input estimation has become a hot research subject with which the full-field structural response can be calculated with a finite element model (FEM). This paper proposes a framework based on the dynamic stiffness theory, to estimate harmonic input, reconstruct responses, and to localize damages from seriously deficient measurements. To begin, Fourier transform converts the dynamic equilibrium equation to an equivalent static one in the frequency domain, which is under-determined since the dimension of measurement vector is far less than the FEM-node number. The principal component analysis has been adopted to “compress” the under-determined equation, and formed an over-determined equation to estimate the unknown input. Then, inverse Fourier transform converts the estimated input in the frequency domain to the time domain. Applying this to the FEM can reconstruct the target responses. If a structure is damaged, the estimated nodal force can localize the damage. To improve the damage-detection accuracy, a multi-measurement-based indicator has been proposed. Numerical simulations have validated that the proposed framework can capably estimate input and reconstruct multi-types of full-field responses, and the damage indicator can localize minor damages even with the existence of noise.

Graphical abstract

Keywords

dynamic stiffness / principal component analysis / response reconstruction / damage localization / under-determined equation

1 Introduction

1.1 Literature review

In structural health monitoring (SHM) of the bridges, the external input, including vehicle and wind loads, is unknown, but important. However, the measurement is severely deficient compared with the number of degrees of freedom (DOFs). Computation of the real load distribution on a bridge from the monitored data, researchers have principally developed two methods of load estimation.

The first method can only estimate the vehicle loads, known as the bridge weight in motion (BWIM). BWIM extracts the vehicular information based on the principle of influence line that is generally obtained by a finite element model (FEM), or in-situ measurements [1]. The optimization algorithms have been adopted to minimize the difference between the measured and computed responses to optimize the vehicle information, including the axle weight, speed, and wheelbase [2,3]. Researchers have also developed the orthotropic BWIM algorithm, influence area method [1], and shear force method [4]. The main difficulties of BWIM include non-constant vehicle speed [5], multi-vehicle problem, lateral distribution [6], and dynamic load identification [7]. The moving force identification can settle the dynamic vehicle load problem [8,9], whereas, the other difficulties are still challenging.

The second method has been developed to apply to general loading distributions. At present, there are two types of expression. The first one [10] rewrites the dynamic equilibrium equation into the modal coordinates, obtains the impulse response matrix in the discrete-time by the Duhamel integral, calculates the modal force acting on the structure, and finally solves the ill-posed problem by regularization method [11] or iteration method [12]. The second type [13–15] adopts the expression of the state-space equation to estimate the input and optimize the estimated force by various approaches [16]. Gillijns and de Moor [17] firstly proposed a method of unbiased determination of system input with mathematical proof [18] when the dimension of the observation vector is equal to the DOFs. Also in complement, Hsieh [19] proposed an input estimation method when the dimension of the observation vector is less than the DOFs.

Since the state space expression has a similar expression with the Kalman filter theory, it has been widely used in input estimation. The extended Kalman filter [20], and robust Kalman filter [21], have been developed to identify both the linear system and the input simultaneously [22]. Besides, the updated Kalman filter applied to non-white noise [23] is proposed for load estimation [24,25] and damage identification [26,27]. The unbiased load estimation method under the modal coordinates and orthogonal decomposition [28] has also been developed.

In structural dynamic analysis, the harmonic load is an important load type but lacks study. It usually results in overly large vibration, such as vortex-induced resonance, and a structure tends to damage more under the harmonic input. Therefore, the focus of this study is the harmonic load.

In the former research of the authors [29,30], a static load estimation and deformation reconstruction method was proposed. This paper extends the method to the dynamic scope. The full-field response is computed using FEM and the estimated input.In addition, the estimated load can localize the structural damages based on the concept of the equivalent damage load.

1.2 Static equivalent damage load

Sun et al. [29,31] proposed the concept of the static equivalent damage load [32,33] to identify damage in beam-like structures. The mechanism is shown in Fig.1. Under identical static load q(x), the deflection curves of the undamaged and damaged beam are w^U and w^D, respectively. The magnitude of w^D is greater than w^U because the stiffness decreases due to damage. The equivalent damage load g(x) is computed from the difference between w^U and w^D and it is a virtual load. The deflection of the undamaged beam under the action g(x) is equal to w^D–w^U. Theoretical deductions demonstrate that the equivalent damage load exists merely in the damaged area, and thus, the equivalent damage load is a good damage indicator.

Inspired by the static equivalent damage load, this paper extends it to the dynamic scope using the dynamic stiffness theory. The theoretical basis of the dynamic stiffness theory is the Fourier transform, by which the dynamic equilibrium equation is transferred to an equivalent static equation in the frequency domain. As the governing equation in the frequency domain is similar to the static one, the static mechanics-based equivalent damage load applies to the dynamic scope as well.

The paper is arranged as follows. Section 2 is the theoretical basis where the basic assumptions, the dynamic stiffness theory, and the solution of the under-determined equation are introduced. Section 3 represents the response reconstruction process using a three-span continuous beam bridge model. Section 4 shows how to localize the damage using the proposed framework. Section 5 is the conclusion. It should be noted that the proposed method merely applies to the linear beam-like structures subjected to harmonic load.

2 Theoretical basis

2.1 Procedures and assumptions of the algorithm

This paper aims to obtain the position, frequency, and amplitude of the harmonic loads using a small number of sensors. With the estimated force and FEM, the full-field response is easy to calculate. The specific steps are listed as follows.

1) Convert the dynamic equilibrium equation into the frequency domain by the Fourier transform.

2) Form an under-determined equation about the unknown input.

3) Adopt the principal component analysis to solve the under-determined equation and estimate the equivalent nodal force in the frequency domain.

4) Use the inverse Fourier transform to transfer the equivalent nodal force back to the time domain.

5) Localize the damage from the estimated force in the time domain.

6) Apply the estimated force to the FEM and calculate the full-field responses.

Basic assumptions of the proposed algorithm include:

1) the position, frequency, and amplitude of the input do not change;

2) the FEM without damage is known;

3) the studied structure is linear.

For different structures, difficulty in satisfying the above assumptions varies. For civil structures, the input is complex and only the artificial exciter satisfies the assumptions. For mechanical and aerospace equipment, the basic assumptions are easier to provide.

2.2 Fourier transform

The Fourier transform is an effective tool to analyze the linear time-invariant system. For completeness, the formulas of the Fourier transform are briefly introduced:

(1)\mathcal{F}[ f(t)]= F(\omega )={\int }_0^{\mathrm{\infty }} f(t){\text{e}}^{-\text{j}\omega t}\text{d}t,

(2){\mathcal{F}}^{-1}[F (\omega )]=f (t)=\frac{1}{2\text{π}}{\int }_{-\text{j}\mathrm{\infty }}^{\text{j}\mathrm{\infty }}F(\omega ){\text{e}}^{\text{j}\omega t}\text{d}\omega ,

where \text{j}={(-1)}^{1/2}. Equations (1) and (2) are rewritten as:

(3)\{\begin{array}{c}\mathcal{F}[f(t)]=F(\omega ),\\ {\mathcal{F}}^{-1}[F(\omega )]=f(t).\end{array}

If f(t) in Eq. (1) is the derivative with respect to time, there is:

(4)\mathcal{F}\left[\frac{\text{d}f(t)}{\text{d}t}\right]=\text{j}\omega F(\omega ),

Equation (4) indicates that the Fourier transform can transfer the differential calculation in the time domain to the multiplication in the frequency domain. It is noted that the Eq. (4) is satisfied only with a zero-initial condition, i.e., f(0) = 0. For engineering problems, it is difficult to provide a zero-initial condition. With the non-zero initial condition, Eq. (4) becomes:

(5)\mathcal{F}\left[\frac{\text{d}f(t)}{\text{d}t}\right]=\text{j}\omega F(\omega )-f(0).

With the non-zero initial condition, the expression of the Fourier transform for the second-order derivative becomes:

(6)\begin{array}{rl}\mathcal{F}\left[\frac{{\text{d}}^{2}f(t)}{\text{d}{t}^2}\right]& ={{\text{e}}^{-\text{j}\omega t}\frac{\text{d}f(t)}{\text{d}t}|}_0^{\mathrm{\infty }}+\text{j}\omega {\int }_0^{\mathrm{\infty }}\frac{\text{d}f(t)}{\text{d}t}{\text{e}}^{-\text{j}\omega t}\text{d}t\\ & =-f^{\prime }(0)+\text{j}\omega [\text{j}\omega F(\text{j}\omega )-f(0)]\\ & =-{\omega }^{2}F(\omega )-\text{j}\omega f(0)-f^{\prime }(0).\end{array}

The displacement, velocity, and acceleration can convert to each other by the Fourier transform. In practice, the fast Fourier transform (FFT) has high computational efficiency and is widely used for numerical computations.

2.3 Dynamic stiffness theory

For a linear mechanical system, the governing equation is:

(7)\mathit{M}\ddot{\mathit{x}}+\mathit{C}\dot{\mathit{x}}+\mathit{K}\mathit{x}=\mathit{p}(t),

where M, C, and K \in {\mathbb{R}}_{n\times n} are the mass, damping, and stiffness matrices, respectively. n is the number of DOFs. p(t) \in {\mathbb{R}}_{n\times 1} is the nodal force vector at time moment t. With the zero-initial condition, the velocity and acceleration after the Fourier transform can be expressed as:

(8)\begin{array}{rl}& \mathit{X}=\mathcal{F}(\mathit{x}),\\ & \mathit{X}=\frac{1}{\text{j}\omega }\mathcal{F}(\dot{\mathit{x}}),\\ & \mathit{X}=\frac{1}{-{\omega }^2}\mathcal{F}(\ddot{\mathit{x}}),\end{array}

where X \in {\mathbb{C}}_{n\times 1} is the displacement vector in the frequency domain. The dynamic equilibrium equation is converted to an equivalent static equilibrium equation:

(9)(-{\omega }^2\mathit{M}+\text{j} \omega \mathit{C}+\mathit{K})\mathit{X}=\mathit{P},

where P = \mathcal{F}[M(t)]. In Eq. (9), {(-\omega }^2\mathit{M}+\text{j}\omega \mathit{C}+\mathit{K})\in {\mathbb{C}}_{n\times n} is the dynamic stiffness D(\omega):

(10)\mathit{D}(\omega )=-{\omega }^2\mathit{M}+\text{j}\omega \mathit{C}+ \mathit{K}.

Under the non-zero initial condition, Equation (8) is not provided. Considering the relationship in Equations (5) and (6), the expressions in Eq. (9) become:

(11)\begin{array}{r}\mathit{M}[-{\omega }^2\mathit{X}-\text{j}\omega \mathit{x}(0)-{\mathit{x}}^{\prime }(0)]+\mathit{C}[\text{j}\omega \mathit{X}-\mathit{x}(0)]+\mathit{K}\mathit{X}=\mathit{P}.\end{array}

Equation (11) has lost the conciseness compared to Eq. (9). Define:

(12)\tilde{\mathit{P}}=\mathit{P}+\text{j}\omega \mathit{M}\mathit{x}(0 )+{\mathit{M}\mathit{x}}^{\prime }(0) +\mathit{C}\mathit{x}(0) .

Equation (11) can be rewrriten as:

(13)(-{\omega }^2\mathit{M}+\text{j} \omega \mathit{C}+\mathit{K})\mathit{X}=\tilde{\mathit{P}}.

Equation (13) is an equivalent static equilibrium equation with the non-zero-initial condition. When merely the velocity or acceleration is available, Equation (13) becomes:

(14)\left(\text{j}\omega \mathit{M}+\mathit{C}+\frac{1}{\text{j}\omega }\mathit{K}\right) {\mathit{X}}_1=\tilde{\mathit{P}},

(15)\left(\mathit{M}+\frac{1}{\text{j}\omega }\mathit{C}+\frac{1}{-{\omega }^2}\mathit{K}\right){\mathit{X}}_2= \tilde{\mathit{P}},

where X₁ and X₂ \in {\mathbb{C}}_{n\times 1} are the velocity and acceleration vectors in the frequency domain, respectively.

2.4 Solution of the under-determined equation

The Fourier transform transfers the dynamic equilibrium equation to an equivalent static one in the frequency domain:

(16)\mathit{D}(\omega )\mathit{X}(\omega )=\mathit{P}(\omega ).

The inverse matrix of the dynamic stiffness D is the dynamic flexibility matrix F \in {\mathbb{C}}_{n\times n}. F is also known as the transfer function in system identification researches. Equation (16) is rewritten as:

(17)\mathit{X}(\omega )=\mathit{F}(\omega )\mathit{P}(\omega ).

X(\omega) is an n \times 1 vector, and whereas, the known elements of X(\omega) are at the DOFs equipped with sensors. In SHM, the dimension of observation is far less than the DOFs. Assuming the number of the measured displacements is N, the rows corresponding to unknown elements of X(\omega) are eliminated. Then, Eq. (17) becomes:

(18){\mathit{X}}_{N\times 1}(\omega )={\mathit{F}}_{N\times n}( \omega ){\mathit{P}}_{n\times 1}(\omega ),

where X is a vector of the knowns and P is the vector of unknowns. Since the dimension of knowns is far less than the unknowns, Eq. (18) is an under-determined equation that cannot be solved directly. The principal component analysis is adopted to calculate the external nodal force P(\omega) for each frequency \omega.

The principal component analysis (PCA) can “condense” the unknowns. The column number n of {\mathit{F}}_{N\times n} is greater than the row number N. Thus, the columns of {\mathit{F}}_{N\times n} are linear correlated. Firstly, the columns of {\mathbf{F}}_{N\times n} are normalized to provide consistency in the magnitudes. Then, the PCA [34,35] is adopted to calculate the maximal independent columns of the dynamic flexural matrix {\mathit{F}}_{N\times n}:

(19){\mathit{R}}_{N\times p}={\mathit{F}}_{N\times n}{\mathit{T}}_{n\times p},

(20){\mathit{R}}_{N\times p}=[ {\mathit{r}}_1,{\mathit{r}}_2,...,{\mathit{r}}_p],

where p is the number of principal components ( p < N ), T is the transferring matrix for computing the principal components, and r_i is the ith principal component column. The matrix R involves the first pth principal components. With Eq. (19), there is:

(21)\mathit{F}= \mathit{R}{\mathit{T}}^{\mathrm{T}}[\mathit{T}{\mathit{T}}^{\mathrm{T}}{]}^{-1}.

Substitute Eq. (21) into Eq. (18):

(22)\mathit{X}= \mathit{R}{\mathit{T}}^{\mathrm{T}}[\mathit{T}{\mathit{T}}^{\mathrm{T}}{]}^{-1}\mathit{P}.

Record T^T[TT^T]^–1P as \mathit{P}{}^{\prime }. There is:

(23){\mathit{X}}_{N\times 1}={\mathit{R}}_{N\times p}\mathit{P}{}_{p\times 1}^{\prime }.

Equation (23) is an over-determined equation when N > p and {\mathit{P}}^{\prime } can be determined by the least square estimation. Finally, the equivalent nodal force \overline{\mathit{P}} is calculated by:

(24)\overline{\mathit{P}}={\mathit{T}\mathit{P}}^{\prime }.

As matrix F is previously normalized, \overline{\mathit{P}} should be scaled back based on the normalization factors. It should be noted that {\overline{\mathit{P}}}_{n\times 1}\left(\omega \right) is the equivalent nodal force calculated from the displacement {\mathit{X}}_{N\times 1} in the frequency domain, and {\overline{\mathit{P}}}_{n\times 1}\left(\omega \right) is not perfectly identical to the real load P\left(\omega \right). However, the structural displacement {\mathit{X}}_{N\times 1} is provided to be identical under \overline{\mathit{P}}\left(\omega \right) and P\left(\omega \right), and the distribution of \overline{\mathit{P}}\left(\omega \right) is similar to that of P\left(\omega \right) for the undamaged structure. Thus, we call \overline{\mathit{P}}\left(\omega \right) as the equivalent nodal force in the frequency domain. The inverse Fourier transform in Eq. (2) is then used to calculate \overline{\mathit{p}} \left(t\right), i.e., the equivalent nodal force in the time domain.

It should be noted that F is a complex matrix and the PCA cannot directly tackle it. The damping matrix C is assumed as zero and F(\omega) becomes a real matrix. Then, the real and imaginary components of P(\omega) can be separately estimated based on the PCA approach.

When the external input is unknown, \overline{\mathit{p}}\left(t\right) can directly indicate the real input’s basic characteristics, including its frequency, amplitude, and position. The complete structural responses, such as deflection, inclination, and strain, can be reconstructed by applying the equivalent nodal force to the FEM. If the velocity and acceleration of the structure are required, the sampling frequency should be increased to provide accuracy.

2.5 Damage localization

If the structure is damaged, its static and dynamic stiffness decreases, and its dynamic response is greater than the intact structure under identical input. When estimating the input, the stiffness matrix is extracted from the intact FEM because the damage detail is unknown, and thus, the intact FEM-estimated input differs from the real input. The estimated equivalent nodal force \overline{\mathit{p}}\left(t\right) involves two components: the real input q(t) and the equivalent damage load g(t):

(25)\overline{\mathit{p}}(t)=\mathit{q}(t) +\mathit{g}(t).

The spatial distribution of q(t) keeps consistent with the real external load p(t). In the damaged area, \overline{\mathit{p}}\left(t\right) differs from p(t) since the equivalent damage load g(t) exists. Thus, the distribution of \overline{\mathit{p}}\left(t\right) can indicate the damage position by observing where g(t) exists. In some cases, such as the in-field experiment, the input position is known and it is easy to detect the existence of g(t). In other conditions, the input position is unknown. The next step is to distinguish g(t) from q(t) with unknown input position.

The deflection, inclination, and strain can reflect the deformation of a structure from varied scales. The inclination is the integration of strain, and the deflection is the integration of inclination. Therefore, the deflection and inclination measurements are global measurements, while the strain merely contains the local stiffness information. The difference in information between the deflection (or inclination) and the strain can be used to localize the damage. The inclination, rather than the deflection, is chosen as the measurement in this paper because it is easier to measure and is free of a reference point.

First, with N inclinometers, the unknown input \overline{\mathit{p}}\left(t\right) is estimated and recorded as \overline{\mathit{p}}\left(t\right)₁.

When identifying the damage, totally N₁ strain gauges have been equipped on the bridge. At each sampling moment, the governing equation is:

(26)\text{ε}(t) =\mathit{S}\mathit{p}(t) ,

where \text{ε}\in {\mathbb{R}}_{N1\times 1} is the vector of the strain measurement, S\in {\mathbb{R}}_{N1\times n} is the influence line matrix of strain, and p is the generalized force vector involving the inertial, damping, and external forces.

In practice, the strain gauge number is smaller than the FEM-DOFs, and therefore, Eq. (26) is an under-determined equation about p(t). Since the form of Eq. (26) is identical to Eq. (18), the proposed PCA-based approach can estimate the unknown generalized p(t) at each sampling moment t. The generalized input estimated by strain is recorded as \mathit{p}(t{)}^{\prime }. By multiplying \mathit{p}(t{)}^{\prime } with the influence line matrix of the inclination (which is extracted from the undamaged FEM), the inclination of the whole bridge can be calculated at each sampling moment. With the strain-estimated inclination, the algorithm in Section 2.4 can work again and the calculated equivalent nodal force is recorded as \overline{\mathit{p}}\left(t\right)₂. It should be noted that \overline{\mathit{p}}\left(t\right)₂ is merely based on the strain measurement. When the strain gauges are not equipped exactly on the damage position (which is the focus of this paper because the damage is easy to find when a strain gauge is exactly situated on it), the damage information is not recorded by the strain gauges. Then, the damage indicator is obtained:

(27)\mathrm{\Delta }\mathit{p}=\overline{\mathit{p}} (t{)}_1-\overline{\mathit{p}} (t{)}_2.

\mathrm{\Delta }p involves the information of both the inclination and strain. Merely the inclination involves the damage information when the strain gauges were not equipped at the damage. This means \overline{\mathit{p}}\left(t\right)₂ merely contains the input information, and whereas, \overline{\mathit{p}} \left(t\right)₁ contains both the input and damage information. The computing formula of \mathrm{\Delta }p eliminates the input information in \overline{\mathit{p}}\left(t\right)₁ and only the damage information is retained in \mathrm{\Delta }p.

3 Response reconstruction

3.1 Introduction of the numerical model

This study adopts a three-span continuous beam bridge as the numerical model (see Fig.2). The bending stiffness of the cross-section is 1.24 × 10⁷ N·m², and the linear density is 75 kg/m. The cross-section is rectangular whose width and height are 1.2 and 0.28 m, respectively. The bearings are all hinged. The bridge is divided into 44 elements with a total of 45 nodes. The element number and node number increase from left to right. Considering merely in-plane deformation, the total DOFs is 130. Fig.3 is the first two modes. The Rayleigh damping is assumed for the bridge, and the first two damping ratios are 1% and 2%, respectively.

This study uses rotations as “measurement”. The deflection and inclination are translational and rotational displacements. However, the inclination has some advantages over deflection. 1) The inclination is the derivative of deflection, and thus, it is more sensitive to damages than deflection. 2) Inclination is easier to measure than deflection because it is free of a reference point. There are fifteen inclinometers equipped on the bridge (see Fig.2).

This section involves different loading conditions to validate the applicability of the algorithm for reconstructing structural responses. It should be noted that when reconstructing the responses, the bridge is undamaged.

3.2 Load estimation and response reconstruction

To excite the responses of more modes, a concentrated harmonic load is applied to node No. 18. The input amplitude and frequency are 1000 N and 15 Hz. The dynamic inclinations at sensor-DOFs are calculated as “measurement”. The sampling frequency is 1000 Hz. For improving the computational efficiency of FFT, the data length is 16384.

The bridge response at the beginning two seconds contains the impulse response component (see Fig.4) and is unstable. Thus, merely the data during the steady response stage (after 2 s) is used to estimate the input. With the estimated input, the corresponding responses can be calculated from the FEM. The reconstructed responses involve both the impulse and steady responses. For the studied model bridge, the vertical load is dominant. Thus, the columns corresponding to the rotational and horizontal DOFs are eliminated in the dynamic flexibility matrix in Eq. (17), and the calculated nodal force is merely in the vertical direction.

Fig.5 is the map of the estimated nodal force vectors \overline{\mathit{p}} \left(t\right). The distribution of \overline{\mathit{p}} \left(t\right) on the bridge is similar to the real load. Values of \overline{\mathit{p}}\left(t\right) are greater near node No. 18 and smaller at other positions. Fig.6 is the distribution of \overline{\mathit{p}}(t= 2.52\mathrm{s}) when it reaches a peak over time. The maximum of \overline{\mathit{p}} in Fig.6 appears near node No. 18 (the real loading position), while the values at other nodes are smaller. Negative forces appear at nodes Nos. 22–26 which is the inherent limitation of the algorithm because the equivalent force \overline{\mathit{p}} is estimated from deficient measurements. The resultant force between node No. 14 and No. 26 is 997 N, which is close to the real load amplitude of 1000 N. Therefore, the bridge responses under the real and estimated loads are identical. One can notice that the estimated force is the smoothness of the real load longitudinally. The erroneous forces have a slight influence on the reconstructed response since they are smaller in value. Fig.7 is the frequency spectrum of the nodal force at node No. 18, which accurately retains the frequency characteristics of real input.

By applying the nodal force \overline{\mathit{p}} \left(t\right) to the FEM, the full-field responses are reconstructed. Fig.4 shows the reconstructed deflections at the bridge center, and Fig.8 is the deflection curve when the deflection reaches the peak in the steady response stage. Obviously, the reconstructed deflection is accurate over both time and space. The minor phase error in Fig.4 can be eliminated by moving the time axis, which does not affect other works in SHM.

3.3 Performance under complex harmonic load

The input can be complex for the in-field bridges when the frequency components at different positions are varied. In this section, a more complex loading condition is adopted. Two concentrated harmonic forces are applied to nodes No. 7 and No. 18 with frequencies of 5 and 15 Hz, respectively. The amplitudes of both of them are 500 N.

Fig.9 to 12 are the estimated nodal force and the reconstructed responses. In Fig.9, the positions of the two concentrated forces are easy to find where two peaks appear near nodes No. 7 and No. 18. Intuitively, one can deduce that the load is applied to these positions. Besides, the estimated force accurately retains the frequency information of the real external load. The frequency of the estimated force accords well with the real forces at nodes No. 7 and No. 18 (see Fig.10). By applying the estimated force to the FEM, the complete responses can be obtained (see Fig.11 and Fig.12) where the maximum errors are 4.2% and 3.6%, respectively. The reconstructed response is accurate over time and space, regardless of the complex input condition.

• Strain reconstruction

Strain is important in SHM because damage always relates to the over-large strain responses. However, the number of strain gauges is seriously deficient because of financial limitations or inaccessibility. The proposed framework can also calculate the dynamic strain response by applying the estimated nodal force \overline{\mathit{p}} \left(t\right) to the FEM. Fig.13 and Fig.14 are the reconstructed strain response, which accords well with the “measurement”, i.e., the FEM calculation. The maximum relative errors in of them are 8.7% and 7.1%, respectively. During the harmonic vibration of a bridge, it is more prone to damage. With the reconstructed full-field deflection, inclination, and strain, a structure can be better assessed.

3.4 Influence of noise

As many factors can affect the in-situ measurement, the monitored data inevitably contains noise. There are two noise sources: 1) uncertainty of the input; 2) measurement noise. No matter which type of noise, the noise is reflected in the measurements. In this section, the dynamic load is evenly acting on the whole bridge (different loading scenario with the former sections), the amplitude of each nodal force is 1000/45 N, and the frequency is 15 Hz. 10% Gaussian noise is added to the FEM-calculated results, and the equivalent nodal force is calculated from the polluted data. Fig.15 is the estimated force map, where the measurement noise affects the calculated equivalent nodal force \overline{\mathit{p}} \left(t\right). However, \overline{\mathit{p}} \left(t\right) accurately retains the frequency characteristics of the real input even though the frequency spectrum is polluted (see Fig.16). By applying \overline{\mathit{p}}\left(t\right) to the FEM, the complete responses can be reconstructed. The reconstructed deflection accords well with the FEM-calculated one in Fig.17. The maximum relative error is 12% that is slightly larger than the noise level. Thus, the proposed approach is robust to measurement noise.

4 Damage identification

4.1 The procedure to localize the damage

In the former section, the bridge is undamaged, where the map of the estimated nodal force \overline{\mathit{p}} \left(t\right) can reflect the position, frequency, and amplitude of the real external load. For the damaged structure, \overline{\mathit{p}} \left(t\right) can indicate the damage position only when the position of the real external load is known, which is easy to implement for in-field tests. On the other hand, the map of \mathrm{\Delta }p can localize the damage without knowing the input position. Compared with \overline{\mathit{p}}\left(t\right), \mathrm{\Delta }p is superior to \overline{\mathit{p}}\left(t\right) when there is no prior information about the input. Thus, \mathrm{\Delta }p is principally discussed in the following.

Many factors influence the performance of the damage indicator \mathrm{\Delta }p, including the damage position, damage severity, loading position, and input frequency. This section adopts numerical simulations to discuss their influence. Fig.18 is the damage and loading conditions. A 10 Hz and 1000 N concentrated harmonic force is applied to node No. 18. The elastic modulus of element No.6 is reduced by 30%.

Fig.19(a) and 19(b) are the computed \overline{\mathit{p}}(t) and \mathrm{\Delta }p, respectively. In Fig.19(a), the estimated nodal force \overline{\mathit{p}}(t) deviates from the real load that merely acts on node No. 18. The estimated forces near node No. 18 correspond to the real external load q(t) in Eq. (25) and they are greater than other positions. However, the nodal forces on the left side span (between nodes No. 1 and No. 13) are great too, where there is no load. This is because the stiffness of the left side span decreases, and its deformation increases. \overline{\mathit{p}}(t) is computed by the intact FEM. The equivalent damage load g(t) appears in the damaged area to provide that the dynamic deformation of the intact FEM keeps consistent with the damaged structure. The absolute equivalent damage load at node No.6 is maximal that can localize the damaged element. It should be noted that the loading position is required when adopting \overline{\mathit{p}}(t) to localize the damage.

Fig.19(b) is the map of \mathrm{\Delta }p, the expression of which is in Eq. (27). Since \mathrm{\Delta }p is calculated from the difference of the inclination-estimated and strain-estimated nodal forces, the nodal force corresponding to the external force q(t) is eliminated in \mathrm{\Delta }p. Compared with \overline{\mathit{p}}(t), the amplitude of \mathrm{\Delta }p decreases near node No. 18, which indicates that merely the damage information is retained in \mathrm{\Delta }p while the loading information is eliminated. The amplitude of \mathrm{\Delta }p is maximum at the damage position. Therefore, \mathrm{\Delta }p can localize the damage without knowing the loading position. In addition, the amplitudes of the applied load and \mathrm{\Delta }p at node No. 6 are 1000 and 234 N, respectively, which indicates \mathrm{\Delta }p is sensitive to damage.

The maps of \overline{\mathit{p}} \left(t\right) and \mathrm{\Delta }p can localize the damage with and without the loading position. Therefore, the following sections only discuss \mathrm{\Delta }p. Varied damage scenarios and loading conditions are discussed in detail.

4.2 Input frequency

The structural response relates to the input frequency. If the input frequency is close to the natural frequency, the response is greater and involves fewer orders of modes. With varied input frequencies, the dominant mode is different. This section discusses the influence of the input frequency on the damage indicator \mathrm{\Delta }p.

The damage and loading scenarios are still those in Fig.18. The only difference is the input frequency. As the first two-order natural frequencies of the bridge are 9.46 and 20.2 Hz, this section set the input frequencies as 7, 9, 11, 15, 20, and 25 Hz, respectively. The input frequencies of 9 and 20 Hz are close to the first two natural frequencies. Using the varied input frequencies, the maps of \mathrm{\Delta }p are calculated. Similar to Fig.19, they are three-dimensional. The peaks of \mathrm{\Delta }p are close to a constant by observing Fig.19. For a better comparison between the maps of \mathrm{\Delta }p, the spatial distribution of \mathrm{\Delta }p is extracted when it reaches a peak over time. Besides, \mathrm{\Delta }p is normalized according to the applied external load (1000 N in this section). The normalized damage indicator \mathrm{\Delta }p is plotted in Fig.20.

The normalized \mathrm{\Delta } p shows that the input frequency has a significant influence. When the input frequency is close to the natural frequencies (9 and 20 Hz), the values of \mathrm{\Delta } p near the damage are great and the absolute value of \mathrm{\Delta }p at element No. 6 is maximum. Thus, \mathrm{\Delta } p can localize the damage with high accuracy. In Fig.20(a), the maximal absolute value of \mathrm{\Delta }p is 0.085, not equal to zero for the non-damaged condition, which indicates that \mathrm{\Delta }p is not zero when the structure is undamaged. Thus, a threshold is required to determine the damage. For the bridge in this paper, the threshold for damage is chosen as 0.1. Besides, the distribution of \mathrm{\Delta } p is smooth along the whole bridge without damage.

With damage, sharp peaks appear in the damaged area. When the input frequency deviates from the natural frequencies within ± 2 Hz, \mathrm{\Delta }p can still localize the damage. However, the sensitivity of \mathrm{\Delta }p decreases. When the input frequency is remote from the natural frequencies, \mathrm{\Delta }p cannot localize the damage. Comparison of \mathrm{\Delta }p with 9 and 20 Hz inputs shows that \mathrm{\Delta }p is more sensitive when the input is 20 Hz. This phenomenon relates to the modal shape. In Fig.3, the magnitude of the second modal shape is greater at element No. 6 (the damaged element). Thus, \mathrm{\Delta }p can better localize the damage when the input frequency is close to the second natural frequency.

This section shows that the proposed algorithm can localize the damage with higher sensitivity when the input frequency is close to the structural natural frequency. The position of the damage also influences the performance of \mathrm{\Delta }p.

4.3 Damage severity

Regarding damage identification, the damage indicators are more sensitive to severe damages. The more minor damage is usually difficult to identify, and thus, it is important to discuss the influence of damage severity. The sensor placement and damage position are identical to the former two sections. The only difference in this section is the damage severity. The damage percentage of element No. 6 ranges from 10% to 25% with a step of 5%. The natural frequency of the damaged structure is smaller than the intact structure. In this section, the maximum elastic modulus reduction of element No. 6 is 25%, and the corresponding first natural frequency becomes 9.44 Hz, which is close to the natural frequency of the undamaged FEM, 9.46 Hz. Section 4.2 shows that \mathrm{\Delta }p can localize the damage better when the input frequency is close to the structural natural frequencies. Thus, the input frequency is still 10 Hz, and the amplitude is 1000 N.

Fig.21 is the normalized \mathrm{\Delta }p with different damage percentages. For 10% damage, \mathrm{\Delta }p cannot localize the damage since normalized \mathrm{\Delta }p is smaller than the threshold 0.1. For 15% to 25% damage, the damage is localized by \mathrm{\Delta }p. \mathrm{\Delta }p in the damaged area has sharp peaks. \mathrm{\Delta }p between node No. 14 and No. 21 is non-zero because \mathrm{\Delta }p is calculated from deficient measurement and it cannot provide 100% accuracy. \mathrm{\Delta }p between node No. 21 and No. 45 is close to zero, indicating no damage or no external load exists in this area. The distribution of \mathrm{\Delta }p shows that the damage is in the left side span. At node No. 6 (the damaged element), the absolute value of \mathrm{\Delta }p is maximum, indicating the absolute value of \mathrm{\Delta }p can accurately localize the damage.

4.4 Damage and loading positions

In the previous simulations, the load and damage are at different spans. Different damage and loading positions are combined in this section to verify the proposed approach. When the damage and the dynamic load are at an identical position, \mathrm{\Delta }q cannot indicate the damage position, because the non-zero values of p(t) and g(t) overlap with each other. When the damage and the dynamic load are near each other but not overlapped, \mathrm{\Delta }p can localize the damage. Many numerical simulations show that the method can identify the damages at different positions. This section only illustrates the results of one damage scenario with varied loading positions. There are three conditions, as listed in Tab.1, where the damage conditions are identical and the difference is merely the input.

The sensitivity of \mathrm{\Delta }p is different in Fig.22 under varied loads. For condition (a), the harmonic load is at node No. 8 on the left side span. The normalized \mathrm{\Delta }p is lower than the threshold 0.1 along the whole bridge, and therefore, the damage cannot be detected in this condition. The reason for this failure is that the response of the whole bridge is small under such loading. For conditions (b) and (c), the loads are acting on the middle span at nodes No. 23 and No. 28, respectively. The corresponding bridge responses are greater than condition (a). The normalized \mathrm{\Delta }p in (b) and (c) both indicate the damaged area between nodes No. 17 and No. 18. The sensitivity of (b) is higher than (c) for two reasons. First, the damage is at element No. 17 where the displacement in modal shape is greater for the first mode (see Fig.3) and the first mode is dominant in condition (b) where the applied force is 10 Hz. The second reason is that the applied force for condition (c) is more remote from the damage than (b).

In practice, the damage position is unknown and the optimal input frequency and position cannot be determined. Then, the averaged \mathrm{\Delta }p can be a damage indicator, as shown in Fig.23. The averaged \mathrm{\Delta }p is calculated by averaging the three normalized \mathrm{\Delta }p in Fig.22. With the averaged \mathrm{\Delta }p, the damage can also be localized at node No. 17.

4.5 Influence of noise

Noise exists in both input and measurement. Section 3.4 demonstrates that the framework is robust to noise when reconstructing responses. When localizing the damage, the influence of noise becomes greater. The difference in structural response due to damage is usually very minor, and therefore, it is important to discuss the algorithm with measurement noise.

The loading and damage conditions in this section are identical to those in Tab.1. When detecting the damages, inclination and strain measurements are used. The inclination measurement contains the damage information while the strain does not. Thus, the inclination measurement is less robust to noise. Adding 1% and 10% root mean square Gaussian noise to the inclination and strain measurements, respectively, the corresponding maps of \mathrm{\Delta }pare calculated. It is found that the maps of \mathrm{\Delta }p cannot localize the damage. To eliminate the influence of noise, the numerical simulation is conducted by ten times with identical loading and noise level (dividing measurements of one test into ten parts works as well). The maps of \mathrm{\Delta }p from the ten “tests” are averaged and shown in Fig.24. With a greater number of tests, the influence of noise can be eliminated. Under loading (b) and (c), the normalized \mathrm{\Delta }p indicates that the damage is at node No.17 where the absolute \mathrm{\Delta }p is maximum and above the threshold of damage. Thus, the method can localize damages under the noise environment. The distributions of \mathrm{\Delta }p are similar in Fig.22 and Fig.24 while the magnitude of \mathrm{\Delta }p is greater in Fig.24, because the noise measurement has higher power and the corresponding estimated force is greater.

The three curves (a), (b), (c) in Fig.24 are also added and averaged in Fig.25. The averaged \mathrm{\Delta }p can localize the damage and it is more stable than a single result when the damage position is unknown.

All the numerical simulations demonstrate that the proposed framework can reconstruct the structural response and localize the damage with high accuracy. For practice, the algorithm requires high sampling frequency sensors to provide calculation accuracy. For reconstructing responses, the accuracy of sensors is not strictly required, but the inclinometers with high accuracy are necessary when localizing the damage. The algorithm is not only applicable to bridges, but also to cables, mechanical machines, and other structures. This study merely adopts a continuous beam bridge to validate it.

5 Conclusions

This paper proposes a dynamic stiffness-based framework for estimating the harmonic input, reconstructing structural responses, and localizing damage(s). The work efforts in the paper can be separated into response reconstruction and damage identification, as they are separately investigated. For an undamaged structure, the estimated nodal force can accurately reflect the spatial and frequency information of the real input. Applying the estimated force to the FEM, the full-field responses, including deflection, rotation, and strain, can be well reconstructed under complex inputs. In addition, the reconstructed response is robust to the measurement noise. Therefore, the proposed method can reconstruct the full-field responses of the in-field structures under harmonic load.

During the damage detection stage, the estimated force-based damage indicator \mathrm{\Delta } p can localize the damage since \mathrm{\Delta } p reaches a maximum at the damage position where there is no external input. \mathrm{\Delta } p can better localize the damage when the input frequency is close to the structural natural frequency and the bridge response is great. The threshold 0.1 is selected for the normalized \mathrm{\Delta }p to determine the existence of damage. Exclusively \mathrm{\Delta }p can detect the damage when the loss of stiffness is more than 10%. The required accuracy of inclination is higher than strain, and because it contains damage information, therefore, the inclination is less robust to noise.

Except for its advantages, the proposed framework still has some limitations. Above all, the harmonic loading condition is relatively critical to conduct since artificial loading is necessary. Importantly, the accuracy of damage detection decreases when the loading condition is inappropriate. Though this problem can be settled by multi-loadings, the whole approach becomes laborious. Lastly, the laboratory or in-field test is absent. The framework feasibility requires further validation.

References

[1] Yang C Q, Yang D, He Y, Wu Z S, Xia Y F, Zhang Y F. Moving load identification of small and medium-sized bridges based on distributed optical fiber sensing. International Journal of Structural Stability and Dynamics, 2016, 16( 4): 1640021

[2] Lydon M, Taylor S E, Robinson D, Mufti A, Brien E J O. Recent developments in bridge weigh in motion (B-WIM). Journal of Civil Structural Health Monitoring, 2016, 6( 1): 69–81

[3] Yu Y, Cai C S, Deng L. State-of-the-art review on bridge weigh-in-motion technology. Advances in Structural Engineering, 2016, 19( 9): 1514–1530

[4] Bao T, Babanajad S K, Taylor T, Ansari F. Generalized method and monitoring technique for shear-strain-based bridge weigh-in-motion. Journal of Bridge Engineering, 2016, 21( 1): 04015029

[5] Lansdell A, Song W, Dixon B. Development and testing of a bridge weigh-in-motion method considering nonconstant vehicle speed. Engineering Structures, 2017, 152 : 709–726

[6] Zhao H, Uddin N, O’Brien E J, Shao X, Zhu P. Identification of vehicular axle weights with a bridge weigh-in-motion system considering transverse distribution of wheel loads. Journal of Bridge Engineering, 2014, 19( 3): 04013008

[7] Chan T H T, Law S S, Yung T H, Yuan X R. An interpretive method for moving force identification. Journal of Sound and Vibration, 1999, 219( 3): 503–524

[8] Yu L, Chan T H T. Moving force identification based on the frequency–time domain method. Journal of Sound and Vibration, 2003, 261( 2): 329–349

[9] Law S S, Chan T H, Zeng Q. Moving force identification: A time domain method. Journal of Sound and Vibration, 1997, 201( 1): 1–22

[10] Amiri A K, Bucher C. A procedure for in situ wind load reconstruction from structural response only based on field testing data. Journal of Wind Engineering and Industrial Aerodynamics, 2017, 167 : 75–86

[11] Kazemi Amiri A, Bucher C. Derivation of a new parametric impulse response matrix utilized for nodal wind load identification by response measurement. Journal of Sound and Vibration, 2015, 344 : 101–113

[12] Law S S, Bu J Q, Zhu X Q. Time-varying wind load identification from structural responses. Engineering Structures, 2005, 27( 10): 1586–1598

[13] Hwang J S, Kareem A, Kim H. Wind load identification using wind tunnel test data by inverse analysis. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99( 1): 18–26

[14] Zhi L, Li Q S, Fang M, Yi J. Identification of wind loads on supertall buildings using Kalman filtering-based inverse method. Journal of Structural Engineering, 2017, 143( 4): 06016004

[15] Zhi L, Fang M, Li Q S. Estimation of wind loads on a tall building by an inverse method. Structural Control and Health Monitoring, 2017, 24( 4): e1908

[16] Li Y, Huang H, Zhang W, Sun L. Structural full-field responses reconstruction by the SVD and pseudo-inverse operator-estimated force with two-degree multi-scale models. Engineering Structures, 2021, 249 : 112986

[17] Gillijns S, de Moor B. Unbiased minimum-variance input and state estimation for linear discrete-time systems with direct feedthrough. Automatica, 2007, 43( 5): 934–937

[18] Fang H, de Callafon R A. On the asymptotic stability of minimum-variance unbiased input and state estimation. Automatica, 2012, 48( 12): 3183–3186

[19] Hsieh C S. Extension of unbiased minimum-variance input and state estimation for systems with unknown inputs. Automatica, 2009, 45( 9): 2149–2153

[20] Pan S, Xiao D, Xing S, Law S S, Du P, Li Y. A general extended Kalman filter for simultaneous estimation of system and unknown inputs. Engineering Structures, 2016, 109 : 85–98

[21] Pan S, Su H, Wang H, Chu J. The study of joint input and state estimation with Kalman filtering. Transactions of the Institute of Measurement and Control, 2011, 33( 8): 901–918

[22] Yong S Z, Zhu M, Frazzoli E. A unified filter for simultaneous input and state estimation of linear discrete-time stochastic systems. Automatica, 2016, 63 : 321–329

[23] Hsieh C S, Chen F C. Optimal solution of the two-stage Kalman estimator. IEEE Transactions on Automatic Control, 1995, 44( 1): 194–199

[24] Niu Y, Fritzen C P, Jung H, Buethe I, Ni Y Q, Wang Y W. Online simultaneous reconstruction of wind load and structural responses-theory and application to Canton Tower. Computer-Aided Civil and Infrastructure Engineering, 2015, 30( 8): 666–681

[25] Nord T S, Lourens E M, Øiseth O, Metrikine A. Model-based force and state estimation in experimental ice-induced vibrations by means of Kalman filtering. Cold Regions Science and Technology, 2015, 111 : 13–26

[26] Hsieh C S, Chen F C. Robust two-stage Kalman filters for systems with unknown. IEEE Transactions on Automatic Control, 2000, 45( 12): 2374–2378

[27] Zhang C D, Xu Y L. Structural damage identification via response reconstruction under unknown excitation. Structural Control and Health Monitoring, 2017, 24( 8): e1953

[28] Zhi L, Li Q S, Fang M. Identification of wind loads and estimation of structural responses of super-tall buildings by an inverse method. Computer-Aided Civil and Infrastructure Engineering, 2016, 31( 12): 966–982

[29] Zhang W, Sun L M, Sun S W. Bridge-deflection estimation through inclinometer data considering structural damages. Journal of Bridge Engineering, 2017, 22( 2): 04016117

[30] Sun L M, Zhang W, Nagarajaiah S. Bridge real-time damage identification method using inclination and strain measurements in the presence of temperature variation. Journal of Bridge Engineering, 2019, 24( 2): 04018111

[31] Li Y, Sun L, Zhang W, Nagarajaiah S. Bridge damage detection from the equivalent damage load by multitype measurements. Structural Control and Health Monitoring, 2021, 28( 5): e2709

[32] Bernal D. Load vectors for damage localization. Journal of Engineering Mechanics, 2002, 128( 1): 7–14

[33] Bernal D. Damage localization from the null space of changes in the transfer matrix. AIAA Journal, 2007, 45( 2): 374–381

[34] Wold S, Esbensen K, Geladi P. Principal component analysis. Chemometrics and Intelligent Laboratory Systems, 1987, 2( 1–3): 37–52

[35] JolliffeI. Principal Component Analysis. 2nd ed. New York: Springer, 2010

RIGHTS & PERMISSIONS

Higher Education Press 2022