# A reduction technique for mistuned bladed disks with superposition of large geometric mistuning and small model uncertainties

Vladislav Ganine, Mathias Legrand, Christophe Pierre, Hannah Michalska

► To cite this version:

Vladislav Ganine, Mathias Legrand, Christophe Pierre, Hannah Michalska. A reduction technique for mistuned bladed disks with superposition of large geometric mistuning and small model uncertainties. 12th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Feb 2008, Honolulu, United States. hal-01393687

HAL Id: hal-01393687

https://hal.science/hal-01393687

Submitted on 12 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
A REDUCTION TECHNIQUE FOR MISTUNED BLADED DISKS WITH
SUPERPOSITION OF LARGE GEOMETRIC MISTUNING AND SMALL
MODEL UNCERTAINTIES

<sup>a</sup>Vladislav Ganine,<sup>*</sup> <sup>a</sup>Mathias Legrand, <sup>a</sup>Christophe Pierre, and <sup>b</sup>Hannah Michalska

<sup>a</sup>Structural Dynamics & Vibration Laboratory, McGill University Montreal QC H3A 2K6, Canada

<sup>b</sup>Department of Electrical, & Computer Engineering, McGill University, Montreal QC H3A 2A7, Canada
vladislav.ganine@mail.mcgill.ca, mathias.legrand@mcgill.ca,
christophe.pierre@mcgill.ca, hannah.michalska@mcgill.ca

ABSTRACT

A new method for vibration analysis of mistuned bladed
disks is presented. The method combines two previously
reported modeling techniques in order to study the effects
of small random parameter variation on geometrically mis-
tuned bladed disks. It is based on the observation that the
nominal projections usually possess a certain degree of ro-
bustness tolerating small perturbations. Hence the subspace
spanned by a sufficient number of compensated tuned sys-
tem normal modes can be used to repeatedly project models
with small random parameter variation. The method is val-
idated numerically on an industrial bladed disk model, both
free and forced responses are compared with full model finite
element analysis. The benchmark cases show acceptable ac-
curacy while retaining low computational cost to build and
evaluate obtained ROM.

INTRODUCTION

Bladed disk assemblies belong to a class of rotationally
periodic systems in which cyclic symmetry is usually ex-
ploited to predict vibrational response because it reduces the
size of a problem to a single sector. However, the cyclic
symmetry is destroyed by blade-to-blade geometrical and/or
mechanical properties variations due to manufacturing tol-
erances and operational wear and tear, which is commonly
referred to as mistuning. The impact of mistuning on a pe-
riodic system sums down to strong mode localization and
significantly larger amplitudes of vibration response unlike

those of the nominal system and as a consequence their inad-
equate prediction when using cyclic symmetry. The stochas-
tic nature of blade to blade properties variations requires
Monte Carlo simulation to quantitatively assess the effect
of mistuning. In this connection model order reduction of
mistuned bladed disks plays an important role to save the
computational effort.

The problem of approximating high order finite element
models by lower order ones received a substantial amount
of interest over the years. Most of the reduction methods
for bladed rotors involve projection into a lower order sub-
space. Methods that follow that general structure [Castanier
et al., 1997, Bladh et al., 2001a, Bladh et al., 2001b, Bladh
et al., 2002, Lim et al., 2003, Lim et al., 2004, Yang and Grif-
fin, 2001] constitute modal truncation, where state transfor-
mation into a realization in which the state variables can be
ranked according to some measure of importance, i.e. nat-
ural frequencies falling into expected excitation frequency
band, followed by truncation of the least important ones.
Unlike others, [Petrov et al., 2002] used Woodbury formula
to calculate the inverse of modified by mistuning symmetric
FRF (Frequency Response Function) matrix.

Beyond the reduction approach the models can also
be classified in a way they introduce mistuning. Here the
substructuring of a system into blade and disk components
is a primary technique to conveniently introduce measur-
able blade to blade properties variations. [Castanier et al.,
1997] applied component mode synthesis technique mod-
ifying blade natural frequencies that appear explicitly in
the synthesized matrices. [Bladh et al., 2001a, Bladh et al.,

<footnote>*Address all correspondence to this author.</footnote>
2001b, Bladh et al., 2002] projected mistuning onto a selected set of tuned blade-alone modes, which allowed to capture nonproportional blade to blade variations in a computationally efficient way. [Lim et al., 2003] extended this approach and obtained a more compact ROM (Reduced Order Model) with CMM (Component Mode Mistuning) technique without component based representation of the reduced system where mistuning is projected directly onto a set of tuned system modes through mode participation factors. Although [Yang and Griffin, 2001] were first to express mistuned system normal modes as a linear combination of tuned ones, their approach required the knowledge of mistuned system stiffness/mass matrices. Thus they only simulated proportional mistuning by perturbing the Young's moduli of individual blades.

Even if using this projection yields accurate low order models, the number of tuned system modes required to capture the effects of large geometric mistuning (e.g. blade geometry change) render them computationally impractical. For that reason [Lim et al., 2004] developed a new SMC technique that generates models of commensurable size in the presence of large geometric mistuning. Later [Sinha, 2007] presented MMDA (Modified Modal Domain Analysis) method to deal with geometric mistuning problem yet his ROM are of significantly higher order - number of tuned system modes kept times number of POD features.

The current paper combines the results of two publications [Lim et al., 2003, Lim et al., 2004] and addresses the issue of how to efficiently carry out a Monte Carlo simulation in order to investigate the small mistuning-dependent performances with a reduced order model retaining the same set of mistuning parameters as those in the original FEM in the presence of deterministic geometric mistuning. We extend the original SMC (Static Mode Compensation) formulation to a non symmetric multiple blade geometry mistuning and demonstrate that CMM and SMC techniques can be effectively combined in a computationally inexpensive way to address the needs of design community.

