Direct numerical simulation of 2D transonic flows around airfoils

Computers & Fluids 88 (2013) 19–37

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Direct numerical simulation of 2D transonic flows around airfoils Tapan K. Sengupta ⇑, Ashish Bhole, N.A. Sreejith Department of Aerospace Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India

a r t i c l e

i n f o

Article history: Received 4 May 2013 Received in revised form 7 August 2013 Accepted 20 August 2013 Available online 28 August 2013 Keywords: Navier–Stokes equation Transonic flow Airfoil aerodynamics Shock capturing Shock-boundary layer interaction

a b s t r a c t High-accuracy, time-accurate compressible Navier–Stokes solvers have been developed for transonic flows. These solvers use optimized upwind compact schemes (OUCS) and four-stage, fourth order explicit Runge–Kutta (RK4) time integration scheme, details of which can be obtained in Sengupta [Sengupta TK. High Accuracy Computing Methods, fluid flows and wave phenomena. UK: Cambridge Univ. Press; 2013]. Although these compact schemes have been developed originally for direct simulation of incompressible flows, it is shown here that the same can be used for compressible flows, with shock-boundary layer interactions clearly captured for flow past NACA 0012 and NLF airfoils. Numerical higher order diffusion terms which are used for incompressible flows, have been replaced here by the pressure-based artificial diffusions proposed by Jameson et al. [Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Ruge–Kutta time stepping schemes. AIAA Paper 1981-1259. AIAA 14th fluid and plasma dynamics conference. Palo Alto, CA; 1981]. Such second and fourth order diffusion terms are used adaptively at selective points, located by the pressure switch. Developed computational methods used here are validated for cases with and without shocks, for which experimental results are available. Apart from surface pressure coefficient, contours of physical quantities are presented to explain the time-accurate results. Presented methods are robust and the results can be gainfully used to study shock formation, drag divergence and buffet onset of flow over airfoils. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Transonic flows represent many complex phenomena, such as creation of shock, contact discontinuities, shock-boundary-layer interaction and buffet onset. Such flows are strongly time-dependent and cannot be fully understood without reference to viscous nature of the fluid, as viscous terms are important in their interaction with the convective process [1]. In real viscous flows, changes in flow field occur in non-monotonic fashion and produce surface fluctuations, even in the absence of unsteady boundary conditions [2]. Therefore, transonic flows must be described by the timedependent, compressible Navier–Stokes equations. Numerical methods for solving Navier–Stokes equations must be capable of predicting all the above flow features accurately. In [3], authors summarize the developments in CFD ranging from solution of potential equation to Reynolds Averaged Navier–Stokes equation (RANS) and use of software packages to model transonic flows. Use of numerical diffusion and turbulence models are common in the simulation of transonic flows. In [4,5] turbulence models have been used, whereas direct numerical simulation (DNS) of transonic flows can be found in [6]. The authors report DNS results to investigate upstream moving pressure waves over a supercritical airfoil. ⇑ Corresponding author. Tel.: +91 512 2597945. E-mail address: [email protected] (T.K. Sengupta). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.08.007

In the literature, various numerical schemes have been proposed for calculation of compressible flows. For example, Allaneu and Jameson [7] proposed kinetic energy preserving scheme; Chiu et al. [8] have used a conservative meshless scheme. Other notable references are with implicit factored scheme [9], essentially non-oscillatory (ENO) shock capturing scheme [10], weighted essentially nonoscillatory (WENO) scheme [11]. In the present work, we use high accuracy, dispersion relation preserving (DRP), optimized upwind compact schemes (OUCS), which were originally developed for incompressible flows [12–14]. DNS of two-dimensional wall bounded turbulent flow is also reported in [15], where OUCS3 scheme for spatial discretization have been used. In [16], symmetrized version of the OUCS4 scheme have been used in parallel computing framework, to solve three-dimensional unsteady compressible Navier–Stokes equations for supersonic flow past a cone–cylinder for M1 = 4 and Re1 = 1.12  106. Parallel computing technique have been developed in this reference by solving the convection of a shielded vortex in subsonic flow. By solving and validating two-dimensional transonic flows around NACA 0012 and SHM-1 airfoils using two OUCS schemes, we show here that the same compact schemes can be used for transonic flows as well. Since the validity of computation is based on available experimental results, knowledge of experimental conditions and possible errors associated with the data should be known quantitatively. In [5,17], problems related to transonic wind tunnel

20

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

Farfield boundary

(a) -0.6

COMPUTATIONAL EXPERIMENTAL

-0.4 -0.2

Cp

0 0.2 0.4

η

0.6

Periodic boudary U∞

0.8

ξ

1 0.2

0.4

0.6

0.8

1

x/c

(b)

NACA0012 airfoil surface

COMPUTATIONAL

-0.8

EXPERIMENTAL

-0.6

No-slip wall boundary

-0.4

Cp

-0.2 0 0.2 0.4 Fig. 1. Computational domain and boundary segments.

0.6 0.8 1

(a)

0.2

0.4

0.8

1

Fig. 3. Comparison of computed and experimental Cp-distribution for flow around NACA 0012 airfoil with Re1 = 3  106 at a = 0.14° and (a) M1 = 0.6 and (b) M1 = 0.758. Computed data in both cases is time-averaged over t = 10–40.

0.1

y/c

0.6

x/c

0.2

0 -0.1 -0.2

0

0.2

0.4

0.6

0.8

COMPUTATIONAL

1

x/c

EXPERIMENTAL

-0.8

(b)

-0.6

0.2

-0.4 -0.2

0

Cp

y/c

0.1

0.2

-0.1 -0.2

0

0.4 0

0.2

0.4

0.6

0.8

1

x/c Fig. 2. Close-up view of orthogonal grid generated around (a) NACA 0012 and (b) SHM-1 airfoils.

0.6 0.8 1 0

