On the Kutta Condition in Compressible Flow over Isolated Airfoils

: This paper presents a novel and accurate method to implement the Kutta condition in solving subsonic (subcritical) inviscid isentropic compressible ﬂow over isolated airfoils using the stream function equation. The proposed method relies on body-ﬁtted grid generation and solving the stream function equation for compressible ﬂows in computational domain using ﬁnite-di ﬀ erence method. An expression is derived for implementing the Kutta condition for the airfoils with both ﬁnite angles and cusped trailing edges. A comparison of the results obtained from the proposed numerical method and the results from experimental and other numerical methods reveals that they are in excellent agreement, which conﬁrms the accuracy and correctness of the proposed method. are run on a PC with Intel Core i5 and 6G RAM. A tolerance of 10 − 7 is used in iterative loops to increase the accuracy of results. The computation time is about 5 min which is high due to the iterative solution of the elliptic grid generation and the stream function equations.


Introduction
Nowadays, computational fluid dynamics complements its experimental and theoretical counterpart. Benefiting from high-speed digital computers, the use of sophisticated numerical methods has made possible the numerical solution of fluid flow problems which were heretofore intractable. Compressible flows over airfoils and wings play a vital role in computational aerodynamics, and demand advanced computational techniques. Since introducing the source and vortex panel methods in the late 1960s [1], they have become the standard tools to numerically solve low-speed flows over bodies of arbitrary shape, and have had extensive applications in flow modeling [2][3][4][5][6]. The panel methods used for the simulation of flows over an airfoil are concerned with the vortex panel strength and circulation quantities; the evaluation of such quantities allows one to calculate the velocity distribution over the airfoil surface, and hence, to determine the pressure coefficients. The Kutta condition (a viscous boundary condition based on physical observation used with inviscid theoretical model) states that the flow leaves the sharp trailing edge of an airfoil smoothly [7]. Various methods have been proposed to impose the Kutta condition [8][9][10]. In panel methods, the Kutta condition is incorporated in the numerical formulation by requiring that the strength of vortex sheet is zero at the airfoil trailing edge. Moreover, the use of panel methods along with the compressibility corrections such as the Prandtl-Glauert method [11] allow one to consider the compressible flows over bodies. The panel methods have been extensively investigated in the aerodynamics literature, so these will not be discussed further here. However, these methods often have trouble with accuracy at the trailing edge of airfoils with zero angle cusped trailing edges [12]. Moreover, the compressibility corrections do not give accurate results for the compressible flows over airfoils of any shape at any angle of attack. For example, the Prandtl-Glauert method is based on the linearized perturbation velocity potential equation, and hence, it is limited to thin airfoils at small angles of attack [4].

Governing Equations
The stream function equation for two-dimensional, irrotational, steady, and isentropic flow of a compressible fluid in non-conservative form is as follows [14][15][16] (c 2 − u 2 )ψ xx + (c 2 − v 2 )ψ yy − 2uvψ xy = 0 (1) where ψ is the stream function, u and v are the components of the velocity vector V, i.e.,V = ui + vj (i and j are the unit vectors in x and y directions, respectively). For a two-dimensional compressible flow, c is the local sound speed [4] c 0 is the stagnation speed of sound, ρ is the density, ρ 0 is the stagnation density, and γ = c p /c v (c v and c p are the specific heats at constant volume and constant pressure, respectively) is the ratio of specific heats of gas (for air at standard conditions, γ = 1.4). Dividing both sides of Equation (1) by c 2 and then substituting the expressions in Equations (2) and (3) into Equation (1) gives the stream function equation as by substituting Equation (4) into Equation (5), we get and from the local isentropic stagnation properties of an ideal gas, we have c is the local Mach number. Equations (6) and (7) should be solved simultaneously to obtain the local values of ψ and ρ in the flow field.

