**2. Reduced-Order Modeling for Dynamic Analysis**

Computational efficiency is of paramount importance in vibration analysis of large-scale, fi‐ nite-element models. Reduced-order modeling (or substructuring) is a common approach to reduce the computational effort. Substructuring methods can be classified in "mathemati‐ cal" and "physical" methods. The "mathematical" substructuring methods include the Au‐ tomatic 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 theo‐ ry to obtain an abstract (mathematical) substructuring using matrix partitioning of the entire finite-element model. The computational savings from the "mathematical" and "physical" methods can be comparable depending on the problem at hand.

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

$$\mathbf{F}\{\mathbf{K} - \phi^2 \mathbf{M}\}\mathbf{d} = \mathbf{F} \tag{1}$$

of the mode shapes and can be therefore, solved separately reducing further the overall

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

Note that for a damped structure with a damping matrix*C*, Equation (1) becomes

adding the modal damping term . For proportional (structural or Rayleigh) damp‐ ing, *C* is a linear combination of *M* and *K* ; i.e. *C* =*α M* + *β K* where *α*and *β*are con‐ stants. 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

To model the dynamics of a complex structure, a finite-element analysis of the entire struc‐ ture 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. Com‐ ponent mode synthesis (CMS) was developed as a practical and efficient approach to model‐ ing 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 sig‐ nificant computational savings, this component-based approach also facilitates distributed design. Components may be modified or redesigned individually without re-doing the en‐

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 "con‐ straint 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 in‐ cluding 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 inter‐ face DOF. If the finite element mesh is sufficiently fine, the constraint-mode DOFs will dom‐

inate 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-Bamp‐ ton approach. For that, a secondary eigenvalue analysis is performed using the constraintmode partitions of the CMS mass and stiffness matrices. The CC modes are the resultant eigenvectors. The basic formulation is described in Sections 2.2.1 and 2.2.2. The interface modes represent more "natural" physical motion at the interface. Because they capture the

**2.2. Substructuring with Component Mode Synthesis (CMS)**

*M d* = *F* and Equation (3) is modified as by

http://dx.doi.org/10.5772/51402

137

computational effort.

*K* + *jωC* −*ω* <sup>2</sup>

straightforward.

tire analysis.

where the displacement **d** is calculated at the forcing frequency*ω*. If the response is required at multiple frequencies, the repeated direct solution of Equation (1) is computationally very expensive and therefore, impractical for large scale finite-element models.

A reduced order model (ROM) is a subspace projection method. Instead of solving the origi‐ nal 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 calcu‐ late the response more efficiently. A set of eigen-frequencies *ω<sup>i</sup>* and corresponding eigenvec‐ tors (mode shapes) *φ<sup>i</sup>* are first obtained. Then, the displacement **d** is approximated in the reduced space formed by the first *n* dominant modes as

$$\mathbf{d} = \ddot{\Phi}\mathbf{U} \tag{2}$$

where is the modal basis and **U** is the vector of principal coordi‐ nates 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 free‐ dom as

$$\left| \hat{\boldsymbol{\Phi}}^{\top} \mathbf{K} \hat{\boldsymbol{\Phi}} - \boldsymbol{\alpha}^{2} \hat{\boldsymbol{\Phi}}^{\top} \mathbf{M} \hat{\boldsymbol{\Phi}} \right| \mathbf{U} = \hat{\boldsymbol{\Phi}}^{\top} \mathbf{F} \tag{3}$$

The response **d** can be recovered by solving Equation (3) for the modal response **U** and pro‐ jecting it back to the physical coordinates using Equation (2). If *ω*max is the maximum excita‐ tion frequency, it is common practice to retain the mode shapes in the frequency range of 0÷2*ω*max. The system modes in the high frequency range can be safely truncated with mini‐ mal 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

*K* + *jωC* −*ω* <sup>2</sup> *M d* = *F* and Equation (3) is modified as by

adding the modal damping term . For proportional (structural or Rayleigh) damp‐ ing, *C* is a linear combination of *M* and *K* ; i.e. *C* =*α M* + *β K* where *α*and *β*are con‐ stants. 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)**

ry to obtain an abstract (mathematical) substructuring using matrix partitioning of the entire finite-element model. The computational savings from the "mathematical" and "physical"

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

where the displacement **d** is calculated at the forcing frequency*ω*. If the response is required at multiple frequencies, the repeated direct solution of Equation (1) is computationally very

A reduced order model (ROM) is a subspace projection method. Instead of solving the origi‐ nal 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 calcu‐

where is the modal basis and **U** is the vector of principal coordi‐ nates 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 free‐

The response **d** can be recovered by solving Equation (3) for the modal response **U** and pro‐ jecting it back to the physical coordinates using Equation (2). If *ω*max is the maximum excita‐ tion frequency, it is common practice to retain the mode shapes in the frequency range of 0÷2*ω*max. The system modes in the high frequency range can be safely truncated with mini‐

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

are first obtained. Then, the displacement **d** is approximated in the

(1)

and corresponding eigenvec‐

(2)

(3)

<sup>2</sup> [ ] **K Md F** - = w

expensive and therefore, impractical for large scale finite-element models.

late the response more efficiently. A set of eigen-frequencies *ω<sup>i</sup>*

reduced space formed by the first *n* dominant modes as

methods can be comparable depending on the problem at hand.

**2.1. Modal Reduction**

136 Advances in Vibration Engineering and Structural Dynamics

tors (mode shapes) *φ<sup>i</sup>*

dom as

mal loss of accuracy.

To model the dynamics of a complex structure, a finite-element analysis of the entire struc‐ ture 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. Com‐ ponent mode synthesis (CMS) was developed as a practical and efficient approach to model‐ ing 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 sig‐ nificant computational savings, this component-based approach also facilitates distributed design. Components may be modified or redesigned individually without re-doing the en‐ tire 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 "con‐ straint 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 in‐ cluding 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 inter‐ face DOF. If the finite element mesh is sufficiently fine, the constraint-mode DOFs will dom‐ inate 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-Bamp‐ ton approach. For that, a secondary eigenvalue analysis is performed using the constraintmode partitions of the CMS mass and stiffness matrices. The CC modes are the resultant eigenvectors. The basic formulation is described in Sections 2.2.1 and 2.2.2. The interface modes represent more "natural" physical motion at the interface. Because they capture the characteristic motion of the interface, they may be truncated as if they were traditional modes of vibration, leading to a highly compact CC-mode-based reduced order model (CCROM). In addition, the CC modes provide a significant insight into the physical mechanisms of vibration transmission between the component structures. This information could be used, for example, to determine the design parameters that have a critical impact on power flow. Figure 1 compares a conventional constraint mode used in Craig-Bampton analysis with an interface mode, for a simple cantilever plate which is subdivided in two substructures.

þ ý ü

*<sup>N</sup>* = *φi*<sup>1</sup> *φi*<sup>2</sup> ⋯ *φin* which are calculated by

*ΩΩ* {*φin*} *for <sup>n</sup>* =1, 2, ... (5)

*<sup>C</sup>*are calculated by enforcing a set of static unit constraints to


**u u Φ Φ** (7)

&& (8)

<sup>W</sup> **u** can be thus represented by the constraint-mode

*<sup>N</sup>* , and constraint (*C*) modes*Φ<sup>i</sup>*

http://dx.doi.org/10.5772/51402

*C*,

139

(9)

W G

*i i*

**f f**

&& (4)

î í ì = þ ý ü

î í ì ú û

W G

*i i*

**u u**

<sup>ù</sup> <sup>ê</sup>

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

*ii ii*

The fixed-interface Craig-Bampton CMS method utilizes two sets of modes to represent the

where *i* denotes the *i th* substructure. The size reduction of the Craig-Bampton method comes

**kk kk**

WG WW GG GW

ë

*ΩΩΦ<sup>i</sup> C*

<sup>1</sup> *<sup>C</sup> i ii*

*i i i C N N i i i i*

*CC CN C CC CN C C i i i ii i i NC NN N NC NN N N i i i ii i i*

é ùì ü é ùì ü ì ü ê úê ú í ý í ýíý + = ë ûî þ ë ûî þ î þ **mm u kk u f mm u kk u f**

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

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

ìü ì ü é ù <sup>º</sup> íý í ý <sup>=</sup> ê ú îþ î þ ë û **u I 0 u u**

G G

*C*

î í ì ú û

W G

*i i*

fixing all interface DOFs and solving the following eigenvalue problem

*<sup>C</sup>* <sup>=</sup>*<sup>Λ</sup> <sup>N</sup> <sup>m</sup><sup>i</sup>*

<sup>G</sup> **u** and *<sup>i</sup>*

&&

W

*ΩΩ* {*φin*} <sup>=</sup>*λ<sup>n</sup> mi*

**u u**

&&

<sup>ù</sup> <sup>ê</sup>

*ii ii*

**mm mm**

from the truncation of the normal modes *Φ<sup>i</sup>*

*ki ΩΩ*−1*Φ<sup>i</sup>*

*ki*

The static constraint modes*Φ<sup>i</sup>*

The original physical DOFs *<sup>i</sup>*

DOFs *<sup>C</sup>* **u***<sup>i</sup>* and the normal-mode DOFs *<sup>N</sup>* **u***<sup>i</sup>* as

the interface DOFs as

expressed as

Equation (8) are

WG WW GG GW

substructure motion; substructure normal (*N*) modes**Φ***<sup>i</sup>*

ë

é

