автореферат диссертации по авиационной и ракетно-космической технике, 05.07.03, диссертация на тему:Исследование перспективных схем интегрированияуравнений Эйлера для трансзвуковых чисел Маха взадачах аэроупругости

РГБ ОД

6 1 ЛНВ 2000

Московский Физико-Технический Институт (государственный университет)

На правах рукописи

УДК 629.735.33

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

уравнений Эйлера для трансзвуковых чисел Маха в

____

задачах аэроутгругости

Азаров Александр Юрьевич

(05.07.03 - Прочность летательных аппаратов)

Автореферат диссертации на соискание ученой степени кандидата технических наук

Жуковский - 1999

Работа выполнена в Центральном Аэрогидродинамическом Институте имени профессора Н.Е. Жуковского.

Научный руководитель — кандидат технических наук П.Г. Карклэ Официальные оппоненты — Микеладзе В.Г., доктор технических наук,

профессор, ЦАГИ.

Баранов Н.И., кандидат технических наук, МАПО "МИГ" Ведущая организация —■ АНПК "ОКБ Сухого"

Защита состоится "____" ............. 1999 г. на заседании Диссертационного совета К 063.91.07 при МФТИ.

Адрес: 140160, Московская обл., г. Жуковский, ул. Гагарина, д. 16 С диссертацией можно ознакомиться в библиотеке ФАЛТ.

Автореферат разослан "____" ............. 1999 г.

Ученый секретарь Диссертационного совета к.ф.м.н. доцент

А.И. Киркинский

0Я4 -046.63-0} <И4б;0

Введение

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

Актуальность

Несмотря на интенсивное развитие вычислительных методов аэродинамики и постоянное совершенствование вычислительной техники, позволяющее < решать все более сложные задачи, до сих пор значительные усилия направлены на создание надежных и быстрых расчетных методов аэроупругости для трансзвукового диапазона скоростей. Авиационная промышленность остро нуждается в таком инструменте, тем более что расчетные режимы полета современных пассажирских самолетов большой вместимости все ближе подходят к скорости звука. Препятствий на пути создания универсального и удобного инструмента много. Физическая картина трансзвукового потока сама по себе достаточно сложна —из-за возникновения в потоке сверхзвуко- у вых областей, и, как следствие, газодинамических разрывов (волн разрежения и слабых ударных волн) на их границах. Эта сложность многократно возрастает в нестационарной задаче, где каждое из различных явлений, таких как

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

Большинство расчетных методов определения нестационарных аэродинаг-мических нагрузок на трансзвуке обладают рядом недостатков. Как правило, затраты времени на создание расчетной сетки для сложной области достаточно велики. Точность расчетных методов существенно различается даже в случае решения одной системы уравнений и может зависеть от навыков и опыта исследователя в подборе различных параметров. Расчетная трансзвуковая аэродинамика чувствительна к ошибкам в алгоритмах и их реализации, что может приводить к неизбежным ошибкам в решении. Значительную трудность представляют потери в точности, обусловленные такими особенностями расчетной нестационарной аэродинамики, как отражения волн, распространяющихся от тела, на внешних границах и даже на границах различных по размеру ячеек.

Представлено в диссертации

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

Соискателем разработаны алгоритмы и создана программа для расчетов

стационарного и нестационарного обтекания тела или конфигурации тел сжимаемым невязким потоком газа на расчетной сетке, состоящей из криволинейных шестигранных блоков по расчетному методу второго порядка аппроксимации Родионова, либо методу Годунова первого порядка аппроксимации. Лично соискателем получены результаты расчетов по методам Родионова и Годунова стационарного обтекания крыла RAE А с профилем RAE 101 для числа Маха 0.8 под углом атаки 2° и числа Маха 0.9 под углом атаки 0°, ста^ ционарного и нестационарного обтекания горизонтального оперения AGARD SMP для расчетного случая M = 0.95, а — 0.4° (для нестационарного случая движение крыла описывалось периодическим изменением угла атаки около положения 0.4° с амплитудой 0.4° и приведенной частотой 0.583). Выполнено сравнительное исследование точности и сходимости новых алгоритмов для стандартных тестовых конфигураций.

Достоверность результатов

Достоверность результатов подтверждается сравнением с результатами, полученными с помощью расчетной программы С.И. Кузьминой [4] и эспери-ментальными и теоретическими данными работы [5].

