Analysis of flow-induced vibrations in turbomachinery by mapping of complex fluid pressures

In this paper we introduce a mapping procedure which facilitates the simulation of flow-induced vibrations in turbomachinery. The transient steady state pressure fluctuations in the flow field (which excite vibrations) are computed in the frequency domain by what are generally referred to as “harmonic CFD” methods where the pressure oscillations are expressed by complex amplitudes. They are mapped using the Fraunhofer software FSIMapper to a structural vibration analysis. A main focus lies in the provision of mapping methods for cyclic symmetric models. The process provides a fast numerical assessment of flow-induced vibrations where the resulting vibration amplitudes can be used for realistic fatigue estimations of flow-excited turbine components. The procedure is applied to a ceramic impeller of a micro gas turbine.


INTRODUCTION
Turbomachinery is deployed in a variety of industrial systems.Flow-induced vibrations can lead to a high noise emission and to blade fatigue which can endanger the integrity of the whole system.The excitation is caused by pressure fluctuations in the flow field generated by interactions between rotating and stationary blade rows [1].Design optimisation for the reduction of product size and weight leads inter alia to a reduction of the distance between rotor blades and stationary guide vanes.This increases the unsteady interactions and thereby the excitation of the already highly loaded blading in the whole flow channel.
Numerical simulations of those excitation forces and vibration responses lead to time-and cost-savings in the prototyping and testing of products since the design process is moved to earlier development phases.However, a classical transient simulation of the steady state flow conditions can be computationally very expensive.A faster simulation approach -the Nonlinear Harmonic method -and its application in a new mapping process are the topic of this paper.

