Investigation of a perturbation-based model updating approach for structural health and event monitoring

Abstract

With the intention to detect structural events, for example, occurrence of damage and assessing their severity under service, structural health and event monitoring proved itself as a reliable evaluation tool. While many algorithmic approaches exist, model-based monitoring offers the prominent advantage to enable extensive analysis of damage states based on a structural digital twin. In order to accurately replicate the dynamical behaviour of a real-world structure, inherent physical properties, for example, mass, stiffness and structural damping, of the digital twin might be corrected from nominal values. This can be achieved by numerical model updating procedures. In this study, the application of a model updating approach is presented to eliminate modal discrepancies. Based on the perturbated modal dynamic residual, an updating equation is formulated. Stiffness and mass correction terms with preserved orthogonality and symmetry conditions are determined by the method of least squares. A case study using computationally generated modal data is conducted, evaluating the numerical updating performance in terms of correction accuracy and reproduction of modal parameters. Finally, the transient structural response to an impulse excitation is compared between updated and reference model. The gathered results prove the method’s suitability as accurate and robust updating procedure, fostering its application in model-based monitoring frameworks.

1. Introduction

With advances in sensor technology and data-driven analysis, the concept of structural health and event monitoring (SHEM) has become of major importance to various engineering disciplines over the last decade in civil, military, aerospace and automotive applications (Balageas et al., 2016; Haidarpour and Tee, 2020). Enabling an almost real-time, full-scale evaluation of the structural integrity without interrupting operation cycles, an embedded monitoring system can be understood as a non-destructive evaluation technique (NDE). This enables condition-based maintenance and eliminates the need to interrupt operations for inspection. In the very special case of aerospace applications, neither manual inspection nor maintenance might be possible at all. However, mission success depends on the structural reliability of spacecraft components and demands early detection of potential issues to prevent catastrophic failures.

According to Beard et al. (2005), SHEM systems consist of three similar components: (1) A permanently attached or embedded sensor/actuator network, (2) a data acquisition system (DAQ) with potential preprocessing capabilities and (3) a data processing unit (DPU). The DPU, considered the core of each SHEM system, is tasked with extracting meaningful features from large, complex datasets to detect singular events, such as sudden impacts or the evolution of fatigue, before they threaten operational safety. However, the success of these predictions relies heavily on the chosen algorithm, making the appropriate selection of the approach crucial, as the high volume and complexity of the data make analysis challenging. Further information on monitoring strategies can be found in Figueiredo et al. (2009), Haynes (2014), Sun and Qin (2014), and Yang et al. (2014). Algorithmic approaches allow distinction between data-driven and model-based methodologies (Fan and Qiao, 2011; Vagnoli et al., 2001; Worden and Manson, 2007). Some authors offer other differentiations, for example, local versus global (Balageas et al., 2016) and linear versus nonlinear (Haidarpour and Tee, 2020), but are beyond the scope of this paper.

Data-driven: Associated procedures use pattern recognition for defining statistical models as baseline projection from data gathered under service (Farrar and Worden, 2012). This allows incorporating environmental conditions and various load cases with ease. However, ensuring high-quality datasets that account for uncertainties and variations is crucial for robust decision thresholds, as outlined in Beck et al. (2001). The preprocessing task of achieving this is challenging, especially when applied on a full-system level. Some applications even demand real-time decisions based on sensor-acquired data which requires high-performance computing (Tibaduiza Burgos et al., 2020).

Model-based: Corresponding approaches rely on a structural model, for example, finite element (FE) model. If such a reliable model is available, load characteristics estimated from sensor signals can be implemented into simulations (Ching, 2006; Giurgiutiu, 2015; Otsuki, 2022). This allows evaluation of possibly induced damages by numerical analysis and inverse techniques (Gomes et al., 2018). However, FE models inherit certain inaccuracies which cause the structural model to deviate from experimentally gathered results (Ereiz et al., 2022). Aside from modelling errors, for example, idealized geometries and simplified contact modelling, inaccuracies stem from physical parameters, that is, nominal mass and stiffness, friction and structural damping that do not match real-world properties. To eliminate those discrepancies, the FE model’s implemented design parameters must be adjusted.

From the operator’s perspective, model-based monitoring allows comprehensive understanding of the system’s structural condition. This offers a superior advantage over data-driven analysis but requires a highly sophisticated model. Implementing such structural model can be challenging as it necessitates a thorough understanding of the system’s physical properties and environmental interactions. If not all relevant information is accessible, the implementation of such a mature model can be problematic even with expertise. In this case, numerical model updating (NMU) techniques can be utilized to improve the accuracy of the model and ultimately enhance the reliability of the monitoring results (Friswell and Mottershead, 1995). NMU has gained broad attention within SHEM applications since it allows deriving high-fidelity baseline models for model-based health monitoring. Once such a digital twin is implemented, NMU allows structural damage identification by various methods which were reviewed from Sharry et al. (2022).

2. Delimitation to existing methods

In this study, mass and stiffness updating is considered. Numerous prior research efforts have been dedicated to this approach. Similar to this study, their updating equations minimize the modal dynamic residual and enforce orthogonality resp. symmetry constraints. Thus, the derivation of the hereinafter investigated updating procedure might be partially similar to existing methods mentioned below. The following section summarizes prior research efforts that are most closely related to the present study.

In the late 1980s, first proposals for individual mass or stiffness updating have been published by Baruch and Itzhack (1978) followed by Berman (1979). They introduced modal based updating for either mass or stiffness matrix while assuming sufficient accuracy for the remainder. Baruch (1982) proposed an early approach to correct the system matrices sequentially by formulating two different updating variations based on the assumption that either system’s mass or stiffness is known with higher confidence. Caesar and Peter (1998) derived a generalized linear least square formulation for individual mass and stiffness updating with included symmetry, orthogonality and additional force equilibrium constraints. However, considering the matrices’ interdependent nature towards the modal dynamic residual, the practitioner is challenged to decide on the more accurate property while limiting the full updating potential as one matrix remains unaltered during each updating step. Consequently, a simultaneous updating methodology for both mass and stiffness should be favoured.

Sarmadi et al. (2016) expanded the mass and stiffness orthogonality based on the dynamic discrepancy theory, establishing two improved updating formulations to determine correction terms as discrepancy matrices. Similar to this study, they transformed the updating formulation into a non-square linear equation system using the Kronecker tensor product (see (21)) and determined the solution by iterative methods. However, their methodology does not facilitate the removal of constrained and sparse matrix entries which becomes a notable drawback when dealing with large-scale systems.