Диссертация является продолжением работы [б] по созданию программы для расчета аэродинамических нагрузок в трансзвуковом диапазоне чисел Маха для применения в задачах аэроупругости. Результаты исследований докладывались на конференциях молодых специалистов МФТИ 1996 г. .[7] и 1997 г. [8]. Результаты работы, проведенной в целях сравнения характеристик методов первого (схема Годунова) и второго (схема Родионова) порядков аппроксимации были доложены на конференции "Современные проблемы аэрокосмической науки" в 1998 г. [9].

Выносится на защиту

IIa защиту выносится разработанный метод расчета аэродинамических ^ нагрузок в трансзвуковом диапазоне скоростей для применения в задачах аэроупругости.

1. Теория. Расчетные схемы

Для интегрирования уравнений Эйлера применяется конечноразностная схема второго порядка аппроксимации A.B. Родионова [2]. Схема относится к схемам типа "предиктор-корректор" и, наряду с рядом других схем второго порядка аппроксимации, базируется на методе Годунова. Повышение порядка аппроксимации достигается за счет замены кусочно-постоянного распределения газодинамических переменных (ГДП) в расчетной области на кусочно-линейное. При определении приращений используется принцип минимальных производных, предложенный В.П. Колганом [10]. Для оценки влияния порядка аппроксимации была также реализована схема первого порядка аппроксимации Годунова [1].

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

2. Реализация расчетного метода

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

Каждый блок сетки имеет собственный, независимый, диапазон индексов i, j, к, представляющих собой логические координаты. Соседние блоки явно стыкуются по границам с точным соответствием количества и расположения узлов. Таким образом устраняется процедура интерполяции из одного блока сетки на границу другого, применяемая в методах с перекрывающимися

сетками. Сам блок характеризуется следующими данными, которые подробно рассмотрены далее:

1) набором интервалов индексов, характеризующих размерность блока и число "заграничных" ячеек для обеспечения граничных условий;

2) поверхностями, задающими грани блока;

3) массивами, определяющими координаты узлов сетки, центров ячеек и их граней, а также объемы ячеек и площади их граней;

4) информацией о связях данного блока с соседними, необходимой для "обмена" геометрическими данными сетки.

Основной объект манипуляций при работе с индексами —■ интервал индексов в трехмерном пространстве, фактически представляющий собой параллелепипед. Для удобства реализации используется два типа интервалов. Один индексирует центры ячеек сетки (на рис. 1 это интервалы АА и ВВ), а другой узлы сетки (на том же рисунке это СС и йИ).

Блок сетки характеризуется в первую очередь областью определения ячеек сетки. Это "внутренняя" область блока, рассчитываемая целиком на основании знания координат узлов границ блока и именно в области определения заданы величины, рассчитываемые в процессе счета. На рис. 1 область определения показана, как прямоугольник с крайними индексами АА. Область данных отмечена как ВВ и включает в себя как область определения, так и "заграничные" ячейки, которые используются для обмена сеточными величинами и переменными счета с соседними блоками сетки и для вычисления граничных условий. Как область определения, так и область данных задаг ны на центрах ячеек сетки. На узлах сетки задана область данных узлов СС, которая просто охватывает область данных ячеек ВВ и область определения узлов БО, связь которой с областью определения ячеек сложнее. Здесь и далее область определения либо область данных обозначает соответствующую область на ячейках сетки; об областях, определенных на узлах сетки, будет сказано особо.

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

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

• координаты узлов сетки;

• нормали к граням ячеек сетки;

• площади граней ячеек;

• центры граней.

И функции, определенные на ячейках блока сетки:

• объемы ячеек блока;

• центры ячеек блока.

Расчет сетки в блоке

При расчете сетки сначала заполняются узлы на границах области определения узлов (см. рис. 1) — т.е. грани блока. Далее рассчитываются координаты узлов сетки внутри области определения ББ по методу трансфинитной интерполяции.

Следующим шагом координаты узлов экстраполируются за область определения йС -— эта процедура призвана заполнить "заграничные" ячейки сетки узлами в случай, если данная грань не отмечена как "свободная" и, вместе

с тем, не имеет соседнего блока. Для фиктивного узла сетки (г]'к), соседний узел которого (¿4- 1,ук) является узлом сетки, простая формула экстраполяции по двум соседним точкам выглядит так:

%*) — 2®(»+1о*) _ к)

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

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

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

Расчет координат узлов грани блока

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

быстродействия узлы грани записываются сразу в сеточную функцию узлов блока в соответствующие координаты {,], к.

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

