Numerical study of characteristic modes and frequencies of flow in high-speed compressors

In this paper, we present the newly developed open-source density-based solver pisoCentralDyMFoam and investigate an application of the Proper Orthogonal Decomposition (POD) algorithms for industrial turbomachinery-related problems. POD is implemented in Apache Spark framework for distributed data processing. This solver is based on hybrid Kurganov-Tadmore/PISO scheme. The research was conducted for geometry close to its real prototype with known resonance frequencies. The solver was previously validated on simple industrial case ERCOFTAC centrifugal pump. The POD coefficient matrices were constructed using data set of the snapshots representing each saved time-step of the whole Navier-Stokes numerical model. Several hundred consecutive snapshots of the static pressure field on the impeller wheel surface as well as the velocity field were used for computation of the POD modes. The eigenvalues determined by POD method corresponds to the kinetic energy contained in each mode. The spatial coefficients represents contribution of each elementary volume to the whole mode and helps to locate the region having influence on the mean flow at specific frequencies after Fast-Fourier-Transform (FFT) applied to time-dependent coefficients of the decomposition. The POD modes were sorted by kinetic energy and the zeroth mode was most energetic representing mean flow with relatively small amplitudes. The described concept was extensively validated using computationally cheap 2D-case and then extended to the highspeed centrifugal pump of the small-scale turbojet. It was found that the third mode of the flow has first peak at 12970 Hz right between 2 construction resonance points at 12000 Hz and 13700 Hz. The third and fourth modes represent pressure fluctuations in the wakes region in the diffuser behind vanes. The demonstrated approach allows engineers to analyze flow dynamics more effectively compare to traditional FFT at certain points or cross-section. In addition, it can be useful for data compression and Reduced-Order Model development.


Introduction
In the turbomachinery-related field, Epureanu et al. [1] elaborated theoretical aspects and application of the POD models itself.The authors investigated the influence of the amount of POD modes on resulting aerodynamic characteristics if the airfoil cascade.They noticed that subsonic and transonic simulations require different degrees of freedom of the POD model to predict correctly flow dynamics for all interblade phase angles.The early work of the ONERA research group [2] was dedicated to extraction of the most energetic modes of the flow in the centrifugal pump to catch rotor-stator interaction using most reliable and robust techniques.Rochuon et al. [2] investigated both Fourier-transform and POD to obtain information about flowfield.They found that POD allows getting first most energetic modes representing 98% of the flow kinetic energy by nine times less harmonics.Clark et al. [3] demonstrated the application of POD-analysis for turbomachinery-related dynamics problems very well in their work about reduced-order models designed to diagnose flow-induced Nonsynchronous vibration (NSV) in turbomachines.The paper contains welldescribed techniques dedicated to obtain POD modes and visualize POD's harmonics as flow field in the blade channel.Also the solution predicted by ROM correlates very well to CFD solution in "blind" test of the two-dimensional Navier-Stokes model of the C1 single stage axial compressor.The authors concluded that POD methods could be successfully used to reduce computational cost of traditional CFD models that can predict NSV correctly.Danaila et al. [4] did POD of the centrifugal pump CFD-solution to investigate rotorstator interaction and used method as concurrent to traditional Fouriertransformation to obtain frequency characteristics of the flow-field.To sum up, they reproduced ONERA results obtained by Rochuon et al. [2].Danaila's group explained numerical effects arising due to interfaces in the fluid domain and pointed reliability of the POD for saving information about unsteady complex flows in the numerical models.Fossati et al. [5] developed the reduced-order model of axial turbomachine.They showed the same advantages of POD as Danaila's group.The numerous publications dedicated to POD analysis shows that this method has various applications in industry and applied science despite its statistical nature.In this paper, the certain problems of relation between physically reasonable simulations and statistical analysis are also investigated.From the engineering point of view, the POD method could be very reliable and useful as a replacement to FFT at certain points or planes.Similar to Rochuon's [2] work we also use Fourier-transform results to prove POD robustness and reliability, but they were obtained by other researchers for the same case.

Research objectives
The aim of this research effort consists of three subtasks.These are:  test POD as potential technique for simplification of the large amount of CFD data by relatively small set of eigenvalues;  compare POD-analysis of the test pump with FFT-analysis in case of amount of harmonics need for pressure-signal reconstruction and obtained information;  determine modes that carry information about dangerous regions of the flow in the sense of high probability of resonance between structure and flow.The first one deals specifically with POD algorithm development and validation, which was carried using various tests.The second one allows us to prove POD advantages relate to Fourier-transform techniques and prove the concept using sophisticated industrial-related problems, taking a simple pump model as an instance.The complexity and specialization of the developed techniques increases significantly at the last stage when all developed software stack is applied to industrial case.Concerning to the third objective we use FFT transformation of the time-dependent POD coefficients to identify frequencies of the modes and then compare them to the resonance frequencies presented in the documentation for the chosen compressor.