é + þ ý ü

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 intro‐ ducing 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 espe‐ cially 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 assump‐ tion. The method is commonly used in CMS algorithms. The finite-element model of the en‐ tire 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 *i th* substructure are then expressed as

Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402 139

$$
\begin{bmatrix}
\mathbf{m}\_i^{\Pi^\Gamma} & \mathbf{m}\_i^{\Pi\Omega} \\
\mathbf{m}\_i^{\Omega\Gamma} & \mathbf{m}\_i^{\Omega\Omega}
\end{bmatrix}
\begin{Bmatrix}
\ddot{\mathbf{u}}\_i^{\Gamma} \\
\ddot{\mathbf{u}}\_i^{\Omega}
\end{Bmatrix} + \begin{bmatrix}
\mathbf{k}\_i^{\Pi^\Gamma} & \mathbf{k}\_i^{\Gamma\Omega} \\
\mathbf{k}\_i^{\Omega\Gamma} & \mathbf{k}\_i^{\Omega\Omega}
\end{bmatrix}
\begin{Bmatrix}
\mathbf{u}\_i^{\Gamma} \\
\mathbf{u}\_i^{\Omega}
\end{Bmatrix} = \begin{Bmatrix}
\mathbf{f}\_i^{\Gamma} \\
\mathbf{f}\_i^{\Omega}
\end{Bmatrix} \tag{4}
$$

The fixed-interface Craig-Bampton CMS method utilizes two sets of modes to represent the substructure motion; substructure normal (*N*) modes**Φ***<sup>i</sup> <sup>N</sup>* , and constraint (*C*) modes*Φ<sup>i</sup> C*, where *i* denotes the *i th* substructure. The size reduction of the Craig-Bampton method comes from the truncation of the normal modes *Φ<sup>i</sup> <sup>N</sup>* = *φi*<sup>1</sup> *φi*<sup>2</sup> ⋯ *φin* which are calculated by fixing all interface DOFs and solving the following eigenvalue problem

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 intro‐ ducing 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 espe‐

This section provides the basics of Craig-Bampton method using the fixed-interface assump‐ tion. The method is commonly used in CMS algorithms. The finite-element model of the en‐ tire 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 *i th*

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

138 Advances in Vibration Engineering and Structural Dynamics

substructure are then expressed as

$$\begin{aligned} \mathbf{k}\_i^{\Omega\Omega}{}^{-1}\boldsymbol{\Phi}\_i^{\mathbb{C}} &= \boldsymbol{\Lambda}^N \boldsymbol{m}\_i^{\Omega\Omega} \boldsymbol{\Phi}\_i^{\mathbb{C}}\\ \mathbf{k}\_i^{\Omega\Omega} \mathbf{l}[\boldsymbol{\varphi}\_{in}] &= \boldsymbol{\lambda}\_n \mathbf{l}[\boldsymbol{m}\_i^{\Omega\Omega}] \{\boldsymbol{\varphi}\_{in}\} \quad \text{for} \quad n = 1, 2, \dots \end{aligned} \tag{5}$$

The static constraint modes*Φ<sup>i</sup> <sup>C</sup>*are calculated by enforcing a set of static unit constraints to the interface DOFs as

$$
\begin{bmatrix}
\mathbf{d}\mathbf{p}^{\rm C} \\
\end{bmatrix} = -\begin{bmatrix}
\mathbf{k}\_{i}^{\alpha\alpha} \\
\end{bmatrix}^{-1} \begin{bmatrix}
\mathbf{k}\_{i}^{\alpha\mathbf{r}^{\rm C}} \\
\end{bmatrix} \tag{6}
$$

The original physical DOFs *<sup>i</sup>* <sup>G</sup> **u** and *<sup>i</sup>* <sup>W</sup> **u** can be thus represented by the constraint-mode DOFs *<sup>C</sup>* **u***<sup>i</sup>* and the normal-mode DOFs *<sup>N</sup>* **u***<sup>i</sup>* as

$$\begin{Bmatrix} \mathbf{u}\_{i}^{\Gamma} \\ \mathbf{u}\_{i}^{\Omega} \end{Bmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0}\mathbf{0}\_{i}^{C} & \mathbf{0}\mathbf{0}\_{i}^{N} \end{bmatrix} \begin{Bmatrix} \mathbf{u}\_{i}^{C} \equiv \mathbf{u}\_{i}^{\Gamma} \\ \mathbf{u}\_{i}^{N} \end{Bmatrix} \tag{7}$$

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

$$
\begin{bmatrix}
\mathbf{m}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{C}} & \mathbf{m}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}} \\
\mathbf{m}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{C}} & \mathbf{m}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}}
\end{bmatrix}
\begin{Bmatrix}
\ddot{\mathbf{u}}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{C}} \\
\ddot{\mathbf{u}}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}}
\end{Bmatrix} + \begin{bmatrix}
\mathbf{k}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{C}} & \mathbf{k}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}} \\
\mathbf{k}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{C}} & \mathbf{k}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}}
\end{bmatrix}
\begin{Bmatrix}
\mathbf{u}\_{\boldsymbol{i}}^{\boldsymbol{\prime}} \\
\mathbf{u}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}}
\end{Bmatrix} = \begin{Bmatrix}
\mathbf{f}\_{\boldsymbol{i}}^{\boldsymbol{C}} \\
\mathbf{f}\_{\boldsymbol{i}}^{\boldsymbol{\prime}\boldsymbol{N}}
\end{Bmatrix} \tag{8}
$$

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

$$\mathbf{m}\_{\ell}^{\rm CC} = \mathbf{m}\_{\ell}^{\rm IT} + \mathbf{m}\_{\ell}^{\rm T\Omega} \boldsymbol{\Phi}\_{\ell}^{\rm C} + \left(\boldsymbol{\Phi}\_{\ell}^{\rm C}\right)^{\rm T} \mathbf{m}\_{\ell}^{\rm DT} + \left(\boldsymbol{\Phi}\_{\ell}^{\rm C}\right)^{\rm T} \mathbf{m}\_{\ell}^{\rm \rm \Omega \Omega} \boldsymbol{\Phi}\_{\ell}^{\rm C} \tag{9}$$

$$\mathbf{m}\_{i}^{\mathsf{CIV}} = \left(\mathbf{m}\_{i}^{\mathsf{NC}}\right)^{T} = \mathbf{m}\_{i}^{\mathsf{T}\Omega} \boldsymbol{\Phi}\_{i}^{\mathcal{N}} + \left(\boldsymbol{\Phi}\_{i}^{\mathcal{C}}\right)^{T} \mathbf{m}\_{i}^{\Omega\Omega} \boldsymbol{\Phi}\_{i}^{\mathcal{N}} \tag{10}$$

$$\mathbf{m}\_{i}^{\rm{NN}} = \left(\boldsymbol{\Phi}\_{i}^{\rm{N}}\right)^{T} \mathbf{m}\_{i}^{\rm{\alpha}\rm{\alpha}} \boldsymbol{\Phi}\_{i}^{\rm{N}} \tag{11}$$

1 2 *ss*

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

is the number of substructures in the global structure. The corresponding synthe‐

L L L

L

*<sup>n</sup>* <sup>=</sup> é ù ë û **d ddd d** <sup>L</sup> (17)

*SS*

(18)

141

http://dx.doi.org/10.5772/51402

*n*

*<sup>T</sup> CMS CT N T N T N T*

1 2

*C CN CN CN*

é ù ê ú

**mmm m mm0 0**

*SS SS*

L L L

L

where the component modal matrices *<sup>N</sup>* **m***i* and *<sup>N</sup>* **k***i* 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 *<sup>C</sup>* **m** and *<sup>C</sup>* **k** , 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 con‐ straint-mode partitions of the CMS matrices may be relatively large. For this reason, we fur‐ ther reduce the CMS matrices by performing a modal analysis on the constraint-mode DOFs

The eigenvectors *ψn* are transformed into the finite-element DOFs for the *i th* component

*CN T N n n*

M M MOM

ë û

**m 00 m**

*SS*

*<sup>k</sup>***¯***<sup>C</sup>ψ<sup>n</sup>* <sup>=</sup>*Λn<sup>m</sup>***¯***<sup>C</sup>ψ<sup>n</sup> for <sup>n</sup>* =1, 2, 3, ... (19)

*CC C C* **Φ ΦβΨ** *i ii* = (20)

*Ψ* = *ψ*<sup>1</sup> *ψ*<sup>2</sup> ⋯ *ψ<sup>n</sup> CC* (21)

*N n*

1 1

*CN T N CMS CN T N*

1

**K 0 0k 0**

*N*

*C*

=

=

structure using the following transformation

*CMS N*

2 2

**M m 0m 0**

2

M M MO M

ë û

**000 k**

é ù ê ú

**k 00 0 0k 0 0**

sized CMS mass and stiffness matrices are as follows

where *n ss*

as follows

where

$$\mathbf{k}\_i^{\rm CC} = \mathbf{k}\_i^{\rm IT} + \mathbf{k}\_i^{\rm IT\Omega} \boldsymbol{\Phi}\_i^C = \mathbf{k}\_i^{\rm IT} - \mathbf{k}\_i^{\rm IT\Omega} \left(\mathbf{k}\_i^{\rm \Omega \Omega}\right)^{-1} \mathbf{k}\_i^{\rm \Omega \Gamma} \tag{12}$$

$$\mathbf{k}\_i^{\text{CN}} = \left(\mathbf{k}\_i^{\text{NC}}\right)^T = \mathbf{0} \tag{13}$$

