Category UNSTEADY AERODYNAMICS, AEROACOUSTICS AND AEROELASTICITY OF TURBOMACHINES

The steady-state solution is found by implicitly time-stepping the unsteady ft>w equations (Eqn. 1) until the steady-state residual has been reduced sev­eral orders of magnitude. Local time-stepping is used to increase the rate of convergence. The traditional method [2, 3] of solving the time-linearised ft>w equations (Eqn. 9) is to add pseudo-time derivatives to the left side of the equations. The equations then have a similar form to the unsteady ft>w equa­tions and hence the same method can be employed to solve them. Instead of the traditional pseudo-time-stepping, the current scheme uses GMRES with preconditioning to solve the linearised flow equations.

The Generalized Minimum RESidual (GMRES) [13] can be used to solve large, sparse and non-symmetric linear systems. GMRES belongs to the class of Krylov based iterative methods. Consider the linear system

Ax = b (12)

where A is a square non-singular n x n complex matrix, and b is a complex vector of length n. Let xo Є Cn be an initial guess for this linear system and ro = b — Ax0 be its corresponding residual. The GMRES algorithm builds an approximation of the solution of Equation 12 under the form

xi = xo + Vi y (13)

where Vi is an orthonormal basis for the Krylov space of dimension l

where y belongs to C1. The vector y is determined so that 2-norm of the residual ||r ||2 is minimal over Ki. The 2-norm is defined as

(15)

For more information about GMRES see Saad [13].

The convergence rate of Krylov subspace methods is highly dependent on the condition number of the matrix A. The smaller the condition number, the better the convergence. Preconditioning can be used to reduce the condition number of the linear system. The current scheme uses ILU(p) [13] to right precondition the linear system. The preconditioner is calculated by performing an incomplete Gaussian elimination. The method is incomplete because the non-zero pattern on the preconditioner is determined by the non-zero pattern of the original matrix and the level of fill p. For inviscid cases, a level of fill p =1 gave the best performance and for viscous cases p = 3. In order to construct the preconditioner, it is necessary to explicitly determine and store the governing matrix of the linearised fbw equation (Eqn. 9).

2. Test Cases

In order to verify that a computational method is accurately solving the given equations for real profiles where no analytical solution exists, it is nec­essary to compare results with previous numerical work. Most of the test cases shown are of Standard Configuration 10 mainly due to the quality, quantity and availability of previous numerical solutions. Unfortunately this is an Euler test case, but many features of the code (correct Euler formulation, correct far-field boundary treatment, and correct treatment of shocks) can be examined by this case. In order to verify the correct formulation of the linearised Navier-Stokes equations and the turbulence model, Standard Configuration 11 is examined.

The governing flow equation is transformed to the frequency domain by assuming that the prescribing motion of grid points is harmonically oscillating about their steady-state positions

x = x + 3?{xeJaJt} (4)

and the resulting unsteady cell-averaged perturbations are also harmonically oscillating about its steady-state values

U = U + 5R{UeJaJt}. (5)

The discretised time derivative of the fbw solution, R (Eqn. 3) is linearised by performing a first-order Taylor expansion about the steady-state flow values, steady-state grid position, and zero grid velocity.

Note that it is assumed that the time-derivative of the fbw solution at the steady-state condition,

R = R(U, x, 0) = 0. (8)

Substituting Eqs. 7, 4, 5 into Equation 6 and dropping the e]wt term gives the following result

AR(U, x + x, 0) + j AR(U, x, ш x) (9)

where

The right hand side of the linearised fbw equation can be easily determined by calculating the change in R from the steady-state due to perturbing the grid by the amplitude of the grid motion and also calculating the change in R due to setting the grid velocity equal to the amplitude of the maximum velocity of the grid motion. The partial derivative of R with respect to the ft>w solution required on the left hand side of the equation is equal to

The Jacobians required are calculated numerically and the Jacobian of the fhx calculation includes the action of MUSCL and the limiter. Note that some of the flix schemes (AUSMDV and the exact Riemann solver) used here are not strictly differentiable for some input ft>w states. However, for the cases to be presented it was possible to obtain accurate solutions using numerically calculated derivatives of these fluxes.

The two-dimensional Reynolds-averaged Navier-Stokes (RANS) equations with the Spalart & Allmaras one-equation turbulence model [4] are used to model the unsteady fbw. The ft>w is assumed to be fully turbulent and no transition model is included. The turbulence model was chosen by the author and previous researchers [2, 3] because of its ability to accurately model the ft>w at a reasonable computational effort [5]. Better turbulence models do exist, but the effort required to solve their equations is significantly greater and prohibits their use for design calculations. The unsteady ft>w equations for a moving control volume can be written