tests, such as wall interference, three-dimensional effects in two-dimensional tests, aeroelastic effects, flow unsteadiness have been discussed which affect the reliability of the data. Corrections for wind tunnel wall effects in experimental data are with uncertainty, and may not be directly used for prediction by CFD. The wall corrections manifest as incremental velocity, angle of incidence, drag and streamline curvature [18] for subsonic wind tunnel. But, when supercritical flow region cannot be considered small with respect to the tunnel dimensions, wall corrections become more complicated as distortions of the flow field [19], which can strongly influence the flow and shock. As discussed in

0.2

0.4

0.6

0.8

1

x/c Fig. 4. Comparison of computed and experimental Cp-distribution for flow around NACA 0012 airfoil with Re1 = 3  106 at a = 0.14° and M1 = 0.779. Computed data is time-averaged over t = 10–45.

[20], flow unsteadiness is noted for velocity, pressure and temperature as fluctuations. Most transonic tunnels have high level of unsteadiness, in particular pressure fluctuations, which have significant adverse effects. In case of free-transition, such noise has influence on transition location [21]. Free stream turbulence can

21

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

also have an influence on attached boundary layers to cause separation [22]. Authors in [5] concluded that available experimental data on transonic airfoils are insufficient for turbulence-model validation. In [23], authors have also concluded that comparison with wind tunnel have been limited to cases with insignificant wall effects, because of the lack of appropriate data. The experimental results for two-dimensional section of NACA 0012 reported in [24], are used here for the validation purpose. Although, the wind tunnel used for experiments was not designed for two-dimensional testing, the geometry of the model with large span-to-chord ratio of 3.43 was useful for two-dimensional testing. Also, the width of the tunnel helped to minimize sidewall boundary layer effects [24] making conditions suitable for two-dimensional test. The correction due to wall-induced lift interference is suggested as incremental change in angle of attack in [25].

E^v ¼ b0; sxx ; sxy ; ðusxx þ v sxy  qx ÞcT F^v ¼ b0; syx ; syy ; ðusyx þ v syy  qy ÞcT In Eq. (1), the variables q, u, v, et, p, T represent nondimensionalized density, Cartesian velocity components, total specific energy, pressure and absolute temperature of the fluid, respectively. The total specific energy et is given by 2 2 et ¼ e þ ðu þ2 v Þ, where e represents specific internal energy and sxx, sxy, syx, syy are the components of the symmetric viscous stress tensor defined as

  1 @u 2l þ kr  V Re1 @x   1 @v ¼ 2l þ kr  V Re1 @y   l @u @ v ¼ syx ¼ þ Re1 @y @x

sxx ¼ syy sxy

The heat conduction terms are given by 2. Governing equations and boundary conditions The governing equations of motion are the conservation laws for mass, momentum and energy given here by the unsteady, compressible Navier–Stokes equation. In nondimensionalized strong conservation form, these can be expressed in the physical (x, y)plane as

b @b @Q E @b F @b E v @ F^v þ þ ¼ þ @t @x @y @x @y where the conservative variables are given by

b ¼ bq; qu; qv ; qet cT Q The convective flux vectors b E and b F are given as

b E ¼ bqu; qu2 þ p; quv ; ðqet þ pÞucT b F ¼ bqv ; quv ; qv 2 þ p; ðqet þ pÞv cT and the viscous flux vectors b E v and b F v are given as

ð1Þ





l @T PrRe1 ðc  1ÞM 21 @x   l @T qy ¼  2 PrRe1 ðc  1ÞM 1 @y

qx ¼ 

All terms in Eq. (1), are nondimensionalized with scale as the chord (L = c) for length, the free stream density (q1) for density, pffiffiffiffiffiffiffiffiffiffiffiffi the free stream velocity (U 1 ¼ M1 cRT 1 ) for velocity, (q1 U 21 ) for the pressure, the free stream temperature (T1) for temperature and (c/U1) for the time, where c denotes the ratio of specific heats, R denotes the universal gas constant and M1 as the free stream Mach number. The non-dimensional parameters, Reynolds number (Re1), Prandtl number (Pr) and Mach number (M1) are defined as,

Re1 ¼

q1 U 1 c lC p U1 ; Pr ¼ ; M1 ¼ l1 j a1

where j denotes the thermal conductivity and a1 denotes the speed of sound at the free stream. The two coefficients of viscosity are related by Stokes hypothesis: 3k + 2l = 0.

Fig. 5. Distribution of mean streamwise velocity at various streamwise locations in wall coordinates, for flow around NACA 0012 airfoil with Re1 = 3  106 at a = 0.14° and M1 = 0.779.

22

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

In Eq. (1), additional information is given by Sutherland’s viscosity law relating the coefficient of viscosity (l), with T as 3=2 l ¼ 1:45T  106 . The closure of system of equations is provided Tþ110 by the ideal gas equation given by p = qRT and it is used to define 2 2 et as, et ¼ cv T þ ðu þ2 v Þ where cv is the specific heat at constant R volume and is given by, cv ¼ c1 . Eq. (1) is transformed to the

(a)

generalized curvilinear coordinates, by introducing transformation [n = n(x, y), g = g(x, y)]. The transformed equation in the strong conservative form can be written as

@Q @E @F @Ev @F v þ þ ¼ þ @t @n @ g @n @g

(f)

γξ

γη

0.05

0 .0 0 1

0.001

y

0.3

ð2Þ

y

0.2

0

min = 0 max = 0.1246

0.1

-0.05 0.2

0.3

0.4

0.5

0.6

0.7

0.8

min = 0 max = 0.0232 1

(b)

1.1

1.2

(g)

Dξ1

Dη1

0.05

0.3

y

5 E -0

y

0.2

0

min = -0.00014 max = 0.00018

0.1

-0.05 0.2

0.3

0.4

0.5

0.6

0.7

(c) y

0.8

0.1

0.5

0.6

0.7

0.8

1.1

1.2

(i)

Dξ3

0.3

min = -0.000064 max = 0.000059 1

Dη3

0.05

0.2

y

y

(d)

0.4

Dη2

05

0

-0.05 0.3

1.2

1E-

0.05

min = -0.000094 max = 0.000078 0.2

