# Nonlinear Cyclic Transient Dynamic Analysis for Bladed Disk Tip Deflection

Tong Hui, Letian Wang, Suryarghya Chakrabarti, Siddharth Ashar, Steven Szczap, Anthony Laude

► To cite this version:

Tong Hui, Letian Wang, Suryarghya Chakrabarti, Siddharth Ashar, Steven Szczap, et al.. Nonlinear Cyclic Transient Dynamic Analysis for Bladed Disk Tip Deflection. 16th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC 2016), Apr 2016, Honolulu, United States. hal-01516498

HAL Id: hal-01516498

https://hal.science/hal-01516498

Submitted on 1 May 2017

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear Cyclic Transient Dynamic Analysis for Bladed Disk Tip Deflection

Tong Hui<fnref n="1" />*, Letian Wang<fnref n="1" />**, Suryarghya Chakrabarti<fnref n="1" />, Siddharth Ashar<fnref n="1" />, Steven Szczap<fnref n="2" />, Anthony Laude<fnref n="2" />

<figure><img src="image_2.png" /></figure>

ISROMAC 2016

International Symposium on Transport Phenomena and Dynamics of Rotating Machinery

Hawaii, Honolulu

April 10-15, 2016

## Abstract

High aspect ratio rotating structures such as fan blade and first stage compressor blades have the potential to undergo large dynamic deflections under unsteady aerodynamic excitations. The presence of nonlinear geometric effects in the vibratory response limits the use of traditional cyclic symmetry in linear harmonic forced response analyses. In this paper a methodology is developed to calculate the nonlinear dynamic response of cyclic structures undergoing large vibratory deformation under a traveling wave excitation using a double-sector model and transient dynamic analysis. It is shown that the results from the double-sector model are in good agreement with the results obtained from a full wheel nonlinear transient dynamic analysis when nonlinear geometric effects are concentrated near the blade tip away from the sector edges where cyclic constraint equations are enforced. The quality of agreement gradually degrades as the amount of nonlinear geometric effects near the cyclic sector boundaries increase.

## Keywords

cyclic symmetry; nonlinear dynamics; large deformation; traveling wave.

## INTRODUCTION

Accurate prediction of blade tip deflection is critical to prevent rub and wear while maintaining high performance for a jet engine fan, compressor and turbine blades. Due to their rotational periodicity, the theory of cyclic symmetry has been widely used to investigate the vibratory characteristics of bladed disks through harmonic forced response analyses.

Comprehensive treatment of the application of linear vibration theory to rotationally periodic structures exists in the literature [1, 2, 3]. While linear vibration theory is used extensively in the industry during the design process through frequency avoidance (Campbell diagrams, Interference diagrams etc.) and linear forced response calculations, there are certain situations where modeling nonlinear phenomena in the structure become critical.

One such situation arises when there is a need to model friction damping between contacting interfaces in the structure (turbine blade shrouds [4], platform dampers [5], combustor seals [6] etc.). A time stepping analysis on the full finite element model is computationally prohibitive in most practical applications. Therefore the literature on modeling rotationally periodic structure with localized friction nonlinearity is focused on efficient model reduction and solution techniques for such systems. Dynamic substructuring methods based on component mode synthesis [7, 8, 9] have been widely used for reducing large finite element models to a small number of nonlinear (master) degrees of freedom (which experience friction) and a small number of modal coordinates which approximate the behavior of the (slave) degrees of freedom. The

computational gains of model order reduction are complemented by replacing the time stepping solution with more efficient frequency domain techniques such as the harmonic balance method [10] or the hybrid frequency-time domain method [11, 12] which work in the frequency domain assuming time periodicity of the nonlinear response.

While the methods discussed above work well for localized nonlinearities, they cannot be applied to systems with material or geometrical nonlinearities [13, 14]. In such cases solution is generally obtained through a time-stepping approach which typically requires factorization of the stiffness matrix at regular intervals during the time integration process. Geometrically nonlinear phenomena may occur in the fan and front stage compressor blades of a jet engine during stall and surge events due to large dynamic deflections. Conventional linear analyses based on linear mode superposition theory can incur significant errors in tip clearance calculations during these events. Theoretically, the presence of such nonlinearities breaks down the block-circulant nature of the full wheel stiffness matrix such that the response of the full wheel structure cannot be accurately captured by a single geometric sector any more. This would mean that a nonlinear transient time domain analysis on the full wheel structure is required to accurately solve the nonlinear forced response problem. The computational cost of such an analysis makes it prohibitive for use in any practical application.

This paper presents an approximate method for calculating max deflection values from a single-sector model using the commercial finite element software ANSYS. The prestress

<footnote n="1">¹GE Global Research Center, Niskayuna 12309, United States</footnote>

<footnote n="2">²GE Aviation, West Chester, OH 45069, United States</footnote>

<footnote>*Corresponding author: tong.hui@ge.com</footnote>

<footnote>**Corresponding author: wangl@ge.com</footnote>
and statically deformed shape is computed using circumferen-
tial harmonic index 0 coupling on the sector boundaries and
performing a static analysis using the appropriate thermal and
structural loads. Thereafter, the sector coupling is replaced
with a non-zero harmonic index coupling according to the
excitation engine order and a transient time dynamic analysis
is performed. The results are compared to a similar analysis
performed on the full wheel structure.

1. THEORETICAL FORMULATION FOR STUC-
TURES WITH CYCLIC SYMMETRY

