Open access peer-reviewed chapter

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

By Zissimos P. Mourelatos, Dimitris Angelis and John Skarakis

Submitted: March 3rd 2012Reviewed: July 8th 2012Published: October 2nd 2012

DOI: 10.5772/51402

Downloaded: 2447

1. 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-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-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-19]. Accurate results and significant computational savings have been reported. All reported studies on the CA method [12-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.

A large number of modes must be calculated and used in a dynamic analysis of a large finite-element model with a high modal density, even if a reduced-order modeling approach (Section 2) is used. In such a case, the implicit assumption of the CA method that the cost of 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.


2. Reduced-Order Modeling for Dynamic Analysis

Computational efficiency is of paramount importance in vibration analysis of large-scale, finite-element models. Reduced-order modeling (or substructuring) is a common approach to reduce the computational effort. Substructuring methods can be classified in “mathematical” and “physical” methods. The “mathematical” substructuring methods include the Automatic Multi-level Substructuring (AMLS) and the Automatic Component Mode Synthesis (ACMS) in NASTRAN. The “physical” substructuring methods include the well known fixed-interface Craig-Bampton method. Both the AMLS and ACMS methods use graph theory 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.

2.1. Modal Reduction

For an undamped structure with stiffness and mass matrices Kand Mrespectively, under the excitation force vectorF, the equations of motion (EOM) for frequency response are


where the displacement dis 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 ωiand corresponding eigenvectors (mode shapes) φiare first obtained. Then, the displacement dis approximated in the reduced space formed by the first ndominant modes as


where is the modal basis and Uis 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


The response dcan be recovered by solving Equation (3) for the modal response Uand projecting it back to the physical coordinates using Equation (2). If ωmaxis the maximum excitation frequency, it is common practice to retain the mode shapes in the frequency range of0÷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 [K+jωCω2M]d=Fand Equation (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.

2.2. 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 [2023]. 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 constraint-mode 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.

Figure 1.

Illustration of interface modes.

2.2.1. Craig-Bampton Fixed Interface CMS

This section provides the basics of Craig-Bampton method using the fixed-interface assumption. The method is commonly used in CMS algorithms. The finite-element model of the entire structure is partitioned into a group of substructures. The DOFs in each substructure are divided into interface (Γ) DOF and interior (Ω) DOF. The equations of motion for the ithsubstructure are then expressed as


The fixed-interface Craig-Bampton CMS method utilizes two sets of modes to represent the substructure motion; substructure normal (N) modesΦiN, and constraint (C) modesΦiC, where idenotes the ithsubstructure. The size reduction of the Craig-Bampton method comes from the truncation of the normal modes ΦiN=[φi1φi2φin]which are calculated by fixing all interface DOFs and solving the following eigenvalue problem


The static constraint modesΦiCare calculated by enforcing a set of static unit constraints to the interface DOFs as


The original physical DOFs and can be thus represented by the constraint-mode DOFs and the normal-mode DOFs as


Using the above Craig-Bampton transformation, the original EOM of Equation (7), can be expressed as


where the superscripts Cand Nare used to indicate partition related to static constraint mode DOFs and fixed-interface normal mode DOFs, respectively. The matrix partitions of Equation (8) are


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 is usually a full matrix. Therefore Equation (9) can be computationally expensive due to the triple-product involving constraint modes. The computational cost of the Craig-Bamtpon method is mostly related to

  1. solving for the normal modes,

  2. solving for the constraint modes, and

  3. the transformation calculation in Equation (9).

2.2.2. Craig-Bampton CMS with Interface Modes

In Craig-Bampton CMS, the matrices from all substructures are assembled into a global CBROM with substructures coupled at interfaces by enforcing displacement compatibility. This synthesis yields the modal displacement vector of the synthesized system to be partitioned as


where nssis the number of substructures in the global structure. The corresponding synthesized CMS mass and stiffness matrices are as follows


where the component modal matrices and are diagonalized and their sizes depend on the number of selected modes for the frequency range of interest. However, the number of constraint-mode DOFs, or the size of matrices and, is equal to the number of DOFs 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


The eigenvectors ψnare transformed into the finite-element DOFs for the ithcomponent structure using the following transformation




is a selected set of nCCinterface eigenvectors which are few compared to the number of the constraint-mode DOFs. The matrix maps the global (system) interface DOFs back to the local (subsystem i) DOF. The vectors in are referred to as the interface modesor characteristic constraint (CC) modes, because they represent the characteristic physical motion associated with the constraint modes. Relatively few CC-mode DOFs are used compared to the number of interface DOFs.

Finally, the CMS matrices are transformed using the CCmodes and the reduced-order CMS matrices are obtained similarly to Equations (18). Now, the unknown displacement vector dROM is partitioned as


where superscript CCindicates the partition associated with the CCmodes. The equations of motion of the reduced order CMS model (ROM) are expressed by


The mass matrix MROM, the stiffness matrix KROM, and the applied force vector fROM, are explicitly written as






2.2.3. Filtration of Constraint Modes

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.

Figure 2.

Illustration of a “filtered” constraint mode.

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


where φpqCis the pthelement of the qthconstraint mode. 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.

2.3. 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.

2.3.1. 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 Hare frequency dependent and complex if damping is present. A bold letter indicates a matrix or vector. According to Equation (30), HAC for example, represents the displacement XA of DOF A due to a unit force FC on DOF C. Sub2 is a finite element (FE) type substructure. Its dynamic behavior is described using the stiffness K, mass Mand damping Bmatrices.

Figure 3.

Two-substructure example and notation.

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 Hin 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 Bform the frequency dependent dynamic matrix Z=K+iωBω2Mwhich 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 Hmatrix of Sub1 with the Zmatrix of Sub2 and calculate (solve for) the system matrix His 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 Bindicates the interior DOFs of Sub2, and subscript Cindicates the connection/common/coupling DOFs between Sub2 and Sub1. Because of displacement compatibility at the interface, XCappears 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 FC at Sub1 and Sub2.

To couple Sub1 and Sub2, compatibility of forces at the interface is applied aswhere the force vector with is obtained from the second row of Equation (30) and the second row of Equation (31) provides the force vector. We thus have


From Equation (31),


where . Substitution of Equation (33) in Equation (32) yields


where . From Equation (34), XC can be expressed in terms of FA , FB and FC as


Substitution of Equation (35) in Equation (33) gives XB in terms of FA , FB and FC as


Solving Equation (30) for and substituting yields


We can now express XA in terms of FA , FB and FC by substituting Equation (35) in Equation (37), as


Based on Equations (38), (36) and (35), XA, XB and XC are expressed in terms of FA , FB and FA as follows


resulting in the following FRF of the assembled system


where S= [ΦΘI] and .

2.4.2. Algorithm for FRF/FRF Coupling

Figure 4 shows the coupling of two FRF type substructures. The equations of motion for Sub1 and Sub2 and are expressed as




Figure 4.

Two FRF type substructures example and notation

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


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

3. Reanalysis Methods for Dynamic Analysis

3.1. 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 K0 and M0, the exact mode shapes in Φ0 are obtained by solving the eigen-problem


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 ΔKand ΔMare 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 asKR=Φ0TKpΦ0and MR=Φ0TMpΦ0and the following reduced eigen-value 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.

3.2. 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 mdesign variables, Zhang [9] suggested that the representative designs include a baselinedesign for which all parameters are at their lower limits plus mdesigns 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.

Figure 5.

Design space of three parameters.

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 Pincludes 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Φ0is the modal matrix composed of the dominant mode shapes of the baseline design, and Φiis the modal matrix of the ithcorner point. The mode shapes of the new design satisfy the following eigenvalue problem,


where Λis a diagonal matrix of the first seigenvalues.

A reduced eigenvalue problem is obtained by pre-multiplying both sides of Equation (50) by PTas


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 KRandMR.

For mdesign variables, (m+1)eigenvalue problems must be solved in order to form the basis Pof Equation (49). Therefore, both the cost of obtaining the modal matrices Φiand the size of matrix Pincrease 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 mcorner points in the design space, and form subspace basisP.

  2. Calculate the reduced stiffness and mass matrices KRand MRfrom Equation (52).

  3. Solve eigenproblem (51) for matrixΘ.

  4. Reconstruct the approximated eigenvectors in Φ˜pusing Equation (48).

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

  1. the cost of performing (m+1)full eigen-analyses to form subspace basis Pin Equation (49), and

  2. the cost of reanalysis of each new design in steps 2 to 4.

The former is the fixed cost of PROM because it does not depend on the numbers of reanalyses and the latter is the variable cost of PROM because it is proportional to the number of reanalyses. The fixed cost is not attributed to the calculation of the response for a particular design. It is simply required to obtain the information needed to apply PROM. The variable cost (cost of reanalysis of a new design in part b) is small compared to the fixed cost. The fixed cost of PROM is proportional to the number of design variables mbecause it consists of the dominant eigenvectors Φ0of the baseline design and the dominant eigenvectors Φi,i=1,...,mof the mcorner design points (see Equation 49). When the size of basis Pincreases 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).

3.3. 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-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 K0and M0are the stiffness and mass matrices of the original (baseline) design, the exact mode shapes Φ0are obtained by solving the eigen-problemK0Φ0=λ0M0Φ0. We want to approximate the mode shapes of a modified design (subscript p) with stiffness and mass matrices Kp=K0+ΔKand Mp=M0+ΔMwhere ΔKand ΔMrepresent large perturbations. The CA method estimates the new eigenvalues λpand eigenvectors Φpwithout performing an exact eigenvalue analysis.

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 Rj=Φp,jΦp,j1to form a subspace basis to approximate the modes of the new design. In order to simplify the calculations, λpK01MpΦp,j1in Equation (53) is replaced with λpK01MpΦ0and Equation (54) becomes Φp,j=λpK01MpΦ0K01ΔKΦp,j1showing that the basis vectors satisfy the following recursive equation


where the first basis vector is assumed to beR1=K01MpΦ0.

The CA method forms a subspace basis


where sis usually between 3 and 6 [16-18, 23] and the mode shapes of the new (Kp,Mp) design are then approximated in the subspace spanned by Rusing the following algorithm:

  • Condense the stiffness and mass matrices as

  • Solve the reduced eigen-problem (using matrices KR and MR) to calculate the eigenvector matrixΘ.

  • Reconstruct the approximate eigenvectors of the new designΦ˜pas


The eigenvalues of the new design are approximated by the eigenvalues λ˜pof the reduced eigen-problem.

The CA method has three main advantages:

  1. it only requires a single matrix decomposition of stiffness matrix K0in Equation (55) to calculate the subspace basisR,

  2. it is accurate because the basis is updated for every new design, and

  3. the eigenvectors of a new design are efficiently approximated by Equation (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 Rdoes 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 Rand 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 mcorner 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 Ris not proportional to the number of parameters mas 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.

3.4. 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-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), the number of columns in Φ0may increase considerably, thereby increasing the cost of the repeated FBS. When the number of dominant modes becomes very large, the cost of performing the calculations in Equation (55) can be dominated by the FBS cost. For example, the vehicle model of Section 4.3 (Figure 7) has 65,000 degrees of freedom and 1050 modes in the frequency range of 0-300 Hz. The cost of one DCMP is 1.1 seconds (using a SUN ULTRA workstation and NASTRAN v2001) and the cost of one FBS is less than 0.1 seconds if Φ0has only one mode. In this case, the total cost is dominated by the DCMP, and the CA method reduces the cost of one reanalysis considerably. However, if Φ0has 1050 modes, the cost of FBS increases to 29 seconds dominating the cost of the DCMP. The CA method can therefore, improve the efficiency only when the number of retained modes is small. Otherwise, the computational savings do not compensate for the loss of accuracy from using K0(stiffness matrix of baseline design) instead of Kp(stiffness matrix of new design). The modified combined approximations (MCA) method of this section addresses this issue.

The MCA method uses a subspace basis Twhose columns are constructed using the recursive process


instead of that of Equation (55). The selection of the appropriate number of basis vectors sis discussed later in this section. The only difference between Equations (55) and (59) is that matrix K0is inverted in the former while matrix Kpis inverted in the latter. The DCMP of Kpmust be repeated for every new design. However, the cost of the repeated DCMP does not significantly increase the overall cost in Equation (59) because the latter is dominated by the FBS cost. The iterative process of Equation (59) provides a continuous mode shape updating of the new design. If the process converges in siterations, the mode shapes Tswill be the exact mode shapesΦp. Equation (55) does not have the same property. The vectors Tiprovide therefore, a more accurate approximation of the exact mode shapes Φpthan the Rivectors of the original CA method. This is an important advantage of MCA.

Because the mode shapes Tiin Equation (59) can quickly converge to the exact mode shapesΦp, for many practical problems only one iteration (i.e. s= 1) may be needed, resulting in


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 Φ0of the baseline design. Because the approximate modes Tiare more accurate than the CA vectors Riin 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-17].

