Olivier ´ezvan

Xiaoshu Xiaoshu Zeng

Roger Ghanem

Bora Gencturk
email: gencturk@usc.edu

Olivier Ezvan

| Xiaoshu Zeng

| Roger Ghanem

Dominant vibration modes for broadband frequency analysis of multiscale structures with numerous local vibration modes

Keywords: Structural dynamics, model reduction, modal analysis, dominant vibration modes, local vibration modes, nuclear fuel assembly

published or not. The documents may come

Dominant vibration modes for broadband frequency analysis of multiscale structures with numerous local vibration modes

INTRODUCTION

This research deals with the dynamics of structures in the framework of finite element (FE) analysis [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. This paper focuses on linear structural dynamics with the assumption of infinitesimal strain. In this case, the structure eigenfrequencies and eigenmodes (vibration modes, or elastic modes) are fundamental to characterize its dynamic behavior [START_REF] Craig | Fundamentals of Structural Dynamics[END_REF]. In the linear case, the structure dynamics can be fully described in the frequency (Fourier) domain by the frequency response functions (FRFs). In the low-frequency (LF) range, the FRFs present sharp and distinct resonance peaks, each one corresponding to a given elastic mode. In the highfrequency (HF) range, the eigenfrequencies are closely spaced, with approximately constant average density. Thus, in this range, the FRFs are devoid of resonance peaks and present rather smooth variations. In the intermediate band, the medium-frequency (MF) band, the modal density is erratic and the FRFs exhibit numerous peaks that overlap [START_REF] Ohayon | Structural acoustics and vibration[END_REF]. Modal analysis, based on the projection of the FE model onto a reduced-order basis (ROB) of elastic modes, is very well adapted to the study of the LF range, since the dynamic behavior can be described by superposition of the contribution of each mode (and since the modal density is low in the LF range). A reduced-order model (ROM) is obtained, for which the coordinates are the coefficients of the modal expansion. In contrast, for HF analysis, the modal approach is not as attractive and statistical energy analysis [START_REF] Lyon | Theory and Application of Statistical Energy Analysis[END_REF] (SEA) is generally preferred. For the MF band, SEA is no longer valid and FE analysis is costly. Various hybrid methods mixing the modal and energetic approaches have been proposed [START_REF] Morand | A modal hybridization method for vibroacoustic studies at medium frequencies[END_REF][START_REF] Ladevèze | A new computational approach for structure vibrations in the medium frequency range[END_REF][START_REF] Langley | A hybrid method for the vibration analysis of complex structural-acoustic systems[END_REF][START_REF] Zhang | The discontinuous enrichment method for elastic wave propagation in the medium-frequency regime[END_REF][START_REF] Ji | A mode-based approach for the mid-frequency vibration analysis of coupled long-and short-wavelength structures[END_REF].

In this research, we are interested in predicting the dynamic behavior of complex multiscale structures. These structures are characterized by the presence of numerous substructures attached to a main body or frame. The substructures are small-scale flexible components with fundamental eigenfrequency in the same range as that of the whole structure. As a consequence, these complex multiscale structures are characterized by the presence of numerous local modes in addition to the few usual global modes. The global modes consist of ensemble deformations of the structure skeleton (main body), whereas the local modes exhibit localized deformations of one or several substructures while the structure as a whole remains still. In such a context, modal analysis in the LF range becomes more expensive, since the modal density is increased. In other words, the modeling of the small-scale complexity results in a high number of modes. Even though the use can be made of a detailed FE model that represents the small-scale structural complexity, the interest might not be the prediction of the local behavior of the small scales, but still be the global dynamic behavior. Being able to predict the local dynamics of the small scales is different and more challenging than being able to predict the global dynamics of the whole structure. It could be argued that the prediction of the global dynamics do not require the fine representation of the small-scale structural complexity, but it is not necessarily true. In addition, the structure skeleton for which one would be interested in predicting the global dynamics could be constituted of a complex assembly of small-scale components. Or, one might be interested in assessing the sensitivity of the global behavior with respect to variations of the small scales. Complex multiscale structures can be encountered in automotive (cars), nuclear (fuel assemblies), naval (complex ship structures) or aerospace (large space structures) engineering. The most common method to tackle the present difficulty is to select the dominant modes based on the mass participation factors [START_REF] Chopra | Dynamics of Structures: Theory and Applications to Earthquake Engineering[END_REF]. This approach is widespread in civil engineering (for seismic analysis). The mass participation factors have been used for removing the insignificant modes from the dynamic analysis of pressurized water reactor (PWR) fuel assemblies [START_REF] Klymyshyn | Fuel Assembly Shaker Test Simulation[END_REF][START_REF] Jiang | Spent nuclear fuel system dynamic stability under normal conditions of transportation[END_REF]. It should also be noted that the presence of numerous local modes makes experimental modal identification difficult. In addition, physical interpretation of the modal data can be very difficult even in a numerical setting. For instance, it is possible that none or a very few of the mode shapes of an automobile exhibit ensemble (global) deformations. This is due to the fact that the global deformations may be coupled with local deformations with very high relative amplitudes. Recently, a general methodology has been developed to precisely address the present problem, thanks to the introduction of a ROB of global modes [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF]. This global ROB is made up of the eigenmodes of a system for which the representation (kinematics) of the kinetic energy is altered, so as to filter the local modes. The methodology has been applied to automotive vehicles [START_REF] Arnoux | Stochastic reduced order computational model of structures having numerous local elastic modes in low frequency dynamics[END_REF] and PWR fuel assemblies [START_REF] Batou | Reduced-order computational model in nonlinear structural dynamics for structures having numerous local elastic modes in the low-frequency range. Application to fuel assemblies[END_REF]. It has also been extended to multilevel reduction [START_REF] Ezvan | Multilevel model reduction for uncertainty quantification in computational structural dynamics[END_REF]. The present research is a further improvement and extension of this method.

The objective of this research is to extract the dominant and representative small-dimension subspace that describes the global structural dynamics behavior of complex multiscale structures. For a fixed accuracy of the associated ROM, two aspects are of interest. First, that the ROM be of the smallest dimension as possible. Second, that its construction be efficient and unaffected by the presence of the numerous local elastic modes. This research enables practical FE analysis of complex structures in a broad frequency range that is characterized by the presence of heterogeneous and high modal density.

In order to build such a small-dimension ROM, an automatic procedure for the selection of the most significant elastic modes is proposed. For that, an error indicator is introduced to define what is intended by an optimal ROM. Since the mode selection procedure still requires the computation of the numerous elastic modes, an alternative methodology for constructing, a priori, the small-dimension ROB is proposed, based on the previous research aforementioned and devoted to the construction of a global ROM. The present developments assume a scale separation between the structure skeleton and the flexible components. The concept of static condensation is used to filter the local vibrations. In order to take into account the inertia of the flexible components and the influence of their local vibrations on the global vibrations, a reduced kinematics is introduced.

In Section 2, modal analysis in linear computational structural dynamics is recalled, needed to introduce the framework of the present research. In Section 3, the automatic mode selection strategy is described and implementation details are given. In Section 4, the methodology for constructing the global ROM is presented. Static condensation is recalled, followed by a general method for constructing a global ROM. Several kinematics based on scale separation are investigated. Then, static condensation is combined with a reduced kinematics for representing the coupling of the global dynamics with the inertia and local dynamics of the small scales. The details for an efficient implementation are described. Finally, in Section 5 the proposed methods are applied to three examples, namely a heterogeneous plate, a simplified nuclear fuel assembly, and a detailed boiling water reactor fuel assembly.

MODAL ANALYSIS THEORY AND DIFFICULTY ENCOUN-TERED FOR MULTISCALE STRUCTURES

In this paper, we are interested in solving the structural dynamic problem in the linear regime by using the FE method, under the assumption of infinitesimal strain. We consider the formulation in the frequency domain and aim at predicting the displacement field for all circular frequency ω belonging to the frequency band of analysis, B ω = 2π×]0, f max ], with f max the upper bound of the frequency band in Hz. The vector U(ω) denotes the discretization of the displacement field over the bounded structural domain Ω and is the solution of the following matrix equation, for all ω in B ω ,

-ω 2 [M] + iω [D] + [K] U(ω) = F(ω) . (1)
In Eq. ( 1), [M], [D], and [K] are the mass, damping, and stiffness matrices, vector F(ω) is the discretization of the applied forces, and i denotes the complex number. The structure is assumed to be fixed on a part of its boundary and also, it is assumed that these three matrices (the system matrices) are real and positive definite. For the applications of interest, the number of nodal degrees of freedom (DOFs), that is denoted N and which is the dimension of the system matrices, can be very large (between a few million and a few hundred million) but the matrices are very sparse (only the DOFs of connected nodes give nonzero extradiagonal entries). As the mesh is refined to capture the structural complexity as well as to represent small-wavelength deformations associated with high frequencies, the cost of solving the linear system of Eq. ( 1), for fixed ω, rapidly increases. This equation has to be solved for a given sampling of frequency points, but possibly also for several distinct applied forces. The use of the structural vibration eigenmodes allows to replace the high-dimensional model by a ROM, which alleviates these two difficulties. For all α = 1, . . . , N , the vibration eigenmode (or elastic mode) ϕ α and associated eigenvalue λ α are the solutions of the generalized eigenvalue problem,

[K]ϕ α = λ α [M]ϕ α . (2)
The eigenvalues verify 0 < λ 1 ≤ λ 2 ≤ . . . ≤ λ N < ∞, the circular eigenfrequencies ω α are such that ω 2 α = λ α , and the eigenfrequencies f α (in Hz) are such that ω α = 2πf α . In practice, only the first n e N elastic modes, associated with the lowest frequencies, are calculated and the ROB is given by the matrix [Φ] = [ϕ 1 . . . ϕ ne ]. The elastic modes are real and can be arbitrarily normalized. In this paper, the unit mass normalization is systematically used, which gives

[Φ] T [M][Φ] = [I ne ] , (3)
in which [I ne ] denotes the identity matrix of dimension n e . The projection of the stiffness matrix is such that

[Φ] T [K][Φ] = [Λ] , (4)
in which [Λ] = diag(λ 1 , . . . , λ ne ) is the diagonal matrix of the first n e eigenvalues. In general, the projection of the damping matrix is not a diagonal matrix. A general practice for lightly-damped structures is to adopt a damping model for which the projected damping matrix is diagonal. In this paper, a modal damping model is considered, such that

[Φ] T [D][Φ] = 2[Ξ][Ω] , (5)
in which [Ω] = [Λ] and where [Ξ] = diag(ξ 1 , . . . , ξ ne ) is the diagonal matrix of the modal damping ratios ξ α > 0 . The classic ROM of modal analysis is obtained by approximating, for all ω ∈ B ω , the solution U(ω) of Eq. ( 1) by the following expansion,

U(ω) ne α=1 q α (ω)ϕ α = [Φ]q(ω) , (6)
in which q(ω) = (q 1 (ω), . . . , q ne (ω)) is a n e -dimensional complex vector of generalized coordinates q α (ω) .

That is, the dynamic response U(ω) is written as a linear combination of the first n e elastic modes. In other words, the dynamic response is approximated into the subspace spanned by the first n e elastic modes. The Galerkin projection of Eq. ( 1) onto the elastic modes leads to the reduced matrix equation,

∀ω ∈ B ω , -ω 2 [I ne ] + 2iω [Ξ][Ω] + [Λ] q(ω) = F(ω) , (7)
in which F(ω) = [Φ]

T F(ω). Since Eq. ( 7) is a diagonal matrix equation, it is solved very efficiently and consequently, the overall cost only pertains to the eigenvalue computation of Eq. [START_REF] Craig | Fundamentals of Structural Dynamics[END_REF]. In general, the modal ROB [Φ] is also used for design optimization, stochastic analysis, robust design, etc. It can also be used for nonlinear dynamical systems (e.g., with geometric nonlinearities) or viscoelastic dynamical systems. That is, the physical model is modified or perturbed with respect to the nominal high-dimensional model from which the ROB is built. In this case, the validity of the ROB has to be assessed for the given range of perturbations, and it may be needed to increase the number of elastic modes to be included in the ROB. In addition, the associated ROM is no longer given by diagonal matrix equations. In this case, the efficiency of the ROM is a function of the number n e of modes kept. The extended use of the modal ROB for such equations that differ from the deterministic linear conservative problem of Eq. ( 2) enables great computational savings. Nevertheless, sometimes the ROB may not remain valid and consequently, it may be needed to recompute it for each perturbation or for a sampling of perturbations of the system. As a consequence of these facts, in general, both the cost to build the ROB (here, Eq. ( 2)) and the dimension of the ROB are key factors for efficiency.

For complex structures exhibiting several structural scales, numerous elastic modes associated with isolated vibrations may be intertwined with the usual global elastic modes. More precisely, the short-wavelength local displacements of the flexible substructures may be mixed with the long-wavelength global displacements of the structure skeleton, yielding numerous additional modes, such that the number of modes in the frequency band of analysis may be very large. Very often, these numerous additional modes are insignificant for representing the global dynamics of the structure. These local modes lead to a high modal density, hence a high dimension n e for the ROB. In addition, as there are more modes to compute, the cost for building the ROB is also increased. In this paper, we thus address the problem of (i) constructing an accurate ROM with small dimension and the problem of (ii) avoiding the calculation of the numerous elastic modes of Eq. (2). As an illustration, Figure 1 depicts a heterogeneous plate with two structural scales, namely 12 flexible panels that are supported by a stiff skeleton. This simple case will be studied in Section 5.1. Figure 2 depicts the deformation shape of a local mode, where it can be seen that only one flexible panel undergoes a bending deformation. Figure 3 depicts the deformation shape of a global mode, where it can be seen that the whole structure undergoes a bending deformation. Figure 4 depicts the deformation shape of a multiscale mode, that is to say that is composed of both global and local deformations. This paper considers the case of lightly-damped linear systems. The high-dimensional FE model can be replaced, without loss of fidelity, by a ROM whose ROB is made up of elastic modes. The elastic modes must cover the frequency band of analysis and a little beyond, since the contribution of each mode is not exactly confined to its eigenfrequency but spread around it. In addition, it may be necessary to include more modes to take into account their contribution to the static response, although for broadband frequency analysis, the static response may converge without the use of additional higher-frequency modes. In the literature, the term "dominant modes" often refers to the elastic modes associated with the lowest eigenfrequencies. In contrast, in the context of the presence of insignificant local elastic modes, dominant elastic modes have, most of the time, been defined based on their modal participation factors.

The modal participation factor Γ i,α for mode ϕ α and direction i is given by

Γ i,α = ϕ T α [M]r i ϕ T α [M]ϕ α , (8)
in which the influence vector r i is composed of ones for the DOFs associated with direction i and of zeros otherwise. The effective mass Γ 2 i,α for mode ϕ α and direction i is such that

m tot = m α=1 Γ 2 i,α , (9)
in which m tot denotes the total mass of the structure. Finally, the participation mass ratio (or fraction) for mode ϕ α and direction i,

Γ 2 i,α /m tot , is such that m α=1 Γ 2 i,α m tot = 1 . (10)
Thereby, the modal participation factors (or equivalently, the effective masses or the participation mass ratios) can be used to rank the elastic modes by order of importance. It should be noted that the participation factors are independent of the applied forces. However, they are dependent on the choice of mode normalization. In addition, the participation factors are zero for antisymmetric modes. Another remark is that very often the cumulative mass fraction (the truncated sum of Eq. ( 10)) practically reaches unity by only including the low frequency modes, which contradicts the fact that a prediction in a broad frequency band most likely requires to include modes located in its HF part. Therefore, in the next section, a new criterion is introduced to define the mode importance.

FIRST METHOD PROPOSED: AUTOMATIC SELECTION OF THE DOMINANT VIBRATION MODES

In this section, an alternative to the participation factors is presented to rank the importance of each mode.

The importance criterion is used to discard the elastic modes that are insignificant for representing the global dynamic response of the structure skeleton. It should be noted that such modes may still be very important to represent the local dynamic response of the flexible substructures. Therefore, the upcoming definition of modal importance is not as general as that of the modal participation factors, as it is goal oriented. Denoting by