Consider a structure with geometrically cyclic periodicity
as shown in Figure 1. The structure has nS geometrically
identical sectors numbered from 0 to nS – 1. The positive
rotation direction (θ) is defined as increasing from sector #0
to nS – 1. The cyclic sector #0 is referred to as the reference
cyclic sector. The cyclic edges are the interfaces between two
adjacent sectors. By definition, the high cyclic edge of the
reference cyclic sector if obtained by a positive rotation of
2π/nS radians of the low cyclic edge. With the definitions
above, the degrees of freedom (DOFs) of the reference cyclic
sector can be partitioned into nI internal DOFs, nL low cyclic
edge DOFs, and nH high cyclic edge DOFs. The total num-
ber of DOFs nT of the structure is therefore calculated as:
nT = nS(nL + nI). Throughout this manuscript the following
conventions are followed:

• sector numbers are written within {} in the superscript

• spatial harmonics are written within [] in the superscript

* low edge, internal and high edge quantities are written
with subscript L, I and H respectively.

• complex quantities are underlined

• full wheel quantities are written with a tilde

<figure><img src="image_3.png" /><figcaption>Figure 1. A rotationally periodic structure.</figcaption></figure>

1.1 Full wheel matrices

Let $\mathbf{x}_{\text{sector}}^{\{0\}}$ be the vector made of the low (L), internal (I), and high (H) DOFs associated with the reference cyclic sector, and let us assume it can be written as:

$$
\mathbf{x}_{\text{sector}}^{\{0\}} = \begin{bmatrix} \mathbf{x}_{\text{I}}^{\{0\}} \\ \mathbf{x}_{\text{L}}^{\{0\}} \\ \mathbf{x}_{\text{H}}^{\{0\}} \end{bmatrix}, \qquad (1)
$$

noting that $\mathbf{x}_H^{(0)}$ is actually $\mathbf{x}_L^{(1)}$. The stiffness (and mass)
matrix of the reference sector can therefore be written as:

$$
\mathbf{K}_{\text{sector}}^{\{0\}} = \begin{bmatrix} \mathbf{K}_{\text{LL}} & \mathbf{K}_{\text{LI}} & \mathbf{K}_{\text{LH}} \\ \mathbf{K}_{\text{LI}}^{\text{T}} & \mathbf{K}_{\text{II}} & \mathbf{K}_{\text{IH}} \\ \mathbf{K}_{\text{LH}}^{\text{T}} & \mathbf{K}_{\text{IH}}^{\text{T}} & \mathbf{K}_{\text{HH}} \end{bmatrix}. \quad (2)
$$

In the interest of being concise, the stiffness matrix is used for
illustrating the theory in this section. The mass matrix follows
identical transformation laws. Let us assume that the DOFs x̃
of the structure are ordered so that:

$$
\tilde{\mathbf{x}} = \begin{bmatrix} \mathbf{x}^{\{0\}} \\ \mathbf{x}^{\{1\}} \\ \vdots \\ \mathbf{x}^{\{n_S - 1\}} \end{bmatrix}, \qquad (3)
$$

where the superscripts indicate the cyclic sector numbers, and
where $\mathbf{x}^{[k]}$ contains the low cyclic edge and internal DOFs of
cyclic sector $k$, as shown below:

$$
\mathbf{x}^{[k]} = \begin{bmatrix} \mathbf{x}_{\mathrm{L}}^{[k]} \\ \mathbf{x}_{\mathrm{I}}^{[k]} \end{bmatrix}. \qquad (4)
$$

Since the structure is cyclic, the stiffness (and mass) matrix of
each cyclic sector must be identical to each other, that is:

$$
\mathbf{K}_{\text{sector}}^{\{0\}} = \mathbf{K}_{\text{sector}}^{\{1\}} = \dots = \mathbf{K}_{\text{sector}}^{\{n_S - 1\}}. \qquad (5)
$$

Taking into account their inherent symmetry properties, the
assembled stiffness (and mass) matrices of the full structure
can therefore expressed as:

$$
\tilde{\mathbf{K}} = \begin{bmatrix} \mathbf{K}_A & \mathbf{K}_B & & & & & \mathbf{K}_B^T \\ \mathbf{K}_B^T & \mathbf{K}_A & \mathbf{K}_B & & & & \\ & \mathbf{K}_B^T & \mathbf{K}_A & \mathbf{K}_B & & & \\ & & \mathbf{K}_B^T & \mathbf{K}_A & \mathbf{K}_B & & \\ & & & \ddots & \ddots & \ddots & \\ & & & & \mathbf{K}_B^T & \mathbf{K}_A & \mathbf{K}_B \\ & \mathbf{K}_B & & & \mathbf{K}_B^T & \mathbf{K}_A & \mathbf{K}_B \end{bmatrix}, \quad (6)
$$

where

$$
\mathbf{K}_A = \begin{bmatrix} \mathbf{K}_{LL} + \mathbf{K}_{HH} & \mathbf{K}_{LI} \\ \mathbf{K}_{LI}^T & \mathbf{K}_{II} \end{bmatrix} \quad \text{and} \quad (7)
$$

$$
\mathbf{K}_B = \begin{bmatrix} \mathbf{K}_{\mathrm{LH}} & \mathbf{0} \\ \mathbf{K}_{\mathrm{IH}} & \mathbf{0} \end{bmatrix}. \qquad (8)
$$
## 1.2 Complex single sector formulation with the Extended Fourier Matrix

The stiffness (and mass) matrix of the full cyclic structure is block circulant and therefore, can be block-diagonalized by the extended Fourier Matrix $\underline{E}$ of size $[n \times (n_L + n_I)]$ defined as:

$$ \underline{E} = \underline{E}_{[n]} \otimes \underline{I}_{[n_L+n_I]}, \quad (9) $$

where $\underline{E}_{[n]}$ is the n-by-n Fourier matrix, and $\underline{I}_{[n_L+n_I]}$ is the $(n_L + n_I)$-by-$(n_L + n_I)$ identity matrix.

$$ \underline{E}^* \tilde{\underline{K}} \underline{E} = \begin{bmatrix} \underline{K}^{[0]} \\ \underline{K}^{[1]} \\ \vdots \\ \underline{K}^{[n_S-1]} \end{bmatrix} \quad (10) $$

Matrix $\underline{K}^{[p]}$ is the stiffness matrix corresponding to circumferential harmonic index $p$. It is complex and Hermitian of the form

$$ \underline{K}^{[p]} = \begin{bmatrix} \underline{K}_{LL} + \underline{K}_{HH} + \underline{K}_{LH} e^{ip\theta} + \underline{K}_{LH}^{\mathrm{T}} e^{-ip\theta} & \underline{K}_{LI} + \underline{K}_{II} e^{-ip\theta} \\ \underline{K}_{LI}^{\mathrm{T}} + \underline{K}_{IH} e^{ip\theta} & \underline{K}_{II} \end{bmatrix}. \quad (11) $$

The mass matrix $\underline{M}^{[p]}$ follows an identical structure. An arbitrary full wheel forcing is transformed by the Extended Fourier matrix into its different spatial harmonic components

$$ \underline{E}^* \tilde{\underline{f}} = \begin{bmatrix} \underline{F}^{[0]} \\ \underline{F}^{[1]} \\ \vdots \\ \underline{F}^{[n_S-1]} \end{bmatrix} \quad (12) $$

If the forcing corresponds to a single spatial harmonic (say p), only the block corresponding to the $p^{th}$ circumferential harmonic ($\underline{F}^{[p]}$) is non-zero. Similarly the full wheel state vector can be transformed using the Extended Fourier matrix as

$$ \underline{E}\tilde{\underline{x}} = \tilde{\underline{U}} = \begin{bmatrix} \underline{U}^{[0]} \\ \underline{U}^{[1]} \\ \vdots \\ \underline{U}^{[n_S-1]} \end{bmatrix} \quad (13) $$

### 1.2.1 Equations of motion

The forced response problem for the full wheel structure (assuming beta damping) can be written in the time domain as

$$ \tilde{\underline{M}}\tilde{\underline{x}} + \beta \tilde{\underline{K}}\tilde{\underline{x}} + \tilde{\underline{K}}\tilde{\underline{x}} = \tilde{\underline{f}}. \quad (14) $$

For a forcing consisting of a single spatial harmonic (say p), transforming the state vector using (13), pre-multiplying by $\underline{E}^*$, and applying (10) and (12) the system of equations reduce to

$$ \underline{M}^{[p]} \tilde{\underline{U}}^{[p]} + \beta \underline{K}^{[p]} \tilde{\underline{U}}^{[p]} + \underline{K}^{[p]} \underline{U}^{[p]} = \underline{F}^{[p]}. \quad (15) $$

## 1.3 Double-sector formulation

Matrix $\underline{K}^{[p]}$ (and $\underline{M}^{[p]}$) can also be represented as a real valued matrix by using the double-sector approach. In this approach the cyclic symmetry problem is solved by using two sectors - one sector (the basic or real sector) stores the real part of the solution while the second sector (duplicate or imaginary sector) stores the imaginary part of the solution (Figure 2). This is true for real forcing in a time domain analysis. In a frequency domain analysis, typically the forcing vector is complex (representing relative time phasing of forces at different locations on the structure.) In this case both sector dofs will be complex. However the double-sector model separates out spatial phasing from time phasing. The basic and duplicate sectors store the spatial phasing due to cyclic symmetry whereas the individual sector dofs store the time phasing information. In other words the double-sector model reduces a cyclic problem to an equivalent non-cyclic system with real system matrices. In the double-sector formulation,

<figure><img src="image_4.png" /><figcaption>Figure 2. Double-sector model.</figcaption></figure>

the full wheel response is expressed in terms of a Fourier series of two quantities $\underline{u}^c$ and $\underline{u}^s$ such that the $k^{th}$ sector quantities can be expressed as

$$ \underline{x}_{\text{sector}}^{[k]} = \sum_{h=0}^{\infty} \left[ \underline{u}^{c,[h]} \cos(kh\theta) - \underline{u}^{s,[h]} \sin(kh\theta) \right], \quad (16) $$

Here $\theta$ is the geometric sector angle $\frac{2\pi}{n_S}$ and $h$ represents the spatial harmonic. Assuming the full wheel response is governed by a single spatial harmonic (say p), the transformed state vectors $\underline{u}^{c,[p]}$ and $\underline{u}^{s,[p]}$ are ordered as follows

$$ \underline{u}^{c,[p]} = \begin{bmatrix} \underline{u}_{L}^{c,[p]} \\ \underline{u}_{L}^{c,[p]} \\ \underline{u}_{H}^{c,[p]} \end{bmatrix} \quad \text{and} \quad \underline{u}^{s,[p]} = \begin{bmatrix} \underline{u}_{L}^{s,[p]} \\ \underline{u}_{L}^{s,[p]} \\ \underline{u}_{H}^{s,[p]} \end{bmatrix}, \quad (17) $$

