Frequency domain damage detection of plate and shell structures by finite element model updating

Abstract

Plate- and shell-like elements are the main constituent of many engineering structures such as ships and airplanes. This paper develops a frequency domain model updating method using power spectral density function and seeks its viability in damage identification of plates and shells. The derivation of the sensitivity equation is exact and measured natural frequencies of few lower modes are used to overcome the drawbacks of incomplete measurement. The method is numerically examined in several damage cases. The robustness of method against measurement and mass modelling error is numerically assessed using Monte-Carlo analysis and numerically simulated error contaminated data. The method proves to be robust against both kind of errors and performs well, both in cases of plate and shell elements.

Keywords:

• Finite-element model updating
• plates and shells
• damage detection
• power spectral density
• measured natural frequency

AMS Subject Classification:

• 34A55
• 65F10
• 65N20

1. Introduction

Plates and shells are widely used for structural applications such as, ships and boats, planes, oil and gas reservoirs, pipelines and bridges. These structures are generally constructed to withstand harsh loading conditions. Most of these structures are intended to carry dangerous and critical fluids. Thus, keeping their integrity in check and monitoring their health is of certain importance. It is promising to detect damages at the earliest possible stage after initiation as it prevents sudden catastrophic structural failure. Among the nondestructive methods for structural inspection, vibration-based methods have known to be capable of diagnosing damage without prior knowledge on its location. In most vibration-based methods, modal data such as natural frequency [1–3], anti-resonance frequencies [4], mode shape [5,6], modal flexibility [7–9], modal strain energy [10–12] as well as Frequency Response Function (FRF) [13–17] have widely been used in model-based and non-model-based damage identification schemes. Modal data (natural frequencies and mode shapes) are inherent structural properties and render great information about structural condition. The accuracy of parameter estimation weakens when the data from fewer lower modes are accessible. In such situations, the frequency domain data are more admissible because they provide data from more frequency points. Availability of data from more frequency points is necessary for robust model updating and solution uniqueness by setting an over-determined equation [18]. Innovative [19,20] or hybrid [21] methods that circumvent the shortcomings of conventional methods are necessary for robust model updating and damage detection using modal data.

The use of vibration-based damage detection methods for plates and shells has recently drawn the attention of researchers [22–25]. Zubaydi et al. [26] used the autocorrelation function to detect damage in stiffened plates. They found the random response metrics suitable for damage detection in plates. Li et al. [11] presented two damage indicator indices namely, bending moment and residual strain mode shape, to locate the damage in plate-like structures. Despite the simplicity of method, some damages at certain locations might not be detected. Wu and Law [27] used incomplete measured mode shape to calculate uniform load surface (ULS) sensitivity for damage detection in plates. They used Chebyshev polynomials of higher order to calculate the sensitivities of ULS. Bayissa and Haritos [28] proposed a non-model-based method for identification and localization of the damage in a plate-like structure using spectral strain energy (SSE). The SSE was derived from moment response power spectral density and curvature power spectral density. They compared their method with modal strain energy (MSE) method proposed by Cornwell et al. [29]. Later Hu et al. [30] extended the use of same method for identifying the damaged regions of cylindrical shells. They numerically tested the damage index for the examples of pipeline structures. Bayissa and Haritos [31] introduced another damage sensitive index based on changes in bending moment power spectral density for plate-like structures. They numerically investigated the method and proved its advantages over changes in natural frequency, modal flexibility and curvature for damage identification. Despite certain advantages of such methods, they need region by region inspection of the whole structure to find the vicinity of damage. Trendafilova and Manoach [32] developed a method based on time series for damage detection and localization in thin plates. Trendafilova et al. [33] applied their method to detect damage in circular plates. Huh et al. [34] introduced a two-dimensional damage index based on vibratory power of plates. The vibratory power of thin plates was calculated based on its response to arbitrary excitations. Wei et al. [20] introduced a two stage damage localization and quantification method based on modal strain energy change ratio (MSECR). They also proposed a novel method to weaken the vicinity effect in localization of damage.

FE model updating using the statistical metrics of structural response such as correlation function and its frequency domain indicator, Power Spectral Density (PSD), have recently been noticed [35,36]. There are a number of advantages for model updating using PSD [35,36]. First, it embraces the data on structural behaviour at wide frequency ranges. Second, it is more sensitive to damage as compared to FRF and modal data. Third, it embraces both auto and cross spectral terms, hence providing more data for model updating. Despite the given advantages of PSD for FE model updating, it has not received enough attention. This is certainly because PSD is a highly sensitive function of structural parameters and the calculation of a reliable linear equation that attributes the changes in PSD to variation of structural parameters is a challenging issue. Hence, some already developed model updating schemes [35–37] resort Taylor series expansion of PSD to develop the sensitivity equation.

This paper presents a model-based damage detection method using PSD function suitable for application to plate- and shell-like structures. A novel sensitivity-based FE model updating method using PSD for damage detection in plates and shells is developed and its challenges are addressed. The sensitivity equation is derived in an exact manner and an approximate formula is used to circumvent the adverse of incomplete measurement. Several damage cases are studied for both kinds of elements and under assumption that only low frequency range data are available. Robustness of method against measurement and mass modelling error is examined by adding random error to simulated PSD of damaged model and mass parameters.

2. Theory

2.1. Development of the sensitivity equation

Power spectral density of structural response subject to multiple stationary excitations is expressed by,(1) $$S_{xx}^{}(\mathit{\omega })=\sum _{i=1}^{Ne}{H}_{}^{\ast }(\mathit{\omega })F_i{H}_{}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(1)

