Dynamics of random coupled structures through the wave finite element method

Purpose – The purpose of this paper is to develop a new formulation using spectral approach, which can predict the wave behavior to uncertain parameters in mid and high frequencies. Design/methodology/approach – The work presented is based on a hybridization of a spectral method called the “ wave finite element (WFE) ” method and a non-intrusive probabilistic approach called the “ polynomial chaos expansion (PCE). ” The WFE formulation for coupled structures is detailed in this paper. The direct connection with the conventional finite element method allows to identify the diffusion relation for a straight waveguide containing a mechanical or geometric discontinuity. Knowing that the uncertainties play a fundamental role in mid and high frequencies, the PCE is applied to identify uncertainty propagation in periodic structures with periodic uncertain parameters. The approach proposed allows the evaluation of the dispersion of kinematic and energetic parameters. Findings – The authors have found that even though this approach was originally designed to deal with uncertainty propagation in structures it can be competitive with its low time consumption. The Latin Hypercube Sampling (LHS) is also employed to minimize CPU time. Originality/value – The approach proposed is quite new and very simple to apply to any periodic structures containing variabilities in its mechanical parameters. The Stochastic Wave Finite Element can predict the dynamic behavior from wave sensitivity of any uncertain media. The approach presented is validated for two different cases: coupled waveguides with and without section modes. The presented results are verified vs Monte Carlo simulations.


