авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 |
-- [ Страница 1 ] --

РОССИЙСКАЯ АКАДЕМИЯ НАУК

СИБИРСКОЕ ОТДЕЛЕНИЕ

ИНСТИТУТ ВЫЧИСЛИТЕЛЬНЫХ ТЕХНОЛОГИЙ

На правах

рукописи

Слюняев Андрей Юрьевич

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

ОСНОВЕ МОДИФИЦИРОВАННОГО МЕТОДА РАСЩЕПЛЕНИЯ

05.13.18 – математическое моделирование, численные методы

и комплексы программ

Диссертация на соискание ученой степени кандидата физико-математических наук Научный руководитель д.ф.-м.н., проф. В.М. Ковеня Новосибирск, 2009 2 ОГЛАВЛЕНИЕ Введение....................................................................................................................... Глава I. Исходная система уравнений и метод решения...................................... 1.1 Исходная система уравнений Эйлера и Навье-Стокса............................ 1.1.1 Системы уравнений в одномерном случае............................................. 1.2.1 Двумерная система уравнений Эйлера и Навье-Стокса....................... 1.3 Преобразование системы координат............................................................. 1.4 Система уравнений в безразмерном виде..................................................... 1.5 Разностные схемы для решения уравнений Эйлера и Навье-Стокса......... 1.5.1 Разностные схемы для решения одномерных уравнений Эйлера и Навье-Стокса....................................................................................................... 1.5.2. Конечно-разностные для решения двумерных уравнений Эйлера и Навье-Стокса....................................................................................................... 1.6 Аппроксимация производных........................................................................ 1.7 Краевые условия на дробных шагах.............................................................. Глава II. Исследование свойств разностной схемы............................................... 2.1 Задача о распаде произвольного разрыва..................................................... 2.2 Задача о квазиодномерном течении газа в канале....................................... 2.3 Оценки области оптимальной скорости сходимости................................... 2.4 Тестирование алгоритма на двумерных задачах.......................................... Глава III. Моделирование течений газа в канале воздухозаборника со вдувом газа с части поверхности.......................



................................................................... 3.1 Постановка задачи........................................................................................... 3.2 Критерий установления и выход решения на автоколебательный режим 3.3 Расчеты с различными числами Маха........................................................... 3.4 Расчеты с различными числами Рейнольдса................................................ Глава IV. Моделирование сверх- и гиперзвукового обтекания элементов летательного аппарата.............................................................................................. 4.1 Постановка задачи........................................................................................... 4.2 Результаты численных расчетов.................................................................... 4.2.1 Расчеты течений при варьировании чисел Маха.................................. 4.2.2 Расчеты течений при варьировании угла атаки..................................... 4.2.3 Течение газа около элементов ЛА в расширенной области................. 4.2.4 Влияние краевых условий на характеристики течения газа в канале Выводы..................................................................................................................... Литература............................................................................................................... Введение Применение аппарата математического моделирования предоставляет современным исследователям возможность предварительного изучения физических процессов в широком диапазоне параметров. Стремление к моделированию явлений, наиболее полно соответствующих реальному физическому процессу, приводит к необходимости использования наиболее сложных математических моделей. В задачах аэродинамики наиболее полной моделью течения вязкого сжимаемого теплопроводного газа является газодинамическая модель, описываемая системой уравнения Навье-Стокса.

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

Разнообразие физико-математических моделей в аэродинамике и соответствующих им типов уравнений обусловило необходимость применения различных численных алгоритмов их решения. Известными численными алгоритмами решения задач аэродинамики являются конечно-разностные, конечно-объемные, конечно-элементные методы и другие (например, метод частиц в ячейках, метод статистических испытаний [2,8,16,28,29,30,42,48,52,55,56,68,98-100]). Наиболее широкое распространение в механике сплошных сред из-за своей универсальности получил метод конечных разностей.





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

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

Явные схемы для решения задач газовой динамики были предложены в работах С.К.Годунова, А.В.Забродина и Г.П.Прокопова [16,17], В.В.Русанова [58]. Развитие этого метода для решения стационарных и нестационарных задач проведено А.Н.Крайко, Н.В.Михайловым, М.Я.Ивановым [16,26,27], В.П.

Колганом, А.П.Косых, А.Н.Минайлосом [43,44]. Явные разностные схемы типа предиктор-корректор рассматривали P.D.Lax, R.W.MacCormak [91,92], G.Moretti [95], P.Kutler и H.Lomax [87,89,90]. Этот класс схем можно рассматривать как разновидность метода расщепления, предложенного Н.Н.Яненко [73]. Другим подходом к решению задач аэродинамики является метод частиц в ячейках предложенный P.H.Harlow, R.A.Gentry, B.I.Daly [72,81].

Некоторые модификации этого метода рассматривались Н.Н.Яненко, Н.Н.Анучиной, В.Е.Петренко, Ю.И. Шокиным [3,51,74]. Простота реализации явных схем обеспечила широкое их применение для целого класса задач аэродинамики. Однако жесткие ограничения на устойчивость таких схем увеличивали время расчета задач при попытках получить более точное решение на измельченной сетке. При таких ограничениях на временной шаг использование явных разностных схем может оказаться неэффективным (особенно в случае решения стационарных задач методом установления).

Неявные разностные схемы предложены в работах Н.Н.Яненко [73], R.Beam, R.F.Warming, H.Macdonald, I.L.Steger, W.R.Briley, R.Kutler и T.H.Pulliam [9,77,78,79,88,93,96,103,105], А.И.Толстых [50,67,69,70]. Они основаны на приближенной факторизации многомерных операторов с последующей линеаризацией нелинейных членов. Реализация таких схем сводится к векторным прогонкам вдоль каждого пространственного направления. В трехмерном случае схемы приближенной факторизации теряют свойство безусловной устойчивости. Однако использование векторных прогонок в случае многомерных задач может оказаться неэффективным ввиду необходимости проводить обращение матриц для нахождения коэффициентов прогонки. Такое ограничение может существенно увеличить абсолютное время решения задачи при увеличении числа узлов в расчетной области.

Неявные разностные схемы, основанные на расщеплении уравнений по физическим процессам и пространственным направлениям, предложены в работах Н.Н.Яненко и В.М.Ковени [39,41,42,76]. Они обладают свойством безусловной устойчивости и реализуются на дробных шагах скалярными прогонками, что делает их экономичными. Многочисленные проведенные расчеты показали их достаточную точность и эффективность на решении стационарных задач. Для численного решения нестационарных задач газовой динамики в лагранжевых координатах в работах А.А.Самарского и Ю.П.Попова предложен класс неявных разностных схем, основанный на интегро интерполяционном методе [53,59,60].