Where, N_e is the number of excitations, F_i is the mapping matrix consists zero and one. S_ffi(ω) is the PSD of ith excitation at ωth frequency. In this equation, S_xx(ω) is the PSD matrix of response in a sense that its diagonal terms are auto spectral density and the non-diagonal ones are cross spectral density terms. The system frequency response function H(ω), is defined as:(2) $$H(\mathit{\omega })={(K-{\mathit{\omega }}^{2}M+i\mathit{\omega }C)}^{-1}=Z{(\mathit{\omega })}^{-1}$$(2)

The transfer function H(ω), is n_d × n_d matrix composed of the system characteristic matrices, mass (M), damping (C) and stiffness (K), and $i=\sqrt{-1}$. In order to develop the sensitivity matrix, the Equation (1(1) $$S_{xx}^{}(\mathit{\omega })=\sum _{i=1}^{Ne}{H}_{}^{\ast }(\mathit{\omega })F_i{H}_{}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(1) ) is rewritten for the intact (denote by I index) and damaged (denote by D index) states of structure as,(3) $$Z_I^{\ast }{S}_{xx}^I(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i{H}_{}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(3) (4) $$Z_D^{\ast }(\mathit{\omega })S_{xx}^D(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i{H}_D^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(4)

Damage will cause δM, δK and δC changes in mass, stiffness and damping matrices, respectively. Hence for the damaged structure, Equation (4(4) $$Z_D^{\ast }(\mathit{\omega })S_{xx}^D(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i{H}_D^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(4) ) is stated in the following format,(5) $$\left((K_I+\mathit{\delta }K)-{\mathit{\omega }}^2(M_I+\mathit{\delta }M)-i\mathit{\omega }(C_I+\mathit{\delta }C)\right)(S_{xx}^D(\mathit{\omega }))=\sum _{i=1}^{Ne}{F}_i{H}_D^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(5)

And by further expansion of Equation (4(4) $$Z_D^{\ast }(\mathit{\omega })S_{xx}^D(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i{H}_D^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(4) ), it will transform as,(6) $$(K_I-{\mathit{\omega }}^2{M}_I-i\mathit{\omega }{C}_I)S_{xx}^D(\mathit{\omega })+(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M-i\mathit{\omega }\mathit{\delta }C)S_{xx}^D(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i{H}_D^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(6)

Upon the subtraction of Equation (2(2) $$H(\mathit{\omega })={(K-{\mathit{\omega }}^{2}M+i\mathit{\omega }C)}^{-1}=Z{(\mathit{\omega })}^{-1}$$(2) ) from (5), it yields,(7) $$Z_I^{\ast }\mathit{\delta }{S}_{xx}(\mathit{\omega })=\sum _{i=1}^{Ne}{F}_i\mathit{\delta }{H}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })-(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M-i\mathit{\omega }\mathit{\delta }C)S_{xx}^D(\mathit{\omega })$$(7)

Finally by substituting Equation (1(1) $$S_{xx}^{}(\mathit{\omega })=\sum _{i=1}^{Ne}{H}_{}^{\ast }(\mathit{\omega })F_i{H}_{}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })$$(1) ) for damaged state, some mathematical manipulation δZ = δK − ω²δM + iωδC the final expression for the change in the power spectral density is derived as,(8) $$\mathit{\delta }{S}_{xx}(\mathit{\omega })=\sum _{i=1}^{Ne}{H}_I^{\ast }(\mathit{\omega })F_i\mathit{\delta }{H}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })-H_I^{\ast }(\mathit{\omega })(\mathit{\delta }{Z}^{\ast })S_{xx}^D(\mathit{\omega })$$(8)

In this equation, δH(ω) is the change in the reacceptance matrix as expressed in the Appendix.(9) $$\mathit{\delta }H(\mathit{\omega })=-H_D(\mathit{\omega })(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M+i\mathit{\omega }\mathit{\delta }C)H_I(\mathit{\omega })$$(9)

By substitution of the transpose of Equation (9(9) $$\mathit{\delta }H(\mathit{\omega })=-H_D(\mathit{\omega })(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M+i\mathit{\omega }\mathit{\delta }C)H_I(\mathit{\omega })$$(9) ) in the (8), the change in the spectral density matrix is expressed as(10) $$\mathit{\delta }{S}_{xx}=-H_I^{\ast }(\mathit{\omega })\left(\sum _{i=1}^{Ne}(F_i{H}_I^T(\mathit{\omega })(\mathit{\delta }Z)+(\mathit{\delta }{Z}^{\ast })H_D^{\ast }(\mathit{\omega })F_i)S_{ffi}(\mathit{\omega })\right)H_D^T(\mathit{\omega })$$(10)

The change of PSD with respect to change in stiffness ($\mathit{\delta }{S}_{xx}^s$), mass ($\mathit{\delta }{S}_{xx}^m$) and damping ($\mathit{\delta }{S}_{xx}^c$) matrix is expressed using the following equations.(11) $$\mathit{\delta }{S}_{xx}^s=-H_I^{\ast }(\mathit{\omega })\left(\sum _{i=1}^{Ne}(F_i{H}_I^T(\mathit{\omega })(\mathit{\delta }K)+(\mathit{\delta }K)H_D^{\ast }(\mathit{\omega })F_i)S_{ffi}(\mathit{\omega })\right)H_D^T(\mathit{\omega })$$(11) (12) $$\mathit{\delta }{S}_{xx}^m={\mathit{\omega }}^2{H}_I^{\ast }(\mathit{\omega })\left(\sum _{i=1}^{Ne}(F_i{H}_I^T(\mathit{\omega })(\mathit{\delta }M)+(\mathit{\delta }M)H_D^{\ast }(\mathit{\omega })F_i)S_{ffi}(\mathit{\omega })\right)H_D^T(\mathit{\omega })$$(12) 13) $$\mathit{\delta }{S}_{xx}^c=i\mathit{\omega }{H}_I^{\ast }(\mathit{\omega })\left(\sum _{i=1}^{Ne}(-F_i{H}_I^T(\mathit{\omega })(\mathit{\delta }C)+(\mathit{\delta }C)H_D^{\ast }(\mathit{\omega })F_i)S_{ffi}(\mathit{\omega })\right)H_D^T(\mathit{\omega })$$13)

