Category UNSTEADY AERODYNAMICS, AEROACOUSTICS AND AEROELASTICITY OF TURBOMACHINES

Dual time stepping interest. Aeroelasticity and fliid-structure coupling computations are usually performed using a very small global time step value. This is especially true when studying low frequency phenomena. Therefore a very large number of iterations is required, and this leads to very expensive computations. In the viscous case, large and very refined meshes are used, and the time step requirements for numerical stability are even more critical. Moreover, moving meshes computations are required, which increases CPU costs.

This is the reason why the use of dual time stepping for Navier-Stokes aeroe – lastic computations becomes very interesting. The physical time step used to describe the unsteady phenomenon is no longer constrained by stability time step values in the smallest cells. At each physical time step, a modified steady problem is solved in a dual pseudo time. Usual convergence acceleration tech­niques such as local time stepping or multi-grid scheme may be used. Dual local time steps are bounded by specific stability requirements.

As far as moving meshes computations are concerned, the dual time step technique helps to reduce the number of remeshing computations. Dual time iterations are performed at a fixed physical time step, that is to say in a fixed mesh. This is much more important in the case of coupled fhid structure com­putations, where the position and velocity of the grid is not prescribed, but derives from the resolution of the coupled equations.

1.1.1 Dual time stepping scheme for moving meshes. Let us consider

du

equation (E) : — + /(u) = 0, where и is a numerical function of time. Writ­ing 2nd order Taylor’s expansions for u(tn) and u(tn-1) at time tn+1 allows us to write equation (E) at time tn+1 as follows :

du* u(t„_i) — 4u(t„) + 3u* ~dt* H 2At

+ f (u*) + o(At)

u(tn+1) is then obtained as the solution at steady state of a differential equation for the variable u*, function of the pseudo-time t*:

Pseudo-time t* is called dual time. The resolution of the unsteady problem is now performed within a system of two time loops. The external one is the physical time loop. The inner one is the dual time loop. The dual time loop is carried out using local time stepping, because it solves a steady state problem. The usual four step Runge-Kutta scheme is used to perform the dual loop res-

Q*(k) = Q

with Q*(0)

At* ■

Vn-1Qn-i – 4VnQn + 3Vn+1Q*(k-1) ]

= Qn, Q*(n+1)

= Q*(4),

olution. Within this loop, the grid is fixed at physical time tn+1. Applying this approach to the conservative fliid equations in moving grid, we obtain:

RVq-1) = ^ F(q-1) NSi – VTV0) and F(q-1) = F(q-1) + Fd(q-1) + D(0)

i=1

In these formulae, At and At* denote the physical and dual time steps, respectively. At steady state of the dual loop, the conservative variables vector Q is obtained at time t(n+1). The stability condition for the dual time scheme depends on the value of the physical time step as follows:

4

At* < CFL x – At 3

This condition is added to the those given by the properties of the Jacobian matrices of the convective and viscous flixes.

Several turbulence models are available. The first one is the mixing length turbulence model of Michel. It gives an expression of the turbulent viscosity as a function of a mixing length depending on the local distance from the wall and on the boundary layer thickness.

The second one is the one-equation model of Spalart-Allmaras (Spalart et al., 1992). In this case, one transport equation for the kinematic turbulent viscosity v is added to the set of mean fbw equations. Several algorithms may be used, with either strong or weak coupling between mean flow and turbulent equations at each time step.

The third model is the Launder-Sharma ke two-equation model. Using this model, the mean flow equations are closed by two additional equations for the turbulence kinetic energy and for its dissipation rate. Low-Reynolds number corrective terms are used.

The time-step is adapted to ensure the stability of the conservative turbulent equations. Numerical 2nd and 4th order viscosity terms are added as well as limiting functions for the turbulent variables.

In the following numerical validations, only the Spalart-Allmaras model will be considered.

In order to reduce the size of the unsteady harmonic response computations, a specific time-space periodicity boundary condition is used for the turboma­chinery numerical simulations. The cascade is supposed to be made of N ge­ometrically and structurally identical blades. Therefore, using this boundary condition, only a single channel needs to be meshed and computed for various inter-blade phase angles, in order to obtain the values of the unsteady aerody­namic forces for the complete cascade.

The chorochronic boundary condition is applied between the upper and lower boundaries of the channel. This condition reads:

2n

F(x, R,d, t) = F(x, R,d + j — ,t + jo)

