Improvement of Airfoils Aerodynamic Efficiency by Thermal Camber Phenomenon at Low Reynolds Number

1.Ferdowsi University of Mashhad – Faculty of Engineering – Mechanical Engineering Department – Mashhad – Iran. Correspondence author: Mohammad Hassan Djavareshkian | Ferdowsi University of Mashhad – Faculty of Engineering – Mechanical Engineering Department | Azadi Square | 9177948974 – Mashhad – Iran | E-mail: javareshkian@um.ac.ir Received: May 27, 2017 | Accepted: Nov. 22, 2017 Section Editor: T John Tharakan ABSTRACT: In this research, viscous, laminar and steady flow around symmetric and non-symmetric airfoils is simulated at Low Reynolds Number (LRN). Navier-Stokes (N-S) equations are discretized by Finite Volume Method (FVM) and are solved by the SIMPLE algorithm in an open source software, namely OpenFOAM. The main objective of this paper is the introduction of the thermal camber phenomenon. This phenomenon is used to improve the aerodynamic performance. Hence, a symmetric airfoil, like NACA0012, with thermal camber is compared with the airfoils with the physical camber, including NACA2412 and NACA4412, to specify which camber type has more effects on the aerodynamic efficiency. Furthermore, various temperatures are tested in order to find the optimum condition. After validation, results indicated that cooling upper surface and heating lower surface, namely thermal camber, generate lift force and improve aerodynamic performance for symmetric airfoils at a 0° Angle of Attack (AOA). These improvements are more than the airfoils with physical camber. Also, when this method is applied to the NACA2412 and NACA4412 airfoils, lift to drag coefficient ratio will increase more than the condition with only cooling or heating the surfaces.


INTRODUCTION
Due to the importance of Micro Aerial Vehicles (MAV), many studies have been done to develop this field of research.The motivation of recent attempts is to improve the aerodynamic performance in Low Reynolds Number (LRN) flows.In this range of Re, the adverse pressure gradient plays a destructive role in aerodynamic performance and consequently causes flow separation over the airfoil.The separation happens either at leading edge or trailing edge.Finally, this separation increases the drag and reduces the lift (Febi Ponwin and Rajkumar 2015).Therefore, many methods, including optimization of airfoil shape, blowing and suction, have been used in those studies to overcome mentioned problems.On the other hand, the theory of using heat transfer based on the viscous boundary layer characteristics can also be utilized to improve the aerodynamic efficiency.Thermal effects on the flow can be clarified by investigation of the airfoils in steady conditions.Norton et al. (1973)  experimental studies to investigate the boundary layer of a NACA0012 airfoil at the Re of 10 5 .Results of wind tunnel tests showed that the drag increases by heating the surfaces.Landrum and Macha (1987) experimentally studied the effects of leading edge heating on two dimensions boundary layer characteristics.They found no alteration in lift force.Mabey (1990) studied the effects of the heat transfer on aerodynamic of different bodies at the Re of 10 6 .He revealed that cooling the surfaces decreases Re and delays the flow separation.By solving the N-S equations at the Re of 10 6 , Longo and Radespiel (1995) found that the maximum lift, the maximum lift to drag coefficient ratio and the lift curve slope of a NACA0012-64 airfoil are reduced by heating the surfaces.However, they showed that at low angles of attack, the heat transfer enhances lift to drag coefficient ratio.Navier-Stokes equations were used by Raghunathan and Mitchell (1995) to simulate the subsonic and transonic flows over a flat plate and a NACA0012 airfoil, respectively, in the presence of heat transfer.Results indicated that in a steady condition, cooling the surfaces thins boundary layer and increases its momentum.By performing numerical study on a DOAL3 airfoil at the Re of 10 6 , Stock (2002) found that the drag approximately doubles in the presence of the heating compared to the adiabatic condition.To summarize, it can be deduced that cooling the surfaces reduces the drag while it has no influence on the lift force of the symmetric airfoils.Kim et al. (2003) revealed that for small-scale airfoils, cooling the upper surface and heating the lower surface enhance the lift and reduce the drag.However, they only focused on symmetric airfoils, and the simulation was accomplished at the Re of 10 4 .Bekka et al. (2009) concluded that thermal effects are more sensible in LRN range.Also, by simulation of flow around a NACA0012 airfoil, they showed that the results are approached to experimental data by using turbulence models.Using CFD methods and solving N-S equations at the Re of 3000, Hinz et al. (2013) investigated the impacts of the heat transfer on aerodynamic performance of a flapping NACA0012 airfoil.The temperature on both surfaces of the airfoil in their simulations was kept hot or cold.
As is well known, one way to increase the lift force is physical camber.However, this method also enhances the drag force.The purpose of this research is to enhance the aerodynamic performance of the airfoil by thermal chamber phenomenon.The effects of thermal camber phenomenon on cambered airfoils have not been investigated yet.Therefore, the main objective of this research is to apply thermal camber on symmetric airfoils as well as cambered airfoils in order to weigh the influence of the thermal camber on the aerodynamic performance against the physical camber.Furthermore, various temperatures are tested to find optimum condition.In this simulation, although the flow velocity is in incompressible flow regime, it is considered as the compressible flow regime due to the temperature variations.Thus, the Compressible N-S equations are discretized by FVM and are solved by the compressible SIMPLE algorithm in an open source software, namely OpenFOAM (Open CFD Limited 2008).

