УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIX 198 8 ~
№ 5
УДК 533.6.011.3/.5: 533.695
РАСЧЕТ СВЕРХЗВУКОВОГО ОБТЕКАНИЯ КОМБИНАЦИИ КРЫЛА И ФЮЗЕЛЯЖА С ВЫДЕЛЕНИЕМ СКАЧКА УПЛОТНЕНИЯ ОТ КРЫЛА
С. М. Босняков, Е. В. Карпов, С. В. Михайлов
Предложен подход и разработана программа на алгоритмическом языке ФОРТРАН-1У, позволяющая рассчитывать поля течения и распределения параметров течения на поверхности крыла и фюзеляжа. Подход заключается в разбиении расчетной области на подобласти, построении расчетной сетки в каждой подобласти с последующим интегрированием полной системы уравнений Эйлера методом С. К. Годунова [1]. Скачки уплотнения от крыла и фюзеляжа выделены явным образом.
Одним из направлений повышения точности расчета пространственных сверхзвуковых течений идеального газа является выделение газодинамических разрывов. В настоящее время традиционным стало выделение головной волны [2] и некоторых внутренних разрывов [3—5]. В данной работе рассмотрена важная задача — расчет сверхзвукового обтекания крыла в компоновке с фюзеляжем с учетом взаимодействия скачка уплотнения от крыла с поверхностью фюзеляжа. В известных авторам работах (см. например, [6—10]) указанный вопрос не рассматривался.
1. Решается задача Коши для замкнутой безразмерной .^-гиперболической системы уравнений Эйлера, дополненной уравнением Бернулли. Система уравнений в правой прямоугольной декартовой системе координат записывается в следующем виде:
(ри)* + (ргОу + (?w)z = 0,
(Р + + (?uv)y + (?uw)2 = О,
(puv)x + (р + pv2)y + (ркда), = 0,
(р uw)x + (pvw)y + (р + р«У2)г = О,
®2 + V* + и2 + р- = В,
' 1 1 х — ip
В = const.
Значения скорости отнесены к критической скорости V*, плотности — к ПЛОТНОСТИ невозмущенного потока р<х>, давления—удвоенному критическому скоростному напору v. На границах возмущенной области вы-
полняются условия: непротекания на стенках и Ренкина—Гюгонио на скачке уплотнения. Решение системы уравнений Эйлера осуществляется с использованием метода сквозного счета С. К. Годунова и имеет первый порядок аппроксимации по шагу расчетной сетки.
Летательный аппарат располагается в декартовой системе коорди-
лат, так, что его строительная горизонталь совпадает е осью ОХ, плоскость симметрии—х плоскостью Z=0, нос фюзеляжа лежит в плоскости Х = 0, начало крыла — в плоскости Х = Хкр (рис. 1). Геометрия летательного аппарата задается в виде координат точек опорных сечений. В промежуточных сечениях координаты определяются линейной интерполяцией.
Возмущенная область представляется в виде совокупности криволинейных четырехугольников, называемых подобластями (рис. 1,6). Подобласти связаны друг с другом через границы. Границей подобласти может являться скачок уплотнения, волна разрежения, поверхность тела. Наряду с газодинамическими разрывами для удобства построения
Рис. 1
сетки в качестве границ используются «эйлеровы поверхности», на которых выполняются условия гладкости всех функций, кроме сеточной (пунктирные линии на рис. 1,6). Число узлов сетки по разные стороны от границы одинаково. Разностная схема записывается в дивергентном виде и является консервативной.
Расчеты проводились в диапазоне чисел Мсо = 2...2,5 углов атаки а = 0...10° при нулевых углах скольжения. Вследствие симметрии течения расчет ограничивался полуплоскостью Z>0.
Начальные данные на носовом конусе рассчитывались с использованием соотношений работы [8]. С этой целью определялись местные углы наклона носового конуса 6Кон и по соотношению е = 6Кл/6Кон= = (1,72 — 2/<?м°с) пересчитывались эффективные углы наклона клиньев бкл- Начальное положение головного скачка уплотнения и параметры за ним рассчитывались с использованием соотношений для косых скачков уплотнения [11] ПО известным значениям Моо и бкл-
Расчетная область на участке Х<Хкр строилась так, что ее внешней границей являлся головной скачок уплотнения, внутренней — поверхность тела, боковыми — отрезки плоскости симметрии Z = 0, заключенные между телом и скачком, см. рис. 1 , а. При ^=^кр структура расчетной области менялась — область разбивалась на шесть подобластей. Скачок уплотнения от крыла выделялся явным образом (рис. 1,6). Его начальное положение и параметры за ним рассчитывались по соотношениям для косых скачков уплотнения [11]. В подобластях 3—6 параметры потока в ячейках новой сетки пересчитывались при помощи линейной интерполяции по результатам расчета течения около изолированного фюзеляжа в сечении Х = Хкр.
В работе сопоставлены два подхода к выделению скачка уплотнения от крыла. В подходе, предложенном авторами, учитывается взаимодействие скачка уплотнешш от крыла с поверхностью фюзеляжа. Для этого точка пересечения скачка с фюзеляжем в каждом сечении Х = const выделяется специально. С этой точкой связывается «эйлерова поверхность», разделяющая подобласти 3 и 4, представляющая собой прямую линию. Ее уравнение записывается в виде:
У = К0 + R ■ sin 0,
Z = Z0 + R • cos 0,
где F0, Z0 — координаты точки взаимодействия скачка уплотнения с фюзеляжем, 0 — угол наклона прямой к оси OZ, R — радиус-вектор,
Y0=f(X), Z0 = ?(X), 6 = <!>(*).
Угол наклона 0 прямой к оси OZ подбирается исходя из условия построения наилучшей сетки. Траектория движения точки пересечения скачка уплотнения с фюзеляжем Y0 = f(X), Z0=cp(X) определяется расчетом. Алгоритм расчета фронта скачка уплотнения построен с использованием принципа Гюйгенса, изложенного в работе [5].
Во втором подходе, который реализован, в частности, в работе [6], скачок уплотнения выделяется как' граница возмущенной области, огибающая тело (рис. 1,в). Эта граница в начальных сечениях крыла не совпадает со скачком уплотнения. Предполагается, что в процессе счета по мере движения по координате X решение устанавливается и на некотором расстоянии от сечения Х=Хкр(ХЭ>^кр) выделенная граница совпадает с истинным положением скачка. Решения для скачка от крыла, полученные в рамках этих двух подходов, естественным образом
сближаются при значениях Х>Хи где координата соответствует условию 20=ф(А'1) =0. При значениях Хкр<Х<Х1 второй подход дает ошибку в решении.
Разбиение в процессе счета расчетной области на подобласти предъявляет дополнительные требования к используемому численному методу — метод не должен давать забросов параметров в решении на границах подобластей. Этому условию удовлетворяет модификация [12, 13] метода С. К. Годунова [1].
Решение, получаемое с помощью метода С. К. Годунова, зависит от формы ячеек расчетной сетки, но во многих случаях эта зависимость мала и ей можно пренебречь. Авторами разработан вариант программы с использованием схемы С. К. Годунова [1]. Расчеты проводились на ЭВМ БЭСМ-6. Время расчета одного режима обтекания составляет приблизительно 2,5—3 часа.
2. На рис. 2 и 3 приведены линии равных значений статического давления р/рас и чисел М, построенные в возмущенном поле у поверхности крыла и фюзеляжа в сечении ^=30. Скачки уплотнения .и границы областей возмущения показаны жирными линиями. Рисунки построены путем монтажа: в правой полуплоскости в области ниже крыла скачок уплотнения выделен по первому подходу, в левой полуплоскости — по второму. В области выше крыла — наоборот: слева волна
Поле Р
Рис. 2
разрежения выделена по первому подходу, справа — по второму. Скачок уплотнения от крыла представляет собой линию, замыкающуюся на поверхность фюзеляжа и кромку крыла, что указывает на сверхзвуковой характер обтекания кромки. Интенсивность скачка уплотнения убывает по мере приближения к поверхности фюзеляжа. Также в значительной степени убывает давление на крыле при движении от кромки к корню. Особенностью течения является значительное уменьшение интенсивности скачка уплотнения при переходе его с боковой поверхности на нижнюю. Увеличению давления на нижней поверхности крыла способствует наличие фюзеляжа. Анализ поля течения в окрестности фюзеляжа (без крыла) показывает, что давление растет по мере удаления от плоскости симметрии (число М соответственно падает). Это приводит к нарушению классической схемы обтекания крыла со сверхзвуковой кромкой — исчезает область постоянства параметров потока в окрестности кромки [14].-
Граница возмущенной области, рассчитанная с использованием второго подхода [6], в рассматриваемом сечении Х = 30 не соответствует истинной форме скачка уплотнения (левая полуплоскость, область ниже крыла на рис. 2 и 3). Реальный скачок уплотнения на участке —9<Z<—4 размыт, на что указывает сгущение изолиний параметров р/роо и М. Ширина области размытия составляет около 40% от диамет-
ра фюзеляжа в рассматриваемом сечении. Форма размытого скачка уплотнения приблизительно соответствует форме выделенного скачка. Можно определить также приблизительные координаты точки взаимодействия размытого скачка с поверхностью фюзеляжа. Уточнение этих данных требует значительного измельчения расчетной сетки. Результаты, полученные на одинаковых сетках в случае специального выделения скачка уплотнения (первый подход) и в рамках второго подхода [6], сильно отличаются (ср. параметры течения в правых и левых полуплоскостях на рис. 2 и 3). Это отличие обусловлено размытием скачка уплотнения в рамках второго подхода. Например, по рис. 2 (левая полуплоскость) можно сделать ошибочный вывод о наличии значительного градиента давления вблизи боковой поверхности фюзеляжа, в то время как анализ правой части графика показывает, что этот эффект вызван размытием скачка уплотнения. Отсюда следует, что, при наличии в потоке сильных внутренних разрывов, их целесообразно выделять.
В области выше крыла на рис. 2 и 3 представлены результаты расчета внешнего фронта волны разрежения в рамках обоих подходов. В случае специального выделения (левая полуплоскость) волна разрежения, также, как и скачок уплотнения, замыкается на верхнюю поверхность фюзеляжа. Выделение волны по второму подходу дает некоторую границу, которая пересекает ось О У в точке У~9,5. В отличие от скачка уплотнения волна разрежения не имеет резко выраженного скачка параметров на внешней границе. Это затрудняет хотя бы приблизительное построение истинного фронта волны при использовании второго подхода.
На рис. 4 приведено распределение статического давления по поверхности крыла и фюзеляжа в сечении Х = 28,2. По оси абсцисс отло-
жена относительная длина 1 = 1/10 (где /0 — длина кривой, описывающей форму фюзеляжа с крылом в данном сечении), по оси ординат — статическое давление. На графике выделены участки, соответствующие крылу и фюзеляжу. Расчет проведен при двух значениях шага расчетной сетки в подобластях. Обозначение 9X17 на рис. 4 означает, что во всех подобластях число узлов по нормали к телу равно 9, а вдоль тела— 17 (17X17—соответственно 17 узлов по нормали и 17 — вдоль). Сопоставление показывает, что на крыле и верхней части фюзеляжа уже на крупной сетке точность расчета достаточно высока. На нижней части фюзеляжа расхождение достигает 10% и обусловлено большим размером подобласти 3. Увеличение числа ячеек в вертикальном направлении улучшает аппроксимацию решения в этой подобласти. Положение скачка уплотнения от крыла на рис. 4 соответствует значению / = 0,21. В данном сечении скачок попадает на границу максимальной кривизны контура сечения фюзеляжа, которая является участком перехода скачка с боковой поверхности на нижнею. Особенности геометрии фюзеляжа обуславливают перестройку течения на этом участке: давление резко уменьшается, скачок уплотнения вырождается в границу волны разрежения. Этот эффект аналогичен вырождению скачка уплотнения в границу волны разрежения сбоку и сверху пластинки конечной ширины при обтекании ее под положительным углом атаки. Расчетное значение разности давлений за и перед скачком (выродившимся в волну) отрицательно и равно Ар/р^~—0,1. Вслед за разрежением следует сжатие и рост давления в 1,5 раза (участок 0,21 </<0,25).
Расчет течения в рамках второго подхода хуже описывает особенности эпюры давления в области взаимодействия скачка уплотнения с фюзеляжем и в корне крыла (см. рис. 4). На верхней поверхности крыла и фюзеляжа эти подходы дают близкие результаты. На рис. 4 нанесены темными кружками также расчетные значения давления на фюзеляже в отсутствие крыла. Видно, что в рассматриваемом сечении интерференция крыла и фюзеляжа приводит к уменьшению коэффициента подъемной СИЛЫ Су фюзеляжа вследствие того, что крыло сопряжено с фюзеляжем в его верхней суживающейся части.
Следует отметить, что по условиям обтекания кромка крыла, около которого рассчитывалось течение — сверхзвуковая. В этом случае решение в области кромки, полученное в рамках уравнений Эйлера, достаточно хорошо соответствует реальному течению. В случае дозвуковой передней кромки расчет дает ошибку на подветренной поверхности крыла. Ошибка обусловлена наличием в реальном течении в области кромки вихря, который не учитывается в расчете.
На рис. 5 приведены значения местных чисел М перед скачком в области взаимодействия скачка уплотнения с фюзеляжем и значения перепада давления в скачке в зависимости от координаты X. Немонотонность кривых является, по-видимому, следствием относительно грубого задания контура тела. Значения чисел М перед скачком меняются незначительно на всем участке расположения крыла. Перепады давления, в отличие от местных чисел М, на участке 24<Х<27 снижаются более чем вдвое. Полученные зависимости использованы для анализа возможности отрыва пограничного слоя. На графиках нанесены значения критического перепада, соответствующие числам М перед скачком, для двумерного взаимодействия скачка уплотнения с турбулентным пограничным слоем [15] и для трехмерного взаимодействия [16]. Сделанная оценка показывает, что на боковой поверхности возможен отрыв пограничного слоя на участке Х<25. Это предположение подтверждено экспериментом по визуализации течения методом «лазерный
Рис. 5
Ти»Ш& ШШ
Рис. 6
нож». Эксперимент был проведен на модели, геометрии которой соответствовала расчетная конфигурация. На рис. 6 представлены фотография и схема течения в сечении Х = 28,2 при числе Моо = 2, угле атаки <х=10. Положения плоскостей «лазерного ножа» располагались дискретно в интервале 18,2<Х<37,2 с интервалом АХ = 2,Ь. Визуализация течения на участке 18,2<Х<28,2 позволила установить, что в области взаимодействия скачка уплотнения от крыла с пограничным слоем фюзеляжа образуется зона отрыва, которая сворачивается в вихрь, рас-
пространяющийся вниз по потоку. Таким образом, реальное течение в этой области имеет достаточно сложную структуру, описание которой необходимо проводить с привлечением специальных алгоритмов или методов решения уравнений Навье-Стокса.
Авторы признательны А. К. Иванюшкину за внимание к работе и полезные обсуждения.
ЛИТЕРАТУРА
1. Годунов С. К. Численное решение многомерных задач газовой динамики. — М.: Наука, 1976.
2. Thomas P. D., Vinokur U. R., Bastianum R„ Conti R. J. (Numerical solution for the three-dimensional Hypersonic flow-f|ield of a hlunt delta body. — AIAiA J., vol. 10, N 7, 1972.
3. Годунове. К., Жданов А. И., Семендяев К. А. Численные методы решения одномерных неустановившихся задач газовой динамики. Доклад на Всесоюзном съезде по теоретической и прикладной механике.— М.: 1960.
4. М о г е 11 i J. Three-dimensional, supersonic, steady flows with any number of imbedded shocks. — AlAlA Paper, N 74-10, 1974.
5. M а к a p о в В. E. К выделению поверхностей разрывов при численном расчете сверхзвуковых конических течений. — Ж. вычисл. матем. и матем. физ., 1982, № 5.
6. 3 а р у б и н А. Г. Реализация конечно-разностного метода Крайко для расчета обтекания комбинации крыла и фюзеляжа сверхзвуковым потоком идеального газа. — Труды ЦАГИ, 1978, вып. 1941.
7. Marconi F., Jaeger L. Computation of high-speed inviscid flows about real configurations:—NASA SP-347, Langly Res Center, part. II, 1975.
8. А у к и н М. К., Т а г и р о в Р. К. Расчет свехзвукового обтекания модели самолета, имеющего сдвоенные двигатели. — Ученые записки ЦАГИ, 1982, т. 13, № 3.
9. Иванов О. В. Численное исследование сверхзвукового обтекания тел сложной формы идеальным газом.— В сб.: Моделирование процессов гидрогазодинамики и энергетики. ИТПМ СОАН, 1985.
10. Лобановский Ю. Н. Расчет обтекания сверхзвуковым потоком невязкого газа крылатых конических тел. — Ученые записки ЦАГИ,
1980, т. 11, № 6.
11. Лойцянский А. Г. Механика жидкости и газа. — М.: Наука,.
1970.
12. Ко л га н В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики. — Ученые записки ЦАГИ. 1972, т. 3, № 6.
13. Тилляева Н. И. Обобщение модифицированной схемы С. К. Годунова на произвольные нерегулярные сетки. — Ученые записки ЦАГИ, 1986, т. 17, № 2.
14. Зарубин А. Г. К расчету обтекания крыльев со сверхзвуковыми острыми кромками. — Ученые записки ЦАГИ, 1980, т. 11, № 6.
15. Абрамович Г. Н. Прикладная газовая" динамики. — М.:
Наука, 1976.
16. Korkegi А. Н. Asimple correlation for incipient turbulent Boundary separation dueto skewed shock wave. — AIAA J., 1973, vol. 11, N 11.
Рукопись поступила 27/IV 1987 г.