A Linearised Hybrid FE-SEA Method for Nonlinear Dynamic Systems Excited by Random and Harmonic Loadings

The present paper proposes a linearised hybrid finite element-statistical energy analysis (FE-SEA) formulation for built-up systems with nonlinear joints and excited by random, as well as harmonic, loadings. The new formulation was validated via an ad-hoc developed stochastic benchmark model. The latter was derived through the combination of the Lagrange-Rayleigh-Ritz method (LRRM) and the Monte Carlo simulation (MCS). Within the build-up plate systems, each plate component was modelled by using the classical Kirchhoff’s thin-plate theory. The linearisation processes were carried out according to the loading-type. In the case of random loading, the statistical linearisation (SL) was employed, while, in the case of harmonic loading, the method of harmonic balance (MHB) was used. To demonstrate the effectiveness of the proposed hybrid FE-SEA formulation, three different case studies, made-up of built-up systems with localized cubic nonlinearities, were considered. Both translational and torsional springs, as joint components, were employed. Four different types of loadings were taken into account: harmonic/random point and distributed loadings. The response of the dynamic systems was investigated in terms of ensemble average of the time-averaged energy.


Introduction
Manufacture uncertainties are widespread in various industrial applications, e.g., aerospace, civil, mechanical, and marine engineering. The structural vibrations which arise from these applications, are normally investigated by using a variety of computational methods. One of the most used is the finite element method (FEM). It relies on a very large number of degrees of freedom (DOFs), which makes the determination of the dynamic system response rather complex. In addition to this, the uncertainties of the complex systems largely decrease the prediction accuracy of high-frequency structural vibration, due to the fact that high-frequency modes are very sensitive to uncertainties. A computational technique that can successfully be used to overcome these shortcomings is the statistical energy analysis (SEA). It is a powerful tool in the analysis of dynamic systems, and above all when it comes to predict the energy transfer within complex system for response in high-frequency range. In more than half a century's development, the SEA has demonstrated its advantages when analysing several engineering applications.
The SEA aims to obtain the average energy level response of an ensemble of dynamic systems which are featured by uncertainties. It is based on the energy equilibrium and the assumption that on the energy response of nonlinear dynamic systems subjected to random loading. More specifically, nonlinear dynamic systems with localized cubic nonlinearities introduced by translational and torsional springs, as joint components, are taken into account. Both harmonic and random loading-types are considered, and both point or distributed loadings are applied. The response of the dynamic systems was examined in terms of ensemble average of the time-averaged energy.

Benchmark Model-Lagrange-Rayleigh-Ritz Method
This section is entirely devoted to the derivation of the governing equations (GEs) for the benchmark model. The latter is used in all of the proposed case studies for validation purpose. The GEs are based on Kirchhoff's thin-plate theory, the LRRM is performed to solve the linearised GEs. Various scenarios accounting for inclined plates and nonlinear translational and rotational springs, as well as several loading-types, are considered. With respect to the solution of the nonlinear GEs, the MHB and the SL are employed as linearisation techniques for systems excited by harmonic and random loadings, respectively. Some further information can be found in a previous article [19].

Built-Up System with Inclined Plate
The built-up system, schematically shown in Figure 1, is excited by an harmonic force P orthogonal to the inclined plate. The two plates are considered simply-supported, and the plate on top is inclined to an angle α with respect to the plate at the bottom. In the system, the out-of-plane motion is only considered. The stretch/compression of the translational spring can be written as ∆L s = w 1 cos(α) − w 2 ; where w 1 and w 2 denote the transverse displacement of plate 1 and plate 2, respectively. Then, the elastic potential energy of the built-up plate system, is given as: The kinetic energy of the system, including the randomly distributed masses on the plates, is given as follows: The potential energy related to the application of the external force assumed to be concentrated and perpendicular to the upper plate can be written as It should be noted that, in addition to the random and harmonic point loadings, both rain-on-the-roof and harmonically distributed loading-types are taken into account in the present investigation. By using the Lagrange equations, the GEs for a built-up plate system with nonlinear springs, inclination angle end harmonic point loading can be written as In Equations (5) and (6), q mn are the time-dependent modal coordinates; ψ mn are the mass-normalized shape functions; m and n represent the Ritz expansion order in the x and y direction, respectively; N 1 , m and N 2 , m are the number of the distributed lumped masses used to randomize plate 1 and plate 2, respectively; k 1 and k 3 are the linear and nonlinear stiffness coefficients of the connection springs; and N s refers to the number of springs introduced as joint elements amongst subsystems.

