This paper is an extended version of a contribution presented
at the Graphiñon 2025 conference.

1. Introduction

The study of aerodynamic phenomena related to shock wave propagation represents a fundamental problem in the field of gas dynamics and computational fluid dynamics. Scenarios where the interaction of shock waves with geometric obstacles leads to the formation of complex wave structures, such as reflected shock waves and shock wave chains, are of particular interest. Such phenomena are critical for understanding and designing a wide range of engineering systems. Numerical simulation of these processes is becoming an indispensable tool for detailed study of flow dynamics, parameter distribution, and shock wave characteristics that are difficult or impossible to investigate experimentally [1–3].

In recent decades, computational experimentation has become an indispensable tool in the arsenal of engineers and researchers working on aerodynamic problems. The ability to model complex flows on a computer allows for detailed analysis that is often inaccessible or prohibitively expensive when relying solely on laboratory or full-scale experiments. Computational experiments enable the investigation of a wide range of parameters, variation of geometry, study of unsteady processes, and acquisition of detailed information about flow parameters at any point within the computational domain. However, simulating high-speed flows with shock waves presents a particularly challenging task for numerical methods. Shock waves are characterized by sharp gradients, requiring the use of high-precision numerical schemes capable of correctly resolving these discontinuities without significant artifacts such as oscillations or artificial smoothing. An incorrect choice of numerical approach can lead to misdetermination of shock wave position, intensity, and distortion of other important flow characteristics. This underscores the importance of careful selection and testing of solvers used for such tasks [4, 5].

Numerical simulation of such processes requires the application of robust and accurate computational methods capable of adequately describing the discontinuous solutions characteristic of shock waves. The OpenFOAM software package [6, 7], being a powerful open-source tool, offers a wide range of solvers, each with its own characteristics in terms of accuracy, stability, and computational efficiency for solving different types of problems. The selection of a suitable solver is critically important for obtaining reliable results, especially in shock wave problems where high resolution and accurate reproduction of shock wave fronts are required. Understanding the advantages and disadvantages of each solver in the context of this problem will allow for recommendations of the most suitable tools for further research in aerodynamics where shock waves occur.

This work is a continuation of a series of the authors' publications. Previously, reference problems considered included: the formation of a shock wave [8], the formation of a two-dimensional expansion wave [9], and the flow past a cone with a spherical nose [10]. Through the application of a generalized computational experiment [11, 12], practitioners gain the ability to confidently navigate the broad spectrum of developed numerical methods. This allows them to select the most accurate and efficient solutions for their calculations. The methodology itself involves studying a problem by discretizing its defining parameters within a specific range, followed by parametric analysis and visualization of multidimensional results.

2. Problem Statement

A high-speed gas flow with Mach number M flows from left to right (Fig. 1).

 

Figure 1. Flow schematic

 

Let's denote the free stream region as zone “0”. The flow encounters a wedge “a” with an angle β and generates an oblique shock wave; the region behind this shock wave is denoted as zone “1”. The flow in zone “1” is parallel to the wedge “a”, and conditions are set by the oblique shock relations, for example, as given in [13]. The oblique shock then strikes a solid wall and reflects from it, creating a new shock wave. The flow behind the reflected shock wave is denoted as zone “2”. Since the flow in zone “1” is parallel to the wedge “a”, it impinges on the solid wall at an angle “a”, as shown by the dashed white line. The flow in zone “2” is parallel to the solid wall, and conditions in zone “2” are determined by the oblique shock relations, with upstream conditions corresponding to those in zone “1”. The reflected shock wave itself reflects off the wedge, creating a chain of shock waves in the channel formed by the wedge and the solid wall. With each shock wave passage and reflection, the flow Mach number decreases. Eventually, the Mach number in some region becomes too low for an oblique shock wave to form, and a final, normal shock wave is generated.

The defining parameters of the problem in terms of the generalized computational experiment are the Mach number M and the wedge angle β. The Mach number varied from 1.8 to 2.0 with a step of 0.1, and the wedge angle β from 5° to 10° with a step of 2.5°.

