Updating of a Nonlinear Finite Element Model Using Discrete-Time Volterra Series

Abstract

In this study, the discrete-time Volterra series are used to update parameters in a nonlinear finite element model. The main idea of the Volterra series is to describe the discrete-time output of a nonlinear system using multidimensional convolutions between the Volterra kernels represented in a Kautz orthogonal basis and the excitations. A metric based on the residue between the experimental and the numerical Volterra kernels is used to identify the parameters of the numerical model. First, the identification of the linear parameters is performed using a metric based only on the first order Volterra kernels. Then the nonlinear parameters are identified through a metric based on the higher-order kernels. The originality of this nonlinear updating method stems from the decoupling of linear and nonlinear parameters and the use of global nonlinear model. In order to put in light the applicability of this technique, this work focus on the identification of the parameters in a nonlinear finite element model of a beam that was preloaded by compression mechanism. This work shows that the updated numerical model was able to represent the behaviour observed in the experimental measurements.

Keywords:
Discrete-time Volterra series; Kautz filters; model updating; nonlinear finite element model; nonlinear identification; sensitivity

1 INTRODUCTION

In many cases, the dynamical behaviour is needed to design a mechanical system. Moreover, the knowledge of the dynamical behaviour of a structure is important to control its integrity and to predict its service life. In the practical systems, nonlinear effects due to large displacements, gaps, jumps, discontinuities, etc. are very common (e.g. (Worden and Tomlinson, 2001)). The nonlinear finite element method is a very powerful tool to study the nonlinear vibrations. Nevertheless, it is not easy to identify the parameters of the nonlinear model due to the uncertainties in the modelling of the experimental system (e.g. the mechanical properties, the boundary conditions). Despite of the important number of investigations, the nonlinear model updating techniques are not mature as the classical linear ones.

A bibliography review on the state of the art nonlinear system identification techniques in structure dynamics is presented by (Noël and Kerschen, 2017). Furthermore, a bibliography review of the nonlinear updating methods is exposed by (Bussetta et al., 2017). Generally, the updating methods are based on the minimisation of an objective function or a metric that represents the difference between the numerical results and the experimental data. This objective function can use data in the frequency domain or in the time domain. The Volterra series can be used to defined this objective function (e.g. (Wu and Kareem, 2014); (Guo et al., 2013)). In the paper of (Shiki et al., 2012) the authors made a numerical comparison of a nonlinear model updating technique based on Volterra series and proper orthogonal decomposition. The last technique showed severe limitations, especially when considering higher levels of excitation of a model with lumped polynomial stiffness while the Volterra series showed a better performance with higher levels of nonlinearity. Generally speaking, the updating methods are applied by assuming linear models of the structure with lumped non-linearities, often assuming some kind of nonlinear springs as being the source of the nonlinear behaviour. Although this kind of model is able to reproduce the nonlinear phenomena (e.g. (Kerschen et al., 2003)), it is actually a rough simplification of the structure.

1.1 Originality of this Work

The originality of this work stems from the updating of nonlinear finite element model of structures using the Volterra series, rather than using simplified linear elements with lumped nonlinearity. The dynamical behaviour of structure is modelled using nonlinear finite element model and experimental data are used to identify numerical parameters. These parameters are split in linear and nonlinear parameters. The value of the linear parameters is identified using linear error indicator which is computed with the linear response of Volterra series. The identification of nonlinear parameters is based on the full Volterra series.

1.2 Outline

This paper is organised in five sections. First, the experimental setup of a preloaded buckled beam is described (section 2). Then, the nonlinear finite element model is presented in the part 3. The section 4 deals with the model updating based on the Volterra kernels expanded in a Kautz orthogonal basis. Finally, some remarks are summarised (section 5).

2 EXPERIMENTAL SETUP

This experimental setup was already used by (Scussel and da Silva, 2016) to experimentally identify the value of the first three Volterra kernels and by (Shiki et al., 2014a) for damage detection. It has been the subject of experimental study to identify cubic stiffness non-linearity of single degree-of-freedom Duffing system (Tang et al., 2015). The aim of this experimental setup is to study the vibration of a preloaded beam. The test structure is composed of a vertical aluminium beam with dimensions of 460x18x2 mm (see Figure 1). The bottom of the beam is clamped and a preload is applied on the top of the beam. The experimental setup has a screw in the top of the beam to apply an axial preload force. The preload was adjusted by trial and error until just before the buckling load of the beam (Tang et al., 2015). This was done to obtain a stronger nonlinear effect since the linear stiffness is reduced by applying the preload (Kovacic and Brennan, 2011). The excitation of the beam is generated by an electrodynamic shaker attached at 65 mm from the bottom of the beam. A load cell is placed between the shaker and the beam to record the exciting force with an acquisition board system. The oscillation of the beam is measured using the velocity captured by a laser vibrometer in the middle of the beam. The sampling frequency is equal to 1024 Hz and 4096 samples are recorded at each experimental test. In previous papers, (Scussel and da Silva, 2016) and (Shiki et al., 2014a) have demonstrated the nonlinear regime of vibration of this test bed using stepped sine test, FRFs and time-frequency plots. More details can be obtained in these references.

