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

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 $$K$$ and $$M$$ respectively, under the excitation force vector $$F$$ , the equations of motion (EOM) for frequency response are

$$[K-{\omega }^{2}M]d=F$$ E1

where the displacement d is calculated at the forcing frequency $$\omega$$ . 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 $${\omega }_i$$ and corresponding eigenvectors (mode shapes) $${\phi }_i$$ are first obtained. Then, the displacement d is approximated in the reduced space formed by the first $$n$$ dominant modes as

E2

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

E3

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 $${\omega }_{\max }$$ is the maximum excitation frequency, it is common practice to retain the mode shapes in the frequency range of $$0÷2{\omega }_{\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 matrix $$C$$ , Equation (1) becomes $$\left[K+j\omega C-{\omega }^{2}M\right]d=F$$ and Equation (3) is modified as by adding the modal damping term . For proportional (structural or Rayleigh) damping, $$\left[C\right]$$ is a linear combination of $$\left[M\right]$$ and $$\left[K\right]$$ ; i.e. $$\left[C\right]=\alpha \left[M\right]+\beta \left[K\right]$$ where $$\alpha$$ and $$\beta$$ 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 [20–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 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 ( $$\Gamma$$ ) DOF and interior ( $$\Omega$$ ) DOF. The equations of motion for the i ^th substructure are then expressed as

E4

The fixed-interface Craig-Bampton CMS method utilizes two sets of modes to represent the substructure motion; substructure normal (N) modes $${\mathbf{\Phi }}_i^N$$ , and constraint (C) modes $${\Phi }_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 $${\Phi }_i^N=\left[\begin{array}{cccc}{\phi }_{i1}& {\phi }_{i2}& \cdots & {\phi }_{in}\end{array}\right]$$ which are calculated by fixing all interface DOFs and solving the following eigenvalue problem

$$\begin{array}{l}{k}_i^{\Omega \Omega }{}^{-1}{\Phi }_i^C={\Lambda }^N{m}_i^{\Omega \Omega }{\Phi }_i^C\\ \left[k_i^{\Omega \Omega }\right]\left\{{\phi }_{in}^{}\right\}={\lambda }_n\left[m_i^{\Omega \Omega }\right]\left\{{\phi }_{in}^{}\right\}forn=1,2,\mathrm{...}\end{array}$$ E5

The static constraint modes $${\Phi }_i^C$$ are calculated by enforcing a set of static unit constraints to the interface DOFs as

E6

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

$$\left\{\begin{array}{c}{u}_i^{\Gamma }\\ u_i^{\Omega }\end{array}\right\}=\left[\begin{array}{cc}I& 0\\ {\Phi }_i^C& {\Phi }_i^N\end{array}\right]\left\{\begin{array}{c}{u}_i^C\equiv u_i^{\Gamma }\\ u_i^N\end{array}\right\}$$ E7

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

$$\left[\begin{array}{cc}{m}_i^{CC}& m_i^{CN}\\ m_i^{NC}& m_i^{NN}\end{array}\right]\left\{\begin{array}{c}{\ddot{u}}_i^C\\ {\ddot{u}}_i^N\end{array}\right\}+\left[\begin{array}{cc}{k}_i^{CC}& k_i^{CN}\\ k_i^{NC}& k_i^{NN}\end{array}\right]\left\{\begin{array}{c}{u}_i^C\\ u_i^N\end{array}\right\}=\left\{\begin{array}{c}{f}_i^C\\ f_i^N\end{array}\right\}$$ E8

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

E9

E10

E11

E12

E13

E14

E15

E16

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

$$d^{CMS}={\left[\begin{array}{ccccc}{d}^C{}^T& d_1^N{}^T& d_2^N{}^T& \cdots & d_{n^{ss}}^N{}^T\end{array}\right]}^T$$ E17

where $$n^{ss}$$ is the number of substructures in the global structure. The corresponding synthesized CMS mass and stiffness matrices are as follows

$$\begin{array}{l}{M}^{CMS}=\left[\begin{array}{ccccc}{\overline{m}}^C& m_1^{CN}& m_2^{CN}& \cdots & m_{n^{SS}}^{CN}\\ m_1^{CN}{}^T& m_1^N& 0& \cdots & 0\\ m_2^{CN}{}^T& 0& m_2^N& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ m_{n^{SS}}^{CN}{}^T& 0& 0& \cdots & m_{n^{SS}}^N\end{array}\right]\\ K^{CMS}=\left[\begin{array}{ccccc}{\overline{k}}^C& 0& 0& \cdots & 0\\ 0& k_1^N& 0& \cdots & 0\\ 0& 0& k_2^N& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & k_{n^{SS}}^N\end{array}\right]\end{array}$$ E18

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

$${\overline{k}}^C{\psi }_n={\Lambda }_n{\overline{m}}^C{\psi }_{n}forn=1,2,3,\mathrm{...}$$ E19

The eigenvectors $${\psi }_n$$ are transformed into the finite-element DOFs for the i ^th component structure using the following transformation

E20

where

$$\Psi =\left[\begin{array}{cccc}{\psi }_1& {\psi }_2& \cdots & {\psi }_{n^{CC}}\end{array}\right]$$ E21

is a selected set of $$n^{CC}$$ interface 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 modes or 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 CC modes and the reduced-order CMS matrices are obtained similarly to Equations (18). Now, the unknown displacement vector d^ROM is partitioned as

$$d^{ROM}={\left[\begin{array}{ccccc}{d}^{CC}{}^T& d_1^N{}^T& d_2^N{}^T& \cdots & d_{n^{ss}}^N{}^T\end{array}\right]}^T$$ E22

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

$$\left[-{\omega }^2{M}^{ROM}+K^{ROM}\right]d^{ROM}=f^{ROM}$$ E23

The mass matrix M^ROM, the stiffness matrix K^ROM, and the applied force vector f^ROM, are explicitly written as

$$M^{ROM}=\left[\begin{array}{ccccc}{\overline{m}}^{CC}& {\overline{m}}_1^{CN}& {\overline{m}}_2^{CN}& \cdots & {\overline{m}}_{n^{SS}}^{CN}\\ {\overline{m}}_1^{CN}{}^T& m_1^N& 0& \cdots & 0\\ {\overline{m}}_2^{CN}{}^T& 0& m_2^N& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\overline{m}}_{n^{SS}}^{CN}{}^T& 0& 0& \cdots & m_{n^{SS}}^N\end{array}\right]$$ E24

$$K^{ROM}=\left[\begin{array}{ccccc}{\overline{k}}^{CC}& 0& 0& \cdots & 0\\ 0& k_1^N& 0& \cdots & 0\\ 0& 0& k_2^N& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & k_{n^{SS}}^N\end{array}\right]$$ E25

$$f^{ROM}={\left[\begin{array}{ccccc}{f}^{CC}{}^T& f_1^N{}^T& f_2^N{}^T& \cdots & f_{n^{SS}}^N{}^T\end{array}\right]}^T$$ E26

where

E27

and

$${\overline{\mathbf{m}}}_i^{CN}={\mathbf{\Psi }}^T{\mathbf{\beta }}_i^{C^T}{\mathbf{m}}_i^{CN}$$ E28

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

$${\phi }_{pq}^C=0if\left|{\phi }_{pq}^C\right|<\epsilon *\underset{p}{\max }\underset{}{\left|{\phi }_{pq}^C\right|}$$ E29

where $${\phi }_{pq}^C$$ is the p ^th element of the q ^th constraint 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 $$\epsilon$$ , 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 $$\epsilon =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 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.

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 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 $$\mathbf{Z}=\mathbf{K}+i\omega \mathbf{B}-{\omega }^2\mathbf{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

$$\left[\begin{array}{c}{X}_A\\ X_C\end{array}\right]=\left[\begin{array}{cc}{H}_{AA}& H_{AC}\\ H_{CA}& H_{CC}\end{array}\right]\left[\begin{array}{c}{F}_A\\ F_C^1\end{array}\right]$$ E30

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

$$\left[\begin{array}{c}{F}_B\\ F_C^2\end{array}\right]=\left[\begin{array}{cc}{Z}_{BB}& Z_{BC}\\ Z_{CB}& Z_{CC}\end{array}\right]\left[\begin{array}{c}{X}_B\\ X_C\end{array}\right]$$ E31

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

E32

From Equation (31),

E33

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

E34

where . From Equation (34), X_C can be expressed in terms of F_{A ,} F_B and F_C as

E35

Substitution of Equation (35) in Equation (33) gives X_B in terms of F_{A ,} F_B and F_C as

E36

Solving Equation (30) for and substituting yields

E37

We can now express X_A in terms of F_{A ,} F_B and F_C by substituting Equation (35) in Equation (37), as

E38

Based on Equations (38), (36) and (35), X_A, X_B and X_C are expressed in terms of F_{A ,} F_B and F_A as follows

E39

resulting in the following FRF of the assembled system

E40

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

$$\left[\begin{array}{c}{X}_A\\ X_C\end{array}\right]=\left[\begin{array}{cc}{H}_{AA}& H_{AC}\\ H_{CA}& H_{CC}^1\end{array}\right]\left[\begin{array}{c}{F}_A\\ F_C^1\end{array}\right]$$ E41

and

$$\left[\begin{array}{c}{X}_B\\ X_C\end{array}\right]=\left[\begin{array}{cc}{H}_{BB}& H_{BC}\\ H_{CB}& H_{CC}^2\end{array}\right]\left[\begin{array}{c}{F}_B\\ F_C^2\end{array}\right]$$ E42

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

$$\left[\begin{array}{c}{X}_A\\ X_B\\ X_C\end{array}\right]=\left(\left[\begin{array}{ccc}{H}_{AA}& H_{AC}& \\ H_{CA}& H_{CC}^1& \\ & & H_{BB}\end{array}\right]+\left[\begin{array}{c}{H}_{AC}\\ H_{CC}^1\\ H_{BC}\end{array}\right]{\left[H_{CC}^1+H_{CC}^2\right]}^{-1}{\left[\begin{array}{c}{H}_{AC}\\ H_{CC}^1\\ H_{BC}\end{array}\right]}^T\right)\left[\begin{array}{c}{F}_A\\ F_C\\ F_B\end{array}\right]$$ E43

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 K₀ and M₀, the exact mode shapes in Φ₀ are obtained by solving the eigen-problem

E44

where Λ₀ is the diagonal matrix of the baseline eigenvalues. A new design (subscript p) has the following stiffness and mass matrices

$$\begin{array}{ccc}{K}_p=K_0+\Delta K& & M_p=M_0+\Delta M\end{array}$$ E45

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 Φ₀ of the initial design. The new stiffness and mass matrices are then condensed as $$K_R={\Phi }_0^T{K}_p{\Phi }_0$$ and $$M_R={\Phi }_0^T{M}_p{\Phi }_0$$ and the following reduced eigen-value problem is solved to calculate the eigen-vector $$\Theta$$

$$K_p\Theta =M_p\Theta {\Lambda }_p$$ E46

The approximate eigenvalues of the new design are given by Λ_p and the approximate eigenvectors Φ_p are

$${\Phi }_p=R\Theta$$ E47

where $$R={\Phi }_0$$ . Thus, the modal response of the modified structure can be easily obtained and the actual response can be recovered using the eigenvectors $${\Phi }_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 $$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.

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

$$\Phi \approx {\tilde{\Phi }}_p=P\Theta$$ E48

where the modal matrix P includes the basis vectors as in Equation (49) and $$\Theta$$ represents the participation factors of these vectors. The columns of P are the dominant mode shapes of the above $$\left(m+1\right)$$ designs,

$$P=\left[\begin{array}{cccc}{\Phi }_0& {\Phi }_1& \dots & {\Phi }_m\end{array}\right]$$ E49

where $${\Phi }_0$$ is the modal matrix composed of the dominant mode shapes of the baseline design, and $${\Phi }_i$$ is the modal matrix of the $$i^{th}$$ corner point. The mode shapes of the new design satisfy the following eigenvalue problem,

$$K{\tilde{\Phi }}_p=M{\tilde{\Phi }}_p\Lambda \iff KP\Theta =MP\Theta \Lambda$$ E50

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 $$P^T$$ as

$$K_R\Theta =M_R\Theta \Lambda$$ E51

where the reduced stiffness and mass matrices are

$$K_R=P^{T}KPand{M}_R=P^{T}MP$$ E52

Thus, the matrix $$\Theta$$ in Equation (48) consists of the eigenvectors of the reduced stiffness and mass matrices $$K_R$$ and $$M_R$$ .

For $$m$$ design variables, $$\left(m+1\right)$$ eigenvalue problems must be solved in order to form the basis $$P$$ of Equation (49). Therefore, both the cost of obtaining the modal matrices $${\Phi }_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 basis $$P$$ .

2. Calculate the reduced stiffness and mass matrices $$K_R$$ and $$M_R$$ from Equation (52).

3. Solve eigenproblem (51) for matrix $$\Theta$$ .

4. Reconstruct the approximated eigenvectors in $${\tilde{\Phi }}_p$$ using 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 $$\left(m+1\right)$$ full eigen-analyses to form subspace basis $$P$$ in 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 m because it consists of the dominant eigenvectors $${\Phi }_0$$ of the baseline design and the dominant eigenvectors $${\Phi }_i,i=1,\mathrm{...},m$$ of the m corner design points (see Equation 49). When the size of basis $$P$$ increases so does the fixed cost because more eigenvalue problems and mode shapes must be calculated to form basis $$P$$ . 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 $$K_0$$ and $$M_0$$ are the stiffness and mass matrices of the original (baseline) design, the exact mode shapes $${\Phi }_0$$ are obtained by solving the eigen-problem $$K_0{\Phi }_0={\lambda }_0{M}_0{\Phi }_0$$ . We want to approximate the mode shapes of a modified design (subscript p) with stiffness and mass matrices $$K_p=K_0+\Delta K$$ and $$M_p=M_0+\Delta M$$ where $$\Delta K$$ and $$\Delta M$$ represent large perturbations. The CA method estimates the new eigenvalues $${\lambda }_p$$ and eigenvectors $${\Phi }_p$$ without performing an exact eigenvalue analysis.

The eigen-problem for the modified design can be expressed as

$${\Phi }_p={\lambda }_p{K}_0^{-1}{M}_p{\Phi }_p-K_0^{-1}\Delta K{\Phi }_p$$ E53

leading to the following recursive equation

$${\Phi }_{p,j}={\lambda }_p{K}_0^{-1}{M}_p{\Phi }_{p,j-1}-K_0^{-1}\Delta K{\Phi }_{p,j-1}$$ E54

which produces a sequence of approximations of the mode shapes $${\Phi }_{p,j},j=1,2,\dots ,s$$ . The CA method uses the changes $$R_j={\Phi }_{p,j}-{\Phi }_{p,j-1}$$ to form a subspace basis to approximate the modes of the new design. In order to simplify the calculations, $${\lambda }_p{K}_0^{-1}{M}_p{\Phi }_{p,j-1}$$ in Equation (53) is replaced with $${\lambda }_p{K}_0^{-1}{M}_p{\Phi }_0$$ and Equation (54) becomes $${\Phi }_{p,j}={\lambda }_p{K}_0^{-1}{M}_p{\Phi }_0-K_0^{-1}\Delta K{\Phi }_{p,j-1}$$ showing that the basis vectors satisfy the following recursive equation

$$R_j=-K_0^{-1}\Delta K{R}_{j-1}j=2,\dots ,s$$ E55

where the first basis vector is assumed to be $$R_1=K_0^{-1}{M}_p{\Phi }_0$$ .

The CA method forms a subspace basis

$$R=\left[\begin{array}{cccc}{R}_1& R_2& \cdots & R_s\end{array}\right]$$ E56

where $$s$$ is usually between 3 and 6 [16-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

$$\begin{array}{ccc}{K}_R=R^T{K}_{p}R& & M_R=R^T{M}_{p}R\end{array}$$ E57

• Solve the reduced eigen-problem (using matrices K_R and M_R) to calculate the eigenvector matrix $$\Theta$$ .

• Reconstruct the approximate eigenvectors of the new design $${\tilde{\Phi }}_p$$ as

$${\tilde{\Phi }}_p=R\Theta$$ E58

The eigenvalues of the new design are approximated by the eigenvalues $${\tilde{\lambda }}_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) to calculate the subspace basis $$R$$ ,

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 $$\Theta$$ 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 basis $$P$$ .

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 $${\Phi }_i,i=1,\cdots ,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.

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 $${\Phi }_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 $${\Phi }_0$$ may 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 $${\Phi }_0$$ has 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 $${\Phi }_0$$ has 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 $$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

$$\begin{array}{l}{T}_1=K_p^{-1}\left(M_p{\Phi }_0\right)\\ \begin{array}{ccc}{T}_i=K_p^{-1}(M_p{T}_{i-1})& & i=2,3,\cdots ,s\end{array}\end{array}$$ E59

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 $$K_p$$ must 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 s iterations, the mode shapes $$T_s$$ will be the exact mode shapes $${\Phi }_p$$ . Equation (55) does not have the same property. The vectors $$T_i$$ provide therefore, a more accurate approximation of the exact mode shapes $${\Phi }_p$$ than the $$R_i$$ vectors of the original CA method. This is an important advantage of MCA.

Because the mode shapes $$T_i$$ in Equation (59) can quickly converge to the exact mode shapes $${\Phi }_p$$ , for many practical problems only one iteration (i.e. s = 1) may be needed, resulting in

$$T=T_1$$ E60

If the convergence is slow, multiple sets of updated mode shapes can be used so that

$$T=\left[\begin{array}{ccccc}{\Phi }_0& T_1& T_2& \cdots & T_s\end{array}\right]$$ E61

For better accuracy, the above basis also includes the mode shapes $${\Phi }_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 $${\Phi }_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 $${\tilde{\Phi }}_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

$$\begin{array}{ccc}{K}_R=T^T{K}_{p}T& & M_R=T^T{M}_{p}T\end{array}$$ E62

• Solve the following reduced eigen-problem to calculate the eigenvalues and the projections of the modes in the reduced space spanned by T

$$\left(K_R-\lambda M_R\right)\Theta =0$$ E63

• Reconstruct the approximate eigenvectors $${\tilde{\Phi }}_p$$ as

$${\tilde{\Phi }}_p=T\Theta$$ E64

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 basis $$T$$ . 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 $${\tilde{\Phi }}_p$$ are orthogonal with respect to the mass and stiffness matrices. However, the cost of estimating the mode shapes $${\tilde{\Phi }}_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

$$T=\left[\begin{array}{cccc}{T}^1& T^2& \cdots & T^k\end{array}\right]$$ E65

where

$$T^i=\left[\begin{array}{ccccc}{\Phi }_0^i& T_1^i& T_2^i& \cdots & T_s^i\end{array}\right]$$ E66

Each group $$T^i$$ contains roughly $$n/k$$ original modes $${\Phi }_0^i$$ from $${\Phi }_0$$ , and their corresponding improved modes. The eigenvectors of the new design are calculated using $$T^i$$ instead of T in Equations (62) to (64). The process is repeated $$k$$ times using a modal basis that is $$1/k$$ of the size of the original modal basis. All $$k$$ groups 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 $$P$$ of 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 $$R$$ or $$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 $$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.

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 $$P$$ of 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 p₀ all parameters are at their lower limit, to obtain the baseline mode shapes $${\Phi }_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

$$\begin{array}{l}{T}_{i,1}=K_i^{-1}\left(M_i{\Phi }_o\right)\\ \begin{array}{ccc}{T}_{i,j+1}=K_i^{-1}(M_i{T}_{i,j})& & j=2,3,\cdots ,s\end{array}\end{array}$$ E67

• Form the subspace basis T as

$$T=\left[\begin{array}{ccccc}{\Phi }_o& T_{0,s}& T_{1,s}& \cdots & T_{m,s}\end{array}\right]$$ E68

where $$m$$ is the total number of parameters.

• Obtain the approximate mode shapes $${\tilde{\Phi }}_p$$ using the subspace projection procedure of Equations (50) through (52) where $$T$$ is used instead of $$P$$ .

The modal basis $${\tilde{\Phi }}_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.

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

$$\left(\left[\begin{array}{cc}{K}_S& -H_{SF}\\ 0& K_F\end{array}\right]-j{\omega }^2\left[\begin{array}{cc}{M}_S& 0\\ {\rho }_0{c}_0^2{H}_{SF}^T& M_F\end{array}\right]\right)\left[\begin{array}{c}{d}_S\\ p_F\end{array}\right]=\left[\begin{array}{c}{f}_b\\ f_q\end{array}\right]$$ E69

where the vibratory displacement $$d_S$$ and the acoustic pressure $$p_F$$ are the primary variables. The finite-element representation of the two domains consists of stiffness and mass matrix pairs $$⟨K_S,M_S⟩$$ and $$⟨K_F,M_F⟩$$ , respectively. The air density and wave speed are $${\rho }_0$$ and $$c_0$$ . The right hand side of Equation (69) denotes the external forces.

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.

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

$$T=\left[\begin{array}{ccccc}{\Phi }_o& T_{0,1}& T_{1,1}& \cdots & T_{5,1}\end{array}\right]$$ E70

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 ( $$P$$ in 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.

Figure 8.

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

Method	Solving for mode shape $${\Phi }_0$$ at baseline design	Solving for mode shapes at 5 corner design points	Total Cost
PROM	180 sec	180*5=900 sec	1080 sec
PROM+MCA	180 sec	30*5=150 sec	330 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)
1	Bumper	15	Radiator mtg.	29	Tire, front right
2	Rails	16	Radiator mtg., mid.	30	Tire, rear left
3	A-arm, low left	17	Fan cover, low	31	Tire, rear right
4	A-arm, low right	18	Fan cover, up	32	Engine outer
5	A-arm, up left	19	Cabin	33	A-arm conn., up left
6	A-arm, up right	20	Cabin mtg. reinf.	34	A-arm conn., up right
7	Tire rim	21	Door, left	35	A-arm conn., low left
8	Engine Oil-box	22	Door, right	36	A-arm conn., low right
9	Fan	23	Bed	37	Glass, left
10	Hood	24	Brake, front left	38	Glass, right
11	Fender, left	25	Brake, front right	39	Glass, rear
12	Fender, right	26	Rail conn., rear	40	Glass, front
13	Wheel house, left	27	Rail mount	41	Rail conn., front
14	Wheel house, right	28	Tire, 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 $$T_1$$ in Equation (59). The subspace basis at each optimization step is thus $$T=\left[\begin{array}{cc}{\Phi }_o& T_1\end{array}\right]$$ . 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.

	k=1	k=21
Eq. (59)	31 sec	31 sec
Eq. (62)	258 sec	50 sec
Eq. (63)	48 sec	6 sec
Eq. (64)	67 sec	10 sec
Total Cost	404 sec	97 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 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.

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 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 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 $$n_c$$ equals 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 to $$n_c$$ . For any matrix multiplication that involves $${\Phi }_i^C$$ , the cost is proportional to $$n_c$$ . For any triple-product that involves $${\Phi }_i^C$$ the cost is proportional to $$n_c^2$$ .

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

$$K^C={\displaystyle \sum k_i^C},M^C={\displaystyle \sum m_i^C}$$ E71a

and a secondary eigenvalue analysis is performed to calculate the interface modes $${\Phi }^{CC}$$ as

$$\left[K^C-{\lambda }^{CC}{M}^C\right]{\Phi }^{CC}=0$$ E71b

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

$$\begin{array}{l}{m}_i^{CC}={\Phi }^{{CCT}^T}{}^{}{m}_i^C{\Phi }^{CC}\\ m_i^{CCN}={\Phi }^{{CCT}^T}{}^{}{m}_i^{CN}\\ k_i^{CC}={\Phi }^{{CCT}^T}{}^{}{k}_i^C{\Phi }^{CC}\end{array}$$ E72

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 $${\Phi }^{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⁵⁵

$${\widehat{\mathbf{\Phi }}}_i^C={\mathbf{\Phi }}_i^C{\mathbf{\Phi }}^{CC}=-{\mathbf{k}}_i^{\Omega {\Omega }^{-1}}({\mathbf{k}}_i^{\Omega \Gamma }{\mathbf{\Phi }}^{CC})$$ E73

$${\mathbf{m}}_i^{CC}={\left({\Phi }^{CC}\right)}^T{\mathbf{m}}_i^{\Gamma \Gamma }{\Phi }^{CC}+{\left({\Phi }^{CC}\right)}^T{\mathbf{m}}_i^{\Gamma \Omega }{\widehat{\Phi }}_i^C+{\left({\widehat{\Phi }}_i^C\right)}^T{\mathbf{m}}_i^{\Omega \Gamma }{\Phi }^{CC}+{\left({\widehat{\Phi }}_i^C\right)}^T{\mathbf{m}}_i^{\Omega \Omega }{\left({\widehat{\Phi }}_i^C\right)}^T$$ E74

$${\mathbf{m}}_i^{CCN}={\left({\Phi }^{CC}\right)}^T{\mathbf{m}}_i^{\Gamma \Omega }{\Phi }_i^N+{\left({\widehat{\Phi }}_i^C\right)}^T{\mathbf{m}}_i^{\Omega \Omega }{\Phi }_i^N$$ E75

$${\mathbf{k}}_i^{CC}={\left({\Phi }^{CC}\right)}^T{\mathbf{k}}_i^{\Gamma \Gamma }{\Phi }^{CC}-{\left({\Phi }^{CC}\right)}^T{k}_i^{\Gamma \Omega }{\widehat{\Phi }}_i^{CC}$$ E76

The following observations can be made:

1. In Equations (74) to (76), the computation involves $${\Phi }^{CC}$$ and $${\widehat{\Phi }}_i^C$$ and does not involve $${\Phi }_i^C$$ . Therefore, the calculation of original constraint modes $${\Phi }_i^C$$ is no longer needed.

2. In Equation (73), the number of columns of matrix $$({\mathbf{k}}_i^{\Omega \Gamma }{\Phi }^{CC})$$ is equal to the number of interface modes $$n_{cc}$$ which is usually smaller than $$n_c$$ . Therefore, the FBS cost of solving for $${\widehat{\Phi }}_i^C$$ is proportional to $$n_{cc}$$ and it is much smaller than the FBS cost of solving for $${\Phi }_i^C$$ .

3. Because both $${\Phi }^{CC}$$ and $${\widehat{\Phi }}_i^C$$ are of size $$n_{cc}$$ the cost of matrix multiplication and triple-product in Equations (74) to (76) are now proportional to $$n_{cc}$$ and $$n_{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 $${\Phi }^{CC}$$ are obtained using the assembled interface partitions of the CBROM formulation. Thus, it is impossible to know $${\Phi }^{CC}$$ before 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 $${\Phi }^{CC}$$ for 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 CBROM or CCROM method 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 CCROM method 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 Time	Normal Modes	Constraint Modes	Multiplication	Other Cost	Total Cost
CCROM	8 sec	61 sec	65 sec	3 sec	137 sec
New Approach	7 sec	2 sec	0.3 sec	2 sec	11 sec
Substructure 2:
CPU Time	Normal Modes	Constraint Modes	Multiplication	Other Cost	Total Cost
CCROM	108 sec	282 sec	927 sec	10 sec	1327 sec
New Approach	110 sec	16 sec	3 sec	10 sec	139 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

                                       $$Fitness=\frac{{\left[\underset{i=1}{\overset{10}{\max }}\left(\mathrm{Re}s{p}_i\right)\right]}_{No\min al}}{\left[\underset{i=1}{\overset{10}{\max }}\left(\mathrm{Re}s{p}_i\right)\right]}\ast \left[1+p\ast \min \left(c,0\right)\right]$$

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 \left(c,0\right)$$ where p = 10 is a penalty value and $$c=1-\frac{Mass}{Mas{s}_{No\min al}}$$ . Thus, c is positive if $$Mass$$ is less than $$Mas{s}_{No\min al}$$ satisfying the constraint and the value of $$1+p*\min \left(c,0\right)$$ is equal to one.

Otherwise, c becomes negative if $$Mass$$ is greater than $$Mas{s}_{No\min al}$$ and the term $$1+p*\min \left(c,0\right)$$ 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.

DesignVariables	Thickness	LowerBound	UpperBound	NominalDesign	OptimalDesign
X1	Door Shell	0.1	1	0.7	0.6638
X2	Front Frame	0.1	1	0.7	0.3084
X3	Rear Frame	0.1	1	0.7	0.2019
X4	Top Panel	0.1	2	0.7	0.2019
X5	Middel Pipe	0.5	4	2.4	1.3722
X6	Engine mount	19.5	370.5	195	110.6
X7	Engine mount	19.5	370.5	195	161.8
X8	Engine mount	19.5	370.5	195	108.8
X9	Engine mount	19.5	370.5	195	132.1
X10	Exhaust support	19.5	370.5	195	213.9
X11	Exhaust support	19.5	370.5	195	218.1
X12	Exhaust support	19.5	370.5	195	345.9
X13	Exhaust support	19.5	370.5	195	286.1
X14	Exhaust support	19.5	370.5	195	118.7
X15	Exhaust support	19.5	370.5	195	233.4
Max Resp.				1*10^-3	0.47*10^-3
Door Mass				55.12	51.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 X₁-X₂-X₃ 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 X₁-X₂-X₃ 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.