Keywords

1 Introduction

In recent years, there has been a fairly frequent situation where it is necessary to calculate the flow of elongated bodies of rotation (EBR) under specific conditions. Such calculations are usually carried out for practical purposes, taking into account all technological features of the body in question. Naturally, for such calculations, there is a desire to apply some of CFD software packages, which have been widely used in recent times. However, when trying to solve practical problems with the help of such packages, there are some difficulties. The catalogs of mathematical models and finite-difference schemes used in such complexes are imperfect. The acceptability of many models for solving complex problems and determining the limits of their applicability are the subject of a separate study. This refers to the problems of flow around the elongated bodies of rotation and the implementation of turbulence simulation methods for them. For a particular class of EBR it is required to carry out a large number of test calculations to show the applicability of the chosen numerical method and the chosen model of turbulence. These methodological studies are often neglected. Therefore, a user of similar software packages encounters the need to specify a variety of variable parameters, which in practice provides an indefinite result.

In this situation, obviously, we need a computational technology that would be a kind of standard for solving the problems of flow around the EBR and would help to regulate the tunable parameters of both numerical methods and models of turbulence in various software packages. In this capacity, it was decided to recreate on the modern level the computational technology developed earlier in the Keldysh Institute of Applied Mathematics. In the late 80’s - early 90’s this computational technology allowed to make mass industrial computing for a flow around EBR with a high degree of reliability. The error of aerodynamic drag coefficients did not exceed 2–3% in comparison with the experimental results. The essence of this technology was that the aerodynamic drag coefficient Cx, was considered as a sum of three components: Cp – coefficient for inviscid flow, Cf - coefficient for viscous friction and Cd – coefficient for near wake pressure. Such an approach was widely used for industrial analysis of aerodynamic properties of EBR and proved to be very effective. The work presented is a part of the general project to create a similar technology [1, 2]. To calculate the friction coefficient, a computational technique [2] is realized. The technique is based on an approximate semi-empirical model which combines the results of experimental studies and the method of effective length. This computational technology is designed to determine the friction coefficient and estimate the characteristics of the boundary layer for EBR. To calculate the aerodynamic characteristics for inviscid flow around the elongated bodies of rotation, it was proposed to use the OpenFOAM software package (Open Source Field Operation and Manipulation CFD Toolbox) [3]. OpenFOAM is actively used in industry and in science. OpenFOAM contains a number of solvers [4,5,6,7] having different computational properties.

Therefore, it is necessary to make methodological calculations that allow to evaluate the effectiveness of these solvers for practical application. This paper presents a comparative analysis of the OpenFOAM solvers accuracy for the problem of inviscid flow around cones with different angles and different flow velocities at zero angle of attack. Tabular solutions [8] are used as an exact solution for comparison. Presented in a tabular form solutions [8] have high accuracy and for many years are used for testing the computational properties of numerical methods. It should be noted that similar comparisons of solvers were carried out in [9, 10]. However, these comparisons do not give full and clear recommendations on the accuracy of solvers.

2 Formulation of the Problem

The statement of the problem is presented in full accordance with [8], where the results of the inviscid flow around cones with different angles at various Mach numbers are considered. We consider the case of a cone placed in a uniform supersonic flow of an ideal gas at zero angle of attack α = 0° with a Mach number of 2 to 7. The body under investigation is a cone with an angle β = 10–35° in steps of 5°. Here angle β is a half of cone angle as shown in Fig. 1. The conditions of the input flow are denoted by the index “∞”, and at the output by the index ξ, since the solution is self-similar and depends on the dimensionless variable. The Euler equations system is used for the calculation. The system is supplemented by the equation of state of an ideal gas.

Fig. 1.
figure 1

Flow scheme.

Full size image

3 OpenFOAM Solvers

For comparison, 4 solvers were selected from the OpenFOAM software package:

RhoCentralFoam is based on a central-upwind scheme, which is a combination of central-difference and upwind schemes [4, 5]. The essence of the central-upwind schemes consists in a special choice of a control volume containing two types of domains: around the boundary points, the first type; around the center point, the second type. The boundaries of the control volumes of the first type are determined by means of local propagation velocities. The advantage of these schemes is that, using the appropriate technique to reduce the numerical viscosity, it is possible to achieve good solvability for discontinuous solutions — shock waves in gas dynamics, and for solutions in which viscous phenomena play a major role.