where F is any function, 9 the azimuthal angle and a the inter-blade phase angle, a is defined by n 2тг – у where 0 < n < N — 1

In order to reduce the storage, the fbw field at the chorochronic boundaries is simply stored at a reduced number of time steps during the cycle. The field at inner time steps is then rebuilt using a specific time interpolation technique.

Numerical techniques have been developed at ONERA for 2D and 3D mesh deformation (Dugeai et al., 2000). They are based on a linear structural anal­ogy, with discrete spring networks or continuous elastic analogy. A finite ele­ment formulation is used, and special features allow the reduction of the size of the problem, especially in the Navier-Stokes case.

In the case of the spring analogy, two different techniques have been de­veloped. The first one is the method proposed by Batina (Batina, 1989), and the second one is an extension in which the 3 components of the displacement vector are coupled.

In the case of the continuous elastic analogy, 8-node hexahedral finite ele­ments are used to discretize the problem of the deformation of a linear elastic medium. The local stiffness matrix is computed using a numerical Gauss in­tegration procedure with a cheap but not exact one Gauss point integration, which leads to Hour-Glass modes terms. A special procedure is used to re­move the singularity of the stiffness matrix, giving satisfactory enough results for the grid deformation purpose.

For both approaches, spring network or elastic material analogy, the static equilibrium of the discretized system leads to the following linear system:

Kii qi Kif qf

where qi and qf are respectively the induced and prescribed displacement vectors. As the stiffness matrix is positive definite, the system is solved us­ing a pre-conditioned conjugated gradient method. The technique has been implemented in the case of multi-block structured grids. The full mesh defor­mation is defined as a sequence of individual block deformations. Additional conditions are set on the boundaries to impose zero or prescribed displacement values, and to get a continuity of the deformations at block interfaces.

A macro-mesh technique is used for large grid sizes, which is often the case in 3D Navier-Stokes computations. The macro-mesh is defined from the original one by packing several cells, typically 2, 3, or 5 cells, in each direc­tion. In the case of Navier-Stokes meshes, the whole boundary layer region is packed, in normal direction, in a single macro-cell. The coarse macro-mesh is then deformed using the structural analogy techniques, and the inner node displacements are finally interpolated in each macro-cell.

In this section, we present the numerical features of the ALE Navier-Stokes code Canari (Vuillot et al., 1993). This 3D code solves Euler and Reynolds – averaged Navier-stokes equations in multi-block structured grids.

1.1 Space and time ALE discretization for the mean flow

Unsteady Navier-Stokes computations have to be performed in a moving grid framework. An ALE (Arbitrary-Lagrangian-Eulerian) numerical scheme has therefore been developed. The spatial discretization is based on a centered finite volume approach. The fliid motion equations are written in a frame, which rotates at circular frequency D. In this frame, the grid is moving at velocity Vg :

where Q = (p, pW, pE)t is the unknown field vector, Fc[Q, Vg] the convective flix, Fd[Q, gradQ] the diffusive flix and T[Q] the source term. W is the relative velocity of the fluid and E the relative total energy in the rotating frame.

The time integration (Jameson et al.,1981) is performed using a Jameson – like four stage Runge-Kutta scheme. Second and fourth order artificial viscos­ity terms are added to the original scheme in order to obtain suitable dissipative
properties. The implicit spectral radius method of Lerat (Lerat et al., 1982) is used to increase the stability domain.

Anne-Sophie Rougeault-Sens and Alain Dugeai

Structural Dynamics and Coupled Systems Department Office National d’Etudes et de Recherches Aerospatiales B. P. 72, 29 avenue de la Division Leclerc, 92322 Chatillon Cedex, France Anne-Sophie. Sens_Rougeault@onera. fr, Alain. Dugeai@onera. fr

Abstract This paper presents ONERA’s recent advances in the experimental and numeri­cal understandings about the aeroelastic stability of aeronautical turbomachiner­ies. Numerical features of a quasi-3D and a 3D Navier-Stokes unsteady aeroelas­tic solver are discussed: turbulence models, grid deformation techniques, spe­cific boundary conditions, dual time stepping. A dynamically coupled fliid – structure numerical scheme is presented. Isolated profile, rectilinear cascades computational results are compared to experimental data. Results of aeroelastic Navier-Stokes computations for 3D fans are shown.

