Open Access
 Issue Mechanics & Industry Volume 19, Number 3, 2018 312 17 https://doi.org/10.1051/meca/2018001 04 October 2018

## 1 Introduction

The aerodynamics study problems on a numerical way is a relatively new research area. Most previous work either theoretical [16], or numerical [3,5,716] or even experimental in wind tunnel [9] on supersonic flows around airfoils are devoted to rounded airfoils at the leading edge, that is to say a development of a detached shock wave at the leading edge. These studies are generally based on the numerical solution of the Euler equations [16] or the equation of potential speed [1,4,5].

The pointed shape of the airfoil at the leading edge gives the possibility of having an attached shock wave, where a numerical technique can be used to evaluate the aerodynamic parameters of the flow. Since the flow is supersonic in the open air, far from the presence of any other obstacle, this technique makes it possible to progressively follow the flow on the airfoil surface as a function of the parameters of the upstream flow.

Given the complexity of the methods used in the supersonic aerodynamics [16], and in particular that presented in this study, called a relaxation shock method, the authors use an analytical technique named by thin-airfoil theory [36] to evaluate approximately the flow parameters, and in particular the calculation of the aerodynamic coefficients. This method gives acceptable results for very small airfoil thicknesses and upstream Mach number.

The first study on the use of the expansion shock method is presented in reference [17]. This method is used for the calorically and thermally perfect gas. They assume in this case that the specific heat at constant pressure CP is constant and does not depend on the temperature. This approach gives acceptable results only if the three parameters M0, T0 and t/C are very small. It is our opinion that the PG model does not depend on the stagnation temperature T0.

The aim of this work is to develop a new mathematical model and to develop in this context a new numerical calculation program to determine the stagnation temperature effect on the supersonic flow around a pointed airfoil, using the equations of an oblique shock and the Prandtl Meyer expansion in the case at high temperature, calorically imperfect and thermally perfect gas, lower than the dissociation threshold of the molecules, in order to correct the PG model, and to determine the aerodynamic characteristics presented by the lift, drag and pitching moment coefficients as a function of upstream Mach number, incidence angle, airfoil shape and airfoil thickness and especially of the stagnation temperature, when they are high. In this case, the specific heat at constant pressure does not remain constant and varies with the increase in temperature, which will be taken into account in our HT model. The stagnation temperature T0 is an important parameter in our model, which will demonstrate considerable corrections to the results given by the PG model. The latter gives good results only if the values of M0, T0, α, t/C are very small. Then, given the current and future applications in supersonic aerodynamics requiring high values of M0, T0 and t/C witch exceeding respectively 2.00, 1000 K and 10.0, and can arrive at 5.00, 3500 K respectively. The PG model falls failing, and the results given by this model becomes very far from the reality, which requires to make corrections to this model, hence the interest of the application of our HT model. We can consider that the HT model becomes a generalization of the PG model. In the other word, the PG model becomes a particular case of our HT model. The PG model falls failing, and the results given by this model becomes very far from the reality, which requires to make corrections to this model, hence the interest of the application of our HT model. We can consider that the HT model becomes a generalization of the PG model, were the application of the HT model is extended to high value of T0, M0, α and t/C. The application is for air. In this case and for the CP function, one finds in references [5,8,9,18,19] a series of a tabulated values for the variation of CP according to the temperature, between 55 K and 3550 K (limit not to have a dissociation of the molecules). In this temperature margin and according to references [15,16,18,19], only the translational, rotational and chemical vibrational energies are present and are included in the total evaluation of the specific heat CP(T) of air. Other sources of energy, such as the molecular dissociation energy and the ionization energy of the atoms, are absent since the temperature margin does not exceed 3500 K. A polynomial interpolation is made to these values, after several tests, in order to find an analytic function with the variation of CP(T) with the temperature. A choice on a 9th degree polynomial is made, giving a maximum error less than 0.01%. More details are found in references [8,9].

## 2 Mathematical model at high temperature

A flow deviation on the airfoil surface may result in a compression or a Prandtl Meyer expansion. If the compression is produced in the flow, a shock wave develops at the beginning of the deviation of the obstacle as shown in Figure 1. Otherwise, expansion waves develop at the beginning of the deviation as Figure 2 shown it. In both cases, there is a perturbation of the flow which results in a change of all the flow parameters through this deviation. It can be shown that, through the shock, only the total temperature [5,8,9,15,16] is conserved. Then T0 = T01 = T02. But the total pressure and density change values. Hence / = / ≠ 1.0.

For the Prandtl Meyer expansion, the temperature and the total pressure will be preserved. We set the conditions at infinity upstream by T0, ρ0 and P0. We assume that the state equation of a perfect gas (P = ρRT) remains valid, with R = 287.102 J/(kg K). For PG model γ = 1.402 is taken [8,9,14].

 Fig. 1Compression of an angle ψ.
 Fig. 2Expansion of an angle ψ.