$$\mathbf{k}\_i^{\rm{NN}} = \left(\boldsymbol{\Phi}\_i^N\right)^T \mathbf{k}\_i^{\rm{\alpha\Omega}} \boldsymbol{\Phi}\_i^N = \mathbf{D}\_i \equiv \boldsymbol{\alpha} \mathbf{k}\_i \tag{14}$$

$$\mathbf{f}\_i^{\mathbf{C}} = \mathbf{f}\_i^{\Gamma} + \left(\boldsymbol{\Phi}\_i^{\mathbf{C}}\right)^T \mathbf{f}\_i^{\mathbf{C}} \tag{15}$$

$$\mathbf{f}\_{i}^{N} = \left(\boldsymbol{\Phi}\_{i}^{N}\right)^{T} \mathbf{f}\_{i}^{\boldsymbol{\Omega}} \tag{16}$$

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


#### *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 *CMS* **d** of the synthesized system to be partitioned as

Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402 141

$$\mathbf{d}^{\rm CMS} = \begin{bmatrix} \mathbf{d}^{\rm CT} & \mathbf{d}\_1^{\rm NT} & \mathbf{d}\_2^{\rm NT} & \cdots & \mathbf{d}\_{n^{\rm ref}}^{\rm N} \end{bmatrix}^T \tag{17}$$

where *n ss* is the number of substructures in the global structure. The corresponding synthe‐ sized CMS mass and stiffness matrices are as follows

$$\mathbf{M}^{CKS} = \begin{bmatrix} \mathbf{m}^C & \mathbf{m}\_1^{CN} & \mathbf{m}\_2^{CN} & \cdots & \mathbf{m}\_{\boldsymbol{\kappa}^{\mathcal{K}}}^{CN} \\ \mathbf{m}\_1^{CN} & \mathbf{m}\_1^{\mathcal{N}} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{m}\_{2^{\mathcal{N}}}^{CN} & \mathbf{0} & \mathbf{m}\_2^{\mathcal{N}} & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{m}\_{\boldsymbol{\kappa}^{\mathcal{K}}}^{CN} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{m}\_{\boldsymbol{\kappa}^{\mathcal{K}}}^{N} \end{bmatrix} \tag{18}$$
 
$$\mathbf{K}^{Class} = \begin{bmatrix} \mathbf{\bar{K}}^C & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{k}\_1^N & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{k}\_2^N & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{k}\_{\boldsymbol{\kappa}^{\mathcal{K}}}^{N} \end{bmatrix} \tag{19}$$

where the component modal matrices *<sup>N</sup>* **m***i* and *<sup>N</sup>* **k***i* 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 *<sup>C</sup>* **m** and *<sup>C</sup>* **k** , 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 con‐ straint-mode partitions of the CMS matrices may be relatively large. For this reason, we fur‐ ther reduce the CMS matrices by performing a modal analysis on the constraint-mode DOFs as follows

$$\overline{\mathbf{k}}^{C}\psi\_{n} = \Lambda\_{n}\overline{\mathbf{m}}^{C}\psi\_{n} \quad \text{for} \quad n = 1, 2, 3, \dots \tag{19}$$

The eigenvectors *ψn* are transformed into the finite-element DOFs for the *i th* component structure using the following transformation

$$\boldsymbol{\Phi}\_{l}^{\text{CC}} = \boldsymbol{\Phi}\_{l}^{\text{C}} \boldsymbol{\mathfrak{f}}\_{l}^{\text{C}} \boldsymbol{\Psi} \tag{20}$$

where

(10)

(11)

(12)

(14)

(15)

(16)

( )

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

Note that constraint mode matrix *<sup>C</sup>* **Φ***i* is usually a full matrix. Therefore Equation (9) can be

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 *CMS* **d** of the synthesized system to be

modes. The computational cost of the Craig-Bamtpon method is mostly related to

system to obtain the eigenvalues and eigenvectors.

140 Advances in Vibration Engineering and Structural Dynamics

**1.** solving for the normal modes,

partitioned as

**2.** solving for the constraint modes, and

*2.2.2. Craig-Bampton CMS with Interface Modes*

**3.** the transformation calculation in Equation (9).

computationally expensive due to the triple-product ( )

*<sup>T</sup> CN NC* **kk 0** *i i* = = (13)

*<sup>T</sup> C C i ii*

WW **Φ mΦ** involving constraint

$$\mathbf{\dot{\Psi}} = \begin{bmatrix} \psi\_1 & \psi\_2 & \cdots & \psi\_n \end{bmatrix} \tag{21}$$

is a selected set of *n CC* interface eigenvectors which are few compared to the number of the constraint-mode DOFs. The matrix *<sup>C</sup>* **β***i* maps the global (system) interface DOFs *<sup>C</sup>* **d** back to the local (subsystem *i*) DOF *<sup>C</sup>* **d***<sup>i</sup>* . The vectors in *CC* **Φ***i* 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

$$\mathbf{d}^{ROM} = \begin{bmatrix} \mathbf{d}^{CCT} & \mathbf{d}\_1^{NT} & \mathbf{d}\_2^{NT} & \cdots & \mathbf{d}\_{\mu^u}^{N/T} \end{bmatrix}^T \tag{22}$$

where

and

*2.2.3. Filtration of Constraint Modes*

**Figure 2.** Illustration of a "filtered" constraint mode.

constraint modes, the following criterion is used

*φpq*

*<sup>C</sup>* =0 *if* <sup>|</sup>*φpq*

the constraint mode density reduces from 100% to 16% if*ε* =0.03.

turbed interface DOF.

where *φpq*

**m¯***i*

*CN* <sup>=</sup>**Ψ***<sup>T</sup>* **<sup>β</sup>***<sup>i</sup>*

*C <sup>T</sup>* **m***i*

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

Figure 2 shows a typical constraint mode for a plate substructure. The non-zero displace‐ ment field (indicated by red color) is usually limited to a very small region close to the per‐

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. Conse‐ quently, the computational cost in Equation (9) will be considerably reduced. To filter the

> *<sup>C</sup>* | <*ε* \* max*p*

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,

*<sup>C</sup>* is the *p th* element of the *q th* constraint mode *<sup>C</sup>* **φ** . If the ratio of an element of the


*<sup>C</sup>* | (29)

(27)

143

*CN* (28)

http://dx.doi.org/10.5772/51402

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[-\phi^2 \mathbf{M}^{ROM} + \mathbf{K}^{ROM}\right] \mathbf{d}^{ROM} = \mathbf{f}^{ROM} \tag{23}$$

The mass matrix **M**ROM, the stiffness matrix **K**ROM, and the applied force vector **f**ROM, are ex‐ plicitly written as

$$\begin{aligned} \mathbf{M}^{ROM} &= \begin{bmatrix} \bar{\mathbf{m}}^{CC} & \bar{\mathbf{m}}\_1^{CN} & \bar{\mathbf{m}}\_2^{CN} & \cdots & \bar{\mathbf{m}}\_{x^{\mathcal{R}}}^{CN} \\ \bar{\mathbf{m}}\_1^{CN} & \mathbf{m}\_1^{N} & \mathbf{0} & \cdots & \mathbf{0} \\ \bar{\mathbf{m}}\_2^{CN} & \mathbf{0} & \mathbf{m}\_2^{N} & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \bar{\mathbf{m}}\_{x^{\mathcal{R}}}^{CN} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{m}\_{x^{\mathcal{R}}}^{N} \end{bmatrix} \\\\ \mathbf{K}^{ROM} &= \begin{bmatrix} \bar{\mathbf{K}}^{CC} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{k}\_1^{N} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{k}\_2^{N} & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{k}\_{x^{\mathcal{R}}}^{N} \end{bmatrix} \end{aligned} \tag{25}$$

$$\mathbf{f}^{ROM} = \begin{bmatrix} \mathbf{f}^{CCT} & \mathbf{f}\_1^{NT} & \mathbf{f}\_2^{NT} & \cdots & \mathbf{f}\_{n^{\rm IN}}^{N^{\rm IN}} \end{bmatrix}^T \tag{26}$$

Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402 143

where

is a selected set of *n CC* interface eigenvectors which are few compared to the number of the constraint-mode DOFs. The matrix *<sup>C</sup>* **β***i* maps the global (system) interface DOFs *<sup>C</sup>* **d** back to

the local (subsystem *i*) DOF *<sup>C</sup>* **d***<sup>i</sup>* . The vectors in *CC* **Φ***i* 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

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

> *<sup>T</sup> ROM CCT N T N T N T <sup>n</sup>* = é ù

<sup>2</sup> *ROM ROM ROM ROM* é ù -+ =

motion of the reduced order CMS model (ROM) are expressed by

w

=

1 2 *ss*

where superscript *CC* indicates the partition associated with the *CC* modes. The equations of

The mass matrix **M**ROM, the stiffness matrix **K**ROM, and the applied force vector **f**ROM, are ex‐

1 2

**m mm m mm0 0**

*CC CN CN CN*

é ù ê ú

*SS SS*

*CN T N n n*

2

ë û

1 2 *SS*

M M MO M

**0 00 k**

*<sup>T</sup> ROM CCT N T N T N T <sup>n</sup>* = é ù

M M MOM

ë û

é ù ê ú

**k 00 0 0k 0 0**

**m 00 m**

1 1

*CC*

*ROM N*

=

*CN T N ROM CN T N*

2 2

1