The above sensitivity relation is solvable using linear least-squares error minimization through an iterative model updating scheme. Moreover, in model updating process, elemental connection and boundary condition type is of certain importance [38]. In this study, it is assumed that the connection between elements and boundary condition is not affected by the damage. In an iterative model updating scheme, change in stiffness, mass and damping matrices are stated as,(14) $$\mathit{\delta }{K}_{}=\sum _{n=1}^{ne}{K}_n{d}_n^k$$(14) (15) $$\mathit{\delta }M=\sum _{n=1}^{ne}{M}_n{d}_n^m$$(15) (16) $$\mathit{\delta }C=\sum _{n=1}^{ne}{C}_n{d}_n^c$$(16)

k_n, m_n, c_n represent the element level stiffness, mass and damping matrices assembled to the global matrix at the nth element DOFs and $d_n^{k,m,c}$ indicates the normalized damage ratio attributed to the nth elements. The final format of the model updating equation will be,(17) $$\mathit{\delta }{S}_{xx}=S_k{d}^k+S_m{d}^m+S_c{d}^c$$(17)

where, S_k, S_m and S_c are the sensitivity of PSD to stiffness, mass and damping parameters.

The developed sensitivity Equation (17(17) $$\mathit{\delta }{S}_{xx}=S_k{d}^k+S_m{d}^m+S_c{d}^c$$(17) ) should be solved by least-squares method for unknown changes in the structural parameters. It is however, an ill-posed problem. It is worth noting, vibration-based damage detection methods using finite-element model updating engender ill-posed inverse problems, unless special care is given to the well-posedness of method [39]. An ill-posed set of linear equations is solvable if proper side constraint is defined to the problem [40]. In this regard, a set of unprejudiced side constraints are defined to the problem as $-1\le d_n^{k,m,c}\le 1$. Hence, the inverse problem for damage detection and model updating turns to be,(18) $$\begin{array}{c}\text{Find}:d^T={\{d_1^{k,m,c},d_2^{k,m,c},...,d_n^{k,m,c},...,d_{ne}^{k,m,c}\}}^T\hfill \\ \text{Minimise}:W(d)=\left[{\left(\mathit{\delta }{S}_{xx}-S_{k,m,c}{d}^{k,m,c}\right)}^T\left(\mathit{\delta }{S}_{xx}-S_{k,m,c}{d}^{k,m,c}\right)\right]\hfill \\ \text{Subject to:}-1\le d_n^{k,m,c}\le 1.\hfill \end{array}$$(18)

The solution to the above problem is reached by MATLAB’s optimization toolbox through constrained linear least-squares method.

In solving a model updating problem using least-square method, some of the equations with largest coefficients may dominate the solution. In a model updating scheme using frequency domain data, it is possible to weight each equation by the, inverse of the standard deviation of data, the norm of the sensitivity equation, or by the inverse of the excitation frequency (ω⁻¹) [13]. In this paper, each row of the sensitivity is divided by its second norm. Upon implementation of the weighting approach the discrepancy between the impacts of equations on the solution vanishes as the norm of all rows becomes unity. The application of this weighting approach also improves observability of less observable design variables in the model updating process.

2.2. Model updating issues

The success of a model updating scheme based on the frequency domain data depends on the excitation frequency ranges, the number and location of measured and excited DOFs, and the accuracy of the sensitivity equation. The exact sensitivity presented in Equation (8(8) $$\mathit{\delta }{S}_{xx}(\mathit{\omega })=\sum _{i=1}^{Ne}{H}_I^{\ast }(\mathit{\omega })F_i\mathit{\delta }{H}^T(\mathit{\omega })S_{ffi}(\mathit{\omega })-H_I^{\ast }(\mathit{\omega })(\mathit{\delta }{Z}^{\ast })S_{xx}^D(\mathit{\omega })$$(8) ) needs the spectral density of damaged structure. Evaluating $S_{xx}^D(\mathit{\omega })$ necessitates complete measurement of FRF of the damaged structure. In order to avoid the need for measured FRF of damaged structure it is approximated as,(19) $$H_D^{\text{approx.}}(\mathit{\omega })=\sum _{j=1}^m\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jD}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jD}{\mathrm{\Omega }}_{jD}\mathit{\omega }}+\sum _{j=m+1}^{nd}\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jI}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jI}{\mathrm{\Omega }}_{jI}\mathit{\omega }}$$(19)

The first term consists of measured natural frequencies Ω_jD, modal damping loss factors ξ_jD of the damaged structure at a m lower modes and the corresponding mode shapes of the analytical model ϕ_jI and Ω_jI. The second term consists of the analytical model natural frequencies, damping loss factors and mode shapes of analytical model corresponding to unmeasured modes. The second term is added in order to decrease the adverse impacts of measurement incompleteness in higher frequency ranges. It also helps sensitivity equation to be more precise during the model updating process by updating the modal characteristics of higher modes.