### 2.1 Oblique shock wave

Figure 1 shows a general diagram of the development of an oblique shock wave at the beginning of the deflection of an obstacle by an angle ψ = |θ2 − θ1| and the envisaged parameters.

For the determination of the parameters (M2, β, T2/T1, ρ2/ρ1, P2/P1, /, ΔS21) through the oblique shock, the HT model presented in references [1013] is used after making a correction to the relation between β, ψ and M1, since the authors used the equation designed for the PG model to constant CP [10], given the difficulty of finding an analytic form, which gives results far enough of reality, and that does not meet the need for HT assumptions This equation is the most interesting in the calculation of the shock parameters, since all the other parameters depend on β, ψ and M1. Then, in the quality of the results, corrections will be found to the results presented in the said references. Another modification made at the CP(T) level used in these references. It has been observed that CP(T) used exhibits a slight discontinuity in the passage of T = 1000 K with a 27% error between the function used and the tabulated values [18,19].

The relationship can be summarized as follows: (1) (2) (3) (4) (5) (6) (7) (8) (9)

Two solutions can be found depending on the value of M2, implying that all physical parameters will admit two solutions. If M2 ≥ 1.00, a weak shock is obtained. If M2 < 1.00, a strong shock is obtained. In general, the weak shock that occurs in nature.

Prior to the determination of the flow parameters by (18), we must determine the angle β corresponding to ψ and M1. Then in this work, we will determine the deviation β with high precision according to the real HT model presented in this work. Since the development of an analytic relation between β, ψ and M1 is quite complicated, we will use the relations of a normal shock wave to HT model [14].

### 2.2 Expansion of Prandtl − Meyer

The situation of the presence of a Prandtl Meyer expansion is presented in Figure 2. In this case we will have a deviation of an angle ψ = |θ2 − θ1|. The flow becomes parallel to the wall after the deviation, and the calculation of the parameters after the expansion takes place after the calculation of the new value of the Prandtl Meyer function by the following relation: (10)

The function PM at HT can be calculated by the following equation [20,21]: (11) where [8,9]: (12) (13) (14)

First, it is necessary to calculate the critical temperature T* corresponding to the Mach number M = 1.00. This temperature depends only on T0. It can be determined numerically by solving the nonlinear equation obtained from relation (12) by replacing M = 1.00 and T = T* using the bipartition algorithm [2224]. We obtain T* < T0. It is calculated once in the problem.

In equation (10), ψ and T1 are known. Then ν(T1) can be calculated from equation (11) by replacing T = T1. The evaluation of the obtained integral is done by the use of the Simpson quadrature with condensation of the nodes [20] or using the Gauss Legendre quadrature of a function having a weight term to accelerate the numerical process [21]. Let us replace the obtained result in (10) to determine the new value ν(T2). Let us replace again this value in (11). In this case we fall into an inverse problem. That is to say, it is necessary to determine the temperature T2 which gives the integral (11) equal to the value given by (10). To determine T2 from (11), a combination of the bipartition algorithm with the Gauss Legendre quadrature was used. It should be noted that T2 < T*. The precision chosen in the calculation is ε = 10−8. In this case, the bipartition algorithm is used 27 times.

Once T2 is determined, it is possible to obtain the corresponding ratio T2/T0 and the Mach number M2 by relation (12) and ρ2/ and P2/ respectively by relations (5) and (6). We will have an increase in Mach number, i.e., M2 > M1. Integration (5) is done by using the Simpson method with nodes condensation [8,9].

## 3 Numerical procedure

The aim is to determine the aerodynamic characteristics summarized by the variation of M, T/T0, ρ/ρ0, P/P0 along the airfoil surface and consequently the determination of the aerodynamic coefficients CD, CL and Cm for different airfoils.

Figures 1 and 2 shows the case of the flow deflection on the extrados. For the intrados, the opposite occurs. To group the problem in a single relation, we used the absolute value for the evaluation of the angle. In the developed program, we first compute the angle ψ =  θ2 − θ1. Then, if we are on the extrados, the expansion occurs if ψ < 0.0 and the compression occurs if ψ > 0.0. The opposite is taken into account on the intrados.

Subdividing the selected airfoil into k nodes on the upper surface and into l nodes in the lower surface as shown in Figure 3. The total number of the nodes on the chosen airfoil surface is equal to n = k + l − 2. The number of regions on the upper surface is equal to k − 1 and it is equal to l − 1 in the lower surface.

A streamline coming from infinity upstream is divided into two parts at the leading edge of the airfoil (point 1). One is through the upper surface and the other through the lower surface, so that the two lines meet again at the trailing edges. We can consider the flow calculation according to the upper part and the lower surface after the other.

Every two juxtaposed panels are connected by a node. At each node, the flow undergoes a deviation by generating either a compression or expansion, or similar continuous across it, Figure 3.

