Nt1310 Unit 3 Lab Report

erical methPage 48

Licentiate Thesis / Nenad Glodic

6 NUMERICAL INVESTIGATIONS 6.1 Numerical Method
Numerical simulations are carried out employing a commercial CFD code (ANSYS CFX v11). The solver is using a full-scale time-marching 3D viscous model. Underlying equations, three dimensional Navier-Stokes equations in their conservation form, are being solved by using a Finite Volume method, where equations are integrated over the finite control volumes. Thereby, the solution domain is subdivided into a finite number of control volumes employing a suitable grid, which defines the control boundaries around a computational node in each control volume center. 6.1.1 Governing equations In fluid dynamics, the fluid flow is governed by the conservation laws for mass, momentum and energy. The basic conservation laws are formulated by using Leibniz-Reynolds transport theorem, which is an integral relation stating that the changes of some intensive property defined over a control volume must be equal to what is lost (or gained) through the boundaries of the volume plus what is created/consumed by sources and sinks inside the control volume. The
threedimensional conservation equations (mass conservation, x-,y- and z-momentum equations and energy equation) for compressible fluid can be given in differential form (Versteeg & Malalasekera, 1995) as follows
      u   0 t Du   p   xx   yx  zx      S Mx Dt x y z Dv  xy   p   yy   zy     S My  Dt x y z Dv  xz  yz   p   zz       S Mz Dt x z z   u xx   u yx   u zx       y z  x    v   v   v    DE xy yy zy    u           kT   S E  Dt y z  x        w xz    w yz   w zz    x y z   

Eq. 6-1 Eq. 6-2 Eq. 6-3 Eq. 6-4

Eq. 6-5

where the source terms S Mx , S My , S Mz include contributions due to body forces only. The specific energy of the fluid E in Eq.6-5 is defined as the sum of internal 1 (thermal) energy i and kinetic energy u 2  v 2  w 2 , while the effects of 2 gravitational potential energy changes are included as energy source term





Licentiate Thesis / Nenad Glodic

Page 49

(gravitational force regarded as a body force). It is common practice to extract the changes of the kinetic energy in order to obtain an equation for internal energy i . This is done by taking the part of energy attributable to the kinetic energy, which can be obtained by multiplying the momentum equations with corresponding velocity components and adding them together, and subtracting it from the above presented energy equation (Eq.6-5)
  xx   Di   p  u    Dt  zy  u u   yx   zx x y v w   xz   yz z x u   xy z w   zz y v v    yy  x y        kT  S i w  z 

Eq.6-6

Among the unknowns in the above presented conservation equations are four thermodynamic variables:  , p, i and T . The relationship between these thermodynamic variables can be obtained through the assumption of thermodynamic equilibrium, where the state of substance in the equilibrium can be described by means of just two state variables. This yields, for a perfect gas, the well-known equations of state: p  RT and i  cv T

Eq. 6-7

In the flow of compressible fluids the equations of state provide the linkage between the energy equation on one hand and mass conservation and momentum equations on the other. The governing equations also contain the viscous stress components  ij , and in a Newtonian fluid the viscous stresses are proportional to rates of deformation or strain rate. For compressible flow, Newton’s law of viscosity involves two constants of proportionality: dynamic viscosity  , to relate stresses to deformation, and second viscosity  , to relate stresses to volumetric deformation. The second viscosity  is very difficult to determine and is often neglected. The shear stress components can thereby be expressed as

  w u  v  xx  2    u    u  yy  2     u y x z  u v   v w   u w   xy   yx       xz   zx       yz   zy      y x   z y  Eq. 6-8   z x      The introduction of the above stress components into momentum equations (Eq. 6-2 - Eq. 6-4) yields the so called Navier-Stokes equations, which after some rewriting can be given as

 xx  2

Du p      u   S Mx Dt x Dv p       v   S My Dt y Dw p       w  S Mz Dt z



Eq. 6-9 Eq. 6-10 Eq. 6-11

Page 50

Licentiate Thesis / Nenad Glodic

Introducing a Newtonian model for viscous stresses in internal energy equation (Eq. 6-6) gives after some rearrangement  Di    p  u    kT     S i , Eq.6-11 Dt where all the effects due to viscous stresses are described by the dissipation function  . Having now a system of seven equations (five partial differential flow equations and two algebraic state equations) with seven unknowns, the system is mathematically closed and can properly model dynamics of time-depended threedimensional fluid flow of a compressible Newtonian fluid. The introduction of a general variable  into the conservative form of the fluid flow equations discussed above, yields a so-called general transport equation describing the transport of a fluid property  through a continuum

      ( u )       S Eq. 6-12 t The left hand side of the equation contains a rate of change term and the convective term, while on the right hand side the diffusive term (where  is diffusion coefficient) and the source term are present. The equation 6-12 is used as a starting point for computational procedures in the finite volume method, where the transport equation is being integrated over a three-dimensional control volume.
