Метод построения оптимальных контуров фюзеляжа и сопла летательного аппарата на режиме крейсерского сверхзвукового полета Текст научной статьи по специальности «Механика и машиностроение»

Том ХЫУ

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

2013

№ 3

УДК 533.6:621.452:629.73

МЕТОД ПОСТРОЕНИЯ ОПТИМАЛЬНЫХ КОНТУРОВ

ФЮЗЕЛЯЖА И СОПЛА ЛЕТАТЕЛЬНОГО АППАРАТА НА РЕЖИМЕ КРЕЙСЕРСКОГО СВЕРХЗВУКОВОГО ПОЛЕТА

А. П. МАЗУРОВ, С. А. ТАКОВИЦКИЙ

Представлен метод построения обводов элементов планера и сопла силовой установки высокоскоростного летательного аппарата. В качестве целевой функции выбрана подъемная сила на режиме крейсерского полета, максимизируемая при заданных ограничениях на внешние габариты и изопериметрическом условии на продольный момент. Использован прямой оптимизационный метод, объединяющий моделирование в рамках системы уравнений Эйлера и ньютоновский алгоритм минимизации. Сопоставлены результаты последовательной и комплексной оптимизации, усиливающей положительные эффекты интеграции. Выполнен анализ интегральных аэродинамических нагрузок для оптимальных конфигураций.

Ключевые слова: оптимизация, целевая функция, носовая часть фюзеляжа, сопло, подъемная сила, балансировка летательного аппарата.

ВВЕДЕНИЕ

Одна из основных проблем создания высокоскоростного летательного аппарата — интеграция планера и силовой установки. Общие вопросы интеграции на режиме крейсерского полета рассмотрены в [1, 2]. Возможно использование положительной аэродинамической интерференции при расположении воздушно-реактивного двигателя под фюзеляжем. Носовая и кормовая части фюзеляжа обеспечивают, соответственно, эффективный забор воздуха и расширение в реактивном сопле с косым срезом. На этих элементах аппарата сосредоточены значительные аэродинамические силы. Наряду с продольной составляющей эти силы имеют нормальную компоненту. При этом носовая и кормовая части расположены на большом расстоянии от центра масс. Эта особенность осложняет продольную балансировку летательного аппарата (ЛА). Известны решения проблемы, связанные с уменьшением неустойчивости посредством размещения балласта в носовой части фюзеляжа и установкой дополнительных носовых рулей [3]. В обоих случаях потери на балансировку оказываются значительными.

Комплексное исследование с учетом аэродинамической интерференции планера и силовой установки требует использования методов многодисциплинарной оптимизации [4]. Целью оптимизации является обеспечение высокой аэродинамической эффективности ЛА в целом, а не совершенство планера и силовой

МАЗУРОВ Анатолий Павлович

кандидат физико-математических наук, ведущий научный сотрудник ЦАГИ

ТАКОВИЦКИИ Сергей Александрович

доктор технических наук, начальник сектора ЦАГИ

установки в отдельности. При этом большое значение имеет корректный выбор целевой функции. В [5] рассмотрены два варианта оптимизации силовой установки с несимметричным соплом. В качестве оптимизируемого функционала принимались тяга и величина суммарной подъемной силы аппарата. Сопла, профилируемые по подъемной силе, получаются менее расширяющимися, что приводит к значительному возрастанию момента сил давления. При введении изопериметри-ческого условия на момент достигается компенсация его прироста при незначительном снижении подъемной силы [6]. Более общее решение получается при совместной оптимизации планера и силовой установки.

В данной работе представлен метод оптимального профилирования аэродинамической формы фюзеляжа и реактивного сопла с учетом уменьшения вклада носовой части фюзеляжа в создание подъемной силы и изменения угла наклона вектора тяги сопла относительно продольной оси.

Задача аэродинамического проектирования может быть сведена к условной максимизации (минимизации) функции многих переменных. Для рассматриваемого летательного аппарата, выполняющего крейсерский полет при заданных значениях числа Маха и высоты, в качестве целевой функции выбрана подъемная сила ЛА. Увеличение подъемной силы, а значит, и веса ЛА позволяет поместить дополнительное топливо и за счет этого увеличить дальность полета. Крейсерский полет требует выполнения условия равенства нулю тангенциальной силы, действующей на ЛА. Другое условие связано с продольной балансировкой и обеспечивается перераспределением аэродинамической нагрузки. Таким образом, можно сформулировать следующее математическое представление задачи:

Здесь У — подъемная сила, равная весу О ЛА и являющаяся целевой функцией; X — сила лобового сопротивления; Рх — компонента тяги по направлению полета ЛА. Изопериметриче-ское условие накладывается на момент тангажа Мг. Рассматриваются случаи с неактивным и активным ограничением на этот момент. В первом случае балансировка достигается размещением балласта в носовой части фюзеляжа, во втором — перераспределением аэродинамической нагрузки на ЛА. Геометрические ограничения определяются внешними габаритами. Вектор независимых переменных состоит из геометрических параметров И, ( =1... N). Наряду с геометрическими параметрами переменной величиной является угол атаки а. Допускается совместное профилирование фюзеляжа и сопла и профилирование этих элементов по отдельности.

Продольное сечение аппарата показано на рис. 1. Нижняя поверхность носовой части фюзеляжа формирует поле течения перед воздухозаборником и в процессе оптимизации не варьируется. Значение балансировочного угла атаки изменяется в достаточно узком диапазоне (примерно от 3 до 6°). Геометрические параметры ступеней торможения потока выбираются таким образом, что при среднем значении угла атаки скачки уплотнения попадают на кромку обечайки воздухозаборника. Длина силовой установки, высота воздухозаборника и камеры сгорания фиксируются.

1. ПОСТАНОВКА ЗАДАЧИ

У (к, а) = тах, X (, а)-Рх (, а) = 0, Мг (иИ■, а) = 0,

И, т:п < И, < И,

1 тш

г тах •

Варьируемый контур фюзеляжа

Варьируемый контур сопла

Рис. 1. Продольное сечение ЛА

В процессе оптимизации профилируются верхняя поверхность фюзеляжа и верхняя стенка сопла. Данный выбор не случаен. Носовая часть фюзеляжа вносит значительный вклад в суммарное лобовое сопротивление. При этом подъемная сила носовой части фюзеляжа превосходит подъемную силу крыла. С другой стороны, импульс струи на входе воздухозаборника в несколько раз превосходит суммарное лобовое сопротивление. Таким образом, возможно использование реактивной силы не только как движущей, но и как управляющей силы.

Расчет течения проводится в двумерной постановке. Аэродинамические характеристики рассчитываются в скоростной системе координат. Условный центр тяжести располагается на расстоянии 0.6 в долях длины фюзеляжа от передней кромки носовой части.

Оптимизируемые образующие фюзеляжа и сопла представляются набором отрезков, соединяющих узловые точки. Точки сгущены к передней кромке фюзеляжа и к горлу сопла и в продольном направлении не перемещаются (xi = const). Параметрами задачи являются ординаты узловых точек hi. Каждая образующая разбивается на 100 отрезков. Число геометрических параметров составляет N = 199. Фиксируются положение точек образующей фюзеляжа, находящихся на передней кромке и в плоскости донного среза, и точки образующей сопла, находящейся в плоскости горла. На геометрические параметры образующей фюзеляжа накладывается габаритное ограничение на максимальную высоту hi < hi max, где hi max = hф, hф — ордината образующей фюзеляжа в плоскости донного среза.

Два ограничения, фиксирующие длину и высоту входного сечения сверхзвукового сопла, учитываются при выборе системы геометрических параметров. Третье ограничение связано с заданием высоты выходного сечения сопла и определяется положением верхней поверхности фюзеляжа. Высота выходного сечения находится в процессе решения задачи. В рассмотренном случае ограничение сопла по степени расширения оказалось не активным. Дополнительное ограничение типа неравенства — условие, накладываемое на кривизну поверхности вблизи входного сечения сопла. Для точек образующей сопла, примыкающих к плоскости горла, задается ограничение на минимальное значение радиуса кривизны. Возможен также излом контура образующей, однако это неблагоприятно по тепловым характеристикам. Вводится окружность, проходящая через первую узловую точку, с центром, расположенным в плоскости входного сечения. Окружность находится вне внутреннего канала. Узловым точкам запрещается располагаться внутри данной окружности. Ограничение имеет следующую математическую формулировку:

hi < hc + R R2 -(xi - xc )2, xi - xc < R, где xc и hc — продольная координата и ордината для точки входного сечения сопла, R — радиус окружности.

2. МЕТОД ОПТИМИЗАЦИИ

Оптимальное профилирование элементов летательного аппарата выполняется эффективным прямым методом [7, 8]. Метод объединяет решение прямой задачи, связанной с расчетом обтекания аппарата и его аэродинамических характеристик, и решение экстремальной задачи, заключающейся в поиске вариаций формы, направленных на повышение аэродинамического совершенства. Оптимизационный процесс является циклическим. Моделирование выполняется в рамках системы уравнений Эйлера, поиск оптимума — на основе ньютоновского алгоритма. Для повышения эффективности исследований на основе сложных моделей течения требуется ускорение сходимости. Простейшие универсальные методы (такие, как метод покоординатного спуска и градиентный метод) не обеспечивают надежного решения плохо обусловленных задач. Даже при небольшом числе варьируемых параметров линейной скорости сходимости оказывается недостаточно. Линии (поверхности) равных значений целевой функции сильно вытянуты и вариации параметров, направленные на увеличение (уменьшение) целевой функции, оказываются малыми при относительно высокой погрешности их вычисления. В результате процесс оптимизации, осложняемый наличием ограничений, прерывается на значительном расстоянии от точки оптимума.

Значительное ускорение сходимости достигается путем упрощения постановки оптимизационной задачи на основе квадратичной аппроксимации целевой функции. Особенности аэродинамических функций исследуются в рамках локального анализа распределения аэродинамической нагрузки по поверхности ЛА, обеспечивая аналитическое представление целевой функции и

функций ограничений. На основе аппроксимаций матрицы Гессе и градиента находятся вариации формы, направленные на улучшение аэродинамических характеристик. Полученные вариации используются в численном расчете. Данный подход обеспечивает высокую скорость сходимости при большом числе независимых переменных.

3. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕЧЕНИЯ И МЕТОД РАСЧЕТА

Обтекание летательного аппарата и поле течения в сопле моделируется в двумерной постановке в рамках системы преобразованных уравнений Эйлера. Характеристики крыла определяются для изолированного варианта. Сопротивление поверхностного трения вычисляется по соотношениям для пластины.

При расчете тяги сопла определяющими газодинамическими параметрами задачи являются: располагаемая степень понижения давления в сопле пср = рос /рх, показатель адиабаты у и числа

Маха на входе в сопло М0 и во внешнем потоке М. Величины пср = рос/рх, у и М0 находятся

в результате термодинамического расчета камеры сгорания двигателя.

Детальный расчет параметров потока по тракту двигателя, начиная от входа в воздухозаборник и заканчивая выходом из сопла, представляет собой сложную многопараметрическую задачу и практически неприемлем для решения оптимизационных задач, в которых одновременно рассматриваются фюзеляж и сопло. В настоящей работе для расчета параметров потока на входе в сопло используется упрощенный подход, основанный на раздельном определении характеристик элементов двигателя (воздухозаборника, камеры сгорания и сопла) с помощью интегральных законов сохранения массы, импульса и энергии. Для учета термодинамических свойств газа используются эмпирические зависимости теплоемкости газа от температуры в виде полинома 7-й степени. Основные соотношения для термодинамических параметров и значения коэффициентов полинома берутся из работы [9].

Для приближенного расчета моментно-тяговых характеристик выходного устройства применяется методика, основанная на численном интегрировании двумерных стационарных уравнений Эйлера с использованием модифицированной разностной схемы Мак-Кормака.

При расчете течения в сопле расчетная область в плоскости х, у огранивается входной свободной границей слева, верхним контуром сопла уь (х) и нижней составной границей ун (х), одна часть которой является твердой стенкой обечайки, а другая — границей струи. Решение задачи находится в преобразованных независимых переменных < и п, связанных с декартовыми координатами х и у преобразованием:

Е = х,

у - Ун (х) Уь (х)- Ун (х)

п =

Уравнения Эйлера, описывающие стационарное течение невязкого газа в плоском сопле, записываются в преобразованных переменных <, П в следующей векторной форме:

ди А ди п

-+ А-= 0,

д< дп

(1)

Г Р 1 Г в Б 01

и = с , А = в в 0

V 5 V 10 0 А V

А = пх + сПу, в = пх +-