1.1

(h)

y

0.2

min = -0.000032 max = 0.000044 1

Dξ2

0.3

6

-0

-0.05 0.2

0.3

0.4

0.5

0.6

0.7

0.8

(e)

min = -0.00012 max = 0.00011 1

1.1

1.2

(j)

Dξ4

2E

0.05

-0

Dη4 5

y

0.3

5

0

0

min = -0.000036 max = 0.000043

1E

0.1

0.1

y

0.2

0

min = -0.00042 max = 0.00051 -0.05 0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

min = -0.000102 max = 0.000143 1

1.1

1.2

Fig. 6. Added diffusion terms for the case of flow over NACA 0012 airfoil shown in Fig. 4: Frames (a) and (f) show contours for values of i,j terms (Eq. (9)) in n- and gdirections, respectively. Frames (b–e) show contours for values of Dn’s (Eq. (7)) added in the numerical solution of Eq. (2), in n-direction. Frames (g–j) show contours for values of Dg’s (Eq. (8)) added in the numerical solution Eq. (2) for the JST model, in g-direction.

23

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

where corresponding state variables and flux vectors are given as

OUCS3

    b =J; E ¼ nx b E þ ny b F =J; F ¼ gx b E þ gy b F =J Q¼Q     Ev ¼ nx E^v þ ny F^v =J; F v ¼ gx E^v þ gy F^v =J;

OUCS2 EXPERIMENTAL

-0.5

The grid metrics nx, ny, gx and gy of the transformation are defined as

nx ¼ Jyg ;

ny ¼ Jxg ;

gx ¼ Jyn ; gy ¼ Jxn 0

J¼

Cp

and the Jacobian of transformation is given by

1 ðxn yg  xg yn Þ 0.5

The viscous shear stress components in the transformed plane are given by



 2 nx un þ gx ug  ny v n þ gy v g 3 Re1 3

 2  l 4 ¼ ny v n þ gy v g  nx v n þ gx ug 3 Re1 3 i l h ¼ syx ¼ n un þ gy ug þ nx v n þ gx v g Re1 y

sxx ¼ syy sxy



l 4

1 0

0.2

0.4

0.6

0.8

1

x/c Fig. 7. Comparison of computed and experimental Cp-distribution for flow around SHM-1 airfoil with Re1 = 13.6  106 at a = 0.27° and M1 = 0.62. Computed data is time-averaged over t = 10–20.

and heat conduction terms in the transformed plane are given by

  l nx T n þ gx T g ; 2 PrRe1 ðc  1ÞM1   l qy ¼  n T þ g T n g y y PrRe1 ðc  1ÞM21

qx ¼ 

On the airfoil surface, no-slip condition is imposed on the velocity components and adiabatic wall condition is imposed on heat conduction terms

u ¼ 0;

v ¼0

and qx ¼ 0;

qy ¼ 0

A typical sketch for the flow past an airfoil is shown in Fig. 1, which shows that the problem is solved in an O-grid topology. This requires introduction of a cut, as marked in the figure, along which periodic boundary condition is implemented for all flow variables. Free stream initial condition is implemented for all the flow properties with respective free stream values. At any time, the far-field boundary condition can be either of the following: a subsonic inflow or a supersonic inflow; a subsonic outflow or a supersonic outflow boundary condition. Inflow and outflow boundary conditions are applied at the boundary, depending upon the sign of normal component of velocity (Vn). Here, we have used an orthogonal grid, and this velocity component is the same as the contravariant component of velocity (gxu + gyv) in the transformed plane. Additionally, the far-field boundary conditions are determined based on eigenvalues of the locally one-dimensional description of the flow, details of which can be found in [26,27]. 3. Numerical procedure A computational domain with outer boundary segment located at approximately 10c from the airfoil is used, to ensure accurate implementation of far-field boundary condition and to reduce the acoustic wave reflections from the boundary. An orthogonal, body-conforming O-grid is generated using hyperbolic grid generation technique [14,28]. For this purpose, one of the two Beltrami equations given below are solved in a finite difference framework

xg ¼ fyn ; yg ¼ fxn along with the orthogonality condition

xn xg þ y n y g ¼ 0 Here f represents the distortion function, which is given in terms of grid scale factors h11 and h22 as

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2g þ y2g h22 f ¼ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h11 x2 þ y 2 n

n

The distance between successive g = constant lines in the wallnormal direction (DSg = h22Dg), is prescribed by the following distribution,



tanh ½bð1  gÞ h22 ¼ H 1  tanh b where b is the clustering parameter and H is the nondimensional distance of the outer boundary in terms of c. In the present work, high accuracy optimized dispersion relation preserving (DRP) upwind schemes and their symmetrized versions, OUCS and S-OUCS schemes [29], are used for evaluating the convective flux terms in g- and n-directions, respectively. The OUCS and S-OUCS schemes are implicit and have superior stability and DRP properties [12,13,30]. These schemes are obtained by optimizing dissipation and dispersion error in spectral plane [14], thereby ensuring higher resolution. These schemes are used in parallel computing framework by domain decomposition technique [16] to solve the problem in smaller domains independently with message passing interface (MPI) used for parallelizing the code. The interior stencil of the optimal upwind OUCS2 scheme [12] 0 used here to evaluate first derivative (d ), is given by 0

0

0

b1 dj1 þ b2 dj þ b3 djþ1 ¼

2 1X ak djþk h k¼2

ð3Þ

a1 a1 b2 a1 where, b1 ¼ b32  12 ; b3 ¼ b32 þ 12 ; a2 ¼  36 þ 72 ; a1 ¼ 7b92 þ a91 and 6 a0 ¼  a41 . The free parameter a1 is the coefficient of @@x6d , which is treated as the free parameter to fix the order of representation [14]. The interior stencil of the OUCS2 is a fifth order accurate upwind scheme for a1 < 0 and a sixth-order accurate central scheme results for a1 = 0, as can be shown by Taylor’s series expansion. The other compact scheme used here is the OUCS3 scheme [12], whose interior stencil is given by