3. Computational Setup

Four solvers participated in the comparison: two standard solvers — rhoCentralFoam and sonicFoam, as well as two custom ones — pisoCentralFoam [14] and QGDFoam [15]. The latter two solvers were developed by teams from the Institute for System Programming of the Russian Academy of Sciences and the Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences.

These solvers differ significantly: rhoCentralFoam uses the Kurganov-Tadmor scheme, a central-upwind Godunov-type scheme [16]; sonicFoam uses the PIMPLE algorithm, which includes a splitting method [17]; pisoCentralFoam is a hybrid method using both the Kurganov-Tadmor scheme and the PIMPLE method [18]; QGDFoam is based on the quasi-gas-dynamic system of equations [19, 20].

The computational domain is divided into cells. The OpenFOAM package requires specification of boundary and initial conditions for solving. At the inlet boundary “inlet”, parameters of the undisturbed incoming flow are specified (pressure P = 101325 Pa, temperature T = 300 K, x-component of velocity U_x varies from 625.05 m/s to 694.5 m/s, y-component of velocity U_y equals 0 m/s). At the outlet boundary “outlet”, boundary conditions of zero derivatives of gas-dynamic functions normal to the boundary are specified. At the wedge boundary “wedge” and at the upper boundary “top”, a zero gradient condition is specified for pressure and temperature, and a “slip” condition is specified for velocity, corresponding to the no-penetration condition for Euler equations. For the front “front” and back “back” boundaries, a special “empty” condition is used. This condition is specified in cases where calculations in the given direction are not performed, since we are solving a two-dimensional problem. The computational domain scheme for a wedge with angle β = 10° is shown in Fig. 2. It should be noted that in the indicated image, for clarity, the mesh is coarser than in actual calculations.

 

Figure 2. Computational domain layout

 

The number of mesh cells depends on the inclination angle, since the number of oblique shock reflections depends on it: 50000 for β = 5°, 32500 for β = 7.5°, and 25000 for β = 10°. Initial conditions correspond to boundary conditions at the inlet face, i.e., incoming flow parameters are used as initial conditions. In the QGDFoam solver, a smoothing coefficient α = 0.1 over the entire computational domain was also specified as an initial condition. Molar mass M = 28.96 and specific heat capacity at constant pressure Cp = 1004 were also specified.

OpenFOAM stands out among other software packages in that simulation management is performed through text files. This method provides significant flexibility: it allows easy automation of calculation launches, adjustment of modeling parameters, and analysis of obtained data.

A unified approach to conducting calculations is extremely important for comparing solvers. It ensures that all of them will be tested under identical conditions, making the assessment of their performance and accuracy objective. When calculation methodologies, meshes, boundary conditions, and physical models are standardized, the obtained results become comparable and trustworthy. This allows researchers to isolate the influence of extraneous factors and focus on the characteristics of a specific solver. Moreover, uniformity in testing helps better understand the advantages and disadvantages of each solver, which, in turn, facilitates the selection of the most appropriate tool for solving a specific engineering problem.

In working with OpenFOAM, we applied the same settings for the fvSchemes and fvSolution configuration files as in the authors' previous works.

4. Experimental Results

Flow patterns are presented in Figs. 3–5 as pressure, density, and velocity magnitude distributions in the computational domain. The presented pressure distribution was obtained using the rhoCentralFoam solver. Solution breakdown was not observed for any of the solvers, which indicates high stabilizing properties of all solvers participating in the study.

 

Figure 3. Steady-state flow pressure field for the rhoCentralFoam solver, β = 5°

 

Figure 4. Steady-state flow density field for the rhoCentralFoam solver, β = 5°

 

Figure 5. Steady-state flow velocity magnitude field for the rhoCentralFoam solver, β = 5°

 

We will construct error estimates relative to the exact solution over the entire computational domain using an L₂ norm analogue. For this purpose, we define the relative error Err for the L₂ norm analogue as follows:

.