Расчеты вязких течений в приближении полных уравнений Навье-Стокса сопряжены со значительными трудностями. Как уже говорилось, особенностью таких течений является наличие узких областей (порядка O(1/ Re), O(1/Re) ), в которых происходит резкое изменение значений искомых функций. Для получения достоверного результата необходимо обеспечить попадание в такую область достаточного количества узлов расчетной сетки. Одним из вариантов решения проблемы является использование адаптирующихся к решению подвижных сеток, такие подходы рассматривались в работах А.И.Толстых [69, 71], I.F.Tompson, I.L. Steger [80,82,101,102,106], Н.Н.Яненко, В.Д.Лисейкина, Н.Т.Данаева и других [18,24,42,46,49,54,61,62,75]. Их применение совместно с использованием схем повышенного порядка аппроксимации как на трехточечном (А.И.Толстых [50,67,70], R.F.Warming, I.L.Steger [78,93,96,103,105]), так и на расширенном шаблоне (например, работы В.В.Русанов [57], Б.В.Балакин [5], P.F.Warming, P.Kutler, H.Lomax [22, 94,104] совместно с применением неявных разностных схем дает возможность получать решения задач аэродинамики в приближении полной системы уравнений Навье-Стокса.

Одним из актуальных и бурно развивающихся направлений вычислительной математики является теория и технология конструирования адаптивных высокоточных разностных методов, обеспечивающих сохранение монотонных свойств решения. Наиболее широко известными и признанными в аэродинамике являются схемы с полной ограниченной вариацией решения, которые могут подстраиваться под решение на различных участках и, таким образом, сохранять высокий порядок аппроксимации, гасить высокочастотные гармоники решения и сохранять монотонность решения. Одной из ключевых работ в области построения нелинейных адаптивных конечно-разностных методов является построенная A.Harten теория TVD-схемы [83]. Эта работа стала основополагающей для дальнейшего развития и строгому обоснованию высокоточных конечно-разностных схем.

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

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

Предложен алгоритм второго порядка аппроксимации по пространственным переменным.

2. Создана программа расчета течений вязкого сжимаемого газа около тел со сложной геометрией.

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

4. Исследованы сверхзвуковые течения газа в канале воздухозаборника с наличием вдува газа с части поверхности канала. Получены основные закономерности течения, изучено влияние чисел Маха и Рейнольдса на характер течения. Выявлены основные закономерности влияния параметров потока на размер отрывных зон.

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

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

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

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

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

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

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

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

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

Глава III посвящена моделированию течений газа в канале воздухозаборника со вдувом газа с части поверхности. Пункт 3.1 содержит постановку задачи. Здесь приводится описание расчетной области, дана постановка краевых условий. Далее приводится информация о выборе начального распределения, скорости сходимости и параметрах разностной схемы, указывается диапазон изменения параметров газа набегающего потока и параметров вдуваемого в канал газа. Вид выбранного критерия установления, а также особенности поведения численного алгоритма при решении задачи методом установления обсуждается в п. 3.2. Пункт 3.3 посвящен описанию решения задачи при варьировании числа Маха набегающего потока, дано описание характера течения, а так же описаны влияние изменения числа Маха на структуру течения, а так же зависимость между изменением числа Маха и размерами зон возвратного течения, образующихся на верхней и нижней стенке канала. Отдельные результаты моделирования течения газа в канале сравниваются с аналогичными исследованиями других авторов и экспериментальными данными. Здесь же дано сравнение численных решений, полученных с первым и вторым порядком точности, обосновывается эффективность применения оператора второго порядка аппроксимации.

Описание влияния числа Рейнольдса на структуру течения дано в п. 3.4.

Основные выводы по проведенному исследованию даны в конце главы.

В главе IV объединены результаты серии расчетов течения сверх- и гиперзвукового газа около элементов модельного летательного аппарата. В п.

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

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

На основе проведенного исследования в п.п. 4.2.4 изучается влияние краевых условий для температуры на структуру течения. Здесь показано, что изменение типа краевого условия для температуры с тепловой изоляции стенки на условие теплоотвода приводит к перестроению структуры течения исчезновению областей отрывного течения в канале воздухозаборника. Краткие выводы по серии проведенных расчетов приведены в конце главы.

Основные результаты диссертации сформулированы в разделе Выводы.

По мере получения результаты диссертации докладывались на:

1. Всероссийской конференции «Исследования и перспективные разработки в авиационной промышленности», Москва, 2. Конференции молодых ученых по математическому моделированию и информационным технологиям, Красноярск, 3. Международном семинаре по вычислительным технологиям Россия Казахстан, Новосибирск, 4. Международной конференции по методам аэрофизических исследований (ICMAR), Новосибирск, 5. Всероссийской конференции по вычислительной математике, Новосибирск, 6. Конференции молодых ученых по математическому моделированию и информационным технологиям, Новосибирск, 7. VII Международной конференции по неравновесным процессам в соплах и струях (NPNJ'2008) Крым, Алушта, 8. Международной конференции по методам аэрофизических исследований (ICMAR), Новосибирск, 9. Объединенном семинаре ИВТ СО РАН, кафедры математического моделирования НГУ и кафедры вычислительных технологий НГТУ «Информационно-вычислительные технологии» под руководством академика Ю.И. Шокина и проф. В.М. Ковени, Новосибирск, Результаты диссертации опубликованы в работах [34-38,63-66,85,86].

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

Работа является составной частью плановых работ Института вычислительных технологий СО РАН «Информационно-вычислительные технологии в задачах поддержки принятия решений». Работа выполнена при частичной финансовой поддержке РФФИ (гранты № 05-01-00146, 08-01-00264).

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

1.1 Исходная система уравнений Эйлера и Навье-Стокса 1.1.1 Системы уравнений в одномерном случае Представим систему уравнений Навье-Стокса вязкого сжимаемого теплопроводного газа в дивергентной форме в виде (см., например, [47]) U W 0, (1.1) t x v v, E (e v ).

где U v, W v p x E T v v( E p) x x Для замыкания системы уравнений (1.1) зададим уравнения состояния p p (, e), e e(T ) cvT (1.2) и зависимости вязкости и теплопроводности, от температуры (T ), (T ).

При 0 система уравнений (1.1) описывает уравнения газовой динамики – уравнения Эйлера. Здесь - плотность, v - скорость газа, p - давление, T v температура газа, E (e ) - полная энергия газа, e cvT - внутренняя удельная энергия газа, - коэффициент динамической вязкости, коэффициент теплопроводности, - удельная теплоемкость газа при cv постоянном объеме.

Наряду с консервативной формой представим исходные уравнения (1.1) в недивергентной и предельно-дивергентной формах f f U A 1 W, A B f F или. (1.3) t t f Отметим эквивалентность записи уравнений Навье-Стокса в дивергентном (1.1) и в недивергентном (1.3) виде. Выбор вектора f задает форму матричных операторов B, что дает возможность рассматривать различные классы разностных схем. Для систем уравнений Эйлера и Навье-Стокса искомыми функциями обычно выбирают плотность, скорость, а также давление или температуру, что определяется заданием краевых условий. Однако вид матричных операторов B для различных компонент вектора f существенно различается [42]. Выберем в качестве искомых функций плотность, скорость и давление газа, тогда для уравнения состояния совершенного газа p ( 1) e матрица B принимает вид x v 0 0 1 v, B 0, F, (1.4) f v a x x x x p v 1 x 0 vk c x x x x p p 1 p, где a 2, c2,k. Уравнение неразрывности записано в cv T cv T cv дивергентном виде, что позволяет упростить вид матричного оператора B.

1.2.1 Двумерная система уравнений Эйлера и Навье-Стокса Система уравнений Навье-Стокса для вязкого сжимаемого теплопроводного газа в дивергентном векторном виде в декартовой системе координат для двумерного случая может быть записана в следующем виде:

U F E G H, (1.5) t x1 x2 x1 x где u v 2 uv u u p U, E( U ), F(U) 2, uv v p v ( E p)u ( E p )v E (1.6) 0 0 x1 x1 x1 x G (U), H (U) x x, x1x u x x v x x qx u x x v x x qx 11 1 12 12 вектор искомых функций, газодинамические и вязкие потоки по направлениям x1 и x2 соответственно. Компоненты тензора вязких напряжений и тепловые потоки задаются соотношениями:

u v v u u v 2 x x (2 ), x2 x2 (2 ), x1x2 ( ), x1 x2 x2 x1 x1 x 3 T T qx1 k, qx2 k.

x1 x Здесь - плотность, u, v - компоненты вектора скорости V (u, v)T в V направлении x1 и x2, p - давление, T - температура газа, E (e ) - внутренняя удельная энергия газа, полная энергия газа, e cvT коэффициент динамической вязкости, k - коэффициент теплопроводности, cv удельная теплоемкость газа при постоянном объеме. При записи уравнений принято соотношение Стокса, что вторая вязкость связана с коэффициентом динамической вязкости соотношением.

Система уравнений Навье-Стокса замыкается уравнением состояния совершенного газа:

p ( 1) e, (1.7) и законами зависимости вязкости и теплопроводности от температуры T 0 ( ), Pr=cp, (1.8) T0 k где для разных газов параметр изменяется в пределах 0.5 1.0, Pr - число Прандтля, c p - удельная теплоемкость газа при постоянном давлении. В предположении, что 0 система уравнений (1.5) переходит к уравнениям газовой динамики.

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

f (, u, v, p )T. (1.9) Тогда система уравнений примет следующий вид:

f f Bj F, (1.10) t j 1 x j где 1 j uj j j 0 uj Bj, j 0 0 uj 0 uj p 1 p j j 0 1 ( x1x1 x1x2 ) x1 x F x1x2 x2 x2, ( ) x x T T ( 1)[ x k x x k x ] 1 1 2 здесь ij - символ Кронекера.

1.3 Преобразование системы координат Как правило, решение задач аэродинамики требуется найти не в прямоугольной области, а в областях более сложной формы с криволинейными границами. Это связано, например, с тем, что образующая несущей поверхности летательного аппарата имеет сложную форму. Решение задач аэродинамики в криволинейных областях конечно-разностными методами сопряжено с большими трудностями при построении разностной сетки и аппроксимации функций на таких сетках, поэтому обычно при решении уравнений Навье-Стокса в непрямоугольных областях проводят замену переменных [например, 12,40,42,97]. Замена переменных позволяет преобразовать криволинейную расчетную область в единичный квадрат, в котором потом удобно строить равномерную сетку. Естественно, что на равномерной сетке удобно проводить конечно-разностную аппроксимацию непрерывных разностных операторов. Равномерной сетке в единичной области будет соответствовать неравномерная сетка в исходной криволинейной области. Кроме того, можно подобрать такие преобразования координат, которые позволят сгустить сетку в исходной области в зоне больших изменений градиента решения и таким образом учесть физические особенности решения и получить более точные численные результаты. Некоторые подходы построения адаптивных сеток для задач аэродинамики изложены в [12,42,46].

Введем невырожденное преобразование координат в виде:

q1 q1 ( x1, x2 ), q2 q2 ( x1, x2 ). (1.11) Тогда система уравнений Навье-Стокса (1.5) в новых переменных (1.11) вновь может быть записана в консервативном виде:

U E F G H (1.12), t q1 q2 q1 q здесь вектор искомых функций, газодинамические и вязкие потоки в старых и новых переменных связаны следующими соотношениями:

U (, u, v, E )T, J U q1 U q 1 uU q1 pq1x1 q2 x1 1 uU q2 pq2 x q1x1 q1x2 q2 x E E F,F E F, J vU q1 pq1x2 J vU q2 pq2 x J J J J ( E p)U q ( E p )U q 1 0 1 q1x1 x1x1 q1x1 x1x q1x1 q1x2 G G H, J q1x1 x1x2 q1x1 x2 x J J (u x x v x x ) q1x (u x x v x x )q1x Tq (1.13) 11 12 1 12 22 0 q q2 x2 x1x q2 x1 G q2 x2 H 1 2 x1 x1x1 H.

q q J J J 2 x1 x1x 2 x2 x2 x (u x x v x x )q2 x (u x x v x x )q2 x Tq 11 12 1 12 22 При записи уравнений введены обозначения:

U q1 uq1x1 vq1x2, U q2 uq2 x1 vq2 x2, T T T T Tq1 k[( q1x1 q2 x1 )q1x1 ( q1x2 q2 x )q1x2 ], q1 q2 q1 q2 T T T T Tq2 k[( q1x1 q2 x1 )q2 x1 ( q1x2 q2 x )q2 x2 ], q1 q2 q1 q2 2 v 4 u 2 v 4 u x x [( q2 x q1x1 q1x2 ) ( q2 x )], 3 q2 3 q1 3 q1 3 q 4 v 2 u 4 v 2 u x x [( q1x2 q1x1 ) ( q2 x2 q2 x )], 3 q1 3 q1 3 q2 3 q2 u v u v x x [ q1x2 q1x1 q2 x2 (1.14) q2 x )], q1 q1 q2 q2 q j q jxi, i, j 1,2, xi J - якобиан преобразования координат ( x1, x2 ) в (q1, q2 ), задаваемый в виде:

q1 q x x1 q q2 x J det J det det 1x q1x1 q2 x2 q1x2 q2 x1. (1.15) q1x q2 x q1 q2 2 x x 2 Переходя к недивергентной форме, представим систему уравнений Навье Стокса в новых переменных в виде:

f f F, Bj (1.16) t j 1 q j q jx q jx vq jx uq jx1 1 q jx uq jx1 vq jx 0 Bj, q jx uq jx1 vq jx 0 uq jx1 vq jx pq jx pq jx 1 0 1 (q1x x x q1x x y q2 x x x q2 x x x ) 1 11 2 1 1 11 2 1 (q q q q ).

F 1 x1 x1 x2 1 x2 x2 x2 2 x1 x1 x2 2 x2 x2 x ( 1)[(u v )q (u v )q T x1 x1 x1 x2 1 x1 x1 x2 x2 x2 1 x2 q (u x1x1 v x1x2 )q2 x1 (u x1x2 v x2 x2 )q2 x2 Tq2 ] Здесь обозначения x1x1, x1x2, x2 x2, Tq1, Tq2 приведены в (1.14).

Заметим, что для системы уравнений Навье-Стокса в новых переменных в недивергентном векторном виде, вектор искомых переменных не изменит своего вида по сравнению с видом (1.9).

Вообще говоря, использовать при численных расчетах коэффициенты q jx i неудобно, так как мы не знаем шага по сетке в исходной области, поэтому можно воспользоваться следующей связью коэффициентов:

q1x1 Jx2 q2, q1x2 Jx1q2, q2 x1 Jq2 q1, q2 x2 Jq1q1, которые могут быть вычислены на сетке в единичном квадрате.

Как и для уравнений (1.5) и (1.10) система уравнений (1.12) замыкается уравнением состояния совершенного газа (1.7) и законами зависимости вязкости и теплопроводности от температуры (1.8).

Отметим связь между дивергентной и недивергентной формами записи уравнений Навье-Стокса. Учитывая, что вектор функций U U (f ), т.е. зависит от вектора функций f, получим:

U U f F E G H ( ).

t f t x1 x2 x1 x U F E G H и A Введем обозначения W ( ). Тогда система f x y x y уравнений (1.5) может быть представлена в предельно-дивергентной форме:

f A 1 W. (1.17) t где матрицы A и A 1 для уравнения состояния (1.7) равны:

1 0 1 0 0 u u 0 0 0, A.

A v 0 v 2 u v u v 2 1 ( 1) u v 2 ( 1)u ( 1)v Из (1.10) следует, что:

f f B j F.

t x j j Как следствие, получим связь между дивергентной и недивергентной формами записи уравнений Навье-Стокса:

f B j F A 1 W. (1.18) x j j Отметим, что для введенной замены переменных (1.11) матрица A изменит свой вид на A J A.

1.4 Система уравнений в безразмерном виде Использование уравнений в размерном виде неудобно для практических приложений, поэтому обычно переходят к безразмерным переменным, осредненным на некоторые их значения. Введем характерный размер L, 0, скорость u0, температуру T0, равные их значениям в плотность x x tu u v T невозмущенном потоке. Тогда x1 1, x2 2, t 0, u, v, T, u0 u0 T L L L p. Здесь величины со знаком есть безразмерные величины. После p 0u подстановки безразмерных переменных, исходные уравнения (1.12) перепишутся в виде (для простоты знак опущен):

U E F 1 G H ( ), (1.19) t q1 q2 Re q1 q компоненты векторов E и F не изменят своего вида по сравнению с (1.13), изменит свой вид только коэффициент при градиенте температуры:

, k ( 1) M 2 Pr 0 L0u0 u здесь Re - число Рейнольдса, Pr c p 0 - число Прандтля, M 0 k0 c число Маха, c p - удельная теплоемкость газа при постоянном давлении, p c - скорость звука. Уравнения (1.16) примут вид:

f f 1 F.

Bj (1.20) t j 1 q j Re Что касается коэффициента при градиенте температуры, то он будет иметь вид:

.

k ( 1) Pr Теперь будем полагать, что рассматривается безразмерная система уравнений Навье-Стокса.

1.5 Разностные схемы для решения уравнений Эйлера и Навье-Стокса 1.5.1 Разностные схемы для решения одномерных уравнений Эйлера и Навье-Стокса Введем в расчетной области разностную сетку k с равномерным шагом по пространству h и шагом по времени. Сеточные функции U h и f h согласно (1.4) и определим во всех узлах вычислительной сетки. Дифференциальный оператор x будем аппроксимировать разностным оператором в узлах сетки с порядком O (h k ), где k - порядок аппроксимации пространственных производных. При аппроксимации разностных операторов конечно k разностными аналогами будем полагать для (несимметричная противопотоковая аппроксимация):

если v 0, то,, если v 0, то,.

При симметричной аппроксимации производных для k 2 полагаем ( ) 2.

Здесь ( I T1 ) h, T1 fl f l 1 - оператор сдвига. Вторые производные B в матричных операторах аппроксимируем симметричными x x операторами со вторым порядком O(h 2 ). Вектор потоков W x в системе уравнений (1.1) также аппроксимируем симметричным оператором Wh со вторым порядком. С учетом введенных аппроксимаций разностный оператор Bh аппроксимирует дифференциальный оператор B из (1.4) с порядком O (h k ) по формулам:

v a 2, b1, b2 k.

Bh 0 v b1 (1.21) 0 v b c При построении разностной схемы аппроксимацию первых производных в операторе Bh будем задавать несимметричной с учетом знака скорости с первым порядком, а в операторе правых частей Whn – симметричной со вторым порядком. Отметим, что при несимметричной аппроксимации первых производных в операторе Bh для получения безусловно устойчивых схем аппроксимацию членов с давлением в уравнении движения необходимо выбирать по сопряженным к конвективным членам формулам (см., [42]).

Известно [15], что симметричная аппроксимация оператора j W j в разностной схеме приводит к осцилляциям решения. В работах [31,33] предложен разностный оператор второго порядка аппроксимации со сглаживающим оператором для недивергентного случая. В данной работе для решения стационарных задач предложена модификация этого оператора для уравнений, записанных в дивергентной форме (1.1):

Wl 1 Wl 1 W 2Wl Wl a sign(v) l Wl, (1.22) 2h 2h Wl 1 2Wl Wl, где | Wl 1 Wl | | Wl Wl 1 | 0, если | W W | | W W | 0.

l 1 l l l Перейдем от непрерывной модели (1.3) к ее конечно-разностному аналогу.

Заменим её конечно-разностной аппроксимацией и введем аппроксимацию операторов их конечно-разностными аналогами. Тогда разностная схема с x весами:

f n1 f n Bn ( f n1 f n ) F n, аппроксимирует систему уравнений (1.3) с порядком O (h m ). Перепишем разностную схему в каноническом виде:

f n1 f n (I B ) A 1 W n, n (1.23) где с учетом соотношения (1.3) правая часть схемы (1.23) аппроксимирована в предельно-дивергентной форме. Схема (1.23) аппроксимирует систему уравнений (1.4) с порядком O( h 2 ), а при установлении установления (при f n1 f n 0 ) аппроксимирует исходные уравнения в выполнении условия дивергентном виде (1.1). с порядком O(h 2 ). В линейном приближении схема (1.23) безусловно устойчива при 0.5 [29]. Матричный оператор Bn, h аппроксимирующий оператор B, перепишется в следующем виде:

u n n 0 (1.24) Bn 0 u n b1n b2n.

0 u n b2n c2n Для реализации схемы (1.23) рассмотрим эквивалентную ей схему в дробных шагах:

n (A 1 )h Whn, n (1.25) (I Bn ) n1 n, (1.26) f n 1 f n n 1. (1.27) Разностная схема (1.23) в дробных шагах (1.25)-(1.27) реализуется векторными прогонками на каждом дробном шаге и требует обращения матриц размерностью 3 3.

Для построения экономичных разностных схем, реализуемых скалярными прогонками, подобно [32] представим матричный оператор Bh в виде суммы двух операторов:

Bh B1 Bh (1.28) h где разностные операторы имеют вид:

0 0 0 v B1 0 0 b 2, B2 0 v b1 0 (1.29) h h 0 0 c 2 v b Для численного решения системы уравнений Навье–Стокса рассмотрим схему приближенной факторизации:

f n1 f n (I B1n )(I Bh n ) (A 1 ) n Whn (1.30) h h или эквивалентную ей схему в дробных шагах:

n (A h1 ) n Whn, (I B1n ) n1/2 n, h (1.31) (I B2 n ) n1 n1 2, h f n1 f n n1, которая аппроксимирует уравнения (1.1) с порядком O( h 2 ) при всех, при установлении схема консервативна, т.е. Whn 0 в силу того, что A 1 0.

Очевидно, что разностная схема (1.31) реализуется скалярными прогонками на втором дробном шаге, т.к. матрица B2 содержит ненулевые h элементы лишь на главной диагонали. Покажем, что на первом дробном шаге схема реализуется также скалярными прогонками. Запишем систему уравнений, решаемую на первом дробном шаге:

1 n n k u b p un, 2n 2 (1.32) 1 1 1 n n n n p ( p n 1 u u n 1 p b2n p ) p.

k k n 2 2 2 n Из первого уравнения выразим u, подставим во второе уравнение системы n (1.32) и получим уравнение, разрешенное относительно p 2 :

n k (I (u b p b )) p p 1 un, n k n n k 2n n k 1 2 которое может быть решено методом скалярных прогонок. После вычисления 1 n n p, значение u вычисляется явно по формуле:

2 1 n n k u b p n 2n 2 u Введение расщепления (1.28) приводит к появлению в разностной схеме (1.30) дополнительных членов вида D 2 2 B1n Bh n – диссипативных членов h второго порядка малости. Поэтому расщепление матричных операторов Bh необходимо выбирать таким образом, чтобы построенная разностная схема (1.30) сохраняла свойство безусловной устойчивости и имела минимальную диссипацию – минимальное число дополнительных членов. В работе [42] было предложено расщепление одномерного оператора Bh из (1.21) в виде расщепления по физическим процессам в форме 0 0 v 2 B1 a 0 b 2, B2 0 v b1 0, h h 0 v b 0 c 0 данное расщепление позволило свести реализацию разностной схемы вида (1.30) на дробных шагах к четырем скалярным прогонкам, а диссипативная матрица содержала дополнительные члены в уравнениях движения и энергии.

Для вектора f (, v, p)T из (1.4) существует единственное расщепление Bh на их сумму в виде (1.29), для которого диссипативная матрица минимальна [29].

При использовании расщепления такого вида дополнительный диссипативный член D возникает лишь в одном уравнении – уравнении энергии.

Спектральным методом показано [29], что для линеаризованных уравнений схема (1.30) обладает свойством безусловной устойчивости при 0.5 O( ).

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

1.5.2. Конечно-разностные для решения двумерных уравнений Эйлера и Навье-Стокса Введем в расчетной области разностную сетку k с шагами по пространству h j и шагом по времени. Сеточные функции U h и f h согласно (1.6) и (1.9) определим во всех узлах расчетной сетки, совпадающие с исходными функциями в узлах сетки. Дифференциальные операторы x j будем аппроксимировать разностным оператором j в узлах сетки с порядком h j max(hij ), i 1, N j, j 1, N, O(h m ), m – порядок аппроксимации j пространственных производных, N – размерность задачи.

Вернемся к системе уравнений (1.20) в недивергентном векторном виде.

Приведем эту систему уравнений:

f f F.

Bj (1.33) t j 1 q j Перейдем от непрерывной модели (1.33) к конечно-разностному аналогу. Для этого производную по времени заменим её конечно-разностной аппроксимацией и введем аппроксимацию операторов их конечно x j разностными аналогами. Тогда разностная схема с весами:

f n1 f n n B j ( f n1 f n ) F n, 1 (1.34) j аппроксимирует систему уравнений (1.33) с порядком O(h m ). Перепишем разностную схему в каноническом виде:

n fn n f 1 n (I B j ) A W, (1.35) j где с учетом соотношения (1.17) правая часть схемы (1.35) аппроксимирована в предельно-дивергентной форме.

Разностная схема (1.35) безусловно устойчива при 0.5 O( ). При f n1 f n 0, разностная схема установлении, т.е. при выполнении условия 1 n принимает вид A W 0, т.е. аппроксимирует исходные уравнения в дивергентном виде, так как матрица A не вырождена. Решение может быть получено с помощью матричной прогонки, требующей большого числа операций - примерно N 3 в каждом узле расчетной сетки, где N j – число узлов j по пространственному направлению. Поэтому, используя схему (1.35) в качестве базовой, построим более экономичную схему, при этом согласуясь с требованиями устойчивости и консервативности разностных схем.

Рассмотрим разностную схему, в которой решение уравнений на дробных шагах может быть получено более экономичными методами. В работе [32] предложено расщепление операторов в форме, когда влияние дополнительных членов минимально. Такие схемы авторы работы [32] назвали схемами с оптимальным расщеплением. Дадим обобщение схемы из [32] на уравнение Навье-Стокса для двумерного случая. В переменных f f (, u, v, p ) (, u, v, p )T существует единственное расщепление матричных операторов:

1 (1.36) B j B j B j, где 0 0 0 1 k 0 q j 0 jx Bj, (1.37) 1 k 0 q j 0 jy 0 pq k (uq jx vq jy ) j b pq jy kj k jx j q jx kj u q jy kj v 0 (uq jx vq jy ) kj b 0 0 Bj 2, (1.38) (uq jx vq jy ) kj b 0 0 0 0 элементы bii содержат повторные производные по пространственному направлению из компонентов тензора вязких напряжений и градиента температуры.

Учитывая расщепление (1.36), приближенно факторизуем исходный 2 n оператор (I B j ) по формуле:

j n 1n 2n (I B j ) (I B j )(I B j ). (1.39) j 1 j Тогда разностная схема примет вид:

n fn 2n f 1 n 1n (I B j )(I B j ) (A ) n W h. (1.40) j Она аппроксимирует исходную систему уравнений с порядком O(h m ), а при установлении – с порядком O(h m ), в линейном приближении при 0. является абсолютно устойчивой для двумерного случая [34].

Для реализации разностной схемы (1.40) рассмотрим эквивалентную ей схему в дробных шагах:

1 n n (A ) n W h, (1.41) h 1n (I B1 ) n1/4 n, 2n (I B1 ) n1/2 n1/4, (1.42) 1n (I B2 ) n3/4 n1/2, 2n (I B2 ) n1 n3/4, f n1 f n n1. (1.43) Очевидно, что разностная схема (1.41) – (1.43) реализуется скалярными прогонками на 2 и 4 дробных шагах, т.к. матрицы B j содержат ненулевые элементы лишь на главной диагонали. Покажем, что на 1 и 3 дробных шагах схема реализуется также скалярными прогонками. Запишем систему уравнений, решаемую на первом дробном шаге:

1 n n k u q1x p un, 4 n 1 n n k v q vn, 4 n 1y 1 p (1.44) 1 1 n n n ( p n q1x p n q1 y k k 4 4 p 1u 1v 1 n n (u q1x v q1 y ) p b ) p.

n n k n n 4 44 p 1 n n Из первых двух уравнений выразим u и v 4, подставим их в третье n уравнение системы (1.44) и получим однородное относительно p уравнение:

1 1 n k k (I ((u n q1x v n q1 y )1 b44 p n (q1x 1 q ))) p q q1 y 1 k n k k n 1y 1 n 1x p (1 un 1 vn ), n k k n которое может быть решено скалярными прогонками. После вычисления p, 1 n n значения u и v вычисляются явно по формулам:

4 1 n n k u q1x p, n 4 n u 1 n n k v q1 y p.

n 4 n v Схема с оптимальным расщеплением является достаточно экономичным алгоритмом (по числу операций на один узел расчетной сетки) для численного решения систем уравнения Навье-Стокса. Однако при обобщении ее на многомерный случай факторизованный оператор содержит дополнительные члены порядка O( 2 ). Представим схему (1.40) в виде:

n fn 2 in f 1 n (I B j ) (A ) n W h, (1.45) j 1 i где матрица есть матрица добавочных членов, возникающих вследствие расщепления и факторизации исходного оператора:

n fn 1n 2 n 1n 2 n 1n 2 n f in in O( 3 ), B B j. (1.46) 2 2 (B1 B1 B2 B2 B B ) j При нахождении стационарного решения методом установления шаг стремятся выбрать как можно более большим, чтобы максимально быстро придти к установлению, но чем больше мы выбираем шаг, тем больше становится и ошибка факторизации на промежуточных шагах. Это может приводить к понижению скорости сходимости к стационарному решению. При решении нестационарных задач влияние ошибки факторизации (1.46), может оказаться значительным для сильно меняющихся по времени процессов.

Выбор вида вектора искомых функций влияет на вид матричного оператора B. Зададимся целью сохранить уравнения неразрывности и движения в предельно-дивергентной форме на дробных шагах. Оставаясь в тех же переменных f f (, u, v, p) (, u, v, p)T, согласно принятой модификации вида уравнений расщепим исходный оператор B j по формуле:

1 (1.47) B j B j B j, где оператор B j не изменит своего вида по сравнению с оператором (1.37):

0 0 0 1 k 0 q j 0 jx Bj, (1.48) 1 k 0 q j 0 jy 0 pq k (uq jx vq jy ) j b pq jy kj k jx j а оператор B j изменит свой вид на следующий:

q jx kj u q jy kj v 0 q jx kj u q jy kj v b 0 0 2j B, (1.49) q jx kj u q jy kj v b 0 0 0 0 элементы bii содержат повторные производные по пространственному направлению из компонентов тензора вязких напряжений и градиента температуры.

Аналогично (1.39), приближенно факторизуем исходный оператор (I Bnj ) по формуле, учитывая расщепление (1.47):

j n 1n 2n (I B j ) (I B j )(I B j ). (1.50) j 1 j Тогда разностная схема примет вид:

n fn 2n f 1n (I B j )(I B j ) (A 1 ) n Whn. (1.51) j Она аппроксимирует исходную систему уравнений с порядком O(h m ), а при установлении – с порядком O(h m ), в линейном приближении при 0. является абсолютно устойчивой для для двумерного случая [34].

Для реализации схемы (1.51) рассмотрим эквивалентную ей схему в дробных шагах:

n (A 1 )h Whn, n (1.52) 1n (I B1 ) n1/4 n, 2n (I B1 ) n1/2 n1/4, (1.53) 1n (I B2 ) n3/4 n1/2, 2n (I B2 ) n1 n3/4, f n1 f n n1. (1.54) В связи с изменением расщепления, теперь на первом и третьем дробном шаге счет ведется относительно невязок скоростей u и v, а на втором и четвертом дробных шагах счет ведется относительно невязок u и v.

Покажем, что для организации счета по разностной схеме (1.51) в дробных шагах (1.52) - (1.54), нам не нужно дополнительно хранить значения u и v.

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

0, u p t q1x 0, t q u 1 p q1x 0, v q p 0, t q 1y t q v 1 p q1 y 0.

t q1 p u v p p(q1x q1 y ) (uq1x vq1 y ) b44 0.

t q1 q1 q далее из уравнений получим:

u 1 p q1x 0, q t u p 0, u u, q1x t q1 t t v v v 1 p 0, q1 y, q t t t v p 0.

q1 y t q1 затем ( u ) n1/ 4 ( u ) n u n1/ 4 u n n, ( v) n1/ 4 ( v) n v n1/ 4 v n n, и отсюда необходимые значения невязок для второго дробного шага:

1 1 1 n n n n u u, v v, n n 4 4 4 далее после расчета на втором дробном шаге, необходимо получить значения 1 n n u и v для третьего дробного шага:

2 1 n n u 2 n 1 v n u n, v 2 n, аналогично для четвертого дробного шага, получим значение:

3 3 3 n n n n u u, v v.

n n 4 4 4 Чтобы получить значения un1, vn1, необходимые для перехода на новый n 1 ый временной слой, необходимо значения u 1, v1, полученные после n n выполнения дробных шагов, разделить на новое значение плотности - n1 в каждом узле вычислительной сетки.

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

1.6 Аппроксимация производных Аппроксимация производных x в разностных схемах проводилась j разностным оператором j с порядком O (h k ), где j – пространственное направление, по следующим формулам:

При k 1 (противопотоковая схема):

(U q j K )i (U q j K )i j (U q j K )i (U q j K )i, если U q j 0;

j hj (U q j K )i 1 (U q j K )i j (U q j K )i (U q j K )i, если U q j 0;

(1.55) j hj (, u1, u2, E p )T, K J pi 1 pi j pi j pi, если U j,i 0;

hj и (1.56) p pi j pi = pi i, если U j,i 0 (j 1,2, i 1, N j ) j hj где U q j и J определены в (1.13) и (1.15) соответственно. На дробных шагах аппроксимация первых производных выполнена по противопотоковой схеме:

l l n n fi l l 2N 2N n n U q j j U q j f i Uqj, если U q j 0;

2N 2N fi j fi hj l l n n fi l l 2N 2N n n иначе U q j j U q j f i Uqj (1.57), 2N 2N fi j fi hj f (, u1, u2 )T, для невязок давления по формуле:

l l n n p pi l l 2N 2N n n ai j ai j ai, если U q j 0;

2N 2N i pi pi hj l l n n p p i l l 2N 2N n n иначе ai j ai j ai, (1.58) 2N 2N i pi pi hj j 1,2, l 1, 2 N, i 1, N j.

При k 2 :

f f j fi ( fi fi ) i 1 i 1, j j 2 2h j (1.59) f f 1 j fi ( j fi j fi ) i 1 i 1, j 1, 2.

2 2h j При повышении порядка аппроксимации до второго разностный оператор не зависит от знака U q j, определенной согласно (1.13), в узле и переходит к известной симметричной аппроксимации дифференциальных операторов со вторым порядком. Как уже упоминалось ранее, среди схем второго порядка аппроксимации нет монотонных схем, поэтому использовать схему (1.59) в чистом виде невозможно. Введем сглаживающий оператор в (1.59) с целью уменьшения нефизических осцилляций в решении. В п. 1.5.1 Главы I предложен разностный оператор второго порядка точности со сглаживающими членами для аппроксимации дифференциальных операторов в дивергентном виде. Т.к. для решения уравнений Навье-Стокса введено преобразование координат, то вычислительный эксперимент показал, что в таком случае применять разностный оператор в виде (1.22) невозможно. Для аппроксимации дифференциальных операторов в дивергентной форме использовалась формула:

Wl 1 Wl 1 W 2Wl Wl a sign(U q j ) 2 l Wl, (1.60) 2h 2h Wl 1 2Wl Wl, где | Wl 1 Wl | | Wl Wl 1 | 0, если | W W | | W W | 0.

l 1 l l l Оператор, определенный в (1.61), аппроксимирует дифференциальный оператор x со вторым порядком. В частном случае при 2 1 оператор j (1.60) имеет первый порядок аппроксимации, т. е. схема переходит в схему с направленными разностями. Введение замены переменных для системы уравнений Навье-Стокса приводит к тому, что в разностных операторах (1.55) – (1.60) направление аппроксимации дифференциальных операторов по противопотоковой схеме и знак монотонизирующего члена определяется знаком компоненты U j. Это условие следует из формы записи системы уравнений Навье-Стокса в дивергентной форме (1.12).

1.7 Краевые условия на дробных шагах При решении задач методом приближенной факторизации возникают вопросы о постановке краевых условий на дробных шагах. Получаемые на дробных шагах системы линейных алгебраических уравнений (СЛАУ) являются результатом аппроксимации последовательности промежуточных задач, которые декомпозирует исходную задачу на отдельные физические процессы. Каждая такая задача требует постановки краевых условий, дополнительно проблема заключается в том, что за счет введенного расщепления один целый шаг по времени разбивается на четыре дробных шага (для двумерного случая), а сами уравнения разрешаются не относительно искомых переменных, а относительно невязок решения, для которых необходимо поставить краевые условия. Для стационарных краевых условий первого рода проблема решается следующим образом. Пусть на некоторой границе выполняется условие u, тогда если решение не изменяется во времени, т.е. выполняется условие u n1 u n, переход на следующий временной слой формально осуществляется по формуле u n1 u n un1, отсюда следует, что un1 0. Для того, чтобы выполнялось условие un1 необходимо, чтобы на всех дробных шагах выполнялось требование unl 0, которое и является краевым условием.

Аналогично способом можно получить краевое условие второго рода.

Вновь примем, что на некоторой границе выполняется условие:

u, xi u n u n. Отсюда следует, что тогда полагаем, что xi xi u n1 u n (u n1 u n ) 0, 0, xi xi xi учитывая формулу перехода на новый временной слой, получаем краевое условие для невязок:

n 0. (1.62) xi Для того чтобы краевое условие (1.62) выполнялось на целом шаге, необходимо nl 0, последнее чтобы на всех дробных шагах выполнялось условие xi выражение и есть краевое условие второго рода для невязок.

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

Глава II. Исследование свойств разностной схемы В качестве тестовых задач для исследования свойств численных алгоритмов были выбраны задача о распаде произвольного разрыва и задача о квазиодномерном течении газа в канале. Эти задачи были выбраны для тестирования в силу следующих причин:

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

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

2.1 Задача о распаде произвольного разрыва Труба постоянного сечения заполнена покоящимся газом и перегорожена заслонкой, по обе стороны которой заданы различные параметры – давление, скорость и плотность. В момент времени t 0 заслонка мгновенно убирается.

Требуется описать изменение во времени параметров газа. Полная постановка задачи о распаде произвольного разрыва дана, например, в [45]. Решение задачи описывается одномерной нестационарной системой уравнений Эйлера, имеющей точное решение [см. 45].

Для численного эксперимента задавались параметры: число узлов расчетной сетки составило N 401 (hx 0.025), длина расчетной области l единиц, решение приведено на момент времени t 2.5. На нулевом дробном шаге производные аппроксимированы разностными операторами как первого, так и второго порядка точности, на дробных шагах производные аппроксимированы разностными операторами с первым порядком точности по противопотоковой схеме.

На рис. 1 и 2 приведены распределения скорости и плотности газа в канале, полученные со вторым порядком точности по схеме с оптимальным расщеплением операторов (1.30). Треугольниками дано численное решение, точное решение – сплошной линией. Как видно из рисунков, численное решение достаточно хорошо совпадает с точным решением. Количественные характеристики максимальной относительной погрешности в процентах даны в таблице 1, наибольшие отклонения численного решения от точного наблюдались в области ударной волны и контактного разрыва.

Рис. 1. Распределение скорости газа Рис. 2. Распределение плотности газа Таблица 1 – Относительная погрешность вычислений. Задача о распаде произвольного разрыва. Второй порядок точности числ точн uчисл uточн pчисл pточн max,% max,% max,% точн uточн pточн 0,35 0,054 0, На рис. 3 изображен сравнительный график распределения скорости газа для расчетов, выполненных с первым и вторым порядком точности (прямоугольники – первый порядок точности, кружочки - второй). Из рис. видно, что применение разностного оператора второго порядка точности приводит, как и ожидалось, к более точному численному решению задачи, а также к уменьшению размазывания решения, которое для численного решения со вторым порядком на ударной волне составляет 4 узла, тогда как для решения, полученного с первым порядком, – 7 узлов. Кроме того, использование специального оператора сглаживания позволяет в полученном численном решении исключить нефизические осцилляции, характерные для схем второго порядка точности и выше.

Рис. 3. Сравнительный график распределения скорости в канале для расчетов, выполненных с первым и вторым порядком точности.

2.2 Задача о квазиодномерном течении газа в канале Пусть имеется плоский канал (рис. 4), заданный по известной формуле A( x), площадь которого сначала монотонно убывает, а потом монотонно возрастает. Будем задавать течение на входе дозвуковым. Тогда течение в канале будет ускоряться и в критическом сечении x0 становится звуковым, т.е.

u c, где c - скорость звука, а за ним сверхзвуковым. Полагая критическим сечением значение x0, в котором достигается самое малое сечение канала, можно найти точное решение задачи, обеспечивающее непрерывный переход течения от дозвукового к сверхзвуковому. Такой канал носит название сопла Лаваля. Если в задаче на входе задавать дозвуковое течение, а на выходе произвольное течение, то решение уравнений квазиодномерной задачи так же существует, но оно будет разрывным [45]. Задача описывается системой квазиодномерных уравнений Эйлера [45]:

( A( x) u ) 0, x A( x) ( A( x) u 2 A( x) p ) p, x x ( A( x)u ( E p )) 0.

x Рис. 4 Сопло Лаваля Решение отыскивалось методом установления, условием выхода алгоритма на стационарный режим являлось выполнение следующего соотношения:

1 | ( u ) n1 ( u ) n |, 104. (2.1) ( u) n Для численного решения задачи были заданы параметры: шаг по пространству hx 0.0025 (количество узлов N x 401), длина расчетной области l 1.0, шаг по времени выбирался из условий наискорейшего установления и устойчивости схемы, параметр 1.0.

На рис. 5–6 приведены распределения плотности и скорости газа, полученные со вторым порядком точности (даны треугольниками). Из рисунков видно, что численное решение достаточно хорошо повторяет точное (дано сплошной линией). На рис. 7 изображен сравнительный график распределений скорости газа в канале, полученных с первым и вторым порядком точности (прямоугольники – первый порядок точности, кружочки - второй). Из рисунка видно, что использование разностного оператора второго порядка точности позволило существенно уменьшить зону размазывания решения на ударной волне – с 14 до 3 узлов.

Рис. 5. Распределение плотности газа Рис. 6. Распределение скорости газа Рис. 7. Сравнительный график распределения скорости газа в канале Относительная погрешность вычислений для численных решений со вторым порядком точности дана в таблице 2. Введение разностного оператора оказало незначительное влияние на скорость сходимости схемы к стационарному решению, так для установления решения со вторым порядком требуется на 70 итераций больше, чем при установлении решения с первым порядком точности.

Таблица 2 – Относительная погрешность вычислений. Задача о квазиодномерном течении газа в канале. Второй порядок точности числ точн uчисл uточн pчисл pточн Количество Число,%,%,% точн uточн pточн итераций до Куранта установления 0,3 0,2 0,1 400 7, Наилучшая сходимость к стационарному решению достигалась при шагах ( u c) p сетки, соответствующим числам Куранта 1 K 8 ( K, где c h – скорость звука). Оптимальный шаг по времени для заданных условий численного эксперимента составил 1.04hx 0.0026.

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


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

На рисунках 8–9 представлены численные решения для плотности, скорости и давления, полученные по нефакторизованной схеме (1.23).

Рис. 8. Распределение плотности газа Рис. 9. Распределение скорости газа В таблице 3 приведена оценка относительной погрешности вычислений в процентах для плотности, скорости и давления газа, число итераций, необходимых схеме для установления по заданному критерию (2.1), а также оптимальное число Куранта, с которым схема позволяет проводить расчет.

Таблица 3 – Относительная погрешность вычислений. Нефакторизованная схема числ точн uчисл uточн pчисл pточн Количество Число,%,%,% точн uточн pточн итераций до Куранта установления 0,1 0,1 0,25 101 280, На рис. 10 изображен сравнительный график распределений скорости газа в канале, полученных с первым и вторым порядком точности (прямоугольники – первый порядок точности, кружочки - второй). Для данного расчета ширина зоны перехода в области скачка уплотнения составила 3 узла.

Рис. 10. Распределение скорости газа, полученное с первым и вторым порядком точности. Область скачка уплотнения Шаг по времени был выбран из условия наискорейшей сходимости и для заданных параметров вычислительного эксперимента составил 40hx 0.1.

2.3 Оценки области оптимальной скорости сходимости На рис. 11 и 12 дано изображение зависимости между шагом по времени и количеством итераций, понадобившимся рассматриваемым схемам (дополнительно дан график для схемы с модифицированным оптимальным расщеплением) для установления. На рис. 11 графики даны для схемы с оптимальным расщеплением (график дан сплошной линией) и для модифицированной схемы с оптимальным расщеплением (изображено штрихпунктирной линией), на рис. 12 изображен график для нефакторизованной схемы.

Рис. 11. Область оптимальной Рис. 12. Область оптимальной сходимости. Схемы с расщеплением сходимости. Нефакторизованная схема Из рисунков видно, что схемы с расщеплением операторов имеют примерно одинаковый оптимальный параметр, область оптимальной сходимости у схемы с оптимальным расщеплением операторов немного шире, чем у модифицированной схемы. Для нефакторизованной схемы область оптимального параметра значительно смещена относительно первых двух схем и находится в окрестности точки 0.1 (число Куранта см. в таб. 2), где достигается минимальное количество итераций. Так же из рис. 12 видно, что график оптимальной сходимости для нефакторизованной схемы имеет немонотонный характер, что возможно связано с нелинейностью решаемой системы дифференциальных уравнений.

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

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

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

1.5.1 Главы I диссипативный член, возникающий вследствие введения расщепления, имеет вид:

f n1 f n D B B 2 2 B1n B2 n (ft O( )).

2 2 1n 2n (2.2) h h h h Если теперь обратиться к записи системы уравнений в недивергентной векторной форме, то запись (2.2) перейдет к виду 2 2 B1n B2 n Bn f n, которая h h h содержит в себе третьи производные по пространству. Тогда исходная система уравнений газовой динамики после введенного расщепления и факторизации изменит свой вид и перейдет к следующему:

f f 3f B C 3 0 (2.3) t x x C 2 2 B1n B2 n Bn. Уравнение (2.3) представляет собой где оператор h h h уравнение Кортевега – де Фриза, в котором за наличие низкоамплитудных знакопеременных колебаний в решении отвечает член с третьей производной по пространству. Таким образом, вводя расщепление исходного оператора, мы исходную систему уравнений заменяем неким аналогом, в котором присутствуют члены, вызывающие наличие колебаний в решении. Проблема в данном случае решается простым осреднением решения за один полный период колебания. Как было показано на рисунках и приведенных количественных оценках точности полученного численного решения такой подход вполне оправдан и позволяет проводить расчеты с достаточной точностью.

В заключение параграфа приведем еще один интересный факт.

Модифицированную разностную схему с оптимальным расщеплением с расщеплением операторов в форме представим в виде:

f n1 f n (A) (I B )(I B ) Whn.

n 1n 2n (2.4) Введем обозначение:

(A) n (I B1n )(I B2 n ), (2.5) и запишем схему в виде универсального алгоритма:

f n1 f n Whn.

Для решения стационарных задач методом установления, т.е. когда решение задачи находится как предельное решение нестационарной задачи со стационарными краевыми условиями, выбор релаксационного оператора может осуществляться из условий наиболее быстрой сходимости [29]. Для разностной схемы вид релаксационного оператора в форме (2.5) обусловлен выбором расщепления исходных операторов и переменных f f (, u, p ) (, u, p )T, относительно которых разрешены уравнения на дробных шагах. Заметим, что схема с оптимальным расщеплением допускает решение стационарных и нестационарных задач. С другой стороны ничто не мешает изменить вид релаксационного оператора (2.5) и подобрать его из условий наискорейшей сходимости. Оставаясь в рамках расщепления (1.48) – (1.49) изменим вид матрицы A на следующий:

1 0 0 1 u A 0, A -1 u 0.

( 1)u u2 u 2 Разностная схема с релаксационным оператором, в который входит измененная матрица A, позволяет решать только стационарную задачу о квазиодномерном течении газа в канале. Изменение релаксационного оператора позволило повысить устойчивость схемы и решать задачу с шагом по времени до 3.6hx (число Куранта K 25.44 ), тогда как модифицированная схема с оптимальным расщеплением позволяет решать эту задачу с шагом по времени 1.04hx (число Куранта K 7.36 ). В таблице 4 приведены данные об относительной ошибке решения, количестве итераций, понадобившемся схеме для установления по условию, и число Куранта, с которым схема позволила проводить расчет. Данные даны для расчетов с первым порядком аппроксимации на сетке в 401 узел.

Таблица 4 – Относительная погрешность вычислений. Модифицированный релаксационный оператор числ точн uчисл uточн pчисл pточн Количество Число,%,%,% точн uточн pточн итераций до Куранта установления 0,9 0,5 1,1 198 25, Из таблицы видно, что по сравнению с разностной схемой, схема с измененным релаксационным оператором является более устойчивой и позволяет проводить расчет с числом Куранта равным 25,44, устанавливается за значительно меньшее число итераций (198 итераций против 364 итераций).

Отметим, что схема с модифицированным релаксационным оператором не дает возможности решать нестационарные задачи, так как оператор уже не соответствует зависимости между дивергентной и недивергентной формами записи уравнений (см. Главу I).

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

2.4 Тестирование алгоритма на двумерных задачах Рассмотрим результаты тестирования алгоритма на расчете двумерных стационарных течений. Отметим, что все стационарные решения получены методом установления по схеме (1.40) с расщеплением матричных операторов по формулам (1.37) – (1.38). Преобразования координат задавалось в виде:

xa y ( x) y1 ( x) q1, q2, (20) ba y2 ( x) y1 ( x) где y1 ( x), y2 ( x) – кусочно-линейные функции. Рассмотрим результаты расчетов, которые были проведены для задачи о течении вязкого сверхзвукового теплопроводного газа в канале типа «воздухозаборник».

Расчетная область представляет собой плоский канал (рис. 13), нижняя стенка которого начинается в точке 0.0, а верхняя снесена по оси абсцисс относительно нижней стенки на 0.5 единицы. Высота канала составляет единицы, длина – 6 единиц, перед каналом выделена область с невозмущенным потоком длиной 0.5 единицы. На твердых стенках канала поставлены условия прилипания газа к стенкам u v 0, а так же условия теплоизоляции стенок канала T 0, n на выходе из канала поставлены мягкие условия f 0.

x На входе в канал задается невозмущенный сверхзвуковой поток с 1.4, u 1.0, v 0, 1.0, p M 2, Re 1000, параметрами, M Pr 0.72, 0.76 под нулевым углом атаки. Расчет выполнен на равномерной стеке 521 161 узел ( hx hy 0.0125 ), итерационный шаг min(hx, hy ), параметр схемы 1.0, установление получено за 2500 итераций.

Рис. 13 Расчетная область На рис. 14 и 15 представлены распределения полей скорости и температуры газа соответственно. Из рис. 14 хорошо видно, что на стенках канала развивается пограничный слой, а в центре канала остается сверхзвуковое ядро течения. Из рис. 15 видно, что наиболее разогретый газ находится около стенок канала, при продвижении газа вдоль границ канала происходит дальнейший прогрев зон, расположенных около пограничных слоев.

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

Рис. 16 Изолинии плотности в канале.

Рис. 17 Изолинии давления в канале По данным, приведенным на рисунках, можно оценить угол падения ударной волны, он составляет 36o11'. По известной зависимости между числом Маха набегающего потока и углом раствора клина для газодинамических невязких течений можно определить теоретический угол отклонения от плоскости для возникающей ударной волны. Для заданных параметров течения он равен 30o. Учитывая, что мы имеем дело с течением вязкого газа и на пластине в результате торможения образовывается пограничный слой, который играет роль клина, угол косого скачка уплотнения изменяется. Оценив угол раствора клина, который равен 5o, угол скачка должен составлять 35o, что близко к полученному численному результату.

С целью исследования точности получаемого решения и оценки его сходимости при дроблении сетки были проведены расчеты на последовательности сгущающихся сеток. Для расчетов была выбрана начальная сетка 521 81, которая последовательно дробилась в 2 и в 4 раза. На рис. 18 и 19 приведены распределения величин тангенциальной и нормальной компонент скорости в продольном сечении вдоль координаты Y 0.25, на рисунках 20, – распределения плотности и давлениях в сечении Y 0.125. Решения, полученные на сетках 81, 161 и 321 узел, отмечены на графиках цифрами 1, 2 и 3 соответственно.

Рис. 18 Распределение продольной компоненты скорости Рис. 19 Распределение нормальной компоненты скорости На всех рисунках можно видеть, что решение 1, полученное на самой грубой сетке смещено относительно решений 2 и 3, полученных на сетках в и 321 узел. Это говорит о том, что решение 1 имеет меньший угол наклона ударной волны, чем решения 2 и 3. Кроме того, решение 3 имеет немного больший угол наклона ударной волны и составляет 36o 21'.

Рис. 20 Распределение плотности в продольном сечении Рис. 21 Распределение давления в продольном сечении Решения 2 и 3 достаточно хорошо совпадают до отметки x 4. (максимальная разница решений составляет 2,5%), а далее происходит небольшое расхождение в количественных характеристиках потока. На рис. представлены распределения давления в поперечном сечении канала X 3.0.

Из приведенных рисунков хорошо видно, что характер течений 2 и 3 совпадает, максимальная разница в количественных характеристиках составляет 8%.

Рис. 22 Распределение давления в поперечном сечении На рис. 23 и 24 приведены распределения числа Маха, полученные на последовательности сеток, в поперечных сечениях X 1.125 и X 2.25.

Распределения числа Маха, приведенные на рис. 23, соответствуют течению, частично прошедшему через первый скачок уплотнения, на рис. 24 – течению, прошедшему через первый скачок уплотнения и частично через второй скачок.

Из рис. 24 видно, что решения 2 и 3 имеют другой угол наклона ударной волны, хотя количественные характеристики течения до и после скачков уплотнения совпадают у всех трех решений. Расхождение в угле наклона скачка уплотнения для решений 2 и 3 составляет не более 0.8%. На рис. наблюдается небольшое расхождение численных решений в окрестности точки Y 1.5, однако максимальное расхождение 2 и 3 решений составляет не более 1%.

Рис. 23 Распределение числа Маха в поперечном сечении Рис. 24 Распределение числа Маха в поперечном сечении Из приведенных выше результатов можно сделать следующие выводы:

1) предложенная схема второго порядка точности позволяет проводить расчеты сложных течений без возникновения у них нефизических осцилляций;