24

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

0

0

0

mj1 dj1 þ dj þ mjþ1 djþ1 ¼

2 1X q djþk h k¼2 k

ð4Þ

f f f 11f where, mj1 ¼ D  60 ; q2 ¼  4F þ 300 ; q1 ¼  2E þ 30 and q0 ¼  150 with D = 0.3793894912;E = 1.57557379 and F = 0.183205192. The 4  free parameter f is the coefficient of @@x4d which is treated as the

Fig. 8. M contours for flow around SHM-1 airfoil with M1 = 0.62, Re1 = 13.6  106 at a = 0.27° computed with scheme used for evaluating convective fluxes: (a) OUCS2 (b) OUCS3.

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

free parameter to fix the order of representation [14]. The interior stencil of OUCS3 scheme is second-order accurate. For both the OUCS2 and OUCS3 schemes, following explicit closure schemes for the boundary and near boundary points have been used with the stencils proposed as [12] 0

d1 ¼ 0

d2 ¼

1 h

3d1 þ 4d2  d3 2h      

2b2 1 8b2 1 8b2 1 2b  d1  þ d2 þ ð4b2 þ 1Þd3  þ d4 þ 2 d5 3 2 6 3 3 3 3

ð5Þ ð6Þ

25

Similar closures are written for j = N and (N  1), with the coefficients on the right hand side changing signs. The value of bj is chosen as 0.02 for j = 2 and 0.09 for j = (N  1), while using Eq. (6) with OUCS2 interior stencil. These values are obtained for optimal behavior in terms of dispersion error, by studying global properties of Eq. (3) along with Eq. (5) and (6) in the spectral plane [12,14]. Similarly, value of bj in Eq. (6) is chosen as 0.025 for j = 2 and 0.09 for j = (N  1), while using with OUCS3 interior stencil. While implementing the symmetrized versions of respective schemes (S-OUCS), the derivatives at all the nodes are calculated using OUCS schemes by traveling from j = 1 to N and then again

Fig. 9. q contours for flow around SHM-1 airfoil with M1 = 0.62, Re1 = 13.6  106 at a = 0.27° computed with scheme used for evaluating convective fluxes: (a) OUCS2 and (b) OUCS3.

26

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

going from N to 1. Symmetrized derivative is taken as the average of the derivatives obtained by following these two directions. By this procedure, symmetrized scheme removes inherent bias due to upwinding, which is present in all OUCS schemes [29]. Procedure for full domain analysis of any numerical scheme for stability, DRP properties, phase error and the error evolution analysis is provided in [29,31]. Artificial diffusion terms due to Jameson, Schmidt and Turkel (JST) [32] are added explicitly to the numerically evaluated convective flux terms, in order to damp higher wavenumber components which arise due to non-linearities and discontinuities in the

solution and to accurately capture discontinuities. The artificial diffusion terms (Dn) and (Dg) added in n- and g-directions, respectively are given by

    Dn ¼ rn J 1 dn2  dn4

ð7Þ

    Dg ¼ rg J 1 dg2  dg4

ð8Þ

where, rn and rg are the spectral radii of the flux Jacobian matrices [26]. These terms are blend of second and fourth diffusion terms, dn2 , dn4 , dg2 , and dg4 added selectively at various points of the computational domain and are defined as

Fig. 10. p Contours for flow around SHM-1 airfoil with M1 = 0.62, Re1 = 13.6  106 at a = 0.27° computed with scheme used for evaluating convective fluxes: (a) OUCS2 and (b) OUCS3.

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

dn2 ¼ 2n rn Dn ðJQ Þ

2 4 2g 4g n

2

dn4 ¼ 4n ðrn Dn Þ ðJQ Þ

n

dg2 ¼ 2g rg Dg ðJQ Þ  2 dg4 ¼ 4g rg Dg ðJQ Þ with 2n , 4n , 2g , and tional domain as

4g

are defined at points i, j in the computa-

27

  ¼ j2 max  iþ1;j ;  i;j ;  i1;j   ¼ j4 max 0; j4  2n   ¼ j2 max  i;jþ1 ;  i;j ;  i;j1   ¼ j4 max 0; j4  2g

where j2 and j4 are the parameters which controls the amount of numerical diffusion to be added. Values of j2 = 0.25 and j4 = 0.15

Fig. 11. Computed M contours for flow around NACA 0012 airfoil with M1 = 0.779, Re1 = 3  106 at a = 0.14°.

28

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

have been used. The terms i,j, are the pressure based switch defined at the point (i, j) as [32]

 i;j ¼

jpiþ1;j  2pi;j þ pi1;j j jpiþ1;j j þ j2pi;j j þ jpi1;j j

ð9Þ

The term i,j has very high value near a shock. Hence, it acts as a pressure switch activating the second difference term, wherever there are discontinuities in pressure. The fourth difference term is

switched to zero near shock. These artificial diffusion terms are added to the numerically evaluated convective flux derivative terms for numerical stability and better dispersion properties. Viscous flux terms represent diffusion phenomenon, which is isotropic. Hence, the derivatives are calculated using second order central discretization (CD2) in self-adjoint form, preserving isotropic nature. Classical RK4 scheme is used for time integration [32]. The parallel computing is performed by domain decomposition technique, developed using above mentioned numerical methods,

Fig. 12. Computed p contours for flow around NACA 0012 airfoil with M1 = 0.779, Re1 = 3  106 at a = 0.14°.

29

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

time = 9.00

time = 14.00

-1

-1

-0.5

-0.5

0

0

0.5

0.5

1

1 0

0.2

0.4

x/c

0.6

0.8

1

0

0.2

0.4

0.2

0.4

0.2

0.4

0.2

0.4

x/c

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

time = 15.00

time = 10.00 -1

-1

-0.5

-0.5

0

0

0.5

0.5 1

1 0

0.2

0.4

0.6

0.8

0

1

x/c

x/c

time = 16.00

time = 11.00 -1

-1

-0.5

-0.5

0

0

0.5

0.5 1

