Mathematical Problems in Engineering

Volume 2018, Article ID 4347650, 15 pages

https://doi.org/10.1155/2018/4347650

## CFD Prediction of Airfoil Drag in Viscous Flow Using the Entropy Generation Method

School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China

Correspondence should be addressed to Jun Wang; nc.ude.tsuh@tsuhjgnaw

Received 15 January 2018; Revised 6 March 2018; Accepted 13 March 2018; Published 16 May 2018

Academic Editor: Eusebio Valero

Copyright © 2018 Wei Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.