[A(ω)] = -ω 2 [M] + iω [D] + [K]
the dynamic stiffness matrix, the generalized coordinates vector q(ω) solution of Eq. ( 7) can be written as

q(ω) = [Φ] T [A(ω)][Φ] -1 [Φ] T F(ω) . (11)
The (n e × n e ) diagonal complex matrix [

H(ω)] = [Φ] T [A(ω)][Φ] -1 is introduced, whose components h α (ω)
for α = 1, . . . , n e are given by

h α (ω) = -ω 2 + 2 iωξ α ω α + ω 2 α -1 . ( 12
)
The linear system (Eq. ( 1)) can be entirely described by the FRF matrix, whose entry i, j is the response (displacement) at the DOF i due to a unit load applied to the DOF j. The displacement induced by a unit load is called compliance. The response to any arbitrary load can be deduced from the FRF matrix, by linearity. From Eqs. [START_REF] Ladevèze | A new computational approach for structure vibrations in the medium frequency range[END_REF] and [START_REF] Klymyshyn | Fuel Assembly Shaker Test Simulation[END_REF] it can be seen that the (N

× N ) complex matrix [Φ][H(ω)][Φ]
T is the order-n e approximation of the FRF matrix. It is assumed that the order n e is sufficient to have reached convergence in frequency band B ω and consequently, that such approximation of the FRF matrix is exact for all ω in B ω . Consequently, entry i, j of the FRF matrix, denoted U ij (ω), can be calculated by

U ij (ω) = ne α=1 U ijα (ω) , (13)
in which U ijα (ω) is the contribution of mode ϕ α to the response (to the compliance), which can be written as

U ijα (ω) = Φ iα h α (ω) Φ jα . (14)
In Eq. ( 14), scalar quantity Φ kα is given by Φ kα = [Φ] kα , for all k = 1, . . . , N and α = 1, . . . , n e . We are interested in analyzing the modulus in dB scale of the FRF, which is given by

u ij (ω) = 20 log 10 |U ij (ω)| . (15)
In the sequel of the paper, by FRF we mean its modulus. The error measure (u, v) associated with the FRFs in dB scale u(ω) and v(ω) is introduced, with the expression

(u, v) = 1 |B ω | Bω (u(ω) -v(ω)) 2 dω . ( 16
)
Error measure is the root mean square error between the FRFs in dB scale u and v over the frequency band of analysis B ω . It thus gives the error level in dB. The error measure has to be used for different combinations of input j and output i, for which the magnitude of the compliance U ij (ω) may vary; the error measure sought should not depend on the magnitude of the response U ij , that is to say it should be normalized. Since the arguments of function are FRFs in dB scale, it can be verified that the error measure is independent of the magnitude of the response, that is to say the error measure is normalized. It should also be noted that, thanks to the dB scale, the error is independent of considering either the displacement, velocity or acceleration (or, respectively, the compliance, mobility or accelerance). The error measure is introduced to quantify the error introduced by the use of a ROM.

Optimal mode selection

Let u

(n) ij (ω) u ij (ω) be the FRF in dB scale, for output observation i and input excitation j, calculated by an optimal ROM of dimension n. The objective of such a ROM is to accurately represent the global dynamics of the structure skeleton, with no interest in representing the local dynamics of the flexible substructures. Therefore, we are only interested in input and output locations that are not support of local vibrations and consequently, that are not primarily sensitive to the local dynamics of the flexible substructures. For instance, such locations can belong to the structure stiff skeleton (or main frame), if it exists. Let Ω G denote a subdomain of Ω that is not the support of local vibrations of flexible substructures, or more generally, a subdomain for which the user is interested in predicting the dynamic response. It should be noted that for some complex structures without scale separation, it may be impossible to identify a subdomain that would not be part of the support of any of the local vibrations. In this paper, we restrain ourselves to the study of complex multiscale structures with scale separation. Let S G denote a given subset (sampling) of nodal DOFs belonging to subdomain Ω G . The requirement for the optimal ROM is that it be of minimal dimension n under the constraint that

∀i, j ∈ S G × S G (u ij , u (n) ij ) < tol , (17)
with tol a prescribed error tolerance. It means that the ROM must be valid for all the input and output locations considered. Therefore, the maximum error over all the combinations of output i and input j must be minimized. In the sequel, for brevity, by error it is meant the maximum error over all the combinations of output i and input j. The search for the optimal ROM (or ROB) is restrained to the search of the optimal ROM for which the ROB is constituted of elastic modes. The problem is thus simplified to finding the optimal selection of elastic modes. Despite the FRF U ij being a linear function of the n e terms Φ iα Φ jα for α = 1, . . . , n e , the error function is not a linear function of said terms. Consequently, the error in removing a given set of modes is not equal to the sum of the individual errors associated with removing separately each of these modes. The optimal ROM of size n e -k is given by removing the set of k modes that induces the smallest error. The number of possible different sets of k modes is extremely large even for a reasonably small n e . It is given by n e !/(k!(n e -k)!). Therefore, the search of the optimal selection in trying all the possibilities is highly intractable.

Tractable mode sorting procedure

The optimal ROM of size n e -1 is given by removing the one mode that induces the smallest error. The number of possibilities is n e . The optimal ROM of size n e -2 is given by removing the one set of two modes that induces the smallest error. The number of possible different sets of two modes is equal to n e (n e -1)/2. In order to decrease the number of cases to be considered, one possible approximation is to restrict the search to the sets of two modes that include the one mode that had been removed for obtaining the optimal ROM of size n e -1. The number of cases to be considered is thus reduced to n e -1. Likewise, a good ROM of size n e -3 (but not necessarily optimal) can be obtained by searching for the mode whose removal, in conjunction with the removal of the two modes of previous step, induces the smallest error. The number of cases to be considered for this is n e -2. Proceeding in the same manner for defining the next modes to be removed, the total number of cases to be considered for defining the k least important modes is k-1

j=0 n e -j = (k/2)(2n e -k + 1). For a significant model reduction, one has k ∼ n e and consequently, the number of cases to be considered is in the order of n 2 e . The present mode selection is optimal if it can be assumed, for all n = 1, . . . , n e -1 , that the optimal ROM of dimension n does not include any mode that is not included in the optimal ROM of dimension n + 1. This mode selection procedure is the same as that introduced in [START_REF] Wang | Component-Centric Reduced Order Modeling of the Dynamic Response of Linear Multibay Structures[END_REF] and devoted to the construction of compact ROBs for the high-fidelity modeling of specific sub-components of a structural dynamic system.

Simplified mode sorting procedure

Instead of sorting the modes in such a sequential manner, one further approximation is to neglect the fact that the error function is not a linear function of the contribution of each mode, and independently define the importance of each mode based on the error induced by its sole removal. For this simplified sorting procedure, the number of cases (number of ROMs) to be considered is only n e . It should be recalled that, for each ROM, said error refers to the maximum error over all the combinations of input load j and output observation i. An algorithm to perform such modal sorting is now presented.

Algorithm for defining the modal importance

• For each mode α = 1, . . . , n e , initialize the importance quantity, α , to zero: α = 0

• For each mode α = 1, . . . , n e , compute h α (ω) defined in Eq. ( 12). This can be stored in a (n e × n ω ) complex matrix, with n ω the number of sampled frequency points.

• Nested loops over the observation outputs i and excitation inputs j, with j ≥ i (from Eqs. [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF] and [START_REF] Arnoux | Stochastic reduced order computational model of structures having numerous local elastic modes in low frequency dynamics[END_REF] it can be deduced that U ij is symmetric):

-For all α = 1, . . . , n e , compute U ijα (ω) defined in Eq. ( 14). This can be stored in a (n e × n ω ) complex matrix.

-Compute U ij (ω) = ne α=1 U ijα (ω) (this is the exact FRF). -Compute u ij (ω) = 20 log 10 |U ij (ω)| (this is the associated modulus in dB).
-For all α = 1, . . . , n e , compute U ijα (ω) = U ij (ω) -U ijα (ω). This can be stored in a (n e × n ω ) complex matrix (these are the approximated FRFs obtained after removal of each mode ϕ α ).

-For all α = 1, . . . , n e , compute u ijα (ω) = 20 log 10 | U ijα (ω)| . This can be stored in a (n e × n ω ) real matrix (these are the associated moduli in dB).

-For all α = 1, . . . , n e , compute ijα = (u ij , u ijα ) with the error measure given in Eq. ( 16).

-For all α = 1, . . . , n e , update the importance quantity α by replacing it with max( α , ijα ).

Modal truncation

The modes are ranked according to the error induced by their individual removal, but the error induced by the combined removal of the k (so-defined) least important modes has not been evaluated. The basis of dominant modes is obtained by picking the n d = n e -k modes with highest importance α , with n d minimum under the constraint that

∀i, j ∈ S G × S G (u ij , u (n d ) ij ) < tol , (18)
with u

(n d ) ij (ω) the FRF in dB scale, for output observation i and input excitation j, calculated by the ROM of the n d dominant modes. The error function (u ij , u

(n d ) ij ) is not necessarily a monotonic function of n d .
As a consequence, one is not guaranteed of having found the smallest n d satisfying Eq. ( 18) (the optimal n d ) without having tried all the n e possibilities. In practice, error function (u ij , u

(n d ) ij )
is roughly monotonic and one may be satisfied with considering an n d deduced from the evaluation of a limited sample.

The presented ROB of dominant elastic modes addresses the problem of dimension reduction (dimension of the associated ROM) but it does not address the problem of the unusually-high cost for constructing this ROM, which is due to the large value of n e . Therefore, in the next section a methodology is presented, which aims at directly constructing the small-dimension ROB, without having to compute the n e elastic modes via Eq. ( 2).

SECOND METHOD PROPOSED: CONSTRUCTION OF A GLOBAL-DISPLACEMENTS BASIS

In order to construct, a priori, the small-dimension ROB, the strategy is to introduce a filtering of the local deformations (or displacements). In other words, the objective is the construction of basis vectors whose deformation shapes consist of ensemble deformations of the structure. These basis vectors, which will be used instead of the elastic modes, are called "global modes" and form a ROB that is called "global-displacements ROB" or "global ROB". The construction of said global modes is based on the introduction of a modified eigenvalue problem with respect to that of Eq. ( 2). The solutions of this eigenvalue problem are intended to be devoid of the numerous local modes. Prior to presenting the methodology for the construction of the global ROB, the concept of static condensation [START_REF] Chopra | Dynamics of Structures: Theory and Applications to Earthquake Engineering[END_REF] is recalled, useful for contextualization.

Static condensation

Often in the literature, the terms Guyan condensation [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF] and static condensation are used interchangeably. In the sequel, the difference between the two will be highlighted. Both are based on splitting the DOFs of the structure into two sets, namely the master DOFs and the slave DOFs. A ROM is then obtained, for which the coordinates are the master DOFs. The slave DOFs are eliminated. The general concept of condensation is now presented. Partitioning the generalized eigenvalue problem of Eq. (2), one gets

K mm K ms K sm K ss ϕ m ϕ s = λ M mm M ms M sm M ss ϕ m ϕ s , (19)
in which subscript m refers to the master DOFs and subscript s refers to the slave DOFs. From the second row of above equation, one gets

[K sm ]ϕ m + [K ss ]ϕ s = λ ([M sm ]ϕ m + [M ss ]ϕ s ) from which it can be deduced that ϕ s = -([K ss ] -λ[M ss ]) -1 ([K sm ] -λ[M sm ]) ϕ m ≡ [t sm (λ)] ϕ m . (20)
In the above equation, the slave DOFs are expressed in terms of the master DOFs thanks to the condensation matrix [t sm (λ)]. As a consequence, an eigenmode ϕ solution of Eq. ( 19) (which only differs from Eq. ( 2) by a reordering of the DOFs) can be written as ϕ = [T (λ)]ϕ m , in which the transformation matrix [T (λ)] is given by

[T (λ)] = I m t sm (λ) , (21)
with [I m ] denoting the identity matrix of dimension equal to the number of master DOFs. Performing the Galerkin projection of Eq. ( 19) onto the column space of matrix [T (λ)], an eigenvalue problem of reduced dimension is obtained,

[K R (λ)]ϕ m = λ[K R (λ)]ϕ m , (22)
in which

[K R (λ)] = [T (λ)] T [K][T (λ)] and [M R (λ)] = [T (λ)] T [M][T (λ)
] are (m × m) symmetric positivedefinite real matrices that depend on λ . Hence, the above eigenvalue problem is nonlinear in λ . Various condensation methods have been introduced [START_REF] O'callahan | A procedure for an improved reduced system (IRS) model[END_REF][START_REF] Xia | A new iterative order reduction (IOR) method for eigensolutions of large structures[END_REF], that consider different approximations for [T (λ)] or schemes for solving Eq. [START_REF] Golub | Matrix Computations[END_REF]. The goal of condensation was to reduce the cost involved in solving Eq. ( 2), with respect to both memory and time. The Guyan condensation is obtained by substituting [T (λ)] with [T (0)]. The induced approximation depends on the choice of the master DOFs and in general, it is only valid at low frequencies. Moreover, for a poor selection of the master DOFs, not only the eigenfrequencies and eigenmodes obtained by Guyan condensation can exhibit large errors, but also some modes can be missed. From Eq. ( 20) it can be seen that the approximated eigenmodes obtained by Guyan condensation satisfy the following relationship,

ϕ s = -[K ss ] -1 [K sm ]ϕ m , (23)
which means that the master DOFs and the slave DOFs are in static equilibrium. The eigenmodes ϕ approximated by Guyan condensation are expressed in the physical coordinates of the FE model by the equality ϕ = [T (0)]ϕ m . The difference between the Guyan condensation and the static condensation is now highlighted. The static condensation is obtained by eliminating the mass of the slave DOFs, that is to say in Eq. ( 19) matrices [M ms ], [M sm ], and [M ss ] are replaced by zero matrices. Replacing these matrices by zero matrices in Eq. ( 20) leads to Eq. ( 23), which means that the static condensation implies the static equilibrium between the master DOFs and the slave DOFs, just like the Guyan condensation. From the first row of Eq. ( 19) and for the case of eliminated slave mass, one gets

[K mm ]ϕ m + [K ms ]ϕ s = λ[M mm ]ϕ m ,
which, substituting ϕ s by its expression given in Eq. ( 23), yields the reduced-dimension eigenvalue problem,

[K G ]ϕ m = λ[M mm ]ϕ m , (24)
in which

[K G ] = [K mm ] -[K ms ][K ss ] -1 [K sm ] (Schur complement
) is the Guyan reduced stiffness matrix, that is to say it can be verified that [K G ] = [K R (0)]. In contrast, the Guyan reduced mass matrix,

[M G ] = [M R (0)], is given by [M G ] = [M mm ] -[M ms ][K ss ] -1 [K sm ] -[K ms ][K ss ] -1 [M sm ] + [K ms ][K ss ] -1 [M ss ][K ss ] -1 [K sm ] = [M mm ].
The eigenvalue problem of Eq. [START_REF] Aktulga | Parallel eigenvalue calculation based on multiple shift-invert Lanczos and contour integral based spectral projection method[END_REF] does not proceed from a Galerkin projection of Eq. ( 2), as the reduced matrices are projected onto different subspaces, associated with the column space of [T (0)] on one hand and with the column space of the indicator matrix [χ m ] on the other hand, the latter being defined as

[χ m ] = I m 0 . (25)
Therefore, unlike the eigenmodes approximated by Guyan condensation, the eigenmodes approximated by static condensation diagonalize the mass matrix for which the slave mass has been eliminated but not the actual mass matrix, [M]. For this reason, such approximate modes (which can be seen as nothing more than arbitrary basis vectors) are denoted by letter ψ rather than ϕ in the sequel. It should be noted that, similarly to the modes approximated by Guyan condensation, these modes approximated by static condensation are expressed in the physical coordinates of the FE model by the equality ψ = [T (0)] ϕ m , due to Eq. [START_REF] Bennighof | An automated multilevel substructuring method for eigenspace computation in linear elastodynamics[END_REF]. Let [Ψ s ] denote the (N × n s ) real matrix whose n s columns are the first n s ≤ m approximated modes ψ by static condensation. Performing the Galerkin projection of Eq. (2) onto the column space of [Ψ s ] yields the reduced-dimension eigenvalue problem,