This paper is organized as follows. The theoretical basis for the reduced order model formulation is presented in Section 2. In Section 3, the extensive numerical studies with an industrial bladed disk FEM are given. In particular, a numerical comparison with an exact solution using full FEM with different geometrical mistuning patterns is carried out (eigenvalues, mode shapes and forced response). The last section of the paper presents a comprehensive assessment of the numerical properties of the method proposed here.

## SELECTED NOMENCLATURE

* $\Lambda^S$ Diagonal matrix of nominal system eigenvalues

* $\mu_L^{syn}$ Reduced mass matrix with large mistuning only

* $\Phi^S$ Nominal modal matrix

* $\Phi_n^B$ Cantilevered-blade normal modes of the n-th blade

* $\Psi_n^{B,M}$ Modes that describe n-th blade motion at the boundary, mass mistuning

* $\omega_c^2$ Eigenvalue shift

* $\otimes$ Kronecker product

* $\mathbf{F}^Q$ Forcing terms due to mistuning

* $\mathbf{G}_{\Gamma_k}^Q$ k-th blade portion of equivalent forcing terms due to mistuning

* $\mathbf{K}$ Nominal stiffness matrix

* $\mathbf{M}$ Nominal mass matrix

* $\mathbf{M}_n^{\delta L}$ n-th blade mass matrix deviation due to large mistuning

* $\mathbf{M}_n^{\delta S}$ n-th blade mass matrix deviation due to small mistuning

* $\mathbf{q}_{\Phi,n}^M$ Participation factors of the cantilevered-blade modes in the n-th blade portion of compensated nominal modes, mass mistuning

* $\mathbf{q}_{\Psi,n}^M$ Participation factors of the boundary modes in the n-th blade portion of compensated nominal modes, mass mistuning

* $\tilde{\mathbf{B}}_{\text{diag}[\cdot]}$ A pseudo-block-diagonal matrix, where the size of $h=1,...,H$

* a block changes according to harmonic $h$

* $\hat{\Phi}$ Mistuned system modal matrix

* $\hat{\lambda}_i$ Mistuned system eigenvalue

* $\hat{\Psi}_{\Gamma}^{S,Q}$ Blade portion of nominal quasi-static modes matrix in cyclic coordinates

* $\mathcal{F}_k$ k-th row of the real-valued Fourier matrix

## THEORY

Consider a finite element model of an undamped structure around a position of equilibrium in frequency domain

$$(-\omega^2\mathbf{M} + \mathbf{K})\mathbf{x}(\omega) = \mathbf{f}(\omega) \quad (1)$$

where $\mathbf{x} = (x_1, ..., x_n)$ is the displacement vector of n degrees of freedom, $\mathbf{f} = (f_1, ..., f_n)$ is a vector of harmonic excitation forces, $\mathbf{M}$ and $\mathbf{K}$ are real symmetric definite mass and stiffness matrices, respectively.

The objective of the modal superposition reduction procedure is to determine an r-th order reduced-space basis, i.e. a set of modeshapes corresponding to contiguous set of natural frequencies within expected excitation frequency band, onto which the state vector can be projected

$$ \begin{gathered} (-\omega^2\mathbf{M}_r + \mathbf{K}_r)\mathbf{x}_r(\omega) = \mathbf{f}_r(\omega) \\ \mathbf{M}_r = \mathbf{\Phi}^{S^T}\mathbf{M}\mathbf{\Phi}^S, \qquad \mathbf{K}_r = \mathbf{\Phi}^{S^T}\mathbf{K}\mathbf{\Phi}^S \quad (2) \\ \mathbf{f}_r(\omega) = \mathbf{\Phi}^{S^T}\mathbf{f}(\omega), \qquad \mathbf{x}_r(\omega) = \mathbf{\Phi}^{S^T}\mathbf{x}(\omega) \end{gathered} $$

Typically, the modal truncation model reduction approach is accompanied by projection of the mistuned matrices onto the same truncated set of nominal system normal modes $\mathbf{\Phi}^S$,
which leads to a mistuned ROM eigenvalue problem:

$$ \mathbf{\Phi}^{ST}(\mathbf{K} + \mathbf{K}^{\delta})\mathbf{\Phi}^{S}y = \mathbf{\Lambda}^{m}\mathbf{\Phi}^{ST}(\mathbf{M} + \mathbf{M}^{\delta})\mathbf{\Phi}^{S}y \quad (3) $$

Although the procedure is relatively straightforward yielding very compact low order models when applied to bladed disk assemblies, it does not extend to the large geometric mistuning problem. Essentially, the precision in that case is limited by the amplitude of variation in mass and stiffness matrices $\mathbf{M}^{\delta}$ and $\mathbf{K}^{\delta}$ due to mistuning. If it is small compared to nominal properties, the reduced model natural frequencies and modeshapes will be close to those of the full system. It is well known [Stewart and Sun, 1990] that the spectrum of the matrix pencil

$$ (\mathbf{\Phi}^{ST}(\mathbf{K} + \mathbf{K}^{\delta})\mathbf{\Phi}^{S}, \mathbf{\Phi}^{ST}(\mathbf{M} + \mathbf{M}^{\delta})\mathbf{\Phi}^{S}) $$

perturbed from $(\mathbf{\Lambda}^S, \mathbf{I})$ is equivalent to that of the matrix pencil

