Identiﬁcation of stochastic loads applied to a non-linear dynamical system using an uncertain computational model

This paper deals with the identiﬁcation of stochastic loads applied to a non-linear dynamical system for which a few experimental responses are available using an uncertain computational model. Uncertainties are induced by the use of a simpliﬁed computational model to predict the responses of the real system. A non-parametric probabilistic approach of both parameter uncertainties and model uncertainties are implemented in the simpliﬁed computational model in order to take into account uncertainites. The level of uncertainites is identiﬁed using the maximum likelihood method. The identiﬁed stochastic simpliﬁed computational model which is obtained is then used to perform the identiﬁcation of the stochastic loads applied to the real non-linear dynamical system. A numerical validation of the complete methodology is presented.


Introduction
This paper is devoted to the identification of a stochastic load applied to a non-linear dynamical system for which a few measurements of its responses are available and for which an uncertain simplified computational model is used. In the dynamical system, the uncertainties are taken into account in the context of the probability theory. Consequently the uncertain simplified computational model is in fact a stochastic simplified computational model for which the input is a stochastic process (stochastic load) and for which the linear operators of the computational model are random. This identification is then performed using the stochastic simplified computational model which allows the responses of the real system to be predicted and then which allows the stochastic loads to be identified in minimizing a certain distance between the experimental responses and the random responses calculated with the stochastic simplified computational model. In fact, the methodology presented is developed in the context of the non-linear dynamical analysis of tube bundles in Pressurized Water Reactor. The stochastic loads applied to the tubes which have to be identified are then induced by the turbulent flow. Since such a non-linear dynamical system is very complex, the computational model developed can not exactly represent the complexity of the system. Consequently, the identification is not performed using a computational model wich has the capability to accurately predict the experimental responses but is performed using a simplified computational model containing model errors. In order to perform a robust identification of the stochastic loads with respect to model uncertainties in the non-linear dynamical computational model, a probabilistic model of uncertainties allowing both parameter uncertainties and model uncertaintes to be taken into account is introduced. The responses of the computational model are then random and the randomness is due to the stochastic loads and is due to the stochasticity of the system. In a first step, the probability model of uncertainties in the computational model is identified using the maximum likelihood method. We then deduce a stochastic computational model which allows a robust identification of stochastic loads to be carried out with respect to uncertainties in the non-linear computational model. The second step is devoted to the stochastic inverse problem consisting in identifying the stochastic loads. From a theoretical and methodological point of view, we then present a complete probabilistic construction and the associated methodology to solve an inverse problem consisting of the identification of a Gaussian stationary stochastic process which is the input of a continuous nonlinear dynamical system with random operators and for which the stochastic output is measured. It should be noted that, if the parametric probabilistic approach is usual to take into account system parameter uncertainties, in the present paper, both the system parameter uncertainties and the model uncertainties are taken into account using a nonparametric probabilistic approach consisting in directly modeling the linear operators of the dynamical system by random operators using the random matrix theory. Section 2 deals with the construction of the mean computational model. In Section 3, the probabilistic model of the stochastic loads is introduced. Section 4 is devoted to the identification of the stochastic load. The last section presents a numerical validation of the methodology proposed.

Mean computational model
Let Ω be the domain of the dynamical system having a non-linear behaviour due to the presence of elastic stops located to several points of the part of the boundary of Ω. The domain Ω is decomposed in two bounded open subdomains of R 3 , the subdomain Ω A and the subdomain Ω B . The subdomain Ω A is constituted of a three dimensional linear viscoelastic medium with instantaneous memory and there are elastic stops located at κ points x 1 , .., x κ in the boundary Γ A of ∂Ω A . In addition, the subsystem occupying the subdomain Ω A is fixed on the part Γ A 0 of its boundary ∂Ω A . The outward unit normal of ∂Ω A is denoted n A . The subdomain Ω B is constituted of a three dimensional linear viscoelastic medium with instantaneous memory, fixed on the part Γ B 0 of its boundary ∂Ω B . The outward unit normal of ∂Ω B is denoted n B . Consequently, each uncoupled subsystem Ω A and Ω B does not have rigid body displacement. These two subsytems are coupled on the common coupling interface Γ C . One then has We are interested in constructing the stationary random responses of the non-linear stochastic dynamical system excited by stationary stochastic processes. Consequently, we will not introduce the initial conditions and we will assume that the time parameter t belongs to R.

