Analysis of One-Dimensional Inviscid and Two-Dimensional Viscous Flows Using Entropy Preserving Method

In this paper, the entropy preserving (EP) scheme (which is introduced recently by Jameson) has been considered deeply and compared with the other artiﬁcial viscosity and upwind schemes. The discretization of the governing equations in the EP scheme is performed in such a way that the entropy is conserved in all those points with no shock. The purpose of this study was to introduce a stable numerical method that enters a minimum artiﬁcial dissipation only in the vicinity of shocks. In this paper, an inviscid one-dimensional ﬂow through a convergent–divergent nozzle and a viscous two-dimensional ﬂow with axial symmetry are considered. It is shown that the EP scheme is more accurate if the number of mesh points is increased; and in contrast to other schemes, there is no limit in increasing the number of points.


List of symbols
Indicating the derivative with respect to the desired parameter ∞ The infinity variable

Introduction
Having studied the history of computational fluid dynamics methods, so many efforts were taken to invest new methods with lower fluctuation and lower numerical error. Weak solutions of conservative laws (solutions with discontinuity like compression wave) with applied initial conditions may do not have unique solution [1]. To eliminate the irrelevant physical solutions, additional conditions must be applied. For nonlinear problems with discontinuity points (like shocks), it is preferred to solve equations conservatively [2]. Based on Lax and Wandroff theory [3], if mesh is generated in such a way that distances between points goes to zero, uncoupled governing equations will predict shock conditions well with low fluctuations. On the other hand, conserving of quantities such as mass and linear momentum in uncoupled equations is of high importance. For example, small errors in conserving mass, which is more evident in shock problems, may cause enormous errors in the solution. Adding the artificial dissipation terms also cast doubt on the mass conservation law. So, an approach that adds less dissipation to the equations is more reliable for enforcing conservation laws as better solutions can be predicted using them.
It seems that if the governing equations estimate the amount of energy correctly, it will be possible to construct a discretization method such that estimate the same amount of energy with a stable solution, while there is no need to add any artificial dissipation as long as the solution remains smooth. Using the energy estimation for the stability of discrete methods for initial value problems has a long history. The energy method is discussed by Richtmyer and Morton [4].
Capturing discontinuities without producing fluctuations still seems to be necessary. Structure of the shock operators for the Burger equations and the equations of gas dynamics have been discussed by Gustafsson and Olsson [5]. Shock capturing methods have been widely studied by many researchers. For example, we may refer to Godunov [6], Boris and book [7], Van Leer [8], Roe [9], Harten [10], Liou [11] and Jameson [12]. They typically added a form of explicit or implicit artificial dissipation terms to the equations, or the discretization of the equations was done in a way that the scheme introduces a strong diffusion property (such as upwind and CUSP schemes). These schemes may reduce the accuracy of the viscous flow simulations considerably. Therefore, the objective is to introduce a way which adds the dissipation term to the equations only in the vicinity of discontinuity points [13].
There are some discussion to define the entropy function h(u) as a function of u. Harten [14] determines conditions that the function h(u) is a convex function such that if h(u) remains finite, the solutions also remain limited. Multiplying the governing equations by ∂h ∂t = w T provides an improved equation in which w T represents the entropy variable. In fact, the change of variable u to w is used to create symmetric variables for smoothing the solution. Entropy variables have been studied by Hughes, Franca and Mallet [15] and also by Gerritsen and Olsson [16].
Analysis of the optimization in terms of entropy production or entropy preserving (EP) has been widely studied in the recent years. There are different types for studying the entropy production. For example, entropy generation minimization method can be cited. Bejan [17] used this method to study the actual devices (irreversible devices). Bejan [18] also has reviewed the evolution of classical thermodynamics to entropy generation minimization. In the recent years, many researches have been done in this theme. For example, it could be mentioned to the works have been done by Ben-Mansour and Sahin [19], Esmail and Mokheimer [20] and Slimi and Saati [21].
The present paper is based on EP method, proposed by Jameson [2]. We briefly introduce EP method and also some artificial dissipation and upwind schemes, which are implemented to compare with EP method. The EP method discretizes equations in such a way that the rate of the total energy variation be the same as their real value. It also fixes its value in cases which total entropy is constant. Governing equations are valid in continuous domain and near the discontinuity points (like shock), energy or entropy difference values must be calculated using specific operators.
In the subsequent section, discretization method of equations is expressed for one-dimensional scalar conservation equation. Then, the numerical comparison between EP method and some artificial dissipation or upwind methods will be presented. To this end, inviscid one-dimensional flow in a convergent-divergent nozzle is considered. And then, viscous two-dimensional flows with axial symmetric on a flat plate is investigated, and finally, a general overview is presented.

