Abstract
Unsteady flow around an oscillating plate cascade and that through a single compressor rotor subject to vibration have been computationally studied, aimed at examining the predictive ability of two low fidelity frequency methods compared with a high fidelity time-domain solution method for aeroelasticity. The computational solutions demonstrate the capabilities of the frequency domain methods compared with the nonlinear time-domain solution method in capturing small perturbations in the unsteady flow. They also show the great advantage of significant CPU time saving by the frequency methods over the nonlinear time method. Comparisons of two different frequency methods, nonlinear harmonic and phase solution method, show that these methods can produce different results due to the differences in numeric and physical conditioning. The results obtained using phase solutions method are in better agreement with the nonlinear time-domain solution. This is because the same numeric and physical conditioning are used in both the nonlinear time-domain method and phase solution frequency domain method.
1. Introduction
The blading aerodynamic design has long been employing steady flow methods as they are highly efficient and robust [1–3]. The steady solution methods can be simply automated for design applications using optimization technique or inverse approaches [4, 5]. Even at a detailed design stage, the solution to the steady flow equations rather than the unsteady flow equations is sought widely for blading aerodynamics predictions. This is because in many aerodynamics design applications the unsteadiness in the flow field particularly at its design operating point is usually small. As a consequence, the time-averaged flow solution is not greatly influenced by the unsteadiness perturbations. In fact, the current aerodynamic blade and airfoil designs have achieved high aerodynamic performance by using the steady flow model. The solution of unsteady flow equations is much more costly, requiring substantial computer resources [6–8].
There are many cases in which the unsteady forces have significant effects that cannot be ignored. For example, the effects on blade vibration, noise generation, fatigue, or failures of blades are very important. As modern multirow turbomachinery components are typically designed with high loading and more compact structural configurations, those unsteady effects are intensified in them. As aeroloading is increased, the blade mechanical integrity, aerodynamic noise, and vibration stress levels need to be carefully examined as the blades will be more vulnerable to flow-induced vibration problems. Recently a great advancement has been made in computing power and numerical methods for unsteady flows. However, within a multidisciplinary design environment there is always a need for developing efficient and fast prediction methods. This is because it is impractical to perform an aeromechanic related design optimization based upon costly unsteady flow solution methods, considering that a design optimization is carried out in an iterative process. As a result, significant research has been devoted in recent years to developing efficient numerical approaches capable of capturing the major unsteady features which are relevant to the engineering problems of interest while reducing the solution time to an acceptable level for use in routine design. One of the earliest methods is the time-linearized harmonic frequency domain methods which have been widely used for turbomachinery aeromechanical applications [9, 10]. In these methods the unsteady flow equation is considered as one steady equation and one perturbation equation. The unsteady perturbation of the flow solution is represented in a Fourier series. The validity of these methods is limited by the linear assumption which often exhibit solution divergence behaviours for highly nonlinear flows.
Recently advanced frequency domain such as nonlinear harmonic and phase solution methodologies have been developed which on the one hand is efficient for an unsteady flow solution for aeromechanics and on the other hand can be used for blading design optimization without much extra effort. One recent method, phase solution method [11–13], provides a straightforward simple method of modeling unsteady perturbations. In these methods the unsteady flow equations with a single periodic unsteadiness are solved at two or three distinctive phases of a period of unsteadiness. Using this approach the same computational efficiency as a conventional time-linearized method can be obtained. While the other time-linearized harmonic methods express the whole flow solution in a Fourier series, this method is based on casting the unsteady flow equations into a set of steady-like equations at a series of phases of a period of unsteadiness. By using this method the same steady flow solution method used for aerodynamic design can be used for aeromechanical design. Therefore the aeroelasticity of a blade can be optimized in a design optimization process to check both the aerodynamics and aeromechanics simultaneously. Recently, Valero et al. [14] developed a concurrent design optimization method on a phase solution method developed by Rahmati et al. [15] in conjunction with a commercial FEA solver for blade structural dynamics.
Validations and verifications are crucial in establishing our confidence in the numerical methods and thus the conclusions drawn upon those numerical results. It is not the aim of this study to apply the frequency domain method for very complex configurations. The application of various frequency domain methods for complex multirow turbomachines has been already reported [16, 17]. This paper focuses on two relatively simple cases, flow around a 3D flat plate oscillating cascade and through a compressor rotor subject to vibration, to demonstrate the validity and compare the effectiveness of three computational methods. The aim is to highlight the importance of the numerical diffusions and physical conditioning on the computational results. Although in the phase solutions method the same numeric and physical conditioning as the nonlinear time-domain method are used, this is not the case in the nonlinear harmonic method.
A brief description of the governing equations and the time and frequency domain methods is given in the next section; this is followed by the prediction of flow over an oscillating plate and flow through a single compressor rotor using various methods. These cases are provided to verify the capability of various numerical methods based on frequency domain method in predicting unsteady flow in comparison with time-domain. It is also used to show that it is important to develop both frequency and time-domain solution methods using the same numeric and physical conditioning. By doing so, when the high fidelity nonlinear time-domain solution is compared with the low fidelity frequency domain, one can consistently find the difference in the modeling fidelity without interference from differences in numeric.
2. Flow Modeling Formulation
2.1. Flow Governing Equations
The unsteady Reynolds averaged Navier-Stokes equations in a cylindrical coordinate system integral form can be described by the following form:, , and are the convective flux vectors and , , and are the viscous flux vectors while is the source term due to rotating effects.
An absolute frame of reference is used for the solution of the flow equations. The standard one-equation Spalart-Allmaras model is implemented to calculate the eddy viscosity for the turbulence model. The equations are discretized in the finite-volume form with a cell-centred variable storage. The structured -mesh is used in all calculations.
The following form can be used to describe the semidiscrete integral form of the equations:In the above equation the summation is over the boundary faces of a mesh cell. The normal projected areas in the three coordinate directions are , , and , respectively. In this formulation , , and are the flux terms, including the moving grid terms, the nonlinear convective fluxes, and viscous terms in the momentum and energy equations. The right-hand side of (2) can be considered together as a residual term. So, the equation can be simply expressed asA 2nd-order cell-centre based finite-volume scheme with a blend of 2nd- and 4th-order numerical dissipations is used for spatial discretization (see Jameson et al. [18]). The time-marching solution is based on a 4-step Runge-Kutta integration in the pseudo time. The temporal change of flow variables at each fractional time step () isThe multigrid and local time stepping techniques are used to accelerate the time-marching solutions [15].
2.2. Time-Domain Direct Method for Unsteady Flow
The time-domain solution is obtained using the basic time-marching solver. It is important to note that a time-consistent multigrid is used for accelerating time-integration for an unsteady flow. A simple implementation of the time-consistent multigrid is to use just two grids. In this implementation, the spatial resolution is governed by the fine mesh while the coarse block is for enlarging the time step in a time-consistent manner [19]. Several levels of intermediate-meshes between the basic fine mesh and the coarse mesh are introduced. If is the levels of intermediate-meshes, the general multigrid formulation for the temporal change of flow variables for a fixed mesh can be written asThe subscripts “” denote the fine mesh while “” represent the intermediate-mesh (level ) and the coarse mesh is represented by “”. is the net flux; is the cell or block volume. and are the allowable time step lengths on the fine mesh and the intermediate-mesh of level (), respectively.
In order to maintain a uniform global time step length for unsteady calculations, the time step length for the coarse block needs to satisfy the time-consistence condition:
2.3. Harmonic Method
For simplicity, the flow governing equation for unsteady flow is rewritten asSince the grid movement is usually prescribed and these terms can be calculated prior to the main flow field solution, in the above equation, the unsteady terms corresponding to the temporal change of the mesh volume are moved to the right-hand side. In a frequency domain method the unsteady flow is assumed to be composed of a steady and a fluctuating part:where is the steady part of the flow variables while is the fluctuating part. The fluctuating part is harmonic in time for many flows of engineering interest. So the unsteady perturbation flow can be expressed as a Fourier series. Since a linear assumption has been made, the behaviour of each Fourier component can be analysed individually and then summed together to form the total solution. Therefore, for a single periodic disturbance, the fluctuating part can be represented as where is the vector of complex amplitudes of perturbations while is the complex conjugate of .
Substituting the relationships (8) into (9) yields the following equation:
Substituting the relationships (10) into (7) and collecting the zeroth- and first-order terms lead to two equations: one an equation for the steady and the other one to a linearized perturbation equation. The steady equation is solved by ignoring the time-dependent term and has the following form:A similar equation can be found for the linearized perturbation equation [20]. So basically in this method, first, a steady flow solution is obtained by solving the steady Navier-Stokes equation. Then, for a given frequency and interblade phase angles, the coefficients of the time-linearized equations are formed from the steady flow solution. Finally the time-linearized perturbation equations are solved. Thus, a time-dependent unsteady problem in the time-domain is effectively transformed to solving two equations in the frequency domain. This linear method can be extended to nonlinear method based on the work of Adamczyk [21] who showed that time-averaging the Navier-Stokes equations resulted in the inclusion of the effect of periodic perturbations on the mean flow through stress terms. The nonlinear harmonic method is similar to RANS equations in which the effect of turbulence is included by introducing the Reynolds stresses. Extra closure models are required to work out these deterministic stress terms similar to turbulence modeling for Reynolds stress terms. In this approach, a significant modification is that a time-averaged value instead of a steady flow value is used as the basis for the harmonic perturbations. So, it is assumed that the flow field is composed of a time-averaged flow with a small perturbation unsteady flow:
Here is the vector of the time-averaged conservative variables and is the vector of perturbation conservative variables. By substituting expressions (12) into (7) and taking a time-averaging, the resultant time-averaged Navier-Stokes equations are obtained. In comparison with the steady form used in linear methods the time-averaging produces extra terms because of the nonlinearity of the flow governing equations. The extra stress terms due to the velocity fluctuations are the result of nonlinearity of the flow governing equations. One assumption made here is that the unsteady perturbation is dominated by first-order terms so the first-order harmonic perturbation equation has the same form as the unsteady perturbation equation in the time-linearized method. Obviously if the time-averaged flow is the same as the steady flow, the above first harmonic perturbation equation reduces to the conventional time-linearized perturbation equation. The time-averaging procedure produces extra unsteady stress terms in the time-averaged equations which are evaluated from unsteady perturbations. While the unsteady perturbations are obtained from solving the harmonic perturbation equation, the coefficients of disturbance equations resulted from the solution of time-averaged equation and this interaction is achieved through a strong coupling procedure. More details regarding the governing equations can be found in [20].
2.4. Phase-Shift Solution Method
In this formulation similar to nonlinear harmonic method, by considering only one harmonic, the flow solution is given byBy substituting (13) into (7) the following equation is obtained:
At three different temporal phases (), (13) can be written as follows:
In (15a)–(15c) the three unknowns , , and are expressed in terms of the flow solutions at the three phases. Substituting these unknowns in terms of flow solutions at three phases into (7) yields the following equations:These equations are very similar to the steady equations. The difference is that there are three sets of equations instead of only one and there is an extra source term in each equation. The equations are instantaneously solved in a similar way to that of the steady RANS equations. The source terms are simply calculated using the flow solution at the previous iteration. One advantage of this phase solution method in comparison with the conventional harmonic solution method or linear harmonic method is that the nonlinearity is automatically retained in the nonlinear convection terms in the discretized equations. This feature leads to much improved solution convergence for harmonic disturbances, particularly for cases with flow separation. More details regarding this method can be found in [16].
2.5. Boundary Conditions
Similar boundary conditions are used in all methods. The slip wall condition is applied and the log-law is applied on solid blade and end wall surfaces to determine the surface shear stress. For inlet and exit to the computational domain boundaries, the same reflective treatment as those for the steady flow solutions can be used by specifying inlet stagnation pressure, stagnation temperature, inlet flow angle, and exit static pressure. Alternatively, the local characteristic disturbance-based nonreflective approach by Giles [22] can be used for inlet and outlet boundary treatments. The wall and far field inlet/exit conditions are applied in exactly the same way for both the time-domain and frequency domain solutions.
For the periodic boundaries, the direct periodic (repeating) condition is applied for the time-domain method. Similarly, a direct periodic condition is applied to a steady flow solution. For the frequency domain solution method using a single passage domain the harmonic components are phase-shifted between the upper and lower periodic boundaries by a given interblade phase angle, (the phase leads to the upper blade vibration relative to its lower neighbouring one).
3. Flow over an Oscillating Plate
A three-dimensional flat plate cascade test case is used for the comparison of the frequency domain and time-domain method.
The geometric parameters and flow conditions are shown in Figure 1 and Table 1. The blades are oscillating in a three-dimensional mode. Each two-dimensional section is subject to torsion mode around its leading edge with the torsion amplitude varying linearly along span from at the lower end to at the upper end.
The present calculations were carried out using both the frequency domain and the time-domain method with a mesh density of in the streamwise, pitchwise, and radial directions, respectively.
Calculated unsteady pressure jumps at each two-dimensional section are presented in the form of :where is the first harmonic pressure jump across the blade; is the torsion amplitude at the tip in radian. is a complex number with its real and imaginary parts, which can be used to work out amplitude and phase angle. The calculations were performed at the spanwise position of various normalized distance measured from the hub section ().
Figures 2–5 show the unsteady pressure on the lower and upper plate. The phase solution analysis results are virtually identical to the time-domain analysis results. Similar results are obtained for the linear and nonlinear harmonic method. Figures 6 and 7 show the chordwise distributions of the real and imaginary parts of the unsteady pressure jump at six spanwise sections for the interblade phase angle of 180 and 0 deg. The numerical results are compared with the analytical results for this case provided by Namba [23, 24] and He and Denton [25]. The linear and nonlinear harmonic results are very similar. This is something within expectation considering that the same numeric is used in these two different fidelity modes and no extra unsteadiness is available in the models. However, compared to nonlinear time-domain method, the phase solution method gives better agreement in particular for interblade phase angle of 180 degrees. It is important to note that in the phase solution formulation no extra numerical method is used for calculations of flow perturbations. So, basically, both phase solution and time-domain solution methods use the same numeric and physical conditioning. So, a comparison between the phase solution method and the time-domain method shows the difference in the modeling fidelity, which is none here, without interference from differences in numeric.
(a) Upper plate
(b) Lower plate
(a) Upper plate
(b) Lower plate
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
A very small discrepancy exists between the calculation and analytical data in all methods around the inlet due to reflection effects at inlet. A typical time-domain solver run for this case takes about 37 hours on a single CPU using a 2 GHz Intel Xeon CPU. The CPU time saved using the frequency domain solution method is quite significant as the solution in all methods takes about thirty to forty minutes using the same computer.
4. DLR Compressor Rotor
Three-dimensional unsteady flow through a transonic compressor from the German Aerospace Center (DLR, Koln, Germany) subject to blade vibration is computed to evaluate the capabilities of the frequency domain and time-domain methods. The rotor rotational speed is 20,260 RPM and it gives a stage pressure ratio of 1.59 at a mass flow rate of 17.15 kg/s. The compressor has 28 rotor blades and 31 stator blades. The detailed information of the compressor and test rig can be found in [3]. In previous studies [9, 11] multistage effects on the aerodamping of this DLR compressor have been investigated using both frequency domain and time-domain methods. In this study, unsteady flow through the rotor row has been computationally studied. The aim is to examine the predictive ability of two low fidelity frequency methods for aeroelasticity when the only source of unsteadiness is the rotor vibration. To make sure that no other source of unsteadiness (e.g., reflections from the downstream outlet) influences the unsteady flow through the compressor, the outlet in the computational domain is placed far away from the trailing edge. Figure 8 shows the 3D computational model of the rotor while Figure 9 shows the blade-to-blade view of the computational mesh at the mid-span section used in current computation. A total mesh point of 162500 is used in a single row.
In the modal analysis to get more realistic mode shape and fundamental vibration frequencies, the chosen material for the rotor blade has a density of 4428.8 kg/m3, Young’s modulus of 1251011 Pa, and Poisson ratio of 0.27. These properties are close to the property of titanium at the ambient temperature.
The unsteady flow through the rotor blade at various Nodal Diameters (ND) at one modal frequency is investigated. The ND of 4 with IBPA of 51.428° at mode 1 vibration frequency of 1350 Hz is examined in this paper. The nonlinear time-domain method is used for validations of two frequency domain methods: nonlinear harmonic and phase solution method. A typical time-domain solver run for this case takes about 3 days on a single CPU with a 2 GHz Intel Xeon CPU. The CPU time required for the frequency domain solution is only 2 hrs using the same computer.
Figure 10 shows the time-averaged pressure values at the mid-span of the blade calculated using time-domain and frequency domain method. Figure 11 shows the unsteady part of pressure values at the mid-span of the rotor. No noticeable difference can be observed in those pressure contours. Detailed comparisons of time-averaged pressure, 1st harmonic pressure amplitude, and phase angle are shown in Figures 12–14. Compared to the time-domain method both nonlinear harmonic and phase solution method predict the same time-averaged values for pressure at the mid-span. The comparisons of the 1st harmonic pressure amplitudes and phases show that the phase solution method predicts the same pressure amplitude and phase at the mid-span as those of the time-domain one. The results obtained using the nonlinear harmonic method, however, are slightly different. This method predicts lower values for the 1st harmonic pressure amplitudes at both suction side and pressure side and a higher value for phase at pressure side. As the only differences between these two frequency domain methods are the formulation and numeric used for predicting the unsteady perturbation, and there is no other source of unsteadiness in the flow, it can be concluded that these differences are because of the numerical diffusion in this model. The phase solution method predicts the same results as those of the time-domain method because in this method the same numeric and physical conditioning are used and no extra formulation is used to take into account the unsteady perturbation. Indeed, in this method, the same equations are solved at three phases. Since there is only one unsteady harmonic in the model, one expects that these two models predict the same results. This case highlights the importance of using the same numeric and physical conditioning for both low fidelity and high fidelity solution methods. When the high fidelity nonlinear time-domain solution is compared with the low fidelity phase solution method, one can consistently find the difference in the modeling fidelity without interference from differences in numeric. However, the same cannot be said for the nonlinear harmonic method due to the difference in numeric and the subsequent differences in the numerical diffusions compared with the time-domain method.
(a) Frequency domain
(b) Time-domain
(a) Frequency domain
(b) Time-domain
5. Conclusion
A linear oscillating plate cascade and a transonic compressor rotor are investigated computationally to assess the validity, merits, and limitations of two low fidelity frequency domains and a high fidelity time-domain solution method. The frequency domain methods can produce one or two orders of magnitude saving in computational effort. On practical relevance and implications, this low computational cost makes the frequency domain method a valuable tool in an aeroelasticity design environment in which a design optimization is carried out in an iterative process. Comparisons of the nonlinear harmonic and phase solution results demonstrated good agreement with the time-domain method over a wide range of spanwise positions. The results using phase solution method are virtually identical to the time-domain analysis results. The outcomes of the present work suggest that in aeroelasticity applications, in which the perturbation magnitudes are usually very small compared to the mean flow, the solution of a nonlinear phase solution method can be used for analysing unsteady flow. By comparing the time-domain solution with the phase solution method, one can consistently find the difference in the modeling fidelity without interference from differences in numeric. This is because in the phase solution method the same numeric and physical conditioning as the high fidelity time-domain are used.
Nomenclature
| : | Torsion amplitude at the lower end |
| : | Torsion amplitude at the upper end |
| : | Stagger angle |
| : | First harmonic pressure jump, Pa |
| : | Complex pressure jump coefficient |
| : | Real part of pressure jump coefficient |
| : | Imaginary part of pressure jump coefficient |
| : | Vibration frequency, Hz |
| : | Radial direction, m |
| : | Spanwise position |
| : | Axial direction, m |
| : | Spanwise direction, m |
| : | Chordwise position |
| : | Phase angle, deg |
| : | Density, kg/m3 |
| : | Steady flow component |
| : | Fluctuating flow component |
| : | Time-averaged flow component |
| : | Unsteady flow variable. |
Disclosure
Present address of M. T. Rahmati is Department of Mechanical, Aerospace, and Civil Engineering, Brunel University London, Uxbridge UB8 3PH, UK.
Competing Interests
The author declares having no competing interests.
Acknowledgments
First and foremost thanks are due to Professor L. He from Oxford University for Providing Namba’s analytical data and valuable discussion and help throughout the course of this research at Oxford. The author would like to thank Siemens Industrial Turbomachinery Ltd., Lincoln, for providing funding for this research at Oxford University. Thanks are also due to Professor M. Namba from Sojo University, Japan, for valuable discussion regarding his analytical results for the three-dimensional flat plate cascade used in this study.