1 0

0.2

0.4

x/c

0.6

0.8

0

1

x/c

time = 17.00

time = 12.00 -1

-1

-0.5

-0.5

0

0

0.5

0.5 1

1 0

0.2

0.4

x/c

0.6

0.8

0

1

x/c

time = 18.00

time = 13.00 -1

-1

-0.5

-0.5

0

0

0.5

0.5 1

1 0

0.2

0.4

x/c

0.6

0.8

1

0

0.2

0.4 x/c

Fig. 13. Comparison of computed instantaneous Cp-distribution with the experimental values for flow around NACA 0012 airfoil with M1 = 0.779, Re1 = 3  106 at a = 0.14°.

30

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

Fig. 14. Computed s contours for flow around NACA 0012 airfoil with M1 = 0.779, Re1 = 3  106 at a = 0.14°.

which require overlap of subdomains. Motivation behind the overlap is to overcome problems due to the nonphysical growth and attenuation of the solution near the interfaces of neighboring subdomains, due to the effects of boundary closures [16] for compact schemes. It is established in [33] that one requires sufficiently large overlap of subdomains to minimize or eliminate the problems of convergence. Selective addition of artificial second and fourth order diffusions [32] help to use smaller number of overlapping points [16]. Here eight overlapping points have been used.

4. Results and discussion Here, we report results for transonic flows past NACA 0012 and SHM-1 airfoils displaying mixed flow nature, i.e., there are regions where the flow speed exceeds sonic speed locally, while in other parts the flow Mach number remains subsonic. As a consequence, one notices discontinuities in the physical variables, leading to formation of shock(s) and contact discontinuities. These flows are inherently unsteady and thus, the usage of unsteady formulation

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

-0.8 -0.6 -0.4

Cp

-0.2 0 0.2 0.4 0.6 COMPUTATIONAL

0.8

EXPERIMENTAL

1 0

0.2

0.4

0.6

0.8

1

x/c Fig. 15. Comparison of computed and experimental Cp-distribution for flow around NACA 0012 airfoil with Re1 = 3  106 at a =  0.14° and M1 = 0.82. Computed data is time-averaged over t = 10–20.

is absolutely necessary. In comparing with experimental results, we will time-average the computed solution taking meaningful time intervals, with justification provided. 4.1. Validation of methodologies for compressible flow calculations and shock capturing For computing transonic flow past NACA 0012 airfoil, we have used 636 points in the azimuthal and 350 points in the wall-normal directions in O-grid topology, for the solution of compressible Navier–Stokes equation for Re1 = 3  106 and three Mach numbers, M1 = 0.6, 0.758, and 0.779. For these flow parameters experimental results exist with which numerical results are validated. The close-up view of grids used for NACA 0012 and SHM-1 airfoils are shown in Fig. 2(a) and (b), respectively. As noted in Section 3, the employed grid generation technique is accurate with desirable numerical properties of orthogonality and various metrics have smooth variations in the domain. For the SHM-1 airfoil, 663 points have been taken in n- and 450 points in g-directions, respectively. This is an airfoil with higher (15%) thickness, yet it has better performance at cruise Mach number (M1 = 0.69), as compared to NACA 0012 airfoil, as can be seen by comparing the drag polar of SHM-1 airfoil in [34] with results from literature on NACA 0012 airfoil. For SHM-1 airfoil, we validate our results for M1 = 0.62 and Re1 = 13.6  106 at an angle of attack, a = 0.27°. The near-wall resolution for NACA 0012 and SHM-1 airfoils used are given by Dsg = 1.89  105c and Dsg = 8.20  106c, respectively. There are many reasons for choosing such values of wall-normal spacing with the used schemes. Employed compact schemes [12–14] have near-spectral accuracy and these values have been found more than adequate, like in any other global methods. Secondly, the flow is not turbulent right from the leading edge of the airfoil for both the cases considered. Specifically, for the NLF airfoil the flow remains laminar up to and beyond the mid-chord, as per the design requirement. For example, at the locations where the flow pressure gradient is zero, these spacing are obtained in wall unit as, Dy+jw = 0.1669 for the NACA 0012 airfoil and Dy+jw = 0.1312 for the SHM1 airfoil. We also emphasize that for downstream locations past the shock location, wall friction value decreases due to flow separation and there Dy+jw values reduce even further for both the airfoils. Such behavior of wall friction with streamwise

31