where the subscripts L, I and H represent the low, internal, and high sector edge quantities respectively. The low and high sector edge DOFs are related by the $p^{th}$ circumferential harmonic constraint equations as follows

$$ \begin{bmatrix} \underline{u}_{H}^{c,[p]} \\ \underline{u}_{H}^{s,[p]} \end{bmatrix} = \begin{bmatrix} \cos(p\theta) & -\sin(p\theta) \\ \sin(p\theta) & \cos(p\theta) \end{bmatrix} \begin{bmatrix} \underline{u}_{L}^{c,[p]} \\ \underline{u}_{L}^{s,[p]} \end{bmatrix}. \quad (18) $$
The corresponding real-valued constraint-eliminated double-sector stiffness matrix can be expressed as

$$ \mathbf{K}_{\mathrm{DS}}^{[p]} = \begin{bmatrix} \mathbf{K}_{\mathrm{B}}^{[p]} & \mathbf{K}_{\mathrm{BD}}^{[p]} \\ \mathbf{K}_{\mathrm{BD}}^{[p]^{\mathrm{T}}} & \mathbf{K}_{\mathrm{B}}^{[p]} \end{bmatrix}, \qquad (19) $$

where,

$$ \mathbf{K}_{\mathrm{B}}^{[p]} = \begin{bmatrix} \mathbf{K}_{\mathrm{LL}} + \mathbf{K}_{\mathrm{LH}}\cos(p\theta) + \mathbf{K}_{\mathrm{LH}}^{\mathrm{T}}\cos(p\theta) + \mathbf{K}_{\mathrm{HH}} & \mathbf{K}_{\mathrm{LI}} + \mathbf{K}_{\mathrm{IH}}\cos(p\theta) \\ \mathbf{K}_{\mathrm{LI}}^{\mathrm{T}} + \mathbf{K}_{\mathrm{IH}}\cos(p\theta) & \mathbf{K}_{\mathrm{II}} \end{bmatrix}, \qquad (20) $$

$$ \mathbf{K}_{\mathrm{BD}}^{[p]} = \begin{bmatrix} -\mathbf{K}_{\mathrm{LH}} \sin(p\theta) + \mathbf{K}_{\mathrm{LH}}^{\mathrm{T}} \sin(p\theta) & \mathbf{K}_{\mathrm{IH}}^{\mathrm{T}} \sin(p\theta) \\ -\mathbf{K}_{\mathrm{IH}} \sin(p\theta) & \mathbf{0} \end{bmatrix}. \qquad (21) $$

It can be seen that the double-sector stiffness matrix is a combination of real and imaginary parts of the complex valued stiffness matrix (11). The double-sector mass matrix can be derived similarly.

Rotating turbomachinery components are typically exposed to travelling wave aerodynamic excitation. This means that each blade experiences an identical dynamic excitation but separated by a predetermined time delay [15]. Similar to deflections (16), the time domain forcing on the k^th sector is expressed as a Fourier series

$$ \mathbf{f}^{[k]} = \sum_{h=0}^{\infty} \left[ \mathbf{f}^{c, [h]} \cos(kh\theta) - \mathbf{f}^{s, [h]} \sin(kh\theta) \right]. \qquad (22) $$

Considering the forcing to contain a single spatial harmonic (say p), the basic and duplicate sector forcing $\mathbf{f}^{c, [p]}$ and $\mathbf{f}^{s, [p]}$ can be related as follows.

Substituting $k = 0$ in Eq. (22) yields $\mathbf{f}^{c, [h]} = \mathbf{f}^{[0]}$. (i.e the basic sector forcing is the same as the reference sector forcing). Next consider a harmonic form of the reference sector (equal to the basic sector) forcing

$$ \mathbf{f}^{[0]}(t) = \mathbf{f}^{c, [h]}(t) = \mathbf{f}_r \cos(\omega t) - \mathbf{f}_i \sin(\omega t). \qquad (23) $$

With this form the first sector forcing (sector following the reference sector) can be written as

$$ \mathbf{f}^{[1]}(t) = \mathbf{f}_r \cos(\omega t + p\theta) - \mathbf{f}_i \sin(\omega t + p\theta). \qquad (24) $$

With some trigonometrical manipulations and comparing the resulting form to (22), it can be shown that the duplicate sector time domain forcing can be written as

$$ \mathbf{f}^{s, [h]}(t) = -\mathbf{f}_i \cos(\omega t) - \mathbf{f}_r \sin(\omega t). \qquad (25) $$

Thus the equations of motion for the double-sector model with single spatial harmonic (p) excitation become:

$$ \mathbf{M}_{\mathrm{DS}}^{[p]} \ddot{\mathbf{u}}_{\mathrm{DS}}^{[p]} + \beta \mathbf{K}_{\mathrm{DS}}^{[p]} \dot{\mathbf{u}}_{\mathrm{DS}}^{[p]} + \mathbf{K}_{\mathrm{DS}}^{[p]} \mathbf{u}_{\mathrm{DS}}^{[p]} = \mathbf{f}_{\mathrm{DS}}^{[p]}. \qquad (26) $$

where

$$ \mathbf{u}_{\mathrm{DS}}^{[p]} = \begin{bmatrix} \mathbf{u}^{c, [p]} \\ \mathbf{u}^{s, [p]} \end{bmatrix}, \qquad (27) $$

