Течение в окрестности точки излома передней кромки тонкого крыла на режиме сильного взаимодействия

Том XЬїї

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

№ 2

УДК 532.526.2

ТЕЧЕНИЕ В ОКРЕСТНОСТИ ТОЧКИ ИЗЛОМА ПЕРЕДНЕЙ КРОМКИ ТОНКОГО КРЫЛА НА РЕЖИМЕ СИЛЬНОГО ВЗАИМОДЕЙСТВИЯ

Г. Н. ДУДИН, А. В. ЛЕДОВСКИЙ

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

Ключевые слова: пограничный слой, крыло, гиперзвуковое течение, сильное вязко-невязкое взаимодействие.

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

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

ДУДИН еоргий Николаевич

доктор физико-латематических наук, главный научный сотрудник ЦАГИ, зав. кафедрой огидромеханики МФТИ

ЛЕДОВСКИИ Алексей Вячеславович

студент МФТИ

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

1. Рассматривается обтекание полубесконечной пластины с изломом передней кромки под нулевым углом атаки гиперзвуковым потоком вязкого газа на режиме сильного вязко-невязкого взаимодействия. Газ предполагается совершенным с постоянным показателем адиабаты у. Температура поверхности считается заранее известной и постоянной по всей поверхности крыла. Задача является в общем случае несимметричной из-за наличия излома передней кромки и угла скольжения. Используется двухслойная схема течения — внешняя область течения невязкого газа, описываемая уравнениями Эйлера, и вязкий ламинарный пограничный слой, описываемый уравнениями Прандтля. На рис. 1 представлена схема крыла и цилиндрическая система координат г, 0, у с центром в точке излома. Угловая координата 9 отсчитывается от биссектрисы О А угла излома. Ось у направлена по нормали к поверхности крыла. Основные геометрические параметры задачи: Р — угол между направлением набегающего потока и биссектрисой угла раствора крыла, 0 — угол между биссектрисой ОА и передней кромкой пластины (полуугол раствора крыла). В цилиндрической системе координат уравнения пограничного слоя для сжимаемого газа имеют вид [1]:

1 5рw i 5pv ри _ q

дг г 50

ду

ди w ди ди w2

и------1----------ь V--------------

дг г 50 ду г

dw w dw dw wu

її-----------1---------------h v---------1---------

dr r 50 dy r

1 dp 1 5

pr 50 p dv

dH wdH дН 1 5

dr r 59 dy p dy

1 dH 1-Pr d

Pr dy 2Pr

dy

H =

y — 0: u — v — w — 0, H-Hw,

Рис. 1. Схема крыла и система координат

у —» со : и = Um cos 9 - Р , w = -Um sin 9 - Р ,

Н = иЦ2.

Компоненты вектора скорости и, w, v направлены вдоль г, 9, у соответственно; р — плотность; ц — динамический коэффициент вязкости; Pr — число Прандтля; H — полная удельная энтальпия; Hw — энтальпия на поверхности крыла; С/да — скорость набегающего потока газа. Давление и температура газа в невозмущенном потоке стремятся к нулю, когда число Маха Мю —> со.

Для решения краевой задачи необходимо знать распределение давления, которое не задано

г

и должно определяться из совместного решения внутренней и внешней задачи. Для простоты используется приближенная формула «касательного клина» [1, 5] в форме, справедливой при условии Мда5 »1, где 5 — характерная безразмерная толщина пограничного слоя:

У + 1 7-7-2

Р = —^~ Роо^оо

58 81П е-р 55е

сое 0-В —---------------------------------------

дг г 50

В соответствии с геометрическими параметрами задачи и характерными оценками параметров течения для пограничного слоя в гиперзвуковом потоке [1, 5] вводятся безразмерные переменные:

г = Ьг°, 0 = 00°, у = Ьду°, Р = ©Р°,

и = и.лі\ м> = и~м>°.