я

Л (»,.?) =

(1)

7(м) = Л(м) + Х>С?)[ы*') - мш]

*=1

Расчет координат узлов ребра блока

Для того, чтобы применить формулу (1) и вычислить координаты узлов грани блока сетки, сначала требуется вычислить координаты узлов ребер блока и заполнить элементы с соответствующими индексами в сеточной функции узлов блока. Очевидно, что ребро блока всегда образуется изменением только лишь одной из координат {,], к, оставляя остальные две неизменными. Этим естественным свойством, вытекающим из того, что форма блока в индексном пространстве образует параллелепипед, удобно пользоваться при реализации алгоритмов вычисления координат.

Выл разработан ряд вариантов заполнения ребра блока сетки координатами узлов.

1) Линейное распределение точек по ребру. Задаются крайние точки о и Ь и координаты узлов ребра распределяются равномерно между ними. Так, для ребра I = [¿о, »1],$ = .70, к = к0:

£ = а 4- *—~-(5—&) »1-го

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

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

Расстояние между узлами будет меньше у нижнего конца ребра г'о, если 5 < 0. Иначе, при д > 0, узлы будут сосредоточиваться к г\.

3) Иногда бывает необходимо сосредоточить узлы сетки у точки, не совпав дающей с концами ребра и, к тому же, узлы могут быть расположены на кривой, имеющей сложное описание в пространстве (например, профиль крыла), заданной не аналитически, а массивом точек. Предлагается алгоритм для реализации такого типа ребра, который действует в два этапа. Сначала на точки, заданные заранее, "натягивается" сплайн, т.е. рассчитываются параметры аналитической кривой, проходящей через заданные точки. Далее сосредоточение узлов проводится уже на сплайне.

Так как классический сплайн является скалярной функцией одного аргумента, возникает вопрос об интерполяции сплайном пространственной кривой, имеющей три компоненты х,у,г. Для этого применяется широко распространенный метод, суть которого в том, что вводится переменная "пути" кривой 0 < ^ 1,» = г0,...,н, используемая далее как независимая переменная. Переменная пути записывается, как

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

¿1

«5 • ■

5; = —, 1 = 10,...,»!

Р)

заданный набор узлов, в зависимости от скалярной переменной в, определенной в диапазоне 0 ^ в ^ 1. Такой диапазон изменения независимой переменной предпочтителен для более простого представления функций сгущения.

Использование функций сгущения для изменения интервалов координат является наиболее часто используемым средством контроля за генерацией сетки. Для реализации была выбрана функция сгущения с гиперболическим тангенсом (диапазоны аргумента и функции были приведены к интервалу [0,1]):

гм - {п+°(тч-с))-ъч-Щ}

ПП> ~ {1 +а(Ш[Ь(1 - с)] - 1Ь[-Ьс])} ^

Для концентрации узлов около точки с параметр а должен быть отрицателен. Чтобы отображение было взаимно однозначным, требуется выполнение условия — аЪ < 1. Параметр с может быть выбран лежащим вне интервала [0,1] для концентрации узлов у края интервала. На рис. 2 изображены графики }{г)) при различных значениях параметров а, Ь и с.

Учет движения сетки

При движении сетки для алгоритма расчета требуется знание ряда сеточных функций, характеризующих переход от геометрии блока сетки на момент времени £ к моменту 4 + т. Индексное пространство блока остается неизменным. Исходными данными для вычисления параметров движения сетки служат два положения блока на оба момента времени. Движение узлов сетки предполагается равномерным от начального положения в момент времени Ь в течение временного интервала г. Все сеточные функции, описывающие движение сетки, определены на узлах блока сетки и описывают поведение граней ячеек. Они включают объем, заметаемый гранью ячейки при перемещении с одного положения в следующее, площадь и нормаль к грани и скорость ее движения.

Расчет течения в блоке

Расчет поля течения в блоке сетки может осуществляться как по методу Годунова, так и по методу Родионова. Если заведомо известно, что сетка, не будет перемещаться, в алгоритмах используются формулы, упрощенные с учетом этого предположения. Так, например, очевидно, что для статичной сетки объем, заметаемый гранью ячейки при движении сетки, равен нулю.