In summary, the proposed MCA method involves four steps in calculating the approximate eigenvectors Φ˜pas follows

  • 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 Φ˜pas


The slightly increased cost of using Equation (61) instead of Equation (60) is usually smaller than the realized savings in steps 2 through 4 of Equations (62) through (64) due to the smaller size of the reduced basisT. The bases of Equations (60) and (61) are smaller in size than the CA basis of Equation (56) for comparable accuracy. The MCA method requires therefore, less computational effort for steps 2 through 4. The computational savings compensate for the increased cost of DCMP for dynamic reanalysis of large finite-element models with a large number of dominant modes.

All mode shapes in Equation (63) must be calculated simultaneously in order to ensure that the approximate mode shapes Φ˜pare orthogonal with respect to the mass and stiffness matrices. However, the cost of estimating the mode shapes Φ˜pusing 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 Tin Equation (61) is divided into kgroups as




Each group Ticontains roughly n/koriginal modes Φ0ifromΦ0, and their corresponding improved modes. The eigenvectors of the new design are calculated using Tiinstead of T in Equations (62) to (64). The process is repeated ktimes using a modal basis that is 1/kof the size of the original modal basis. All kgroups of eigenvectors are then collected together to calculate the frequency response of the new design. As demonstrated in Section 4.3.2, this reduces the cost considerably with minimal loss of accuracy.

