Category UNSTEADY AERODYNAMICS, AEROACOUSTICS AND AEROELASTICITY OF TURBOMACHINES

The objective is to show the variations of amplitudes and phase lead to­wards bump motion of both the shock wave movement and the unsteady static pressure relatively to the reduced frequencies characterizing this experimental study.

2. Description of the experimental set-up

The test facility features a straight rectangular cross section. The oscillating model used in the here presented study is of 2D prismatic shape and has been investigated as non-vibrating in previous studies (Bron et al. [2001] ,Bron et al. [2003]), from where extensive baseline data are available. In order to intro­duce capabilities for the planned fliid-structure tests, a fhxible version of the model was built. Figure 1 shows the way the generic model oscillates in the test section and presents the optical access offered by this test facility. The ft>w

entering the test section can be set to different operating conditions character­ized by different inlet Mach number, Reynolds number and reduced frequency (Table 1). The generic model is molded of polyurethane, at defined elasticity (E=36.106MPa) and hardness (80 shore), by vulcanization over a steel metal bed. As shown in Figure 2, it includes a fully integrated mechanical actuator allowing smooth surface deformations. This oscillating mechanism actuates the fhxible model (bump) in a first bending controlled mode shape. While the highest point located at 57% of the chord vibrates in a sinusoidal motion of 0.5mm amplitude, the two edges of the chord stay fixed. A 1D laser sensor measures the model movement through the optical glass top window in one direction with a bandwidth of 20kHz and a resolution of +/-0.01mm. Time-

resolved pressure measurements are performed on the oscillating surface using pressure taps and Kulite fast response transducers. To achieve this, Tefbn tubes are directly moulded in the 2D fhxible generic model and plugged to the Kulite transducers mounted with the long line probe technique far from the os­cillating measured surface (Schaffer and Miatt [1985], see Table 2). These fast response transducers deliver signals with delays and large damping but exempt of resonance effect. The delays, damping, tubes vibrations and tubes elonga­tions have been carefully calibrated. All components of this test facility are fully described in Allegret-Bourdon et al. [2002]. The test section offers op­tical access from three sides (Figure 1). While the instantaneous model shape is scanned using the geometry measurement system through the top window, Schlieren measurement can be performed using the access through two sides windows. A high-speed video camera produces the Schlieren videos with a sampling frequency of 8kHz.

Figure 2. Cut view of the generic model (bump)

Davy Allegret-Bourdon

Chair of Heat and Power Technology

Royal Institute of Technology, Stockholm, Sweden

davy@energy. kth. se

Torsten H. Fransson

Chair ofHeat and Power Technology

Royal Institute of Technology, Stockholm, Sweden

fransson@energykth. se

Abstract A ffexible generic model has been developed at the Chair of Heat and Power Technology in order to perform ffitter experiments in a more fundamental fash­ion. It is made of engineered fexible material and oscillate in a controlled way at non-uniform amplitude and variable frequencies. Time-resolved measurements of the unsteady surface pressures, the instantaneous model geometry as well as unsteady Schlieren visualizations are performed in order to study the shock wave motion and the aerodynamic load acting over this ffexible generic bump. The model oscillates at reduced frequencies from 0.015 to 0.294 at transonic fow condition. The mode shapes of such a fexible bump strongly depends on the excitation frequency of the generic model. Schlieren pictures are obtained for an operating point characterized by an inlet Mach number of 0.63. More­over, the presented results demonstrate that 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 taps located after the shock wave formation, the phase of pressure fuctuations towards bump local motion presents the same decreasing trend.

Keywords: Fluid-structure interaction, Schlieren, Long line probe, shock wave movement,

Unsteady static pressure, first bending mode shape, Flexible generic model

409

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

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

Nomenclature

1. Introduction

Structure oscillating phenomena occur in many industrial applications in the field of energy technology. Under certain conditions a curved shape located in an uniform fbw, such as a blade, an airfoil or the surface of a nozzle, can enter into a self-excited vibration known as flitter. Under flitter condition, aerody­namic loads can rapidly increase the amplitude of vibration of a structure until its failure.

Experiments on controlled vibrating models have been performed by sev­eral researchers to investigate such aerodynamic loads and observe the shock wave motion related to it. Kobayashi et al. [1994] studied this relationship on an annular blade row oscillating in torsional mode with interblade phase angle. They drew the conclusion that at an inlet Mach number of 1.19 the shock wave movement significantly changed between the reduced frequencies