$$ \mathbf{f}_{\mathrm{DS}}^{[p]} = \begin{bmatrix} \mathbf{f}_r \cos(\omega t) - \mathbf{f}_i \sin(\omega t) \\ -\mathbf{f}_i \cos(\omega t) - \mathbf{f}_r \sin(\omega t) \end{bmatrix} \qquad (28) $$

This system of equations is used to solve the dynamics of the double-sector model, and for linear systems, the full wheel deflections are obtained by using (16).

## 2. LARGE-DEFORMATION EFFECTS

The large deformation theory can be explained by combining the effects of rotation and deformation [16]. Consider a body undergoing deformation from an initial (reference) condition to a final (current) condition. In a deformation process described by a displacement field $\mathbf{u}(\mathbf{X})$, the point P moves to position p with coordinate x and in the spatial coordinate frame where $\mathbf{x} = \mathbf{X} + \mathbf{u}(\mathbf{X})$ and the point Q moves to position q whose coordinate $\mathbf{x} + d\mathbf{x}$ can be written as

$$ \mathbf{x} + d\mathbf{x} = \mathbf{X} + d\mathbf{X} + \mathbf{u}(\mathbf{X} + d\mathbf{X}). \qquad (29) $$

The deformation gradient is defined as:

$$ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \mathbf{I} + \frac{\partial \mathbf{u}}{\partial \mathbf{X}}. \qquad (30) $$

Commercial solvers such as ANSYS break the deformation gradient into a rotation component and a shape change component using the polar decomposition theorem

$$ \mathbf{F} = \mathbf{R}\mathbf{U}, \qquad (31) $$

where $\mathbf{R}$ is a proper orthogonal tensor, i.e., $\mathbf{R}^{-1} = \mathbf{R}^T$ and det($\mathbf{R}$) = 1, representing the average rotation of the material at a point; $\mathbf{U}$ denotes the right stretch tensor which is a symmetric tensor, i.e., $\mathbf{U} = \mathbf{U}^T$

Implementation is done following an incremental solution method [16]. The net effect is that the stiffness matrix of the system becomes state dependent. As a result, there is sector to sector variation in the stiffness matrix and the transformation (10) does not create a block diagonal matrix any longer which prevents us from decoupling the equations of motion for each spatial harmonic.

Large dynamic deflections may typically occur in high aspect ratio blades such as the fan or the front stage compressor blades. The large deformation effect is predominantly observed in the regions of highest displacements (typically blade tips). However, in the portion of the structure where there is sector to sector mating (typically the disk), the deformations are much smaller and hence large deformation effects are also small. In such a situation the double-sector cyclic constraint equations can still be used to describe the time domain behavior of a single sector. However, expansion of the displacements to the full wheel following (16) will be incorrect for the dofs which do undergo large deformation. For a traveling wave excitation each sector undergoes the identical response with phase shifted in time. Thus, for the case of a traveling wave excitation it is sufficient to track the response of a single sector using cyclic coupling on the sector edges through constraint equations with the assumption that dofs belonging to the sector edges encounter little geometric nonlinearity.
<figure><img src="image_5.png" /><figcaption>Figure 3. Deformation of a body from a reference condition to current condition [17].</figcaption></figure>

# 3. FINITE ELEMENT IMPLEMENTATION

## 3.1 Time integration

In a large deformation problem, the force-deflection relationship is nonlinear and cannot be expressed as a constant stiffness matrix. Instead the equations of motion are expressed as

$$M\ddot{u}(t) + C\dot{u}(t) + f^i(u(t)) = f^e(t), \quad (32)$$

where $f^i$ represents the state-dependent internal restoring force vector and $f^e(t)$ is the time dependent external forcing. For nonlinear large deformation problems, commercial FEA softwares (such as ANSYS) typically follow an implicit and incremental solution scheme [16]. First the equations of motion are expressed at time $t_{n+1}$

$$M\ddot{u}_{n+1}^{(k)} + C\dot{u}_{n+1}^{(k)} + f^i(u_{n+1}^{(k)}) = f_{n+1}^e, \quad (33)$$

where the superscript (k) represents the Newton-Raphson iteration number and

$$u_{n+1}^{(k)} = u_{n+1}^{(k-1)} + \Delta u^{(k)}. \quad (34)$$

Following the Newmark time integration scheme, $u_{n+1}^{(k)}$ and $\dot{u}_{n+1}^{(k)}$ are written as

$$\dot{u}_{n+1}^{(k)} = \dot{u}_n + [(1 - \delta)\ddot{u}_n + \delta\ddot{u}_{n+1}^{(k)}]\Delta t \quad (35)$$

$$u_{n+1}^{(k)} = u_n + \dot{u}_n \Delta t + \left[ \left( \frac{1}{2} - \alpha \right) \ddot{u}_n + \alpha \ddot{u}_{n+1}^{(k)} \right] \Delta t^2, \quad (36)$$

where $\alpha$ and $\delta$ are scalar parameters for Newmark's time integration scheme. Introducing the residual vector $R_{n+1}$ as

$$R_{n+1}(u_{n+1}^{(k)}) = f_{n+1}^e - M\ddot{u}_{n+1}^{(k)} - C\dot{u}_{n+1}^{(k)} - f^i(u_{n+1}^{(k)}), \quad (37)$$

the solution to the nonlinear equation $R_{n+1}(u_{n+1}^{(k)}) = 0$ is found using the Newton-Raphson method, which yields

$$R_{n+1}(u_{n+1}^{(k)}) + \frac{\partial R_{n+1}(u_{n+1}^{(k)})}{\partial u_{n+1}^{(k)}} \Delta u_{n+1}^{(k)} = 0. \quad (38)$$

