# Computational Fluid Dynamic Analysis of a Vibrating Turbine Blade

Authors: Osama N. Alshroof; Gareth L. Forbes; Nader Sawalhi; Robert B. Randall; Guan H. Yeoh
Journal: International Journal of Rotating Machinery
 (2012)

---

## Abstract

This study presents the numerical fluid-structure interaction (FSI) modelling of a vibrating turbine blade using the commercial software ANSYS-12.1. The study has two major aims: (i) discussion of the current state of the art of modelling FSI in gas turbine engines and (ii) development of a “tuned” one-way FSI model of a vibrating turbine blade to investigate the correlation between the pressure at the turbine casing surface and the vibrating blade motion. Firstly, the feasibility of the complete FSI coupled two-way, three-dimensional modelling of a turbine blade undergoing vibration using current commercial software is discussed. Various modelling simplifications, which reduce the full coupling between the fluid and structural domains, are then presented. The one-way FSI model of the vibrating turbine blade is introduced, which has the computational efficiency of a moving boundary CFD model. This one-way FSI model includes the corrected motion of the vibrating turbine blade under given engine flow conditions. This one-way FSI model is used to interrogate the pressure around a vibrating gas turbine blade. The results obtained show that the pressure distribution at the casing surface does not differ significantly, in its general form, from the pressure at the vibrating rotor blade tip.


---

## Body

## 1. Introduction

Blade vibration is an inherent characteristic in the operation of gas turbines, due to the passing of the rotor blades through the trailing wakes from the upstream stator blade stage. This unavoidable blade vibration needs to be taken into account in the design of these engines and as such the measurement of blade vibration is needed. Due to the high operating temperatures, pressures and the rotating environment within a turbine engine, the need for noncontact measuring techniques is apparent. A recent research program has revealed that the gas turbine internal casing pressure, and resulting casing vibration, can be used to determine some characteristics of the internal blade vibration [1, 2]. This work, however, relied on the assumptions that the pressure at the casing surface was not significantly different from that at the blade tip and that the pressure profile developed around the rotor blades followed the blade oscillatory motion.Investigation of the pressure at the casing surface and at the tip of a vibrating rotor blade is undertaken within this paper through a “tuned” three-dimensional one-way fluid-structure interaction (FSI) model of a vibrating turbine blade stage using the commercial numerical code ANSYS (ANSYS Mechanical and ANSYS CFX). This paper also introduces and discusses the limitations on numerical modelling of complex structures, such as gas turbine engines, undergoing FSI, and the various simplifications currently employed to model this behaviour.


## 2. When Is FSI Modelling Required?

Most often problems which involve the interaction between fluids and structures can quite sensibly be modelled as  uncoupled  problems within their separate domains. In these uncoupled problems it is assumed that one domain is driven by the other, with the driven domain having no feedback effect on the driving domain. This results in what is called one-way fluid-structure interaction [3]. The choice of which is the driving and driven domain, respectively, is made from holistic knowledge of the problem being investigated. For instance, the acoustic excitation of air being driven by a speaker can, most often, sensibly be modelled without the need to look at the feedback effects of the driven air on the structural deflection of the load speaker surface. That is therefore an example of a classic FSI problem where the coupling effect of the interaction between the fluid and structural domains only needs to be taken into account in one direction, namely, the effect of the structure on the fluid. Another example, conversely, where the fluid is the driving domain and the structure is being driven, with little feedback effect between the domains, is the dynamic wind loading on many civil structures. The dynamic nature of the flow, and gusting of winds, provides variable dynamic driving loads on the structures’ surfaces. However, the motion and displacement of the structure rarely have any influence on the wind loads which drive this motion. Thus the structure can be analysed separately, or uncoupled from the wind, under the given loads of interest. This, however, is not always the case for wind loads on civil structures, as the infamous Tacoma Narrows Bridge collapse is thought to be a result of FSI [4]. The Tacoma Narrows Bridge is believed to have collapsed due to the dynamic motion of the bridge under wind loads which introduced aeroelastic flutter in the bridge deck. The motion of the bridge caused by the wind created structural deflections which further enhanced the wind forces driving its motion, so that the air and the structure were interacting in a feedback system. It is when conditions such as these arise that modelling of the two-way FSI needs to be undertaken. This is often the case in situations where the amplitude of the structural deflection is large and/or the structural system is highly compliant.Analytical modelling of complex structures has limited scope, so that numerical analysis is often required. The need for numerical solutions to most complex engineering problems has driven the development and wide application of Finite Element Analysis (FEA) and Computational Fluid Dynamics (CFD) numerical simulation in the respective structural and fluid domains. One way to numerically model coupled FSI problems is to couple these two computational domains together with independent solvers for each domain, passing boundary conditions from each domain to the other at each computational step [5]. This setup of independent numerical solvers for each domain is the so-called partitioned solver technique; this is in contrast to a monolithic FSI solver procedure where the coupled boundary conditions are solved together in both domains at each time step [6]. It should be noted, however, that a monolithic solver still often requires iteration of the obtained solution at each time step if the system is nonlinear, which is the case for strongly coupled domains.Although two-way fluid-structure coupling is starting to be more widely available in commercial CFD/FEA packages, it still has not reached a level of integration and computational efficiency to enable it to be used for the analysis of complex structure and flow interactions without significant computational and time expense.