where U is the vector of conserved variables, # is the control volume, S is the control surface, H is the fhx of conserved variables at the control surface, t is time, and S is the vector of source terms due to the turbulence model.

The ft>w domain is discretised by a multi-block structured grid, and the ft>w equations are discretised using a cell-centered finite-volume scheme. The discretised ft>w equations can be written

d 4

-[U + = s#,

i=1

or

where U is the vector of the cell-averaged conserved variables, # is the volume of the cell, Si is the area of a cell interface between two cells, and Hi is the average flix through a cell interface. The inviscid flix is determined by re­constructing the ft>w states at both sides (typically called left and right) of the cell interface using MUSCL interpolation with a modified van Albada limiter

. The flix at the interface can be estimated from the left and right states us­ing one of the following upwind schemes: Roe’s approximate Riemann solver

, AUSMDV [9], an exact Riemann solver [10], and EFM [11]. The spatial derivatives required for the viscous fluxes are calculated by an application of the divergence theorem [12].

The boundary conditions applied at the wall, on the periodic boundaries, and at the far-field boundaries are the same as applied by the previous researchers [2, 3]. It should be stated that the exact (for uniform ft>w) non-refecting boundaries described by Giles [6] has been implemented at the far-field bound­aries. The implementation and verification of this boundary condition required significant effort.

Paul Petrie-Repar

Institute ofAeroelasticity, German Aerospace Center Bunsenstrasse 10, 37073 Gottingen, Germany.

Email: Paul. Petrie-Repar@dlr. de

Abstract A time-linearised Navier-Stokes flaw solver is presented. The turbulence model of Spalart & Allmaras has been included in the fOw model. The main feature of the current method is that GMRES with preconditioning is used to solve the linearised equations in a robust and efficient manner. The method is verified by the presentation of solutions of test cases from Standard Configuration 10 and 11. The use of various upwind schemes is also examined.

Keywords: linearised Navier-Stokes, Spalart and Allmaras turbulence model, GMRES, pre­

conditioning

1. Introduction

Blade flitter in turbomachinery is highly undesirable as it can lead to dra­matic blade failure. In order to design blading that will be free of flitter, the design engineer requires a tool that is capable of accurately predicting blade flitter. Unsteady fl»w simulation is an essential component of flitter analy­sis. For transonic fl>w or separated fl>w, it is necessary to consider viscous effects. Time-marching of the Navier-Stokes equations [1] is one method of simulating the unsteady flow, however this is computationally expensive and hence has a limited use in real design. An alternative to the time-marching method is the so called time-linearised method [2]. Here the unsteady fl)w is assumed to be small harmonic linear perturbations about the steady flow. This assumption allows the equations to be transformed to the frequency domain. The resulting linear equations require significantly less computational effort to solve (1 to 2 orders of magnitude faster than time-marching methods). Here, an efficient and robust method for solving the time-linearised Navier-Stokes equations is presented. The formulation of the linearised equations is similar to the methods presented by Clark & Hall [2] and Sbardella & Imregun [3].

437

K. C. Hall et al. (eds.),

Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines, 437-448. © 2006 Springer. Printed in the Netherlands.

The main difference between the current method and the previous methods is the use of GMRES with preconditioning instead of pseudo-time-stepping to solve the linearised equations. This was found to be a very efficient and robust method.

5.3.1 Unsteady Navier-Stokes response to harmonic motion. Pre­scribed harmonic motion Navier-Stokes simulations have been performed for a 3D wide chord fan. This fan is made up with 22 swept blades. The maxi­mum radius of the fan is about 0.9 m. A Navier-Stokes grid of moderate size has been built in order to run Spalart steady and unsteady computations. It is made up with 6 blocks, and its total number of nodes is 397044. The first grid layer thickness at the wall is about 5.e-06 m. A view of the grid and of its multi-block topology is given in the next figure.

A steady computation is initially per­formed for an upstream absolute Mach number of 0.5. The rotating speed of the compressor for this computation is 4066.4 RPM. 4000 iterations were run using the Spalart-Allmaras model at a CFL value of 5. The computation time was about 12 hours on a single Itanium2 900MHz pro­cessor. Here is shown the quadratic resid­ual convergence history of the conserva­tive variables for one block. Figure 6 Convergence history