Keywords: fbid-structure coupling, aeroelasticity, turbomachinery

1. Introduction

For several years, ONERA has been interested in aeronautical turbomachin­ery aeroelasticity studies. The goal of this research has been to improve the ex­perimental and numerical knowledge about the aeroelastic stability and forced response of aeronautical turbomachineries.

One of the main challenges in this matter concerns the prediction of the aeroelastic stability of fans, especially in the case of the transonic regime. In this case, the dynamic behavior of the boundary layer needs to be accu­rately predicted using RANS numerical modeling with transport equations tur­bulence models. Numerical simulations have to be performed in a deforming grid framework using an Arbitrary Lagrangian Eulerian formulation.

In order to perform validations of the developed numerical tools, several unsteady data bases were built first for an isolated profile, and then for a rec­tilinear cascade. Theses databases have been extensively used to conduct nu­merical unsteady Navier-Stokes aeroelastic validations.

Another point concerns computational time reduction. Unsteady aeroelastic Navier-Stokes computations are extremely time-consuming due to the small

423

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

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

time-step needed to keep the numerical scheme stable in the small boundary – layer cells. This is the reason why the numerical technique of dual time­stepping has been implemented in the various unsteady Navier-Stokes codes used at ONERA. This technique allows one to reduce the time of aeroelastic Navier-Stokes computations in such a way that simulations that would have been unaffordable using global time stepping are now possible.

The last purpose of this paper is to show some results for direct fliid – structure coupling simulations. A coupled scheme using Newmark’s time dis­cretization has been developed and implemented in our aeroelastic Navier – Stokes codes (Girodroux et al.,2003). Coupled time domain simulations have been performed in the case of a compressor fan blade.

This paper presents some of the unsteady aerodynamic numerical develop­ments and results of the experimental campaigns. Some results of the valida­tion processes of the 2.5D and 3D aeroelastic Navier-Stokes codes will be de­tailed. An example of a dynamically coupled 3D Navier-Stokes Arid-structure computation will be given.

The unsteady pressure fluctuations are measured along the bump and the corresponding unsteady pressure coefficient and phase leads towards bump motion are deduced for five chosen pressure taps. The amplitudes of the unsteady pressures fluctuations shift significantly at the reduced frequency k=0.221 for the pressure taps located 20% upstream and downstream of the

200 400 600 800 10OO

200 400 600 800 10OO

1 (Hz) 1 — 200Hz

600 lOOO 1500 2000

f – 200Hz

500 lOOO 1500 2000

f — 200Hz

500 10OO 1500 2000

f (Hz)

Figure 6. Time-variant and power spectra of static pressure, shock wave movement and bump top motion at 10Hz, 75Hz and 200Hz perturbation frequencies

bump axial chord as shown in Figure 8. Moreover the unsteady pressure co­efficients remain stable and range between 2 and 4 for the three pressure taps located within 40% to 80% of the bump axial chord. The phase lead towards bump motion of the static pressure flictuations range between 90Deg. and 180Deg. for the pressure taps located before the bump max height, and be­tween -180Deg. and 90Deg. for the pressure taps located after the max bump height. At the pressure tap located close to the shock wave mean location (67% of the bump chord) and at y/H=0.25, the phase leads towards bump motion follow the same decreasing trend. In comparison with the shock wave motion phase variation, a global decrease in phase close to 270Deg. is observed for the pressure taps located after the shock wave.

3. Conclusion

Phase relations among oscillatory bump motion, shock wave movement and unsteady pressure fluctuations are investigated in the case of a flexible generic model controlled-oscillated in bending mode shapes at an inlet Mach number of 0.63, over a range of reduced frequencies from 0.015 to 0.294. The follow­ing conclusions are drawn:

• The mode shapes of such a fhxible bump strongly depends on the exci­tation frequency of the generic model.

The phase of shock wave movement towards bump local motion shows a decreasing trend for the third bending mode shapes at reduced frequency higher than k=0.074.

At the pressure tap located after the shock wave formation (67% of the bump chord), the phase of pressure flictuations towards bump local mo­tion presents the same decreasing trend as for the shock wave movement analysis.

For those same pressure taps, lower and stable pressure coefficients are also observed.

Acknowledgements

