Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods

Finite element analysis (FEA) is a well-established numerical simulation method for struc‐ tural dynamics. It serves as the main computational tool for Noise, Vibration and Harshness (NVH) analysis in the low-frequency range. Because of developments in numerical methods and advances in computer software and hardware, FEA can now handle much more com‐ plex models far more efficiently than even a few years ago. However, the demand for com‐ putational capabilities increases in step with or even beyond the pace of these improvements. For example, automotive companies are constructing more detailed models with millions of degrees of freedom (DOFs) to study vibro-acoustic problems in higher fre‐ quency ranges. Although these tasks can be performed with FEA, the computational cost can be prohibitive even for high-end workstations with the most advanced software.


Introduction
Finite element analysis (FEA) is a well-established numerical simulation method for structural dynamics. It serves as the main computational tool for Noise, Vibration and Harshness (NVH) analysis in the low-frequency range. Because of developments in numerical methods and advances in computer software and hardware, FEA can now handle much more complex models far more efficiently than even a few years ago. However, the demand for computational capabilities increases in step with or even beyond the pace of these improvements. For example, automotive companies are constructing more detailed models with millions of degrees of freedom (DOFs) to study vibro-acoustic problems in higher frequency ranges. Although these tasks can be performed with FEA, the computational cost can be prohibitive even for high-end workstations with the most advanced software.
For large finite element (FE) models, a modal reduction is commonly used to obtain the system response. An eigenanalysis is performed using the system stiffness and mass matrices and a smaller in size modal model is formed which is solved more efficiently for the response. The computational cost is also reduced using substructuring (superelement analysis). Modal reduction is applied to each substructure to obtain the component modes and the system level response is obtained using Component Mode Synthesis (CMS).
When design changes are involved, the FEA analysis must be repeated many times in order to obtain the optimum design. Furthermore in probabilistic analysis where parameter uncertainties are present, the FEA analysis must be repeated for a large number of sample points. In such cases, the computational cost is even higher, if not prohibitive. Reanalysis methods are intended to analyze efficiently structures that are modified due to various changes. They estimate the structural response after such changes without solving the complete set of modified analysis equations. Several reviews have been published on reanalysis methods [1][2][3] which are usually based on local and global approximations. Local approximations are very efficient but they are effective only for small structural changes. Global approximations are preferable for large changes, but they are usually computationally expensive especially for cases with many design parameters. The well-known Rayleigh-Ritz reanalysis procedure [4,5] belongs to the category of local approximation methods. The mode shapes of a nominal design are used to form a Ritz basis for predicting the response in a small parametric zone around the nominal design point. However, it is incapable of capturing relatively large design changes.
A parametric reduced-order modeling (PROM) method, developed by Balmes [6,7], expands on the Rayleigh-Ritz method by using the mode shapes from a few selected design points to predict the response throughout the design space. The PROM method belongs to the category of local approximation methods and can handle relatively larger parameter changes because it uses multiple design points. An improved component-based PROM method has been introduced by Zhang et al. [8,9] for design changes at the component level. The PROM method using a 'parametric' approach has been successfully applied to design optimization and probabilistic analysis of vehicle structures. However, the 'parametric' approach is only applicable to problems where the mass and stiffness matrices can be approximated by a polynomial function of the design parameters and their powers. A new 'parametric' approach using Kriging interpolation [10] has been recently proposed [11]. It improves efficiency by interpolating the reduced system matrices without needing to assume a polynomial relationship of the system matrices with respect to the design parameters as in [6,7].
The Combined Approximations (CA) method [12][13][14] combines the strengths of both local and global approximations and can be accurate even for large design changes. It uses a combination of binomial series (local) approximations, called Neumann expansion approximations, and reduced basis (global) approximations. The CA method is developed for linear static reanalysis and eigen-problem reanalysis [15][16][17][18][19]. Accurate results and significant computational savings have been reported. All reported studies on the CA method [12][13][14][15][16][17][18][19] use relatively simple frame or truss systems for static or dynamics analysis with a relatively small number of DOF and/or modes. For these problems, the computational efficiency was improved by a factor of 5 to 10. Such an improvement is beneficial for a single design change evaluation or even for gradient-based design optimization where only a limited number of reanalyses (e.g. less than 50) is performed. However, the computational efficiency of the CA method may not be adequate in simulation-based (e.g. Monte-Carlo) probabilistic dynamic analysis of large finite-element models where reanalysis must be performed hundreds or thousands of times in order to estimate accurately the reliability of a design. solving a linear system is dominated by the cost of matrix decomposition way no be longer valid (see Section 3.4) and the computational savings from using the CA method may not be substantial. For this reason, we developed a modified combined approximation (MCA) and integrated it with the PROM method to improve accuracy and computational efficiency. The computational savings can be substantial for problems with a large number of design parameters. Examples in this Chapter demonstrate the benefits of this reanalysis methodology.
The Chapter presents methodologies 1. for accurate and efficient vibration analysis methods of large-scale, finite-element models, 2. for efficient and yet accurate reanalysis methods for dynamic response and optimization, and 3. for efficient design optimization methods to optimize structures for best vibratory response.
The optimization is able to handle a large number of design variables and identify local and global optima. It is organized as follows. Section 2 presents an overview of reduced-order modeling and substructuring methods including modal reduction and component mode synthesis (CMS). Improvements to the CMS method are presented using interface modes and filtration of constraint modes. The section also overviews two Frequency Response Function (FRF) substructuring methods where two substructures, represented by FRFs or FE models, are assembled to form an efficient reduced-order model to calculate the dynamic response. Section 3 presents four reanalysis methods: the CDH/VAO method, the Parametric Reduced Order Modeling (PROM) method, the Combined Approximation (CA) method, and the Modified Combined Approximation (MCA) method. It also points out their strong and weak points in terms of efficiency and accuracy. Section 4 demonstrates how the reanalysis methods are used in vibration and optimization of large-scale structures. It also presents a new reanalysis method in Craig-Bampton substructuring with interface modes which is very useful for problems with many interface DOFs where the FRF substructuring methods cannot be used. Section 5 presents a vibration and optimization case study of a large-scale vehicle model demonstrating the value of reduced-order modeling and reanalysis methods in practice. Finally, Section 6 summarizes and concludes.
ry to obtain an abstract (mathematical) substructuring using matrix partitioning of the entire finite-element model. The computational savings from the "mathematical" and "physical" methods can be comparable depending on the problem at hand.