The solver
The pisoCentralDyMFoam solver is an extension of the newly developed pisoCentralFoam [6] with an option of the mesh motion handling in the numerical model.The OpenFOAM dynamicFvMesh abstract class generalizes methods to handle mesh transformation at the level of points and cells.In our work we only use mesh revolution around one axis of rotation but the mesh motion could be more complex.In the connection with turbomachinery problems such dynamic mesh motion solvers usually applicable to complex transient flows.The possible applications could include rotor precession modelling and tip clearance variation due to thermo-mechanical expansion very similar to Amirante's et al. work [7].The dynamic mesh modelling deals with the sliding interface approach which is opposite to frozen rotor when stationary problem is simulated.The main advantage of the hybrid Kurganov-Tadmor solver is in the possibility of simulation compressible flows in wide range of Mach numbers (0<M<6 [6]).The numerous amount of turbomachines equipped with multistage compressors have different operating conditions varying from relatively low velocities at the inlet to transonic and supersonic at the outlet of the last compressor stage.Thus, using one universal open-source solver makes typical engineering workflow even simpler and closer to commercial computational codes.

The test case
The ERCOFTAC centrifugal pump is a basic well-described test case for CFD-code validation.Petit et al. [8] carried out successful testing of their codes specially developed for turbomachinery applications.Combès [9] carried out the initial publication of the described test case at a Turbomachinery Flow Prediction ERCOFTAC Workshop in 1999.The experimental results were obtained at the test rig built by M. Ubaldi et al. [10].We used two-dimensional version of the initial case due to its simplicity and possibility to include as basic solver validation test.The geometry of the model is identical to description given in [9,10] and all dimensions are shown on the Figure 1.The test-case mesh was generated using ICEM-HEXA, consists of 94 000 cells, with average Y+ value about 35.For sliding mesh interface, we used Arbitrary Mesh Interface (AMI).Also we used slightly different boundary conditions: the average static pressure at the inlet and the mass flow rate at the outlet.Petit and Nilsson [11] used constant velocity at the inlet, which is equivalent to mass flow, and constant static pressure at the outlet.The same as in [8,11], the tip clearance which is equal to 1% (0.4 mm) of the blade span was not modelled in our simplified 2D verification simulation.Also we used slightly different boundary conditions: the average static pressure at the inlet and the mass flow at the outlet.Petit and Nilsson [11] used constant velocity at the inlet, which is equivalent to mass flow, and constant static pressure at the outlet.Based on the experimental data [10] we supposed that the inlet turbulent intensity is assumed to be 5% with a viscosity ratio of 10.The static pressure at the inlet is constant and equal to 101325 Pa (not given explicitly by [10]).Subsequently, the outlet boundary conditions for all variables except velocity are given a zero gradient condition, but for velocity are equal to flow rate.Mass flow was calculated using M. Ubaldi [10] pump performance data by the simple formula: The simulations were conducted using variable time-step adjusted by the maximum Courant number (equal to 0.5).The typical value for the time-step during simulation was about 5e-6 s.After several rotations the final residuals was equal to values: RMS U x = 1.56629e-11,RMS U y = 1.36025e-11,RMS p = 7.74527e-10, Continuity errors: sum local = 4.22966e-12, global = 1.10594e-12, cumulative = -3.41295e-09,RMS ω = 7.6337e-12, RMS k = 7.98248e-12.

Validation results
The used solver with k -ω SST turbulence model demonstrates good corresponding to both experimental and simulation results obtained by Petit and Nilsson using k − ε turbulence model.The calculated velocity profiles are presented in the Figures 1, 2. The difference between solvers could be explained by different types of solutions (2D in our and 3D in [11]) including numerical schemes effects and different turbulence treatment.

