# Model updating of locally non-linear systems based on multi-harmonic extended constitutive relation error

Inge Isasa, Aurélien Hot, Scott Cogan, Emeline Sadoulet-Reboul

► To cite this version:

Inge Isasa, Aurélien Hot, Scott Cogan, Emeline Sadoulet-Reboul. Model updating of locally non-linear systems based on multi-harmonic extended constitutive relation error. Mechanical Systems and Signal Processing, 2011, 25 (7), pp.2413-2425. 10.1016/j.ymssp.2011.03.010. hal-02300225

HAL Id: hal-02300225

https://hal.science/hal-02300225

Submitted on 22 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Model updating of locally non-linear systems
based on multi-harmonic extended constitutive relation error

I. Isasa<sup>a,*</sup>, A. Hot<sup>b</sup>, S. Cogan<sup>b</sup>, E. Sadoulet-Reboul<sup>b</sup>

<sup>a</sup> Ikerlan - IK4, Department of Mechanical Engineering, J.M. Arizmendiarreta 2, 20500 Arrasate-Mondragon, Spain

<sup>b</sup> FEMTO-ST Institute, Structural Dynamics Research Group, Rue de l'Epitaphe 24, 25000 Besançon, France

Improving the fidelity of numerical simulations using available test data is an important activity in the overall process of model verification and validation. While model updating or calibration of linear elastodynamic behaviors has been extensively studied for both academic and industrial applications over the past three decades, methodol-ogies capable of treating non-linear dynamics remain relatively immature. The authors propose a novel strategy for updating an important subclass of non-linear models characterized by globally linear stiffness and damping behaviors in the presence of local non-linear effects. The approach combines two well-known methods for structural dynamic analysis. The first is the multi-harmonic balance (MHB) method for solving the non-linear equations of motion of a mechanical system under periodic excitation. This approach has the advantage of being much faster than time domain integration procedures while allowing a wide range of non-linear effects to be taken into account. The second method is the extended constitutive relation error (ECRE) that has been used in the past for error localization and updating of linear elastody-namic models. The proposed updating strategy will be illustrated using academic examples.

# 1. Introduction

Non-linear phenomena are commonplace in mechanical systems containing mechanisms, joints and contact interfaces [1]. Engineers often simplify the behavior of complex structural models by considering them to be linear for dynamic analysis, thus neglecting non-linear effects due to large displacements, contact, clearance and impact phenomena, among others.

The following paper is devoted to the revision of non-linear models in the field of structural dynamics based on measured responses. During the past two decades, linear model updating has been extensively studied to improve the accuracy of simulations [2]. Non-linear model updating techniques on the other hand have received much less attention. Both time domain or frequency domain approaches can be found in the literature. In the time domain, the restoring force surface method (RFS) and proper orthogonal decomposition (POD) are described in detail in the overview paper by Kerschen et al. [3] with complete references to the literature. More recently, Gondhalekar et al. have

E-mail addresses: iisasa@orona.es (I. Isasa), aurelien.hot@univ-fcomte.fr (A. Hot), scott.cogan@univ-fcomte.fr (S. Cogan), emeline.sadoulet-reboul@univ-fcomte.fr (E. Sadoulet-Reboul).

<footnote>* Corresponding author. Tel.: +34 943 33 65 26; fax: +34 943 33 55 88.</footnote>
## Nomenclature

### Physical system

$M$ mass matrix

$K$ stiffness matrix

$C$ damping matrix

$p(t)$ vector of external forces in the time domain

$q(t)$ displacement time responses

$\dot{q}(t)$ velocity time responses

$\ddot{q}(t)$ acceleration time responses

$f_{NL}(q(t), \dot{q}(t))$ non-linear forces in the time domain

### Multi-harmonic system

$K$ multi-harmonic stiffness matrix

$K_R$ reduced multi-harmonic stiffness matrix

$Z(\omega)$ multi-harmonic dynamic stiffness matrix

$Q$ vector of harmonic coefficients

$Q_{\omega}^e$ vector of experimentally identified harmonic coefficients

$F(\omega)$ multi-harmonic non-linear forces

$P(\omega)$ multi-harmonic external force

$Q_{\omega}, V_{\omega}$ two admissible harmonic displacement fields

$m_j$ harmonic indices

$Q_0$ static term of the Fourier series

$Q_j^c$ jth cosine term of the Fourier series

$Q_j^s$ jth sine term of the Fourier series

$\omega$ angular frequency of the external harmonic excitation

$T$ period

$D$ measure of distance between two displacement fields

$E$ measure of the multi-harmonic modeling and test-analysis error

$r_{\omega}$ multi-harmonic residual displacement vector

$H$ transformation matrix between the numerical and experimental local reference frames

$\alpha$ scalar expressing the relative confidence in measurement data

$g$ objective function

$\gamma$ vector of Lagrange multipliers

### Other symbols

$N_p$ number of points per period

$N$ number of physical model degrees of freedom

$n$ number of harmonics

$N_e$ number of degrees of freedom measured

### Abbreviations

$ECRE$ extended constitutive relation error

$FE$ finite element

$MHB$ multi-harmonic balance

$POD$ proper orthogonal decomposition

$RFS$ restoring force surface

$dofs$ degrees of freedom

