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

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

Том XXIII 1992 № 3

УДК 533.6.011.55 : 629.7.025.1 629.782.015.3 : 533.6.011.6

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

А. В. Гончар

Путем численного моделирования в рамках модели параболизованных уравнений Навье — Стокса исследуются сверхзвуковые течения вязкого теплопроводного газа около затупленных треугольных крыльев и комбинации «корпус — крыло». Исследован эффект возникновения на наветренной плоскости затупленного треугольного крыла пиков в тепловых потоках, известный из эксперимента. Эффект носит чисто трехмерный гиперзвуковой характер. Показано, что он обусловлен взаимодействием вязкого пограничного слоя с энтропийным. Расчетное распределение тепловых потоков в зоне эффекта хорошо согласуется с экспериментальным.

Считается, что решение уравнений Навье — Стокса (УНС) может дать полное описан-ие течений газа около гиперзвуковых летательных аппаратов при условии, что все важные масштабы явлений хорошо разрешаются расчетной сеткой. Поскольку прямое решение УНС требует больших вычислительных мощностей в 1 — 10 флопс, то на практике широко используют различные упрощения.

Параболизованные уравнения Навье — Стокса (ПУНС) являются эффективной альтернативой к полным уравнениям Навье — Стокса для сверх- и гиперзвуковых течений с большими числами Рейнольдса. Основное предположение, учитываемое ПУНС, состоит в незначительности влияния возмущений, идущих против течения (слабая эллиптичность). Это предположение разумно для большей части поля течения из-за доминирующей роли гиперболического невязкого потока и параболического пограничного слоя. Как и УНС, модель ПУНС автоматически учитывает эффекты вязко-невязкого взаимодействия, что позволяет моделировать течения с отрывами потока около тел сложной формы.

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

Исследован эффект возникновения на наветренной плоскости треугольного крыла пиков в тепловых потоках, обнаруженный ранее в экспериментах.

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

1. Описание метода. В данной работе ПУНС для ламинарных течений вязкого теплопроводного совершенного газа следуют из уравнений Навье — Стокса, записанных относительно произвольной криволинейной системы координат, связанной с телом, в результате следующих упрощений:

нестационарные члены равны нулю;

производные в касательных к телу направлениях в членах с вязкостью незначительны по сравнению с производными в нормальном к поверхности тела направлении;

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

Согласно проведенному анализу получающиеся уравнения имеют относительно продольного направления неполно параболический тип. Задача Коши вдоль этого пространственного направления корректна.

На поверхности тела задаются граничные условия прилипания (1?=0) и фиксируется температура поверхности Тш или тепловой поток Внешняя граница расчетной области совпадает с головным скачком уплотнения, на ней ставятся условия Ренкина — Гюгонио. На боковых плоских срезах ставятся условия зеркальной симметрии. Коническая поверхность начальных данных располагается в сверхзвуковой части сферического носового затупления. Начальные данные для ПУНС определяются интерполяцией осесимметричного решения для задачи обтекания сферы. Решение этой задачи в рамках уравнений Навье — Стокса методом установления [1] занимает около двух минут на ЭВМ ЕС-1066 и предшествует расчету с помощью ПУНС.

Расчетная сетка на каждой маршевой поверхности строится численно с помощью встроенного генератора, использующего решение уравнений в частных производных эллиптического типа для сетки. Осуществляется автоматическая адаптация сетки в нормальном к поверхности тела направлении в соответствии с градиентами полей скорости и плотности в газодинамическом решении. Около 40% точек от общего их числа М/ в нормальном направлении формируют равномерное распределение, на фоне которого происходит сгущение 40% точек в зонах повышенных градиентов скорости, а 20% — в зонах повышенных градиентов плотности. Такой подход обеспечивает адекватное разрешение расчетной сеткой пограничных и вихревых слоев, не имея о них никакой априорной информации. На каждом маршевом шаге газодинамическое решение и расчетная сетка полностью согласованы.

Для численной аппроксимации ПУНС использована неявная схема. Производные в членах с вязкостью аппроксимированы со вторым порядком центральными разностями. В маршевом направлении используются односторонние разности второго порядка. Для невязких членов применены ТУБ аппроксимации третьего порядка по технологии схем типа Годунова с приближенным решателем задачи Римана, предложенным в работе [2]. Проведено их обобщение на неравномерные сетки. Параллельно на конкурентной основе используется вариант программы с центральными конечными разностями и нелинейными членами искусственной диссипации, обладающий вторым порядком аппроксимации.

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

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