6.1.2 RANS and Turbulence modelling The majority of flows of engineering significance are turbulent. Turbulence consists of fluctuations in the flow field in time and space. It is a complex, threedimensional, unsteady process consisting of a wide range of length and time scales. Turbulence occurs when the inertia forces in the fluid become significant compared to viscous forces, and is characterized by a high Reynolds Number. In principle, the Navier-Stokes equations describe both laminar and turbulent flows without the need for additional information. However, turbulent flows at realistic Reynolds numbers span a large range of turbulent length and time scales, and would generally involve length scales much smaller than the smallest finite volume mesh, which can be practically used in a numerical analysis. Although being possible, the direct numerical solution of the Navier–Stokes equations for turbulent flow is extremely difficult and expensive in computational efforts. To enable the effects of turbulence to be predicted, time-averaged equations such as the Reynolds-averaged Navier-Stokes equations (RANS), supplemented with turbulence models (such as the k - or k -model), are used. Due to the averaging procedure, information from the full Navier-Stokes equations is lost and it is supplied back into the code by the turbulence model. The idea behind the equations is Reynolds decomposition, whereby an instantaneous quantity is decomposed into its time-averaged and fluctuating quantities. The Reynolds decomposition defines flow property  as the sum of a steady mean (time

Licentiate Thesis / Nenad Glodic

Page 51

average) component  and a time varying fluctuating component  t  with zero mean value. Hence Eq. 6-13  (t )     (t ) The mean or averaged component is given by: t  t 1   (t )dt t  t

Eq. 6-14

where t is a time scale that is large relative to the turbulent fluctuations, but small relative to the time scale to which the equations are solved. For compressible flows, the averaging is actually weighted by density (Favreaveraging). For simplicity, it is assumed that density fluctuations are negligible, but mean density variations are not. Substituting the averaged quantities into the instantaneous transport equations results in Reynolds averaged flow equations for turbulent compressible flows: Continuity: ~       U   0     t

Eq. 6-15

Momentum: ~   U P ~~ ~    u2   uv   uw       UU       U         S Mx Eq. 6-16   t x x y z     ~   V P ~~ ~    uv   v2   vw      VU       V         S My Eq. 6-17   t y x y z     ~   W P ~~ ~   uw   vw   w2       WU       W         S Mz Eq.6-18   t z x y z    

   

















 
 









Scalar transport equation: ~    ~~ ~    u    v    w          U            S   z  y x t 





Eq. 6-19

In the presented equations the “overbar” indicates a time-averaged variable and the “tilde” indicates a density-weighted or Favre-averaged variable (Versteeg & Malalasekera, 1995) Extra terms appear in the Reynolds averaged flow equations due to interactions between various turbulent fluctuations. These additional terms (   vi v j ) are socalled Reynolds stresses, arising from the non-linear convective term in the unaveraged equations. On the energy equation an equivalent term is found, the Reynolds flux. In order to be able to compute turbulent flows with RANS equations, it is necessary to develop turbulence models to approximate Reynolds stresses and close the RANS equation system. The RANS turbulence models are classified on the basis of the additional transport equations that need to be solved together with

Page 52

Licentiate Thesis / Nenad Glodic

the RANS flow equations. Most widely used are the two-equation models (such as k   or k   model), consisting of two additional partial differential equations for different turbulence quantities. The model used in the present work is the standard k   model, where the gradient diffusion hypothesis is used to relate the Reynolds stresses to the mean velocity gradients and the turbulent viscosity. The model delivers an additional equation for turbulent kinetic energy k and one for the dissipation rate  . The equations are coupled via eddy viscosity with the momentum equations. The model is based on the eddy viscosity concept stating that the effective viscosity is achieved as the sum of dynamic and turbulent viscosity, which describes the increase of the viscosity due to turbulent fluctuations. The turbulent viscosity is modeled as the product of a turbulent velocity and turbulent length scale. The turbulence velocity scale  and length scale l are defined by using turbulent kinetic energy k and eddy dissipation  as follows:
2

 Applying dimensional analysis, the eddy viscosity  t can be specified as

 k

,

l

k3

Eq. 6-20

t  C l  C

k2



Eq. 6-21

where C is a dimensionless constant. The standard k   model (Launder and Spalding, 1974) uses the following transport equations for kinetic energy and eddy dissipation:

    k     kU     t k   2 t Sij .Sij   t  k         2    U     t    C1 Sij .Sij  C2  t k k   





