CC BY
5УДК 519.63, 532.511 DOI: 10.34759/trd-2021-119-04
Бессеточный алгоритм расчёта сверхзвуковых течений невязкого газа
Способин А.В.
Московский авиационный институт (национальный исследовательский университет), Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия
e-mail: spise@inbox.ru
Статья поступила 10.07.2021
Аннотация
Изложен бессеточный метод численного решения системы уравнений движения невязкого сжимаемого газа в двумерной и трёхмерной постановке. Описан способ построения набора узлов расчётной области для заданной геометрии. На примере решения модельной задачи численного моделирования сверхзвукового течения газа вокруг сферы продемонстрирована применимость рассматриваемого подхода к задачам обтекания. Проведённый сравнительный анализ показал хорошее согласование полученных результатов моделирования с другими вычислительными методами.
Ключевые слова: численное моделирование, бессеточный метод, уравнения Эйлера, сверхзвуковое обтекание тел.
Введение
В работе рассматривается решение актуальной задачи численного моделирования обтекания затуплённого тела сверхзвуковым невязким потоком газа [1,
2] в трёхмерном пространстве. Подобная задача может быть решена довольно широким спектром известных методов, в частности, для её решения в двумерной постановке успешно использовался метод погруженной границы с фиктивными ячейками на декартовых сетках [3, 4], посредством которого в рамках серии работ [5-7] было проведено исследование эволюции течения в ударном слое при движении крупной частицы вдоль оси симметрии. Продолжение указанных исследований для асимметричных течений потребовало решения задачи моделирования течения газа в трёхмерной постановке. Выбор бессеточного метода решения уравнений газовой динамики [8-11] обусловлен меньшими требования к вычислительным ресурсам, прежде всего, объему оперативной памяти вследствие анизотропии решётки по сравнению с реализованными в предыдущих работах [4-7] адаптивными методами на декартовых сетках, а также необходимостью моделирования движения частицы посредством перемещения вычислительных сеток, что представляет сложность для конечнообъемных методов на криволинейных сетках, обладающих, что следует отметить, как правило, большей точностью по сравнению с бессеточными методами.
Математическая модель и численный метод
Система дифференциальных уравнений Эйлера описывает течение невязкого сжимаемого газа в трёхмерном пространстве
да д¥ (д) | дС (д) < дН (д)
+
+
д? дх ду д2
q
Г ри л ( Рv Л Г рw л
ри 2 ри + р рш> рuw
Pv , Р = рuv , с = 2 рv + р , н = рvw
рw рuw рvw рм/1 + р
кРе у ч риН у ч р^ у ч РwH у
где р - плотность, р - давление, и, V, w - компоненты скорости газа, у - показатель адиабаты.
Пространственная дискретизация уравнений выполняется согласно методу наименьших квадратов [10-12]. Для каждого узла г, облака окружающих его соседних узлов ] е С и заданной в пространстве некоторой функции р( х, у, 2 ):
л др (р] = р + Дху —
+ Ду,
др
V
ду
Л др
дг
+ 0
(*2)
1
Дхг] = Х] - X, ДУ] = У] - Уг, Д2Ч = 2] - 2г, ДР] = Р] - Р ,
Л2
С
х ^
]еСг
V
Л Л др
Др] -Дху дх
-Ду
др
]
ду
Л др
- Дг„ —
и дг
^ min.
г у
Аппроксимация производных функции р по пространству осуществляется посредством решения системы уравнений
др
дх
X адрч-
др
]еСг
ду
= X в]Др
г ]еСг
др
X П] р
]еСг
в
ц
Е ®уАх| Е ^]АхуАУу Е ®уАхуАгу
Е Е
Е Е
Е АхуАУу Е АУ(! Е ®уАУуАгу
Е ауАхуАгу Е ауАУуАгу Е
Е
№
Е
-1 1- -|
Е ®у Аху
Е ®у АУу
Е ®у Агу
_г _
Полученные для функции р соотношения используются для аппроксимации частных производных газодинамических величин, таких как плотность, давление, температура, проекций скорости газа на координатные оси, компонентов потоков Р, С, Н и т.д. В соответствии с полученными приближениями уравнения газовой динамики могут быть записаны в полу-дискретной форме:
д?
1+ 2 Е \_ау (Ру - Р ) + ву (Су - С Н (Ну - Н )
= 0
здесь р = Р (аг-), Сг- = С (аг-), Нг- = Н (дг-), множитель 2 обусловлен вычислением потоков Ру, Су, Ну в середине отрезка, соединяющего узлы г и у. Согласно методу Хартена-Лакса-ван Лира [13], используются выражения:
Р = ГУ
р (а ь)
Чр (а ь)- Чр (а к) + (а к - д z):
як -ЯК
Р(ая)
о <Ль ЧЬ < 0 <ЧЯ , 0 <ЧЯ
а ь
ау ,хг < ху
4+, х г > х у у
а я
ау ,х г < х у
У
а- ,хг > ху
Собственные значения XL и Xr рассчитываются по Roe-усреднению вектора
физических переменных: Xl = u* - c*, Xr = u* + c*,
P =4plpr ,
*
u =
= ul4p~l + ur4pr 4pl+4pr
v
vl4pl + vRy[PR
4pl+4pr '
w
wL
4pL + wRy[PR *_ eLyfPL + eRy/PR
4pl +4Pr ' 4pl +4pr
*
c =
* 2 *2 * 2 * u + v + w e--
r(r~ 1)
Вычисления потоков С] и Н] вдоль осей у и г аналогичны приведённым выше.
Для реконструкции векторов консервативных переменных q±j = q () к компонентам вектора Т = (р, u, v, w, e) применяется схема MUSCL с ограничителем van Albada [11]:
Vij = Щ+-
(1 - ksi) Aj + (1 + ksi)(vj - Щ) , A- = 2Arj • V щ - (pj - щ),
Vj =pj-T
(1 - kSj) A j + (1 + kSj )(щ - щ)], Aij = 2Arj ^Vj - (щ -щ),
Sj = max
2A ij (Vj -Vi ) + £ A-+(vj -Vj)2 + s
Sj = max
2A + (vj -Vi ) + £ 2
A+2+(vj -Vi)
'у
С \
ху хг
Уу - У К 2 у - 2 У
, Vу
п
Г ду \ /
дх п
ду
дУ п
ду
дг п У V
Е атп (ут уп )
т<ЕС„
Е Ртп (Ут-Уп )
Е Ттп (ут уп)
1 -П
значения параметров к = —, б = 10 .
Интегрирование системы дифференциальных уравнений по времени выполняется явным методом Рунге-Кутты третьего порядка:
+я (аг ) = о, я (аг ) = 2 Е[а (% - Р ) + ву (С у - С г ) + Уу (Ну - Н )""
С (л\\
а «=а п-а? я (а п), а (2 } = | а п + \ а ?)-\а? я I а (
(1)
/
к У
а п+1 = 3 а п + 2 а (2}- | а? я
' ( 2 а/
к У
Шаг по времени при явном интегрировании определяется согласно критерию Куранта: СБЬ
А?
тах
/
ах| ч >ч=е [у+У+уу+У а+вУ,сть=05..
Описанный метод может быть легко адаптирован для численного решения уравнений Эйлера в двумерной постановке, в этом случае параметры газа по оси г
др
принимаются тождественно равными нулю w = 0, у у = 0
д2
= 0, а метод
наименьших квадратов применятся для аппроксимации производных функции р по координатам x и у :
1
а
V
вг
У
Е ®у ^У Е ®у ДхУ ДУу
Е
Е ®у АхгV ДУу Е ®у Ауу
Е
Е Дхг]
^Сг
Е ^ дуу
Решения осесимметричной задачи может быть выполнено в цилиндрической системе координат, в этом случае система уравнений Эйлера принимает вид [14]:
ддув , дс2Р(д)+1 Кс = 0,
Ы дх ду г
д 2 в =
Г ри л Г ру Л Г р л
ри N С) 2 ри + р , С 2В - риу 2 , Кс - риу 2
ру риу ру + р ру
ч риН у руН у
здесь х выступает в качестве оси симметрии, у ей ортогональна, г - расстояние от точки до оси симметрии.
Выполнена программная реализация описанных алгоритмов в двумерной и трёхмерной постановке на языке программирования С++ с применением технологии распараллеливания вычисления ОрепМР [15, 16].
Вычислительный эксперимент
В качестве модельной решалась задача сверхзвукового обтекания сферы радиусом = 55 мм потоком воздуха. На входной границе расчётной области задаётся граничное условие первого рода, определяющее невозмущённый сверхзвуковой поток при числе Маха М^ = 3. На выходной границе области заданы
однородные условия Неймана.
Поверхность сферы рассматривается как твёрдая стенка с условиями
др др рТ дуТ
непротекания: —= 0, = , —т- = 0, уп = 0, где уп, уТ - нормальная и
дп дп дп
касательная компоненты скорости газа на поверхности. Реализация граничных условий
др также основана на применении метода наименьших квадратов [17]:
дп
др
Епуру дп
р = ^-, Пу = аупх + впУ + Ууп2
Е п
у^СЬ
У
где пх, пу , п2 - компоненты вектора внешней нормали п в узле I на границе
поверхности, Сь - множество его соседних узлов, не принадлежащих границе.
Во внешней области расчёта распределение узлов представляет собой прямоугольную декартову сетку с равномерным шагом. Вблизи обтекаемой поверхности узлы располагаются на поверхности концентрических сфер, расстояние между которыми увеличивается линейно по мере удаления от поверхности тела пропорционально радиусу (см. рис. 1). Узлы на каждой из сфер распределяются
относительно равномерно согласно следующему алгоритму [18]. Для заданного числа
узлов N на сфере единичного радиуса вычисляются площадь элемента a =
4п
N
количество окружностей Мд
п \[а
и угловой шаг (д
п
м
по «широте» между
в
ними. Для каждой окружности с широтой угла в
( 1Л п т +-
I 2 у
т
м
в диапазоне
в
0 < т < Мд -1 определяется количество принадлежащих ей узлов
М.
(рт
2пЛв sinв
т
а
координаты каждого из которых определяются выражениями:
гX^
У
V у тп
Г^пвт Ртп^
^пвт ^П(тп V с^вт у
2пп
(тп =-, 0 < п < М
М,
рт
1. Пример распределения узлов
рт
приведён на рис. 2.
В области стыковки сеток при слишком близком расположении соседних узлов или плохой обусловленности матрицы аппроксимации производных осуществляется коррекция положения узлов декартовой сетки в направлении внешней нормали к обтекаемой поверхности. Расположение расчётных узлов в сечении области
плоскостью 0Ху, проходящей через центр обтекаемой сферы приведено на рис. 1.
Рис. 1. Распределение узлов в центральном сечении расчётной области
Рис. 2. Расположение узлов на сферической поверхности
На рис. 3 представлены полученные в вычислительных экспериментах теневые Шлирен-картины течения газа [19]. Штриховой кривой отмечено положение головной ударной волны, которое определяется приближённо-аналитическим выражением [20]:
х
(у )=
я-а + я« (м х
(м I -1)
1 +
у
я«2 (м I -1)
Яс -1.143 • Я • ехр
г \
0.54
,1.2
, А-0.143•Я, -ехр
г \ 3.24
чМ2 у
(М ю-1)1
Можно отметить большую степень схожести теневых картин течения, полученных различными вычислительными методами, и хорошее соответствие между
аналитическим и численным пространственным положением скачка уплотнения.
а б
Рис. 3. Теневая картина течения (Шлирен): бессеточный метод (а), метод погруженной границы на декартовых сетках (б).
Рис. 4. Поле давления вблизи поверхности сферы при М ^ = 3: бессеточный метод (а), метод погруженной границы на декартовых сетках (б).
Рис. 5. Поле плотности газа вблизи поверхности сферы при М ^ = 3: бессеточный метод (а), метод погруженной границы на декартовых сетках (б).
Рис. 6. Поле числа Маха вблизи поверхности сферы при М ^ = 3: бессеточный метод (а), метод погруженной границы на декартовых сетках (б).
На рис. 4-6 показаны поле давления, плотности и числа Маха, полученные рассмотренной выше реализацией бессеточного метода и ранее неоднократно апробированным двумерным методом погруженной границы с фиктивными ячейками
на декартовых сетках в осесимметричной постановке, также основанном на схеме
Хартена-Лакса-ван Лира [3, 4]. Результаты решения задачи обтекания сферы двумя различными методами очень близки, метод, использующий декартовы сетки, демонстрирует несколько большую диссипацию вблизи головного скачка уплотнения.
121 10 8
теБЫеББ 1_иЫтОу. Риаапоу
\
\
\
30 60
а
а
90
4.5 4 3.5 3
0^2.5 О. 2 1.5 1
0.5 0,
те£|п1е£5 ииЫтоу, Иизапоу
\
V
30 60
а
90
Рис. 7. Распределение давления (а) и плотности (б) газа на поверхности сферы.
На рис 7 приведено сравнение давления и плотности газа на поверхности сферы, полученных бессеточным методом, с эталонными значениями [21].
В целом получено хорошее согласование результатов расчётов бессеточным методом с эталонными данными и результатами расчётов методом на декартовых сетках, что свидетельствует о применимости описанного метода к рассматриваемому кругу задач, а именно - численному моделированию обтекания тел сверхзвуковым
газовым потоком.
Заключение
Описанная математическая модель и разработанное программное обеспечение представляют собой этап в реализации проекта в области численного исследования обтекания тел запылёнными потоками. Следующие ступени предполагают моделирование течений вязкого газа путём решения уравнений Навье-Стокса, определение тепловых потоков от несущей фазы к обтекаемой поверхности, использование графических процессоров и вычислительных кластеров для расширения масштаба решаемых задач.
Работа выполнена в рамках государственного задания Минобрнауки России,
номер темы FSFF-2020-0013.
Библиографический список
1. Ревизников Д.Л., Сухарев Т.Ю. Гиперзвуковое обтекание затупленных тел в условиях атмосферы Земли и Марса. Сравнительный анализ математических моделей // Тепловые процессы в технике. 2018. Т. 10. № 1 - 2. С. 5 - 15.
2. Быков Л.В., Никитин П.В., Пашков О.А. Математическое моделирование процессов обтекания затупленного тела высокоскоростным потоком // Труды МАИ. 2014. № 78. URL: http://trudymai.ru/published.php?ID=53445
3. Винников В.В., Ревизников Д.Л. Метод погруженной границы для расчета сверхзвукового обтекания затупленных тел на прямоугольных сетках // Труды МАИ. 2007. № 27. URL: http://trudymai.ru/published.php?ID=34018
4. Винников В.В., Ревизников Д.Л., Способин А.В. Двухфазный ударный слой при обтекании тел сверхзвуковым запыленным потоком // Математическое моделирование. 2009. Т. 21. № 12. С. 89 - 102.
5. Ревизников Д.Л., Способин А.В., Иванов И.Э. Изменение структуры течения под воздействием высокоинерционной частицы при обтекании тела сверхзвуковым гетерогенным потоком // Теплофизика высоких температур. 2018. Т. 56. № 6. С. 908 -913. DOI: 10.31857/S004036440003569-9
6. Reviznikov D.L., Sposobin A.V., Ivanov I.E. Oscillatory flow regimes resulting from the shock layer - particle interaction // High Temperature, 2020, vol. 57, no. 2, pp. 278 - 283. DOI: 10.1134/S0018151X20020169
7. Способин А.В., Ревизников Д.Л., Иванов И.Э., Крюков И.А. Колебания давления и теплового потока, индуцированные газодинамическим взаимодействием высокоинерционной частицы с ударным слоем // Известия высших учебных заведений. Авиационная техника. 2020. № 4. C. 108 - 115.
8. Толстых А.И., Широбоков Д.А. Бессеточный метод на основе радиальных базисных функций // Журнал вычислительной математики и и математической физики. 2005. Т. 45. № 8. С. 1498 - 1505.
9. Vasilyev A.N., Kolbin I.S., Reviznikov D.L. Meshfree computational algorithms based on normalized radial basis functions // International Symposium on Neural Networks, 2016, DOI:10.1007/978-3-319-40663-3_67
10. Sattarzadeh S., Jahangirian A. 3D implicit mesh-less method for compressible flow calculations // Scientia Iranica, 2012, vol. 19, no. 3, pp. 503 - 512. DOI: 10.1016/j.scient.2012.04.013
11. Y. Wang, X. Cai, M. Zhang, X. Ma, D. Ren, J. Tan The study of the three-Dimensional meshless solver based on AUSM+-up and MUSCL scheme // International Conference on Electromechanical Control Technology and Transportation, 2015. DOI: 10.2991/icectt-15.2015.52
12. Дринфельд Г.И. Интерполирование и способ наименьших квадратов. - Киев, Вища школа, 1984. - 103 с.
13. Harten A., Lax P.D. and van Leer B. On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws // SIAM Review, 1983, vol. 25, no. 1, pp. 35 -61.
14. Молчанов А.М. Математическое моделирование задач газодинамики и тепломассообмена. - М.: Изд-во МАИ, 2013. - 206 с.
15. Антонов А.С. Параллельное программирование с использованием технологии OpenMP. М.: Изд-во МГУ, 2009. - 79 с.
16. Малявко А.А. Параллельное программирование на основе технологий OpenMP, CUDA, OpenCL, MPI. - М.: Изд-во Юрайт, 2021. - 135 с.
17. Hashemi M.Y., Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations // Journal of Computational and Applied Mathematics, 2011, no. 235, pp. 4687 - 4700. D01:10.1016/j.cam.2010.08.002
18. Deserno M. How to Generate Equidistributed Points on the Surface of a Sphere, 2004. URL: https://www.cmu.edu/biolphys/deserno/pdf/sphere equi.pdf
19. Бодрышев В.В., Абашев В.М., Тарасенко О.С., Миролюбова Т.И. Интенсивность изображения, как количественная характеристика параметров газового потока // Труды МАИ. 2016. № 88. URL: http://trudymai.ru/published.php?ID=70428
20. Billig F.S. Shock-wave shapes around spherical-and cylindrical-nosed bodies // Journal of Spacecraft and Rockets, 1967, vol. 4, issue 6, pp. 822 - 823.
21. Любимов А.Н., Русанов В.В. Течение газа около тупых тел. - М.: Наука, 1970. ч.1. - 287 с.
Meshless algorithm for supersonic inviscid gas flows calculating
Sposobin A.V.
Moscow Aviation Institute (National Research University), MAI, 4, Volokolamskoe shosse, Moscow, A-80, GSP-3, 125993, Russia
e-mail: spise@inbox.ru
Abstract
The presented work regards the problem solution of numerical modeling of blunted body flow-around by the supersonic inviscid gas flow [1, 2]. Such problem may be solved by rather wide specter of the well-known methods. Cartesian-grid-based immersed boundary method was successfully applied to solve this problem in the two-dimensional setting [3, 4]. It was employed to study the flow evolution in the shock layer while large-size particle movement along the symmetry axis was performed within the framework of a series of works [5-7]. The above said studies continuation while a particle movement along various trajectories required the problem solution of the gas flow modeling in the 3D setting.
Unlike the finite volume method, the meshless method does not divide the computational domain into mesh cells, but uses the cloud of points to represent the values of functions. The least square method is being used for spatial approximation of function derivatives. The article describes the Harten-Lax-van Leer method in combination with the MUSCL scheme and Van Albada limiter as applied to Euler system solution in 3D space. The least square method is being used to approximate Neuman boundary conditions on the boundary of the surface [17].
The meshless method selection was stipulated by the possibility of gas flow modeling in the domains with complex geometry. This method is less fastidious to the computer memory compared with those on the Cartesian meshes, employed in the previous works [47].
The meshless method for Euler equations solving was adapted also to the 2D space for the flat and axisymmetric cases [14].
The software implementation of the described method is accomplished in the C++ programming language and employs the OpenMP parallelization technology [15, 16].
The series of computational experiments on the bodies flow-around by the supersonic flow was performed to test the method operation accuracy. The article describes in detail the problem of modeling the sphere flow-around by the inviscid gas flow at the incident flow Mach number of 3. The algorithm for the nodes set generation in the 3D space for the given problem is described. A good agreement with solutions obtained by the finite volume methods on Cartesian meshes [3] and reference data [21] is demonstrated.
The next stage in the development of the meshless method is supposed to be solution of the Navier-Stokes equations in order to simulate supersonic viscous flows. Keywords: numerical simulation, meshless method, Euler equation, supersonic flows-around of bodies.
References
1. Reviznikov D.L., Sukharev T.Yu. Teplovye protsessy v tekhnike, 2018, vol. 10, no. 1 - 2, pp. 5 - 15.
2. Bykov L.V., Nikitin P.V., Pashkov O.A. Trudy MAI, 2014, no. 78. URL: http://trudymai.ru/eng/published.php?ID=53445
3. Vinnikov V.V., Reviznikov D.L. Trudy MAI, 2007, no. 27. URL: http://trudymai. ru/eng/publi shed.php? I D=34018
4. Vinnikov V.V., Reviznikov D.L., Sposobin A.V. Matematicheskoe modelirovanie, 2009, vol. 21, no. 12, pp. 89 - 102.
5. Reviznikov D.L., Sposobin A.V., Ivanov I.E. Teplofizika vysokikh temperature, 2018, vol. 56, no. 6, pp. 908 - 913. DOI: 10.31857/S004036440003569-9
6. Reviznikov D.L., Sposobin A.V., Ivanov I.E. Oscillatory flow regimes resulting from the shock layer - particle interaction, High Temperature, 2020, vol. 57, no. 2, pp. 278 - 283. DOI: 10.1134/S0018151X20020169
7. Sposobin A.V., Reviznikov D.L., Ivanov I.E., Kryukov I.A. Izvestiya vysshikh uchebnykh zavedenii. Aviatsionnaya tekhnika, 2020, no. 4, pp. 108 - 115.
8. Tolstykh A.I., Shirobokov D.A. Zhurnal vychislitel'noi matematiki i i matematicheskoi fiziki, 2005, vol. 45, no. 8, pp. 1498 - 1505.
9. Vasilyev A.N., Kolbin I.S., Reviznikov D.L. Meshfree computational algorithms based on normalized radial basis functions, International Symposium on Neural Networks, 2016, DOI: 10.1007/978-3-319-40663-3 67
10. Sattarzadeh S., Jahangirian A. 3D implicit mesh-less method for compressible flow calculations, Scientia Iranica, 2012, vol. 19, no. 3, pp. 503 - 512. DOI: 10.1016/j.scient.2012.04.013
11. Y. Wang, X. Cai, M. Zhang, X. Ma, D. Ren, J. Tan The study of the three-Dimensional meshless solver based on AUSM+-up and MUSCL scheme, International Conference on Electromechanical Control Technology and Transportation, 2015. DOI: 10.2991/icectt-15.2015.52
12. Drinfel'd G.I. Interpolirovanie i sposob naimen'shikh kvadratov (Interpolation and the method of least squares), Kiev, Vishcha shkola, 1984, 103 p.
13. Harten A., Lax P.D. and van Leer B. On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws, SIAM Review, 1983, vol. 25, no. 1, pp. 35 - 61.
14. Molchanov A.M. Matematicheskoe modelirovanie zadach gazodinamiki i teplomassoobmena (Mathematical Modeling of Gas Dynamics and Heat and Mass Transfer Problems), Moscow, Izd-vo MAI, 2013, 206 p.
15. Antonov A.S. Parallel'noe programmirovanie s ispol'zovaniem tekhnologii OpenMP (Parallel programming using OpenMP technology), Moscow, Izd-vo MGU, 2009, 79 p.
16. Malyavko A.A. Parallel'noe programmirovanie na osnove tekhnologii OpenMP, CUDA, OpenCL, MPI (Parallelprogramming based on OpenMP, MPI, CUDA technology), Moscow, Izd-vo Yurait, 2021, 135 p.
17. Hashemi M.Y., Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations, Journal of Computational and Applied Mathematics, 2011, no. 235, pp. 4687 - 4700. DOI: 10.1016/j.cam.2010.08.002
18. Deserno M. How to Generate Equidistributed Points on the Surface of a Sphere, 2004. URL: https://www.cmu.edu/biolphys/deserno/pdf/sphere equi.pdf
19. Bodryshev V.V., Abashev V.M., Tarasenko O.S., Mirolyubova T.I. Trudy MAI, 2016, no. 88. URL: http://trudymai.ru/eng/published.php?ID=70428
20. Billig F.S. Shock-wave shapes around spherical-and cylindrical-nosed bodies, Journal of Spacecraft and Rockets, 1967, vol. 4, issue 6, pp. 822 - 823.
21. Lyubimov A.N., Rusanov V.V. Techenie gaza okolo tupykh tel (Gas Flow Near Blunt Bodies), Moscow, Nauka, 1970, Part I, 287 p.