distance is also reported in [35] and is due to shock boundary layer interaction. Thus the obtained results truly show adequacy of these choices. All the computations performed here using OUCS3-RK4 scheme (unless otherwise stated) are with a small time step of Dt = 106, which is necessary for better DRP properties. However, the main reason for the choice rests on the requirement of neutral stability for the chosen space–time discretization scheme. For the detailed analysis of error propagation one is referred to [13], where it is shown that for DNS there is no other recourse but to take neutrally stable scheme. Such requirement also restricts time step for explicit time integration schemes, such as the one used here. Numerical properties of the chosen OUCS3-RK4 discretization, the numerical amplification factor, phase speed and group velocities are also provided in [29,31] and not repeated here. We emphasize that these numerical schemes have been generated with the purpose of having excellent DRP properties, measured by calibrating the methods in capturing discontinuities, without dissipation and dispersion errors. Thus, the developed methods produce accurate results for incompressible and compressible flows, without requiring any flow-specific fixes. Presented results would testify this assertion that one does not require different numerical methods for compressible and incompressible flows, while utilizing near-spectral resolution of compact schemes, with excellent DRP properties. In Fig. 3(a), comparison between computed and experimental [24] pressure coefficient (Cp) distribution is shown for flow over NACA 0012 airfoil for M1 = 0.6,Re1 = 3  106 and a = 0.14°. The experimental correction for angle of attack (Da = 1.55CL) is implemented in the code, which is due to wall-induced lift interference effects [24], where CL denotes the sectional normal-force coefficient. Fig. 3(b) shows comparison between computed and experimental [24] Cp distribution for flow over NACA 0012 airfoil for higher Mach number case of M1 = 0.758, Re1 = 3  106, and a = 0.14°. One notices very good match for the increased M1 case, between numerical and experimental Cp distribution. To further validate the performance of the developed method in capturing the shock, the solver is used for simulating flow over NACA 0012 airfoil for M1 = 0.779 and corresponding results are shown in Fig. 4. The results once again shows good agreement between computed and experimental [24] Cp-distribution for the same Reynolds number and angle of attack. In Fig. 4, one notices that the shock has been captured with good accuracy, with the pre-shock dip in Cp noted upstream of the shock. In performing the reported calculations here, we have not used any models for transition and/or turbulence. The near-spectral accuracy of compact schemes with good dispersion properties enables one to capture shock with good resolution. To understand the distinctive feature of the developed methodology to simultaneously capture the laminar, transitional and turbulent flows, we show the velocity profile of the computed flow past NACA 0012 airfoil in Fig. 5 in wall unit, for the case of Fig. 4. Such a representation is preferred for turbulent flows [36], where one notices the presence of a viscous sublayer at the wall and this is followed by the inertial log law region. A comprehensive description of the flow profile is given by the Spalding law [36]. In Fig. 5, one can clearly notice pairs of profiles representing laminar, transitional and turbulent flows for this Reynolds number of 3  106. We note that the flow only becomes turbulent near the trailing edge with the distinct match shown by the plotted time-averaged profile with that given by the Spalding law. To further elucidate the fact that selectively added diffusion does not affect the solution accuracy, apart from controlling numerical instability and dispersion. It is the magnitudes of added terms proportional to j2 and j4 which are important and these are shown in Fig. 6, for the case of flow past NACA 0012 airfoil shown in Fig. 4. In each frame, the maximum and minimum values of the

32

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

Fig. 16. Computed M contours for flow around NACA 0012 airfoil with M1 = 0.82, Re1 = 3  106 at a = 0.14°.

quantities involved are also marked. It is readily evident that the second diffusion term is only added near the shock and only a limited region of the flow field is affected by these added terms, showing the reason why it should be viewed as performing DNS. To further test capabilities of the developed procedure in solving compressible flows, we report results for computed flow over natural laminar flow (NLF) SHM-1 airfoil for M1 = 0.62 and Re1 = 13.6  106 at a = 0.27° in Fig. 7. These simulations are performed with convective flux derivatives evaluated using both OUCS2 and OUCS3 schemes and results compared in Fig. 7, with experimental and computed Cp-distribution [34]. Although, the

OUCS2 scheme has fifth order accuracy, as compared to second order accuracy of the OUCS3 scheme, identical results are obtained with both the schemes in capturing flow features correctly. This is due to both the schemes having comparable spectral resolution [13]. Figs. 8(a), Figs. 9(a) and 10(a) show computed M,p and q-contours, respectively, at different time instants for flow around SHM1 airfoil using the OUCS2 scheme, while Figs. 8(b), 9(b) and 10(b) show the corresponding computed M, p, and q-contours computed with the OUCS3 scheme at the same time instants. One can notice identical flow features captured by both the schemes. For example, the Mach contours are positioned identically for these two

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

33

Fig. 17. Computed p contours for flow around NACA 0012 airfoil with M1 = 0.82, Re1 = 3  106 at a = 0.14°.

methods. Similar observations can be made with respect to pressure and density contours, as obtained by both these methods. 4.2. Computing strong shock cases We have established the correctness of the adopted procedures with respect to Cp-distribution for both sub- and super-critical

cases. It is worth noting that such accuracy of results are demonstrated without the use of any models for transition and turbulence for compressible flows. Validation of numerical with experimental results, not only establishes the credibility of the adopted DRP compact schemes, but also establishes the consistency of compact schemes by comparing detailed results in the domain by using a fifth order (OUCS2) and a second order (OUCS3) method. Having

34

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

Fig. 18. Computed q contours for flow around NACA 0012 airfoil with M1 = 0.82, Re1 = 3  106 at a = 0.14°.

established accuracy and consistency, in the following we compute cases for higher Mach numbers representing transonic flows with shock waves creating entropy. In Figs. 11 and 12, computed M and p-contours are shown, respectively, at different time instants for the flow around NACA 0012 airfoil for M1 = 0.779, Re1 = 3  106 and at a = 0.14°. The flow encounters a favorable pressure gradient region on the top

surface near the leading edge of the airfoil, up to approximately 30% of the airfoil chord (this location is referred to as xcr, a location downstream of the shock formation region). The flow accelerates creating a pocket of supersonic flow over the top surface. Beyond xcr, the flow tends to recover pressure and flow retards as a result of adverse pressure gradient. This is observed as a sharp rise in Cp, as shown in Fig. 4 near x = 0.30c. Pressure contours in Fig. 12

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

35

Fig. 19. (a) Computed M contours showing subsonic, sonic and supersonic region. (b) Computed s contours showing entropy generation in the region of the shock for M1 = 0.82 and Re1 = 3  106 at a = 0.14°.

indicate the shock formation, which is also noted in Fig. 10 from the Mach contours. 4.3. Unsteady flow behavior of compressible flows We have already noted in Figs. 8–10, that even when shock waves have not formed, the flow over the airfoil displays time dependent behavior, as clearly captured in the contour plots of

physical variables. When strong shocks form, such unsteadiness increases due to physical dispersion of higher wavenumber components. For a DRP scheme, such convection of high wavenumbers are naturally captured, as noted in Fig. 13, for flow past NACA 0012 airfoil for M1 = 0.779, where snapshots of the Cp-distribution on the airfoil surface is shown at the indicated time instants. For this flow, from the near vicinity of the shock, one notices convection of pressure pulses in the downstream direction. It is associated

36

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

Table 1 Computed and experimental [24] values of CL and CD for viscous flow past NACA 0012 airfoil. Sr. no.

1 2 3

M1

Re1

0.6 0.779 0.82

3  106 3  106 3  106

a (in °)

0.14 0.14 0.14

Computed

Experimental

CL

CD