3.5. 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 Pof Equation (49) (see Section 3.2). Therefore, the CA/MCA method is more suitable, when the number of reanalyses is comparable to the number of design parameters. This is usually true in gradient-based design optimization. The CA and MCA methods can become expensive however, when many reanalyses are needed because, for each reanalysis, they require a new basis Ror T(see Equations 56and 61, respectively) and new condensed mass and stiffness matrices in Equations (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 Pdoes 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 MethodPROM Method
Overhead CostNoneRequired cost to constructP.
Cost proportional to the number of design parametersm.
Basis VectorVariable basisR/T.
Size proportional tonands.
Constant basisP.
Size proportional tonandm.
Reanalysis CostModerate
Relatively small size ofR/T.
Must recalculateR/Tat every new design.
High if no parametric relationship exists due to the condensation of large size and denseP.
Low if a parametric relationship exists.
Best ApplicationSmall 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.

Table 1.

Comparison of the CA/MCA and PROM methods.

4. 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.

4.1. 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.

Figure 6.

Flowchart for mca-enhanced nastran optimization.

4.2. Integration of MCA and PROM Methods

The PROM method requires exact calculation of the mode shapes for all designs corresponding to the corner points of the parameter space in order to calculate the subspace basis Pof Equation (49). The required computational effort can be prohibitive for a large number of parameters (optimization design variables). This effort can be reduced substantially if the modes of each corner point are approximated by the MCA method. In this case, an exact eigen-solution is required only for the baseline design. The following steps describe an algorithm to integrate the MCA and PROM methods:

  • Perform an exact eigen-analysis at the baseline design point p0 all parameters are at their lower limit, to obtain the baseline mode shapesΦ0.

  • Use the MCA method at design point pi with all parameters at their low limit except the ithparameter which is set at its upper limit. Obtain approximate mode shapes for the ithcorner point using the following recursive process

  • Form the subspace basis Tas


where mis the total number of parameters.

  • Obtain the approximate mode shapes Φ˜pusing the subspace projection procedure of Equations (50) through (52) where Tis used instead ofP.

The modal basis Φ˜pcan 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.

4.3. 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 vibro-acoustic 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.

Figure 7.

FE model of a pickup truck.

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 dSand the acoustic pressure pFare the primary variables. The finite-element representation of the two domains consists of stiffness and mass matrix pairs KS,MSandKF,MF, respectively. The air density and wave speed are ρ0andc0. The right hand side of Equation (69) denotes the external forces.

The spatial coupling matrix HSFindicates 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 (fqin 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.

4.3.1. 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 Tin 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.

Figure 8.

Comparison of sound pressure at driver’s ear between combined MCA and PROM method and NASTRAN.

MethodSolving for mode shapeΦ0at baseline designSolving for mode shapes at 5 corner design pointsTotal Cost
PROM180 sec180*5=900 sec1080 sec
PROM+MCA180 sec30*5=150 sec330 sec

Table 2.

CPU time to construct reduced basis.

4.3.2. 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.

Prm.#Description(thickness of)Prm.#Description(thickness of)Prm.#Description(thickness of)
1Bumper15Radiator mtg.29Tire, front right
2Rails16Radiator mtg., mid.30Tire, rear left
3A-arm, low left17Fan cover, low31Tire, rear right
4A-arm, low right18Fan cover, up32Engine outer
5A-arm, up left19Cabin33A-arm conn., up left
6A-arm, up right20Cabin mtg. reinf.34A-arm conn., up right
7Tire rim21Door, left35A-arm conn., low left
8Engine Oil-box22Door, right36A-arm conn., low right
9Fan23Bed37Glass, left
10Hood24Brake, front left38Glass, right
11Fender, left25Brake, front right39Glass, rear
12Fender, right26Rail conn., rear40Glass, front
13Wheel house, left27Rail mount41Rail conn., front
14Wheel house, right28Tire, front left

Table 3.

Description of design parameters.

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 T1in Equation (59). The subspace basis at each optimization step is thusT=[ΦoT1]. 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.

Eq. (59)31 sec31 sec
Eq. (62)258 sec50 sec
Eq. (63)48 sec6 sec
Eq. (64)67 sec10 sec
Total Cost404 sec97 sec

Table 4.

CPU time of the MCA method.

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.

Figure 9.

Comparison of sound pressure at driver’s ear between initial and optimal designs.

Figure 10.

Percent increase of optimal design parameters relative to baseline design parameters.

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 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.

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 Rin 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.

Figure 13.

Comparison of sound pressure at driver’s ear between direct nastran and mca.

4.4. 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.

4.4.1. 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 modesthat 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 partitions of the CMS mass and stiffness matrices. The CC modes are the resultant eigenvectors. Details are provided in Sections 2.2.1 and 2.2.2.

The number of constraint modes ncequals to the number of interface DOF. For many FE models of large structures, the number of interface DOF can be rather large. The calculation of constraint modes in Equation (6) involves a decomposition step and a FBS step. The cost of FBS is proportional tonc. For any matrix multiplication that involvesΦiC, the cost is proportional tonc. For any triple-product that involves ΦiCthe cost is proportional tonc2.

The matrices from all substructures are assembled into a global CBROM with substructures coupled at interfaces by enforcing displacement compatibility. If kiCand miCare the component (substructure) matrices, the global matrices KCand MCare assembled as


and a secondary eigenvalue analysis is performed to calculate the interface modesΦCCas


The matrices in Equations (9), (10) and (12) are then reduced as


where the matrices miCC , miCCN and kiCC are of much smaller size than matrices miC , miCN and kiC .

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 Kand Mmatrices. 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 ΦCCwere known before hand, the calculations in Equations (6), (9), (10) and (12) and Equation (72) could be performed much more efficiently as follows55


The following observations can be made:

  1. In Equations (74) to (76), the computation involves ΦCCand Φ^iCand does not involveΦiC. Therefore, the calculation of original constraint modes ΦiCis no longer needed.

  2. In Equation (73), the number of columns of matrix (kiΩΓΦCC)is equal to the number of interface modes nccwhich is usually smaller thannc. Therefore, the FBS cost of solving for Φ^iCis proportional tonccand it is much smaller than the FBS cost of solving forΦiC.

  3. Because both ΦCCand Φ^iCare of size nccthe cost of matrix multiplication and triple-product in Equations (74) to (76) are now proportional to nccandncc2. Therefore, the cost is much smaller than the corresponding cost in Equations (9), (10) and (12).

In this CCROMmethod which is based on CBROM, the interface modes ΦCCare obtained using the assembled interface partitions of the CBROM formulation. Thus, it is impossible to know ΦCCbefore hand for a new design. For this reason, Equations (73) to (76) can not be theoretically implemented to improve efficiency. For this reason, we propose a reanalysis approach where the calculated interface modesΦCCfor original (baseline) design can be used as an approximation of the new interface modes at any modified design. In this case, Equations (73) to (76) are applied to improve the computational efficiency.

4.4.2. A Car Door Example

The car door model of Figures 14 and 15 is used to demonstrate the proposed reanalysis method for substructuring with Craig-Bampton method using interface modes. It has 25,800 nodes and 25,300 elements and is divided into two substructures. The first substructure includes the outer door shell and a bar attached to it. The second substructure includes the rest of the door. There are 293 nodes (1758 DOFs) on the interface. Therefore, the CBROMor CCROMmethod must calculate 1758 constraint modes according to Equation (6) for both substructures. The 1758 constraint modes are involved in matrix multiplication or triple-products in Equations (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 (substructure 1) and inner door (substructure 2) are doubled. To provide baseline numbers, the CCROMmethod is used on the new design to solve for the system natural frequencies. The new reanalysis approach is used on the new design to calculate approximate natural frequencies which are then compared with the baseline numbers. The interface modes calculated at the original design are used as an initial guess for the interface modes of the new design.

Figure 14.

Outside and inside views of car door model.

Figure 15.

Substructure 1 (outer door shell) and substructure 2 (rest of door).

Figure 16.

Interface nodes indicated by white dots.

Figure 17.

Comparison of natural frequencies between original CCROM method and CCROM with reanalysis for the car door example.

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.

Substructure 1:
CPU TimeNormal ModesConstraint ModesMultiplicationOther CostTotal Cost
CCROM8 sec61 sec65 sec3 sec137 sec
New Approach7 sec2 sec0.3 sec2 sec11 sec
Substructure 2:
CPU TimeNormal ModesConstraint ModesMultiplicationOther CostTotal Cost
CCROM108 sec282 sec927 sec10 sec1327 sec
New Approach110 sec16 sec3 sec10 sec139 sec

Table 5.

Summary of computational cost for the car door example.

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.

5. 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.

Figure 18.

Details of “BETA” car model.

Figure 19.

Ten response locations on two front doors.

We form an optimization problem in terms of the maximum vibratory displacement at any location of the outer shell of the two front doors by minimizing the maximum displacement among ten locations of the two front doors (Figure 19) due to a hypothetical engine excitation in the vertical (up-down) direction. The engine is represented by a lumped mass connected rigidly to the engine mounts (Figure 20). The powertrain-exhaust model has about 1.3 million DOFs and is composed of 29 PSHELL components and 12 PSOLID components. 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).

Figure 20.

Description of the fifteen design variables.

Figure 21.

Description of the five design variables on the doors.

Fifteen design variables are chosen; five structural elements of each door (thickness of door shell, front frame, rear frame, top panel, middle pipe), vertical stiffness of each of the four engine mounts, and vertical stiffness of each of the six exhaust system supports. All design variables are schematically indicated in Figures 20 and 21.

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 gradient-based optimizer (fmincon in MATLAB) using the best estimate of the optimal point from the GA as initial point. This ensures a rapid convergence to the final optimum because although all GA optimizers can move quickly to the vicinity of the final optimum, they have a very slow convergence rate in pinpointing the final optimum.

FRF substructuring is used to assemble all components of the vehicle (body, doors, and engine-exhaust) into a small reduced-order model. This keeps the computational cost of each dynamic analysis low (4 minutes per analysis). A modal model is created only once for the body subsystem and then used to generate an FRF representation. This modal model does not change during the optimization because the chosen design variables are not associated with the body. However, the modal models of the doors change during the optimization. The final model for the entire vehicle is created by assembling the FRF models of each component. The FRF assembly operation is repeated during optimization because the FRF models of the two doors keep changing.

The Niching GA optimizer maximizes a fitness function by modifying all design variables. 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=1MassMassNominal. Thus, c is positive if Massis less than MassNominalsatisfying the constraint and the value of 1+p*min(c,0)is equal to one.

Otherwise, c becomes negative if Massis greater thanMassNominaland 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).

Figure 22.

Comparison of optimal and initial designs.

Figure 23.

Vehicle local mode at 105 Hz indicating door deformation.

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 units to the final 51.92 units.

X1Door Shell0.110.70.6638
X2Front Frame0.110.70.3084
X3Rear Frame0.110.70.2019
X4Top Panel0.120.70.2019
X5Middel Pipe0.542.41.3722
X6Engine mount19.5370.5195110.6
X7Engine mount19.5370.5195161.8
X8Engine mount19.5370.5195108.8
X9Engine mount19.5370.5195132.1
X10Exhaust support19.5370.5195213.9
X11Exhaust support19.5370.5195218.1
X12Exhaust support19.5370.5195345.9
X13Exhaust support19.5370.5195286.1
X14Exhaust support19.5370.5195118.7
X15Exhaust support19.5370.5195233.4
Max Resp.1*10-30.47*10-3
Door Mass55.1251.92

Table 6.

Summary of optimal design.

Figure 24 shows the actual function evaluations (design points where the vehicle dynamic response was calculated) in the X1-X2-X3 space and indicates the vicinity of the optimal design point. The GA optimizer needed only 359 function evaluations and used a population size of 5*(15+1) = 80 and a maximum of 10 generations. The population size and the number of allowed generations were kept at a minimum in order to locate the vicinity of the optimum quickly without “wasting” valuable computational effort.


Figure 24.

Function evaluations of the Niching GA in the X1-X2-X3 space.

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.


6. 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.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zissimos P. Mourelatos, Dimitris Angelis and John Skarakis (October 2nd 2012). Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods, Advances in Vibration Engineering and Structural Dynamics, Francisco Beltran-Carbajal, IntechOpen, DOI: 10.5772/51402. Available from:

chapter statistics

2447total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Vibration Analysis of Cracked Beams Using the Finite Element Method

By A. S. Bouboulas, S. K. Georgantzinos and N. K. Anifantis

Related Book

First chapter

Adaptive Tuned Vibration Absorbers: Design Principles, Concepts and Physical Implementation

By Philip Bonello

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us