у = их 50-4°, Н = иІН°/2,

<?2, 0 гт2с'2 0 Ос?

Р = Роо5 Р , р = ртию8 р , ц = , Ъе=ЬЪЪе.

Здесь Ь — некоторый характерный линейный размер; ц0 — динамический коэффициент вязкости при температуре торможения Т0.

Выполняя преобразование Дородницына [6], переходим к новым переменным:

Х = |р °сіу0, Уд =р%° +

о , * о дГк

' ~7Г—о+0М —

г° 50° дг

о '

В результате система уравнений (1) с учетом (3) принимает вид:

0

дг°и° дм>° дг°у^

аои

дХ

= 0,

©її

0 0 д о

ди ди

----------1----------------Ь V?

з 0 0 дг.0 о

дг г 50

о

..02

Г' _ 0 др° + д

(

ц°р°

дХ

о дм>° м>° дм?0 дм>° м>°и°

&и -------7Г + —7,------+ ------+ -------7Г~

..0 ..0 0 0 ..О

дги Г 50

дХ

гО

1 др° д - +

г

рV Э0° ах

оо ^ М- Р

о А

дХ

Л 0 дН° м>° 5Я° 5Яи 5 I о о

&и —— + —-----------— + Ух------= —ц р

дг

о „о

г° 50°

дХ дХ

1 5Я° 1-Рг д

Рг дХ

Рг

ах

Я° = -^—^— + и02 +м>02,

У _ 1 р

л /л 0 0 /~\ т т 0 т т О

л = 0: и =у§=м/ = 0, Н =Нм;,

X —> оо: и

о

сое

0

м/

-эш

0

Я —И.

Для динамического коэффициента вязкости предполагается степенная зависимость вязкости

О 0 02 02 ^

от температуры \х = Н —и — м> . Используя преобразование (3), для толщины вытеснения пограничного слоя получаем:

У 1

5" = ^ [ я°-и02 -и-02 ах.

2УР°

2

0

о

Распределение индуцированного давления (2) в переменных (3) принимает вид:

I 2

О У + 1 I

р =-------ІСОЄ

2

0 0° - р°

55° 1 .

п :8111

дг г°&

0 0° - р°

55

50

о

о

Для численного решения краевой задачи (4) — (6) необходимо учесть особенности поведения функций течения вблизи передних кромок и точки излома. Так как рассматривается обтекание полубесконечной пластины, то можно ввести автомодельные переменные [3]:

. * 0 * 0 3/4 0/~. О ах

X = г X, V =у, г +г &и —-

дг°

О 0 _1/2 * пО 0 0 _1/2 * пО Л * еО 0 3/4 аО

р = г р Ь , р = г р 0 , А , Ое = г ое0

0 о * о *

Предполагая, что функции не зависят от радиальной координаты г , т. с. и —и 0 . X ,

0 И« 0 И« 0 И« 0 И«

м/ =м/ 0 , X , Н =Н 0 , X , р = р X ,5=5 X , получаем из (4) двумерную систему уравнений:

ду»' ду* 5 *

—- +---------+ -&и =0,

дв° дХ* 4

* ди * ди „02 Р 5 м> —- + V-------------------м> = +-

50 дХ* 2р* дХ*

о * ди

ц р

дХ*

* дм> * ды оо I др 5

—- + у -----+ 0»т=--------^-7- +---

50° 5Х* р* 50° 5Х*

” “Г(Г + у

* * 5Я* 5

50 дХ дХ

1 дН* 1-Рг

о * ды

Iі Р —;

ч дк ,

*2 . *2 о и + м?

Рг дХ* Рг

зх*

Данная система является автомодельной для полубесконечного крыла и описывает течение в точке излома крыла г° = 0 .

Далее для упрощения предполагается, что вязкость линейно зависит от температуры со = 1 . На кромках крыла значение давления стремится к бесконечности, поэтому для численного решения краевой задачи необходимо учесть особенности поведения давления р* и толщины

вытеснения пограничного слоя 8* в окрестности передних кромок при 0° =+1. Для этого вводятся новые переменные:

Л *

2у

.У“1

4\-

л02

1-

л02

3/4

V* =У

1^1 + ,_еИ

I 2у р* 50°

р'=р% 1-0О2~1/2, 8:= 1-002 3/4А.

В результате преобразований (9) система уравнений (8) принимает вид:

5к и-0©0 1-

\02

5г| 2р Р*

(я, 0 ^

ды 5 о —- + -0м 50° 4

1-

\02 f а..О

Л

WO^_0WO2 V 50° ,

+v^ = 0lz^lzl я°+Цм^Л

I Эт1 У

дц р* 4у

1-

s02

„00 w ——+ 0М w

90

y-1 o 02 02

-v*-------=-----------H -u -w

dr\ 2y

( qO

0 1-0 дрл

+

02 Л я Z' д.,,0 Л

A pl 90°

+ -

V аЧ /

1-002 0 дН° 9Я° 9

-w

+ V*

N

1 9Я° 1-Рг'

^02+w023í

N = -и02-w02

Рг дг\ Рг 9г|

^^1

г) = 0: u°=v*=w°= О, Я°=Я°,

о

г| —> со: и —> eos

© 0° -р°

w —> — sin

© 0° -(3°

Выражения для толщины вытеснения (5) и давления (6) с учетом (7) и (9) принимают окончательный вид:

А = ——- Í Н° -и02 -w02 dr\ P*i 2у j

A=I±lJ 1-002 —A eos 2 4

0 0o -р°

02

На кромках крыла множитель 1 — 0 обращается в ноль, и система уравнений в частных

производных (10) вырождается в системы обыкновенных дифференциальных уравнений. Их интегрирование позволяет сформулировать краевые условия на передних кромках, которые в дальнейшем используются для решения двумерной задачи.

2. Система уравнений пограничного слоя (10) решается конечно-разностным методом [1]. При этом для аппроксимации исходных уравнений применяются конечные разности второго порядка точности с учетом направления течения. Сначала решаются уравнения на кромках крыла

при 0° =±1. Затем полученные профили скоростей и энтальпии используются в качестве граничных условий для двумерной задачи. В расчетах применяется следующая итерационная процедура. Пусть в начале п-й итерации задано распределение давления рп 0° . Тогда система уравнений пограничного слоя (10), записанная в линеаризованном разностном виде, решается методом прогонки. Рассчитанные при этом поля скоростей и энтальпии после проведения их релаксации позволяют определить толщину вытеснения А” 0° (11), которая вычислялась с помощью

формулы Симпсона четвертого порядка точности, и рассчитать новое распределение давления рп0 0° по формуле (12). Для определения давления на следующей итерации используется соотношение:

1 .

н—sin 0

0 0°-Э°

30°Д_ i.002 5А

2

90

„и+1 Г\0

„и г\0

reí р  Ар

и п0

где п — номер итерации; reí р — коэффициент релаксации давления; Ар" 0° — поправка к давлению, которая определяется из решения дифференциального уравнения второго порядка [4]:

2

i2Apn 0°

-aApn 0o = a pn 0o -p0n 0o

Apn 0°=±1 =0,

где а — положительная константа, которая в расчетах в зависимости от определяющих параметров краевой задачи изменяется от 0.5 до 2.

В конечно-разностной записи уравнение (13) принимает вид:

А0

где к — номер узла сетки по координате 0 {к = 1, 2,...); А0 — шаг сетки. Для решения уравнения (14) используется метод прогонки:

&Рк-1 = ^Рк^к-1 + В к-1, Арт = 0,

где Ак_х и Вк_х — коэффициенты прогонки, определяемые по формулам:

Л -

1

2 + аА02 - А

Вк =

а Рк-Рок А0 ~вк-1

lk-1

2 + аД02 - А,

Aq-B0 - 0.

-к-\

В отличие от [3] в данном методе для определения давления на очередной итерации используются заданные и полученные значения давления во всех точках расчетной области, включая граничные. Применение этого подхода позволяет выполнять расчеты на достаточно мелких сетках с высокой скоростью, что было почти невозможно при использовании «локальной» релаксации [3]. В случае [3] коэффициенты релаксации, в зависимости от варианта расчета, приходится уменьшать до 0.001, тогда как с использованием дифференциального уравнения (13) достаточно использовать коэффициенты релаксации порядка 0.5—0.8, что позволяет значительно ускорить расчеты.

Для верификации решения задачи, полученного методом прогонки, использовался метод установления с явной схемой по времени [7]. В данном подходе решается соответствующая нестационарная задача, которая при t —> со сходится к стационарному решению. Результаты, полученные этим способом, хорошо согласуются с результатами при использовании прогонки. Хотя этот метод проще в реализации, суммарное время расчетов оказывается в несколько раз больше из-за ограничения на шаг по времени, и поэтому для получения окончательных результатов используется метод прогонки.

Представленные в работе результаты получены на равномерной сетке со следующими параметрами: число узлов в вертикальном направлении — 400 с шагом сетки 0.02; число узлов по угловой координате — 50 с шагом 0.04. Расчеты показали, что использование сетки с числом узлов в 2 раза больше по каждому направлению приводит к изменению результатов менее, чем на 0.1%. Это свидетельствует о том, что выбранная сетка является достаточно подробной для данной задачи. Сходимость решения проверялась по давлению (наиболее медленно сходящейся

величине) согласно условию рт 11 - рт Jрт < reí р  10 5, где т — номер итерации.

3. Для численного решения используются преобразованные уравнения (10), поэтому удобно представить результаты в переменных:

р* Г ди ' 2 + W ' dw*' 2 р* дН* * _ р*

1-е»2 314 i { дц] { J ? W н i-©02 3^4 ‘ Vl-e02

8,= 1-

л02

3/4

А,

где хи --коэффициент поверхностного напряжения трения, а Хц — тепловой поток.

2

\ =

Для исследования влияния формы крыла в плане были проведены расчеты при температурном факторе Нп. =0.5, показателе адиабаты у = 1.4 и числе Прандтля Рг = 0.72 для следующих геометрических параметров крыла: 1) 0 = 45°, Р = 15°; 2) 0 = 67.5°, р = 22.5°; 3) 0 = 112.5°, Р = 22.5°; 4) 0 = 105°, р = 35°; 5) 0 = 135°, Р = 15°. На рис. 2 представлены распределения коэффициента напряжения трения по угловой координате 9. На крыльях с углом ©>90° (кривые 3—5) возникают локальные области повышенного трения с максимумом в направлении набегающего потока (в области линии стекания). На треугольных крыльях (кривые 1, 2) такое не наблюдается. Коэффициент напряжения трения выше вблизи кромки с меньшей стреловидностью. Распределение тепловых потоков по крылу качественно совпадает с распределением коэффициента поверхностного трения (рис. 3). При этом максимум в области локального повышения тепловых потоков более выраженный, чем для трения. Так для случая 0 = 135°, Р = 15° (кривая 5) тепловые потоки в этой области увеличиваются приблизительно в 2 раза. На рис. 4 представлено распределение давления по координате 9 для данных крыльев. Для треугольных крыльев (кривые 1, 2) давление в 2—3 раза выше, чем для крыльев с изломом (кривые 3—5). В области линии стекания для треугольных крыльев наблюдается увеличение давления (см. рис. 4) и толщины вытеснения (рис. 5). Для крыльев с изломом передней кромки (кривые 3—5) заметно уменьшение толщины вытеснения (см. рис. 5) в области повышенного трения, а на давление эта область практически не оказывает влияния.