## 3. Current FSI Modelling of Gas Turbines

Fluid flow and structural motion within a gas turbine represent a highly complex system to be numerically modelled. Gas turbines consist of multiple rotor and stator blade stages, as shown in Figure1, which depicts a single stator and rotor blade stage, rotating relative to each other and producing complex flow and forces around both the stator and rotor blades. Full investigation of the flow and structural motion in such an engine requires a three-dimensional two-way coupled FSI model to be developed. Full three-dimensional two-way coupled modelling of a gas turbine stage with rotor blade vibration is theoretically possible with today’s commercial numerical packages; however, the computational expense and setting up of such models are prohibitive and do not appear to be within feasible limits for a number of years yet. The limitations of commercial software on the two-way FSI modelling have also been noted by other authors [7, 8]. For instance, a three-dimensional single stator and rotor blade sector with two-way coupling on the rotor blade surface using ANSYS Workbench 12.1 may take months for a converged solution. As another example, the model within this paper required approximately 4 weeks to obtain a converged solution. Two-way coupling, on average, can take an order of magnitude longer in computational time when compared to a one-way model. One-way modelling is no more onerous than a pure CFD model, except for the more complicated flow that can develop due to the structural motion and the transient solution method which is often required.Figure 1
Schematic view of gas turbine stage, showing the hub, the external casing, and the blades (40 stator blades and 40 rotor blades).Adding to this computational time is the significant increase in computational expense which arises when modelling a turbine stage, as the entire stage often needs to be included in the model. In order to obtain a complete analysis, the number of blade sectors required to be modelled should be the lowest common ratio of the number of rotor to stator blades. The need to model this increased domain can be explained with the aid of Figure2. In Figure 2(a) a two-dimensional blade cascade is shown for a stage with equal numbers of rotor and stator blade numbers (note that three sectors have been shown for clarity). This allows, with appropriate boundary conditions, a single sector to be modelled. Flow from the stator side of the rotor/stator junction can flow into the rotor side (highlighted in yellow). In Figure 2(b) a noninteger pitch ratio is shown. In this case, continuity of the flow from the stator to the rotor is not provided. Therefore in order to fully model a turbine stage an equal width stator and rotor domain need to be created. In the case of most real turbine stage designs, noninteger pitch numbers are desirable which can result in a large computational domain. For instance, a stator of 13 blades and a rotor of 17 blades would require the entire stage to be modelled for the continuity of the flow between the stator/rotor. However, the most recent version of ANSYS includes multirow methods, which reduce the computational domain needed to allow for noninteger pitch ratios without the need for modelling the whole domain [8]. It is unclear yet whether these multirow methods result in any loss of accuracy.Figure 2
Modelling of rotor/stator domains requires the lowest common multiple of rotor/stator blades to be modelled for continuity of flow out of one domain and into another. (a) Integer pitch ratio, all flow out of stator domain flows back into rotor domain (b) noninteger pitch ratio, has unequal flow from stator to rotor domains.In addition to often requiring a large computational domain for complete analysis of gas turbine engines, the unsteady nature of the flow within the turbine stage also requires the use of a transient solution procedure. For instance, in reference to Figure1, the rotor stage moves relative to the stator stage and causes oscillating flow and forces on the rotor blades. Therefore large computational domains and transient solution procedures have meant that a number of simplifying FSI modelling techniques for gas turbine engines have arisen.To the best of the authors’ knowledge, a complete three-dimensional transient two-way FSI of a turbine engine including blade vibration has not been achieved. Some notable simplified models which have been developed recently are(i)
two-way transient three-dimensional FSI model of an axial hydro turbine; this does not include the preceding stator stage, so it only considers static deformation and does not include blade vibration [7];(ii)
pseudo two-dimensional two-way transient FSI model of a gas turbine rotor blade vibration [9];(iii)
one-way transient three-dimensional FSI model of a gas turbine blade undergoing vibration [10];(iv)
significant numerical and experimental investigation of the one-way transient three-dimensional FSI model of a turbine blade undergoing vibration and the effect on blade flutter and casing tip gap [11].The use of a harmonic balance approach to the solution of the fluid dynamics is by far the most efficient for steady state harmonic response forms [12]. Harmonic balance methods, however, do not seem to be available in commercial software packages at this stage. The increase in efficiency is reduced as full transient flow or motion should be modelled. Neglecting the full transient flow effects may not produce the actual blade motion, but as the dominant pressure pulses on the turbine blades during operation are from harmonic orders of the rotational speed, the dominant motion effects should still be present.Figure3 shows schematically the increase in computational expense for different FSI modelling approaches, with each employing simplifying assumptions from the full three-dimensional two-way FSI model. One-way FSI models give the option of reduced complexity by enforcing the vibration on the blade surface and solving the moving boundary condition CFD model of a fluid sector with a vibrating rotor blade. However, these models do not take into account the preceding stator blade row effects. Two-dimensional models give an obvious reduction in complexity by limiting the computational domain, but this eliminates the three-dimensional behaviour of the flow. This simplification does not allow resolution of the complex three-dimensional flow behaviour in the gap between the tip of the blade and casing wall which is the aim of this paper.Figure 3
Schematic of different FSI numerical modelling procedures for gas turbines. The simplified tuned model used within this work lies between 3D One-way FSI and 2D Two-way FSI in computational expense for the same mesh density.The model presented within this paper is a “tuned” three-dimensional one-way FSI model which enforces the rotor blade motion but also includes the preceding stator blade domain and “tunes” the blade motion to take into account the dynamic pressure on the rotor blade due to the stator passing period. This model was developed to give the computational efficiency of a one-way model and to include the effects of the two-way interaction of the blade motion and the driving pressure fluctuations on the rotor blade. Once this model is developed, the investigation of the pressure distribution at the rotor blade tip and the casing surface is undertaken.