Расчёт начинается с выбора величины маршевого шага, который обычно соответствует числу Куранта, равному 1 — 2 для зон сверхзвукового течения. Задается распределение точек по контуру тела в новом маршевом сечении. По наклону скачка определяется его положение на новом шаге. Начальным приближением к решению в новом маршевом сечении служит решение с предыдущего шага. По этим данным строится расчетная сетка. Решаются линеаризованные уравнения газодинамики и уточняются форма скачка и текущее приближение к газодинамическому решению. С учетом новых данных производится перерасчет сетки. Переход к новому маршевому шагу оказывается устойчивым уже после двух итераций. Критерием окончания итерационного процесса является сходимость во всей расчетной области плотности, давления и скорости с относительной погрешностью 10~3. Контролируется также итерационная сходимость в распределении тепловых потоков qw и коэффициентов трения С]. Реально число итераций на маршевом шаге лежит в пределах 3—6, и, как правило, уменьшается с уменьшением маршевого шага.

Чтобы иметь информацию о контуре тела в каждом маршевом сечении, было решено вписывать произвольное тело набором плоских сечений, ортогональных некоторой трехмерной кривой F(t) — криволинейной оси тела. В каждом сечении t — const контур тела представляет собой составную плоскую кривую, каждый сегмент которой является одним из простых объектов: линейным отрезком, дугой эллипса, дугой кубической кривой и т. д. При задании геометрии тела, создается подпрограмма, которая должна для любого значения параметра t выдавать параметры составной кривой: количество сегментов, параметры каждого сегмента. Для тел простой формы эти параметры можно задать явно. В более сложных случаях используется библиотека подпрограмм численного определения точек пересечения, касания между основными типами кривых. Помимо определяющих сегмент параметров, вычисляются или задаются и производные этих параметров относительно /, что позволяет аналитически определить не только координаты точек на составной кривой, но и нормали к поверхности тела.

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

2. Об оценке точности и достоверности расчетов. Ниже будут представлены результаты расчетов течений вязкого газа около затупленных треугольных крыльев и комбинации «корпус — крыло». Но прежде несколько слов об оценке точности и достоверности получаемых результатов.

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

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

и

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

Данная программа в ее невязком варианте тестировалась на равномерных сетках в задаче обтекания затупленных конусов. Был доказан второй порядок аппроксимации. Результаты расчетов на 30 калибров на сетках с 17 и 33 точками по нормали к телу отличались по всем параметрам менее чем на 0,5%.

Расчеты течений вязкого газа методом ПУНС проводились на неравномерных сетках. Разворот течения с параметрами М^, = 6, Ие» = 3,55-104, Тт/То = 0,7 на сфере от 0 — 54° до 0 = 90° исследовался на последовательности сеток с 26, 30, 40, 60, 90 узловыми точками по нормали к телу и с разными законами сгущения. Было установлено, ^то 30 — 40 точек достаточно для расчета с погрешностью в 1%.

Использование адаптирующихся к газодинамическому решению сеток при М/ = 40 гарантирует, что не менее 16 точек создадут равномерное фоновое распределение узлов, способное, как было показано, неплохо разрешить крупномасштабные особенности течения. Около 16 точек всегда попадут в пограничный слой, будь он присоединенный или оторвавшийся. Еще около 8 точек будут отслеживать градиенты плотности по мере их развития. Ясно, что такого количества точек вполне достаточно для схемы второго порядка, чтобы обеспечить решение с погрешностью на уровне 1%.

С разрешением особенностей течения в азимутальном направлении дело обстоит сложнее, поскольку адаптация сетки в этом направлении не производится. Для затупленных треугольных крыльев стреловидностью 65 — 75° равномерное распределение МК = 60-т- 120 точек по контуру тела приемлемо при расчетах на длину 10 — 20 калибров. Такие расчеты выявили основные особенности течения, что позволило в дальнейшем проводить целенаправленное сгущение узловых точек. Как правило, не менее 20 точек всегда присутствует на кромке крыла. Обеспечивается плавное изменение шага сетки при переходе с кромки на крыло. К зоне на наветренной плоскости крыла, где были выявлены пики теплового потока, было дополнительно привязано около 20 узловых точек, которые не уходили из нее с ростом размаха крыла. Переход в процессе марша ПУНС с равномерного распределения на неравномерное происходил плавно, по мере того как равномерный шаг сетки оказывался больше, чем было предписано по закону сгущения.