2) предложенная разностная схема устойчива в широком диапазоне итерационных параметров и позволила проводить расчет течения с параметром min(hx, hy ), что соответствует числам Куранта в диапазоне 1.5 K 2.5 ;

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

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

Глава III. Моделирование течений газа в канале воздухозаборника со вдувом газа с части поверхности Данная глава посвящена численному исследованию стационарного плоского течения газа в канале воздухозаборника. В качестве модели для описания течения используется полная система уравнений Навье – Стокса вязкого сжимаемого теплопроводного газа. Для исследуемого диапазона параметров потока течение предполагалось ламинарным. Основное внимание уделено изучению структуры течения, образующейся внутри канала в результате вдува газа, а так же исследованию влияния параметров газа в набегающем потоке на характер течения.

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

Решение полной системы уравнений Навье-Стокса вязкого сжимаемого теплопроводного газа, принятой для описания течения, отыскивается в расчетной области, представленной на рис. 25. Расчетная область R {0.5 x1 6, 0 x2 2}, внутри представляет собой прямоугольник которого находится канал. Нижняя стенка канала отнесена от границы расчетной области на расстояние 0.5 и начинается в точке x1 0.0, а верхняя стенка смещена относительно нижней стенки вправо и начинается в точке x1 0.5.

Рис. 25. Расчетная область Стационарное решение задачи находилось в приближении уравнений Навье-Стокса в преобразованных координатах (1.11) по разностной схеме с оптимальным расщеплением операторов (1.40) методом установления.