SonicFoam is based on the PISO algorithm (Pressure Implicit with Splitting of Operator) [6]. The basic idea of the PISO method is that two difference equations are used to calculate the pressure for the correction of the pressure field obtained from discrete analogs of the equations of moments and continuity. This approach is due to the fact that the velocities corrected by the first correction may not satisfy the continuity equation, therefore, a second corrector is introduced which allows us to calculate the velocities and pressures satisfying the linearized equations of momentum and continuity.

RhoPimpleFoam is based on the PIMPLE algorithm, which is a combination of the PISO and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms. An external loop is added to the PISO algorithm, thanks to which the method becomes iterative and allows to count with the Courant number greater than 1.

PisoCentralFoam is a combination of a Kurganov-Tadmor scheme [4] with the PISO algorithm [7].

For all solvers the calculations were carried out using the OpenFOAM version 2.3.0. Solver sonicFoam in the standard version does not support dynamic time step change, so the necessary corrections have been made to the code of solver. Also the calculations were made for pimpleCentralFoam solver. This solver exists only for OpenFOAM version 3.0.1 and higher. The results for this solver were similar to the results of pisoCentralFoam, so it was decided not to include these results in the tables below.

4 Computations and Results

4.1 Mesh Generation, Initial and Boundary Conditions

Figure 2 shows the computational domain. On the upper boundary indicated as “top”, the zero gradient condition for the gas dynamic functions, is specified. The same conditions are set on the right border, denoted by “outlet”. On the left border, designated as “inlet”, the parameters of the oncoming flow are set: pressure P = 101325 Pa, temperature T = 300 K, speed U from 694.5 m/s (Mach number = 2) to 2430.75 m/s (Mach number = 7). On the boundary of the cone (“cone”) for pressure and temperature, the condition of zero gradient is given, for the speed is given the condition “slip”, corresponding to the non-flow condition for the Euler equations. To model the axisymmetric geometry in the OpenFoam package, a special “wedge” condition is used for the front (“front”) and back (“back”) borders. The OpenFoam package also introduces a special “empty” boundary condition. This condition is specified in cases when calculations in a given direction are not carried out. In our case, this condition is used for the “axis” border.

Fig. 2.
figure 2

Computational domain and boundaries.

Full size image

The number of grid cells is 13200. The initial conditions correspond to the boundary conditions on the inlet edge, that is, the initial conditions are used for the parameters of the oncoming stream. The molar mass M = 28.96 kg/mol and the specific heat at constant pressure C p  = 1004 were also set.

To estimate the effect of the grid partition on the accuracy of calculations, the calculations were carried out on three grids, denoted as coarse, fine, finest. The number of cells: coarse – 3000, fine – 12000, finest – 48000.

4.2 Parameters of Solvers

In the OpenFOAM package, there are two options for approximating differential operators: directly in the solver’s code or using the fvSchemes and fvSolution configuration files. In order for the comparison to be correct, we used the same parameters where possible. In the fvSchemes file: ddtSchemes – Euler, gradSchemes – Gauss linear, divSchemes – Gauss linear, laplacianSchemes – Gauss linear corrected, interpolationSchemes– vanLeer. In the fvSolution file: solver – smoothSolver, smoother symGaussSeidel, tolerance – 1e−09, nCorrectors – 2, nNonOrthogonalCorrectors – 1.

4.3 Calculation of the Axisymmetric Flow

Figure 3 presents the steady-state flow field for pressure when using the solver rhoCentralFoam. The figure indicates that, as a result of the establishment, a qualitative picture of the flow is obtained that corresponds to the known solutions [8].

Fig. 3.
figure 3

Pressure field for steady flow.

Full size image

Tables 1, 2, 3, 4, 5, 6, 7, 8 and 9 show the results of calculations in the form of an analog of the L2 norm:

Table 1. Deviation from the exact solution for coarse grid
Full size table
Table 2. Deviation from the exact solution for fine grid
Full size table
Table 3. Deviation from the exact solution for finest grid
Full size table
Table 4. Deviation from the exact solution, U = 2M
Full size table
Table 5. Deviation from the exact solution, U = 3M
Full size table
Table 6. Deviation from the exact solution, U = 4M
Full size table
Table 7. Deviation from the exact solution, U = 5M
Full size table
Table 8. Deviation from the exact solution, U = 6M
Full size table
Table 9. Deviation from the exact solution, U = 7M
Full size table
$$ {{\sqrt {\sum\limits_{m} {\left| {y_{m} - y_{m}^{exact} } \right|^{2} } V_{m} } } \mathord{\left/ {\vphantom {{\sqrt {\sum\limits_{m} {\left| {y_{m} - y_{m}^{exact} } \right|^{2} } V_{m} } } {\sqrt {\sum\limits_{m} {\left| {y_{m}^{exact} } \right|^{2} V_{m} } } }}} \right. \kern-0pt} {\sqrt {\sum\limits_{m} {\left| {y_{m}^{exact} } \right|^{2} V_{m} } } }} $$
(1)

In Tables 1, 2 and 3 y m is the velocity U x , U y , pressure p and density ρ in the cell, V m is the cell volume for the cone angle β = 20° and the Mach number M = 2. These tables show the grid convergence for the variant considered. Grid convergence for all variants was considered similarly.

In Tables 4, 5, 6, 7, 8 and 9 ym is the pressure p in the cell, V m is the cell volume for the cone angle β = 10–35° in steps of 5° and the Mach numbers M = 2–7. The minimum values are highlighted in bold. The symbol “x” in the tables means that at a given speed and given cone angle, the solver became unstable. The values of y exactm are obtained by interpolating tabular values from [8] into grid cells. It should be noted that the authors of the tables [8] indicate the admissibility of interpolation for all parameters and table values. Further we will use abbreviations for solvers. rCF (rhoCentralFoam), pCF (pisoCentralFoam), sF (sonicFoam), rPF (rhoPimpleFoam).

Figure 4 presents a diagram of the deviation from the exact solution in the analogue of the L 2 norm for the pressure for all used solvers by the example of the problem of flow past a cone with a cone angle β = 20° with Mach number M = 2. The smallest deviation from the exact solution is shown by the solver rhoCentralFoam, the maximum deviation is shown by the solver rhoPimpleFoam.

Fig. 4.
figure 4

Deviation from the exact solution for pressure M = 2, β = 20.

Full size image

Figure 5 shows the change in the deviation from the exact solution in the analogue of the L 2 norm for the pressure for all solvers, depending on the cone angle for a fixed Mach number M = 2. The smallest deviation from the exact solution is shown by the solver rhoCentralFoam, the largest deviation with an increase in the cone angle is shown by the sonicFoam solver.

Fig. 5.
figure 5

Change in deviation from the exact solution for pressure depending on the cone angle for all solvers, M = 2.

Full size image

Figure 6 shows the dependence of the deviation on the exact solution in the analog of the L 2 norm for the pressure for the solver rhoCentralFoam with the variation of the cone angle and the initial velocity. An increase in the Mach number of the oncoming stream has the greatest effect on the increase in the deviation of the numerical result from the exact solution.

Fig. 6.
figure 6

Change in deviation from the exact solution for pressure depending on the cone angle and the velocity for the solver rhoCentralFoam.

Full size image

5 Conclusion

Using well-known problem of a supersonic inviscid flow around a cone at zero angle of attack we compared four OpenFoam solvers with the exact solution. Grid convergence was shown for all solvers. According to the results obtained, the solver rhoCentralFoam has minimal error in almost all cases. The only drawback of rhoCentralFoam is the appearance of oscillations near the surface at the head of the cone. Solver pisoCentrlFoam is in second place in accuracy, however, when using this solver, the appearance of oscillations is not observed. The methodical research can serve as a basis for selecting the OpenFoam solver for calculating the inviscid supersonic flow around the elongated bodies of rotation. The results of solvers comparison can also be useful for developers of OpenFoam software content. The results obtained made it possible to get a general idea of the calculation errors for all solvers.

In further studies it is proposed to make a similar comparison of solvers for the problem of flow around a cone with a variation of the angle of attack. It is also proposed to investigate the matrix of mutual errors for solutions obtained by different solvers by constructing elastic maps.