Random Loading and Statistical Linearization
Consider the situation schematically shown in Figure 1 with the random loading acting on plate 1 but its inclination angle is equal to zero (α = 0). The external force is vertically downward and concentrated, assumed to be white noise with mean value µ P and standard deviation σ P . The linearised equations of motion can be written as where M 0 and C 0 are mass and damping matrices; q is the generalized displacement; P represents generalized force; K eq is the equivalent stiffness matrix. To solve the above equations, the calculation of the equivalent stiffness is fundamental. In this respect, the SL can be successfully used.
With respect to the single-degree system, described by the differential equation mẍ + cẋ + k 1 x + k 3 x 3 = F, where m, c, k denote the mass, damping, and stiffness of this system; x is the displacement; and the external force F follows normal distribution F ∼ N(µ, σ 2 ), statistical linearisation can be performed and the equation is transformed as mẍ + cẋ + k eq x = F, where k eq is referred as equivalent stiffness; x 2 denotes the expectation of x 2 [23].
For the built-up system in Figure 1, the stretch or the compression of the translational spring ∆ can be given as thence, It should be born in mind that the linearised stiffness matrix derived by using MHB and SL, and considering a cubic nonlinearity, are different. For the former, it can be written as k eq = k 1 + 3 4 k 3 x 2 , while, for the latter, it is given in Equation (8). In addition, since Equation (7) is in time domain, it is necessary to obtain the energy in a time-averaged form (see Ref. [19]).

The Linearized Hybrid FE-SEA Formulation
The present section provides a overview of the hybrid FE-SEA formulation accounting for nonlinearities in the point joints of dynamic system. The hybrid FE-SEA formulation firstly requires the identification of those components, within the system, which are assumed to behave statistically. These components are modelled as SEA subsystems. The remaining components are deemed to be deterministic and are modelled by FE method. The relationship between the SEA and the FE subsystems is considered to satisfy the following conditions [5]: where q is the general displacement vector of FE parts under the frequency of ω; f represents the external forces vector exerted to the FE components; f k rev is the forces vector resulting from the reverberant field in k-th subsystem; D d corresponds to the dynamic stiffness matrix of the deterministic components; and D k dir is the the dynamic stiffness matrix arising from k-th direct field. Considering the diffuse field reciprocity relation between direct fields and reverberant fields [6], the energy equilibrium equation for each subsystem and the cross spectral matrix S qq is given as [5]: where In Equation (13), η j is the loss factor of j-th subsystem; η d,j corresponds to the power dissipation in j-th master system; η jk is the coupling loss factor; n j is the modal density; E j is the ensemble average energy of j-th subsystem; P in,j and P ext in,j represent the power input from the loadings to subsystems and to master systems, respectively. In Equation (14), S f f denotes the cross spectral matrix of external forces to master systems. Usually, Equations (13) and (14) are used to obtain the response of subsystems and FE components. To solve Equation (13), P ext in,j , η jk and η d,j can be calculated by Equations (15)- (17). Then, the responses of deterministic components are obtained using Equation (14). As far as the localized nonlinearity is concerned, it is treated exactly in the same way of the LRRM benchmark model. In the previous authors' article [19], it is illustrated the linearisation of the MHB corresponding to harmonic loading; this section similarly applies the SL for the localized cubic nonlinearities (translational and/or torsional springs) under random loading.