$$ (\mathbf{\Phi}^{ST}(\mathbf{K} + \mathbf{K}^{\delta})\mathbf{\Phi}^{S} + O(\Delta^2), \mathbf{\Phi}^{ST}(\mathbf{M} + \mathbf{M}^{\delta})\mathbf{\Phi}^{S} + O(\Delta^2)) $$

up to the terms in $\Delta^2 = \|\mathbf{K}^{\delta}\|^2 + \|\mathbf{M}^{\delta}\|^2$. To address this issue, [Yang and Griffin, 2001] discuss the possibility of using the modes of a geometrically mistuned system as nominal modes to project additional random mistuning. Clearly, for very large systems that will be computationally infeasible, given that one must repeatedly solve the generalized eigenvalue problem for a full non symmetrical system. Here we present an alternative approach that does not require full system analysis: by using SMC method we calculate the compensated tuned system normal modes that accurately approximate those of a full system featuring large geometric mistuning, upon which we project small random mistuning using CMM.

## Algorithm Formulation

The combined model reduction algorithm can be summarized as follows:

1. Use SMC method to obtain a set of compensated basis vectors $(\mathbf{\Phi}_{\Gamma}^{\delta} - \mathbf{\Psi}^{S,Q}\mathbf{G}_{\Gamma}^{Q})$.

2. Project large mistuning $\mathbf{K}^{\delta L}$ and $\mathbf{M}^{\delta L}$ onto the subspace spanned by the basis vectors.

3. Calculate mode participation factors of CMM method using $(\mathbf{\Phi}_{\Gamma}^{\delta} - \mathbf{\Psi}^{S,Q}\mathbf{G}_{\Gamma}^{Q})$.

4. Project small mistuning $\mathbf{K}^{\delta S}$ and $\mathbf{M}^{\delta S}$.

5. Assemble total reduced order model mass and stiffness matrices $\boldsymbol{\mu}_{L}^{\text{syn}} + \boldsymbol{\mu}_{S}^{\text{syn}}$ and $\mathbf{k}_{L}^{\text{syn}} + \mathbf{k}_{S}^{\text{syn}}$.

6. Solve the reduced system.

## Calculation of the compensated basis vectors

[Lim et al., 2004] used Modal Acceleration Method to modify tuned system normal modes in such a way so they are closer to the mistuned system modes. In so doing one needs much fewer compensated modes to approximately span the eigenspace of geometrically mistuned system

$$ \mathbf{\Phi}^{S} - ((\mathbf{K} + \mathbf{K}^{\delta L}) - \omega_{c}^{2}(\mathbf{M} + \mathbf{M}^{\delta L}))^{-1}\mathbf{F}(\omega) = \sum_{i=1}^{N} \frac{(\omega^{2} - \omega_{c}^{2}) \hat{\mathbf{\Phi}}_{i}^{T} \mathbf{F}(\omega)}{\hat{\lambda}_{i} - \omega_{c}^{2}} \hat{\mathbf{\Phi}}_{i} \quad (4) $$

Notice that $\omega$ correspond to nominal system natural frequencies $\mathbf{\Lambda}^S$ and $\mathbf{F}(\omega)$ is an equivalent to geometric mistuning forcing excitation matrix defined as

$$ \mathbf{F}^{Q} = (-[\mathbf{\Lambda}^{S} - \omega_{c}^{2}\mathbf{I}][\mathbf{M} + \mathbf{M}^{\delta L}] + [\mathbf{K} + \mathbf{K}^{\delta L}] - \omega_{c}^{2}[\mathbf{M} + \mathbf{M}^{\delta L}])\mathbf{\Phi}^{S} = \left\{ \begin{array}{c} \mathbf{0} \\ (-[\mathbf{\Lambda}^{S} - \omega_{c}^{2}\mathbf{I}]\mathbf{M}^{\delta L} + \mathbf{K}^{\delta L} - \omega_{c}^{2}\mathbf{M}^{\delta L})\mathbf{\Phi}_{\Gamma}^{S} \end{array} \right\} \quad (5) $$

Rather than performing full matrix inversion, the quasi-static compensation terms can be obtained much more efficiently by taking advantage of the number of zero terms in the perturbation matrices $\mathbf{M}^{\delta L}$ and $\mathbf{K}^{\delta L}$

$$ (\mathbf{K} + \mathbf{K}^{\delta L} - \omega_c^2 (\mathbf{M} + \mathbf{M}^{\delta L}))^{-1} = [\mathbf{K} - \omega_c^2 \mathbf{M}]^{-1} (\mathbf{I} + [\mathbf{K}^{\delta L} - \omega_c^2 \mathbf{M}^{\delta L}][\mathbf{K} - \omega_c^2 \mathbf{M}^{\delta L}]^{-1})^{-1} \quad (6) $$

Then the equivalent force can be expressed as

$$ \mathbf{G}_{\Gamma}^{Q} = \left\{ \begin{array}{c} \mathbf{0} \\ (\mathbf{I} + [\mathbf{K}^{\delta L} - \omega_c^2 \mathbf{M}^{\delta L}]\mathbf{\Psi}_{\Gamma}^{S,Q})^{-1} \mathbf{F}_{\Gamma}^{Q} \end{array} \right\} \quad (7) $$

In order to efficiently calculate the quasi-static modes $\mathbf{\Psi}^{S,Q}$, we use the cyclic symmetry domain, where $[\mathbf{K} - \omega_c^2 \mathbf{M}]^{-1}$ is inverted blockwise

$$ \tilde{\mathbf{\Psi}}_{\Gamma}^{S,Q} = \tilde{\mathbf{B}}\mathrm{diag} \left[ (\tilde{\mathbf{K}}_h - \omega_c^2 \tilde{\mathbf{M}}_h)^{-1} \right] \quad (8) $$