Mean boundary value problems
2.1.1 Mean boundary value problems for the linear subsystem Ω B Let x = (x 1 , x 2 , x 3 ) be the cartesian coordinates and u B (x, t) be the displacement field for the linear subsystem Ω B at time t. The external prescribed volumetric and surface forces fields applied to Ω B and to its bound- is the linearized strain tensor. The fourth-order tensors a B (x) and b B (x) verify the usual properties of symmetry and positiveness ( [4]). Then, the displacement field u B (t) verifies, for all t ∈ R and 3 hal-00686141, version 1 -7 Apr 2012 for i = 1, 2, 3, the mean boundary value problem in which a dot means the partial time derivative and a double dot means the double partial time derivative, where f B coupl = (f B coupl,1 , ..., f B coupl,3 ) is the forces induced by the subsystem Ω A on Ω B via the coupling interface Γ C . One has used the classical convention for summations over repeated Latin indices. The parameter ρ B (x) is the mass density for the subsystem Ω B .

Mean boundary value problems for the non-linear subsys-
tem Ω A Let u A (x, t) be the displacement field for the non-linear subsystem Ω A at time t. The external prescribed volumetric and surface forces fields applied to Ω A and to its boundary Since Ω A is occupied by a linear viscoelastic material with instataneous memory, the stress tensor σ A (x) is written as is the linearized strain tensor. The fourth-order tensors a A (x) and b A (x) verify as above, the usual properties of symmetry and positiveness. Then, the displacement field u A (t) verifies, for all t ∈ R and for i = 1, 2, 3, the mean boundary value problem in which f A coupl = (f A coupl,1 , ..., f A coupl,3 ) is the forces induced by the subsystem Ω B on Ω A via the coupling interface Γ C . The forces −f NL,k (u(x k , t)) represent the actions exerted by the elastic stop located at point x k on the subsystem Ω A and δ 0 (x − x k ) is the surface Dirac measure such that, for all continuous function g defined on Γ A , one has Γ A δ 0 (x − x k )g(x)ds(x) = g(x k ). The parameter ρ A (x) is the mass density for the subsystem Ω A .

Interface conditions for the coupling of Ω
The coupling conditions on Γ C are written as f A coupl + f B coupl = 0 on Γ C .

Mean finite element model
The mean finite element method ( [16]) is applied to the variational formulation of the boundary value problems defined in subsection 2.1.

Mean finite element model for subsystem Ω B
The R n B vector U B (t) of the n B DOF of the subsystem Ω B is written as c vector function of the n B c coupling DOF on the interface. From Eq. 1, it can be deduced that the mean finite element model of subsystem Ω B is written as

Mean finite element model for subsystem Ω A
Similarly to Subsection 2.2.1, the R n A vector U A (t) of the n A DOF of the subsystem Ω A is written as it can be deduced that the mean finite element model of subsystem Ω A is written as are respectively the positive-definite symmetric real positive (n A ×n A ) mass, damping and stiffness matrices. The R n A vectors F A (t), F A coupl (t) and F NL (U A (t)) of the external forces, of the coupling forces and of the non-linear forces are written as

Interface conditions for the coupling of Ω A with Ω B
The finite element discretization of the interface conditions defined by Eq. 3 yields

Reduced mean matrix model
The continuous linear subsytem Ω B (linear dynamical subsystem) contains elastic modes in the frequency band of analysis. In addition, the computational model of the continuous linear subsytem Ω B is uncertain (presence of both the system parameter uncertainties and the model uncertainties). As we have explained in Section 1, these uncertainties are taken into account using the nonparametric approach of uncertainties which requires a reduced matrix order model (see refs. [11], [12], [14], [13]). Since, we have to represent the effects of this substructure on the nonlinear substructure Ω A through the coupling interface, it is natural to use the Craig Bampton method [1] in order to reduce the finite element model of subsytem Ω B . Finally to reduce the computational cost of the coupled system, subsystem Ω A is also reduced with the same technique.