**K 0 0k 0**

*N*

**M m 0m 0**

ë û **d d dd d** L (22)

ë û **M Kd f** (23)

*SS*

(24)

(25)

*n*

L L L

L

L L L

L

*SS*

ë û **f f ff f** L (26)

*N n*

the number of interface DOFs.

142 Advances in Vibration Engineering and Structural Dynamics

**d**ROM is partitioned as

plicitly written as

$$
\bar{\mathbf{m}}^{\mathbb{C}\mathbb{C}} = \Psi^{\mathbb{T}} \bar{\mathbf{m}}^{\mathbb{C}} \Psi, \; \bar{\mathbf{k}}^{\mathbb{C}\mathbb{C}} = \Psi^{\mathbb{T}} \bar{\mathbf{k}}^{\mathbb{C}} \Psi \tag{27}
$$

and

$$\overline{\mathbf{m}}\_{i}^{\text{CN}} = \mathbf{W}^{T} \mathbf{g}\_{i}^{\text{C}} \overset{\text{T}}{\mathbf{m}}\_{i}^{\text{CN}} \tag{28}$$

#### *2.2.3. Filtration of Constraint Modes*

Figure 2 shows a typical constraint mode for a plate substructure. The non-zero displace‐ ment field (indicated by red color) is usually limited to a very small region close to the per‐ turbed 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. Conse‐ quently, the computational cost in Equation (9) will be considerably reduced. To filter the constraint modes, the following criterion is used

$$\left| \left| \boldsymbol{\varrho} \right|\_{pq}^{\mathbb{C}} = 0 \quad \text{if} \quad \left| \left| \boldsymbol{\varrho} \right|\_{pq}^{\mathbb{C}} \right| < \varepsilon \text{ } \text{max} \left| \left| \boldsymbol{\varrho} \right|\_{pq}^{\mathbb{C}} \right| \text{ } \tag{29}$$

where *φpq <sup>C</sup>* is the *p th* element of the *q th* constraint mode *<sup>C</sup>* **φ** . 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 re‐ duced-order model of small size can be developed using an FRF representation of each sub‐ structure. 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 funda‐ mentals of FRF substructuring for both FRF-FE and FRF-FRF coupling.

The procedure to assemble the **H** matrix of Sub1 with the **Z** matrix of Sub2 and calculate

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

*A AA AC A C CA CC C* é ù é ùé ù <sup>=</sup> ê ú ê úê ú ë û ë ûë û **X HH F**

*B BB BC B C CB CC C* é ù é ùé ù <sup>=</sup> ê ú ê úê ú ë û ë ûë û **F ZZ X**

where subscript A indicates the interior (excitation plus response) DOFs of Sub1, and sub‐ script C indicates the connection/common/coupling DOFs between Sub1 and Sub2. The

where subscript *B* indicates the interior DOFs of Sub2, and subscript *C* indicates the connec‐ tion/common/coupling DOFs between Sub2 and Sub1. Because of displacement compatibili‐ ty at the interface, *XC*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

To couple Sub1 and Sub2, compatibility of forces at the interface is applied as 1 2 **FFF** *C CC* = +

the second row of Equation (30) and the second row of Equation (31) provides the force vec‐

1

**X HH F** (30)

**F ZZ X** (31)

*CC CA*


http://dx.doi.org/10.5772/51402

145

(32)

(33)

(34)

(solve for) the system matrix **H** is described below.

The equations of motion for Sub1 are expressed as

equations of motion for Sub2 are expressed as

interface forces **F**C at Sub1 and Sub2.

tor <sup>2</sup> **F ZX ZX** *C CC C CB B* = + . We thus have

From Equation (31),

2

where the force vector with <sup>1</sup>

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

## *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 ma‐ trices which are denoted by **H** (see Equation 30 for notation). The elements of **H** are frequen‐ cy 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**<sup>A</sup> of DOF A due to a unit force **F**C on DOF C. Sub2 is a finite element (FE) type substructure. Its dynamic be‐ havior 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 fre‐ quency domain while in the latter a modal model is first developed and then used to calcu‐ late 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 excita‐ tion and response DOFs. However in order to assemble two substructures, FRFs are also cal‐ culated between interface DOFs and excitation/response DOFs. The Sub1 FRF matrix **H** in Equation (30) is thus partitioned into interior (A) DOFs and interface/coupling (C) DOFs. The interior DOFs include all excitation and all response DOFs (Figure 3).

The second substructure Sub2 is expressed in FE format. The system FE matrices **K**, **M**, and **B** form the frequency dependent dynamic matrix **Z**=**K** + *iω***B**−*ω* <sup>2</sup> **M** which is then parti‐ tioned 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

**2.3. Frequency Response Function (FRF) Substructuring and Assembling**

mentals of FRF substructuring for both FRF-FE and FRF-FRF coupling.

havior is described using the stiffness **K**, mass **M** and damping **B** matrices.

The interior DOFs include all excitation and all response DOFs (Figure 3).

**B** form the frequency dependent dynamic matrix **Z**=**K** + *iω***B**−*ω* <sup>2</sup>

*2.3.1. Algorithm for FRF/FE Coupling*

144 Advances in Vibration Engineering and Structural Dynamics

**Figure 3.** Two-substructure example and notation.

If the number of interface nodes (or DOFs) between connected substructures is small, a re‐ duced-order model of small size can be developed using an FRF representation of each sub‐ structure. 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 funda‐

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 ma‐ trices which are denoted by **H** (see Equation 30 for notation). The elements of **H** are frequen‐ cy 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**<sup>A</sup> of DOF A due to a unit force **F**C on DOF C. Sub2 is a finite element (FE) type substructure. Its dynamic be‐

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 fre‐ quency domain while in the latter a modal model is first developed and then used to calcu‐ late 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 excita‐ tion and response DOFs. However in order to assemble two substructures, FRFs are also cal‐ culated 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 second substructure Sub2 is expressed in FE format. The system FE matrices **K**, **M**, and

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

**M** which is then parti‐

$$
\begin{bmatrix} \mathbf{X}\_{\mathcal{A}} \\ \mathbf{X}\_{\mathcal{C}} \end{bmatrix} = \begin{bmatrix} \mathbf{H}\_{\mathcal{A}d} & \mathbf{H}\_{\mathcal{A}\mathcal{C}} \\ \mathbf{H}\_{\mathcal{C}d} & \mathbf{H}\_{\mathcal{C}\mathcal{C}} \end{bmatrix} \begin{bmatrix} \mathbf{F}\_{\mathcal{A}} \\ \mathbf{F}\_{\mathcal{C}}^{\mathrm{I}} \end{bmatrix} \tag{30}
$$

where subscript A indicates the interior (excitation plus response) DOFs of Sub1, and sub‐ script C indicates the connection/common/coupling DOFs between Sub1 and Sub2. The equations of motion for Sub2 are expressed as

$$
\begin{bmatrix}
\mathbf{F}\_{\mathcal{B}} \\
\mathbf{F}\_{\mathcal{C}}^{2}
\end{bmatrix} = \begin{bmatrix}
\mathbf{Z}\_{\mathcal{B}\mathcal{B}} & \mathbf{Z}\_{\mathcal{B}\mathcal{C}} \\
\mathbf{Z}\_{\mathcal{C}\mathcal{B}} & \mathbf{Z}\_{\mathcal{C}\mathcal{C}}
\end{bmatrix} \begin{bmatrix}
\mathbf{X}\_{\mathcal{B}} \\
\mathbf{X}\_{\mathcal{C}}
\end{bmatrix} \tag{31}
$$

where subscript *B* indicates the interior DOFs of Sub2, and subscript *C* indicates the connec‐ tion/common/coupling DOFs between Sub2 and Sub1. Because of displacement compatibili‐ ty at the interface, *XC*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 1 2 **FFF** *C CC* = + where the force vector with <sup>1</sup> *CC CA* - **ΦHH** = is obtained from the second row of Equation (30) and the second row of Equation (31) provides the force vec‐ tor <sup>2</sup> **F ZX ZX** *C CC C CB B* = + . We thus have

$$\mathbf{F}\_{\text{C}} - \left(\mathbf{H}\_{\text{CC}}^{-1} + \mathbf{Z}\_{\text{CC}}\right)\mathbf{X}\_{\text{C}} - \Phi \mathbf{F}\_{\text{A}} + \mathbf{Z}\_{\text{CB}} \mathbf{X}\_{\text{B}} \tag{32}$$

From Equation (31),

$$\mathbf{X}\_B = \mathbf{Z}\_{BB}^{-1} \mathbf{F}\_B - \mathbf{Z}\_{BB}^{-1} \mathbf{Z}\_{BC} \mathbf{X}\_C = \mathbf{Z}\_{BB}^{-1} \mathbf{F}\_B + \boldsymbol{\Theta}^T \mathbf{X}\_C \tag{33}$$

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

$$\begin{aligned} \mathbf{F}\_{\mathcal{C}} &= \left( \mathbf{H}\_{\mathcal{C}\mathcal{C}}^{-1} + \mathbf{Z}\_{\mathcal{C}\mathcal{C}} + \mathbf{Z}\_{\mathcal{C}\mathcal{B}} \boldsymbol{\Theta}^{T} \right) \mathbf{X}\_{\mathcal{C}} - \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{A}} - \boldsymbol{\Theta} \mathbf{F}\_{\mathcal{B}} \\ &- \mathbf{R} \mathbf{X}\_{\mathcal{C}} - \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{A}} - \boldsymbol{\Theta} \mathbf{F}\_{\mathcal{B}} \end{aligned} \tag{34}$$