Combination of (33) to (38) yields

$$\left[ \frac{1}{\alpha \Delta t^2} M + \frac{\delta}{\alpha \Delta t} C + K_{n+1}^t (u_{n+1}^{(k)}) \right] \Delta u_{n+1}^{(k)} = R_{n+1} (u_{n+1}^{(k)}) \quad (39)$$

where

$$K_{n+1}^t (u_{n+1}^{(k)}) = \frac{\partial f_{n+1}^i (u_{n+1}^{(k)})}{\partial u_{n+1}^{(k)}}, \quad (40)$$

is the tangent stiffness matrix. During time integration, (39) is solved iteratively at each time step till the $R_{n+1}(u_{n+1}^{(k)})$ and $\Delta u_{n+1}^{(k)}$ are below specified tolerances.

## 3.2 Inclusion of centrifugal effects

Rotating turbomachinery components are subject to a combination of static loads (centrifugal load, steady pressure loads etc.) with zero spatial harmonic and vibratory loads with non-zero spatial harmonic (circumferential flow distortions, nozzle wakes, structurally transmitted synchronous drivers). The static loads induce static structural deformations with zero spatial harmonic while the vibratory loads induce non-zero spatial harmonic deformations. Therefore, the total deformation of a cyclic structure is a combination of zero and non-zero spatial harmonics. Moreover, the static stresses on the structure cause stress stiffening which affects the dynamic behavior of the structure.

This poses a difficulty in direct application of the double-sector model as the double-sector system matrices are derived for a specific spatial harmonic. For linear systems this is tackled using a two-step process where first a static analysis is performed using zero spatial harmonic matrices and forcing. Then the prestressed stiffness matrices are used to calculate the vibratory response of the structure using a (mode superposition or full) harmonic analysis. This process needs to be modified for a structure where dynamic deflections can be large enough to require modeling of geometrically nonlinear effects. The modified process is as follows:

1. Static analysis for the bladed disk subject to steady (spatial harmonic 0) loads: Constraint equations are applied on the cyclic sector edges using Eq. (18) with $p=0$ and the resulting system of equations are solved. Note that substituting $p=0$ in (20) and (21), the off-diagonal blocks in the stiffness matrix are zero. This yields a solution $u^{c,[0]} = u^{s,[0]}$ and for the sector edges $u_H^{c,[0]} = u_H^{s,[0]} = u_L^{c,[0]} = u_L^{s,[0]} = u_0$. At the end of the step, the nodal reaction forces corresponding to all nodes of the structure are stored along with the nodal displacements for the high and low edge nodes.
<figure><img src="image_6.png" /><figcaption>Figure 4. Flow chart of static analysis.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Flow chart of dynamic analysis.</figcaption></figure>

2. Dynamic analysis for the bladed disk subject to steady loading and dynamic excitation: The reaction force at each node output from the previous step is first applied on each node, and then the constraint equations are modified for the spatial harmonic corresponding to the excitation (say p) as follows:

$$ \begin{bmatrix} \mathbf{u}_H^{c,[p]} - \mathbf{u}_0 \\ \mathbf{u}_H^{s,[p]} - \mathbf{u}_0 \end{bmatrix} = \begin{bmatrix} \cos(p\theta) & -\sin(p\theta) \\ \sin(p\theta) & \cos(p\theta) \end{bmatrix} \begin{bmatrix} \mathbf{u}_L^{c,[p]} - \mathbf{u}_0 \\ \mathbf{u}_L^{s,[p]} - \mathbf{u}_0 \end{bmatrix}. \quad (41) $$

The ANSYS based implementation for static and dynamic analyses is summarized in Fig. 4 and Fig. 5 respectively. The purpose of static analysis is to obtain the reaction force for spatial harmonic zero condition. In the first step of the transient dynamic analysis, a quasi-static analysis is performed (by switching off dynamic effects) after applying the reaction forces obtained from static analysis. Then harmonic zero constraint equations are then switched with the modified constraint equations for the excitation spatial harmonic according to Eq. (41). A nonzero spatial harmonic traveling wave excitation is then applied, dynamic effects are switched on and conventional time integration is performed following a scheme outlined in section 3.1 till steady state behavior is observed. It should be noted that, gyroscopic effects which can be important in rotor dynamics are negligible in the present problem since axis of the structure is assumed to be perfectly fixed.

# 4. NUMERICAL EXAMPLES

The double-sector model based approach is further examined and applied on a bladed-disk model as shown in Fig. 6. The full wheel consists of total 30 sectors each sector containing one blade. Eight-noded brick elements are used for meshing. In the analysis, the inner rim of the disk is fully constrained.

## 4.1 Modal analysis

The cyclic boundary condition is first examined through the modal analysis. The frequencies for harmonic index: 3, 5, 10 and 15 are calculated using the sector and full wheel models respectively and compared in Table 1. It is shown that both the models give identical frequencies for all the modes at each nodal diameter. The mode shapes in the full wheel model are presented in Fig. 7. The number of nodal diameters is clearly shown in the full wheel mode shapes. The same mode shapes can be obtained through mode shape expansion of the sector model in ANSYS.

## 4.2 Dynamic analysis