CL

CD

0.0161 0.00989 0.0406

0.0022 0.0057 0.0238

0.016 0.013 NA

0.006 0.007 NA

with the sharp pressure discontinuity which is responsible for wide-spectrum of physical variables. However, we want to emphasize that such unsteadiness for high wavenumbers is not numerical. It is well known that discontinuities are sites of q-waves (spurious upstream propagating waves due to extreme dispersion), as explained in [37] and in Fig. 13, no such upstream propagating waves are noted. This type of unsteady downstream propagating waves are attributes of adverse pressure gradient in the flow, due to creation of inflection point in the velocity profile. The time average of Cp-distribution shows very good match with experimental value, indicating the fact that such pressure measurements are not instantaneous, but averaged over a finite time interval. Similar kind of Cp-distribution results are also reported in [38], which show unsteadiness in instantaneous values, whereas time-averaged values show match with experimental values. For the time accurate solution of Navier–Stokes equation, definitive initial conditions are not known. Thus, a uniform flow condition is assumed at t = 0 and the results at early time are not relevant physically. In most of the cases, it takes about t ’ 10 for the flow to develop into a meaningful viscous flow state. Hence, in the reported cases, we show time averages taken from time series from t = 10 to t = 45. 4.4. Creation of rotational effects In the case of NACA 0012 airfoil at small angles of attack, two different streams of flow are created which originate from the top and bottom surface of the airfoil and reach the trailing edge with different values of entropy. Such dissimilar flows will create a mixing layer, which is unstable and roll into vortices of different length scales, as shown in Fig. 14, where entropy contours are plotted, with the values calculated from

s  s0 ¼ cv ln

p

qc

The fact that entropy remains constant for any flow which has uniform total enthalpy, also indicates that such flows are irrotational. Whenever vorticity is created, as near the trailing edge for the mixing layer due to flow instability, the flow creates an entropy gradient across any streamline, which is a consequence of Crocco’s theorem given by

~¼ x

T dS u dn

The relationship between vorticity, total enthalpy and entropy gradient is given by Crocco’s theorem [39] given as

! ! @V ~ ¼ rh0 þ T rS þ V  x @t 4.5. Strong shock case and creation of entropy gradient We also report here, results obtained by numerical simulation performed for a case of M1 = 0.82, Re1 = 3  106 and a = 0.14° for flow over NACA 0012 airfoil. The value of M1 here is much higher than the drag divergence Mach number (MDD) for NACA

0012 airfoil. For this case the flow gradient is expected to be much higher and an orthogonal grid around NACA 0012 airfoil is generated with 2658 points in n- and 360 points in g-directions, respectively, to capture such gradients. Near-wall resolution of the grid is given by Dsg = 8.70  106c. Very strong shocks are observed on upper and lower surfaces of the airfoil for this case. In Fig. 15, computed Cp-distribution is compared with the experimental value [24], showing excellent match with the sharp drop in the Cp-distribution near the shock region. Contours of computed M and p in the vicinity of the airfoil for this case are shown in Figs. 16 and 17, respectively, for same time instants. Since the oncoming flow has higher value of M1, the flow accelerates near the leading edge to higher value of M than that observed in the case of M1 = 0.779. For some time instants, as at t = 15.0, the maximum M reaches a high value of 1.6872. The computed q-contours at the corresponding times are presented in Fig. 18. In Fig. 19(a) and (b), M and scontours are plotted, respectively, for the case which display the presence of very strong shock. Moreover, the strength of the shock is evident from the creation of additional entropy, aligned with the location of the shock. As usual, one notices the creation of entropy in the shear layer at the aft region of the airfoil and in the near wake. 4.6. Lift and drag calculation The coefficient of lift (CL) and drag (CD) are computed by integrating the static pressure and shear stress over the airfoil surface. Computed results are compared with available experimental values. For flow over NACA 0012 airfoil with M1 = 0.6 and Re1 = 3  106 at a = 0.14° experimental values of CL and CD are provided in [24] as 0.016 and 0.006, respectively. The computed values of CL and CD for the same case are noted as 0.0161 and 0.0022, respectively, if the data are time-averaged in the range 10 6 t 6 55. For the flow around NACA 0012 airfoil for M1 = 0.779, Re1 = 3  106 and a = 0.14° experimental values of CL and CD are provided in [24] as 0.013 and 0.007, respectively. The present computed values of CL and CD for the same case are noted as 0.00989 and 0.0057, respectively, if time-averaged in the range 10 6 t 6 45. The time-averaged CL and CD values obtained by computations for the flow around NACA 0012 airfoil for M1 = 0.82; Re1 = 3  106 and a = 0.14° are 0.0406 and 0.0238, respectively, for 10 6 t 6 25. For this case, experimental value for CL is given as 0.008 in [24]. All these computed values are consolidated in Table 1 for ease of understanding. 5. Conclusion High-resolution, high accuracy 2D compressible Navier–Stokes solvers without any transition and/or turbulence models are developed and applied successfully to study the transonic flow past NACA 0012 and SHM-1 airfoils. It is observed that, the transonic flow phenomena are captured correctly, including shock formation, post drag divergence unsteady behavior. The computed pressure distributions are noted to be in good agreement with the experimental results available in [24,34]. Although the compact schemes used here were originally developed for direct simulation of incompressible flows, it is shown here that same schemes can directly be used for compressible flows. Moreover, presented computations do not require any models for transition and turbulence, except the JST terms adding adaptively second and fourth order diffusion at selective points with the help of pressure switch, as was originally proposed by Jameson et al. [32]. Hence, these compact schemes can be used to study phenomena that occur in transonic flows such as shocks, drag divergence and shock-boundary layer interaction.

T.K. Sengupta et al. / Computers & Fluids 88 (2013) 19–37