Introduction
Wave propagation in structures is one of the important subject treated in this last decade. Wave characteristics (wavenumber, velocities, deformed shapes) can be identified by different methods such as spectral finite element (SFE) and the wave finite element method (WFE).
The wave finite element is a spectral method developed for phononic and photonic media and for wave propagation in crystal lattices (Brillouin, 1946). Zhong and Williams (1995) extended this approach for wave propagation in periodic structures which is composed by identical sub-elements. This formulation, with direct connection with the conventional finite element method, allows to predict wave characteristics (propagation constant, deformed shapes, etc.). This spectral study starts from the discretization of one sub-element with the conventional FEM, then a spectral eigenvalue problem is formulated regarding the periodicity of the media. The general formulation is proposed by Mead (1973). The WFE is then extended for different problems: Mace et al. (2005) and Mencik and Ichchou (2005) used the WFE to identify the characteristics of free structural vibration. In addition, forced dynamic response with harmonic load (Duhamel et al., 2006) and general load (Renno and Mace, 2010) was investigated.
The WFE formulation is based on an principal assumption: the periodicity of the structure. In reality, structures can contain some irregularities or section changes; in general manner, every discontinuity of mechanical impedance can be modeled as a mechanical discontinuity. In this case, the WFE must be modified to predict the behavior of waves regarding these discontinuities. Many researchers studied wave scattering and tried to express the expression of the diffusion coefficients using analytical procedure (Mace, 1984;Mei et al., 2006). Mencik and Ichchou (2005) developed a Diffusion Matrix Model (DMM) which is based on the hybridization of the conventional finite element method and the WFE. This formulation allows us to obtain the reflection and transmission coefficients of waves encountering mechanical discontinuities. This formulation is then used and validated by different researchers Zhou et al., 2009). Bareille et al. (2012) applied the DMM to identify the size and location of defects on metallic pipelines. The DMM approach allows building a numerical cartography regarding the different sizes of the defect. This cartography will be used regarding experimental reflection coefficients.
Knowing that uncertainties play an important role in high frequencies, the WFE which uses deterministic inputs, will be penalized for structures with uncertain parameters. For this purpose, the exploration of probabilistic and non-probabilistic approaches is needed. These approaches take into account the variability of physical parameters which can modify the mechanical behavior of structures in high-frequency range. In order to identify the level of variability in non-deterministic systems (Blanchard et al., 2009), different mathematical approaches are used such as Fuzzy theory (Moens and Hanss, 2011;Giannini and Hanss, 2008), interval arithmetic (Dessombz et al., 2001), Bayesian theory (Goller and Schueller, 2011), parametric (Ghanem and Spanos, 1991) and non-parametric approaches (Soize, 2000(Soize, , 2005. The aim of these stochastic methods is to reduce the gap among the simulated and the experimental results. The polynomial chaos expansion (PCE) is an efficient probabilistic tool to quantify the uncertainty propagation in structures. This theory is based on the expansion of a second order process into orthogonal polynomes. This approach allows the separation of the deterministic and the stochastic component of the original process. Two different approaches are often used to evaluate the stochastic modes; intrusive and non-intrusive approaches. The intrusive one is based on Galerkin projection which allows the expansion of the stochastic equation into deterministic ones, so the formulation is modified for each studied problem. The non-intrusive approach is more efficient because it is based on some iterations of the deterministic code. Finally, an adequate post-traitement is then used to quantify the uncertainty propagation in the structure studied.
The parametric probabilistic approach consists in developing the stochastic random variables on a series by decoupling the mean and the fluctuation parts. The WFE is already extended to uncertain structures using a parametric probabilistic approach. Ichchou et al. (2011) uses the first-order perturbation to express variability as an other dimension of the spectral problem. The presented work offered an explicit expressions of statistics of wave characteristics (propagation constant, mode shapes, group velocity). Ben Souf et al. (2013a) used the WFE and the Stochastic Finite Element Method (SFEM) to express the variabilities' effects on the coupling loss factors for two periodic waveguides connected through a junction with mechanical and geometrical variability. This study begins by the evaluation of the dispersion of the kinematic and energetic diffusion coefficients. This study is then extended to predict the forced response of stochastic coupled structures (Ben Souf et al., 2013b).
In previous works, the first-order perturbation is used to quantify the uncertainty effects on the wave propagation in periodic structures. The approach developed proves its efficiency by expressing kinematic and energetic parameters for uncertain and periodic structures. Although, it was limited for low dispersion and for few number of uncertain parameters. The main novelty in this paper consist in the hybridization of the WFE (spectral and forced response formulations) and the Generalized Polynomial Chaos Expansion (GPCE). The main assumption in this work, that the uncertainties are periodic as well as the structure (The mechanical parameters variabilities are periodic because the deterministic spectral formulation is based on the periodicity of media). The presented formulation tries to identify the wave sensitivity on the uncertain inputs. The use of a non-intrusive approach allows to identify the stochastic modes with no changes in the deterministic code. This paper is composed by three different parts: the first one is a general presentation of the deterministic WFE for periodic media. Section 2 deals with the forced response of coupled structures and wave behavior regarding to mechanical discontinuities. Section 3 presents the Generalized PCE. Finally, different studied cases are presented in the last section and validated vs Monte Carlo simulations.

The WFE
The spectral WFE is a numerical method which helps to identify wave characteristics in periodic structures. It is applicable for media composed by identical sub-elements connected along a specific direction. The studied structure is supposed to be linear, dissipative and elastic ( Figure 1).
This study begins with the discretisation of one sub-element of length d of the structure to express the dynamic stiffness matrix is using the extracted mass M and stiffness K matrices from the FEM model and condensed into the left (denoted by L) Figure 1. An illustration of a periodic waveguide 3 and right (denoted by R) sides as: where q and F R represent the displacement and the forces vectors, respectively. Let's define two state vectors in the left and the right side of each element as: These vectors are related by a 2n × 2n symplectic matrix S as: The dynamic of the global system is formulated from the description of the waves propagating along the xaxis. Using the periodicity and the Bloch's theorem, a spectral eigenvalue problem is formulated as (Mencik and Ichchou, 2005): To face out numerical problems, this eigenvalue problem can be solved using the approach proposed by Zhong and Williams (1995). The resolution of Equation (4) leads to identify the propagation constant and the mode shapes. Since the structure is dissipative, the waves can be classified as an incident 9m inc i 9o 1 ÀÁ and reflected ones 9m ref i 94 1 ÀÁ . Consequently, the wave basis Φ composed by the eigenvectors, can be expressed as follows: where F inc q , F inc F , F ref q and F ref F are n × n matrices and q and F refers to the displacement and the force component of the eigenvectors.
Finally, the dynamics of the global system is formulated from the projection of the kinematic variables for the kth element on the wave basis as: represents the amplitude of waves in the kth element.

Forced response for coupled structures through a spectral approach
The dynamic behavior of two finite structures connected through a coupling element is presented in this section. The formulation presented is a solution for loss of periodicity of the structure. The first step consists on expressing wave behavior near junction, the formulation for the forced response of the whole media is then offered.

Wave behavior regarding discontinuities
This part deals with prediction of wave behavior regarding mechanical and geometric singularities. In general manner any mechanical impedance in structures, causes a reflection and transmission waves from generated ones. The presented formulation is based on the hybridization of the conventional finite element method and the Wave finite element. The formulation developed called the Diffusion Matrix Model (DMM) (Mencik and Ichchou, 2005;Ichchou et al., 2009;Renno and Mace, 2012) which model the junction as coupling stiffness and helps to quantify the reflection and transmission coefficients. Bareille et al. (2012) and Zhou et al. (2009) used the DMM formulation to identify the location and the size of discontinuities in pipelines regarding experimental results. In this paper, the deterministic formulation is shortly presented.
Let us define two semi-infinite structures connected through a coupling element. This junction could be a crack, an elbow, etc. In a general manner, it is considered as discontinuity of mechanical impedance caused by mechanical or geometrical changes ( Figure 2).
The DMM formulation begins by the expression of the dynamical equation, condensed into the right and the left sides, of the coupling element and the sub-elements directly connected to the junction which can be written as: where matrix K stands for the dynamical stiffness of the coupling element q c 1 ; F c 1 ÀÁ and ðq c 2 ; F c 2 Þ represent the displacements and the forces applied at the dof 's of the coupling element on surfaces Γ 1 and Γ 2 , respectively. The coupling conditions can then be expressed on the interface in the matrix form as: Assuming that the that the left and the right cross-sections of the given subsystem contains the same number of degrees of freedom, Equations (8) and (7) lead to the relation: Guide 1 Guide 2 Coupling element d 1 d c d 2 Figure 2. Two E-B beam connected through a stochastic coupling element In addition, the projection of the state vectors into the wave basis Φ (1) and Φ (2) , respectively, for the waveguide 1 and 2, allows to appear the amplitude of waves as: Indeed, using the Equations (9) and (8) we can write a relationship between the amplitude of waves incident to and reflected by the junction as follows: Finally, the diffusion relation can be expressed: where: where + is the pseudo-inverse.

Forced response and energy densities in periodic structures
This section deals with the forced response of finite coupled structures to harmonic force excitation using a spectral approach. The approach used is based on wave approach. The dynamics of the global system is described through o projection of kinematic variables on wave basis constructed using the WFE.
Let us consider two waveguides with different lengths L 1 and L 2 and composed, respectively, by N 1 and N 2 sub-elements with two different origins of coordinates.
The first step consist on performing the WFE for two sub-elements for each waveguide to identify different waves. Then the DMM is used to express the relation between the reflected and incident waves in waveguide 1 and waveguide 2) as: where I and II refer to the first and second waveguide, respectively. In order to evaluate kinematic variables in each node in the periodic media, this step consist in identifying the amplitude of waves in the first sub-element for two waveguides by applying the general wave-mode expansion for boundary conditions (forced-clamped coupled structures Figure 3).
Hence, for the first waveguide, the boundary and transfer conditions on x 1 ¼ 0 and on x 1 ¼ L 1 , respectively, can be expressed using wave mode amplitude as: Equivalently for the second waveguide, the boundary and transfer conditions on x 2 ¼ 0 and on x 2 ¼ L 2 , respectively, can be expressed using wave mode amplitude as: C o u p li n g e le m e n t X2 Z2 Y2 Figure 3. Two coupled periodic waveguides According to the Bloch's theorem, the stochastic wave amplitude on kth element can be obtained from the amplitude of waves on the first sub-element as: Finally, the boundary conditions can be express in matrix form as: The direct inverse of the presented matrix is impossible regarding its diagonal components, so an appropriate scaling is needed by decoupling the ill-conditioned matrix into two matrices (well conditioned and diagonal one). This equation leads to identify the wave amplitude of waves in x 1 ¼ 0 and x 2 ¼ 0.
Using Equation (20), the amplitude of waves can be calculated for each node and the displacements can be easily expressed for kth element as: Assuming that the length of the sub-elements is sufficiently small, different energy densities can be approximated for discrete structures: for the kinematic energy density, and the strain energy density can be expressed as: where q L and q R represent the displacements in the left and the right sides on the sub-element k.

GPCE
The PCE is an efficient tool to describe uncertainty propagation in mechanical systems. The theory, developed by Wiener (1938), helps to expand any second order process u (with finite variance) in a series of orthogonal polynomes as: where H p ðx i 1 ; ...; x i p Þ represents orthogonal polynomes (chaos polynomes) of order p.
In a compact form, Equation (25) can be expressed as: where x ¼½x i 1 ; ÁÁÁ; x i p T , and M denoting the number of the uncertain parameters.
The PCE is used to develop the Stochastic Spectral Finite Element Method (SSFEM) by Ghanem and Spanos (1991), used for instance, in fluid analysis (Najm, 2009), etc. Since in most applications the stochastic input variables are not normal, Xiu and Karniadakis (2002) proposed a generalized form of Hermite PCE using other orthogonal polynomes in terms of non-Gaussian random variables called Wiener-Askey. Table I resumes usual random variables and their orthogonal polynomials.
When the input parameters have not a non-Gaussian behavior, the parameterization of the problem is quite difficult. Rosenblatt (1952) proposed simple transformations of non-Gaussian distributions to Gaussian ones.
Let us define Y ¼ [Y 1 , …, Y n ] T and U ¼ [U 1 ,…, U n ] T a random vector with dependent variables and normal vector, respectively. The probabilistic transformation U ¼ Tr(Y) can be obtained using:

Intrusive/non-intrusive approach
Expressing the output variables using the PCE requires the estimation of the coefficients to quantify the uncertainty propagation on the studied structure. Two different methods are often used: intrusive and non-intrusive methods.
The intrusive approach is based on Galerkin projection of the stochastic equation using the the orthogonality of the polynomial chaos basis. The original problem is reduced to a linear (N × P) system where N is the number of dofs and P represents the number of terms in the PCE. This method needs a specific implementation for each studied problem. For complex systems, this method requires an important computational time and memory. In order to face out these problems, non-intrusive schemes will be used.
4.1.1 Non-intrusive methods. The non-intrusive methods allow us the evaluation of the polynomial chaos coefficients using the deterministic code without any specific modifications. In addition, a pre-processing affecting the random input variables is necessary. Two different approaches are proposed: the projection and the regression approaches.
To illustrate these approaches, let us assume that a second order random variable is presented onto a generalized PCE as follows: The projection approach consists in using of the orthogonality of the Polynomial Chaos basis. By premultiplying Equation (28) par C i fxðyÞg M k¼1 ÀÁ and taking the expectation, the coefficients s j can be obtained as: To evaluate s j , the numerator must be calculated in numerical way (denominator is well known). The numerator can be expressed in integral form as: Distribution Transformation Uniform (a, b) a þðbÀaÞ 0:5 þ0:5erf ðx= ffiffi ffi 2 p Þ ÀÁ Normal (μ, σ) μ + σξ Lognormal (μ, σ) exp(μ + σξ) Gamma (a, b) ab x ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 9a þ 1 À 1 9a ÀÁ q 3 Exponential (λ) À 1 l log 1 2 þ 1 2 erf x= ffiffi ffi 2 p ÀÁ ÀÁ Table II.

Random variables and transformations
In order to evaluate this integral, different methods are used such as Monte Carlo simulations or multidimensional quadrature schemes.
The regression approach is based on the minimization, in the sense of least square measure, the difference ΔS between the infinite development of stochastic parameter and the truncated one regarding the coefficients {s i , i ¼ 1, …, P − 1}, and: So, the minimization is obtained as: In explicit form, Equation (32) can be written as: Finally a system of equations can be solved: (34) In matrix form, the Equation (34) can be expressed as: This formulation, developed by Berveiller (2005), allows the identification of the coefficients by evaluating the deterministic term S (i) ¼ S(X (i) ) where X (i) represents a realization of random input parameters. The direct inversion can cause some numerical problems regarding the ill-conditioned matrix. Isukapalli (1999) proposed to build an experimental design based on the roots of orthogonal polynomes. Berveiller (2005) proved that n ≃ (M−1)P is an optimal size of the experimental design to obtain accurate coefficients. Finally, an appropriate post-treatment leads to identify statistical moments as (Berveiller, 2005): d ijk s i s j s k

Validation cases
This section deals with numerical validations of the formulation proposed to identify uncertainty propagation through a spectral method. Two different cases are tested: t h ef i r s to n et r e a tt h ec a s eo ft w oi d e n t i c a l waveguides with classical propagating waves. The second case studied is based on two waveguides with hollow rectangular sections which are the seat of a lot of propagating waves (classical and sectional ones). For two cases, the uncertainties is supposed to be periodic as well as the studied structure. The structure studied is supposed to be excited in x 1 ¼ 0 m and clamped in x 2 ¼ 0m. The first step consist on evaluating the propagating waves on the structures studied. The direct application of the deterministic WFE proves that waveguides are the seat of only classical waves in the whole studied frequency range.

Two connected waveguides with rectangular cross section
Assuming that the uncertainties are homogenous in the media, the second step consist in evaluating the effect of uncertain parameters on the spectral behavior regarding the coupled waveguides. The Table III summarizes the used random parameters and their distributions. This concerns five uncertain parameters (excitation force and mechanical properties).
Let's define n realizations of ffx ðiÞ k g M k¼1 ; i¼ 1; ÁÁÁ; ng. Using the iso-probabilistic transformations, we can obtain the n realizations of the input random variables {X (1) , …, X (n) }. For every realization, the deterministic WFE is used to evaluate the output variables to identify the coefficients of the PCE series. The regression approach is used to quantify the uncertainty effects on kinematic and energetic parameters. In order to reduce time consuming, the Latin Hypercube Sampling (LHS) is used which represents a type of stratified Monte Carlo sampling.
The structure studied contains a mechanical discontinuity which is modeled as a junction between two waveguides. When waves propagate, the impedance discontinuity generates reflected and transmitted waves. Figure 5 represents the reflection and the transmission coefficients for propagating waves in waveguides (1: flexural 1, 2: flexural 2, 3: torsional, 4: longitudinal). The mean value, calculated by WFE (--) and MC (*), represents the evolution of propagating waves regarding the coupling element. The yellow areas quantify the level of dispersion of wavenumbers regarding random inputs.
From the knowledge of the wave behavior through a junction, the global behavior is then observed. The Figures 6 and 7 represent the mean of the displacement evaluated in x 1 ¼ 0.25 m and x 2 ¼ 0.25 m. Figure 6 represents the mean and the envelop (min-max) constructed using the PCE. The use of the Monte Carlo simulations presents a good agreement with the proposed approach using the fourth order expansion. The colored envelop is larger in mid and high frequencies which proves the important effect of uncertainties of mechanical behavior in this frequency range. Figure 7 represents the standard deviation of displacement evaluated in x 1 ¼ 0.25 m and x 2 ¼ 0.25 m. It's represents the dispersion of kinematic output variables due to uncertain input parameters. The convergence of the approach proposed is studied here. The third order PCE is first tested (dashed line). We remark that even the frequency increases, the results of the third order expansion does not converge to the reference results, so, the use of the fourth order expansion is more accurate to predict the variance of kinematic variables especially in high-frequency range.  Table III. Random variables 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000      As presented in the theoretical part, it is interesting to evaluate the variance of energetic parameters in the periodic structure and their dispersion regarding uncertain inputs. Assuming that the energy densities weakly varies in one layer of the structure (length of sub-element is small regard to wavelength), so the stochastic energy densities can easily expressed using discrete kinematic variables. Figure 8 represents the mean and the standard deviation of the kinematic energy density, respectively, in x 1 ¼ 0.25 (m). For these two cases, the fourth order expansion can predict the energetic behavior of the coupled structures. These results are validated by Monte Carlo simulations using 3,000 samples.

Two coupled structures with hollow section
This second validation case deals with a waveguide with hollow section which contains a mechanical discontinuity. This structure will be modeled as two periodic waveguides connected through a junction. The mean difference vs the previous case, this type of media has, even classical modes, a lot of sectional modes. The prediction of the random behavior of this type of modes, which can exchange energy with conventional ones, represents a numerical challenge (Figure 9). 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000 0   For this purpose, let us define two identical hollow waveguides connected through a coupling element. The studied structure is discretized using SHELL63 element. The waveguides mechanical characteristics are: density ρ ¼ 7,800 kg/m 3 , Young's modulus E ¼ 2 × 10 11 Pa, Poisson's ratio ν ¼ 0.3, structural damping η ¼ 1 percent, dimensions of cross section area: 3 · 10 −2 m × 2 · 10 −2 m and the thickness ¼ 1 mm. The first step consists in identifying propagating waves, so the deterministic WFE is applied for one layer. Figure 10 represents the wavenumbers of different propagating waves. Below 2150 Hz, we have only classical waves. Going higher in frequency, the successive triggering of sectional waves modifies the wave behavior. Figure 10 emphasizes the energy exchanged between waves in the propagating phase. For example, waves 6 and 7 (shear modes) interact with torsional wave (10) and modify its behavior beyond 2150 Hz (normally linear). Table IV resumes the different waves and their deformed shapes ( Figure 11). Now, let us consider that the structure studied contains some random variables (Table V). These uncertain parameters affect directly the kinematic and energetic behavior of the structure studied. The same process, as previous case, is used to quantify the variance of output parameters. Figure 12 represents the mean and the standard deviation of the displacement in x 1 ¼ 0.25 m. The yellow envelop which represents the effect of random input parameters, is larger in high frequencies. The standard deviation, calculated by the formulation presented, can predict correctly the stochastic behavior of the structure and coincide with the reference results (Monte Carlo simulations) in all frequency range. Figures 13 and 14 represents the statistics of kinematic and potential energy densities (mean and standard deviation) evaluated in x 1 ¼ 0.25 m, respectively. The presented formulation proves its efficiency (using the fourth order expansion) regarding the Monte Carlo simulations, while, some fluctuations can be identified for the standard deviation. This behavior, in high frequency, can be explained by the number of the random inputs and the level of uncertainties and their effects in this frequency range.

Conclusion
In this paper, the uncertainty propagation in periodic coupled structure is presented and validated. This uncertainty is quantified through a spectral method using a constructed wave basis to describe the dynamic behavior of structures. Hybridization between the WFE and a parametric probabilistic approach was presented. It has been applied for two different cases: connected waveguides without sectional modes and connected structures with section modes. In the two cases, the formulation presented proves a good agreement with reference method.
In our case, the non-intrusive approach was used, so no changes in the deterministic code are required. Just an adequate post-processing allows us to express the variance of  random output variables. The presented formulation presents an alternative stochastic approach with low time and memory consuming. The presented formulation can predict the parametric uncertainties on structures. The PDF of stochastic inputs can be identified using experiments. But, the formulation presented can not take into account the model uncertainties. In this case, an other approach, such as non-parametric approach, can be used. As a future work, this formulation can be coupled with other dynamics problem such as control process, defects detection, vibro-acoustics and energetic methods for high-frequency vibrations.