where . From Equation (34), **X**<sup>C</sup> can be expressed in terms of **F**<sup>A</sup> , **F**B and **F**<sup>C</sup> as

$$\mathbf{X}\_{\circ} = \mathbf{R}^{-1} \boldsymbol{\Phi} \mathbf{F}\_{\circ} + \mathbf{R}^{-1} \boldsymbol{\Theta} \mathbf{F}\_{\circ} + \mathbf{R}^{-1} \mathbf{F}\_{\circ} \tag{35}$$

(40)

147

http://dx.doi.org/10.5772/51402

where **S** = [**Φ Θ I**] and .

Figure 4 shows the coupling of two FRF type substructures. The equations of motion for

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

*A AA AC A C CA CC C* é ù é ùé ù <sup>=</sup> ê ú ê úê ú ë û ë ûë û **X HH F**

*B BB BC B C CB CC C* é ù é ùé ù <sup>=</sup> ê ú ê úê ú ë û ë ûë û **X HH F**

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

1 1

2 2

**X HH F** (41)

**X HH F** (42)

*2.4.2. Algorithm for FRF/FRF Coupling*

Sub1 and Sub2 and are expressed as

**Figure 4.** Two FRF type substructures example and notation

be re-arranged in matrix form as

and

Substitution of Equation (35) in Equation (33) gives **X**B in terms of **F**<sup>A</sup> , **F**B and **F**C as

$$\begin{split} \mathbf{X}\_{B} &= \mathbf{Z}\_{BB}^{-1} \mathbf{F}\_{B} + \boldsymbol{\Theta}^{T} \left( \mathbf{R}^{-1} \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{A}} + \mathbf{R}^{-1} \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{B}} + \mathbf{R}^{-1} \mathbf{F}\_{\mathcal{C}} \right) \\ &= \boldsymbol{\Theta}^{T} \mathbf{R}^{-1} \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{A}} + \left( \mathbf{Z}\_{BB}^{-1} + \boldsymbol{\Theta}^{T} \mathbf{R}^{-1} \boldsymbol{\Theta} \right) \mathbf{F}\_{B} + \boldsymbol{\Theta}^{T} \mathbf{R}^{-1} \mathbf{F}\_{\mathcal{C}} \end{split} \tag{36}$$

Solving Equation (30) for and substituting yields

$$\begin{aligned} \mathbf{X}\_{\mathcal{A}} &= \mathbf{H}\_{\mathcal{A}\mathcal{A}} \mathbf{F}\_{\mathcal{A}} + \mathbf{H}\_{\mathcal{A}\mathcal{C}} \left( \mathbf{H}\_{\mathcal{C}\mathcal{C}}^{-1} \mathbf{X}\_{\mathcal{C}} - \boldsymbol{\Phi} \mathbf{F}\_{\mathcal{A}} \right) \\ &- \left( \mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}} \boldsymbol{\Phi} \right) \mathbf{F}\_{\mathcal{A}} + \boldsymbol{\Phi}^{T} \mathbf{X}\_{\mathcal{C}} \end{aligned} \tag{37}$$

We can now express **X**A in terms of **F**<sup>A</sup> , **F**B and **F**C by substituting Equation (35) in Equation (37), as

$$\begin{aligned} \mathbf{X}\_{A} &= \left(\mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}}\boldsymbol{\Phi}\right) \mathbf{F}\_{A} + \boldsymbol{\Phi}^{T} \mathbf{X}\_{C} \\ &= \left(\mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}}\boldsymbol{\Phi}\right) \mathbf{F}\_{A} + \boldsymbol{\Phi}^{T} \left(\mathbf{R}^{-1}\boldsymbol{\Phi}\mathbf{F}\_{A} + \mathbf{R}^{-1}\boldsymbol{\Phi}\mathbf{F}\_{B} + \mathbf{R}^{-1}\mathbf{F}\_{C}\right) \\ &= \left(\mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}}\boldsymbol{\Phi} + \boldsymbol{\Phi}^{T}\mathbf{R}^{-1}\boldsymbol{\Phi}\right) \mathbf{F}\_{A} + \boldsymbol{\Phi}^{T}\mathbf{R}^{-1}\boldsymbol{\Phi}\mathbf{F}\_{B} + \boldsymbol{\Phi}^{T}\mathbf{R}^{-1}\mathbf{F}\_{C} \end{aligned} \tag{38}$$

Based on Equations (38), (36) and (35), **X**A, **X**B and **X**<sup>C</sup> are expressed in terms of **F**<sup>A</sup> , **F**B and **F**<sup>A</sup> as follows

$$
\begin{bmatrix} \mathbf{X}\_{\mathcal{A}} \\ \mathbf{X}\_{\mathcal{B}} \\ \mathbf{X}\_{\mathcal{C}} \end{bmatrix} = \begin{bmatrix} \left[ \mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}} \boldsymbol{\Phi} \mathbf{I} \\\\ \mathbf{Z}\_{\mathcal{B}\mathcal{B}}^{-1} \\\\ \mathbf{0} \end{bmatrix} + \begin{bmatrix} \boldsymbol{\Phi}^{T} \\\\ \boldsymbol{\Phi}^{T} \\\\ \mathbf{I} \end{bmatrix} \mathbf{R}^{-1} \begin{bmatrix} \boldsymbol{\Phi} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \end{bmatrix} \begin{bmatrix} \mathbf{F}\_{\mathcal{A}} \\\\ \mathbf{F}\_{\mathcal{B}} \\\\ \mathbf{F}\_{\mathcal{C}} \end{bmatrix} \tag{39}
$$

resulting in the following FRF of the assembled system

Vibration and Optimization Analysis of Large-Scale Structures using Reduced-Order Models and Reanalysis Methods http://dx.doi.org/10.5772/51402 147

$$\begin{aligned} \mathbf{H}\_{\mathcal{g}z} &= \begin{bmatrix} \mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}}\boldsymbol{\Phi} \\ & \mathbf{Z}\_{\mathcal{B}\mathcal{B}}^{-1} \\ & \mathbf{0} \end{bmatrix} + \begin{bmatrix} \boldsymbol{\Phi}^{\mathrm{T}} \\ \boldsymbol{\Phi}^{\mathrm{T}} \\ \mathbf{I} \end{bmatrix} \mathbf{R}^{-1} [\boldsymbol{\Phi} \quad \boldsymbol{\Phi} \quad \mathbf{I}] \\ &= \operatorname{diag} \Big( \mathbf{H}\_{\mathcal{A}\mathcal{A}} - \mathbf{H}\_{\mathcal{A}\mathcal{C}} \boldsymbol{\Phi}, \quad \mathbf{Z}\_{\mathcal{B}\mathcal{B}}^{-1}, \quad \mathbf{0} \Big) + \mathbf{S}^{\mathrm{T}} \mathbf{R}^{-1} \mathbf{S} \\ &= \mathbf{H} \mathbf{K} + \mathbf{S}^{\mathrm{T}} \mathbf{R}^{-1} \mathbf{S} \end{aligned} \tag{40}$$

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

$$
\begin{bmatrix}
\mathbf{X}\_A \\
\mathbf{X}\_C \\
\end{bmatrix} = \begin{bmatrix}
\mathbf{H}\_{\mathcal{A}} & \mathbf{H}\_{\mathcal{A}} \\
\mathbf{H}\_{\mathcal{C}\mathcal{A}} & \mathbf{H}\_{\mathcal{C}\mathcal{C}}^\dagger \\
\end{bmatrix} \begin{bmatrix}
\mathbf{F}\_{\mathcal{A}} \\
\mathbf{F}\_{\mathcal{C}}^\dagger \\
\end{bmatrix} \tag{41}
$$

and

(35)

(36)

(37)

(38)

(39)

where . From Equation (34), **X**<sup>C</sup> can be expressed in terms of **F**<sup>A</sup> , **F**B and **F**<sup>C</sup>

Substitution of Equation (35) in Equation (33) gives **X**B in terms of **F**<sup>A</sup> , **F**B and **F**C as

Solving Equation (30) for and substituting yields

We can now express **X**A in terms of **F**<sup>A</sup> , **F**B and **F**C by substituting Equation (35) in Equation

Based on Equations (38), (36) and (35), **X**A, **X**B and **X**<sup>C</sup> are expressed in terms of **F**<sup>A</sup> , **F**B and **F**<sup>A</sup>

resulting in the following FRF of the assembled system

as

146 Advances in Vibration Engineering and Structural Dynamics

(37), as

as follows

$$
\begin{bmatrix} \mathbf{X}\_{B} \\ \mathbf{X}\_{C} \end{bmatrix} = \begin{bmatrix} \mathbf{H}\_{BB} & \mathbf{H}\_{BC} \\ \mathbf{H}\_{CB} & \mathbf{H}\_{CC}^{2} \end{bmatrix} \begin{bmatrix} \mathbf{F}\_{B} \\ \mathbf{F}\_{C}^{2} \end{bmatrix} \tag{42}
$$

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