proposed a strategy combining the RFS method with model reduction [4]. In the frequency domain, Böswald and Link [5] have developed a methodology based on the first-order harmonic balanced method to obtain a suitable representation of non-linear effects and they have applied their approach to update non-linear joint parameters in complex structural assemblies. Another frequency domain method is investigated by Puel [6] where the extended constitutive relation error (ECRE) for linear dissipative systems is generalized to non-linear model updating with a first order harmonic balance approximation [7]. Non-linear identification methods can also be characterized as being either direct or modal [8–11]. Direct methods lead to an estimation of a physical stiffness or damping model property while modal methods project the non-linear effects into modal space. The latter represents a useful compromise when the structural responses on the non-linear dofs are not available. The identification methods referred to above as well as the method that will be proposed here belong to the class of direct methods. Direct methods have an important advantage over modal approaches in that they potentially provide information that is useful under different structural boundary and loading conditions.

With the exception of the RFS method, the existing methods for non-linear updating are based on some form of linearization and this naturally limits their application to relatively weak non-linear effects. As for the RFS approach, the major weakness of the formulation in physical coordinates (as opposed to the modal approach) lies in the fact that it requires that the structural responses be measured on all model degrees of freedom where significant non-linear effects are present.

In this paper, a novel methodology is presented which effectively combines the multi-harmonic balance method for calculating the periodic response of a non-linear system and the extended constitutive relation error method for establishing a well-behaved metric for modeling errors and test-analysis errors. The ECRE approach can be used as a benchmark criterion for comparing different candidate updated models on an absolute scale. Used as an error localization indicator, it only points to the subdomains of the model which are responsible for the test-analysis distances and gives no explicit clues as to what model characteristics (mesh refinement, element formulation, mechanical properties, etc.) are erroneous. In practice, this is determined by trial and error.

The proposed approach is neither based on any linear approximations nor does it require the observation of all non-linear degrees of freedom. The method presented in this paper lies in the group of direct methods, where multi degrees of freedom systems can be studied. An academic beam example with simulated experimental data will be used to illustrate the advantages and limitations of the methodology.
## 2. Mathematical formulation

### 2.1. Equations of motion

The equations of motion of a discrete linear structure can be written:

$$M\ddot{q}(t) + C\dot{q}(t) + Kq(t) = p(t) \quad (1)$$

where $K, M, C \in \mathbb{R}^{N,N}$ are, respectively, the symmetric stiffness, mass and damping matrices, with the stiffness matrix assumed to be non-negative definite; $p(t) \in \mathbb{R}^{N,1}$ is a vector of external forces; $q(t) \in \mathbb{R}^{N,1}$ is the vector of time responses on the $N$ degrees of freedom (dofs).

The equation of motion of a non-linear structure can be written in the same way as a linear structure with the addition of a non-linear term, $f_{NL}(q(t), \dot{q}(t)) \in \mathbb{R}^{N,1}$, which can depend on the system displacements and velocities:

$$M\ddot{q}(t) + C\dot{q}(t) + Kq(t) + f_{NL}(q(t), \dot{q}(t)) = p(t) \quad (2)$$

The origin of these non-linear forces can be quite diverse, including:

* Some large displacement systems for example the classical pendulum.

* Material non-linearities including locally plastic or viscoplastic behaviors, and so on.

* Local interface non-linearities including Hertz contact, dry friction, intermittent contact or clearance phenomena.

The response of a non-linear system can be qualitatively very different from a linear one. In a linear system the steady-state response to a periodic excitation is at the same frequency as the excitation force once the transient term vanishes in time and is independent of the initial conditions. The periodic response of a non-linear system, when it exists, generally exhibits primary and secondary resonances and can depend on the initial conditions [12].

Although transient behavior may be important, the study of periodic solutions and their stability remains essential to capturing the behavior of a vibrating system. Non-linear time domain simulations are extremely burdensome especially when they are used to calculate the steady-state response of large-order models. The multi-harmonic balance method is based on a Fourier series approximation and was developed with the objective of solving for the periodic response of non-linear systems more efficiently.

### 2.2. Multi-harmonic balance method

The MHB method is a frequency domain approach developed to solve equation (2) for a periodic excitation. Many extensions to the first-order harmonic balance approach to include higher harmonics were developed in the 1980s, for example [13] or [14]. The developments presented in this article are based on the formulation proposed by Cardona et al. [15] and more recently applied to complex industrial structures with contact effects by Petrov et al., for example [16].

The equilibrium equation of a non-linear system of $N$ degrees of freedom is given by Eq. (2). Expressing the vector of time responses $q(t)$ as a Fourier series yields:

$$q(t) = Q_0 + \sum_{j=1}^{n} (Q_j^c \cos m_j \omega t + Q_j^s \sin m_j \omega t) \quad (3)$$

where

* $Q_0$ represents the constant or static contribution;

* $Q_j^c$ and $Q_j^s$ are, respectively, the $j$th cosine and sine coefficients of the Fourier series;

* $m_j$ expresses the harmonic of the excitation frequency $\omega$.

Introducing this expression into Eq. (2) yields:

$$K \left( Q_0 + \sum_{j=1}^{n} Q_j^c \cos m_j \omega t + Q_j^s \sin m_j \omega t \right) + C \left( \sum_{j=1}^{n} -m_j \omega Q_j^c \sin m_j \omega t + m_j \omega Q_j^s \cos m_j \omega t \right) + M \left( \sum_{j=1}^{n} -(m_j \omega)^2 Q_j^c \cos m_j \omega t - (m_j \omega)^2 Q_j^s \sin m_j \omega t \right) + f_{NL}(q(t), \dot{q}(t)) - p(t) = 0 \quad (4)$$