Modal Reduction
For an undamped structure with stiffness and mass matrices K and M respectively, under the excitation force vectorF , the equations of motion (EOM) for frequency response are where the displacement d is calculated at the forcing frequencyω. If the response is required at multiple frequencies, the repeated direct solution of Equation (1) is computationally very expensive and therefore, impractical for large scale finite-element models.
A reduced order model (ROM) is a subspace projection method. Instead of solving the original response equations, it is assumed that the solution can be approximated in a subspace spanned by the dominant mode shapes. A modal response approach can be used to calculate the response more efficiently. A set of eigen-frequencies ω i and corresponding eigenvectors (mode shapes) φ i are first obtained. Then, the displacement d is approximated in the reduced space formed by the first n dominant modes as (2) where is the modal basis and U is the vector of principal coordinates or modal degrees of freedom (DOF). Using the approximation of Equation (2), the EOM of Equation (1) can be transformed from the original physical to the modal degrees of freedom as (3) The response d can be recovered by solving Equation (3) for the modal response U and projecting it back to the physical coordinates using Equation (2). If ω max is the maximum excitation frequency, it is common practice to retain the mode shapes in the frequency range of 0 ÷ 2ω max . The system modes in the high frequency range can be safely truncated with minimal loss of accuracy.
Due to the modal truncation, the size of the ROM is reduced considerably, compared to the original model. However, the size increases with the maximum excitation frequency. An added benefit of the ROM is that Equations (3) are decoupled because of the orthogonality of the mode shapes and can be therefore, solved separately reducing further the overall computational effort.
Note that for a damped structure with a damping matrixC, Equation (1) becomes (3) is modified as by adding the modal damping term . For proportional (structural or Rayleigh) damping, C is a linear combination of M and K ; i.e. C = α M + β K where αand βare constants. In this case, the reduced Equations (3) of the modal model are decoupled. Otherwise, they are not. In this Chapter for simplicity, we present all theoretical concepts for undamped systems. However for forced vibrations of damped systems, the addition of damping is straightforward.

Substructuring with Component Mode Synthesis (CMS)
To model the dynamics of a complex structure, a finite-element analysis of the entire structure can be very expensive, and sometimes infeasible, due to computer hardware and/or software constraints. This is especially true in the mid-frequency range, where a fine finite element mesh must be used in order to capture the shorter wavelengths of vibration. Component mode synthesis (CMS) was developed as a practical and efficient approach to modeling and analyzing the dynamics of a structure in such circumstances [20][21][22][23]. The structure is partitioned in component structures and the dynamics are described by the normal modes of the individual components and a set of modes that couple all component. Besides the significant computational savings, this component-based approach also facilitates distributed design. Components may be modified or redesigned individually without re-doing the entire analysis.
One of the most accurate and widely-used CMS methods is the Craig-Bampton method [22] where the component normal modes are calculated with the interface between connected component structures held fixed. Attachment at the interface is achieved by a set of "constraint modes." A constraint mode shape is the static deflection induced in the structure by applying a unit displacement to one interface DOF while all other interface DOF are held fixed. The motion at the interface is thus completely described by the constraint modes. The Craig-Bampton reduced-order model (CBROM) results in great model size reduction by including only component normal modes within a certain frequency range. However, there is no size reduction for constraint modes because CBROM must have one DOF for each interface DOF. If the finite element mesh is sufficiently fine, the constraint-mode DOFs will dominate the CBROM mass and stiffness matrices, and increase the computational cost.
We address this problem by using interface modes (also called characteristic constraint -CC-modes) in order to reduce the number of retained interface DOFs of the Craig-Bampton approach. For that, a secondary eigenvalue analysis is performed using the constraintmode partitions of the CMS mass and stiffness matrices. The CC modes are the resultant eigenvectors. The basic formulation is described in Sections 2.2.1 and 2.2.2. The interface modes represent more "natural" physical motion at the interface. Because they capture the characteristic motion of the interface, they may be truncated as if they were traditional modes of vibration, leading to a highly compact CC-mode-based reduced order model (CCROM).
In addition, the CC modes provide a significant insight into the physical mechanisms of vibration transmission between the component structures. This information could be used, for example, to determine the design parameters that have a critical impact on power flow. Figure 1 compares a conventional constraint mode used in Craig-Bampton analysis with an interface mode, for a simple cantilever plate which is subdivided in two substructures.
It should be noted that the calculation of the CC modes is essentially a secondary modal analysis. Therefore, the benefits are the same as those of a traditional modal analysis. For instance, refining the finite-element mesh increases the accuracy of a CCROM without introducing any additional degrees of freedom. The ability of the CC mode approach to produce CCROM whose size does not depend on the original level of discretization makes it especially well suited for finite-element based analysis of mid-frequency vibration.

