Calculation of transonic flows around profiles with blunt and angled leading edges

Vietnam Journal of Mechanics, VAST, Vol. 38, No. 1 (2016), pp. 1 – 13 DOI:10.15625/0866-7136/38/1/4177 CALCULATION OF TRANSONIC FLOWS AROUND PROFILES WITH BLUNT AND ANGLED LEADING EDGES Hoang Thi Bich Ngoc∗, Nguyen Manh Hung 1Hanoi University of Science and Technology, Vietnam ∗E-mail: hoangthibichngoc@yahoo.com Received December 07, 2014 Abstract. Transonic flow is a mixed flow of subsonic and supersonic regions. Because of this mixture, the solution of transonic flow problems is obtained only when solving the differential equations of motion with special treatments for the transition from subsonic re- gion to supersonic region and vice versa. We built codes solving the full potential equation and Euler equations by applying the finite difference method and finite volume method, and also associated with software Fluent to consider the viscous effects. The analysis of results calculated for cases of transonic flow over profiles with blunt and angled leading edges shows more clearly the physical nature of the gas - solid interaction at leading edges in the mixed flow and the optimal application of each profile in transonic flows. Keywords: Transonic flow, finite volumes, finite differences, blunt LE, angled LE. 1. INTRODUCTION In order to solve transonic flows (with free flow Mach numbers M∞ ≥ 0.7), it is necessary to use the equations of compressible flow. In the assumption of potential flow, differential equations of compressible flow are Euler equations and full potential equation. For incompressible flows (with Mach numbers M∞ < 0.3), by considering the constant density, the Euler equations and the full potential equation are reduced to the Laplace equation of potential (elliptic form). Whatever the method, the transonic prob- lem needs to treat the transition from subsonic flow zones to supersonic flow zones and vice versa. Motion differential equations in this transition zone change back and forth from elliptic form to hyperbolic form. The appearance of a supersonic region can cause shock waves. The method of Euler equations permits to calculate cases with relatively strong shock waves, but the algorithm is complex, the computer memory is required large and the running time is long. The method of full potential equation (FPE) allows calculating cases with not too strong shock waves with a maximal local Mach number M < 1.5. The FPE method saves the memory and running time. In fact, when flow hav- ing strong shock wave, it may appears a strong interaction between the shock wave and c© 2016 Vietnam Academy of Science and Technology 2 Hoang Thi Bich Ngoc, Nguyen Manh Hung the boundary layer. That makes a relatively large difference of results calculated by the viscous flow theory and the inviscid flow theory. In these cases, it is necessary to use the viscous method. Viscous flow calculations have be done here by using Fluent software with stringent verification operations to ensure accuracy in the application domain. 2. METHODS OF RESOLUTION 2.1. Euler equations and method of solving Differential equations of inviscid flow in the general case are Euler’s equation [1]. ∂ ∂t ∣∣∣∣∣∣∣∣ ρ ρu ρv ρE ∣∣∣∣∣∣∣∣+ ∂ ∂x ∣∣∣∣∣∣∣∣ ρu ρu2 + p ρuv ρuH ∣∣∣∣∣∣∣∣+ ∂ ∂y ∣∣∣∣∣∣∣∣ ρv ρuv ρv2 + p ρvH ∣∣∣∣∣∣∣∣ = ∣∣∣∣∣∣∣∣ 0 fex fey W f ∣∣∣∣∣∣∣∣ , (1) where u, v are velocity components; fex, fey are external force components; ρ is the den- sity; p is pressure, p = (γ− 1)ρ [E− (u2 + v2)/2]; E and H are the internal energy and enthalpy H = E + p/ρ; E = c2/γ(γ− 1) + (u2 + v2)/2 with c is the speed of sound and γ is the ratio of specific heats. For subsonic flows, entropy remains constant and uniform in the whole flow do- main. The numerical dissipation generated by the numerical scheme will mimic in some way the physical dissipation. Seen from the aspect of the discretization, the computer cannot distinguish physical dissipation from numerical dissipation. However for tran- sonic flows, the transition from subsonic flow regions to supersonic flow regions and vice versa, and the increase of entropy require suitable numerical schemes in order to avoid differences caused by numerical dissipation. When physical dissipations in the left-hand side and in the right-hand side at a point are different, it is necessary to use a scheme depending on time in the steady flow. Therefore, time integration methods for space discretized equations are used in this case. The relationship between the time step and the space interval is expressed by the Courant number. We can write Eq. (1) under the following general form ∂U ∂t + ∂ f ∂x + ∂g ∂y = Q. (2) The use of coupled time scheme allows treating a discontinue transition from left to right, when one side the flow is subsonic and the other side it is supersonic, and vice versa [2]. The integral time is done by the method of Runge-Kutta 4. Eqs. (1) and (2) are solving by a finite volume method. Conservative integral equation for a volume element Ωj is written ∫ Ωj ∂U ∂t dΩ+ ∫ Ωj ~∇.~FdΩ = ∫ Ωj QdΩ, (3) where ~F is a vector of components f and g. Calculation of transonic flows around profiles with blunt and angled leading edges 3 (a) (b) Fig. 1. a) Cell vertex (finite volume mesh); b) Generation of triangular mesh for Naca 4412 Consider a mesh having triangular cells (Fig. 1). Using Green’s formula for the domain Ωj leads to the following discretized scheme Ωj ∂Uj ∂t + 1 2 6 ∑ side=1 [( fk + fk+1)(yk+1 − yk)− (gk + gk+1)(xk+1 − xk)] = 0, where k and k + 1 are the starting and ending points of each side and the point k = 7 coincides with k = 1. The summation extends over all the nodes by assembling terms such as ( f1 + f2)(y2 − y1) + ( f2 + f3)(y3 − y2) + ... = f2(y3 − y1) + ...− y2( f3 − f1) + ... we obtain the alternative formulation Ωj ∂Uj ∂t + 1 2 6 ∑ k=1 [ fk(yk+1 − yk−1)− gk(xk+1 − xk−1)] = 0. (4) With a scheme of 6 triangles associated to point j: y0 = y6; y7 = y1; x0 = x6; x7 = x1 (and similar for f , g). The index ′ j′ is for the control volume and the index ′k′ is for the triangle next to the control triangle. f jk and gjk are determined by the scheme of flux splitting. The first order upwind scheme for a numerical flux is written: f ∗(1)i+1/2 = f ∗(Ui, Ui+1), and the second order upwind scheme is: f ∗(2) i+1/2 = f ∗(ULi+1/2, U R i+1/2). ULi+1/2 = Ui + (1−δ)(Ui−Ui−1) + (1+δ)(Ui+1−Ui) 4 , URi+1/2 = Ui+1 − (1+δ)(Ui+1−Ui) + (1−δ)(Ui+2−Ui+1) 4 , (5) 4 Hoang Thi Bich Ngoc, Nguyen Manh Hung where δ indicates the wave property, getting the values 0, 1,−1. The indexes ′L′ and ′R′ are for the left and right sides at the considering position. In summary, the second order linear scheme with time (split flux scheme) for flow having different rules in left and right sides are written as follows f ∗(2)i+1/2,j = f +(ULi+1/2,j) + f −(URi+1/2,j), g∗(2)i,j+1/2 = g +(ULi,j+1/2) + g −(URi,j+1/2). (6) Four steps of the Runge-Kutta method with time are written U(1)j = U n j − ∆t α1Rnj Ωj , U(2)j = U n j − ∆t α2R (1) j Ωj , U(3)j = U n j − ∆t α3R (2) j Ωj , U(n+1)j = U n j − ∆t α4R (3) j Ωj , (7) where α1 = 1/4; α1 = 1/6; α3 = 3/8; α4 = 1. 2.2. Full potential equation and method of solving The equation for the full potential Φ is written as [3]( a2 −Φ2x ) Φxx + ( a2 −Φ2y ) Φyy − 2ΦxΦyΦxy = 0, (8) where a is speed of sound; derivatives Φx = ∂Φ/∂x; Φy = ∂Φ/∂y; Φxx = ∂2Φ/∂x2; Φyy = ∂2Φ/∂y2; Φxy = ∂2Φ/∂x∂y. Eq. (8) is solved by finite difference method. The transition from subsonic region to supersonic region is done by the transfer from centered finite difference scheme to the backward scheme (Fig. 2). (a) (b) Fig. 2. a) Finite difference scheme for supersonic regions; b) Generation of grid by solving Laplace equation (profile f10) Calculation of transonic flows around profiles with blunt and angled leading edges 5 Derivatives with centered scheme in subsonic region are written as( ∂Φ ∂x ) k,j = Φk,j+1 −Φk,j−1 xj+1 − xj−1 ,( ∂2Φ ∂x2 ) k,j = (Φx)k,j+1/2 − (Φx)k,j−1/2 xj+1/2 − xj−1/2 . (9) Derivatives with backward scheme in supersonic region are written as( ∂Φ ∂x ) k,j = Φk,j −Φk,j−2 xj+1 − xj−1 ,( ∂2Φ ∂x2 ) k,j = (Φx)k,j − (Φx)k,j−1 xj+1/2 − xj−1/2 . (10) Boundary conditions for the Euler equations (1) and the full potential equation (6) are slid conditions on the profile contour, Joukowski condition at the trailing edge and conditions of small disturbance potential for the far-field boundary. 3. ANALYSIS OF RESULTS 3.1. Comparison of results Numerical results calculated from the programming codes by solving the Euler equations and the full potential equations were compared with each other and compared with experimental results to verify the accuracy of the built programs. Results calculated by using Fluent software are also compared with experimental results and the numerical results to ensure the correctness of Fluent software application operations. The following graphs represent aerodynamic characteristics of profile about the distribution of pressure coefficients cp of upper and lower sides on the dimensionless chord x/C (C is the profile chord, cp = p− p∞ 1 2 ρV∞ with p∞ and V∞ are the pressure and the velocity of the free flow in experimental conditions of zero altitude). Fig. 3 shows numerical results on pressure coefficient calculated from the program- ming codes by the Euler method and the FPE method for Naca 0012, the angle of inci- dence α = 0◦ and the free Mach number M∞ = 0.8 in comparison with experimental results [4]. This case is symmetric, so the results on the upper an on the lower are over- lapping. The comparison shows that similarities of the experimental result and the nu- merical results calculated by the Euler code and the FPE code. There are some small differences at the discontinue position where occurs the shock wave. At this position, numerical results depend on the method of resolution, and it is difficult to determine them exactly by experiment. Fig. 4 presents results on coefficient of pressure and for the case of Naca 0018, α = 2.65◦ and M∞ = 0.7 calculated by Euler code and by Fluent software in comparison with experimental results [4]. In this case, the shock wave on upper side is relatively strong. Values of pressure coefficient represented in Fig. 5 are calculated by Fluent soft- ware in comparison with experimental results for Naca 0012, α = 2.98◦, M∞ = 0.8. The 6 Hoang Thi Bich Ngoc, Nguyen Manh Hung Fig. 3. Pressure coefficient - Comparison of present numerical results (Euler method and FPE method) with experimental results [4] (Naca 0012, α = 0◦, M∞ = 0.8) Fig. 4. Pressure coefficient - Comparison of nu- merical results by Euler method, by Fluent software with experimental results [4] (Naca 0018, α = 2.65◦, M∞ = 0.7) Fig. 5. Pressure coefficient - Comparison of results of Fluent with experiment [4] (Naca 0012, α = 2.98◦, M∞ = 0.8) two cases in Fig. 4 and Fig. 5 are flows having strong shock wave with maximum local Mach number greater than 1.4. These comparisons show similarities between numerical results and experimental results that verify the accuracy of the built codes, as well as the exploitation of Fluent software, and allow applications to the study. 3.2. Cases of subsonic-transonic flow 0.65 < M∞ < 1 Consider the profile Naca 0012 with blunt leading edge (blunt LE) and the para- bolic profile f8 with angled leading edge (angled LE) for an example of calculation. The two profiles Naca 0012 and f8 have maximum thickness around 12% of chord. Fig. 6 presents results of Mach contour and Mach field calculated for Naca 0012 and f8 with incidence angle α = 0◦ and free Mach number M∞ = 0.8. At this flow regime, on the pro- file appear big supersonic regions ended by shock waves with maximum Mach number Calculation of transonic flows around profiles with blunt and angled leading edges 7 Mmax = 1.24 for Naca 0012 and Mmax = 1.3 for f8. However, the position of supersonic regions is very different for two profiles Naca 0012 and f8. The distributions of Mach number and pressure coefficient on the profile contour are shown in Fig. 7. Fig. 6. Mach Contour and Mach field for Naca 0012 and f8 (α = 0◦, M∞ = 0.8) Fig. 7. Mach number and pressure coefficient on profile contour for Naca 0012 and f8 (α = 0◦, M∞ = 0.8) With the angle of incidence α = 0◦, the profiles Naca 0012 and f8 do not create the lift. Consider the drag for these symmetrical cases. In Fig. 8 are presented drag coefficients with free Mach number calculated for Naca 0012 and f8 by using inviscid method (Euler code) and the viscous method (Fluent). It is observed that for the range of transonic flow with M∞ < 1, drag coefficients of the two profiles Naca 0012 and f8 are not very different, and there are not great differences between inviscid and viscous calculations. 8 Hoang Thi Bich Ngoc, Nguyen Manh Hung Fig. 8. Drag coefficient for Naca 0012 and f8 (α = 0◦) Consider the case of the incidence angle α = 4◦. The useful aerodynamic force is the lift. Fig. 9 shows results on the lift coefficient and drag coefficient calculated for Naca 0012 and f8. It might be seen that in the range of transonic flow with M∞ < 1, lift coefficients and drag coefficient of the profile Naca 0012 are respectively greater than those of the parabolic profile f8. Fig. 9. Lift coefficient and drag coefficient for Naca 0012 and f8 (α = 4◦) A special point for the range of free Mach numbers 0.7 < M∞ < 0.9, results calcu- lated by viscous theory and inviscid theory are very different for the two profiles Naca 0012 and f8. In order to explain the effect, see results on Mach lines calculated for flow M∞ = 0.86 in Fig. 10. With the free Mach number M∞ = 0.86, differences between in- viscid and viscous results are relatively largest. In inviscid calculations, the supersonic region on the upper profile is much larger than that on the lower profile. This large dif- ference much shortens in viscous calculations. It is shown that differences in the position Calculation of transonic flows around profiles with blunt and angled leading edges 9 and dimension of supersonic regions appeared on Naca and parabolic profiles are caused by the attack of flow on blunt LE and angled LE. Fig. 10. Mach line - Supersonic zone shock wave by viscous and invicid calculations (Naca 0012 and f8, α = 4◦, M∞ = 0.86) Consider the ratio of lift coefficient on drag coefficient to assess the quality of aero- dynamic profiles that is shown in Fig. 11. With free Mach number M∞ < 0.86, the aero- dynamic quality of Naca 0012 (blunt LE) is much better than that of f8 (angled LE). How- ever, with 1 > M∞ > 0.86, drag coefficients increase and aerodynamic qualities decrease. This is the range of transonic flows having too low efficiency, it is often not used. Fig. 11. Ratio of lift coefficient on drag coefficient for Naca 0012 and f8 α = 4◦ 10 Hoang Thi Bich Ngoc, Nguyen Manh Hung 3.3. Cases of supersonic-transonic flow 1 < M∞ < 1.5 With free Mach number M∞ > 1, shock waves appear right at the surrounding leading edge. Supersonic flows can produce oblique shock waves attached the angled leading edge in conditions of the free Mach number being high enough and the pointed angle being small enough. For free Mach numbers 1 < M∞ < 1.5, it is possible to do not produce an oblique shock wave at angled leading edge, and these are supersonic transonic flows. For these flows, shock waves are detached from the angled leading edge and form a curved shock wave located in front of leading edge. The detached curved shock wave for a profile with angled LE (as f8) causes a big energy loss. That is similar to a profile with blunt LE (as Naca 0012). This phenomenon is illustrated in Fig. 12 with the graphs of drag coefficient calculated for Naca 0012 and f8 at incidence angles α = 0◦ and α = 4◦. It can be seen that only when free Mach numbers M∞ > 1.5, losses of flow through the profile f8 (angled LE) are smaller than losses for the profile Naca 0012 (blunt LE). Note that in these supersonic transonic flows, drag coefficients calculated by inviscid and viscous theories are not much different. This is explained by the absence of interactions between boundary layer and shock wave, and losses of boundary layer are too small by comparison with losses of shock waves. Fig. 12. Drag coefficient for Naca 0012 and f8 (α = 0◦, α = 4◦, M∞ > 1) The formation of detached curved shock wave for a profile with angled LE with free Mach number M∞ = 1.3 (i.e. 1 < M∞ < 1.5) can be seen in Fig. 13, in which on the left side are shown streamlines and shock waves for the parabolic profile f8 (α = 0◦), and on the right side presented the limiting line for oblique shock waves of wedge [5]. The parabolic profile f8 has angled leading edge of θ = 13.86◦ that corresponds to the critical Mach number M∞min = 1.55. So, with M∞ = 1.3 < M∞min a detached curved shock wave appears in front of leading edge. With the incidence angle α = 4◦, the flow becomes asymmetric between the upper and the lower sides. Due to this asymmetry, at the trailing edge, oblique shock wave on the upper side much stronger than oblique shock wave on lower side. Results in Fig. 14 illustrate this comment. At the free Mach number M∞ = 1.3, subsonic regions formed Calculation of transonic flows around profiles with blunt and angled leading edges 11 Fig. 13. Streamlines, detached curved shock wave for f8 and limiting line for oblique shock waves of wedge Fig. 14. Curved shok waves and subsonic regions for f8 and Naca 0012 (M∞ = 1.3, α = 4◦) around the trailing edge are very different for the blunt LE profile Naca 0012 and the angled LE profile f8. Fig. 15 presents comparisons between results of aerodynamic quality for the profile f8 and the profile Naca 0012 with α = 4◦, 1 ≤ M∞ ≤ 2. It is observed that with M∞ < 1.3, the aerodynamic quality of profiles f8 and Naca 0012 is equivalent. Only with free Mach numbers M∞ > 1.6 (supersonic flow), the aerodynamic quality of profile with angle leading edge is increased by the shift from detached shock waves into attached shock wave. Obviously, the aerodynamic quality of profile Naca 0012 with blunt leading edge is too bad in supersonic flows M∞ > 1.6 due to a strong detached shock wave in front of leading edge. 12 Hoang Thi Bich Ngoc, Nguyen Manh Hung Fig. 15. Ratio CL/CD for f8 and Naca 0012 (α = 4◦, M∞ > 1) 4. CONCLUSIONS In fact, transonic flows are used for large civil aircraft and other flight objects. The implementation of computational methods and deep understandings of physical nature of the transonic flow are always needed [6, 7]. The above analysis for subsonic-transonic flow and supersonic-transonic flow allows drawing the following conclusions. It is clear that transonic flows are very complex problems. Mathematically, the differential equations of transonic motion take the mixture between elliptic, parabolic and hyperbolic types. Physically, transonic flows are very sensitive to any changes in the aerodynamics, meaning that only a very small change in geometry or dynamics, it can lead to very large changes in the physical nature of the phenomenon, and of course a very large change in aerodynamic characteristics. For transonic flow problems, it is difficult to draw any analytic rule. The problems require solving the differential equations for each specific case. This report refers to the transonic flow problems, the use of profiles having blunt and angled leading edges with advantages and disadvantages and practical applications. In order to draw conclusions on this issue, it is necessary to split the Mach range 0.65 < M∞ < 1.5 into two domains 0.65 < M∞ < 1 (subsonic-transonic flow) and 1 < M∞ < 1.5 (supersonic-transonic flow). - With the supersonic-transonic flow domain (M∞ > 1), the profile with blunt leading edge is not usable due to the formation of detached shock wave in front of the leading edge that causes too large losses of energy. However, if using the profile with angled leading edge, it is only effective when forming oblique shock waves attached the leading edge. But with free Mach numbers 1 < M∞ < 1.5, shock waves are detached the nose, even using an angled nose. So, the range of free Mach numbers 1 < M∞ < 1.5 should not be used in practice. - With the subsonic-transonic flow domain (0.65 < M∞ < 1), using the profile with blunt leading edge is more effective. However, supersonic flow regions formed on the profile are very different in terms of strength and size, very different on the upper side and the lower side, and they can produce normal shock waves ended on the profile or oblique shock waves at the trailing edge. The aerodynamic characteristics of the gas - solid interaction change very sensitively. The range of free Mach numbers 0.9 < M∞ < 1 Calculation of transonic flows around profiles with blunt and angled leading edges 13 has very low efficiency, should not be used. The range of free Mach numbers 0.7 < M∞ < 0.9 has very high performance if choosing optimal modes. - Especially, for transonic flows with free Mach numbers 0.7 < M∞ < 0.9, nor- mal shock waves often appear on the profile. This is also the mode that has a sensitive interaction between the shock wave and the boundary layer. The interaction can make a separation of boundary layer and causes large variations of aerodynamic character- istics [8]. So, inviscid and viscous calculations give different results. For this range of free Mach number, it is necessary to use the viscous theory and especially experimental results. These transonic flows are extremely sensitive to very small changes in profile ge- ometry, incidence angle and free Mach number. Only a very small change of one among these parameters can lead to very large changes in the aerodynamic characteristics and they can be transformed from positive to negative and vice versa. Therefore, it requires a high accuracy in the calculation. REFERENCES [1] C. Hirsch. Numerical computation of internal and external flows: Fundamentals of computational fluid dynamics, Vol. 1. Elsevier Pub. B. H., Great Britain, (2007). [2] C. Hirsch. Numerical computation of internal and external flows: Computaional methods for invicid and viscous flows, Vol. 2. John Wiley & Sons Pub., England, (1994). [3] J. D. Cole and L. P. Cook. Transonic aerodynamics. Elsevier, Amsterdam, (2012). [4] F. W. Riegels. Aerofoil sections. Butterworths Pub., London, (1961). [5] H. T. B. Ngoc and N. M. Hung. Calculating shock wave angles and drags of supersonic flows through cones and wedges. In Proceedings of the 14th Asia Congress of Fluid Mechanics, (2013), pp. 196–201. [6] W. Hassan and M. Picasso. An anisotropic adaptive finite element algorithm for transonic viscous flows around a wing. Computers & Fluids, 111, (2015), pp. 33–45. [7] S. V. S. A. Hema Sai Chand, K. Giridhar, T. Keerthi Goud, B. Vamshi Bhargav, and P. Srini- vas Rao. Transonic shockwave/boundary layer interactions on naca 5 series -24112. Interna- tional Journal of Current Engineering and Technology, (Special Issue-2, 2014), pp. 629–634. [8] H. T. B. Ngoc and N. M. Hung. Study of separation phenomenon in transonic flows produced by interaction between shock wave and boundary layer. Vietnam Journal of Mechanics, 33, (3), (2011), pp. 170–181.