NONLINEAR HARMONIC METHOD
A possibility of approximating the transient steady state behaviour of the flow field is a finite Fourier decomposition into periodic oscillations as developed by He [2], [3].This approach is used in the Nonlinear Harmonic (NLH) method implemented in Numeca's turbine simulation software FINE/Turbo [4].It decomposes the state variable u(r, t) at position r and time t into the time-averaged variable u − (r) and periodic fluctuations around this mean value.
The periodic oscillations are split into K frequency terms which correspond to multiples of the blade passing frequency ω, also called harmonics (see eqn (1)). (1) The factor u ~k (r) ∈  denotes the complex amplitude of the k-th harmonic corresponding to frequency ω k = k ⋅ ω.The complex amplitude can be reformulated as magnitude and phase lag of the oscillation (using Euler's formula).Since the time dependent state variables u are real numbers, u ~k is complex conjugated to u ~k .So the number of free variables is K + 1.
In the NLH method this decomposition is used for time-averaging the unsteady Navier-Stokes equations in order to compute an approximation of the steady state transient solution of the turbomachinery problem.This is analogous to Reynolds averaging, except that the periodic fluctuations are assumed to predominate over the turbulent fluctuations [4].As in the concept of turbulence modelling, the periodic fluctuations contribute additional terms to the timeaveraged Navier-Stokes equation, referred to as deterministic stresses.For the model closure, a transport equation for the unsteady perturbations is obtained by retaining the first-order terms in the basic unsteady flow equations.Casting this first order linearised equation into the frequency domain gives the remaining equations (besides the turbulence model) to close the problem [4].The reader is referred Vilmin [4] for detailed explanation of the NLH method.
As shown by Membera [5] this approach is much faster than the classical transient simulation since the computationally expensive calculation of the initial transient response is avoided.Also, computer storage is minimized since the transient information is represented by a small amount of complex flow data, which needs not to be saved for each time-step.

MAPPING OF COMPLEX PRESSURES
For the computation of flow-induced vibrations the pressure excitations are transferred (mapped) to the structural mesh using the Fraunhofer SCAI tool FSIMapper.The excitations are provided by the NLH method in terms of the time-averaged pressure pand a certain number of harmonic pressures p ~k ∈ , k = 1... K.The basic procedure for simulating the flow-induced vibrations is shown schematically in Figure 1.
In a first step the harmonic turbine flow simulation is performed.FINE/Turbo saves the results in the CFD General Notation System (*.cgns) file, which FSIMapper is able to read.
The harmonic data i.e. the time-averaged pressure pand the complex amplitudes p ~k, k = 1... K of the computed harmonics are read by FSIMapper and are mapped to the structural target mesh.The two meshes do not necessarily need to coincide but represent approximately the same geometric shape.The FSIMapper-algorithms are able to handle different mesh densities or element formulations [6].
FSIMapper exports for each harmonic a loading file which contains the corresponding complex excitation forces on the target mesh.Moreover, it exports a file with the timeaveraged load.The data contained in these files is used in the harmonic structural simulation (stage three in Figure 1).
The first step of the structural vibration analysis is the computation of the static deformations x − at the time-averaged pressure loading.This serves as the base state for the subsequent frequency response steps.In each of those steps, the mapped real and imaginary load data are included and the response (complex deformation x ~k, k = 1 ... K) at the corresponding frequency is calculated.Here also resulting stress, strain, etc. responses are available.
The system loaded by the overall pressure fluctuation given by (2) responds by the linear superposition (inverse Fourier decomposition) of the single responses to the time-averaged pressure and the harmonic pressure fluctuations, as given by eqn (3). ( This is possible due to the linearity of the equations of motion which are solved in the frequency response analyses.

PERIODIC MODELS AND NODAL DIAMETERS
In turbomachinery simulation, periodic models are often used in order to reduce computation times.Usually, the mapping algorithms need a more-or-less coarse match between the source and target meshes.
In order to offer flexible modelling of source and target mesh the mapping procedure presented here provides the possibility of mapping between different periodic sections which in fact represent the same full model.
Figure 2 shows schematically two different cyclic symmetric meshes (black lines) for the use in a data mapping.The mapping algorithm uses the periodicity information to map the data from the source mesh to the -at first glance non-matching -target mesh by "rotating" the data to be present on the virtual full source model (grey lines).
Modal shapes, dynamic excitations, deformations, etc. of periodic systems such as bladed disks usually are described using the concept of nodal diameters (also known as circumferential wave number or cyclic symmetry mode) by Wildheim [7], [8], [9].The nodal diameter ND of a deformation or excitation shape is equal to the number of "inflexion lines" across the disk.The maximum valid nodal diameter ND max is given by eqn (4) where n denotes the number of periodic sections.
(4) A shape of nodal diameter y ∈ [0, ND max ] is abbreviated by NDy. Figure 3 shows a schematic example of a six bladed disk where ND max = 3.In the ND0 shape all blades are excited exactly in phase (shown in Figure 3a).For an even number of blades (as it is the case here) the shape with ND = ND max is characterized by a phase lag of p, i.e. each blade is excited in the opposite direction to its neighbours (Figure 3d).The ND1 resp.ND2 shape is characterized by two resp.four direction changes in the blade row, see Figure 3b and Figure 3c.
To provide the complex data from the periodic model (s = 1) to the remaining sectors s = 2 ... n, eqn (5) is applied. ( Here σ ND is referred to as forward or backward inter blade phase angle which is defined by the excitation shape (respective to the rotation sense).They are closely tied to the nodal diameter ND by the following equations.(6) Eqn (5) corresponds to a rotation in the complex plane: the amplitude of the complex number stays the same and its phase is augmented by σ ND ⋅ (s -1).
If the excitation has ND0 shape the data can be simply copied from one section to another.This is because the matrix from eqn (5) becomes the unit matrix.For an even number of blades and nodal diameter equals ND max the data on the full model can be created by copying the data and applying a sign-change alternatingly.
(1) In both cases the matrix from eqn ( 5) is diagonal, so real and imaginary part are decoupled and can be treated independently.
For vector quantities the derived values from eqn (5) have to be additionally rotated by the angular pitch around the cyclic symmetry axis.

DERIVING EXCITATION AND RESPONDING SHAPE
The specific shape of an excitation, i.e. nodal diameter and forward or backward mode, can be derived by the periodicity of the components in the turbomaehinery system [7], [8], [9].A component of periodicity m leads to excitation frequencies which are multiples of the blade passing frequency ω = m ⋅ Ω, Ω representing the relative rotation speed of the two components towards each other.These multiples are called k-th harmonic frequencies The k-th harmonic associated with an m-periodic part is called engine order E = k ⋅ m and excites the nodal diameter ND in forward mode if there exists a positive integer a with (7) The forward inter blade phase angle is used in eqn (5).If there exists a positive integer b with (8) the excitation has a backward nodal diameter shape ND.So the backward inter blade phase angle is used in eqn (5).
These two conditions correspond to the ZZENF (Zig Zag shaped Excitation line in the Nodal diameter versus Frequency) diagram developed by Wildheim [8] which is shown in Figure 4 as nodal diameter versus engine order diagram.It shows which engine order E causes which excitation shape on an n-periodic part.In the left example diagram shown in Figure 4 the engine orders 2, 4, 8, etc. will lead to an excitation shape of nodal diameter 2 on a six bladed disk.The engine orders 2 and 8 excite in forward mode, engine order 4 in backward mode.
In summary, with the described tools and formulas we are able to derive the excitation shape (as forward or backward nodal diameter) from the turbomachinery configuration (eqns (7) and ( 8)).This information is used to provide the complex data (e.g.exciting pressures) from one periodic section to the corresponding full model (eqn (5)).Thus FSIMapper is enabled to map data between cyclic symmetric meshes which in fact represent the same full model.The resulting deformation shape is the same as the excitation shape.

RESULTS AND DISCUSSION
The presented mapping procedure has been applied during an internal Fraunhofer research project called MAVO TurboKeramik.The goal of this project is to develop a fully ceramic rotor of a micro gas turbine which withstands the operational static and dynamic loading.
The considered turbine configuration is shown in Figure 5.The hot exhaust gas from the combustion chamber is directed radially through the stator which includes 13 blades.The 14bladed rotor is turned by the relaxing gas.The cooled off fluid leaves the turbine in the axial direction in order to co-generate heat.
After the assessment of the static loading, which comprises pre-tension, thermal gradients and centrifugal load, the dynamic excitation and response is considered.
The NLH method in FINE/Turbo computes the average pressure and the first three harmonics on the periodic model (see Figure 6) which together approximate the transient pressure behaviour.The pressure amplitudes decrease with increasing frequency (or harmonic index) whereas the complexity of the pressure distribution is increasing.
For a flow-induced structural vibration analysis these data have to be transferred to the structural mesh.The periodic target model in Abaqus exhibits a different geometry, cf. Figure 7 (right).One section comprises parts of two neighbouring blades with planar cutting faces.In order to provide the periodic data to the full model the excitation shapes of the pressure need to be determined using the ZZENF diagram.The results are shown in Table 1.
Figure 7 visualizes the cyclic symmetric mapping using eqn (5) of the first harmonic complex pressure from the source mesh to the target mesh.The data on the periodic source model (Figure 7 left) exhibits a "ND1-periodicity" and is "rotated" continuously over the periodic boundaries (Figure 7 middle).These data are mapped to the periodic target model (Figure 7 right).
The mapping process uses simultaneously the real and imaginary parts of the data, since the matrix of eqn ( 5) is dense and couples them.An independent mapping of real and imaginary part is only possible for ND0 and ND7.
The complex data differ on the periodic source and target model since they correspond to different blades.This difference is simply induced by the phase lag between the blades resulting from the backward inter blade phase angle.The corresponding amplitude of the pressure excitation is the same.
In the mapping process, include files in Abaqus input format are created to define the complex loading for the frequency response analyses.For each considered harmonic a frequency response step is defined which uses the complex loading files.Abaqus uses the information about the excitation shape in order to build the periodic boundary conditions.
The time-averaged pressure at frequency 0Hz results in a deformation shown in Figure 8a.The maximal deformation caused by the mean pressure is located in the blade tip at trailing edge (near outlet).
In the frequency response simulations a damping has been assumed.Since the magnitude of responses is closely tied to its value, the maximal amplitudes have been scaled to 1.
The deformation magnitudes for each harmonic are shown in Figure 8.Each colour range is scaled to its maximal amplitude.
For the two blade tips at leading and trailing edge the displacement response frequency spectrum is shown in Figure 9. Since the computed results belong to two neighbouring blades, first, the simulation results of the trailing edge blade tip are turned by the angular pitch and then the dedicated inter blade phase angle is added via eqn (5).
The time-averaged pressure is plotted at 0Hz and puts a constant contribution to the oscillation.For each deformation direction the first harmonic has the biggest influence.The responses to the higher harmonics converge to zero.
The inverse Fourier transformation of the spectra shown in Figure 9 leads to the transient steady state behaviour of the blade tips, cf. Figure 10a.As already seen in the spectra, frequency 20800Hz dominates the oscillation in all three degrees of freedom.A coordinate transformation leads to the transient steady state behaviour in cylindrical coordinates, i.e. in radial (r), circumferential (φ) and in axial (z) direction, see Figure 10b.
The resulting amplitudes of stress or strain cycles can be used for a fatigue analysis [10], [11].Figure 11 shows hot spots of high stress oscillation magnitudes at frequency ω 1 = 20.8kHz (i.e. at first harmonic) where failure-probability is highest.The final durability  assessment uses fatigue principles and material properties in order to estimate the lifetime of the dynamically loaded part.

CONCLUSION
This paper presents a mapping procedure for the analysis of flow-induced vibrations in turbomachinery applications.Pressure excitations calculated by harmonic CFD (Computational Fluid Dynamics) methods are mapped to the structural model where a frequency response analysis is performed.In this way the influence of pressure fluctuations in the turbine flow to vibrations of the structure are estimated.The approach is much faster than the classical transient computational procedure, where the deflection is calculated time stepwise until steady state conditions are reached.In a post-processing step a life-time assessment can be performed which takes flow-induced vibrations into account.The results do not encounter aerodynamic damping since the influence of the structural vibration to the flow is not considered.Here a complete coupling of the complex quantities is planned as future work using the vendor-neutral coupling interface MpCCI developed at Fraunhofer SCAI.
The mapping method is designed for the use of periodic models which are often applied in turbomachinery simulations.The algorithm is able to map data between periodic models while not matching geometrically but represent the same equivalent full models.This feature reduces modelling and simulation effort by using the periodicity information of the data.Also, it allows the mapping between periodic and full models.
The procedure was applied to the ceramic impeller of a micro gas turbine.Blade vibration responses to the pressure fluctuations are calculated and hot spots with high fatigue probability are located.
The presented procedure can be easily transferred to other application areas such as electromagnetic induced vibrations in motors or generators.A result of the frequency response analyses can also be the sound pressure level for acoustic assessments.
For source (resp.to the mapping) simulation codes which do not provide harmonic but transient analyses the presented procedure is also applicable.After a transient simulation on the source model, the data of time-steps which build the steady state behaviour are converted using a Fourier transformation to frequency dependent complex amplitudes.They correspond to the harmonic amplitudes used in the mapping method, where the procedure can be continued.This approach has already been implemented in FSIMapper for electromagnetic and CFD applications.

Figure 1 :
Figure 1: Scheme of mapping procedure.The time-averaged pressure p − and the complex amplitudes p ~k, k = 1... K as result of the harmonic Computational Fluid Dynamics (CFD) simulation are transferred via the software FSIMapper to the target structural mesh.The resulting files contain the boundary conditions of the structure in the target code format.They are used in the harmonic structural analysis in order to compute the time-averaged displacement x − and the complex displacements x ~k, k = 1... K for each considered harmonic excitation

Figure 2 :
Figure 2: Mapping of data between different periodic sections (black lines) which represent the same full model (grey lines).The data is mapped from the left source model to the right target model by "rotating" the data to be present on the virtual full source model

Figure 3 :
Figure 3: Mode shapes of a six bladed disk for all valid nodal diameters ND

Figure 4 :
Figure 4: Derived Zig Zag shaped Excitation line in the Nodal diameter versus Frequency (ZZENF) diagram for a six (left) and nine (right) bladed disk.Using the diagram the excitation shape (nodal diameter in forward (black lines) or backward (grey lines) mode) on a n-periodic component caused by a certain engine order can be determined

Figure 5 :
Figure 5: Turbomachinery configuration in the Fraunhofer project MAVO TurboKeramik.The hot exhaust gas is directed radially through the 13-bladed stator to the 14-bladed rotor.The outflow leaves the rotor in the axial direction

Figure 6 :
Figure 6: Mean pressure amplitude and the first three harmonic pressure amplitudes (colouring uses different ranges) as a result of the Nonlinear Harmonic method in FINE/Turbo ω 1 = 20.8kHzω 2 = 41.6kHzω 3 = 62.4kHz

Figure 7 :
Figure7: Cyclic symmetry mapping procedure demonstrated for the first harmonic pressure with FSIMapper.The data on the periodic source model (left) is "rotated" using the excitation shape information to build the full source mesh (middle).These data are mapped to the periodic target model (right)

Figure 8 :Figure 9 :
Figure 8: Magnitude of resulting displacements for the time-averaged pressure and the first three pressure harmonics.Colour range scaled to the particular maximum.Simulation result of frequency response analysis in Abaqus ω 1 = 20.8kHzω 2 = 41.6kHzω 3 = 62.4kHz

Figure 10 :Figure 11 :
Figure 10: Transient steady state deformation as linear superposition of the timeaveraged and harmonic responses

Table 1 :
Excitation shapes of the first three harmonics derived by the ZZENF diagram