Craig-Bampton Fixed Interface CMS
This section provides the basics of Craig-Bampton method using the fixed-interface assump- The fixed-interface Craig-Bampton CMS method utilizes two sets of modes to represent the substructure motion; substructure normal (N) modesΦ i N , and constraint (C) modesΦ i C , where i denotes the i th substructure. The size reduction of the Craig-Bampton method comes from the truncation of the normal modes Φ i N = φ i1 φ i2 ⋯ φ in which are calculated by fixing all interface DOFs and solving the following eigenvalue problem The static constraint modesΦ i C are calculated by enforcing a set of static unit constraints to the interface DOFs as The original physical DOFs i G u and i W u can be thus represented by the constraint-mode DOFs C i u and the normal-mode DOFs N i u as Using the above Craig-Bampton transformation, the original EOM of Equation (7), can be expressed as where the superscripts C and N are used to indicate partition related to static constraint mode DOFs and fixed-interface normal mode DOFs, respectively. The matrix partitions of Equation (8) are (9) The matrices of each substructure are then assembled by applying displacement continuity and force balance along the interface to obtain the EOM of the reduced system. A secondary modal analysis is finally carried out using the mass and stiffness matrices of the reduced system to obtain the eigenvalues and eigenvectors.
Note that constraint mode matrix C i Φ is usually a full matrix. Therefore Equation (9) can be computationally expensive due to the triple-product ( )  (17) where n ss is the number of substructures in the global structure. The corresponding synthesized CMS mass and stiffness matrices are as follows  of the interfaces between components and is therefore, determined by the finite-element mesh. If the mesh is fine in the interface regions, or if there are many substructures, the constraint-mode partitions of the CMS matrices may be relatively large. For this reason, we further reduce the CMS matrices by performing a modal analysis on the constraint-mode DOFs as follows k C ψ n = Λ n m C ψ n for n = 1, 2, 3, ...
The eigenvectors ψ n are transformed into the finite-element DOFs for the i th component structure using the following transformation where is a selected set of n CC interface eigenvectors which are few compared to the number of the Finally, the CMS matrices are transformed using the CC modes and the reduced-order CMS matrices are obtained similarly to Equations (18). Now, the unknown displacement vector d ROM is partitioned as where superscript CC indicates the partition associated with the CC modes. The equations of motion of the reduced order CMS model (ROM) are expressed by The mass matrix M ROM , the stiffness matrix K ROM , and the applied force vector f ROM , are explicitly written as  where (27) and Figure 2 shows a typical constraint mode for a plate substructure. The non-zero displacement field (indicated by red color) is usually limited to a very small region close to the perturbed interface DOF. If the small-displacement part of the constraint mode shape is explicitly replaced by zero, the density of the resulting "filtered" constraint mode will be significantly reduced. Consequently, the computational cost in Equation (9) will be considerably reduced. To filter the constraint modes, the following criterion is used

Filtration of Constraint Modes
where φ pq C is the p th element of the q th constraint mode C φ . If the ratio of an element of the constraint mode vector to the maximum value in the vector is smaller than a defined smallε, the element of the constraint mode is truncated to zero. For the constraint mode of Figure 4, the constraint mode density reduces from 100% to 16% ifε = 0.03.

Frequency Response Function (FRF) Substructuring and Assembling
If the number of interface nodes (or DOFs) between connected substructures is small, a reduced-order model of small size can be developed using an FRF representation of each substructure. This is known as FRF substructuring. The FRF representation can be easily obtained from a finite element (FE) model or even experimentally. If the FE model of one substructure is very small (e.g. a vehicle suspension model), it can be easily coupled directly to another substructure which is represented using FRFs. This section provides the fundamentals of FRF substructuring for both FRF-FE and FRF-FRF coupling.

Algorithm for FRF/FE Coupling
The numerical algorithm is explained using the two-substructure example of Figure 3. Sub1 is an FRF type substructure, meaning that its dynamic behavior is described using FRF matrices which are denoted by H (see Equation 30 for notation). The elements of H are frequency dependent and complex if damping is present. A bold letter indicates a matrix or vector. According to Equation (30), H AC for example, represents the displacement X A of DOF A due to a unit force F C on DOF C. Sub2 is a finite element (FE) type substructure. Its dynamic behavior is described using the stiffness K, mass M and damping B matrices. The FRF matrix of Sub1 can be calculated by either a direct frequency response method or a modal response method. In the former case, the original FE equations are used in the frequency domain while in the latter a modal model is first developed and then used to calculate the FRF matrix. The size of the FRF matrix is small and depends on the number of DOFs of the excitation, response and interface DOFs. Usually FRFs are calculated between excitation and response DOFs. However in order to assemble two substructures, FRFs are also calculated between interface DOFs and excitation/response DOFs. The Sub1 FRF matrix H in Equation (30) is thus partitioned into interior (A) DOFs and interface/coupling (C) DOFs.
The interior DOFs include all excitation and all response DOFs ( Figure 3).
The second substructure Sub2 is expressed in FE format. The system FE matrices K, M, and B form the frequency dependent dynamic matrix Z = K + iωB − ω 2 M which is then partitioned according to the interior (B) and interface (C) DOFs. The interface DOFs for Sub1 and Sub2 have the same node IDs so that they can be assembled to obtain the system FRF matrix.
The procedure to assemble the H matrix of Sub1 with the Z matrix of Sub2 and calculate (solve for) the system matrix H is described below.
The equations of motion for Sub1 are expressed as where subscript A indicates the interior (excitation plus response) DOFs of Sub1, and subscript C indicates the connection/common/coupling DOFs between Sub1 and Sub2. The equations of motion for Sub2 are expressed as where subscript B indicates the interior DOFs of Sub2, and subscript C indicates the connection/common/coupling DOFs between Sub2 and Sub1. Because of displacement compatibility at the interface, X C appears on the left-hand side of Equation (30) for Sub1 and on the right-hand side of Equation (31) for Sub2. Superscripts 1 and 2 are used to differentiate the interface forces F C at Sub1 and Sub2.
To couple Sub1 and Sub2, compatibility of forces at the interface is applied as  where S = [Φ Θ I] and . Figure 4 shows the coupling of two FRF type substructures. The equations of motion for Sub1 and Sub2 and are expressed as To couple Sub1 and Sub2, we enforce displacement compatibility at the interface and apply the interface force relationship . In this case, the assembled system equations can be re-arranged in matrix form as Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402

Algorithm for FRF/FRF Coupling
The FRF/FRF coupling is a special case of the FRF/FE coupling.

CDH/VAO Method
The CDH/VAO method, developed by CDH AG for vibro-acoustic analysis, is a Rayleigh-Ritz type of approximation. If the stiffness and mass matrices of the baseline design structure are K 0 and M 0 , the exact mode shapes in Φ 0 are obtained by solving the eigen-problem (44) where Λ 0 is the diagonal matrix of the baseline eigenvalues. A new design (subscript p) has the following stiffness and mass matrices For a modest design change where ΔK and ΔM are small, it is assumed that the change in mode shapes is small and the new response can be therefore, captured in the sub-space spanned by the mode shapes Φ 0 of the initial design. The new stiffness and mass matrices are then condensed asK R = Φ 0 T K p Φ 0 and M R = Φ 0 T M p Φ 0 and the following reduced eigenvalue problem is solved to calculate the eigen-vector Θ The approximate eigenvalues of the new design are given by Λ p and the approximate eigenvectors Φ p are whereR = Φ 0 . Thus, the modal response of the modified structure can be easily obtained and the actual response can be recovered using the eigenvectorsΦ p .