Все расчеты, в пределах 20 — 30 калибров дублировались на различных сетках: менялись как число узловых точек, так и функции* их сгущения. Для схемы с центральными разностями выяснялось влияние на решение уровня искусственной диссипации. Исследовалось влияние критерия сходимости итераций, параметров регуляризации. Было установлено, что общая картина течения инвариантна относительно всех этих изменений. Для задачи обтекания конфигурации «корпус — крыло» потоком газа с параметрами Мос=14, Ие^ = 10\ Та/То = 0,3 и размерами сетки МК = 120, М1 = 40 на длину 50 калибров был сделан вывод о погрешности полученного решения на наветренной стороне крыла в 1%. Расчеты на наиболее грубых сетках с числом узловых точек МК = 60, Аи = 30 имеют погрешность около 5%. На подветренной стороне ошибки, особенно в уровнях тепловых потоков, могут быть и больше.

Размер маршевого шага для представленных задач варьировался в пределах от 0,05 калибра вблизи носового затупления до 0,5 калибра вдали от него. Число маршевых шагов составляло несколько сотен. Число расчетных точек полной трехмерной сетки 120 X 40 X 300 достигает миллиона. Характерное время расчета одного варианта составляет 15—30 часов на ЕС-1066.

3. Результаты расчетов.

Расчет 1. Исследовалось течение газа с параметрами = 6,8, Ие«, к = = 1,3-105, Тж = 60,2 К, Тш = 317 К, у = 1,4, Рг = 0,72 около затупленного

Рис. 2

в носовой части по сфере, а на кромке по цилиндру треугольного крыла стреловидностью 70° и длиной 10 калибров, обтекаемого под углом атаки а= 10°. Расчет был проведен на сетке, содержащей МК = 60 и МУ = 40 узловых точек. На рис. 1 ,а,б,в представлены соответственно поля скорости, плотности и расчетная сетка в сечении X//? = 9,45. На рис. 2 приведено распределение давления по контуру тела в сечении Х/Я = 9. По оси абсцисс отложена длина дуги контура тела, отсчитываемая от наветренной образующей крыла. По оси ординат отложено давление на поверхности тела, отнесенное к Роо и\о- Сплошная линия представляет результат расчета ПУНС, светлые квадраты — данные эксперимента [3]. Результаты расчета и эксперимента хорошо согласуются между собой.

Расчет 2. Исследовалось течение газа с параметрами М,„ = 14, Ре^ д = = 1 • 10\ = 32 К, Тт = 300 К, у = 1,4, Рг = 0,72 около комбинации «кор-

пус— крыло», обтекаемой под углом атаки а =10° к плоскости крыла. Корпус представлял собой затупленный по сфере конус с полууглом раствора 5°. Треугольное крыло стреловидностью 75° состыковано с корпусом по схеме низкоплан. Место стыковки сглажено с радиусом, равным половине радиуса конуса для данного сечения. Радиус носового затупления в два раза больше радиуса затупления кромки крыла. Расчетная сетка содержала МК = 120 и М] = 40 точек. Маршевые поверхности, за исключением носовой части, были плоскими и ортогональными оси конуса.

На рис. 3 представлены поля плотности в ударном слое на наветренной части крыла для значений продольной координаты X« 16, 21, 54. Видно появление характерного прогиба в энтропийном слое в окрестности четырех калибров от плоскости симметрии тела. Этот прогиб формируется в месте взаимодействия головного скачка уплотнения с кромкой крыла. Здесь находится локальный максимум в давлении и тепловых потоках и минимум в отходе головного скачка от кромки. Особенность в энтропийном слое, возникшая на кромке, сносится в дальнейшем течением на плоскости крыла. Ее удаление от плоскости симметрии остается неизменным с ростом X и равным примерно четырем калибрам.