## 4. Rotor Blade Vibration

As stated earlier, vibration of gas turbine blades is an inherent characteristic of their use and is driven predominantly by the periodic traversing of high and low pressure regions from the preceding stator blade row. This process can be explained, with the aid of Figure1, by conceptualising how rotor blades would need to pass through the flow of the stator rows as the turbine rotates. This blade vibration has been enforced on the rotor blade in the following numerical modelling, forcing the rotor blade section to move  transverse to the turbine’s angle of attack.A three-dimensional gas turbine was modelled as shown in Figure1, with 40 blades on a 100 mm diameter hub rotating with angular velocity, ω, of 500 rad/s. The axial gap between the stator blades and the rotor blades is 10 mm and various radial gaps between the rotor blade tips and casing surface, of 0.5 mm, 1 mm, and 1.5 mm, were tested. In practice the smallest gap possible between the blade tip and casing surface is desirable so the sensitivity to the size of this gap was investigated. An equal number of stator and rotor blades were used to reduce the computational domain (as explained in Section 3), such that the enforced motion on the rotor blades was given the stator passing frequency of 20000 rad/s  with the maximum blade tip peak to peak amplitude of 0.15 mm. In order to reduce the computational time needed to solve this problem, a three-dimensional wedge encompassing one stator blade and another three-dimensional wedge encompassing one rotor blade, as shown in Figure 4(a), were simulated using the boundary conditions outlined in the next section.Figure 4
(a) Three-dimensional blade model showing blade vibration direction. (b) Three-dimensional model showing polylines where the pressure is observed. (c) Zoomed polylines. (d) Polylines and vertical line across blade tip/casing gap where the pressure is observed.


## 5. Boundary Conditions

The following boundary conditions were employed to obtain the time-dependent solution (Table1).Table 1
Boundary conditions for the FSI model as shown in Figure2.
Boundary
Description
Inlet
Inlet velocity = [44.38i − 61.08j + 0k] m/s, in which i, j, and k are the unit vectors for x,y, and z coordinates as shown in Figure 3.
Outlet
Average relative static pressure = 0.
Casing surface
Nonslip wall
Hub
Nonslip wall
Stator blade
Nonslip wall
Rotor blade
Nonslip wall with forced vibration added to represent the response to forces from the periodic traversing of high and low pressure regions from the preceding stator blade row. The displacement of the rotor blade is presented as a sinusoidal displacement as a function of time and position in thez-direction as follows:δ(t,z)=45(-A×sin(20000t-tf))×(|z|-R0)2, in which there are 40 stator blades, tf = 1.2868 s is a phase offset calculated from the applied pressure, and R0 = 0.1 m is the radius of the hub
The two sides of the wedge
Defined using rotational periodicity
The interface between the two domains
Defined as transient rotor stator [13]It should also be noted that the rotor domain is rotating at 500 rad/s, as shown in Figure1, whereas the casing surface is defined as a stationary nonslip wall.


## 6. Rotor Blade Structural Analysis and Property Settings

The forces present on the rotor blade were first calculated under the flow conditions it would experience without any blade vibration. It is the resulting motion of the rotor blade, calculated as a result of these forces, which is used in the later enforced motion of the CFD model. The response of the rotor blade under these forces is discussed in this section. The material properties of the blade were altered to obtain the desired maximum blade tip displacement of 0.15 mm. What is undertaken within this section is a classic decoupled structural analysis under the incoming flow conditions.A judicious choice of damping is also made to suppress higher modes of vibration for two reasons. Firstly, the first mode of vibration is often the most dominant and thus the most important to consider. Secondly, limiting the motion to only one mode of vibration allows for simpler motion to be observed when investigating the flow between the blade tip and casing surface (this avoids the necessity of low-pass filtering the calculated results).


### 6.1. Blade Geometry

The blade under consideration is shown in Figure5; it has a length of 150 mm and a chord length of 50 mm with an aerodynamic profile of NACA4506.Figure 5
150 mm long rotor blade model.


### 6.2. Natural Frequencies and Mode Shapes