Then by sequentially pre-multiplying Eq. (4) by the harmonic functions $(1, \cos m_1 \omega t, \sin m_1 \omega t, \dots, \cos m_n \omega t, \sin m_n \omega t)$ and integrating over the period $T = 2\pi/\omega$ and regrouping the resulting equations for each harmonic in the Fourier expansion, the following frequency domain expression can be obtained:

$$Z(\omega)Q + F(Q, \omega) - P = 0 \quad (5)$$
where $Q = \{Q_0, Q_1, Q_2, \dots, Q_{2n-1}, Q_{2n}\}$ is the vector of harmonic coefficients with $Q_i \in \mathbb{R}^{N,1}$. The matrix $\mathcal{Z} \in \mathbb{R}^{(2n+1)N, (2n+1)N}$ is given by

$$ \mathcal{Z} = \begin{bmatrix} K & 0 & 0 & \dots & 0 & 0 \\ 0 & K-(m_1\omega)^2M & m_1\omega C & \dots & 0 & 0 \\ 0 & -m_1\omega C & K-(m_1\omega)^2M & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & K-(m_n\omega)^2M & m_n\omega C \\ 0 & 0 & 0 & \dots & -m_n\omega C & K-(m_n\omega)^2M \end{bmatrix} \quad (6) $$

and the vectors $\mathcal{F}, \mathcal{P} \in \mathbb{R}^{(2n+1)N,1}$ corresponding, respectively, to the non-linear forces and the external excitations are given by

$$ \mathcal{F} = \left\{ \begin{array}{l} \int_0^T f_{NL}(q(t), \dot{q}(t)) \, dt \\ \frac{\omega}{\pi} \int_0^T f_{NL}(q(t), \dot{q}(t)) \cos \omega t \, dt \\ \frac{\omega}{\pi} \int_0^T f_{NL}(q(t), \dot{q}(t)) \sin \omega t \, dt \\ \vdots \\ \frac{\omega}{\pi} \int_0^T f_{NL}(q(t), \dot{q}(t)) \cos m_n \omega t \, dt \\ \frac{\omega}{\pi} \int_0^T f_{NL}(q(t), \dot{q}(t)) \sin m_n \omega t \, dt \end{array} \right\} \quad (7) $$

and

$$ \mathcal{P} = \left\{ \begin{array}{l} \int_0^T p(t) \, dt \\ \frac{\omega}{\pi} \int_0^T p(t) \cos \omega t \, dt \\ \frac{\omega}{\pi} \int_0^T p(t) \sin \omega t \, dt \\ \vdots \\ \frac{\omega}{\pi} \int_0^T p(t) \cos m_n \omega t \, dt \\ \frac{\omega}{\pi} \int_0^T p(t) \sin m_n \omega t \, dt \end{array} \right\} \quad (8) $$

The following remarks apply to the aforementioned equations:

* The Fourier series expansion that enables the transformation from non-linear time domain into linearized frequency domain can be observed in Eqs. (7) and (8). It can be seen that each harmonic of the periodic excitation yields a corresponding sine- and cosine-amplitude for the excitation for $\mathcal{P}$ but also for the non-linear force $\mathcal{F}$.

* Eq. (5) is still non-linear, even after linearization by Fourier series expansion that comes along with MHB. Thus, Eq. (5) must be solved iteratively for the unknown terms $Q$ and $\mathcal{F}(Q, \omega)$. A Newton iteration scheme can be used for this. In [17], a predictor-corrector continuation scheme has been applied to solve this equation.

* Model reduction (e.g. Guyan reduction) can be used effectively for the linear system matrices in order to reduce the computational burden for very large models.

* The number of harmonics taken increases the size of the system to be solved and thus increases time for numerical calculations.

## 2.3. Extended constitutive relation error

The constitutive relation error was initially proposed by Ladevèze et al. in the early 1980s as an error estimator for finite element models [18]. An extended version for use in model updating was described mid of the 1990s, for example [19], taking into account both modeling error and test-analysis errors for linear elastodynamic behaviors. An example of a discrete formulation of the approach for dissipative linear structures can be found in [20]. The basic philosophy of the ECRE methodology consists in dividing the relations of interest (constitutive behavior laws, equations of motion, measured displacements, initial conditions, etc.) into two groups: the reliable and the less reliable quantities. The solution to the problem is sought so as to satisfy the reliable equations exactly while minimizing the errors in the less reliable equations. For the sake of simplicity and better interpretation of the method, the following equations are derived for non-linear elastodynamic systems which contain only non-linear stiffness errors. Extensions to non-linear dissipative effects as well as combined errors in both linear and non-linear properties can be formulated in an analogous manner.
Let $Q_\omega$ and $V_\omega$ be two admissible displacement fields of Eq. (5) and $D^2(Q_\omega, V_\omega)$ a measure of distance between these two vectors such that

$$D^2(Q_\omega, V_\omega) = \|Q_\omega - V_\omega\|_K^2 = (Q_\omega - V_\omega)^T K (Q_\omega - V_\omega) \quad (9)$$