Application of the Proper Orthogonal Decomposition for turbomachinery cases
The vast majority of papers dealing with POD applications in turbomachinery (for instance [2], [3], [12]) uses traditional derivation of the basic equations as it was initially described by Lumley and then developed by him with co-authors in latest work [13].The ground work [13] also discuss in detail the physical meaning of the evaluated modes of the experimentally observed or simulated flow quantities.Using POD we are seeking for approximation: where u(x,t) usually implies function (i.e.velocity) in certain spatial location, represents spatial modes and ) ( t a i its time-dependent coefficients.Since POD is purely statistical method and do not account dimensionality (and physics of the phenomena), using it for problems with mesh motion or geometrical deformations is quite problematic and needs to be carefully considered.Present well-described works about POD applications in the turbomachinery do not include elaboration of the rotating regions (Fjällman J. et al. [12]), process only section in the stator part (ONERA work [2]) or dealing with "frozen rotor" approach ([3], [4], [14]).For nonstationary, full-dimensional "sliding mesh" approach, the proper application of the POD stays questionable.Closely connected research works about deforming meshes use two different methods to handle mesh moving or deformation.Anttonen et al.
[17] interpolated time-dependent computational fluid domain on the separate static mesh used to evaluate POD modes.In our work, we evaluate flow modes in the indexed space (see Figure 3b) of the elementary fluid volumes due to simplicity of the implementation in the CFD code and reliability when cylindrical coordinate system, which is typical for turbomachinery, is used.Thus, the equation (1) for the velocity component could be rewritten as: where n -volume number in the CFD mesh hierarchy.Consequently in the each moment of time (snapshot), a i (t N ) represents the contribution of the whole mode the in elementary volume, -time-invariant spatial contribution to the kinetic energy of the incompressible flow.According to [13] for compressible flow we should introduce non-dimensionalized sum of the flow parameters.Across the whole snapshot collection the relative position of the different volumes in the index list does not change (mesh transformation handles only on the point level which is the lowest in the hierarchy), so main limitation (points does not change its relative position in index space) of the POD is fulfilled [13].It needs to be pointed that POD -algorithm itself was tested on various simple cases including vibrating string with various initial and boundary conditions and Korteweg-de Vries equation using known soliton solution as reference.In the present work, we validate and demonstrate POD capabilities using only Ubaldi [10] test pump results.

POD analysis of the ERCOFTAC compressor
Relatively low computational cost of the 2D test case is good advantage for the developing of various algorithms to meet the technology needs.That is the main reason why Ubaldi [10] centrifugal pump was chosen as "proof-of-concept" case for demonstrations of the POD applications in the turbomachinery field.The studying of the flow phenomena using simple 2D representation of the ERCOFTAC compressor have an additional property such as simplicity of the visualization and analysis.The calculated lowest mode 0 for radial and tangential velocities is shown in Figure 4. Higher 6th mode is shown in the Figure 5.The characteristic of the spatial distribution of the spatial coefficients across the domain on the mode 0 (Figure 4) is very different compare to 6th mode where we can observe separation on the static and rotating regions.The wakes behind the impeller blades are almost uniform and independent from the relative position of the stator and rotor.We can say that mode 0 represents mean flow distributed across all computational domain and various turbulent effects in the wakes regions (for example, wake breaks when impeller blade passes near stator one) are captured by higher modes.The comparison between modes of the radial and tangential components also tells that tangential mode has larger absolute contribution to wakes region than radial (Figure 5).All these facts say that POD as statistical method could carry physically reasonable information about flow in compressor if we chose proper coordinate basis (using orthogonal coordinate for mode decomposition is potentially erroneous or needs much more modes and makes further analysis very hard).

POD versus FFT
In the original work about flow analysis, in the chosen test pump [11] the authors done FFT analysis using two probes at specific points and 20 frequencies for signal reconstruction.It can be shown that POD could catch all detected peak frequencies using fewer modes than frequencies in FFT, which has correspondence to the results of the ONERA group ([2]).In addition to the demonstrated visualization capabilities, the FFT analysis of the modal time-coefficients allows engineer to catch characteristic frequencies of the flow and reconstruct certain value of the quantity in the specified location.In the Figure 6 the value versus time, graphs were plotted to demonstrate time-dependent properties of the found POD modes.It can be seen, that mode 0 has constant value while higher modes have not.However, it is only due to the scale of the signal versus time graph: the FFT analysis of the mode 0 time coefficient shows that mode 0 has small peaks in the specific frequencies but with different amplitude from the other higher modes (Figure 7)., 7), on the phase-space diagram, such coefficients form concentric circles and the distribution of spatial mode coefficients is similar except little relative rotation.Spatial concentration of the non-zero modal coefficients for mode 3-4 is located in the rotor-stator clearance (very similar to mode 5 in Figure 5).These modes catch frequencies 400, 800, 1200 Hz, which is the frequency of the diffuser blades.Both radial and tangential components have these frequencies (Figure 7).There is no significant difference between velocity and pressure POD (see Figure 8), if we have an intention to describe frequencies of the flow.The pressure scalar field is much simpler to decompose and analyze.Thus for complex 3-dimensional flows we can evaluate only pressure field modes if we do not need to capture and store separation and wake -effects in the velocity vector field.

