UDC 533.661.2
DOI: 10.18698/0236-3941-2019-6-4-16
USING THE MAXIMUM PRESSURE PRINCIPLE FOR VERIFICATION OF CALCULATION OF STATIONARY SUBSONIC FLOW
V.A. Anikin 1 V.V. Vyshinsky 2 O.A. Pashkov 1 E.V. Streltsov 2
[email protected] [email protected] [email protected] [email protected]
1 LCC "VR-Technologies", Moscow, Russian Federation
2 Zhukovsky Central Aerohydrodynamic Institute TsAGI, Zhukovsky, Moscow Region, Russian Federation
Abstract
The principle of maximum pressure for subsonic stationary three-dimensional vortex flows of an ideal gas (author Sizykh G.B., 2018) is applied to verify the calculation method and its implementation on a specific computer technology. The four criteria for solution's verification are proposed. The method for obtaining flow parameters is based on solving of discrete analogs of the Navier — Stokes system of equations on three-dimensional non-structured computational meshes. For example, there was consider the vortex tear-off flow around the fuselage of a helicopter with an empennage and landing gear at obviously insufficient computing resources. Conclusions of the feasibility of applying the author's criteria for evaluation of a particular calculation and for estimation of reliability of the results have been made
Keywords
Nonlinear partial differential equations, boundary value problems, Euler equations, Reynolds averaged Navier — Stokes equations, subsonic vortex flows, subsonic maximum pressure principle
Received 31.10.2018 © Author(s), 2019
Introduction. The methods of computational aerohydromechanics within the framework of grid approaches realize the solution of initial-boundary value problems for nonlinear partial differential equations on a specific computing technique. The correctness of the boundary value problems being solved (with the exception of the simplest ones) has no proof, and the convergence of the numerical solution to the solution of the original boundary value problem is proved only in the linear case (Lax theorem). In most cases, the available computing power is not enough for full verification, that is, evidence that the solution obtained depends only on the physical parameters of the problem and does not depend on the parameters of the numerical scheme (dimensions of the
computational domain for external aerodynamics problems, mesh condensation, number of iterations). Therefore, any additional way of checking the adequacy of the calculation results to the solution of the initial boundary value problem, as well as the possibility of identifying weaknesses in its numerical implementation, is of value.
Verification methods can be conditionally divided into two types. The first is a comparison of the obtained solution with a known exact solution or with a solution obtained by another, well-established numerical method. For example, the calculation method within the framework of the Neumann boundary value problem for the complete equation with respect to the velocity potential can be verified using exact solutions for the jumpless flow around the aerodynamic profile with a local supersonic zone [1]. The second type is verification of the implementation of known patterns. For example, when solving problems in the potential approximation, it is possible to verify by checking the condition that the velocity circulation is equal to zero along the contours of all cells of the computational grid and the mass flow through their boundaries; in the case of adiabatic flows, the entropy is constant for each computational node.
In this paper, as an additional verification of numerical methods for calculating the flow of a viscous gas, we use the consequences of the maximum pressure principle [2]. In the most complex viscous gas flows, there are zones in which a simpler boundary-value problem can be solved, for example, laminar sections in a turbulent flow, flow areas where viscosity or compressibility can be neglected, and stationary zones. If there are zones in the flow field that are simulated with high accuracy in the framework of the boundary value problem for Euler equations, then the consequences of these equations must also be fulfilled.
The principle of maximum pressure and criteria for independent verification. The first theorems on extreme pressure values resulting from the equations of fluid motion include Rowland's theorems [3] for an ideal income-pressible fluid. Poincare [4] for an incompressible fluid in a gravitational field created by the fluid itself, and Hamel [5] for viscous incompressible fluid [6]. In these theorems, statements are made about the minimum or maximum pressure depending on the sign of the parameter Q, the expression for which includes only the first spatial derivatives of the velocity components, which is convenient for verifying numerical solutions because it reduces the requirements on the order of approximation accuracy and on the fineness of the grid step in determining sign of the parameter Q. These theorems relate to an income-pressible fluid and are not applicable to gas.
In [7], the maximum principle was obtained for the quantity J dp jp in barotropic gas flows, the calculation of which uses the values of the second derivatives of the velocity components. Compared with the conditions of Hamel's theorem [4] (for an incompressible fluid), these quantities contain the "additional" term d I/dt + vWI, where I = divv, which contains the second spatial derivatives of the velocity components. In addition, the theorems [7] were proved only for a barotropic gas p = f (P), that substantially narrows the class of flows under consideration.
In this paper [2], the principle of maximum pressure in stationary subsonic vortex flows of an ideal gas is obtained. The conditions of this maximum principle include only the parameter Q (more precisely, only its sign). In this case, the flow can be non-barotropic, i.e., the entropy function, remaining constant along the streamline, canbe different on different streamlines. Accordingly, energy along the streamline is conserved. The maximum principle [2] differs from the maximum principle obtained by Schiffman [8] in that the Schiffman principle does not apply to pressure but to velocity and is proved only for barotropic and vortex-free flows [9], while that obtained in [2] the principle of maximum pressure can also be used to verify the calculations of non-barotropic vortex gas flows.
The interfaces of many software systems (see, for example, [10-13]) contain the possibility of constructing parameter level surfaces Q:
Q = 0.5 (Qy : Qj - Sy : Sy ). Where the sign : means the convolution of tensors,
Qij = 0.5{duijdxj -duyIdxi) and Sy = 0.5{du^jdxy +duyldxi)
are antisymmetric vorticity tensor and symmetric strain rate tensor.
We consider the stationary subsonic flow of an ideal gas in a certain bounded closed region G, which is described by the Euler equations written in the form [4]:
Ox v = - - Vp - 1W2 + g;
P 2
Q = rot v; div (pv) = 0,
where p is density; p is pressure; v is velocity; V = |v| is module of velocity; g is gravitational acceleration. The pressure p and the density p are related by the relation p = opk, where k is the adiabatic exponent; a is the entropy function.
It can be different on different streamlines, but is maintained along each specific streamline
v-V(p"kp ) = 0.
The flow is assumed to be subsonic, i.e., it is believed that the velocity V at each point in the field G is less than the local speed of sound -yjkp/p :
M = v/^kp/p < 1.
In the field G all hydrodynamic parameters are assumed to be twice continuously differentiable coordinate functions.
It was proved in [2] that under these conditions the following principle of pressure maximum is satisfied: if at all points of the domain G, including boundaries Q < 0, then at the internal point G the pressure cannot have a local minimum; if Q > 0, then at the internal point G the pressure cannot have a local maximum. It is not only about the impossibility of achieving strict extrema, but also the impossibility of achieving non-stringent extrema at the internal point G.
The presence in the expression for Q of only the first derivatives of the velocity components makes it possible to determine the sign of Q even for numerical solutions obtained by low-order approximation methods. Moreover, the fact of the presence of a local minimum or a local maximum of pressure in a certain region of the flow is determined by the usual "enumeration" of values at the nodes of the computational grid.
Based on the principle of maximum pressure, G.B. Sizykh for the verification of numerical methods of a wide class and their specific implementations, the following criteria are proposed.
Criterion 1. The positivity of Q in the nuclei of vortices. Since there is a rarefaction in the core of the vortex, a local minimum of pressure (otherwise the vortex would not exist), then Q < 0 indicates a violation of the principle of maximum pressure.
Criterion 2. The positivity of Q in the regions of flow separation. Since the separation region is characterized by low pressure (source of bottom resistance), Q < 0 indicates a violation of the principle of maximum pressure. Therefore, it should be checked that the condition Q > 0 is satisfied inside these zones.
Criterion 3. The absence of sharp oscillations of Q in the flow region. The presence of oscillations indicates a sharp anisotropy of the flow and, as a consequence, incorrect linear models of gradient diffusion, which are the basis for the usual closure models. In addition, the introduction of coherent structures into equilibrium turbulence violates the ergodicity hypothesis, on which the
averaging of the Navier — Stokes equations by Reynolds is based. Sharp oscillations indicate an excess of the measure of this violation in a specific implementation of the numerical method.
In zones far from the "sources" of vorticity, the flow can be vortex, but the viscosity of the gas can be neglected and its motion described by Euler equations. In such zones, the maximum pressure principle obtained in [2] is valid and can be used to verify calculations of viscous gas flows. Such verification was carried out during the numerical simulation of the flow around the layout of the helicopter fuselage + plumage + landing gear without a rotor.
The solution of the rotor translational motion problem requires additional justification, since it is solved in a non-inertial coordinate system that rotates with the rotor, which deprives the ergodicity system on which the averaging of the Navier — Stokes Reynolds equations (RANS) is based. In addition, in this coordinate system there is no isotropy, on which linear gradient diffusion models are based.
Calculation method. The boundary-value problem for RANS is solved with a two-parameter SST closure model combining k-& (in the flow field) and ks (near the body) of the model. In models of this class, the closure of the governing equations is carried out using two transport equations. The "exact" expression for k and the "semi-empirical" differential equation for the turbulence dissipation rate s are used. The quantity s plays a decisive role in the structure of turbulence, since it is associated with the flow of energy of turbulent pulsations along the cascade of vortices (energy sink) and characterizes the frequency of turbulent pulsations (o = s / k). To find correlations ui uy (the bar means Reynolds averaging) second-order partial differential equations are solved. Since k and s in some cases numerical implementation can take on very small or too large values ("poor conditionally"), in k-s the model can arise the problems with accuracy and stability of calculations. In order to eliminate them, a class of k-& models used in the composite SST circuit model has been introduced. The hypothesis of isotropic linear gradient diffusion with all the disadvantages arising from this is used.
Reynolds averaging is based on a modern understanding of turbulence as a combination of weakly correlating vortex movements. Despite the intermittency and random distribution over a wide range of spatial and temporal scales, turbulence actually consists of local vortices that interact with each other during movement.
Part of the energy spectrum is resolved explicitly, while the other is taken into account using approximate modeling. When calculating the parameters of
turbulent flows using two-parameter turbulence models during interpolation, at the initial stages of the calculation, k and/or s and/or ra can take negative values, as a result of which the solution is corrected so that they remain positive. Strong nonlinear relationship between equations can lead to instability of the multigrid method.
The presence of regions containing many fragments Q < 0 in the zone of formation of vortex structures indicates a violation of local isotropy, which reduces the accuracy of closure models and indicates insufficient power of the computational grid.
Calculation results. The calculations of the fuselage + vertical tail (VT) + + horizontal tail (HT) + landing gear (Fig. 1) were performed using the
Fig. 1. The appearance of the layout of the helicopter and the position of the points Pmin and Pmax, starting with the hundredth iteration (the front edge of the landing gear and the rear edge of VT)
EWT-TsAGI computer code [14] using a second-order approximation scheme in space. A stationary solution is constructed using the establishment method. Detailed thickening in the area of the streamlined body location made it possible to describe in more detail the flow, especially the area of stern separation and the near wake, which significantly affect the field of bevels of the stream in the area of the streamlined body (Fig. 2).
The number of nodes of the computational grid (9 853 000) and the dimensions of the computational domain (distance to the boundaries of 15-20 characteristic sizes of the problem) are obviously insufficient to solve the problem of this level of complexity. This makes the additional method of verification and identification of weaknesses in the numerical implementation of
Fig. 2. Field of numbers M in the plane of symmetry (a = 8°)
the method particularly influential, affecting both the accuracy of determining the integral characteristics (primarily aerodynamic drag) and the modeling of individual flow fragments.
The calculations were performed at the free-stream Mach number M = 0.3, static pressure Pstat = 101 325 Pa and the temperature 288.15 K atmosphere. The characteristic area for dimensionlessness when calculating the integral characteristics of 1 m2.
In Fig. 3-5 show the results of the calculation of the flow around at zero values of the angles of attack (a = 0) and slip (P = 0). It is assumed that the symmetry of the boundary conditions corresponds to a symmetric flow field. In Fig. 3 shows half the surface of the fuselage and the surface Q = 0, visualizing vortices. Inside the vortices, the pressure is reduced. Criterion 1 is fulfilled in these zones only with a sufficient density of the computational grids and a sufficiently large size of the computational domain.
Fig. 3. The investigated layout and surfaces Q = 0, visualizing vortex structures
CP »10°
oja -0 2$ -ОВД
Fig. 4. Pressure coefficient Cp and streamlines in the plane of symmetry
Macft Number
I040
0 30
020
Fig. 5. The field of Mach numbers in the plane of symmetry
Starting with the hundredth iteration, the position of the minimum and maximum pressure points is stabilized (see Fig. 1). After the convergence of the solution Pmin - Pstat = -25 009.5 Pa, Pmin - Pstat = 7423.15 Pa. At all iterations, the points of pressure extremum lie on the boundaries of the flow (on the surface of the streamlined body), which corresponds to the maximum principle and speaks in favor of the correctness of the problem being solved and the implemented algorithm for its solution (Criterion 4).
The point of the maximum value of the characteristic Q is located on the front edge of the chassis Q max — 1.66 • 1012 s-2, the pressure at this point is APq max — Pqmax - Pstat = -9410.92 Pa. The point of the minimum value of the characteristic Q is located at the trailing edge of the horizontal tail Qmin = = -2.96 • 1011 s-2, APqmin = Pqmin - Pstat = -4107.52 Pa.
The pattern of streamlines in the plane of symmetry (Fig. 4) indicates the presence of a gap in the lower part of the aft fuselage region, where, by virtue of a local decrease in pressure, the criterion must be fulfilled 2.
In Fig. 6, a-c show surfaces Q = 0 at 1, 3, 5, 20, 1000 iterations of the solution. The table shows the minimum and maximum values of the parameter Q in the flow field for these solutions and the pressure values at these points. One can see the generation of vortices at the boundaries of the grid condensation, as well as vortices arising when the initial field is specified (Fig. 6, a).
с
Fig. 6. Surfaces Q = 0 at iterations 1 and 3 (a), 5 and 20 (b) and at the thousandth iteration (converged solution) (c)
Values of characteristics Q and pressure at iterations
Niter Qmax> s Qmim s APqmax, Pa APQmin, Pa
1 1.55E + 12 -1.72E + 11 91 192.5 91 192.5
3 5.95E + 12 -6.56E + 11 99 460.67 96 731.68
5 5.60E + 12 -5.44E + 11 98 490.7 99 641.43
20 1.97E + 12 -7.50E + 11 87 411.2 98 118.47
1000 1.64E + 12 - 88 664.4 93 990.66
Strong diffusion of vortices in the near wake (Fig. 6, b) indicates the insufficient size of the computational domain, which in particular reduces the accuracy of determining the integral characteristics. The converged solution (Fig. 6, c) demonstrates the destruction of the vortex wake at a distance of the order of three lengths of the streamlined body, which cannot but affect the accuracy of determining the forces and moments and indicates the need to increase the distance to the exit boundary and/or power of the grid.
In Fig. 7 shows the cross sections x = const of the surfaces Q = 0 for the last solution. The flat contours Q = 0 are filled with pressure fields (left) and fields Q > 0 (right). It can be seen that inside the vortex structures, where the pressure is lowered, criterion 1 is satisfied. The appearance of subdomains inside the vortex structures and in the area of fodder separation, where Q < 0 (Fig. 8), indicates a violation of criteria 1 and 2, which indicates insufficient accuracy of the solution in the fodder area, and due to the ellipticity of the problem, then in the entire flow field.
Pressure Iso Clip 10 Pressure Iso Clip 10
2875 00 —375 00
4500.00 1250.00 ' -2000.00 100.00 75.00 50.00 25.00 0.00
Fig. 7. The pressure fields (left) and the Q-characteristic (right) inside the contours formed by the intersection of the surfaces Q = 0 and the planes x = const
Fig. 8. The fields of the parameter Q inside the contours formed by the intersection of the surfaces Q = 0 and the planes x = const
Conclusion. Along with existing verification methods, it is proposed to use Criteria 1-3 based on the principle of maximum pressure proposed in [2] as additional criteria for the accuracy of methods for calculating stationary flows, according to which surfaces divide the flow region into subdomains where, in which there can be no local pressure maximum, and where there can be no local minimum.
Criterion 1. The positivity of Q in the nuclei of vortices. Criterion 2. The positivity of Q in the regions of flow separation. Criterion 3. The absence of sharp oscillations of Q in the flow region. Criterion 4 (optional). The points of pressure extremum lie on the boundaries of the flow. In this paper, the authors show the application of these criteria to assess the quality of the calculation (the adequacy of the size of the computational domain, the density of the grid, the number of iterations) of a particular implementation of the method using available computer technology.
In order to increase the accuracy of calculating the integral characteristics in this case, it is necessary to increase the size of the computational domain, first of all, behind the streamlined body and/or increase the number of nodes of the computational grid, since the fast numerical diffusion of vortices in the near wake substantially affects the field of the bevels of the flow in the layout field.
Validation of the calculation results according to the data of a pipe experiment, as shown by the authors' experience, makes it possible to achieve a matching drag in the calculation with the experiment, due to an increase in the
flow energy of turbulence. The degree of turbulence in the calculation may exceed the degree of turbulence of the flow at the entrance to the working part declared in the passport of the wind tunnel in the absence of a model. However, the course of the curves along the angle of attack, both for the drag coefficient and for the lifting force and longitudinal moment, can be corrected in the right direction (to improve coordination with experience) only by increasing the resolution (calculation quality) of the near vortex wake, for which, in particular, intended criteria, obtained on the basis of the principle of maximum pressure G.B. Sizykh.
Translated by K. Zykova
REFERENCES
[1] Vyshinsky V.V. Exact solutions generation of shock-free flow of symmetrical profile with local supersonic region. Uchenye zapiski TsAGI, 1975, vol. 6, no. 3, pp. 1-8 (in Russ.).
[2] Vyshinsky V.V., Sizykh G.B. The verification of the calculation of stationary subsonic flows and the presentation of results. Math. Models Comput. Simul, 2019, vol. 11, no. 1, pp. 97-106. DOI: 10.1134/S2070048219010162
[3] Rowland H. On the motion of a perfect incompressible fluid when no solid bodies are present. Am. J. Math., 1880, vol. 3, no. 3, pp. 226-268. DOI: 10.2307/2369424
[4] Lamb H. Hydrodynamics. Cambridge, University Press, 1916.
[5] Hamel G. Ein allgemeiner Satz uber den Druck bei der Bewegung volumbestandiger Flussigkeiten. Monatsh. f. Mathematik und Physik., 1936, vol. 43, no. 1, pp. 345-363. DOI: 10.1007/BF01707614
[6] Serrin J. Mathematical principles of classical fluid mechanics. Berlin, Springer-Verlag, 1959.
[7] Truesdell C. Two measures of vorticity. Indiana Univ. Math. J., 1953, vol. 2, no. 2, pp. 173-217. DOI: 10.1512/iumj.1953.2.52009
[8] Shiffman M. On the existence of subsonic flows of a compressible fluid. Proc. Natl. Acad. Sci. USA, 1952, vol. 38, no. 5, pp. 434-438. DOI: 10.1073/pnas.38.5.434
[9] Bers L. Mathematical aspects of subsonic and transonik gas dynamics. Wiley, Chapman and Hall, 1958.
[10] Jeong J., Hussain F. On the identification of a vortex. J. Fluid Mech., 1995, vol. 285, pp. 69-94. DOI: 10.1017/S0022112095000462
[11] Dubief Y., Delcayre F. On coherent-vortex identification in turbulence. J. Turbul., 2000, vol. 1, art. 11. DOI: 10.1088/1468-5248/1/1/011
[12] Lesieur M., Begou P., Briand E., et al. Coherent-vortex dynamics in large-eddy simulations of turbulence. J. Turbul., 2003, vol. 4, art. 16.
DOI: 10.1088/1468-5248/4/1/016
[13] Cala C.E., Fernandes E.C., Heitor M.V., et al. Coherent structures in unsteady swirling jet flow. Exp. Fluids, 2006, vol. 40, no. 2, pp. 267-276.
DOI: 10.1007/s00348-005-0066-9
[14] Bosnyakov S.M., Vlasenko V.V., Engulatova M.F., et al. Programmnyy kompleks dlya sozdaniya geometrii LA, sozdaniya mnogoblochnoy 3kh mernoy raschetnoy setki, polucheniya poley techeniya pri pomoshchi resheniya sistemy uravneniy Eylera i sistemy uravneniy Nav'ye — Stoksa, osrednennykh po vremeni i obrabotki rezul'tatov rascheta (EWT). Svidetel'stvo o gosudarstvennoy registratsii programmy dlya EVM № 2008610227 [Software package for formation of aircraft geometry, multiblock 3D grid and flow fields by solving system of time-averaged Euler and Navier — Stokes equations and by processing calculation data (EWT). State certificate of software registration no. 2008610227].
Anikin V.A. — Dr. Sc. (Eng.), Deputy Director General for Science, LCC "VR-Technologies" (Krasnopresnenskaya naberezhnaya 12, Moscow, 123610 Russian Federation).
Vyshinsky V.V. — Dr. Sc. (Eng.), Professor, Chief Researcher, Department of Aerodynamics of Airplanes and Rockets, Zhukovsky Central Aerohydrodynamic Institute TsAGI (Zhukovskogo ul. 1, Zhukovsky, Moscow Region, 140180 Russian Federation).
Pashkov O.A. — Cand. Sc. (Eng.), Leading Design Engineer, Department of Helicopter Aerodynamics, LCC "VR-Technologies" (Krasnopresnenskaya naberezhnaya 12, Moscow, 123610 Russian Federation).
Streltsov E.V. — Researcher, Zhukovsky Central Aerohydrodynamic Institute TsAGI (Zhukovskogo ul. 1, Zhukovsky, Moscow Region, 140180 Russian Federation).
Please cite this article as:
Anikin V.A., Vyshinsky V.V., Pashkov O.A., et al. Using the maximum pressure principle for verification of calculation of stationary subsonic flow. Herald of the Bauman Moscow State Technical University, Series Mechanical Engineering, 2019, no. 6, pp. 4-16. DOI: 10.18698/0236-3941-2019-6-4-16