The outlet boundary condition prescribes the value of the output pressure on the hub. For this computation, the mass fbw of the compressor is about 458 kg/s. A shock occurs near the tip, either on the suction side or on the pressure side. The maximum Mach number is about 1.5. Figure 7 shows the Mach contours on the suction and pressure sides.

Dynamic Navier-Stokes fluid – structure coupling. Navier-Stokes dy­namic fliid-structure coupling computa­tions have been performed for an ad­vanced wide-chord swept blade fan. The linear structural model was made of 10 modes. The numerical coupled scheme as described in a previous section has been used. The computation was performed on a 6 block grid gathering 372036 nodes. 320 physical iterations have been per­formed for a total simulated time of about 0.1 s. 100 dual time steps have been run at each physical time step, leading to a total computation CPU time of about 150 hours onItanium2. We present in Figs. 10 and 11 the time history of the generalized coordi­nates and that of the mechanical energy of the blade, for specific operating point and initial conditions.

The blade is clearly aeroelastically stable, which can be more precisely char­acterized through the processing of the generalized coordinates time histories, in order to extract frequencies and damping for this operating point.

4. Conclusion

A Navier-Stokes numerical tool has been developed for the computation of unsteady turbomachinery applications. An Arbitrary Lagrangian Eulerian formulation has been developed, and the dual time stepping acceleration tech­nique has been implemented in the 3D code. The basic scheme has also been modified in order to allow moving meshes computations. Static and dynamic fhid-structure coupling schemes have also been developed in the case of a modal structural model. Some results of the validation processes of the 2.5D

and 3D aeroelastic Navier-Stokes codes have been presented. An example of a dynamically coupled 3D Navier-Stokes Arid-structure computation has been given. We intend to go on with 3D developments in order to be able to per­form fully 3D Navier-Stokes unsteady turbomachinery computations for more complex configurations.

References

Batina, J. T. (1989). Unsteady Euler airfoil solutions using unstructured dynamics meshes. 27th Aerospace sciences meeting, AIAA Paper 89-0115.

Dugeai, A., Madec, A., and Sens, A. S. (2000). Numerical unsteady aerodynamics for tur­bomachinery aeroelasticity. In P., Ferrand and Aubert, S., editors, Proceedings of the 9th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines, pages 830-840. Lyon, PUG.

Girodroux-Lavigne, P., and Dugeai, A. (2003) Transonic aeroelastic computations using Navier – Stokes equations. International Forum on Aeroelasticity and Structural Dynamics, Amster­dam, June 4-6.

Jameson, A., Schmidt, W. , and Turkel, S. (1981). Numerical Solution of the Euler Equation by Finite Volume Methods using Runge-Kutta Time Stepping schemes. 14th Fluid and Plasma Dynamics Conference, Palo Alto (CA), USA, AIAA Paper 81-1259.

Leconte, P., David, F., Monnier, J.-C., Gilliot, A. (2001). Various measurement techniques in a blown-down wind-tunnel to assess the unsteady aeroelastic behavior of compressor blades 2001IFASD, June 05-07.

Lerat, A., Sides, J. and Daru, V. (1982) An Implicit Finite Volume Method for Solving the Euler Equations Lectures notes in Physics, vol 170, pp 343-349.

Spalart P., and Allmaras, S. (1992). One Equation Turbulence Model for Separated Turbulent Flows 30th Aerospace Science Meeting, AIAA Paper 92-0439, Reno (NV).

Vuillot, A.-M., Couailler, V., and Liamis N. (1993). 3D Turbomachinery Euler and Navier – Stokes Calculation with Multidomain Cell-Centerd Approach. AIAA/SAE/ASME/ASEE 29th Joint propulsion conference and exhibit, Monterey (CA), USA, AIAA Paper 93-2573.

We next present Navier-Stokes computations dealing with a high subsonic unsteady fbw over a PGRC profile cascade. The upstream Mach number is 0.9. The fbw angle of attack is 12 degrees. The computation was performed for a total pressure of 159881 Pa and a total temperature of 285.16 K. For this computation, a 3 domain HCH grid was designed for a single channel using a total of 16433 nodes. Continuity boundary conditions were used in the steady case at channel interfaces. We present in Fig. 3 the steady isomach lines map obtained by PIV technique and by the computation, over 2 channels.

The unsteady computations have been performed for 5 inter-blade phase angles. The dual timestep technique has been used with a CFL number of 4.