The aim of the initial finite element model of one rotor blade is to tune the structural properties of the blade (mass densityρ, modulus of elasticity E) so the first natural frequency of the blade is placed roughly at -/+ 20% from the excitation frequency (blade pass frequency) as shown in Figure 6. Normally the separation has to be even larger than 20% to ensure that the system is not excessively excited. For the purpose of this study the 20% limit was selected to drive the system close to its first natural frequency and obtain reasonable blade deflections.Figure 6
Nondimensional amplitude for the first mode with different driving frequencies. Shown are the chosen stiffness-driven limit 20% below the natural frequency and the mass-driven limit 20% above.To simulate the connection of the blade to the rotor, the rotor blade FEA model was fixed at one end. Thus the model in its simplest configuration represents a cantilever beam with an aerofoil cross-section. The analytical expression to calculate the natural frequencies of a cantilever beam is given in (1) [14]. The analytical expression of (1) allows observation of the main properties which control the behaviour of the blade natural frequency and thus how it reacts to external excitations:(1)ωi=(βil)2EIρAl4,
where,  i: natural frequency Index (i=1,2,3,…);  ωi: natural frequency (rad/s);  βil: boundary condition constants; namely, β1l=1.875104,β2l=4.694091,β3l=7.854757,andβ4l=10.995541;  E: Young’s modulus of elasticity (Pa);  I: second moment of area of the beam (blade) cross-section about the y-axis (m4);  ρ: mass density (kg/m3);  A: cross-sectional area (m2);  l: length (m).For a rotor with 40 blades, and a speed of 500 rad/s (79.58 Hz, 4774.68 rpm), the blade pass frequency, BPF (rotor speed times the number of blades) would be 20,000 rad/s (3183.1 Hz).The selection of the first natural frequency at 20% from the BPF means that the 1st natural frequency would be either 2546.5 (20% below) for a stiffness controlled (driven) system or 3819.7 Hz (20% above) for a mass controlled (driven) system. This also means that the vibration response will be in phase with the exciting force for the stiffness controlled system (180 degrees out of phase if mass was controlled). Depending on the engine of interest and the dominant engine order, both stiffness and mass controlled blade vibrations can be found in turbine engines.To achieve a first natural frequency around 3820 Hz, the mass density of the material (ρ) and the Young’s modulus of elasticity (E) were adjusted. The modulus of elasticity per mass density, which is often referred to as the specific stiffness or specific modulus needs to be optimized. In aerospace applications the specific modulus is required to be high indicating a minimum structure weight (hollow structure) with minimum deflections (higher modulus of elasticity or stiffness). To achieve a high specific stiffness structure with a first natural frequency around 3820 Hz, the mass density (ρ) was selected as 9 kg/m3, while the modulus of elasticity (E) was increased to 3 × 1011 Pa. The first four natural frequencies of the blade are listed in Table 2, while corresponding mode shapes of the first three natural frequencies as shown in Figure 7. The most important of these is the first mode (pure bending). The second mode (mixed torsion and bending) is roughly 6-times higher, which agrees with rough calculations (β2l/β1l)2, based on (1).Table 2
The first four natural frequencies of one rotor blade.
Mode
Frequency (Hz)
(
1
)
3831.2
(
2
)
23112
(
3
)
25322
(
4
)
48583Figure 7
The first 3 mode shapes of the blade: top (pure bending), middle, a mixture of bending and torsion, and bottom, pure torsion.


### 6.3. Damping

Modal damping, proportional to the mass and/or stiffness matrix, was added to the finite element model. To suppress the second and higher modes and retain the first mode, the damping factor (ψ), defined in (2) [15], was adjusted, in the material properties, to achieve this purpose:
(2)ψ=2ζiωi,
where, ζi is the damping ratio of mode i; ωi is the natural frequency of mode i.For a 2% damping in the first mode, (ψ) has been calculated to be 1.67×10-6 s. This value of (ψ) will correspond to a damping ratio of around 24% (very heavily damped) at the second mode. This gets even higher for higher modes, thus suppressing their effect. The suppression of the response of higher frequency modes was employed to avoid numerical problems in the simulation. The natural frequencies were outside the range directly excited by blade passage, and responses would be filtered to exclude higher frequency responses in any case, as has been done in Figures 14 and 15. Additionally the blade motion is related to the phase modulation of the blade pass signal, and this has been arranged to be close to the first mode, so that higher order modes are well outside the frequency band demodulated.


### 6.4. Deriving Blade Forced Vibration Equation (Amplitude and Phase)