Clearly, in physical domain we only need those partitions of $\mathbf{\Psi}^{S,Q}$ that correspond to mistuned DOF of each blade

$$ \mathbf{\Psi}_{k,m}^{S,Q} = (\mathcal{F}_k \otimes \mathbf{I}) \tilde{\mathbf{\Psi}}_{\Gamma}^{S,Q} (\mathcal{F}_m^T \otimes \mathbf{I}) \quad (9) $$

Tuned blades are numbered $k = 1, \dots, N$, mistuned $m = 1, \dots, K$, where $N$ and $K$ denoting the total number and the number of mistuned blades correspondingly. Note that $\mathcal{F}_k$ is referred to $k$-th row of real-valued Fourier matrix.
## Large mistuning projection

The reduced order model mass and stiffness matrices due to geometric mistuning are given by

$$ 
\begin{aligned}
\boldsymbol{\mu}_{L}^{\text{syn}} &= (\boldsymbol{\Phi}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q})^{T} (\mathbf{M} + \mathbf{M}^{\delta_{L}}) (\boldsymbol{\Phi}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q}) \\
&= \mathbf{I} + \boldsymbol{\mu}_{1}^{\text{syn}} - \boldsymbol{\mu}_{2}^{\text{syn}} - \boldsymbol{\mu}_{2}^{\text{syn}T} + \boldsymbol{\mu}_{3}^{\text{syn}} \\
\boldsymbol{\kappa}_{L}^{\text{syn}} &= (\boldsymbol{\Phi}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q})^{T} (\mathbf{K} + \mathbf{K}^{\delta_{L}}) (\boldsymbol{\Phi}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{k}^{Q}) \\
&= \mathbf{I} + \boldsymbol{\kappa}_{1}^{\text{syn}} - \boldsymbol{\kappa}_{2}^{\text{syn}} - \boldsymbol{\kappa}_{2}^{\text{syn}T} + \boldsymbol{\kappa}_{3}^{\text{syn}}
\end{aligned}
\quad (10) $$

Those products can be taken on a sector-by-sector basis i.e.,

$$ 
\begin{aligned}
\boldsymbol{\mu}_{1}^{\text{syn}} &= \sum_{k=1}^{K} \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{M}_k^{\delta_L} \boldsymbol{\Phi}_{\Gamma_k}^{S} \\
\boldsymbol{\mu}_{2}^{\text{syn}} &= \sum_{k=1}^{K} \left( (\boldsymbol{\Lambda}^{S} - \omega_c^2 \mathbf{I})^{-1} \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{G}_{\Gamma_k}^{Q} + \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{M}_k^{\delta_L} \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right) \\
\boldsymbol{\mu}_{3}^{\text{syn}} &= \sum_{k=1}^{K} \left( \left[ \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right]^T \mathbf{M}_k^{\delta_L} \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right. \\
&\qquad \left. + \sum_{m=1}^{K} \mathbf{G}_{\Gamma_m}^{Q,T} \left[ \boldsymbol{\Psi}^{S,Q,T} \mathbf{M} \boldsymbol{\Psi}^{S,Q} \right]_{m,k} \mathbf{G}_{\Gamma_k}^{Q} \right) \\
\boldsymbol{\kappa}_{1}^{\text{syn}} &= \sum_{k=1}^{K} \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{K}_k^{\delta_L} \boldsymbol{\Phi}_{\Gamma_k}^{S} \\
\boldsymbol{\kappa}_{2}^{\text{syn}} &= \sum_{k=1}^{K} \left( \boldsymbol{\Lambda}^{S} (\boldsymbol{\Lambda}^{S} - \omega_c^2 \mathbf{I})^{-1} \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{G}_{\Gamma_k}^{Q} + \boldsymbol{\Phi}_{\Gamma_k}^{S,T} \mathbf{K}_k^{\delta_L} \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right) \\
\boldsymbol{\kappa}_{3}^{\text{syn}} &= \sum_{k=1}^{K} \left( \mathbf{G}_{\Gamma_k}^{Q} \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right. \\
&\qquad \left. + \omega_c^2 \left[ \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right]^T \mathbf{M}_k^{\delta_L} \sum_{m=1}^{K} \boldsymbol{\Psi}_{k,m}^{S,Q} \mathbf{G}_{\Gamma_k}^{Q} \right. \\
&\qquad \left. + \sum_{m=1}^{K} \mathbf{G}_{\Gamma_m}^{Q,T} \left[ \boldsymbol{\Psi}^{S,Q,T} \mathbf{M} \boldsymbol{\Psi}^{S,Q} \right]_{m,k} \mathbf{G}_{\Gamma_k}^{Q} \right)
\end{aligned}
$$

Again, $\boldsymbol{\Psi}^{S,Q,T} \mathbf{M} \boldsymbol{\Psi}^{S,Q}$ terms are calculated using cyclic symmetry in the similar vein

$$ \tilde{\boldsymbol{\Psi}}_{\Gamma}^{S,Q,T} \tilde{\mathbf{M}} \tilde{\boldsymbol{\Psi}}_{\Gamma}^{S,Q} = \tilde{\text{Bdiag}} [\tilde{\boldsymbol{\Psi}}_{\Gamma,h}^{S,Q,T} \tilde{\mathbf{M}}_{h} \tilde{\boldsymbol{\Psi}}_{\Gamma,h}^{S,Q}] \quad (11) $$