Исследования влияния температурного фактора проводились для крыла с параметрами

0 = 105°, Р = 35° при Нн, =0.1, 0.3, 0.5, 0.7, 0.9, 1

5

0.15 -

0.1 -0.05 -

0 —*—*———'——*——*—‘——*—*—*—1—*——*——1 0

-1 -0.5 0 0.5 1

Рис. 2. Распределение коэффициента хи по угловой координате для крыла с Нм, = 0.5 и геометрическими параметрами:

1 — 0 = 45°. Р = 15°; 2 — 0 = 67.5°. Р = 22.5°; 3 — 0 = 112.5°. Р = 22.5°; 4 — в = 105°. Р = 35°; 5 — 0 = 135°. Р = 15°

Рис. 3. Распределение коэффициента тн по угловой координате 9 с теми же параметрами, что на рис. 2

1

0.8

0.6

0.4

0.2

1 \ 1 \ » \ 1 \ 1 \ 1 х 1 \ \ \ % \ \ 1 \ %

ч ^ ^ 1 -% ч ч ч ч ч ч ^ ^ т / 1 \ / 1 \ / 1 х / » X / 1 X ✓ 1 —^ / 1 / ! 9 / /

^ /• . / Л

; г 1 V

Ч \ з ч **’««1млтш: Ч 4 V У

