Next Article in Journal
Experimental Observations on Flow Characteristics around a Low-Aspect-Ratio Wall-Mounted Circular and Square Cylinder
Next Article in Special Issue
High Speed Flows
Previous Article in Journal
The Effects of Buoyancy on Laminar Heat Transfer Rates to Supercritical CO2 in Vertical Upward Flows
Previous Article in Special Issue
Experimental and CFD Investigation of Directional Stability of a Box-Wing Aircraft Concept
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Implementation of Flux Limiters in Simulation of External Aerodynamic Problems on Unstructured Meshes

1
Federal State-Funded Educational Institution of Higher Education, Nizhniy Novgorod State Technical University, 603999 Nizhniy Novgorod, Russia
2
Russian Federal Nuclear Center, All-Russia Research Institute of Experimental Physics, Federal State Unitary Enterprise, 607188 Sarov, Russia
3
Moscow Aviation Institute, 125993 Moscow, Russia
4
Federal State-Funded Higher Education Institution, Baltic State Technological University, 199406 St. Petersburg, Russia
*
Author to whom correspondence should be addressed.
Fluids 2023, 8(1), 31; https://doi.org/10.3390/fluids8010031
Submission received: 9 November 2022 / Revised: 20 December 2022 / Accepted: 6 January 2023 / Published: 15 January 2023
(This article belongs to the Special Issue High Speed Flows)

Abstract

:
The study is dedicated to the peculiarities of implementing the flux limiter of the flow quantity gradient when solving 3D aerodynamic problems using the system of Navier–Stokes equations on unstructured meshes. The paper describes discretisation of the system of Navier–Stokes equations on a finite-volume method and a mathematical model including Spalart–Allmaras turbulence model and the Advection Upstream Splitting Method (AUSM+) computational scheme for convective fluxes that use a second-order approximation scheme for reconstruction of the solution on a facet. A solution of problems with shock wave structures is considered, where, to prevent oscillations at discontinuous solutions, the order of accuracy is reduced due to the implementation of the limiter function of the gradient. In particular, the Venkatakrishnan limiter was chosen. The study analyses this limiter as it impacts the accuracy of the results and monotonicity of the solution. It is shown that, when the limiter is used in a classical formulation, when the operation threshold is based on the characteristic size of the cell of the mesh, it facilitates suppression of non-physical oscillations in the solution and the upgrade of its monotonicity. However, when computing on unstructured meshes, the Venkatakrishnan limiter in this setup can result in the occurrence of the areas of its accidental activation, and that influences the accuracy of the produced result. The Venkatakrishnan limiter is proposed for unstructured meshes, where the formulation of the operation threshold is proposed based on the gas dynamics parameters of the flow. The proposed option of the function is characterized by the absence of parasite regions of accidental activation and ensures its operation only in the region of high gradients. Monotonicity properties, as compared to the classical formulation, are preserved. Constants of operation thresholds are compared for both options using the example of numerical solution of the problem with shock wave processes on different meshes. Recommendations regarding optimum values of these quantities are provided. Problems with a supersonic flow in a channel with a wedge and transonic flow over NACA0012 airfoil were selected for the examination of the limiter functions applicability. The computation was carried out using unstructured meshes consisting of tetrahedrons, truncated hexahedrons, and polyhedrons. The region of accidental activation of the Venkatakrishnan limiter in a classical formulation, and the absence of such regions in case a modified option of the limiter function, is implemented. The analysis of the flow field around a NACA0012 indicates that the proposed improved implementation of the Venkatakrishnan limiter enables an increase in the accuracy of the solution.

1. Introduction

Unstructured meshes are most preferable currently in computations of aerodynamic flows for industrial applications [1,2]. Implementation of these meshes requires adaptation of the numerical methods being used, for example, in the part of numerical schemes construction, rated quantities approximation, and computation of flows and gradients.
One of specific features of the computations of aerodynamic flows is the possible availability of high-gradient areas in the flow, e.g., a shock wave. Robustness of the implemented algorithm depends on the stability of the discretisation scheme of convective fluxes when simulating such flows [3,4,5,6]. It is known that it is strictly recommended to implement schemes of an upgraded order of accuracy in the computations of the flows with shock waves and rarefaction waves, as the schemes of a lower order of accuracy tend to cause considerable “smearing” of the solution. However, they are more reliable in the way of fail-safe features than higher-order schemes.
It is possible to improve robustness of a high accuracy order scheme through the implementation of a scalar limiter of the gradient of the flow quantities. A gradient limiter is basically used in upwind schemes of the second order of accuracy. It allows preventing non-physical oscillations in the solution that are characteristic in computations of shock wave flows. A similar technique of limiter implementation is used in hydrodynamics [7]. Gradient limiters there prevent the value produced in gradient reconstruction of the quantity at the facet of the cell from violating the limits of its minimum and maximum in the cells neighboring the given one. The Venkatakrishnan limiter [8] is one of the most popular options in practice.
The experience in the implementation of limiters is long and initially appeared regarding structured meshes and the meshes with cells of a regular geometric shape [7,8]. It is evident, for example, from the summands as a part of the limiter expression, that they are related to the particular cell size of the computational mesh. In the case of unstructured meshes, it is difficult to find a characteristic size of a cell in the form of a random polyhedron. In this case, it appears logical and reasonable to set the operation threshold based solely on the flow quantities.
This work studies the implementation of the Venkatakrishnan limiter. Unstructured meshes consisting of tetrahedrons, truncated hexahedrons, and polyhedrons are used for the computation. A modified option of the Venkatakrishnan limiter is offered based on the results produced, and its advantages are shown in comparison to the initial option of the function when solving the problems with shock wave processes that occur at supersonic flow around the wedge and transonic flow around NACA0012 airfoil.
The paper is organized in three main sections. The section of basic equations describes the mathematical model used to simulate gas flow. The next section describes implementation of the flow limiter and solution of practical problems. The main conclusions are given in the final part of the paper.

2. Governing Equations

