CC BY
0Russian Journal of Nonlinear Dynamics, 2022, vol. 18, no. 3, pp. 349-365. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd220302
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 76F65
Direct Numerical Simulation of Aerodynamic Flows Based on Integration of the Navier - Stokes Equations
The results of the theoretical solution of aerodynamic problems based on direct numerical simulation by integrating the Navier-Stokes equations without involving additional models and empirical constants are shown. Modern approaches to the theoretical study of high-speed flows are determined. The advantages, problems, development trends and scientific directions of research on various approaches are revealed. The advantages and disadvantages of the direct numerical simulation are analyzed. The velocity vectors of laminar and transient flows in a rectangular channel with a sudden expansion at the inlet are presented in different planes. The convergence of the method is studied when the computational domain is quantized in space. It is discovered that fast relaminarization is characteristic of transitional flows. A mathematical model for calculating bottom drag is presented. The numerical results are compared with the data of physical experiments and the results of other methods. It is shown that the results of simulation based on DNS are not inferior in accuracy to RANS and LES results. The results of a parametric study of a transonic flow around a profile are presented. The high-speed buffet onset is investigated. The distribution surfaces of the velocity pulsation energy generation are shown. The frequency of self-oscillations is determined on the basis of spectral analysis.
Keywords: direct numerical simulation, Navier-Stokes equations, transient flows, base drag, baffet onset
1. Introduction
The existence of high differences in the form of movement of high-speed flows and slow flows was noticed long ago. But the beginning of the theory explaining the differences in the regimes
Received April 18, 2022 Accepted July 12, 2022
Alexey M. Lipanov aml35@yandex.ru
Keldysh Institute of Applied Mathematics Miusskaya pl. 4, Moscow, 125047 Russia
Sergey A. Karskanov ser@udman.ru
Udmurt Federal Research Center UB RAS ul. T. Baramzinoi 34, Izhevsk, 426067 Russia
A. M. Lipanov, S. A. Karskanov
of the flows was laid in the works of Osborne Reynolds at the end of the nineteenth century [1-3]. Further, for more than a hundred years, many mathematicians and mechanics have been searching for forms of mathematical description of high-speed turbulent flows. Such interest and amount of efforts are explained by the exceptional complexity of the physical processes occurring in a turbulent flow and the fact that such a flow regime is often implemented in many technical systems and applications.
On the other hand, unlike many physical phenomena for which there is no clear mathematical description, the flows of Newtonian media follow the classical Navier-Stokes equations, which has been repeatedly confirmed experimentally. From this point of view, the problem of mathematical modeling of viscous medium flows seems to be solved, however, the nonlinearity of the equations and the increasing requirements for computer technology with an increase in the Reynolds number leave the problem relevant to this day.
Even laminar flows, which, it would seem, have long been studied theoretically, are investigated on the basis of assumptions that greatly simplify the system of Navier-Stokes equations. Exact solutions for Poiseuille flows, between two parallel plates, and Couette, when one of the plates is stationary and the other moves at a constant speed, are obtained under the essential assumption of the incompressibility of the medium [4-6]. In addition, for all other more complex exact solutions, it is assumed that the medium flows are stationary. However, the prevailing majority of flows that have to be dealt with in practice are not idealized. And it is not possible to solve the Navier-Stokes equations analytically because of their nonlinearity.
At present, there are three approaches to the numerical simulation of high-speed viscous (more often turbulent) flows. The most traditional approach is based on solving equations arising from the application of Reynolds-averaged Navier-Stokes (RANS) equations and the use of various turbulence models. This method is the least demanding on computing resources. There are many varieties of turbulence models that give good results in comparison with the experiment. Semi-empirical models resolve the boundary layer well, and it is possible to obtain fluctuating characteristics in the boundary layer with acceptable accuracy. In the 1970s-1980s, great hopes were associated with this approach, international conferences were held, which received the unofficial name of "Olympiad of turbulence models" [7]. The purpose of the work seemed clear — to create a universal semi-empirical model of turbulence. However, time passed, and a universal model capable of resolving any aerodynamic flows did not appear. Moreover, it became clear that the construction of such a universal semi-empirical model of turbulence for calculating separated flows is an extremely difficult and perhaps even unsolvable problem [8].
Currently, the most popular RANS model with one additional equation is the Spalart-Allmaras (SA) model [9]. The model works well for supersonic and transonic flows around airfoils and flows in the boundary layer. However, it is not suitable for calculating shear flows and damped turbulence. In this case, models with two additional equations are used. The SST (Shear Stress Transport) model of Menter [10], which allows one to calculate more complex separated flows, has gained the most popularity.
The method of Large Eddy Simulation (LES) involves an accurate calculation of the momentum and heat transfer only by large, energetically important structures. Structures larger than the filter size are resolved exactly. Smaller eddies are simulated based on various subgrid models [11]. This method is especially attractive for the simulation of thermoconvective flows, since it adequately reproduces the time evolution of large-scale vortex structures that determine convection. However, modeling of turbulent flows in the presence of solid boundaries based on the LES method is accompanied by requirements for the grid resolution of near-wall regions in which
relatively small structures are present. This imposes significant limitations on the applicability of this approach.
It should be noted that the desire to overcome the limitations of the RANS and LES methods has led to the emergence of a hybrid approach called the Detached Eddy Simulation (DES) [12]. Taking all the "best" from RANS and LES, the hybrid method yielded excellent results. In regions of attached flow (with small vortices), RANS is used; in separated regions (with large vortices), LES is used. It is with the help of this approach that it is possible to solve well the widespread problems (with complex volumetric geometry) of the flow around vehicles or their individual parts. However, the DES method is not without its drawbacks. Still, the model contains nonuniversal empirical constants, which limits the applicability of the method. For example, the high natural viscosity of the retarded flow region can be artificially suppressed. Weakly unstable flows, when the development of unsteady structures occurs slowly, are also not amenable to empiricism. Nevertheless, work in the field of development and improvement of the DES approach is being actively carried out, and the results of this work seem to be very promising [13].
Separately, it is worth mentioning the tasks of aeroacoustics. It turns out that high-speed perturbed flow influences not only the body with which it directly interacts, but also a distant object by means of acoustic vibrations. The most common example of this phenomenon is the loud noise from aircraft engines. The task of reducing the noise generated by turbulence is extremely relevant. Mathematical models for describing the noise impact involve not only the solution of gas-dynamic equations, but also the subsequent finding of an aeroacoustic solution in the far wake. We can say that this is a difficult complex task. The methods of mathematical modeling in aeroacoustics are described in detail in [14]. One thing is obvious, that in order to obtain an accurate spectrum of acoustic oscillations, the solution of a nonstationary gas-dynamic problem must be as complete, correct and adequate as possible.
Finally, the third approach to the simulation of high-speed viscous flows described by the Navier-Stokes equations is associated with the direct integration of these equations without involving additional models and empirical constants. The method of direct numerical simulation (DNS) is coupled with huge computational efforts. However, the development of computational aerodynamics methods, along with an increase in computer technology performance, made this approach quite viable. Of course, the range of problems solved by the DNS method is very limited. According to forecasts, the widespread practical use of DNS in solving complex problems of aerodynamics may begin only at the end of this century [15]. The research activities of the authors of this paper involve using the DNS method and carrying out numerical experiments based on this approach. On the basis of DNS, we solved some problems of aerodynamics (the results are presented below) and compared the data of numerical and natural experiments.
The purpose of this work is to determine the possibilities of direct numerical simulation at the present stage of computer technology development and computational algorithms using the results of solving several problems of aerodynamics, and to show the advantages and disadvantages of the method of integrating the Navier-Stokes equations by direct methods. It should be noted that by direct numerical simulation we mean the integration of the Navier-Stokes equations without the use of empiricism and additional models.
2. 3D-modeling of subsonic gas flows in a rectangular channel with sudden expansion
2.1. Simulation of transition flows
The laminar-turbulent transition has always been an important problem of hydromechanics, because the operation of a device depends on the flow regime. A feature of transitional flows is that pulsating zones of unsteadiness in a laminar flow appear abruptly. For such flows, it is difficult to choose a turbulence model, since in some area the flow is three-dimensional and unsteady, and somewhere it continues to be laminar. These specificities are mainly characteristic of flows with moderate Reynolds numbers of the order 103.
Let's write the classical system of three-dimensional Navier-Stokes equations in Cartesian coordinates in a dimensionless form:
dW OA <9B dC dt dx dy dz
0.
(2.1)
Here
A
B
C
/
W = (p, pux, pUy, puz, pE) pux
puxux + ~~ LITT
T
\
pUX Uy nxy
pux uz — nxz
\UX (pE + fcipO Qx ~ ~~ Uy^xy ~ UZ^XZ J
(
puy
puy ux 11
\
yx
pUyUy + ¿2 - n yy
puyuz nyz
\Uy (pE + + Qy — UxHyx — UyTlyy — Uzn,yzJ
(
puz pu z u^x n^z puzuy — n
\
puzllz +
zy
-n.
\uz (pE fcM^) Qz ux^zx uy^zy uz^zz)
(2.2)
(2.3)
(2.4)
(2.5)
The viscous stress tensor elements have the form
n
1 (duy dux
xy
n
yx
Re V dx
+
dy
TT 1 ( duz dux zx Re V dx dz
(2.6)
(2.7)
(2.8)
1 ( duy 2
n -n - — (^ + ^ ] (210)
The heat flux vector components are calculated as
1 dT ,
= "77—T^FTO (2-12)
(k - 1)Pr M2 dx' 1 dT
% ~ ~ (k — l)Pr M2 ~dy' (2'13)
, =__I_ÉL (2 M)
qz (k — l)Pr M2 dz ' '
Divergence is calculated using the formula
dux duy duz /„-..x
divu = —^ + —^ + —2.15 dx dy dz
The system (2.1) is closed by the Mendeleev-Clapeyron equation of state
- = RT. (2.16)
P
Thus, it can be seen that the system is based on the classical fundamental conservation laws and the equation of state.
The Reynolds (Re), Mach (M) and Prandtl (Pr) numbers are dimensionless complexes, on which the type of the flow and its characteristics depend.
The problem of subsonic (M = 0.6) viscous gas flowing through a channel with a sudden expansion at the inlet was solved. The gas temperature was assumed to be nearly constant and equal to T = 293 K. The speed of sound was assumed to be c* = 343 m/s. The scale for velocities was u* = 0.6 • c* = 205.8 m/s. For linear values, the gap height h = 1 mm served as the scale.
The channel geometry is shown in Fig. 1. The gap was one third of the channel height. The side walls were located at a distance of 4 units from the central plane Z = 0.
The numerical experiment began with a series of simulations with a low Reynolds number corresponding to laminar flows. The flow at the initial time leaned upwards; as a result, the final steady-state flow pattern turned out to be asymmetric, but stationary (Fig. 2). The figure shows streamlines and flow velocity vectors in the central plane. The upper separation zone is much smaller than the lower one. The characteristic Reynolds number was calculated from the gap height.
As the Reynolds number increased, the flow ceased to settle to a steady state, and the large lower separation zone began to break up into several fragmenting and merging vortices, periodically shearing downstream. The instantaneous velocity vectors in different planes for Re = = 400 are shown in Figs. 3. Figure 3a shows that the lower zone is destroyed, while the upper one retains its configuration. Several interacting vortices are observed in the plane X = 10, the flow has a complex structure. However, under the action of physical viscosity, the flow rapidly
Y
X
3 h
h\
Fig. 1. Channel with sudden expansion
Fig. 2. Steady-state laminar flow, Re = 120
y
(b) (c)
Fig. 3. Flow velocity vectors at Re = 400 in planes: (a) XY, Z = 0, (b) YZ, X = 10, (c) YZ, X = 25
relaminarizes, as evidenced by the data in Fig. 3c (vortex structures in the transverse plane at the end of the channel are almost absent). Thus, due to the stratified geometry, the flow becomes unsteady and three-dimensional in some segments, however, the vortex pulsations quickly decay.
Direct numerical simulation operates laminar and transient flows well at moderate Reynolds numbers and simple volumetric geometry. The flows shown above were calculated on grids of about 1 million nodes.
2.2. Simulation of subsonic developed turbulent flows
It is well known that, with an increase in the Reynolds number, nonstationary pulsating processes in the flow only intensify. Figure 4 shows a developed turbulent flow with the Reynolds
number Re = 0.5 • 104 in a three-dimensional channel with a sudden expansion at the inlet. Here, the velocity equal to 300 m/s was taken as the scale. The gap height was 1.5 mm.
Fig. 4. Flow velocity vectors in a rectangular channel with expansion at Re = 0.5 • 104 in the central plane
To minimize the numerical viscosity (so that it does not suppress the physical one), in addition to a powerful grid, it is necessary to use high-order approximation algorithms for derivatives.
A numerical experiment was carried out, consisting in the following: on computational grids of different thicknesses with different orders of accuracy, the flow in a rectangular channel 20x3x1 was calculated and the maximum value of the vortex vector modulus was determined. The grids considered are shown in Table 1.
Table 1. The grids under consideration
Item Number of nodes in x, y, z Grid power (mln. nodes)
1 400 x 45 x 15 0.27
2 600 x 63 x 21 0.8
3 800 x 93 x 31 2.3
4 1600 x 183 x 61 17.8
5 1920 x 273 x 91 47.7
6 2400 x 363 x 111 96.7
The calculation of spatial partial derivatives was carried out with the second, fourth, eighth and tenth orders of accuracy based on central difference schemes. In time, the equations integrated with the third order in all cases. Figure 5 shows graphs of changes in the maximum values of the vortex vector modulus for different grids and different approximation orders.
It can be seen that only the calculation with the eighth and tenth approximation orders on the 5th and 6th grids give the same results. In addition, the coordinate of the maximum value lies in the vicinity of the same point. It can be concluded that the calculation of flows with 104 order Reynolds numbers must be carried out on grids with a dimensionless step of at least 0.01 with the eighth order of accuracy. Nevertheless, even in this case, to estimate the near-wall parameters, additional verification of the results obtained will be required, depending on the problem statement and the parameter under study. However, it is quite possible to estimate the acoustic characteristics of the flow. Of course, the calculation on a grid of 100 million nodes of complex nonlinear equations takes considerable time even on a supercomputer with the possibility of parallel computing. But the results obtained have an undeniable value, which consists in the fact that they give a complete picture of the flow with the possibility of finding any spectral characteristics.
L s. ^ V ------- ' 1
/ / / /•V >
/ S
/ / / -- 2nd — 4th
1 / / < 8th -- 10th -1—
grid number
Fig. 5. The maximum values of the vortex vector modulus
3. Axisymmetric simulation of supersonic flow behind a cylinder
When it comes to 2D-modeling in relation to the DNS, many critics of the approach immediately declare that this is impossible, since turbulence is purely three-dimensional. We do not argue with this, turbulence, as a physical phenomenon, is a three-dimensional process. Indeed, in the axisymmetric formulation of the problem, turbulence modeling does not occur. We repeat once again that by "direct numerical simulation" we mean modeling based on the solution of the Navier-Stokes equations, without involving any empiricism, and not only calculations of turbulent flows, in which the movements of all scales present in the flow are resolved. Speaking about direct numerical modeling based on the Navier-Stokes equations, we, in fact, turn to the study of solutions to these equations. It is obvious that in practice the axisymmetric flow is not realizable, but nothing prevents the study of the equations (this is the advantage of numerical simulation). We believe that such a formulation is not without meaning, since we can say that, despite the two-dimensional calculation, the axisymmetric formulation does not violate the notions of the three-dimensionality of space.
Based on the axisymmetric model of the Navier-Stokes equations, we calculated the flow behind a circular cylinder. These studies are important when it comes to investigating and evaluating base drag, which is up to half of the total aerodynamic drag of an object moving at high speed and arising due to the formation of a low pressure area behind the object.
The dimensionless Navier-Stokes equations in an axisymmetric formulation in a cylindrical coordinate system have the form
dW dA dB D
dt dx dr r
0.
Here
W = (p, pux, pur, pE)
T
(3.1)
(3.2)
A
pux
pUxUx + ¿.jyj PU x Ur - n
-nx
\UX [pE + ¿.jyp ) + Qx "UTIITT UrTir
(3.3)
\
B
pur
pUrUr I I ^ r r
\ur (ypE ) Qr uxTixr urUrrJ
(
D
pur
pur Ur - nrr + ne
\
\ur (ypE kjyp ) Qr uxTixr urTirr J
j-, _ ux + ur
n.
2 1
k(k - 1)M2 p
Re V^ dx 3
nr
1 ( dur dux
Re V dx
+
dr
1
Re V ^ dr 3 U
Ildd = — fe--divu
dd Re V r 3
1
Qx = -
dT
(k - 1)Pr M2 dx '
1 dT
(k - l)Pr M2 ~dr'
du^™ du^™ u« divu= —- + - + —.
dx dr r
(3.4)
(3.5)
(3.6)
(3.7)
(3.8)
(3.9)
(3.10)
(3.11)
(3.12)
(3.13)
One more point should be noted: in contrast to the Cartesian planar formulation, the 2D axisymmetric formulation allows one to resolve the azimuthal stress (Eq. (3.10)).
The shape of the ogival body under consideration was close to the shape of a projectile. The computational domain (Fig. 6) was covered with a powerful rectangular spatial grid with mesh thickening near the streamlined walls. The grid capacity was 26 million nodes.
r 4
y
Computational domain
0 6=1 2 3 4 5
Fig. 6. Computational domain behind the streamlined body
1
p
The supersonic flow behind a circular cylinder was simulated with the Reynolds number Re = = 1.65 • 106 and the Mach number M = 2.46. An environment with a pressure of one atmosphere and a temperature of 293 K was considered. The scale for the velocity was u* = 844 m/s. The scale for linear quantities was the radius of the cylinder, equal to 3.25 cm.
The density shaded map is shown in Fig. 7. It can be seen that the flow consists of a set of interacting shock waves.
J_I_I_L
-1-1-1-1-1-1-1-
1 2 3 4 5 6 7 x 8
Fig. 7. Shaded density map behind a streamlined cylinder, Re = 1.65 • 106, M = 2.46
To verify the grid, the maximum and minimum values of vorticity were estimated (Fig. 8a), and the distance from the wall to the first computational node along the bottom was calculated in near-wall coordinates (Fig. 8b).
-80
>> -100
V • p-H
tH g -120
S -140
• p-H Ö
s -160
-180
10 20 30 40 50 60 0.0 0.2 0.4 0.6
number of nodes, millions
(a) (b)
Fig. 8. Convergence validation: (a) vorticity, (b) dimensionless distance from the first node to the wall
It was established [16] that the viscous sublayer is located near the surface at a distance of y+ < 5. For turbulent flows, the value y+ ~ 2 is considered to ensure proper resolution of
the viscous sublayer [17]. The data of Fig. 8 show that convergence in mesh grinding and the conditions for resolution of the viscous sublayer are satisfied.
To determine the base drag, it is necessary to have a distribution of the values of the base pressure coefficient along the surface, which is calculated by the formula
Cp —
2(^-1
\Po_
kM2
(3.14)
where pb is the pressure on the wall and p0 is the pressure in the environment.
Figure 9 compares the results of DNS modeling with experiment data and other models results. As expected, the DES method gave the best results in comparison with the experiment. It can be seen that the DNS modeling results are not inferior in accuracy to the RANS and LES methods results. This is especially noticeable if we evaluate the overall base drag coefficient (given in Table 2), defined as
Ch — -
Cp(r) dS.
(3.15)
0.00
CL
-0.05
-0.10
-0.15
• Experiment
-DNS
--RANS (k-e)
--DES
........LES
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig. 9. Distribution of the base pressure coefficient: experiment [18], DNS (the present work), RANS, DES, LES [19]
b
S
Table 2. Integrated base drag coefficient
Experiment [18] DNS RANS DES LES
0.103 100 0.117 0.092 0.137
As expected, the DNS method allows the largest deviations near the axis of symmetry. Refining the grid near the axis partially solves the problem: the base drag coefficient will be calculated more accurately. Especially considering that the closer the point is to the axis, the smaller the contribution to the overall integral. However, such parameters as the point of attachment of the vortex behind the streamlined body or the length of the separation zone will be simulated with some error. Although the parameter necessary for practice is calculated quite adequately, accurately and acceptably. This, in our opinion, is the indisputable value of these results, since they can serve as a basis for proper fitting of the coefficients and functions of turbulence models.
4. 2D simulation of the buffet onset
The buffet phenomenon became known in the 1930s. Buffet is associated with the occurrence of self-oscillations of aircraft structural elements [20]. At first, buffet was considered exclusively as a phenomenon that occurs as a result of an object falling into a wake, or into a turbulence zone. But, besides buffet, in the form of "shaking" due to atmospheric turbulence, there is also speed buffeting. At certain parameters, the wing begins to flow around the transonic flow, which leads to the occurrence of a wave crisis. Due to periodically occurring shock waves and flow separations, vibrations appear. The undesirable consequences of high-speed buffet are associated not only with the problems of aircraft controllability and crew well-being due to shaking, but also with the fatigue strength of the aircraft elements. High-speed buffet is not caused from outside, but is determined by the characteristic features of the flow around the aircraft itself. This type of buffet is extremely dangerous and difficult to reproduce on physical models.
Obviously, at present it is impossible to calculate the flow around an aircraft or even its separate part using DNS. Modern computing power does not allow this. However, it is quite possible to formulate some problems in a flat setting, when one of the linear dimensions of the object under study significantly prevails over the others. For example, consider the profile of a wing. Of course, despite the fact that the flow characteristics in the profile plane can be an order of magnitude greater than the characteristics in other planes, it is categorically impossible to speak of a reliable physical correspondence in 2D modeling. Such a mathematical calculation can help to qualitatively study the process, to identify trends that will really be relevant in practice. In this regard, solving the problem in 2D makes sense, especially considering that the 3D version is impossible.
So, the flow around the NACA0012 airfoil was simulated (Fig. 10). The system of equations (2.1) was written in a two-dimensional formulation by eliminating the z coordinate. The computational grid with mesh thickening near the streamlined object consisted of about 17 million nodes. Here, the scale for the velocity was u* = 240 m/s, the length of the profile l* = 18.2 cm served as the scale for linear quantities.
The results of numerical simulation of the flow around the NACA0012 airfoil with characteristic numbers M = 0.734, Re = 1.0 • 106, on backwind side were compared with the experimental data of [21]. The pressure coefficients were compared (Fig. 11). The simulation results are in good agreement with experiment. On the backwind side, the discrepancies are more significant, which is explained by the presence of flow separations. But taking into account that the experi-
y
x
Fig. 10. NACA0012 Airfoil
ment reproduces a three-dimensional flow, and the simulation is two-dimensional, the results can be considered quite satisfactory.
The authors of [22] showed that a lambda-shaped shock wave is formed on the backwind side. Behind the leg of the lambda-shaped shock, the flow is separated from the profile surface. But further, the separated boundary layer, together with the main flow, is blocked by the main shock wave, behind which a zone of low pressure and low velocities is formed. All this leads to significant pressure fluctuations on the backwind side of the airfoil (especially in the area of its middle).
Also in [22], the changes in the coefficients over time are analyzed. At an angle of attack of 5 degrees, the fluctuations in the coefficients are insignificant. As a increases, the oscillation amplitude increases. Starting from an angle of attack of 6 degrees, the level at which fluctuations of the coefficients occur changes. At 7 and 8 degrees, this is especially pronounced. Moreover, the fluctuations of the level, relative to which the pulsations occur, are quite regular. The periods of increase and decrease of the levels, which have the same length in time, are clearly distinguished, and the amplitude of the fluctuations of the levels themselves is determined. We can talk about the presence of self-oscillations of levels, at which time pulsations of the coefficient Cp occur, with much lower frequencies than the local pulsations themselves.
Based on the frequency analysis, it was determined that at 8 degrees angle of attack, the spectral density of the coefficient CL is concentrated at a frequency of 12.9 Hz. In [23], when solving a similar problem with similar parameters using the RANS method with the SST (Shear Stress Transport) turbulence model, the self-oscillation frequency of 14.3 Hz was obtained. It can be said that the results are consistent.
Figure 12 shows the pulsation energy generation (Pk) surfaces for various angles of attack. Pk is determined from the expression
Pk = ~ « ' <) ^ " « ' <> (|f + " « ' <> |f • (4-1)
The data of Fig. 12 indicate that at the 4 degrees angle of attack unsteady processes on the upper surface are intensified, and the turbulence parameters are nonzero. The gradient maps visualize a shock wave followed by velocity pulsations along the airfoil surface. At an angle of attack equal to 6 degrees, the zone of high generation of pulsation energy increases even more and begins already before the shock wave. It can be said that at an angle of 8 degrees, nonstationary oscillatory processes occur along the entire backwind surface, both before and after the shock
(c)
Fig. 12. Velocity pulsation energy generation surfaces: (a) a = 4°, (b) a = 6°, (c) a = 8°
wave. The shock wave only intensifies the oscillations, forming a zone of increased values of pulsation characteristics. Thus, we can conclude that, at these angles of attack, high-speed buffet arises, which consists in a high intensity of oscillatory processes.
Thus, the flat formulation of the problem and the solution of the two-dimensional Navier-Stokes equations by the DNS method are not without meaning. Of course, one cannot fully rely on 2D modeling data, especially quantitative data. However, when solving mathematical
equations, it is quite possible to identify physical trends, obtain an indirect confirmation of hypotheses and to identify the direction of further research.
5. On boundary conditions, computational schemes, and parallel processing
Definitely, the essential parts of mathematical modeling are the correct and adequate setting of boundary and initial conditions, the choice of efficient computational algorithms with the possibility of implementing these algorithms on the available computing power. Since this article is of an overview and analytical type, we give no exact mathematical formulations of certain conditions or schemes used by us, however, in this section we will briefly indicate the sources and essence of the methods.
Laminar and transient subsonic flows, as well as turbulent subsonic flows, were calculated by us on the basis of high-order central-difference spatial schemes [24]. For the time integration of such flows, the method based on the Taylor formula has proven itself well, which we used. The absence of discontinuities in the solution and differentiability in the entire computational domain considered make, in this case, high-order approximation central-difference schemes advantageous from the point of view of saving computer time.
As boundary conditions in solving the problems of simulation the internal subsonic flows we used the ideas formulated in the works of A. T. Fedorchenko [25], K.Tompson [26] and H. Kreiss [27] for accurate derivation of disturbances.
When calculating the external transonic or supersonic flow around bodies, high-accuracy algorithms are needed that are capable of reproducing shock waves. Therefore, to calculate spatial partial derivatives, we used WENO (Weighted Essentially Non-Oscillatory) schemes based on adaptive templates and automatic smoothness analysis [28]. Integration over time was carried out using the TDV (Total Variation Diminishing) Runge-Kutta scheme [29]. Thus, if the solution in the area considered is smooth, then it was approximated with a high order on a wide stencil. In the presence of discontinuities or high gradients, the template narrowed down to two points.
Nonreflecting boundary conditions for open boundaries were formulated by L. V. Dorodnit-syn [30]. These conditions were used by us when calculating the base drag and the flow around the airfoil.
In all problems, the gas was assumed to be viscous and the walls to be solid; therefore, sticking and impermeability conditions were set on solid surfaces. Additional boundary conditions for the boundary layer near the walls were not specified.
The DNS method requires huge machine costs, so its implementation demands a multiprocessor system with the ability to parallelize the computing process. Even the most modern personal computer is unlikely to handle the DNS. The calculation results presented in this paper were obtained on the cluster of the Interdepartmental Supercomputing Center of the Russian Academy of Sciences and the Uran supercomputer (N. N. Krasovsky Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences).
6. Conclusion
The results of solving problems presented in the paper, obtained on the basis of direct numerical simulation of the Navier-Stokes equations, show the consistency of the method, despite
the skepticism expressed by many scientists in hydromechanics regarding it. Undoubtedly, the advantage in speed and the ability to solve real practical problems make methods based on turbulence models more attractive. However, such methods limit, in a sense, the possibilities of computer simulation as an independent research tool, since no one can vouch for the reliability of the developed technique in the case of flow conditions different from those under which the model was calibrated. In this regard, direct numerical simulation comes in handy, since it is free from any empiricism and can serve as a tool for calibrating models. The comparisons of simulation results shown in the paper testify to this. In addition, there is a class of flows that the DNS method handles better than others. Obtaining a full range of information about the flow during DNS modeling can also make it possible to obtain its acoustic characteristics.
Conflict of interest
The authors declare that they have no conflicts of interest.
References
[1] Reynolds, O., An Experimental Investigation of the Circumstances Which Determine whether the Motion of Water Shall Be Direct or Sinuous, and of the Law of Resistance in Parallel Channels, Philos. Trans. R. Soc. Lond, 1883, vol. 174, pp. 935-982.
[2] Reynolds, O., On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion, Philos. Trans. R. Soc. Lond., 1895, vol. 186, pp. 123-161.
[3] Bradshaw, P., An Introduction to Turbulence and Its Measurement, Oxford: Pergamon, 1971.
[4] Kochin, N.E., Kibel, I. A., and Roze, N. V., Theoretical Hydrodynamics, New York: Wiley, 1964.
[5] Monin, A. S. and Yaglom, A.M. Statistical Fluid Mechanics: Mechanics of Turbulence: Vol. 1, Cambridge, Mass.: MIT Press, 1987.
[6] Landau, L.D. and Lifshitz, E. M., Course of Theoretical Physics: In 10 Vols.: Vol. 6. Fluid Mechanics, 2nd ed., Oxford: Butterworth-Heinemann, 2003.
[7] Garbaruk, A., Strelets, M. Kh., and Shur, M. L., Turbulence Modeling in Complex Flow Calculations, St. Petersburg: SPbPU, 2012 (Russian).
[8] Spalart, P.R., RANS Modeling into a Second Century, Int. J. Comut. Fluid Dyn, 2009, vol. 23, no. 4, pp. 291-293.
[9] Xiao, H. and Cinnella, P., Quantification of Model Uncertainty in RANS Simulations: A Review, Prog. Aerosp. Sci, 2019, vol. 108, pp. 1-31.
[10] Gorman, J., Bhattacharyya, S., Cheng, L., and Abraham, J., Turbulence Models Commonly Used in CFD, in Computational Fluid Dynamics, S. Bhattacharyya (Ed.), London: IntechOpen, 2021, 19 pp.
[11] Berselli, L.C., Galdi, G.P., Iliescu, T., and Layton, W. J., Mathematical Analysis for the Rational Large Eddy Simulation Model, Math. Models Methods Appl. Sci., 2002, vol. 12, no. 8, pp. 1131-1152.
[12] Spalart, P. R., Jou, W. H., Strelets, M. Kh., and Allmaras, S. R., Comments on the Feasibility of LES for Wings, and on a Hybrid RANS/LES Approach, in Advances in DNS/LES: Direct Numerical Simulation and Large Eddy Simulation: Proc. of the 1st AFOSR Internat. Conf. on DND/LES (Ruston, La, 1997), pp. 137-148.
[13] Spalart, P. R. Hybrid RANS-LES Methods, in Advanced Approaches in Turbulence: Theory, Modeling, Simulation, and Data Analysis for Turbulent Flows, P. Durbin (Ed.), Amsterdam: Elsevier, 2021, pp. 133-159.
[14] Acoustic Interactions in Gas Flows, V. N. Emel'yanov, K. N. Volkov (Eds.), Moscow: Fizmatlit, 2021 (Russian).
[15] Spalart, P. R., Strategies for Turbulence Modeling and Simulations, Int. J. Heat Fluid Flow, 2000, vol. 21, no. 3, pp. 252-263.
[16] Pope, S.B., Turbulent Flows, Cambridge: Cambridge Univ. Press, 2000.
[17] Menter, F.R., Kuntz, M., and Langtry, R., Ten Years of Industrial Experience with the SST Turbulence Model, in Turbulence, Heat and Mass Transfer 4: Proc. of the 4th Internat. Symp. on Turbulence, Heat and Mass Transfer (Antalya, Turkey, Oct 2003), K. Hanjalic, Y. Nagano, M. J. Tummers (Eds.), New York: Begell, 2003, pp. 625-632.
[18] Herrin, J. L. and Dutton, J. C., Supersonic Base Flow Experiments in the Near Wake of a Cylindrical Afterbody, AIAA J., 1994, vol. 32, no. 1, pp. 77-83.
[19] Simon, F., Deck, S., and Guillen, P., Reynolds-Averaged Navier-Stokes/Large-Eddy Simulations of Supersonic Base Flow, AIAA J., 2006, vol. 44, no. 11, pp. 2578-2590.
[20] Abdrashitov, G.G., On the Issue of Buffeting the Tail Unit, Moscow: TsAGI, 1939 (Russian). See also: Abdrashitov, G., Tail Buffeting, NACA-TM-1041 (Feb 1943).
[21] McDevitt, J.B. and Okuno, A.F., Static and Dynamic Pressure Measurements on a NACA 0012 Airfoil in the Ames High Reynolds Number Facility, NASA-TP-2485 (1985).
[22] Lipanov, A. M., Karskanov, S. A., Chernyshev, S. L., and Lipatov, I. I., Theoretical Investigation of Conditions for the Appearance of High-Speed Bufting, Vestn. Udmurtsk. Univ. Mat. Mekh. Komp. Nauki, 2019, vol. 29, no. 3, pp. 382-395 (Russian).
[23] Lipatov, 1.1., Fam, T.V., and Prikhod'ko, A. A., Numerical Simulations of Buffeting Appearence, Trudy Mosk. Fiz.-Tekhn. Inst, 2014, vol. 6, no. 2, pp. 122-132 (Russian).
[24] Lipanov, A.M., Theoretical Fluid Mechanics of Newtonian Media, Moscow: Nauka, 2011 (Russian).
[25] Fedorchenko, A. T., Numerical Investigation of Unsteady Subsonic Viscous Gas Flows in a Plane Channel with a Sudden Expansion, Fluid Dyn., 1988, vol. 23, no. 4, pp. 509-517; see also: Izv. Akad. Nauk SSSR Mekh. Zidk. Gaza, 1988, no. 4, pp. 32-41.
[26] Thompson, K.W., Time-Dependent Boundary Conditions for Hyperbolic Systems: 2, J. Comput. Phys., 1990, vol. 89, no. 2, pp. 439-461.
[27] Kreiss, H.-O. and Lorenz, J., Initial-Boundary Value Problems and the Navier - Stokes Equations, Pure Appl. Math., vol. 136, Boston, Mass.: Acad. Press, 1989.
[28] Liu, X.-D., Osher, S., and Chan, T., Weighted Essentially Non-Oscillatory Schemes, J. Comput. Phys, 1994, vol. 115, no. 1, pp. 200-212.
[29] Gottlieb, S. and Shu, C.-W., Total Variation Diminishing Runge-Kutta Schemes, Math. Comp., 1998, vol. 67, no. 221, pp. 73-85.
[30] Dorodnicyn, L. V., Nonreflecting Boundary Conditions and Numerical Simulation of External Flows, Comput. Math. Math. Phys., 2011, vol. 51, no. 1, pp. 143-159; see also: Zh. Vychisl. Mat. Mat. Fiz, 2011, vol. 51, no. 1, pp. 152-169.
