The need to organize and conduct a
comparative assessment of the
numerical
methods accuracy is growing significantly today. The reason for this is the
emergence of a large number of new numerical methods, solvers developed on
their basis, integrated into various software packages. At the same time,
solver developers are far from always able to organize full testing and
comparative assessment of accuracy on problems that have a reference solution.
It should be remembered that, ultimately, both software packages and the
solvers implemented in them are targeted at the end user who uses these
software tools to solve practical problems. In this case, the end user must
rely on something, and here comparative estimates of the solvers accuracy can
help him. In order to provide the most comprehensive and complete comparison,
it is necessary to carry out parametric studies. They allow comparative
assessment not for one specific problem, but for a class of problems.
Parametric numerical research is now being carried out very actively in various
fields [1-4].
A generalized computational experiment
will help to organize such numerical studies. The theory and methods of
generalized computational experiment are being actively developed at the Keldysh
Institute of Applied Mathematics of the RAS. A generalized computational
experiment is a computational technology based on the synthesis of solutions to
problems of mathematical modeling, the use of parallel technologies and methods
of visual analytics. The construction of such an experiment makes it possible
to obtain a numerical solution for a class of problems determined by variations
in the defining parameters in certain ranges. The region of the defining
parameters is subdivided into a grid. At each point of such a partition, the
problem of mathematical modeling is solved using parallel technologies. The
result of such a calculation is multidimensional data, with the help of which
valuable functionals are constructed as functions of several variables. The
obtained data and functionals are explored using scientific visualization and
visual analytics methods. The theoretical description of methods, tools and
examples of constructing a generalized computational experiment is described in
detail in the cycle of works [5-11]. With regard to the tasks of comparative
assessment of the accuracy of solvers, the purpose of the study is to obtain
comparison results in a visual form and in an analytical form.
This work considers an example of the
application of the generalized computational experiment to the problem of
accuracy comparative estimation for several solvers of the open software
package OpenFOAM [12]. The classical problem of a rarefaction wave formation in
a flow is used as the basic problem. The class of problems is determined by the
variation in the given ranges of two defining parameters - the freestream Mach
number M and the flow deflection angle β. This problem has an exact
solution described in [13-15]. This solution is used as a reference one. The
deviation from the exact solution calculated in the L1 and L2 norms is
considered in the role of the objective functional. For each solver, a
parametric problem is solved by varying the Mach number and flow deflection
angle. The results obtained represent the dependence of the objective functional
on the defining parameters for each solver. The results are presented in visual
form as surfaces and in analytical form as polynomials. The results obtained
make it possible to obtain a fairly complete idea of the accuracy of each of
the solvers participating in the comparison in such an important class of
problems as the rarefaction wave.
This work continues a large cycle of works
devoted to the development of methods for constructing a generalized
computational experiment and its application to practical problems in
computational gas dynamics [5-11,17-22]. In particular, the application of a
generalized computational experiment to problems of comparative estimation of
numerical methods accuracy is considered using examples of typical classes of
problems.
A number of characteristic directions
presented in previous works can be distinguished:
- the basic principles and methods of
constructing a generalized computational experiment for problems of
computational gas dynamics are presented in [5-10];
- a description of the main tasks of data
visualization arising in the construction of a generalized computational
experiment is presented in [11, 16];
- a description of the construction and
implementation of a generalized computational experiment for problems of
comparative assessment of numerical methods
accuracy
is presented in [17-21]. Papers [17-19] consider the problem of comparing the
accuracy of various OpenFOAM solvers when flowing around cones at an angle of
attack. This task has four defining parameters. Variations of the Mach number,
the cone half-angle, the angle of attack, and the choice of the solver are
considered. Results in a similar setting for the problem of of an oblique shock
wave formation are presented in [20, 21];
- a description of the construction of
analytical dependencies for visual images is presented in [22].
All the approaches, methods and software
tools developed at the previous stage are used in this work for a comparative
analysis of the OpenFOAM solvers accuracy for the problem of a rarefaction wave
formation.
This
work is devoted to the accuracy comparative assessment of OpenFOAM solvers when
calculating the rarefaction wave. The comparison of accuracy was carried out on
the basis of constructing a generalized computational experiment. The
well-known classical problem of two-dimensional rarefaction wave formation in
an inviscid gas flow around a flat plate at an angle of attack was considered
as the basic problem. A scheme of such a flow is shown in Figure 1. A
supersonic flow of an inviscid compressible gas runs onto a flat half-plate at
an angle of attack β, as shown in the figure. A rarefaction wave fan is
formed at the end of the plate.
Fig. 1. Flow
scheme. The area of calculations.
This
problem is widely known as the Prandtl-Mayer flow. This problem has an exact
solution, the description of which can be found in [13-15]. The exact solution
in our case acts as a standard one.
The
problem is solved within the framework of a two-dimensional system of Euler
equations for an inviscid gas. At the inlet boundary, the parameters of the
unperturbed free stream are set at the Mach number M and a certain value of
β. On the part of the lower boundary corresponding to the flat plate, the
no-slip condition is set. At the outlet boundary, the boundary conditions are
set for the zero derivatives of the gas-dynamic functions along the normal to
the boundary. At the upper boundary for the velocity components, the boundary
conditions are set similarly to the conditions for the input boundary. For the
rest of the gas dynamic functions of the upper boundary, the conditions are set
similarly to the conditions for the outlet boundary.
Parametric
studies were carried out for a Mach number equal to 2–4, with a step of 0.5,
and an angle β = 6°, 10°, 15°, 20°. Four solvers of the OpenFOAM open
software package [24-29] -
rhoCentralFoam, pisoCentralFoam, sonicFoam, and
QGDFoam
- took part in the comparative assessment of accuracy. For all
these solvers, some actions were taken to unify the calculations when setting
the boundary and initial conditions, similar to the procedures described in
[17-21].
This
section contains a short description of the solvers of the open source software
package OpenFOAM, participating in the study on comparative assessment of their
accuracy. The following solvers were used in the calculations:
1)
rhoCentralFoam (rCF) - based on the central-upwind scheme, which is a
combination of the central-difference and upwind schemes [23]. The original
Kurganov-Tadmor scheme [24] is a central difference scheme based on the
Lax-Friedrichs scheme. It is characteristic of Godunov-type schemes that the
calculation of the flow through the lateral faces of the cell requires the
solution of the Riemann problem on the decay of an arbitrary discontinuity. The
schemes with a central difference calculate the required values
at the next time step without solving the Riemann problem; they
integrate according to the Riemann fan. Consequently, there is no need to
search for the Riemann integral and expansion of the solution in terms of
characteristics, thus methods of this type are less computationally expensive.
However, the OpenFOAM software package uses the Kurganov – Noel – Petrova
scheme [25], which is a central upstream scheme. The scheme is central, since
it does not require the solution of the Riemann problem, and upstream, since
the information taken against the flow rate is used to estimate the width of
the Riemann fan. This more accurate estimate makes the schemes less dissipative
than the original Kurganov-Tadmor scheme.
2)
sonicFoam (sF) - based on the PISO (Pressure Implicit with Splitting of
Operator) algorithm [26]. The essence of splitting methods is to separate the
influence of various terms on the change in momentum in the computational cell.
The contribution of the pressure gradient is decoupled from the contribution of
spatial transport and viscous terms. When implementing the method, the predicted
and corrected values of the velocity and pressure of the flow are introduced.
Then the equations for pressure correction and momentum balance are solved
iteratively.
3)
pisoCentralFoam [27] is a hybrid method, which is a combination of a
central-uupwind scheme with the PISO algorithm. To smoothly switch between the
implicit version of the Kurganov-Tadmor method and the splitting method,
depending on the local flow parameters, a linear mixing of the flow expressions
from both methods is used. This solver is not included in the standard set of
solvers in the OpenFOAM software package. It was created by an independent team
of developers at the Ivannikov Institute for System Programming of the RAS.
4)
The QGDFoam solver [28, 29] is based on the system of quasi-gasdynamic
equations [30-32], developed by a group of researchers led by B.N.
Chetverushkin. A significant difference between the QGD approach and the
Navier-Stokes theory is the procedure of space-time averaging to determine the
main gas-dynamic quantities. Additional smoothing in time is the reason for the
appearance of additional dissipative terms in the QGD system. These dissipative
terms have the form of the second spatial derivatives and are proportional to a
small parameter having the dimension of time. The presence of a controllable
parameter with dissipative terms makes it possible to successfully suppress
unwanted oscillations in the numerical simulation of problems with
discontinuities. Additional terms appear to be efficient regularizers. This
solver is also not included in the standard set of solvers in the OpenFOAM
software package. It was developed at Keldysh Institute of Applied Mathematics
of the RAS under the leadership of professor T.G. Elizarova.
This section presents the results of the
calculations. Figure 2 shows the pressure field in the computational domain for
the exact solution calculated in accordance with [13-15]. Figure 3 shows a
typical pressure field for all solvers. The figure shows the pressure field
calculated using the QGDFoam solver. Note that all calculations for the QGDFoam
solver were carried out at the value of the parameter α, which controls
the dissipative properties of the solver, equal to α = 0.1. The figures
show that the calculated pressure field presents a more blurry picture compared
to the pressure field for an exact solution.
Fig. 2. Pressure
field for an exact solution.
Fig. 3. Pressure
field calculated using the QGDFoam solver.
Figure 4 below shows a
comparison of the pressure calculated for all solvers with the exact solution
along the horizontal line ÀÀ1 (Fig. 1). The horizontal line is located at a
height of 0.1 from the plate. Figure 5 shows the same comparison in close-up.
Fig. 4. Comparison
of solutions along the horizontal line AA1 for all solvers with exact solution
Fig. 5. Comparison
of solutions along the horizontal line AA1 for all solvers with the exact
solution - close-up.
Figures 4 and 5 show that
the sonicFoam solver deviates most from the exact solution. The rhoCentralFoam,
pisoCentralFoam and QGDFoam solvers show very close results. It should be noted
that the curve for pisoCentralFoam shows constant noticeable fluctuations,
which indirectly indicates the presence of unwanted non-destructive
oscillations in the solution.
Thus, we have a
comparative assessment of the accuracy of the solvers along the line. However,
this is clearly not enough. Let us construct estimates of the deviation from
the exact solution for the entire computational domain in the L1 and L2 norms.
To do this, calculate the relative error Err for the L1 and L2 standards as
follows:
Here
ym
is the pressure
p,
Sm
is the cell area. The
ymexact
values are obtained from the exact solution of the problem of the rarefaction
wave formation [13-15]. The parameters of the incident flow vary as follows:
the angle of flow deflection β = 6°, 10°, 15°, 20°, Mach number M = 2,
2.5, 3, 3.5, 4. Next, a generalized computational experiment is implemented,
based on the results of which we calculate the error for of each solver for
each combination (M, β). This enables us to construct for both norms the
error in the form of a function of two variables Err (X1, Y1), where X1 = M and
Y1 = β. These designations will be used in all subsequent figures.
Next, consider a typical
picture resulting from the construction of the error Err (X1, Y1) in different
norms. For example, consider the results for the QGDFoam solver.
Figure 6 shows the
deviation from the exact solution in the L1 norm for the QGDFoam solver in the
ranges of the governing parameters. Similarly, Figure 7 shows the results for
the L2 norm. By the appearance of the surfaces shown in these figures, it can
be argued that the error, in comparison with the exact solution, increases both
with an increase in the Mach number and with an increase in the angle of flow
deflection.
Fig. 6. The
deviation from the exact solution in the L1 norm for the QGDFoam solver in the
ranges of variation of the governing parameters
Fig. 7. The
deviation from the exact solution in the L2 norm for the QGDFoam solver in the
ranges of variation of the governing parameters
Now let's look at
comparing all four solvers at different error norms. Figure 8 shows the
deviations from the exact solution in the L1 norm for all solvers in the ranges
of the governing parameters. Similarly, Figure 9 shows deviations from the
exact solution in the L2 norm for all solvers.
Fig. 8. The
deviation from the exact solution in the L1 norm for all solvers in the ranges
of the governing parameters.
Fig. 9. The
deviation from the exact solution in the L2 norm for all solvers in the ranges of
the governing parameters.
As can be seen from Figures 8 and 9,
qualitatively similar results were obtained for both norms. The rhoCentralFoam
solver provides the smallest deviation from the exact solution. The next in
order to deviate from the exact solution is the QGDFoam solver. The
pisoCentralFoam solver gives slightly worse results. It should be noted that in
the selected ranges of variation of the Mach number and the angle of
deflection, the error surfaces for the QGDFoam and pisoCentralFoam solvers are
located close to each other. The greatest deviation from the exact solution is
observed in the results for SonicFoam solver.
Further, all the obtained surfaces can be
presented in an analytical form according to the method proposed and described
in [22].
According to [22], to approximate curved surfaces we
will use second-order polynomials, where the error for the surface under
consideration can be represented as a function of the following form:
Err
= AX1 + BY1 + CX12
+ DY12
+ EX1Y1 + F
Here also X1 is the Mach number M, Y1 is the flow
deflection angle β, Err
is the error of comparison with the
exact solution in the L1 or L2 norm. Coefficients A, B, C, D. E, F are
calculated for a specific surface.
Here is an example of a similar analytical dependence
for a surface corresponding to the results for solver QGDFoam in the L2 norm
(Fig. 9).
Approximating the required surface by a polynomial of
the second order using the least squares method, we obtain the following values
for the coefficients
A
= 0.0007426759754917806
B
= 0.0005021159520976077
C
= 0.00020442857142857106
D
= - 0.000012584191705984598
E
= - 0.0000167566591422123
F
= 0.0038062569399619252
For a more general comparative assessment,
we will calculate the Error Index (EI), similarly to that proposed in [33].
According to [33], the Error Index (EI) is an average value for each error
surface. The results for each solver in the L2 norm in accordance with Figure 9
are presented in Table 1.
Table 1.
Error Index values for the problem of rarefaction wave formation
|
Solver
|
rCF
|
QGDF
|
pCF
|
sF
|
|
Error
Index
|
0.00894
|
0.01134
|
0.01182
|
0.02285
|
Table 1 shows that the values of the error
index EI fully correspond to the relative positions of the numerical results
presented in Figure 9.
Thus, the results obtained in the form of
visual representations of error surfaces, their analytical representations and
calculated error indices allow the user of the above solvers to get a full
understanding of their accuracy in the class of problems of a rarefaction wave
formation.
The paper considers the processing and
visualization of the results of a generalized computational experiment using a
specific example. As an example, we used the results of a generalized
computational experiment on the comparative assessment of the accuracy of four
solvers of the OpenFOAM open software package. As a basic problem, the problem
of the formation of a rarefaction wave when a supersonic flow of an inviscid
gas flows around a plate at an angle of attack is used. The space of the
defining parameters is set by the variation of two parameters in the selected
ranges - the Mach number and the angle of flow deflection. For each solver, a
discrete solution is obtained in the form of a dependence of the error on the
defining parameters. A visual representation of the solution is shown, the
example of analytical form of the solution is given as 2nd order polynomial.On
the whole, the results obtained make it possible to estimate the degree of
accuracy of each solver for the class of problems under consideration.
The study was supported by a grant from
the Russian Science Foundation ¹ 19-11-00169,
https://rscf.ru/project/19-11-00169/
[1]
E.
V. Myshenkov, V. I. Myshenkov, “Numerical simulation of efflux from
Znamenskii's nozzle: Parametric investigations”, TVT, 37:1 (1999), 142–149;
High Temperature, 37:1 (1999), 136–143
[2]
K.
N. Volkov, V. N. Emelyanov, M. S. Yakovchuk, “Multiparameter optimization of
operating control by the trust vector based on the jet injection into the
supersonic part of a nozzle”, Num. Meth. Prog., 19:2 (2018), 158–172.
[3]
Rahim
M.F.A., Jaafar A.A., Mamat R., Taha Z. (2020) Parametric Study of CNG-DI Engine
Operational Parameters by Using Analytical Vehicle Model. In: Osman Zahid M.,
Abd. Aziz R., Yusoff A., Mat Yahya N., Abdul Aziz F., Yazid Abu M. (eds)
iMEC-APCOMS 2019. iMEC-APCOMS 2019. Lecture Notes in Mechanical Engineering.
Springer, Singapore.
https://doi.org/10.1007/978-981-15-0950-6_93
[4]
Zangeneh,
Rozie (2021) Parametric Study of Separation and Reattachment in Transonic
Airfoil Flows, AIAA Journal, DOI: 10.2514/1.J060520
[5]
Bondarev
A.E. Analysis of Space-Time Flow Structures by Optimization and Visualization
Methods // Transactions on Computational Science XIX, LNCS 7870, pp. 158-168.
[6]
Bondarev
A.E, Galaktionov V.A. Parametric Optimizing Analysis of Unsteady Structures and
Visualization of Multidimensional Data // International Journal of Modeling,
Simulation and Scientific Computing, 2013, V.04, N supp01, 13 p., DOI
10.1142/S1793962313410043
[7]
Bondarev
A.E. Multidimensional data analysis in CFD problems / Scientific Visualization.
V.6, ¹ 5, p.59-66, 2014.
[8]
Bondarev
A.E., Galaktionov V.A. Multidimensional data analysis and visualization for
time-dependent CFD problems // Programming and Computer Software, 2015, Vol.
41, No. 5, pp. 247–252, DOI:10.1134/S0361768815050023
[9]
Bondarev
A.E. On the Construction of the Generalized Numerical Experiment in Fluid Dynamics
// Mathematica Montisnigri, Vol. XLII, 2018, p. 52-64.
[10]
A.E.
Bondarev, V.A. Galaktionov. Generalized Computational Experiment and Visual
Analysis of Multidimensional Data (2019). Scientific Visualization 11.4: 102 -
114, DOI: 10.26583/sv.11.4.09
[11]
A.E.
Bondarev . On visualization problems in a generalized computational experiment
(2019). Scientific Visualization 11.2: 156 - 162, DOI: 10.26583/sv.11.2.12
[12]
OpenFOAM Foundation. [Online]. Available: http://www.openfoam.org
[13]
Loycyanskiy L.G. 1987 Mechanics of liquid and gas (Mehanika zhidkosti I
gaza) (Moscow: Nauka) p 840 [In Russian]
[14]
Abramovich G.N. 1976 Applied Gas Dynamics (Prikladnaya gazovaya dinamika)
(Moscow: Nauka) p 888 [In Russian]
[15]
Bondarev E.N., Dubasov V.T., Ryzhov Y.A., Svirschevsky S.B. and Semenchikov N.V.
1993 Aerigidromechanics (Aerogidromehanika). (Moscow: Mashinostroenie) p
608 [In Russian]
[16]
Bondarev A.E., Galaktionov V.A., Chechetkin V. M. Analysis of the Development
Concepts and Methods of Visual Data Representation in Computational Physics /
Computational Mathematics and Mathematical Physics, 2011, Vol. 51, No. 4, pp.
624–636.
[17]
Alexander E. Bondarev and Artem E. Kuvshinnikov. Analysis of the Accuracy of
OpenFOAM Solvers for the Problem of Supersonic Flow Around a Cone // LNCS
10862, pp. 221–230, 2018.
DOI:10.1007/978-3-319-93713-7_18
[18]
Bondarev A., Kuvshinnikov A. Comparative Estimation of QGDFoam
Solver Accuracy for Inviscid Flow Around a Cone // IEEE The Proceedings of the
2018 Ivannikov ISPRAS Open Conference (ISPRAS-2018). P. 82-87,
DOI:10.1109/ISPRAS.2018.00019
[19]
A.E. Bondarev, A.E. Kuvshinnikov. Analysis of the behavior of OpenFOAM solvers
for 3D problem of supersonic flow around a cone at an angle of attack // CEUR
Workshop Proceedings, V. 2763, 2020, p. 48-51. DOI:
10.30987/conferencearticle_5fce2771320ef0.90086903
[20]
Alekseev A.K., Bondarev A.E., Kuvshinnikov A.E. Comparative analysis of
the accuracy of OpenFOAM solvers for the oblique shock wave problem //
Matematica Montisnigri, 2019, vol.
XLV, p. 95-105 DOI:
10.20948/mathmontis-2019-45-8
[21]
Bondarev A., Kuvshinnikov A. Parametric Study of the Accuracy of OpenFOAM
Solvers for the Oblique Shock Wave Problem // IEEE The Proceedings of the 2019
Ivannikov ISPRAS Open Conference (ISPRAS-2019) 2019. P. 108-112,
DOI:10.1109/ISPRAS47671.2019.00023
[22]
A.E. Bondarev. Processing of Visual Results of a Generalized Computational
Experiment for the Problem of Supersonic Flow Around a Cone at an Angle of
Attack (2021). Scientific Visualization 13.2: 104 - 116, DOI:
10.26583/sv.13.2.08
[23]
C. Greenshields, H. Wellerr, L. Gasparini, J. Reese “Implementation of
semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite
volume framework, for high-speed viscous flows”, Int. J. Numer. Meth. Fluids,
63(1), 1–21 (2010). https://doi.org/10.1002/fld.2069
[24]
A. Kurganov, E. Tadmor “New high-resolution central schemes for nonlinear
conservation laws and convection-diffusion equations”, J. Comput. Phys.,
160(1), 241–282 (2000). https://doi.org/10.1006/jcph.2000.6459
[25]
A. Kurganov, S. Noelle, G. Petrova, “Semidiscrete central-upwind schemes for
hyperbolic conservation laws and Hamilton–Jacobi equations”, SIAM J. on Sci.
Comp., 23(3), 707–740 (2001). https://doi.org/10.1137/S1064827500373413
[26]
R. Issa, “Solution of the implicit discretized fluid flow equations by operator
splitting”, J. Comput. Phys., 62(1), 40–65 (1986).
https://doi.org/10.1016/0021-9991(86)90099-9
[27]
M. Kraposhin, A. Bovtrikova, S. Strijhak, “Adaptation of Kurganov-Tadmor
numerical scheme for applying in combination with the PISO method in numerical
simulation of flows in a wide range of Mach numbers”, Procedia Computer
Science, 66, 43–52 (2015). doi:10.1016/j.procs.2015.11.007.
[28]
M. V. Kraposhin, D. A. Ryazanov, E. V. Smimova, T. G. Elizarova, M. A.
Istomina, "Development of OpenFOAM Solver for Compressible Viscous Flows
Simulation Using Quasi-Gas Dynamic Equations", Ivannikov ISPRAS Open
Conference, 117-123 (2017).
[29]
M. A. Istomina, "About realization of one-dimensional quasi-gas dynamic
algorithm in the open program OpenFOAM complex", Preprint IPM No. 1
(Moscow: KIAM), (2018).
[30]
B. N. Chetverushkin, T. N. Elizarova “Kinetic algorithms for calculating gas
dynamic flows,” U.S.S.R. Comput. Math. Math. Phys., 5 (5), 164–169 (1985).
[31]
B. N. Chetverushkin, Kinetic schemes and Quasi-Gas Dynamic system of equations,
CIMNE, Barcelona (2008).
[32]
T. G. Elizarova, Quasi-Gas Dynamic Equations, Springer (2009).
[33]
Aleksey Alekseev, Alexander Bondarev, Vladimir Galaktionov, Artem Kuvshinnikov,
Lev Shapiro. On Applying of Generalized Computational Experiment to Numerical
Methods Verification // CEUR Workshop Proceedings, V. 2744, 2020, p.
paper19-1 — paper19-12, DOI: 10.51130/graphicon-2020-2-3-19