Iterative methods, as the one presented in this study, often incorporate the modal dynamic residual or the frequency response function (FRF) as part of their objective function. Arora (2011) investigated the achieved updating accuracy in terms of the obtained frequency response function (FRF) when applying direct or iterative updating approaches. Dong and Wang (2016) evaluated the frequency domain updating performance of two different updating procedures where one approach involves the minimization of the modal dynamic residual. They also compared different iterative solvers to find suitable updating parameters for the system’s stiffness matrix. Li et al. (2018) investigated the structural identification process for stiffness correction, utilizing the modal dynamic residual approach which was minimized by LSQR (least squares method). In their study, they also compared the performance of different gradient-based solvers. Based on this, Li and Wang (2019) also addressed the problem of unfeasible computational efforts when dealing with large datasets. Similar to this study, they exploited sparsity to significantly reduce the size of the optimization problem. However, FRF-based updating formulations are constrained by their inability to target individual modes or eigenfrequencies for updating. Instead, they update the entire frequency spectrum captured by the FRF. The lack of granularity in selecting particular modes or eigenfrequencies limits practitioners seeking to fine-tune specific dynamic characteristics of their structures.

Further studies investigated the application of stochastic and metaheuristic optimization such as evolutionary algorithms (EA) for model updating. Gomes et al. (2018) provided a comprehensive review of vibration-based inverse techniques for damage detection and identification which included FE model updating by genetic algorithms (GA), swarm optimization, etc. They also summarized objective functions which were assembled from modal parameters such as natural frequencies, mode shapes or FRF of the reviewed studies. Jafarkhani and Sami (2011) evaluated the updating performance of stochastic optimization algorithms for finite element model updating to monitor dispersed civil infrastructures. The damage detection results revealed fair agreement between the qualitative reported visual inspections and identified damages through model updating, proving the potential of EA for model-based health monitoring. Another real-world application was demonstrated by Sun et al. (2020), using a two-phased updating technique to assess the mechanical state of an in-service cable-stayed bridge by GA which significantly improved FE model accuracy. However, as those metaheuristic approaches can be powerful for global optimization and exploring complex solution spaces, they demand high computational costs when dealing with large parameter spaces. As in direct mass and stiffness matrix updating every non-sparse matrix entry can be considered an optimization parameter, the parameter space expands rapidly and makes evolutionary algorithms suffer from slow convergence. On the contrary and as proposed in this study, gradient-based methods leverage the gradient information to iteratively refine the algorithmic parameters in a more adaptive manner which often leads to faster convergence.

This study presents a methodology for simultaneously updating of mass matrix M and stiffness matrix K of a structural FE model by correction terms ΔM and ΔK to reproduce set of modal parameters (specific eigenfrequencies ω and mode shapes ϕ). The updating procedure is derived by perturbating the modal dynamic residual which is minimized by a gradient-based solver. The updating equations account for orthogonality and symmetry constraints, exploit sparsity patterns to deal with large datasets and allow to precisely target specific modal sets within the updating objective. Ultimately, a straightforward formulation allows convenient implementation for the practitioner while its high correction accuracy encourages confident application in real-world scenarios.

3. Updating methodology

In structural dynamics, the behaviour of a structure subjected to dynamic loads can be described by second-order ordinary differential equations as follows:

$M\ddot{u}+C\dot{u}+Ku=F\left(t\right).$

(1)

with

$M\in {\mathbb{R}}^{n\times n}\hspace{0.17em}\text{the}\hspace{0.17em}\text{mass}\hspace{0.17em}\text{matrix},$

(2)

$D\in {\mathbb{R}}^{n\times n}\hspace{0.17em}\text{the}\hspace{0.17em}\text{damping}\hspace{0.17em}\text{matrix},$

(3)

$C\in {\mathbb{R}}^{n\times n}\hspace{0.17em}\text{the}\hspace{0.17em}\text{stiffness}\hspace{0.17em}\text{matrix},$

(4)

$\ddot{u}\in {\mathbb{R}}^n\hspace{0.17em}\text{the}\hspace{0.17em}\text{acceleration}\hspace{0.17em}\text{vector},$

(5)

$\dot{u}\in {\mathbb{R}}^n\hspace{0.17em}\text{the}\hspace{0.17em}\text{velocity}\hspace{0.17em}\text{vector},$

(6)

$u\in {\mathbb{R}}^n\hspace{0.17em}\text{the}\hspace{0.17em}\text{displacement}\hspace{0.17em}\text{vector},$

(7)

$F\left(t\right)\in {\mathbb{R}}^n\hspace{0.17em}\text{the}\hspace{0.17em}\text{time}\hspace{0.17em}\text{dependent}\hspace{0.17em}\text{external}\hspace{0.17em}\text{force}\hspace{0.17em}\text{vector}.$

(8)

To examine the structure’s natural vibrational characteristics, (1) is further simplified by neglecting damping forces and external loading. Assuming the stationary solution $u\left(t\right)={\varphi }_i{e}^{{\omega }_{i}t}$ with ${\varphi }_i\in {\mathbb{R}}^n$ , the eigenvector to eigenfrequency ${\omega }_i\mathbb{\in }\mathbb{C}$ , the generalized eigenvalue problem is derived:

$\left(K-{\omega }^{2}M\right)\varphi =\widehat{R}$

(9)

with $\widehat{R}\in {\mathbb{R}}^n$ , the modal dynamic residual can be derived from (1). For suitable eigenvector–eigenfrequency pairs (ϕ _i , ω _i ), the inertial and elastic forces of (9) preserve dynamic equilibrium. This is essentially quantified by the modal dynamic residual, where a small residual indicates sufficient agreement to the underlying differential equation. Conversely, a substantial residual indicates a mismatch between the (possibly experimentally gathered) modal solutions and inherent mass resp. stiffness properties. As the ultimate goal of NMU is to update mass matrix M and stiffness matrix K, the following unconstrained optimization problem summarizes the updating objective:

$\underset{K,M}{\min }{\Vert K{\varphi }_i-{\omega }_i^{2}M\varphi \Vert }_2.$

(10)

with ${\Vert \cdot \Vert }_2=\sqrt{(}{\sum }_i^n\left|{\widehat{r}}_i\right|)$ the Euclidean norm.

To introduce the updating procedure, subscripts indicate the affiliation to different models: Firstly, a reference model (RM) shall be labelled by the subscript R. The RM’s structural properties are known in advance and serve as benchmark to evaluate the updating performance. Secondly, an initial model (IM) with yet unknown parameters is identified by subscript I. The purpose of the latter is to reproduce the RM’s modal parameters based on correction terms ΔK and ΔM as a solution to the updating equation.

With the assembled mass and stiffness matrices from the IM M _I , $K_I\in {\mathbb{R}}^{n\times n}$ with n DOFs, the modal dynamic residual (9) can be formulated as matrix equation:

${\Vert \left(K_I{\mathrm{\Phi }}_I-M_I{\mathrm{\Phi }}_I{\mathrm{\Lambda }}_I\right)\Vert }_2=0$

(11)