References [1] Moulden TH. Fundamentals of transonic flow. Canada: Wiley; 1984. [2] Maybey DG. Physical phenomena associated with unsteady transonic flows. In: Unsteady transonic aerodynamics, progress in astronautics and aeronautics, AIAA series. Washington, DC, USA: AIAA; 1989. [3] Jameson A, Ou K. 50 Years of transonic aircraft design. Prog Aero Sci 2011;47(5):308–18. [4] Visbal MR. Dynamic stall of a constant-rate pitching airfoil. J Aircraft 1990;27:400–7. [5] Garbaruk A, Shur M, Strelets M, Spalart P. Numerical study of wind tunnel wall effects on transonic airfoil flow. AIAA J 2003;41:1046–54. [6] Alshabu A, Olivier H, Klioutchnikov I. Investigation of upstream moving pressure waves on a supercritical airfoil. Aero Sci Tech 2006;10:465–73. [7] Allaneu Y, Jameson A. Direct numerical simulations of a two-dimensional viscous flow in a shocktube using a kinetic energy preserving scheme. AIAA Paper 2009-3797; 2009. [8] Chiu EK, Wang Q, Jameson A. A conservative meshless scheme: general order formulation and application to Euler equations. AIAA Paper 2011-651; 2011. [9] Beam RM, Warming RF. An implicit factored scheme for the compressible Navier–Stokes equations. AIAA J 1978;16(4):393–402. [10] Chakravarthy S, Harten A, Osher S. Essentially non-oscillatory shock capturing schemes of uniformly very high accuracy. AIAA Paper 86-0339; 1986. [11] Liu TCX, Osher S. Weighted essentially non-oscillatory schemes. J Comp Phys 1994;115:200–12. [12] Sengupta TK, Ganeriwal G, De S. Analysis of central and upwind compact schemes. J Comp Phys 2003;192:677–94. [13] Sengupta TK. High accuracy computing methods: fluid flows and wave phenomena. UK: Cambridge Univ. Press; 2013. [14] Sengupta TK. Fundamentals of computational fluid dynamics. Hyderabad, India: Universities Press; 2004. [15] Sengupta TK, Bhaumik S, Bhumkar Y. Direct numerical simulation of twodimensional wall-bounded turbulent flows from receptivity stage. Phys Rev E 2012;85(2):026308. [16] Sengupta TK, Dipankar A, Kameswara Rao A. A new compact scheme for parallel computing using domain decomposition. J Comp Phys 2007;220:654–77. [17] Binion TW. Limitations of available data. AGARD-AR-138, May 1979. [18] Garner HC, Rogers E, Acum W, Maskell E. Subsonic wind tunnel wall corrections. AGARDograph AG-109; 1966. [19] Jacocks JL. An investigation of the aerodynamic characteristics of ventilated test section walls for transonic wind tunnel. Ph.D. dissertation. Univ. Tennessee; 1976. [20] Maybey DG. Some remarks on the design of transonic tunnels with low levels of flow unsteadiness. NASA CR-2722; 1976.

37

[21] Pate SR, Schueler CJ. Radiated aerodynamic noise effects on boundary layer transition in supersonic and hypersonic wind tunnels. AIAA J 1968;7:450–7. [22] Otto H. Systematical investigation of the influence of wind tunnel turbulence on the results of force measurements. AGARD CP-174; 1976. [23] Mercer JE, Geller EW, Johnson ML, Jameson A. Transonic flow calculations for a wing in a wind tunnel. AIAA J 1981;18(9):707–11. [24] Harris CD. Two-dimensional aerodynamic characteristics of the NACA0012 airfoil in the Langley 8-foot transonic pressure tunnel. NASA TM-81927; 1981. [25] Barnwell RW. A similarity rule for sidewall-boundary-layer effect in twodimensional wind tunnels. AIAA Paper 79-0108; 1979. [26] Hoffmann KA, Chiang ST. Computational fluid dynamics, Engineering education system, vol. II. Wichita, Kansas, USA; 2000. [27] Sreejith NA. Effect of plasma for transonic flow past an airfoil. M. Tech. Thesis, Dept. Aero. Engg. Indian Institute of Technology, Kanpur, India; 2011. [28] Nair MT, Sengupta TK. Orthogonal grid generation for Navier–Stokes computations. Int J Num Meth Fluids 1998;28:215–24. [29] Dipankar A, Sengupta TK. Symmetrized compact schemes for receptivity study of 2D channel flow. J Comp Phys 2006;215:245–73. [30] Sengupta TK, Ganeriwal G, Dipankar A. High accuracy compact schemes and Gibbs phenomena. J Sci Comput 2004;21(3):253–68. [31] Sengupta TK, Dipankar A, Sagaut P. Error dynamics: beyond von Neumann analysis. J Comp Phys 2007;226:1211–8. [32] Jameson A, Schmidt W, Turkel E. Numerical simulation of the Euler equations by finite volume methods using Runge–Kutta time-stepping schemes. AIAA Paper 1981-1259. AIAA 14th fluid and plasma dynamics conference. Palo Alto, California; 1981. [33] Dolean V, Lanteri S, Nataf FC. Optimized interface conditions for domain decomposition methods in fluid dynamics. Int J Num Meth Fluids 2002;40:1539–50. [34] Fujino M, Yoshizaki Y, Kawamura Y. Natural-laminar-flow airfoil development for a lightweight business jet. J Aircraft 2003;40(4):609–15. [35] Pirozzoli S, Bernardini M, Grasso F. Direct numerical simulation of transonic shock/boundary layer interaction under conditions of incipient separation. J Fluid Mech 2010;57:361–93. [36] Panton RL. Incompressible flow. USA: John Wiley & Sons Inc.; 2005. [37] Sengupta TK, Bhumkar Y, Rajpoot M, Saurabh S. Purious waves in discrete computation of wave phenomena and flow problems. App Math Comp 2012;218(18):9035–65. [38] Hermes V, Klioutchnikov I, Olivier H. Numerical investigation of unsteady wave phenomena for transonic airfoil flow. Aero Sci Tech 2013;25:224–33. [39] Hirsch C. Numerical computation of internal and external flows, Fundamentals of numerical discretization, vol. 1. New York, USA: Wiley-Interscience Publication; 1994.