The accuracy of the approximated $S_{xx}^D(\mathit{\omega })$, calculated by Equation (19(19) $$H_D^{\text{approx.}}(\mathit{\omega })=\sum _{j=1}^m\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jD}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jD}{\mathrm{\Omega }}_{jD}\mathit{\omega }}+\sum _{j=m+1}^{nd}\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jI}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jI}{\mathrm{\Omega }}_{jI}\mathit{\omega }}$$(19) ) is illustrated by comparing it with its exact counterpart for a typical damage case in Figure 1. The plot presented in Figure 1 proves the ability of Equation (19(19) $$H_D^{\text{approx.}}(\mathit{\omega })=\sum _{j=1}^m\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jD}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jD}{\mathrm{\Omega }}_{jD}\mathit{\omega }}+\sum _{j=m+1}^{nd}\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jI}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jI}{\mathrm{\Omega }}_{jI}\mathit{\omega }}$$(19) ) in presenting a well-approximate power spectral density of damaged model. It is also evident that the approximation is more precise near measured resonances of damaged model. Due to incorporating modal data of damaged model in the first summation, the approximate PSD is well clear from the PSD of intact model up to the mth natural frequency. Since the approximate formula is controlled by the denominator, so the accuracy of Equation (19(19) $$H_D^{\text{approx.}}(\mathit{\omega })=\sum _{j=1}^m\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jD}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jD}{\mathrm{\Omega }}_{jD}\mathit{\omega }}+\sum _{j=m+1}^{nd}\frac{{\mathit{\varphi }}_{jI}{\mathit{\varphi }}_{jI}^T}{{\mathrm{\Omega }}_{jI}^2-{\mathit{\omega }}^2+2i{\mathit{\xi }}_{jI}{\mathrm{\Omega }}_{jI}\mathit{\omega }}$$(19) ) reduces by marching from resonances to anti-resonances.

Figure 1. Approximate damaged PSD and exact damaged PSD.

The appropriate selection of the excitation frequency ranges will result in the smoothness of the optimization process, convergence rate of model updating scheme, reducing the effect of noise on PSD residual. In this study, the frequency ranges in vicinity of resonances are adopted as excitation frequency range. In order to avoid unsmooth behaviour of optimization process, excitation frequencies between natural frequencies of the intact and damaged structure are excluded [14,17,41]. Between natural frequencies of intact and damaged model, the residual between intact and damaged PSD reduces and then increases. Therefore, a smooth reduction of the residual is not possible in successive iterations [41]. At the selected excitation frequency, high accuracy of sensitivity equation and significant difference between the PSD of damaged and intact model decrease negative effect of noise on δS_xx. In other word, the frequency ranges where the measurement noise overshadows the residual should be avoided. Furthermore, the frequency ranges with notable effect of damping, on the response of intact model, like at the resonances and anti-resonances should be avoided in order to remove the need for an accurate damping model [42]. The frequency range of interest will be illustrated in next section.

3. Numerical investigations

3.1. Finite element formulation

According to Bathe [44], the displacement field within a shell element is defined by the following displacement field functions.(20) $$u(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })u_k+\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k(-{\mathit{\alpha }}_i{V}_{jk}^x+{\mathit{\alpha }}_j{V}_{ik}^x)$$(20) (21) $$v(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })v_k+\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k(-{\mathit{\alpha }}_i{V}_{jk}^y+{\mathit{\alpha }}_j{V}_{ik}^y)$$(21) (22) $$w(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })w_k+\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k(-{\mathit{\alpha }}_i{V}_{jk}^z+{\mathit{\alpha }}_j{V}_{ik}^z)$$(22)

Figure 2 shows a schematic of shell element. The given shape function is suitable for modelling thin to moderately thick shells.

Figure 2. Schematic of shell element.

In these equations, nn is the number of nodes, ξ, η, μ are the natural coordinates, u, v, w are the displacement in x, y, z directions and u_k, v_k, w_k are the nodal displacements. Moreover, N_k(ξ, η) is the isoparametric shape function adopted based on the number of nodes within element, t_k is the thickness of element at kth node and $V_{nk}^{}$ is the perpendicular unit vector to the shell mid-surface at nodes. $V_{ik}^{}$ and $V_{jk}^{}$ are two unit vectors perpendicular to $V_{nk}^{}$. Moreover, α_i and α_j are the rotation about these two vectors, respectively.

3.1.1. Stiffness and mass matrix of plate element

In the case of a plate element the unit vectors are V_i = (1, 0, 0), V_j = (0, 1, 0) and V_n = (0, 0, 1). In other words, they are coincident to X–Y–Z axes. Moreover, α_i = ψ_y and α_j = ψ_x, where ψ_x and ψ_y are the rotations about y and x axis. In case of plate elements in-plane displacements are inadmissible and therefore u_k = v_k = 0, hence making the plate element displacement field become,(23) $$u(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k({\mathit{\psi }}_x^k)$$(23) (24) $$v(\mathit{\xi },\mathit{\eta },\mathit{\mu })=-\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k({\mathit{\psi }}_y^k)$$(24) (25) $$w(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })w_k$$(25)