Time-variant signals and corresponding power spectra of pressure fhctua – tions, bump top motions and shock wave movements are shown in Figure 6 for three bump oscillatory frequencies (10Hz, 75Hz and 200Hz). Both pres­sure fhctuation and shock wave motion signals seem to follow the shape of the sinusoidal signal generated by the bump displacement at the oscillatory fre­quencies of 10Hz, 75Hz and 200Hz. At these three excitation frequencies, the pressure fhctuation and shock wave motion power spectra show the same clear fundamental harmonic. The bump top location movement power spectra con­tains one supplementary higher harmonic component that is not shown here. It does not exist in the power spectra of the pressure flictuation and shock wave motion signals. It is interpreted as being linked to external mechanical vibra­tions coming from the oscillation drive train and the wind tunnel. All three oscillations seem to be of a sinusoidal type after ensemble averaging posttreat­ment.

2.2 Schlieren visualization results

Figure 7 characterized the measured oscillations of the shock wave up to k=0.294. The mean location of the shock stays the same for all excitation frequencies. Moreover one can notice that the amplitude of the shock wave oscillations increases slightly from 0.015 to 0.294. The first bending mode shape at k=0.015 is characterized by a phase lag towards bump motion close to 315Deg., and the phases range between 30Deg. and 90Deg. for the second bending mode shape from k=0.03 to k=0.074. The phase decreases signif­icantly from 270Deg. to almost 0Deg. at reduced frequencies higher than k=0.089 for what has been considered as a third bending mode shape.

At this operating condition, the generic bump is controlled-oscillated in bending mode shapes at frequency between 10 and 200Hz. For each oscillating frequency, the synchronized data of the bump motion, shock wave movement and static pressure flictuations are acquired. The shock wave motion is mea­sured at one vertical location corresponding to 15mm (y/H = 0.25) over the top of the bump neutral position (it is symbolized by the white dashed arrows in Figure 3). Figure 5b shows successive pictures of this shock wave oscillat­ing at 10Hz oscillation frequency. A reference line indicates the mean location of the shock wave (67% of the bump chord). From t’=0 to t’=0.250, the shock wave moves through its mean position in an upstream direction. From t’=0.500 to t’=0.750, the shock wave moves again through its mean position in a down­stream direction. Due to the sinusoidal oscillation of the bump, the shock wave

stays a longer time in the two extreme positions (upstream and downstream) and crosses quickly its mean position during one period at 10Hz bump os­cillatory frequency. Figure 5a shows in the same way one oscillation of the vertical part of the shock wave at 200Hz excitation frequency. These pictures demonstrate a movement close to be sinusoidal.

2.1 Description of the operating condition

In these experiments, inlet and outlet time averaged isentropic Mach num­bers are set and a time averaged "lambda" shock wave is generated over the generic model surface at 67% (+/-1%) of the bump chord. Figure 3 shows a typical shape of the shock wave created during those experiments. To de­fine this operating condition, the stagnation pressure and the stagnation tem­perature of the fbw are measured (Pt = 159kPa at Tt = 305K) at ten chords upstream and the corresponding isentropic inlet Mach number is cal­culated (Misoi = 0.63). The downstream static pressure is measured on the ground wall and allows calculation of a downstream isentropic Mach number MiSo2 = 0.61 at two chords after the generic model. Figure 4 shows the chord wise distribution of local static pressures for the same operating condition. The generic model acts as a contraction of the channel. Miso2 decreases until 10% of the chord and then increases until 50% of the chord where the fbw speed is maximal. Then Miso2 decreases through the shock wave formation. Because of the manufacturing method, the pressure taps are not exactly perpendicular to the surface and thus do not measure the exact static pressure profile as well as the unsteady pressure ffictuations.

Figure 4 describes the way the generic model is oscillating. A regular repar­tition of the amplitudes along the bump half chords shows a maximal defor­mation at x/cax = 0.47. Due to its fbxible nature, a first bending mode shape

Schlieren picture of the shock wave created in the test section (Miso1 = 0.63, 0.61) and isentropic Mach number profile at upper and lower bump positions

at k=0.015 changes in a second bending mode shape at k=0.074, and reaches a third bending mode shape at higher reduced frequencies. At the mean shock wave location x/cax = 0.67, the local geometry presents a phase towards bump top motion. This phase is 20Deg. at k=0.015, 45Deg. at k=0.037, 120Deg. at k=0.074, -180Deg. fromk=0.11 to k=0.147, -45Deg. at k=0.221 and -90Deg. at k=0.294.