[Ψ s ] T [K][Ψ s ] s = λ [Ψ s ] T [M][Ψ s ] s , (26)
from which the eigenmodes ϕ approximated by static condensation, but associated with the "true" system given by [K] and [M], are obtained by ϕ = [Ψ s ] s . In addition, it should be noted that the eigenvalues of Eq. ( 26) differ from the eigenvalues of Eq. [START_REF] Aktulga | Parallel eigenvalue calculation based on multiple shift-invert Lanczos and contour integral based spectral projection method[END_REF]. Such an approximation of the elastic modes in a given subspace, as done in Eq. ( 26), is often referred to as "Rayleigh-Ritz analysis".

As mentioned earlier, condensation was introduced to decrease the time and memory cost associated with the eigenvalue computation of Eq. ( 2). In such a context, the number m of master DOFs is small, and the Guyan stiffness matrix [K G ], despite being dense (unlike the stiffness matrix [K]), can be assembled and stored at a low cost. However, our objective is not precisely to reduce the cost associated with solving Eq. ( 2), but to introduce a filtering of the insignificant local modes. Therefore, in the sequel it will be considered to introduce a well-chosen selection of master DOFs so that the local modes are missed. More precisely, by dividing the structure into two parts, one of which being referred to as the flexible part and that is the support of local displacements intended to be filtered, the slave DOFs are defined as all the DOFs of said flexible part and the master DOFs as all the DOFs of the complementary stiff part. In such a context, dimension m is such that m ∼ N rather than m N .

General methodology for constructing the global modes Let S r denote a given subspace of R N (i.e., of the Euclidean space) that is spanned by the columns of a (N × r) ROB, [B]. Here, subscript r stands for "reduced", as in "reduced kinematics" and the quantity r denotes the dimension of subspace S r .

The introduction of such a reduced kinematics (reduced representation) is intended to imply a filtering of the insignificant local displacements. It should be noted that for the Guyan and static condensations, the representation of the elastic energy is kept exact. For most model reduction techniques in structural dynamics, the representation of the elastic energy is kept exact whereas the representation of the kinetic energy is approximated. For instance, the well-known Craig-Bampton substructuring technique [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] considers a Galerkin projection of Eq. ( 2) onto the subspace spanned by the static constraint modes and the component eigenmodes, and is such that the representation of the elastic energy is kept exact in the case that all the component eigenmodes are removed (in practice, the component eigenmodes are truncated with respect to the frequency). Therefore, in order to construct the global-displacements basis, a general framework is presented, that allows the use of any arbitrary reduced representation (reduced kinematics) for the kinetic energy, while keeping an exact representation for the elastic energy.

The Euclidean space is equipped with an inner product such that, for all vectors u and v in R N , their product u, v W is given by

u, v W = v T [W]u , (27)
where [W] is a (N × N ) symmetric positive-definite real matrix that will be referred to as the weight matrix (it serves as a metric matrix). In the sequel, [W] can also be a positive-semidefinite matrix, in which case, the inner product is said to be degenerate. For all u ∈ R N , let u r denote the orthogonal projection of u onto subspace S r . It is given by u r = [P]u , in which the orthogonal-projection matrix [P] is a (N × N ) real matrix that can be written as

[P] = [B] [B] T [W][B] -1 [B] T [W] . (28)
On the other hand, the kinetic energy E k ( u) associated with any time-dependent velocity vector u is given by

E k ( u) = 1/2 uT [M] u. The kinetic energy E r k ( u) = E k ( ur ) associated with the orthogonal projection ur = [P] u of u is given by E r k ( u) = 1/2 ( ur ) T [M] ur = 1/2 uT [M r ] u, in which the (N × N ) low-rank mass matrix [M r ] is written as [M r ] = [P] T [M][P] . (29)
Matrix [M r ] is a (N × N ) symmetric positive-semidefinite real matrix of rank r. The following eigenvalue problem is introduced, which differs from Eq. ( 2), and for which the eigenvectors {ψ α } α=1,...,r are not the elastic modes,

[K]ψ α = σ α [M r ]ψ α . (30)
In Eq. ( 30), σ α is the eigenvalue associated with the eigenvector ψ α . This eigenvalue problem yields only r finite eigenvalues, which are such that 0 < σ 1 ≤ σ 2 ≤ . . . ≤ σ r < ∞ . Similarly to the static condensation presented in the beginning of this section, above eigenvalue problem is not obtained by a Galerkin projection of Eq. ( 2) and consequently, an additional Rayleigh-Ritz step is necessary to obtain the global modes.

Denoting as [Ψ(ν)] the (N × ν) real matrix such that [Ψ(ν)] = [ψ 1 . . . ψ ν ] , the following small-dimension eigenvalue problem is introduced,

[Ψ(ν)] T [K][Ψ(ν)] r α = λ g α [Ψ(ν)] T [M][Ψ(ν)] r α , (31)
from which, for α = 1, . . . , ν , the global mode ϕ g α with associated eigenvalue λ g α is obtained as ϕ g α = [Ψ(ν)] r α . Due to the orthogonality of the eigenvectors r α and taking into account the unit mass normalization, it can be deduced that the global modes ϕ g α are orthogonal with respect to both the mass and stiffness matrices.

Introducing [Φ g ] = [ϕ g 1 . . . ϕ g ng ] the global ROB and [Λ g ] = diag(λ g 1 , .
. . , λ g ng ) the diagonal matrix of the first n g ≤ ν global eigenvalues, the orthogonality reads For each column of matrix [B], the nodal DOFs with a nonzero entry have nonzero coupling terms in the low-rank mass matrix [M r ], by construction. For instance, if a fully-populated vector of ones is included in [B], then matrix [M r ] is fully populated. The use of global shape functions, that is to say deformation shapes whose support is a large part of the structure, can be useful to construct accurate global modes while filtering the insignificant local modes. However, in this case, (N × N ) matrix [M r ] is thus dense. Consequently, in general, it may not be possible to assemble and store [M r ]. Nevertheless, the solution of the eigenvalue problem of Eq. (30) only requires products of [M r ] with right-hand side vectors, hence the assembly and storage of [M r ] is not necessary. However, in the case of dense [M r ], it may no longer be possible to perform a Sturm sequence check (inertia count, based on the LDL decomposition [START_REF] Golub | Matrix Computations[END_REF] of [K] -µ[M r ] for a given µ > 0) to predict the number of eigenvalues in the desired frequency band. As implied earlier, the eigenvalues σ α in Eq. ( 30) and λ g α in Eq. (31) are different; the FRF resonance peaks predicted by the global ROM are located at the eigenfrequencies f g α = √ λ g α /2π if the small shift induced by damping is neglected, whereas there is no such meaning for the eigenvalues σ α . Dimension n g of the global ROM can then be chosen, as usual, with respect to a cutoff frequency related to the upper bound 2πf max of the frequency band of analysis B ω . In contrast, parameter ν may not be chosen in the same manner, as the eigenvalues σ α may differ too much from their counterpart λ g α . Truncation ν can be seen as a parameter that controls the richness of the subspace in which the global modes ϕ g α belong. It is possible that by increasing the value of ν, the modal density of the global modes ϕ g α increase significantly; the fact that the eigenvectors ψ α associated with high eigenvalues σ α may increase the modal density is because their contributions combine to, possibly, generate LF modes, among which there would be many unwanted local modes. The more the eigenvectors ψ α would differ from the elastic modes ϕ α , the more they would combine to generate LF local modes. It has been verified numerically that

[Φ g ] T [M][Φ g ] = [I ng ] , [Φ g ] T [K][Φ g ] = [Λ g ] . (32)
∀α = 1, . . . , ν λ α ≤ λ g α ≤ σ α . (33)
In particular, it means that the global eigenvalues λ g α are smaller than the eigenvalues σ α for which truncation ν has to be defined. It has been observed that in general the amount of difference between λ g α and σ α is related to the ability of the reduced kinematics to represent the global dynamics in the frequency range around f g α . In case the reduced kinematics is satisfactory, the global ROB is not very sensitive to ν, provided that ν is large enough so that σ ν > (2πf cut ) 2 , with f cut ∼ f max , the cutoff frequency from which n g is deduced. The sensitivity with respect to ν can be described by the modal density of the global ROB (i.e., of the eigenvalues λ g α ) with respect to ν. The modal density can be estimated by inertia counts performed at given points of the frequency band. To this end, the LDL decomposition of [Ψ(ν)]

T [K][Ψ(ν)] -µ[Ψ(ν)] T [M][Ψ(ν)] for a given µ > 0 is very cheap because ν N .
A clarification of the meaning of the global modes ϕ g α is now given. For that, let us first recall that Eq. ( 1) is the discretization of the weak form associated with the boundary value problem that is aimed to be solved. The boundary value problem is the mathematical model, governed by partial differential equations, that describes the structure dynamics. The mathematically-exact elastic modes are the solutions of the spectral problem associated with the weak form of the boundary value problem. The discretization Eq. ( 1) of the weak form is obtained by Galerkin projection onto a finite dimensional subspace that is spanned by the FE shape functions (here, the dimension of this subspace is denoted by N , which is the number of DOFs). The term Galerkin projection refers to the fact that both the solution and the test function associated with the weak form are approximated into the same subspace, which yields symmetric system matrices. Hence, the elastic modes ϕ α are the best approximations to the elastic modes, belonging to the FE subspace (the subspace spanned by the FE shape functions). Likewise, it can be seen that the global modes ϕ g α are the best approximations to the elastic modes, belonging to the subspace spanned by the ν vectors ψ α . As a conclusion, terminology "mode" for designating ϕ g α is justified, as only the approximation subspace differs. It should be noted that, since the subspace spanned by the ν vectors ψ α is a subspace of the FE subspace, the global modes ϕ g α are no better than the elastic modes ϕ α in approximating the mathematically-exact elastic modes. In contrast to the global modes ϕ g α , the ν vectors ψ α are not obtained by Galerkin projection of the spectral problem, as different subspaces are used to build [K] and [M r ]. Terminology "mode" is thus not employed for designating ψ α . It is recalled that in the beginning of this section related to static condensation, regarding notation, the same distinction had been made between vectors ψ and ϕ.

In the next sections, a description of a few relevant examples of choice for the reduced kinematics is given. Each reduced kinematics is defined by both matrix [B] (for which the column space is subspace S r ) and weight matrix [W] (which defines the inner product for the orthogonal projection of Eq. ( 28)). For each case, an implementation of the product of matrix [M r ] with a right-hand side vector is provided, which is the building block for solving the eigenvalue problem of Eq. (30). The assembly of matrices [P] and [M r ] is unnecessary.

Reduced kinematics based on scale separation

The assumption is made that there is a clear separation between flexible substructures responsible for the presence of the local vibrations and a stiff structure skeleton (main frame). Therefore, the subdomains Ω 1 and Ω 2 of the structure domain Ω are introduced, such that

Ω = Ω 1 ∪ Ω 2 , Ω 1 ∩ Ω 2 = ∅ , (34)
in which the subdomain Ω 1 is the stiff part and where subdomain Ω 2 is the flexible part, given by the union of all the flexible substructures. In order to completely filter the local vibrations of the flexible substructures, the indicator matrix (similarly to Eq. ( 25)) can be used for [B], that is

[B] = [χ 1 ] = I N1 0 , (35)
in which N 1 denotes the number of nodal DOFs in Ω 1 and assuming that the system matrices are partitioned accordingly. In this case, dimension r = N 1 of subspace S r is of the same order as N , that is to say it can be millions.

Exact in the stiff part and zero in the flexible part Let us consider [W] = [M] as the choice for the inner product. In this case, it can be shown that, for all u in R N , the difference E k ( u) -E r k ( u) is minimized. Using the aforementioned [B] and [W], it can be shown that the orthogonal-projection matrix [P] is given by

[P] = I N1 M -1 11 M 12 0 0 , (36)
from which it can be deduced the expression for [M r ] as,

[M r ] = M 11 M 12 M 21 M 21 M -1 11 M 12 . ( 37
)
As the number of nodal DOFs N 2 = N -N 1 in Ω 2 can be millions and since in general [M r 22 ] is not as sparse as [M 22 ], it may not be practical to build it. For a given vector u ∈ R N , product p = [M r ] u is given by

p = [M r ] u = M 11 u1 + M 12 u2 M 21 u1 + M 21 M -1 11 M 12 u2 , (38)
which requires solving a linear system of dimension N 1 . It should be noted that the orthogonal projection ur = [P] u is such that ur 2 = 0 but ur 1 = u1 . Figure 5 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. ( 36), where it can be seen that the flexible part has zero deformation and that the original deformation of the stiff part is kept unchanged everywhere except at the interface with the flexible part.

Mass elimination, or static condensation Using the matrix [B] given in Eq. ( 35) and considering the positive-semidefinite matrix, Figure 6 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. ( 40), where it can be seen that the flexible part has zero deformation and that the original deformation of the stiff part is kept unchanged. It should be noted that in spite of the discontinuity of the projection at the interface between Ω 1 and Ω 2 (as can be seen in Figure 6), the solutions ψ α of Eq. (30) satisfy the static equilibrium at the interface between Ω 1 and Ω 2 (Eq. ( 23)) and hence, continuity. The difference between Fig. 5 and Fig. 6 is due to considering a different weight matrix, namely the mass matrix [M] or the matrix defined in Eq. (39), respectively.

[W] = M 11 0 0 0 , (39)
Exact in the stiff part and flexible part in static equilibrium, or Guyan condensation The Guyan transformation matrix, [T (0)], is defined through Eqs. ( 20) and ( 21) and is such that the DOFs of Ω 1 and Ω 2 are in static equilibrium, as stated in Eq. [START_REF] Bennighof | An automated multilevel substructuring method for eigenspace computation in linear elastodynamics[END_REF]. Using the Guyan transformation matrix instead of the indicator matrix for [B] as well as using the weight matrix [W] defined in Eq. (39) leads to ur 1 = u1 and ur

2 = -[K 22 ] -1 [K 21 ]
u1 and it can be deduced that

[P] = I N1 0 -K -1 22 K 21 0 , [M r ] = M G 0 0 0 , p = [M r ] u = M G u1 0 , (41)
in which the N 1 -dimensional real vector p 1 = M G u1 can be written as

p 1 = M 11 u1 -M 12 K -1 22 K 21 u1 -K 12 K -1 22 M 21 u1 + K 12 K -1 22 M 22 K -1 22 K 21 u1 . ( 42
)
By grouping the identical terms that right-multiply [K 22 ] -1 and then the identical terms that left-multiply

[K 22 ] -1
, it can be seen that Eq. ( 42) involves solving two linear systems of dimension N 2 (with the same matrix, [K 22 ], and two right-hand sides). In the eigenvalue problem of Eq. (30), using the mass matrix [M r ] given in Eq. ( 41) is equivalent to solving the Guyan eigenvalue problem, corresponding to Eq. ( 22) with λ = 0, and then obtaining the modal displacements of Ω 2 by static equilibrium (as in Eq. ( 23)). Thereby, the present approximation is equivalent to Guyan condensation, and its implementation is such that the assembly of the dense (N 1 × N 1 ) real matrices [K G ] and [M G ] is circumvented. Figure 7 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. ( 41), where it can be seen that the deformation of the stiff part is kept unchanged and that the flexible part deforms accordingly so as to satisfy the static equilibrium with the stiff part. It should be noted that for Guyan condensation the Rayleigh-Ritz step given by Eq. ( 31) is unnecessary (the ν eigenvectors r α form the canonical basis of R ν ) because the eigenvectors ψ α belong to subspace S r (i.e. [P]ψ α = ψ α ).