of 0.102 and 0.136, and that after this range of reduced frequencies the un­steady force damps the blade oscillation. Fujimoto et al. [1997] studied this unsteady fliid structure interaction on a transonic compressor cascade oscillat­ing in a controlled pitching angle vibration. They noticed that although the am­plitude of the shock wave displacement did not change much within the range of this experiment, the phase lag relative to the blade oscillation increases up to almost 90r as the blade oscillation reduced frequency increases to 0.284. Later, Hirano et al. [2000] performed other experimental campaigns on this transonic compressor cascade oscillating in a controlled pitching angle vibration. They conclude that the shock wave movement has a large effect on the amplitude and the phase angle of unsteady pressures on the blade surfaces; the amplitude of unsteady pressure becomes large upstream of the shock wave but decreases rapidly downstream; the phase angle across the shock wave changes largely for the surfaces facing the fl»w passages adjacent to the oscillating blade, the amplitude of shock wave movement increases following the increase of the re­duced frequency, and the phase angle relative to the blade displacement lags almost linearly as the reduced frequency increases.

In such kind of experiments, a driving system is creating an artificial oscil­lation of the rigid structure, whose amplitude and frequency can be controlled. The compressor blade of Lehr and Bolcs [2000], for example, is made oscil­lating in a controlled plunging mode by a hydraulic excitation system. The high-speed pitching vibrator of Hirano et al. [2000] is able to reach a 500Hz frequency of a 2D mode shape controlled oscillation in a linear cascade. In most of the cases, the vibrating structures are designed in metal to be close to real applications. Thus, large amplitudes of vibration at high oscillation fre­quencies prompt the failure of the structures. Moreover, recent research has presented a 2D blade harmonically driven in a 3D mode shape controlled vi­bration such as in Queune et al. [2000].

To date, this kind of flitter experimental investigations have been limited to stiff models made of metal, which oscillate in a pitching mode. Rather than studying the complex geometry of a turbomachine and specific industrial applications, the here presented generic experiments are voluntarily not taking into account inertial effects, radial geometry, numerous blades or 3D aspect of the flow occurring in industrial applications. Thus a generic oscillating flexible model is studied in order to reach a better understanding of the physics of the flutter phenomenon under transonic operating conditions.

A coupling calculation is performed on the same blade at a higher rotation speed for two inter-blade phase angles, o0 = 0° and or = 360°/22. The only changes in the aerodynamic conditions are: rotation speed ^ = 4516.8 rpm, pressure ratio P2/Рц = 1.12. For the inlet ft>w, the Mach number is 0.5, the velocity is 166 m/s and the mass ft>w is 487.95 kg/s. Like in the previous case, the coupled simulations are performed under transonic conditions.

A first coupling calculation is performed for an inter-blade phase angle oо = 0°. The projection basis is formed by the first two bending modes and the

first torsion modes whose frequencies are respectively 101.7 Hz (1F), 212.8 Hz (2F) and 448.9 Hz (1T). The generalized unsteady aerodynamic forces are computed at the frequencies given in Table 3.

Figure 4 shows the flitter diagrams obtained with the double scanning and the smoothing methods. Both methods give similar results : the blade is stable in the upstream infinite velocity range from 100 m/s to 250 m/s.

A time simulation is performed for an upstream infinite velocity = 166 m/s with the minimum state smoothing method. The aerodynamic forces are approximated using 6 poles with a relative error of 1.13% . The time incre­ment is determined to have 60 time steps per period (third mode). Figure 5 shows the time evolution of the three generalized coordinates q i, q2 and q3. Table 4 shows a Fourier analysis of those signals. The resulting frequencies and damping factors are similar to those obtained from the computation in the frequency domain (double scanning method).

w» Bow (Vys)

Figure 6 Case 2, crі – Flutter diagram

Figure 5 Case 2, <ro – Time evolution ol q,

A coupling calculation is now performed for an inter-blade phase angle a1 = 1 x 360°/22. The method is then tested for a non-zero phase angle and for a complex projection basis. The latter is formed by the first two bending modes and the first torsion mode whose frequencies are respectively 105 Hz (1F), 191.8 Hz (2F) and 278.3 Hz (1T). The generalized unsteady aerodynamic forces are computed for the frequencies given in Table 3.

Figure 6 shows the flitter diagram resulting from a computation using the double scanning method. As in the zero phase angle case, the blade is stable in the upstream infinite frequency range from 100 m/s to 250 m/s.

A time calculation is performed for an upstream infinite velocity = 166 m/s. The time increment is determined to have 60 time steps per period (third mode). The smoothing method uses 8 poles and gives a relative error of 1.18%. As for a zero phase angle, the frequencies and the damping factors resulting from the computations in the time domain (smoothing method) and in the frequency domain (double scanning method) are similar (Table 5).

A non-linear friction force (model of Coulomb) is applied on the blade for the phase angle a0. The coupled simulation is performed by projecting the motion equation of the structure on a Craig and Bampton basis which is com­posed of six normal modes and one static constraint modes, and by using the minimum state smoothing of the generalized aerodynamics forces in the time domain. Figure 7 shows the time evolution of the generalized coordinates q 1 (first normal mode) and q7 (static mode, on which the friction force is applied) for different intensities of the friction force: 0 N (linear case), 50 N and 500 N. We can see the infflence of the intensity of the friction force on the response.

2. Conclusion