$$
\begin{bmatrix} \mathbf{X}\_A \\ \mathbf{X}\_B \\ \mathbf{X}\_C \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} \mathbf{H}\_{\mathcal{A}} & \mathbf{H}\_{\mathcal{A}C} \\ \mathbf{H}\_{\mathcal{A}} & \mathbf{H}\_{\mathcal{C}C}^1 \\ & & \mathbf{H}\_{\mathcal{B}B} \end{bmatrix} + \begin{bmatrix} \mathbf{H}\_{\mathcal{A}C} \\ \mathbf{H}\_{\mathcal{C}C}^1 \\ \mathbf{H}\_{\mathcal{B}C} \end{bmatrix} \begin{bmatrix} \mathbf{H}\_{\mathcal{C}C}^1 + \mathbf{H}\_{\mathcal{C}C}^2 \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{H}\_{\mathcal{A}C} \\ \mathbf{H}\_{\mathcal{C}C}^1 \\ \mathbf{H}\_{\mathcal{B}C} \end{bmatrix}^T \begin{bmatrix} \mathbf{F}\_A \\ \mathbf{F}\_C \\ \mathbf{F}\_B \end{bmatrix} \tag{43}
$$

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

$$\mathbf{K}\_{\diamond} \boldsymbol{\Phi}\_{\diamond} = \mathbf{M}\_{\diamond} \boldsymbol{\Phi}\_{\diamond} \boldsymbol{\Lambda}\_{\diamond} \tag{44}$$

**3.2. Parametric Reduced-Order Modeling (PROM) Method**

sults in a more accurate PROM algorithm.

**Figure 5.** Design space of three parameters.

corner points as

the above (*m* + 1) designs,

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 de‐ sign parameters. For a structure with *m* design variables, Zhang [9] suggested that the repre‐ sentative 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 re‐

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

http://dx.doi.org/10.5772/51402

149

The mode shapes of a new design are approximated in the space of the mode shapes of the

where the modal matrix **P** includes the basis vectors as in Equation (49) and *Θ* represents the participation factors of these vectors. The columns of P are the dominant mode shapes of

*<sup>p</sup>* =*PΘ* (48)

*P* = *Φ*<sup>0</sup> *Φ*<sup>1</sup> … *Φ<sup>m</sup>* (49)

*<sup>Φ</sup>* <sup>≈</sup>*Φ***˜**

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

$$\mathbf{K}\_p = \mathbf{K}\_0 + \Delta \mathbf{K} \qquad \mathbf{M}\_p = \mathbf{M}\_0 + \Delta \mathbf{M} \tag{45}$$

For a modest design change where Δ**K** and Δ**M** are small, it is assumed that the change in mode shapes is small and the new response can be therefore, captured in the sub-space spanned by the mode shapes **Φ**0 of the initial design. The new stiffness and mass matrices are then condensed as*K<sup>R</sup>* =*Φ*<sup>0</sup> *<sup>T</sup> <sup>K</sup> <sup>p</sup>Φ*0and *M<sup>R</sup>* <sup>=</sup>*Φ*<sup>0</sup> *<sup>T</sup> <sup>M</sup> <sup>p</sup>Φ*0and the following reduced eigenvalue problem is solved to calculate the eigen-vector *Θ*

$$\mathbf{K}\_p \Theta = \mathbf{M}\_p \Theta \Lambda\_p \tag{46}$$

The approximate eigenvalues of the new design are given by **Λ**<sup>p</sup> and the approximate eigen‐ vectors **Φ**p are

$$
\boldsymbol{\mathfrak{OP}}\_p = \mathbf{R}\boldsymbol{\Theta} \tag{47}
$$

where*R* =*Φ*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**

<sup>1</sup> <sup>1</sup> 1 12 1

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 struc‐ ture are **K**0 and **M**0, 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

For a modest design change where Δ**K** and Δ**M** are small, it is assumed that the change in mode shapes is small and the new response can be therefore, captured in the sub-space spanned by the mode shapes **Φ**0 of the initial design. The new stiffness and mass matrices

The approximate eigenvalues of the new design are given by **Λ**<sup>p</sup> and the approximate eigen‐

where*R* =*Φ*0. Thus, the modal response of the modified structure can be easily obtained and

*<sup>T</sup> <sup>K</sup> <sup>p</sup>Φ*0and *M<sup>R</sup>* <sup>=</sup>*Φ*<sup>0</sup>

**KK K MM M** *p p* = +D = +D 0 0 (45)

*<sup>T</sup> <sup>M</sup> <sup>p</sup>Φ*0and the following reduced eigen-

*K <sup>p</sup>Θ* =*M <sup>p</sup>ΘΛ<sup>p</sup>* (46)

*Φ<sup>p</sup>* = *RΘ* (47)

*A AA AC AC AC A B CA CC CC CC CC CC C C BB BC BC B*

æ ö éù é ùé ù é ù é ù ç ÷ êú ê úê ú ê ú ê ú = ++ é ù êú ê úê ú ê ú ê ú ë û

**X HH H H F X HH H HH H F X H H H F**

ëû ë ûë û ë û ë û è ø

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

**3. Reanalysis Methods for Dynamic Analysis**

êú ê

148 Advances in Vibration Engineering and Structural Dynamics

**3.1. CDH/VAO Method**

the following stiffness and mass matrices

are then condensed as*K<sup>R</sup>* =*Φ*<sup>0</sup>

vectors **Φ**p are

value problem is solved to calculate the eigen-vector *Θ*

the actual response can be recovered using the eigenvectors*Φp*.

*T*

(43)

(44)


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 de‐ sign parameters. For a structure with *m* design variables, Zhang [9] suggested that the repre‐ sentative 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 re‐ sults 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

$$\mathbf{D} \approx \widetilde{\mathbf{O}}\_p = \mathbf{P} \Theta \tag{48}$$

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

$$P = \begin{bmatrix} \Phi\_0 & \Phi\_1 & \dots & \Phi\_m \end{bmatrix} \tag{49}$$

where*Φ*0 is the modal matrix composed of the dominant mode shapes of the baseline de‐ sign, and *Φ<sup>i</sup>* is the modal matrix of the *<sup>i</sup> th* corner point. The mode shapes of the new design satisfy the following eigenvalue problem,

$$\mathbf{K}\,\widetilde{\mathbf{O}}\,\widetilde{\mathbf{O}}\,\_{p} = \mathbf{M}\,\widetilde{\mathbf{O}}\,\_{p}\Lambda \Leftrightarrow \mathbf{K}\mathbf{P}\Theta = \mathbf{M}\mathbf{P}\Theta\,\Lambda\tag{50}$$

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 <sup>T</sup>* as

$$\mathbf{K}\_R \Theta = \mathbf{M}\_R \Theta \Lambda \tag{51}$$

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 *Φ*0 of the baseline design and the dominant eigenvectors

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

, *i* =1, ..., *m* of the *m* corner design points (see Equation 49). When the size of basis *P* in‐ creases 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 ap‐ plied 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).

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 meth‐ od, when the number of reanalyses is less than or comparable to the number of design pa‐

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 origi‐ nal (baseline) design, the exact mode shapes *Φ*0 are obtained by solving the eigen-problem *K*0*Φ*<sup>0</sup> =*λ*0*M*0*Φ*0. We want to approximate the mode shapes of a modified design (subscript *p*) with stiffness and mass matrices *K <sup>p</sup>* = *K*<sup>0</sup> + *ΔK* and *M <sup>p</sup>* =*M*<sup>0</sup> + *ΔM* where *ΔK* and *ΔM* represent large perturbations. The CA method estimates the new eigenvalues *λp* and eigen‐

**3.3. Combined Approximations (CA) Method**

rameters, such as in gradient-based design optimization.

vectors *Φp* without performing an exact eigenvalue analysis.

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

*Φ<sup>p</sup>* =*λpK*<sup>0</sup>

*Φp*, *<sup>j</sup>* =*λpK*<sup>0</sup>

leading to the following recursive equation

(53) is replaced with *λpK*<sup>0</sup>

−1

−1

which produces a sequence of approximations of the mode shapes*Φp*, *<sup>j</sup>*

modes of the new design. In order to simplify the calculations, *λpK*<sup>0</sup>

*M <sup>p</sup>Φ<sup>p</sup>* − *K*<sup>0</sup>

*M <sup>p</sup>Φp*, *<sup>j</sup>*−<sup>1</sup> − *K*<sup>0</sup>

−1

CA method uses the changes *R <sup>j</sup>* =*Φp*, *<sup>j</sup>* −*Φp*, *<sup>j</sup>*−1 to form a subspace basis to approximate the

−1

−1

*ΔKΦ<sup>p</sup>* (53)

http://dx.doi.org/10.5772/51402

151

*ΔKΦp*, *<sup>j</sup>*−<sup>1</sup> (54)

−1

*M <sup>p</sup>Φ*0 and Equation (54) becomes

, *j* =1, 2, …, *s*. The

*M <sup>p</sup>Φp*, *<sup>j</sup>*−1in Equation

*Φi*

where the reduced stiffness and mass matrices are

$$\mathbf{K}\_R = \mathbf{P}^T \mathbf{K} \mathbf{P} \quad \text{and} \quad \mathbf{M}\_R = \mathbf{P}^T \mathbf{M} \mathbf{P} \tag{52}$$

Thus, the matrix *Θ* in Equation (48) consists of the eigenvectors of the reduced stiffness and mass matrices *KR* and*MR*.

For *m* design variables, (*m* + 1)eigenvalue problems must be solved in order to form the ba‐ sis *P* of Equation (49). Therefore, both the cost of obtaining the modal matrices *Φ<sup>i</sup>* 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:


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