Here y_m is the value of the investigated quantity (pressure, density and velocity magnitude), V_m is the cell volume. The values y_m^exact are obtained from the analytical solution of the problem [1, 3]. The comparative accuracy analysis involved the solvers sonicFoam, QGDFoam, rhoCentralFoam, and pisoCentralFoam. Hereinafter in the tables, abbreviated designations are used for the solvers: rCF (rhoCentralFoam), pCF (pisoCentralFoam), sF (sonicFoam), QGDF (QGDFoam). The values of deviation from the exact solution over the entire computational domain are calculated for pressure p, density ρ, and velocity magnitude U and are presented in Tables 1–3. The smallest values in each row are highlighted in bold.

Table 1. Errors for m=1.8

Table 2. Errors for m=1.9

Table 3. Errors for m=2.0

For table analysis, let us visualize the data as error surfaces. The result is presented in Figs. 6–8.

 

Figure 6. Error surfaces of the considered solvers for pressure

 

Figure 7. Error surfaces of the considered solvers for density

 

Figure 8. Error surfaces of the considered solvers for velocity magnitude

 

Analysis of the presented tables allows drawing a number of important conclusions regarding their accuracy and dependence on input parameters. First and foremost, it should be noted that the rhoCentralFoam solver consistently demonstrates the smallest error values for both pressure and density calculations in all considered scenarios without exception, which unambiguously indicates its superiority in accuracy compared to other investigated solvers. This superiority manifests in absolute error values, but not in relative stability to changes in problem conditions. Considering the dependence of solver accuracy on wedge angle, a clear general trend can be traced: with increasing wedge angle, numerical solution errors for all solvers tend to increase. This is a natural phenomenon due to flow physics: increasing wedge angle leads to higher flow parameter gradients, formation of more intense shock waves, and, consequently, to problem complexity for numerical methods. The sonicFoam solver, however, demonstrates the smallest absolute error increase with increasing wedge angle. For example, at M = 1.8, transitioning from 5° to 10° wedge angle for pressure leads to an error increase of 0.006294 for rhoCentralFoam (from 0.017524 to 0.023818), while for sonicFoam this increase is 0.002243 (from 0.030166 to 0.032409), and for QGDFoam — 0.007186 (from 0.019833 to 0.027019). This emphasizes its greater stability to geometric changes. A similar pattern is observed for density calculations, where sonicFoam also demonstrates the smallest error growth.

The influence of Mach number on solver accuracy is also significant. As can be seen from the tables, increasing Mach number leads to increased errors for all solvers. This is due to the increasing role of flow compressibility, which requires more accurate numerical schemes for adequate description of phenomena such as shock waves, especially at high speeds. The rhoCentralFoam solver, despite the growth of absolute error values with increasing Mach number, as in the case with wedge angle, maintains its accuracy advantage.

Comparing the remaining solvers, it can be noted that pisoCentralFoam shows results very close to rhoCentralFoam, but typically with somewhat higher errors. It also demonstrates a similar tendency to increase errors with increasing Mach number and wedge angle. QGDFoam occupies an intermediate position; its results are often comparable to pisoCentralFoam, but in some cases it may be inferior to it in accuracy. The sonicFoam solver systematically demonstrates the highest errors among all investigated solvers. It is also important to note that errors in velocity magnitude calculations are generally lower, followed by errors in density calculations, and maximum errors occur in pressure calculations, for all solvers and under all conditions. This may be related to the magnitude of parameter changes (for M = 2 and β= 10°, pressure increases by a factor of 2.8, density by 2.08, and velocity magnitude decreases by a factor of 1.55).

Thus, based on detailed analysis of the presented data, it can be concluded that the rCF solver is the most reliable and accurate choice for solving the wedge flow problem in the considered parameter range, due to its consistently low error and lower sensitivity to increases in Mach number and wedge angle.

5. Conclusion