Reduced mean matrix model for subsystem Ω B
The following change of coordinates is introduced , in which [Φ B ] is the (n B p × N B ) real matrix whose columns are the N B first elastic modes for the subsystem Ω B with a fixed coupling interface. Those

Reduced mean matrix model for subsystem Ω A
Using the same reduction method and introducing the elastic modes of the linear subsystem Ω A with fixed interface and without elastic stops, the R n A qvector q A (t) = (y A (t), U A c (t)) verifies the following matrix equation in which the matrices [

Transient dynamical response of the reduced non-linear computational model
Let n U = n A p + n B p + n c be the total number of DOF for the non-linear ) of the mean non-linear computational model is written as Let n q = N A + N B + n c . Then, using the coupling conditions defined by Eq. (6), the R nq -vector q(t) = (y A (t), y B (t), U c (t)) is solution of the reduced non-linear dynamical system with

Stochastic non-linear computational model including system uncertainties and identification
In this part, firstly the non-parametric probabilistic approach will be used to take into account both data uncertainties and model uncertainties in the reduced mean computational model of the linear subsystem Ω B of the computational model. This approach which has recently been introduced consists in replacing the mass, damping and stiffness matrices of reduced mean computational model by random matrices for which the probability distributions are explicitly given by the theory and for which a generator of independant realizations is known. Such an approach has been validated for many cases. For the details concerning the non-parametric probabilistic approach, one refers the reader, for instance, to references ( [11], [12], [14], [13]). In such an approach, the levels of uncertainties for the mass, damping and stiffness random matrices are defined by the dispersion parameters which are defined below. Secondly, these dispersion parameters will be identified using the maximum likelihood method. Finally, the stochastic non-linear computational model will be introduced and deduced from Section 2.2.3. It should be noted that (1) only the linear subsystem Ω B is assumed to be uncertain and (2) the mean non-linear subsystem Ω A is representative and consequently, that both data uncertainties and model uncertainties can be neglected. If such an assumption was not verified, then the non-parametric probabilistic approach of uncertainties could always be implemented without any difficulties in this non-linear subsystem (see for instance [2], [7]).

Stochastic linear subsystem Ω B modeling uncertainties
Therefore The dispersion parameter δ B P allows the level of uncertainites of the random matrix [P B ] to be controlled. It can be found in (see [11], [14]) an algebraic representation of random matrix [P B ] which allows independent realizations to be explicitly constructed in order to solve the random equations by the Monte Carlo method. For each random matrix, this random generator depends only on the mean value, on the dimension of the matrix and on the dispersion parameter. Such an approach is used in this paper.

Identification of the dispersion parameters
As explained in Section 3.1, the probability distributions of the random matrices (and then of the random generators) depend on the vector δ = (δ B M , δ B D , δ B K ) of the dispersion parameters which is identified using the measurements. The observation of the stochastic computational model is defined introducing the n B q × n B q random complex dynamic stiffness matrix [A B (ω)] of the linear subsystem Ω B written as Then the random condensed dynamical stiffness matrix [Z B (ω)] of the linear subsystem Ω B on the coupling interface is such that Taking into account the properties of the probabilistic model , it can be shown that, for all ω fixed in B, the random matrix [Z B (ω)] is invertible almost surely and the random variable J(δ) defined by is calculated using experimental data. The method used to identify vector δ is the maximum likelihood method (see for instance [8]) for the random variable J(δ) for which J exp is one realization. We then have to solve the following optimization problem in which δ opt is the identified value of δ.

Random transient dynamical response of the stochastic non-linear computational model
Using the probabilistic model defined in Section 3.1, the deterministic Eqs. (11) to (15) give the following stochastic non-linear computational model in which, for all fixed t, the vector-valued random variable Q(t) verifies

Identification of stochastic loads
The transient load is modelled by a stochastic process {F(t), t ∈ R}. Since all the degrees of freedom of the computational model are not excited by this stochastic load, we then introduce the usual projection operator Proj in order to extract the vector F(t) = Proj(F(t)) of the non zero random components of the random vector F(t). This equation can easily be inversed and yields F(t) = Lift(F(t)).

