Sensitivity Analysis Using a Reduced Finite Element Model for Structural Damage Identification

Sensitivity analysis is widely used in engineering fields, such as structural damage identification, model correction, and vibration control. In general, the existing sensitivity calculation formulas are derived from the complete finite element model, which requires a large amount of calculation for large-scale structures. In view of this, a fast sensitivity analysis algorithm based on the reduced finite element model is proposed in this paper. The basic idea of the proposed sensitivity analysis algorithm is to use a model reduction technique to avoid the complex calculation required in solving eigenvalues and eigenvectors by the complete model. Compared with the existing sensitivity calculation formulas, the proposed approach may increase efficiency, with a small loss of accuracy of sensitivity analysis. Using the fast sensitivity analysis, the linear equations for structural damage identification can be established to solve the desired elemental damage parameters. Moreover, a feedback-generalized inverse algorithm is proposed in this work in order to improve the calculation accuracy of damage identification. The core principle of this feedback operation is to reduce the number of unknowns, step by step, according to the generalized inverse solution. Numerical and experimental examples show that the fast sensitivity analysis based on the reduced model can obtain almost the same results as those obtained by the complete model for low eigenvalues and eigenvectors. The feedback-generalized inverse algorithm can effectively overcome the ill-posed problem of the linear equations and obtain accurate results of damage identification under data noise interference. The proposed method may be a very promising tool for sensitivity analysis and damage identification based on the reduced finite element model.


Introduction
The derivatives of the eigenvalue and eigenvector with respect to the system parameters are usually referred to as the sensitivity coefficients [1][2][3]. Sensitivity analysis is often used in the following areas: finite element model (FEM) updating [4,5], structural optimization design [6,7], vibration control [8,9], damage identification [10,11], and so on. Therefore, there has been much interest in efficient computation methods for sensitivity analysis. The classical sensitivity analysis approaches are the modal superposition method [12] and Nelson's method [13]. The key idea of the mode superposition method is to express the eigenvector sensitivity as a linear combination of all modes. This leads to great limitations in the application of mode superposition method, since only a few modes are available for large-scale structures. In view of this, several improved modal superposition algorithms [14][15][16][17] are proposed by using only the partial lower modes. Using these methods, sensitivity analysis for the damping systems [18,19] and complex-valued systems [20,21] are studied by some researchers. Generally, the calculation formula of this kind of method is complex, especially for the higher-order sensitivity analysis. Relatively speaking, Nelson's method has an obvious advantage that only the eigenvector of interest is required. Using the inverse iteration technique, Lin et al. [22,23] further improved the computation efficiency of Nelson's method. Adhikari and Friswell [24] extended Nelson's method for sensitivity analysis of nonviscously damped systems. Wu et al. [25] improved Nelson's method for sensitivity analysis with distinct and repeated eigenvalues by reducing the condition number of the coefficient matrix. Using Nelson's method, Guedria et al. [26] performed the second-order sensitivity analysis for symmetric and asymmetric damped systems with distinct eigenvalues. Wang and Dai [27,28] further studied the eigensensitivities for symmetric and asymmetric damped systems with repeated eigenvalues. Ruiz et al. [29] presented a general framework for sensitivity analysis when tracking specific mode shapes selected beforehand. Lin and Ng [30] used Nelson's method to analyze the eigenpair sensitivities of fractional vibration systems. Recently, Yang and Peng [31] proposed an exact method for calculating the eigenvector sensitivities, which can be seen as an improvement of Nelson's method due to its simple operation.
However, the above sensitivity calculation formulas are generally derived from the complete FEM of the structure, which require a large amount of calculation for large-scale structures. On the other hand, only a few degrees of freedom (DOFs) of the structure can be measured in engineering practice due to the limitation of modal measurement technology. As a result, the reduced FEM [32][33][34][35] is often used in structural mechanics analysis. It is very necessary to study the fast sensitivity analysis approach based on the reduced FEM. To this end, a new sensitivity analysis algorithm is proposed in this paper by using the reduced FEM. Subsequently, structural damage identification based on the proposed sensitivity analysis is studied. A feedback-generalized inverse algorithm is also proposed for solving the damage identification equations more accurately.
The presentation of this work is organized as follows. In Section 2, the common formulas for calculating the eigenvalue and eigenvector sensitivities are briefly reviewed. In Section 3, the new algorithm for fast sensitivity analysis is proposed by using the IRS-z reduced model. The application of this algorithm in structural damage identification is illustrated in Section 4. Numerical and experimental examples are presented in Section 5 to verify the proposed method. Finally, the conclusions of this work are summarized in Section 6.