The former is the fixed cost of PROM because it does not depend on the numbers of reanaly‐ ses 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 *Φ*0 of the baseline design and the dominant eigenvectors *Φi* , *i* =1, ..., *m* of the *m* corner design points (see Equation 49). When the size of basis *P* in‐ creases 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 ap‐ plied 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**

where*Φ*0 is the modal matrix composed of the dominant mode shapes of the baseline de‐ sign, and *Φ<sup>i</sup>* is the modal matrix of the *<sup>i</sup> th* corner point. The mode shapes of the new design

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

Thus, the matrix *Θ* in Equation (48) consists of the eigenvectors of the reduced stiffness and

For *m* design variables, (*m* + 1)eigenvalue problems must be solved in order to form the ba‐

size of matrix *P* increase linearly with *m*. The PROM approach uses the following algorithm

**1.** Calculate the mode shapes of the baseline design and the designs corresponding to the

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

**1.** the cost of performing (*m* + 1) full eigen-analyses to form subspace basis *P*in Equation

The former is the fixed cost of PROM because it does not depend on the numbers of reanaly‐ ses 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

**2.** Calculate the reduced stiffness and mass matrices *KR* and *MR* from Equation (52).

sis *P* of Equation (49). Therefore, both the cost of obtaining the modal matrices *Φ<sup>i</sup>*

*m* corner points in the design space, and form subspace basis*P*.

than the cost of a full analysis. The computational cost of PROM consists of

*<sup>p</sup>Λ*⇔*KPΘ* =*MPΘ Λ* (50)

*KRΘ* =*MRΘΛ* (51)

*<sup>p</sup>* using Equation (48).

and the

*<sup>K</sup><sup>R</sup>* <sup>=</sup>*<sup>P</sup> <sup>T</sup> KP* and *<sup>M</sup><sup>R</sup>* <sup>=</sup>*<sup>P</sup> <sup>T</sup> MP* (52)

satisfy the following eigenvalue problem,

150 Advances in Vibration Engineering and Structural Dynamics

*P <sup>T</sup>* as

*KΦ***˜**

where **Λ** is a diagonal matrix of the first *s* eigenvalues.

where the reduced stiffness and mass matrices are

to compute the mode shapes of a new design:

**3.** Solve eigenproblem (51) for matrix*Θ*.

**4.** Reconstruct the approximated eigenvectors in *Φ***˜**

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

mass matrices *KR* and*MR*.

(49), and

*<sup>p</sup>* <sup>=</sup>*<sup>M</sup> <sup>Φ</sup>***˜**

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 meth‐ od, when the number of reanalyses is less than or comparable to the number of design pa‐ rameters, 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 origi‐ nal (baseline) design, the exact mode shapes *Φ*0 are obtained by solving the eigen-problem *K*0*Φ*<sup>0</sup> =*λ*0*M*0*Φ*0. We want to approximate the mode shapes of a modified design (subscript *p*) with stiffness and mass matrices *K <sup>p</sup>* = *K*<sup>0</sup> + *ΔK* and *M <sup>p</sup>* =*M*<sup>0</sup> + *ΔM* where *ΔK* and *ΔM* represent large perturbations. The CA method estimates the new eigenvalues *λp* and eigen‐ vectors *Φp* without performing an exact eigenvalue analysis.

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

$$
\boldsymbol{\Phi}\_p = \lambda\_p \mathbf{K}\_0^{-1} \mathbf{M}\_p \boldsymbol{\Phi}\_p - \mathbf{K}\_0^{-1} \boldsymbol{\Delta} \mathbf{K} \boldsymbol{\Phi}\_p \tag{53}
$$

leading to the following recursive equation

$$
\boldsymbol{\upPhi}\_{p,j} = \boldsymbol{\uplambda}\_p \mathbf{K}\_0^{-1} \boldsymbol{\upmu}\_p \boldsymbol{\upPhi}\_{p,j-1} - \mathbf{K}\_0^{-1} \boldsymbol{\upmu} \mathbf{K} \, \boldsymbol{\upPhi}\_{p,j-1} \tag{54}
$$

which produces a sequence of approximations of the mode shapes*Φp*, *<sup>j</sup>* , *j* =1, 2, …, *s*. The CA method uses the changes *R <sup>j</sup>* =*Φp*, *<sup>j</sup>* −*Φp*, *<sup>j</sup>*−1 to form a subspace basis to approximate the modes of the new design. In order to simplify the calculations, *λpK*<sup>0</sup> −1 *M <sup>p</sup>Φp*, *<sup>j</sup>*−1in Equation (53) is replaced with *λpK*<sup>0</sup> −1 *M <sup>p</sup>Φ*0 and Equation (54) becomes *Φp*, *<sup>j</sup>* =*λpK*<sup>0</sup> −1 *M <sup>p</sup>Φ*<sup>0</sup> − *K*<sup>0</sup> −1 *ΔKΦp*, *<sup>j</sup>*−1 showing that the basis vectors satisfy the following recur‐ sive equation

$$\mathbf{R}\_{j} = -\mathbf{K}\_{0}^{-1} \Delta \mathbf{K} \mathbf{R}\_{j-1} \quad j = 2, \ldots, s \tag{55}$$

suitable for these problems because the subspace basis *R* does not change for every new de‐ sign point. Note that steps 1 and 3 (Equations 57and 58) are similar to steps 2 and 4 of

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

http://dx.doi.org/10.5772/51402

153

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

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 re‐ sponse 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 frequen‐ cy 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 devel‐ oped 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