Для хранения газодинамических параметров (ГДП), подлежащих расчету, создается сеточная функция, определенная на ячейках блока сетки, в каждой ячейке хранится вектор £/. Фиктивные ячейки используются для задания граничных условий различных типов. Значения ГДП {/ в фиктивных ячейках вычисляются на основании величин и в соседних ячейках, принадлежащих области определения блока.

В дополнение к сеточной функции и в методе Родионова требуются массивы для хранения приращений 511 к граням ячейки и оценочных величины {/'+т. На эти сеточные функции, как и на II, требуется постановка соответствующих граничных условий.

Объединение блоков сетки в расчетную область

Для объединения блоков сетки в единую расчетную область требуется состыковать границы смежных блоков. На практике это означает, что размерности границ блоков должны совпадать и все сеточные функции, участвующие в процессе расчета, должны "обмениваться" данными со своими соседями из смежных блоков. Это касается как сеточных функций методов Годунова и Родионова и, 511 и 0, так и, что не менее важно, собственных сеточных функций блоков сетки, поскольку геометрические данные сетки для корректности расчета должны быть верными в фиктивных ячейках.

Алгоритмы расчетного метода (Годунова либо Родионова) должны использовать один шаг по времени, наименьший среди ограничений рассчитанных в блоках. Абсолютная ошибка Я на расчетной области равна максимуму вычисленных величин Л в блоках.

3. Результаты численных исследований

В целях сравнения характеристик методов Годунова и Родионова автором была проведена модификация расчетной программы С.И. Кузьминой таким образом, чтобы можно было рассчитывать течение существующей программой, но по схеме Родионова. Целью данной работы [9] было получение предварительных результатов, позволяющих сравнить в интересующем нас трансзвуковом диапазоне скоростей такие характеристики методов Годунова и Родионова, как точность расчета особенностей течения, сходимости и времени счета и сделать вывод о целесообразности применения схемы Родионова.

Программа С.И. Кузьминой [4] предназначена для расчета трехмерного нестационарного течения невязкого идеального газа около крыла по схеме Годунова решения уравнений Эйлера. В используемой регулярной сетке передняя, задняя и боковые грани ячеек всегда параллельны плоскостям ху, yz. Только нижняя и верхняя грани могут быть не параллельны плоскости xz над и под крылом. В этом случае считается, что все узлы, расположенные в этой области, смещаются вверх или вниз на величину у поверхности крыла. Столь простая организация сетки позволяет достаточно компактно представить вычисление потоков в разностной схеме, тем самым значительно ускоряя процедуру интегрирования. Одновременно значительно снижаются требования к памяти вычислительной системы на хранение геометрических характеристик сетки — требуется знать лишь шаги сетки в направлениях осей и углы наклона граней над и под крылом.

Малые колебания крыла или его части (например, элерона) вводятся как добавки к углам наклона граней ячеек в соответствующей области над и под крылом. Можно заметить, что разбивка сетки, не учитывающая геометрию крыла в плане, приводит к возникновению дробных ячеек у передней и задней кромок. Как следствие, это требует специальной обработки в алгоритме численного интегрирования.

По исходной и модифицированной программам были выполнены расчеты стационарного обтекания крыла RAE А с профилем RAE 101 и проведено сравнение с экспериментальными и расчетными результатами, полученными

в работе [5]. Исследовалось распределение коэффициента давления по профилю на 15% размаха крыла, графики сходимости, т.е. зависимости ошибки от времени процесса t, а также времена счета и количество шагов до достижения заданной точности. (

Как видно из сравнения с расчетами по теории потенциала малых возмущений и экспериментальными данными [5], метод Родионова позволяет достичь большей точности в целом нежели схема Годунова. Кроме того, метод Родионова намного точнее показывает положение слабого скачка уплотнения (а именно они и возникают на околозвуковых скоростях), что особенно заметно на грубых сетках — это демонстрируют расчеты обтекания при числе Маха потока равном 0.9.

4. Расчеты стационарного обтекания

Представлены результаты расчетов по программе автора и их сравнение с программой TRAN С.И. Кузьминой [4] и экспериментальными и теоретическими результатами из работы [5]. Рассчитывались следующие крылья: крыло RAE А с профилем RAE 101 и горизонтальное оперение AGARD SMP. Далее в автореферате приведены только результаты расчетов крыла RAE101A.

В задачах аэроупругости для нахождения напряженно-деформированного состояния крыла, как балки, достаточно знать распределение аэродинамических силы и момента по размаху крыла. Поэтому здесь приводятся графики распределения аэродинамических силы Cv и момента Мг относительно передней кромки крыла. Приняты обозначения "TRAN (Годунов)" для расчетов по программе [4] и "Годунов" и "Родионов" для результатов, полученных программой автора методами Годунова и Родионова соответственно.

Используемая в расчетах по программе автора сетка Н-типа состоит из 12 блоков, окружающих крыло. Так, блоки 0 и 1 находятся над и под крылом и являются, соответственно, прилегающими к поверхностям крыла. Остальные блоки дополняют расчетную область по направлению х перед и за крылом и в направлении г. На гранях блоков, соседствующих со свободной границей или с т&аом ставятся граничные условия. Так, на границе z = 0 (блоков с нулевого по пятый) ставится граничное условие отражения, на внешних гра-

ницах всех блоков условие на бесконечности, и на поверхности крыла (блоки О и 1) — условие непротекания. Варьируются такие параметры, как количество ячеек в блоках, размеры блоков относительно корневой хорды крыла и параметры сгущения по определенным направлениям, либо отдельных ребер. Приведенные результаты были получены на сетке, покрывающей 22x24 ячеек по крылу.

Для сравнения в программе TRAN генерироваласть сетка с таким же количеством ячеек.

Расчеты проводились для случаев M — 0.9, а = 0° и M = 0.8, а = 2°.

На рис. 3 и рис. 4 показаны кривые распределения коэффициента давления по хорде крыла на различных процентах размаха крыла. Видно, что уже начиная с 50% размаха крыла программа TRAN значительно теряет точность из-за постоянного уменьшения количества ячеек по хорде крыла. На рис. 5 показаны зависимости интегральных характеристик Су, Мг и положения центров давления (приведенных к местной хорде) от координаты z по размаху крыла.

Кривые на рис. 6 демонстрируют распределение давления по хорде крыла RAE101А при числе Маха 0.9 и нулевом угле атаки. На этих графиках можно отметить разницу в "размазывании" ударной волны методами Годунова и Родионова при различном числе ячеек. Для одинакового количества ячеек у корня, равного 22, зона скачка на 15% размаха крыла в методе [4] составляет 18% хорды, а у метода Родионова — 10%. На 25% размаха зона скачка в методе [4] вырастает до более 20%, в то время как в схеме Родионова около 10%. На 50% размаха крыла по программе TRAN скачок размазывается и теряет ярко выраженную форму, в то время как по методу Родионова зона скачка так же занимает не более 10% хорды.

5. Расчет нестационарного обтекания

Результаты расчетов нестационарного обтекания горизонтального оперения AGARD SMP для перидических колебаний небольшой амплитуды покат заны на рис. 7 и рис. 8. Расчеты были проведены по модифицированной программе TRAN. Для сравнения приведены результаты расчетов из работы [11]

и по исходной программе TRAN.

Рассчитывалось нестационарное обтекание горизонтального оперения, колеблющегося по периодическому закону аа{ 1 — cosuit). Число Маха набегающего потока 0.95. Приведенная частота равна 0.583. Амплитуда угла атаки составляла 0.4°. Приведены реальная и мнимая компоненты коэффициента давления на размахе 69%.

Выводы

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

2) Выполнено сравнительное исследование точности и сходимости новых алгоритмов для стандартных тестовых конфигураций. Проведено сравнение результатов расчетов стационарного обтекания крыла RAE101A и горизонтального хвостового оперения AGARD SMP с программой TRAN С.И. Кузьминой [4] и зарубежными экспериментальными и расчетными данными. Также даны графики сравнения интегральных ха/-рактеристик обтекания, важных для задач аэроупругости — распределений аэродинамических силы, момента и положения центра давления по размаху крыла.

3) Проведены расчеты нестационарного обтекания горизонтального оперения AGARD SMP по схеме Родионова и дано сравнение результатов с программой TRAN и зарубежными расчетными данными.

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

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

6) Время, затрачиваемое на шаг по времени в одной ячейке сетки у разработанной программы по схеме Родионова приблизительно в 5 раз больше, чем у программы [4]. Причиной этому — общая трактовка ячейки сетки и метод более высокого порядка, требующий большего числа операций.