To expansion either on the upper or lower surfaces, it will be an increase in Mach number, that is to say that M2 > M1. As against a compression will decrease the Mach number M2 < M1. Since it was considered that the shock is weak (real case).

Note that these parameters are constant along the segment for the shock and expansion. The flow properties in a region on the right are function of the flow parameters in the region of the left. Since these are known. There may be a detached shock if the angle ψ exceed the angle ψmax.

The jump in entropy is zero for an expansion [16]. For a compression, the entropy is different to zero. The relation (8) gives the variation of the entropy between the passages from one adjacent segment to another. To determine the total flow leap around the airfoil, all the local entropy leaps must be summed.

We may encounter a cases where the flow on the upper and lower surfaces has no deviation, i.e., ψ = (θ2 − θ1) = 0, θ2 = θ1. In this case, the flow properties remain unchanged.

To determine the position (x, y) of a node number i on the upper surface of coordinates (, ), or the lower surface of coordinates (,) we choose a reference coordinate and a small step. The step can be constant or variable. Typically the reference coordinates is selected at the leading edge of the airfoil. By replacing the value of in the airfoil equation, is easily to determine the value of its ordinate . Similarly for the lower airfoil surface, there will by the following relations: (15) (16)

For the determination of the angle deviation θi, made by a line segment connected between the points i and i + 1 with the horizontal, in order to determine the flow angle deflection ψ, the following relationship is used: (17) with: i = 1, 2,...,k for extrados, and i = 1, 2,...,l for intrados.

The region number 0 is the upstream infinity provided free data.

The flow in the region number 1 of the upper and the lower surface is given by the region number 0 of the upstream infinity.

The region number i of the upper or lower surface is limited by the nodes number i on the left and the node number i + 1 on the right.

The pressures PEi and PIj are broken down into two components each one, as shown in Figure 4. The vertical component that represents the lift force is LEi of the upper, and LIj on the lower surface. Similarly, the horizontal components that represent the drag are DEi on the upper and DIj on the lower surface. The pitching moment of the two forces exerted on the extrados is named by mEi and the moment of the two forces exerted on the intrados is mIj. The pitching moment is calculated with respect to the point p, Figure 4. For the applications, p = O is taken. Then, on the upper surface and per unit of depth we have (18) (19) (20) with i = 1, 2,…,k − 1

On the lower surface and per unit of depth we have (21) (22) (23) with j = 1, 2,…,l − 1

The number of nodes k on the upper surface is not necessarily equal to the number of the nodes l on the upper surface. For applications, the leading edge is placed at the point O of the reference of calculation.

The total lift and drag across the airfoil are respectively considered as the sum of the forces on all segments (regions) of the upper and lower surfaces. (24) (25) (26)

The pitching moment is considered positive when it rotates counter clockwise, Figure 4.

The aerodynamic coefficients are obtained as follows: (27) (28) (29) where (30)

The surface of reference S is considered to be the airfoil chord per unit of depth.

By varying the parameters M0, t/C, α and T0 as well as the airfoil shape, it is possible to find all the possible parameters and in particular the effect of T0 when it is high. A comparison between the PG model and the HT model is performed.

 Fig. 3Airfoil discretization in some panels.
 Fig. 4Presentation of the forces on the segment.

## 4 Applications

Our applications are limited by three types of airfoils which are symmetrical lozenge, Symmetrical curved airfoil and unsymmetrical curved airfoil. The calculation of the pitching moment is made to the leading edge.

### 4.1 Symmetrical lozenge

The lozenge is shown in Figure 5. We can meet this type of airfoil for building applications wings of supersonic aircraft. The geometry is given by the maximum thickness t. It is chosen in the middle of the chord. With 0 ≤ x/C ≤ 1.

We discretize the airfoil into k = l = 3 nodes. Two regions in upper and lower surface will be sufficient to describe the flow and to have the exact solution.

 Fig. 5Symmetrical lozenge

### 4.2 Symmetrical curved airfoil

The form of this type of airfoil is shown in Figure 6. The equation of the extrados and the intrados is chosen as a polynomial of 3rd degree. Its equation is given by: (31) where 0 ≤ x/C ≤ 1. The maximum thickness t of this airfoil is at a distance x/C = 1/3 from the leading edge.

 Fig. 6Symmetrical curved airfoil.

### 4.3 Curved non symetrical airfoil

The shape of the upwards curved airfoil is shown in Figure 7. The upper and lower surface are selected from the parabolic equation witch respect 3 chosen conditions. We encountered this type of airfoil applications for blade of compressor. The equation of the extrados and the intrados are respectively chosen by: (32) With 0 ≤ x/C ≤ 1. In this case, the maximum thickness t of this airfoil lies in the middle of the chord. This type of airfoil is referred to as a skeletal airfoil. The value of tE is chosen arbitrarily greater than 0. While t/C = 0.03 for the applications.

 Fig. 7Curved non symetrical airfoil.