Parametric Reduced-Order Modeling (PROM) Method
The PROM method approximates the mode shapes of a new design in the subspace spanned by the dominant mode shapes of some representative designs, which are selected such that the formed basis captures the dynamic characteristics in each dimension of the parameter space. Balmes et al. [6,7] suggested that these representative designs should correspond to the middle points on the faces of a box in the parameter space representing the range of design parameters. For a structure with m design variables, Zhang [9] suggested that the representative designs include a baseline design for which all parameters are at their lower limits plus m designs obtained by perturbing the design variables from their lower limits to their upper limits, one at a time. The points representing these designs in the space of the design variables are called corner points (see Figure 5). This selection of representative designs results in a more accurate PROM algorithm. The mode shapes of a new design are approximated in the space of the mode shapes of the corner points as where the modal matrix P includes the basis vectors as in Equation (49) and Θ represents the participation factors of these vectors. The columns of P are the dominant mode shapes of the above ( m + 1 ) designs, where Λ is a diagonal matrix of the first s eigenvalues.
A reduced eigenvalue problem is obtained by pre-multiplying both sides of Equation (50) by where the reduced stiffness and mass matrices are Thus, the matrix Θ in Equation (48) consists of the eigenvectors of the reduced stiffness and mass matrices K R andM R .
For m design variables, ( m + 1 ) eigenvalue problems must be solved in order to form the basis P of Equation (49). Therefore, both the cost of obtaining the modal matrices Φ i and the size of matrix P increase linearly with m. The PROM approach uses the following algorithm to compute the mode shapes of a new design:

1.
Calculate the mode shapes of the baseline design and the designs corresponding to the m corner points in the design space, and form subspace basisP.

2.
Calculate the reduced stiffness and mass matrices K R and M R from Equation (52).
Step 1 is performed only once. A reanalysis requires only steps 2 to 4. For a small number of mode shapes and a small number of design variables, the cost of steps 2 to 4 is much smaller than the cost of a full analysis. The computational cost of PROM consists of ). When the size of basis P increases so does the fixed cost because more eigenvalue problems and mode shapes must be calculated to form basisP. The PROM method results in significant cost savings when applied to problems that involve few design variables (less than 10) and require many analyses (e.g. Monte Carlo simulation or gradient-free optimization using genetic algorithms).

Combined Approximations (CA) Method
The PROM method requires an eigenvalue analysis for multiple designs (corner points) to form a basis for approximating the eigenvectors at other designs. It is efficient only when the number of design parameters is relatively small. On the contrary, the CA method of this section does not have such a restriction because the reanalysis cost is not proportional to the number of design parameters. The CA method is thus more suitable than the PROM method, when the number of reanalyses is less than or comparable to the number of design parameters, such as in gradient-based design optimization.
The fundamentals of the combined approximations (CA) method [15][16][17][18][19] are given below. A subspace basis is formed through a recursive process for calculating the natural frequencies and mode shapes of a system. If K 0 and M 0 are the stiffness and mass matrices of the original (baseline) design, the exact mode shapes Φ 0 are obtained by solving the eigen-problem The eigen-problem for the modified design can be expressed as leading to the following recursive equation which produces a sequence of approximations of the mode shapesΦ p, j , j = 1, 2, … , s. The CA method uses the changes R j = Φ p, j − Φ p, j−1 to form a subspace basis to approximate the modes of the new design. In order to simplify the calculations, showing that the basis vectors satisfy the following recursive equation where the first basis vector is assumed to beR 1 The CA method forms a subspace basis [ ] where s is usually between 3 and 6 [16][17][18]23] and the mode shapes of the new (K p ,M p ) design are then approximated in the subspace spanned by R using the following algorithm: • Condense the stiffness and mass matrices as • Solve the reduced eigen-problem (using matrices K R and M R ) to calculate the eigenvector matrixΘ.
• Reconstruct the approximate eigenvectors of the new designΦ p as The eigenvalues of the new design are approximated by the eigenvalues λ p of the reduced eigen-problem.
The CA method has three main advantages: 1. it only requires a single matrix decomposition of stiffness matrix K 0 in Equation (55) (58) where the eigenvectors Θ correspond to a much smaller reduced eigen-problem.
However for a large number of reanalyses, the computational cost can increase substantially because a new basis and the condensed mass and stiffness matrices in Equation (57) must be calculated for every reanalysis. Examples where many analyses are needed are optimization problems in which a Genetic Algorithm is employed to search for the optimum, and probabilistic analysis problems using Monte-Carlo simulation. The PROM method can be more suitable for these problems because the subspace basis R does not change for every new design point. Note that steps 1 and 3 (Equations 57and 58) are similar to steps 2 and 4 of PROM. CA uses basis R and PROM uses basisP.
The CA method is more efficient than PROM for design problems where few reanalyses are required for two reasons. First, it does not require calculation of the eigenvectors Φ i , i = 1, ⋯ , m, of the m corner design points, and second the cost of matrix condensation of Equation (57) is much lower than that of Equation (52), because the size (number of columns) of basis R is not proportional to the number of parameters m as in basis P . For problems with a large number of design parameters, the PROM approach is efficient only when a 'parametric' relationship is established [7] because a large overhead cost, proportional to the number of design parameters, is required. In contrast, the CA method does not require such an overhead cost because the reanalysis cost is not proportional to the number of design parameters. The CA method is therefore, more suitable than PROM, if the number of reanalyses is less than or comparable to the number of design parameters.