where ${\mathrm{\Phi }}_I=\left[{\varphi }_{I,1},\dots ,{\varphi }_{I,n}\right]\in {\mathbb{R}}^{n\times n}$ is the modal matrix and Λ _I the spectral matrix with eigenvalues $\left[{\lambda }_1,\dots ,{\lambda }_n\right]\in {\mathbb{R}}^{nxn}$ , that is, eigenfrequencies $\left[{\omega }_1^2,\dots ,{\omega }_n^2\right]$ on the main diagonal. The updating approach is based on perturbating the modal dynamic residual (11) by introducing correction terms ΔM, $\mathrm{\Delta }K\in {\mathbb{R}}^{n\times n}$ , $\mathrm{\Delta }\mathrm{\Phi }\in {\mathbb{R}}^{n\times m}$ and $\mathrm{\Delta }\mathrm{\Lambda }\in {\mathbb{R}}^{m\times m}$ with m the number of included modes and eigenfrequencies to the updating process. The correction terms set out the link between RM and IM:

$\begin{array}{ll}{M}_R=\mathrm{\Delta }M+M_I& K_R=\mathrm{\Delta }K+K_I\\ {\mathrm{\Phi }}_R=\mathrm{\Delta }\mathrm{\Phi }+{\mathrm{\Phi }}_I& {\mathrm{\Lambda }}_R=\mathrm{\Delta }\mathrm{\Lambda }+{\mathrm{\Lambda }}_I\end{array}$

(12)

Usually, m ≪ n holds true since the considered set of modal parameters represents only a small subset to the solution space of (11). In a real-world updating scenario, only M _R and K _R remain unknown since Φ _R and Λ _R can be determined by experimental modal analysis (EMA) of the examined structure. Including all correction terms (12) into the modal dynamic residual (11) yields different updating formulations: Either solely updating of stiffness by ΔK or simultaneous updating of both mass and stiffness by ΔM and ΔK. In the following, both approaches are introduced in Section 3.1 and Section 3.2, and the reader can directly skip to the desired chapter based on their updating preference. Isolated updating of M might be possible as well following the procedure in Section 3.1. However, several authors suggest to neglect isolated mass matrix updating since assembling from nominal material properties achieves sufficient accuracy (Halevi et al., 2003).

3.1. Isolated updating of stiffness matrix

For the purpose of isolated updating of stiffness matrix K, mass correction ΔM is disregarded from the modal dynamic residual (11):

${\Vert \left(K_I+\mathrm{\Delta }K\right)\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)-\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)\left({\mathrm{\Lambda }}_I+\mathrm{\Delta }\Lambda \right)\Vert }_2=0$

(13)

Consideration of (11) allows reducing to (14):

${\Vert K_I\mathrm{\Delta }\mathrm{\Phi }+{\mathrm{\Phi }}_I\mathrm{\Delta }K+\mathrm{\Delta }K\mathrm{\Delta }\mathrm{\Phi }-{M_1\mathrm{\Phi }}_I\mathrm{\Delta }\Lambda -M_I\mathrm{\Delta }\mathrm{\Phi }\mathrm{\Delta }\Lambda -M_I\mathrm{\Delta }\mathrm{\Phi }{\mathrm{\Lambda }}_I\Vert }_2=0$

(14)

Finally, separation of the correction term ΔK leads to the following equation:

$\Vert \mathrm{\Delta }K\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)-(M_I{\mathrm{\Phi }}_I\mathrm{\Delta }\Lambda -K_I\mathrm{\Delta }\mathrm{\Phi }+\hspace{0.17em}{M}_I\mathrm{\Delta }\mathrm{\Phi }\mathrm{\Delta }\Lambda +M_I\mathrm{\Delta }\mathrm{\Phi }{\mathrm{\Lambda }}_I)=0$

(15)

At this point, matrix norm ‖ ⋅‖₂ is neglected for improved readability even though mapping matrix–vector terms to scalar quantities is only applicable when applying matrix norms. Finally, (15) can be reshaped to a linear equation system:

$AX=B$

(16)

with

$\begin{array}{l}A=\mathrm{\Delta }K\in {\mathbb{R}}^{n\times n}\end{array}$

(17)

$\begin{array}{l}X={\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\in {\mathbb{R}}^{n\times m}\end{array}$

(18)

$B=M\mathrm{\Delta }\Lambda {\mathrm{\Phi }}_I+M{\mathrm{\Lambda }}_I\mathrm{\Delta }\mathrm{\Phi }I+M\mathrm{\Delta }\Lambda \mathrm{\Delta }{\mathrm{\Phi }}_I-K_I\mathrm{\Delta }{\mathrm{\Phi }}_I\in {\mathbb{R}}^{n\times m}$

(19)

To solve for matrix A resp. ΔK, rearrangement of (16) is required. This can be achieved by the following expression:

$\underset{C}{\underbrace{\left(I\otimes X^T\right)}}A=B$

(20)

with I the identity matrix and ⊗ the Kronecker tensor product:

$\left(\widehat{A}\otimes \widehat{B}\right)=\left[\begin{array}{lll}{\widehat{a}}_{11}\widehat{B}& \cdots & {\widehat{a}}_{1m}\widehat{B}\\ ⋮& \ddots & ⋮\\ {\widehat{a}}_{n1}\widehat{B}& \cdots & {\widehat{a}}_{nm}\widehat{B}\end{array}\right]$

(21)

Due to the incomplete solution set Φ _I , Λ _I to the modal dynamic residual (11), $C\in {\mathbb{R}}^{nm\times m}$ remains highly unsquare. In a more detailed view, equation (20) can be stated as

$\left[\begin{array}{ccccc}{\underset{¯}{X}}_1& 0& \cdots & \cdots & 0\\ 0& {\underset{¯}{X}}_2& 0& \cdots & 0\\ & & \ddots & & \\ 0& \cdots & 0& {\underset{¯}{X}}_{m-1}& 0\\ 0& \cdots & \cdots & 0& {\underset{¯}{X}}_m\end{array}\right]\left[\begin{array}{l}{\underset{¯}{A}}_1\\ {\underset{¯}{A}}_2\\ ⋮\\ {\underset{¯}{A}}_{n-1}\\ {\underset{¯}{A}}_n\end{array}\right]=\left[\begin{array}{l}{\underset{¯}{B}}_1\\ {\underset{¯}{B}}_2\\ ⋮\\ {\underset{¯}{B}}_{n-1}\\ {\underset{¯}{B}}_n\end{array}\right]$

(22)

with

$\begin{array}{l}{\underset{¯}{X}}_i=I\otimes \left[X_{1,i},\dots ,X_{n,i}\right]=I\otimes {\left[{\varphi }_{I,i}+\mathrm{\Delta }{\varphi }_i\right]}^T,\end{array}$

(23)

${\underset{¯}{B}}_i={\left[B_{1,i}\cdots b_{n,i}\right]}^T=[M{\varphi }_{I,i}\mathrm{\Delta }{\omega }_i+M\mathrm{\Delta }{\varphi }_i{\omega }_i+M\mathrm{\Delta }{\varphi }_i\mathrm{\Delta }{\omega }_i-K_I\mathrm{\Delta }{\varphi }_i]$

(24)

and $I\in {\mathbb{R}}^{n\times n}$ . The terms ${\underset{\_}{A}}_i$ of A are basically the corresponding rows i of ΔK, which can be expressed as

${\underset{\_}{A}}_i={\left[A_{i,1},\dots ,A_{i,n}\right]}^T={\left[\mathrm{\Delta }{K}_{1,i},\dots ,\mathrm{\Delta }{K}_{n,i}\right]}^T$

(25)

To clarify, terms A and B of (20) are described by (25) and (24) while C is assembled by (23) when isolated updating of stiffness matrix K is foreseen.

3.2. Simultaneous updating of mass matrix and stiffness matrix

This section introduces the updating equation for simultaneous updating of stiffness matrix K and mass matrix M based on the same perturbation approach as in Section 3.1. As an entry point, the modal dynamic residual (11) with correction terms including ΔM and ΔK is considered:

${\Vert \left(K_I+\mathrm{\Delta }K\right)\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)-\left(M_I+\mathrm{\Delta }M\right)\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)\left({\mathrm{\Lambda }}_I+\mathrm{\Delta }\Lambda \right)\Vert }_2=0$