СП у

Б =

УП у

1 -

(а 2/и 2 ) 1 -(а 2/и 2 )

(

е=-М>

уи

1 -(а V и 2 )

П >

Р = 1пР,

Ро

V

с = —,

и

5 - 5о

5 =-

где а — скорость звука; и, V — проекции вектора скорости на оси декартовых координат х, у; Ро, р — полное и статическое давление; 5о, 5 — энтропия на входе в сопло и в струе; сЛ! — удельная теплоемкость при постоянном объеме.

Для численного решения уравнения (1) применяется модифицированная разностная схема Мак-Кормака, в которой матрица коэффициентов уравнения движения расщепляется на две матрицы в соответствии со знаками собственных значений. Такой подход используется, например, в работе [10], где показано преимущество схемы с ориентированными разностями по сравнению с симметричной схемой. Чтобы придать исходной схеме Мак-Кормака [11] свойство монотонности, матрица А в уравнении (1) представляется в виде суммы

А = А++А-;

где А+ и А — матрицы, имеющие соответственно только положительные Я+ и только отрицательные Я- собственные значения.

Для численного решения уравнения (1) модифицированная схема Мак-Кормака записывается в следующем виде:

Ок+1 = О* А+5О+ + А-5 Ор),

] ] ] + ] /'

ик+1 = 0.5

О

к+1

и-Дп(а+5О+ + А -5Ос

где Д^, Дп — шаги разностной сетки \Ъ,к = Д^(к -1), Пу = Дп(7 -1), к = 1...КЕ, у = 1... ^ соответственно в направлениях ^ и п.

Разности, входящие в эту систему уравнений, имеют вид [7]:

5О+Р = 2Оу - 30у_1 - О]_2, 5О_Р = О]+1 - О]

с

V

5О+ = О у - О у -1, 5 О- = 3Оу+1 - 2Оу - О у+2.

Применение ориентированных, с учетом знаков собственных значений, разностей делает схему монотонной без введения дополнительных сглаживающих операторов. Расчет производится маршевым методом вниз по течению от входного сечения сопла до его среза.

Чтобы упростить задачу и не рассчитывать внешнее поле течения при определении границы струи, распределение статического давления на ней вычисляется приближенно в зависимости от угла наклона струи с использованием линейной теории сверхзвуковых течений. При этом давление на кромке обечайки принимается равным давлению за косым скачком уплотнения, который возникает при взаимодействии струи с внешним потоком.

Приведенная методика используется также при расчете параметров течения на верхней поверхности фюзеляжа ЛА. Параметры потока на нижней поверхности фюзеляжа определяются по теории плоских сверхзвуковых течений. В результате совместного расчета течения в сопле и около фюзеляжа вычисляются значения суммарных аэродинамических сил и момента, действующих на ЛА. При этом донное давление на срезе фюзеляжа принимается равным давлению в невозмущенном потоке. Оптимальные обводы фюзеляжа и сопла находятся при помощи итерационного процесса, в котором каждая итерация состоит из двух этапов. На первом этапе решается вариационная задача, в которой находятся вариации контуров сопла и верхней поверхности фюзеляжа. На втором этапе по этим вариациям задается новая геометрия сопла и фюзеляжа и

проводится численный расчет, в результате которого находится новое распределение давления на поверхности ЛА и сопла, которое затем используется на первом этапе следующей итерации.

4. РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

В результате численной оптимизации построены аэродинамические конфигурации элементов ЛА для четырех вариантов постановки задачи. В первом случае фюзеляж и сопло проектируются совместно при неактивном ограничении на момент тангажа. Предполагается, что балансировка аппарата обеспечивается добавлением балласта. Второй случай представляет комплексную оптимизацию образующих фюзеляжа и сопла при изопериметрическом ограничении на момент. Кроме этого, по отдельности определяются формы фюзеляжа и сопла, обеспечивающие продольную балансировку летательного аппарата. При этом форма элемента, не участвующего в оптимизации, задается по результатам первой постановки задачи. Для трех последних вариантов постановки задачи ограничение на момент тангажа является активным и размещение дополнительного балласта не требуется. Балансировка достигается перераспределением аэродинамической нагрузки по поверхности летательного аппарата в результате изменения аэродинамической формы и угла атаки.