Review of Sensitivity Calculation Formulas
In this section, the common formulas for calculating the eigenvalue and eigenvector sensitivities are briefly reviewed. Assuming that K and M denote the n × n stiffness and mass matrices of structural FEM, the free vibration modes of the structure can be obtained by solving the following generalized eigenvalue equation [12][13][14][15][16][17][18][19][20][21] as: where λ j and ϕ j denote the jth eigenvalue and eigenvector, j = 1, 2, · · · , n. Clearly, K,M,λ j , and ϕ j are the functions of the physical parameters p i , such as elastic modulus, crosssectional area, and density. By taking a partial derivative of Equation (1) with respect to p i , one has: where ∂λ j ∂p i and ∂ϕ j ∂p i are the first-order eigenvalue and eigenvector sensitivities, respectively. From Equation (1), one can obtain: Using Equations (4) and (2), the eigenvalue sensitivity can be obtained by multiplying Equation (3) by ϕ T j as: For eigenvector sensitivity, Fox and Kapoor [12] proposed the calculation formula based on the modal superposition as: where: It is known that Equation (6) is inefficient, since all the eigenvectors are needed in computation. In contrast, Nelson's method [13] has the great advantage that only the eigenvector of interest is required. Its basic idea is to convert the eigenvector sensitivity into the sum of the general solution and the special solution. Then, the general solution and the special solution are solved, respectively. Based on Nelson's method, Yang and Peng [31] present a direct calculation formula for eigenvector sensitivity as: Note that it is still necessary to further improve the computational efficiency of Equation (10) for large-scale structures. For a structure with more than thousands of DOFs, the computational burden of Equation (10) mainly lies in the solving process of eigenvalue λ j and eigenvector ϕ j from Equation (1). Therefore, it is very meaningful to study a fast calculation method that can reduce the above calculation burden. To this end, a new sensitivity analysis method based on reduced FEM is proposed in the next section.

The Proposed Algorithm for Fast Sensitivity Calculation
As stated before, the existing sensitivity calculation formulas are generally derived from the complete FEM of the structure, which leads to the large amount of calculation. On the other hand, only a few measurement points can be arranged on the structure in engineering practice due to the limitation of modal measurement technology. In view of this, a new algorithm for fast sensitivity analysis is proposed in this work by using the reduced FEM.
Model reduction is a commonly used algorithm to fast estimate some low eigenvalues and eigenvectors of structures by reducing the order of the original structural model to a smaller one. Guyan method [32,33] is the earliest method for FEM reduction. By partitioning the DOFs into the master DOF and slave DOF, Equation (1) where the subscripts "m" and "s" represent the number of the master and slave DOFs, m + s = n. Note that the master DOFs are associated with the measured DOFs in structural dynamic testing. From Equation (13), the Guyan reduced model can be derived as [12,13]: where K r and M r are the reduced stiffness and mass matrices, and T 0 is the DOF transformation matrix. It is clear that the dimensions of K r and M r are both m × m, which just match with the measured DOFs. By solving Equation (14), λ j and ϕ m j can be obtained as the approximate solutions of the low eigenvalues and eigenvectors of the complete FEM. Clearly, the calculation time of solving Equation (14) is far less than that of solving Equation (1) due to m << n. Based on Guyan reduction, the improved reduced system (IRS) reduced model is proposed in [34] by considering the first-order inertia item to improve the computation accuracy as: Furthermore, the IRS-2 reduced model is proposed in [35] by considering the first two inertia items as: As an extension for the IRS reduced method, IRS-z reduced model is proposed in this paper by considering the first z inertia items as follows: Using the IRS-z reduced model, the fast sensitivity calculation formulas can be derived as follows. Using ϕ j = T z ϕ m zj , Equations (5) and (10) can be rewritten as: Equations (32) and (33) are the proposed formulas for quickly computing the sensitivities of eigenvalues and eigenvectors. As will be shown in the numerical examples, the proposed formulas may increase efficiency, with a small loss of accuracy of sensitivity analysis compared with the existing sensitivity calculation formulas.
Finally, it is necessary to discuss how to determine the suitable value of z in the IRS-z reduced model. Note that the inertias in Equation (31) are approximated by Guyan reduced FEM. It is known that the eigenvalues calculated by the Guyan reduced model are generally greater than the exact values. As a result, the accuracy of the IRS-z reduced model first increases, and then decreases with increasing z. Therefore, a criterion to determine the suitable value of z is needed. In this work, the sum of squares of the calculated eigenvalues for each reduced FEM with different z is used as the criterion. That is: The reduced FEM with a minimum value of Θ z is associated with the most suitable value of z.