The cost of calculating the subspace basis in Equation (55) consists of one matrix decomposi‐ tion (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 frequen‐ cy range of interest increases, more modes are needed to predict the structural response ac‐ curately. In such a case, although a single DCMP is needed in Equation (55), the number of columns in *Φ*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 calcu‐ lations 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 *Φ*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 *Φ*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

, *i* =1, ⋯, *m*, of the *m* corner design points, and second the cost of matrix condensation of Equation (57) is much lower than that of Equation (52), because the size (number of col‐ umns) of basis *R* is not proportional to the number of parameters *m* as in basis **P** . For prob‐ lems 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 pa‐ rameters. The CA method is therefore, more suitable than PROM, if the number of reanaly‐

PROM. CA uses basis *R* and PROM uses basis*P*.

ses is less than or comparable to the number of design parameters.

**3.4. Modified Combined Approximations (MCA) Method**

*Φi*

size models.

where the first basis vector is assumed to be*R*<sup>1</sup> = *K*<sup>0</sup> −1 *M <sup>p</sup>Φ*0.

The CA method forms a subspace basis

$$\mathbf{R} = \begin{bmatrix} \mathbf{R}\_1 & \mathbf{R}\_2 & \cdots & \mathbf{R}\_s \end{bmatrix} \tag{56}$$

where *s* is usually between 3 and 6 [16-18, 23] and the mode shapes of the new (*K <sup>p</sup>*,*M <sup>p</sup>*) de‐ sign are then approximated in the subspace spanned by *R* using the following algorithm:

**•** Condense the stiffness and mass matrices as

$$\mathbf{K}\_{\boldsymbol{\aleph}} = \mathbf{R}^{\boldsymbol{\aleph}} \mathbf{K}\_{\boldsymbol{\rho}} \mathbf{R} \qquad \quad \mathbf{M}\_{\boldsymbol{\aleph}} = \mathbf{R}^{\boldsymbol{\aleph}} \mathbf{M}\_{\boldsymbol{\rho}} \mathbf{R} \tag{57}$$


$$
\widetilde{\boldsymbol{\Phi}}\_p = \mathbf{R}\boldsymbol{\Theta} \tag{58}
$$

The eigenvalues of the new design are approximated by the eigenvalues *λ*˜ *<sup>p</sup>* of the reduced eigen-problem.

The CA method has three main advantages:


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 proba‐ bilistic 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 de‐ sign 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 *Φi* , *i* =1, ⋯, *m*, of the *m* corner design points, and second the cost of matrix condensation of Equation (57) is much lower than that of Equation (52), because the size (number of col‐ umns) of basis *R* is not proportional to the number of parameters *m* as in basis **P** . For prob‐ lems 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 pa‐ rameters. The CA method is therefore, more suitable than PROM, if the number of reanaly‐ ses is less than or comparable to the number of design parameters.

## **3.4. Modified Combined Approximations (MCA) Method**

*Φp*, *<sup>j</sup>* =*λpK*<sup>0</sup>

sive equation

matrix*Θ*.

eigen-problem.

−1

*M <sup>p</sup>Φ*<sup>0</sup> − *K*<sup>0</sup>

−1

152 Advances in Vibration Engineering and Structural Dynamics

where the first basis vector is assumed to be*R*<sup>1</sup> = *K*<sup>0</sup>

**•** Condense the stiffness and mass matrices as

The CA method has three main advantages:

calculate the subspace basis*R*,

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

*Φ***˜**

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

The CA method forms a subspace basis

1

*ΔKΦp*, *<sup>j</sup>*−1 showing that the basis vectors satisfy the following recur‐

0 1 2, , *j j j s* - **R K KR** =- D = - K (55)

**RRR R** = [ 1 2 L *<sup>s</sup>* ] (56)

*T T* **K RKR M RMR** *Rp R p* = = (57)

*<sup>p</sup>* as

*<sup>p</sup>* = *RΘ* (58)

−1 *M <sup>p</sup>Φ*0.

where *s* is usually between 3 and 6 [16-18, 23] and the mode shapes of the new (*K <sup>p</sup>*,*M <sup>p</sup>*) de‐ sign are then approximated in the subspace spanned by *R* using the following algorithm:

**•** Solve the reduced eigen-problem (using matrices **K**R and **M**R) to calculate the eigenvector

The eigenvalues of the new design are approximated by the eigenvalues *λ*˜ *<sup>p</sup>* of the reduced

**1.** it only requires a single matrix decomposition of stiffness matrix *K*0 in Equation (55) to

**3.** the eigenvectors of a new design are efficiently approximated by Equation (58) where

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 proba‐ bilistic analysis problems using Monte-Carlo simulation. The PROM method can be more

the eigenvectors *Θ* correspond to a much smaller reduced eigen-problem.

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 re‐ sponse 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 frequen‐ cy 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 devel‐ oped 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 decomposi‐ tion (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 frequen‐ cy range of interest increases, more modes are needed to predict the structural response ac‐ curately. In such a case, although a single DCMP is needed in Equation (55), the number of columns in *Φ*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 calcu‐ lations 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 *Φ*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 *Φ*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*<sup>0</sup> (stiffness matrix of base‐ line design) instead of *K <sup>p</sup>* (stiffness matrix of new design). The modified combined approxi‐ mations (MCA) method of this section addresses this issue.

The MCA method uses a subspace basis **T** whose columns are constructed using the recur‐ sive process

$$\begin{aligned} T\_1 &= \mathbb{K}\_p^{-1} \{ \mathbf{M}\_p \: \Phi\_0 \} \\ T\_i &= \mathbb{K}\_p^{-1} \{ \mathbf{M}\_p \: T\_{i-1} \} \qquad i = \text{2, 3, } \cdots \text{, s} \end{aligned} \tag{59}$$

**•** Calculate basis T using Equation (60) or Equation (61).

tions of the modes in the reduced space spanned by **T**

**•** Reconstruct the approximate eigenvectors *Φ***˜**

els with a large number of dominant modes.

Equation (61) is divided into *k* groups as

where

trices. However, the cost of estimating the mode shapes *Φ***˜**

*T i* = *Φ*<sup>0</sup> *<sup>i</sup> <sup>T</sup>*<sup>1</sup> *<sup>i</sup> <sup>T</sup>*<sup>2</sup>

the approximate mode shapes *Φ***˜**

**•** Calculate the condensed stiffness and mass matrices **K**R and **M**R as

*T T* **K TKT M TMT** *Rp R p* = = (62)

(*K<sup>R</sup>* −*λMR*)*Θ* =0 (63)

*<sup>p</sup>* are orthogonal with respect to the mass and stiffness ma‐

1 2 *<sup>k</sup>* <sup>=</sup> é ù ë û **T TT T** <sup>L</sup> (65)

*<sup>p</sup>*using Equations (62) to (64) may

http://dx.doi.org/10.5772/51402

155

*<sup>i</sup>* (66)

*<sup>p</sup>* =*TΘ* (64)

**•** Solve the following reduced eigen-problem to calculate the eigenvalues and the projec‐

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

*<sup>p</sup>* as

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 com‐ pensate for the increased cost of DCMP for dynamic reanalysis of large finite-element mod‐

All mode shapes in Equation (63) must be calculated simultaneously in order to ensure that

increase quickly (quadratically) with the number of modes, and as a result, the MCA meth‐ od 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 frequen‐ cy 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

*<sup>i</sup>* <sup>⋯</sup> *<sup>T</sup><sup>s</sup>*

*Φ***˜**

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 <sup>p</sup>* is inverted in the latter. The DCMP of *K <sup>p</sup>* 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 up‐ dating of the new design. If the process converges in *s* iterations, the mode shapes *Ts* will be the exact mode shapes*Φp*. Equation (55) does not have the same property. The vectors *T<sup>i</sup>* provide therefore, a more accurate approximation of the exact mode shapes *Φ<sup>p</sup>* than the *R<sup>i</sup>* vectors of the original CA method. This is an important advantage of MCA.

Because the mode shapes *T<sup>i</sup>* in 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

$$\mathbf{T} = \mathbf{T}\_1 \tag{60}$$

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

$$T = \begin{bmatrix} \Phi\_0 & T\_1 & T\_2 & \cdots & T\_s \end{bmatrix} \tag{61}$$

For better accuracy, the above basis also includes the mode shapes *Φ*0 of the baseline design. Because the approximate modes *T<sup>i</sup>* are more accurate than the CA vectors *R<sup>i</sup>* in approximat‐ ing 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 *Φ***˜** *<sup>p</sup>*as follows


$$\mathbf{K}\_{\mathcal{R}} = \mathbf{T}^{\top} \mathbf{K}\_{\rho} \mathbf{T} \qquad \quad \mathbf{M}\_{\mathcal{R}} = \mathbf{T}^{\top} \mathbf{M}\_{\rho} \mathbf{T} \tag{62}$$

**•** Solve the following reduced eigen-problem to calculate the eigenvalues and the projec‐ tions of the modes in the reduced space spanned by **T**

$$(\mathbf{K}\_R - \lambda \mathbf{M}\_R)\Theta = 0\tag{63}$$

**•** Reconstruct the approximate eigenvectors *Φ***˜** *<sup>p</sup>* as

$$
\widetilde{\boldsymbol{\Theta}}\_p = \mathbf{T} \boldsymbol{\Theta} \tag{64}
$$

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 com‐ pensate for the increased cost of DCMP for dynamic reanalysis of large finite-element mod‐ els 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 *Φ***˜** *<sup>p</sup>* are orthogonal with respect to the mass and stiffness ma‐ trices. However, the cost of estimating the mode shapes *Φ***˜** *<sup>p</sup>*using Equations (62) to (64) may increase quickly (quadratically) with the number of modes, and as a result, the MCA meth‐ od 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 frequen‐ cy 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

$$\mathbf{T} = \begin{bmatrix} \mathbf{T}^1 & \mathbf{T}^2 & \cdots & \mathbf{T}^k \end{bmatrix} \tag{65}$$

where

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*<sup>0</sup> (stiffness matrix of base‐ line design) instead of *K <sup>p</sup>* (stiffness matrix of new design). The modified combined approxi‐

The MCA method uses a subspace basis **T** whose columns are constructed using the recur‐

(*M <sup>p</sup>Ti*−1) *i* =2, 3, ⋯, *s*

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 <sup>p</sup>* is inverted in the latter. The DCMP of *K <sup>p</sup>* 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 up‐ dating of the new design. If the process converges in *s* iterations, the mode shapes *Ts* will be the exact mode shapes*Φp*. Equation (55) does not have the same property. The vectors *T<sup>i</sup>* provide therefore, a more accurate approximation of the exact mode shapes *Φ<sup>p</sup>* than the *R<sup>i</sup>*

in Equation (59) can quickly converge to the exact mode shapes

*T* = *Φ*<sup>0</sup> *T*<sup>1</sup> *T*<sup>2</sup> ⋯ *T<sup>s</sup>* (61)

**T T**= <sup>1</sup> (60)

(59)

in approximat‐

mations (MCA) method of this section addresses this issue.

154 Advances in Vibration Engineering and Structural Dynamics

*T*<sup>1</sup> = *K <sup>p</sup>* −1

*T<sup>i</sup>* = *K <sup>p</sup>* −1 (*M <sup>p</sup>Φ*0)

vectors of the original CA method. This is an important advantage of MCA.

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

Because the approximate modes *T<sup>i</sup>* are more accurate than the CA vectors *R<sup>i</sup>*

For better accuracy, the above basis also includes the mode shapes *Φ*0 of the baseline design.

ing 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

In summary, the proposed MCA method involves four steps in calculating the approximate

sive process

Because the mode shapes *T<sup>i</sup>*

[13-17].

eigenvectors *Φ***˜**

*<sup>p</sup>*as follows

$$\mathbf{T}^i = \begin{bmatrix} \Phi\_0^i & T\_1^i & T\_2^i & \cdots & T\_s^i \end{bmatrix} \tag{66}$$

Each group *T <sup>i</sup>* contains roughly *n* / *k* original modes *Φ*<sup>0</sup> *<sup>i</sup>* from*Φ*0, and their corresponding improved modes. The eigenvectors of the new design are calculated using *T <sup>i</sup>* 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.

**4. Reanalysis Methods in Dynamic Analysis and Optimization**

tically applicable.

**4.1. Integration of MCA Method in Optimization**

MCA and PROM methods is more suitable.

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. Depend‐ ing on the problem and the type of analysis, a particular reanalysis method may be prefer‐ red 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 ac‐ curacy 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 prac‐

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

http://dx.doi.org/10.5772/51402

157

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 probabil‐ istic analysis and gradient-free (e.g. genetic algorithms) optimization, a combination of the

Figure 6 shows a flowchart of the optimization process for modal frequency response prob‐ lems. The DMAP (Direct Matrix Abstraction Program) capabilities in NASTRAN have been used to integrate the MCA method and the NASTRAN modal dynamic response and opti‐ mization (SOL 200). The highlighted boxes indicate modifications to the NASTRAN opti‐ mizer. 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 practi‐ cal 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

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 re‐ analysis method can be inserted into the procedure as shown in Figure 6 to provide an ap‐ proximate 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

was reduced in half compared with the existing NASTRAN approach.

design variables and the number of expected iterations.

## **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 gradi‐ ent-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 employ‐ ing a Genetic Algorithm for example, and in simulation-based probabilistic analysis prob‐ lems 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.


**Table 1.** Comparison of the CA/MCA and PROM methods.