Образующие верхней поверхности фюзеляжа показаны на рис. 2. При неактивном ограничении на продольный момент фюзеляж имеет клиновидную верхнюю поверхность. В этом случае обеспечивается уменьшение сопротивления и увеличение подъемной силы. Задание условия самобалансировки летательного аппарата изменяет форму фюзеляжа. Верхняя поверхность становится выпуклой. При этом становится активным габаритное ограничение — появляется горизонтальный участок в кормовой части фюзеляжа. В случае отдельного профилирования фюзеляжа кривизна оптимальной образующей возрастает и увеличивается длина плоского элемента поверхности.

Оптимальные образующие стенки сверхзвукового сопла представлены на рис. 3. Наибольшую степень расширения имеет сопло, спроектированное при отсутствии ограничения по моменту тангажа. Требование самобалансировки аппарата приводит к увеличению кривизны направляющей и уменьшению высоты сопла в плоскости донного среза. Наиболее сильно данное свойство проявляется при оптимизации сопла при фиксированной форме фюзеляжа. Для всех рассмотренных вариантов задачи ограничение по габаритному размеру оказывается неактивным.

Сопоставление полей течения около оптимизируемых элементов ЛА для четырех вариантов постановки задачи представлено на рис. 4. Показаны линии равных значений числа Маха в плоскости симметрии с шагом ДМ = 0.15 и начальным значением М = 1.1. Течение около фюзеляжа

0 0.25 0.5 0,75 * I

Рис. 2. Оптимальные образующие верхней поверхности фюзеляжа:

1 — оптимизация фюзеляжа с активным ограничением на момент тангажа; 2 — оптимизация фюзеляжа и сопла с активным ограничением на момент тангажа; 3 — оптимизация фюзеляжа и сопла с неактивным ограничением на момент тангажа

V

1 - - — - _ _]