Modified Combined Approximations (MCA) Method
In the literature, the accuracy and efficiency of the CA method has been mostly tested on problems involving structures with up to few thousands of DOFs, such as frames or trusses [12][13][14][15][16][17][18][19]. We have tested the CA method using, among others, the structural dynamics response of a medium size (65,000 DOFs) finite-element model (Figure 7 of Section 4.3). Due to its high modal density, there were more than two hundred dominant modes in the frequency range of zero to 50 Hz. It was observed that the computational savings of the CA method, using the recursive process of Equation (55), were not substantial. For this reason, we developed a modified combined approximations method (MCA) by modifying the recursive process of Equation (55) which is much more efficient than the original CA method for large size models.
The cost of calculating the subspace basis in Equation (55) consists of one matrix decomposition (DCMP) and one forward-backward substitution (FBS). The DCMP cost is only related to the size and density of the symmetric stiffness matrix, while the FBS cost depends on both the size and density of the stiffness matrix and the number of columns ofΦ 0 . As the frequency range of interest increases, more modes are needed to predict the structural response accurately. In such a case, although a single DCMP is needed in Equation (55) In this case, the total cost is dominated by the DCMP, and the CA method reduces the cost of one reanalysis considerably. However, if Φ 0 has 1050 modes, the cost of FBS increases to efficiency only when the number of retained modes is small. Otherwise, the computational savings do not compensate for the loss of accuracy from using K 0 (stiffness matrix of baseline design) instead of K p (stiffness matrix of new design). The modified combined approximations (MCA) method of this section addresses this issue.
The MCA method uses a subspace basis T whose columns are constructed using the recursive process instead of that of Equation (55). The selection of the appropriate number of basis vectors s is discussed later in this section. The only difference between Equations (55) and (59) is that matrix K 0 is inverted in the former while matrix K p is inverted in the latter. The DCMP of If the convergence is slow, multiple sets of updated mode shapes can be used so that For better accuracy, the above basis also includes the mode shapes Φ 0 of the baseline design.
Because the approximate modes T i are more accurate than the CA vectors R i in approximating the exact mode shapesΦ p , MCA can achieve similar accuracy to the CA method using fewer modes. The example of Section 4.3 demonstrates that the MCA method achieves good accuracy with only 1 basis vector whereas the CA method requires 3 to 6 basis vectors [13][14][15][16][17].
In summary, the proposed MCA method involves four steps in calculating the approximate eigenvectors Φ p as follows

• Calculate basis T using Equation (60) or Equation (61).
• Calculate the condensed stiffness and mass matrices K R and M R as • Solve the following reduced eigen-problem to calculate the eigenvalues and the projections of the modes in the reduced space spanned by T • Reconstruct the approximate eigenvectors Φ p as The slightly increased cost of using Equation (61)  All mode shapes in Equation (63) must be calculated simultaneously in order to ensure that the approximate mode shapes Φ p are orthogonal with respect to the mass and stiffness matrices. However, the cost of estimating the mode shapes Φ p using Equations (62) to (64) may increase quickly (quadratically) with the number of modes, and as a result, the MCA method may become more expensive than a direct eigen-solution when the number of dominant modes exceeds a certain limit. One way to circumvent this problem is to divide the frequency response into smaller frequency bands and calculate the frequency response in each band separately instead of solving for the frequency response in one step. The modal basis T in Equation (61) is divided into k groups as

Comparison of CA/MCA and PROM Methods
As we have discussed, a large overhead cost which is proportional to the number of design parameters is required before the PROM reanalysis is carried out. However, the CA/MCA method does not require this overhead cost because it does not need the basis P of Equation  (57) and (62). This is the case in gradient-free optimization problems employing a Genetic Algorithm for example, and in simulation-based probabilistic analysis problems employing the Monte-Carlo method. For these problems, the PROM method is more suitable because the subspace basis P does not change for every new design point. Table 1 summarizes the main characteristics, advantages and application areas of the CA/MCA and PROM methods.

CA/MCA Method PROM Method
Overhead Cost None Required cost to construct P.
Cost proportional to the number of design parameters m.

Basis Vector Variable basis R/T .
Size proportional to n and s.

Constant basis P.
Size proportional to n and m.

Reanalysis Cost Moderate
Relatively small size of R/T .

Must recalculate R/T at every new design.
High if no parametric relationship exists due to the condensation of large size and dense P.
Low if a parametric relationship exists.
Best Application Small number of reanalyses compared to the number of design parameters.
Evaluation of few design alternatives and gradient-based optimization.
Very large number of reanalyses.
Gradient-free optimization (e.g. genetic algorithms) and probabilistic analysis.

Reanalysis Methods in Dynamic Analysis and Optimization
The reanalysis methods of Section 3 can be used in different dynamic analyses such as modal or direct frequency response and free or forced vibration in time domain. Depending on the problem and the type of analysis, a particular reanalysis method may be preferred considering how many times it will be performed and how many design parameters will be allowed to change. This section demonstrates the computational efficiency and accuracy of reanalysis methods in dynamic analysis and optimization. It also introduces a new reanalysis method in Craig-Bampton substructuring with interface modes which is very useful for problems with many interface DOFs where FRF substructuring is not practically applicable.

Integration of MCA Method in Optimization
We have mentioned that the MCA method provides a good balance between accuracy and efficiency for problems that require a moderate number of reanalyses, as in gradient-based optimization. For problems where a large number of reanalyses is necessary, as in probabilistic analysis and gradient-free (e.g. genetic algorithms) optimization, a combination of the MCA and PROM methods is more suitable. Figure 6 shows a flowchart of the optimization process for modal frequency response problems. The DMAP (Direct Matrix Abstraction Program) capabilities in NASTRAN have been used to integrate the MCA method and the NASTRAN modal dynamic response and optimization (SOL 200). The highlighted boxes indicate modifications to the NASTRAN optimizer. Starting from the original design, the code first calculates the design parameter sensitivities in order to establish a local search direction and determine an improved design along the local direction. At the improved design, an eigen-solution is obtained to calculate a modal model and the corresponding modal response. The dynamic response at certain physical DOFs is then recovered from the modal response. At this point, a convergence check is performed to decide if the optimal design is obtained. If not, further iterations are needed and the above procedure is repeated. Many iterations are usually needed for practical problems to obtain the final optimal design. Section 4.3.2 demonstrates how this process was used to optimize the vibro-acoustic behavior of a 65,000 DOF, finite-element model of a truck. Using the MCA method, the computational cost of the entire optimization process was reduced in half compared with the existing NASTRAN approach.
As for a stand alone modal frequency response, the eigen-solution accounts for a large part of the overall optimization cost for vibratory problems where a modal model is used. A reanalysis method can be inserted into the procedure as shown in Figure 6 to provide an approximate eigen-solution saving therefore, substantial computational cost. Other reanalysis methods such as the CDH/VAO, CA or PROM can also be used depending on the number of design variables and the number of expected iterations.