This work confirmed the possibility of successful application of OpenFOAM for modeling oblique shock wave chains. Comparative evaluation of solvers showed that rhoCentralFoam is the most preferable choice for achieving high accuracy, although other investigated solvers also solve the posed problem. The research results contribute to understanding the applicability of CFD methods in aerodynamic calculations. The obtained data have important practical significance for engineers and researchers working in the field of aerodynamics. This study will help them make informed decisions when choosing numerical methods for their tasks, and opens the way to further development and application of OpenFOAM in advanced areas of science and technology.

Acknowledgments

Calculations were performed on the hybrid supercomputer K-100 installed in the Supercomputer Centre of Collective Usage of KIAM RAS.

References

1. Anderson J. D. Modern Compressible Flow: With Historical Perspective. McGraw-Hill Education, 2021. 778 p.

2. Anderson J. D. Fundamentals of Aerodynamics. McGraw-Hill Education, 2016. 1152 p.

3. Chernuy G. G. Introduction to Hypersonic Flow. Academic Press, 1961. 276 p.

4. Toro E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Heidelberg: Springer Berlin, 2009. 724 p.

5. Ferziger J. H., Peric M., Street R. L. Computational Methods for Fluid Dynamics. Springer Cham, 2009. 724 p.

6. OpenFOAM Foundation: [Online]. URL: http://www.openfoam.org (Accessed: 02.07.2025).

7. Jasak H. OpenFOAM: Open source CFD in research and industry // Int. J. Nav. Archit. Ocean Eng. 2009. Vol. 1. P. 89–94.

8. Alekseev A. K., Bondarev A. E., Kuvshinnikov A. E. Comparative analysis of the accuracy of OpenFOAM solvers for the oblique shock wave problem // Mathematica Montisnigri, 2019, Vol. XLV. P. 95-105.

9. Bondarev A. E., Kuvshinnikov A. E. Analysis and Visualization of the Computational Experiments Results on the Comparative Assessment of OpenFOAM Solvers Accuracy for a Rarefaction Wave Problem // Scientific Visualization. 2021. Vol. 13. ¹ 3. P. 34–46.

10. Bondarev A. E., Kuvshinnikov A.E. Integrating Scientific Visualization in the Assessment of OpenFOAM Solvers for the Flow Around a Spherically Blunted Cone // Scientific Visualization. 2024. V. 16. ¹ 4. P. 25–36.

11. Bondarev A. E. On the Construction of the Generalized Numerical Experiment in Fluid Dynamics // Mathematica Montisnigri. 2018. Vol. XLII. P. 52–64.

12. Bondarev A. E., Galaktionov V.A. Generalized Computational Experiment and Visual Analysis of Multidimensional Data // Scientific Visualization. 2019. V. 11. ¹ 4. P. 102–114.

13. Deich M. E. Tehnicheskaja gazodinamika. Moscow–Leningrad: Gosjenergoizdat, 1961. 670 p. [In Russian]

14. United collection of hybrid Central solvers — one-phase, two-phase and multicomponent versions: [Online]. URL: https://github.com/unicfdlab/hybridCentralSolvers (Accessed: 10.07.2025).

15. OpenFOAM framework for simulation of fluid flows using regularized (QGD/QHD) equations: [Online]. URL: https://github.com/unicfdlab/QGDsolver (Accessed: 10.07.2025).

16. Kurganov A., Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations // J. Comput. Phys. 2000. Vol. 160. ¹ 1. P. 241–282.

17. Issa R. Solution of the implicit discretized fluid flow equations by operator splitting // J. Comput. Phys. 1986. Vol. 62. ¹ 1. P. 40–65.

18. Kraposhin M. V., Banholzer M., Pfitzner M., Marchevsky I. K. A hybrid pressure-based solver for nonideal single-phase fluid flows at all speeds // Int. J. Numer. Meth. Fluids. 2018. Vol. 88. ¹ 2. P. 79–99.

19. Istomina M.A. About realization of one-dimensional quasi-gas dynamic algorithm in the open program OpenFOAM complex // KIAM Preprint ¹ 1, Moscow, 2018.

20. Chetverushkin B.N. Kinetic Schemes and Quasi-Gas Dynamic System of Equations. CIMNE Barcelona, Spain, 2008, 298 p.