Introduction of Some Artificial Dissipation and Upwind Schemes
In this section, several artificial dissipation and upwind schemes will be explained briefly. Scalar Artificial Dissipation Scheme (SCDS): In 1981, Jameson et al. [22] proposed an artificial method that has been very effective in practice and became popular as scalar artificial dissipation method. This method is based on the supposition that fourth-order artificial dissipation is added in all studied region to prevent nonlinear instability. Next, to prevent the fluctuations near shocks, second-order artificial dissipation term is added locally. However, this method was really successful on its time, but it adds high value of artificial viscosity to the equations.
Matrix Artificial Dissipation Method (MADS): Swanson and Turkel [23] suggested another method based on so-called scalar one, named the matrix artificial dissipation. In this method, a matrix is used to estimate the necessary dissipation coefficient in the system of governing equations instead of single value. However, using this method, significant improvement in the accuracy of results in the vicinity of compression waves is obtainable, but the time needed to estimate the dissipation term is increased significantly. The purpose of this scheme is to apply appropriate value of dissipation in each flow equation such that the dissipation is reduced desirably.
CUSP scheme: This scheme is based on adding dissipations by separating the pressure terms in the flow flux relations. Jameson tried to reduce the complexity of calculations and the required time, while acceptable solution can be obtainable [12].
Roe Upwind Scheme: In this approach, initial parameters on the cell surface are calculated by combining saved information in two sides of cell surfaces by density-based averaging. This approach was introduced by Roe in 1981 [9]. EP method will be discussed in the next section.

Introducing the Entropy Preserving Method
This section will introduce the EP method and shows how to discretize the equations. It should be noted that governing equations were written in references [1] and [2]. In this section, derivation of the available equations is shown in more details comparing to the cited references.

One-Dimensional Scalar Conservative Law
The approach used in this article is the EP scheme. For onedimensional scalar conservation law is: which becomes discretized as following: The rate of energy changes which is predicted by discrete equation will be the same as the exact solution only if the amount of flux is calculated using the following equations [1] In the above equations, u is the state variable, f is flux, and t and x are the time and the distance between two adjacent points, respectively. The index j denotes the studied point and θ is a variable between zero and one which is an integration point. For viscous flows, it is seen that shock waves could be removed using more refined mesh. In this case, the cell Reynolds number must be less than 2 [24]. Finally, the best time step method should be used. The amount of energy or any other quantity must be preserved using an implicit time step method such as Crank-Nicolson [25]. To calculate the state vector, we can use the average between the start and the end of one time step as: For validity of Eq. (5), internal marching in each interval or each time step must be implemented. To avoid complicated programming and reduce the costs, Shu's total variation diminishing scheme [26] has been used. For the following time step equation in which R(u) is the residual quantity, the following multistage time step method is used to solve the problem.
in which u (0) indicates the amount of u parameter, in first time step. u (1) and u (2) are its value in the second and third steps of the mentioned method, and u (3) indicates the u parameter value at the end of time step. The general form of the governing equation is shown in Eq. (1). If F defines in such a way that (8) in which subscript u denotes derivation of function with respect to u. It could be shown [1] that the rate of change in total energy is calculated as Suppose that there would be a G(u) function in such a way that its derivative with respect to u equals to f function.
F and G can be obtained as follows [2]: Assuming that Eq. (1) is discretized on a uniform mesh based on the Eq. (2). If there is a G function so as If G u j+ 1 2 can be obtained using the following equation, and G(u j ) is denoted by G j , the rate of total energy variation is calculated as follows [1]: which is the same as what obtains from Eq. (9). The above conditions are valid only if Eq. (13) is valid. The essential condition for conserving this equation is to calculate f flux in the form of Eqs. (3) and (4). It will subsequently be shown how Eqs. (3) and (4) satisfy Eq. (13).
As it can be interpreted by Eq. (4) if θ value be zero and one, the termû(θ ) equals to u j and u j+1 , respectively. So it could be written as Alternatively, using the chain rule, it can be written as Using Eqs. (15) and (16) In the above equations, indices denote derivatives with respect to the studied variable. By considering Eq. (4) and substituting the u θ , the Eq. (16) gets the following form By comparing Eqs. (17) and (19) Considering Eq. (12), the term G u is substituted by its equivalent term. Then, the above integral is substituted with its equivalent from Eq. (3).
which are the same conditions as Eq. (13). Equations (3) and (4) are essential equations forming the EP method. In fact, a method is called EP if these formulas are used for computing the flux. To calculate the integration (3), the Lobatto Rule is used which uses the boundary points and n − 2 internal points for exact solution of polynomials of order 2n − 3. If n = 3, Simpson three points formula is obtained [2].