## 5 Error of PG model compared to HT model

For each parameter, the error given by the PG model compared to our HT model can be calculated by the following relation, for aim to compare the two models: (33)

The results were divided into eight parts. In order to obtain graphical results, we have used a discretization of k = l = 1000 points on the extrados and on the intrados. For the tabulated results, we have used a discretization of k = l = 8000 points. This discretization is chosen in such a way that there will be 5 exact digits of the solution. Figures 921 contain four curves. Curves 1–4 respectively represent the variation of the parameter at HT for T0 = 3000 K, T0 = 2000 K, T0 = 1000 K and the case PG for γ = 1.402. Figure 912 contain 8 curves. 4 curves in continuous line are to represent the variation of the selected parameter on the extrados of the airfoil with HT. The 4 dashed curves represent the variation of the selected parameter on the airfoil intrados.

In this study three airfoils have been considered as Figures 57 shown them. The third airfoil (Fig. 7) is characterized by two parameters. While the 1st and 2nd airfoils are characterized by one parameter. One can even consider airfoils with several parameters.

The results for the PG model can be found in references [17]. They are presented for comparison with the HT model.

### 6.1 Typical example

In this example, we chose three very interesting airfoils in aerodynamics. The aim is to present the effect of T0 on CD, CL, Cm and ΔS21 in a numerical way as well as the calculation of the error committed by the PG model with respect to the HT model for each value of T0.

Table 1 shows the effect of T0 on CD, CL, Cm and ΔS21 for the symmetrical lozenge of Figure 5 when α = 2.0 °C, M0 = 4.00 and t/C = 0.1, followed by the given errors of PG model on these coefficients with respect to the HT model as represented in Table 2.

Table 3 shows the effect of T0 on CD, CL, Cm and ΔS21 for the symmetric curved airfoil of Figure 6 when α = 2.0 °C, M0 = 4.00 and t/C = 0.1, followed by the given errors of PG model on these coefficients with respect to HT model as represented by Table 4.

Table 5 shows the effect of T0 on CD, CL, Cm and ΔS21 for the non-symmetric curved airfoil of Figure 7 when α = 2.0 °C, M0 = 4.00 and t/C = 0.2, followed by the given errors of PG Model on these coefficients with respect to HT model as represented by Table 6. For this third airfoil, tE/C = 0.03 was taken for the application.

We clearly notice the T0 effect on all aerodynamic parameters for the selected three airfoils. The difference between the PG and HT models increases with increasing of T0. The maximum error is noticed on the coefficient CL which can arrive at 35% when T0 = 3000 K and which can arrive at 39.38% when T0 = 3500 K for M0 = 4.00 and α = 2.00. It is noted that the shape of the airfoil also influences the difference between the two models, although the parameters α, M0 and t/C are the same for the three airfoils. The PG model determines the aerodynamic parameters with excess. It is designed to solve low T0 problems. Then, if T0 increases, the performance of the flow will be degraded with respect to the PG model.

Table 1

Effect of T0 on CD, CL, Cm and ΔS21/R for the lozenge when α = 2.0 °C, M0 = 4.00 and t/C = 0.1.

Table 2

Error of PG model compared to HT model according to Table 1.

Table 3

Effect of T0 on CD, CL, Cm and ΔS21/R for the symmetrical curved airfoil when α = 2.0 °C, M0 = 4.00 and t/C = 0.1.

Table 4

Error of PG model compared to HT model according to Table 3.

Table 5

Effect of T0 on CD, CL, Cm and ΔS21/R for the non-symmetrical curved airfoil when α = 2.0 °C, M0 = 4.00 and t/C = 0.03 and tE/C = 0.2.

Table 6

Error of PG model compared to HT model according to Table 5.

### 6.2 Variation of the flow deviation along airfoil surface

Figure 8 shows the variation in the deviation of the airfoil wall which also represents the variation of the flow deviation along the surface of the selected three airfoils. This deviation enters into the calculation of the value of ψ determining the type of the flow, whether compression or expansion, given by this deviation.

 Fig. 8Variation of the flow angle deviation along the extrados and the intrados of the selected airfoils.

### 6.3 Variation of the parameters along the airfoil surface

Figures 912 show the effect of T0 on the variation of M, T/T0, ρ/ρ0 and P/P0 along the surface of the extrados and the intrados of the three selected airfoils. In Figures 11 and 12 the ratios ρ/ρ0 and P/P0 are calculated with respect to the total local conditions of the segments. We clearly notice the effect of T0 on all thermodynamic parameters with a degradation of M and T/T0 when T0 decrease, and increase of ρ/ρ0 and P/P0 when T0 increases gradually. This influence necessarily gives the influence of T0 on all the aerodynamic coefficients CD, CL and Cm.