Numerical Results
This section provides both validation and assessment of the linearised hybrid FE-SEA formulation regarding to both harmonic and random excitations. Three case studies are addressed; the first case study is made up of a three-plate build-up system accounts for both harmonic and random point loadings, rain-on-the-roof loading-type and inclination angle of the driven plate; the second case study focuses on the harmonically distributed excitation on a four-plate dynamic system; the third case study investigates a different four-plate dynamic system loaded by a white noise random point loading. The case study is based on the simulation of built-up plate systems with linear and nonlinear translational and/or torsional springs. In all of the addressed case studies, the thin plate is homogeneous, isotropic, and linear elastic with Young's modulus E = 70 GPa; Poisson's ratio ν = 0.3; and density ρ = 2700 kg/m 3 . The plates' damping loss factors, modal densities, and sizes, including a and b (plate's sides), as well as h (thickness), are given in the Table 1. The linear and nonlinear translational spring elastic coefficients are given as k l = 2 × 10 5 and k nl = 2 × 10 15 , respectively. Those for torsional spring are given as k θl = 10 3 and k θnl = 10 12 . With regard to the LRRM+MCS, the lumped masses are used to break the system symmetries inducing the Gaussian Orthogonal Ensemble (GOE) [21,33]. This paper applies 20 masses with 2.0% mass rate to the bare plate according to the authors' previous paper [21]. Besides, 50 Monte Carlo samples are utilized to calculate the ensemble average energy response. The schematic figure of first case study is shown in Figure 2. It includes three plates one of which is an with inclined angle α and excited by various loading types perpendicular to plate 1 middle surface. We consider several situations for this case studies for the purpose of finding what factors could effect the energy response of each subsystem within the built-up system. In the first situation, we explore the energy cascade through subsystem by changing the inclination angles. The translational springs are set to be nonlinear, while the torsional ones are linear. Different inclination angles, e.g., 0 • , 30 • , 60 • , and 80 • , are applied to plate 1, the external excitation applied to plate 1 is considered to be an harmonic point load. Figure 3 depicts the ensemble average energy for both linear and linearised FE-SEA and LRRM+MCS analysis with α = 0. In the figure, the average energy of LRRM+MCS analysis fluctuates dramatically in lower-frequency range but tends to keep stable and close to the response obtained by using the hybrid FE-SEA formulation in higher-frequency range. This is because lower-frequency modes are hard to be randomized by uncertainties, due to the fact that the energy response in low-frequency range is mainly influenced by resonant modes. The higher-frequency modes are, instead, more affected by the randomization induced by the lumped masses, which leads to the mixing and veering of the modes and then to the Rayleigh distribution. It is noted that the energy responses obtained via LRRM+MCS analysis compare well with those computed through the hybrid FE-SEA method for both linear and linearised formulation. This can also be seen in Figure 4, which shows the linear and linearised FE-SEA and LRRM+MCS analysis of plate 2 and plate 3 at different inclination angles. In this figure, it should be noted that, when the angle varies from 0 • to 30 • , the average energy response just slightly decreases, whilst, from 60 • to 80 • , a significant reduction is observed.     The second scenario of the first case study explores the influence by the spring position. The inclination angle is set to be zero, and the translational springs are linear and torsional springs are nonlinear. Three conditions in the spring position are considered: (i) the centre of plate 1 [coordinate (0.5a,0.5b)]; (ii) a remote position from the centre [coordinate (0.05a, 0.5b)]; and (iii) a random position in every MCS sample. The loading is still set as the harmonic point excitation. Results are shown in Figure 5 for both linear and linearised FE-SEA and LRRM+MCS analysis. As expected, all the energy responses calculated by means of the hybrid FE-SEA approach do not change prominently. This is because the FE-SEA method randomizes the spring position and estimates the average results; in other words, no specific positions of the joints are required. It can also be noted an excellent match between the linearised FE-SEA formulation and the benchmark model. Next, investigation within the case study 1 considers different values of the springs stiffness coefficients. As the exploration to translational spring stiffness coefficients has been made in the previous work [21], this investigation focus on the torsional spring. We separately increase the linear and nonlinear stiffness coefficients of the torsional springs, while keeping the harmonic point loading as external excitation. Figure 6a is obtained by increasing the linear coefficients from 10 to 10 2 , 10 3 , and 10 4 , respectively. The energy response by linearised FE-SEA which matches the benchmark model increases very smoothly with the rise of linear coefficients. In Figure 6b, the increase of the nonlinear stiffness coefficient from 10 8 to 10 9 , 10 10 , and 10 11 generates an energy level rise in different frequency range. Smaller nonlinear stiffness coefficient values influence the energy level in lower-frequency range, while the larger ones effect the higher-frequency range. A similar trend was obtained in a previous work for focused on translational springs [21].
The built-up system in Figure 2 is now considered to be excited by a white noise loading on plate 1. The statistical linearisation is the one used to derive the linearised FE-SEA formulation. For the white noise point excitation, Figure 7 presents both the linear and the linearised results computed by using both the hybrid FE-SEA method and the benchmark model. A good match can be seen between two different analyses. We also considered another situation: rain-on-the-roof excitation on plate 1. The energy responses can be found in Figure 8. Besides the good agreement between the linearised hybrid FE-SEA model and LRRM+MCS formulation, the results yielded by the benchmark model are smoother than those with point random loading, and the MCS sample cloud of rain-on-the-roof is thinner than those of the point random load, due to that fact that rain-on-the-roof is evenly distributed on the surface of the plate, and it can help realize better randomization for the modes.
In Figure 9, the system energy responses with both point random load and rain-on-the-roof as external excitations for both linear and linearised hybrid FE-SEA formulation are depicted. It can be noted that the energy level of plate 2 around 4000 rad/s excited by rain-on-the-roof remains steady, while those evaluated by using the point random load show some oscillations.