$$ (\boldsymbol{\Psi}^{S,Q,T} \mathbf{M} \boldsymbol{\Psi}^{S,Q})_{k,m} = (\mathcal{F}_k \otimes \mathbf{I})(\tilde{\boldsymbol{\Psi}}_{\Gamma}^{S,Q,T} \tilde{\mathbf{M}} \tilde{\boldsymbol{\Psi}}_{\Gamma}^{S,Q})(\mathcal{F}_m^T \otimes \mathbf{I}) \quad (12) $$

## Small mistuning projection

Following calculation of basis vectors and contribution of large geometric mistuning to the reduced matrices, we can obtain the total reduced mass and stiffness matrices $\boldsymbol{\mu}^{\text{syn}}$ and

$\boldsymbol{\kappa}^{\text{syn}}$ by projecting small random blade-to-blade variations with CMM method [Lim et al., 2003]

$$ 
\begin{aligned}
\boldsymbol{\mu}^{\text{syn}} &= \boldsymbol{\mu}_{L}^{\text{syn}} + (\boldsymbol{\Phi}_{\Gamma}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q})^{T} \mathbf{M}^{\delta_{S}} (\boldsymbol{\Phi}_{\Gamma}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q}) \\
&= \boldsymbol{\mu}_{L}^{\text{syn}} + \sum_{n=1}^{N} \mathbf{q}_{n}^{M} \mathbf{U}_{n}^{M} \mathbf{M}_{n}^{\delta_{S}} \mathbf{U}_{n}^{M} \mathbf{q}_{n}^{M} \\
\boldsymbol{\kappa}^{\text{syn}} &= \boldsymbol{\kappa}_{L}^{\text{syn}} + (\boldsymbol{\Phi}_{\Gamma}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q})^{T} \mathbf{K}^{\delta_{S}} (\boldsymbol{\Phi}_{\Gamma}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q}) \\
&= \boldsymbol{\kappa}_{L}^{\text{syn}} + \sum_{n=1}^{N} \mathbf{q}_{n}^{K} \mathbf{U}_{n}^{K} \mathbf{K}_{n}^{\delta_{S}} \mathbf{U}_{n}^{K} \mathbf{q}_{n}^{K}
\end{aligned}
\quad (13) $$

where the blade portion of compensated normal modes is expressed as

$$ (\boldsymbol{\Phi}_{\Gamma}^{S} - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{\Gamma}^{Q}) = \tilde{\text{Bdiag}} [\mathbf{U}_n] \mathbf{q} \quad (14) $$

Note that matrices $\mathbf{U}_n$ and modal participation factors $\mathbf{q}_n$, which represent the motion of n-th blade in the coordinates of basis vectors and described by cantilevered-blade normal modes together with boundary modes, are calculated for both mass and stiffness mistuning

$$ \mathbf{U}_n^M = \begin{bmatrix} \mathbf{\Phi}_n^B & \mathbf{\Psi}_n^{B,M} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} = \begin{bmatrix} \mathbf{\Phi}_n^B - (\mathbf{M}_{n,\Gamma_{ii}} + \mathbf{M}_{n_{ib}}^{\delta_L})^{-1}(\mathbf{M}_{n,\Gamma_{ib}} + \mathbf{M}_{n_{ib}}^{\delta_L}) \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \quad (15) $$

$$ \mathbf{U}_n^K = \begin{bmatrix} \mathbf{\Phi}_n^B & \mathbf{\Psi}_n^{B,K} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} = \begin{bmatrix} \mathbf{\Phi}_n^B - (\mathbf{K}_{n,\Gamma_{ii}} + \mathbf{K}_{n_{ib}}^{\delta_L})^{-1}(\mathbf{K}_{n,\Gamma_{ib}} + \mathbf{K}_{n_{ib}}^{\delta_L}) \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \quad (16) $$

where the subscripts i and b indicate that each n-th blade component is partitioned into internal and boundary DOFs.

$$ \mathbf{q}_n^M = \begin{bmatrix} \mathbf{q}_{\Phi,n}^M \\ \mathbf{q}_{\Psi,n}^M \end{bmatrix} = \begin{bmatrix} \boldsymbol{\Phi}_n^{B,T} (\mathbf{M}_{n,\Gamma} + \mathbf{M}_n^{\delta_L}) (\boldsymbol{\Phi}_{n,\Gamma}^S - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{n,\Gamma}^Q) \\ (\boldsymbol{\Phi}_{n,\Gamma}^S - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{n,\Gamma}^Q)_b \end{bmatrix} \quad (17) $$

$$ \mathbf{q}_n^K = \begin{bmatrix} \mathbf{q}_{\Phi,n}^K \\ \mathbf{q}_{\Psi,n}^K \end{bmatrix} = \begin{bmatrix} \boldsymbol{\Lambda}_n^{B^{-1}} \boldsymbol{\Phi}_n^{B^T} (\mathbf{K}_{n,\Gamma} + \mathbf{K}_n^{\delta_L}) (\boldsymbol{\Phi}_{n,\Gamma}^S - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{n,\Gamma}^Q) \\ (\boldsymbol{\Phi}_{n,\Gamma}^S - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{n,\Gamma}^Q)_b \end{bmatrix} \quad (18) $$

subscript b denotes blade boundary DOF partition of compensated modes $(\boldsymbol{\Phi}_{n,\Gamma}^S - \boldsymbol{\Psi}^{S,Q} \mathbf{G}_{n,\Gamma}^Q)$. Unlike CMM based on tuned system nominal modes, the cantilevered-blade normal and boundary modes of each blade with rogue geometry must also be precalculated in addition to those of a tuned blade.
## SMC Convergence

The CMM concept extends readily to the combined large plus small mistuning case; however, it is important to accurately approximate the eigenspace of the system undergoing large mistuning before any small mistuning projection. Therefore, it is useful to revisit the classical MAM formulation given in [Akgun, 1993]. If the mid-frequency truncation scheme with eigenvalue shifting is applied, the FRF can be expressed as