The rotor blade will be excited at the blade pass frequency (3183.1 Hz). The forced vibration equation for a cantilever beam along its lengthz and its change in time t, δ(t,z), is given in (3) [14]. The direction axial direction z for the blade of interest is shown in Figure 5. For use in this model the deflection of the blade along its length z has been simplified to a quadratic equation, as shown in (4). Application of (4) will result in a blade deflection shape which will look similar to the deflection shape shown in the top image in Figure 7. This deflection of the rotor blade is zero at the rotor hub and maximum at the blade tip whilst oscillating around its equilibrium position in time. The spatial distribution of the deflection is given by the “mode shape equation” and the osculation in time is given by the “time varying” section highlighted in (4):
(3)δ(t,z)=C1︸Constant×[sinβ1z-sinhβ1z-α1(cosβ1z-coshβ1z)]︷Mode  shape  equation×(-0.00015×sin(500×40×t-1.2868))︷Time  varying,
where,
(4)α1=(sinβ1l+sinhβ1lcosβ1l+coshβ1l),β1l=1.875104,δ(t,z)=1(0.25-0.1)2︸Scaling  factor×(z-0.1)2︷Mode  shape  equation×(-0.00015︸Max  deflection×sin(500×40︸BPF×t-1.2868︸Phase))︷Time  varying.
The terms given in (4) are described as follows.(i)
t: simulation time steps defined in the CFD solver.(ii)
z: rotor blade longitudinal axis. Note that z values are defined as 0.1≤z≤0.25. The limits of the blade z given such that zero deflection is obtained at the attachment to the rotor (0.1 m from the centre of the rotor) and maximum deflection at the tip (0.25 m). The scaling factor was incorporated to give a maximum deflection of 0.15 mm.(iii)
The rotor blade is already deflected under the average static pressure, that is, the blade starting position assumes that it has already deflected under the static pressure.(iv)
The phase of (3) was derived from running the CFD model on its own and observing the averaged pressure on one-third of the free end of the blade. The period for the pressure variation is 3183.1 Hz (BPF). The phase angle was derived as 1.2 rad and was used in the forced vibration equation to ensure that the response (deflection) is in phase with the pressure (excitation). However, as the model was run with the forced equation, it is the deflection of the blade that creates the pressure variation rather than the other way around (this is the limiting assumption of any one-way FSI model).(v)
Amplitude of the blade vibration was arbitrarily selected as±0.15 mm. This value corresponds to the blade response under a tip pressure of around 1000 Pa. The harmonic response (total deflection) under a constant 1000 Pa pressure applied to one surface is shown in Figure 8 (100–30000 Hz, dB—log scale) and zoomed in Figure 9 (100–5000 Hz linear). Figure 8, which displays the blade response over a wide range, shows that the dominant natural frequency in the model comes from the first natural frequency. In Figure 9, the maximum blade defection at the tip, when the blade is excited at the BPF, is shown to be 0.15 mm. The static displacement is 0.05 mm, while the resonant deflection will be about 0.6 mm.Figure 8
Harmonic response total maximum deflection Decibels (re1.3×10-6) for a 1000 Pa pressure acting on one surface of the blade.Figure 9
Harmonic response linear total maximum deflection for a 1000 Pa pressure acting on one surface of the blade.


## 6.1. Blade Geometry

The blade under consideration is shown in Figure5; it has a length of 150 mm and a chord length of 50 mm with an aerodynamic profile of NACA4506.Figure 5
150 mm long rotor blade model.


## 6.2. Natural Frequencies and Mode Shapes

The aim of the initial finite element model of one rotor blade is to tune the structural properties of the blade (mass densityρ, modulus of elasticity E) so the first natural frequency of the blade is placed roughly at -/+ 20% from the excitation frequency (blade pass frequency) as shown in Figure 6. Normally the separation has to be even larger than 20% to ensure that the system is not excessively excited. For the purpose of this study the 20% limit was selected to drive the system close to its first natural frequency and obtain reasonable blade deflections.Figure 6
Nondimensional amplitude for the first mode with different driving frequencies. Shown are the chosen stiffness-driven limit 20% below the natural frequency and the mass-driven limit 20% above.To simulate the connection of the blade to the rotor, the rotor blade FEA model was fixed at one end. Thus the model in its simplest configuration represents a cantilever beam with an aerofoil cross-section. The analytical expression to calculate the natural frequencies of a cantilever beam is given in (1) [14]. The analytical expression of (1) allows observation of the main properties which control the behaviour of the blade natural frequency and thus how it reacts to external excitations:(1)ωi=(βil)2EIρAl4,
where,  i: natural frequency Index (i=1,2,3,…);  ωi: natural frequency (rad/s);  βil: boundary condition constants; namely, β1l=1.875104,β2l=4.694091,β3l=7.854757,andβ4l=10.995541;  E: Young’s modulus of elasticity (Pa);  I: second moment of area of the beam (blade) cross-section about the y-axis (m4);  ρ: mass density (kg/m3);  A: cross-sectional area (m2);  l: length (m).For a rotor with 40 blades, and a speed of 500 rad/s (79.58 Hz, 4774.68 rpm), the blade pass frequency, BPF (rotor speed times the number of blades) would be 20,000 rad/s (3183.1 Hz).The selection of the first natural frequency at 20% from the BPF means that the 1st natural frequency would be either 2546.5 (20% below) for a stiffness controlled (driven) system or 3819.7 Hz (20% above) for a mass controlled (driven) system. This also means that the vibration response will be in phase with the exciting force for the stiffness controlled system (180 degrees out of phase if mass was controlled). Depending on the engine of interest and the dominant engine order, both stiffness and mass controlled blade vibrations can be found in turbine engines.To achieve a first natural frequency around 3820 Hz, the mass density of the material (ρ) and the Young’s modulus of elasticity (E) were adjusted. The modulus of elasticity per mass density, which is often referred to as the specific stiffness or specific modulus needs to be optimized. In aerospace applications the specific modulus is required to be high indicating a minimum structure weight (hollow structure) with minimum deflections (higher modulus of elasticity or stiffness). To achieve a high specific stiffness structure with a first natural frequency around 3820 Hz, the mass density (ρ) was selected as 9 kg/m3, while the modulus of elasticity (E) was increased to 3 × 1011 Pa. The first four natural frequencies of the blade are listed in Table 2, while corresponding mode shapes of the first three natural frequencies as shown in Figure 7. The most important of these is the first mode (pure bending). The second mode (mixed torsion and bending) is roughly 6-times higher, which agrees with rough calculations (β2l/β1l)2, based on (1).Table 2
The first four natural frequencies of one rotor blade.
Mode
Frequency (Hz)
(
1
)
3831.2
(
2
)
23112
(
3
)
25322
(
4
)
48583Figure 7
The first 3 mode shapes of the blade: top (pure bending), middle, a mixture of bending and torsion, and bottom, pure torsion.