Figure 1
Experimental setup.

3 NONLINEAR FINITE ELEMENT MODEL

The model of the preloaded buckled beam uses the classical updated Lagrangian formulation of the 3D nonlinear finite element method. In each time step the displacement field is computed according to the momentum equation (see (Wriggers, 2008); (Belytschko et al., 2014)):

$\mathrm{\rho }\frac{{\partial }^2\mathbf{u}}{\partial {\mathrm{t}}^2}=\mathit{}\nabla .\mathbf{\sigma }+\mathbf{}\mathrm{\rho }\mathbf{b}$ (1)

Where u is the displacement, ρ is the mass density, σ is the Cauchy stress tensor and b is the specific body forces. The material is modelled as visco-elastic using the Kelvin-Voigt law, then the value of the stress is given by:

$\mathbf{\sigma }=\mathit{H}:\mathbf{\epsilon }+\mathrm{}\mathrm{\eta }\mathrm{}\widehat{\mathit{D}}$ (2)

where H is the fourth-order tensor of the Hooke's law, ε the strain tensor, $\widehat{D}$ the deviatoric parts of the strain-rate tensor and η the coefficient of viscosity. Note that the evolution of the Cauchy stress tensor, σ , is computed thanks to the Jaumann's objective time derivative. The mesh used in this model is composed of linear hexahedral elements. The displacement fields is computed at each node of the elements. The deviatoric part of the Cauchy stress is evaluated using a 2x2x2 Gauss quadrature scheme. To overcome the locking phenomenon, the pressure is considered constant over the element and computed only at a central quadrature point. This is a so called selective reduced integration. After integration of the weak form of the equation (1) over the mesh, the system of equations becomes:

$\mathbf{M}\mathrm{}\ddot{\mathit{u}}+\mathrm{}\underset{{\mathit{F}}_{int}(\mathit{u},\dot{\mathit{u}})}{\underbrace{\mathit{K}\left(\mathit{u},\dot{\mathit{u}}\right)\mathit{}\mathit{u}}}={\mathbf{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}$ (3)

where M is the mass matrix and _^Fint and F _ext are the internal and external force vectors. K is the stiffness matrix. The dynamical behaviour is computed thanks to the generalised-α method proposed by (Chung and Hulbert, 1993). With this time integration scheme, the unknowns are computed at the time t + ∆t due to the value at the time t. The equilibrium equation is rewritten pondering inertia forces by the parameter _^αM , and internal and external forces by the parameter _^αF:

$\begin{array}{cc}\left(1-{\mathrm{\alpha }}_{\mathrm{M}}\right)\mathbf{M}\mathrm{}{\ddot{\mathit{u}}}^{t+\mathrm{\Delta }t}& +\mathrm{}{\mathrm{\alpha }}_{\mathrm{M}}\mathbf{M}\mathrm{}{\ddot{\mathit{u}}}^t\\ & +\mathit{}\left(1-{\alpha }_F\right){\mathit{F}}_{int}^{t+\mathrm{\Delta }t}+{\alpha }_F{\mathit{F}}_{int}^t\\ & =\mathit{}(1{-{\alpha }_F)\mathit{F}}_{ext}^{t+\mathrm{\Delta }t}+{\alpha }_F{\mathit{F}}_{ext}^t\end{array}$ (4)

Moreover, the relations between the displacements, the velocities and the accelerations are:

${\mathit{u}}^{t+\mathrm{\Delta }t}={\mathit{u}}^t+\Delta t\mathbf{}\mathbf{}{\dot{\mathit{u}}}^t\mathbf{}+\mathrm{\Delta }{\mathrm{t}}^2\left(\left(\frac{1}{2}-\mathrm{\beta }\right)\mathrm{}{\ddot{\mathit{u}}}^{\mathit{t}}+\mathrm{\beta }\mathrm{}{\ddot{\mathit{u}}}^{\mathit{t}+\mathbf{\Delta }\mathit{t}}\right)$ (5)

${\dot{\mathit{u}}}^{t+\mathrm{\Delta }t}={\dot{\mathit{u}}}^t+\Delta t\mathbf{}\left(\right(1-\mathrm{\gamma })\mathbf{}{\ddot{\mathit{u}}}^{\mathit{t}}\mathbf{}+\mathrm{\gamma }\mathrm{}{\ddot{\mathit{u}}}^{\mathit{t}+\mathbf{\Delta }\mathit{t}}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}$ (6)

Second order accuracy, unconditional stability and optimal numeric dissipation of high frequencies is obtained for:

${\mathrm{\alpha }}_{\mathrm{M}}=\frac{2\mathrm{}{\mathrm{\rho }}_{\mathrm{\infty }}-\mathrm{}1}{{\mathrm{\rho }}_{\mathrm{\infty }}+\mathrm{}1}$ (7)

${\mathrm{\alpha }}_{\mathrm{F}}=\frac{{\mathrm{\rho }}_{\mathrm{\infty }}}{{\mathrm{\rho }}_{\mathrm{\infty }}+\mathrm{}1}$ (8)

$\mathrm{\beta }=\frac{1}{{{(\mathrm{\rho }}_{\mathrm{\infty }}+\mathrm{}1)}^2}$ (9)

$\gamma =\frac{1}{2}\frac{3-{\mathrm{\rho }}_{\mathrm{\infty }}}{{\mathrm{\rho }}_{\mathrm{\infty }}+\mathrm{}1}$ (10)

where ρ_∞ is approximately the percentage of numerical damping for the highest frequency of the structure.

4 MODEL UPDATING

Firstly, the test structure is submitted to a low level of amplitude force generated by the shaker in order to observe the linear response of the beam. Next, a larger amplitude force is used to study the nonlinear behaviour of the structure. The input used is a chirp signal sweeping up the frequency range from 20 to 50 Hz during the first half of the test (2 seconds). This frequency range comprises the first vibration mode which was considerably affected by the nonlinearity of the structure. During the second half of the test the shaker is turned off and the beam vibrates freely until it reaches the rest position. The low level of amplitude force and the high level one are respectively used to identify the first Volterra kernel and the higher-order Volterra kernels (i.e. the second and the third Volterra kernels). Both input signals are used to update the 3D nonlinear finite element model. Finally, a middle level of amplitude force generated by the chirp input signal as well as a low, a middle and a high level of a 35 Hz sinusoidal input signal are used to show that essentially the output signal from the updated nonlinear finite element model and the experimental one are similar.

Figure 2 presents the scheme of the numerical model. This numerical model use the homemade code Metafor of the University of Liege (see the official website of Metafor 2017). Like the experimental results, the solution of the numerical model is computed to the time of 4 seconds. In absence of any indication, the value of the numerical parameters are the following ones. The time step is 1/1024 seconds, the percentage of numerical damping for the highest frequency, ρ_∞, is 0.01. The mesh is composed of 9 600 hexahedral elements. The length, the width and the depth are respectively meshed with 300, 8 and 4 elements. The aluminium beam is modelled as visco-elastic material (Young's modulus: 58 GPa, Poisson's ratio: 0.33, density: 2 700 kg m^-3, coefficient of viscosity: η = 3 MPa s). Unfortunately, an important unknown of this numerical model is the boundary condition on the top of the beam. Since in the used experimental setup it is not possible to identify the value of the preloaded force, it is assumed that the boundary conditions on the top of the beam are modelled with an constant force, F, a linear spring - with a stiffness noted k_s - and a viscous damper - with a damping coefficient noted c (see Figure 2).

Figure 2
Scheme of the numerical model.

The constant force - F - is incorporated into the numerical model by a constant pressure apply on the boundary elements. The behaviour of the linear spring as well as the behaviour of the viscous damper are respectively modelled by linear spring and viscous damper elements connected to each node of this boundary. The value of the stiffness - k^e _s - and the damping coefficient - c^e - of these elements are divide by two if the connected node is on the edges and by four if it is on the corners. The stiffness and the damping coefficient of these elements are defined such as: ∑k^e _s ⁼ k_s and ∑c^e =c.

The updating process will be dealt with these three parameters - F, k_s and c - as well as the coefficient of viscosity - η. These parameters can be split into the linear and the nonlinear parameters. The linear parameter - η - can be identified with only the linear behaviour of the structure. The nonlinear parameters - F, k_s and c - have to be identified thanks to nonlinear vibration of the structure.

4.1 Volterra Model

The Volterra kernels are used to update the 3D nonlinear finite element model. The Volterra series are based on the fact that the output of a system (noted y) can be given by the sum of linear and nonlinear relations with the input (noted u). The output is described by:

$\mathrm{y}\left(\mathrm{k}\right)=\mathrm{}{\mathrm{y}}_1\left(\mathrm{k}\right)+\mathrm{}{\mathrm{y}}_2\left(\mathrm{k}\right)+\mathrm{}{\mathrm{y}}_3\left(\mathrm{k}\right)+\cdots$ (11)

where y₁ =(k) is the linear term and y_η =(k) with η>1 are the nonlinear terms. Each term of the equation (11) can be written using the discrete-time Volterra series expansion:

${\mathrm{y}}_{\mathrm{\eta }}\left(\mathrm{k}\right)=\sum _{{\mathrm{n}}_1=0}^{{\mathrm{N}}_1}\cdots \sum _{{\mathrm{n}}_{\mathrm{\eta }}=0}^{{\mathrm{N}}_{\mathrm{\eta }}}{\mathit{H}}_{\mathrm{\eta }}\left({\mathrm{n}}_1,\mathrm{}\cdots ,\mathrm{}{\mathrm{n}}_{\mathrm{\eta }}\right)\prod _{\mathrm{i}=1}^{\mathrm{\eta }}\mathrm{u}(\mathrm{k}-{\mathrm{n}}_{\mathrm{i}})$ (12)

where H _η is the η-order Volterra kernel and N_j with 1 ≤ j ≤ η are the size of this one. Problems of overparametrisation can appear due to the high numbers of unknowns (N_1×...N_j×...Nη) needed to represent the discrete Volterra kernels, H η. To overcome the ill-posed and convergence problems, the Volterra kernel can be expanded in orthonormal basis. It means that each Volterra kernel can be approximated by it projection over an orthonormal basis. The η-order Volterra kernel is written as a combination of orthonormal functions ψ _η,ij(n_j) (see (Marmarelis, 2004)):

${\mathit{H}}_{\mathrm{\eta }}\left({\mathrm{n}}_1,\mathrm{}\cdots ,\mathrm{}{\mathrm{n}}_{\mathrm{\eta }}\right)\approx \sum _{{\mathrm{i}}_1=0}^{{\mathrm{J}}_1}\cdots \sum _{{\mathrm{i}}_{\mathrm{\eta }}=0}^{{\mathrm{J}}_{\mathrm{\eta }}}{\mathit{B}}_{\mathrm{\eta }}\left({\mathrm{i}}_1,\mathrm{}\cdots ,\mathrm{}{\mathrm{i}}_{\mathrm{\eta }}\right)\prod _{\mathrm{j}=1}^{\mathrm{\eta }}{\psi }_{\mathrm{\eta },\mathrm{}{\mathrm{i}}_{\mathrm{j}}}\left({\mathrm{n}}_{\mathrm{j}}\right)$ (13)

where B _η is the projection of the η-order Volterra kernel over this orthonormal basis. Obviously, the number of terms of the Volterra kernels, J_j, is smaller than the number of terms with the direct computation (N_j in equation (12)).

In the case of nonlinear oscillating system, the orthonormal Kautz functions can be used (Kautz, 1954). These functions can be represented in the frequency domain by:

$\begin{array}{ccc}{\mathrm{\Psi }}_{\mathrm{\eta },\mathrm{}2\mathrm{g}}\left(\mathrm{z}\right)& =& \frac{\sqrt{(1-{\mathrm{d}}_{\mathrm{\eta }}^2)(1-{\mathrm{b}}_{\mathrm{\eta }}^2)\mathrm{z}}}{{\mathrm{z}}^2+{\mathrm{b}}_{\mathrm{\eta }}({\mathrm{d}}_{\mathrm{\eta }}-1)\mathrm{z}-{\mathrm{d}}_{\mathrm{\eta }}}{\left[\frac{-{\mathrm{d}}_{\mathrm{\eta }}{\mathrm{z}}^2+{\mathrm{b}}_{\mathrm{\eta }}\left({\mathrm{d}}_{\mathrm{\eta }}-1\right)\mathrm{z}+1}{{\mathrm{z}}^2+{\mathrm{b}}_{\mathrm{\eta }}({\mathrm{d}}_{\mathrm{\eta }}-1)\mathrm{z}-{\mathrm{d}}_{\mathrm{\eta }}}\right]}^{g-1}\\ {\mathrm{\Psi }}_{\mathrm{\eta },\mathrm{}2\mathrm{g}-1}\left(\mathrm{z}\right)& =& \frac{\mathrm{z}-{\mathrm{b}}_{\mathrm{\eta }}}{\sqrt{1-{\mathrm{b}}_{\mathrm{\eta }}^2}}{\mathrm{\Psi }}_{\mathrm{\eta },\mathrm{}2\mathrm{g}}\left(\mathrm{z}\right)\end{array}$ (14)

where z is the complex variable, g∈ ⟦1, ..., J_η/2⟧ and b _η and d _η are parameters defined thanks to the pair of conjugate poles of the Kautz function (Z _η and $\stackrel{-}{Z_n}$). These poles are linked to the sampling frequency of the time-domain signal (F_s) and the continuous pole, S _η (i.e. the natural frequency ω _η and damping ratio ζ _η).

$\begin{array}{ccc}{\mathrm{b}}_{\mathrm{\eta }}& =& \frac{{\mathrm{Z}}_{\mathrm{\eta }}+{\stackrel{-}{\mathrm{Z}}}_{\mathrm{\eta }}}{1+{\mathrm{Z}}_{\mathrm{\eta }}{\stackrel{-}{\mathrm{Z}}}_{\mathrm{\eta }}}\\ {\mathrm{d}}_{\mathrm{\eta }}& =& -{\mathrm{Z}}_{\mathrm{\eta }}{\stackrel{-}{\mathrm{Z}}}_{\mathrm{\eta }}\\ \begin{array}{c}{\mathrm{Z}}_{\mathrm{\eta }}\\ {\mathrm{S}}_{\mathrm{\eta }}\end{array}& \begin{array}{c}=\\ =\end{array}& \begin{array}{c}{e}^{{\mathrm{S}}_{\mathrm{\eta }}/F_s}\\ -{\mathrm{\zeta }}_{\mathrm{\eta }}{\mathrm{\omega }}_{\mathrm{\eta }}+\mathrm{}\mathrm{j}{\mathrm{\omega }}_{\mathrm{\eta }}\sqrt{1-{\mathrm{\zeta }}_{\mathrm{\eta }}^2}\end{array}\end{array}$ (15)

Generally, the values describing the pair of conjugate poles of the Kautz functions (i.e. ω _η and ζ _η) are close to the values representing the linearised system (i.e. the natural frequency and damping ratio). The value in the time domain of ψ _η,ij is obtained by an inverse Z-transform of ψ _η,ij.

Finally, the value of y_η can be written as:

${\mathrm{y}}_{\mathrm{\eta }}\left(\mathrm{k}\right)=\sum _{{\mathrm{i}}_1=0}^{{\mathrm{J}}_1}\cdots \sum _{{\mathrm{i}}_{\mathrm{\eta }}={\mathrm{i}}_{\mathrm{\eta }-1}}^{{\mathrm{J}}_{\mathrm{\eta }}}{\mathit{B}}_{\mathrm{\eta }}\left({\mathrm{i}}_1,\mathrm{}\cdots ,\mathrm{}{\mathrm{i}}_{\mathrm{\eta }}\right)\prod _{\mathrm{j}=1}^{\mathrm{\eta }}{\mathrm{l}}_{\mathrm{\eta },{\mathrm{i}}_{\mathrm{j}}}\left(\mathrm{k}\right)$ (16)

where l _η,ij is the input signal filtered by the orthonormal Kautz function ψ _η,ij (i.e. the convolution between the input signal and the Kautz function):

$l_{\mathrm{\eta },{\mathrm{i}}_{\mathrm{j}}}\left(\mathrm{k}\right)=\sum _{{\mathrm{n}}_{\mathrm{j}}=0}^{{\mathrm{N}}_{\mathrm{j}}}{\psi }_{\mathrm{\eta },{\mathrm{i}}_{\mathrm{j}}}\left({\mathrm{n}}_{\mathrm{j}}\right)\mathrm{u}(\mathrm{k}-{\mathrm{n}}_{\mathrm{j}})$ (17)

Since the Volterra model is linear with respect to the parameters, the relationship between the input filtered by the orthonormal Kautz functions, the output and the Volterra kernels can be rewritten as a simple matrix equation:

$\mathbf{y}\mathbf{}=\mathbf{\Gamma }\mathbf{}\mathbf{\Phi }$ (18)

where Γ is the matrix of the input signal filtered by the orthonormal Kautz functions and Φ is the vector with the Volterra kernels considered in the truncation of the model. The Volterra kernels can be evaluated by the classical least squares method:

$\mathbf{\Phi }={\left({\mathbf{\Gamma }}^{\mathrm{T}}\mathbf{\Gamma }\right)}^{-1}{\mathbf{\Gamma }}^{\mathrm{T}}\mathbf{y}$ (19)

For more information about the evaluation of Volterra kernels see (Shiki, 2016); (Shiki et al., 2014b).

4.2 Identification of the First Volterra Kernel - Updating of Linear System

The experimental data acquired during the test with the low level of chirp signal is used to identify the value of the parameters - the coefficient of the viscous damper c, the preloaded force F, the spring stiffness k_s and the material damping coefficient η. This level of amplitude force corresponds to a voltage of 0.01 V applied in the signal generator. Figure 3 shows the value of the exciting force recorded during the test and applied in the numerical model. Due to the low level of amplitude force, the response of the beam is linear. Then, the metric based on the first Volterra kernel is used (i.e. using the output corresponding to the first Volterra kernel, see equation (16)):

${\mathrm{e}}_{\mathrm{l}\mathrm{i}\mathrm{n}}(\mathrm{F},\mathrm{}{\mathrm{k}}_{\mathrm{s}},\mathrm{}\mathrm{c},\mathrm{}\mathrm{\eta })=\frac{‖{\mathrm{y}}_1^{\mathrm{r}\mathrm{e}\mathrm{f}}-{\mathrm{y}}_1^{\mathrm{m}\mathrm{o}\mathrm{d}}(\mathrm{F},\mathrm{}{\mathrm{k}}_{\mathrm{s}},\mathrm{}\mathrm{c},\mathrm{}\mathrm{\eta })‖}{‖{\mathrm{y}}_1^{\mathrm{r}\mathrm{e}\mathrm{f}}‖}$ (20)

where ^ref and ^mod refer respectively to the experimental data and the numerical model.

Figure 3
Value of the exciting force recorded during the test and applied in the numerical model (low level chirp input).

Figure 4 presents the evolution of the linear error, e_lin, versus the value of the preloaded force and the value of spring stiffness for different values of the coefficient of the viscous damper. This figure allows us to observe that the value of the preloaded force and the value of spring stiffness are linked. In addition, the value of the coefficient of the viscous damper, c, has only a weak influence over the value of the linear error. The minimal value of e_lin versus the value of the coefficient of the viscous damper and the corresponding value of the preloaded force and the spring stiffness are presented in Table 1. The linear error is minimised by: c = 10 N s mm^-1, F = 74 N and k_s = 0 N/mm.

Figure 4
Value of the linear error, e_lin, (i.e. based on the first order Volterra kernel) versus the value of the preloaded force and the value of spring stiffness, for the low level chirp input.

Table 1
The minimal value of the linear error, e_lin and the corresponding value of the preloaded force, F and the spring stiffness, k_s versus the value of the coefficient of the viscous damper, c, for the low level chirp input.

c (in N s/mm)	elin (in %)	F (in N)	ks (in N/mm)
0	2.27	80.25	390
5	1.76	82.75	540
10	1.59	74.00	0
15	1.61	77.00	180
20	1.64	79.50	330
25	1.67	81.50	450

Figure 5 shows the evolution of the velocity measured by the laser vibrometer and the velocity computed at the same position by the updated numerical model (c = 10 N s mm^-1, F = 74 N and k_s = 0 N mm^-1). The output of the updated model is near to the output of the experimental setup. Moreover, the output signal (Figure 5) is quasi identical to the linear response of the Volterra model (Figure 6) as far as the signal of the experimental setup and the numerical results are considered. The frequency response function (FRF) and the magnitude-squared coherence of both experimental setup and numerical model are presented in Figure 7. In this figure, the output of the updated model and the output of the experimental setup are very similar. Consequentially, it is proved that with this input level, the system behaviour is essentially linear.

Figure 5
Response to the low level chirp excitation of the experimental setup and computed with the numerical model (F = 74 N, k_s = 0 N mm^-1 and c = 10 N s mm^-1).

Figure 6
Linear response of the Volterra model of system (y₁) under the low level chirp excitation (F = 74 N, k_s = 0 N mm^-1 and c = 10 N s mm^-1).

Figure 7
Frequency response to the low level chirp excitation of the experimental setup and computed with the numerical model (F = 74 N, k_s = 0 N mm^-1 and c = 10 N s mm^-1).

Finally, we considered that the value of the material damping coefficient is independent from the value of the others parameters. In this way, the value of the material damping coefficient is updated independently. Figure 8 presents the value of the linear error versus the value of the material damping coefficient. The linear metric is minimised for the value of the material damping coefficient of 3 MPa s.

Figure 8
Value of the error based on the first Volterra kernel, e_lin, versus the value of the material damping coefficient, for the low level chirp input.

4.3 Identification of the Second and the Third Volterra Kernels - Updating of Nonlinear System

The system is excited by higher level chirp signal to put in light its nonlinear behaviour. This signal corresponds to 0.1 V applied in the signal generator. Figure 9 presents the evolution of the value of the exciting force.

Figure 9
Value of the exciting force recorded during the test and applied in the numerical model (high level chirp input).

Another metric based on the second and the third Volterra kernels is defined to update the nonlinear model (i.e. using the output corresponding to the second and the third Volterra kernels):

${\mathrm{e}}_{\mathrm{n}\mathrm{l}\mathrm{i}\mathrm{n}}(\mathrm{F},\mathrm{}{\mathrm{k}}_{\mathrm{s}},\mathrm{}\mathrm{c})=\frac{‖{\mathrm{y}}_{2+3}^{\mathrm{r}\mathrm{e}\mathrm{f}}-{\mathrm{y}}_{2+3}^{\mathrm{m}\mathrm{o}\mathrm{d}}(\mathrm{F},\mathrm{}{\mathrm{k}}_{\mathrm{s}},\mathrm{}\mathrm{c})‖}{‖{\mathrm{y}}_{2+3}^{\mathrm{r}\mathrm{e}\mathrm{f}}‖}$ (21)

with y*₂₊₃ = y*₂ +y*₃. The minimum of this nonlinear error is different from the minimum of the linear error, ^_elin . The difference between both error functions can be explained by the fact that the value of these parameters - F, k_s and c - have an influence on the linearised response as well as the nonlinear one. A new objective function is built to take into consideration the nonlinear behaviour of the structure as well as linear part:

${\mathrm{e}}_{\mathrm{l}\mathrm{i}\mathrm{n}+\mathrm{n}\mathrm{l}\mathrm{i}\mathrm{n}}={\mathrm{e}}_{\mathrm{l}\mathrm{i}\mathrm{n}}+{\mathrm{e}}_{\mathrm{n}\mathrm{l}\mathrm{i}\mathrm{n}}$ (22)

The minimum value of e_{lin +nlin} as well as the corresponding value of F and k_s are presented in Table 2 versus c. The minimum value is obtained for: c = 15 N s mm^-1, F = 77 N and k_s = 180 N/mm.

Table 2
The minimal value of e_{lin +nli n} and the corresponding value of the preloaded force, F and the spring stiffness, k_s versus the value of the coefficient of the viscous damper, c, for the high level chirp input.

c (in N s/mm)	elin+nlin (in %)	F (in N)	ks (in N/mm)
0	55.92	80.75	420
5	18.44	82.75	540
10	16.71	83.00	540
15	16.67	77.00	180
20	16.77	74.25	0
25	19.16	74.25	0

Figure 10 presents the value of the velocity recorded by the laser vibrometer and the velocity computed by the updated numerical model. The output of the updated model is very similar to the output of the experimental setup. In addition, with this model updated thanks to e_{lin +nlin}, the linear error, e_lin computed under the excitation of the low level chirp signal is close to the minimal value (see Table 1). The FRF and the magnitude-squared coherence (Figure 11) let us see that the numerical solution and the experimental data are close in the frequency domain. Moreover, the structure exhibit clearly a nonlinear behaviour. It can be verified that the components of both Volterra models (base on experimental data and using the updated numerical model) are also close (see Figure 12). The nonlinear part of the Volterra model (y₂ and y₃) is higher than the linear one (y₁). This is a clear demonstration that under the high level chirp excitation the system exhibits a nonlinear behaviour.

Figure 10
Response to the high level chirp excitation of the experimental setup and computed with the numerical model (F = 77 N, k_s = 180 N mm^-1 and c = 15 N s mm^-1).

Figure 11
Frequency response to the high level chirp excitation of the experimental setup and computed with the numerical model (F = 77 N, k_s = 180 N mm^-1 and c = 15 N s mm^-1). The FRF distortions of the plot indicate a nonlinear behaviour of the system.

Figure 12
Output of the Volterra Model computed with the experimental data and to the updated numerical model (F = 77 N, k_s = 180 N mm^-1 and c = 15 N s mm^-1), under the high level chirp excitation.

To put in a nutshell, the updating of the nonlinear finite element model is effected in two steps. First, the all parameters is evaluated thanks to the minimisation of the linear error (objective function using the first Volterra kernel). Then, only the nonlinear parameters are corrected by the minimisation of an objective function considering the three first Volterra kernels. The Table 3 resume the value of the parameters identified versus the updating process and the objective function.

Table 3
Value of the parameters identified versus the updating process and the objective function.

updating process	linear updating	nonlinear updating
objective function	elin	elin+nlin
parameters	η (in MPa s)	c (in N s/mm)	F (in N)	ks (in N/mm)
value	3	15	77.00	180

4.4 Validation of the Nonlinear Finite Element Model

To validate the updated nonlinear finite element model, a medium level chirp input is used. This signal is generated by 0.05 V applied in the signal generator. Figure 13 (a) presents the value of the exciting force. Figure 13 (b) allows us to verify that the velocity computed by the updated nonlinear model is near to the one recorded by the laser vibrometer.

Figure 13
Exciting force generated by the medium level chirp input and normal velocity in the middle of the beam recorded with the experimental setup and computed thanks to the numerical model.

Another kind of signal is used to validate the updated nonlinear model. Three levels of 35 Hz sinusoidal signal are used - low, medium and high level. These input signals correspond respectively to 0.01, 0.05 and 0.1 V applied in the signal generator. The power spectral density (PSD) of both experimental and numerical data are presented in Figure 14. These results show that the nonlinear finite element model is able to represent all harmonics in this range of frequencies (between 0 and 140 Hz). The comparisons of PSDs of experimental data and PSDs of numerical results validate the updated nonlinear model. In the same way, this results prove that the model updating method based on the Volterra kernels can be used for nonlinear finite element models with large number of degree of freedom.

Figure 14
Power spectral density (PSD) of the experimental setup and computed thanks to the numerical model.

5 FINAL REMARKS

In this paper, a model updating method using the Volterra series is proposed for nonlinear finite element model with large number of degree of freedom. One of the main advantages of the Volterra series is that the relationship between the input ant the output can be split into a linear response and a nonlinear one. This model updating method make use of the Volterra series to simplify the procedure. Consequently, this updating method is split in two parts. First, the linearised system (i.e. the nonlinear system under low amplitude of vibrations) is used to update the numerical model using the first order Volterra kernel. Then, the nonlinear parameters are updated using a high amplitude of vibrations and the higher-order Volterra kernels. This process is used to update nonlinear finite element model of preloaded buckled beam thanks to data of experimental setup. Chirp input signals with low and high level of vibration are used to update the nonlinear finite element model. The updated nonlinear model is validated with a medium level chirp input and with three levels of sinusoidal signal. These examples show that the updated nonlinear finite element model essentially reproduce the output of the experimental setup in the frequency range next to the first vibration mode which was strongly affected by the nonlinearity. In the way of reducing the computation time, the future study will be focused on the computation of the first Volterra kernel using the linearised numerical finite element model.

Acknowledgements

The authors gratefully acknowledge the financial support for this research from “Conselho Nacional de Desenvolvimento Científico e Tecnológico” (CNPq); grant number 313661/2014-6 and grant number 307520/2016-1 and FAPESP; grant number 12/09135-3. The second author is thankful to FAPESP for his scholarship grant number 13/25148-0.

References

• Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K., (2014). Nonlinear Finite Elements for Continua and Structures. New York: Wiley.
• Bussetta, P., Shiki, S.B., da Silva, S., (2017). Nonlinear updating method: a review. Journal of the Brazilian Society of Mechanical Sciences and Engineering.
• Chung, J., Hulbert, G., (1993). A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method. Journal of Applied Mechanics; 60(2):371-375. doi: 10.1115/1.2900803.
  » https://doi.org/10.1115/1.2900803
• Guo, Y., Guo, L.Z., Billings, S.A., Coca, D., Lang, Z.Q., (2013). Volterra series approximation of a class of nonlinear dynamical systems using the adomian decomposition method. Nonlinear Dynamics; 74(1):359-371. doi: 10.1007/s11071-013-0975-8.
  » https://doi.org/10.1007/s11071-013-0975-8
• Kautz, W.H., (1954). Transient synthesis in the time domain. Transactions of the IRE Professional Group on Circuit Theory; 3:29-39. doi: 10.1109/TCT.1954.1083588.
  » https://doi.org/10.1109/TCT.1954.1083588
• Kerschen, G., Lenaerts, V., Golinval, J.C., (2003). Identification of a continuous structure with a geometrical non-linearity. Part I: Conditioned reverse path method. Journal of Sound and Vibration; 262(4):889-906. doi: 10.1016/S0022-460X(02)01151-3.
  » https://doi.org/10.1016/S0022-460X(02)01151-3
• Kovacic, I., Brennan, M.J., editors, (2011). The Duffing equation: Nonlinear oscillators and their behaviour. John Wiley & Sons. doi: 10.1002/9780470977859.
  » https://doi.org/10.1002/9780470977859
• Marmarelis, V.Z., (2004). Nonlinear dynamic modeling of physiological systems. John Wiley & Sons.
• Noël, J., Kerschen, G., (2017). Nonlinear system identification in structural dynamics 10 more years of progress. Mechanical Systems and Signal Processing; 83:2-35. doi: 10.1016/j.ymssp.2016.07.020.
  » https://doi.org/10.1016/j.ymssp.2016.07.020
• Official website of Metafor, (2017). URL: http://metafor.ltas.ulg.ac.be/dokuwiki
  » http://metafor.ltas.ulg.ac.be/dokuwiki
• Scussel, O., da Silva, S., (2016). Output-only identification of nonlinear systems via volterra series. Journal of Vibration and Acoustics; 138(4):041012. doi: 10.1115/1.4033458.
  » https://doi.org/10.1115/1.4033458
• Shiki, S.B., (2016). Application of Volterra series in nonlinear mechanical system identification and in structural health monitoring problems. Ph.D. thesis; Universidade Estadual Paulista - UNESP. URL: http://hdl.handle.net/11449/137761
  » http://hdl.handle.net/11449/137761
• Shiki, S. B., Lopes Jr, V. and da Silva, S., (2012). Model updating of the non-linear vibrating structures through Volterra series and proper orthogonal decomposition. Proceedings Of International Conference On Noise And Vibration Engineering (isma2012) / International Conference On Uncertainty In Structural Dynamics (usd2012); 2199-2212. URL: http://hdl.handle.net/11449/117716
  » http://hdl.handle.net/11449/117716
• Shiki, S.B., Hansen, C., da Silva, S., (2014a). Nonlinear features identified by Volterra series for damage detection in a buckled beam. In: International Conference on Structural Nonlinear Dynamics and Diagnosis.
• Shiki, S.B., Lopes, V., da Silva, S., (2014b). Identification of nonlinear structures using discrete-time Volterra series. Journal of the Brazilian Society of Mechanical Sciences and Engineering; 36(3):523-532. doi: 10.1007/s40430-013-0088-9.
  » https://doi.org/10.1007/s40430-013-0088-9
• Tang, B., Brennan, M.J., Lopes, J.V., da Silva, S., Ramlan, R., (2015). Using nonlinear jumps to estimate cubic stiffness nonlinearity: An experimental study. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science; (x):xxx-xxx. doi: 10.1177/0954406215606746.
  » https://doi.org/10.1177/0954406215606746
• Worden, K., Tomlinson, G.R., (2001). Nonlinearity in Structural Dynamics - Detection, Identification and Modelling. University of Sheffield: Institute of Physics Publishing.
• Wriggers, P., (2008). Nonlinear Finite Element Methods. Springer Berlin Heidelberg. doi: 10.1007/978-3-540-71001-1.
  » https://doi.org/10.1007/978-3-540-71001-1
• Wu, T., Kareem, A., (2014). Simulation of nonlinear bridge aerodynamics: A sparse third-order Volterra model. Journal of Sound and Vibration; 333(1):178-188. doi: 10.1016/j.jsv.2013.09.003.
  » https://doi.org/10.1016/j.jsv.2013.09.003

Publication Dates

• Publication in this collection
  Aug 2017

History

• Received
  21 Mar 2017
• Reviewed
  21 Apr 2017
• Accepted
  25 Apr 2017