$$x(\omega) \approx \left( \sum_{j=1}^{J} \left( (\omega^2 - \omega_c^2)(\mathbf{K} - \omega_c^2 \mathbf{M})^{-1} \mathbf{M} \right)^{j-1} (\mathbf{K} - \omega_c^2 \mathbf{M})^{-1} f(\alpha_i) \right) + \sum_{i=N_{\text{low}}}^{N_{\text{high}}} \left( \frac{\omega^2 - \omega_c^2}{\lambda_i - \omega_c^2} \right)^J \frac{\mathbf{\Phi}_i^T \mathbf{f}(\omega)}{\lambda_i - \omega_c^2} \mathbf{\Phi}_i \quad (19)$$

where J is a positive integer denoting the level of acceleration, $N_{low}$ and $N_{high}$ are indexes of the first and last kept modes, and the frequency shift is usually taken as $\omega_c^2 = (\lambda_{N_{low}} - \lambda_{N_{high}})/2$, must satisfy $\omega_c^2 \neq \lambda_i$ and, in our case, $\omega_c^2 \neq \omega$. Clearly, the FRF $x(\omega)$ obtained from (19) reach the exact values when the number of modes included increases. Although the accuracy of the FRF may also be improved by increasing the level of acceleration J, it seems impractical computationally. Since we do not know how many modes are enough to compute the FRFs to the prescribed accuracy, the iterative scheme becomes necessary. Also, a proper selection of $\omega_c^2$ closer to some harder to approximate $\lambda_i$ could help to accelerate the convergence to that eigenpair.

## NUMERICAL STUDIES

The theory presented in the last section was implemented in FORTRAN code using high performance Intel MKL numerical library. The program incorporates mistuned system modes calculation and forced response analysis. In this section the performance of the algorithm will be illustrated using an industrial bladed disk model. First the natural frequencies and mode shapes of geometrically mistuned disc are analyzed. Second a realistic small mistuning is added on top of geometrical mistuning and forced response is calculated. Both mistuned modes and dynamic response are also computed with a full finite element model using the commercially available ANSYS FE code, which serves as a benchmark to assess the accuracy of ROM solution.

### Large Mistuning

Consider the finite element model of an industrial 29 blade rotor depicted in Fig. 1. The model is built with standard linear eight node brick elements totaling 126846 DOF. In order to provide a simple test case the bladed disk was mistuned by significantly changing geometry of three arbitrarily selected blades with the assumption that mass density and Young's modulus are unchanged. We construct ROM

<figure><img src="image_2.png" /><figcaption>Figure 1. Finite element model of mistuned bladed disk.</figcaption></figure>

for two frequency ranges 33 – 36 KHz and 14.5 – 16.5 KHz. The first band corresponds to 3T torsion and 2S stripe mode families, and the second to 2T/2F mixed torsion - flexural bending family of modes, which is shown in Fig. 2. The frequency bands of interest are marked by two horizontal lines.

<figure><img src="image_3.png" /><figcaption>Figure 2. Natural frequencies versus nodal diameters for the tuned FEM</figcaption></figure>

The selection of 2S and 3T families is a good candidate for a benchmark test because the amount of interaction between them would typically require to include both in ROM. 2T/2F mode family is relatively isolated, therefore it would provide an alternative test case for ROM.

To build the reduced order model, the following finite
element models are required:

- A single sector 4372 DOF model of the tuned system on which we run ANSYS modal analysis in cyclic symmetry and obtain system normal modes with harmonic blocks of mass and stiffness matrices in real-cyclic coordinates.

- Tuned blade 2496 DOF model on which we run ANSYS modal analysis to obtain the cantilevered-blade modes (used later when small mistuning is projected) and blade alone mass and stiffness matrices.

- For each blade with geometrical mistuning we run ANSYS modal analysis on 2496 DOF model to obtain the cantilevered-blade modes (used later to project small mistuning) and mistuned blade mass and stiffness matrices.

Fig. 3 and Fig. 4 depict the mistuned frequencies errors and MAC ratio between mistuned mode shapes predicted by ROM and ANSYS for the 33-36 KHz frequency band calculated with $\omega_c = 32400$ KHz. 71 natural frequencies that exist in that region are chosen for comparison with full FEM results. Note that even for the smallest model of 71 DOF

<figure><img src="image_4.png" /><figcaption>Figure 3. Mistuned natural frequencies errors for 33 – 36 KHz region</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Comparison of MAC ratio between FEM and ROM mistuned mode shapes for 33 – 36 KHz region</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Highly localized modes.</figcaption></figure>

the natural frequency error for all approximated modes is below 0.04% and MAC value above 0.995, which demonstrates that 71 DOF ROM is sufficiently accurate to capture the dynamics of mistuned system in the range 33 – 36 KHz. The models of greater size are shown to demonstrate the convergence towards FEM results. Here the modal basis was gradually augmented by including $3F$ and $4F/1R$ families of modes, which resulted in ROM of the size up to 136 DOF. Problematic for the proper development of an iterative algorithm is the fact that in some cases the asymmetric inclusion of mode families might potentially increase the error of approximation for the modes closer to the opposite side the studied frequency band, which is the case for the models of 102 and 105 DOF.

It should also be pointed out that the modes with higher error of approximation and slower convergence rate are those strongly localized to a single blade. Some of them that appear in 33 – 36 KHz frequency range are depicted in Fig. 5. In particular note 41-st mode located in the middle of the frequency range and localized to 7-th blade. The errors at the end of the analyzed frequency range are typical for any modal projection approach with modal truncation [Yang and Griffin, 2001] because the modes in the middle are more accurately represented than those at the end.