На рис. 4 представлено распределение теплового потока по поверхности тела в сечении Х//?«51, По оси абсцисс отложено значение координаты У вдоль размаха крыла. Значение У = 0 соответствует плоскости симметрии тела. По оси ординат отложено значение теплового потока к поверхности тела, отнесенное к РжИ'п- Сплошной линией представлен результат расчета ПУНС,

•м

z

-з

+ ♦

Рис. 6

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

На рис. 5 представлена величина А отхода внешней границы энтропийного слоя от поверхности крыла в зависимости от координаты У при Х/к= 54,43. Внешняя граница определена как место, где скорость потока равна 0,9 от скорости набегающего потока иВ области четырех калибров от плоскости симметрии наблюдается локальный минимум А, соответствующий максимуму В тепловом потоке. Положим для оценки <7а,~(1/Й) и выберем коэффициент пропорциональности таким, чтобы оценочный тепловой поток совпал с определенным с помощью ПУНС в одной точке У = 10. Полученное таким образом оценочное распределение qw представлено штриховой линией на рис. 4. Как видно из рисунка, эта оценка столь же ^хорошо согласуется с экспериментальными данными на всей поверхности крыла, как и та, что получена с помощью ПУНС. Таким образом, объяснение пика в тепловых потоках заключается в форме вихревого слоя.

Толщина вихревого слоя Л на наветренной поверхности треугольного крыла вблизи плоскости симметрии соответствует толщине этого слоя в момент

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

На рис. 6 представлены результаты расчета теплового потока </«, на поверхности тела, отличающегося радиусом носового затупления. Теперь он уменьшен в два раза и равен радиусу затупления кромки крыла. Расчетная сетка содержала МК = 60 и MJ = 35 точек. Система координат, связанная с крылом, и единица измерения длины, равная двум радиусам затупления кромки крыла, остались прежними. Возможно прямое сравнение с результатами предыдущего расчета. Первая особенность нового расчета заключается в появлении пиков фо,, начиная с меньших X. Вторая и главная особенность — это уменьшение именно в два раза расстояния от плоскости симметрии до пиков фш, что вполне согласуется с предложенным выше объяснением эффекта появления пиков.

Расчет 3. В третьем расчете исследовалось течение газа с параметрами М00 = 6, Ие«, я = 3,55-10 , (Та/Тж) = 5,74, 7=1,4, Рг = 0,70, (ц/ц.) = = (7’/7'<х>)0' около такой же комбинации «корпус — крыло», что и в расчете 2, обтекаемой под углом атаки а= 10° к плоскости крыла. Расчетная сетка

Г-14,10М-6,0 ;Ж1=3,55-10Ч1 Т„=5,7Ч

Рис. 7

Рис. 8

содержала МК = 120 и М/ = 40 точек. Расчет был выполнен на длину 86 калибров.

На рис. 7 приведена расчетная сетка в сечении X//? = 14. По сгущению линий сетки можно судить о положении пограничного слоя, оторвавшегося от кромки крыла, и месте его присоединения к боковой поверхности конуса.

На рис. 8 представлено поле плотности в сечении X//? = 86. На наветренной плоскости крыла вблизи плоскости симметрия находится неплохо разрешенный сеткой вихревой слой, имеющий цилиндрическую форму. На подветренной стороне виден скачок уплотнения, сформированный при взаимодействии поперечного течения с корпусом тела, а также сошедший с кромки крыла пограничный слой. Как и в расчете 2, наличие вихревого слоя на наветренной плоскости крыла ведет к появлению локальных пиков в тепловых потоках, но меньшей амплитуды. Более резким является переход к минимальному значению на плоскости симметрий.

ЛИТЕРАТУРА

t

1. Головачев Ю. П., Леонтьев Н. В. Циркуляционное течение у лобовой поверхности сферы, обтекаемой сверхзвуковым потоком типа следа// Изв. АН СССР МЖГ,—1985, № 3.

• 2. R о е P. L. Approximate Riemann Solvers, parameter vectors, and difference schemes//J. Comp. Phys.— 1981. Vol. 43.

3. Tannehill J. C., Venkatapathy E:, Rakich J. V. Numerical solution of supersonic viscous flow over blunt delta wings//A!AA Paper.—1981, N 49.

Рукопись поступила 17/VII 1990 г.