0.75 0.8 0.85 0.9 0.95 * [

Рис. 3. Оптимальные образующие сопла:

1 — оптимизация фюзеляжа и сопла с неактивным ограничением на момент тангажа; 2 — оптимизация фюзеляжа и сопла с активным ограничением на момент тангажа; 3 — оптимизация сопла с активным ограничением на момент тангажа

Рис. 4. Распределение чисел Маха около оптимизируемых элементов ЛА

зависит как от формы поверхности, так и от угла атаки ЛА. Для первого и четвертого вариантов постановки задачи давление на подветренной стороне практически постоянно. Отличие заключается в уменьшении угла атаки при оптимизации сопла. Оптимизация фюзеляжа при активном условии на продольный момент (второй и третий варианты) приводит к увеличению нагрузки в носовой части верхней поверхности фюзеляжа. Влияние угла атаки на течение в области сопла более слабое, чем для фюзеляжа. Это подтверждается сопоставлением полей распределения числа М для первого и третьего вариантов компоновки, имеющих одинаковую форму верхней стенки сопла и обтекаемых при разных значениях угла атаки.

Наибольшие значения целевой функции — подъемной силы ЛА — получены при оптимизации без ограничения на продольный момент. В этом случае для балансировки аппарата необходимо поместить в носовой части фюзеляжа балласт, вес которого составляет ~20% подъемной силы. При этом полетное значение угла атаки составляет а = 5.9°. При задании условия на самобалансировку аппарата получены более выпуклые аэродинамические формы контура фюзеляжа. В результате, при фиксированном угле атаки, увеличивается лобовое сопротивление планера и уменьшается тяга двигателя. Балансировочные потери проявляются в уменьшении подъемной силы на ~20%. Таким образом, оба варианта аэродинамической компоновки близки по весовым характеристикам размещаемого топлива и полезной нагрузки. Сопоставление летательных аппаратов по весам с выделением веса балласта показано на рис. 5. Преимущество самобалансирующегося ЛА проявляется в увеличении полезного объема и уменьшении полетного угла атаки до а = 4.5°, что приводит к уменьшению секундного расхода топлива примерно на 8%.

Оптимизация фюзеляжа и сопла по отдельности оказывается существенно менее эффективной. Оба этих варианта проигрывают по интегральным характеристикам рассмотренным выше компоновкам. Наименьшее значение подъемной силы получается при профилировании фюзеляжа: на 22% меньше, чем в случае комплексной оптимизации. При этом угол атаки незначительно

I 2 3

Рис. 5. Относительный вес ЛА:

1 — оптимизация фюзеляжа и сопла с неактивным ограничением на момент тангажа; 2 — оптимизация фюзеляжа и сопла с активным ограничением на момент тангажа; 3 — оптимизация фюзеляжа с активным ограничением на момент тангажа; 4 — оптимизация сопла с активным ограничением на момент тангажа

Рис. 6. Сопоставление по критерию m = (G - G6ajl)/qt (обозначения те же, что на рис. 5)

увеличивается, до а = 4.6°. Построение сопла при фиксированной форме фюзеляжа приводит к уменьшению подъемной силы на 17% при сравнении с комплексной оптимизацией. В этом варианте получено наименьшее значение полетного угла атаки а = 2.9°, что объясняется увеличением вклада реактивного сопла в создание подъемной силы.

При сопоставлении полезного веса найденных оптимальных аэродинамических компоновок варианты комплексной оптимизации при активном и неактивном ограничении на момент тангажа имеют близкие характеристики. В то же время, компоновка, самобалансирующаяся в продольном движении, обеспечивает более эффективный расход топлива и, соответственно, имеет наибольшую дальность полета. В качестве примера на рис. 6 дано сопоставление компоновок летательного аппарата по критерию т = (О - Gбал)/, где G — вес летательного аппарата, Gбал — вес балласта, д( — секундный расход топлива в двигателе. В этом случае проявляется эффективность оптимизации сопла. Изолированное профилирование сопла обеспечивает характеристики, близкие к результатам, полученным комплексной оптимизацией при неактивном ограничении на продольный момент.

По результатам исследования можно заключить следующее. Традиционный подход, основанный на раздельном рассмотрении внешней и внутренней аэродинамики, не обеспечивает корректного решения задач аэродинамического проектирования высокоскоростных летательных аппаратов. Представлен метод комплексной оптимизации элементов планера и силовой установки. Эффективность метода продемонстрирована на примере тестовой компоновки. Отмечена важность корректного выбора целевой функции.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (10-01-00208а).

ЛИТЕРАТУРА

1. Гусев В. Н. Интеграция планера гиперзвукового летательного аппарата с воздушно-реактивным двигателем // Ученые записки ЦАГИ. 1991. Т. XXII, № 5, с. 3 —11.

2. Guse v V. N. Optimal aerodynamic shapes of hypersonic vehicle with airbreathing engine // Development in high-speed-vehicle propulsion systems. V. 165, Progress in astronautics and aeronautics, AIAA, 1996, p. 17—49.

3. Heppenheimer T. A. Facing the heat barrier: a history of hypersonics // NASA history series. NASA SP-2007-4232. — Washington, DC, 2007.

4. Bowcutt K. Multidisciplinary optimization of airbreathing hypersonic vehicles // J. of Propulsion and Power. 2001. 17(6), p. 1184—1190.

5. Бафталовский С. В., Крайко А. Н., Макаров В. Е., Тилляева Н. И. Оптимизация силовой установки гиперзвукового летательного аппарата с прямоточным воздушно-реактивным двигателем // Изв. РАН. МЖГ. 1997. № 4, с. 127—135.

6. Миско Г. Ю. Построение оптимального сопла гиперзвукового летательного аппарата при заданных габаритах и моменте // Изв. РАН. МЖГ. 1999. № 1, с. 118—124.

7. Таковицкий С. А. Оптимальные формы треугольного крыла в сверхзвуковых потоках идеального и вязкого газов // Ученые записки ЦАГИ. 1998. Т. XXIX, № 3—4, с. 61—69.

8. Takovitskii S. A. Direct optimization method and aerodynamic shape design at supersonic flight conditions // Proc. 27th Intern. Congr. of Aeronaut. Sci. (ICAS). — Nice, France, 2010, p. 110.1 — 110.8.

9. Янкин В. И. Система программ для расчета характеристик ВРД на ЭЦВМ. — М.: Машиностроение, 1976, 168 с.

10. Daywitt J. E., Szostowski D. J., Anderson D. A. Split coefficient / locally monotonic scheme for multishocked supersonic flow // AIAA Paper 82-0287, 1982.

11. MacCormack R. W. The effect of viscosity in hypervelocity impact cratering // AIAA Paper 69-354, 1969.

Рукопись поступила 2/III2012 г.