Литература

[1] C.K. Годунов, A.B. Забродин, М.Я. Иванов, А.Н. Крайко, Г.П. Прокопов. Численное решение многомерных задач газовой динамики. Москва, Наука, 1975.

[2] A.B. Родионов. Повышение порядка аппроксимации схемы С.К. Годунова. Журнал вычислительной математики и математической физики, 27(12):1853—1860, 1987.

[3] W.J. Gordon and С.А. Hall. Construction of curvilinear coordinate systems and application to mesh generation. International Journal for Numerical Methods in Engineering, 7(4), 1973.

[4] С.И. Кузьмина. Расчетные исследования трансзвукового флаттера самолета. Ученые записки ЦАГИ, 20(6), 1989.

[5] P.A. Newman and Б.В. Klunker. Computation of transonic flow about finite lifting wings. AI A A Journal, 10(7), 1972.

[6] А.Ю. Азаров. Расчет трансзвукового обтекания профиля в задаче о флаттере по теории потенциала малых возмущений с использованием метода итераций Ньютона. Дипломная работа, Московский Физико-Технический Институт, 1996.

[7] А.Ю. Азаров. Расчет трансзвукового обтекания профиля в задаче о флаттере. Конференция молодых специалистов МФТИ, Жуковский, Ноябрь 1996.