Exact in the stiff part and flexible part in static equilibrium, with minimum residual kinetic energy It is recalled that for the case [W] = [M], the residual kinetic energy, E k ( u) -E r k ( u), is minimum, for all u in R N . More precisely, for all v in subspace S r , it can be shown that

E k ( u) -E k ( ur ) = E k ( u -ur ) ≤ E k ( u -v) . (43)
The expressions for matrices

t G ] = -[K 22 ] -1 [K 21 ]
, the matrices [P] and [M r ] can be defined by their blocks, such that

[P 11 ] = [M G ] -1 ([M 11 ] + [t G ] T [M 21 ]) , [P 12 ] = [M G ] -1 ([M 12 ] + [t G ] T [M 22 ]) , [P 21 ] = [t G ][P 11 ] , [P 22 ] = [t G ][P 12 ] , [M r 11 ] = ([M 11 ] + [M 12 ][t G ])[M G ] -1 ([M 11 ] + [t G ] T [M 21 ]) , [M r 12 ] = ([M 11 ] + [M 12 ][t G ])[M G ] -1 ([M 12 ] + [t G ] T [M 22 ]) , [M r 21 ] = ([M 21 ] + [M 22 ][t G ])[M G ] -1 ([M 11 ] + [t G ] T [M 21 ]) , [M r 22 ] = ([M 21 ] + [M 22 ][t G ])[M G ] -1 ([M 12 ] + [t G ] T [M 22 ]) . (44)
Since [M r ] = [M][P], the product p = [M r ] u can be deduced from the expressions of ur 1 and ur 2 , which are now given. Introducing the N 1 -dimensional real vector b 1 , such that

b 1 = [M 11 ] u1 + [M 12 ] u2 -[K 12 ][K 22 ] -1 ([M 21 ] u1 + [M 22 ] u2 ) , (45)
it can be deduced that ur

1 = [M G ] -1 b 1 and ur 2 = [t G ] ur 1 .
It can be seen that in this case the computation of the orthogonal projection ur of a N -dimensional real vector u is very expensive. Not only because of the computation of ur 2 from ur 1 and of b 1 , which each involve solving a linear system with matrix [K 22 ], but more importantly because of solving for ur 1 the linear system,

[M G ] ur 1 = b 1 . (46)
The (N 1 × N 1 ) positive-definite real matrix [M G ] is dense. In order to solve Eq. ( 46), it is advised to use a preconditioned conjugate gradient (PCG) algorithm with sparse matrix [M 11 ] as a preconditioner. Each iteration of the PCG algorithm is based on the right-product of [M G ] with a given vector of R N1 . Such product is detailed in Eq. ( 42) and followed by a cost analysis. The PCG algorithm does not require the assembly of dense matrix [M G ]. Figure 8 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. ( 44), where it can be seen that the static equilibrium between the flexible and stiff parts is satisfied, that the local displacements are reproduced, and that the deformation of the stiff part is changed accordingly so as to satisfy the static equilibrium with the flexible part. Thus, the choice [W] = [M] for the weight matrix, which implies an orthogonal projection that minimizes the residual kinetic energy, may not be suited to the filtering of local displacements. It can be shown that, although the flexible and stiff parts are in static equilibrium, the eigenvectors ψ α do not belong to subspace S r and consequently, the Rayleigh-Ritz step given by Eq. ( 31) remains necessary. The difference between Fig. 7 and Fig. 8 is due to considering a different weight matrix, namely the matrix defined in Eq. ( 39) or the mass matrix [M], respectively.

Reduced kinematics without scale separation

There are some multiscale structures exhibiting numerous local vibration modes but for which it is not possible to attribute a structure subdomain that would be the support of said local vibration modes. This can be the case for very complex structures such as cars. This paper is a continuation of previous work [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF], whose primary concern was the robust and efficient model reduction of car FE models [START_REF] Arnoux | Stochastic reduced order computational model of structures having numerous local elastic modes in low frequency dynamics[END_REF]. For extracting the dominant and meaningful global modes, the strategy was to introduce a spatial filtering of small wavelengths by using a coarse (reduced) kinematics associated with macro-elements. These macro-elements are defined by the introduction of a domain partition, such that

Ω = Ns i=1 Ω i , ∀i = j Ω i ∩ Ω j = ∅ , (47)
in which N s denotes the number of subdomains (macro-elements). An orthogonal-projection operator, h r , was then introduced [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF], such that

{h r (u)} (x) = Ns i=1 1 Ωi (x) 1 m i Ωi ρ(x )u(x )dx , (48)
in which u(x) and ρ(x) are the displacement field and the mass density at point x, quantity m i = Ωi ρ(x)dx is the mass of subdomain i, and 1 Ωi is the indicator function such that 1 Ωi (x) = 1 if x ∈ Ω i and 1 Ωi (x) = 0 otherwise. Then, the discretization [H r ] of h r was introduced to obtain the low-rank mass matrix

[M r ] = [H r ] T [M][H r
] as in Eq. ( 29), with [H r ] a special case of [P]. Finally, the global eigenvectors were also obtained using Eq. (30). Operator h r is an orthogonal projector onto the subspace of functions that are constant in each subdomain. Orthogonal-projection matrix [H r ] can be constructed [START_REF] Ezvan | Multilevel model reduction for uncertainty quantification in computational structural dynamics[END_REF] by using the matrix [B] such that

[B] =       S 1 0 • • • 0 0 S 2 . . . . . . . . . . . . . . . 0 0 • • • 0 S Ns       , ( 49
)
in which submatrix [S i ] associated with subdomain Ω i is given by

[S i ] =      T i (X αi,1 ) T i (X αi,2 ) . . . T i (X α i,N i )      , ( 50
)
in which [T i ] is a nodal transformation matrix of dimension (6 × r i ), with r i the dimension of the kinematics of subdomain Ω i (r i = 3 for a constant field in Ω i ). In Eq. ( 50), quantity α i,j denotes the index of node j of subdomain Ω i , X αi,j is its coordinate vector, and N i is the number of nodes in subdomain Ω i . In fact, for a displacement field that is constant within each subdomain, the nodal transformation matrix [T i ] is independent of space and also, it is independent of the subdomain. Consequently, for the case of a piecewise constant kinematics, the notation for the nodal transformation matrix can be simplified. For all subdomain i = 1, . . . , N s and all node j = 1, . . . , N i , it can be rewritten as

[T i (X αi,j )] = [T ], in which matrix [T ] is given by [T ] =         1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0         . ( 51
)
In Eqs. (49), (50), and (51) it has been assumed that the DOFs are sorted by subdomains, and also by nodes with the DOF ordering T X, T Y , T Z, RX, RY , RZ, with T X the translation DOF following X, RY the rotation DOF around Y , etc. In this case of a piecewise constant approximation, one has r = Ns i=1 r i = 3N s . Because of the definition of operator h r , its discretization [H r ] = [P] is obtained following the construction for [P] given in Section 4 with reduced basis [B] defined through Eqs. ( 49), (50), and (51) but also with the following choice of weight matrix, [W] = [M]. This piecewise reduced kinematics cannot represent the local fluctuations within each subdomain and consequently, the filtering of local displacements is controlled by the domain decomposition. In some cases, it can be useful to introduce enriched kinematics within each subdomain. A piecewise rigid-body motion (under infinitesimal strain assumption) is given by considering the nodal transformation matrix [T i ] such that

[T i (X α )] =         1 0 0 0 Z α -Y α 0 1 0 -Z α 0 X α 0 0 1 Y α -X α 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         , (52)
with X α = (X α , Y α , Z α ) the coordinate vector of the node α of subdomain Ω i . In this case of a piecewise rigid-body motion approximation, one has r = 6N s . A piecewise linear approximation can be obtained with the following nodal transformation matrix,

[T i (X α )] = 1[T ] X α [T ] Y α [T ] Z α [T ] , (53)
which is such that the approximation is independent for each direction. Similarly, a piecewise polynomial approximation can be obtained with the following nodal transformation matrix,

[T i (X α )] = p 1 (X α )[T ] p 2 (X α )[T ] . . . p ni (X α )[T ] , (54)
in which n i denotes the number of polynomials considered for the reduced kinematics of subdomain Ω i and where p k with k = 1, . . . , n i is a polynomial function of space. Figure 9 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by a piecewise constant approximation, where it can be seen that the local displacements are filtered and that the global deformation shape is reproduced very roughly. Figure 10 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by a piecewise rigid body approximation, where it can be seen that the local displacements are filtered and that the global deformation shape is approximately reproduced. Figure 11 depicts the deformation shape of the orthogonal projection of the multiscale mode of Fig. 4 given by a degree-5 polynomial approximation over the whole structure domain, where it can be seen that the local displacements are filtered and that the global deformation shape seems to be reproduced fairly accurately.

Figure 9: Orthogonal projection of the multiscale mode of Fig. 4 given by a piecewise constant approximation Figure 10: Orthogonal projection of the multiscale mode of Fig. 4 given by a piecewise rigid body approximation Figure 11: Orthogonal projection of the multiscale mode of Fig. 4 given by a degree-5 polynomial approximation over the whole structure domain

Combined kinematics

As already stated, in this paper we restrain ourselves to the case of multiscale structures with scale separation. Nevertheless, the ingredients of both strategies presented in Sections 4.1 and 4.2 can be useful to this case.

One drawback of the strategy presented in Section 4.2 is that the approximation of the global deformations is dependent upon the local fluctuations, in case these do not cancel out in each subdomain. In practice, the amplitude of the local deformations can be so large that they dominate the resulting averaged displacement of each subdomain. The strategy presented in Section 4.1 alleviates this problem by affecting zero displacement to the spatial support of local displacements and, in addition, it considers an exact representation of the global deformation instead of a piecewise coarse one (another possibility is, instead of affecting zero displacement to the spatial support of local displacements, to consider the spatial support of local displacements to have small or zero weight for the orthogonal projection). One drawback of the strategy presented in Section 4.1 is that it may eliminate the inertia of the flexible substructures, which can be significant. The static condensation method given by Eq. ( 40) completely neglects the inertia of the flexible substructures and consequently, the obtained global modes do not take into account the coupled resonant vibrations of the main frame and the substructures. In addition, the neglected mass of the substructures may induce a shifting of the global eigenfrequencies. The combination of the two strategies presented in Sections 4.1 and 4.2 is now investigated. In particular, we are interested in combining the static condensation method with the introduction of a reduced kinematics for the substructures, so that (i) the shifting of the global eigenfrequencies be corrected, (ii) the resonant global vibrations arisen by coupling with the resonant local vibrations be reproduced, and (iii) the local vibrations of the substructures that give rise to the numerous local elastic modes be not represented. To this end, the following reduced basis is considered, 30) and can be deduced from the expressions of ur 1 and ur 2 , which are now given. To do so, the following quantities are sequentially introduced,

[B] = I N1 0 0 b , (55)
• the (N 1 × r 2 ) real matrix [V 0 ] = [M 12 ][b], • the (N 1 × N 1 ) upper-triangular real matrix [U 1 ] from the Cholesky decomposition [M 11 ] = [U 1 ] T [U 1 ],
• the (

N 1 × r 2 ) real matrix [V 3 ] = [U 1 ] -T [V 0 ] (computed by forward substitutions), • the (r 2 × r 2 ) real matrix [M 2 ] = [b] T [M 22 ][b],
• the (

r 2 × r 2 ) real matrix [H] = [M 2 ] -[V 3 ] T [V 3 ] -1 , • the (N 1 × r 2 ) real matrix [V 1 ] = [U 1 ] -1 [V 3 ][H] (computed by backward substitutions), • the (N 2 × r 2 ) real matrix [V 2 ] = [M 22 ][b][H] -[M 21 ][V 1 ].
Then, it can be shown that, for all u in R N ,

ur 1 = u1 + v( u2 ) + [V 1 ] [b] T ([M 21 ]v( u2 ) -[M 22 ] u2 ) , ur 2 = [b] [V 2 ] T u2 , (56)
in which the real vector v of dimension N 1 is given by

v( u2 ) = [U 1 ] -1 [U 1 ] -T [M 12 ] u2 , (57)
and can be computed by forward and backward substitutions. In Eq. ( 56), the computation of v given in Eq. (57) has to be done for each u, whereas matrices [V 1 ] and [V 2 ] can be precomputed.

Simplification of the weight matrix

I N1 0 0 P 2 , [M r ] = M 11 M 12 P 2 P T 2 M 21 M 22 P 2 , ( 58
)
in which the (N 2 × N 2 ) real matrix [P 2 ] of rank r 2 is the orthogonal-projection matrix onto the column space of [b] (the reduced basis for the flexible substructures). Its expression is given by

[P 2 ] = [b] [b] T [M 22 ][b] -1 [b] T [M 22 ] . (59)
In order to speed up the computation of product p

= [M r ] u, the (N 2 × r 2 ) real matrix [b] is orthonormalized with respect to [M 22 ]. Let [U 2 ] be the (N 2 × N 2 ) upper-triangular real matrix from the Cholesky decompo- sition [M 22 ] = [U 2 ] T [U 2 ]
. The QR decomposition of the (N 2 × r 2 ) real matrix given by the product

[U 2 ][b] is introduced, such that [U 2 ][b] = [q][r] , [q] T [q] = [I r2 ] . (60)
In Eq. ( 60), [q] is a (N 2 × r 2 ) real matrix and [r] is a (r 2 × r 2 ) positive-definite upper-triangular real matrix

(provided that [b] is nonsingular). The orthonormalization [ b] of [b] satisfying [ b] T [M 22 ][ b] = [I r2 ] is obtained as [ b] = [U 2 ] -1 [q] , (61)
which can be computed by backward substitutions. As the column spaces of [b] and [ b] coincide, orthogonalprojection matrix [P 2 ] can be rewritten as

[P 2 ] = [ b][ b] T [M 22 ] . (62)
The (N 1 × r 2 ) and (N 2 × r 2 ) real matrices [Z 1 ] and [Z 2 ] are introduced, such that

[Z 1 ] = [M 12 ][ b], [Z 2 ] = [M 22 ][ b]. (63)
Then it can be deduced that, for all u in R N , product p = [M r ] u is given by

p =   M 11 u1 + [Z 1 ] [Z 2 ] T u2 [Z 2 ] [Z 1 ] T u1 + [Z 2 ] T u2   . ( 64
)
Lumped mass approximation A common practice in structural dynamics is the use of a diagonallylumped approximation for the mass matrix [M]. In general, for a fine mesh, such an approximation induces a negligible error. A diagonally-lumped mass matrix can be obtained by considering FE shape functions that are constant by element. Such an approximation is small compared to that of replacing consistent mass ] u (it is recalled that this is the building block to solve the eigenvalue problem of Eq. (30), hence the denomination online). It can be deduced that

p 2 = [Z 2 ] [Z 2 ] T u2 . (65)
Thanks to the diagonally-lumped approximation, the Cholesky decomposition

[M 22 ] = [U 2 ] T [U 2 ] is straight- forward (it is diagonal). In addition, matrix [Z 2 ] can be directly computed, without constructing [ b], as [Z 2 ] = [U 2 ] T [q]
. Moreover, thanks to the block structure of [b] as shown in Eq. ( 49) and since matrix [U 2 ] is diagonal, the QR decomposition of Eq. ( 60) can be performed blockwise. Furthermore, the sparsity inherent to [b] is preserved, such that the products in Eq. (65) are faster. In fact, introducing the notation

Ω = Ω 1 ∪ Ω 2 , Ω 1 ∩ Ω 2 = ∅ , Ω 2 = Ns-1 i=1 Ω 2,i , ∀i = j Ω 2,i ∩ Ω 2,j = ∅ , (66)
which highlights the decomposition of the flexible part (union of the flexible substructures) into several subdomains (with number N s -1, in which N s is the total number of subdomains), the product p 2 = [M 22 ] ur 2 can be written as

p 2 =      p 2,1 p 2,2 . . . p 2,Ns-1      , ( 67
)
in which, for all i = 1, . . . , N s -1,