5

_____I_____________I_______________I______________I____

_____I____________________________I____

___I___________I___________I___________1____

Л

-1

-0.5

0.5

Рис. 4. Распределение индуцированного давления с теми же параметрами,

что на рис. 2

Рис. 5. Распределение толщины вытеснения пограничного слоя на крыльях с теми же параметрами, что на рис. 2

Рис. 6. Распределение коэффициента тм для крыла с геометрическими параметрами © = 105°, (3 = 35° при:

1 — Н - 0.1; 2 — Я,„ = 0.3; 5 — Я = 0.5; 4 — Я,„ = 0.7; 5 — Я,„ = 0.9;

\У ~ \У ~ И'“ ~ И'"

0 I I .......................................... , I . . , . I 0

-1 -0.5 0 0.5 1

Рис. 7. Распределение коэффициента тя для крыла с геометрическими параметрами 0 = 105°, (3 = 35° при:

1 — Я,„ =0.1; 2 — Я,„ = 0.3; 3 — Я,„ = 0.5; 4 — Я,„ = 0.7; 5 — Я,„ = 0.8

-1 -0.5 0 0.5 1

Рис. 8. Распределение толщины вытеснения (сплошные кривые) и изолинии М = 1 (штриховые) для крыла с геометрическими параметрами 0 = 105°, Р = 35° при:

1 — Я,„ =0.1; 2 — Я,„ = 0.3; 5 — Я,„ = 0.5; 4 — Я,„ = 0.7; 5 — Я,„ = 0.9

р

О 0.2 0.4 0.6 0.8 1

Рис. 10. Профили скорости (сплошные кривые) и энтальпии (штриховые) в направлении набегающего потока для крыла с геометрическими параметрами 0= 105°,

(3 = 35° при:

1 — Я,„ = 0.1; 2 — Я,„ = 0.5; 5 — Я,„ = 0.9

1 -0.5 0 0.5 1

Рис. 12. Распределение толщины вытеснения на крыле при: 1 — у = 1.4;2 — у = 1.2; 5 — у = 1.1

-0.5 0 0.5 1

Рис. 13. Распределение давления по угловой координате 0 при:

1 — у= 1.4; 2 — у= 1.2; 3— у= 1.1

Рис. 14. Распределение тепловых потоков для крьша с Нм, = 0.5 и геометрическими параметрами © = 105°, (3 = 35° при:

1 — Рг = 0.72; 2 — Рг= 1

при изменении Ик от 0.1 до 0.9 возрастает почти в 2 раза, что вызывает повышение давления также примерно в 2 раза (рис. 9). На рис. 8 также представлены изолинии числа Маха М= 1, которые определяют размер области дозвукового течения в пограничном слое. При повышении температуры поверхности крыла увеличивается размер области дозвукового течения внутри пограничного слоя. Так для крыла с температурным фактором Н№ =0.1 (кривая 1) область дозвукового течения составляет примерно 40% от толщины пограничного слоя, а при Нк =0.9 (кривая 5) — более 60%. Кроме того, для крыльев с Нм> = 1 характерны более наполненные профили скорости и энтальпии (рис. 10).

Рассмотрено влияние показателя адиабаты у = 1.1, 1.2, 1.4 на течение в пограничном слое для крыла с параметрами 0 = 105°, Р = 35°, Нп. — 0.5 и числе Прандтля Рг = 0.72. При уменьшении у от 1.4 до 1.1 наблюдается уменьшение коэффициента трения примерно на 40% (рис. 11). Теплообмен, распределение которого качественно совпадает с трением, также уменьшается при снижении показателя адиабаты примерно на треть. При этом пик трения и тепловых потоков немного сдвигается при изменении показателя адиабаты. Толщина вытеснения пограничного слоя уменьшается примерно на 20% (рис. 12). При этом положение минимума толщины вытеснения зависит от показателя адиабаты, как и положение максимума тепловых потоков. Из-за меньшей толщины вытеснения давление уменьшается почти на 30% (рис. 13).

Исследование влияния числа Прандтля проведено при обтекании крыла с параметрами 0 = 105°, Р = 35°, Нп. =0.5 и значением у= 1.4. Изменение числа Прандтля от 0.72 до 1 повышает тепловые потоки на крыле примерно в 1.5 раза (рис. 14). При этом локальное повышение трения становится немного более выраженным. На распределение напряжения трения, давления и толщины вытеснения изменение числа Прандтля почти не оказывает влияния при Н№ =0.5. Для других температурных факторов число Прандтля может оказывать большее влияние на течение.

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

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

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

Изменение показателя адиабаты от 1.4 до 1.1 приводит к уменьшению трения, тепловых потоков, толщины вытеснения и давления примерно на 20—40%. Увеличение числа Прандтля от 0.72 до 1 приводит к увеличению тепловых потоков в 1.5 раза. Для крыла с температурным фактором Н№ = 0.5 остальные характеристики течения почти не зависят от числа Прандтля.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 10-01-00173-а) и ФЦП ННПКИР ГК № 02.740.11.0154.

ЛИТЕРАТУРА

1. Башкин В. А., Дудин Г. Н. Теория гиперзвуковых вязких течений. — М.:

МФТИ, 2006, 328 с.

2. Н ей л ан д В. Я., Боголепов В. В., Дудин Г. Н., Липатов И. И. Асимптотическая теория сверхзвуковых течений вязкого газа. — М.: Физматлит, 2003, 456 с.

3. Дудин Г. Н., Нейланд В. Я. Теплообмен в окрестности точки излома передней кромки пластины при гиперзвуковом полете // Изв. АН СССР. МЖГ. 1980. № 3, с. 40—45.

4. Дудин Г. Н., Лыжин Д. О. Об одном методе расчета режима сильного вязкого взаимодействия на треугольном крыле // Изв. АН СССР. МЖГ. 1983. № 4, с. 119—124.

5. Хейз У. Д., Пробстин Р. Ф. Теория гиперзвуковых течений. — М.: ИЛ, 1962,

607 с.

6. Лойцянский Л. Г. Механика жидкости и газа. — 7-е изд. — М.: Дрофа, 2003,

840 с.

7. Годунов С. К., Рябенький В. С. Разностные схемы (введение в теорию). — 2-е изд. — М.: Наука, 1977, 440 с.

Рукопись поступила 14/ІІІ2010 г.