Two fliid-structure coupling methods for a rotating bladed disk system based on the cyclic symmetry properties and the projection on the complex modes are presented. The double scanning method is more robust but it is only applica­ble for solving the flutter equation in the frequency domain. The minimum state smoothing method can be used for both frequency and time domain sim­ulations but its implementation is more complex and it requires more attention and know-how from the user. The main advantage of the proposed indirect cou­pling method in the time domain is that it is generally less time consuming than the direct coupling method since the aerodynamic calculations are performed only once at the beginning of the simulation and not at each time step. There­fore, the aerodynamic forces do not need to be computed again if the modes that generate them are unchanged and they can be re-used for other computa­tions, for example when the applied external forces change or when we want to study the influence of a friction damper on the stability of the bladed-disk. However, the proposed indirect coupling method is more restrictive and less accurate since it is based on the assumptions of linearized aerodynamics and harmonic motion, which is not the case of the direct coupling method.

The good agreement between the results of the different methods shows that the extension of the proposed coupling methods from aircraft applications to turbomachinary is feasible. However, further validation tests should be carried out under various conditions such as other rotation speeds, phase angles, pres­sure ratios, viscous fluid etc. The future work consists in numerical simulations of an unstable configuration and the introduction of the blade mistuning.

References

Craig, Jr., R. R. and Bampton, M. C. C. (1968). Coupling of substructures for dynamic analysis. AIAA Journal, 6(7):1313-1319.

Dat, R. and Meurzec, J.-L. (1969). Sur les calculs de fbttementpar lamethode dite du “balayage en frequence reduite”. La Recherche Aerospatiale, (133):41-43.

Dugeai, A., Madec, A., and Sens, A. S. (2000). Numerical unsteady aerodynamics for turbo­machinery aeroelasticity. In Proceedings of the 9th International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines, Lyon, pages 830-840. Karpel, M. (1982). Design for active flitter suppression and gust alleviation using state-space aeroelastic modeling. AIAA Journal of Aircraft, 19(3):221-227.

Poirion, F. (1995). Modelisation temporelle des systemes aeroservoelastiques. Application a l’etude des effets des retards. La Recherche Aerospatiale, (2):103—114.

Sayma, A. I., Vahdati, M., and Imregun, M. (2000). An integrated nonlinear approach for turbo­machinery forced response prediction. Part I: Formulation. Journal of Fluids and Structures, 14:87—101.

Tran, D.-M., Labaste, C., and Liauzun, C. (2003). Methods of fliid-structure coupling in fre­quency and time domain using linearized aerodynamics for turbomachinery. Journal of Flu­ids and Structures, 17(8):1161-1180.

The simulations were performed for the following aerodynamic conditions: rotation speed U = 4066.4 rpm, upstream total temperature Tц = 288 K, upstream total pressure Рц = 101325 Pa, pressure ratio P2/Рц = 1.05, phase angle do = 0 °. The steady aerodynamic simulation gives for the inflow a Mach number of 0.5, a velocity of 166 m/s and a mass ft>w of 462.92 kg/s. Figure 2 shows a shock whose position depends on the applied pressure ratio. The coupled simulations are then performed under transonic conditions.

For the coupling computation, only the first two bending modes of the blade whose frequencies are 97.4 Hz (1F) and 204.3 Hz (2F) are retained to form the projection basis. Unsteady aerodynamic simulations are performed using the mode shapes as oscillating motion shape for the frequencies given in Table 1.

In order to perform a time simulation, a minimum state smoothing of the generalized aerodynamic forces is performed with a relative error of 0.197%. The time increment is determined to have 60 time steps per period (second mode). An initial velocity is applied to all generalized coordinates. Figure 3 shows the time evolution of both generalized coordinates computed using three

methods : the first method is the indirect coupling in the time domain (smooth­ing method), the second one is a direct coupling numerical simulation using a grid deformation technique, and the last one is a direct coupling numerical simulation using a blowing condition. The results from the three methods are similar, although the aeroelastic damping computed with the direct simulation using the grid deformation technique is quite a lot smaller than those com­puted with the other methods. The good agreement between the results of the smoothing method and those of the direct coupling method with a blowing con­dition can be explained by the fact that a blowing condition was also used in the smoothing method. Table 2 gives the frequencies and the damping factors computed using the three methods. Moreover, the results from the smooth­ing method in the time domain is similar to the one from the double scanning method in the frequency domain. The time simulation has been performed for an upstream infinite velocity = 166 m/s.

The previously described coupling methods have been applied to a numeri­cal model of a compressor disk composed of 22 large chord blades, for differ­ent rotation speeds and phase angles. The structural finite element model of a reference sector with one blade has 19 539 degrees of freedom. The eigenfre – quencies and modes in vacuum are computed by using the cyclic symmetry and
by taking into account the stiffening and softening rotational effects. The aero­dynamic computations were performed on a two-block structured grid, each block having 61x18x30 points with 60 points on the profile. Figure 1 shows the mesh of one channel of the embedding Arid. For the coupling calculation, the blade is assumed not to have any structural damping. In addition to the indirect coupling methods described in this paper, and for comparison, we also use the direct coupling method in the time domain (Case 1).