At each physical timestep, the first component of the fOw field must de­crease of two order. To describe a cycle we need 64 time-steps. When the inter-blade phase angle is different from zero, we need to run about 20 peri­ods ( 6 periods are enough for the zero inter-blade phase angle).We present the unsteady results, for a pitching motion of the central blade at a frequency of about 300 Hz, compared to experimental data. A rather good agreement with the experimental distribution can be noticed.

We present first some unsteady Navier-Stokes results obtained with the 2.5D solver for an isolated profile and the rectilinear cascade.

3.1 Isolated PFSU profile

The aerodynamic conditions for this 2D computation are an upstream Mach number value of 0.75, a total pressure of 1108121 Pa, and a total temperature of 299.8 Kelvin. The steady angle-of-attack of the ft>w is 3 degrees. The chord of the profile is 0.3 meter. The profile is moving in pitch at a frequency of 40Hz, with an amplitude of 0.25 degree. Navier-Stokes steady and un­steady computations were run using the turbulence model of Michel and that of Spalart-Allmaras. We used a 300×100 C-like mesh.

The unsteady run needed the computation of 4 periods of about 50000 iter­ations using a uniform time step at a maximum CFL value of 10. The previous sketch shows the mean and unsteady pressure coefficients on the profile (Fig. 1) and a comparison of steady velocity profiles for three axial positions (Fig. 2). Spalart-Allmaras results and available experimental data compare fairly well. Algebraic turbulent model seems to be insufficient to get a good description of the leading edge separated zone.

The first test campaign concerned the isolated PFSU profile for which a large amount of well-documented data has been obtained during the experi­ments located at Onera Modane in the S3 transonic wind tunnel. Steady and unsteady measurements have been performed for inlet Mach numbers of 0.5 to 0.75, with various static incidence angles (0 to 5 degrees) and pitching move­ments of the profile at a frequency of 40 Hertz. Beside this, the LDV (Laser Doppler Velocimetry) technique was used to assess the velocity profile in the vicinity of the profile, in order to get a proper description of the separated zone, which occurs at the leading edge.

The second test campaign aimed at improving the knowledge of the aerody­namic fbw field around the central blade of a straight blade cascade (Leconte et al., 2001). This test took place in ONERA R4 blow-down wind-tunnel for both steady and unsteady configurations. Conventional measurement techniques such as steady and unsteady pressure recording on the surfaces of the blades were used. To acquire a knowledge of the flow velocity field in the channels contiguous to the central blade of the cascade, PIV (Particle Image Velocime – try) technique was used. The test matrix featured various Mach numbers (sub – sonic, transonic and supersonic), cascade angle-of-attack, plunging and pitch­ing movements of the central blade.

Using the variable Q = (q, q)1 gives a first order in time differential equa­tion:

AQ + BQ = F * A=(J£) B=(“-1) F=(£)

At time tn+2 , a second order in time discretization is obtained using New- mark’s scheme.

With aerodynamic forces being calculated at time /" 1 ^ , a mechanical cou­pling iteration has to be performed in order to equilibrate generalized coordi­nate at time tn+1. The following numerical scheme is implemented:

■ Modal mesh deformations computation TIMELOOP: Physical time loop

Generalized coordinates estimate at t”+1 MECALOOP: Coupled fluid – structure equilibrium loop

Grid velocity computation DUALLOOP: Dual time loop

RKLOOP: Runge – Kutta loop ■ Unsteady Aerodynamic dual step END RKLOOP END DUALLOOP

Generalized coordinates and metric updating at tn+1 END MECALOOP if convergence criterium is reached

END TIMELOOP

Direct dynamic coupling methodology for moving meshes depends on the structural model. In the case of a linear or weakly non-linear model, the de­formation of the structure may be described using a modal basis. The grid deformation and velocity are then given by a linear combination of modal grid deformations and velocities at any time. The grid motion is interpolated. In the strongly non-linear case, structural deformations and velocities have to be computed at each time step of the coupled system from a finite element model. The grid’s motion has to be fully computed at each time step as well.

Assuming now that we use a linear structural model, the dynamic behavior of the structure is properly given by a linear combination of modal deforma­tions. We may write at any time t the vector 1І of the displacement at node M as:

Ft (M, t) = ^ qi (t) hi (M)

i

where qi (t) stands for the instantaneous generalized coordinates. Writing the mechanical principle in the modal basis gives :

Mq(t) + Dq(t) + Kq(t) = Fa (q, <q, t)

where M, D, K and Fa respectively stand for the mass, damping and stiffness matrices and the instantaneous aerodynamic forces.