Transformation
The solution of the governing PDE is based on transformation of the physical domain (x, y) and the governing equations into the regular computational domain (ξ, η) (see Figure 1). Therefore, the derivatives such as ψ x , ψ y , ψ xx , ψ yy , and ψ xy in the stream function equation, Equation (1), should be transformed from the physical domain (x, y) to the computational domain (ξ, η) [17][18][19]. This transformation can be stated as ξ ≡ ξ(x, y) η ≡ η(x, y) (9) and the inverse transformation is given as below.
x ≡ x(ξ, η) (10) y ≡ y(ξ, η) Since the stream function equation involves first and second derivatives, relationships are needed to transform such derivatives from the (x, y) system to the (ξ, η) one. In order to do this, the Jacobian of the transformation is needed, which is given below As will be shown, the transformation relations involve the Jacobian in denominator. Hence, it cannot be zero. Since we deal with the stream function equation, it is necessary to find relationships for the transformation of the first and second derivatives of the variable ψ with respect to the variables x and y. By using the chain rule, it can be concluded that By interchanging x and ξ, and y and η, the following relationships can also be derived By solving Equations (15) and (16) for ∂ψ ∂x and ∂ψ ∂y , we finally obtain where J = x ξ y η − x η y ξ is Jacobian of the transformation. By comparing Equations (13) and (17), and (14) and (18), it can be shown that To transform terms in the stream function equation (Equation (6)), the second order derivatives are needed. Therefore, which result in the following expressions for the second order derivatives where , , f x y ψ ≡ .

Boundary Conditions
Conditions at infinity: Far away from the airfoil surface (toward infinity), in all directions, the flow approaches the free stream conditions. The known free stream conditions are the velocity ∞ V , the pressure ∞ p , the density ∞ ρ , and the temperature ∞ T . Thus, the free stream Mach number is

Boundary Conditions
Conditions at infinity: Far away from the airfoil surface (toward infinity), in all directions, the flow approaches the free stream conditions. The known free stream conditions are the velocity V ∞ , the pressure p ∞ , the density ρ ∞ , and the temperature T ∞ . Thus, the free stream Mach number is where R is the specific gas constant, which is a different value for different gases. For air at standard conditions, R air = 287J/(kg.K), γ air = 1.4, ρ ∞ = 1.23kg/m 3 , and p ∞ = 1.01 × 10 5 N/m 2 . The critical Mach number M critical is that free stream Mach number at which sonic flow is first achieved on the airfoil surface. Condition on the airfoil surface: The relevant boundary condition at the airfoil surface for the inviscid flow is the no-penetration boundary condition. Thus, the velocity vector must be tangential to the surface. This wall boundary condition can be expressed by where s is tangent to the surface.

Grid Generation
Here, an O-grid is initially generated around the airfoil using elliptic grid generation method [17] and then the stream function equation is solved in the computational domain. The size of mesh is M × N where M is the number of nodes on the branch cut (a horizontal line connecting the trailing edge to outer boundary) and N is the number of nodes on the airfoil surface (and hence, the outer boundary), as shown in Figure 1. The physical domain before and after meshing using an O-grid of size 110 × 101 is shown in Figure 2 (using a NACA 0012 airfoil). Moreover, the magnified view of different parts of domain is shown in Figure 3 to highlight the orthogonality of the gridlines at airfoil surface and outer boundary. [17] and then the stream function equation is solved in the computational domain. The size of mesh is × M N where M is the number of nodes on the branch cut (a horizontal line connecting the trailing edge to outer boundary) and N is the number of nodes on the airfoil surface (and hence, the outer boundary), as shown in Figure 1. The physical domain before and after meshing using an O-grid of size × 110 101 is shown in Figure 2 (using a NACA 0012 airfoil). Moreover, the magnified view of different parts of domain is shown in Figure 3 to highlight the orthogonality of the gridlines at airfoil surface and outer boundary.

Kutta Condition
In inviscid flows, because the flow cannot penetrate the surface, the velocity vector must be tangential to the surface. In other words, the component of velocity normal to the surface must be zero and only the tangential velocity component must be considered. The unit tangent vector on the airfoil surface can be expressed as where ( ) ξ n and ( ) η n are the outward-pointing unit normal vector to a airfoil surface in ξ and η directions, respectively, and k is the unit vector in z direction. At the airfoil surface, corresponding to surface 3 S in Figure 4, we have