The matrix An (Фп, i u) has been computed for nK reduced frequencies к1, …,кПк with the assumption of harmonic motion. For arbitrary motions like those defined by Eq. (6), it is necessary to extend the values of the aerody­namic coefficient matrix to an area of the complex plane containing the imagi­nary axis, i. e. to determine An (Фп, p) for p = iu(1 + ia) with a = 0.

The minimum state smoothing method Karpel, 1982 consists in modelling the generalized aerodynamic forces by using a rational approximation and aux­iliary state variables :

K*

gn

c*n 0

En – I

M* 0

gn

An0, An1, An2, Dn and En are matrices with dimensions (mn x mn) for An0, A„i and An 2, (mn x np) for Dn and (np x mn) for En. Rn = diag(ri,…,Гпр) is a diagonal matrix where np is the degree of the denominator of the rational function or the number of poles and r * < 0 are the poles. These matrices are computed by using a method of least squares minimization Poirion, 1995.

Mg-[3] [K*n + Gn (p)] – M;

*—1 c*

gn gn

– ^gn 4" 2 Poo 4(зо Ang, Сgn Cgn + 2 Poo ^ ^oo ^-n 1 >

Mgn = Mgn + ^Poo C[4] A„2 and G„ ip) = v£ {p c/O D„ [ip c/Vqq) I – Rn ]—1En. The matrices Gn (p) and H(p) are complex and depend on V^. The nonlinear eigenvalue problem Eq. (15) is solved for nV increasing velocities V1 ,…,Vnv by using an iterative process based on the method of successive approximations for finding a fixed point of a function Tran et al., 2003.

In Eq. (14), a vector of np auxiliary state variables zn have been defined. They satisfy a system of first-order differential equations in the time domain :

Using Eq. (14), the flitter equation (8) can be written under the form of a non-linear eigenvalue system of dimension 2 mn :

where K*n, C*n and M*n are the same matrices as in Eq. (15). This second – order system is solved using the Newmark numerical integration scheme.

For motions defined by Eq. (6) in the frequency domain, the aerodynamic coefficient matrix An (Фп, p) defined in Eq. (11) depends only on the quotient pc/V^ and can be written as An (Фп, pc/Vrx) where c is a reference length, for exemple the blade chord. Let’s define the reduced frequency к = u c/V^.

Using Eq. (11) and denoting An(Фп, pc/Vrx) and A^(Фп, pc/V! X>) the real and imaginary parts of A n (Фп, pc/Vrx:), the flitter equation (8) can be written

under the form of a non-linear eigenvalue system of dimension 2 m n :

orH(pc/V^) x = px, with K*„(pc/F00)=KS„ + y<X)V20k’n{&n, pclV00) and C*„ (p с/) = Cs„ + cp^ V^A" (Ф„, p с/) / (p с/).

If the damping factor is small (|a| ^ 1), we have pc/Vn ~ i к and the matrices A(Фп, pc/Vn), An(Фп, pc/Vn) and H(pc/Vn) can be approxi­mated by An (Фп, iк), An (Фп, iк) and H(i к). Eq. (12) becomes :

or H(i/c)x = px, with K*„(i/t) = Ks„ + ^Роо^А^Фп. ік), C*n(i/c) = cs„ + і cPqq А"(Ф„, in)/к and к = Srm(p) c/V^. H(i/t) is a real ma­trix which depends on An(Фп, iк), An(Фп, iк) and Vn. The solutions of the approximated eigensystem (13) should be determined for nV increasing velocities V^ ,…,Vcny and the frequencies should satisfy u = 3m(p) = кVcю/c. For each velocity V^, the following eigensytems are solved for i = 1,…, 2 mn and j = 0,1,2,…, until the convergence on к is obtained:

H(i «i, j ) xi, j + 1 pi, j +1 Xi, j + 1, кІ, І U*,j c/Vrx ^m(pi, j ) c/Vn,

where (p*-, xj j) are the i-th complex eigenvalue and eigenvector obtained at the j-th iteration. The matrix H(i к*-) is computed by interpolating A, n (Фп, i к) and An (Фп, i к)/к from the tabulated values. It is remarked that the double scanning method is not valid for small velocities since the latter lead to very large reduced frequencies and the extrapolation of H(i к*-) will be incorrect.

The unsteady aerodynamic forces are computed from a basis of mn real mode shapes Ф of the reference sector, for an oscillation frequency и and an inter-blade phase angle an. By expressing the displacements of the reference sector as a linear combination of the modes Ф and by assuming that the struc­tural motion is harmonic and that all the sectors have the same motion with a constant phase angle an between two adjacent sectors, the generalized aerody­namic forces generated by the displacement un (t) are written, by linearity:

fagn (un (t), un (t)) Fagn (^, i t) fin •

Fagn (Ф, i u, t) is the time-dependant aerodynamic coefficient matrix whose (i, j)-term is obtained by projecting the unsteady aerodynamic force generated by the harmonic motion of the j-th mode on the i-th mode:

Fagn, ij (Ф, iu, t) = – / [Pn (M, , iu, t) – Ps (M)]tt(M). (M) dS

JM ЄЕ

(11)

where Pn is the unsteady pressure generated by the j-th mode ^j, Ps the steady pressure, Фі is the displacement vector of the i-th mode, it is the unit external normal vector to the surface S and dS is an elementary surface of S. Let’s introduce the aerodynamic coefficient matrix An (Ф, iu, t) obtained from the integral in Eq. (10) with Pn and Ps replaced by the associated pres­sure coefficient Cp = (P — P00)/(^ Poo V^>), where P^, Рос and are the pressure, the density and the velocity of the upstream unperturbed fhid.

By performing a Fourier analysis of Fagn (Ф, i u, t) and by keeping only the first harmonic term, we have:

Fagn ^, iu, t) – Fagn (ф iи)Є

The generalized aerodynamic forces generated by un (t) become:

fagn (un (t)> un (t)) – Fagn ( Ф, i u) fin еІ

Fagn (Ф, i u) and An (Ф, i u) are complex, asymmetric square matrices of di­mension mn. They are computed for иш oscillation frequencies.

In the flitter equation (8), the aerodynamic coefficient matrix should be de­termined from the mn complex modes Фп. Using the linearity hypothesis, the aerodynamic coefficient matrix Fagn (Фп, i u) generated by Фп can be ex­tracted from the aerodynamic coefficient matrix computed from the basis of the 2 mn real vectors formed by the real and imaginary parts of Фп.

The unsteady aerodynamic forces are obtained solving the Euler equations for an ideal gas using an aerodynamic code called CANARI and developed for years at ONERA Dugeai et al., 2000. Because of the cyclic symmetry of the flow, a chorochronic boundary condition is applied to the simulated channel :