where $K \in \mathbb{R}^{(2n+1)N, (2n+1)N}$ is the multi-harmonic stiffness matrix corresponding to the linear system defined by

$$K = \begin{bmatrix} K & 0 & 0 & \cdots & 0 & 0 \\ 0 & K & 0 & \cdots & 0 & 0 \\ 0 & 0 & K & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & K & 0 \\ 0 & 0 & 0 & \cdots & 0 & K \end{bmatrix} \quad (10)$$

A multi-harmonic ECRE can be defined for non-linear stiffness errors in the following way:

$$E_{\omega}^{2} = r_{\omega}^{T} K r_{\omega} + \alpha (H Q_{\omega} - Q_{\omega}^{e})^{T} K_{R} (H Q_{\omega} - Q_{\omega}^{e}) \quad (11)$$

where

* $r_\omega = Q_\omega - V_\omega$ with $Q_\omega \in \mathbb{R}^{(2n+1)N, 1}$ and $V_\omega \in \mathbb{R}^{(2n+1)N, 1}$ two admissible displacement fields for multi-harmonic equation of motion, Eq. (5). Admissible vectors are displacements fields which satisfy the boundary conditions and geometric constraints of the system. In the present case, the vector $Q_\omega$ represents the experimentally identified harmonic coefficients expanded to either the total number of model dofs or to the total number of reduced model dofs if a preliminary model reduction was employed.

* $Q_\omega^e$ is the vector of harmonic coefficients in the Fourier series expansion of the experimental response. They are obtained directly from the experimentally observed time responses via the fast Fourier transform [21]:

$$Q_j^e = \frac{1}{N_p} \sum_{k=0}^{N_p-1} q(k) \cos \left( \frac{2\pi k j}{N_p} \right) \quad (12)$$

$$Q_j^s = \frac{1}{N_p} \sum_{k=0}^{N_p-1} q(k) \sin \left( \frac{2\pi k j}{N_p} \right) \quad (13)$$

where $N_p$ is the number of points per period.

* $H \in \mathbb{R}^{(2n+1)N, N}$ is a projection matrix allowing the model responses $Q_\omega$ to be projected onto the set of $N_e$ measurement directions so as to account for the limited number of measurement dofs and any differences in local reference frames between the FE model and the experimental model.

* $K_R \in \mathbb{R}^{(2n+1)N, (2n+1)N}$ is the multi-harmonic stiffness matrix of the linear system reduced to the measurement degrees of freedom. In practice, the Guyan stiffness matrix is generally used.

* $\alpha$ is a real positive scalar allowing the relative confidence in the identified harmonic coefficients to be taken into account.

Eq. (11) is composed of two terms. The first term is a measure of the modeling error whereas the second term is a measure of the distance between the experimentally identified harmonic coefficients and those predicted by the model. Both of these terms correspond to the less reliable quantities in the present ECRE formulation. The reliable quantities correspond to the equilibrium equations of the system expressed by Eq. (5). The ECRE approach can be used as a benchmark criterion for comparing different candidate updated models on an absolute scale and thus can be used in some cases (if the differences are significant) to indicate if it is the Young's modulus that is in error or, for example, a problem of mesh refinement. In Eq. (11) $r_\omega$ and $Q_\omega$ are the results of the analysis. The residual vector, $r_\omega$, will be used to evaluate the modeling error, whereas $Q_\omega$ will be used to evaluate test-analysis error. The remaining terms of Eq. (11) are known.

Therefore, in order to evaluate $r_\omega$ and $Q_\omega$ in this case, the following minimization problem has to be solved:

$$\begin{cases} \text{Minimize} & E_{\omega}^{2} = r_{\omega}^{T} K r_{\omega} + \alpha \|H Q_{\omega} - Q_{\omega}^{e}\|_{K_{R}}^{2} \\ \text{Under the constraint} & K r_{\omega} = Z(\omega) Q_{\omega} + \mathcal{F} - P \end{cases} \quad (14)$$

Combining the minimization problem and the constraint in a single minimization problem yields

$$\min g = r_{\omega}^{T} K r_{\omega} + \alpha(H Q_{\omega} - Q_{\omega}^{e})^{T} K_{R} (H Q_{\omega} - Q_{\omega}^{e}) + \gamma^{T} (K r_{\omega} - Z(\omega) Q_{\omega} - \mathcal{F} + P) \quad (15)$$

where $g$ is the objective function and $\gamma \in \mathbb{R}^{(2n+1)N, 1}$ is a vector of Lagrange multipliers.

The stationarity conditions require:

$$\frac{\partial g}{\partial r_{\omega}} = 0 \Rightarrow K(2r_{\omega} + \gamma) = 0$$
$$ \frac{\partial g}{\partial Q_{\omega}} = 0 \Rightarrow -Z(\omega)\gamma - \frac{\partial F}{\partial Q_{\omega}}\gamma + 2\alpha H^T K_R (H Q_{\omega} - Q_{\omega}^e) = 0 $$

$$ \frac{\partial g}{\partial \gamma} = 0 \Rightarrow \kappa r_{\omega} - Z(\omega) Q_{\omega} - F + P = 0 \qquad (16) $$

Eliminating $\gamma$ and regrouping the equations yields the following non-linear matrix equation:

$$ \begin{bmatrix} Z(\omega) + \frac{\partial F}{\partial Q_{\omega}} & \alpha H^T K_R H \\ K & -Z(\omega) \end{bmatrix} \begin{Bmatrix} r_{\omega} \\ Q_{\omega} \end{Bmatrix} + \begin{Bmatrix} 0 \\ -F \end{Bmatrix} = \begin{Bmatrix} \alpha H^T K_R Q_{\omega}^e \\ -P \end{Bmatrix} \qquad (17) $$

The following remarks apply to the aforementioned equations:

* Eq. (17) requires the solution of a non-linear system of order $2N(2n+1)$. It can be solved with a classical Newton-Raphson iterative procedure.

* The solution of Eq. (17) comprises two unknown vectors. First, the residual vector $r_{\omega}$ represents the displacement field resulting from the unbalanced forces in the multi-harmonic equations of motion and provides the basis for calculating the modeling error. Second, the response vector $Q_{\omega}$ represents the experimental multi-harmonic response expanded onto all model dofs and provides a means for evaluating the test-analysis distances.

* Given the vectors $r_{\omega}$ and $Q_{\omega}$, the total MBH-ECRE error equation (11) for the point in model space defined by the nominal linear system matrices and the nominal non-linear model used to estimate the multi-harmonic non-linear forces can now be evaluated.

* The model updating problem simply consists in minimizing the total MBH-ECRE error over the space defined by coefficients of the non-linear model.

* Uncertainty in the experimentally identified quantities can be taken into account via a weighting factor. The weighting factor $\alpha$ provides a means of expressing the analysts relative confidence between the measurement data and the nominal finite element model. In practice, $\alpha$ is chosen in such a way as to obtain a *posteriori* the test-analysis errors that are anticipated in the measured displacements as a function of the testing conditions.

### 3. Illustrative academic example

The proposed methodology will be illustrated on a simulated academic example based on the COST action F3 project benchmark [22]. The model consists in a clamped linear beam attached to a thinner beam at one end. The main beam has a length of 0.7 m and a thickness of 0.014 m, whereas the thin beam has a length of 0.04 m with a thickness of 0.0005 m. Both beams have a width of 0.014 m and the material of both of them is steel with a Young's modulus of 210 GPa and a Poisson ratio of 0.3.

The structure is excited at node number 3 (see Fig. 1) with a frequency increasing stepped sine excitation having an amplitude of 2 N. This amplitude level was chosen based on the results of [3] in order to insure a large enough deflection for non-linear effects to come into play.

As stated in [3,23], the non-linear behavior appears mainly in the first mode (30.76 Hz). A grounded cubic stiffness is introduced at node number 8 as an equivalent modeling to the true geometric stiffness non-linearity. Node number 13 is placed in the same location as node 8. Both nodes are connected by a rotational stiffness such that both nodes got the same translation displacement but rotations are different. Moreover, in this example the influence of this cubic non-linearity is studied for the first mode, for different harmonics and for different degrees on non-linearity. The nominal value of the non-linear coefficient was chosen to be $6.1 \times 10^9$ N/m [23].

The FRF is calculated between the excitation point (node 3) and the response point (node 8) and plotted in Fig. 2(a) in order to visualize the distortion resulting from the non-linear effects. Fig. 2(b) displays the FFT of the time response of

<figure><img src="image_2.png" /><figcaption>Fig. 1. CostF3 beam.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Frequency domain responses of the beam to a constant force of 2 N. (a) The upper branch of the FRF between the excitation node 3 and the response node 8. (b) FFT of node 8 time response for a 32 Hz excitation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Upper branch of the FRF between the excitation point (node 3) and the response point (node 8) with the studied three frequencies.</figcaption></figure>

node 8 to a 32 Hz sine excitation. A peak at the fundamental frequency is observed as well as the 3rd (96 Hz) and the 5th (160 Hz) harmonics. The 7th harmonic (224 Hz) is also present but barely visible.

In order to apply the non-linear model updating procedure based on the MBH-ECRE approach the values of the experimental vector of harmonic coefficients contained in $Q_\omega^e$ are needed. In this paper they are obtained numerically using a Newmark algorithm developed based on [24] followed by a FFT analysis. In cases where a constant excitation frequency is used a Fourier series expansion or harmonic curve fitting can be used.

In order to illustrate the advantages and limitations of the proposed non-linear updating strategy, it will be applied in three different simulated test configurations. Three excitation frequencies will be investigated corresponding to different response levels and thus different degrees of non-linearity, see Fig. 3. The objective here is simply to examine the shape of the error expressed by Eq. (11) as a function of a single non-linear model parameter. To simplify the interpretation of the results, the experimental harmonic coefficients have been generated based on the nominal non-linear model. As such, in what follows it is expected to see a minimum in the MHB-ECRE curve at a value of the correction coefficient that multiplies the non-linear stiffness ($K_{NL}$) equal to 1. The number of harmonics taken into account in the MHB-ECRE procedure will also be studied here.

## 3.1. Test case 1: verification of the implemented algorithm

The objective of this first configuration is simply to verify the implemented MHB-ECRE algorithm. In this case, all model degrees of freedom (dofs) have been measured, that is to say, all 21 dofs (10 translations and 11 rotations) of the beam in
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Complete model and 21 dofs measured: MHB-ECRE results. —, 28 Hz; --, 30 Hz; . ., 32 Hz. (a) Fundamental + H3 + H5 + H7, (b) fundamental + H3 + H5, (c) fundamental + H3 and (d) fundamental only.</figcaption></figure>