The variation of these parameters is related to the deviation of the wall, shown in Figure 8, along the airfoil surface. As the wall deviation takes two values on the extrados and on the intrados of lozenge, one will consequently have 2 values of these parameters on the extrados and on the intrados. While for the other two airfoils there is a continuous variation along the surface of the airfoils. Thus the airfoil shape, in addition to T0, affects the variation of the parameters on the airfoil surface.

 Fig. 9Effect of T0 on the variation of M along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00.
 Fig. 10Effect of T0 on the variation of T/T0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00.
 Fig. 11Effect of T0 on the variation of ρ/ρ0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00.
 Fig. 12Effect of T0 on the variation of P/P0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00.

### 6.4 Effect of α to HT on the aerodynamic coefficients for fixed M0 and t/C

Figures 1315 represent the effect of T0 respectively on the variation of CD, CL and Cm as a function of α for the three airfoils when M0 = 4.00. For the symmetrical airfoils, we have found a symmetry of the variation of CD, CL and Cm. While for the third airfoil, this symmetry is not present, since the airfoil is not symmetrical. We clearly notice the effect of T0 on the three parameters. To determine the actual value of CD, CL and Cm, the value obtained from Figures 1315 must be divided by 103.

In Figures 1315, the range of variation of the incidence angle α is taken [−10 °C, 10 °C]. Generally beyond this interval, we will have the phenomenon of stall.

It is very important to determine an angle of incidence which makes it possible to find CL = 0, called angle of zero lift, and the incidence angle making it possible to find Cm = 0, called zero moment angle. For symmetrical airfoils, this angle is equal to 0.0. Whereas for the non-symmetric airfoils, force will have a value other than zero and which depends on T0. Table 7 shows the effect of T0 on the zero lift angle for the third non-symmetric airfoil and Table 8 shows the effect of T0 on the zero moment angle for the third airfoil for some values of M0 when t/C = 0.03 and tE/C = 0.1. We notice the effect of T0, M0, t/C, the shape of airfoil and tE/C on these two remarkable angles.

It is clear that for an incidence angle α in the vicinity of the angle of zero lift, or the angle of zero moment, of one degree, there is no influence of T0 on these parameters.

 Fig. 13Effect of T0 on the variation of CD as a function of α.
 Fig. 14Effect of T0 on the variation of CL as a function of α.
 Fig. 15Effect of T0 on the variation of Cm as a function of α.
Table 7

Effect of T0 the zero lift angle for the third non-symmetric airfoil in function of M0.

Table 8

Effect of T0 the zero moment angle for the non-symmetric airfoil in function of M0.

### 6.5 Effect of M0 at HT on the aerodynamic coefficients for fixed α and t/C

Figures 1618 represent the effect of T0 respectively on the variation of the aerodynamic coefficients CD, CL and Cm as a function of M0 for the three selected airfoils when α = 2.0 °C. The presentation is done on a Logarithmic scale because of the fact that the small values are grouped in a Figure with the large values of the aerodynamic coefficients. We clearly notice the effect of T0 on these coefficients. This effect becomes important as M0 increases gradually. Then for small values of M0 up to about M0 < 2.00, the difference between the PG and HT models is not significant, which shows that the results given by the PG model are acceptable. But if M0 increases, independently of T0, the corrections made by the HT model are necessary, which shows the interest of the HT model for the large values of M0 > 2.00. Generally this limit depends on the error considered in the calculation. It should also be noted that if M0 decreases, there will be the appearance of a detached shock wave which develops at the leading edge of the airfoil. Commentaries remain valid for any values of α, t/C and tE/C.

The variation of M0 in Figures 1618 starts from up to 5.0. We note the effect of T0 on the minimum value of the Mach number upstream , that can have the flow to limit the attached shock with the detached shock, respectively for the three selected airfoils. Forcing this limit depends on T0, α, t/C and tE/C. The limit given by the PG model is considered as an attached shock for the HT model, since, that is to say one can have a Mach number M0 <  for HT with the possibility of finding an attached shock. It is noted that the margin of finding a detached shock is greater for the third non-symmetric airfoil than the two other. Then not only depends on α, t/C, but also depends on the shape of the airfoil.

 Fig. 16Effect of T0 on the variation of CD as a function of M0.
 Fig. 17Effect of T0 on the variation of CL as a function of M0.
 Fig. 18Effect of T0 on the variation of Cm as a function of M0.

### 6.6 Effect of t/C at HT on the aerodynamic coefficients for fixed M0 and α.