Results and Discussion
In this section, results obtained from EP method is presented and compared with previous upwind and artificial dissipation schemes. First, results for an inviscid one-dimensional flow in a convergent-divergent nozzle has been considered, and then, the viscous flow over a flat plate both for two-dimensional and axisymmetric flows has been considered. It should be noted that in all cases, EP method adds minimum artificial dissipation term to the equations.  Figure 1 shows the studied geometry. This case is 1D flow, and a uniform grid is considered to solve the flow field. Number of grid points varies for each section and is mentioned at the beginning of these sections. The convergent-divergent geometry of nuzzle is as follows: In the above equation, y max = 0.75 and the value of y min are determined by considering the area ratio. If area of inlet, exit and throat face are shown by A 1 , A 2 and A * , respectively, velocity could be controlled by the ratio of back pressure to the total pressure or A 2 /A * ratio. If the ratio of back pressure to inlet pressure be 0.75 and A 2 /A * be 1.5, Mach number in the nozzle would be about 1.6 and a normal shock in the nozzle would happen. Figure 2 shows the changes of Mach number in the nozzle obtained using EP method and other schemes near the shock.
It can be observed that EP method has low fluctuation before the shock. First-order Roe method shows that the shock wave is too thick. The third-order Roe method shows the shock is sharper without considerable fluctuation. CUSP scheme performance is acceptable in determining the place of shock and its thickness. However, EP method shows the shock place better.
Theory of gas dynamics leads to the shock place in x = 0.7184. The Mach number before and after the shock is 1.61 and 0.6655, respectively. Table 1 shows how shock place and Mach numbers can be predicted before and after the shock for different schemes.
This table shows that the predicted shock thickness obtained by using the start and the end place of the shock is less for CUSP method. Furthermore, the shock place and Mach numbers are predicted better by EP method.
As could be seen in Table 1, EP method predicts the shock place with more precision compared with other schemes, and shock thickness in this method is better compared with scalar and first-order Roe methods. It shows that calculations obtained by EP method are more precise before the shock, and this is the property of this method. Figure 3 compares the results obtained by different methods before the shock. As mentioned earlier, EP method is more precise before the shock.   Fig. 3 Mach number before the shock for pressure ratio of 0.75 and area ratio of 1.5 (n = 100) Figure 4 shows the effect of number of mesh points on the results of EP method (back pressure to inlet pressure ratio is 0.75 and A 2 /A * is 1.5). As it was expected, the finer the mesh, the higher accuracy is obtained. On the other hand, with other methods by increasing the number of points, the solution becomes unstable. For instance, in Fig 5, this effect is shown for scalar method (back pressure to inlet pressure ratio is 0.7 and A 2 /A * is 1.5). As could be seen, for EP method for n = 6,000, the solution is converged yet, but in scalar method points more than 800, the solution is not converged. (To show this point, the convergence history diagram for scalar method is plotted in Fig. 6). So in Fig 5, results were presented for mesh points up to this number of points. Besides, by considering Figs. 4 and 5, more fluctuations are observed for increased number of points in scalar method compared with EP method.
By considering the first law of thermodynamics and continuity equation, in an adiabatic steady flow, the total enthalpy and mass flow rate are constant. But in numerical methods, the value of mentioned parameters changes as a result of numerical errors. In subsequent sections, these parameter variations will be discussed. First, consider the case that shock does not happen in the nozzle. In this case, the rate of back pressure to inlet pressure is 0.93 and A 2 /A * is considered to be 1.5. Figure 7 shows the total enthalpy changes for different methods. It could be observed in Fig. 7 that differences in total enthalpy is less in EP and CUSP methods compared with the others, and this property could be considered as a positive point for them. In EP method, the enthalpy difference is the least in the final part of the chart. Figure 8 shows the mass flow rate difference with respect to position for pressure ratio of 0.6 and area ratio of 3 in case of occurrence of an intense shock. As higher mass flow rate indicates higher value of added artificial dissipation in the mentioned method [27], it could be inferred that the scalar method adds the most amount of artificial dissipation. In addition, the fluctuations visible in this method shows the inappropriate amount of added artificial dissipation. However, EP method shows the least mass flow rate differences compared with other methods. Figure 9 compares two different methods with low number of points (n = 50). Fluctuations of the solution in EP method are less, compared with scalar method. Figure 10 does the same comparison for the case with very high number of grid points. It can be seen that in scalar method, when the number of points exceeds n = 1,200, the solution does not converge. However, in EP method, not only it is possible to increase the number of points up to more than 6,000, but still it does not diverge and results in much more accurate predictions.
Finally, as was shown in the previous figures, EP method predicts reasonable results for low number of points. Further, by increasing the number of points, it also tends to lead to more precise predictions. It seems EP method is a more appropriate model of flux calculations for real turbulent flows. In this section, a comparison between EP and scalar method for both low and very high number of grid points is illustrated. In conclusion, for inviscid one-dimensional flow, EP method is able to give results with higher precision in the shock domain by increasing the number of points and choosing an appropriate artificial dissipation term.