Eq. 6-22 Eq. 6-23





where

S ij

is the mean component of the rate of deformation, while

C  , C1 , C 2 ,  k ,   are adjustable constants, usually having values that are arrived

at by comprehensive fitting for a wide range of turbulent flows. At high Reynolds numbers the viscous sublayer of a boundary layer becomes very thin, so that an increasing mesh density near the wall is necessary. In situation like this the standard k   model avoids the need to integrate the model equations all the way to the wall surface by making use of universal behaviour of the flow in this region. In order to model the near-wall flow properly without resolving the turbulence equations in the boundary layer, the wall functions are being used. Thus a lot of computational effort is being saved. These empirical functions relate the local wall shear stress to the mean velocity relaying on the existence of a logarithmic region in the velocity profile. Wall functions need the first calculation node to be placed on the log-law region. At low Reynolds numbers the log-law is not valid so modifications to the k   model equations (Eq. 6-21 through 6-23) needs to be performed (reviewed in Patel et al. ,1985).

Licentiate Thesis / Nenad Glodic

Page 53

6.2 Description of Numerical Model
The simulated domain contains a sector of seven freestanding blades where for the unsteady simulations a harmonic sinusoidal motion of the center blade in various rigid-body modes has been employed. The motion of the oscillating blade is described by a set of equations and imposed as moving mesh boundary condition. The geometry of the modeled sector is depicted in Figure 6.1.

Figure 6.1: Modeled sector cascade The effects of having a finite cascade of blades on implementation of the influence coefficient technique has been previously addressed by Vogt (2005) where a comparison between results obtained from simulations in the influence coefficient domain and results from travelling wave mode simulations were compared. The outcome of that investigation showed that the employed method of simulating in the influence coefficient domain in circumferentially limited sector cascade is numerically correct and the observed differences were negligible. In present investigations the boundary conditions on the sidewalls of the modeled sector were pre-set to rotationally periodic. Based on the findings presented by Vogt (2005), where a high accurate third order Euler approximation is compared with a lower first order discretization scheme, it was decided that a second order backward Euler approximation should be used in the present investigations. Although it was demonstrated that low order numerical schemes can generally be used for these kinds of flow situations, the choice was made to use higher order scheme due to the fact that the first order scheme may lead to minor differences in response magnitude thanks to an increased numerical damping behaviour. For the discretization of the convective terms, the High Resolution advection scheme is selected. The temperature throughout the flow is predicted with the Total Energy heat transfer model. It models the transport of enthalpy and includes kinetic effects and is especially suitable and recommended for compressible flows. All simulations are conducted using a standard k-ε turbulence model with wall functions. In order to be able to relate the current numerical results with the ones obtained by Vogt (2005), the same grid was also employed in the present work. The mesh was created in VOLMOP grid generator, a part of an in-house code from Volvo Aero. The mesh topology consists of an H-O multiblock structure with an O-gird around

Page 54

Licentiate Thesis / Nenad Glodic

the blades and H-grid in the blade-to-blade passages. Mårtensson and Vogt (2005) carried out an extensive mesh sensitivity study where several different mesh topologies and resolutions were tested (Figure 6.2). The unsteady simulations were performed using an in-house CFD solver employing a linearized Euler method for prediction of aeroelastic properties. The simulations were conducted in the travelling wave mode with phase-lagged periodic boundary conditions. The conclusion was that unsteady results depend both on proper mesh distribution and mesh resolution i.e. having a too coarse mesh and disadvantageous mesh distribution might seriously affect the correctness of the unsteady simulation results. The differences in unsteady results between medium and fine meshes were however found to be small. Hence, the medium mesh was chosen as a good compromise between accuracy and computational effort. The current simulation domain consists of 507276 hexahedral volume elements with 540498 nodes.

Coarse

Medium

Fine

Figure 6.2: Different investigated mesh types; Mårtensson (2005) In order to optimize simulation effort, tip clearance has not been modelled. Results from previous investigations (Vogt et al., 2007) have shown that models without tip clearance feature a more favourable convergence while the prediction accuracy at midspan is not affected considerably. The investigation carried out by Glodic et al. (2012) confirmed these findings and showed that model without tip clearance was equally capable to predict the unsteady response at midspan as the model with tip clearance included. At the inlet of the domain the total pressure profile boundary condition has been applied. To account for the viscous effects from the hub and shroud, a total pressure distribution at the inlet includes a turbulent boundary layer and is generated by specifying the values of total pressure, static pressure and the boundary layer thickness expressed as a percentage of the channel height. 5% was retained to be suitable considering the characteristics of the test case. At the outlet the average static pressure value was specified. Turbulence intensity was pre-set to medium level (5%). Both the Tu level as well as the BL thickness has been measured. BL thickness of 5% has been verified, while the Tu level was measured to 2%. However, the simulation results obtained with 5% Tu level did not change when the level was set to 2%. Hub, shroud and blades are defined as smooth walls with no-slip condition. The side walls are rotationally periodic interfaces. The flow regime in the domain is predefined as subsonic.

