# Identification of Stochastic Loads Applied to

Authors: C. Soize; A. Batou
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

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


---

## Body

## 1. Introduction

This paper is
devoted to the identification of a stochastic load applied to a nonlinear
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, 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 nonlinear 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 nonlinear dynamical system is very complex, the computational model
developed cannot exactly represent the complexity of the system. Consequently,
the identification is not performed using a computational model which 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 nonlinear 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 nonlinear 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.Section2 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.


## 2. Mean Computational Model

LetΩ be the domain of the dynamical system having a
nonlinear behavior 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 ℝ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 x1,…,xκ in the boundary ΓA of ∂ΩA.
In addition, the subsystem occupying the subdomain ΩA is fixed on the part Γ0A of its boundary ∂ΩA.
The outward unit normal of ∂ΩA is denoted by
nA.
The subdomain ΩB is constituted of a three-dimensional linear
viscoelastic medium with instantaneous memory, fixed on the part Γ0B of its boundary ∂ΩB.
The outward unit normal of ∂ΩB is denoted by
nB.
Consequently, each uncoupled subsystem ΩA and ΩB does not have rigid body displacement. These
two subsystems are coupled on the common coupling interface ΓC. One then has ∂ΩA=Γ0A∪ΓA∪ΓC and ∂ΩB=Γ0B∪ΓB∪ΓC.
We are interested in constructing the stationary random responses of the
nonlinear 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 ℝ.


### 2.1. Mean Boundary Value Problems

#### 2.1.1. Mean Boundary Value Problems for the Linear SubsystemΩB

