CC BY
75
Software package to calculate the aerodynamic characteristics of aircrafts
V.N. Koterov <vkoterov@yandex.ru> V.M. Krivtsov <vlkrivtsov@yandex.ru> V.I. Zubov <vladimir.zubov@mail.ru> Dorodnicyn Computing Centre FRC CSC RAS, 40, Vavilov st., Moscow, 119333, Russia. Moscow Institute of Physics and Technology (State University), Institutskii per. 9, Dolgoprudnyi, Moscow oblast, 141700, Russia
Abstract. The software package to calculate parameters of three-dimensional steady and unsteady gas flows in complex devices is presented. The mathematical flow model used in package is based on the Reynolds-averaged Navier Stokes equations for a two-component equilibrium turbulent medium and a two-parameter semiempirical turbulence model. A numerical implementation of the software package to modelling three complex flows in modern devices is described.
Keywords: software package; mathematical simulation; gas flows; turbulence; exhaust jet; Navier Stokes equations; aircraft nozzle.
DOI: 10.15514/ISPRAS-2017-29(6)-17
For citation: Koterov V.N, Krivtsov V.M., Zubov V.I. Software package to calculate the aerodynamic characteristics of aircrafts. Trudy ISP RAN/Proc. ISP RAS, vol. 29, issue 6, 2017, pp. 271-288. DOI: 10.15514/ISPRAS-2017-29(6)-17
1. Introduction
The computation of characteristics of steady gas flows around prospective flight vehicles is a complicated task combining external and internal aerodynamics problems, which can be addressed with state-of-the-art computers. Formally, an aircraft can be represented as a streamlined body equipped with an internal gasdynamic nozzle. The ambient gas flow interacts heavily with the nozzle flow. Additionally, it should be noted that the flow characteristics in an extended region surrounding the nozzle are not known beforehand. The difficulties arising in the numerical simulation of a flow of this type are as follows:
• the flow is essentially three-dimensional and turbulent and, in the general case, can be of mixed type (sub- and supersonic);
• the computational domain has to be large, which is required for the flow characteristics near the aircraft to be reliably determined;
• relatively narrow boundary layers developing on the aircraft surface have to be resolved;
• the moving medium is multicomponent. In the general case, this is a mixture of fuel combustion products and the ambient air flowing past the aircraft;
• additionally, some difficulties are usually associated with the representation of numerically computed three-dimensional gasdynamic fields in a convenient form.
The goal of this work is to demonstrate developed mathematical model of the flow, numerical algorithm and a software package in order to estimate the influence exerted on the aircraft characteristics by user-specified input parameters.
2. Mathematical formulation of the problem
The moving gas was simulated as an equilibrium two-component mixture of ambient air (fraction 1) and fossil fuel combustion products (fraction 2). The flow of an equilibrium gas mixture was described by the three-dimensional unsteady Reynolds-averaged compressible Navier-Stokes equations closed with a model of eddy viscosity and thermal conductivity (RANS model) (see, for example, [1, 2]).
2.1 Basic Equations
In Cartesian coordinates, these equations can be written as
Here, t is time; x1 = x, x2 = y and x3 = z are the Cartesian coordinates; V=|| u ,||T is the flow velocity with the components ui, (i = 1, 2 ,3); 5,j - Kronecker Delta; p, p, T, h,
and Cp are the medium's density, pressure, temperature, specific enthalpy, and heat capacity at constant pressure, respectively; Ro is the universal gas constant; Cpi, mi, Cp2, and m2 are the respective heat capacities and molar masses of the combustion gases and air; pa is the reduced air density in the mixture; у is the mass fraction of the air in the mixture; and ks are the effective viscosity and heat conductivity of the medium; ц and цт are the molecular and eddy viscosities of the mixture; Pr and PrT are the laminar and turbulent Prandtl numbers; D2 is the effective diffusivity; and Sc and ScT are the laminar and turbulent Schmidt numbers. The eddy viscosity and heat conductivity were computed using Coacley's д-ю model [3]. In Cartesian coordinates, the equations of the q-ю model can be written as
Here, q and ra are the "pseudovelocity" and "pseudovortitity"; f(n) is the near-wall function introduced to correctly describe the flow parameters in the laminar sublayer developing at the walls; n is the normal distance from a given point to the nearest surface; S is a dissipative function.
2.2 Boundary Conditions
The conditions on the boundaries of the flow region must take into account the following data:
• The state of the unperturbed air flow far away from the aircraft (at the outer boundary of the computational domain and in the outlet cross section);
• The heat transfer conditions on the solid walls;
• The state of the flow in the inlet cross section of the nozzle (if present).
In the case of supersonic gas flows, the flow state away from the aircraft was specified by three parameters: the free-stream velocity V№ the static pressure pm and the static temperature Tw.
The boundary conditions on the solid walls were set as follows:
ut = 0, 0=1,2,3), T = Tw (or dT/dn = 0 ), q = 0, dra/dn = 0,
where Tw is a given wall temperature and n is the normal to the wall. 3. Computational Algorithm
The numerical solution to the boundary value problem stated above was obtained using finite differences.
3.1 Numerical Grid
The first step in the finite-difference implementation was the generation of a uniform three-dimensional mesh in the computational domain. This problem was an important inherent element of the numerical algorithm and was comparable in complexity with the other algorithmic parts. The accuracy of the results was determined to a large extent by the quality of the mesh used. The grid algorithm generates three-dimensional boundary-conforming structured meshes with hexahedral cells. Degenerate (pentahedral) cells could be used on the axis of symmetry of the domain (if any).
There are numerous commercial and open-access grid generators for constructing boundary-conforming hexahedral structured meshes (see [4-6]). As a rule, algorithms for generating such grids involve two principal stages:
• The construction of a parametric representation of the surfaces bounding the computational domain (including grid node arrangement on these surfaces);
• The generation of a three-dimensional structured mesh inside the domain. In domains of complex geometry, these generators can sometimes produce degenerate meshes with self-intersecting cells. It should also be kept in mind that, even in the case of nondegenerate meshes, a fairly accurate solution for internal and external turbulent gas flows can be obtained if the equations and the boundary conditions are well approximated, which imposes rather severe restrictions on the quality of these meshes.
The term "quality of a numerical grid" cannot be defined rigorously. In this case, by such we mean the following collection of intuitive and, generally speaking, quite incompatible requirements:
• The mapping to a parametric parallelepiped used in structured mesh generation must be continuously differentiable and its Jacobian must not vanish.
• To resolve the boundary layers, the mesh must be refined toward the solid walls. Moreover, the transverse (to the flow direction) mesh size near these surfaces can be 10-4 and less of the characteristic length of the problem.
• The grid lines near the walls must be nearly orthogonal to them. Wherever possible, the mesh cells must be similar to parallelepipeds.
• The numerical grid must be quasi-uniform in the sense that, as the number
of grid nodes tends to infinity, the difference between the characteristic lengths of neighboring cells tends to zero faster than the lengths themselves. Anyway, the ratio of the characteristic lengths of neighboring cells must not be too large. The numerical results presented below were obtained on grids generated by the algorithms described in [6].
3.2 Conservative Approximation of Spatial Operators
The approximations of the spatial operators used in the algorithm are described as applied to the Navier-Stokes equations given in Section 2. These approximations are underlain by the following identities, which are often used as definitions of the corresponding operators:
where n is the inward normal to the boundary r of the domain Q, |Q| is the volume of Q, and g and p are arbitrary vector and scalar fields. Note also that such approximations lead to skew-symmetric matrices approximating first spatial derivatives.
To explain what was said above, we consider the case of flat structured meshes consisting of quadrilaterals. Cells in structured meshes are numbered by indices. Each face is adjacent to two cells, whose indices differ in one position. The faces are numbered by half-integer indices: (i -1/2, j ) , (i + 1/2, j), (i, j -1/2), and (i, j+1/2). Assume that all the sought quantities g and f are associated with cells (in the smooth case, with barycenters, i.e., the centers of mass of cells). To use previous formulas, the sought quantities have to be extended to faces. To recover the values 9 of the function 9 on faces, we introduce two extension procedures, M+ and M_. The input data are the values of 9 on both sides of a face. For the face between the cells (i, j) and (i, j+1), the procedures are defined as
For the face between the cells (i,j) and (i + 1, j), they are defined in a similar manner. Note that these extension procedures ensure only the first order of accuracy.
For simplicity, we consider the approximation of the derivative
(Pij+i/2 - Л/+У+1/2(ф,>Ф;.;+|) - фu, фiJ+l/2 - M':J+4\9,j,9u+1) - (p/J+|.
where nx is the projection of the inward normal onto the x axis.
To use the above formula, we need to know the values of f on the faces of the current cell. Assume that, on a structured quadrilateral mesh, the face extensions are produced by the procedure M+. Then the above derivative is approximated as
,+ (f*Snx)i_UXj + (f+Snx)M/2J + (f+Snx)kH/2 + (f+Snx)u+1/2
where S is the edge length and |Qy| is the cell area. For each cell (i, j), we have
Consider the scalar product of the derivative ff and an arbitrary compactly supported function 9:
In this sum, we group the terms involving fj. There are four of them: two correspond to the cell (i, j) and a single term corresponds to each of the cells (i +1, j) and (i, j +1); specifically,
In the terms corresponding to the cell (i, j), the sum of edge lengths multiplied by the corresponding components of the normals is replaced using the above equality. As a result, we obtain the relation
where the derivative q> , of cp is given by the same formula as for ff but with cp extended to the faces with the help of the procedure M_.
The first y-derivatives fy of the function f can be approximated in a similar fashion. Using the basis vectors ex and ey , we can construct the operations grad±f div±g, and curl±g.
To determine the order of accuracy of the proposed formulas, we consider a uniform rectangular mesh with edge lengths .V, and S,-. The formula for computing the derivative fx iJ involves fh]]/ and fLh\. which can be represented as
After substituting these representations into the formula for fx,ij- it is easy to see that the truncation error in the approximation of ( flex by ff is first order.
Note that the leading terms of the truncation errors in f j and fx are identical in absolute value but opposite in sign. Therefore, the derivative с flex is approximated with second order accuracy by the combination 0.5(fx,i,j + fx,¡4) in a cell and by
0.5(f*jj + fXA, 1 j and 0.5(fXIJ + fxij,,) on the corresponding edge. Tliis idea was used in [7].
An alternative approach to an increase in the order of accuracy is to extrapolate or interpolate function values from the cells adjacent to the given face. Next, the values on cell boundaries can be obtained using exact (iterative) or approximate Riemann solvers (see [8]). This approach was implemented in a software code.
3.3 Approximation of Viscous Operators
The second-order operators related to viscosity are approximated by applying the variation principle. Specifically, the viscous operators in momentum equations (see Section 2) can be obtained by varying a special functional with a symmetric nonnegative quadratic part
/(V) = j[/(V) - 2(V ■ ®)]dxdydz
with respect to the components u, v and w of the velocity V; here,
is a dissipation function, 0 is the vector 0=(div Tx, div x„ div Tz), and t is the viscous stress tensor. The first variation of functional I(V) is computed in the class of functions satisfying given boundary conditions.
In each mesh cell c, the components of the velocity gradient are determined as described for the pressure gradient. Then the approximate value Jh (V) of functional I(V) is calculated as
A(V) = Z{°-5[/+ + - 2<V . 0)}|Qj,
where the sum is taken over all cells of the computational domain. |fl.| is the volume of the c-th cell, and I are approximations of /(V) with first derivatives determined by the procedures M+ andM_, respectively. Note that
Relying on this inequality, a simple algorithm can be constructed in which the viscous terms in the considered equations are taken into account in the form of correction terms at the next iteration (see [9]).
3.4 Implicit Scheme and Iterative Algorithm
Governing Equations are numerically solved using a two-level fully implicit difference scheme of the form
where 9 is a grid vector function involving all unknown functions, k is the time level index, At is the time step, H is an operator containing first and second difference derivatives with respect to spatial variables, and A-1 is an operator improving the convergence of the iterations (in the simplest case, A = E, where E is the identity operator). The norm of H(9k) is a measure indicating the proximity of the solution to a steady state.
We use conservative (flux) approximations of the spatial operators. Conservative variables are associated with mesh cells, while flux variables, with cell boundaries. The resulting system of nonlinear algebraic equations is solved using the following simple iterative procedure at each grid node:
Here, 5 is the iteration number and 4 > 0 is an iteration parameter.
Note that one step of this iterative procedure with 4=1 is equivalent to computation
based on an explicit scheme.
At every time step, the iterations are continued until || <^(<P("))|/| (*p"o>) || < e' -where e is the prescribed accuracy of these inner iterations. It turns out that the iterative error introduced at the k-th time step is not accumulated at the subsequent steps but decays. Accordingly, in the case of nonstationary problems, we can use e ~ 0.1, which ensures the stability of the computation. In the case of a linear operator H and A=E, necessary conditions for the convergence of iterative process were obtained in [7].
Consider the general case in more detail. Let the discrete approximation of the spatial derivatives be represented as
Here, the first term approximates the first derivatives (convective terms of the equations) and the second term approximates the viscous terms. Both terms are nonlinear functions of their arguments. To analyze them, we consider their linear approximation of the form
where the coefficients of the corresponding matrices T and V are calculated using the values of the unknowns from the preceding iteration step. Let the operator A be given by
where ||K'|| < |l "| and the matrix A is easily invertible (see [13]). Then we have the estimate
This estimate is used to analyze the convergence of the iterative process mentioned. Indeed, the matrix T approximates the first derivatives of the unknown quantities. It is determined by the approximation of only the convective terms. As a result, the norm of the operator A-1H' to be inverted can be substantially reduced in the case of fine grids (see also [10]).
4. Numerical Results
Several interesting problems were numerically solved by applying the algorithm described above and developed software package. Some of the results are presented below.
4.1 Gas Flows around Aircraft with Allowance for the Flow/Exhaust Jet Interaction
Below are the numerical results obtained for the complex gas flow around a prospective flight vehicle. An external view of this vehicle is shown in Fig. 1.
Fig. 1. External view of vehicle.
The vehicle is a blunt cylinder of radius R=0.3 m with a spherical cap of radius R/10 (the distance along the axis of symmetry from the cylinder to the outer boundary of the sphere is 1.9 R). The cylinder and the cap are adjusted by a spline surface of revolution around of the cylinder's axis of symmetry. The nozzle radius is 80% of the radius of the vehicle rear end. The tail of the vehicle is equipped with rudders. It is assumed that fuel combustion products are exhausted from the nozzle of the vehicle. The exhaust jet has a high temperature. In the general case, the physical characteristics of the jet differ noticeably from those of the ambient air. The task is to determine the flow parameters near the vehicle and some of its integral characteristics.
Fig. 2. Mach number distribution.
H - TfiS' ii
■ ■■WOT
■ ■ VCWP11 -.csiliiv
■ ■ WEI«
Fig. 3. Temperature distribution.
Fig. 4. Fraction distribution.
Fig. 2 shows the Mach number distributions in a longitudinal cross section of the computational domain containing a rudder. As expected, the Mach number decreases near the body surface, the rudders, and in the gas mixing zone. Inspection of the figure reveals an external shock wave determined by the overall vehicle shape, an internal shock wave arising due to the interaction of the outer airflow with the combustion products, and a rarefaction wave caused by the airflow past the convex part of the body. Figs. 3 and 4 present the temperature and fraction distributions near a rudder. The temperature near body surface, the rudders, and in the mixing zone increases and becomes close to that of the combustion products exhausted from the nozzle. The fraction distribution near a rudder shows that the combustion products nearly do not propagate upstream toward the rudders.
Fig. 5. Mach number distribution.
Fig. 6. Temperature distribution.
Fig. 7. Fraction distribution.
Figs. 5-7 illustrate the distributions in a transverse cross section of the computational domain near the rear end. As in the case of longitudinal cross sections, the Mach number near the body surface and the rudders is observed to decrease. The temperature near the body and the rudders increases and becomes close to that of the combustion products exhausted from the nozzle. The fraction distribution shows that the combustion products do not propagate upstream toward the rudders.
An analysis of the results suggests that complex flows of a mixture of widely different gases can be fairly well described by applying the algorithm developed. This algorithm well reproduces external and internal shock waves, rarefaction waves, and gas mixing zones. A priori information on the features of a flow (obtained, for example, via computations on a coarse grid) can be used to improve the quality of the computed flow by the refining mesh in areas of strong variations in the flow parameters.
4.2 Three-Dimensional Turbulent Gas Flows in Complex Nozzle Systems
The computation of gas flow parameters in nozzles is a classical numerical problem in internal aerodynamics. Another classical problem is that of external aerodynamics consisting in the computation of gas flows past various bodies. In some cases, problems combining the features of flows of both types have to be considered in practice.
An example of such problems combining internal and external aerodynamics is the computation of the steady gas flow in a no axisymmetric supersonic ejector nozzle system. Fig. 8 displays an example of a nozzle system of this kind. The system consists of an internal gasdynamic (primary) nozzle and a surrounding ejector contour. Such a primary nozzle is placed, for example, at the exit of a jet combustor with the combustion gases being exhausted through it. The ejector contour is a
Fig. 8. Ejector nozzle.
system of air inlets designed so that the air stream from the surrounding space mixes with the hot combustion gas flow at the exit from the primary nozzle. The gas flows in the different regions of the ejector nozzle strongly interact with each other. It should be mentioned that the flow characteristics in the rather long external-flow region is not known in advance.
Fig. 9. Temperature distribution in the middle section.
Fig. 10. Streamlines in a vertical plane through the interwall spacing.
Some results of test computations imitating a cold flow are presented. In this case, air is drawn in the primary nozzle and air flows past the system. Figs. 9-11 display some of the computed gasdynamic fields. Fig. 9 shows the distribution of the static temperature (°K) in the middle section x0y of the nozzle. The behavior of some streamlines for gas flows is shown in Figs. 1°, 11.
The figures demonstrate a complex flow pattern, various vortex regions inside the nozzle (Fig. 1°), and a vortex sheet forming behind the nozzle (Fig. 11). For more details of this work and comparison with results of other authors see [11].
4.3 Simulation of gas flow in cooled axial turbines
Numerical simulations of the working flow in gas turbines are the complex problem of internal aerodynamics. The flow under consideration occurs in the regions of complex shapes, representing the channels between rotor and stator blades moving relative to each other. In addition, from surfaces of the turbine's blades can usually be blowing the cooling air.
Software package based on algorithm described above are considered all flows within all channels between the blades of each blade row (stators and rotors) as the same, that the process of generation of long wave axial disturbances is ignored. The governing equations of gas flow within the interblade channels of stators are written with the use of a fixed system of coordinates and within the interblade channels of rotors, these equations are written with the use of a rotating system of coordinates. The numerical solution for the fixed and rotating regions are joint at the expense of requirements based on the continuity of the mass, momentum end energy fluxes. The calculation of these fluxes are carried out with the use of axial average values. This axial average allows to avoid the high frequency pulsation of the flow and to reach its average stationary state with the use of a numerical iterative procedure. Some results of computations are shown on Fis. 12-14.
Fig. 11. Streamlines for the flow behind the nozzle exit..
Fig. 12. Numerical grid in middle radial and axial plans.
Fig. 13. Calculated temperature in first stage stator of turbine in middle radial and axial plans.
Fig. 14. Calculated Mach number in first stage rotor of turbine in middle radial and axial plans.
5. Conclusion
The submitted results show the ability to use the software package for solving wide range of stationary and non-stationary aerodynamic problems
References
[1]. Nigmatulin R.I. Dynamics of Multiphase Media. Part 1. Nauka, Moscow, 1987 (in Russian)
[2]. Lapin Yu.V. and Strelets M.Kh. Internal Flows of Gas Mixtures. Nauka, Moscow, 1989 (in Russian).
[3]. Coacley T.J., Turbulence Modeling Methods for the Compressible Navier-Stokes Equations. AIAA Paper, no. 83-1693, 1983.
[4]. Koterov V.N. Three-dimensional grid generation in multistage axial turbines based on the variational barrier method. Comput. Math. Math. Phys., vol. 45, no. 8, pp. 13251333, 2005.
[5]. Official home of The Open Source Computational Fluid Dynamics (CFD) Toolbox. Available at: http://www.openfoam.com/
[6]. Garanzha V.A., Kudryavtseva L.N., and Utyzhnikov S.V. Untangling and optimization of spatial meshes. J. Comput. Appl. Math., no. 269, pp. 24-41, 2014.
[7]. Zubov V. I., Inyakin V. A., Koterov V. N., and Krivtsov V. M. Numerical simulation of three-dimensional turbulent gas flows in complex nozzle systems. Comput. Math. Math. Phys. vol. 45, no. 10, pp. 1802-1814, 2005.
[8]. Rodionov A.V. Monotonic scheme of the second order of approximation for the continuous calculation of nonequilibrium flows. USSR Comput. Math. Math. Phys. vol. 27, no. 2, pp. 175-180, 1987.
[9]. Krivtsov V.M. On a numerical scheme for solving the Navier—Stokes equations. USSR Comput. Math. Math. Phys., vol. 26, no. 3, pp. 172-178, 1986.
[10]. Saad Y. Itertive methods for sparse linear systems. 2nd edition with corrections. — SIAM, 2003.
[11]. Zabarko D.A., Zubov V.I., Kotenev V.P., Krivtsov V.M., Polezhaev Yu.A. Numerical Simulation of Gas Flows around Aircraft with Allowance for the Flow/Exhaust Jet Interaction. Comput. Math. Math. Phys., vol. 55, no. 4, pp. 677-689, 2015.
Программный пакет для расчета аэродинамических характеристик летательных аппаратов
В.Н. Котеров <vkoterov@yandex.ru> В.М. Кривцов <vlkrivtsov@yandex.ru> В.И. Зубов <vladimir.zubov@mail.ru> Вычислительный центр им. А.А. Дородницына РАН Федерального исследовательского центра «Информатика и управление» РАН, 119333, Россия, Москва, ул. Вавилова, 40 Московский физико-технический институт, 141700, Россия, Московская облазь, г. Долгопрудный, Институтский пер., 9
Аннотация. Представлен программный пакет для расчета параметров трехмерных стационарных и нестационарных потоков газа в сложных устройствах. Модель математического потока, используемая в пакете, основана на усредненных по Рейнольдсу уравнениях Навье-Стокса для двухкомпонентной равновесной турбулентной среды и двухпараметрической полуэмпирической модели турбулентности. Описана численная реализация программного пакета для моделирования трех сложных потоков в современных устройствах.
Ключевые слова: программный пакет; математическое моделирование; потоки газа; турбулентность; реактивная струя; уравнения Навье-Стокса; летательный аппарат
DOI: 10.15514/ISPRAS-2017-29(6)-17
Для цитирования: Котеров В.Н., Кривцов В.М., Зубов В.И. Программный пакет для расчета аэродинамических характеристик летательных аппаратов. Труды ИСП РАН, том 29, вып. 6, 2017 г., стр. 271-288 (на английском языке). Б01: 10.15514Л8РРЛ8-2017-29(6)-17
Список литературы
[1]. Нигматулин Р. И. Динамика многофазных сред. Часть I. М.: Наука, 1987, 464 стр.
[2]. Лапин Ю.В. Стрелец М.Х. Внутренние течения газовых смесей. М.: Наука, 1989, 368 стр.
[3]. Coacley T.J., Turbulence Modeling Methods for the Compressible Navier-Stokes Equations. AIAA Paper, no. 83-1693, 1983.
[4]. Koterov V.N. Three-dimensional grid generation in multistage axial turbines based on the variational barrier method. Comput. Math. Math. Phys., vol. 45, no. 8, pp. 13251333, 2005.
[5]. Official home of The Open Source Computational Fluid Dynamics (CFD) Toolbox. http://www.openfoam.com/
[6]. Garanzha V.A., Kudryavtseva L.N., and Utyzhnikov S.V. Untangling and optimization of spatial meshes. J. Comput. Appl. Math., no. 269, pp. 24-41, 2014.
[7]. Zubov V. I., Inyakin V. A., Koterov V. N., and Krivtsov V. M. Numerical simulation of three-dimensional turbulent gas flows in complex nozzle systems. Comput. Math. Math. Phys. vol. 45, no. 10, pp. 1802-1814, 2005.
[8]. Rodionov A.V. Monotonic scheme of the second order of approximation for the continuous calculation of nonequilibrium flows. USSR Comput. Math. Math. Phys. vol. 27, no. 2, pp. 175-180, 1987.
[9]. Krivtsov V.M. On a numerical scheme for solving the Navier—Stokes equations. USSR Comput. Math. Math. Phys., vol. 26, no. 3, pp. 172-178, 1986.
[10]. Saad Y. Itertive methods for sparse linear systems. 2nd edition with corrections. — SIAM, 2003.
[11]. Zabarko D.A., Zubov V.I., Kotenev V.P., Krivtsov V.M., Polezhaev Yu.A. Numerical Simulation of Gas Flows around Aircraft with Allowance for the Flow/Exhaust Jet Interaction. Comput. Math. Math. Phys., vol. 55, no. 4, pp. 677-689, 2015.