Fig. 1 are observed. Fig. 4 plots the results of the MHB-ECRE updating procedure. In the case where the fundamental, the 3rd, 5th and 7th harmonics are taken into account, in Fig. 4(a), the MHB-ECRE curves are, as expected, minimum for a correction coefficient equal to 1. In Fig. 4(b) and (c), where the 7th and 7th + 5th harmonics are, respectively, removed from the MHB-ECRE calculation, the results still give a good estimation of the non-linear parameter. In the case where only the fundamental contribution is retained, in Fig. 4(d), the procedure is still accurate for frequencies 32 and 30 Hz, whereas for 28 Hz the minimum value is slightly overestimated. One can explain this slight shift by the fact that, at this last frequency, the response amplitude is lower than at the two others and thus the non-linear behavior is less influent. However, in this case, this lack of non-linear information contained in the fundamental term is counterbalanced by the one contained the third harmonic. It can also be noted that all three curves are convex, which is an advantage in finding the minimum of the MHB-ECRE function.

Fig. 5 shows the two terms of Eq. (11) (modeling error and test-analysis error) as a function of the non-linear correction coefficient for a frequency of 32 Hz and a weighting factor $\alpha = 0.5$. In practice, the relative values of these two terms will depend on the pertinence of the updating parameters, the level of measurement noise and the selected value of $\alpha$. In the present case, the relatively small value of the test-analysis distances is due to the absence of measurement noise and the fact that the correct model parameter was chosen for the study.

## 3.2. Test case 2: impact of a reduced set of measurement dofs

The second case aims at illustrating the impact of observing only a subset of a model dofs. In the present case, only four translations are assumed to be measured corresponding to nodes 3, 4, 6 and 8. Model reduction has been performed based
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Modeling (—) and test-analysis (---) errors plotted separately as a function of the non-linear correction coefficient (32 Hz).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Four dofs reduced model and four dofs measured: MHB-ECRE results. —, 28 Hz; ---, 30 Hz; . ., 32 Hz. (a) Fundamental only, (b) fundamental + H3 + H5 + H7.</figcaption></figure>

on the static Guyan procedure [25]. It is important to highlight that the non-linear dofs must be retained as master dofs when reducing the model even if they have not been observed experimentally. Indeed, they are required in order for the unknown non-linear forces to appear explicitly in the equation. Finally, the model dofs corresponding to the non-linear cubic spring is assumed to be included in the set of observed dofs. The comparison between the results of the MHB-ECRE for the three different excitation frequencies, taking into account the fundamental contribution only and the fundamental plus the three first odd harmonics, is plotted in Fig. 6(a) and (b), respectively. The curves are still convex in both cases. However, in Fig. 6(b), the non-linear parameters are now underestimated although the MHB-ECRE curve for the highest excitation frequency (corresponding to the highest response amplitude and thus the largest non-linear effect) still has a minimum in the vicinity of 1. These shifts are due to the fact that the Guyan reduction is no longer an exact representation of the dynamics of the linear system. However, as the non-linear effects increase in magnitude, this discrepancy becomes less and less important. Moreover, a compensation effect between errors due to model reduction and the loss of information due to harmonic truncation can be observed in the results at 28 Hz. Indeed, the model reduction tends to shift the minimum to the left while the harmonic truncation tends to shift the minimum to the right.

The same study is now carried out on a reduced two dofs system. The Guyan reduction technique is still used and the two master dofs are the translations at nodes 3 and 8. The results plotted in Fig. 7 lead to an even larger shift in the minimums of the curves. It is thus important to understand the shift in the MHB-ECRE curves observed between a complete model and a reduced model. Indeed, an indicator is clearly required to quantify the relative impacts of model reduction errors and non-linear effects.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Two dofs reduced model and two dofs measured: MHB-ECRE results. —, 28 Hz; - - -, 30 Hz; . . ., 32 Hz. (a) Fundamental only, (b) fundamental + H3 + H5 + H7.</figcaption></figure>

<table><caption>**Table 1**
Eigenfrequencies and relative errors of the different studied models.</caption><thead><tr><th></th><th>Complete model (Hz)</th><th>Four dofs reduced model, relative error</th><th>Two dofs reduced model, relative error</th></tr></thead><tbody><tr><td>Mode 1</td><td>30.76</td><td>30.74 Hz, 0.02%</td><td>31.08 Hz, 1%</td></tr><tr><td>Mode 2</td><td>150.62</td><td>151.35 Hz, 0.5%</td><td>161.72 Hz, 7.4%</td></tr></tbody></table>

One way to quantify the errors due to model reduction is to check the eigenfrequencies accuracy. The eigenfrequencies of the three different models are summarized in Table 1. Small differences can be noted for the four dofs reduced model: 0.02% and 0.5% of relative error for, respectively, the first and second modes. However, for the two dofs reduced model, the errors are more important: 1% for the first mode and 7.4% for the second. Moreover, the 5th harmonic of the three excitation frequencies are 140, 150 and 160 Hz, which is close to the second eigenfrequency. The poor accuracy of this reduced model for the first mode and even more on the second mode may explain the large shift observed in Fig. 7.

Another way to understand the errors due to model reduction is to quantify the ratio between the residual error $\mathcal{R}$ of the reduced equilibrium equation (18) and the effective non-linear force $\mathcal{F}$:

$$ \mathcal{R} = \mathcal{Z}^r(\omega)\mathcal{Q}^r + \mathcal{F}^r(\mathcal{Q}^r) - \mathcal{P}^r \qquad (18) $$

where $\mathcal{Z}^r \in \mathbb{R}^{(2n+1)N_r,(2n+1)N_r}$ and $\mathcal{R}, \mathcal{Q}^r, \mathcal{F}^r, \mathcal{P}^r \in \mathbb{R}^{(2n+1)N_r,1}$ and $N_r$ is the number of dofs retained in the reduced model. Fig. 8 plots the ratio $|\mathcal{R}|/|\mathcal{F}^r|$ calculated for the three frequencies, for the nominal value of the non-linear parameter and taking into account all the harmonics. These results show that for a two dofs reduced model, the non-linear information included in $\mathcal{F}$ is the same order of magnitude as the residual error and thus is not sufficient to have a good estimation of the non-linear parameter. However, for a four dofs reduced model, and even more with 21 dofs (complete model), the ratio tends to 0. The non-linear force is now more significant and, as shown in test cases 1 and 2, the non-linear coefficient can be effectively estimated.

In this case, the two dofs static reduction is clearly insufficient and either more master dofs must be included or an alternative model reduction technique must be used that takes into account the dynamic behavior of the slave structure, such as the Craig-Bampton [26] or Petersmann [27] methods.

In order to improve the previous results, the matrices are now reduced using the Petersmann dynamic condensation method with the same number of master dofs. Fig. 9 plots the comparison between results using a Petersmann's dynamic reduction method and the Guyan's static one's for a four dofs and a two dofs reduction. The results clearly show that the dynamic reduction method is, in this case, much more accurate. For a four dofs reduced model, Fig. 9(a), the non-linear parameter is precisely estimated and for a two dofs reduced model, Fig. 9(b), the results are more accurate. These results are actually in agreement with the previous remarks.

### 3.3. Test case 3: impact of a lack of measurements on non-linear dofs

This third case aims at illustrating a very important characteristic of the proposed updating strategy, namely that it is not necessary to experimentally observe the model degrees of freedom corresponding to the location of the non-linear
<figure><img src="image_9.png" /><figcaption>Fig. 8. Influence of model reduction on results accuracy.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between Guyan static reduction and Petersmann dynamic reduction. ——, 28 Hz Guyan; ---, 30 Hz Guyan; ×, 32 Hz Guyan; —, 28 Hz Petersmann; —, 30 Hz Petersmann; ×, 32 Hz Petersmann.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Four dofs reduced model and three dofs measured: MHB-ECRE results. —, 28 Hz; —, 30 Hz; ., 32 Hz. (a) Fundamental only, (b) fundamental + H3 + H5 + H7.</figcaption></figure>
physics (the translational displacement at node 8 in this example). In this test case, the four dofs reduced system using Petersmann method is studied, but it is assumed that the displacement of node 8 is no longer available. Results taking into account only the fundamental and all the harmonics are, respectively, plotted in Fig. 10(a) and (b). They are qualitatively similar to those of test case 1: no parameter estimation error if all harmonics are included in the MHB-ECRE procedure and a slight shift for the low level response, at 28 Hz, if only the fundamental is retained. That means that the measurement of the non-linear dofs is not required for an accurate estimation of the non-linear coefficient.

## 4. Discussion

The proposed MBH-ECRE approach provides a valuable measure of the distance between the simulated and experimental identified responses of a locally non-linear elastodynamic structure and includes the weighted contributions of two terms: an implicit modeling error and an explicit model response error on the measured dofs.

The potential advantages of this methodology with respect to existing non-linear identification methods are

* No locally linearized model is assumed thus allowing strongly non-linear systems to be analyzed by including the dominant sub- and super-harmonic contributions.

* The approach does not implicitly require to experimentally observe the structural displacements at the locations of the non-linear physics. From a practical point of view, this is clearly very important since non-linear effects often occur at inaccessible locations and, more generally speaking, the rotational dofs are very difficult to observe experimentally.

* The form of the MBH-ECRE objective function is globally convex, thus facilitating the use of efficient local optimization algorithms in an automated updating process.

* Model reduction techniques can be used very effectively to reduce the associated linear system model to the measured dofs thus reducing calculation costs.

* The model responses do not need to be re-evaluated at every updating iteration as time domain responses obtained from the numerical integration of the non-linear equation of motion. For a given set of updating parameters, the objective function is minimized based on iteratively solving Eq. (17). It should be noted that $F$ in Eq. (17) depends on the response amplitudes of the non-linear system. Thus, the computational burden of the method is dependent on the number of harmonics considered but is still less in comparison to performing a numerical integration in the time domain.

* Uncertainty in the experimentally identified quantities can be taken into account via a weighting factor that expresses the degree of confidence in the test data.

## 5. Conclusions

This paper presents a novel non-linear model updating approach that combines two well-known strategies for structural dynamic analysis, namely the multi-harmonic balance method for calculating the periodic response of a non-linear system and the extended constitutive relation error method for establishing a well-behaved metric for modeling and test-analysis errors. The proposed updating strategy has been illustrated using simulated data based on the COST-F3 beam benchmark to evaluate the impact of a reduced set of measurement dofs and to study the influence of uncertainty in the identified experimental quantities.