Fig. 6 and Fig. 7 depict the natural frequency errors and MAC ratio in the 14.5-16.5 KHz frequency band calculated with $\omega_c = 15390$ KHz. 29 modes are chosen for comparison with full FEM results. The error in the frequencies for the smallest model is less than 0.01 percent and the worst MAC ratio is 0.9995. Again the maximum error was observed for the 29-th mode, which is a highly localized to 7-th blade mode. Clearly, the error has been significantly reduced by
<figure><img src="image_7.png" /><figcaption>Figure 6. Mistuned natural frequencies errors for 14.5 – 16.4 KHz region</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Comparison of MAC ratio between FEM and ROM mistuned mode shapes for 14.5 – 16.5 KHz region</figcaption></figure>

increasing the number of modes included in the basis. The important point to note is that SMC better approximates the lower frequency modes and those of belonging to an isolated family, which results in very compact models of the number of blades order.

## Combined Large and Small Mistuning

In this section SMC and CMM techniques are used to compute the dynamic response of the geometrically mistuned bladed disk system, Fig. 1, with the nominal Young's modulus of the *n*-the blade mistuned as

$$E_n = E_0(1 + \delta_n^e) \quad (20)$$

where $E_0$ is the nominal Young's modulus and $\delta_n^e$ is a non-dimensional mistuning value. The specific pattern used in this case is shown in Tab. 1 Although the proportional blade to blade stiffness variation introduced by changing Young's modulus of the blades is a very rough way to model random small mistuning phenomena, it has been chosen for simplicity to validate the accuracy of the developed method. In this study the structural damping coefficient is set to 0.006, the modal participation factors are calculated using

<table><caption>Table 1. Eigenvalue mistuning pattern</caption><thead><tr><th>Blade</th><th>δ<sub>n</sub><sup>e</sup></th><th>Blade</th><th>δ<sub>n</sub><sup>e</sup></th></tr></thead><tbody><tr><td>1</td><td>0.05704</td><td>16</td><td>0.04934</td></tr><tr><td>2</td><td>0.01207</td><td>17</td><td>0.04479</td></tr><tr><td>3</td><td>0.04670</td><td>18</td><td>0.03030</td></tr><tr><td>4</td><td>-0.01502</td><td>19</td><td>0.00242</td></tr><tr><td>5</td><td>0.05969</td><td>20</td><td>0.01734</td></tr><tr><td>6</td><td>-0.03324</td><td>21</td><td>0.02919</td></tr><tr><td>7</td><td>-0.00078</td><td>22</td><td>-0.00328</td></tr><tr><td>8</td><td>-0.01688</td><td>23</td><td>0.00086</td></tr><tr><td>9</td><td>0.00242</td><td>24</td><td>-0.03654</td></tr><tr><td>10</td><td>-0.02747</td><td>26</td><td>-0.03631</td></tr><tr><td>11</td><td>-0.03631</td><td>26</td><td>-0.01665</td></tr><tr><td>12</td><td>-0.03570</td><td>27</td><td>0.00783</td></tr><tr><td>13</td><td>-0.03631</td><td>28</td><td>-0.01169</td></tr><tr><td>14</td><td>-0.03631</td><td>29</td><td>-0.01332</td></tr><tr><td>15</td><td>0.00242</td><td></td><td></td></tr></tbody></table>

30 cantilevered-blade modes, with only a third of them dominant. Despite the large changes in geometry the aerodynamic mistuning is neglected; the forcing function applied corresponds to the 3 and 5 nodal diameter excitation for 71 DOF ROM in the 33 – 36 KHz region, and 2 and 3 engine order excitation for 29 DOF ROM in the 14.6 – 16.6 KHz region. An arbitrary load is applied at a tip of each blade. Figures 8 and 9 show the Euclidean displacement norms

<figure><img src="image_9.png" /><figcaption>Figure 8. Forced response in 33 – 36 KHz region due to engine order 5 excitation</figcaption></figure>

for maximum responding blade versus frequency of excitation for 33 – 36 KHz region EO 5 and EO 3 excitation respectively. The results for 14.5 – 16.5 KHz region under EO 3 and EO 2 excitation are depicted in Fig. 10 and Fig. 11. ANSYS results are marked as FEM L+S and those of re-
<figure><img src="image_10.png" /><figcaption>Figure 9. Forced response in 33 – 36 KHz region due to engine order 3 excitation</figcaption></figure>

duced model as ROM L+S. The responses of the tuned sys-
tem denoted as FEM T and geometrically mistuned system
denoted as ROM L are presented for comparison. Observe

<figure><img src="image_11.png" /><figcaption>Figure 10. Forced response in 14.5 – 16.5 KHz region due to engine order 3 excitation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Forced response in 14.5 – 16.5 KHz region due to engine order 2 excitation</figcaption></figure>

that the ROM in both frequency bands provide an accurate
prediction of the mistuned systems forced response to engine
order excitation in the tuned system veering regions as com-
pared to the results computed by ANSYS. Even at the natural
frequencies of slower to converge highly localized modes 71
and 29 DOF ROM are able to match FEM maximum forced
response.

DISCUSSION