Construction of the stochastic loadF(t)
The stochastic load is modelled by a R m -valued Gaussian stationary centred second order stochastic process {F(t), t ∈ R} defined on a probability space (Θ ′ , T ′ , P ′ ) different from the probability space (Θ, T , P). In addition, it is assumed that the stochastic process is mean square continous on R, physically realizable (causal) and for which its matrix-valued autocorrelation function τ → [RF(τ )] = E{F(t + τ )F(t) T } is integrable on R. This stochastic process is then completely defined by its matrix-valued spectral density function [SF(ω)] = (2π) −1 R e −iωτ [RF(τ )] dτ which is a continuous and integrable function on R and which is in values in the set of all the positive (m × m) hermitian matrices. In addition, we will assume that for all ω in R, the matrix [SF(ω)] is with values in the set M + m (C) of all the positive definite (m × m) hermitian matrices. Since the stochastic process is assumed to be physically realizable, the matrix valued spectral density function must satisfy the following usual inequality ( [6], [10]) The numerical simulation of independent realizations {F(t, θ ′ ), t ∈ R} for θ ′ ∈ Θ ′ ( trajectories) can easily be generated by using adapted algorithms (see for instance [9], [5]).

Stochastic equation for simulation of responses
We have to identify the stochastic processF in presence of uncertainties in the linear subsystem Ω B . This identification consists in identifying the matrix-valued spectral density function [SF(ω)] which completely describes the stochastic process. This stochastic inverse problem is formulated as a stochastic optimization problem. Such an identification is performed using the stochastic equation deduced from Eqs. (23) to (26) with (15) in which the deterministic load F(t) is replaced by the stochastic load F(t). We then have to construct the R n U -valued stationary solution U s (t) = (U A s (t), U B s (t), U c s (t)) (corresponding to U(t)) which is written as in which subindex s is relative to the stationary solution and where the stationary stochastic process {Q s (t), t ∈ R} satisfies the stochastic equation in whichQ s (t) andQ s (t) are the mean-square first and second derivative of Q s (t). For the identification of [SF], for all t fixed, we introduce the R µ -valued random variable Z s (t) = (Z s,1 (t), ..., Z s,µ (t)) which represents the observations of the stochastic computational model made up of components of the vector-valued random response U s (t). Thus there exists a projection

Identification of the stochastic loads
The identification [SF] is performed in introducing a parametric representation of this function which is rewritten as in which C r ⊂ R νr is the admissible set of the parameter r with values in R νr where ν r is the number of unknown scalar parameters which have to be identified and where (ω, r) → [S(ω, r)] is a given function from R×R νr into M + m (C).
Therefore, the identification of the stochastic load {F(t), t ∈ R)} consists in identifying the R νr -valued vector r. Let {Z exp s (t) = (Z exp s,1 (t), ..., Z exp s,µ (t)), t ∈ R} be the R µ -valued stationary stochastic process which is measured for the manufactured real system and corresponding to the observation stochastic process {Z s (t), t ∈ R}. The matrix-valued spectral density function , ω ∈ R} of this stochastic process is estimated using the periodogram method. Then, the parameter r is estimated in minimizing the distance F dω between the matrix-valued spectral density function calculated with the stochastic computational model and the experimental matrix-valued spectral density function. We then have to solve the following optimization problem in which r opt is the identified value of the vector r.

Application
In this section a numerical simulation of a simple example is presented in order to validate the methodology developed in this paper.

Data for the experimental model
The measurement are generated by an experimental model which is made up of one linear subsystem and one non-linear subsystem. The linear subsystem is made up of four parallel beams fixed at their ends. The non-linear subsytem is made up of a beam also fixed at its ends, parallel to the other beams and with one transversal symmetric elastic stop (two identical stops, see Fig. 1). The five beams are linked by three transversal grids, each grid being modelled by four transverval springs (see Fig. 1). Therefore, the coupling interface between the two subsytems is composed of three points located in the neutral fiber of the beam of the non-linear subsystem. Each beam is modelled by eight Euler beam finite elements of equal lengthsfor which the DOF of the two nodes at the ends of the beam are locked. The twelve springs defining the the three tranversal grids are modelled by twelve spring elements. The two elastic stops are modeled by two springs. We are only interested in the y-direction displacements of the beam of the non-linear subsystem (see Fig. 1). Consequently, each beam has 14 DOF (y-translation and z-rotation). The total number of the free DOF for the linear subsystem is then 59 and the total number of the free DOF for the non-linear subsystem is then 14. The beam of the non-linear subsystem is exited by seven transversal