Kutta Condition
In inviscid flows, because the flow cannot penetrate the surface, the velocity vector must be tangential to the surface. In other words, the component of velocity normal to the surface must be zero and only the tangential velocity component must be considered. The unit tangent vector on the airfoil surface can be expressed as where n (ξ) and n (η) are the outward-pointing unit normal vector to a airfoil surface in ξ and η directions, respectively, and k is the unit vector in z direction. At the airfoil surface, corresponding to surface S 3 in Figure 4, we have The velocity component tangential to the airfoil surface ( 3 S ) is For compressible flows, the velocity components of u and v can be expressed in terms of stream function ψ as and where 0 ρ is a constant. On the airfoil surface, = =0 From the transformation relationships (Equation (19)), we have where α = x η 2 + y η 2 . Using Equation (35), we get τ (ξ) The velocity component tangential to the airfoil surface (S 3 ) is For compressible flows, the velocity components of u and v can be expressed in terms of stream function ψ as and where ρ 0 is a constant. On the airfoil surface, ψ s = ψ η = 0 where s is the distance measured along the airfoil surface because the airfoil contour is a streamline of the flow. Thus Equations (41) and (42) become By substituting Equations (43) and (44) into Equation (40), we get Based on the Kutta condition, the velocity at the trailing edge of the airfoil must be the same when approached from the upstream direction along the upper and lower airfoil surfaces [20] (see Figure 5). Thus, For the finite angle trailing edge, we have (The above expression is based on forward finite-difference coefficient with first-order accuracy. The forward finite-difference coefficient with second-order accuracy can be used if more accurate Kutta condition implementation is needed.) Therefore, the wall boundary condition for the finite angle trailing edge becomes because the stream function on the airfoil surface is constant. For the cusped trailing edge, we have The points (1, 1) and (1, N) as well as the points (2, 1) and (2, N) are the same points (see Figure 5). Therefore, ψ 1,1 = ψ 1,N and ψ 2, On the other hand, we know that Thus, the only possibility to satisfy Equation (49) is that which this relation again results in the same expression as for the finite angle trailing edge case. So, the general expression to implement the Kutta condition based on the stream function for the 2D compressible flow is which is the same as one for the incompressible flow [21] Fluids 2019, 4, x 10 of 27 which this relation again results in the same expression as for the finite angle trailing edge case. So, the general expression to implement the Kutta condition based on the stream function for the 2D compressible flow is which is the same as one for the incompressible flow [21] (a) (b)

Computation Procedure
According to the mapping scheme adopted in Figure 1, there are four sections where the nodal value of the flow variables , 1. Inside the domain to calculate the variables , 2. On the airfoil surface to calculate the variables 1, 3. At the outer boundary (far-field) to calculate the variables 4. On the branch cut to calculate the variables ,1 The known free-stream variables are from the local isentropic stagnation properties of an ideal gas, we have

Computation Procedure
According to the mapping scheme adopted in Figure 1, there are four sections where the nodal value of the flow variables f i,j ( f ≡ ψ, ρ, p, u, v, V, c, M, . . .) should be calculated.

1.
Inside the domain to calculate the variables On the airfoil surface to calculate the variables f 1, j ( j = 1, . . . , N).