Figures 1921 represent the effect of T0 respectively on the variation of CD, CL and Cm as a function of the thickness t/C of the three airfoils when M0 = 4.00 and α = 2.00. For the non-symmetrical airfoil, the curvature of the extrados was varied without touching the airfoil thickness taken at t/C = 0.03. We note the effect of T0 on these parameters, despite the small thickness, which requires the use of the HT model for corrections despite at low t/C. It is also noted that if T0 decreases, the HT model becomes non-comparable with the PG model. Then when T0 < 240 K, the PG model gives acceptable results. It should also be noted that the shape of the airfoil will still influences the variation of these coefficients.

It is noted that if t/C increases, it will be possible to find a detached shock. If t/C ≤ (t/C)max, there will necessarily be an attached shock, hence a solution can be given by the developed program. If t/C > (t/C)max, we will have the appearance of a detached shock wave.

 Fig. 19Effect of T0 on the variation of CD as a function of t/C.
 Fig. 20Effect of T0 on the variation of CL versus t/C.
 Fig. 21Effect of T0 on the variation of Cm versus t/C.

### 6.7 Variation of the aerodynamic coefficients as a function of T0 for fixed α, M0 and t/C

Figures 2224 represent the variation of the aerodynamic coefficients of the three selected airfoils as a function of T0 and the comparison with the PG model. The latter does not depend on T0, where it is represented by a horizontal straight line along the entire interval. The CD, CL and Cm variation is chosen for M0 = 4.00 and α = 2.00 °C. We can note the effect of T0 on these parameters. Hence the need to use the HT model for possible corrections. One notices again when T0 is small, the HT model becomes confounded with the PG model. We can go to about 240 K. Then when T0 < 240 K, the PG model gives acceptable results. It is also noted that the shape of the airfoil still influences the variation of these coefficients since the variation of the aerodynamic coefficients of the three airfoils is not the same. We note that the HT model degrades the CD, CL and Cm values due to their decrease with increasing of T0, which does not work with the physical and real behavior of the flow.

 Fig. 22Variation of CD as a function of T0.
 Fig. 23Variation of CL as a function of T0.
 Fig. 24Variation of Cm as a function of T0.

### 6.8 Variation of the error caused by the PG model compared to HT model as a function of M0

Figures 2527 represent the variation of the relative error caused by the use of the PG model with respect to the HT model on the aerodynamic coefficients CD, CL and Cm of the three selected airfoils. The application is made for the temperatures T0 = 1000 K, 2000 K and 3000 K and for α = 2.00 °C and t/C = 0.1. It can be seen that the error increases considerably and can reach 55% when T0 = 3000 K and M0 = 5.00. This error is noticed for the coefficients CL and Cm for the symmetrical curved airfoil. This value also increases with the airfoil shape, t/C and α. This value of the maximum error give the obligation to use the HT model for possible corrections to the results given by the PG model when T0 is high and begins to exceed the 240 K approximately.

 Fig. 25Variation of the relative error caused by the PG model compared to HT model on the coefficient CD versus M0.
 Fig. 26Variation of the relative error caused by the PG model compared to HT model on the coefficient CL versus M0.
 Fig. 27Variation of the relative error caused by the PG model compared to HT model on the coefficient Cm versus M0.

## 7 Conclusion

The conclusions that can be deducted are:

A detached shock occurs for lower M0.

Detached shock still occurs when the angle ψ exceeds a certain maximum angle.

The computational accuracy for the PG and HT models depends on the discretization, which results in the choice of k and l for the curved shapes. More k and l will be high, we will have a good accuracy.

The PG model gives acceptable results for small values of M0 < 2.00, T0 < 240 K and t/C < 1.0 approximately. On the other hand, when M0 or T0 or t/C increases, the PG model gives results which are different from the real case, hence the need for the HT model.

The developed numerical program can process any gas found in nature. In this case, we must add the variation of the specific heat CP(T) and the constant R of the gas with the calculation of H(T).

The convergence of the results requires an additional calculation time for the HT model compared to the PG model for the same accuracy.

A condensation of the nodes in the Simpson quadrature is used to integrate the function with high precision in a reduced time.

For M0 given, there is a maximum deflection of the airfoil to avoid the detached shock. This limit also depends on T0.

For each airfoil geometry, there is a minimum value of the upstream Mach number to avoid the detached shock. This limit depends on T0, α and the gas used.

The T0 is an essential parameter of our HT model. The PG model results do not depend on T0.

The stagnation temperature T0 degrades the CD, CL and Cm parameters compared to the PG models. This difference increases with the increase in T0.

The maximum error is noticed on the coefficient CL which can arrive at 39.38% when T0 = 3500 K for the selected airfoils when M0 = 4.00 and α = 2.00 and can reach 59% when T0 = 3500 K for M0 = 5.00 and α = 2.00.

For t/C < 1.0 or M0 < 2.00 or T0 < 240 K, the PG model can be used to evaluate the flow independently of the values of M0 and T0. But if t/C > 1.0 or M0 > 2.00 or T0 > 240 K, the corrections given by the HT model become necessary to evaluate the flow parameters accurately.