The physical and mathematical model to describe 3D flows is realized in the Russian software package called LOGOS. The LOGOS software package is designed for computational hydro- and aero-dynamics problems on parallel systems [9,10,11,12].
Non-stationary 3D turbulent flows of a viscous thermally conductive gas are described with Reynolds-averaged Navier–Stokes equations [13,14]. In a conservative form in Cartesian coordinates, the system of equations has the following form (averaging signs are skipped):
{ ρ t + ( ρ u ) = 0 , ( ρ u ) t + ( ρ u u ) = p + ( τ μ + τ t ) , ( ρ E ) t + ( ρ u h ) = [ u ( τ μ + τ t ) ( q μ + q t ) ] .
Here, ρ is density; u is the vector of the flow velocity with components u, v, w; p is pressure; E = C v T + 0.5 ( u 2 + v 2 + w 2 ) is total energy; h = C p T + 0.5 ( u 2 + v 2 + w 2 ) is total enthalpy; τ μ and τ t are molecular and turbulent components of the tensor of tangential stresses, respectively; q μ and q t are molecular and turbulent heat flux, respectively; T is temperature; C v = ( C p T R / m ) is specific heat capacity at constant volume; C p is specific heat capacity at constant pressure; R is a universal gas constant; and m is a molar mass of the gas.
The values of the molecular component of the tangential stress tensor of the Newtonian medium meet the rheological Newton's law, and the components of the density vector of the heat flow are connected with the local temperature gradient by Fourier’s law [13,14]:
τ μ = 2 μ ( T ) ( S 1 3 I u ) ,
S = 1 2 ( u + [ u ] t ) ,
q μ = λ ( T ) T .
The dynamic viscosity, μ ( T ) , and heat conductivity, λ ( T ) , are found from the Sutherland formula as a function of the flow temperature [14,15].
μ = μ 0 ( T T 0 ) 0.5 T 0 + T s T + T s ,
λ = λ 0 ( T T 0 ) 0.5 T 0 + T s T + T s ,
where μ 0 and λ 0 are dynamic viscosity and heat conductivity at temperature T 0 ; T s is the Sutherland constant.
System of Equation (1) is open due to the unknown connection of some of the basic variables of this system with averaged parameters of the flow. This connection can be established with additional ratios that, in a general case, are called turbulent models. Linear differential models of turbulence use empirical ratios for the turbulent viscosity factor. Here, the Spalart–Allmaras model [16,17] proved to be effective.
In the Spalart–Allmaras model, a single transport equation is considered. The transport equation is written with respect to the modified kinematic turbulent viscosity ν ˜ .
ρ ν ˜ t + ρ u j ν ˜ x j = 1 σ { x j [ ( μ + ρ ν ˜ ) ν ˜ x j ] + c b 2 ρ ( ν ˜ x j ) 2 } + P v D v
The generation and dissipation terms in Equation (7) are the source terms and they are expressed
D v = ( c w 1 f w c b 1 κ 2 f t 2 ) ( ν ˜ d ) 2
P v = c b 1 ρ S ˜ ν ˜ c b 1 ρ f t 2 S ˜ ν ˜
where d is the closet distance to the rigid wall, k is the von Karman constant.
The other parameters in the transport equation of turbulent viscosity can be found as follows:
S ˜ = Ω + f v 2 ν ˜ κ 2 d 2 .
Here, Ω is the rate of vorticity tensor
Ω = ( 2 Ω i j Ω i j ) 1 / 2
Ω i j = 1 2 ( u j x i u i x j )
f w = g ( 1 + C w 3 6 g 6 + C w 3 6 ) 1 / 6
g = r + C w 2 ( r 6 r )
r = ν ˜ S ˜ κ 2 d 2
C w 1 = c b 1 κ 2 + ( 1 + c b 1 ) σ
where the f t 2 function provides the suppression of the transition from the laminar flow calculations in the boundary layer to the turbulent flow calculations, and is expressed as
f t 2 = C t 3 exp ( C t 4 χ 2 ) .
Effective turbulent viscosity of the model is given next expression:
μ t = ρ ν ˜ f v 1
f v 1 = χ 3 χ 3 + C v 1 3 ,
χ = ν ˜ t ν .
Empirical constants of the model are as follows: σ = 2 3 , κ = 0.41 , c b 1 = 0.1355 , c b 2 = 0.622 , C w 2 = 0.3 , C w 3 = 2 , C v 1 = 7.1 , and C t 3 = 1.2 , C t 4 = 0.5 .

3. Numerical Method

The system above is approximated with the finite element method and it uses an integral formulation of the basic conservation laws. Discrete analogs of summands are written for the reference volume by summation over facets.
A finite volume method is based on integration of initial differential equations by the reference volume. Reference volumes (cells of the mesh) are arbitrary polyhedrons that cover the domain without gaps and overlaps. Each polyhedron is limited with a random number of facets. The vertices of the facets are the nodes of the mesh. A general view of the cell is given in Figure 1.
A system of the Navier–Stokes equation is written in a vector form for numerical solution using a finite volume method:
d d t Δ V W d V + Δ P ( F G ) d S = Δ V H   d V ,
where W is a vector of conservative variables, F and G are vectors of convective and diffusion fluxes, and H is a source term
W = ( ρ ρ u ρ v ρ w ρ E ) ,   F = ( ρ u n ρ u u n + p n x ρ v u n + p n y ρ w u n + p n z ρ H u n + p u n ) ,   G = ( 0 τ n x τ n y τ n z τ u + q ) ,
where un is a normal component of the velocity, q is a heat flux, and τij are components of the tensor of viscous stresses.
The full description of the way of approximation of the system of Equation (22) is given in [18,19].

4. Flux Limiters

When solving a problem, computation accuracy for the convective flows is of great importance. Schemes on the basis of the Riemann problem solution became quite popular in the case of aerodynamic problems. Such schemes comprise Advection Upstream Splitting Method (AUSM)–family schemes [20,21,22,23,24], and in particular AUSM+ scheme [21,22,23,24].
According to [25], a convective flow through the facet is computed in AUSM+ as follows:
F f = c f ( M ¯ L + U L + M ¯ R U R ) + ( P L + | α = 3 16 P L + P R | α = 3 16 P R )
where c f is the sound velocity at the facet; U L and U R are vectors of primitive variables on the left and on the right facet f ; P L and P R are vectors of pressure P = P { 0 , n x , n y , n z , 0 } T on the left and on the right facet f , M ¯ L + , M ¯ R , P L + | α = 3 16 , and P R | α = 3 16 are parameters of the scheme.
  • If M L + + M R 0 , where M L and M R are Mach numbers on the left and on the right of the facet, then
    M ¯ L + = M L + + M R [ ( 1 ω ) ( 1 + f R ) + f R f L ] .
  • If M L + + M R < 0 , then
    M ¯ R = M R + M L + [ ( 1 ω ) ( 1 + f L ) + f L f R ] .
Parameter ω is set by the function that depends on cubes of relations of pressure, and takes a minimum value in the larger part of the domain, except for the areas with a high gradient of pressure, such as the areas of shock waves and discontinuities
ω ( p L , p R ) = 1 min ( p L p R , p R p L ) 3 .
Parameter f L , R also takes a minimum value, except for the areas with oscillations of the solution
f L , R = { ( p L , R p s 1 ) 0 min ( 1 , min ( p 1 , L , p 1 , R , p 2 , L , p 2 , R ) min ( p L , p R ) ) 2 ,   P L + p L + P R p R 0 .
Second-order polynomials are used to find parameters on the facet:
M ± = { ± 1 4 ( M ± 1 ) 2 ,   | M | 1 , 1 2 ( M ± | M | ) ,   | M | > 1 , , P α ± = { 1 4 ( M ± 1 ) 2 ( 2 M ) ± α M ( M 2 1 ) , | M | 1 , 1 2 ( 1 ± sign ( M ) , | M | > 1 .
In order to compute convective fluxes, reconstruction of the solution is conducted. It lies in the definition of parameters on the left and on the right of facet f . When solving flow problems, reconstruction of the solution is performed with regard to primitive variables Q, conservative variables W, and with regard to acoustic invariants. For the first order of approximation, the values from the center of a respective cell (Figure 2) are taken as parameters on the right and on the left from the facet:
ϕ f = ϕ P , ϕ f + = ϕ E .
Reconstruction of the solution of the second-order approximation is usually taken to find the value on the facet [26,27].
ϕ f = ϕ P + α f ( Δ R P f ϕ P ) ϕ f + = ϕ E + α f + ( Δ R E f ϕ E )
where ϕ f and ϕ f + are values of the variable on the left and right facet, ϕ P and ϕ E are the values of the variable in the center of cells E and P (Figure 2), R P f and Δ R E f are the distance from the center of cells E and P to the center of the facet, ϕ E and ϕ P are the values of the gradient in cells E and P, and α f and α f + are limiters designed to prevent oscillations at discontinuous solutions. Implementation of the limiter function allows for controlling the gradient value on unstructured grids (decreasing the gradient value multiplying it by value α f 1 ); it is used to restore the value on the left and on the right from the facet [28,29].

5. Implementation of Flux Limiters

Implementation of so-called limiter functions is necessary to preserve monotonicity property of the solution in the areas with high gradients. In fact, introduction of the limiter of the scheme comes to “smoothening” of false maximum extremums of the values in the flow.
For example, in [28], as a result of limiter implementation, they managed to produce a smooth solution of trans-sonic flow without oscillations even on irregular meshes. The limiter function should be equal to zero in case of strong ruptures to produce the scheme of the first order that would guarantee monotonicity property, but in the areas of “monotonous” flow, the limiter function takes the value of a unit, and reconstruction of values at the facet is not limited. Transition from the limited value to the unlimited one should be smooth; only in this case would you expect upgraded convergence. Implementation of limiters is described in [29,30,31,32,33,34,35].
Correct behaviour of the limiter function is especially important when it is used in engineering codes to solve industrially specific problems on unstructured meshes.
More than a dozen different limiter functions were made available and published; they were reviewed in [36], for example. The most often used ones are:
  • Barth–Jespersen limiter [37]: its general view is α f = min ( r f , 1 ) ;
  • Van Albada–Leer [7]: its general view is α f = min ( r f 2 + r f r f 2 + 1 , 1 ) ;
  • Venkatakrishnan [38]: its general view is α f = min ( r f 2 + 2 r f r f 2 + r f + 2 , 1 ) .
Figure 3 shows a Sweby diagram [29] that demonstrates dependence of α f values on factor r f .
The α f = 0 values correspond to the scheme of the first order of accuracy, and α f = 1 to the scheme of the second order of accuracy. You see that that “strict” limiter function is the Venkatakrishnan limiter that makes transition to the scheme of the second order the latest of all (at r f = 2 ), which will provide high stability and ensure good monotonous property [7,8,38].
The Venkatakrishnan limiter controls the value of gradient ϕ E in cell E according to the expression:
α E = { 1 Δ 2 [ ( Δ l , max 2 + ε 2 ) Δ 2 + 2 Δ 2 2 Δ l , max Δ l , max 2 + 2 Δ 2 2 + Δ l , max Δ 2 + ε 2 ] , Δ 2 > 0 1 Δ 2 [ ( Δ l , min 2 + ε 2 ) Δ 2 + 2 Δ 2 2 Δ l , min Δ l , min 2 + 2 Δ 2 2 + Δ l , min Δ 2 + ε 2 ] , Δ 2 < 0 1 , Δ 2 = 0
Δ l , max = ϕ max ϕ E ,   Δ l , min = ϕ min ϕ E , Δ 2 = 1 2 ( ϕ f Δ R E f ) .
ε 2 = ( K Δ h ) 3
where ϕ max and ϕ min are maximum and minimum values in all neighboring cells, including the values in cell E itself, and Δ R E f is the distance from the center of cell E (or P) to the center of the facet. Parameter ε2 controls the value of the limiter, where K is a constant (a normalizing coefficient), Δ h is a characteristic size of the cell.
In (31), ε2 is a symbolic operation threshold of the limiting function. Oscillations lower than this threshold are allowed in the solution and are not considered by the limiter. A zero value of parameter ε2 means that the limiter is active even in the near-constant regions, when a very high value of parameter ε2 means practically no limit. Such modification of the limiting function makes it possible to protect from random operations and reach improvement of convergence and stable solution on unstructured meshes.
Let us mark
y = Δ l , max Δ 2 o r ( y = Δ l , min Δ 2 )
and write the function from expression (31) as follows [8,37]:
α E = y 2 + 2 y + ε y 2 + y + 2 + ε
where Δ l , max or Δ l , min from (32) is increment of the solution.
With quantity ε exceeding the increment value of the solution, or with the unlimited quantity ε, the limiter takes the value of a unity, i.e., the value of the gradient (in the expression to find the value per facet) is not limited to anything. Where the increment of the solution exceeds quantity ε (e.g., in the region with large gradients or at a small value of ε), the solution itself determines the value of the limiting function (summand y 2 + 2 y y 2 + y + 2 of expression (34)); in this way, it sets the degree of limitation.

6. Results and Discussion

Several CFD benchmark problems are considered to validate the robustness of the flux limiters and their parameters.

6.1. Flow in a Channel with Wedge

The operation of the limiter (at various values of parameter K) is considered using the example of the problem of supersonic flow around the wedge, where the parameters of the incident flow are as follows: Mach number is 2, the pressure is 101,325 Pa, and the temperature is 300 K [39]. Structured computational mesh is used for simulation with the number of cells at 95,000. General view of the computation domain and mesh are shown in Figure 4.
At supersonic flow in the channel with a wedge, an attached compression shock is formed, which results in the formation of the shock wave structure of the flow in the channel. Figure 5 shows the origination of the compression shock, its development and reflection from the walls of the channel, and its interaction with the system (finitary spread) of rarefaction waves [39].
Each shock wave is characterized with its front (the surface where flow quantities have a jump in the development, whereas outside the front they change continuously). From the practical point of view, the implemented numerical scheme should provide stability and preserve monotonicity of the solution in all regions of the flow, including the shock wave fronts.
Let us study the section of the first series of compression shocks to explore the solution as a function of the limiter in Figure 5. Look at the plot in Figure 6, of the distribution of the full pressure along the line in Figure 4 in the specified region produced with and without the limiter (Venkatakrishnan limiter) for different values of parameter K in expression (33).
According to the plot, at K = 0.1 and K = 1 the produced result complies with the computation result without a limiter, i.e., at such values of K, the solution is fully unlimited (complete computed value of the gradient is taken at the reconstruction step). The solution produced at K = 0.01 is characterized with local maximums that have a non-physical (fictitious) origination character that can cause a considerable error in the solution, as it acquires extremely non-monotonic behaviour.
The solution has the best accuracy (as compared to the analytical solution) and monotonicity properties at K = 0 and K = 0.001. However, K = 0 means probability of switching on the limiter in the entire domain, i.e., it has a random reaction character in Figure 7, and this is also intolerable.
The solution at K = 0.001 has a minimum amplitude of oscillations and, in general, is characterized with a good monotonicity property. However, due to a very small value of parameter K, the limiter is practically always activated (random cases of operation on other mesh models can be revealed), and the picture of the produced solution actually reflects the property of this limiter that depends on the size of the mesh cell.
Computation of parameter ε by formula (33) is related to the characteristic cell size only, and the function depends only on the geometric parameters of the mesh. In this case, the limiter operation depends only on the parameters of the computation model and increases the probability of the limiter’s reaction in the regions of the local refinement mesh model. At the same time, in this formulation parameter, ε is not related in any way to the flow quantities. This also increases the probability of random operation of the limiter function for a particular flow quantity. The characteristic cell size is valid only in the case of cells of a regular geometric shape; whereas, in the case of a cell with a shape of a random polynomial, this value does not have a clear definition.
Modification of expression (33) to compute ε is necessary to use the limiter on cells of a random shape. The idea is to make parameter ε the function of the flow quantities, i.e., the function of that quantity, for which the limiter is implemented. Another option to define ε is
ε ˜ = K ϕ
where K is an operation threshold of function and ϕ is a flow variable. The value of the limiter is computed on the basis of the flow variable itself, for which the limiter is computed (when pressure is computed, ϕ shows pressure, and when density is computed, it shows density; same for velocity).
Introduction of a flow variable as a parameter when computing ε ˜ disconnects the limiter from geometric parameters of the mesh cells. In the current formulation parameter, ε ˜ has a physical sense. Changing constant K, the value of oscillations of computational flow quantities filtered with the limiter is found. For example, K = 0.01 means that the operation threshold of the limiter is equal to 1% from flow quantity ϕ , i.e., the limiter is switched on when the solution becomes oscillating and the increment of the solution in the cell is higher than 1% of the solution in the cell.
Let us analyse the implementation of the modified formula for the problem of a supersonic flow in the channel with a wedge in Figure 8, and look at the case when K = 0.01 in (35).
Distribution of total pressure in case K = 0.01 from (35) has a comparable solution with the one produced in case K = 0.001 in (33), the property of monotonicity. However, if expression (35) is implemented, the limiter function has a clear physical interpretation. In this case, the region of activation of the limiter is characterized by the absence of random reactions in Figure 9.
To study the implementation on unstructured meshes, 3D meshes (one-cell thick) were generated for this geometry consisting of polynomials, tetrahedrons, and truncated hexahedrons. The behaviour of the limiters is shown in Figure 10. Note that, when K = 0.001 in (33), multiple areas of random reactions are observed in all options of mesh models.
So, having analysed quantity ε as it influences the behaviour of the solution, it can be concluded that implementation of the Venkatakrishnan limiter, where ε is computed by expression (33) at K = 0.001, contributes to the monotonicity of the convergence process of the solution and reduces the number of non-physical oscillations. However, ε depends on geometric sizes of the cell of the computation mesh. This could be the reason for the possible random character of the reaction of the function under consideration (revealed on unstructured meshes) and introduces a numerical error into the solution. To eliminate this problem, option (35) was suggested at K = 0.01, which has comparable monotonicity properties. However, it is characterized with no random character of the reaction of the limiter on random unstructured meshes. Form (35) provides dependence of quantity ε on the intensity of the flow, allowing for more precise definition of the activation threshold of the limiter.

6.2. Flow around Airfoil

Let us consider the limiter as it influences the accuracy of the produced solution on the example of transonic flow around NACA0012 airfoil [40,41]. The mesh is one cell thick in Figure 11. The total number of mesh elements is 731,000.
According to the research on the components of the drag of aerodynamic airfoils [13], the basic constituent of the drag for a medium-thick airfoil at small attack angles is friction resistance. The airfoil drag (pressure drag) due to the incomplete restoration of pressure in the tail part of the airfoil is 20–30% of the total profile drag (at small attack angles of 0–3° for symmetric and slightly bended airfoils with an average relative thickness of 11–15%).
As [13,39] show, the accuracy of the pressure drag computation depends directly on the level of so-called numerical viscosity (approximation errors that work as additional dissipation, resulting in the loss of complete pressure in the flow and growth of the resistance of the object under investigation). It can be caused by some peculiarities of the numerical method. For example, the selection of the limiter of the gradient increment could reduce the accuracy order greatly in the regions with large gradients of gas parameters. As a result of the incorrect operation of the limiter, not enough rarefaction can be observed in the region of minimum values of the pressure coefficient on the upper surface of the airfoil.
Let us consider this problem with the following boundary conditions. Parameters of the incident flow at the external boundary of the computational domain correspond to the values: the pressure is 46,066.2 Pa, the temperature is 248 K, the Mach number is 0.7, and the angle of attack is 1.489° and 3.046° [39,40,41]. The surface of the airfoil is considered to be a no-penetrated wall; a symmetric boundary condition is set at the sides.
At the set parameters of the airfoil overflow, the formation of the tear-off zone takes place on its leeside. A curvilinear shock wave is formed near the surface of the airfoil. It corresponds to the normal that transforms a supersonic flow into the subsonic flow. In case the angle of attack is 1.489°, a low-intensity compression shock is formed with smeared boundaries of the transfer of the supersonic flow into the subsonic flow in Figure 12.
When increasing the angle of attack to 3.046°, a more intensive compression shock is formed above the airfoil. It has a more expressed front and creates a large area of decreased pressure. Mach number and pressure contours are presented in Figure 13.
No differences in the solution with different options of the limiter are observed in the general structure of the flow (in the shock wave generation above the airfoil). Cognominal fields correspond to each other. However, the differences in the solution are evident if integral characteristics are estimated, for example, the value of the drag force coefficient, C d r a g , as compared to the experimental data in Table 1, [39,40,41].
Maximum error in the solution for both angles of attack is produced in the case of computation No 2. The main contribution into the error in the resistance computation is due to the pressure force component that happens to be overestimated by more than 10%. One of the most probable reasons for these results could be the work of the gradients increment limiter, which is of a random reaction character in all the computation domains revealed in this mesh model in Figure 14.
As for computation No 3, the limiter in this case was active only in the region of the shock wave above the airfoil in Figure 14, which is completely true for form ἐ (35) dependent on the flow quantities.
As a result, a solution was produced that has a minimum deviation from the experimental data (2.3% and 9.1%). It is worth noting that the error in the computation grows with the increase in the angle of attack, and it is approaching the critical values. It is connected to the formation of a more intensive shock wave above the surface of the airfoil and a more complicated operation of the limiter in this region. In particular, it reduces the accuracy order of the scheme most precisely in the regions with large gradients of flow quantities and not to introduce additional numerical viscosity. So, evident advantages of implementation of the limiter in combination with parameter ἐ at K = 0.01 are shown using this problem.

6.3. Flow around Bullet

A modified option of the limiter can also be applied to compute supersonic flow around the “168 Grain Sierra International Bullet” [42]. An unstructured mesh of truncated hexahedrons with a general number of elements of 386 thousand was generated for the computation of this problem in Figure 15.
Regions of the flow, where the Mach number is close to 1, are characterized with complexity due to the formation of compression shocks and possible separation of the flow that influence aerodynamic properties of the object considerably. The results produced here can also be generalized for the solids geometrically similar to the object under consideration.
The geometry of the bullet, the properties of the flow, and experimental aerodynamic parameters are taken from the work of Author [42]. Let us take the problem with the following boundary conditions. Parameters of the striking flow at the external boundary of the domain correspond to the following quantities: the pressure is 101,325 Pa, the temperature is 288.15 K, the angle of attack is 0°, and the Mach number is 1.05 and 1.6 [42]. The surface of the bullet is taken as a solid wall, and a symmetric boundary condition is set on the side boundary.
In case Mach number is equal to 1.05 and 1.6, the flow is characterized by the presence of the head shock wave in Figure 16. Distribution of flow quantities near the surface of the object influences aerodynamic properties considerably. That is why it is necessary to provide correct operation of the limiter in this region to get a high-quality solution. There are no differences in the solutions with different variants of the limiter function in the general structure of the flow in Figure 16 and Figure 17.
Let us estimate the value of the drag force coefficient, C d r a g , in comparison with experimental data in Table 2, [42]. The largest deviation of the results from the experimental data are observed in the computation with limiter ε (33) at K = 0.001 for both modes. So, the modified form of the limiter allows for producing a more accurate solution in case the flow with shock wave processes is considered on an unstructured mesh.

7. Conclusions

The study analyses implementation of the system of Navier–Stokes equations to simulate the problems of aerodynamics. A mathematical model is supplemented with equations of the Spalart–Allmaras turbulence model and AUSM+ scheme to compute convective fluxes that use a second-order approximation scheme for reconstruction of the solution on the facet. The Venkatakrishnan limiter is implemented to prevent the generation of false oscillations of the solution when computing shock wave processes.
The work describes the research on the influence of this limiter on the behaviour of the numerical solution of aerodynamic problems. It shows that implementation of the gradient limiter helps to improve monotonicity of the convergence when simulating the problems with shock waves and local compression shocks. Venkatakrishnan limiter implementation was studied on structured and unstructured meshes. It was found out that the initial option of the Venkatakrishnan limiter on unstructured meshes can produce the regions of its random activation. To eliminate this phenomenon, it was suggested to modify the summand that controls the limiter operation threshold. The suggested option is based on the flow quantities, and it showed correct behaviour on unstructured meshes that helps to improve the accuracy of the produced solution.
The proposed option of the limiter allows for producing a more accurate solution for the problems of trans- and supersonic flows, and that was demonstrated on the example of overflowing an airfoil and a bullet.
Computations on hypersonic flows are planned, as well the study of their peculiarities.

Author Contributions

Conceptualization, A.S.K. and D.Y.S.; data curation, A.V.S. and R.N.Z.; formal analysis, A.S.K., R.N.Z. and K.N.V.; investigation, A.S.K., A.V.S., D.Y.S. and R.N.Z.; methodology, A.S.K. and R.N.Z.; software, A.V.S. and R.N.Z.; supervision, A.S.K.; validation, A.V.S.; visualization, A.V.S. and K.N.V.; writing—original draft, A.S.K. and D.Y.S. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by the program for the creation and development of the World-Class Research Center “Supersonic” for 2020–2025 funded by the Ministry of Science and Higher Education of the Russian Federation (Grant agreement of 20 April 2022 № 075-15-2022-309).

Acknowledgments

The results have been obtained with financial support from the Science & Universities National Project under the Young Scientists Lab Program of the RF Ministry of Education and Science # FSWE-2021-0009 (Research Topic: Development of CFD methods, models and algorithms to simulate liquids and gases in natural and industrial environments under normal and critical conditions on petascale supercomputers) and the Council of the Grants of the President of the Russian Federation for state support of Leading Scientific Schools of the Russian Federation (Grant No. NSH-70.2022.1.5).

Conflicts of Interest

The authors declare no conflict of interest. All authors have read and agreed to the published version of the manuscript.

References

  1. Kozelkov, A.S.; Krutyakova, O.L.; Kurulin, V.V.; Lashkin, S.V.; Tyatyushkina, E.S. Application of numerical schemes with singling out the boundary layer for the computation of turbulent flows using eddy-resolving approaches on unstructured grids. Comput. Math. Math. Phys. 2017, 57, 1036–1047. [Google Scholar] [CrossRef]
  2. Kozelkov, A.S.; Lashkin, S.V.; Efremov, V.R.; Volkov, K.N.; Tsibereva, Y.A.; Tarasova, N.V. An implicit algorithm of solving Navier–Stokes equations to simulate flows in anisotropic porous media. Comput. Fluids 2018, 160, 164–174. [Google Scholar] [CrossRef] [Green Version]
  3. Kumar, V.; Sharma, A.; Singh, R.K. Central upwind scheme based immersed boundary method for compressible flows around complex geometries. Comput. Fluids 2020, 196, 104349. [Google Scholar] [CrossRef]
  4. Pirozzoli, S. Numerical methods for high-speed flows. Annu. Rev. Fluid Mech. 2011, 43, 163–194. [Google Scholar] [CrossRef]
  5. Du, X.; Corre, C.; Lerat, A. A third-order finite-volume residual-based scheme for the 2D Euler equations on unstructured grids. Comput. Phys. 2011, 230, 4201–4215. [Google Scholar] [CrossRef]
  6. Zhu, W.; Xiao, Z.; Fu, S. Studies of flight-velocity effects on near-field and intermittent properties of a subsonic jet. Comput. Fluids 2020, 196, 104351. [Google Scholar] [CrossRef]
  7. Van Albada, G.D.; Van Leer, B.; Roberts, W.W., Jr. A comparative study of computational methods in cosmic gas dynamics. Astron. Astrophys. 1982, 108, 76–84. [Google Scholar]
  8. Venkatakrishnan, V. Convergence to steady state solutions of the Euler equations on unstructured grids with limiters. J. Comput. Phys. 1995, 118, 120–130. [Google Scholar] [CrossRef]
  9. Kozelkov, A.S.; Kurulin, V.V.; Lashkin, S.V.; Shagaliev, R.M.; Yalozo, A.V. Investigation of Supercomputer Capabilities for the Scalable Numerical Simulation of Computational Fluid Dynamics Problems in Industrial Applications. Comput. Math. Math. Phys. 2016, 56, 1506–1516. [Google Scholar] [CrossRef]
  10. Kozelkov, A.S.; Pogosyan, M.A.; Strelets, D.Y.; Tarasova, N.V. Application of mathematical modeling to solve the emergency water landing task in the interests of passenger aircraft certification. AS 2021, 4, 75–89. [Google Scholar] [CrossRef]
  11. Kozelkov, A.S. The Numerical Technique for the Landslide Tsunami Simulations Based on Navier-Stokes Equations. J. Appl. Mech. Tech. Phys. 2017, 58, 1192–1210. [Google Scholar] [CrossRef]
  12. Tyatyushkina, E.S.; Kozelkov, A.S.; Kurkin, A.A.; Pelinovsky, E.N.; Kurulin, V.V.; Plygunova, K.S.; Utkin, D.A. Verification of the LOGOS Software Package for Tsunami Simulations. Geosciences 2020, 10, 385. [Google Scholar] [CrossRef]
  13. Landau, L.D.; Lifshits, V.M. Theoretical physics. In Hydrodynamics; Nauka: Moscow, Russia, 1988; Volume IV, 736p. [Google Scholar]
  14. Loytsyanskiy, L.G. Mechanics of Fluids and Gases; Nauka: Moscow, Russia, 1979; 904p. [Google Scholar]
  15. Shlikhting, G. Theory of a Boundary Layer; Nauka: Moscow, Russia, 1974; 712p. [Google Scholar]
  16. Spalart, P.R.; Allmaras, S.R. A One-Equation Turbulence Model for Aerodynamic Flows; No. 0439; AIAA Paper: Reno, NV, USA, 1992. [Google Scholar]
  17. Shur, M.; Strelets, M.; Travin, A.; Spalart, P.R. Turbulence modeling in rotating and curved channels: Assessment of the Spalart—Shur correction term. AIAA J. 2000, 38, 784–792. [Google Scholar] [CrossRef]
  18. Deryugin, Y.N.; Zhuchkov, R.N.; Zelenskiy, D.K.; Kozelkov, A.S.; Sarazov, A.V.; Kudimov, N.F.; Lipnickiy, Y.M.; Panasenko, A.V.; Safronov, A.V. Validation Results for the LOGOS Multifunction Software Package in Solving Problems of Aerodynamics and Gas Dynamics for the Lift-Off and Injection of Launch Vehicles. Math. Model. Comput. Simul. 2015, 7, 144–153. [Google Scholar] [CrossRef]
  19. Kozelkov, A.S.; Struchkov, A.V.; Strelets, D.Y. Two Methods to Improve the Efficiency of Supersonic Flow Simulation on Unstructured Grids. Fluids 2022, 7, 136. [Google Scholar] [CrossRef]
  20. Roe, P.L. Characteristic Based Schemes for the Euler Equations. Annu. Rev. Fluid Mech. 1986, 18, 337–365. [Google Scholar] [CrossRef]
  21. Liou, M.-S. A Sequel to AUSM: AUSM+. J. Comput. Phys. 1996, 129, 364–382. [Google Scholar] [CrossRef]
  22. Rodionov, A.V. Artificial viscosity to cure the carbuncle phenomenon: The three dimensional case. J. Comput. Phys. 2018, 361, 50–55. [Google Scholar] [CrossRef]
  23. Kotov, D.V.; Surzhikov, S.T. Calculation of Viscous and Inviscid Gas Flows on Unstructured Grids Using the AUSM Scheme. Fluid Dyn. 2011, 46, 809–825. [Google Scholar] [CrossRef]
  24. Shingo, M. Performance of all-speed AUSM-family schemes for DNS of low Mach number turbulent channel flow. Comput. Fluids 2014, 91, 130–143. [Google Scholar]
  25. Jan, V.; Bart, M.; Erik, D. Blended AUSM+ Method for All Speeds and All Grid Aspect Ratios. AIAA J. 2001, 39, 2278–2282. [Google Scholar]
  26. Kim, K.H.; Ch, K.; Rho, O.-H. Methods for the accurate computations of hypersonic flows. I AUSMPW+ scheme. J. Comput. Phys. 2001, 174, 38–80. [Google Scholar] [CrossRef]
  27. Ferziger, J.H.; Peric, M. Computational Method for Fluid Dynamics; Springer: New York, NY, USA, 2002. [Google Scholar]
  28. Jasak, H. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flow. Ph.D. Thesis, Department of Mechanical Engineering, Imperial College of Science, London, UK, 1996. [Google Scholar]
  29. Sweby, P.K. High resolution using flux limiters for hyperbolic conservation laws. SIAM J. Numer. Anal. 1984, 21, 995–1011. [Google Scholar] [CrossRef]
  30. Darwish, M.S.; Moukalled, F. TVD schemes for unstructured grids. Int. J. Heat Mass Transf. 2003, 46, 599–611. [Google Scholar] [CrossRef]
  31. Balcázar, N.; Jofre, L.; Lehmkuhl, O.; Castro, J.; Rigola, J. A finite-volume/level-set method for simulating two-phase flows on unstructured grids. Int. J. Multiph. Flow 2014, 64, 55–72. [Google Scholar] [CrossRef]
  32. Li, L.-X.; Liao, H.-S.; Qi, L.-J. An improved r-factor algorithm for TVD schemes. Int. J. Heat Mass Transf. 2008, 51, 610–617. [Google Scholar] [CrossRef]
  33. Hou, J.; Simons, F.; Hinkelmann, R. Improved total variation diminishing schemes for advection simulation on arbitrary grids. Int. J. Numer. Meth. Fluids 2012, 70, 359–382. [Google Scholar] [CrossRef]
  34. Denner, F.; van Wachem, B. TVD differencing on three-dimensional unstructured meshes with monotonicity-preserving correction of mesh skewness. J. Comput. Phys. 2015, 298, 466–479. [Google Scholar] [CrossRef] [Green Version]
  35. Balcázar, N.; Antepara, O.; Rigola, J.; Oliva, A. A level-set model for mass transfer in bubbly flows. Int. J. Heat Mass Transf. 2019, 138, 335–356. [Google Scholar] [CrossRef]
  36. Gradient Computation. Available online: http://www.cfd-online.com/Wiki/Gradient_computation (accessed on 25 August 2019).
  37. Barth, T.J.; Jespersen, T.J. The Design and Application of Upwind Schemes on Unstructured Grids; No. 89-0366; AIAA Paper: Reno, NV, USA, 1989. [Google Scholar]
  38. Venkatakrishnan, V. On the Accuracy of Limiters and Convergence to Steady State Solution; No. 93-0880; AIAA Paper: Reno, NV, USA, 1993. [Google Scholar]
  39. Povkh, I.L. Technical Hydromechanics; Mashinostroenie: Leningrad, Russia, 1976; 504p. [Google Scholar]
  40. Kolesnikov, G.A.; Markov, V.K.; Mikhaykyuk, A.A. Aerodynamics of Flying Vehicles; Mashinostroeniye: Moscow, Russia, 1993; 544p. [Google Scholar]
  41. Pilipenko, A.A.; Polevoy, O.B.; Prihodko, A.A. Numerical simulation of Mach number and attack angle as they influence the modes of trans-sonic turbulent overflowing of aerodynamic airfoils. Tech. Notes TsAGI 2012, XLIII, 1–110. [Google Scholar]
  42. Chang, P. Separated Flow; Mir: Moscow, Russia, 1973; Volume 3, 335p. [Google Scholar]
Figure 1. Control volume.
Figure 1. Control volume.
Fluids 08 00031 g001
Figure 2. Reconstruction of values.
Figure 2. Reconstruction of values.
Fluids 08 00031 g002
Figure 3. Sweby diagram (–Barth-Jespersen, ▬ ▬–Van Albada–Leer, ▬▬–Venkatakrishnan).
Figure 3. Sweby diagram (–Barth-Jespersen, ▬ ▬–Van Albada–Leer, ▬▬–Venkatakrishnan).
Fluids 08 00031 g003
Figure 4. General view of the computation domain and mesh.
Figure 4. General view of the computation domain and mesh.
Fluids 08 00031 g004
Figure 5. Contours of total pressure.
Figure 5. Contours of total pressure.
Fluids 08 00031 g005
Figure 6. The distribution of full pressure.
Figure 6. The distribution of full pressure.
Fluids 08 00031 g006
Figure 7. Random behaviour of the limiter (31) at K = 0.
Figure 7. Random behaviour of the limiter (31) at K = 0.
Fluids 08 00031 g007
Figure 8. The distribution of total pressure.
Figure 8. The distribution of total pressure.
Fluids 08 00031 g008
Figure 9. The region of activation of the limiter when K = 0.01 in (35).
Figure 9. The region of activation of the limiter when K = 0.01 in (35).
Fluids 08 00031 g009
Figure 10. Activation area of the limiter. Fragments (a,c,e) correspond to the limiter (31) when K = 0.001 in (33). Fragments (b,d,f) correspond to the limiter (31) when K = 0.01 in (35).
Figure 10. Activation area of the limiter. Fragments (a,c,e) correspond to the limiter (31) when K = 0.001 in (33). Fragments (b,d,f) correspond to the limiter (31) when K = 0.01 in (35).
Fluids 08 00031 g010
Figure 11. General view of the computation domain and mesh.
Figure 11. General view of the computation domain and mesh.
Fluids 08 00031 g011
Figure 12. Contours of Mach number (ac) and contours of pressure (df): computation without a limiter (a,d), a computation with a limiter function at K = 0.001 for (33) (b,e), and computation with a limiter at K = 0.01 for (35) (c,f).
Figure 12. Contours of Mach number (ac) and contours of pressure (df): computation without a limiter (a,d), a computation with a limiter function at K = 0.001 for (33) (b,e), and computation with a limiter at K = 0.01 for (35) (c,f).
Fluids 08 00031 g012
Figure 13. Contours of Mach number (ac) and contours of pressure (df): computation without a limiter (a,d), a computation with a limiter function at K = 0.001 for (33) (b,e), and computation with a limiter at K = 0.01 for (35) (c,f).
Figure 13. Contours of Mach number (ac) and contours of pressure (df): computation without a limiter (a,d), a computation with a limiter function at K = 0.001 for (33) (b,e), and computation with a limiter at K = 0.01 for (35) (c,f).
Fluids 08 00031 g013
Figure 14. The region of the limiter activation (ε (33) at K = 0.001—on the left, ἐ (35) at K = 0.01—on the right).
Figure 14. The region of the limiter activation (ε (33) at K = 0.001—on the left, ἐ (35) at K = 0.01—on the right).
Fluids 08 00031 g014
Figure 15. General view of the computation domain and mesh.
Figure 15. General view of the computation domain and mesh.
Fluids 08 00031 g015
Figure 16. Contours of Mach number (a,b) and contours of pressure (c,d): computation with a limiter at K = 0.001 for (33) (a,c), and computation with a limiter at K = 0.01 for (35) (b,d).
Figure 16. Contours of Mach number (a,b) and contours of pressure (c,d): computation with a limiter at K = 0.001 for (33) (a,c), and computation with a limiter at K = 0.01 for (35) (b,d).
Fluids 08 00031 g016
Figure 17. Contours of Mach number (a,b) and contours of pressure (c,d): computation with a limiter at K = 0.001 for (33) (a,c), and computation with a limiter at K = 0.01 for (35) (b,d).
Figure 17. Contours of Mach number (a,b) and contours of pressure (c,d): computation with a limiter at K = 0.001 for (33) (a,c), and computation with a limiter at K = 0.01 for (35) (b,d).
Fluids 08 00031 g017
Table 1. Drag coefficients.
Table 1. Drag coefficients.
NoComputationα = 1.489°α = 3.046°
C d r a g Δ C d r a g ,   % C d r a g Δ C d r a g ,   %
0Experiment0.008190.01267
1Without limiter0.008483.60.0142312.3
2ε (33) at K = 0.0010.008503.70.0142512.5
3ἐ (35) at K = 0.010.008382.30.013829.1
Table 2. Drag coefficient.
Table 2. Drag coefficient.
NoComputationM = 1.05M = 1.6
C d r a g Δ C d r a g ,   % C d r a g Δ C d r a g ,   %
1Experiment0.4490.385
2ε (33) at K = 0.0010.45892.20.42199.6
3ἐ (35) at K = 0.010.45240.70.41217.1
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Struchkov, A.V.; Kozelkov, A.S.; Zhuchkov, R.N.; Volkov, K.N.; Strelets, D.Y. Implementation of Flux Limiters in Simulation of External Aerodynamic Problems on Unstructured Meshes. Fluids 2023, 8, 31. https://doi.org/10.3390/fluids8010031

AMA Style

Struchkov AV, Kozelkov AS, Zhuchkov RN, Volkov KN, Strelets DY. Implementation of Flux Limiters in Simulation of External Aerodynamic Problems on Unstructured Meshes. Fluids. 2023; 8(1):31. https://doi.org/10.3390/fluids8010031

Chicago/Turabian Style

Struchkov, A. V., A. S. Kozelkov, R. N. Zhuchkov, K. N. Volkov, and D. Yu. Strelets. 2023. "Implementation of Flux Limiters in Simulation of External Aerodynamic Problems on Unstructured Meshes" Fluids 8, no. 1: 31. https://doi.org/10.3390/fluids8010031

APA Style

Struchkov, A. V., Kozelkov, A. S., Zhuchkov, R. N., Volkov, K. N., & Strelets, D. Y. (2023). Implementation of Flux Limiters in Simulation of External Aerodynamic Problems on Unstructured Meshes. Fluids, 8(1), 31. https://doi.org/10.3390/fluids8010031

Article Metrics

Back to TopTop