______ УЧЕНЫЕ ЗАПИСКИ Ц АГ И
Том XVIII 1987
№ 2
УДК 629.735.33.015.3 : 533.695.12
ЧИСЛЕННЫЙ АЛГОРИТМ РАСЧЕТА СВЕРХЗВУКОВОГО НЕВЯЗКОГО ОБТЕКАНИЯ КОМБИНАЦИЙ КРЫЛА
С КОРПУСОМ
Ю. И. Лобановский
Рассматривается способ построения численного алгоритма, позволяющего в рамках уравнений Эйлера получать поля течений и интегральные аэродинамические характеристики основных схем планера сверх- и гипер-звуковых летательных аппаратов. Полученные решения сравниваются с другими численными и экспериментальными результатами.
1. Для исследования аэродинамики сверх- и гиперзвуковых летательных аппаратов необходимо применение численных методов, позволяющих получать подробную и многообразную информацию об изучаемых объектах. Для этих целей в настоящее время активно развиваются численные алгоритмы решения уравнений Эйлера, все более усложняющиеся с приближением формы поверхности рассчитываемых конфигураций к реальным компоновкам, увеличением точности результатов и выявлением локальных особенностей течения [1—8]. Создаются пакеты программ, включающие несколько последовательно работающих алгоритмов [9], что позволяет рассматривать тонкие эффекты, возникающие в поле течения [10].
Однако важнейшей задачей гиперзвуковой вычислительной аэродинамики является создание методов, позволяющих достаточно быстро и с достаточной степенью точности получать интегральные аэродинамические характеристики компоновок летательных аппаратов, а также поля течения в основных возмущенных зонах, что необходимо для анализа основных закономерностей гиперзвукового обтекания, а также для поиска оптимальных или близких к ним конфигураций. В попытках решения этой задачи предлагается, например, численный алгоритм приближенного решения уравнений Эйлера в рамках гиперзвуковой теории малых возмущений [11]. Для достижения указанной цели необходимо, чтобы алгоритм был работоспособен в широком диапазоне режимов течения и применим к различным, в том числе и нетрадиционным компоновкам, например, с различными носовыми наплывами. При этом требования к точности решения в отдельных локальных зонах течения могут быть менее жесткие, чем для пакетов алгоритмов.
Взаимодействие скачков уплотнения, создаваемых различными элементами компоновки, например, пересечение головного скачка и скач-
ка уплотнения, возникшего на передней кромке крыла, может привести к трудностям в работе универсального алгоритма с их выделением. А такой режим обтекания является типичным для практически интересных компоновок гиперзвуковых летательных аппаратов при Мсо = 4 ... 8. Дальнейшее усложнение геометрии компоновок, особенно в их носовой части, что наиболее перспективно для улучшения их аэродинамических характеристик, может увеличить эти трудности. Следовательно, несмотря на всю привлекательность выделения скачков уплотнения, для решения поставленной задачи представляется целесообразным применять принцип сквозного счета, который при повышении требований к точности определения локальных параметров во всем поле течения может быть заменен принципом выделения плавающих скачков [7].
Другой важнейшей особенностью численного алгоритма является способ формирования расчетной области из той физической области, в которой первоначально сформулирована краевая задача [12]. Уже традиционным стал способ выбора отображающей функции или системы функций, преобразующей исходную физическую область в расчетную в виде прямоугольного параллелепипеда, где и строится равномерная прямоугольная счетная сетка. По мере перехода от простых линейных отображений [13] к конформным [8] усложняется геометрия рассчитываемых конфигураций и возрастают возможности численных алгоритмов, однако в случае рассмотрения крылатых летательных аппаратов отображающие функции при таком способе отображения приводят к сильным деформациям физической области, и равным счетным ячейкам соответствуют очень сильно различающиеся по размерам ячейки физического пространства (см., например, [8]). На передней кромке бесконечно тонкого крыла такие преобразования имеют особенность [5].
Так как обычно используемые в рассматриваемом классе численных алгоритмов явные конечно-разностные схемы являются условно устойчивыми, максимально возможный шаг по гиперболической переменной, зависящий от поперечных размеров наименьших счетных ячеек, мал и в большей части поля течения далек от оптимального. Поэтому при таком способе формирования расчетной области количество счетных шагов по гиперболической переменной для конфигурации с развитыми тонкими поверхностями велико, а, следовательно, велико и время счета. Переход к абсолютно устойчивым (при малых возмущениях) неявным схемам [13] увеличивается сложность алгоритма, а вследствие резкого увеличения числа арифметических операций при расчете одной точки поля течения делает проблематичным заметное сокращение времени счета.
Таким образом, так как для решения поставленной задачи необходимо, чтобы время, затрачиваемое на единичный расчет, было сравнительно небольшим, нужно отказаться от отображений физического пространства на расчетную область, приводящих к сильным его деформациям. Следовательно, расчетная область в поперечном сечении должна представлять собой не четырехугольник, а явно содержать в себе образы основных элементов компоновки, аналогично тому, как это было осуществлено в алгоритме расчета обтекания крылатых конических конфигураций [14].
При этом можно уменьшить степень деформации физического пространства на порядок, соответственно снизив при фиксированном числе счетных точек в поперечном сечении на порядок затраты времени на расчет компоновки. Зоны потери аппроксимации конечно-разностной схемы в окрестности передних кромок крыла и ребер корпуса, составляющие 1—2 счетные ячейки, увеличиваются пропорционально росту
размера ячеек [2, 14], а точность решения в остальных областях поля течения даже несколько возрастает, вследствие существенного увеличения размера шага по гиперболической переменной и приближения его к местным оптимумам [15].
Большая точность и меньшая степень размазывания газодинамических разрывов в методах сквозного счета при использовании конечно-разностных схем второго порядка по сравнению со схемами первого порядка аппроксимации являются существенными преимуществами, которые естественно использовать при построении численного алгоритма. Неблагоприятные дисперсионные характеристики однородных схем второго порядка аппроксимации, которые могут привести к нефизическим колебаниям численного решения в окрестности газодинамических разрывов, значительно улучшаются введением неоднородного оператора, позволяющего монотонизировать решение [15].
2. Реализация численного алгоритма расчета сверхзвукового невязкого обтекания комбинаций крыла с корпусом, отвечающего сформулированным в п. 1 требованиям, осуществляется путем обобщения алгоритма расчета обтекания крылатых конических тел [14].
Пусть в цилиндрической системе координат поверхность корпуса летательного аппарата быть может за исключением одной линии описывается гладкой функцией г1п(х, <р), поверхность крыла нулевой толщины— частью плоскости Гу,{х, <р0), внешняя граница расчетной области, находящаяся вне зоны возмущенного течения, — гладкой функцией гех(х, ф) (рис. 1 ,а). В силу симметричности течения относительно вертикальной плоскости расчеты проводятся в полупространстве, ограниченном плоскостью симметрии.
и
Наряду с цилиндрической системой координат х, г, ф вводится неортогональная криволинейная система координат х, s(x, г, <р), <р, где
f>)^ <’> Фо— полярный угол плоскости крыла; х0 — абсцисса расположения наибольшей толщины корпуса при ф = ф0.
Если рассматривается наиболее важная с точки зрения практических приложений несимметричная комбинация крыла и корпуса (схемы низко- и высокоплана), то функция Г{П(х, ф) определена только при Ф<Фо, а при ф>фо она доопределяется следующим образом:
!'ех (•*, у) [гех (х0 , Уо) rin (.X, Уо — 0) — гех (.X, у0) гin (х0, у0 -Гех (х, Уо) Vex С*о > <Ро) — гт (хй , Уо — 0;J
■0)]
При <р><р0 Г1П(Х> ?)<0, и функция г1п(х, <р) имеет разрыв при ? — ?о- Однако, если
д гщ (х, уо + 0) а гт (х, уо — 0)
д у д у
=0,
что означает ортогональность поверхности корпуса плоскости крыла на линии их стыка, преобразование (1) является всюду непрерывным и гладким.
В системе криволинейных координат х, я, ф поверхность корпуса компоновки и внешняя граница расчетной области являются поверхностями 5 = сопз1:, а поверхность крыла — плоскостью ф = сопз1:. При х=х0 и ф = фо поперечные деформации физического пространства при его отображении отсутствуют, в других точках они остаются невелики.
Система дифференциальных уравнений газовой динамики стационарного течения совершенного невязкого газа в безразмерной дивергентной форме, записанная в неортогональной системе координат х, я, ф, когда за характерные величины приняты значения параметров невозмущенного потока У», рос, роо, имеет следующий вид:
л
дЕ
дх
л
dF
л
дй
л
Е-
Е\
ds
/» = Р[1+/(1
с д s
+ ~л7~ + "XT + Н — 0,
dy
U1
■V
®2)];
д s
Е + P-F + —G-
'or ' а у
л
G — G]
-^(~)e — ~(~)f-~(~) G;
д s \ д х ) д s \ д г / д s \ д у /
О
д s дх ри k р + р и2
р UV р UW
рда
р UW
pvw k р + р W* 1
d_l’ds д s Р и р UV kp -f р V2 р VW
0 0
kp — р Z02 pVW
-М2 •
* 100»
(2)
и, V, 1Ю — компоненты вектора скорости в цилиндрической системе координат; р — давление; р — плотность; Моо — число Маха невозмущенного потока; х — отношение удельных теплоемкостей.
Решение задачи осуществляется маршевым методом сквозного счета по гиперболической переменной х. В исходном сечении, совпадающим с острым носком конфигурации, задается невозмущенный поток. Для численного решения системы (2) используется конечно-разностная однородная схема Мак-Кормака второго порядка аппроксимации [16], дополненная неоднородным оператором монотонизации типа [15].
Рассчитываемая компоновка может быть крылом, корпусом, поперечное сечение которого линейно отображается на круг, корпусом с плоским низом или верхом, сечение которого отображается на полукруг, и их комбинацией. Для каждого из этих случаев преобразование (1) создает соответствующий образ компоновки в расчетной области. Общий вид сечения расчетной области при фиксированном значении х для низко- или высокоплана показан в плоскости 5, ц> на рис. 1,6. Область представляет собой многоугольный прямоугольник АВСйЕЕ с разрезом Ев, на который отображается крыло. Корпус отображается на прямоугольник РРЕБ. Узлы счетной сетки, лежащие слева от линии Ж2, которая определяется из условия г = О, являются мнимыми, так как соответствуют отрицательным значениям полярного радиуса г цилиндрической системы координат, и расчеты в них не производятся. При изменении переменной л: граница РС} смещается, возникают новые действительные расчетные точки, и при х = х0 линии О А и РС} совпадают, так что все счетные точки становятся действительными.
Если рассматривается изолированное крыло, отрезки ОА, РС} и ОЕ сливаются. Когда рассчитывается конфигурация в схеме среднеплана, граница ЛF становится продолжением линии БЕ, и, наконец, если рассматривается изолированный корпус, разрез ЕС исчезает. При изменении переменной х вид расчетной области, вообще говоря, меняется.
Вся расчетная область покрывается равномерной прямоугольной счетной сеткой. На линии Еб расчетные точки являются двойными и содержат информацию о параметрах потока как на наветренной, так и на подветренной сторонах крыла [2]. На линии ВС (внешняя граница) задаются параметры набегающего невозмущенного потока. Для моделирования условий непротекания на твердых поверхностях (РЕ, Ей, БЕ) и в плоскости симметрии (СИ и ВС?) вводятся дополнительные ряды счетных точек и используется принцип отражения [1].
Принцип отражения аппроксимирует условия непротекания со вторым порядком точности только на плоских поверхностях, поэтому для увеличения точности выполнения граничных условий на всех твердых границах используется метод элементарных волн сжатия — разрежения [1] в качестве дополнительного шага коррекции численного решения. Особенно сильно эта коррекция влияет на решение на первых двухтрех шагах по гиперболической переменной_х, когда на поверхности компоновки происходит резкая перестройка от невозмущенного течения к течению, удовлетворяющему системе уравнений (2) и граничным условиям.
Параметры поля в узлах счетной сетки, являющихся вершинами нестандартных по размеру ячеек, рассекаемых линией (рис. 1,6) во избежание возникновения неустойчивости численного решения [16] по полярному радиусу г, определяются с помощью интерполяции нормальной компоненты вектора скорости У„ по ближайшим узлам стандартных ячеек и линии ЯС1, на которой Уп = 0, и применением вслед за этим метода элементарных волн сжатия—разрежения.
Для удовлетворения условий устойчивости [16] по полярному углу Ф при малых 5, как и в [17] с уменьшением параметра в происходит последовательное удвоение счетных ячеек по переменной ф, что ограничивает уменьшение их линейных размеров в этом направлении.
В непосредственной окрестности передней кромки крыла, в области очень резкого изменения параметров потока по полярному углу ф, конечно-разностная модель граничных условий, как и конечно-разностная система уравнений, перестают аппроксимировать решение дифференциальной краевой задачи. Для того чтобы локализовать погрешности численного решения в узкой окрестности передней кромки крыла и уменьшить их величину, используя в соответствии с принципом отражения линейность нормальной компоненты скорости по нормальному расстоянию от поверхности крыла, величину этой скорости на слое ф0±Дф определим как
1'.(?. + 'Ч)-,!'^1Мл±2а?). (3)
Значения нормальной компоненты скорости из (3) используются для определения параметров течения в фиктивных отраженных точках при Ф=ф0::ьЛф. Вне рассматриваемой области течения подстановка (3) с точностью до членов второго порядка аппроксимации является тождеством, а в окрестности передней кромки крыла обеспечивает ликвидацию нежелательных явлений, связанных с отсутствием аппроксимации.
Если передняя кромка крыла не совпадает с узлом счетной сетки, решение в ближайшем к ней внутреннем (не лежащем на крыле) узле сетки является линейной комбинацией решения с граничными условиями на поверхности крыла и решения во внутренних точках поля течения с весовыми коэффициентами, соответственно пропорциональными частям стороны ячейки, занимаемым крылом и свободным от него. Эта процедура была отработана в алгоритме расчета конических течений.
На ребре корпуса, образованном пересечением его боковой и нижней поверхностей, как и на передней кромке крыла, расчетные точки являются двойными. В этих парах точек вычисляются нормали к соответствующим участкам поверхности корпуса. Скалярные параметры поля течения определяются на ребре корпуса путем интерполяции данных из ближайших узлов счетной сетки, а компоненты вектора скорости вычисляются по его модулю и нормалям для выполнения условий непротекания.
Все эти кратко рассмотренные процедуры, применяемые в нестандартных узлах счетной сетки, по существу представляют собой интерполяцию решения не более чем на удвоенном по сравнению с размером счетной ячейки интервале, дополненную в некоторых случаях использованием метода элементарных волн сжатия — расширения, что сохраняет второй порядок аппроксимации решения по крайней мере там, где аппроксимация существует. Согласование результатов расчета конической конфигурации по описываемому алгоритму и алгоритму [14], в котором нет нестандартных узлов, показано ниже (см. рис. 2,6).
Выше была сформулирована конечно-разностная краевая задача для компоновки с плоским крылом нулевой толщины. Введя дополнительную деформацию по полярному углу ф, можно рассматривать в строгой постановке и крылья с толщиной. Однако показано [18, 19, 12], что при численных расчетах в рамках уравнений Эйлера тонких крыльев сверх- и гиперзвуковых компоновок с высокими аэродинамическими характеристиками можно при МооТ'С! с приемлемой точностью сносить граничные условия с поверхности крыла на некоторую плоскость так,
как это обычно делается в методах линейной теории (т — характерный наклон поверхности крыла, вызванный его толщиной). Вследствие больших простоты и гибкости этого способа и его удобства в применении при исследовании влияния различных геометрических параметров крыла на аэродинамические характеристики компоновки, именно он и был реализован в описываемом алгоритме автором с М. Е. Нестеровым. В плоскости ср = (ро как на верхней, так и на нижней поверхности определяются все необходимые для формирования граничных условий поверхностные производные неплоского и обладающего толщиной крыла.
3. Сравнение распределений давления вдоль поверхности оживаль-но-цилиндрического тела вращения, полученных с помощью предлагаемого алгоритма (точки), с численными (линия) и экспериментальными (треугольники) данными из [9] при Мос = 4,05 и угле атаки а = 0 показано на рис. 2, а. Продольная координата х отнесена к длине ожи-вальной носовой части с удлинением л=4. Число Рейнольдса в эксперименте, вычисленное по максимальному диаметру тела вращения, составляет Рвоо = 2,8 * 106. Согласование полученных данных с результатами расчета по специализированной программе для осесимметричных течений хорошее, с экспериментом — удовлетворительное. Наибольшее расхождение расчетов и эксперимента наблюдается там, где экспериментальные точки имеют наибольший разброс.
На рис. 2, б сравниваются результаты расчета обтекания треугольной пластины с углом стреловидности %=79° при М<х, = 6 и а = 4°, полученные при помощи алгоритмов расчета конических конфигураций [14] (линии) и рассматриваемого (точки). Распределения давления по размаху крыла г хорошо согласуются между собой как на подветренной (1), так и на наветренной (2) его сторонах.
На рис. 3 показаны компоновки гиперзвуковых летательных аппаратов, интегральные характеристики которых представлены на рис. 4— 7. Все рассматриваемые конфигурации являются комбинациями крыла с корпусом в схеме низкоплана. Угол стреловидности передней кромки крыла компоновки 1 %=70°, относительная толщина параболического с плоским низом профиля крыла с = 3%. Корпус образован полови-
Компонсдка 1
2
/
Рис. 3
ной степенного тела с показателем 3/4, постепенно переходящего в полуцилиндр, а также клином, присоединенным к описываемому телу снизу.
Углы стреловидности передней кромки крыла компоновки 2 составляют до излома 80°, а после — 60°. Профиль крыла такой же, как и у конфигурации 1. Образующая носовой части корпуса является параболой, поперечное сечение его — полуэллипс с соотношением вертикальной и горизонтальной полуосей 1:0,61. Крыло заклинено по отношению к корпусу на угол 1,5°.
Компоновка 3 отличается от компоновки 2 отсутствием вертикального оперения и заклинения крыла, а компоновка 4 — кроме того, и крылом нулевой толщины. Отношение вертикальной и горизонтальной полуосей полуэллипса, образующего поперечное сечение корпуса компоновки 5 — 0,96: 1. Наплыв крыла почти постоянной ширины начинается у этой конфигурации от самой вершины корпуса, причем в стыке этого наплыва с корпусом имеются щели, длина которых составляет около половины длины наплыва. Крыло этой компоновки также нулевой толщины. Объемы корпуса конфигураций 2—5 одинаковы.
На рис. 4, 5 приведено сравнение интегральных характеристик компоновки 1 при М<х> = 4 и 7 и компоновки 2 при Моо=6, полученных численно и в экспериментах, проведенных в ЦАГИ. Определение суммарного коэффициента сопротивления трения Ср, составляющего в рассматриваемых условиях (числа Рейнольдса, определенные по длине модели, Кеоо~5 • 106-ь35-106) значительную часть полного коэффициента сопротивления сх компоновок, проводилось в соответствии с работой [20]. Согласование расчетных и экспериментальных данных для обеих компоновок во всем рассмотренном диапазоне параметров хорошее.
Сравнение интегральных волновых (с^ = 0) характеристик компоновок 3—5 при Моо = 6 показано на рис. 6, 7. Из численных данных следует, что влияние толщины крыла на зависимости су(а) компоновок рассматриваемого типа незначительно (кривые 3 и 4 на рис. 6). Носовой наплыв заметно увеличивает производную с*. Волновые поляры
(рис. 7) компоновок с нулевой толщиной крыла (кривая 4) и с крылом с = 3% (кривая.?) проходят практически эквидистантно. Наплыв и уменьшение относительной высоты корпуса уменьшают волновое сопротивление компоновки при фиксированных положительных значениях подъемной силы (кривые 4 и 5).
Время расчета рассмотренных компоновок на сетке с максимальным числом счетных точек в поперечном сечении 30X64 соответственно по координатам г и ф при использовании процессора с быстродействием -1,3- 10е операций в секунду составляет —45 минут. При исследовании групп конфигураций, носовая часть которых одинакова (например, при сравнении компоновок с различной толщиной крыла), до появления различий в геометрических размерах при маршевом движении по переменной х для всех них можно проводить единый расчет, начиная индивидуальные вычисления только с тех значений х, где начинаются эти различия. Такой подход приводит к сокращению суммарного времени расчета групп компоновок.
Проведенное совместно с В. В. Коваленко сравнение рассмотренного алгоритма и численного алгоритма [8], в котором применяются выделение головного скачка уплотнения и конформные отображения для формирования счетной области, показало, что при использовании по-
2 - «Ученые записки» № 2
17
следнего на счетной сетке с числом точек в 1,5 раза меньшим по обеим координатам число шагов по гиперболической переменной в описываемом алгоритме меньше в 12—15 раз, а процессорное время при пересчете на одинаковую производительность процессора меньше в 2,5 раза (при одинаковом числе счетных точек в сечении времена различались бы в 8—9 раз). При этом рассогласование интегральных аэродинамических характеристик компоновки 2 при Мех, = 6, полученных расчетом по сравниваемым алгоритмам, вызванное отчасти различной аппроксимацией формы поверхности компоновки, составляет не более 3—4%.
ЛИТЕРАТУРА
1. Kutler P., Lomax Н., Warming R. F. Computation of space shuttle flow using noncentered finite difference schemes. — AIAA Paper N 72—193, 1972.
2. Б а з ж и н А. П., Челышева И. Ф. О численном решении задачи обтекания плоского треугольного крыла сверхзвуковым потоком газа под малыми углами атаки.—Ученые записки ЦАГИ, 1974, т. 5, №5.
3. Косых А. П., М и н а й л о с А. Н. Расчет сверхзвукового течения у несущих тел и крыльев методом сквозного счета. — Труды ЦАГИ,
1977, вып. 1809.
4. 3 а р у б и н А. Г. Реализация конечно-разностного метода Край-
ко для расчета обтекания комбинации крыла и фюзеляжа сверхзвуковым
потоком идеального газа. — Труды ЦАГИ, 1978, вып. 1941.
5. К1 о р f е г G. Н., Nielsen J. N. Euler solution for wing and
wing-body combination at supersonic speeds with leading edge separa-
tion. — AIAA-80-0126, 1980.
6. А у к и н М. К., Тагиров Р. К. Метод расчета сверхзвукового обтекания летательного аппарата при наличии воздухозаборников, крыльев и оперения.—В кн.: Численные методы механики сплошной среды, Новосибирск, т. 11, № 6, 1980.
7. Ямамото И., Карасима К. Расчет трехмерных невязких сверхзвуковых течений методом плавающего скачка. — Ракетная техника и космонавтика, 1982, т. 20, № 2.
8. Коваленко В. В., Минайлос А. Н. Расчет невязкого сверхзвукового течения около комбинации крыло—корпус. — Труды ЦАГИ,
1984, вып. 2251,
9. Волков В. Ф., Жибинов С. Б., Колобов Б. П., Чернышева Р. Т., Я н е н к о Н. Н. Об опыте применения пакетной технологии для решения задач внешнего обтекания. — В кн.: Задачи аэродинамики тел пространственной конфигурации. — СО АН СССР, ИТПМ, Новосибирск, 1982.
10. Боровой В. Я., Голубинский А. А., Нерсесов Г. Г., П о х в а л и н с к ий С. М., Струминская И. В., Юмашев В. Л. Особенности гиперзвукового обтекания толстого крыла переменной стреловидности с затупленными кромками. — Ученые записки ЦАГИ, 1985, т. 16, № 1.
11. Воеводенко Н. В., Широносов В. А. Расчет гиперзвукового обтекания трехмерных тел с использованием закона плоских сечений. ■—Труды ЦАГИ, 1985, вып. 2262.
12. Wardlow А. В., Priolo F. J., Solomon J. М., В а 1-t а с i s F. P. Inviscid multiple zone calculations for supersonic tactical missiles. — AIAA Paper, 1984, 84—2099.
13. Бабенко К. И., Воскресенский Г. П., Любимов Л. И., Русанов В. В. Пространственное обтекание гладких тел идеальным газом. — М.: Наука, 1964.
14. Лобановский Ю. И. Расчет обтекания сверхзвуковым потоком невязкого газа крылатых конических тел. —• Ученые записки ЦАГИ,
1980, т. 11, № 6.
15. Лобановский Ю. И. О монотонизации конечно-разностных решений в методах сквозного счета. — Ж. вычисл. матем. и матем. физ.,
1979, т. 19, № 4.
16. Мае С or mack R. W. The effect of viscosity in hypervelocity impact cratering. — AIAA Paper N 69—354, 1969.
17. Иванов М. Я.,. Край ко А. Н. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. Ч. I. — Ж. вычисл. матем. и матем. физ., 1972, т. 2, № 2.
18. Хейз У. Д., Пробе тин Р. Ф. Теория гиперзвуковых течений. М.: Изд. иностр. лит-ры, 1962.
19. Wardlow А. В., Solomon J. М., Baltacis F. P. Supersonic inviscid flowfield computations of missile type bodies. — AIAA Paper 80—0271, 1980.
20. Гарбузов В. М., Колина Н. П., Пятнова А. И. Расчет коэффициентов сопротивления трения и теплоотдачи пластины и острого конуса, обтекаемого сверхзвуковым потоком, при турбулентном течении в пограничном слое. — Труды ЦАГИ, вып. 1881, 1977.
Рукопись поступила 4/XII 1985