Application in Structural Damage Identification
The fast sensitivity analysis method proposed above can be used in many engineering fields, such as structural damage identification, model updating, vibration control, etc. In this section, the application of the above sensitivity analysis in damage identification is illustrated as follows.
Generally, structural damage only leads to the reduction in its stiffness matrix. Correspondingly, the stiffness matrix change can be expressed as the sum of the elemental stiffness matrix multiplied by a damage parameter as: where N represents the total number of elements in FEM, K i denotes the i-th elemental stiffness matrix, and x i denotes the corresponding damage parameter (0 ≤ x i ≤ 1). In particular, x i = 0 represents that the i-th element is not damaged and x i = 1 denotes that the i-th element is completely damaged. Using Taylor series expansion, the changes of eigenvalue and eigenvector before and after damage can be approximated as: where ∂λ j ∂α i and ∂ϕ j ∂α i can be computed by the fast sensitivity analysis algorithm proposed above. Combining Equation (36) with (37), the linear equations of damage identification can be obtained as: Commonly, the generalized inverse is used to solve Equation (38) to obtain structural damage parameters as: where x gi is the generalized inverse (GI) solution, A + is the Moore-Penrose generalized inverse [35,36]. For damage identification equations, the condition number of the coefficient matrix A may be very large. For this reason, the calculation results obtained by Equation (42) may be seriously distorted when the data used in the calculation contain random noise [37,38]. In view of this, this paper proposes a feedback-generalized inverse to compute damage parameters more accurately. To this end, a diagonal matrix with small parameters is firstly added to the coefficient matrix of Equation (38) as: The purpose of the above operation is to improve the ill condition of the linear equations. The diagonal elements ε ii in the matrix [ε] are calculated by: where a max is the maximum value of all diagonal elements of A, and cond(A) denotes the condition number of A. From Equation (43), the first solution of x can be obtained as: According to x gi1 , less than 0.05 of the calculated values of damage parameters can be determined to relate with the undamaged elements. From this feedback, Equation (43) can be further reduced by removing the related coefficients associated with undamaged elements as: where x * is the reduced damage parameter vector, and A * is the corresponding coefficient matrix. From Equation (49), the feedback-generalized inverse (FGI) solution of x can be obtained as: Finally, structural damages are evaluated according to x f gi . This feedback process can be repeated until satisfactory results are obtained. In summary, a flowchart, as shown in Figure 1, is given to illustrate the whole technique more clearly.

Numerical Example
The proposed method is verified firstly by the beam structure, as shown in Figure 2.
The basic physical parameters of the structure are: the cross-sectional area Using the reduced model, Table 1 gives the sum of squares of the first seven eigenvalues calculated by different z . For comparison, the calculation results using the complete model are also given in Table 1. It can be seen from Table 1 that the is the smallest. This means that the reduced model with z = 3 is the best condensation model and should be used in the following sensitivity calculation. Tables 2 and 3 present the calculation results of frequency and eigenvector sensitivities based on the complete and IRS-3 reduced model, respectively. It can be seen from Table 2 that the first five eigenvalue sensitivities obtained by the reduced model are very close to those obtained by the complete model. From Table 3, the eigenvector sensitivity obtained by the reduced model is exactly the same as that obtained by the complete model. The computation times by the complete and reduced models are 0.125 and 0.109 s, respectively. These results show that the proposed sensitivity calculation formulas based on the reduced FEM are reasonable and effective. Next, the damage identification method is verified by assuming that the elastic modulus of element 9 in the structure is reduced to 80%. Moreover, 5% random noise is added to the data to simulate the measurement error. Using the first two eigenpair sensitivities, the first calculation result calculated by the generalized inverse is given in Figure 4. It can be found from Figure 4 that unit 9 is the most likely damage unit position. Subsequently, Figure 5 presents the final calculation result based on the feedback. One can see from Figure 5 that unit 9 is the only damage unit, and the calculated damage parameters are very close to the assumed value (0.2). It is clear that the calculated