Licentiate Thesis / Nenad Glodic

Page 55

6.3 Simulation approach
The simulations using ANSYS-CFX were performed as follows: first a steady state solution for a specific operating point was obtained.
This solution is then used as initial value for the transient simulation. In the unsteady simulations, the motion of the center blade was prescribed by a set of equations for the requested mode. The correctness of the imposed motion is confirmed by an analytical model. A timemarching solution was acquired spanning typically 3 oscillation periods. The flow could be regarded as time periodic already after the second oscillation cycle (criterion used: >0.5% phase-locked difference). A time trace of unsteady pressure coefficient at a chosen point on the blade surface is presented in Figure 6.3. One oscillation period was resolved by 20 time steps and with three iteration loops per time
step.
Time history at point span 0.5 arc −0.2 on blade 0 −3.4 osc.period 1 osc.period 2 osc. period 3 osc.period 4

−3.42 pressure coefficient Cp, −

−3.44

−3.46

−3.48

−3.5

−3.52

0

0.2

0.4 0.6 normalized time, t/T

0.8

1

Figure 6.3: Time history of unsteady response on the surface of the oscillating blade ANSYS CFX offers a possibility of applying mesh deformation on particular mesh regions if the motion of geometry domain in this region is known. In this case the motion of the blade is determined in terms of node displacements relative to the initial mesh. The nodes of the mesh are displaced in accordance with specified values defined by expressions written in CFX Expression Language (CEL) code. Mesh deformation of the other sub-domains or boundaries are either unspecified (as in case of mesh on the hub and shroud) or constrained (all remaining subdomains or boundaries). The oscillation of the blade 0 is assumed to be a 3D rigid body movement where the blade can oscillate in three orthogonal modes (axial bending, circumferential bending and torsion) or in combined bending-torsion modes. Rigid body movement is an idealization where deformations are ignored. Accordingly, the blade movements are described as ordinary rigid body motion via vector analysis by rotating the blade around the corresponding axes by the angle φ. The timedependent angle  is defined as

 t    max  cos

 2  t    T 

(6.24)

Page 56

Licentiate Thesis / Nenad Glodic

This angle of rotation has in current simulations never exceeded value of 0.25 degrees i.e. maximum displacement during the blade oscillations was considerably less than 1 % of the blade chord. Hence, it is assumed that the principle of linear superposition could be applied. Simulation data post-processing includes a reduction of time dependent data to complex unsteady pressures and calculation of stability data. Time-resolved blade surface pressure data is exported form ANSYS CFX into Matlab, where transformation into harmonic complex pressure data is performed. 6.3.1 Aerodynamic mistuning simulations The variation of blade-to-blade stagger angle is introduced in the same manner as the blade oscillation. The stagger angle of an individual blade is changed by rotating the blade around the torsion axis with a targeted de-stagger angle. The rotation is applied in a first time step of transient simulation. The blade oscillation is thereafter introduced first after the flow field around the de-staggered blade has stabilized and a steady state is reached (typically after a time period corresponding to two oscillation periods). The blade is thereafter oscillated for three oscillation periods as the solution is converged already after the second period of oscillation. The applied de-stagger angles were in range of -2.5° to +2.5° degrees. In order to assess the perturbed influence coefficients matrices described in section 2.6, in addition to the nominal case simulation, at least three more simulations per each investigated stagger have been performed: one simulation where de-stagger is applied on the oscillating blade and two additional simulations where each of the neighbouring blades was de-staggered. Together with the nominal geometry simulation, it is in overall 4 simulations needed if the aerodynamic influence coefficients on blades -1, 0 and +1 are considered. The resulting perturbed matrix is of tri-diagonal character where diagonal terms are the influence of the blades on themselves and the off-diagonal terms are containing the influence of the neighbouring blades. Aeroelastic stability analyses of a randomly mistuned blade rows are in this case constrained by a requirement to have sufficient spacing between the de-staggered blades to assume that the influence of the two de-staggered blades is not affecting each other. Consequently, between two de-staggered blades there must be at least four nominal blades to avoid interference between perturbation influences.