Letx=(x1,x2,x3) be the cartesian coordinates and uB(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 boundary ΓB are denoted by fvolB(x,t) and fsurfB(x,t).
The stress tensor σB(x) is written as σijB(x)=aijkhB(x)εkhB(uB)+bijkhB(x)εkhB(u˙B) where εkhB(uB)=(∂ukB/∂xh+∂ukB/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aB(x) and bB(x) verify the usual properties of symmetry and
positiveness [1].
Then, the displacement field uB(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.1)ρBu¨iB−∂σijB∂xj=fvol,iBin⁢ΩB,σijBnjB=fsurf,iBon⁢ΓB,σijBnjB=fcoupl,iBon⁢ΓC,uiB=0on⁢Γ0B,in which a dot means the partial
time derivative and a double dot means the double partial time derivative,
where fcouplB=(fcoupl,1B,…,fcoupl,3B) 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.


#### 2.1.2. Mean Boundary Value Problems for the Nonlinear SubsystemΩA

LetuA(x,t) be the displacement field for the nonlinear
subsystem ΩA at time t.
The external prescribed volumetric and surface forces fields applied to ΩA and to its boundary ΓA are denoted by fvolA(x,t) and fsurfA(x,t).
Since ΩA is occupied by a linear viscoelastic material
with instantaneous memory, the stress tensor σA(x) is written as σijA(x)=aijkhA(x)εkhA(uA)+bijkhA(x)εkhA(u˙A) where εkhA(uA)=(∂ukA/∂xh+∂ukA/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aA(x) and bA(x) verify, as above, the usual properties of
symmetry and positiveness. Then, the displacement field uA(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.2)ρAu¨iA−∂σijA∂xj=fvol,iAin⁢ΩA,σijAnjA=fsurf,iA−∑k=1κfiNL,k(u(xk,t))δ0(x−xk)on⁢ΓA,σijAnjA=fcoupl,iAon⁢ΓC,uiA=0on⁢Γ0A,in which fcouplA=(fcoupl,1A,…,fcoupl,3A) is the forces induced by the subsystem ΩB on ΩA via the coupling interface ΓC.
The forces—fNL,k(u(xk,t)) represent the actions exerted by the elastic
stop located at point xk on the subsystem ΩA and δ0(x−xk) is the surface Dirac measure such that, for
all continuous function g defined on ΓA,
one has ∫ΓAδ0(x−xk)g(x)ds(x)=g(xk).
The parameter ρA(x) is the mass density for the subsystem ΩA.


#### 2.1.3. Interface Conditions for the Coupling ofΩA with ΩB

The coupling
conditions onΓC are written as(2.3)uA=uBon⁢ΓC,fcouplA+fcouplB=0on⁢ΓC.


### 2.2. Mean Finite Element Model

The mean
finite element method [2]
is applied to the variational formulation of the boundary value problems
defined in Section 2.1.


#### 2.2.1. Mean Finite Element Model for SubsystemΩB

TheℝnB vector 𝕌B(t) of the nB DOF of the subsystem ΩB is written as 𝕌B(t)=(𝕌pB(t),𝕌cB(t)),
where 𝕌pB(t) is the ℝnpB vector of the npB internal DOF and where 𝕌cB(t) is the ℝncB vector function of the ncB coupling DOF on the interface. From (2.1), it
can be deduced that the mean finite element model of subsystem ΩB is written as(2.4)[M¯B]U¨B(t)+[D¯B]U˙B(t)+[K¯B]UB(t)=FB(t)+FcouplB(t),in which [𝕄¯B],[𝔻¯B],
and [𝕂¯B] are, respectively, the positive-definite
symmetric real positive (nB×nB) mass, damping, and stiffness matrices. The ℝnB vectors 𝔽B(t) and 𝔽couplB(t) of the external forces and of the coupling
forces are written as 𝔽B(t)=(𝔽pB(t),0) and 𝔽couplB(t)=(0,𝔽cB(t)).


#### 2.2.2. Mean Finite Element Model for SubsystemΩA

Similarly to Section2.2.1, the ℝnA vector 𝕌A(t) of the nA DOF of the subsystem ΩA is written as 𝕌A(t)=(𝕌pA(t),𝕌cA(t)),
where 𝕌pA(t) is the ℝnpA vector of the npA internal DOF and where 𝕌cA(t) is the ℝncA-valued function of the ncA coupling DOF. From (2.2), it can be deduced
that the mean finite element model of subsystem ΩA is written as(2.5)[M¯A]U¨A(t)+[D¯A]U˙A(t)+[K¯A]UA(t)+FNL(UA(t))=FA(t)+FcouplA(t),in which [𝕄¯A],[𝔻¯A],
and [𝕂¯A] are, respectively, the positive-definite
symmetric real positive (nA×nA) mass, damping, and stiffness matrices. The ℝnA vectors 𝔽A(t), 𝔽couplA(t),
and 𝔽NL(𝕌A(t)) of the external forces, of the coupling forces
and of the nonlinear forces, are written as 𝔽A(t)=(𝔽pA(t),0),𝔽couplA(t)=(0,𝔽cA(t)),
and 𝔽NL(𝕌A(t))=(𝔽pNL(𝕌A(t)),0).


#### 2.2.3. Interface Conditions for the Coupling ofΩA with ΩB

The finite
element discretization of the interface conditions defined by (2.3)
yields(2.6)UcA(t)=UcB(t)on⁢ΓC,FcouplA(t)+FcouplB(t)=0on⁢ΓC.


### 2.3. Reduced Mean Matrix Model

The continuous linear subsystemΩB (linear dynamical subsystem) contains elastic
modes in the frequency band of analysis. In addition, the computational model
of the continuous linear subsystem Ω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
[3–6]). 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 [7] 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.


#### 2.3.1. Reduced Mean Matrix Model for SubsystemΩB

The following change of coordinates is introduced:(2.7)[UpB(t)UcB(t)]=[HB][yB(t)UcB(t)],[HB]=[[ΦB][SB][0][I]],in which [ΦB] is the (npB×NB) real matrix whose columns are the NB first elastic modes for the subsystem ΩB with a fixed coupling interface. Those modes (ϕ1,…,ϕNB) are associated with the NB first eigenvalues 0<ω12≤⋯≤ωNB2 such that(2.8)[K¯ppB]ϕB=ω2[M¯ppB]ϕB,where [𝕂¯ppB] and [𝕄¯ppB] are the internal DOF blocks of the matrices [𝕂¯B] and [𝕄¯B],
where [SB]=[𝕂¯ppB]−1[𝕂¯pcB] is an (npB×ncB) matrix, where [I] is the (ncB×ncB) unity matrix and where yB(t) is an ℝNB-vector. Let nqB=NB+ncB.
Then, the ℝnqB vector qB(t)=(yB(t),𝕌cB(t)) is a solution of the reduced mean
computational model(2.9)[M¯B]q¨B(t)+[D¯B]q˙B(t)+[K¯B]qB(t)=[HB]TFB(t)+[HB]TFcouplB(t),in which the matrices [M¯B]=[HB]T[𝕄¯B][HB],[D¯B]=[HB]T[𝔻¯B][HB],
and [K¯B]=[HB]T[𝕂¯B][HB] are positive-definite symmetric real (nqB×nqB) matrices.


#### 2.3.2. 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 ℝnqA-vector qA(t)=(yA(t),𝕌cA(t)) verifies the following matrix equation:(2.10)[M¯A]q¨A(t)+[D¯A]q˙A(t)+[K¯A]qA(t)+[HA]TFNL([HA]qA(t))=[HA]TFA(t)+[HA]TFcouplA(t),in which the matrices [M¯A]=[HA]T[𝕄¯A][HA],[D¯A]=[HA]T[𝔻¯A][HA],
and [K¯A]=[HA]T[𝕂¯A][HA] are positive-definite symmetric real (nqA×nqA) matrices.


#### 2.3.3. Transient Dynamical Response of the Reduced Nonlinear Computational Model

Letn𝕌=npA+npB+nc be the total number of DOF for the nonlinear
computational model. The ℝn𝕌-vector 𝕌(t)=(𝕌pA(t),𝕌pB(t),𝕌c(t)) of the mean nonlinear computational model is
written as(2.11)[U(t)]=[H][q(t)],in which the matrix [H] is constructed by the assemblage of [HA] and [HB].
Let nq=NA+NB+nc.
Then, using the coupling conditions defined by (2.6), the ℝnq-vector q(t)=(yA(t),yB(t),𝕌c(t)) is a solution of the reduced nonlinear
dynamical system(2.12)[M¯]q¨(t)+[D¯]q˙(t)+[K¯]q(t)+ℱNL(q(t))=ℱ(t),with(2.13)[M¯]=[M¯yyA0M¯ycA0M¯yyBM¯ycBM¯cyAM¯cyBM¯ccA+M¯ccB],[D¯]=[D¯yyA0D¯ycA0D¯yyBD¯ycBD¯cyAD¯cyBD¯ccA+D¯ccB],(2.14)[K¯]=[K¯yyA0K¯ycA0K¯yyBK¯ycBK¯cyAK¯cyBK¯ccA+K¯ccB],(2.15)ℱNL(q(t))=[[ΦA]TFpNL([HA]qA(t))0[SA]TFpNL([HA]qA(t))],(2.16)ℱ(t)=[[ΦA]TFpA(t)[ΦB]TFpB(t)[SA]TFpA(t)+[SB]TFpB(t)].


## 2.1. Mean Boundary Value Problems

### 2.1.1. Mean Boundary Value Problems for the Linear SubsystemΩB

Letx=(x1,x2,x3) be the cartesian coordinates and uB(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 boundary ΓB are denoted by fvolB(x,t) and fsurfB(x,t).
The stress tensor σB(x) is written as σijB(x)=aijkhB(x)εkhB(uB)+bijkhB(x)εkhB(u˙B) where εkhB(uB)=(∂ukB/∂xh+∂ukB/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aB(x) and bB(x) verify the usual properties of symmetry and
positiveness [1].
Then, the displacement field uB(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.1)ρBu¨iB−∂σijB∂xj=fvol,iBin⁢ΩB,σijBnjB=fsurf,iBon⁢ΓB,σijBnjB=fcoupl,iBon⁢ΓC,uiB=0on⁢Γ0B,in which a dot means the partial
time derivative and a double dot means the double partial time derivative,
where fcouplB=(fcoupl,1B,…,fcoupl,3B) 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.


### 2.1.2. Mean Boundary Value Problems for the Nonlinear SubsystemΩA

LetuA(x,t) be the displacement field for the nonlinear
subsystem ΩA at time t.
The external prescribed volumetric and surface forces fields applied to ΩA and to its boundary ΓA are denoted by fvolA(x,t) and fsurfA(x,t).
Since ΩA is occupied by a linear viscoelastic material
with instantaneous memory, the stress tensor σA(x) is written as σijA(x)=aijkhA(x)εkhA(uA)+bijkhA(x)εkhA(u˙A) where εkhA(uA)=(∂ukA/∂xh+∂ukA/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aA(x) and bA(x) verify, as above, the usual properties of
symmetry and positiveness. Then, the displacement field uA(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.2)ρAu¨iA−∂σijA∂xj=fvol,iAin⁢ΩA,σijAnjA=fsurf,iA−∑k=1κfiNL,k(u(xk,t))δ0(x−xk)on⁢ΓA,σijAnjA=fcoupl,iAon⁢ΓC,uiA=0on⁢Γ0A,in which fcouplA=(fcoupl,1A,…,fcoupl,3A) is the forces induced by the subsystem ΩB on ΩA via the coupling interface ΓC.
The forces—fNL,k(u(xk,t)) represent the actions exerted by the elastic
stop located at point xk on the subsystem ΩA and δ0(x−xk) is the surface Dirac measure such that, for
all continuous function g defined on ΓA,
one has ∫ΓAδ0(x−xk)g(x)ds(x)=g(xk).
The parameter ρA(x) is the mass density for the subsystem ΩA.


### 2.1.3. Interface Conditions for the Coupling ofΩA with ΩB

The coupling
conditions onΓC are written as(2.3)uA=uBon⁢ΓC,fcouplA+fcouplB=0on⁢ΓC.


## 2.1.1. Mean Boundary Value Problems for the Linear SubsystemΩB

Letx=(x1,x2,x3) be the cartesian coordinates and uB(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 boundary ΓB are denoted by fvolB(x,t) and fsurfB(x,t).
The stress tensor σB(x) is written as σijB(x)=aijkhB(x)εkhB(uB)+bijkhB(x)εkhB(u˙B) where εkhB(uB)=(∂ukB/∂xh+∂ukB/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aB(x) and bB(x) verify the usual properties of symmetry and
positiveness [1].
Then, the displacement field uB(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.1)ρBu¨iB−∂σijB∂xj=fvol,iBin⁢ΩB,σijBnjB=fsurf,iBon⁢ΓB,σijBnjB=fcoupl,iBon⁢ΓC,uiB=0on⁢Γ0B,in which a dot means the partial
time derivative and a double dot means the double partial time derivative,
where fcouplB=(fcoupl,1B,…,fcoupl,3B) 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.


## 2.1.2. Mean Boundary Value Problems for the Nonlinear SubsystemΩA

LetuA(x,t) be the displacement field for the nonlinear
subsystem ΩA at time t.
The external prescribed volumetric and surface forces fields applied to ΩA and to its boundary ΓA are denoted by fvolA(x,t) and fsurfA(x,t).
Since ΩA is occupied by a linear viscoelastic material
with instantaneous memory, the stress tensor σA(x) is written as σijA(x)=aijkhA(x)εkhA(uA)+bijkhA(x)εkhA(u˙A) where εkhA(uA)=(∂ukA/∂xh+∂ukA/∂xh)/2 is the linearized strain tensor. The
fourth-order tensors aA(x) and bA(x) verify, as above, the usual properties of
symmetry and positiveness. Then, the displacement field uA(t) verifies, for all t∈ℝ and for i=1,2,3,
the mean boundary value problem(2.2)ρAu¨iA−∂σijA∂xj=fvol,iAin⁢ΩA,σijAnjA=fsurf,iA−∑k=1κfiNL,k(u(xk,t))δ0(x−xk)on⁢ΓA,σijAnjA=fcoupl,iAon⁢ΓC,uiA=0on⁢Γ0A,in which fcouplA=(fcoupl,1A,…,fcoupl,3A) is the forces induced by the subsystem ΩB on ΩA via the coupling interface ΓC.
The forces—fNL,k(u(xk,t)) represent the actions exerted by the elastic
stop located at point xk on the subsystem ΩA and δ0(x−xk) is the surface Dirac measure such that, for
all continuous function g defined on ΓA,
one has ∫ΓAδ0(x−xk)g(x)ds(x)=g(xk).
The parameter ρA(x) is the mass density for the subsystem ΩA.


## 2.1.3. Interface Conditions for the Coupling ofΩA with ΩB

The coupling
conditions onΓC are written as(2.3)uA=uBon⁢ΓC,fcouplA+fcouplB=0on⁢ΓC.


## 2.2. Mean Finite Element Model

The mean
finite element method [2]
is applied to the variational formulation of the boundary value problems
defined in Section 2.1.


### 2.2.1. Mean Finite Element Model for SubsystemΩB

TheℝnB vector 𝕌B(t) of the nB DOF of the subsystem ΩB is written as 𝕌B(t)=(𝕌pB(t),𝕌cB(t)),
where 𝕌pB(t) is the ℝnpB vector of the npB internal DOF and where 𝕌cB(t) is the ℝncB vector function of the ncB coupling DOF on the interface. From (2.1), it
can be deduced that the mean finite element model of subsystem ΩB is written as(2.4)[M¯B]U¨B(t)+[D¯B]U˙B(t)+[K¯B]UB(t)=FB(t)+FcouplB(t),in which [𝕄¯B],[𝔻¯B],
and [𝕂¯B] are, respectively, the positive-definite
symmetric real positive (nB×nB) mass, damping, and stiffness matrices. The ℝnB vectors 𝔽B(t) and 𝔽couplB(t) of the external forces and of the coupling
forces are written as 𝔽B(t)=(𝔽pB(t),0) and 𝔽couplB(t)=(0,𝔽cB(t)).


### 2.2.2. Mean Finite Element Model for SubsystemΩA

Similarly to Section2.2.1, the ℝnA vector 𝕌A(t) of the nA DOF of the subsystem ΩA is written as 𝕌A(t)=(𝕌pA(t),𝕌cA(t)),
where 𝕌pA(t) is the ℝnpA vector of the npA internal DOF and where 𝕌cA(t) is the ℝncA-valued function of the ncA coupling DOF. From (2.2), it can be deduced
that the mean finite element model of subsystem ΩA is written as(2.5)[M¯A]U¨A(t)+[D¯A]U˙A(t)+[K¯A]UA(t)+FNL(UA(t))=FA(t)+FcouplA(t),in which [𝕄¯A],[𝔻¯A],
and [𝕂¯A] are, respectively, the positive-definite
symmetric real positive (nA×nA) mass, damping, and stiffness matrices. The ℝnA vectors 𝔽A(t), 𝔽couplA(t),
and 𝔽NL(𝕌A(t)) of the external forces, of the coupling forces
and of the nonlinear forces, are written as 𝔽A(t)=(𝔽pA(t),0),𝔽couplA(t)=(0,𝔽cA(t)),
and 𝔽NL(𝕌A(t))=(𝔽pNL(𝕌A(t)),0).


### 2.2.3. Interface Conditions for the Coupling ofΩA with ΩB

The finite
element discretization of the interface conditions defined by (2.3)
yields(2.6)UcA(t)=UcB(t)on⁢ΓC,FcouplA(t)+FcouplB(t)=0on⁢ΓC.


## 2.2.1. Mean Finite Element Model for SubsystemΩB

TheℝnB vector 𝕌B(t) of the nB DOF of the subsystem ΩB is written as 𝕌B(t)=(𝕌pB(t),𝕌cB(t)),
where 𝕌pB(t) is the ℝnpB vector of the npB internal DOF and where 𝕌cB(t) is the ℝncB vector function of the ncB coupling DOF on the interface. From (2.1), it
can be deduced that the mean finite element model of subsystem ΩB is written as(2.4)[M¯B]U¨B(t)+[D¯B]U˙B(t)+[K¯B]UB(t)=FB(t)+FcouplB(t),in which [𝕄¯B],[𝔻¯B],
and [𝕂¯B] are, respectively, the positive-definite
symmetric real positive (nB×nB) mass, damping, and stiffness matrices. The ℝnB vectors 𝔽B(t) and 𝔽couplB(t) of the external forces and of the coupling
forces are written as 𝔽B(t)=(𝔽pB(t),0) and 𝔽couplB(t)=(0,𝔽cB(t)).


## 2.2.2. Mean Finite Element Model for SubsystemΩA

Similarly to Section2.2.1, the ℝnA vector 𝕌A(t) of the nA DOF of the subsystem ΩA is written as 𝕌A(t)=(𝕌pA(t),𝕌cA(t)),
where 𝕌pA(t) is the ℝnpA vector of the npA internal DOF and where 𝕌cA(t) is the ℝncA-valued function of the ncA coupling DOF. From (2.2), it can be deduced
that the mean finite element model of subsystem ΩA is written as(2.5)[M¯A]U¨A(t)+[D¯A]U˙A(t)+[K¯A]UA(t)+FNL(UA(t))=FA(t)+FcouplA(t),in which [𝕄¯A],[𝔻¯A],
and [𝕂¯A] are, respectively, the positive-definite
symmetric real positive (nA×nA) mass, damping, and stiffness matrices. The ℝnA vectors 𝔽A(t), 𝔽couplA(t),
and 𝔽NL(𝕌A(t)) of the external forces, of the coupling forces
and of the nonlinear forces, are written as 𝔽A(t)=(𝔽pA(t),0),𝔽couplA(t)=(0,𝔽cA(t)),
and 𝔽NL(𝕌A(t))=(𝔽pNL(𝕌A(t)),0).


## 2.2.3. Interface Conditions for the Coupling ofΩA with ΩB

The finite
element discretization of the interface conditions defined by (2.3)
yields(2.6)UcA(t)=UcB(t)on⁢ΓC,FcouplA(t)+FcouplB(t)=0on⁢ΓC.


## 2.3. Reduced Mean Matrix Model

The continuous linear subsystemΩB (linear dynamical subsystem) contains elastic
modes in the frequency band of analysis. In addition, the computational model
of the continuous linear subsystem Ω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
[3–6]). 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 [7] 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.


### 2.3.1. Reduced Mean Matrix Model for SubsystemΩB

The following change of coordinates is introduced:(2.7)[UpB(t)UcB(t)]=[HB][yB(t)UcB(t)],[HB]=[[ΦB][SB][0][I]],in which [ΦB] is the (npB×NB) real matrix whose columns are the NB first elastic modes for the subsystem ΩB with a fixed coupling interface. Those modes (ϕ1,…,ϕNB) are associated with the NB first eigenvalues 0<ω12≤⋯≤ωNB2 such that(2.8)[K¯ppB]ϕB=ω2[M¯ppB]ϕB,where [𝕂¯ppB] and [𝕄¯ppB] are the internal DOF blocks of the matrices [𝕂¯B] and [𝕄¯B],
where [SB]=[𝕂¯ppB]−1[𝕂¯pcB] is an (npB×ncB) matrix, where [I] is the (ncB×ncB) unity matrix and where yB(t) is an ℝNB-vector. Let nqB=NB+ncB.
Then, the ℝnqB vector qB(t)=(yB(t),𝕌cB(t)) is a solution of the reduced mean
computational model(2.9)[M¯B]q¨B(t)+[D¯B]q˙B(t)+[K¯B]qB(t)=[HB]TFB(t)+[HB]TFcouplB(t),in which the matrices [M¯B]=[HB]T[𝕄¯B][HB],[D¯B]=[HB]T[𝔻¯B][HB],
and [K¯B]=[HB]T[𝕂¯B][HB] are positive-definite symmetric real (nqB×nqB) matrices.


### 2.3.2. 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 ℝnqA-vector qA(t)=(yA(t),𝕌cA(t)) verifies the following matrix equation:(2.10)[M¯A]q¨A(t)+[D¯A]q˙A(t)+[K¯A]qA(t)+[HA]TFNL([HA]qA(t))=[HA]TFA(t)+[HA]TFcouplA(t),in which the matrices [M¯A]=[HA]T[𝕄¯A][HA],[D¯A]=[HA]T[𝔻¯A][HA],
and [K¯A]=[HA]T[𝕂¯A][HA] are positive-definite symmetric real (nqA×nqA) matrices.


### 2.3.3. Transient Dynamical Response of the Reduced Nonlinear Computational Model

Letn𝕌=npA+npB+nc be the total number of DOF for the nonlinear
computational model. The ℝn𝕌-vector 𝕌(t)=(𝕌pA(t),𝕌pB(t),𝕌c(t)) of the mean nonlinear computational model is
written as(2.11)[U(t)]=[H][q(t)],in which the matrix [H] is constructed by the assemblage of [HA] and [HB].
Let nq=NA+NB+nc.
Then, using the coupling conditions defined by (2.6), the ℝnq-vector q(t)=(yA(t),yB(t),𝕌c(t)) is a solution of the reduced nonlinear
dynamical system(2.12)[M¯]q¨(t)+[D¯]q˙(t)+[K¯]q(t)+ℱNL(q(t))=ℱ(t),with(2.13)[M¯]=[M¯yyA0M¯ycA0M¯yyBM¯ycBM¯cyAM¯cyBM¯ccA+M¯ccB],[D¯]=[D¯yyA0D¯ycA0D¯yyBD¯ycBD¯cyAD¯cyBD¯ccA+D¯ccB],(2.14)[K¯]=[K¯yyA0K¯ycA0K¯yyBK¯ycBK¯cyAK¯cyBK¯ccA+K¯ccB],(2.15)ℱNL(q(t))=[[ΦA]TFpNL([HA]qA(t))0[SA]TFpNL([HA]qA(t))],(2.16)ℱ(t)=[[ΦA]TFpA(t)[ΦB]TFpB(t)[SA]TFpA(t)+[SB]TFpB(t)].


## 2.3.1. Reduced Mean Matrix Model for SubsystemΩB

The following change of coordinates is introduced:(2.7)[UpB(t)UcB(t)]=[HB][yB(t)UcB(t)],[HB]=[[ΦB][SB][0][I]],in which [ΦB] is the (npB×NB) real matrix whose columns are the NB first elastic modes for the subsystem ΩB with a fixed coupling interface. Those modes (ϕ1,…,ϕNB) are associated with the NB first eigenvalues 0<ω12≤⋯≤ωNB2 such that(2.8)[K¯ppB]ϕB=ω2[M¯ppB]ϕB,where [𝕂¯ppB] and [𝕄¯ppB] are the internal DOF blocks of the matrices [𝕂¯B] and [𝕄¯B],
where [SB]=[𝕂¯ppB]−1[𝕂¯pcB] is an (npB×ncB) matrix, where [I] is the (ncB×ncB) unity matrix and where yB(t) is an ℝNB-vector. Let nqB=NB+ncB.
Then, the ℝnqB vector qB(t)=(yB(t),𝕌cB(t)) is a solution of the reduced mean
computational model(2.9)[M¯B]q¨B(t)+[D¯B]q˙B(t)+[K¯B]qB(t)=[HB]TFB(t)+[HB]TFcouplB(t),in which the matrices [M¯B]=[HB]T[𝕄¯B][HB],[D¯B]=[HB]T[𝔻¯B][HB],
and [K¯B]=[HB]T[𝕂¯B][HB] are positive-definite symmetric real (nqB×nqB) matrices.


## 2.3.2. 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 ℝnqA-vector qA(t)=(yA(t),𝕌cA(t)) verifies the following matrix equation:(2.10)[M¯A]q¨A(t)+[D¯A]q˙A(t)+[K¯A]qA(t)+[HA]TFNL([HA]qA(t))=[HA]TFA(t)+[HA]TFcouplA(t),in which the matrices [M¯A]=[HA]T[𝕄¯A][HA],[D¯A]=[HA]T[𝔻¯A][HA],
and [K¯A]=[HA]T[𝕂¯A][HA] are positive-definite symmetric real (nqA×nqA) matrices.


## 2.3.3. Transient Dynamical Response of the Reduced Nonlinear Computational Model

Letn𝕌=npA+npB+nc be the total number of DOF for the nonlinear
computational model. The ℝn𝕌-vector 𝕌(t)=(𝕌pA(t),𝕌pB(t),𝕌c(t)) of the mean nonlinear computational model is
written as(2.11)[U(t)]=[H][q(t)],in which the matrix [H] is constructed by the assemblage of [HA] and [HB].
Let nq=NA+NB+nc.
Then, using the coupling conditions defined by (2.6), the ℝnq-vector q(t)=(yA(t),yB(t),𝕌c(t)) is a solution of the reduced nonlinear
dynamical system(2.12)[M¯]q¨(t)+[D¯]q˙(t)+[K¯]q(t)+ℱNL(q(t))=ℱ(t),with(2.13)[M¯]=[M¯yyA0M¯ycA0M¯yyBM¯ycBM¯cyAM¯cyBM¯ccA+M¯ccB],[D¯]=[D¯yyA0D¯ycA0D¯yyBD¯ycBD¯cyAD¯cyBD¯ccA+D¯ccB],(2.14)[K¯]=[K¯yyA0K¯ycA0K¯yyBK¯ycBK¯cyAK¯cyBK¯ccA+K¯ccB],(2.15)ℱNL(q(t))=[[ΦA]TFpNL([HA]qA(t))0[SA]TFpNL([HA]qA(t))],(2.16)ℱ(t)=[[ΦA]TFpA(t)[ΦB]TFpB(t)[SA]TFpA(t)+[SB]TFpB(t)].


## 3. Stochastic Nonlinear Computational Model Including System Uncertainties and Identification

In this part, firstly the nonparametric 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 independent realizations is known. Such an approach has
been validated for many cases. For the details concerning the nonparametric probabilistic
approach, one refers the reader, for instance, to [3–6]. 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 nonlinear computational model will be
introduced and deduced from Section 2.2.3. It should be noted that (2.1) only the
linear subsystem ΩB is assumed to be uncertain and (2.2) the mean
nonlinear 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 nonparametric probabilistic approach of
uncertainties could always be implemented without any difficulties in this
nonlinear subsystem (see for instance [8, 9]).


### 3.1. Stochastic Linear SubsystemΩB Modeling Uncertainties

Therefore, the matrices[M¯B],[D¯B],
and [K¯B] of the reduced mean computational model are
replaced by the random matrices [MB],[DB],
and [KB] defined on a probability space (Θ,𝒯,𝒫) and such that(3.1)∀θ∈Θ,⁢[MB(θ)],[DB(θ)],[KB(θ)]∈MnqB+(ℝ),E{[MB]}=[M¯B],E{[DB]}=[D¯B],E{[KB]}=[K¯B],E{∥[MB]−1∥F2}<+∞,E{∥[DB]−1∥F2}<+∞,E{∥[KB]−1∥F2}<+∞,in which 𝕄n+(ℝ) is the set of all the positive-definite
symmetric (n×n) matrices, where E{⋅} is the mathematical expectation and where ∥⋅∥F is the Frobenius norm such that ∥A∥F2=tr{[A]∗[A]} with [A]∗=[A¯]T, [A¯] is the conjugate of [A] and tr is the trace for matrices. Let [PB] be the random matrix denoting [MB],[DB],
or [KB].
The probability distribution of the random matrix [PB] depends on the dispersion parameter δPB related to the coefficient of variation δPB of the random matrix [PB] by the equation(3.2)(δPB)2=E{∥[PB]−[P¯B]∥F2}∥[P¯B]∥F2=(δPB)2n+1(1+(tr[PB])2tr([PB]2)).The dispersion parameter δPB allows the level of uncertainites of the
random matrix [PB] to be controlled. It can be
found in [3, 5] an algebraic
representation of random
matrix [PB] 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.


### 3.2. Identification of the Dispersion Parameters

As explained in Section3.1, the probability
distributions of the random matrices (and then of the random generators) depend
on the vector δ=(δMB,δDB,δKB) of the dispersion parameters which is
identified using the measurements. The observation of the stochastic
computational model is defined introducing the nqB×nqB random complex dynamic stiffness matrix [AB(ω)] of the linear subsystem ΩB written as(3.3)[AB(ω)]=−ω2[MB]+iω[DB]+[KB].Then the random condensed
dynamical stiffness matrix [ZB(ω)] of the linear subsystem ΩB on the coupling interface is such that [ZB(ω)]=[AccB(ω)]−[AcyB(ω)][AyyB(ω)]−1[AycB(ω)].
Taking into account the properties of the probabilistic model, it can be shown
that, for all ω fixed in ℬ,
the random matrix [ZB(ω)] is invertible almost surely and the random
variable J(δ) defined by(3.4)J(δ)=∫ℬ∥[ZB(ω)]−1∥F2dωexists and has a finite mean
value. This random variable gives a measure over ℬ of the dynamical effects of subsystem ΩB on the subsystem ΩA at the coupling interface. It should be noted
that the random variable J(δ) depends on δ because the probability distributions of the
random matrices [MB],[DB],
and [KB] depend on δ.
Let x↦pJ(x,δ) be the probabilty density function of the
random variable J(δ) with respect to dx.
For any x fixed in [0,+∞[ and for any value of the vector δ belonging to the admissible set 𝒞ad of the dispersion parameters, the value pJ(x,δ) of the probability density function is
estimated by using the above probabilistic model and the Monte Carlo
simulation. The corresponding deterministic experimental value Jexp of J(δ) is calculated using experimental data. The
method used to identify vector δ is the maximum likelihood method (see for
instance [10]) for the
random variable J(δ) for which Jexp is one realization. We then have to solve the
following optimization problem:(3.5)δopt=arg⁡max⁡δ∈Cad(pJ(Jexp;δ)),in which δopt is the identified value of δ.


### 3.3. Random Transient Dynamical Response of the Stochastic Nonlinear Computational Model

Using the
probabilistic model defined in Section3.1, the deterministic (2.11) to (2.15)
give the following stochastic nonlinear computational model:(3.6)[U(t)]=[H][Q(t)],in which, for all fixed t,
the vector-valued random variable Q(t) verifies(3.7)[M]Q¨(t)+[D]Q˙(t)+[K]Q(t)+ℱNL(Q(t))=ℱ(t),and where the random matrices [M],[D],
and [K] are written as(3.8)[M]=[M¯yyA0M¯ycA0MyyBMycBM¯cyAMcyBM¯ccA+MccB],[D]=[D¯yyA0D¯ycA0DyyBDycBD¯cyADcyBD¯ccA+DccB],[K]=[K¯yyA0K¯ycA0KyyBKycBK¯cyAKcyBK¯ccA+KccB].


## 3.1. Stochastic Linear SubsystemΩB Modeling Uncertainties

Therefore, the matrices[M¯B],[D¯B],
and [K¯B] of the reduced mean computational model are
replaced by the random matrices [MB],[DB],
and [KB] defined on a probability space (Θ,𝒯,𝒫) and such that(3.1)∀θ∈Θ,⁢[MB(θ)],[DB(θ)],[KB(θ)]∈MnqB+(ℝ),E{[MB]}=[M¯B],E{[DB]}=[D¯B],E{[KB]}=[K¯B],E{∥[MB]−1∥F2}<+∞,E{∥[DB]−1∥F2}<+∞,E{∥[KB]−1∥F2}<+∞,in which 𝕄n+(ℝ) is the set of all the positive-definite
symmetric (n×n) matrices, where E{⋅} is the mathematical expectation and where ∥⋅∥F is the Frobenius norm such that ∥A∥F2=tr{[A]∗[A]} with [A]∗=[A¯]T, [A¯] is the conjugate of [A] and tr is the trace for matrices. Let [PB] be the random matrix denoting [MB],[DB],
or [KB].
The probability distribution of the random matrix [PB] depends on the dispersion parameter δPB related to the coefficient of variation δPB of the random matrix [PB] by the equation(3.2)(δPB)2=E{∥[PB]−[P¯B]∥F2}∥[P¯B]∥F2=(δPB)2n+1(1+(tr[PB])2tr([PB]2)).The dispersion parameter δPB allows the level of uncertainites of the
random matrix [PB] to be controlled. It can be
found in [3, 5] an algebraic
representation of random
matrix [PB] 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.


## 3.2. Identification of the Dispersion Parameters

As explained in Section3.1, the probability
distributions of the random matrices (and then of the random generators) depend
on the vector δ=(δMB,δDB,δKB) of the dispersion parameters which is
identified using the measurements. The observation of the stochastic
computational model is defined introducing the nqB×nqB random complex dynamic stiffness matrix [AB(ω)] of the linear subsystem ΩB written as(3.3)[AB(ω)]=−ω2[MB]+iω[DB]+[KB].Then the random condensed
dynamical stiffness matrix [ZB(ω)] of the linear subsystem ΩB on the coupling interface is such that [ZB(ω)]=[AccB(ω)]−[AcyB(ω)][AyyB(ω)]−1[AycB(ω)].
Taking into account the properties of the probabilistic model, it can be shown
that, for all ω fixed in ℬ,
the random matrix [ZB(ω)] is invertible almost surely and the random
variable J(δ) defined by(3.4)J(δ)=∫ℬ∥[ZB(ω)]−1∥F2dωexists and has a finite mean
value. This random variable gives a measure over ℬ of the dynamical effects of subsystem ΩB on the subsystem ΩA at the coupling interface. It should be noted
that the random variable J(δ) depends on δ because the probability distributions of the
random matrices [MB],[DB],
and [KB] depend on δ.
Let x↦pJ(x,δ) be the probabilty density function of the
random variable J(δ) with respect to dx.
For any x fixed in [0,+∞[ and for any value of the vector δ belonging to the admissible set 𝒞ad of the dispersion parameters, the value pJ(x,δ) of the probability density function is
estimated by using the above probabilistic model and the Monte Carlo
simulation. The corresponding deterministic experimental value Jexp of J(δ) is calculated using experimental data. The
method used to identify vector δ is the maximum likelihood method (see for
instance [10]) for the
random variable J(δ) for which Jexp is one realization. We then have to solve the
following optimization problem:(3.5)δopt=arg⁡max⁡δ∈Cad(pJ(Jexp;δ)),in which δopt is the identified value of δ.


## 3.3. Random Transient Dynamical Response of the Stochastic Nonlinear Computational Model

Using the
probabilistic model defined in Section3.1, the deterministic (2.11) to (2.15)
give the following stochastic nonlinear computational model:(3.6)[U(t)]=[H][Q(t)],in which, for all fixed t,
the vector-valued random variable Q(t) verifies(3.7)[M]Q¨(t)+[D]Q˙(t)+[K]Q(t)+ℱNL(Q(t))=ℱ(t),and where the random matrices [M],[D],
and [K] are written as(3.8)[M]=[M¯yyA0M¯ycA0MyyBMycBM¯cyAMcyBM¯ccA+MccB],[D]=[D¯yyA0D¯ycA0DyyBDycBD¯cyADcyBD¯ccA+DccB],[K]=[K¯yyA0K¯ycA0KyyBKycBK¯cyAKcyBK¯ccA+KccB].


## 4. Identification of Stochastic Loads

The transient load𝔽(t) defined by 𝔽(t)=(𝔽pA(t),𝔽pB(t),0) corresponding to the displacement vector 𝕌(t)=(𝕌pA(t),𝕌pB(t),𝕌cB(t)) is modeled by a stochastic process {F(t),t∈ℝ}.
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 nonzero random components of the random
vector F(t).
This equation can easily be inversed and yields F(t)=Lift(F(t)).


### 4.1. Construction of the Stochastic LoadF(t)

The stochastic load is modeled by anℝm-valued Gaussian stationary centred
second-order stochastic process {F(t),t∈ℝ} defined on a probability space (Θ′,𝒯′,𝒫′) different from the probability space (Θ,𝒯,𝒫). In addition, it is assumed that the stochastic process is mean square
continuous on ℝ,
physically realizable (causal) and for which its matrix-valued autocorrelation
function τ↦[RF(τ)]=E{F(t+τ)F(t)T} is integrable on ℝ.
This stochastic process is then completely defined by its matrix-valued
spectral density function [SF(ω)]=(2π)−1∫ℝe−iωτ[RF(τ)]dτ which is a continuous and integrable function
on ℝ 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 ℝ,
the matrix [SF(ω)] is with values in the set 𝕄m+(ℂ) 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 [11, 12]:(4.1)∫ℝlog⁡(det⁡[S¯F(ω)])1+ω2>−∞.The numerical simulation of
independent realizations {F(t,θ′),t∈ℝ} for θ′∈Θ′ (trajectories) can easily be generated by
using adapted algorithms (see, e.g., [13, 14]).


### 4.2. 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 (3.6) to (3.8) with (2.15) in which the deterministic
load 𝔽(t) is replaced by the stochastic load F(t).
We then have to construct the ℝn𝕌-valued stationary solution Us(t)=(UsA(t),UsB(t),Usc(t)) (corresponding to U(t)) which is written as(4.2)[Us(t)]=[H][Qs(t)],in which subindex s is relative to the stationary solution and
where the stationary stochastic process {Qs(t),t∈ℝ} satisfies the stochastic
equation(4.3)[M]Q¨s(t)+[D]Q˙s(t)+[K]Qs(t)+ℱNL(Qs(t))=[H]TF(t),in which Q˙s(t) and Q¨s(t) are the mean-square first and second
derivative of Qs(t).
For the identification of [SF],
for all t fixed, we introduce the ℝμ-valued random variable Zs(t)=(Zs,1(t),…,Zs,μ(t)) which represents the observations of the
stochastic computational model made up of components of the vector-valued
random response Us(t).
Thus there exists a projection Proj′ from ℝn𝕌 into ℝμ such that Zs(t)=Proj′(Us(t)).
For all θ in Θ,
the stationary stochastic process {Zs(t,θ),t∈ℝ} is such that Zs(t,θ)=Proj′([H]Qs(t,θ)),
where the stationary stochastic process {Qs(t,θ),t∈ℝ} is such that(4.4)[M(θ)]Q¨s(t,θ)+[D(θ)]Q˙s(t,θ)+[K(θ)]Qs(t,θ)+ℱNL(Qs(t,θ))=[H]TLift(F(t;[SF])),where {F(t;[SF]),t∈ℝ} is a stochastic process defined in Section
4.1. Section 4.3 is devoted to the identification of [SF].
In order to perform this identification, we need to introduce an observation
relative to the stochastic equation and which is useful to construct the cost
function. For all θ∈Θ,
the matrix-valued spectral density function {[SZs(ω,θ)],ω∈ℝ} can be estimated. Generating νθ independent realizations of the random
matrices [M],[D],
and [K],
the matrix-valued spectral density function [SZs] is estimated by the Monte Carlo simulation method.
For all ω∈ℝ,
one has(4.5)[SZs(ω)]=1νθ∑i=1νθ[SZs(ω,θi)].


### 4.3. Identification of the Stochastic Loads

The identification[SF] is performed in introducing a parametric
representation of this function which is rewritten as(4.6)[SF(ω)]=[S(ω,r)],ω∈ℝ⁢,⁢r∈Cr,in which 𝒞r⊂ℝνr is the admissible set of the parameter r with values in ℝν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 into 𝕄m+(ℂ).
Therefore, the identification of the stochastic load {F(t),t∈ℝ)} consists in identifying the ℝνr-valued vector r.
Let {Zsexp(t)=(Zs,1exp(t),…,Zs,μexp(t)),t∈ℝ} be the ℝμ-valued stationary stochastic process which is
measured for the manufactured real system and corresponding to the observation
stochastic process {Zs(t),t∈ℝ}.
The matrix-valued spectral density function {[SZsexp(ω)],ω∈ℝ} of this stochastic process is estimated using
the periodogram method. Then, the parameter r is estimated in minimizing the distance D(r)=∫ℬ∥[SZs(ω,r)]−[SZsexp(ω)]∥F2dω 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(4.7)ropt=arg⁢⁡min⁡r∈Cr⁢D(r),in which ropt is the identified value of the vector r.


## 4.1. Construction of the Stochastic LoadF(t)

The stochastic load is modeled by anℝm-valued Gaussian stationary centred
second-order stochastic process {F(t),t∈ℝ} defined on a probability space (Θ′,𝒯′,𝒫′) different from the probability space (Θ,𝒯,𝒫). In addition, it is assumed that the stochastic process is mean square
continuous on ℝ,
physically realizable (causal) and for which its matrix-valued autocorrelation
function τ↦[RF(τ)]=E{F(t+τ)F(t)T} is integrable on ℝ.
This stochastic process is then completely defined by its matrix-valued
spectral density function [SF(ω)]=(2π)−1∫ℝe−iωτ[RF(τ)]dτ which is a continuous and integrable function
on ℝ 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 ℝ,
the matrix [SF(ω)] is with values in the set 𝕄m+(ℂ) 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 [11, 12]:(4.1)∫ℝlog⁡(det⁡[S¯F(ω)])1+ω2>−∞.The numerical simulation of
independent realizations {F(t,θ′),t∈ℝ} for θ′∈Θ′ (trajectories) can easily be generated by
using adapted algorithms (see, e.g., [13, 14]).


## 4.2. 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 (3.6) to (3.8) with (2.15) in which the deterministic
load 𝔽(t) is replaced by the stochastic load F(t).
We then have to construct the ℝn𝕌-valued stationary solution Us(t)=(UsA(t),UsB(t),Usc(t)) (corresponding to U(t)) which is written as(4.2)[Us(t)]=[H][Qs(t)],in which subindex s is relative to the stationary solution and
where the stationary stochastic process {Qs(t),t∈ℝ} satisfies the stochastic
equation(4.3)[M]Q¨s(t)+[D]Q˙s(t)+[K]Qs(t)+ℱNL(Qs(t))=[H]TF(t),in which Q˙s(t) and Q¨s(t) are the mean-square first and second
derivative of Qs(t).
For the identification of [SF],
for all t fixed, we introduce the ℝμ-valued random variable Zs(t)=(Zs,1(t),…,Zs,μ(t)) which represents the observations of the
stochastic computational model made up of components of the vector-valued
random response Us(t).
Thus there exists a projection Proj′ from ℝn𝕌 into ℝμ such that Zs(t)=Proj′(Us(t)).
For all θ in Θ,
the stationary stochastic process {Zs(t,θ),t∈ℝ} is such that Zs(t,θ)=Proj′([H]Qs(t,θ)),
where the stationary stochastic process {Qs(t,θ),t∈ℝ} is such that(4.4)[M(θ)]Q¨s(t,θ)+[D(θ)]Q˙s(t,θ)+[K(θ)]Qs(t,θ)+ℱNL(Qs(t,θ))=[H]TLift(F(t;[SF])),where {F(t;[SF]),t∈ℝ} is a stochastic process defined in Section
4.1. Section 4.3 is devoted to the identification of [SF].
In order to perform this identification, we need to introduce an observation
relative to the stochastic equation and which is useful to construct the cost
function. For all θ∈Θ,
the matrix-valued spectral density function {[SZs(ω,θ)],ω∈ℝ} can be estimated. Generating νθ independent realizations of the random
matrices [M],[D],
and [K],
the matrix-valued spectral density function [SZs] is estimated by the Monte Carlo simulation method.
For all ω∈ℝ,
one has(4.5)[SZs(ω)]=1νθ∑i=1νθ[SZs(ω,θi)].


## 4.3. Identification of the Stochastic Loads

The identification[SF] is performed in introducing a parametric
representation of this function which is rewritten as(4.6)[SF(ω)]=[S(ω,r)],ω∈ℝ⁢,⁢r∈Cr,in which 𝒞r⊂ℝνr is the admissible set of the parameter r with values in ℝν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 into 𝕄m+(ℂ).
Therefore, the identification of the stochastic load {F(t),t∈ℝ)} consists in identifying the ℝνr-valued vector r.
Let {Zsexp(t)=(Zs,1exp(t),…,Zs,μexp(t)),t∈ℝ} be the ℝμ-valued stationary stochastic process which is
measured for the manufactured real system and corresponding to the observation
stochastic process {Zs(t),t∈ℝ}.
The matrix-valued spectral density function {[SZsexp(ω)],ω∈ℝ} of this stochastic process is estimated using
the periodogram method. Then, the parameter r is estimated in minimizing the distance D(r)=∫ℬ∥[SZs(ω,r)]−[SZsexp(ω)]∥F2dω 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(4.7)ropt=arg⁢⁡min⁡r∈Cr⁢D(r),in which ropt is the identified value of the vector r.


## 5. Application

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


### 5.1. Data for the Experimental Model

The
measurements are generated by an experimental model
which is made up of one linear subsystem and one nonlinear subsystem. The
linear subsystem is made up of four parallel beams fixed at their ends. The
nonlinear subsystem 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 Figure1). The five beams are linked by three transversal grids, each
grid being modeled by four transversal springs (see
Figure 1). Therefore, the coupling interface between the two subsystems is
composed of three points located in the neutral fiber of the beam of the
nonlinear subsystem. Each beam is modeled by eight Euler beam finite elements
of equal lengths for which the DOF of the two nodes at the ends of the beam are
locked. The twelve springs defining the three transversal grids are modeled 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
nonlinear subsystem (see Figure 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 nonlinear subsystem is then 14. The beam of
the nonlinear subsystem is exited by seven transversal forces applied following
the y-direction. The vectors of these seven nonzero
components are denoted by fexp.
The stochastic process {fexp(t),t∈ℝ} is a second-order centred stationary Gaussian
stochastic process for which its matrix-valued spectral density function [Sfexp(ω)] is such that (2.1) for all i in {1,…,7}, [Sfexp(ω)]ii is a constant on the frequency band of
analysis ℬ=2π×[−100,100]rad/s,
and (2.2) for all i and j in {1,…,7}, |[Sfexp(ω)]ij|2=γij(ω)[Sfexp(ω)]ii[Sfexp(ω)]jj where γij(ω)=exp(−|xi−xj|/λ) in which |xi−xj| is the distance between the two excited points
and the value of λ is equal to the quarter of the beam length. In
the frequency band of analysis ℬ,
there are 21 eigenfrequencies for the linearized coupled system made up of the
linear subsystem coupled with the linear beam of the nonlinear subsystem
(nonlinear subsystem without the stops), for which the first three
eigenfrequencies are 5.78 Hz, 15.9 Hz, and 31.1 Hz.(a) Experimental model: 3D view. (b)
Transversal view. (c) Tranversal view in the plane of one grid: the 6 diagonal
lines represent the 12 springs.
(a)(b)(c)


### 5.2. Data for the Mean Computational Model

This part is
devoted to the construction of a simplified mean computational model for the
nonlinear dynamical system described in Section5.1. This simplified mean
computational model will be used to identify the stochastic loads. It consists
in modeling (2.1) the four beams and the three transversal grids of the linear
subsystem of the experimental model by a unique equivalent linear Euler beam
and by three equivalent springs (see Figure 2) and (2.2) the linear beam with
elastic stops of the nonlinear subsystem of the experimental model by a linear
beam with two springs for the elastic stops. The 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.Figure 2
Mean model.


### 5.3. 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 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 Pobs 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 Pobs (see Figure 3) is estimated using the
periodogram method. It can be seen that the prevision given by the mean
simplified computational model is good enough in the frequency band [0,50]Hz.
Nevertheless, there are significant differences in the frequency band [50,100]Hz induced by model uncertainties. This is the
reason why the model uncertainties are taken into account in order to extend
the domain of validity of the simplified mean computational model in the
frequency band [50,100]Hz in order to perform a robust identification of
the stochastic loads.For pointPobs,
power spectral density function (PSD) for (a) the y-displacement and (b) the z-rotation. Comparison between the experimental
model (thin line) and the simplified mean computational model (thick
line).
(a)(b)


### 5.4. System Uncertainties Modeling and Dispersion Parameter Identification

The nonparametric probabilistic approach of model
uncertainties introduced in Section3.1 is used for stiffness part of the
linear subsystem of the simplified mean computational model. We then have to
identify the dispersion parameter δ=(δKB).
Note that the identification procedure which is proposed is independent of the
stochastic loads. The estimation of the probability density function in (3.5) is carried out with 200 realizations for the Monte Carlo simulation.
Figure 4 shows the likelihood function calculated using (3.5) with 𝒞ad=[0,22/34].
The maximum is reached for δopt=0.45.Figure 4
Graph of functionδ↦pJ(Jexp;δ).


### 5.5. Case of an Unknown Stochastic Load and Its Identification

In this
section, the responses of the experimental model are given (those constructed
in Section5.3) and the stochastic load F(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 load F(t) introduced in Section 4.1. We have then chosen
to model F(t) as {F(t)=(T(t),M(t)),t∈ℝ} in which T(t) is a y-force and M(t) is a z-moment applied to the middle of the beam of
the nonlinear subsystem (see Figure 5). This force and
this moment are independent second-order centred stationary Gaussian stochastic
processes. So, they are both completely defined by their power spectral density
functions [ST(ω)] and [SM(ω)].
The matrix-valued spectral density function of the stochastic process {F(t),t∈ℝ} is then defined by(5.1)[SF(ω)]=[ST(ω)00SM(ω)],ω∈ℝ⁢.It is assumed that the function ω↦[SF(ω)] is constant in the frequency band of analysis ℬ and is such that (4.1) is verified. The
experimental stochastic process {Zsexp(t),t∈ℝ} defined in Section 4.3 is composed of μ=7 stochastic y-displacements. Taking into account (4.6),
the function ω↦[SF(ω)], which is a constant diagonal hermitian matrix,
can then be rewritten for all ω in ℬ as(5.2)[SF(ω)]=[S(ω,r)]=[r100r2],ω∈ℬ,⁢r∈Cr,in which the admissible set 𝒞r={r=(r1,r2);r1>0,r2>0}.
This vector r is identified using the trial method to solve
the optimization problem defined by (4.7). 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↦log10(D(r)) which allows the optimal value ropt to be determined. The confidence
region associated with a probability level Pc=.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 Figure 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 Pobs.Figure 5
Definition of the stochastic load.Figure 6
Graph of the cost function(r1,r2)↦log10(D(r1,r2)).For pointPobs,
power spectral density function (PSD) for (a) the stochastic y-displacement and (b) the stochastic z-rotation: upper and lower envelopes and mean
response (mid thin line); experimental model (thick line).
(a)(b)


## 5.1. Data for the Experimental Model

The
measurements are generated by an experimental model
which is made up of one linear subsystem and one nonlinear subsystem. The
linear subsystem is made up of four parallel beams fixed at their ends. The
nonlinear subsystem 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 Figure1). The five beams are linked by three transversal grids, each
grid being modeled by four transversal springs (see
Figure 1). Therefore, the coupling interface between the two subsystems is
composed of three points located in the neutral fiber of the beam of the
nonlinear subsystem. Each beam is modeled by eight Euler beam finite elements
of equal lengths for which the DOF of the two nodes at the ends of the beam are
locked. The twelve springs defining the three transversal grids are modeled 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
nonlinear subsystem (see Figure 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 nonlinear subsystem is then 14. The beam of
the nonlinear subsystem is exited by seven transversal forces applied following
the y-direction. The vectors of these seven nonzero
components are denoted by fexp.
The stochastic process {fexp(t),t∈ℝ} is a second-order centred stationary Gaussian
stochastic process for which its matrix-valued spectral density function [Sfexp(ω)] is such that (2.1) for all i in {1,…,7}, [Sfexp(ω)]ii is a constant on the frequency band of
analysis ℬ=2π×[−100,100]rad/s,
and (2.2) for all i and j in {1,…,7}, |[Sfexp(ω)]ij|2=γij(ω)[Sfexp(ω)]ii[Sfexp(ω)]jj where γij(ω)=exp(−|xi−xj|/λ) in which |xi−xj| is the distance between the two excited points
and the value of λ is equal to the quarter of the beam length. In
the frequency band of analysis ℬ,
there are 21 eigenfrequencies for the linearized coupled system made up of the
linear subsystem coupled with the linear beam of the nonlinear subsystem
(nonlinear subsystem without the stops), for which the first three
eigenfrequencies are 5.78 Hz, 15.9 Hz, and 31.1 Hz.(a) Experimental model: 3D view. (b)
Transversal view. (c) Tranversal view in the plane of one grid: the 6 diagonal
lines represent the 12 springs.
(a)(b)(c)


## 5.2. Data for the Mean Computational Model

This part is
devoted to the construction of a simplified mean computational model for the
nonlinear dynamical system described in Section5.1. This simplified mean
computational model will be used to identify the stochastic loads. It consists
in modeling (2.1) the four beams and the three transversal grids of the linear
subsystem of the experimental model by a unique equivalent linear Euler beam
and by three equivalent springs (see Figure 2) and (2.2) the linear beam with
elastic stops of the nonlinear subsystem of the experimental model by a linear
beam with two springs for the elastic stops. The 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.Figure 2
Mean model.


## 5.3. 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 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 Pobs 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 Pobs (see Figure 3) is estimated using the
periodogram method. It can be seen that the prevision given by the mean
simplified computational model is good enough in the frequency band [0,50]Hz.
Nevertheless, there are significant differences in the frequency band [50,100]Hz induced by model uncertainties. This is the
reason why the model uncertainties are taken into account in order to extend
the domain of validity of the simplified mean computational model in the
frequency band [50,100]Hz in order to perform a robust identification of
the stochastic loads.For pointPobs,
power spectral density function (PSD) for (a) the y-displacement and (b) the z-rotation. Comparison between the experimental
model (thin line) and the simplified mean computational model (thick
line).
(a)(b)


## 5.4. System Uncertainties Modeling and Dispersion Parameter Identification

The nonparametric probabilistic approach of model
uncertainties introduced in Section3.1 is used for stiffness part of the
linear subsystem of the simplified mean computational model. We then have to
identify the dispersion parameter δ=(δKB).
Note that the identification procedure which is proposed is independent of the
stochastic loads. The estimation of the probability density function in (3.5) is carried out with 200 realizations for the Monte Carlo simulation.
Figure 4 shows the likelihood function calculated using (3.5) with 𝒞ad=[0,22/34].
The maximum is reached for δopt=0.45.Figure 4
Graph of functionδ↦pJ(Jexp;δ).


## 5.5. Case of an Unknown Stochastic Load and Its Identification

In this
section, the responses of the experimental model are given (those constructed
in Section5.3) and the stochastic load F(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 load F(t) introduced in Section 4.1. We have then chosen
to model F(t) as {F(t)=(T(t),M(t)),t∈ℝ} in which T(t) is a y-force and M(t) is a z-moment applied to the middle of the beam of
the nonlinear subsystem (see Figure 5). This force and
this moment are independent second-order centred stationary Gaussian stochastic
processes. So, they are both completely defined by their power spectral density
functions [ST(ω)] and [SM(ω)].
The matrix-valued spectral density function of the stochastic process {F(t),t∈ℝ} is then defined by(5.1)[SF(ω)]=[ST(ω)00SM(ω)],ω∈ℝ⁢.It is assumed that the function ω↦[SF(ω)] is constant in the frequency band of analysis ℬ and is such that (4.1) is verified. The
experimental stochastic process {Zsexp(t),t∈ℝ} defined in Section 4.3 is composed of μ=7 stochastic y-displacements. Taking into account (4.6),
the function ω↦[SF(ω)], which is a constant diagonal hermitian matrix,
can then be rewritten for all ω in ℬ as(5.2)[SF(ω)]=[S(ω,r)]=[r100r2],ω∈ℬ,⁢r∈Cr,in which the admissible set 𝒞r={r=(r1,r2);r1>0,r2>0}.
This vector r is identified using the trial method to solve
the optimization problem defined by (4.7). 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↦log10(D(r)) which allows the optimal value ropt to be determined. The confidence
region associated with a probability level Pc=.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 Figure 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 Pobs.Figure 5
Definition of the stochastic load.Figure 6
Graph of the cost function(r1,r2)↦log10(D(r1,r2)).For pointPobs,
power spectral density function (PSD) for (a) the stochastic y-displacement and (b) the stochastic z-rotation: upper and lower envelopes and mean
response (mid thin line); experimental model (thick line).
(a)(b)


## 6. Conclusions

We have presented a methodology and its validation to
perform the identification of a stochastic loads applied to a complex nonlinear
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 nonlinear dynamical system used for this validation is
representative of real industrial systems and then validates the methodology
proposed.



---

*Source: 181548-2008-07-15.xml*