POD analysis of the industrial centrifugal pump
The present research work was conducted using model of the small-scale jet engine compressor developed in the early 90's by Andreas Funke and Klaus Wittig (Figure 9).The whole engine is also used as test case for open-source finite-element software CalculiX [18].The aluminium impeller of the selected compressor is designed for operation with maximum circumferential tip velocity 600 m/s, which is achieved for the diameter of 44 mm at 130000 revolutions per minute.Figure 9b shows the simulation domain including static and rotating regions of the mesh.
For the given case, initial and boundary conditions were set up as described:  undisturbed fluid in the domain at the starting point;  mass flow and static temperature at the inlet;  the fixed static pressure at the outlet 400 kPa and zero gradient for all others variables;  slip condition for velocity at the walls, zero gradient for pressure, adiabatic for temperature and wall-functions for turbulent parameters (so LES actually transforms to DES);  various models for turbulence treatment: k-ω SST and LES with Smagorinsky sub-scale model and cube-root filter.

CFD simulation results
Due to limited space we would not show various figures with simulated flow quantities and elaborate only integral results to control quality of the whole centrifugal pump model.Mesh convergence tests were performed to prove efficiency of the new hybrid density-based central difference solver in relation with standard pressure-based open-source sonicDyMFoam solver.Various meshes were generated using the same geometry and block structure as it was shown in Figure 9b except that the edge discretization variated from case to case.
For the hybrid solver equipped with 1.5 mill.mesh the calculated pressure ratio was 3.17 at N=108 000 rev/min which is close to the documentation value 3.30 [18] for the given mass flow 4,2 kg/s at N=110000 rev/min.In such conditions the error was -4%.The evaluated shaft power was 73.4 kW.
The convergence dynamic for the pressure-based solver is shown on the Table 1.We conclude that it has less performance than density-based central-difference solver using only 1.5 mil mesh and the last is significantly closer to the experimental results.

POD modes and flow dynamics
For the investigated centrifugal pump the simulation results are presented excluding repeating neighboring modes.The FFT-transform of the time coefficients of the odd modes from 0th to 9th is presented in the Figure 10.The visualization of the selected odd modes is presented in the Figures 11-13.The same as in simple two-dimensional validation case, three-dimensional compressor has main 0th mode with significantly lower amplitudes of the timecoefficients.Figure 12 demonstrates distribution of the spatial coefficients across the domain in the impeller and diffuser parts consequently.For the main 0th mode at the Figure 12a we can see that the coefficients are quite uniform and do not depend on what part of the domain, rotating or stating, they belong.Specific wake pressure fluctuations behind vanes were captured by third mode (and its neighboring mode 4th, which is not showed here).Combining it with the FFT results of the third mode (Figure 10) we can clearly say that characteristic frequency (12970 Hz peak) of this wake could potentially coincide with resonance frequencies 12000 Hz and 13700 Hz ([19]) in case of changing operating conditions.The cross-section data (Figures 12, 13

Conclusion
The Proper Orthogonal Decomposition could be reliable and effective tool for analysis of the transient properties of the flow.Evaluated modes help to locate specific regions of the simulation domain, which are connected to its own frequencies.
With the help of the FFT these frequencies could be easily identified.
The application of the method to the industrial case demonstrated successful identification of the flow instabilities with certain frequency near construction resonance.It is the main advantage of the POD relate to FFT analysis in which researcher have to specify the region of interest before processing it.Approaching simulation of the large-scale problems, POD can save more information for further analysis and be an effective basis for ROM.

Fig. 4 .
Fig. 4. Mode 0 for radial (left) and tangential (right) components of the flow velocity in the ERCOFTAC centrifugal pump (cylindrical coordinate system)

Figure 5
Figure 5 Mode 5 for radial (left) and tangential (right) components of the flow velocity in the ERCOFTAC centrifugal pump (cylindrical coordinate system)

Fig. 6 .
Fig. 6.The time coefficient for the first six modes of the radial flow component, obtained during one revolution of the impeller

Fig. 8 .
Fig. 8.The FFT-transform of the time coefficients of the first ten odd modes of the pressure field.The colored vertical lines are the impeller (blue), diffuser (green) frequencies, and its harmonics.The light pink vertical lines are the rotor harmonics Fig. 9. a) The compressor geometry [18]; b) simulation domain (1.5e6 cells mesh size)

Fig. 10 .
Fig. 10.The FFT-transform results of the time coefficient of the first odd modes from 0th to ninth.The colored vertical lines are the impeller (blue), diffuser (green) frequencies, and its harmonics.The light pink vertical lines are the rotor harmonics ) obtained at approximately half of the channel height in the outlet of the impeller demonstrates the same results as three-dimensional data in the Figure11.The distribution of the spatial coefficients in the cross-section similar to pressure distribution at the 0th mode and various non-stationary effects are localized by higher modes.All these effects are concentrated in the clearance regions and depend of the time-coefficient amplitude and spatial coefficient contribute more or less at certain moments of time(Figures 12b, 13).
Fig. 12. a) The original pressure field (left) and mode 0 (right) of the high-speed centrifugal pump in the diffuser cross-section; b) Pressure field mode 1 (left) and mode 3 (right) of the high-speed centrifugal pump in the diffuser cross-section