Abstract
A new aerodynamic force of drag prediction approach was developed to compute the airfoil drag via entropy generation rate in the flow field. According to the momentum balance, entropy generation and its relationship to drag were derived for viscous flow. Model equations for the calculation of the local entropy generation in turbulent flows were presented by extending the RANS procedure to the entropy balance equation. The accuracy of algorithm and programs was assessed by simulating the pressure coefficient distribution and dragging coefficient of different airfoils under different Reynolds number at different attack angle. Numerical data shows that the total entropy generation rate in the flow field and the drag coefficient of the airfoil can be related by linear equation, which indicates that the total drag could be resolved into entropy generation based on its physical mechanism of energy loss.
1. Introduction
Accurate prediction of the aerodynamic force is a critical requirement in aircraft design. CFD methods have been widely applied in aerodynamic design and optimization of aircraft. However, even for airfoils with attached flow at relatively low angles of attack, the predicted drag based on integration of the surface pressure and skin-friction distributions can be off by more than 100% even though the computed surface pressure and skin friction are in good agreement with the experimental data [1]. Therefore, designers might need a calculation method which could resolve the drag according to its mechanism of production.
There are two common approaches to predict total drag force of airfoils or wings, a standard surface integration method, and a wake integration method. The surface integration method relies on calculations of pressures and skin friction over a series of flat surfaces. However, the surface integration was met with difficulties especially for the complex configuration due to the need to approximate the curved surfaces of the body with flat faces which can be affected by significant errors introduced by the “numerical viscosity” and “discretization” error of the numerical solution [2]. This has led various researchers to look at the experimental wake integral method, which is derived from the momentum equation of the governing equations of fluid mechanics and to attempt to apply them to CFD computations [3, 4]. The main advantage of this technique is that no detailed information on the surface geometry of the configuration is required and also has the drag decomposition capability into wave, profile, and induced drag component. Zhu et al. [5] applied the wake integration method to predict drags of some examples including airfoil, a variety of wings, and wing-body combination. They introduced the cutoff parameters method to reduce the computational time required for integrating the wake cross plane. In the method, cells in the wake cross plane which contain small proportions of vorticity and entropy are not included in the integration. Numerical results were compared with those of traditional surface integration method, showing that the predicting drag values with the wake integration method are closer to the experimental data.
Oswatitsch [6] derived a far-field formula of the entropy drag considering first-order effects, in which the drag is expressed as the flux of a function only dependent on entropy variations. In [7, 8], Oswatitsch’s formula is used for computing the entropy drag in RANS solutions by limiting the far-field flux computation to a box enclosing the aircraft. By applying the divergence theorem of Gauss, the surface integral in Oswatitsch’s formula can be replaced with a volume integral. Therefore, the integrand can be set to zero a priori in regions where it is known that physical entropy variations should be zero, thus, removing spurious contributions to drag. In a viscous domain, all drag on an airfoil is eventually realized as entropy generation, so we can establish a relationship between total entropy generation of the flow field and drag on the airfoil. Li and Stewart reported on 2D validation studies for predicting drag by means of calculating the generation both for airfoils and for viscous pipe flows [9]. Stewart [10] and Monsch et al. [11] calculated the drag force of a 3D wing by performing a volume integration within the numerical domain and a surface integral over the outlet of the domain.
It is pointed out that a representation of losses in terms of entropy generation offers significant insight into the flow and thermal transport phenomena over the airfoil and provides an effective tool for drag prediction. Entropy analysis is a method to evaluate a process based on the second law of thermodynamics. It is basically calculating entropy generation in a system and its surroundings and using it as a proxy for the evaluation of the energy loss [12]. Mahmud and Fraser [13] analyzed the mechanism of entropy generation and its distribution through fluid flows in basic channel configurations including two fixed plates and one fixed and one moving plates by considering simplified or approximate analytical expressions for temperature and velocity distributions and derived analytically general expressions for the number of entropy generation and Bejan number. A detailed review of entropy and its significance in CFD is presented by Kock and Herwig [14].
In this paper entropy generation and its relationship to drag were derived for viscous flow. we present model equations for the calculation of the local entropy generation rates in turbulent flows by extending the RANS procedure to the entropy balance equation. This equation serves to identify the entropy generation sources, without need to solve the equation itself. The main objective of this paper is to demonstrate the viability of entropy-based drag calculation method and to compare the consistency of predicting the drag of single-element airfoils using surface integration, wake integration, and entropy generation integration.
2. Alternative Methods for Airfoil Drag Prediction
The conservation law of momentum to the control volume enclosing the configuration makes it possible to predict the total aerodynamic drag force by an integration of stresses on the aircraft configuration (surface integration) or by an integration of momentum flux on a closed surface far from the configuration (wake integration). We consider a steady flow with freestream velocity , pressure , and temperature around the configuration; the only external force acting on the body is due to the fluid. Neglecting body forces, the fundamental formula for the total aerodynamic force acting on the configuration in the Cartesian reference system can be derived as follows [1]:where represents the entire control surface, denotes the outward unit vector normal to , denotes the velocity vector, represents the fluid density, and represents the shear stress tensor. Note that consists of , and with specifies the aircraft surface and the closed surface far from the configuration. Consequently, (1) can be written asThe momentum flux term in the first integral is zero in the general case of solid body surface. Thus (2) reduces towhere the integral on the left-hand side represents the total aerodynamic drag force by an integration of stresses on the aircraft configuration which can be written as
The far-field expression of drag prediction is given by the right-hand side of (3):
The wake expression for drag is derived from (5) by moving the inlet and side faces of the control volume to infinity and the following wake integration for the drag is obtained:where is perpendicular to the freestream direction. The third integral on the right-hand side of (6) is often neglected if an integration plane sufficiently far down stream [4]. Considering the conservation of the mass flux we have
In a two-dimensional flow, as we move the boundary towards infinity (), the pressure term vanishes as proved by Paparone and Tognaccini [4]:and in the far-field drag expressions the viscous stresses are usually neglected for a high Reynolds number flow; (5) can be simplified to
For ideal gas, the module of the velocity can be expressed in terms of variations of total enthalpy (), entropy (), and static pressure () [4]:where is the ratio between the specific heats of the fluid and is the gas constant for air. Hence, the drag can be expressed as the flux of across :
Equation (10) can be expanded in Taylor’s series ignoring second-order term: the coefficients of the series expansion depend on and the freestream Mach number (); they are
After substitution of expression (12) into (10), the contributions of pressure, entropy, and total enthalpy are isolated. Moreover, all terms with vanish on . Hence, the far-field drag expression around a body in two-dimensional flows becomes
The term depending on is negligible in the case of power-off condition and the far-field drag expression is represented by the entropy:
Finally, (15) can be expressed in volume integral form by applying Gauss’s theorem in the domain :The differential balance equation of the entropy leads to
According to Gouy-Stodola theorem, the relationship between the exergy destruction and the entropy generation is defined bywhere is the absolute temperature of the environment in Kelvin [15]. As a result of this theorem, the amount of irreversibility is directly proportional to the entropy generation and is responsible for the inequality sign in the second law of thermodynamics. Based on (17), the drag and entropy generation rate can be related by a linear balance equation. This is essential because drag can be directly estimated from entropy generation without other effects and the balance is correct only when the numerical volume completely encloses the entropy changes caused by the airfoil. Equation (17) provides a method to predict total drag force of airfoil directly by performing a volume integration of entropy generation rate within the numerical domain.
3. The CFD Procedure for Entropy Generation
3.1. The Direct Method of Calculating Entropy Generation Rate
The entropy is a state variable and the transport equation for entropy per unit volume in Cartesian coordinates can be expressed aswhere stands for the three directions and denotes the velocity component in the direction [14].
There are basically two methods how entropy generation can be determined [16]. In the direct method, the entropy generation rate is connected with the entropy transport equation:Equation (19) consists of two terms. The first term denotes the entropy generation due to viscous dissipation () while the second term denotes the entropy generation due to heat conduction (). The entropy generation terms are calculated in the postprocessing phase of a CFD calculation based on (20). That means, they are determined by using the known field quantities velocity and temperature. Integration of these field quantities over the whole flow domain results in the overall entropy generation rate.
In the indirect method, the entropy generation is calculated by equating it to the rest of (18):Since one is interested in the total entropy generation rate of the flow field, (20) must be integrated over the entire flow domain. This corresponds to the fact that the global balance can be cast into the following form:
Obviously the direct method is superior and should be applied in complex flow situations. And, there is one more important advantage of this method [16]: from the direct method we get the information of how the overall entropy generation is distributed, an information which the indirect method cannot provide. It may, however, help to understand the physics of the complex process and be important in finding ways to reduce the overall entropy generation in a technical device.
When the turbulent flow is considered, the derivation of the entropy generation rate is carried out in terms of the RANS equations which splits the velocities and temperature into time-mean and fluctuating components; that is, . Note that average entropy generation rate per unit volume in turbulent flow can be expressed as follows:the eddy viscosity-type assumption is adopted which is consistent with the RANS turbulence model there [12]:where and and are laminar and turbulent viscosity, respectively. And the effective thermal conductivity is also replaced by , where is the specific heat at constant pressure and is the turbulent Prandtl number.
3.2. Flow Field Simulation
In the flow field calculation, an in-house code was used. The RANS based turbulence models are used in conjunction with the Navier-Stokes equations for viscous flow simulations. The numerical simulations discussed herein use the general steady viscous transport equations in conservative form which can be casted into the following compact notation form:where is the general transport variables, is the nominal diffusion coefficient, and is the source term which can be expressed as a linear function of [17]. Thus can be written as follows:where stands for the constant part of , while is the coefficient of .
We discrete the transport equations by FVM (Finite Volume Method) on a nonorthogonal collocated grid that all transport variables are stored at cell centers and the integration and discretization about the control volume yieldswhere the summation is over the faces of the control volume.
The deferred correction method [18] is used to discrete the convection term while it can be expressed as where is the mass flow rate (defined to be positive if flow is leaving ) through the face.
The diffusion term at the face iswhere is the area vector. The form of this term in nonorthogonal grid schemes can be decomposed into orthogonal and nonorthogonal terms [19]. Thuswhere at the face is taken to be the average of the derivatives at the two adjacent cells. The first term on the right-hand side of (29) represents the primary gradient and the discretization is equivalent to a second-order central different representation which is treated implicitly and leads to a stencil that includes all neighboring cells while the second term is the secondary or cross-diffusion term which is treated explicitly.
For evaluating , the cell derivative of , Gauss’s theorem is adopted and the reconstruction gradient is estimated as where the face value can be linear reconstructed from the cell neighbors of the face:
For steady flows, strongly implicit procedure (SIP) [20] iteration technique is adopted to solve the algebraic equations and to accelerate the rate of convergence. Since pressure and velocity components are stored at cell centers, computing face mass flow rate by averaging the cell velocity is prone to checker boarding. In this paper, momentum interpolation method (MIM) [21] is adopted to overcome this.
4. Numerical Examples and Discussion
4.1. Entropy Generation Calculation Validation
So far, we have shown how the total entropy generation in turbulent flows can be calculated in a postprocess. In this section, entropy generation calculation has been carried out using a turbulent flow over NACA standard series airfoils to validate the eddy viscosity model used for entropy generation calculations in the present paper, details of which are provided below. We discussed the relationship between the total entropy generation in the flow field and the drag coefficient at various angle-of-attack under different Reynolds number with different turbulence models.
The nonorthogonal structured mesh is generated by algebraic method. As shown in Figure 1, a grid is used and 120 grid points are distributed over the airfoil. The top and bottom far-field boundaries are located at 12.5 chord lengths from the airfoil. The upstream velocity inlet boundary is 12.5 chord length away from the airfoil trailing edge while the downstream outflow boundary is located 21 chord lengths away from the airfoil. The upstream boundary is set to be velocity inlet and the downstream boundary is set to be outflow while the velocity profile normal to the exit plane is adjusted to satisfy the principle of global conservation of mass flow rate.
(a)
(b)
4.1.1. Turbulence Model Comparison
Shuja et al. [22] reported a dependency between the various turbulence models used and the entropy generation estimate on an impinging jet flow. This is a consequence of the differing estimates of the effective viscosity used in the different models. In this paper, five turbulence models, -, RNG -, -, SST -, and S-A, which are based on eddy viscosity-type assumptions are used to model the flow around NACA0012 airfoil at under different turbulence models in order to discuss the effect of the turbulence models on entropy generation calculation.
The surface pressure distribution at different angle-of-attack compared with the experimental data [23, 24] is plotted in Figure 2. In general, all models performed quite well in the simulation yielding predictions which are in excellent agreement with measurements, as shown in Figure 2. Thus, we thought the numerical methods and discretization scheme for flow simulation are justified.
(a)
(b)
(c)
Contours of entropy generation rate around the NACA0012 airfoil for angle-of-attack , , and at under different turbulence models can be seen in Figure 3. Because the order of magnitude of the entropy generation rate is , this paper performs a logarithmic process to clearly show the source of the entropy generation. For all turbulence models, it can be concluded from the present results that the bad regions that most of the entropy generated are the front, the near wall, the turbulent wake, and the recirculation zones. The high entropy generation in the near wall region is due to the large velocity gradients that develop to establishing a smooth velocity transition from the wall values while the turbulent wake region is due to high eddy viscosity.
(a)
(b)
(c)
Figure 4 shows the total (i.e., integral) entropy generation rate (on the left) of the flow field and drag coefficient (on the right) curves of the airfoil for the five turbulence models under the attack angle from to . The value of the entropy generation rate and drag coefficient is different due to the different turbulence models for the calculation of the turbulence viscosity coefficient. The increase of angle-of-attack had a direct effect to the total entropy generation in the flow over both airfoils for all turbulence models. It can be seen from the figure that the variation of the drag coefficient and entropy generation rate calculated by the corresponding turbulence model has the same trend.
(a)
(b)
Figure 5 shows the relationship between the entropy generation rates in the flow field and the drag coefficient of the airfoil under different turbulence models. Because of turbulent modeling, the evaluation of shear stress in CFD is different which result in differences in drag calculation. It can be seen that the entropy generation rate is linear with the drag coefficient for all turbulence models with almost uniform slope.
4.1.2. Reynolds Number Effect
To study the relationship between the Reynolds number effect and the entropy generation, we simulated the flow over NACA0012 airfoil under different Reynolds number with turbulence model. Contours of entropy generation rate around the NACA0012 airfoil for angle-of-attack , , and at different Reynolds number can be seen in Figure 6.
(a)
(b)
(c)
Figure 7 shows the axial variation of entropy generation rate along the NACA0012 airfoil for angle-of-attack , , and under different Reynolds number. When Reynolds number is increased from to , the entropy generation rate increased correspondingly. Note that when the angle-of-attack is , the recirculation bubble occurred, which results in value of entropy generation dropping to very low values. The results show good agreement with [25].
(a)
(b)
(c)
Figure 8 shows the effect of Reynolds number on the total (i.e., integral) entropy generation of the flow field (on the left) and drag coefficient (on the right) by surface integration. The value of entropy generation rate increases with Reynolds number for the whole range of angle-of-attack. It also can be seen from the figure that the variation of the drag coefficient and entropy generation rate calculated under the corresponding Reynolds number have the same trend.
(a)
(b)
We consider a dimensionless version of (17) for the expression of drag in order to highlight the relations between the drag, total entropy generation rate, and Reynolds number:where the superscript “” means nondimensional form and average entropy generation rate can be expressed in nondimensional form as follows:
Equation (33) and Figure 9 clearly reflect the relationship between the total entropy generation in the flow field and the drag coefficient of the airfoil under different Reynolds number for NACA0012 airfoil. The entropy generation rate is linear with the drag coefficient under all Reynolds number while the slope decreases with Reynolds number.
4.1.3. Airfoil Model Testing
Moreover, in order to study the relationship between the airfoil configuration and the entropy generation, we simulated the flow over different standard series NACA airfoils under with turbulence model. Figure 10 shows a comparison between the total entropy generation and drag coefficient for five NACA airfoils. The increase of angle-of-attack has a direct effect to the total entropy generation in the flow over both airfoils. However, as shown in Figure 10, each airfoil has a different entropy generation signature for different angle-of-attack. NACA0012 airfoil yields less entropy generation than NACA0024 at low angle-of-attack, while such trend is reversed when higher angles were considered. It also can be seen from the figure that the variation of the drag coefficient and entropy generation rate calculated under the corresponding Reynolds number have the same trend.
(a)
(b)
Figure 11 shows the relationship between the entropy generation in the flow field and the drag coefficient of the five NACA airfoils. The entropy generation rate is linear with the drag coefficient for all airfoils with almost uniform slope.
In this section, the entropy generation and its relationship to drag is validated by numerical simulating the total entropy generation rate and drag coefficient of different airfoils under different Reynolds number at different attack angle. Therefore, we are confident to say that entropy generation in turbulent flow using the effective viscosity leads to a good approximation of aerodynamic drag of airfoil.
4.2. Drag Prediction Comparison
Next we chose NACA0012 airfoil as the benchmark airfoil as it is well documented in the open literature. By comparing drag from alternative numerical methods, we are able to verify and validate the computational procedure of airfoil drag based on entropy generation method. The free stream conditions are set to be at different angle-of-attack. The mesh is a single-block O-type grid which is shown in Figures 12(a) and 2(b). In order to cut an exit plane perpendicular to the freestream direction, in this paper, we adopted a multivariate interpolation using radial basis functions (RBFs) [26]. It provides a direct mapping between the control points, the surface geometry, and the locations of grid points in the CFD volume mesh. The deformed mesh of NACA0012 is plotted in Figure 12(c) when the angle-of-attack is .
(a)
(b)
(c)
In order to ensure that the numerical model is free from numerical diffusion and artificial viscosity errors, several grids are tested to estimate the number of grid elements required to establish a grid independent solution. Table 1 shows the specifications of different grids used in such test. Figure 13 plots the nondimensional normal distance from the first grid point to the wall along the airfoil. Noting that is reduced with the mesh refinement increase, the grid convergence study is summarized in Table 1 where drag coefficient values (direct integration of pressure and friction forces at the airfoil surface), as well as the entropy generation rate, are presented. All the drag coefficients are given in drag counts (one drag count is 1/10000 of drag coefficient). There was a logical decrease of the pressure drag and skin friction drag value from the coarse to the medium fine grid. This is due to a better discretization of the computational domain, which leads to a more accurate solution. For the entropy generation, an obviously increase is observed from the coarse to the fine grid. Figure 14 shows the pressure coefficient distribution on the upper and lower surfaces of the airfoil as computed by the four grids. In general, the results are very consistent and the medium-fine gird is chosen to conduct the analysis presented hereafter.
Using the obtained velocity field, drag force on the airfoil from entropy generation is calculated and compared with surface and wake integration methods. The drag can be expressed in the following forms:where is aircraft surface, is the wake line, and denotes the numerical domain.
Table 2 gives the comparison of drag values of NACA0012 airfoil between calculated and experimental data. From the table it can be seen that the calculated drag values by wake integration method increase when the integration is done near the trailing edge ( means the downstream position of wake plane after the wing). Nevertheless, with the wake line placed beyond , the predicted drag values remains nearly constant.
Figure 15 gives the comparison of drag values obtained from different methods with the measured drag. It depicts that the calculation drag values from entropy generation are much lower than both the surface integration and wake integration methods. The wake integration is performed with the wake-line placed at four chord lengths downstream of the airfoil trailing edge. At , the drag coefficient from surface integration and wake integration is 94.2 and 92.1 drag counts, which is 18.8% and 15.1% higher than the experimental value, respectively. This is due to a well known phenomenon of the over prediction of drag, present in all CFD solvers. However, if we calculate the drag coefficient from entropy generation, it brings down the number to 89.1 drag counts. Although it still overpredicts the drag, it is closer to the experimental data and much better than the traditional surface integration and wake integration methods. The test case verifies that drag can be estimated from entropy generation with very high accuracies.
5. Conclusions and Extensions
Traditionally, surface integration of the pressure and stress tensor on the body surface of aircraft, which is called surface integration, is used for the drag prediction in CFD computations. But it is pointed out that the drag computed by the near-field method has inaccuracy relating to the numerical diffusion and error. An advanced drag prediction method (wake integration) based on the momentum conservation theorem around aircraft is watched with keen interests which can compute the drag components from the surface integration on the wake plane of the downstream of aircraft. The method has the drag decomposition capability into wave, profile, and induced drag component. However, the method is closely related to the wake cross surface. In this paper, a drag prediction method based on entropy generation was integrated into a RANS-based CFD solver and an approach was developed to compute the airfoil drag via entropy generation rate in the flow field which is a volume integral method derived from far-field method by applying the divergence theorem of Gauss. Present paper has been devoted to the analysis of some RANS calculations around two-dimensional airfoils, at subsonic freestream conditions. The main objective of this paper is to compare the consistency of predicting the drag of single-element airfoils using surface integration, wake integration, and entropy generation integration. Overall, entropy generation integration shows potential as a simpler method than surface integration and wake integration for calculating drag. The main advantage of this technique is that no detailed information on the surface geometry of the configuration is needed. The results show that drag prediction using CFD and entropy generation integration is possible within engineering accuracy and that the proposed method also have the drag visualization capability in the flow field for designers to take measures to minimize drag. The main conclusions are summarized as follows:(1)The drag and entropy generation in 2D domain can be related by a linear balance equation, so the drag of the airfoil can be directly estimated from total entropy generation without other effects.(2)Entropy generation consists of two parts: heat transfer and viscous dissipation. In turbulent flow, the fluctuating velocity and temperature contribution can be accounted for by using the effective viscosity and effective thermal conductivity.(3)In the 2D turbulent flow of airfoil, most of the entropy is generated in the front region, the near wall region, and the turbulent wake region. The total entropy generation of flow field increases with the angle-of-attack and Reynolds number.(4)Drag prediction using CFD and entropy generation integration is possible within engineering accuracy. Future work should concentrate on 3D airfoil drag and compressible turbulent flow to demonstrate the universality of entropy-based method in CFD prediction of airfoil drag.
Conflicts of Interest
The authors declare that they have no conflicts of interest.