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

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.

Latin American Journal of Solids and Structures 14 (2017) 1183-1199 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 nonlinearities, 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.

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.

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).

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-offreedom 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.

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)): Latin American Journal of Solids and Structures 14 (2017) 1183-1199 Where is the displacement, ρ is the mass density, is the Cauchy stress tensor and 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: where is the fourth-order tensor of the Hooke's law, the strain tensor, 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: where is the mass matrix and and are the internal and external force vectors.
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 Δ due to the value at the time .The equilibrium equation is rewritten pondering inertia forces by the parameter , and internal and external forces by the parameter : Moreover, the relations between the displacements, the velocities and the accelerations are: Second order accuracy, unconditional stability and optimal numeric dissipation of high frequencies is obtained for: Latin American Journal of Solids and Structures 14 (2017) 1183-1199 where ρ is approximately the percentage of numerical damping for the highest frequency of the structure.

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 signal1 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 ks -and a viscous damper -with a damping coefficient noted c (see figure 2).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 -and the damping coefficient -c -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 k and ∑c c.
The updating process will be dealt with these three parameters -F, ks 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, ks and c -have to be identified thanks to nonlinear vibration of the structure.

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: 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: where is the -order Volterra kernel and N with 1 j η are the size of this one.Problems of overparametrisation can appear due to the high numbers of unknowns (N ⋯ N ⋯ N ) needed to represent the discrete Volterra kernels, .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 , n (see Marmarelis (2004)): where is the projection of the -order Volterra kernel over this orthonormal basis.Obviously, the number of terms of the Volterra kernels, J , is smaller than the number of terms with the direct computation (N in equation ( 12)).
Latin American Journal of Solids and Structures 14 (2017) 1183-1199 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: where z is the complex variable, g ∈ 1, ⋯ , J /2 and and are parameters defined thanks to the pair of conjugate poles of the Kautz function ( and ̅ ).These poles are linked to the sampling frequency of the time-domain signal (Fs) and the continuous pole, (i.e. the natural frequency and damping ratio ).
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 , is obtained by an inverse Z-transform of Ψ , .
Finally, the value of y can be written as: where , is the input signal filtered by the orthonormal Kautz function , (i.e. the convolution between the input signal and the Kautz function): 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: (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: Latin American Journal of Solids and Structures 14 (2017) 1183-1199 (19) For more information about the evaluation of Volterra kernels see Shiki (2016); Shiki et al. (2014b).

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 ks 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)):    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 ks = 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.
Latin American Journal of Solids and Structures 14 (2017) 1183-1199  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.

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.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): with y * y * y * .The minimum of this nonlinear error is different from the minimum of the linear error, . The difference between both error functions can be explained by the fact that the value of these parameters -F, ks 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: The minimum value of e as well as the corresponding value of F and ks are presented in table 2 versus c.The minimum value is obtained for: c = 15 N s mm -1 , F = 77 N and ks = 180 N/mm., the linear error, e 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.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.

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.

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.oidal 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.

Figure 2 :
Figure 2: Scheme of the numerical model.
the experimental data and the numerical model.

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

Figure 4 Figure 4 :
Figure 4 presents the evolution of the linear error, e , 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 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 ks = 0 N/mm.

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

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

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

Figure 10
Figure10presents 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 , the linear error, e computed under the excitation of the low level chirp signal is close to the minimal value (see table1).The FRF and the magnitude-squared coherence (figure11) 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 figure12).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 :
Figure 10: Response to the high level chirp excitation of the experimental setup and computed with the numerical model (F = 77 N, ks = 180 N mm -1 and c = 15 N s mm -1 ).

Figure 11 :Figure 12 :
Figure 11: Frequency response to the high level chirp excitation of the experimental setup and computed with the numerical model (F = 77 N, ks = 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 13 :Figure 14 :
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.
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 sinus-Latin American Journal of Solids and Structures 14 (2017) 1183-1199

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

Table 2 :
The minimal value of e and the corresponding value of the preloaded force, F and the spring stiffness, ks versus the value of the coefficient of the viscous damper, c, for the high level chirp input.

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