Numerical Example
The proposed method is verified firstly by the beam structure, as shown in Figure 2. The basic physical parameters of the structure are: the cross-sectional area A = 1.8 × 10 −4 m 2 , the elastic modulus E = 71 GPa, the density ρ = 2.2 × 10 3 Kg/m 3  Using the reduced model, Table 1 gives the sum of squares of the first seven eigenvalues calculated by different z. For comparison, the calculation results using the complete model are also given in Table 1. It can be seen from Table 1 that the Θ z with z = 3 is the smallest. This means that the reduced model with z = 3 is the best condensation model and should be used in the following sensitivity calculation. Tables 2 and 3 present the calculation results of frequency and eigenvector sensitivities based on the complete and IRS-3 reduced model, respectively. It can be seen from Table 2 that the first five eigenvalue sensitivities obtained by the reduced model are very close to those obtained by the complete model. From Table 3, the eigenvector sensitivity obtained by the reduced model is exactly the same as that obtained by the complete model. The computation times by the complete and reduced models are 0.125 and 0.109 s, respectively. These results show that the proposed sensitivity calculation formulas based on the reduced FEM are reasonable and effective. Next, the damage identification method is verified by assuming that the elastic modulus of element 9 in the structure is reduced to 80%. Moreover, 5% random noise is added to the data to simulate the measurement error. Using the first two eigenpair sensitivities, the first calculation result calculated by the generalized inverse is given in Figure 4. It can be found from Figure 4 that unit 9 is the most likely damage unit position. Subsequently, Figure 5 presents the final calculation result based on the feedback. One can see from Figure 5 that unit 9 is the only damage unit, and the calculated damage parameters are very close to the assumed value (0.2). It is clear that the calculated results obtained by feedback have a higher accuracy than the calculated results without feedback.  1 2 3  6  5  4  10 11 12  9  8  7  19 20 21  24  22  16 17 18  15  14  13  25 26 27  30  29  28  31   31  28 29 30  27  26  25  13 14 15  18  17  16  22 23 24  21  20  19  7 8 9  12  11  10  4 5 6  3  2  1 23 32 Figure 2. A beam structure.
(d) Figure 3. Locations of the master DOFs (denoted by the black circle points).

Element Number
Complete model

Experimental Example
The proposed method is verified again by using the experimental data of a CFRP composite cantilever beam from [39]. The schematic diagram and FEM division diagram of the structure are shown in Figure 6. The three-dimensional size of the structure is 200 mm × 20.253 mm × 1.7 mm. The basic physical parameters of the material are: elastic modulus E = 1.33 × 10 5 MPa and density ρ = 1.376 × 10 3 kg/m 3 . The structure is evenly divided into 11 beam elements, and the complete FEM has 66 degrees of freedom. From dynamic testing, the first four natural frequencies of the structure in the undamaged and damaged states are measured as shown in Table 4.