Viscous Two-Dimensional Flow Over a Flat Plate
In this section, viscous flow in two-dimensional Cartesian and axisymmetric coordinate has been studied, and then, a comparison between EP method and previous upwind and artificial dissipation methods has been presented. The geometry considered herein is a flat plate. In this section, the variation of dimensionless velocity with respect to dimensionless vertical position for section of plate in x = 3.4 is studied. Two turbulent models, namely K-L and Cebeci-Smith, are used. The airflow is assumed to be with Mach number 0.4 and Reynolds number 1,360,000 at one unit of dimension. The amount of specific heat capacity is C P = 1,004. J kg.K and Prandtl number is Pr = 0.72. Using the Sutherland relation [2], molecular viscosity of air is calculated as follows: where T ∞ is the free stream temperature. Figure 11 shows the studied geometry and mesh generation. Boundary conditions are defined in the following form. In the inlet, velocity, pressure and density are considered equal to the velocity, pressure and density of free stream. Values of velocity and density in the exit cross section are considered to be the same as their values in the points next to them in the domain, and pressure values in the exit cross section are equaled to the pressure of free stream. For the points at the bottom, values of density and pressure are set to be equal to their upper points, and velocity in the vertical direction of these points is considered to be zero. Horizontal velocity over the flat plate (from x = 0 to the end) is considered to be zero. The generated mesh to study the results is 150×100 points. In this section, turbulent flow with different turbulence models was studied to compare EP method with the previous upwind and artificial dissipation scheme. Figures 12, 13  and EP method was compared with upwind and artificial dissipation scheme and theoretical results [28].
By focusing on Figs. 12, 13, 14, 15, it is obvious that scalar method is less useful to detect the velocity domain compared with other methods. CUSP method has discontinuity in velocity distribution in some regions. In conclusion, EP method presents more accurate results compared with other methods Finally, two points should be mentioned. As was stated before, EP method presents more accurate results by increasing the number of points. In this method, there is no limitation for this increase. Since this method is practically designed for computational grid with high number of points, in this paper, the independence of the numerical results from the mesh generation is not the attention case. Numerical results show that flux calculation in EP method is done in such a way that presents more accurate results.

Conclusion
In this paper, EP method was investigated and the relations were re-derived in more details. EP method uses Eqs. (3) and (4).
Jameson [2] proved that by this little difference in flux calculation, better results could be obtained. In this research work, it is found that using this scheme, although no artificial dissipation term is added to the equations, it causes stable solutions with low fluctuations. But its lower convergence rate is considered as its deficiency. It must be noted that in EP method, like most of previous methods, dissipation terms must be added to the domains with high gradients of pressure (like shock); but unlike the previous methods, it is not necessary to add these terms in the domains without high pressure gradients. Obtained results in this article shows that EP method presents very good results for points without shock even with very small number of points. Of course, better results are obtainable by refining the mesh in comparison with other methods. The constraint for the number of grid points is usually very less for EP method, compared with other schemes, which become unstable very earlier.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.