The PG model use limit in terms of maximum values of M0, t/C, and T0 is set to the required accuracy.

The presentation of the results on the choice of three airfoils differs. The developed program can do the calculation for any airfoil shape.

Not only do the parameters α, T0, M0, t/C influence the aerodynamic coefficients, the airfoil geometry also influences.

The flow around an airfoil is characterized by the generation of a shock wave at the leading edge in addition to the possibility of having a progressive shock through the surface of a certain airfoil, like the third airfoil.

The flow around an airfoil is characterized by an entropy jump due to developed shock on the airfoil.

It can be considered that the work carried out can be considered as a numercial wind tunnel. It allows to numerically validate the new HT model with the old existing PG model that is experimentally validated.

As a perspective, this problem can be studied for the airfoils having a slat and ailerons to have a more maneuverability for the variation of the aerodynamic coefficients.

## Nomenclature

θ: Deviation angle of an airfoil segment

ψ: Flow angle deviation

M: Mach number

β: Shock wave deviation

µ: Mach angle

γ: Specific heats ratio

R: Thermodynamic constant of air

CP: Specific heat at constant pressure

ν: Prandtl-Meyer function

α: Angle of incidence of the airfoil

D: Drag force

L: Lift force

m: Pitchnig moment

q0: Dynamic pressure at upstream infinity

S: Reference surface

P: Pressure

T: Temperature

ρ: Density

C: Airfoil chord

Cm: Pitching moment coefficient

CD: Drag coefficient

CL: Lift coefficient

t: Maximum thickness of the airfoil

x, y: Position of the point on the airfoil.

ε: Error of computation

fE: Equation of the airfoil extrados

fI: Equation of the airfoil Intrados

k: Nodes number on the extrados

l: Nodes number on the intrados

n: Total number of node on the airfoil

PG: Perfect gas

HT: High temperature

PM: Prandtl Meyer

ΔS: Total variation of entropy

## Subscripts

0: Upstream condition

1: Upstream state at the panel under consideration

2: Considered panel

p: Point for the calculation of the pitching moment

## Acknowledgments

The authors acknowledges Khaoula, AbdelGhani Amine, Ritadj and Assil Zebbiche and Mouza Ouahiba for granting time to prepare this manuscript.

## References

1. J. Moran, Introduction to theoretical and computational aerodynamics, John Wiley & Sons, 1984 [Google Scholar]
2. I.H. Abbott, A.E. VonDoenhoff, Theory of wing sections including summary of airfoil, Dover Publication, Inc., New York, 1959 [Google Scholar]
3. G. Emanuel, Gasdynamic: theory and application, AIAA Educational Series, New York, 1986 [CrossRef] [Google Scholar]
4. J.D.Jr. Anderson, Fundamentals of aerodynamics, 2nd Edn, McGraw-Hill Book Company, New York, USA, 1988 [Google Scholar]
5. J.D.Jr. Anderson, Modern compressible flow with historical perspective, 2nd Edn, McGraw-Hill Book Company, New York, USA, 1982 [Google Scholar]
6. J.D.Jr. Anderson, Hypersonic and high temperature gas dynamics, McGraw-Hill Book Company, New York, 1989 [Google Scholar]
7. M.J. Zucro, J.D. Hoffman, Gas dynamics, Wiley, New York, Vol. 1 and 2 , 1976 [Google Scholar]
8. T. Zebbiche, Z. Youbi, Effect of stagnation temperature on the supersonic flow parameters with application for air in nozzles, Aeronaut. J. 111 (2007) 31–40 [CrossRef] [Google Scholar]
9. T. Zebbiche, Z. Youbi, Supersonic flow parameters at high temperature. application for air in nozzles, German Aerospace Congress 2005, DGLR-2005-0256, Friendrichshafen, Germany, 2005 [Google Scholar]
10. E.T. Kenneth, Computation of thermally perfect properties of oblique shcok waves, NASA CR- 4749, 1996 [Google Scholar]
11. E.T. Kenneth, Computation of thermally perfect oblique shcok waves properties, AIAA-97-0868, 35th Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings, 1997 [Google Scholar]
12. K.E. Tatum, Computation of thermally perfect oblique shock wave properties, 35th Aerospace Sciences Meeting and Exhibit, Reno, AIAA-97-0868, 1997 [Google Scholar]
13. S.M. Yahya, Gas tables for compressible flow calculation, Fifth edition, New Age International Publishers, New Delhi, 2006 [Google Scholar]
14. T. Zebbiche, Effect of stagnation temperature on the normal shock wave, Int. J. Aeronaut. Space Sci. 10 (2009) 1–14 [CrossRef] [Google Scholar]
15. C.R. Peterson, P.G. Hill, Mechanics and thermodynamics of propulsion, Addition-Wesley Publishing Company Inc., New York, USA, 1965 [Google Scholar]
16. G.P. Sutton, O. Biblarz, Rocket propulsion elements, 8ème Édn, John Wiley and Sons, 2010 [Google Scholar]
17. T. Zebbiche, M. Salhi, M. Boun-jad, Numerical computation of supersonic flow around a pointed airfoils, J. Comput. Methods Sci. Eng. 12 (2012) 213–233. [Google Scholar]
18. G.J. Van Wylen, Fundamentals of classical thermodynamics, John Wiley and Sons, Inc., 1973 [Google Scholar]
19. B.J. McBride, S. Gordon, M.A. Reno, Thermodynamic Data for Fifty Reference Elements, NASA TP- 3287, 1993 [Google Scholar]
20. T. Zebbiche, Stagnation temperature effect on the Prandtl Meyer function, AIAA J. 45 (2007) 952–954 [CrossRef] [Google Scholar]
21. T. Zebbiche, M. Boun-jad, Numerical quadrature for the Prandtl–Meyer function at high temperature with application for air, Thermophys. Aeromech. 19 (2012) 381–384 [CrossRef] [Google Scholar]
22. A. Raltson, A. Rabinowitz, A First Course in Numerical Analysis, McGraw Hill Book Company, 1985 [Google Scholar]
23. A. Iserles, A first course in the numerical analysis of differential equations, Cambridge University Press, 1996 [Google Scholar]
24. B. Démidovitch, I. Maron, Eléments de calcul numérique, Editions MIR, Moscow, Russia, 1987 [Google Scholar]