Future work will seek to develop a decision-making indicator to weigh the impact of the different sources of uncertainty (e.g. measurement noise, lack-of-knowledge in the associated linear model, form of non-linearity) on the capacity of the MBH-ECRE algorithm to correctly localize and update the unknown non-linear coefficients. The methodology will also be applied to a more complex model based on measurement data.

## Acknowledgments

The work presented in this paper has been carried out with the generous support of Orona EIC (Hernani, Spain), the Centre National des Etudes Spatiales (Toulouse, France), and Thales Alenia Space (Cannes La Bocca, France).

## References

[1] A.H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods, Wiley Series in Nonlinear Science, NY, USA, 1995.

[2] J. Mottershead, M. Friswell, Model updating in structural dynamics: a survey, Journal of Sound and Vibration 167 (2) (1993) 347-375.

[3] G. Kerschen, K. Worden, A.F. Vakakis, J. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mechanical Systems and Signal Processing 20 (3) (2006) 505-592.

[4] A.C. Gondhalekar, E.P. Petrov, M. Imregun, Parameter identification in nonlinear dynamic systems via genetic algorithm optimization, Journal of Computational Nonlinear Dynamics 4 (4) doi:10.1115/1.3187213.

[5] M. Böswald, M. Link, Identification of non-linear joint parameters by using frequency response residuals, in: Proc. of the International Modal Analysis Conference IMAC-XXIII, Orlando, Florida, 2005.
[6] G. Puel, Mise en évidence et recalage des non-linéarités locales en dynamique des structures, Master's Thesis, Université Paris VI, France, 2001.

[7] P. Ladevèze, D. Nedjar, M. Reynier, Updating of finite element models using vibration test, AIAA Journal 32 (7) (1994) 1485-1491.

[8] S. Masri, T. Caughey, A nonparametric identification technique for nonlinear dynamic problems, Journal of Applied Mechanics 46 (1979) 433-447.

[9] E. Crawley, A. Aubert, Identification of nonlinear structural elements by force-state mapping, AIAA Journal 24 (1) (1986) 155-162.

[10] M. Platten, J. Wright, J. Cooper, M. Sarmast, Identification of multi-degree of freedom non-linear systems using an extension of force appropriation, in: Proc. of IMAC XX, Los Angeles, USA, 2002.

[11] J. Wright, M. Platten, J. Cooper, M. Sarmast, Identification of multi-degree of freedom weakly non-linear systems using a model based in modal space, in: Proc. of International Conference on Structural System Identification, Kassel, Germany, 2001.

[12] K. Worden, G. Tomlinson, Nonlinearity in Structural Dynamics: Detection, Identification, and Modeling, Institute of Physics Publishing (IoP), Bristol, Philadelphia, 2000.

[13] C. Pierre, A. Ferri, E.H. Dowell, Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method, ASME Journal of Applied Mechanics 52 (4) (1985) 958-964.

[14] T. Cameron, J. Griffin, An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems, ASME Journal of Applied Mechanics 56 (1) (1989) 149-154.

[15] A. Cardona, T. Coune, A. Lerusse, M. Geradin, A multiharmonic method for non-linear vibration analysis, International Journal for Numerical Methods in Engineering 37 (1994) 1593-1608.

[16] E. Petrov, D. Ewins, Effects of damping and varying contact area at blade-disk joints in forced response analysis of bladed disk assemblies, Journal of Turbomachinery 128 (2006) 403-410.

[17] J. Ferreira, A. Serpa, Application of the arc-length method in nonlinear frequency response, Journal of Sound and Vibration 284 (2005) 133-149.

[18] P. Ladevèze, D. Leguillon, Error estimate procedure in the finite element method and application, SIAM Journal of Numerical Analysis 20 (3) (1983) 485-509.

[19] P. Ladevèze, M. Reynier, N. Maia, Error on the constitutive relation in dynamics: theory and application for model updating, in: Inverse Problems in Engineering, Balkema, Rotterdam 1994, pp. 251-256.

[20] A. Deraemaeker, P. Ladevèze, P. Leconte, Reduced bases for model updating in structural dynamics based on constitutive relation error, Computer Methods in Applied Mechanics and Engineering 191 (2002) 2427-2444.

[21] M. Berthillier, C. Dupont, R. Mondal, J. Barrau, Blades forced response analysis with friction dampers, ASME Journal of Vibration and Acoustics 120 (1998) 468-474.

[22] F. Thouverez, Presentation of the ECL benchmark, Mechanical Systems and Signal Processing 17 (1) (2003) 195-202.

[23] V. Lenaerts, G. Kerschen, J. Golinval, ECL benchmark: application of the proper orthogonal decomposition, Mechanical Systems and Signal Processing 17 (1) (2003) 237-242.

[24] M. Géradin, D. Rixen, Mechanical Vibrations: Theory and Applications to Structural Dynamics, second ed., John Wiley & Sons, 1997.

[25] R. Guyan, Reduction of stiffness and mass matrices, AIAA Journal 3 (2) (1965) 380.

[26] R. Craig, M. Bampton, Coupling of substructures for dynamics analysis, AIAA Journal 6 (7) (1968) 1313-1319.

[27] N. Petersmann, Calculation of eigenvalues using substructures and dynamic condensation, in: Proc. of the Second International Conference on Recent Advances in Structural Dynamics, Southampton, UK, 1984.