F(r,9 + к (3,z, t) = F(r,9,z, t + for n

F is any fbwfield variable, r is the rotation radius, and 9 the azimuthal an­gle. In a first step, a steady state is computed depending on the rotation speed, on the pressure ratio and on the far-field total temperature, total pressure, and velocity. In a second step, unsteady simulations are performed by forcing an oscillating blade motion at different frequencies. These simulations depend on the steady flowfield previously computed and used as initial conditions, on the inter-blade phase angle, and on the forced motion shape and frequency. A blowing condition is then used to simulate the blade motion. Once a pseudo­steady oscillating state has been reached (no transient effect), a Fourier trans­form is performed over the pressure to get the unsteady aerodynamic forces.

A structure with cyclic symmetry is composed of N identical sectors S0, Si,…, SN – i which close up on themselves to form a circular system. The whole structure is obtained by N — 1 repeated rotations of a reference sector S0 through the angle в = 2 n/N. Each sector is limited by a left and a right frontier L and Lr with the adjacent sectors. The fliid surrounding the structure is also assumed to have the same cyclic symmetry while the applied external forces can vary arbitrarily from one sector to another sector.

Using the cyclic symmetry properties, the motion equation of the structure comes down to N motion equations of the reference sector S0, in terms of the travelling wave coordinates and with the appropriate second members and boundary conditions. Only sector S0 has then to be modelized and, after a finite element discretization, the following reduced matrix systems will be solved to obtain the vector of the travelling wave coordinates un = un (S0,t) of sector S0, for each phase number n = 0,…,N — 1 and with an = пв:

Kun + Clin + Miin = fan (un, un) + fn + rn, (1)

1 N—1 . N—1

fan = ^ E fa(^)e_ifc<Jre and = ^ E (2)

k=0 k=0

K is the stiffness matrix of sector S0, including the centrifugal stress stiffening and the spin softening effects, C is the damping and gyroscopic effect matrix and M is the mass matrix. fa (Sk) is the vector of the aerodynamic forces
applied to sector Sk and which depends on the displacements and the velocities of sector Sk. f(Sk) is the vector of the other external forces applied to sector Sk, including the centrifugal forces. rn is the vector of the interface reactions applied to the frontiers of So with the adjacent sectors, it does not intervene in the solutions of Eqs. (1-3) and it is only present due to the contraints Eq. (3). The cyclic symmetry boundary conditions Eq. (3) are expressed in the cylindrical coordinate system. It is remarked that, the aerodynamic forces fan applied to the coordinates un in Eq. (1) depend only on un and not on the other travelling wave coordinates um for m = n, and they are equal to the physical aerodynamic forces fa (S0, un, Un) induced by un on sector S0.

For each phase number n, the travelling wave coordinates are expressed as a linear combination of the first mn complex modes Фп of the undamped, rotating structure in vacuum (solutions of Eqs. (1-3) without C, fn and fan):

un = фп qn• (4)

Introducing Eq. (4) in Eq. (1) and premultiplying by 4Ф„, we obtain the reduced coupled system in the mn complex modal coordinates qn (t):

Kgn qn + Cgn qn + Mgn fin fagn (Фп qn, Фп qn) + fgn, (5)

un(t) = un ept and qn(t) = qn ept with p = iш (1 + ia), (6)

where ш is the unknown aeroelastic eigenfrequency (ш > 0) and a is the un­known aeroelastic damping factor (a Є ). Using the hypothesis of linearity, the generalized aerodynamic forces are written as :

fagn (Фп qn, Фп qn) Fagn (Фп ,p) qn e • (7)

Substituting Eqs. (6) and (7) in Eq. (5), we obtain the flitter equation :

[ Kgn + P Cgn + P2 Mgn – Fagn (Фп, P) ] qn = 0, (8)

which is a complex, nonlinear eigenvalue system in which the aerodynamic coefficient matrix Fagn (Фп, p) depends on the complex modes Фп and the unknown complex eigenvalue p.

When structural nonlinearities such as friction or free-play are to be intro­duced in the reduced coupled system Eq. (5), we need to keep some physical coordinates among the generalized coordinates qn. These "nonlinear” coordi­nates un|NL can be for example the displacements of the reference sector nodes located at the junction between the blade and the disk, where friction dampers are introduced. For this aim, the Craig and Bampton Craig and Bampton, 1968 projection basis Qn = [Фп, Фсп ] is used instead of the eigenmodes :

un = фп n + ^cn un|NL = Qn qn • (9)

Фп are the first mn complex normal modes of So with the cyclic symmetry conditions and with un|NL fixed; Фсп are the complex static constraint modes of S0 with the cyclic symmetry conditions and a unit displacement on one coordinate of un|NL, while the remaining coordinates of un|NL are fixed; and n is the vector of the modal generalized coordinates. By projecting Eq. (1) on Qn, we obtain a complex reduced coupled system similar to Eq. (5) whose size is the number of eigenmodes mn plus the number of the nonlinear coordinates.

Duc-Minh Tran and Cedric Liauzun

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 minh. tran_duc@onera. fr, cedric. liauzun@onera. fr

Abstract Two methods of fUid-structure coupling for turbomachinery are presented, the first one in the frequency domain and the second in both frequency and time domains, with the assumptions of linearized aerodynamics and cyclic symmetry.

Keywords: fliid-structure coupling, aeroelasticity, turbomachinery

1. Introduction

This paper is concerned with the coupled fliid-structure dynamic analysis of turbomachinery. The structure consists of a rotating bladed disk submitted to the unsteady aerodynamic forces exerted by the fliid, which are themselves generated by the structural motion. The structure and the fluid are assumed to have a perfect circumferential cyclic symmetry, so that the classical reduc­tion of the analysis to only one reference sector can be applied. The study of the structure comes down to that of the reference sector by applying the ap­propriate boundary conditions for each phase number. The displacements of the reference sector in the travelling wave coordinates are expressed as a lin­ear combination of the complex modes or Craig and Bampton’s basis and the motion equations are projected on these bases to obtain a reduced system.

In the coupling methods proposed here Tran et al., 2003, the unsteady aero­dynamic forces are assumed to depend linearly on the structural displacements and velocities and they are expressed in terms of those induced by the modes. The mode-induced aerodynamic forces are computed only once at the begin­ning of the simulation by using an aerodynamic code (solving the Euler equa­tions) with the assumption of harmonic motion of the modes, for an inter-blade phase angle and a number of reduced oscillation frequencies.

397

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

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

In the frequency domain, the projection of the mode-induced aerodynamic forces on the modes provides a complex matrix of aerodynamic coefficients whose product with the modal coordinates gives the generalized aerodynamic forces, leading to a nonlinear eigenvalue system. This flitter equation is solved by using two well-known iterative techniques already used for aircrafts, which are the double scanning method Dat and Meurzec, 1969 (also called p — k method) and Karpel’s minimum state smoothing method Karpel, 1982.

In the time domain, Karpel’s minimum state smoothing of the aerodynamic coefficient matrix is used to obtain a time-domain approximation of the gen­eralized aerodynamic forces by means of auxiliary state variables. Structural nonlinearities such as friction or free-play are taken into account. This nu­merical strategy is applied to a compressor blade for different configurations (phase numbers, rotation speeds). The results are compared between the pro­posed methods and to those of the direct coupling method Sayma et al., 2000 where the structural and the fluid equations are solved alternately at each time step and the assumptions of linearity and harmonic motion are not necessary.

The flow past a circular rod is an interesting academic study about vortex shedding. It will be used to pursue a deeper analysis, taking advantage of the extensive experimental database. Three URANS k — и turbulence models will be compared, in order to estimate the infhence of this parameter. Also, a large-eddy simulation (LES) has been carried-out. It enables an interesting comparison and creates an outlook.

4 Description of the confi guration. In the present study, the ft>w past the rod is characterized by a Reynolds number of Red = 48,000 [the infbw velocity is U ж = 72 m/s and the rod diameter: d = 0.01 m]. Accord­ing to the experiments, the vortex shedding occurs at the Strouhal number: St = fd/Ucv ~ 0.2, and the transition to turbulence takes place in the shear layer, just after the separation (sub-critical regime). As a consequence, the wake is fully turbulent, though dominated by the von Karman vortices. This ft>w topology is sufficiently complex to be compared to practical applications.

The results analyzed in this paper are extracted from a larger configuration previously used to study the emission of broadband noise. This configuration involves a symmetric airfoil (chord length: c = 10 d) located in the wake of a rod, one chord length downstream. We are here concerned with the vortex shedding, so only the flow past the rod will be analyzed. The studies have shown that the airfoil has nearly no influence on the rod because of the distance between the two bodies. Information about the study on the whole rod-airfoil configuration can be found in references Jacob et al., 2002 and Boudet, 2003.

5 Numerical characteristics. Figure 6 presents the overall rod – airfoil computational domain and a view of the mesh near the rod. The URANS mesh is 2D and involves 197 points around the circumference, with the first grid line located at У + < 8 from the wall. A 3D mesh is not necessary be­cause URANS only features the mean components of the ft>w, which are bi­dimensional in this case.

The LES mesh is the same, but duplicated along span to cover 3 diameters (31 points), in order to simulate the tri-dimensionality induced by the turbu­lence. In the turbulent region, the resolution is characterized by y + < 3 per-

Figure 6. Rod test-case. Left: Rod-airfoil computational domain (7 blocks). R/ght: mesh near the rod (1pt/2).

pendicularly to the wall, x+ < 50 tangentially to the curvature of the wall, and z+ < 200 along the span.

Three URANS computations were carried-out, using respectively three k — и models: the linear model of Wilcox, 1993b, the low-Reynolds number model of Wilcox, 1993a, and the non-linear model of Shih et al., 1995. A no-slip condition is imposed at the walls, and a non-refecting condition reduces the spurious confinement at the outer boundaries. The computation was led to convergence for 10 aerodynamic cycles, then it was saved 100 times per cycle for 18 cycles.

The LES computation used the auto-adaptive model of Casalino et al., 2003, to estimate the subgrid scale viscosity. As for URANS, no-slip and non­refecting conditions are imposed at the relevant boundaries, and a slip condi­tion is applied to the limiting plane in the spanwise direction. The computation was led to convergence for 6 aerodynamic cycles (initial field: URANS), then it was recorded 100 times per cycle for 18 cycles.

4.0.6 Results. The vortex shedding occurs naturally for all the computa­tions. However, as mentioned in table 2, URANS overestimates the frequency: the Strouhal number is St = 0.24 for the three k — и models, whereas the experimental value is St = 0.20. LES performs better, it predicts the experi­mental value.

Table 2 also presents the integral forces on the rod: the mean drag (< CD >), the fhctuating drag (C D), and the fhctuating lift (C L). The lin­ear URANS model underestimates these forces. As a comparison, the low – Reynolds number model performs slightly better. The non-linear model is fi-

nally better than the two other URANS computations. It performs as well as LES, they both exhibit results in agreement with the experimental ranges.

Next, Fig. 7 presents the mean pressure and the mean friction on the wall. Concerning the mean pressure, we first notice a discrepancy in the recirculation region, represented by the plateau near the angle 180 deg. This turbulent recir­culating region is the most sensitive. However, the non-linear URANS and the LES exhibit better results. The region near the separation is more interesting. There, all the URANS computations fail, exhibiting similar curves. They em­phasize the minimum, and this can be related to a delaying of the separation. As a comparison, LES achieves to predict the experimental curve.

The mean friction confirms the preceding comments. The separation, lo­cated by Cf = 0, is delayed by the URANS computations, whereas the LES features the experimental location.

The overestimate of the Strouhal number by the URANS can be related to the delaying of the separation. The mean recirculation bubble is thiner than

in the experiments, and the vortex shedding frequency is consequently higher. Also, the good estimate of the shedding frequency by LES is consistent with the simulation of the separation, located in agreement with the measurements.

Figure 8 presents the mean and flictuating velocities in a section perpen­dicular to the wake, at x/d = 7.5 downstream of the rod axis. This location is influenced by the airfoil, but computations and experiments are consistent, carried-out for the same configuration (rod + airfoil). We notice similar results for the three URANS computations, but level differences appear. The URANS wake is thiner than in the experiment, in agreement with the delaying of the separation. Also, the two peaks on each side of the axis show that the von Karman vortices are convected along the same path from one cycle to another. This is not true for the experiment and the LES, because they take the turbulent dispersion of the trajectories into account.

Finally, Fig. 9 shows views of the instantaneous vortex identification func­tion, Г2, designed by Graftieaux et al., 2001. The experiment (PIV), the linear URANS and the LES are presented. We notice a good agreement between LES and PIV, many small vortical scales agglomerates to shape the major eddies. On the contrary, URANS only exhibits large von Karman vortices, with a small inter-vortices distance due to the higher frequency.

Figure 10 presents the same quantity, but reconstructed using only the mean field and the two major POD modes (POD: Proper Orthogonal Decomposi­tion). The agreement is good between the three images. It shows that URANS is able to feature the major modes, in agreement with PIV and LES.

.

2. Conclusion

The present studies focused on the abilities of the numerical approaches to predict vortex shedding. Firstly, various U-RANS results has been presented on the VKI test-case. Then, an academic study was pursued on an isolated rod undergoing a sub-critical vortex street. One observes that RANS is able to predict the main features of the fbw, qualitatively. However, differences are noticeable between the various turbulence models used. As an outlook, LES performs better because of its direct description of the flow, including the diversity of the scales. For example, boundary layer separation has an important influence on the vortex shedding, and it is accurately featured by LES.

Unsteady RANS simulations have been achieved for the VKI turbine case. A coarse grid (18 000 nodes, referenced as stp2) was firstly used, but in order to resolve the shocks appropriately, a finer grid has been also designed (72 000 nodes, referenced as stp1). The sonic region and the Von Karman street can be seen on the instantaneous map of the computed Mach number contours on Fig. 1-left.

The time averaged isentropic Mach number distributions along the pressure and the suction sides are shown on Fig. 2. The computational results (lin­ear model of Wilcox, 1993b, and non linear model of Craft et al., 1996) are

very close to the experimental results. The prediction is very good on the pressure side, and a good agreement is noticeable on the suction side where the numerical distribution is above the experimental points, probably a time – integration effect of the experimental data. The unsteadiness due to the vortex shedding modify mainly the area located downstream of the nearly sonic area (x/c > 0.5), on the suction side.

The stabilization on one specific frequency is difficult to obtain with the non­linear k — и model. It seems that unsteady shock boundary layer interactions was not completely stabilized. All simulated Strouhal numbers are quite close to the experimental value of 0.219, as shown in Table 1.

Results are better on the coarse grid than on the fine grid. Consequently, one needs to be careful when analyzing the numerical results. k — и models give a frequency a little too high, which suggests a too small value of the turbulent kinetic energy.

Figure 1-right shows the density iso-contours for the experiment and for the linear and non-linear k — и models. The vortex locations and the shock wave structures are close, but the intensity of the numerical results seems higher than the experimental results. The vortex location obtained by experiment is well compared to linear and non-linear k — и model (Fig. 3-left). For the numerical results, the distance between the lower and upper vortex is a little too wide suggesting that the amplitude of the flictuations is too strong near the trailing edge.

The density at the vortex center is compared to the extrema of VKI exper­imental data for linear and non-linear k — и model (Fig. 3-right). The exper­imental results show a symmetric vortex shedding. The linear and non-linear k — и models increase the deficit of density which shows that the fhctuation level is over-predicted as was suggested on Fig. 1-right. For the linear model, the modification is non-symmetric and affects mainly the pressure side of the trailing edge. As for the non-linear k — и model, it predicts correctly the lower vortex intensity, but over-estimates the upper vortex intensity.

Time averaged static pressure at different angular locations around the trail­ing edge is presented in Fig. 4-left. With the linear model, the agreement is reasonable on the pressure side, whereas more discrepancies are visible on the suction side. Otherwise, the non-linear model gives better results on the suc­tion side except at 80 degrees where the shock wave appears. It seems that the computation is not completely stabilized (locally). It will be noticed that the experimental values are almost symmetric, the numerical solution appearing more influenced by the small supersonic regions upstream of the trailing edge. More investigation will be necessary on spatial schemes and limiters. An in­teresting point is to observe the comparison between the time averaged value

and the steady value. The higher level on steady static pressure confirms that the vortex shedding increases the losses.

The excessive amplitude of the deficit given by the linear k — и model is also visible looking the instantaneous static pressure flictuations due to vor­tex shedding in Fig. 4-right. The curve shapes are similar, but the amplitude is higher for the finer grid, particularly on the suction side where the increase of flrctuations represents 40% of the total flrctuation amplitude. This is con­sistent with a sharper description of local gradients. Two decompositions of raw experimental data were done at VKI using phase-locked Fourier analysis (solid line) or wavelet analysis (dots). The amplitude measured using wavelet analysis is in better agreement with the numerical results. This could be ex­plained by a slight random frequency variation during the experimental vortex shedding to which the Fourier analysis is more sensitive.

Total temperature and total pressure flictuations were recorded 2.5 diame­ters downstream the trailing edge (Fig. 5). Once again, the curve shapes are similar, particularly on the pressure side, the amplitude being slightly over pre­dicted.