In the dynamic analysis, the traveling wave excitation together with the steady rpm loading is applied on blades of the full wheel and sector models. In the full wheel model, the traveling wave excitation is applied in terms of pressure as shown in Fig. 8, with the appropriate inter-blade phase angle. The double-sector model is excited by the forcing as shown in Eq. (28). In the simulation, we consider 3000 rpm spinning and spatial harmonic $p = 3$ excitation. The system is excited at its first natural frequency (227.90Hz). In order to prevent
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Full wheel mode shapes for $ND = 3$ and $ND = 5$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** FEA full wheel and sector models (dimensions in inches).</figcaption></figure>

unbounded vibration, structural damping with the damping ratio, $10^{-3}$, is applied. The displacements in axial ($u_a$) and radial direction ($u_r$) are calculated at the same node (green dot in Fig. 8) in the sector and full wheel models for comparison. The axial and radial displacements are plotted in Figs 9 and 10.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Traveling wave excitation on the full wheel and sector models.</figcaption></figure>

It is shown that the sector model calculated displacements

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Full wheel vs. sector for $u_a$ with $\times 1$ scaled force.</figcaption></figure>

are in good agreement with the full wheel model. Due to the structural damping, the vibration is controlled after a few cycles. For this forcing level the nonlinear geometric effects on the blade tips are noticeable but on the disk sector edge nodes are negligible. Thus the sector model approach is able to predict the nonlinear blade tip deflections accurately. Next, the dynamic excitation is scaled up $\times 10$ such that large deformation effects on the disk cyclic edges are not negligible (Figure 11). Thus the differences between the sector model and full wheel blade tip deflection traces are more evident (Figs. 12 and 13). However, the results are still much more accurate than performing a linear analysis. The sector model underestimates the $u_r$ by around 10% when compared to the full wheel model, whereas the error between the linear and nonlinear full wheel analysis can be as large as 400% (Figure 14). The nonlinearity in the response is evident from the frequency doubling observed in Figure 13. The double-sector model is able to reproduce similar frequency doubling behavior.

The nonlinear effects on blade deflection are investigated further. Specifically, the radial displacement at the tip of blades is of the greatest interest, since it is a crucial factor for clearance between blade and casing. Insufficient clearance
<table><caption>**Table 1.** Frequency comparison of sector and full wheel models.</caption><thead><tr><th rowspan="2">mode number</th><th colspan="2">ND=3</th><th colspan="2">ND=5</th><th colspan="2">ND=10</th><th colspan="2">ND=15</th></tr><tr><th>sector</th><th>full wheel</th><th>sector</th><th>full wheel</th><th>sector</th><th>full wheel</th><th>sector</th><th>full wheel</th></tr></thead><tbody><tr><td>1</td><td>199</td><td>199</td><td>204</td><td>204</td><td>225</td><td>225</td><td>236</td><td>236</td></tr><tr><td>2</td><td>912</td><td>912</td><td>939</td><td>939</td><td>1076</td><td>1076</td><td>1074</td><td>1074</td></tr><tr><td>3</td><td>1073</td><td>1073</td><td>1076</td><td>1076</td><td>1087</td><td>1087</td><td>1255</td><td>1255</td></tr><tr><td>4</td><td>1923</td><td>1923</td><td>1885</td><td>1885</td><td>1730</td><td>1730</td><td>1576</td><td>1576</td></tr><tr><td>5</td><td>2276</td><td>2276</td><td>2313</td><td>2313</td><td>2483</td><td>2483</td><td>2815</td><td>2815</td></tr><tr><td>6</td><td>4752</td><td>4752</td><td>4501</td><td>4501</td><td>3563</td><td>3563</td><td>2882</td><td>2882</td></tr><tr><td>7</td><td>4987</td><td>4987</td><td>5215</td><td>5215</td><td>5408</td><td>5408</td><td>5445</td><td>5445</td></tr><tr><td>8</td><td>5514</td><td>5514</td><td>5439</td><td>5439</td><td>5470</td><td>5470</td><td>5483</td><td>5483</td></tr><tr><td>9</td><td>5550</td><td>5550</td><td>5504</td><td>5504</td><td>5498</td><td>5498</td><td>5529</td><td>5529</td></tr><tr><td>10</td><td>6954</td><td>6954</td><td>6747</td><td>6747</td><td>6482</td><td>6482</td><td>6390</td><td>6390</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Full wheel vs. sector for $u_r$ with ×1 scaled force.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of radial displacements on the sector cyclic edge obtained from linear and nonlinear analysis. Percentages show difference between linear and nonlinear analysis.</figcaption></figure>

may result in severe rubbing of blades, whereas too much clearance will compromise stage efficiency. The radial deflection of one point at the blade tip is tracked at two positions A and B (Figure 14).

Figure 14 summarizes the linear and nonlinear analysis deflections of the observation point at position A and B versus force scaling factor. The deflections in plot are normalized by the force scaling factor. It is apparent that linear and nonlinear calculations suggest completely different radial deflections especially when force scaling factor increases beyond 5. The following conclusions are obtained: first, linear calculation predicts the normalized radial deflection with little variation when scaling factor increases beyond 5. However, the nonlinear calculation indicates nonlinear dependency of the radial deflection on the scaling factor within the whole range. Second, the nonlinear calculation of radial deflection shows significant antisymmetric displacements for position A and B, whereas linear estimation wrongly suggests approximately symmetric displacements for the deflection at A (-0.0045) and B (0.0045). These findings support the necessity of using nonlinear deformation based calculation in predicting radial deflection of bladed disks.

Although the model used for this study was a heavily simplified planar geometry, it justified the use of a double-sector model to get fairly accurate nonlinear response prediction as