[8] А.Ю. Азаров. О представлении конструкции ЛА и потока единой динамической системой. Конференция молодых специалистов МФТИ, Жуковский, Ноябрь 1997.

[9] А.Ю. Азаров. Использование схемы второго порядка аппроксимации для определения аэродинамических нагрузок при расчете на флаттер в трансзвуковом диапазоне скоростей. Конференция "Современные проблемы аэрокосмической науки", Жуковский, Май 1998.

[10] В.П. Колган. Применение принципа минимальных значений производной к построению конечноразностных схем для расчета разрывных решений газовой динамики. Ученые записки ЦАГИ, Ш(6):68-77,1972.

[11] M.J. Knott. Transonic aeroelastic calculations in both the time and frequency domains. AGARD Report, CP-507, 1992.

Рисунки

Рис. 1: Интервалы индексов блока сетки

• --0 4.6-2,0-0-•/ «--0 2.Ь»4,0-0 5 —Л

; |

1 | 1 !

а 02 0.4 о в а.в 1

ч

Рис. 2: График функции сгущения при различных значениях параметров

1 та. 11 (Г{ •щ

г? V.

_ _ -- — *

!

(а) 15% размаха. Верхняя поверхность.

'¿г'1 т 35». НШ1 т

* 1

/

1

(Ь) 15% размаха. Нижняя поверхность.

(с) 25% размаха. Верхняя поверхность.

(А) 25% размаха. Нижняя поверхность.

1к рщ е! ==

/

„ № Г1»« =гг

У

■54

Г

(е) 50% размаха.

верхность.

{{) 50% размаха. Нижняя поверхность.

Рис. 3: Расчет обтекания крыла ЯАЕ101А. Распределение коэффициента даг вления. Сетка 22 х 24. М„ = 0.8, а = 2°.

ер}

/

ч

ч

А

/

1

УЯ "7

/

Ь^гН

<

/

/

(а) 75% размаха. Верхняя поверхность.

(Ь) 75% размаха. Нижняя поверхность.

те, "тай —.

/

Г ч

1

>

*

//

I — — — -

1

(с) 100% размаха. Верхняя поверхность.

(а) 100% размаха. Нижняя поверхность.

Рис. 4: Расчет обтекания крыла ИАЕ101А. Распределение коэффициента дат вления. Продолжение. Сетка 22 х 24. Моа = 0.8, а = 2°.

ч» Ш:

— N N

N

\

•

ч

\

(а) Эпюра подъемной силы

(Ъ) Распределение аэродинамического момента

1А

\\

\

•.

Ч

V /

"V

(с) Ось центров давления

Рис. 5: Расчет обтекания крыла RA.E1.01A. Интегральные характеристики по размаху крыла. Сетка 22 х 24. М,= 0.8, а = 2°.

(а) 15% размаха крыла

ТК| Иг

-С

у

г — —

/

(Ь) 25% размаха крыла

"Л | -

1 I I

1

Г - — Цс

V

/ !

] 1

(с) 50% размаха крыла (с!) 75% размаха крыла

Рис. б: Расчет обтекания крыла Г1АЕ101А. Распределение коэффициента давления. Сетка 22 х 24. М.= 0.9, а = 0°.

Рис. 7: Горизонтальное оперение AGARD SMP. Реальная компонента коэффициента давления на нижней поверхности. М«, = 0.95, а = 0.4°, к = 0.583

Рис. 8: Горизонтальное оперение AGARD SMP. Мнимая компонента коэффициента давления на нижней поверхности. Мао = 0.95, а = 0.4°, к = 0.583