Cite this article as: R. Takhnouni, T. Zebbiche, A. Allali, Stagnation temperature effect on the supersonic flow around pointed airfoils with application for air, Mechanics & Industry 19, 312 (2018)

## All Tables

Table 1

Effect of T0 on CD, CL, Cm and ΔS21/R for the lozenge when α = 2.0 °C, M0 = 4.00 and t/C = 0.1.

Table 2

Error of PG model compared to HT model according to Table 1.

Table 3

Effect of T0 on CD, CL, Cm and ΔS21/R for the symmetrical curved airfoil when α = 2.0 °C, M0 = 4.00 and t/C = 0.1.

Table 4

Error of PG model compared to HT model according to Table 3.

Table 5

Effect of T0 on CD, CL, Cm and ΔS21/R for the non-symmetrical curved airfoil when α = 2.0 °C, M0 = 4.00 and t/C = 0.03 and tE/C = 0.2.

Table 6

Error of PG model compared to HT model according to Table 5.

Table 7

Effect of T0 the zero lift angle for the third non-symmetric airfoil in function of M0.

Table 8

Effect of T0 the zero moment angle for the non-symmetric airfoil in function of M0.

## All Figures

 Fig. 1Compression of an angle ψ. In the text
 Fig. 2Expansion of an angle ψ. In the text
 Fig. 3Airfoil discretization in some panels. In the text
 Fig. 4Presentation of the forces on the segment. In the text
 Fig. 5Symmetrical lozenge In the text
 Fig. 6Symmetrical curved airfoil. In the text
 Fig. 7Curved non symetrical airfoil. In the text
 Fig. 8Variation of the flow angle deviation along the extrados and the intrados of the selected airfoils. In the text
 Fig. 9Effect of T0 on the variation of M along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00. In the text
 Fig. 10Effect of T0 on the variation of T/T0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00. In the text
 Fig. 11Effect of T0 on the variation of ρ/ρ0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00. In the text
 Fig. 12Effect of T0 on the variation of P/P0 along the surface of the three selected airfoils when α = 2.0 °C and M0 = 4.00. In the text
 Fig. 13Effect of T0 on the variation of CD as a function of α. In the text
 Fig. 14Effect of T0 on the variation of CL as a function of α. In the text
 Fig. 15Effect of T0 on the variation of Cm as a function of α. In the text
 Fig. 16Effect of T0 on the variation of CD as a function of M0. In the text
 Fig. 17Effect of T0 on the variation of CL as a function of M0. In the text
 Fig. 18Effect of T0 on the variation of Cm as a function of M0. In the text
 Fig. 19Effect of T0 on the variation of CD as a function of t/C. In the text
 Fig. 20Effect of T0 on the variation of CL versus t/C. In the text
 Fig. 21Effect of T0 on the variation of Cm versus t/C. In the text
 Fig. 22Variation of CD as a function of T0. In the text
 Fig. 23Variation of CL as a function of T0. In the text
 Fig. 24Variation of Cm as a function of T0. In the text
 Fig. 25Variation of the relative error caused by the PG model compared to HT model on the coefficient CD versus M0. In the text
 Fig. 26Variation of the relative error caused by the PG model compared to HT model on the coefficient CL versus M0. In the text
 Fig. 27Variation of the relative error caused by the PG model compared to HT model on the coefficient Cm versus M0. In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.