Преобразование координат задавалось в виде x1 a x1 ( x1 ) q1, q2 2 2, ba x2 ( x1 ) x1 ( x1 ) здесь x1 ( x1 ), x2 ( x1 ) – уравнения нижней и верхней стенок канала;

a, b – левая и правая границы расчетной области.

Система уравнений Навье-Стокса дополняется следующими граничными условиями: на границе области при x1 0.5 задавался невозмущенный u0 1, v 0, 0 1, p сверхзвуковой поток с параметрами, M параллельный оси x1, на линиях x2 0, 0.5 x1 0.0 и x2 2, 0.5 x1 0.5 – условия симметрии течения:

v p u 0, n n n где n - вектор нормали к границе. На твердых стенках канала поставлены условия прилипания газа к стенкам:

u v 0, а так же условия теплоизоляции стенок канала T 0.

n Так как при решении задачи заранее не известно, как ведет себя поток на выходе из канала, то обычно на такой границе ставят так называемые «мягкие условия»:

f 0, x которые предполагают экстраполяцию значений искомых функций на границу. На нижней стенке канала задавался источник подачи газа, начинающийся в точке x 1.0 шириной 0.05 единицы. Скорость вдува газа M1 1.0, число Маха набегающего потока полагалась звуковой, т.е.

варьировалось в пределах M 0 1.5 2.5, число Рейнольдса – в пределах Re 103 104, значение адиабаты, число Прандтля и показатель степени были выбраны следующими 1.4, Pr 0.72, 0.76.

Параметры разностной схемы и итерационный шаг для данного расчета были выбраны следующие: 1.0, 0.4min(hx, hy ). Выбор значений параметров основан на следующих размышлениях: выбирая релаксационный параметр, который определяет степень «неявности»

схемы, равным единице, мы полагаем, что схема является абсолютно неявной, т.е. f f n1 f n f n1. Выбор параметра является отдельной оптимизационной задачей. Как было показано в Главе II п. 2.3, даже в одномерном случае у всех рассматриваемых схем, в т.ч. у нефакторизованной схемы, есть область оптимальной скорости сходимости, несмотря на их абсолютную устойчивость. Однако абсолютная устойчивость схем была показана для линеаризованной схемы, где не учитывалось влияние нелинейности исходной системы уравнений. Путем методических расчетов, начиная с параметра min(hx, hy ), была определена область оптимального значения параметра при заданном параметре, при котором достигается минимальное количество итераций, за которое решение устанавливается.



Pages:   || 2 |
 

Похожие работы:





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.