p 2,i = [Z 2,i ] [Z 2,i ] T u2,i . (68)
In Eq. (67), p 2,i is the subvector of p 2 associated with the DOFs of subdomain Ω 2,i (and it is assumed that the DOF ordering corresponds to the domain decomposition). In Eq. (68), u2,i is the subvector of u2 associated with the DOFs of subdomain Ω 2,i and [Z 2,i ] is the matrix block of [Z 2 ] associated with the DOFs of subdomain Ω 2,i and with the r 2,i shape functions that constitute the reduced kinematics of Ω 2,i , with r 2 =

Ns-1 i=1 r 2,i . As an example, if the reduced kinematics for subdomain Ω 2,i is that of a rigid-body motion, then real matrix [Z 2,i ] is of dimension (N 2,i × r 2,i ) with r 2,i = 6 and N 2,i the number of nodal DOFs in Ω 2,i .

Compatibility with state-of-the-art eigenvalue solvers The ability of the eigenvalue problem of Eq. (30) to be solved by state-of-the-art algorithms is now discussed. The difficulty lies in the fact that low-rank mass matrix [M r ] is not as sparse as [M] can be. In the sequel, a diagonally-lumped approximation of [M] is assumed.

The automated multilevel substructuring (AMLS, see [START_REF] Bennighof | An automated multilevel substructuring method for eigenspace computation in linear elastodynamics[END_REF]) solver is very efficient in providing for approximate modes, whose computation is based on Craig-Bampton substructuring technique along with a recursive domain partitioning with several levels. The approximate modes belong to a subspace spanned by characteristic constraint modes and by components modes of the numerous substructures and, in general, the approximation is satisfactory for the purpose of FRF analysis. The substructure decomposition is based on the sparsity of the mass and stiffness matrices. The substructures are connected to each other solely by the interface DOFs. The interior DOFs of each substructure are uncoupled with respect to the DOFs of the other substructures. The component modes are the eigenmodes of said component (substructure) with fixed boundary condition at its interface. When dealing with Eq. (30), the limited sparsity of matrix [M r ] is a drawback to Craig-Bampton substructuring. With respect to the domain decomposition introduced in Eq. (66) and devoted to the filtering of local displacements, in general the substructure decomposition cannot further divide the subdomains Ω 2,i , as all the DOFs of Ω 2,i may be coupled with each other (not only the neighboring DOFs associated with the FE node connectivity may be coupled with each other, but all the DOFs of Ω 2,i ). As a conclusion, the speedup attainable by AMLS solver is limited, in this case, to that attainable with the domain decomposition that is devoted to the filtering of the local modes.

The multiple shift-invert Lanczos (MSIL, see [START_REF] Aktulga | Parallel eigenvalue calculation based on multiple shift-invert Lanczos and contour integral based spectral projection method[END_REF]) algorithm is very efficient in solving Eq. ( 2) for many eigenpairs (λ, ϕ). It is based on a division of the frequency band (spectrum slicing) and the eigenpairs are computed independently for each band. For each band, based on the introduction of an eigenvalue shift µ, the eigenpairs closest to shift µ are computed by Lanczos algorithm. The elementary operation to solve Eq. ( 2) without a shift is to compute p = [K]

-1 [M] u for many right-hand sides u ∈ R N . This can be done efficiently (through forward and backward substitutions) by introducing the Cholesky decomposition

[K] = [U ] T [U ] such that p = [U ] -1 [U ] -T [M] u. Or, alternatively, by introducing the LDL decomposition [K] = [U] T [D][U] such that p = [U] -1 [D] -1 [U] -T [M] u, in which [D]
is a block diagonal matrix with block size 1 or 2 and [U] is the right-product of an upper-triangular matrix with a permutation matrix. The elementary operation to solve Eq. ( 2) with a shift µ is to compute p = ([K] -µ[M])

-1 [M] u for many right-hand sides u ∈ R N . This can be done efficiently (through forward and backward substitutions) by introducing the LDL decomposition

[K] -µ[M] = [U(µ)] T [D(µ)][U(µ)] such that p = [U(µ)] -1 [D(µ)] -1 [U(µ)] -T [M] u, in which [D(µ)
] is a block diagonal matrix with block size 1 or 2 and [U(µ)] is the right-product of an upper-triangular matrix with a permutation matrix. In spite of the rather dense structure of low-rank mass matrix [M r ], an efficient method is now presented to compute p = ([K] -µ[M r ])

-1 [M r ] u for many right-hand sides u ∈ R N , which is the building block for solving Eq. ( 30). The computation of product p = [M r ] u has been detailed previously in this section (see for instance Eqs. ( 65) and ( 68)). Therefore, the computation of p

= ([K] -µ[M r ]) -1 u for a given u ∈ R N is now described. Introducing the sparse matrix [K µ ] such that [K µ ] = [K] -µ M 11 0 0 0 , (69)
as well as the (N × r 2 ) real matrix [Z] with expression

[Z] = 0 Z 2 , ( 70
)
in which matrix [Z 2 ] is defined in Eq. (63), expression p = ([K] -µ[M r ]) -1 u can be rewritten as

p = [K µ ] -µ[Z][Z] T -1 u . (71)
Equation ( 71) can be seen as involving the inverse of a large sparse matrix with a low-rank perturbation (with rank r 2 N ). Thanks to Woodbury matrix identity, it can be rewritten as

p = [K µ ] -1 u -[K µ ] -1 [Z] - 1 µ [I r2 ] + [Z] T [K µ ] -1 [Z] -1 [Z] T [K µ ] -1 u . ( 72
)
The following quantities are sequentially introduced,

• the (N × N ) block diagonal real matrix with block size 1 or 2, [D µ ], as well as the (N × N ) real matrix [U µ ], right-product of an upper-triangular matrix with a permutation matrix; these matrices are obtained by LDL decomposition of

[K µ ] such that [K µ ] = [U µ ] T [D µ ][U µ ],
• the (N × r 2 ) real matrix [S] = [U µ ] -T [Z] (computed by forward substitutions),

• the (r 2 × r 2 ) real matrix [H µ ] = 1 µ [I r2 ] -[S] T [D µ ] -1 [S] -1 , • the (N × r 2 ) real matrix [Y ] = [U µ ] -1 [D µ ] -1 [S][H µ ] (computed by backward substitutions).
Then it can be shown that, for all u in R N , product p

= ([K] -µ[M r ]) -1 u is given by p = w( u) + [Y ] [Z] T w( u) , (73)
in which the real vector w of dimension N is given by

w( u) = [U µ ] -1 [D µ ] -1 [U µ ] -T u , (74)
and can be computed by forward and backward substitutions. In Eq. ( 73), the computation of w given in Eq. (74) has to be done for each u, whereas matrices [Y ] and [Z] can be precomputed. As a conclusion, it is still very efficient to use MSIL for solving Eq. ( 30), at the expense of additional matrix-vector multiplications involving the thin matrices [Y ], [Z], and [Z 2,i ] for i = 1, . . . , N s -1. However, the filtering of the local modes thanks to having considered [M r ] instead of [M] leads to a decrease in the computational cost, since the number of eigenpairs to be calculated is reduced.

APPLICATIONS

This section applies the above described methods to three structures. The first one is a simple plate structure (also shown earlier) that is composed of two well-identified scales, the second one is a simplified model of a nuclear fuel assembly, and the third one is a detailed model of a nuclear fuel assembly. The purpose of the first application -the simple plate -is to assess the behavior of the proposed methods in a controlled setting. Four subsections are devoted to this simple structure. The first one deals with the exact same heterogeneous plate as in [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF] but considers a broader frequency band, and for this case, the effectiveness of the static condensation introduced in Section 4.1 is assessed. The second one considers the same frequency band as in [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF], and for this case, the effectiveness of a few kinematics introduced in Section 4.2 is assessed.

In subsections three and four, the effectiveness of static condensation is assessed for the case of an increased mass for the subscale (the flexible part). For the fuel assembly applications, the effectiveness of a combined kinematics, as introduced in Section 4.3, is assessed.

Heterogeneous plate with two structural scales

Analysis in a broad frequency range

In this first application, we deal with a simple multiscale structure, namely the heterogeneous plate depicted in Fig. 1. The structure is the same as in [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF], in which all the details and a descriptive drawing are given. The details are repeated here for completeness. The plate dimensions L x and L y along directions x and y are L x = 26 m and L y = 20 m, respectively. There are 12 panels (that are referred to as "flexible panels") supported by a main frame (that is referred to as "stiff skeleton"). Each flexible panel is an isotropic and homogeneous thin square plate with sides of 4 m and with a constant thickness h 1 = 0.002 m. The stiff skeleton is constituted of isotropic and homogeneous thin plates. For the stiff skeleton, the thickness is h = 0.17 m, the Young modulus is E = 2.1 × 10 12 N/m 2 , the Poisson ratio is ν = 0.29, and the mass density is ρ = 9800 kg/m 3 . For the flexible panels, the Poisson ratio is ν = 0.29 and the mass density is ρ 1 = 7800 kg/m 3 . Going first from the left to the right and then from the bottom to the top, the Young moduli of the flexible panels are 1.31 A domain decomposition is performed as in Eq. (34), with Ω 1 the stiff skeleton (the stiff part) and Ω 2 the 12 flexible panels (the flexible part). Each subdomain corresponds to a set of nodes (and to the corresponding set of DOFs), and the subdomains do not share any nodes. The interface nodes between the stiff skeleton and the flexible panels are considered to belong to the flexible part. There are N 1 = 30,723 DOFs in the stiff part and N 2 = 22,500 DOFs in the flexible part. We are interested in calculating the dynamic response of the stiff part under excitation of said stiff part, in the frequency band B ω = 2π×]0, 100] rad/s. Consequently, a subdomain Ω G is introduced (as in Section 3.1), defined as the restriction of the frame Ω 1 to its central part of width 0.5 m, so that its dynamic response be less sensitive to the panels local dynamics. A subset S G of nodal DOFs belonging to Ω G is considered, for which the constructed ROMs are intended to be valid. Thus, the constructed ROMs are not intended to accurately predict the dynamic response of the flexible panels. Subset S G is constituted of 22 TZ DOFs, among which 20 are associated with nodes randomly selected in Ω G , and where the remaining 2 are associated with the nodes P 1 : (x 1 = 10 m, y 1 = 7 m) and P 2 : (x 2 = 19 m, y 2 = 7 m). The random selection is based on introducing a homogeneous 5 × 4 domain partition; for each subdomain, one node is randomly selected (with equiprobability) among the nodes of the restriction to Ω G of said subdomain.

× k 1 N/m 2 , 1.47 × k 1 N/m 2 , 1.54 × k 1 N/m 2 , 1.74 × k 1 N/m 2 , 1.47 × k 1 N/m 2 , 1.71 × k 1 N/m 2 , 0.49 × k 1 N/m 2 , 1.09 × k 1 N/m 2 , 1.74 × k 1 N/m 2 , 1.34 × k 1 N/m 2 , 1.68 × k 1 N/m 2 ,
In the sequel, one particular FRF is plotted for comparing the several ROMs (the same as in [START_REF] Soize | Stochastic reduced-order model in low-frequency dynamics in presence of numerous local elastic modes[END_REF]). It is referred to as the control FRF, and is given by the modulus in dB of the accelerance along the vertical direction of point P 2 under vertical point excitation at P 1 , for all ω in B ω . Each ROM will also be characterized by its maximum error among all the combinations of output observation i = 1, . . . , 22 and input excitation j = 1, . . . , 22 as described in Section 3 (thanks to symmetry of the FRF matrix, there are only 253 different combinations). There are 961 elastic modes in frequency band B ω and 1156 elastic modes below 120 Hz. The reference model is obtained by replacing the FE model (Eq. ( 1)) by the classic ROM given by Eqs. ( 6) and [START_REF] Langley | A hybrid method for the vibration analysis of complex structural-acoustic systems[END_REF] with n e = 1156. Five different ROMs are compared. First, the classic ROM of the first n e elastic modes. Second, the ROM of the n p elastic modes with largest participation factor along the vertical direction (see Eq. ( 8)). Third, the ROM of the n d dominant elastic modes ϕ α with highest importance α as defined in the algorithm for the simplified mode sorting procedure of Section 3.3. Fourth, the dominant ROM of n D elastic modes given by the tractable mode sorting procedure of Section 3.2. Fifth, as introduced in Section 4.1 and specified by Eqs. (39) and (40), the global ROM obtained by static condensation with all the DOFs of Ω 1 (the stiff part) as the master DOFs. For each ROM, the FRFs are computed thanks to diagonal matrix equations, similarly to Eq. ( 7). For each ROM, a modal damping model is used, for which each damping ratio is given by ξ = 0.01. The very same damping modeling is used throughout the rest of this paper. Figure 12 depicts the mode importance as a function of the dominant mode rank. It is recalled that the mode importance α is the error (maximum error among the 253 FRFs) resulting from removing the elastic mode ϕ α from the reference model of the n e = 1156 elastic modes. Said error is defined through Eq. ( 16) for which the meaning is the error level in dB. For instance, it can be seen in Fig 12 that removing one dominant mode whose rank is greater than 100 induces an error level less than 1 dB. Figure 13 depicts the global ROM dimension n g as a function of the truncation parameter ν used for its construction (see Eq. ( 31)), given that n g be deduced from the cutoff frequency f cut = 120 Hz such that f g ng ≤ f cut . For each ν, dimension n g is efficiently calculated by Sturm sequence check, based on the LDL decomposition of the matrix [Ψ(ν)]

T [K][Ψ(ν)] -µ cut [Ψ(ν)] T [M][Ψ(ν)] with µ cut = (2πf cut ) 2 .
It can be seen that function ν → n g (ν) reaches a plateau starting from a given point (ν = 70), for which n g = 56. The global ROM of dimension n g = 56 is constructed with ν = 301. For this ROM, the error (n g ) = max ij (u ij , u (ng) ij ) with n g = 56 is given by (n g ) = 2.7 dB (where u (ng) ij is the FRF in dB scale for output observation i and input excitation j calculated by the global ROM of dimension n g ).

For comparison, the error (n D ) associated with the dominant ROM of n D elastic modes given by the tractable mode sorting procedure is, for n D = 56, given by (n D ) = 3.9 dB. The error (n d ) associated with the ROM of the n d dominant elastic modes given by the simplified mode sorting procedure is, for n d = 56, given by (n d ) = 5.1 dB. The error (n p ) associated with the ROM of the n p elastic modes with largest participation factor is, for n p = 56, given by (n p ) = 24.8 dB. Finally, the error (n e ) associated with the ROM of the n e elastic modes with lowest eigenfrequency is, for n e = 56, given by (n e ) = 59.0 dB. Figure 14 Several observations can be made based on Fig. 14. First, basis truncation solely based on eigenfrequency gives an error that slowly decreases with a roughly constant rate. Second, a mode selection based on the mass participation factors leads to large errors, which means that some important modes have a low mass participation factor. Finally, the proposed methods allow to obtain accurate ROMs for which the dimension is significantly reduced. For this application, the global ROM based on static condensation is better than any of the other ROMs considered, as it is Pareto dominant. It can also be seen that the dominant ROM of the tractable mode sorting procedure is Pareto dominant with respect to the ROM of the dominant elastic modes of the simplified mode sorting procedure. The extent of the error introduced by the proposed smalldimension ROMs is illustrated in Figs. 15 and16, which depict the resulting approximated control FRF in comparison to the exact one computed with the ROM of the n e = 1156 elastic modes. It can be seen in Fig. 16 that an error of 1 dB is quite small. Figures 15 and16 show that it is possible to capture the structure's global dynamics with a reduced subspace. An advantage of the global ROM is that it does not require the computation of the n e = 1156 elastic modes. For this simple structure, the computation of the global ROB took a few seconds whereas the computation of the elastic modes took a few minutes. It should be noted that, given the very small thickness h 1 = 0.002 m of the flexible panels with respect to the thickness h = 0.17 m of the stiff skeleton, the total mass of the flexible part is very small in comparison to that of the stiff part. It can be verified that the flexible part represents roughly 0.55% of the total mass. In such a case, it seems reasonable that a static condensation of the flexible part, obtained by mass elimination of said flexible part, lead to a small error. It has been verified that all the global modes ϕ g α effectively have a global deformation shape. Thus, the contribution of the local displacements, which could be defined as belonging to the orthogonal complement of the global ROB, is very small for this structure.