In developing ROM there is an initial computational
overhead of calculating tuned system quasi-static modes
**Ψ**<sup>S,Q</sup> that are used as input to SMC reduction algorithm.
Since the they need to be calculated only once for the se-
lected frequency band and mistuned DOF distribution, the
real computational advantage in using SMC occurs when
simulating a large number of geometrically mistuned bladed
disks either to analyze physical blade damage phenomena
or for design optimization purposes. Thus we must solve
*N*<sub>harmonics</sub> times a symmetric sparse system of the order of
2*N*<sub>DOF</sub> of a basic sector. Additionally there is a small cost
associated with calculating **Ψ**<sup>S,Q<sup>T</sup></sup> **M****Ψ**<sup>S,Q</sup> matrix product for
each harmonic. For the example presented here with 4372
DOF per sector model the time of calculation on a 2.66 GHz
XEON Quad-Core with 4 GB of RAM is of the order of 8
minutes, however this time is dominated by disk I/O oper-
ations and is amenable to further optimization. The deriva-
tion of compensated basis vectors incurs the solution of a
fully populated system in equation (7) of the order of total
number of geometrically mistuned DOF. If a large number of
DOF need to be mistuned, then the SMC reduction approach
will be prohibitively expensive. The projection of large mis-
tuning is carried out in an efficient, potentially paralleliz-
able way, because the perturbation matrices corresponding
to the geometric mistuning of blades have block-diagonal
structure. The total time of calculation of 71 geometrically
mistuned system mode shapes for the test case model with
3 mistuned blades, 594 mistuned DOF per blade, is 45 sec-
onds. The same task with FEM requires ~ 5 minutes of AN-
SYS analysis running on the same platform, which demon-
strates the potential of a substantial analysis cost reduction.

Depending on the applied geometric mistuning and the
frequency band of interest we can obtain a very compact
ROM of the order of number of blades. If system dynam-
ics cannot be accurately captured by a reasonable number of
modes while the smallest possible model is desired, it might
be more efficient to first solve the reduced system to find
mode shapes of a geometrically mistuned system and then
use the approximated modes in CMM projection.

Throughout the algorithm, there are several decisions to
be made. First, the centering frequency ω_c must be carefully
selected not to cross the tuned system and potential geomet-
rically mistuned system natural frequencies. Second, the
number of cantilevered-blade normal modes to accurately
represent small mistuning have to be determined given the
tuned and mistuned blades geometry. Finally, we have to
know how many compensated normal modes should be re-
tained in the reduced order model based on the expected er-
ror of approximation. Work is currently underway to develop
an inexpensive iterative approach to better select the projec-
tion basis as well as to assess the perturbation bounds.

CONCLUSION

Motivated by the practical needs of the design commu-
nity to assess the effect of random mistuning on geometri-
cally mistuned bladed disks a new reduced order modeling
method has been presented. By the use of two previously
reported techniques, SMC and CMM, compact ROM were
efficiently calculated that accurately capture the mistuned
system dynamics with a very low number of states over a
broad range of excitation frequencies and different geome-
try mistuning patterns. It was demonstrated with numerical
examples that this method can accurately predict mistuned
frequencies, mistuned mode shapes and forced response of
a bladed disk with non symmetric geometric mistuning and
simultaneous small perturbations in both mass and stiffness
matrices of blades. The required amount of computation is
significantly lower than that of the projection method using
nominal geometrically mistuned system modes calculated by
FEM analysis. In this light the developed numerical tool
is very promising for future investigations of interaction of
large deterministic and small random mistuning, as well as
intentional mistuning.

REFERENCES

[Akgun, 1993] Akgun, M. (1993). A new family of mode-
superposition methods for response calculations. *Journal*
*of Sound and Vibration*.

[Bladh et al., 2001a] Bladh, R., Castanier, M., and Pierre, C. (2001a). Component-mode-based reduced order modeling techniques for mistuned bladed disks - part 1: Theoretical models. *ASME Journal of Engineering for Gas Turbines and Power*.

[Bladh et al., 2001b] Bladh, R., Castanier, M., and Pierre, C. (2001b). Component-mode-based reduced order modeling techniques for mistuned bladed disks - part 2: Ap-plicaiton. *ASME Journal of Engineering for Gas Turbines and Power*.

[Bladh et al., 2002] Bladh, R., Castanier, M., Pierre, C., and Kruse, M. (2002). Dynamic response predictions for a mistuned industrial turbomachinery rotor using reduced order modeling. *ASME Journal of Engineering for Gas Turbines and Power*.

[Castanier et al., 1997] Castanier, M., Ottarsson, G., and Pierre, C. (1997). A reduced-order modeling technique for mistuned bladed disks. *Journal of Vibration and Acustics*.

[Lim et al., 2003] Lim, S., Bladh, R., Castanier, M., and

Pierre, C. (2003). A compact, generalized com-
ponent mode mistuning representation for modeling
bladed disk vibration. In *Proceedings of the 45-th
AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dy-
namics and Material Conference, paper 2003-1545*,
Palm Springs, California, USA.

[Lim et al., 2004] Lim, S., Castanier, M., and Pierre, C. (2004). Vibration modeling of bladed disks subject to geometric mistuning and design changes. In *Proceedings of the 45-th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material Conference, paper 2004-1686*, Palm Springs, California, USA.

[Petrov et al., 2002] Petrov, E., Sanliturk, K., and Ewins, D. (2002). A new method for dynamic analysis of mistuned bladed disks based on the exact relationship between tuned and mistuned systems. *ASME Journal of Engineering for Gas Turbines and Power*.

[Sinha, 2007] Sinha, A. (2007). Reduced-order model of a bladed rotor with geometric mistuning.

[Stewart and Sun, 1990] Stewart, G. and Sun, J. (1990). *Matrux Perturbation Theory*. Academic Press, San Diego CA.

[Yang and Griffin, 2001] Yang, M. and Griffin, J. (2001). A reduced order model of mistuning using a subset of nominal system modes. *ASME Journal of Engineering for Gas Turbines and Power*.