4.
On the branch cut to calculate the variables f i,1 (i = 2, . . . , M − 1). We know that f i,N = f i,1 .
The known free-stream variables are M ∞ , p ∞ , ρ ∞ , T ∞ , and the angle of attack α. Thus, we can write from the local isentropic stagnation properties of an ideal gas, we have so the local total density ρ 0 M,j can be computed as In a similar fashion, we can write We know that if the general flow field is isentropic throughout, then the local total density ρ 0 , the local total pressure p 0 , and the local total temperature T 0 are constant values throughout the flow. Thus Also, on the outer boundary (far-field) we have If the number of nodes on the outer boundary (the side BG in Figure 1) is N, then, as shown in Figure 6, the node number of two points at top (point E) and bottom (point F) of the circular outer boundary would be (M, N+3 4 ) and (M, 3N+1 4 ), respectively.   Then we can express the magnitude of the stream function on the outer boundary ψ M,j in terms of far-field velocity V M,j as (by considering the equal magnitudes but opposite in sign for stream functions at top and bottom points of E and F, that is, Therefore, Now the magnitude of the stream function on the outer boundary can be calculated as a (based on the relation In Equation (69), we should note that the two points E and F have the same x-coordinates and The iterative process to obtain the value of ψ i,j and ρ i,j may be initiated by assuming This initial assumption implies that the magnitude of the speed of sound is infinite ( c = ∞ ⇒ 1 c 2 = 0 ) and Equation (5) reduces to Laplace's equation The first iteration is concerned with the solution of Laplace's equation which is described thoroughly in [21] and we do not elaborate further on it. After solving the Laplace's equation to find the value of the stream function at each node, ψ i,j , in the flow field, the values of the velocity components of u i,j and v i,j can be computed from Equations (2) and (3). Then, the value of the speed of sound at each node, c i,j ,can be calculated from Equation (4). Finally, the value of the density at each node, ρ i,j , can be determined from Equation (7). If the flow is compressible, ρ i, j ρ 0 i, j and from the second iteration onwards, instead of Laplace's equation, the stream function given in Equation (6) must be solved by having the density ρ i,j and the stream function ψ i,j from the last iteration (i.e., iteration 1) to get new values of the stream function ψ i,j . The iterative process (solving Equations (6) and (7) simultaneously) repeated until successive iterations produce a sufficiently small change in density ρ i,j and the stream function ψ i,j . Equation (6) is solved by initially discretizing the partial differential terms in Equations (21)-(26) using relations given in Equations (27)-(31) and then substituting the terms into Equation (6), and finally, solving the obtained algebraic equation using an algebraic solver (such as Maple) to get an explicit algebraic expression in terms of ψ i,j . As stated before, the Kutta condition is implemented as ψ 1,1 = ψ 2,1 (73) and the wall boundary condition is implemented as and since the branch cut (AB or HG in Figure 1) is inside the flow field, the same procedure employed to obtain an algebraic expression for ψ i,j can also be used to obtain an expression for the stream function on the branch cut, ψ i,1 . However, some changes are needed in the terms discretized by the finite-difference method (Equations (27)-(31)) as follows (see Figure 7) where f ≡ x, y, ψ.
In order to solve the elliptic grid generation equations (for x and y ) and the stream function equation, the iterative method Successive Over Relaxation (SOR) is used due to its high convergence rate where ≡ , , f x y ψ . In these relation for SOR, k is iteration number, and ω is relaxation factor which has a value of < < 1 2 ω . A value of = 1.8 ω is chosen in this study and the code FOILcom. We should define the stopping criteria for convergence of solution of the stream function equation (Equation (6) and density equation (Equation (7)). These two equations constitute the system of equations to be solved simultaneously. The stopping criteria are defined as follows In order to solve the elliptic grid generation equations (for x and y) and the stream function equation, the iterative method Successive Over Relaxation (SOR) is used due to its high convergence rate where f ≡ x, y, ψ. In these relation for SOR, k is iteration number, and ω is relaxation factor which has a value of 1 < ω < 2. A value of ω = 1.8 is chosen in this study and the code FOILcom.
We should define the stopping criteria for convergence of solution of the stream function equation (Equation (6) and density equation (Equation (7)). These two equations constitute the system of equations to be solved simultaneously. The stopping criteria are defined as follows where k is iteration number. A value of λ ψ = λ ρ = 10 −4 can be considered to get sufficiently accurate results. The other parameters of interest can also be computed after obtaining the density ρ i,j and the stream function ψ i,j . The components of the velocity at each node, u i,j and v i,j , can be computed from Then, the velocity at each node, V i,j , is given by and the local (nodal) Mach number can be obtained by where the local speed of sound,c i,j , is calculated from Equation (4). And finally, the local pressure, p i,j , is calculated from Now the values of drag and lift forces on the airfoil can be computed as follows where A and N are axial and normal forces, respectively ( Figure 8) [4]. We can write the drag and lift forces as and the drag and the lift coefficients may be expressed as Furthermore, the flowchart of the computational procedure is presented in Figure 9. The theoretical and computational details presented here for compressible flows are implemented in the freely available code FOILcom. Interested readers can refer to the code and use it by changing the airfoil shape, the free stream subcritical Mach number, and the angle of attack. The freely available code FOILincom(DOI: 10.13140/RG.2.2.21727.15524) also takes advantage of the same computational procedure but only for incompressible flow. Furthermore, the flowchart of the computational procedure is presented in Figure 9.  The theoretical and computational details presented here for compressible flows are implemented in the freely available code FOILcom. Interested readers can refer to the code and use it by changing the airfoil shape, the free stream subcritical Mach number, and the angle of attack. The freely available code FOILincom(DOI: 10.13140/RG.2.2.21727.15524) also takes advantage of the same computational procedure but only for incompressible flow.

Results
A few test cases are given to reveal the accuracy and robustness of the numerical scheme. Here, the results from FOILcom are compared with experimental and numerical results (alternative numerical schemes) to show its accuracy and efficiency.
Test case 1: Critical Mach number for the NACA 0012 airfoil at zero angle of attack. The airfoil surface pressure coefficient distribution using FOILcom is shown in Figure 10. The critical Mach number for the NACA 0012 airfoil at zero angle of attack is obtained as 0.722 using FOILcom with a grid of size 110 × 101 and x M,1 = 10 m. The critical Much number in references [4,22] is obtained as 0.725 (the experimental data in [22] is investigated in [4]). A comparison of results is shown in Figure 11 which reveals an excellent agreement.  [4,22] is obtained as 0.725 (the experimental data in [22] is investigated in [4]). A comparison of results is shown in Figure 11 which reveals an excellent agreement.   (Figures 12 and 13). The results from FOILcom ( Figure 12) and both finite volume and meshless methods (the results from finite volume and meshless methods in [23] are presented in one plot due to an excellent agreement between the results) are compared and depicted in Figure 13.   [4,22] is obtained as 0.725 (the experimental data in [22] is investigated in [4]). A comparison of results is shown in Figure 11 which reveals an excellent agreement.   (Figures 12 and 13). The results from FOILcom ( Figure 12) and both finite volume and meshless methods (the results from finite volume and meshless methods in [23] are presented in one plot due to an excellent agreement between the results) are compared and depicted in Figure 13. Refs. Figure 11. The pressure coefficient distribution. Comparison of FOILcom result and the one from References [4,22].
Test case 2: The airfoil surface pressure coefficient distribution for the NACA 0012 using a grid of size 110 × 101, M ∞ = 0.5, and α = 3 • (Figures 12 and 13). The results from FOILcom ( Figure 12) and both finite volume and meshless methods (the results from finite volume and meshless methods in [23] are presented in one plot due to an excellent agreement between the results) are compared and depicted in Figure 13.   (Figures 14 and 15). In this test case, the result from FOILcom is compared to the one from the code FLO42. There is excellent agreement between the results.    (Figures 14 and 15). In this test case, the result from FOILcom is compared to the one from the code FLO42. There is excellent agreement between the results. Test case 3: The airfoil surface pressure coefficient distribution for the NACA 0012 using a grid of size 110 × 101, M ∞ = 0.63, and α = 2 • (Figures 14 and 15). In this test case, the result from FOILcom is compared to the one from the code FLO42. There is excellent agreement between the results.      The results are shown in Figure 18. In this test case, four different grid sizes are used to obtain the airfoil surface pressure coefficient distribution. As shown in Figure 18, the distribution can be obtained with reasonable accuracy using even a course grid ( 30 41 × ). The effect of grid size on the drag and lift coefficients is depicted in Figure 19. Moreover, the drag and lift coefficients using the   The results are shown in Figure 18. In this test case, four different grid sizes are used to obtain the airfoil surface pressure coefficient distribution. As shown in Figure 18, the distribution can be obtained with reasonable accuracy using even a course grid ( 30 41 × ). The effect of grid size on the drag and lift coefficients is depicted in Figure 19. Moreover, the drag and lift coefficients using the The results are shown in Figure 18. In this test case, four different grid sizes are used to obtain the airfoil surface pressure coefficient distribution. As shown in Figure 18, the distribution can be obtained with reasonable accuracy using even a course grid (30 × 41). The effect of grid size on the drag and lift coefficients is depicted in Figure 19. Moreover, the drag and lift coefficients using the four different grid sizes are compared to the ones from XFOIL ( Figure 20) and are given in Table 1. As can be seen, the results from FOILcom are in excellent agreement with the result from XFOIL. four different grid sizes are compared to the ones from XFOIL ( Figure 20) and are given in Table 1.
As can be seen, the results from FOILcom are in excellent agreement with the result from XFOIL.  four different grid sizes are compared to the ones from XFOIL ( Figure 20) and are given in Table 1.
As can be seen, the results from FOILcom are in excellent agreement with the result from XFOIL.    Test case 6: The airfoil surface pressure coefficient distribution for NACA 64-012airfoil and M ∞ = 0.3, α = 6 • . The mesh size used in FOILcom is 50 × 141. The computation time is 46 s. The results from two codes FOILcom and FOILincom along with the Karman-Tsien compressibility correction are compared and depicted in Figure 21. Moreover, the convergence history of λ ψ and λ ρ are shown in Figure 22. As can be seen, in this test case five iterations are needed to satisfy the stopping criteria (λ ψ = λ ρ = 10 −4 ).  . In this test case, a thick airfoil, NACA 2240, is used to reveal the accuracy of the proposed numerical method in dealing with thick airfoils at high angles of attack. The results from FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL are compared and depicted in Figure 23 for the angle of attack 3 α =  and Figure 24 for the angle of attack 6 α =  . Furthermore, the results from XFOIL are given in Figure 25 for the angle of attack 3 α =  and Figure  26 for the angle of attack 6 α =  . The drag and lift coefficients using both codes for two different angles of attack are given in Table 2   . In this test case, a thick airfoil, NACA 2240, is used to reveal the accuracy of the proposed numerical method in dealing with thick airfoils at high angles of attack. The results from FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL are compared and depicted in Figure 23 for the angle of attack 3 α =  and Figure 24 for the angle of attack 6 α =  . Furthermore, the results from XFOIL are given in Figure 25 for the angle of attack 3 α =  and Figure  26 for the angle of attack 6 α =  . The drag and lift coefficients using both codes for two different angles of attack are given in Table 2 Figure 22. Convergence history of λ ψ and λ ρ .

Test case 7:
The airfoil surface pressure coefficient distribution for NACA 2240 airfoil and M ∞ = 0.3 at two different angles of attack α = 3 • and α = 6 • . The mesh size used in FOILcom is 120 × 141. In this test case, a thick airfoil, NACA 2240, is used to reveal the accuracy of the proposed numerical method in dealing with thick airfoils at high angles of attack. The results from FOILcom, FOILincom along with the Karman-Tsien compressibility correction, and XFOIL are compared and depicted in Figure 23 for the angle of attack α = 3 • and Figure 24 for the angle of attack α = 6 • . Furthermore, the results from XFOIL are given in Figure 25 for the angle of attack α = 3 • and Figure 26 for the angle of attack α = 6 • . The drag and lift coefficients using both codes for two different angles of attack are given in Table 2.

Conclusions
The solution of 2D steady, irrotational, subsonic (subcritical) compressible flow over isolated airfoils using the stream function equation and a novel method to implement the Kutta condition have been presented. The numerical scheme takes advantage of transformation of the flow solver and the boundary conditions from the physical domain to the computational domain. The physical domain was meshed by an O-grid elliptic grid generation method and the transformed flow solver is discretized by finite-difference method, a method chosen for its simplicity and ease of

Conclusions
The solution of 2D steady, irrotational, subsonic (subcritical) compressible flow over isolated airfoils using the stream function equation and a novel method to implement the Kutta condition have been presented. The numerical scheme takes advantage of transformation of the flow solver and the boundary conditions from the physical domain to the computational domain. The physical domain was meshed by an O-grid elliptic grid generation method and the transformed flow solver is discretized by finite-difference method, a method chosen for its simplicity and ease of

Conclusions
The solution of 2D steady, irrotational, subsonic (subcritical) compressible flow over isolated airfoils using the stream function equation and a novel method to implement the Kutta condition have been presented. The numerical scheme takes advantage of transformation of the flow solver and the boundary conditions from the physical domain to the computational domain. The physical domain was meshed by an O-grid elliptic grid generation method and the transformed flow solver is discretized by finite-difference method, a method chosen for its simplicity and ease of implementation. The numerical scheme is exempt from considering the panels and the quantities such as the vortex panel strength and circulation used in the panel method. An accurate Kutta condition scheme is proposed and implemented into the computational loop by an exact derived expression for the stream function at the airfoil trailing edge. The exact expression is general, and encompasses both the finite-angle and cusped trailing edges. Through several test cases, the proposed numerical scheme was validated by results from the experimental and the other numerical methods. The obtained results revealed that the proposed algorithm is very accurate and robust.