It should be noted that the use of the Guyan condensation instead of the static condensation for constructing the global ROM led to nearly identical results, which are not shown for brevity. The same can be said for the use of the global ROM for which the modified mass matrix is defined in Eq. (37) (and for which the results are not presented for brevity). Despite the fact that Guyan condensation might seem more suitable than static condensation, especially when comparing Fig. 6 to Fig. 7, the latter has been chosen for several reasons. First, it is more efficient than Guyan condensation, as can be seen when analyzing Eq. (42). Second, although for Guyan condensation the kinetic energy of the flexible part is nonzero, said kinetic energy is associated with the static deformation of the flexible part, which can be very different from its exact dynamic deformation. Rather than improving the results, introducing such a nonzero kinetic energy for the flexible part is likely to let some local modes be reproduced (as has been observed by the authors when applying the methodology to other structures). To sum up, static condensation is preferred, so as to efficiently and effectively obtain a ROB that is devoid of local modes. This ROB might not be the best basis, especially in the case of strong coupling between the global displacements and the local displacements. As stated in Section 4.3, the static condensation does not take into account the coupled resonant vibrations of the stiff part with the flexible part and, in addition, the neglected mass of the flexible part may induce a shifting of the global eigenfrequencies. In Fig. 16, a small shift of the resonance peaks can be observed. Moreover, it can be observed that the amplitude of the said resonance peaks is overestimated. This apparent loss of damping is a phenomenon described by the fuzzy structure theory [START_REF] Soize | A model and numerical method in the medium frequency range for vibroacoustic predictions using the theory of structural fuzzy[END_REF], and is related to the fact that the mechanical energy absorbed by the flexible part is not represented by the global ROB (and consequently, in the global ROM).

Analysis in a low-frequency range

In this section, the same analysis is performed, but by considering the LF band B ω = 2π×]0, 11] rad/s as frequency band of analysis. The purpose of this analysis is to investigate other possible kinematics for the construction of the global ROB. In particular, attention will be given to reduced kinematics without scale separation. Similarly to the previous section, the reference model is obtained by using the first n e = 1156 elastic modes. Fig. 17 depicts the mode importance as a function of the dominant mode rank. It is seen that the dominant subspace is concentrated in an even smaller number of modes, which is expected since the frequency band is narrower. Similarly to previous section, the dimension n g of the global ROM by static condensation of the flexible part reaches a plateau, starting from ν = 8, and for which n g = 8. Said global ROM is constructed with ν = 301 and its error is (n g ) = 10.9 dB. In comparison, its error for the control FRF is only 1.8 dB. The big error (maximum error among the 253 FRFs) is due to some other FRF for which the antiresonances are not converged (the convergence of the static contributions is very slow). The error (n D ) associated with the dominant ROM of n D elastic modes given by the tractable mode sorting procedure is, for n D = 8, given by (n D ) = 12.0 dB (1.8 dB for the control FRF). The error (n d ) associated with the ROM of the n d dominant elastic modes given by the simplified mode sorting procedure is, for n d = 8, given by (n d ) = 13.7 dB (2.8 dB for the control FRF). The error (n p ) associated with the ROM of the n p elastic modes with largest participation factor is, for n p = 8, given by (n p ) = 22.9 dB (12.5 dB for the control FRF). Finally, the error (n e ) associated with the ROM of the n e elastic modes with lowest eigenfrequency is, for n e = 8, given by (n e ) = 38.1 dB (19.0 dB for the control FRF).

Figure 18 depicts the error function (n) with respect to the ROM dimension n, for each of these ROMs. Similarly to previous section, it can be seen that the error of the global ROM by static condensation is on the trajectory of the error functions of the ROMs by automatic mode selection. It can also be seen that the mass participation criterion is not satisfactory at all, which indicates that some preponderant elastic modes are discarded. Compared to previous section, the convergence toward a small error (less than a few decibels) is slow and is due to the slow convergence of the static contributions. There are 74 elastic modes in B ω and 100 elastic modes below 13.2 Hz. The classic ROM of the first n e = 100 elastic modes has an error (n e ) = 10.9 dB (1.7 dB for the control FRF). Similarly to Fig. 13, Figure 19 depicts the function ν → n g (ν), but this time for several global ROMs and with n g deduced from the cutoff frequency f cut = 13.2 Hz for which f g ng ≤ f cut . The first global ROM is based on a piecewise constant kinematics (see Eq. ( 51)), corresponding to a homogeneous 5 × 4 domain decomposition that can be deduced from Fig. 9. The second global ROM is based on a piecewise rigid body kinematics (see Eq. ( 52)), with the same domain decomposition. The third global ROM is not based on any domain decomposition, as the reduced kinematics is based on polynomials (see Eq. ( 54)) with the whole structure as support. The fourth global ROM is constructed with the same reduced kinematics as the third global ROM, except that the weight matrix [W] is not [W] = [M] but is given by Eq. (39); indices 1 and 2 are associated with the same subdomains Ω 1 and Ω 2 as in previous section, which means that this fourth global ROM exploits the scale separation. This weight matrix [W] is such that the associated orthogonal projection (through matrix [P]) is independent of the flexible part deformation. Hence, it will be referred to as "discriminatory weight matrix". The polynomial basis is based on the numerical orthogonalization of a monomial basis. All the monomials of x and y with degree less or equal to 5 are included, which yields r = 21 as dimension of the reduced kinematics. The orthogonalization is mandatory because the polynomials rapidly become numerically linearly dependent. The implementation of the global ROM for the fourth case is not straightforward, as the inner product associated with [W] is degenerate. For the polynomial reduced kinematics, the columns of matrix [B] are the r monomials.

Let [ B] denote the orthonormalization of [B] such that [ B] T [W][ B] = [I r ]. ( 75
)
Similarly to Eqs. ( 60) and (61), matrix [ B] of the fourth global ROM can be obtained based on the QR decomposition such that

U 1 0 0 0 B 1 B 2 = [Q][R], ( 76
)
in which the (N 1 × N 1 ) positive-definite real sparse matrix [U 1 ] is obtained from the Cholesky decomposition

[M 11 ] = [U 1 ] T [U 1 ], the (N × r) real matrix [Q] is such that [Q] T [Q] = [I r ],
and where [R] is a (r × r) positive-definite upper-triangular real matrix. From Eq. (76) it can be deduced that the (N 2 × r) submatrix [Q 2 ] of [Q] associated with the DOFs of Ω 2 is zero. Consequently, the QR decomposition of Eq. ( 76) can be obtained by performing the following QR decomposition,

[U 1 ][B 1 ] = [Q 1 ][R 1 ], (77)
in which [Q 1 ] coincides with the (N 1 × r) submatrix of [Q] associated with the DOFs of Ω 1 and where

[R 1 ] = [R]. From Eq. (76) it can be deduced that [Q] T [Q] = [R] -T [B] T [W][B][R] -1 . Hence matrix [ B] can be obtained as [ B] = [B][R 1 ] -1 .
It should be noted that this construction differs from that described in Section 4.3. It can be observed in Fig. 19 that the polynomial kinematics with discriminatory weight matrix gives a number of modes n g equal to that obtained with static condensation of the flexible part, and such that n g converge with respect to truncation ν. In contrast, the other kinematics, which are not based on scale separation, do not exhibit the same behavior. The dimension n g is strongly dependent on ν, which indicates that the addition of higher-frequency eigenvectors ψ leads to the reconstitution of local modes with low eigenfrequency. This phenomenon is put into evidence through Figs. [START_REF] Xia | A new iterative order reduction (IOR) method for eigensolutions of large structures[END_REF] to 23, which depict the functions α → √ σ α /2π and α → f g α (ν) for all ν = 1, . . . , r for each global ROM and in comparison with the same functions associated with the global ROM by static condensation. Figure 20 It can be observed in Fig. 20 that the relative difference between σ α and λ α is very important for the global ROM by piecewise constant kinematics, whereas this is not the case for the global ROM by static condensation. Moreover, it can be seen that for the global ROM by piecewise constant kinematics the eigenfrequencies are sensitive to truncation parameter ν. For the piecewise rigid body kinematics, the same remarks can be done (see Fig. 21); in addition, for most truncations ν, the modal density is low in the neighborhood of the cutoff frequency (there is a jump of function α → f g α (ν)), which is symptomatic of the fact that the mixing of HF eigenvectors ψ leads to the reconstitution of LF local modes ϕ g (such a jump means that the modal density is higher in low frequencies than in high frequencies). For the polynomial kinematics (see Fig. 22), the eigenfrequencies are also quite sensitive to truncation parameter ν, which puts into evidence the apparition of local modes induced by the addition of HF eigenvectors ψ. However, an interesting result can be observed, namely that the eigenvalues σ α are very close to the eigenvalues λ g α obtained by static condensation. Nevertheless, the orthogonalization by Rayleigh-Ritz analysis reveals local modes in [Φ g ], for which the number increases with ν. For the polynomial kinematics with discriminatory weight (see Fig. 23), it can be observed that the global eigenfrequencies f g α coincide very closely with those of static condensation. Second, for this kinematics, the eigenvalues σ α and λ g α are very close, far more than it is the case for static condensation. To interpret this, it is recalled that Guyan and static condensations yield very close FRFs (hence very close eigenvalues λ g α although not necessarily very close eigenvalues σ α ), and that no Rayleigh-Ritz step is necessary for Guyan condensation since λ g α = σ α for this case (as explained in Section 4.1). It can be deduced that the orthogonal projection given by the polynomial kinematics and defined with the discriminatory weight matrix [W] is a very good approximation to the orthogonal projection associated with Guyan condensation. It should be noted that the deformation shapes in Figs. ted lines) and α → f g α (ν) for ν = 1, . . . , r (solid lines) for static condensation (in gray) and for polynomial kinematics with discriminatory weight matrix (in black).

Because of the important sensitivity of these global ROMs with respect to truncation parameter ν, their error (n g ) is a function of ν (with n g deduced from frequency truncation). Figure 24 depicts the error function ν → (n g (ν)), for the four global ROMs. It can be seen that except for the polynomial kinematics with discriminatory weight matrix, the error is bigger and also, it keeps fluctuating with the increase of ν. As expected, the errors, given by the global ROMs of static condensation on one hand and of polynomial kinematics with discriminatory weight matrix on the other hand, coincide. The global ROM with piecewise constant kinematics is constructed with ν = 20 and n g = 14 (n g is deduced from frequency truncation with respect to f cut = 13.2 Hz) and exhibits an error of (n g ) = 16.7 dB (6.7 dB for the control FRF). The global ROM with piecewise rigid body kinematics is constructed with ν = 60 and n g = 44 and exhibits an error of (n g ) = 12.8 dB (2.4 dB for the control FRF). The global ROM with degree-5 polynomial kinematics is constructed with ν = 21 and n g = 16 and exhibits an error of (n g ) = 17.4 dB (4.9 dB for the control FRF).

The global ROM with degree-5 polynomial kinematics with discriminatory weight is constructed with ν = 21 and n g = 8 and exhibits an error of (n g ) = 10.9 dB (1.9 dB for the control FRF). It is recalled that the ROM error (maximum error among the 253 FRFs) is perturbed by the slow convergence of the static contributions.