NUMERICAL SIMULATION
Fluid properties are generally not constant with temperature alteration.Property variations with temperature are shown in Fig. 1.Therefore, the following relations are used in these studies (Eqs.1-3) (Eucken 1940;Incropera and Dewitt 1990;Sutherland 1893): (1) (2) (3) equations are discretized by FVM and are solved by the compressible SIMPLE alg open source software, namely OpenFOAM (Open CFD Limited 2008).

NUMERICAL SIMULATION
Fluid properties are generally not constant with temperature alteratio variations with temperature are shown in Fig. 1.Therefore, the following relations these studies (Eqs.1-3) (Eucken 1940;Incropera and Dewitt 1990;Sutherland 1893) where: kT, µ, cv and R denote thermal conductivity, dynamic viscosity, specific hea volume and particular gas constant, respectively.Equation 2 is evaluated by Suthe The governing equations, namely the conservation of mass, momentum and enthalpy are given in the Cartesian tensor form by (Eqs.4-6): where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent var Newtonian and fully laminar viscous tensor of stress, τij, is given by (Eq where: δ is Kronecker symbol.Governing equations in the form of ve coordinates are shown as follows (Eq.8):
where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor and u is subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent variables.F Newtonian and fully laminar viscous tensor of stress, τij, is given by (Eq.7) (Hinz where: δ is Kronecker symbol.Governing equations in the form of vector and coordinates are shown as follows (Eq.8):
where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor and u is v subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent variables.Flow Newtonian and fully laminar viscous tensor of stress, τij, is given by (Eq.7) (Hinz e where: δ is Kronecker symbol.Governing equations in the form of vector and in coordinates are shown as follows (Eq.8):
where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor and u is v subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent variables.Flow Newtonian and fully laminar viscous tensor of stress, τij, is given by (Eq.7) (Hinz e where: δ is Kronecker symbol.Governing equations in the form of vector and in coordinates are shown as follows (Eq.8):
where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor and u is subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent variables.F Newtonian and fully laminar viscous tensor of stress, τij, is given by (Eq.7) (Hin where: δ is Kronecker symbol.Governing equations in the form of vector and coordinates are shown as follows (Eq.8):
where: ρ is density; p is pressure; h is enthalpy; τ is stress tensor and u is velocity.Also, subscripts of i and j denote coordinate direction.
The stress tensor is usually expressed by basic dependent variables.Flow is assumed Newtonian and fully laminar viscous tensor of stress, τ ij , is given by (Eq.7) (Hinz et al. 2013): where: δ is Kronecker symbol.Governing equations in the form of vector and independent of coordinates are shown as follows (Eq.8): where: φ is variable and can be u i , u j and h.Equation 8 is discretized by FVM and is solved by pressure based algorithm.
A Preconditioned Conjugate Gradient (PCG) solver is used in this paper to solve pressure equation systems.Also, rhoSIMPLEFOAM solver is considered in order to simulate the present study.

RESULTS
As shown in Fig. 2, a C-type mesh is used to model the flow around the airfoil.The computational area that is depicted in Fig. 3 is far enough from the airfoil where normal and back distances from the airfoil are set to 20 c and 25 c, respectively.The boundary conditions of non-slip velocity, zero-gradient pressure and isothermal temperature are considered on the airfoil surfaces.At the inlet, the free stream velocity and the zero-gradient pressure are chosen as the physical and numerical boundary conditions, respectively.In other words, pressure is determined by extrapolating inside of the computational area.Conversely, the boundary conditions of the pressure outlet and the zero-gradient velocity are used for the outlet.Velocity is computed by extrapolating inside of the computational area.Also, the slip boundary condition is used for the top and bottom boundaries.Velocity, pressure and temperature in these boundaries are the same as the first cells of the boundary.Mesh independency is done in Reynolds number, Re, of 3000, free stream temperature of 300 K and airfoil surface temperatures, T S , of 300 K, 350 K and 400 K. Various mesh numbers of 2000, 3000, 6000, 18000 and 25000 are considered to study the mesh independency.Drag coefficients, C d , of these mesh numbers are shown in Fig. 4. By raising mesh numbers from 18000 to 25000, although computational cost increases, the accuracy is almost constant.Therefore, the mesh number of 18000 is employed.
In order to validate current results, the drag coefficient of the present research is compared with that of Hinz et al. (2013).Geometrical characteristics and flow conditions of their research are listed in Table 1. Figure 5 indicates high accuracy of the present simulation in comparison with the mentioned numerical data.This is because of using proper mesh with almost orthogonal cells.

xx/xx 05/13
Sutherland and Eucken's relations are used to compute the dynamic viscosity and the thermal conductivity, respectively.According to Fig. 1, the dynamic viscosity increases with temperature increment.For Newtonian fluids, the viscous shear stress is proportional to dynamic viscosity and velocity gradient normal to the wall.As a result, the drag coefficient increases when the dynamic viscosity enhances.Drag enhancement increment is not desired from the aerodynamic point of view.Variation of the drag coefficient with surface temperatures between 200 K and 400 K are illustrated in Fig. 6.As shown in Fig. 6, temperature increment from 300 K to 400 K causes dynamic viscosity to be increased by 24.6%.However, dynamic viscosity is reduced by 28.2% as a result of reduction in temperature from 300 K to 200 K.Therefore, variation of drag in cooling is more than of heating and the diagram has become slightly non-linear.According to Kim et al. (2003) and this research results, cooling the upper surface, T u , and heating the lower surface, T l , can significantly enhance the aerodynamic efficiency.In other words, lift can be produced in    Hinz et al. (2013) for NACA0012 airfoil at Re = 3000.
xx/xx 06/13 aforementioned condition even at a 0° AOA for a symmetric airfoil like NACA0012.This influence is introduced as thermal camber phenomenon.Hence, the upper and lower surface temperatures of a NACA0012 airfoil are set to 200 K and 400 K, respectively.
Cooling the upper surface causes acceleration of the flow and subsequently decreases the pressure.This procedure is reversed for the lower surface.Consequently, pressure difference is created in order to generate lift force.Velocity and pressure contours for a NACA0012 airfoil at the surface temperature being equal to the free stream temperature are reported in Figs.     and lower surface temperature of 400 K. Enhancement and decrement of velocity, respectively, in upper and lower surfaces can be obviously seen in these figures.In the other words, symmetric pressure distribution will be converted to a non-symmetric pressure distribution in order to create positive lift.To investigate this point in more detail, pressure coefficient diagram is plotted in Fig. 11.According to this figure, generated lift in this condition is equal to 0.88.Therefore, it can be deduced that a NACA0012 airfoil, which is a symmetric type, behaves like a non-symmetric airfoil and produces lift at a 0° AOA.
xx/xx 08/13 Herein, two airfoils with physical camber, including NACA2412 and NACA4412 airfoils, are analyzed at a 0° AOA for the surface temperature of the airfoils being equal to the free stream temperature.To validate the simulations of non-symmetric airfoils, the results of present simulations are compared with those of Abbott and Von Doenhoff (1959).Simulation is conducted in high Re.It is worth mentioning that using of turbulence model is unavoidable in this condition.Therefore, two turbulence models, Spalart-Allmaras (S-A) and SST k -w, are considered to do simulation.Validation results are listed in Table 2.It is clear that the obtained results have a good agreement with those of Abbott and Von Doenhoff (1959) especially by using SST k -ω and S -A models for NACA2412 and NACA4412 airfoils, respectively.Static pressure contour of NACA2412 and NACA4412 airfoils are shown, respectively, in Figs. 12 and 13 at Re = 3000, c = 4 cm and for the surface temperature of the airfoils being equal to the free stream temperature.According to these figures, pressure has the same behavior with that of NACA0012 airfoil in Fig. 10.The pressure coefficients for three airfoils are presented in Fig. 14.The upper and lower surface temperatures of NACA0012 airfoil are set to 200 K and 400 K, respectively.However, the two other airfoils are at a surface temperature similar to the free stream temperature.As shown in Fig. 14, NACA0012 airfoil with the help of thermal camber creates pressure difference in order    to produce positive lift.Drag and lift coefficients of these airfoils are shown in Table 3 at the conditions of Fig. 14.As shown in this Table, the lift coefficient, C l , is reduced to almost 0.2 of the validation result because the Re is low in this research and this point causes large growth of laminar boundary layer in comparison with the condition with high Re.In high Re, the boundary layer will be attached strongly to the airfoil due to the great momentum of the flow but in LRN flow, high growth of laminar boundary layer on the airfoil can reduce the influence of the thermal camber.It can be deduced from Table 3 that in this Re, the NACA0012 airfoil can produce the lift coefficient almost equal to the NACA4412 airfoil when the upper and lower surfaces are cooled and heated to 200 K and 400 K, respectively.Furthermore, it has less drag coefficient than the NACA4412 airfoil.Consequently, the lift to drag coefficient ratio of NACA0012 airfoil is more than that of NACA4412 airfoil.In addition, it can be concluded that thermal camber has more influence than physical camber on the aerodynamic coefficients.Additionally, when the upper and lower surfaces of NACA0012 airfoil are cooled and heated, respectively, to 250 K and 350 K, the lift coefficient of the NACA0012 airfoil can be reached to almost the NACA2412 airfoil.Also in this condition, the lift to drag coefficient ratio of the NACA0012 airfoil is more than that of the NACA2412 airfoil.
Herein, the NACA2412 and NACA4412 airfoils are analyzed in various conditions of cooling, heating and cooling-heating.Lift and drag coefficient variations of the NACA2412 and NACA4412 airfoils are shown in Figs. 15 and 16, respectively.Also, the lift to drag coefficient ratio of NACA2412 and NACA4412 airfoils are represented in Fig. 17.According to these figures, the drag coefficient has same trends with the NACA0012 airfoil.But cooling or heating the NACA0012 airfoil surfaces only affects the drag coefficient and, due to the symmetric geometry of the airfoil, the lift coefficient does not change considerably.However, cooling both surfaces of the non-symmetric airfoils leads to lift production because cooling has more influence on the upper surface of the non-symmetric airfoil than the other surface.These effects enhance the flow velocity on the upper surface and cause more pressure reduction on this surface.Consequently, more lift is produced than the condition with the surface temperature of the airfoil being  equal to the free stream temperature.The reversed condition happens for heating.It seems that in order to reach the maximum aerodynamic efficiency, cooling the upper surface and heating the lower surface have better results.These outcomes for two airfoils are listed in Table 4.A similar trend as the NACA0012 airfoil in    4. In other words, by choosing specific temperature to cool and heat, respectively, the upper and lower surfaces of the NACA2412 airfoil, its efficiency can be raised to the value of NACA4412 airfoil.For instance, the same lift to drag coefficient ratio with the NACA4412 airfoil can be achieved by the NACA2412 airfoil with the upper surface temperature of 250 K and the lower surface temperature of 350 K.As explained before, this is because of the acceleration of the velocity at the upper surface and reduction of the lower surface velocity.This can cause pressure difference, which causes the lift force to be increased.This pressure difference is more than the condition when both surface of the airfoil is cooled.Additionally, according to Table 4 and Fig. 17, decreasing temperature of both NACA2412 airfoil surfaces from free stream temperature to 200 K causes a 36% increase in the lift to drag coefficient ratio.However, when the upper and lower surface temperatures are considered 250 K and 350 K, this improvement reaches 89%.Also, this value rises to 196% when the upper and lower surface temperatures are 200 K and 400 K, respectively.This improvement for the NACA4412 airfoil is less than the NACA2412 airfoil.For instance, when both surfaces of the NACA4412 airfoil are cooled from 300 K to 200 K, the lift to drag coefficient ratio increases by 33%, that is 3% less than of the NACA2412 airfoil.When the upper and lower surface temperatures of the NACA4412 airfoil are set, respectively, to 200 K and 400 K,improvement of the lift to drag coefficient ratio is 97%, which is 99% less than the NACA2412 airfoil in the same condition.Generally, it can be concluded that cooling the upper surface and heating the lower surface have more advantages than equal cooling of both surfaces.Pressure coefficient diagram for the NACA2412 airfoil is demonstrated in Fig. 18.Increment and decrement of the pressure, respectively, on lower and upper surfaces by thermal effects, can be seen in this figure.Figure 19 represents the static pressure contour of the NACA2412 airfoil with the upper surface temperature of 200 K and lower surface temperature of 400 K.In comparison with Fig. 12, increment and decrement of the pressure, respectively, on lower and upper surfaces can be justified.Also, in comparison with Fig. 13, the aerodynamic improvement of the NACA2412 airfoil against the NACA4412 airfoil without thermal effects is clear.

Figure 4 .
Figure 4. Mesh independency for NACA0012 airfoil at Re = 3000 and various cell numbers.

Figure 5 .
Figure 5.Comparison between drag coefficient of present simulation and that ofHinz et al. (2013) for NACA0012 airfoil at Re = 3000.

Figure 7 .
Figure 7. Velocity contour around NACA0012 airfoil at Re = 3000 and surface temperature of 300 K.

Figure 6 .
Figure 6.Variation of drag coefficient with surface temperature for NACA0012 airfoil at Re = 3000.

Figure 10 .
Figure 10.Pressure contour around NACA0012 airfoil at Re = 3000 and T u = 200 K and T l = 400 K.

Figure 12 .
Figure 12.Pressure contour of NACA2412 airfoil in free stream temperature.

Figure 16 .
Figure 16.Variation of lift and drag coefficient with surface temperature for NACA4412 airfoil.

Figure 17 .
Figure 17.Variation of lift to drag coefficient ratio with surface temperature for NACA2412 and NACA4412 airfoils.
200 K T l = 400 K performed analytical and

Table 2 .
Validation of NACA2412 and NACA4412 airfoils with experimental data of Abbott and Von Doenhoff (1959).
Table 3 is seen for the NACA2412 and NACA4412 airfoils in d Figure 15.Variation of lift and drag coefficients with surface.

Table 4 .
Lift and drag coefficients for NACA2412 and NACA4412 airfoils at Re = 3000.