## 6.3. Damping

Modal damping, proportional to the mass and/or stiffness matrix, was added to the finite element model. To suppress the second and higher modes and retain the first mode, the damping factor (ψ), defined in (2) [15], was adjusted, in the material properties, to achieve this purpose:
(2)ψ=2ζiωi,
where, ζi is the damping ratio of mode i; ωi is the natural frequency of mode i.For a 2% damping in the first mode, (ψ) has been calculated to be 1.67×10-6 s. This value of (ψ) will correspond to a damping ratio of around 24% (very heavily damped) at the second mode. This gets even higher for higher modes, thus suppressing their effect. The suppression of the response of higher frequency modes was employed to avoid numerical problems in the simulation. The natural frequencies were outside the range directly excited by blade passage, and responses would be filtered to exclude higher frequency responses in any case, as has been done in Figures 14 and 15. Additionally the blade motion is related to the phase modulation of the blade pass signal, and this has been arranged to be close to the first mode, so that higher order modes are well outside the frequency band demodulated.


## 6.4. Deriving Blade Forced Vibration Equation (Amplitude and Phase)

The rotor blade will be excited at the blade pass frequency (3183.1 Hz). The forced vibration equation for a cantilever beam along its lengthz and its change in time t, δ(t,z), is given in (3) [14]. The direction axial direction z for the blade of interest is shown in Figure 5. For use in this model the deflection of the blade along its length z has been simplified to a quadratic equation, as shown in (4). Application of (4) will result in a blade deflection shape which will look similar to the deflection shape shown in the top image in Figure 7. This deflection of the rotor blade is zero at the rotor hub and maximum at the blade tip whilst oscillating around its equilibrium position in time. The spatial distribution of the deflection is given by the “mode shape equation” and the osculation in time is given by the “time varying” section highlighted in (4):
(3)δ(t,z)=C1︸Constant×[sinβ1z-sinhβ1z-α1(cosβ1z-coshβ1z)]︷Mode  shape  equation×(-0.00015×sin(500×40×t-1.2868))︷Time  varying,
where,
(4)α1=(sinβ1l+sinhβ1lcosβ1l+coshβ1l),β1l=1.875104,δ(t,z)=1(0.25-0.1)2︸Scaling  factor×(z-0.1)2︷Mode  shape  equation×(-0.00015︸Max  deflection×sin(500×40︸BPF×t-1.2868︸Phase))︷Time  varying.
The terms given in (4) are described as follows.(i)
t: simulation time steps defined in the CFD solver.(ii)
z: rotor blade longitudinal axis. Note that z values are defined as 0.1≤z≤0.25. The limits of the blade z given such that zero deflection is obtained at the attachment to the rotor (0.1 m from the centre of the rotor) and maximum deflection at the tip (0.25 m). The scaling factor was incorporated to give a maximum deflection of 0.15 mm.(iii)
The rotor blade is already deflected under the average static pressure, that is, the blade starting position assumes that it has already deflected under the static pressure.(iv)
The phase of (3) was derived from running the CFD model on its own and observing the averaged pressure on one-third of the free end of the blade. The period for the pressure variation is 3183.1 Hz (BPF). The phase angle was derived as 1.2 rad and was used in the forced vibration equation to ensure that the response (deflection) is in phase with the pressure (excitation). However, as the model was run with the forced equation, it is the deflection of the blade that creates the pressure variation rather than the other way around (this is the limiting assumption of any one-way FSI model).(v)
Amplitude of the blade vibration was arbitrarily selected as±0.15 mm. This value corresponds to the blade response under a tip pressure of around 1000 Pa. The harmonic response (total deflection) under a constant 1000 Pa pressure applied to one surface is shown in Figure 8 (100–30000 Hz, dB—log scale) and zoomed in Figure 9 (100–5000 Hz linear). Figure 8, which displays the blade response over a wide range, shows that the dominant natural frequency in the model comes from the first natural frequency. In Figure 9, the maximum blade defection at the tip, when the blade is excited at the BPF, is shown to be 0.15 mm. The static displacement is 0.05 mm, while the resonant deflection will be about 0.6 mm.Figure 8
Harmonic response total maximum deflection Decibels (re1.3×10-6) for a 1000 Pa pressure acting on one surface of the blade.Figure 9
Harmonic response linear total maximum deflection for a 1000 Pa pressure acting on one surface of the blade.


## 7. CFD Solution Method

The air flow is considered to be compressible. A transient solver with a second-order backward Euler differencing scheme is used for the transient term and the SST turbulence model with a high-resolution scheme is employed for the advection terms in the transport equations governing mass, momentum, and energy. At each time step, internal iterations are carried out until the residuals for mass, momentum, and energy residuals fall below the criterion 10−6. A time step of 1 × 10−6 s is used.


### 7.1. Details of the Grid