Equations ((23(23) $$u(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\frac{\mathit{\mu }}{2}\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })t_k({\mathit{\psi }}_x^k)$$(23) )–(25(25) $$w(\mathit{\xi },\mathit{\eta },\mathit{\mu })=\sum _{k=1}^{nn}{N}_k(\mathit{\xi },\mathit{\eta })w_k$$(25) )) are the well-known Mindlin–Reissner plate elements shape functions. So, plate element has one transitional (w_k) and two rotational (${\mathit{\psi }}_x^k$, ${\mathit{\psi }}_y^k$) on each node. A total of twelve DOFs are attributed to each element. For plate elements, the normal (ε_xx, ε_yy) and shearing (γ_xy, γ_yz, γ_xz) strains are the admissible components. Normal strain components as well as in-plane shear components (ε_xx, ε_yy, γ_xy) constitute the bending portion of the stiffness matrix and (γ_yz, γ_xz) form the shearing portion of the stiffness matrix. So, the stiffness matrix is represented in natural coordinate system as,(26) $$K_e^p=\stackrel{1}{\underset{-1}{\int \int \int }{B}_b^T}{D}_b{B}_b\left|J\right|d\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }+\stackrel{1}{\mathit{\kappa }\underset{-1}{\int \int \int }{B}_s^T}{D}_s{B}_s\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }d\mathit{\mu }$$(26) (27) $$D_b=\frac{E}{(1-{\mathit{\upsilon }}^2)}\left[\begin{array}{ccc}1& \mathit{\upsilon }& 0\\ \mathit{\upsilon }& 1& 0\\ 0& 0& 0.5(1-\mathit{\upsilon })\end{array}\right]$$(27) (28) $$D_s=\frac{E}{2(1+\mathit{\upsilon })}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$(28) (29) $$M_e=\stackrel{1}{\underset{-1}{\int \int \int }{N}^T}MN\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }$$(29) (30) $$M=\mathit{\rho }\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$(30)

In above equations, E is modulus of elasticity, υ is the poison ratio, κ is the shear correction factor depending on the boundary condition and $\left|J\right|$ is the determinant of Jacobin matrix. B_b is the (3*12) bending strain mapping matrix and B_s is the shear strain mapping matrix. N is displacement field to DOF mapping matrix. In case of plate element, the integration with respect to μ is held explicitly and the integration with respect to the other variables is held numerically using Gaussian quadrature integration scheme. In order to avoid the problem of shear locking; the shear portion of stiffness matrix is integrated using less integration points. It is worth noting here, the thickness is considered as uniform in all the nodes in this study.

3.1.2. Stiffness and Mass matrix of shell element

In case of shell element, the element behaviour is influenced by membrane effect in addition to bending and shearing effects. For a shell element, there are 5 degrees of freedom at each node, three transitional (u_k, v_k and w_k) and two rotational (α_i and α_j). Hence, a 4-node shell element has 20 DOFs. The stiffness matrix of a shell element is introduced as,(31) $$K_e^{sh}=\stackrel{1}{\underset{-1}{\int \int \int }{B}_b^T}{D}_b{B}_b\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }+\stackrel{1}{\underset{-1}{\int \int \int }{B}_m^T}{D}_m{B}_m\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }+\stackrel{1}{\underset{-1}{\int \int \int }{B}_s^T}{D}_s{B}_s\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }$$(31) (32) $$M_e^{sh}=\stackrel{1}{\underset{-1}{\int \int \int }{N}^T}MN\left|J\right|\text{d}\mathit{\xi }\text{d}\mathit{\eta }\text{d}\mathit{\mu }$$(32) (33) $$M=\mathit{\rho }\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$(33)

In this case the integral could not be held explicitly on μ, so triple integral is held numerically to calculate the three above integrals. The first integral is bending, second integral is membrane and the third integral is shear portion of $K_e^{sh}$ matrix. The bending portion is integrated with 3 * 3 * 3 Gaussian integration scheme and the other integrals are held 2 * 2 * 2 to avoid membrane and shear locking. Further detail about the factors introduced above could be found in Refs. [43–45].

3.2. Measurement and excitation condition

Each excitation at each frequency point will lead to the auto-spectral density (ASD) of the response at each sensor location and cross spectral density (CSD) of the response between the sensors. The data collected by each sensor will lead to Auto-spectrum with itself and Cross-spectrum with the data collected by the other sensors in frequency domain. So, considering n_m sensors, the size of PSD matrix is n_m × n_m. The data on the main diagonal are the ASDs of each sensor with itself. The out of diagonal element of the matrix are the CSDs between the sensors. Since it is a symmetric matrix, the data from the lower triangle are just accounted as,(34) $$S_{xx}^D=\left[\begin{array}{ccccc}{A}_{1,1}& & & & \\ C_{1,2}& A_{2,2}& & Sym.& \\ ⋮& ⋮& \ddots & & \\ C_{1,nm-1}& C_{2,nm-1}& \cdots & A_{nm-1,nm-1}& \\ C_{1,nm}& C_{2,nm}& \cdots & C_{nm-1,nm}& A_{nm,nm}\end{array}\right].$$(34)

Therefore, with n_m measured responses, the total number of data points generated at each frequency because of an excitation will be(35) $$n_{\text{data}}=\frac{n_m(n_m+1)}{2}$$(35)

Regarding the above-mentioned points, the methods leads to n_f × n_e × n_data equations to be solved by least-square method. Moreover n_f and n_e are the number of frequency points and excitation points, respectively. This is one of the major advantages of using PSD for FE model updating rather than FRF. In other words, FRF methods give n_m equations for a given excitation at a certain frequency point.

3.3. Damage detection results

3.3.1. Plate example

The method is applied numerically to the model of a flat plate with four simply supported edges (SSSS) as depicted in Figure 3.

Figure 3. Finite element model of plate.