Integration of MCA and PROM Methods
The • Perform an exact eigen-analysis at the baseline design point p 0 all parameters are at their lower limit, to obtain the baseline mode shapesΦ 0 .
• Use the MCA method at design point p i with all parameters at their low limit except the i th parameter which is set at its upper limit. Obtain approximate mode shapes for the i th corner point using the following recursive process • Form the subspace basis T as

Advances in Vibration Engineering and Structural Dynamics
where m is the total number of parameters.
• Obtain the approximate mode shapes Φ p using the subspace projection procedure of Equations (50) through (52) where T is used instead ofP.
The modal basis Φ p can be subsequently used in a modal dynamic response solution. Only step 4 is repeated in reanalysis. The computational savings can be substantial especially for problems where many reanalyses are needed.

Combined MCA and PROM Methods: Vibro-Acoustic Response of a Vehicle
The pickup truck vehicle model with 65,000 DOFs of Figure 7 is used in this section to demonstrate the advantages of the combined MCA and PROM method in optimizing the vibroacoustic response of a vehicle. The model has 78 components and roughly 11,000 nodes and elements. The example is performed on a SUN ULTRA workstation using NASTRAN v2001. The MCA and PROM methods have been implemented in NASTRAN DMAP. The sound pressure level at the driver's ear location is calculated using a vibro-acoustic analysis. The structural forced vibration response due to unit harmonic forces in x, y, and z directions at the engine mount locations, is coupled with an interior acoustic analysis. The first and second eigen-frequencies of the acoustic volume inside the cabin are 95.9 Hz and 128.3 Hz. The sound pressure level is calculated in the 80 to 140 Hz frequency range. The structure and fluid domains are coupled through boundary conditions ensuring continuity of vibratory displacement and acoustic pressure. A finite-element formulation of the coupled undamped problem yields the following system equations of motion [24].
where the vibratory displacement d S and the acoustic pressure p F are the primary variables. The spatial coupling matrix H SF indicates coupling between the two domains which is usually referred to as "two-way coupling." Due to this coupling, the combined structural-acoustic system of equations is not symmetric. If the acoustic effect on the structural response is small, the coupling term can be omitted, resulting in the so-called "one-way coupling," where the structural response is first calculated and then used as input ( f q in Equation 69) to solve for the acoustic response. The coupled structure-acoustic system can be solved either by a direct method, or more efficiently by a modal response method which can be applied to both the structural and acoustic domains.

Combined MCA and PROM Methods
To demonstrate the computational effectiveness and accuracy of integrating MCA in PROM, a reanalysis was performed for a modified design where five plate thickness parameters vary; chassis and its cross links, cabin, truck bed, left door, and right door. All parameters were increased by 100% from their baseline values. The sound pressure at the driver's ear was calculated using "two-way" coupling. A structural modal frequency response was used. The acoustic response was calculated using a direct method because the size of the acoustic model is relatively small. For the structural analysis, 1050 modes were retained in the 0 to 300 Hz frequency range. The combined MCA and PROM approach was compared against the NASTRAN direct solution for a modified design where all five parameters were at their upper limits. Only one iteration was used in Equation (59) in order to get the set of once-updated mode shapes for each corner design point. The subspace basis, which includes information for all five design parameters, is therefore, represented by The maximum error in natural frequencies as predicted by the combined MCA and PROM method and NASTRAN, is less than 0.45% in the entire frequency range. Figure 8 indicates that the sound pressures calculated by both methods are almost identical. The computational effort for the MCA method to obtain approximate mode shapes at each corner design point is about 30 seconds. In contrast, it takes about 180 seconds for an exact eigen-solution using NASTRAN. The computational cost to construct the reduced basis (Pin PROM and T in PROM+MCA) is compared in Table 2. The total cost was reduced from 1080 seconds to 330 seconds. The computational saving is more significant if the number of design parameter increases.