As the blade deforms, the fluid mesh needs to follow the displacement of the blade; therefore the mesh has to be expanded and compressed. As shown in Figure10, hexahedral elements are created around the blade as well as for the rest of the computational domain. In order to adequately resolve the complex flow around the blade, an increased number of elements were placed around the surface of the blade. A total of 1.7 million elements was generated for the whole domain.Figure 10
View of a structured grid at (a) horizontal plane crossing the blade tip and (b) cross-section plane located at the blade centre.


## 7.1. Details of the Grid

As the blade deforms, the fluid mesh needs to follow the displacement of the blade; therefore the mesh has to be expanded and compressed. As shown in Figure10, hexahedral elements are created around the blade as well as for the rest of the computational domain. In order to adequately resolve the complex flow around the blade, an increased number of elements were placed around the surface of the blade. A total of 1.7 million elements was generated for the whole domain.Figure 10
View of a structured grid at (a) horizontal plane crossing the blade tip and (b) cross-section plane located at the blade centre.


## 8. Investigation of the Pressure between the Blade Tip and Casing Surface

Now that the computational model has been described and set up, the pressure between the vibrating rotor blade tip and the casing surface will now be investigated for a range of gaps between the blade tip and casing surface, namely, 0.5, 1, and 1.5 mm. The pressure between the blade tip and casing surface is being investigated to find if there is any significant difference between the pressure at the blade tip and that at the casing surface and whether the pressure profile around the blade follows its motion as it is oscillating around its equilibrium position.During rotor rotation and the subsequent enforced rotor blade oscillation, a high pressure region, as expected, is maintained underneath the blade aerofoil and a low pressure region is maintained above the blade aerofoil. This is also observed at the blade tip; see Figure11(I) which shows the pressure contours on a surface located at a radius equal to the blade tip, for the case of 0.5 mm gap, at maximum-minimum-mean values of blade deflection. Interestingly, the shape of the pressure contours at the blade tip cross-section surface, Figure 11(I), is similar to the pressure contours at the casing wall, Figure 11(II); however, the magnitude of the pressure varies slightly.Pressure contours of surface located at (I) blade tip radius and at (II) the casing surface, for the case of 0.5 mm gap, when the blade is at (a) maximum deflection, (b) mean deflection, and (c) minimum deflection.
(a)
(b)
(c)As a result of the low and high pressure regions around the blade, the flow escapes from the high pressure region to the low pressure region through the gap between the blade tip and the casing. This causes a blade tip vortex, which has been very well documented by previous researchers [11]. The development of this tip vortex can be seen in Figure 12(b), and is seen more clearly for each blade gap in Figure 13. It can be noted that the tip vortex becomes more developed for increasing tip gap. One aim of this study was to check that this vortex did not greatly affect the relationship between the pressure at the blade tip and the casing inner wall, at least for small tip clearances.Pressure contours of surface located at (a) midway between blade tip and hub and (b) blade tip radios, for the case of 0.5 mm gap between the blade tip and casing surface.
(a)
(b)Figure 13
Streams of fluid crossing the gap between the top of the blade and the casing for the case of (a) 0.5 mm gap, (b) 1.0 mm gap and (c) 1.5 mm gap. Pressure distribution in the longitudinal plane is also shown.Pressure across the three polylines, shown in Figure4(d), when the blade deflection is (a) at a maximum, (b) at a mean, and (c) at a minimum, for the case of 0.5 mm gap. (d), (e), and (f) are the first spatial harmonic of the pressure across the same polylines for the respective blade displacements. Note the angle given is relative to the centre of the aerofoil.
(a)
(b)
(c)
(d)
(e)
(f)First spatial harmonic of the pressure across each of the three polylines, shown in Figure4(d), when the blade deflection goes through its range of motion from maximum deflection, through mean deflection and to its minimum deflection position. (a) Polyline at tip of blade (b) polyline between surface and blade tip, and (c) polyline at casing surface.
(a)
(b)
(c)The angle of attack of the stator and rotor blades was selected to maintain smooth flow around the blade, with no separation, as shown in Figure12(a). This smooth motion of the flow continued along the length of the rotor blade, excluding the top 5 mm of the blade tip, where separation of the flow was observed, as can be seen in Figure 12(b).Now, as seen in Figure11, the pressure appears to have a very similar form on the casing surface as it does at the blade tip. Whether or not the pressure profile around the rotor blade follows the blade motion as it vibrates is however not easily investigated by the use of Figure 11. This will now be investigated in more depth by looking at the pressure along three lines across the blade sector width during the blade’s oscillation. These three lines were located at the blade tip, midway between the tip and at the casing surface; see Figure 4.The pressure profile across the width of the rotor blade sector, over these lines, for three instances in time when the blade tip deflection is at its maximum, at its mean, and at its minimum deflection can be seen in Figures14(a), 14(b), and 14(c), respectively. Figures 14(d), 14(e), and 14(f) then display the first harmonic, in space, of the pressure across the polylines for the corresponding displacements. Figures 14(d)–14(f) are simply the result of low pass filtering Figures 14(a)–14(c). Although the magnitude of the pressure distributions in Figures 14(a), 14(b), and 14(c) is not the same for the three pressure lines, the pressure distributions are generally of a similar shape and only differ at the top and bottom surface of the blade. Comparison of the dominant first harmonic of the pressure distribution shows an almost identical match between the pressures at the casing surface and that at the blade tip, shown in Figures 14(d), 14(e), and 14(f). It is definitely evident from observation of Figure 14 that the pressure distribution at the casing surface across the blade passage is very similar to that measured at a distance corresponding to the blade tip.Figure14 compared the pressure distributions across the rotor blade sector for different tip/casing gaps at a single instant during the blade oscillation. Figure 15 looks at the change in pressure distribution across the rotor blade sector at a single position, for example, blade tip during the blades range of motion. From observation of Figure 15 it can be seen that there is little phase shift between these first harmonics of pressure distribution; however, there is a relatively significant change in magnitude. Hence the pressure distribution is amplitude modulated by a pressure oscillation at a frequency corresponding to the blade passing period. This amplitude modulation of the pressure is something which should be considered in further analysis of the internal pressure signal of a gas turbine engine for possible use as a measure of the blade vibration.The relationship of the pressure across the gap between the blade tip and casing surface is illustrated in Figure16. See Figure 4 for the radial line for which the pressure distribution is plotted for the three blade deflections. The general shape of the pressure distribution along the line can be seen to be similar across the blade tip gap. Although the pressure distribution is similar in shape across the gap for the three blade deflection positions, the magnitude can be seen to differ. It is thought that this magnitude difference is a result of the different tip leakage (or flow across from one side of the blade to the other across the gap) resulting from the rotor and stator being at different angular offsets from each other at these times.Pressure distributions along a line located at a radial direction from the blade tip to the casing, for the case of 0.5 mm gap, at (a) maximum, (b) mean value and (c) minimum blade displacement.
(a)
(b)
(c)Figures17 and 18 show the harmonic nature of the pressure at both the casing surface and the blade tip. In Figure 17 the pressure is observed over two stator passing periods at the casing surface as if a probe were inserted in the casing. The rotor blades will therefore be rotating past this point, while also oscillating back and forth about their equilibrium position. For a stationary probe located at the casing surface, the pressure around the blade will be modulated by both the vibrating motion and the motion of the rotation of the blade stage. This will cause higher harmonics of the stator passing frequency [1], which will be out of phase with each other. The pressure traces against time for the three different gaps between the blade tip and casing surface all seem to show this trend of higher harmonics of stator passing frequency. The first and second harmonics of the frequency corresponding to the stator passing period, 0.0031 s, are clearly visible.Pressure at a fixed point located at the casing wall for (a) 0.5 mm, (b) 1.0 mm, and (c) 1.5 mm gap.
(a)
(b)
(c)Pressure at the two points located at the mean value of the oscillated blade tip and at a corresponding point on the casing wall, in which both points are rotating with the same angular velocity as the rotor, for (a) 0.5 mm, (b) 1.0 mm, and (c) 1.5 mm gap.
(a)
(b)
(c)Figure18 shows the pressure at both the blade tip and the casing surface, as if a probe were placed in the fluid which followed the blade as it rotated. The pressure traces which follow the blade as it rotates also tend to vary harmonically with the stator passing period, dominated by the first harmonic and higher harmonics that are all generally in phase with each other. Once again the resulting pressures at the blade tip and casing surface show similar trends. The results for the 1.5 mm gap, however, do vary more significantly, which is believed to be due to the greater leakage of flow around the blade tip.