For each of these four global ROMs, the obtained control FRF is plotted in comparison to that of the reference of the first n e = 1156 elastic modes (see Figs. 27 through 30). In addition, Figure 25 depicts the control FRF calculated by the dominant ROM of n D = 8 elastic modes selected by the tractable procedure, and Figure 26 depicts the control FRF calculated by the global ROM of static condensation and for which n g = 8. Large errors in the predictions of the control FRF can be observed in Figs. 27 and 29, in particular with the presence of spurious resonance peaks, in conjunction with the inability to predict the clustered peaks in the band 2π × [8, 9] rad/s. The global ROM based on rigid body kinematics being composed of a large number of modes (n g = 44) as constructed with ν = 60, it is able to roughly represent the control FRF but also presents numerous small spurious resonance peaks as well as significant errors (see Fig. 28). These anomalies can be explained by the fact that spurious global deformation shapes can be obtained by least-square fit of deformation shapes dominated by local deformations (the orthogonal projection given by matrix [P] comes down to a least-square approximation). It is observed that the global ROM based on the discriminatory weight matrix (for defining [P]) is very satisfactory and devoid of spurious resonances (see Fig. 30). Frequency ( Frequency (

-70

Accelerance (dB) exact weighted polynomial

Figure 30: Control FRF calculated with the global ROM of dimension n g = 8 obtained by polynomial kinematics with discriminatory weight matrix and for which = 1.9 dB

Heterogeneous plate with increased panel mass

The analysis presented in Section 5.1.1 showed that the global ROM based on static condensation of the flexible panels gave a very good accuracy for a low dimension. The validity of the global ROM by static condensation is now assessed for the case of a flexible part with increased mass. For this, the material properties of the flexible panels are changed, all other things being equal. The material properties are changed such that the flexible panels represent about 5.4% of the total mass (instead of 0.55%). The eigenfrequencies of a rectangular thin plate with simply supported edges are proportional [START_REF] Craig | Fundamentals of Structural Dynamics[END_REF] to h E/ρ with h the thickness, E the Young modulus, and ρ the mass density. In order to keep similar eigenfrequencies for the flexible panels and consequently, similar modal density, parameters h 1 , k 1 , and ρ 1 are replaced by h 2 = 0.0167 m, k 2 = 18 × 10 9 , and ρ 2 = 9800 kg/m 3 . For this structure, there are 961 elastic modes in frequency band B ω and 1157 elastic modes below 120 Hz, which is very similar to the previous case. Figure 31 depicts the mode importance as a function of the dominant mode rank, in comparison to the same curve that was obtained for previous structure.

It can be seen that the dominant modes with rank around 100 have an importance roughly twice as much as their counterpart of the previous structure with low panel mass; the dominant contributions are less concentrated. Similarly to Section 5.1.1, Figure 32 depicts the function ν → n g (ν) for the global ROM given by static condensation of the flexible part. Dimension n g reaches a plateau, for which n g = 60. The global ROM of dimension n g = 60 is constructed with ν = 301. For this ROM, the error is (n g ) = 9.6 dB. For comparison, the error (n D ) associated with the dominant ROM of n D elastic modes given by the tractable mode sorting procedure is, for n D = 60, given by (n D ) = 7.1 dB. The extent of the error introduced by the proposed small-dimension ROMs is illustrated in Figs. 34 and35, which depict the resulting approximated control FRF in comparison to the exact one computed with the ROM of the n e = 1157 elastic modes. It can be seen that it is still possible to capture the structure global behavior through a reduced subspace. Nevertheless, the approximation is far less satisfactory than for the previous structure with low panel mass. It can also be observed that the FRF presents more fluctuations and peaks. Such a richer dynamics would require more dominant modes (increased n D ) to attain the same accuracy as before. For this structure, the contribution of the local displacements is more important than previously. As for the previous structure but in an amplified manner, it can be observed in Fig. 35 that the prediction of the global ROM is tarnished by a frequency shift that increases with frequency as well as by resonances with overestimated amplitudes. In contrast, the major flaw of the dominant ROM by automatic selection lies rather in the fact that it underestimates the dynamic response. An advantage of the ROM by mode selection is that the error can be decreased at will by the addition of more modes, whereas this is not the case for the global ROM by static condensation. In this last analysis conducted on the heterogeneous plate, the properties of the flexible panels are changed once again, this time such that the total mass of the flexible part represent about 54.4%. In order to keep similar eigenfrequencies for the flexible panels and consequently, similar modal density, parameters h 1 , k 1 , and ρ 1 are replaced by h 3 = 0.12 m, k 3 = 1 × 10 9 , and ρ 3 = 28,300 kg/m 3 . For this structure, there are 1010 elastic modes in frequency band B ω and 1212 elastic modes below 120 Hz, which is quite similar to the previous cases. The panels are practically as thick as the main frame, are very heavy and very soft (very high mass density and very low Young moduli). It should be noted that for this structure the eigenvalues σ α and λ g α differ largely (there can be a more than 50% difference). It is recalled that for Guyan condensation the equality σ α = λ g α holds. It has been verified that for this structure the results by Guyan and static condensations are very close, in particular their respective global eigenfrequencies f g α . Consequently, in comparison to the ROM by polynomial kinematics analyzed in Section 5.1.2, the Rayleigh-Ritz step corrects the eigenvalues, instead of the contrary. It can be seen that for both ROMs the error is very significant. The FRF exhibits an even more enriched dynamic content (compared to the previous structures with lower panel mass), which fails to be represented by the small-dimension ROMs. The following conclusions can be drawn regarding this simple example of a heterogeneous plate. First, the complexity of the dynamic response increases with the mass of the flexible panels. Second, the contribution of the local modes increases with the mass of the flexible panels.

Finally, for a fixed accuracy, the number of dominant modes to be retained increases with the mass of the flexible panels. It should be noted that there is no guarantee that the tractable mode sorting procedure (defined in Section 3.2) give the optimal mode selection. In addition, the optimal ROM of elastic modes may not be the optimal ROM (for instance, this was not the case for the application of Section 5.1.1). Therefore, it is unsure whether there exists a small-dimension subspace that would be able to capture the main dynamic content of the global response of this structure, let alone the fast fluctuations of said response with respect to the frequency.

Simplified nuclear fuel assembly

The motivation of the present research is the vibration analysis of nuclear fuel assemblies, for the safety analysis and damage identification of spent nuclear fuel packages. A model reduction is proposed, in the spirit of that performed in [START_REF] Batou | Reduced-order computational model in nonlinear structural dynamics for structures having numerous local elastic modes in the low-frequency range. Application to fuel assemblies[END_REF] for PWR fuel assemblies.

In this section, a simplified nuclear fuel assembly is considered. It is composed of a 10 × 10 regular array of fuel rods held together by ten regularly-spaced tie plates (the fuel rods are clamped to the plates). The fuel rods are 4,124.6 mm long and with varying diameter that is mainly of 10. For this structure, it can be considered that the stiff part is the ten tie plates. In contrast to the heterogeneous plate of Section 5.1, this simplified fuel assembly does not have any connected stiff skeleton that would be the support of the flexible part. Rather, the 100 fuel rods together constitute the skeleton. This beamlike structure still exhibits global modes, that is, beam-like modes. Each fuel rod is divided into 9 parts, separated by the tie plates. Each of the 900 parts is the support of local vibrations. Therefore, the fuel rods are the support of simultaneously the global deformations and the local deformations. Since the fuel rods are all identical, the local deformations of the 900 parts occur at the same frequency. A local elastic mode is thus constituted of a concurrent deformation of all the fuel rod parts, with the fuel assembly skeleton remaining still, as is the case for the local mode plotted in Fig. 41. For both of them, the same point load is applied to the edge of one intermediate tie plate in the x direction (z is the longitudinal direction), for all ω in B ω . The load is not applied at the center of the edge, so that there be torsion of the fuel assembly. The control FRF 1 is the accelerance following x at one opposite edge of another intermediate spacer grid (not in the center of the edge). The control FRF 2 is the accelerance at the same location but following direction y. We are interested in predicting the global response of the fuel assembly. Therefore, subdomain Ω G is defined as the union of the ten tie plates. A subset S G of 240 DOFs is selected, based on a random selection of 80 nodes belonging to Ω G -for each node, the 3 translations are considered. Similarly to Section 5.1, the random selection is based on the introduction of a homogeneous domain partition. For each subdomain, one node is randomly selected (with equiprobability) among those that belong to its restriction to Ω G .

The simplified mode sorting procedure of Section 3.3 is used to construct a ROM of dominant elastic modes. Concerning the construction of the global ROM, three cases are investigated. For all of them, a domain decomposition is introduced. The fuel assembly is cut into 100 transverse slices of equal thickness (similarly to the strategy employed in [START_REF] Batou | Reduced-order computational model in nonlinear structural dynamics for structures having numerous local elastic modes in the low-frequency range. Application to fuel assemblies[END_REF]). The first global ROM is constructed by considering a piecewise constant kinematics for the kinetic energy. The second one is constructed by considering a piecewise rigid body kinematics. The third global ROM is based on a combined kinematics, as introduced in Section 4.3. The combined kinematics is based on the decomposition into 100 slices, but restricted to the fuel rods; the kinematics of the tie plates is kept exact. The reduced kinematics for the fuel rods is that of a rigid body motion by slice. It should be noted that the corresponding ROM by static condensation (zero displacement for the fuel rods and discriminatory weight matrix) did not give satisfactory results, which can be explained by the fact that the fuel rods represent about 94.8% of the structure total mass. The global ROM by piecewise constant kinematics is constructed with a truncation ν = 92 such that √ σ ν /2π ≤ 1500 < √ σ ν+1 /2π, which yields a dimension n g = 76 deduced by frequency truncation with respect to f cut = 1200 Hz. The error of this global ROM in predicting the control FRF 1 is = 5.6 dB (see Fig. 42). The overall prediction is satisfactory, despite missing a few secondary resonance peaks in the band ]0, 100] Hz and the noticeable discrepancy in the band ]800, 1000] Hz. It should be noted that similarly to the plate example of Section 5.1, the prediction is characterized by a small frequency shift to the left as well as by slightly overestimated resonances. The error of this global ROM in predicting the control FRF 2 is = 51.0 dB (see Fig. 46). This global ROM is unable to predict the control FRF 2. This is due to the fact that the piecewise constant kinematics is unable to represent fuel assembly torsion. The global ROM by piecewise rigid body kinematics is constructed with a truncation ν = 120 such that √ σ ν /2π ≤ 1500 < √ σ ν+1 /2π, which yields a dimension n g = 113 deduced by frequency truncation with respect to f cut = 1200 Hz. The error of this global ROM in predicting the control FRF 1 is = 4.1 dB (see Fig. 43). The prediction is nearly perfect in the frequency band ]0, 400] Hz, it is quite satisfactory in the band ]400, 800] Hz although tarnished by a small shift, whereas the prediction in the band ]800, 1000] Hz is not good as tarnished by a large frequency shift. Contrary to the global ROM by piecewise constant kinematics, it can be seen in Fig. 46 that the prediction of control FRF 2 remains as satisfactory as that of control FRF 1, with an error = 3.7 dB. Then, for the combined reduced kinematics, it can be seen in Figs. 44 and 48 that the prediction is similar to that of piecewise rigid body kinematics. Therefore, the remaining discrepancy (in particular, in the frequency band ]800, 1000] Hz) is due to the reduced kinematics of the fuel rods. Finally, the prediction given by the ROM of the 113 dominant elastic modes is shown in Figs. 45 and 49 with = 1.4 dB for FRF 1 and = 1.3 dB for FRF 2. The ROM gives an excellent prediction everywhere, except for low amplitudes. Concerning the numerical cost, the computation of the global ROB took about three minutes whereas the computation of the elastic modes took about one hour and forty minutes. To conclude, it is possible to extract a small-dimension subspace that allows an excellent representation of the fuel assembly global dynamics, in spite of the preponderant mass of the flexible part (the fuel rods). A satisfactory extraction of the relevant dynamics can be performed through the introduction of a rather simple reduced kinematics, although it is not necessarily as good as an automatic mode selection.

Boiling water reactor fuel assembly

This last application considers the vibration analysis of a detailed nuclear fuel assembly. The fuel assembly model is based on the General Electric GE14 design. This fuel assembly is used in boiling water nuclear reactors (BWRs). According to General Electric (GE), it is one of the most popular, widely used fuel bundle designs in the world. [. . . ] Today, GE14 is the workhorse of the global BWR fleet, fueling reactors throughout Europe and North America. [. . . ] To date, more than 34,000 GE14 bundles have been installed in BWRs globally.

The fuel assembly is composed of a 10 × 10 regular array of fuel rods held together by eight regularly-spaced intermediate spacer grids and by upper and lower tie plates. Eight fuel rods are removed to let space to two hollow water rods (with outer diameter 24.13 mm and thickness 0.762 mm) whose main function is neutron moderation. The water rods also support the spacer grids. Among the 92 fuel rods, there are 14 partlength fuel rods with a length of 2,481.2 mm. The other fuel rods are identical to those of the simplified fuel assembly of Section 5.2. The fuel assembly is composed of several components that are held together by elastic connections. The very long and thin flexible fuel rods are held into position by the spacer grids thanks to springs and dimples. The springs and dimples of each spacer grid cell are modeled by linear elastic translational springs in the transverse plane. There is a square tube, called channel, that surrounds the array of rods throughout the whole length. The surrounding channel is in sliding fit with the upper and lower tie plates and sits onto the lower tie plate structure. There are finger springs between the channel and the lower tie plate, and the channel is attached to a post of the upper tie plate structure. In the FE model, it is assumed that the channel is clamped to the lower tie plate and to the post and there are springs between the channel and the upper tie plate (in the transverse plane). Concerning the connection between the eight spacer grids and the surrounding channel, there are protrusions to the spacer grids, in the form of bubble-like projections, that maintain the spacer grids into position with respect to the channel box. Again, there are springs to model this connection (in the transverse plane).

The fuel assembly length is of 4,517.7 mm and its width of 140.2 mm. The fuel assembly mass is of 279.2 kg, among which 82.5% belongs to the fuel rods. Figure 50 shows three local views of the fuel assembly model. From the top image, it can be seen the upper tie plate structure, which includes a handle as well as posts, and which is surrounded by the channel. The middle image depicts the fuel assembly at the elevation of the bottommost spacer grid. The thin faces of the channel have been removed and it is noted that the remaining part of the channel, the corners, is thicker. It can be seen that the water rod diameter is smaller (13.32 mm) in the small bottom region of the fuel assembly (this is also the case in the top region). It can also be seen that the spacer grids are made up of perpendicular thin plates, which constitutes a great simplification in comparison to the original designs (the spacer grids are designed for optimized flow under reactor operation). Through each spacer grid square cell goes a fuel rod, which is represented by a light-gray line. In the bottom image of Fig. 50, the lower tie plate bottom structure is represented. It can be seen that the whole fuel assembly weight is supported by a slender bottom nose piece. Similarly to Section 5.2, the fuel rods are modeled with beam elements. Except the upper and lower tie plate structures, most of the rest of the structure is modeled using thin shell elements. The FE model has N = 1,912,506 DOFs. We are interested in the vibration analysis in the broad frequency band B ω = 2π×]0, 1000] rad/s. Concerning boundary conditions, the extremity of the bottom nose piece is fully constrained and, for each channel face there are 2 nodes located close to the top and close to the curved channel corners and for which the normal translation DOFs are constrained. These boundary conditions aim at representing the situation for which the fuel assembly is placed inside a basket slot in a nuclear canister and cask, in vertical configuration. The constraints at the top of the fuel assembly channel represent the spacer buttons and centering springs that hold the fuel assembly in place within its basket slot. There are 5946 elastic modes in B ω and 6605 elastic modes below 1200 Hz. The reference ROM is obtained by using the first n e = 6605 elastic modes as a projection basis. Figure 51 depicts the fuel assembly, a global mode, and a local mode. These modes are very similar to those depicted in Fig. 41 for the simplified fuel assembly.

In the sequel, four additional ROMs are constructed. Their ability to describe the fuel assembly global dynamics is assessed. To do so, two control FRFs are considered. For both of them, the same point load is applied in the normal direction (x) to one node of the upper part of the stiff bottom structure, and for which the position is not centered. For FRF 1, the accelerance following the same direction (x) is evaluated at one node located in one of the adjacent (perpendicular) faces. For FRF 2, the accelerance is evaluated at the same node, but following the normal direction (y). This work is part of an effort toward damage identification of the canister internals (in particular, for fuel integrity assessment) by a nonintrusive approach. In a loaded nuclear canister in vertical configuration, the fuel assemblies rest on the canister bottom plate through the fuel assembly bottom nose piece. As the canister is completely sealed for radiation protection, vibration excitation and monitoring is only possible over the outer surface of the canister. Therefore, in this paper, we are not interested in predicting the local behavior of the fuel assembly. For instance, we are not interested in predicting the response of the fuel rods. In Section 5.2, we were not interested in predicting the response of the fuel rods either. However, the effect of the fuel rods on the global dynamics was very important as the fuel rods represented about 94.8% of the total mass. Nevertheless, the effect of the local vibrations of the fuel rods were not very important, as a ROM with great reduction could be obtained. For fuel integrity assessment, it would thus be needed to have a ROM that is able to take into account the effect of the small scales reflected on the dynamics of the big scale (the canister scale). In this section, we are thus interested in predicting the global dynamics of the detailed fuel assembly, with no interest in predicting the dynamics of the thin channel or that of the fuel rods, of the spacer grids, etc. This is why we considered control FRFs with input and output locations belonging to the stiff bottom structure. The subdomain Ω G for which we would like our ROM to be valid is defined as the union of the stiff bottom structure and the upper tie plate structure. The subset S G of entries of the FRF matrix that we aim at predicting is based on a random selection of 80 nodes belonging to Ω G . For each node, the 3 translations are considered and consequently, S G comprises 240 DOFs. Similarly to Section 5.2, the random selection is based on the introduction of a homogeneous domain partition and, for each subdomain, one node is randomly selected (with equiprobability) among those that belong to its restriction to Ω G .

The simplified mode sorting procedure of Section 3.3 is used to construct a ROM of dominant elastic modes. Concerning the construction of the global ROM, three cases are investigated, similarly to Section 5.2. For all of them, a domain decomposition is introduced, namely the fuel assembly is cut into 100 transverse slices of equal thickness. The first global ROM is constructed by considering a piecewise constant kinematics for the kinetic energy. The second one is constructed by considering a piecewise rigid body kinematics. The third global ROM is based on a combined kinematics, for which the flexible part is the fuel rods (similarly to Section 5.2). That is to say, the kinematics of the water rods, channel, spacer grids and stiff lower and upper tie plate structures is kept exact. Again, the reduced kinematics for the fuel rods is that of a rigid body motion by slice. It should be noted that the corresponding ROM by static condensation did not give satisfactory results (the fuel rods account for about 82.5% of the structure total mass).

The global ROM by piecewise constant kinematics is constructed with a truncation ν = 185 such that √ σ ν /2π ≤ 5000 < √ σ ν+1 /2π, which yields a dimension n g = 151 deduced by frequency truncation with respect to f cut = 1200 Hz. The error of this global ROM in predicting the control FRF 1 is = 6.1 dB (see Fig. 52). It can be seen that compared to the simplified fuel assembly of Section 5.2, the dynamic content of the response is more important. The accuracy of the global ROM by piecewise constant kinematics is reasonable, although it misses many resonance peaks in both the low-and mid-frequency parts in addition to not having the correct amplitude level in the HF part of the band. For control FRF 2, the error is far more important, especially in the LF part (see Fig. 56). The error is = 11.6 dB. The global ROM by piecewise rigid body kinematics is constructed with a truncation ν = 297 such that √ σ ν /2π ≤ 5000 < √ σ ν+1 /2π, which yields a dimension n g = 217 deduced by frequency truncation with respect to f cut = 1200 Hz. The error of this global ROM in predicting the control FRF 1 is = 4.1 dB (see Fig. 53). It can be seen that the overall prediction is improved, in particular that of the HF part of the band. The error of this global ROM in predicting the control FRF 2 is = 8.2 dB (see Fig. 57). Like the global ROM by piecewise constant kinematics, the prediction of FRF 2 is less satisfactory. Nevertheless, in comparison, the prediction in the LF part is far better. The global ROM by combined kinematics is constructed with a truncation ν = 815 such that √ σ ν /2π ≤ 1500 < √ σ ν+1 /2π, which yields a dimension n g = 681 deduced by frequency truncation with respect to f cut = 1200 Hz. The error of this global ROM in predicting the control FRF 1 is = 2.6 dB (see Fig. 54). The prediction is greatly improved, although there are some noticeable discrepancies in the frequency band 2π×]150, 200] rad/s. The error of this global ROM in predicting the control FRF 2 is = 6.6 dB (see Fig. 58).