Optimization using MCA Method
The goal here is to minimize the sound pressure at the driver's ear. A total of 41 design parameters are used representing the thickness of all vehicle components modeled with plate elements. All thicknesses are allowed to change by 100% from their baseline values. Table 3 describes all design parameters. At the initial point of the optimization process, all parameters are at their low bound.   Because of the large number of design parameters, the combined MCA and PROM approach Section 4.3.1 is not computationally efficient because the size of the PROM basis is very large (see Equation 70). For this reason, we use the MCA reanalysis method and demonstrate its capability to handle a large number of parameters. It approximates the mode shapes at intermediate design points using only T 1 in Equation (59). The subspace basis at each optimization step is thusT = Φ o T 1 . Because 1050 modes exist in the frequency range of 0 to 300 Hz of the initial design, the size of the MCA modal basis is 2*1050 = 2100.  The cost of solving for 1050 modes directly from NASTRAN is 180 seconds (see Table 2). In the MCA method, the cost of solving the linear system of equations in Equation (59) is 31 seconds, and the additional combined cost of Equations (62) to (64) is 373 seconds, resulting in a total cost of 404 seconds (see Table 4). To reduce this cost, the 1050 modes are divided into 21 groups and the modes in each group are obtained separately as explained in the last paragraph of Section 3.4. This reduces the cost of Equations (62) to (64) to 66 seconds for a total cost of 97 seconds, which is about half the cost of the direct NASTRAN method.  The gradient-based optimizer in NASTRAN (SOL 200) using the optimization process of Figure 6 needed three iterations to calculate the optimal design. Figure 9 compares the Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402 163 sound pressure at the driver's ear between the optimal and initial designs. Figure 10 shows the percentage increase of optimal values relative to the initial values for all 41 design parameters. In the frequency range of 80-140Hz, the maximum sound pressure is slightly reduced from 7.9E-7 to 7.2E-7 Pascal. Most parameters are minimally changed. The largest increase is 20% for the rail mount thickness (parameter #27). Figure 11. Comparison of sound pressure at driver's ear between baseline and optimal designs with 20 initial populations and 4 generations.

Figure 12.
Comparison of sound pressure at driver's ear between baseline and optimal designs with 100 initial populations and 6 generations.

Advances in Vibration Engineering and Structural Dynamics
The Sequential Quadratic Programming (SQP) algorithm of NASTRAN can only find a local optimum. To obtain a more significant design improvement, two additional studies were performed using a Genetic Algorithm with the MCA method. The first study used 20 initial populations and 4 generations, and the second study used 100 initial populations and 6 generations. Figures 11and 12 show that the number of initial populations and the number of generations, affect the optimization results. While a higher number of initial populations and generations results in a slightly better result, both studies produced a much better optimum than the SQP algorithm. In the case of 100 initial populations and 6 generations, the sound pressure is reduced from 7.9E-7 Pascal to 2.0E-7 Pascal, which is equivalent to about 15 dB in sound pressure level (SPL).
To verify the accuracy of the MCA approximation, the sound pressure response at the optimal design from MCA+GA with 100 initial populations and 6 generations was evaluated by both direct NASTRAN and MCA. Figure 13 shows that the MCA method is very accurate. For a similar to MCA accuracy, the original CA method needed three sets of mode shapes to form the subspace basis, requiring 90 seconds to solve the linear equations. The much larger mode basis R in CA increases the computational cost to calculate the triple matrix products of Equation (57). Therefore for large scale, finite-element models with a high modal density, the proposed MCA method can be more efficient compared to either a complete NASTRAN analysis or the original CA method.

Reanalysis in Craig-Bampton Substructuring with Interface Modes
The FRF substructuring of Section 2.3 couples two structures using FRF information between the coupling (interface) DOFs, and the excitation and/or response DOFs. Although this approach is very efficient, it is practical only if we have a few coupling DOFs; e.g. connection of a vehicle suspension to chassis or connection of the exhaust system to body through a few hangers. If the physical substructures have interfaces with many DOFs, a different reduced-order modeling (ROM) approach must be used such as the Craig-Bampton ROM of Section 2.2.1. The Craig-Bampton ROM can be large however, if the number of retained interface DOFs is large. We address this problem by performing a secondary eigenvalue analysis which yields the so-called interface modes (see Section 2.2.2). The following section describes a reanalysis methodology for physical substructuring with Craig-Bampton ROMs using interface modes. We show that its accuracy is very good and the computational savings are substantial.

Craig-Bampton with Interface Modes and Reanalysis
In the Craig-Bampton CMS method (Craig-Bampton reduced-order model or CBROM), the mass and stiffness matrices of each substructure are partitioned into interface sub-matrices, interior (omitted DOF) sub-matrices, and their coupling sub-matrices. The dynamics of a structure are then described by the normal modes of its individual components, plus a set of modes called constraint modes that couple the components. In CBROM, there is no size reduction for constraint modes since all of them are kept in the reduced equations. If the finite element mesh is sufficiently fine, the constraint-mode DOFs will dominate the size of CBROM mass and stiffness matrices and result in a large computational cost. This issue is addressed by using interface modes (also called characteristic constraint -CC-modes). For that, a secondary eigenvalue analysis is performed using the constraint-mode The matrices from all substructures are assembled into a global CBROM with substructures coupled at interfaces by enforcing displacement compatibility. If k i C and m i C are the component (substructure) matrices, the global matrices K C and M C are assembled as and a secondary eigenvalue analysis is performed to calculate the interface modes Φ CC as The matrices in Equations (9), (10) and (12) are then reduced as where the matrices m i CC , m i CCN and k i CC are of much smaller size than matrices m i C , m i CN and k i C .
The interface modes reduce the interface size producing a smaller reduced order model (ROM) compared with the traditional Craig-Bampton ROM (CBROM). However, they are calculated from the assembled interface K and M matrices. Thus, the calculation of constraint modes and all matrix multiplications related to constraint modes are still necessary. The interface mode method reduces the size of ROM but it does not reduce the computational cost related to the constraint modes.
If the interface modes Φ CC were known before hand, the calculations in Equations (6), (9), (10) and (12) and Equation (72) could be performed much more efficiently as follows 55 The following observations can be made:

1.
In Equations (74) to (76), the computation involves Φ CC and Φ i C and does not involve Φ i C . Therefore, the calculation of original constraint modes Φ i C is no longer needed. (73), the number of columns of matrix (k i ΩΓ Φ CC ) is equal to the number of interface modes n cc which is usually smaller thann c . Therefore, the FBS cost of solving for Φ i C is proportional ton cc and it is much smaller than the FBS cost of solving forΦ i C .

3.
Because both Φ CC and Φ i C are of size n cc the cost of matrix multiplication and tripleproduct in Equations (74) to (76) are now proportional to n cc andn cc 2 . Therefore, the cost is much smaller than the corresponding cost in Equations (9), (10) and (12).
In this CCROM method which is based on CBROM, the interface modes Φ CC are obtained using the assembled interface partitions of the CBROM formulation. are applied to improve the computational efficiency.

A Car Door Example
The car door model of Figures 14 and 15 (9), (10) and (12). Figure 16 shows the interface nodes.
For the initial design using the CCROM method, 52 interface modes are calculated below 600 Hz. A modified design is created where the shell thicknesses for the outer door (substruc-    Figure 17 compares the natural frequencies of the new (modified) design between the original CCROM (Craig-Bampton with Interface modes) method and the new approach where reanalysis is used in CCROM to approximate the interface modes. We observe that the natural frequencies of the modified design are very different from those of the original design. Also, the accuracy of the proposed reanalysis method is excellent. The frequencies for the modified design calculated by the original CCROM and the proposed new approach are almost identical. The percentage error of the new approach versus the original CCROM approach is less than 1% on average. The computation cost is summarized in Table 5.

Advances in Vibration Engineering and Structural Dynamics
In the new approach to reduce the cost related to constraint modes, the total remaining cost is dominated by the cost of calculating the normal modes for each substructure. For example, the calculation of the normal modes for Substructure 2 took 110 seconds out of a total of 139 seconds (see Table 5). It should be noted that the normal modes cost can be further reduced by applying another reanalysis method such as CDH/VAO, CA or MCA to approximate the normal modes. Therefore, the overall cost of substructuring based on Craig-Bampton with interface modes, can be drastically reduced by using the proposed reanalysis to approximate the constraint modes and a CDH/VAO or MCA reanalysis to approximate the normal modes at a new design.

Optimization of a Vehicle Model
A detailed optimization study is presented using a large-scale FE model of a vehicle. For simplicity, we call it "BETA" car model. It is composed of approximately 7.1 million DOFs and 1.1 million elements. Figure 18 shows all modeling details. There are also some RBE2 and PBUSH elements which are used as connectors. The maximum displacement at each of the ten door locations is observed in the y direction (lateral direction -perpendicular to the door plane). The optimization problem is stated as follows: The optimal value of each of the fifteen design variables is calculated in order to minimize the maximum response among the ten locations on the doors while the mass of the vehicle remains less or equal to the mass of the initial (nominal) vehicle. The response is calculated in the 100 Hz to 200 Hz frequency range and a 3% structural damping is used. The optimization problem is numerically very challenging because of 1. the many local optima and 2. the computational cost of each dynamic analysis.
The former was handled by using a hybrid optimization algorithm which first explores the entire design space using a Niching Genetic Algorithm (GA) [25] and then switches to a gra- A proper fitness function which minimizes the maximum response among the ten door locations while satisfying the vehicle mass constraint is chosen as follows The ratio of the nominal maximum response over the actual maximum response is used so that the fitness value increases when the actual response is reduced. This ratio is multiplied by 1 + p * min ( c, 0 ) where p = 10 is a penalty value andc = 1 − Mass Mass Nominal . Thus, c is positive if Mass is less than Mass Nominal satisfying the constraint and the value of 1 + p * min ( c, 0 ) is equal to one.
Otherwise, c becomes negative if Mass is greater thanMass Nominal and the term 1 + p * min ( c, 0 ) assumes a large negative value which reduces the fitness value considerably. As a result, the GA optimizer always satisfies the mass constraint while maximizing the value of the fitness function. Figure 22 summarizes the optimization results by comparing the maximum door response between the optimal and initial designs. The optimizer determined that the maximum response occurs at location 9 (center of left front door of Figure 19) at approximately 105 Hz. Figure 23 shows that this represents a vehicle local mode involving motion of the doors only. At the optimal design the maximum response was reduced from the initial 10 -3 m to 0.47*10 -3 m ( Table 6).   Table 6 compares the value of each design variable between the initial (nominal) and final optimal designs. It also indicates that all designed variables were allowed to vary within a lower and upper bound. The values of the five door design variables changed considerably between the initial and optimal designs. This is expected because the optimizer tried to suppress the local door mode. The stiffness of the four engine mounts and the six exhaust supports also changed. Although we intuitively expect the stiffness of the engine mounts to change but not the stiffness of the exhaust supports, this is not the case in this example. Table 6 also indicates that at the optimum we not only reduced the maximum response from 10 -3 m to 0.47*10 -3 m but the vehicle mass was also reduced from the initial 55. 12   Considering that the computational cost for each function evaluation was 4 minutes, the total computational time was (359 function evaluations) * (4 minutes per evaluation) = 1436 minutes or 23.9 hours. This is acceptable considering the size and type of performed analysis. The computational cost, in terms of number of function evaluations, was kept low by coupling the Niching GA with a Lazy Learning metamodeling technique [26,27]. The latter estimates the value of the fitness function from existing values at close by designs without calculating the actual response. It uses an error measure to figure out if the estimation is accurate. The error is small if enough previous designs, for which the fitness value was evaluated, are close to the new design. In this case, the metamodel estimates the current fitness value without running an actual dynamic response. If the error is large, an actual response is calculated and the fitness value of this new design is added to the "pool" of previous designs the Lazy Learning metamodeling technique will use downstream.

Conclusions and Future Work
Reduced-order models and reanalysis methodologies were presented for accurate and efficient vibration analysis of large-scale, finite element models, and for efficient design optimization of structures for best vibratory response. The optimization is able to handle a large number of design variables and identify local and global optima.
For large FE models, it is common to solve for the system response through modal reduction in order to improve computational efficiency. An eigenanalysis is performed using the system stiffness and mass matrices and a modal model is formed which is then solved for the response. The computational cost can be also reduced using substructuring (or reduced-order modeling) methods. A modal reduction is applied to each substructure to obtain the component modes and the system level response is then obtained using component mode synthesis. In optimization of dynamic systems involving design changes (e.g. thicknesses, material properties, etc) the FEA analysis must be repeated many times in order to obtain the optimum design. Also in probabilistic analysis where parameter uncertainties are present, the FEA analysis must be repeated for a large number of sample points. In such cases, the computational cost is very high, if not prohibitive.
To drastically reduce the computational cost without compromising accuracy beyond an acceptable level, we developed and used various reanalysis methods in conjunction with reduced-order modeling, in optimization of vibratory systems. Reanalysis methods are intended to efficiently calculate the structural response of a modified structure without solving the complete set of modified analysis equations. We presented a variety of reanalysis methods including the CDH/VAO method, the Combined Approximations (CA) and Modified Combined Approximations (MCA) method, and the Parametric Reduced-Order Modeling (PROM) method. Their advantages and limitations were fully described and demonstrated with practical examples.
Future work will concentrate on developing reanalysis methodologies for shape and topology optimization of vibratory systems and extend the presented work in optimization under uncertainty where efficient deterministic reanalysis methods will be combined with efficient probabilistic reanalysis methods.