Data for the mean computational model
This part is devoted to the construction of a simplified mean computational model for the non-linear dynamical system described in Section 5.1. This simplified mean computational model will be used to identify the stochastic loads. It consists in modelling (1)  section of the equivalent beam for the linear subsystem is arbitrarily chosen and its Young's modulus and its mass density are identified so that the three first eigenfrequencies of this mean computational model are the same that the three first eigenfrequencies of the experimental model. Note that only the three first eigenfrequencies are correctly fitted and consequently, there are model uncertainties in this simplified mean computational model which are taken into account as explained in Section 3. It should be noted that the objective of this paper is not to construct an accurate mean computational model in order to exactly represent the experimental model, but to test the validity of the use of a simplified mean computational model in order to represent a much more complex system. After identification, the first three eigenfrequencies of the simplified mean computational model are 5.74 Hz , 15.3 Hz and 30.8 Hz which have to be compared to the first three eigenfrequencies 5.78 Hz , 15.9 Hz and 31.1 Hz of the experimental model.

Comparison between the dynamical responses of
the experimental model and of the mean computational model for the same given stochastic load.
In this section, it is assumed that the stochastic load is given and is the same for the experimental model and for the simplified mean computational model. Then, for the two models, the stationary stochastic responses are calculated in the time interval [0, 220] s using an explicit Euler integration scheme. Let P obs be the point of the non-linear subsystem located at the impact point of the elastic stops. The power spectral density functions of the stochastic y-displacement and of the stochastic z-rotation in point P obs (see Fig. 3

System uncertainties modelling and dispersion parameter identification.
The non-parametric probabilistic approach of model uncertainties introduced in Section 3.1 is used for stiffness part of the linear subsystem of the simplified mean computational model. We then have to identify the dispersion parameter δ = (δ B K ). Note that the identification procedure which is proposed is independent of the stochastic loads. The estimation of the probability density function in Eq. (22) is carried out with 200 realizations for the Monte Carlo simulation. Figure 4 shows the likelihood function calculated using Eq.

Case of an unknown stochastic load and its identification.
In this section, the responses of the experimental model are given (those constructed in Section 5.3) and the stochastic loadF(t) is assumed to be unknown and has to be identified using the uncertain simplified computational model, that is to say the stochastic simplified computational model for which the dispersion parameter has been identified in Section 5.4. We begin defining a model as simple as possible for the stochastic loadF(t) introduced in Section 4.1. We have then chosen to modelF(t) as {F(t) = (T (t), M(t)), t ∈ R} in which T (t) is a y-force and M(t) is a z-moment applied to the middle of the beam of the non-linear subsystem (see Fig. 5). This force and this in which the admissible set C r = {r = (r 1 , r 2 ); r 1 > 0, r 2 > 0}. This vector r is identified using the trial method to solve the optimization problem defined by Eq. (33). Such a method consists in calculating the cost function D(r) for 100 values of the vector r. Figure 6 shows the graph of the function r → log 10 (D(r)) which allows the optimal value r opt to be determined. The confidence region associated with a probability level P c = 0.95 of the reponse of the stochastic simplified computational model on which the identified stochatic load is applied can then be estimated. The comparison between the experimental responses with the responses constructed with the stochastic simplified computational model is given in Fig. 7. This figure displays the confidence region of the power spectral density function of the stochastic y-displacement and the stochastic z-rotation for point P obs .  Figure 6: Graph of the cost function (r 1 , r 2 ) → log 10 (D(r 1 , r 2 )).

Conclusions
We have presented a methodology and its validation to perform the identification of a stochastic loads applied to a complex non-linear dynamical system for which a few measurements of its responses are available. To carry out this identification, a simplified computational model of the real system is introduced. Since such a simplified computational model induces model uncertainties, a probabilistic model of these uncertainties is introduced in the simplified computational model. The identification of the stochastic loads is then performed using this stochastic computational model which takes into account model uncertainties and consequently, we have validated a method to perform a robust identification with respect to model uncertainties. It should be noted that the non-linear dynamical system used for this validation is representative of real industrial systems and then validates the methodology proposed.