Conclusions

An automatic mode selection procedure has been presented which allows to get rid of the numerous insignificant local modes that may arise in the vibration analysis of complex multiscale structures. It has been shown effective in capturing a small-dimension dominant subspace of the global vibration response of the detailed model of a nuclear fuel assembly. An alternative reduced-order model based on the direct computation of a small-dimension basis of global modes has been presented. The construction of this global reduced-order model is based on the introduction of a general framework that exploits the structural scale separation in order to filter the local deformation modes. The proposed framework allows the local inertia and vibrations to influence the dominant global vibrations via the introduction of a reduced representation (kinematics) for the kinetic energy. The proposed method has been proven effective in capturing the dominant global subspace for simple structures. For the complex fuel assembly, the use of a simple reduced kinematics allowed to obtain a decent approximation of the dominant global subspace. The introduction of a more sophisticated kinematics that better takes into account the influence of the massive fuel rods onto the fuel assembly dynamics would be required to reach the same level of accuracy as the automatic selection. The construction of the global reduced-order model is more efficient in that it does not imply the computation of the numerous local modes. An efficient numerical implementation has been given and analyzed, which can easily be integrated into existing codes.

Figure 1 :Figure 2 :

 12 Figure 1: Heterogeneous plate with two structural scales (flexible panels in light blue, stiff skeleton in dark red)

  Thus, the global ROM can substitute the classic ROM of modal analysis presented in Section 2 by simply replacing [Φ] by [Φ g ] and [Λ] by [Λ g ]. It should be noted that the construction of the global modes depends on the value of ν. The choice for the value of truncation parameter ν is now discussed.

  instead of [W] = [M] as the weight matrix leads to ur 1 = u1 and ur 2 = 0, and it can be deduced that [P] = I N1 0 0 0 , [M r ] = M 11 0 0 0 , p = [M r ] u = M 11 u1 0 , (40) which corresponds to the definition of static condensation for the case where the master DOFs are all the DOFs belonging to Ω 1 .

Figure 5 :

 5 Figure 5: Orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. (36).

Figure 6 :

 6 Figure 6: Orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. (40).

  [P] and [M r ] and for product [M r ] u are now given, for the case of weight matrix [W] = [M] and reduced basis [B] = [T (0)]. It can be verified that the equality [W] = [M] implies low-rank mass matrix [M r ] to be such that [M r ] = [M][P]. Introducing the Guyan condensation matrix, [

Figure 7 :

 7 Figure 7: Orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. (41).

Figure 8 :

 8 Figure 8: Orthogonal projection of the multiscale mode of Fig. 4 given by the matrix [P] defined in Eq. (44).

  in which the (N 2 × r 2 ) real matrix[b] provides for a reduced representation of dimension r 2 for the DOFs of subdomain Ω 2 (substructures) that is independent of the exact representation of the DOFs of subdomain Ω 1 (main frame). The choice of [b] is intended to represent the global dynamics of the substructures (which emerges by coupling with the structure skeleton), while preventing the representation of the local dynamics of said substructures. The construction of matrix [b] follows the strategy presented in Section 4.2 for the case of structures without scale separation, as subdomain Ω 2 does not exhibit any stiff main frame. For the case [W] = [M], the expressions for [P] and [M r ] = [M][P] are quite long and are not given. The expression of product p = [M r ] u is sufficient to solve Eq. (

  In order to remove the cost associated with Eq. (57), one solution is to simplify weight matrix [W] by considering [W 12 ] = [0] and [W 21 ] = [0]. The implication of this simplification is small because the coupling terms of the mass matrix are not preponderant. By doing so, the expressions for [P] and for [M r ] = [P] T [M][P] = [M][P] are greatly simplified, namely [P] =

  matrix [M] by low-rank mass matrix [M r ]. For the case [W] = [M] but considering a diagonally-lumped approximation of [M], the expressions for [P] and [M r ] = [M][P] are quite simplified. They are the same as in Eq. (58), except that the extradiagonal blocks of [M r ] are zero. For all u in R N , product p = [M r ] u is thus given by p = (p 1 , p 2 ) T with p 1 = [M 11 ] u1 and p 2 = [M 22 ] ur 2 , in which ur 2 = [P 2 ] u2 . Similarly to before, the use of the orthonormalization [ b] of [b] allows great computational savings for the online computation of p = [M r

Figure 12 :Figure 13 :

 1213 Figure 12: Plot of the sorted modal importance

Figure 14 :

 14 Figure14: Plot of the error function (n) with respect to the ROM dimension n, for several ROMs: the ROM of the first n = n e elastic modes (dashed line), the ROM of the n = n p elastic modes with largest participation factor (dotted line), the ROM of the n = n d dominant modes given by the simplified mode sorting procedure (black line), and the dominant ROM of n = n D modes given by the tractable mode sorting procedure (gray line); plot of the point (n g = 56, (n g ) = 2.7) for the global ROM given by static condensation (black dot).

Figure 15 :Figure 16 :

 1516 Figure 15: Control FRF calculated with the dominant ROM of dimension n D = 56/1156 given by the tractable mode sorting procedure and for which = 1.8 dB

Figure 17 :

 17 Figure 17: Plot of the sorted modal importance

Figure 18 :

 18 Figure18: Plot of the error function (n) with respect to the ROM dimension n, for several ROMs: the ROM of the first n = n e elastic modes (dashed line), the ROM of the n = n p elastic modes with largest participation factor (dotted line), the ROM of the n = n d dominant modes given by the simplified mode sorting procedure (black line), and the dominant ROM of n = n D modes given by the tractable mode sorting procedure (gray line); plot of the point (n g = 8, (n g ) = 10.9) for the global ROM given by static condensation (black dot).

Figure 19 :

 19 Figure 19: Plot of function ν → n g (ν) for several global ROMs: piecewise rigid body kinematics (black line), polynomial kinematics (black dashed line), piecewise constant kinematics (gray line), polynomial kinematics with discriminatory weight matrix (gray dashed line), and static condensation (dotted line).

  compares the first global ROM (obtained by piecewise constant kinematics) with the global ROM by static condensation. Figure 21 compares the second global ROM (obtained by piecewise rigid body kinematics) with the global ROM by static condensation. Figure 22 compares the third global ROM (obtained by degree-5 polynomial kinematics) with the global ROM by static condensation. Figure 23 compares the fourth global ROM (obtained by degree-5 polynomial kinematics with discriminatory weight matrix [W]) with the global ROM by static condensation.

  Fig.23), it can be observed that the global eigenfrequencies f g α coincide very closely with those of static condensation. Second, for this kinematics, the eigenvalues σ α and λ g α are very close, far more than it is the case for static condensation. To interpret this, it is recalled that Guyan and static condensations yield very close FRFs (hence very close eigenvalues λ g α although not necessarily very close eigenvalues σ α ), and that no Rayleigh-Ritz step is necessary for Guyan condensation since λ g α = σ α for this case (as explained in Section 4.1). It can be deduced that the orthogonal projection given by the polynomial kinematics and defined with the discriminatory weight matrix [W] is a very good approximation to the orthogonal projection associated with Guyan condensation. It should be noted that the deformation shapes in Figs. 7 and 11 look the same, although the orthogonal projection was performed using [W] = [M].

Figure 20 :

 20 Figure 20: Plot of the functions α → √ σ α /2π (dotted lines) and α → f g α (ν) for ν = 1, . . . , r (solid lines) for static condensation (in gray) and for piecewise constant kinematics (in black).

Figure 21 :

 21 Figure 21: Plot of the functions α → √ σ α /2π (dotted lines) and α → f g α (ν) for ν = 1, . . . , r (solid lines) for static condensation (in gray) and for piecewise rigid body kinematics (in black).

Figure 22 :

 22 Figure 22: Plot of the functions α → √ σ α /2π (dotted lines) and α → f g α (ν) for ν = 1, . . . , r (solid lines) for static condensation (in gray) and for polynomial kinematics (in black).

Figure 23 :

 23 Figure 23: Plot of the functions α → √ σ α /2π (dot-

Figure 24 :

 24 Figure24: Plot of the error function (n g (ν)) with respect to the truncation parameter ν, for several global ROMs: piecewise rigid body kinematics (black line), polynomial kinematics (black dashed line), piecewise constant kinematics (gray line), and polynomial kinematics with discriminatory weight (gray dashed line); plot of the point (ν = 8, (n g (ν)) = 10.9) for the global ROM given by static condensation (black dot).

Figure 25 :Figure 26 :

 2526 Figure 25: Control FRF calculated with the dominant ROM of dimension n D = 8 given by the tractable mode sorting procedure and for which = 1.8 dB

Figure 27 :Figure 28 :

 2728 Figure 27: Control FRF calculated with the global ROM of dimension n g = 14 obtained by piecewise constant kinematics and for which = 6.7 dB

Figure 29 :

 29 Figure 29: Control FRF calculated with the global ROM of dimension n g = 16 obtained by polynomial kinematics and for which = 4.9 dB

Figure 31 :Figure 32 :Figure 33 :

 313233 Figure 31: Plot of the sorted modal importance (black line) juxtaposed with that of Fig. 12 (gray line)

Figure 34 :Figure 35 :

 3435 Figure 34: Control FRF calculated with the dominant ROM of dimension n D = 60/1157 given by the tractable mode sorting procedure and for which = 6.2 dB

Figure 36 depicts the

  Figure36depicts the mode importance as a function of the dominant mode rank, in comparison to the same curves that were obtained for previous structures. This figure clearly shows that the important contributions are far less concentrated than for previous structures. It can be seen that the removal of a dominant mode of rank beyond 500 still induces an error level of about 1 dB. It is anticipated that the local vibrations of the flexible panels have a very important contribution to the global dynamics of the stiff skeleton. Figure37depicts the function ν → n g (ν) for the global ROM given by static condensation of the flexible part. Dimension n g reaches a plateau, but not as quickly as for previous structures with lower panel mass. For ν = 121, . . . , 264 the ROM dimension is n g = 83 whereas it is n g = 84 beyond that. It has been verified that the global modes ϕ g effectively exhibit global deformation shapes.

Figure 36 :Figure 37 :Figure 38 :

 363738 Figure 36: Plot of the sorted modal importance (black line) juxtaposed with those of Fig. 31 (gray lines)

Figure 38

 38 Figure38depicts the error function (n) with respect to the ROM dimension n, for the 5 different ROMs aforementioned. In the present case, the global ROM by static condensation is not Pareto dominant. Nevertheless, it is still on the trajectory of the error functions of the ROMs by automatic mode selection. The dominant ROM of the tractable mode sorting procedure is Pareto dominant. The extent of the error introduced by the proposed small-dimension ROMs is illustrated in Figs.39 and 40, which depict the resulting approximated control FRF in comparison to the exact one computed with the ROM of the n e = 1212 elastic modes. It can be seen that for both ROMs the error is very significant. The FRF exhibits an even more enriched dynamic content (compared to the previous structures with lower panel mass), which fails to be represented by the small-dimension ROMs. The following conclusions can be drawn regarding this simple example of a heterogeneous plate. First, the complexity of the dynamic response increases with the mass of the flexible panels. Second, the contribution of the local modes increases with the mass of the flexible panels. Finally, for a fixed accuracy, the number of dominant modes to be retained increases with the mass of the flexible panels. It should be noted that there is no guarantee that the tractable mode sorting procedure (defined in Section 3.2) give the optimal mode selection. In addition, the optimal ROM of elastic modes may not be the optimal ROM (for instance, this was not the case for the application of Section 5.1.1). Therefore, it is unsure whether there exists a small-dimension subspace that would be able to capture the main dynamic content of the global response of this structure, let alone the fast fluctuations of said response with respect to the frequency.

Figure 39 :Figure 40 :

 3940 Figure 39: Control FRF calculated with the dominant ROM of dimension n D = 84/1212 given by the tractable mode sorting procedure and for which = 10.2 dB

  3 mm. The spacing between two adjacent rods (rod pitch) is of 13.0 mm. Each tie plate is square with side of 129.6 mm and thickness of 10 mm. The fuel rods are mainly made up of uranium dioxide (Young modulus of 230 GPa, Poisson ratio of 0.32, mass density of 10,960 kg/m 3 ). The intermediate tie plates are made up of Zircaloy (Young modulus of 88 GPa, Poisson ratio of 0.34, mass density of 6500 kg/m 3 ), whereas the upper and lower tie plates are made up of steel (Young modulus of 200 GPa, Poisson ratio of 0.29, mass density of 7800 kg/m 3 ). The total mass of the fuel rods is of 262.0 kg whereas the total mass of the structure is of 276.4 kg. The tie plates are modeled with 10,240 Belytschko-Tsay shell elements and the fuel rods with 51,500 Belytschko-Schwer tubular beam elements. The FE model has N = 368,892 DOFs. Figure 41 depicts the fuel assembly, a global mode, and a local mode.

Figure 41 :

 41 Figure 41: Fuel assembly: undeformed configuration (left), a global mode (center), a local mode as seen from two points of view (right)

Figure 42 : 5 Figure 43 :

 42543 Figure 42: Control FRF 1 calculated with the global ROM of dimension n g = 76/5441 obtained by piecewise constant kinematics (constant by slice) and for which = 5.6 dB

Figure 44 :Figure 45 :

 4445 Figure 44: Control FRF 1 calculated with the global ROM of dimension n g = 113/5441 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 4.2 dB

Figure 46 :Figure 47 :Figure 48 :Figure 49 :

 46474849 Figure 46: Control FRF 2 calculated with the global ROM of dimension n g = 76/5441 obtained by piecewise constant kinematics (constant by slice) and for which = 51.0 dB

Figure 50 :

 50 Figure 50: Fuel assembly details: view of the top (top), and views of the midspan (middle) and of the bottom (bottom) without the channel thin faces

Figure 51 :

 51 Figure 51: Fuel assembly: undeformed configuration with and without the channel thin faces (left) and, without the channel thin faces, a global mode (center) and a local mode as seen from two points of view (right)

Figure 52 :

 52 Figure 52: Control FRF 1 calculated with the global ROM of dimension n g = 151/6605 obtained by piecewise constant kinematics (constant by slice) and for which = 6.1 dB

Figure 53 :

 53 Figure 53: Control FRF 1 calculated with the global ROM of dimension n g = 217/6605 obtained by piecewise rigid body kinematics (rigid body by slice) and for which = 4.1 dB

Figure 54 :Figure 55 :

 5455 Figure 54: Control FRF 1 calculated with the global ROM of dimension n g = 681/6605 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 2.6 dB

  and 0.61 × k 1 N/m 2 , with k 1 = 10 12 . This variability in the Young modulus of the flexible panels is such that there is variability in the eigenfrequencies of the local modes of said panels. The FE model is constituted of 4-node thin shell elements for which each node has 3 DOFs, namely TZ, RX, and RY (with z the direction perpendicular to the plate). As boundary conditions, the four corner nodes are fully constrained. The FE model has N = 53,223 DOFs.

Acknowledgement

The financial support for this project was provided by the United States Department of Energy through the Nuclear Energy University Program under the Contract No. DE-NE0008529. The findings presented herein are those of the authors and do not necessarily reflect the views of the sponsor.

Again, the error for FRF 2 is more important. This can be explained by the fact that FRF 2 exhibits more complexity than FRF 1 (i.e., exhibits more frequency fluctuations). To conclude, it is possible to extract a small-dimension subspace that allows a good representation of the fuel assembly global dynamics, in spite of the preponderant mass of the flexible part (the fuel rods). The relevant dynamics can be decently approximated through the introduction of a rather simple reduced kinematics, although it is not as satisfactory as an automatic mode selection. The remaining discrepancy is due to the reduced kinematics of the flexible part (the fuel rods). A better global ROM would require the introduction of a suitable kinematics (or suitable shape functions) for representing the global displacements of the flexible part that arise from its connection and coupling with the skeleton.