The plate is 1 m × 1 m and its finite-element model consists of 25 elements (5 divisions both in x and y direction). Numbers of nodes with active transitional DOFs are shown in Figure 3. The material of the model is considered to be steel with Young modulus of 200 Gpa, Poisson ratio of 0.3 and mass density of 7800 kgm⁻³. The uniform thickness of 10 mm is considered for the plate. It is assumed that the structure is excited at nodes 6, 7, 10, 11 and structural responses are measured at 8 measurement location. These points are selected on the transitional DOFs of 2, 3, 5, 9, 8, 12, 14, and 15 nodes. No excitation and measurement is considered for the rotational DOFs.

In order to assess the robustness of the proposed method against measurement errors, 15% random error is added to the simulated PSD data of the damaged model. The simulated PSD of the damaged model is contaminated with noise by adding a certain per cent of normally distributed random error multiplied by the PSD of damaged model as,(36) $$S_{xx}^{Dn}=S_{xx}^D+Ep.\text{Noise}.S_{xx}^D$$(36)

where, Ep is the error percent. Noise is a zero mean normally distributed random error vector. $S_{xx}^D$ is the simulated PSD of damaged model. $S_{xx}^{Dn}$ is the simulated PSD of damaged model contaminated with noise. An example of simulated noisy PSD of damaged model is given in Figure 4. The frequency ranges used for model updating according to Section 2.2, are also marked on this figure. In such ranges, the model updating process is least affected by the measurement noise and the adopted model for damping of intact model.

Figure 4. Simulated noisy PSD of the damaged model and frequency range of interest.

Fifty sets of the simulated error contaminated PSD data are used to conduct Monte-Carlo simulation. The COV of the predicted parameters is also used to investigate the robustness of the method. The lower COVs of the estimated damage parameters proves the robustness of the method against measurement error. In order to investigate the performance of the proposed method for plate elements, several damage cases with different severities are investigated. The studied damage cases are introduced in the Table 1.

Table 1. Damage cases studied for plate example.

Case number	Stiffness reduction
Case P-1	Damaged elements	12	13	17	18
	Damage per cent (%)	25	25	25	25
Case P-2	Damaged elements	8	14	17	23
	Damage per cent (%)	40	40	50	50
Case P-3	Damaged elements	5	7	9	13	17	23
	Damage per cent (%)	10	20	5	10	20	15

In this study, in order to probe the method in practical condition, lower frequency range is considered for model updating. In this regard, the frequency range of interest is selected around the first six natural frequencies. The natural frequencies of intact and damaged plate models are summarized in Table 2. Moreover, the excitation frequency ranges used for model updating are given in the last column of this table.

Table 2. Natural frequencies and excitation frequency ranges (Hz)for considered damage case.

Case	Natural frequencies						Excitation frequency ranges
	1st	2nd	3rd	4th	5th	6th
Intact	50.31	142.58	142.58	231.093	382.22	382.22	–
Case P-1	48.28	137.73	139.38	226.31	365.94	365.53	43–47; 52–56
							144–148; 221–225
							232–236; 384–388
Case P-2	48.80	131.65	135.69	222.66	358.34	370.61	43–47; 52–56
							144–148; 232–238
							383–387
Case P-3	49.27	138.35	140.52	225.13	374.38	378	44–48; 51.5–57
							144–148; 220–224
							232–236; 383–388

The model updating results of three damage scenarios are presented in Figures 5–7. Parts (a) of these figures present the mean of the parameter estimation results, while part (b) shows the COV of the estimated parameters.

Figure 5. P-1 damage case.

Figure 6. P-2 damage case.

Figure 7. P-3 damage case.

Close values of the predicted and actual damage parameters; demonstrate that the method is capable of localization and quantification of damage irrelevant to the location of damaged elements. The norm of the residual ||δS_xx − S_k . d^k||₂ for the studied damage cases is plotted in Figure 8. As this figure shows, the method exhibits a fast convergence rate. It is worth noting that due to presence of noise the residual cannot reach to zero in the model updating process.

Figure 8. The residual of PSD changes at each iteration for the studied damage cases.

3.3.1.1. Mass modelling errors

Mass matrix is the discrete representation of the distributed structural mass [45]. There is a possibility that the developed mass matrix is not the perfect representation of the structural inertia. On, the other hand, quality of mass parameters used for modelling the inertia properties of the structure is critical to the accuracy of the method. In order to simulate inaccuracy in mass matrix 10% random error in mass parameters of the model is considered and P-2 and P-3 damage case are reinvestigated. Predicted damage parameters are given in Figures 9 and 10. By comparing the results in Figures 9 and 10 with Figures 5 and 7, it is understood that incorporation of mass parameter errors will not affect the robustness of the method.

Figure 9. P-1 damage case with 10% mass modelling error.

Figure 10. P-2 damage case with 10% mass modelling error.

3.3.2. Shell example

Shells are widely used in the construction of various engineering structures. So, in this section the performance of method is investigated in such cases. A half cylinder with 3.2 m length and 1 m radius as depicted in Figure 11 is considered for numerical investigation.

Figure 11. Geometry of cylindrical shell example.

All edges of the model are considered to be clamped. The material properties of steel same as plate example are considered here with the thickness of 10 mm. According to Figure 12, model is discretized by 32, 4-node elements. With the considered boundary conditions it has 21 active nodes. Having 5 DOFs at each node, the model has 105 active DOFs. The transitional DOFs normal to the element surface on nodes 8, 10, 12 and14 are considered as excitation points. Eighteen considered measurement locations are selected among the transitional DOFs at three directions. Transitions of nodes 1, 9, 13 and 2 in X direction, Transitions of nodes 5, 11 and 17 in Y direction and Transitions of nodes 1, 3, 5, 7, 9, 11, 13, 15, 17, 19 and 21 in Z direction are the selected measurement locations. In order to investigate the performance of the proposed method for shell elements, several damage cases with different severities are investigated. The studied damage cases are introduced in the Table 3.