## 9. Conclusions

The work presented in this paper aimed to present the current state of the art of FSI numerical modelling of gas turbine engines and to develop a computationally efficient one-way FSI model which can be developed within the current commercial ANSYS software package. With this numerical model the pressure differences between the vibrating rotor blade tip and the casing surface were then investigated with two aims:(i)
to show if the pressure at the blade tip was significantly different from that at the casing surface;(ii)
to investigate if the pressure profile around the rotor blade followed its motion as it vibrates about its equilibrium position.A three-dimensional CFD “tuned” one-way FSI model of a single stator and rotor stage, with vibrating rotor blades, was constructed. This model enforced the structural motion of the rotor blade. The rotor blade motion was calculated from the pressure forces which would be present if there was no blade oscillation but still including the rotor/stator relative rotation. This allows the correct amplitude and phase of the structural motion to be used in the one-way CFD model. Setting up an FSI model in this manner gives the reduced computational expense of a one-way FSI model, but allows the correct deflection of the structure to be applied to the fluid domain boundaries.Although the pressure at the casing surface did indeed differ from that at the rotor blade tip, the results presented in Figures11, 14, 15, and 18 show generally that the pressure at the casing surface and rotor blade tip have similar trends and are representative of each other and are almost identical for the dominant first spatial harmonic of the pressure distribution. The difference which exists between these two locations is believed to be as a result of the flow which leaks over the blade tip and the blade tip vortices which develop. From the results presented it is also concluded that the pressure profile around the rotor blade tends to follow its motion. It is also shown that the internal pressure of the turbine is amplitude modulated by the blade passing period. However, due to both the integer number of stator to rotor blades and the relatively small dynamic blade deflection, further investigation into the importance of these two effects would need to be undertaken to make a more conclusive assessment.



---

*Source: 246031-2012-11-12.xml*