long as cyclic edge dofs do not experience significant nonlinear effects. Aircraft engine components which would experience large vibration phenomena, such as first few compressor stages, are usually free standing blades with the cyclic edge confined to the disk. In these components the blades are thin and long and do undergo large vibrational motion but the disk is bulkier which restricts the cyclic edge deformation to be close to the linear regime. The method presented here is well-suited to model such cases.

## 5. CONCLUSIONS

This paper presented a double-sector model based solution for simulating the nonlinear dynamic response of bladed disks subject to large vibrating deformation. The finite element implementation of the double-sector model is done in the commercial software ANSYS. The proposed method is capable of including both the static spatial harmonic zero deformations (from steady loads such as wheel spin) and nonzero spatial harmonic deformations due to unsteady aerodynamic excitation. The methodology works under the assumption that large deformation effects are considerable only near the blade tips and are small near the sector cyclic edges.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Full wheel vs. sector for $u_a$ with force ×10 scaled force.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Full wheel vs. sector for $u_r$ with ×10 scaled force.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Radial displacements at position A and B for different force scaling factors.</figcaption></figure>

The performance of the proposed double-sector model is verified against the full wheel model through modal analysis and nonlinear transient dynamic analysis. In the modal analysis, frequencies for the first 10 modes with different nodal diameters are examined. Both models gave identical results. In the transient dynamic analysis, radial and axial displacements are calculated using the double-sector model and compared with the full wheel model results. It is found that there is excellent agreement between the sector and full wheel models when nonlinear deformation effects on the cyclic sector edges are small. The degree of agreement gradually degrades with the increase in dynamic excitation as large deflection effects on the cyclic sector edges become more prominent. Tip deflection of the blade along with increasing excitation force was investigated further. The effects on large deformation are discussed through the asymmetric vibratory response observed at the blade tips. It is shown that when dynamic deflections during the vibration cycle are large, it is essential to model nonlinear large deformation effects as it can have a significant impact on the radial tip deflection and hence clearance calculation.

This proposed double-sector based approach has been applied to GE Aviation compressor dynamic tip clearance analysis.

## ACKNOWLEDGMENTS

The authors are grateful to GE Aviation for providing the financial support for this project, and to GE Global Research Center and GE Aviation for giving permission to publish this work.

## REFERENCES

[1] D. L. Thomas. Dynamics of rotationally periodic structures. *International Journal for Numerical Methods in Engineering*, 14(1):81–102, 1979.

[2] F. W. Williams. An algorithm for exact eigenvalue calculations for rotationally periodic structures. *International journal for numerical methods in engineering*, 23(4):609–622, 1986.

[3] J. Wildheim. Vibrations of rotating circumferentially periodic structures. *The Quarterly Journal of Mechanics and Applied Mathematics*, 34(2):213–229.

[4] B-D Yang and C-H Menq. Modeling of friction contact and its application to the design of shroud contact.
*Journal of Engineering for Gas Turbines and Power,*
119(4):958-963, 1997.

[5] BD Yang and CH Menq. Characterization of contact kinematics and application to the design of wedge dampers in turbomachinery blading: part 1 stick-slip contact kinematics. *Journal of engineering for gas turbines and power*, 120(2):410-417, 1998.

[6] Chulwoo Jung, Bogdan I Epureanu, Sanghum Baik, and Marcus B Huffman. Nonlinear reduced order models for the structural dynamics of combustor systems with prestress and friction. *Journal of Computational and Nonlinear Dynamics*, 10(1):011009, 2015.

[7] Walter C Hurty. Dynamic analysis of structural systems using component modes. *AIAA journal*, 3(4):678-685, 1965.

[8] Mervyn CC Bampton and Roy R CRAIG, JR. Coupling of substructures for dynamic analyses. *AIAA Journal*, 6(7):1313-1319, 1968.

[9] S Rubin. Improved component-mode representation for structural dynamic analysis. *AIAA Journal*, 13(8):995-1006, 1975.

[10] C-H Menq, JH Griffin, and J Bielak. The influence of a variable normal load on the forced vibration of a frictionally damped structure. *Journal of Engineering for Gas Turbines and Power*, 108(2):300-305, 1986.

[11] O Poudou and C Pierre. *Blades forced response analysis with friction dampers*. PhD thesis, Ph. D. thesis (University of Michigan), 2007.

[12] Jérôme Guillen. *Studies of the dynamics of dry-friction-damped blade assemblies*. PhD thesis, University of Michigan-Ann Arbor, 1999.

[13] A Kaveh, H Rahami, and P Pezeshky. Geometrically nonlinear analysis of circulant structures using an efficient eigensolution method. *Acta Mechanica*, 223(10):2167-2182, 2012.

[14] Evangéline Capiez-Lernout, Christian Soize, and Moustapha Mbaye. Computational geometrically nonlinear vibration analysis of uncertain mistuned bladed disks. In *ASME Turbo Expo 2014: Turbine Technical Conference and Exposition*, pages V07BT33A001-V07BT33A001. American Society of Mechanical Engineers, 2014.

[15] Sanghum Baik and Mehmet Dede. A forced response method for annular combustors excited by traveling acoustic pressure waves. In *ASME Turbo Expo 2008: Power for Land, Sea, and Air*, pages 345-353. American Society of Mechanical Engineers, 2008.

[16] ANSYS Release. 15.0. ANSYS Theory Reference, 2015.

[17] Suryarghya Chakrabarti. *Modeling of 3D magnetostrategic systems with application to Galfenol and Terfenol-D transducers*. PhD thesis, The Ohio State University, 2011.