Experimental Example
The proposed method is verified again by using the experimental data of a CFRP composite cantilever beam from [39]. The schematic diagram and FEM division diagram of the structure are shown in Figure 6. The three-dimensional size of the structure is 200mm 20.253mm 1.7mm   . The basic physical parameters of the material are: elastic modulus 5 1.33 10 E   MPa and density 3 1.376 10    kg/m 3 . The structure is evenly divided into 11 beam elements, and the complete FEM has 66 degrees of freedom. From dynamic testing, the first four natural frequencies of the structure in the undamaged and damaged states are measured as shown in Table 4.  The eleven vertical DOFs of nodes 1-11 are chosen as the master DOFs for constructing the reduced FEM. Since only the first four vibration frequencies are measured in modal testing, the damage identification equation is established only by the eigenvalue sensitivities using the reduced FEM. Table 5 gives the sum of squares of the first seven eigenvalues calculated by the reduced FEM with different z . It can be seen from Table 5 that the z  with z = 2 is the smallest. This means that the reduced model with z = 2 can be used in the following sensitivity calculation. Table 6 shows the first four eigenvalue sensitivities obtained by the complete and reduced model. It can be seen from Table 6 that the calculation results obtained by the reduced model are very close to those obtained by the complete model. It was found again that the proposed sensitivity calculation formula based on reduced FEM has high accuracy. In the case where element 1 is damaged, the first calculation result using the measured frequency data before and after damage is given in Figure 7. It can be found from Figure 7 that unit 1 is the most likely damage unit position. Subsequently, Figure 8 presents the final calculation result based on the feedback. It is clear that Figure 8 indicates that unit 1 is the only damage unit, and the calculated damage parameter is very close to the real damage extent (0.3). In the case where elements 5 and 8 are damaged, the first calculation result using the generalized inverse is  The eleven vertical DOFs of nodes 1-11 are chosen as the master DOFs for constructing the reduced FEM. Since only the first four vibration frequencies are measured in modal testing, the damage identification equation is established only by the eigenvalue sensitivities using the reduced FEM. Table 5 gives the sum of squares of the first seven eigenvalues calculated by the reduced FEM with different z. It can be seen from Table 5 that the Θ z with z = 2 is the smallest. This means that the reduced model with z = 2 can be used in the following sensitivity calculation. Table 6 shows the first four eigenvalue sensitivities obtained by the complete and reduced model. It can be seen from Table 6 that the calculation results obtained by the reduced model are very close to those obtained by the complete model. It was found again that the proposed sensitivity calculation formula based on reduced FEM has high accuracy. In the case where element 1 is damaged, the first calculation result using the measured frequency data before and after damage is given in Figure 7. It can be found from Figure 7 that unit 1 is the most likely damage unit position. Subsequently, Figure 8 presents the final calculation result based on the feedback. It is clear that Figure 8 indicates that unit 1 is the only damage unit, and the calculated damage parameter is very close to the real damage extent (0.3). In the case where elements 5 and 8 are damaged, the first calculation result using the generalized inverse is given in Figure 9. Using the feedback process, Figure 10 presents the final calculation results. It is clear that Figure 10 indicates that elements 5 and 8 are the true damage units, and the calculated damage parameters are close to the real damage extents (0.2 and 0.3). This shows that the calculation accuracy after feedback is obviously improved. given in Figure 9. Using the feedback process, Figure 10 presents the final calculation results. It is clear that Figure 10 indicates that elements 5 and 8 are the true damage units, and the calculated damage parameters are close to the real damage extents (0.2 and 0.3). This shows that the calculation accuracy after feedback is obviously improved.

Conclusions
In this work, a fast sensitivity analysis technique is firstly proposed by using the reduced FEM, which can effectively solve the mismatch problem between the number of complete DOFs and that of measured DOFs in practice. It was found that the sensitivity calculation values obtained by the reduced FEM are almost the same as those obtained by the complete FEM for low-order vibration modes. Subsequently, the linear equations of structural damage identification are established by using the proposed sensitivity analysis method with the reduced FEM. A feedback-generalized inverse algorithm is then proposed for solving the damage identification equation more accurately. It was found that the ill-conditioned problem of the linear equations can be overcome to a certain extent and more accurate damage assessment results can be obtained. Numerical and experimental examples verified the feasibility of the proposed method. The proposed method may be a very promising tool for sensitivity analysis and damage identification based on reduced FEM. It is worth noting that the proposed method is based on the assumption of linear elastic vibration. It is also assumed that structural damage only causes the change in stiffness matrix and does not cause the change in mass matrix. Sensitivity analysis of nonlinear vibration and corresponding damage identification methods can be further studied in the future.