Figure 12. Finite-element model of shell.

Table 3. Damage cases studied for shell example.

Case	Stiffness reduction
Case Sh-1	Damaged elements	7	10	15	22	25	31
	Damage per cent (%)	10	15	10	10	10	15
Case Sh-2	Damaged elements	12	18	22
	Damage per cent (%)	60	40	50
Case Sh-3	Damaged elements	1	7	19	22	30
	Damage per cent (%)	25	20	20	30	15

The first six lowest natural frequencies are assumed as measured and excitation frequency ranges for model updating are considered around them. The natural frequencies of the shell model in intact and damaged cases, as well as excitation frequency ranges for model updating are listed in Table 4.

Table 4. Natural frequencies and excitation frequency ranges (Hz) for considered damage cases.

Case	Natural frequencies						Excitation frequency ranges
	1st	2nd	3rd	4th	5th	6th
Intact	445.069	573.27	659.37	672.56	787.82	892.56	–
Case Sh-1	442.1647	566.37	650.61	663.68	778.56	886.77	437–441; 446–450
							561–565; 575–579
							661–664; 773–777
							789–793;894 898
Case Sh-2	443.18	561.98	637.42	656.60	774.03	857.13	427–432; 446–450
							556–561; 575–579
							789–793
Case Sh-3	440.83	560.53	648.023	659.61	771.27	880.52	436–440; 446–450
							555–559; 575–579
							661–665; 766–770
							790–794

The model updating results are presented in Figure 13 through 15. The results show robust localization and quantification of structural damages in case of shell structure.

Figure 13. Sh-1 damage case.

Figure 14. Sh-2 damage case.

Figure 15. Sh-3 damage case.

3.3.2.1. Mass modelling errors

Assuming 10% random error in mass parameter of the model, Sh-2 and Sh-3 damage case are reinvestigated. Presented results in Figures 16 and 17 show that the method can withstand the adverse effects of the above-mentioned errors.

Figure 16. Sh-2 damage case with 10% mass modelling error.

Figure 17. Sh-3 damage case with 10% mass modelling error.

3.3.3. Zone damage and distributed damage

In this section, three cases of zone damage (damage in a portion of elements) and deterioration damage are investigated for plate and shell model. Occurrence of such damage cases is very probable in the highly corrosive environments such as marine and offshore. These damage cases are produced by adding up to 10% random changes to the stiffness parameters. For better illustration and readability of results, linear bar chart representation is adopted here. The method proves to be successful in such cases as presented in Figures 18–20.

Figure 18. Deterioration damage case of plate.

Figure 19. Deterioration damage case of shell.

Figure 20. Zone damage case of shell.

4. Effect of frequency ranges on the observability of elements

By using frequency domain representation of structural responses by PSD or FRF adequate data are available to perform model updating in wide frequency ranges. This is an advantage in case of shell- and plate-like structures. Generally it is indicated that higher frequency ranges data are more sensitive to structural damages and are recommended for model updating. However, for flexural elements like plates and shells more insight must be attended to this issue. Because, By selecting the structural response at higher frequency ranges, the number of low amplitude vibrating DOFs will increase, which leads to low sensitivity of the selected metric to change of structural parameters. In order to demonstrate this, the studied damage cases are reconsidered here.

Among the frequency ranges introduced in Table 2 for each damage case, four frequency ranges are considered to as given in the Table 5 to investigate the effect of frequency ranges. For the plate model with SSSS boundary condition as it is given in Table 2, second and third natural frequencies, as well as fifth and sixth natural frequencies are coincident in the same frequency. For each damage case, the first frequency range is in the vicinity of first natural frequency. Second range is in the vicinity of second and third natural frequencies. Third range is around the forth natural frequency and forth range is near fifth and sixth natural frequencies. The mode shape corresponding to the selected excitation frequency ranges and of norm of sensitivity equation columns are plotted in the Figures 21–23. Norm of the sensitivity matrix’s columns is evaluated as a measure of element observability.

Table 5. Study frequency ranges for damage cases.

Frequency range	Case P-1	Case P-2	Case P-3	Distributed damage
Range 1	52–56	52–56	51.5–57	52–56
Range 2	144–148	144–148	144–148	146–150
Range 3	232–236	232–238	232–236	226–230
Range 4	384–388	383–387	383–388	384–388

Figure 21. The effect of frequency range on the observability of elements.

Figure 22. The effect of frequency range on the observability of elements.

Figure 23. The effect of frequency range on the observability of elements.

From these figures, it is discerned that, design variables relevant to elements at higher curvature, exhibit more observability as compared to other elements. Furthermore, elements located on the null point of the corresponding mode shape, represent the minimum observability in related excitation frequency ranges. The other elements are moderately observable. In brief as the number of null points in a mode shape increases, the number of low observable elements in the frequency ranges of its vicinity will increase too. In order to compare the observability of the elements regarding the selected frequency range, the norms of the sensitivities matrix before and after weighting are plotted simultaneously in a bar chart plot of Figure 24 for P-1 damage case. This comparison shows that as marching from low frequency ranges toward the high frequency ranges, the sensitivity of the metric to change of structural parameters drastically reduces. So, it is advantageous to consider low frequency ranges for model updating of plates.