Case Study 2: Distributed Loading
This case study focuses on the harmonic distributed loading on the built-up plate system. The four-plate system with both translational and torsional springs, schematically shown in Figure 10, with localized nonlinearity in the spring set 1, is investigated. The distributed loading is set to excite plate 1 orthogonally to the middle surface. Different loading areas are applied in order to explore its effects on the energy response. The loading area varies from the small value of 0.2 × 0.2 to the larger of 0.6 × 0.6, 0.8 × 0.8, 1 × 1, where the case 1 × 1 means that the distributed harmonic load area equals the plate surface. The average energy responses related to this latest case are depicted in Figure 11. The energy response of the linearised analysis increases comparing to those of linear analysis for the reason of cubic harden stiffness. In Figure 12, energy response of plate 2 and plate 3 for different loading area on plate 1 is shown. It can be observed that a larger gap between the energy responses with loading area 0.2 × 0.2 and 0.6 × 0.6 occurs.

Case Study 3: Four-Plate Built-Up System
To further test the linearised FE-SEA formulation towards random loading, a more complex case study consisting of four-plate built-up system, shown in Figure 13, is addressed. All the plate parameters can be found in Table 1. The white noise point load in the figure orthogonally excites plate 1. Four cases are considered: (i) all spring sets are linear; (ii) only the first spring set is nonlinear; (iii) only the second spring set contains nonlinearity; (iv) only the third spring set is nonlinear. The energy responses from both linear and linearised FE-SEA and LRRM+MCS analysis are shown in Figure 14.
Comparing the linear energy response given in Figure 14a with those of the second case shown in Figure 14b, the energy response of plate 2 significantly increases as the nonlinearity is applied, while those of the other plate subsystems are only slightly affected. A very similar result can be observed by comparing Figure 14a with Figure 14c, where only the energy level of plate 3 ramps up significantly due to the nonlinearity in second spring set. However, the nonlinearity existing in third spring set changes the energy response in a very different manner with respect to the previous two cases. Figure 14d demonstrates that: (1) the energy level of both plate 3 and 4 steps up remarkably due to the nonlinearity; (2) a cross of the curve of energy level of plate 2 and plate 3 can occur. Moreover, an excellent match between linear and linearised FE-SEA method and LRRM+MCS analysis is presented in all the faced cases.

Conclusions
The present article proposes a linearised hybrid FE-SEA formulation for the dynamic response of build-up systems featured by nonlinear joints and subjected to both harmonic and random excitations. The formulation was validated by developing a benchmark model based on the combination of both the Lagrange-Rayleigh-Ritz method and the Monte Carlo Simulation technique. Within the framework of the benchmark model, each plate subsystem of the dynamic system is modelled by using Kirchhoff's thin plate theory. The two different linearisation procedures are used according to the external excitation type. More specifically, in the case of harmonic excitation, the method of harmonic balance was employed; in the case of random excitation, the statistical linearisation was used. Various case studies were examined to both validate and assess the new hybrid FE-SEA formulation. From all the analyses carried out, the following main conclusions can been drawn:

•
The plate inclination angle within the built-up systems slightly affects the energy response for small values; on the contrary, its effect tends to be prominent for inclination angle close to 90 • .

•
The springs' position-acting as joint components-, as expected, does not affect the energy response.

•
Larger values of the cubic nonlinear stiffness coefficients of the torsional springs increase the energy level in a wider frequency range affecting also the higher frequency.

•
Comparing the random point load with the rain-on-the-roof excitation can realize better randomization from the perspective of the LRRM, namely the energy level of the rain-on-the-roof tends to be closer to that of FE-SEA formulation. In all of the addressed case studies, the MHB and the SL, employed in both the hybrid FE-SEA formulation and the benchmark model, turned out to be highly effective in the linearisation process of built-up systems with localized nonlinearity.