(26)

which is reduced with respect to (11) while separation of correction terms ΔM and ΔK allows rearranging to

$\Vert \mathrm{\Delta }K\left({\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi }\right)-\mathrm{\Delta }M({\mathrm{\Phi }}_I\mathrm{\Delta }\Lambda +{\mathrm{\Phi }}_I{\mathrm{\Lambda }}_I+\mathrm{\Delta }\mathrm{\Phi }\mathrm{\Delta }\Lambda +\mathrm{\Delta }{\mathrm{\Lambda }}_I)-(M_I{\mathrm{\Phi }}_I\mathrm{\Delta }\Lambda -K_I\mathrm{\Delta }\mathrm{\Phi }+M_I\mathrm{\Delta }\mathrm{\Phi }\mathrm{\Delta }\Lambda +M_I\mathrm{\Delta }\mathrm{\Phi }{\mathrm{\Lambda }}_I){\Vert }_2=0$

(27)

Similar to Section 3.1, equation (27) can be rearranged to a non-square linear equation system:

$AX=B$

(28)

with

$\begin{array}{l}A=\left[\mathrm{\Delta }K,\mathrm{\Delta }M\right]\in {\mathbb{R}}^{n\times 2n}\end{array}$

(29)

$X=[{\mathrm{\Phi }}_I+\mathrm{\Delta }\mathrm{\Phi };{\mathrm{\Phi }}_I\mathrm{\Delta }\Lambda +{\mathrm{\Phi }}_I{\mathrm{\Lambda }}_I+\mathrm{\Delta }\mathrm{\Phi }\mathrm{\Delta }\Lambda +\mathrm{\Delta }{\mathrm{\Lambda }}_I]\in {\mathbb{R}}^{2n\times m}$

(30)

$B=M\mathrm{\Delta }\Lambda {\mathrm{\Phi }}_I+M{\mathrm{\Lambda }}_I\mathrm{\Delta }{\mathrm{\Phi }}_I+M\mathrm{\Delta }\Lambda \mathrm{\Delta }{\mathrm{\Phi }}_I-K_I\mathrm{\Delta }{\mathrm{\Phi }}_I\in {\mathbb{R}}^{n\times m}$

(31)

Again, a rearrangement is required to solve for A as stated in (20). This time, A consists of vectors ${\underset{\_}{A}}_i$ which are the corresponding stacked rows i of ΔK and ΔM:

${\underset{¯}{A}}_i={\left[A_{i,1},\dots ,A_{i,2n}\right]}^T={\left[\mathrm{\Delta }{K}_{i,1},\dots ,\mathrm{\Delta }{K}_{i,n},\mathrm{\Delta }{M}_{i,1},\dots ,\mathrm{\Delta }{M}_{i,n}\right]}^T$

(32)

Since ${\underset{\_}{A}}_i$ now includes both mass and stiffness correction, ${\underset{\_}{X}}_i$ must be expanded to account for all correction terms of (29) as stated in (33).

$\begin{array}{l}{\underset{¯}{X}}_i=I\otimes \left[X_{1,i},\dots ,X_{n,i}\right]=[I\otimes {\left[{\varphi }_{I,i}+\mathrm{\Delta }{\varphi }_i\right]}^T,\\ I\otimes {\left[{\varphi }_{I,i}\mathrm{\Delta }{\omega }_i+{\varphi }_{I,i}{\omega }_{I,i}+\mathrm{\Delta }{\varphi }_i\mathrm{\Delta }{\omega }_i+\mathrm{\Delta }{\omega }_i\right]}^T],\end{array}$

(33)

${\underset{¯}{B}}_i={\left[B_{1,i}\cdots b_{n,i}\right]}^T=\left[M{\varphi }_{I,i}\mathrm{\Delta }{\omega }_i+M\mathrm{\Delta }{\varphi }_i{\omega }_i+M\mathrm{\Delta }{\varphi }_i\mathrm{\Delta }{\omega }_i-K_I\mathrm{\Delta }{\varphi }_i\right]$

(34)

and $I\in {\mathbb{R}}^{n\times n}$ . For simultaneous updating of mass and stiffness, terms A, B and C are defined by (32), (34) and (33).

3.3. Including matrix symmetry conditions

So far, the updating equations do not comply with symmetry conditions. For symmetric correction terms, the symmetry constraints

$\mathrm{\Delta }{K}^T=\mathrm{\Delta }K\hspace{1em}\text{and}\hspace{1em}\mathrm{\Delta }{M}^T=\mathrm{\Delta }M$

(35)

must be fulfilled. Symmetry constraints are included by extending (20) to

$\left[\begin{array}{c}C\\ {\underset{\_}{S}}_{i,j}\end{array}\right]\left[A\right]=\left[\begin{array}{l}B\\ 0\end{array}\right]$

(36)

with

${\underset{¯}{S}}_{i,j}=\left[{\delta }_{1,1},\dots ,{\delta }_{i,j},\dots ,{\delta }_{j,i},\dots ,{\delta }_{n,n}\right]\in {\mathbb{R}}^{n^2}$

(37)

The vector ${\underset{¯}{S}}_{i,j}$ indicates two off-diagonal matrix elements A_i,j and A_j,i to be symmetric by utilizing a modified Kronecker delta:

${\delta }_{i,j},{\delta }_{j,i}=\{\begin{array}{ll}1,& \text{if}\hspace{0.17em}i\ne j\\ 0,& \text{else}\end{array}$

(38)

The modification denotes the Kronecker delta’s inverted properties which is equal to 1 − δ_i,j. By assigning i, j to all non-zero matrix elements for each ${\underset{\_}{S}}_{i,j}$ , ΔK and ΔM maintain symmetry. When extending (36) by ${\underset{\_}{S}}_{i,j}$ , mass and stiffness correction must be considered individually.

3.4. Including mass and stiffness orthogonality

A similar approach can be adopted to achieve mass and stiffness orthogonality regarding the modal set Φ _R , Λ _R . For mass orthonormal modes, the orthogonality constraints

${\varphi }_{R,i}^{T}M{\varphi }_{R,j}=\{\begin{array}{ll}1& \text{if}\hspace{0.17em}i=j\\ 0& \text{else}\end{array}$

(39)

${\varphi }_{R,i}^{T}K{\varphi }_{R,j}=\{\begin{array}{ll}{\omega }_{R,i}^2& \text{if}\hspace{0.17em}i=j\\ 0& \text{else}\end{array}$

(40)

must be fulfilled. Substituting M with (M _I + ΔM) resp. K with (K _I + ΔK) allows deriving the following expression:

${\varphi }_{R,i}^T\mathrm{\Delta }M{\varphi }_{R,j}=\{\begin{array}{ll}1-{\varphi }_{R,k}^T{M}_I{\varphi }_{R,k}& \text{if}\hspace{0.17em}i=j\\ 0& \text{else}\end{array}$

(41)

${\varphi }_{R,i}^T\mathrm{\Delta }K{\varphi }_{R,j}=\{\begin{array}{ll}{\omega }_{R,k}^2-{\varphi }_{R,k}^T{K}_I{\varphi }_{R,k}& \text{if}\hspace{0.17em}i=j\\ 0& \text{else}\end{array}$

(42)

Since ΔM and ΔK are square matrices, (41) and (42) can be included into (20) by

$\left[\begin{array}{l}C\\ {\underset{¯}{S}}_{i,j}\\ {\underset{¯}{Q}}_k\end{array}\right]\left[A\right]=\left[\begin{array}{l}B\\ 0\\ q_k\end{array}\right]$

(43)

with

${\underset{\_}{Q}}_{R,k}=X_{R,k}^T\otimes X_{R,k}^T\in {\mathbb{R}}^{n^2}$

(44)

The symmetry constraints (36) are already included in (44) for the sake of completeness. To demand orthogonality of pairing modes, q _k can be set to

$q_k=\{\begin{array}{ll}{\omega }_{R,k}^2-{\varphi }_{R,k}^T{K}_I{\varphi }_{R,k}& \mathrm{w}.\mathrm{r}.\mathrm{t}.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }K\\ 1-{\varphi }_{R,k}^T{M}_I{\varphi }_{R,k}& \mathrm{w}.\mathrm{r}.\mathrm{t}.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }M\end{array}$

(45)

while q _k = 0 is set for non-pairing modes. Again, mass and stiffness orthogonality must be considered individually.

3.5. Additional constraints

Further constraints can be considered by excluding specific matrix elements from the updating procedure. This applies to, for example, fixed single point constraints in terms of translation and rotation or zero-elements of M and K whose altering would generate non-existent stiffness and mass coupling. Therefore, all corresponding rows of A and columns of terms C, ${\underset{\_}{S}}_{i,j}$ and ${\underset{\_}{Q}}_k$ shall be removed from (20), (37) and (44). This allows significant reduction of (20) since M and K are usually sparse, band-structured matrices.

3.6. Solving as optimization problem

Finally, the updating equation allows determining correction terms ΔM and ΔK to given mass and stiffness matrices M _I and K _I such that the modal subset Λ _R , Φ _R is reproduced. This can be summarized in the unconstrained optimization problem:

$\underset{\mathrm{\Delta }K,\mathrm{\Delta }M}{\min }{\Vert \left(K_I+\mathrm{\Delta }K\right){\mathrm{\Phi }}_R-\left(M_I+\mathrm{\Delta }M\right){\mathrm{\Phi }}_R{\mathrm{\Lambda }}_R\Vert }_2$

(46)

Under consideration of symmetry (35) and orthogonality constraints (39) and (40), the non-square linear equation system (20) optionally with (36) and (43) is solved instead.

$\underset{\mathrm{\Delta }K,\mathrm{\Delta }M}{\min }{\Vert CA-B\Vert }_2$

(47)

4. Numerical investigation

In order to evaluate the updating performance and to demonstrate the approach’s suitability as SHEM tool, a numerical investigation is conducted. In this study, a rectangular cantilever plate made of aluminium is utilized as specimen with geometric dimensions summarized in Table 1. The updating objective is to reproduce the first six modes Φ _R and eigenfrequencies Λ _R by updating the mass matrix M _I and stiffness matrix K _I . Since the updating examines a virtual test case, Φ _R and Λ _R are computed numerically. In a real-world scenario, Φ _R and Λ _R might be determined by experimental modal analysis without altering the overall updating procedure. The results of the simulative study are summarized in Section 4.2.

Table 1.

Model geometries.

	Dimensions
	Height (mm)	Width (mm)	Depth (mm)
Reference model	420	420	1
Initial model	420	420	1

4.1. Correlation metrics

In this study, different modal and time domain correlation tools are used to identify discrepancies in the dynamic structural response. Those are listed in the following with the introduced notation for IM and RM:

Modal Assurance Criterion: The Modal Assurance Criterion (MAC) allows evaluating the level of correlation between sets of pairing mode shapes. A high correlation is indicated by an MAC value of 1 while a lack of correlation results in an MAC value of 0. The MAC is formulated by:

${\text{MAC}}_{k,l}=\frac{{\left|{\varphi }_{I,k}^T{\varphi }_{R,l}\right|}^2}{\left({\varphi }_{I,k}^T{\varphi }_{I,k}\right)\left({\varphi }_{R,l}^T{\varphi }_{R,l}\right)}$

(48)

According to Pástor et al. (2012), MAC is sensitive to nodal DOFs with large modal displacements and tends to neglect minor mode shape deviations which might contribute to a decrease of correlation. Therefore, it should be mainly used for global mode shape comparison instead of local correspondence assessment.

Coordinate Modal Assurance Criterion: The Coordinate Modal Assurance Criterion (CoMAC) identifies DOFs that diminish the degree of mode shape correspondence. This allows to monitor their individual correlation drift throughout the updating process. The CoMAC is computed for every nodal DOF over the set of mode pairs with high correlation (pairing modes are determined by MAC). Therefore, it allows evaluating local mode shape correspondence but is not able to summarize global correlation quality. The CoMAC is formulated as:

${\text{CoMAC}}_q=\frac{{\sum }_{r=1}^m{\left|{\varphi }_{I,r}^{\left(q\right)}{\varphi }_{R,r}^{\left(q\right)}\right|}^2}{{\sum }_{r=1}^m{\left({\varphi }_{I,r}^{\left(q\right)}\right)}^2{\sum }_{r=1}^m{\left({\varphi }_{I,r}^{\left(q\right)}\right)}^2}$

(49)

where q is the particular nodal DOF to be examined. To illustrate regions with lack of correlation, the CoMAC will be projected onto the mode shapes as 1 − CoMAC _q in the following (see Figure 3).

Test Response Assurance Criteria: The Test Response Assurance Criteria (TRAC) allows correlating time domain data between different structural dynamical responses. Similar to MAC, a TRAC value of 1 indicates high correlation while low similarity results in a value of 0.

$\text{TRAC}=\frac{{\left(x_I^T{x}_R\right)}^2}{\left(x_I^T{x}_I\right)\left(x_R^T{x}_R\right)}$

(50)

According to Reyes et al. (2017), TRAC can be understood as a regression analysis with the TRAC value being the R² coefficient.

Percentage Error ϵ: Even though not a common correlation tool, the percentage error can be utilized to compare the correlation of individual matrix elements which can be useful for evaluating the local updating accuracy. The percentage error for both mass and stiffness elements is stated here for the sake of completeness:

${ϵ}_{k,ij}=\left|1-k_{I,ij}^{*}/k_{R,ij}\right|\cdot 100$

(51)

${ϵ}_{m,ij}=\left|1-m_{I,ij}^{*}/m_{R,ij}\right|\cdot 100$

(52)

Compared to more commonly used matrix norms, such as the L₂-Norm, the percentage error provides an insight into local error distributions. Matrix norms provide a measure of the overall correction accuracy of a matrix but may not accurately reflect local error accumulations or magnitude of entries within the matrix. This can result in a misleading assessment of the correction accuracy if relying solely on matrix norms for error evaluation.

4.2. Simulative setup

For the numerical investigation of the study, the structural time and frequency response data are computed using the general-purpose finite element analysis (FEA) environment ABAQUS FEA, developed by Dassault Systèmes SE (3DS). The finite element model is created twice, each with different material properties. The first model, referred to as the reference model (RM), inherits the assumed-to-be-correct material properties and correctly assembled mass and stiffness matrices Φ _R and Λ _R . The modal properties of the RM (Φ _R and Λ _R ) serve as the reference. The second model – initial model (IM) – is assigned strongly distinct uniform mass and stiffness properties to emphasize their updating capability. Consequently, the system matrices M _I and K _I and modal properties Λ _I and Φ _I differ from the RM. The assigned homogeneous material properties of RM and IM are summarized in Table 2.

Table 2.

Model’s material properties.

	Material properties
	ρ (g/mm3)	E (GPa)	μ [-]
Reference model	2.7e-2	70	0.3
Initial model	3.5e-2	60	0.3

Modelling is realized by solid second-order hexahedral elements as standard volume elements of this study. Despite the fact that thin-walled structures may be more efficiently modelled with structural plate or shell elements, standard elements are used to emphasize the updating procedure’s versatility independent of element type constraints. It is noteworthy to mention that different element types, global mesh refinement, etc., will introduce variations in mass and stiffness matrix entry’s magnitude. However, the updating formulation remains unaltered since core assumptions, that is, matrix symmetry and orthogonality, persist despite the choice of modelling. As this study is intended to investigate the updating performance on a general simplified structural model with uniform material and mesh properties, effects of advanced modelling techniques, for example, non-uniform material/mesh partitioning, are not covered in the numerical analysis.

The described numerical setup and boundary conditions, resulting from the clamping at the specimen’s bottom edges, are illustrated in Figure 1. Additionally, the resulting mode shapes Φ _R and eigenfrequencies ω _R of the RM are displayed in Figure 2. The updated IM’s matrices and modal properties are labelled by (*) as $M_I^{*}$ , $K_I^{*}$ , and ${\mathrm{\Phi }}_I^{*}$ , ${\mathrm{\Lambda }}_I^{*}$ . The numerical study allows extensive preliminary analysis of the updating performance before approaching real-world scenarios. Three key aspects are formulated as scientific questions and investigated in this study:

Question 1: Can the desired modal subset Λ _R , Φ _R be reproduced from the IM’s updated mass and stiffness matrices?

Question 2: Can a sufficient correction term accuracy of ΔK, ΔM w.r.t orthogonality and symmetry constraints be achieved?

Question 3: Can the updated IM be utilized to compute a transient structural response caused, for example, from an impact event which is equivalent to the RM?

Figure 1.

Description of the numerical setup. Subfigure (a): Illustration of cantilever plate’s FE mesh with dimensions, boundary conditions, measurement point X(t) and excitation point F(t). Subfigure (b): Cantilever plate’s depth.

Figure 2.

Normalized mode shapes Φ _R and eigenfrequencies ω _R of the RM from 1st to 6th mode. For the purpose of illustration, the modal displacements are normalized and scaled.

In the following, Section 4.3, Section 4.4 and Section 4.5 are dedicated to answer those questions. In order to examine Question 1, the modal subset Φ _R , Λ _R of RM and updated IM are compared. To evaluate their modal correlation, both MAC and CoMAC are used. Additionally, preservation of mass and stiffness orthogonality after updating are investigated. To evaluate Question 2, the Percentage Error of mass matrix and stiffness matrix elements are analyzed by the empirical cumulative distribution function. Regarding Question 3, an impulse response of the updated IM’s is compared to the RM. The response is measured as nodal acceleration from measurement point X(t). To induce vibrations, an impact event is simplified by an impulse function F(t) applied to an excitation point as nodal force. The modelling of the impulse excitation function is described in Section 4.5. TRAC is used to examine the time domain data correlation. Additionally, the power spectral density (PSD) is computed to allow comparison in the frequency domain.

4.3. Question 1: Reconstruction of modal parameters

The computed modal parameters from RM, IM and updated IM are summarized in Table 3. To allow comparison between different eigenvalue deviations, the individual percentage error completes the table.

Table 3.

Comparison of eigenfrequencies.

		Modes
Before updating		ϕ 1	ϕ 2	ϕ 3	ϕ 4	ϕ 5	ϕ 6
Eigenfreq.	Reference model (Hz)	23.5	35.4	73.9	93.8	119.7	136.6
	Initial model (Hz)	17.8	28.5	67.8	88.5	113.2	129.3
	Percentage error (%)	24.3	19.5	8.3	5.7	5.4	5.3
After updating
Eigenfreq.	Reference model (Hz)	23.5	35.4	73.9	93.8	119.7	136.6
	Updated model (Hz)	23.1	35.7	74.2	93.9	118.4	136.7
	Percentage error (%)	1.7	0.8	0.5	0.1	1.1	0.1

Before updating, considerable discrepancies appear throughout the examined eigenfrequency spectrum, resulting in percentage errors ranging from approx. 5−25%. The lower modes, such as the 1st and 2nd modes, tend to have larger eigenfrequency deviations with errors of 19.5−24.3%. In contrast, 4th to 6th modes remain with similar errors at 5.3−5.7%. After the updating, the eigenfrequencies of IM and RM closely match across the considered frequency spectrum. The diminished eigenfrequency deviations result in percentage errors below 2%. The highest deviation of 1.7% is found again in the fundamental mode. The tendency of enforcing high-order modes to be more accurate might be explained by (33) and (34) since the included eigenfrequency difference Δω _i acts as natural weighting factor for high-frequency content. Although the purpose of this unintentional weighting might be questioned, it proves to have a meaningful consequence since higher modes are more prominent in the dynamical acceleration response which are usually measured by acceleration signals. Lower modes, on the other hand, dominantly participate in the overall dynamical response. Therefore, this eigenfrequency weighting gives balanced attention to both high and low frequency content. Depending on the included modes within the updating procedure, this balancing will eventually suffer and demand additional compensation weights. The magnitude of this compensation might be customized for the individual application.

The correlation of the updated IM’s mode shapes ${\mathrm{\Phi }}_I^{*}$ to the RM’s mode shapes Φ _R is quantified by the MAC and depicted in Figure 3, Subfigure A. Pairing modes show high correlation (MAC ≥0.96) while sufficient distinction between the remaining modal subset is given (MAC ≤0.05). Severe mode shape distortions don’t become apparent based on the overall agreement of pairing modes. Additionally, the existence of spurious modes as solely numerical solutions from (47) to (11) with no physical meaning must be mentioned (Carvalho et al., 2007). Those physical unmeaningful modes might be generated throughout the updating procedure, for example, within eigenfrequency gaps, and can be recognized by shifting former pairing modes apart from the main diagonal. Since all pairing modes remained on the main diagonal, it can be concluded that the updating procedure does not produce spurious modes in the considered frequency spectrum of the updated IM.

Figure 3.

Mode shape correlations from updated IM’s mode shapes ${\varphi }_I^{*}$ and RM’s mode shapes ϕ _R . Subfigure (a): Modal Assurance Criterion (MAC) from 1st–6th mode. Subfigure (b): Projection of subtracting 1 − CoMAC for DOFs from pairing 1st–6th modes of ϕ _I − ϕ _R and ${\varphi }_I^{*}-{\varphi }_R$ . High correlation has been retained throughout the updating process for MAC and CoMAC.

The results of subtracting 1 − CoMAC and projecting onto the mode shapes before and after updating are compared in Figure 3, Subfigure B. It can be seen that correlation properties of DOFs in X-direction are mainly affected. Before the updating process, most uncorrelated nodes are distributed around a pole in the plate’s left centre half. Afterwards, nodal correlation properties evolve into a more axisymmetric shape and split into separated poles in the plate’s upper region. For both Y- and Z-direction, CoMAC values remain almost unchanged. The maximum value of 1 − CoMAC doesn’t exceed 1e-2 for X, Y- and Z-direction which indicates high nodal mode shape correspondence. Therefore, no severe decrease of nodal correlation is apparent.

Finally, mass and stiffness orthogonality for Φ _R to $M_I^{*}$ and $K_I^{*}$ as described in Section 3.4 must be preserved. The computed values of ${\varphi }_R^T{M}_I^{*}{\varphi }_R$ and ${\varphi }_R^T{K}_I^{*}{\varphi }_R$ for all considered modes are shown in Table 4 and Table 5. The latter was square rooted and divided by 2π to allow comparison with eigenfrequencies Λ _R in Table 3. In both cases, the updating procedure retained mass and stiffness orthogonality. The diagonal terms of the updated mass matrix $M_I^{*}$ consists of ones with near-zero off-diagonal terms in the range of |7.1e-8| to |1.1e-10|, indicating sufficient mass matrix orthogonality. The same applies to the updated stiffness matrix $K_I^{*}$ which is nearly equivalent to Λ _R with near-zero off-diagonal values of |5.1e-4| to |1.0e-5| for non-pairing modes.

Table 4.

Mass orthogonality evaluation.

Orthogonality of ${\varphi }_R^T{M}_I^{*}{\varphi }_R$
	ϕ R,1	ϕ R,2	ϕ R,3	ϕ R,4	ϕ R,5	ϕ R,6
ϕ R,1	1	−7.9e-9	7.1e-8	−5.4e-9	1.3e-9	−7.6e-10
ϕ R,2	−2.6e-8	1	1.1e-8	1.9e-10	1.1e-10	1.3e-10
ϕ R,3	1.5e-8	1.3e-9	1	1.6e-10	−1.4e-9	−6.3e-8
ϕ R,4	−2.6e-9	1.1e-9	1.4e-10	1	−2.4e-8	−6.9e-8
ϕ R,5	1.2e-9	6.9e-9	−7.8e-10	−2.7e-9	1	−2.2e-9
ϕ R,5	−8.6e-10	2.1e-10	−3.2e-10	−3.2e-8	1.1e-8	1

Table 5.

Stiffness orthogonality evaluation.

Orthogonality of ${\left({\varphi }_R^T{K}_I^{*}{\varphi }_R\right)}^{1/2}/2\pi$
	ϕ R,1	ϕ R,2	ϕ R,3	ϕ R,4	ϕ R,5	ϕ R,6
ϕ R,1	23.13	4.4e-4	8.3e-5	1.8e-4	1.5e-4	4.6e-4
ϕ R,2	1.1e-5	35.67	1.1e-4	1.4e-5	1.1e-4	9.2e-5
ϕ R,3	8.5e-5	1.2e-4	74.23	1.0e-5	9.9e-5	9.2e-5
ϕ R,4	1.9e-4	2.4e-4	9.1e-5	93.85	7.6e-5	1.2e-4
ϕ R,5	1.3e-4	1.4e-4	1.2e-4	8.1e-5	118.41	5.1e-4
ϕ R,5	4.2e-4	3.4e-5	9.5e-5	8.0e-5	1.6e-4	136.72

4.4. Question 2: Accuracy performance of correction terms

Special attention must be placed at the deviations of updated IM’s $M_I^{*}$ and $K_I^{*}$ compared to RM’s M _R and K _R in order to evaluate the credibility of the solution. Therefore, the percentage errors (ϵ_M,ij and ϵ_K,ij) are calculated to evaluate the local updating accuracy. As an example, high ϵ _ij (low accuracy) on isolated nodes surrounded by nodes with low ϵ _ij (high accuracy) can be interpreted as numerically induced abrupt changes in element’s mass and/or stiffness properties. This degrades the credibility of ΔM and ΔK but still might be compensated in the modal dynamic residual (11) by a sufficient amount of nodes with low ϵ _ij . Therefore, the absence of high singular errors should be verified since the updated IM’s dynamical behaviour will be dominated by that erroneous node. This investigation might not be possible in real-world scenarios if now RM is available. However, as this study investigates a numerical study, this supplementary analysis will prove that the updating does not generate artificial singularities within material characteristics. Obviously, this is the procedure’s fundamental requirement as model-based SHEM tool: Sudden changes in material properties must reliably indicate structural changes due to, for example, impact events and shall not be traced back to numerical artifacts.

In Figure 4, the percentage error distributions and their empirical cumulative distribution function (ECDF) for mass and stiffness are displayed. It can be seen that the majority of percentage errors remains in the lowest error bin of observed values. In fact, the number of low-order errors exceeds high error terms by several magnitudes of order, resulting in an ECDF ≥0.99 after the first histogram bin of ϵ _M = 1.75e-6 resp. ϵ _K = 2.5e-5. The maximum error does not exceed ϵ _M = 5e-5 resp. ϵ _K = 7 × 10⁻⁴. As a conclusion, the updating procedure determined accurate correction terms ΔK and ΔM.

Figure 4.

Percentage Error ϵ _ij evaluation from updated K, M. Subfigure (a) and Subfigure (b): Distribution of ϵ_K,ij and ϵ_M,ij for all non-zero entries. Additionally, the empirical cumulative distribution function is computed indicating low percentage errors ϵ _ij , that is, high updating accuracy.

4.5. Question 3: Transient structural response

Finally, the updated IM’s capability to reproduce time domain data is investigated. This is accomplished by applying a nodal impulse excitation F(t) function with force amplitude of $\widehat{F}=10$ N to the updated IM. The impulse excitation F(t) is defined in the following:

$\widehat{F}\cdot {\delta }_p\left(t\right)=\{\begin{array}{ll}\widehat{F},& \text{if} t=t_p\\ 0,& \text{else}\end{array}$

Both impact location and measurement location are depicted in Figure 1. The impact site is situated at the centre of the plate, 140 mm away from the unclamped top edge, while the nodal acceleration $\ddot{X}\left(t\right)$ response is recorded 140 mm below. To assess the impulse response, the free vibrations of IM, updated IM and RM are computed for a duration of t = 1 s with the impact event taking place at t _p = 0.1 s.

The acquired time domain data is shown in Figure 5, Subfigure A and Subfigure B. In upper Subfigure A, it can be seen that the dynamic response of the IM lacks from lower acceleration amplitudes compared to the RM. This was expected, since the IM was implemented with higher density and lower stiffness (see Table 2). The correlation of nodal acceleration at the measurement point is quantified by TRAC = 0.37. Therefore, the computed impulse response of IM and RM only corresponds to a minor extent. In lower Subfigure A, the nodal acceleration of the updated IM can be seen. The updated IM and RM reveal a high correlation in their dynamic response, resulting in TRAC = 0.79. It can be seen that the time domain correlation has improved significantly but not full correspondence of TRAC → 1 is achieved. This is due to the fact that only a small subset of six modes was incorporated into the updating procedure. However, this improvement is expected to be even higher when the included modal subset is expanded.

Figure 5.

Comparison of the transient dynamic response to an impulse excitation F(t) measured at the measurement point X(t). Subfigure (a): Nodal acceleration from IM, RM and updated IM* and RM. Subfigure (b): Power spectral density (PSD) of nodal acceleration from IM, RM and updated IM* and RM. For computing the PSD, only acceleration data after impulse excitation t > t _p is considered.

The computed power spectral densities (PSDs) before and after updating are displayed in Subfigure B. In both PSD plots, the relevant frequency spectrum up to the 6th mode is highlighted. As visualized in upper Subfigure B and discussed in Section 4.3, the eigenfrequencies – localized at the spectral power peaks – of the IM strongly differ from the RM. Additionally, the distributed spectral energy shows considerable differences at most of the eigenfrequencies. This corresponds to the observed lower acceleration amplitudes of the IM in Subfigure A. On the contrary, the comparison of PSD from updated IM against RM in lower Subfigure B reveals a high correlation in both eigenfrequencies and spectral power curve. Small deviations in the PSD can be observed at the anti-resonances. However, these discrepancies are of minor relevance for the generated structural response. Consequently, both time and frequency domain data of the updated IM show high correlation to the RM and can be utilized for model-based SHEM frameworks.

5. Conclusion

Structural health and event monitoring (SHEM) has proven maturity in a broad field of industrial applications to assess the structural integrity in the aftermath of singular events, environmental changes and material degradation. To allow model-based SHEM approaches, a high-fidelity structural model is necessary. This can be achieved by numerical model updating (NMU) techniques. In this study, a perturbation-based NMU procedure is investigated to reproduce a modal subset. Balancing the modal dynamic residual and incorporating orthogonality and symmetry constraints allows to derive a non-square linear equation system. Minimizing its residual by any suitable optimization method, for example, LSQR, allows determination of the correction terms. The approach allows including modes of arbitrary number in any frequency regime.

An FE model of an aluminium cantilever plate is utilized for evaluating the updating performance. The model is implemented twice with different material properties as reference model (RM) and initial model (IM). The IM is updated on behalf of the first six modes. The updating procedure determined highly accurate correction terms, allowing to reproduce the RM’s eigenfrequencies within an error margin of 2%. The mode shape correlation was evaluated by the Modal Assurance Criterion (MAC) and the Coordinate Modal Assurance Criterion (CoMAC). Both MAC and CoMAC proved high mode shape correlation between updated IM’s and RM’s modes. Lastly, a transient response to a nodal impulse excitation was computed for updated IM and RM. Time and frequency domain characteristics were evaluated by Test Response Assurance Criteria (TRAC) and Power Spectral Density (PSD). Again, the updated IM was able to reproduce the RM’s dynamic structural response in both time and frequency domain.

While the obtained results encourage to apply the updating procedure to real-world SHEM scenarios, further questions of more practical nature might arise which have not been covered in this study:

• The generation of spurious modes can’t be ruled out. This might happen if neighbouring updating modes are far separated, increasing the chance of spurious mode appearance in the frequency gap. Additionally, modes above the considered frequency spectrum might be altered as well. To suppress spurious mode generation, all included modes should be densely packed with regard to their eigenfrequencies. A subsequent correlation check, for example, by MAC to identify spurious modes seems reasonable.

• The updating formulation does not allow handling non-nominal constraints resp. boundary conditions as the related DOFs are fully removed from the updating procedure (see Section 3.5). Including those constraint DOFs can be considered fully unlocking of restrictions resp. supports and might enable physically unmeaningful solutions. Additionally, real-world constraints might be subject to more complex phenomena and cannot be sufficiently represented by altering local inertia and/or elastic properties of the structural model.

• As more complex structural models resp. geometries require varied element types, material and/or mesh partitioning, mass and stiffness matrices sparsity pattern might be affected and matrix entries with discrepancies in magnitude arise. The latter eventually impacts the condition number of matrix updating equations, causing both the eigenvalue problem and the updating formulation to become sensitive to perturbations in the input data and compromise numerical stability. As this becomes particularly critical when relying on experimentally gathered modal data under less controllable environments, inaccuracies from the measurement setup resp. sensor noise potentially amplify numerical instability of the updating procedure. This might lead to erroneous results that do not accurately represent the actual structural behaviour. To mitigate effects of ill-conditioning, regularization and preconditioning techniques should be employed to strengthen updating robustness which will be examined in upcoming studies.

• Ultimately, application to real-world data might be more burdensome as the extraction and projection of experimentally gathered mode shapes onto IM’s nodal DOFs requires dimensional reduction and/or mode shape interpolation techniques. The effects on updating robustness of such techniques must be evaluated in further studies.