Figure 24. Comparison of the norm of sensitivity at four frequency ranges.

The effect of weighting approach on the observability of design variables in the study frequency ranges is the other matter discussed here. Norm of the weighted sensitivity is also plotted in Figure 24(b). The weighting approach does not affect the sensitivity of an element with respect to other elements in a frequency range. However, the significant difference between the sensitivity norms of elements before weighting is eliminated after weighting (Figure 24(a) and (b)). Thus, a balanced effect of the frequency ranges on the results is expected by application of weighting. The same study is done in the case of shell model and the same observations are made. Yet, the results are not presented to avoid repeated discussions.

5. Discussion

In order to calculate the sensitivity of PSD it is also possible to adopt first order Taylor series expansion. The Derivation-based sensitivity of PSD is given in the following equation.(37) $$\mathit{\delta }{S}_{xx}(\mathit{\omega })=\sum _{i=1}^{Ne}-(H_I^{\ast }(\mathit{\omega })F_i{H}_I^T(\mathit{\omega })\mathit{\delta }Z{H}_I^T(\mathit{\omega })+H_I^{\ast }(\mathit{\omega })(\mathit{\delta }{Z}^{\ast })H_I^{\ast }(\mathit{\omega })F_i{H}_I(\mathit{\omega }))S_{ffi}(\mathit{\omega })$$(37)

A comparison case of slight damage in elements 8 and 12 of plate model is considered to highlight the merits of exact sensitivity equation. The results presented in Figure 25, proves that with the same condition in term of excitation frequency range, measurement and excitation location and upper and lower bounds the exact method exhibits more precise solution of the comparison case in the noise-free situation.

Figure 25. Comparison case in noise-free condition.

The robustness of both sensitivity equations is also tested in presence of 5% random noise in simulated PSD of damaged model. The damage detection by Monte-Carlo analysis of simulated noisy sets proves that the exact method leads to more precise results as compared to derivation-based method. The Closeness Index (CI) of the comparison case in both conditions is given in the Table 6. According to definition [15], the closer the CI is to 1 the more precise the solution is. It is worth noting here, the derivation-based method does not converge to a solution in case of all the studied cases before and in presence of higher level of noise. The comparison of the mean parameter estimation results is given in the Figure 26.

Table 6. Closeness index for comparison case.

	Noise-free	5% Noise
Derivation-based sensitivity	0.7936	0.5932
Proposed sensitivity	0.9697	0.7522

Figure 26. Comparison case with 5% noise.

6. Conclusion

Most structures comprised of shell and plate elements should withstand harsh loading situation in their lifetime. Using PSD as metric, this paper presents a FE model updating method and investigates its performance in case of plate and shell models. The sensitivity equation is derived in an exact manner and an approximate formula is used to overcome incomplete measurement adversities. Several damage cases are studied for both kinds of elements in low frequency range and under assumption that only the first few natural frequency of damaged models is available. Damage cases test the performance of method against different level and location of damage. The method is successfully applied to both models considering the adverse effects of incomplete measurement both in spatial coordinate and modal parameters. The amount of equations that are provided enables the method to work at low frequency ranges. The method converges to high quality results using six lowest natural frequencies and selecting excitation frequency ranges up to this level. The method proves to be robust against measurement and mass modelling error.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix 1

The change of the FRF due to the damage in the structural elements is formulated in the following manner. The equation of motion of a MDOF structural system is,(A.1) $$F(\mathit{\omega })=(K-{\mathit{\omega }}^{2}M+i\mathit{\omega }C)X(\mathit{\omega })$$(A.1)

Due to the damage the equation of motion becomes,(A.2) $$F(\mathit{\omega })=(\underset{K_D}{\underbrace{(K_I+\mathit{\delta }K)}}-{\mathit{\omega }}^2\underset{M_D}{\underbrace{(M_I+\mathit{\delta }M)}}+i\mathit{\omega }(\underset{C_D}{\underbrace{C_I+\mathit{\delta }C}}))(\underset{X_D}{\underbrace{X(\mathit{\omega })+\mathit{\delta }X(\mathit{\omega })}})$$(A.2)

Expansion of the Equation (A.2(A.2) $$F(\mathit{\omega })=(\underset{K_D}{\underbrace{(K_I+\mathit{\delta }K)}}-{\mathit{\omega }}^2\underset{M_D}{\underbrace{(M_I+\mathit{\delta }M)}}+i\mathit{\omega }(\underset{C_D}{\underbrace{C_I+\mathit{\delta }C}}))(\underset{X_D}{\underbrace{X(\mathit{\omega })+\mathit{\delta }X(\mathit{\omega })}})$$(A.2) ) and subtracting Equation (A.1(A.1) $$F(\mathit{\omega })=(K-{\mathit{\omega }}^{2}M+i\mathit{\omega }C)X(\mathit{\omega })$$(A.1) ) from it yields,(A.3) $$\mathit{\delta }X(\mathit{\omega })=-H_D(\mathit{\omega })(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M+i\mathit{\omega }\mathit{\delta }C)X(\mathit{\omega })$$(A.3)

This by assuming the same loading condition for the intact and damaged structure will transform as, which expresses the change of FRF due to damage.(A.4) $$\mathit{\delta }H(\mathit{\omega })=-H_D(\mathit{\omega })(\mathit{\delta }K-{\mathit{\omega }}^2\mathit{\delta }M+i\mathit{\omega }\mathit{\delta }C)H(\mathit{\omega })$$(A.4)