WWW.DISS.SELUK.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА
(Авторефераты, диссертации, методички, учебные программы, монографии)

 

Pages:     || 2 |

«ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ГАЗА НА ОСНОВЕ МОДИФИЦИРОВАННОГО МЕТОДА РАСЩЕПЛЕНИЯ ...»

-- [ Страница 1 ] --

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

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

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

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

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

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

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

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

Научный руководитель д.ф.-м.н., проф. В.М. Ковеня Новосибирск,

ОГЛАВЛЕНИЕ

Введение

Глава 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, где приведены результаты тестирования алгоритма на сходимость численного решения на последовательности сеток.

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

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

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

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

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

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

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

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

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

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

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

вычислительных технологий СО РАН «Информационно-вычислительные технологии в задачах поддержки принятия решений».

Работа выполнена при частичной финансовой поддержке РФФИ (гранты № 05-01-00146, 08-01-00264).

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

1.1 Исходная система уравнений Эйлера и Навье-Стокса 1.1.1 Системы уравнений в одномерном случае теплопроводного газа в дивергентной форме в виде (см., например, [47]) Для замыкания системы уравнений (1.1) зададим уравнения состояния и зависимости вязкости и теплопроводности, от температуры (T ), (T ).

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

Наряду с консервативной формой представим исходные уравнения (1.1) в недивергентной и предельно-дивергентной формах Отметим эквивалентность записи уравнений Навье-Стокса в дивергентном (1.1) и в недивергентном (1.3) виде. Выбор вектора f задает форму матричных операторов B, что дает возможность рассматривать различные классы разностных схем. Для систем уравнений Эйлера и Навье-Стокса искомыми функциями обычно выбирают плотность, скорость, а также давление или температуру, что определяется заданием краевых условий. Однако вид матричных операторов B для различных компонент вектора f существенно различается [42]. Выберем в качестве искомых функций плотность, скорость и давление газа, тогда для уравнения состояния совершенного газа p ( 1) e матрица B принимает вид дивергентном виде, что позволяет упростить вид матричного оператора B.

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

где вектор искомых функций, газодинамические и вязкие потоки по направлениям x1 и x2 соответственно. Компоненты тензора вязких напряжений и тепловые потоки задаются соотношениями:

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

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

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

функций. Выберем вектор искомых функций в виде:

Тогда система уравнений примет следующий вид:

где здесь ij - символ Кронекера.

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

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

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

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

При записи уравнений введены обозначения:

J - якобиан преобразования координат ( x1, x2 ) в (q1, q2 ), задаваемый в виде:

Переходя к недивергентной форме, представим систему уравнений НавьеСтокса в новых переменных в виде:

Здесь обозначения x1x1, x1x2, x2 x2, Tq1, Tq2 приведены в (1.14).

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

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

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

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

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

уравнений (1.5) может быть представлена в предельно-дивергентной форме:

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

Из (1.10) следует, что:

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

Отметим, что для введенной замены переменных (1.11) матрица A изменит свой вид на A J A.

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

компоненты векторов E и F не изменят своего вида по сравнению с (1.13), изменит свой вид только коэффициент при градиенте температуры:

число Маха, c p - удельная теплоемкость газа при постоянном давлении, c - скорость звука. Уравнения (1.16) примут вид:

Что касается коэффициента при градиенте температуры, то он будет иметь вид:

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

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

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

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

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

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

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

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

где с учетом соотношения (1.3) правая часть схемы (1.23) аппроксимирована в уравнений (1.4) с порядком O( h 2 ), а при установлении установления (при выполнении условия дивергентном виде (1.1). с порядком O(h 2 ). В линейном приближении схема аппроксимирующий оператор B, перепишется в следующем виде:

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

Разностная схема (1.23) в дробных шагах (1.25)-(1.27) реализуется векторными прогонками на каждом дробном шаге и требует обращения матриц размерностью 3 3.

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

где разностные операторы имеют вид:

Для численного решения системы уравнений Навье–Стокса рассмотрим схему приближенной факторизации:

или эквивалентную ей схему в дробных шагах:

которая аппроксимирует уравнения (1.1) с порядком O( h 2 ) при всех, при установлении схема консервативна, т.е. Whn 0 в силу того, что A 1 0.

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

(1.32) и получим уравнение, разрешенное относительно p 2 :

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

Введение расщепления (1.28) приводит к появлению в разностной схеме (1.30) дополнительных членов вида D 2 2 B1n Bh n – диссипативных членов второго порядка малости. Поэтому расщепление матричных операторов Bh необходимо выбирать таким образом, чтобы построенная разностная схема (1.30) сохраняла свойство безусловной устойчивости и имела минимальную диссипацию – минимальное число дополнительных членов. В работе [42] было расщепления по физическим процессам в форме данное расщепление позволило свести реализацию разностной схемы вида (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) определим во всех узлах расчетной сетки, совпадающие с исходными функциями в узлах сетки. Дифференциальные операторы будем аппроксимировать разностным оператором j в узлах сетки с порядком пространственных производных, N – размерность задачи.

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

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

Перейдем от непрерывной модели (1.33) к конечно-разностному аналогу. Для разностными аналогами. Тогда разностная схема с весами:

аппроксимирует систему уравнений (1.33) с порядком O(h m ). Перепишем разностную схему в каноническом виде:

где с учетом соотношения (1.17) правая часть схемы (1.35) аппроксимирована в предельно-дивергентной форме.

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

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

где элементы bii содержат повторные производные по пространственному направлению из компонентов тензора вязких напряжений и градиента температуры.

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

Тогда разностная схема примет вид:

Она аппроксимирует исходную систему уравнений с порядком O(h m ), а при установлении – с порядком O(h m ), в линейном приближении при 0. является абсолютно устойчивой для двумерного случая [34].

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

Очевидно, что разностная схема (1.41) – (1.43) реализуется скалярными прогонками на 2 и 4 дробных шагах, т.к. матрицы B j содержат ненулевые элементы лишь на главной диагонали. Покажем, что на 1 и 3 дробных шагах схема реализуется также скалярными прогонками. Запишем систему уравнений, решаемую на первом дробном шаге:

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

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

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

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

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

а оператор B j изменит свой вид на следующий:

элементы bii содержат повторные производные по пространственному направлению из компонентов тензора вязких напряжений и градиента температуры.

(I Bnj ) по формуле, учитывая расщепление (1.47):

Тогда разностная схема примет вид:

Она аппроксимирует исходную систему уравнений с порядком O(h m ), а при установлении – с порядком O(h m ), в линейном приближении при 0. является абсолютно устойчивой для для двумерного случая [34].

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

В связи с изменением расщепления, теперь на первом и третьем дробном шаге счет ведется относительно невязок скоростей u и v, а на втором и четвертом дробных шагах счет ведется относительно невязок u и v.

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

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

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

затем и отсюда необходимые значения невязок для второго дробного шага:

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

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

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

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

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

где U q j и J определены в (1.13) и (1.15) соответственно. На дробных шагах аппроксимация первых производных выполнена по противопотоковой схеме:

для невязок давления по формуле:

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

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

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

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

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

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

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

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

Глава 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 – Относительная погрешность вычислений. Задача о распаде произвольного разрыва. Второй порядок точности На рис. 3 изображен сравнительный график распределения скорости газа для расчетов, выполненных с первым и вторым порядком точности (прямоугольники – первый порядок точности, кружочки - второй). Из рис. видно, что применение разностного оператора второго порядка точности приводит, как и ожидалось, к более точному численному решению задачи, а также к уменьшению размазывания решения, которое для численного решения со вторым порядком на ударной волне составляет 4 узла, тогда как для решения, полученного с первым порядком, – 7 узлов. Кроме того, использование специального оператора сглаживания позволяет в полученном численном решении исключить нефизические осцилляции, характерные для схем второго порядка точности и выше.

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

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

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

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

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

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

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

квазиодномерном течении газа в канале. Второй порядок точности Наилучшая сходимость к стационарному решению достигалась при шагах – скорость звука). Оптимальный шаг по времени для заданных условий численного эксперимента составил 1.04hx 0.0026.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

расщеплением операторов в форме представим в виде:

Введем обозначение:

и запишем схему в виде универсального алгоритма:

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

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

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

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

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

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

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

Расчетная область представляет собой плоский канал (рис. 13), нижняя стенка которого начинается в точке 0.0, а верхняя снесена по оси абсцисс относительно нижней стенки на 0.5 единицы. Высота канала составляет единицы, длина – 6 единиц, перед каналом выделена область с невозмущенным потоком длиной 0.5 единицы. На твердых стенках канала поставлены условия прилипания газа к стенкам а так же условия теплоизоляции стенок канала на выходе из канала поставлены мягкие условия На входе в канал задается невозмущенный сверхзвуковой поток с Pr 0.72, 0.76 под нулевым углом атаки. Расчет выполнен на равномерной стеке 521 161 узел ( hx hy 0.0125 ), итерационный шаг min(hx, hy ), параметр схемы 1.0, установление получено за 2500 итераций.

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

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

По данным, приведенным на рисунках, можно оценить угол падения ударной волны, он составляет 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. Расчетная область представляет собой прямоугольник которого находится канал. Нижняя стенка канала отнесена от границы расчетной области на расстояние 0.5 и начинается в точке x1 0.0, а верхняя стенка смещена относительно нижней стенки вправо и начинается в точке x1 0.5.

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

Преобразование координат задавалось в виде здесь x1 ( x1 ), x2 ( x1 ) – уравнения нижней и верхней стенок канала; a, b – левая и правая границы расчетной области.

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

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

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

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

варьировалось в пределах 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 ), была определена область оптимального значения параметра при заданном параметре, при котором достигается минимальное количество итераций, за которое решение устанавливается.

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

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

Аналогичные явления наблюдались и другими авторами (например, [40]).

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

Расчеты были выполнены на сетке, содержащей 521 161 узел. Выбор такой сетки был обусловлен проведением расчетов как на более грубых, так и на более подробных стеках. Расчеты на последовательности сеток по аналогии с Главой II п. 2.3 показали, что на выбранной сетке при сгущении узлов в 2 раза, решение отличается не более чем на 1.5% для всех искомых функций.

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

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

В рамках данного исследования использовались оба варианта начального распределения. Установление к стационарному решению достигалось за 8500 итераций при использовании в качестве начального распределения не согласованного распределения параметров потока (при x1 0.5, 0 x задавались параметры набегающего потока, а в расчетной области задавалась нулевая скорость и постоянные значения плотности и температуры). При задании в качестве начального приближения решения задачи при других параметрах (например, при другом числе Рейнольдса или числе Маха) количество итераций сокращалось в 1,5 – 2 раза до 4500 5500 итераций.

3.3 Расчеты с различными числами Маха Перейдем к рассмотрению результатов численного моделирования. В первой серии расчетов изучено влияние изменения числа Маха набегающего потока на картину течения при фиксированном числе Рейнольдса Re 104.

Дадим описание типовой картины течения газа в канале. На рис. приведены распределение поля скорости и изолиний давления газа для числа Маха набегающего потока M 2.0. При попадании невозмущенного потока газа в канал на его стенках образуются головные скачки уплотнения. При продвижении потока газа внутрь канала на его стенках развивается пограничный слой. Наличие источника вдува на нижней стенке канала существенно влияет на общую картину течения. В отличие от результатов расчетов, приведенных в предыдущей главе, наличие источника вдува приводит к тому, что перед ним формируется область повышенного давления и как следствие формируется сильный скачок уплотнения. Этот возникший скачок уплотнения распространяется от нижней к верхней стенке и взаимодействует с пограничным слоем. Отражаясь от пограничного слоя, скачок уплотнения (ударная волна) вызывает образование системы из волны разрежения и скачков уплотнения. В зоне взаимодействия пограничного слоя и скачка уплотнения в зависимости от интенсивности последнего может возникать отрыв пограничного слоя от стенки канала. Влияние интенсивности набегающего потока на образование отрыва пограничного слоя на верхней стенке будет продемонстрировано ниже.

Рис. 26. Распределение поля скорости (а) и изолинии давления газа (б) в Перед источником вдува на нижней стенке формируется область сложного течения. На отдалении от источника вдува газа перед ним образуется зона отрыва пограничного слоя, которая вызвана следующими факторами: интенсивно выдуваемый газ вызывает образование сильного скачка уплотнения, этот скачок взаимодействуя с пограничным слоем на нижней стенке и вызывает его отрыв. Непосредственно перед источником вдува газа образуется зона вихревого движения газа. Направление движения газа в этих двух зонах вихревого течения различается. В области перед источником вдува направление течения против часовой стрелки, направление движения газа в области отрыва пограничного слоя по часовой стрелке. Эта особенность движения газа подтверждается так же расчетами других авторов [6] и результатами экспериментальных данных [1,13,14]. За источником вдува образуется протяженная зона отрыва пограничного слоя. За зоной отрыва пограничного слоя возникает скачок уплотнения, который, пересекая канал, взаимодействует с верхней стенкой, отражаясь от нее в зависимости от параметров набегающего потока либо уходит за границы вычислительной области, либо успевает отразиться от пограничного слоя на нижней стенке канала.

Рис. 27 Теплеровская фотография течения газа и схематичное изображение течения газа около струйного препятствия на пластине. Данные На рис. 27 представлена теплеровская фотография обтекания струи и схематичное изображение течения, полученного в результате экспериментального исследования обтекания струйного препятствия на пластине [14]. На пластине на некотором отдалении от ее начала, так чтобы успел образоваться пограничный слой, в плоскости ее симметрии расположено круглое отверстие, через которое выдувается газ с числом Маха M 1.0. На пластину набегает невозмущенный сверхзвуковой поток с числом Маха M 3.0. На рисунке представлена теплеровская картина течения для p0 j P0 j / P 128, где P0 j - полное давление в струе, P статическое давление в набегающем потоке. Основное внимание авторами было уделено изучению структуры возникающего явления, а так же выявлению зависимости между структурой течения и изменением полного давления в струе. На основе полученных данных авторы работы [14] сделали вывод о характере течения и графически представили его в виде схемы, изображенной на рис. 27б. Из сравнения рис. 26 и 27 можно отметить, что наблюдается хорошее совпадение качественной картины течений. Однако, скачок уплотнения 2 на рис. 27а слабо виден на рис. 26, что объясняется тем, что для данных расчетов число Маха набегающего потока было выбрано равным M 0 2. Как будет показано ниже, при увеличении числа Маха набегающего перед источником вдува газа четка видна сложная -образная система скачков уплотнения.

Рис. 28. Области отрыва пограничного слоя на нижней (а) и верхней (б) На рис. 28 представлены области около источника вдува газа и область в районе взаимодействия головного скачка уплотнения и пограничного слоя на верхней стенке канала в увеличенном масштабе. Как хорошо видно из рис.

28, поток, проходя через скачки уплотнения, изменяет направление движения, что хорошо согласуется как с теоретическим представлениями о характере движения потока газа через косой скачок уплотнения, так и с экспериментальными данными [14].

Отдельно остановимся на применении схемы первого и второго порядка аппроксимации. На рис. 29 приведены сравнительные графики распределения продольной компоненты скорости и плотности газа в продольных и поперечных сечениях. Расчет выполнен для параметров Re 104, =1.4, Pr 0.72, =0.76, M 2.0, M 1 1.0, на сетке 521 161 узел.

Как хорошо видно из рисунков, в тех подобластях, где отсутствуют особенности течения (скачки уплотнения, отрыв пограничного слоя) количественные и качественные характеристики течения в приведенных сечениях совпадают. Однако, для сечения по линии x 2.5 из рис. 29а видно, что расчет с первым порядком размазывает структуру скачка уплотнения, при расчете со вторым порядком видно, что скачок в действительности имеет более сложную структуру. Таким образом, использование для расчета схемы первого порядка аппроксимации не позволяет детально описать структуру скачков уплотнения. Аналогичные выводы можно сделать и для рис. 29б. Для двух приведенных продольных сечений видно, что при расчете со вторым порядком точности удается четко выделить три скачка уплотнения (область 1.8 x 3.0 ), тогда как расчет с первым порядком аппроксимации позволяет увидеть только два скачка (тройная структура пиков плотности, график красного цвета и двойная структура пиков на графике зеленого цвета). Для сечения по координате Y 0.5 видно, что расчет с первым порядком точности в области 3.0 x 4. не отражает прохождение потока через скачок уплотнения, тогда как для расчета со вторым порядком аппроксимации это хорошо видно (голубой график – первый порядок точности, синий график - второй). В целом из приведенных рисунков можно сделать вывод о том, что использование разностных операторов второго порядка аппроксимации позволяет выделять структуры скачков уплотнения, которые не позволяет получать расчет с первым порядком аппроксимации, а так же уточнять количественные характеристики в областях прохождения потока через скачки уплотнения.

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

Рис. 29. Сравнительные графики продольной компоненты скорости (а) и плотности газа, выполненные с первым и вторым порядком аппроксимации, Перейдем к подробному рассмотрению результатов расчетов с варьированием числа Маха набегающего потока и анализу его влияния на характер течения в канале. На последующих рисунках представлены поля скорости и изолинии давления газа для двух предельных случаев рассматриваемого диапазона чисел Маха - 1.5 и 2.5. Из рисунков хорошо видно, что изменение числа Маха в набегающем потоке существенно влияет на картину течения в канале. Как следует из теории ударных волн наклон головного скачка уплотнения для случая M 1.5 (рис. 30, 32) существенно больше, чем для числа Маха 2.5 (рис. 31 и рис. 33). В первом случае головной скачок уплотнения, возникший около источника вдува, падает на верхнюю стенку и отражается от нее (точнее от пограничного слоя), затем отражается от нижней стенки канала. Во втором случае скачок уплотнения вызывает отрыв пограничного слоя, возникают два скачка: один инициирован отрывом пограничного слоя, а второй отражением от пограничного слоя. Далее они также взаимодействуют с потоком в канале и в силу того, что они имеют значительно меньший угол наклона, то для рассматриваемой расчетной области не успевают взаимодействовать с нижней стенкой.

От значений чисел Маха набегающего потока существенно зависят размеры областей отрыва пограничного слоя на нижней и верхней стенках канала. При малых числах Маха интенсивность головного скачка недостаточна для отрыва пограничного слоя у верхней стенки и отрыв возникает лишь при M 1.5. На рис. 34 и 35 представлены распределения коэффициентов поверхностного трения для нижней и верхней стенок канала это подтверждающие. Как и ожидалось, наименьшие размеры отрывной зоны в окрестности вдува наблюдаются для меньшего числа Маха набегающего потока M 1.5, наибольшие – для M 2.5 (рис. 34).

Рис. 34. Коэффициент трения на нижней стенке канала Рис. 35. Коэффициент трения на верхней стенке канала На рис. 36 и 37 представлены распределения полей скорости и области отрыва пограничного слоя около источника вдува для чисел Маха 1.5 и 2. соответственно. Из рисунков видно, что размеры отрывных зон за источником вдува газа значительно отличаются, что связано с интенсивностью потока газа после прохождения головных скачков. При оценке размеров отрывной зоны для верхней стенки канала можно сказать, что для числа Маха 1.5 интенсивности ударной волны не хватает для отрыва пограничного слоя (рис. 35, красная линия, помечена цифрой 1). Начиная со следующего исследуемого числа Маха M 1.75 и для всех последующих значений, наблюдается отрыв пограничного слоя. Наибольшая ширина отрывной зоны, как и для нижней стенки, наблюдается для максимального числа Маха M 2.5.

Рис. 38. Распределение чисел Маха в сечении Y=0. Рис. 39. Распределение чисел Маха в сечении X=1. На рис. 38 – 39 представлены распределения чисел Маха в продольном и поперечном сечении канала. Отметим, что поток газа после прохождения через головной ударный скачок на уровне Y=0.5 затормаживается до примерно одинаковой скорости независимо от скорости набегающего потока (отличие скоростей потоков составляет не более 6%). В поперечном сечении x 1.5 профили скорости подобны для различных чисел Маха.

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

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

Напомним, что приведенные результаты соответствуют числу Рейнольдса 104. Влияние чисел Re на течение исследовано ниже.

3.4 Расчеты с различными числами Рейнольдса В следующей серии расчетов было проведено моделирование течения газа в канале воздухозаборника со вдувом газа с части поверхности для различных чисел Рейнольдса набегающего потока при фиксированном числе M 2.0. Число Re последовательно задавалось равным 103, 5 103, 104.

Расчет выполнен для параметров =1.4, Pr 0.72, =0.76, M 0 2.0, M 1 1.0, на сетке 521 161 узел.

На рис. 42 – 43 представлены поля скоростей и линии равного уровня давления для всех исследуемых чисел Рейнольдса. Из рисунков видно, что с уменьшением числа Re толщина пограничного слоя на стенках возрастает и от носков нижней и верхней стенок возникает все более сильные скачки уплотнения. Одновременно с уменьшением числа Re происходит увеличение толщины скачков уплотнения, что связано с увеличением вязкости газа. В окрестности нижней стенки первый головной скачок и пограничный слой взаимодействуют с головным скачком, возникшим в области вдува. За отрывной зоной в области вдува возникает висячий скачок, который взаимодействует с отраженными скачками. Дальнейшее их взаимодействие приводит к сложной картине течения, которая, конечно, зависит от чисел Re.

Отрывные зоны у верхней стенки канала, возникшие после падения головного скачка, приводят к перестройке течения и возникновению вторичных отраженных скачков.

Рис. 43. Поле скоростей при различных Re Рис. 44. Изолинии давления при различных Re Видно (рис. 43в и 44в), что с уменьшением чисел Re толщина пограничного слоя и область с дозвуковыми скоростями возрастают и скачки уплотнения отражаются уже не от стенок канала, а от развитого пограничного слоя. Увеличение числа Re также приводит к уменьшению зоны отрыва на верхней стенке канала и к увеличению зоны отрыва пограничного слоя на нижней стенке канала. Об этом факте также можно судить по рис. 45. Скорость течения (красная линия, помечена цифрой 1) значительно падает и после прохождения через скачок уплотнения сохраняет только примерно 5% от своей первоначальной величины. Это и говорит о том, что рассматриваемая область находится в пограничном слое. Совсем другую картину можно наблюдать для максимального числа Рейнольдса.

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

Рис. 45. Распределение x-компоненты скорости в сечении Y=0. На рис. 46 и 47 приведены картины течения, около источника вдува для чисел Рейнольдса 104 и 103. Существенных различий в структурах течений нет: около источника вдува газа формируются зоны отрывного течения, на верхней стенке канала, в области взаимодействия скачка уплотнения и пограничного слоя также происходит отрыв пограничного слоя. Однако из приведенных рисунков можно сделать вывод и о том, что уменьшение числа Рейнольдса приводит к формированию более мощных структур вихревого течения (сравнительные размеры зоны отрывного течения на верхней стенке канала), а также в случае числа Рейнольдса 103 образованию перед источником вдува газа системы из двух мощных зон вихревого течения. Как уже говорилось, образование таких зон подтверждается расчетами других авторов и данными эксперимента. Направление движения газа в этих зонах является особенностью течения и также подтверждается в работах [6,14].

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

Рис. 47. Отрыв пограничного слоя и вихрь на верхней стенке канала, Re В качестве заключения к результатам моделирования течения газа в канале с источником вдува с варьированием чисел Маха и Рейнольдса в невозмущенном потоке можно сделать следующие выводы:

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

Достоверность расчетов подтверждена расчетами экспериментальными данными.

2. Получены зависимости угла наклона головного скачка уплотнения и возникновения зон отрыва пограничного слоя на стенках канала от числа Маха набегающего потока.

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

Глава IV. Моделирование сверх- и гиперзвукового обтекания элементов летательного аппарата Данная глава посвящена численному исследованию плоского стационарного течения вязкого сжимаемого теплопроводного газа вблизи элементов модельного гиперзвукового летательного аппарата. Численное исследование задачи проводилось в рамках системы уравнений НавьеСтокса. Поставленная начально-краевая задача решалась методом установления по схеме с оптимальным расщеплением операторов. Основное внимание уделено изучению структуры течения, образующегося около элементов ЛА, влиянию определяющих параметров газа в невозмущенном потоке на характер течения, а также влиянию краевых условий на характер течения. Течение предполагалось ламинарным.

4.1 Постановка задачи.

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

Рис. 48. Модель гиперзвукового летательного аппарата (а) и расчетная область (б): 1, 2, 3 – точки первого, второго, третьего излома соответственно; ВЗ – На тело (рис. 48б, представлено поперечное сечение ЛА в плоскости симметрии, верхняя часть на рисунке не приведена), движущееся со сверхзвуковой скоростью, слева набегает однородный поток газа под различными углами атаки. Поток, взаимодействуя с носовой частью летательного аппарата, попадает в канал воздухозаборника, проходя который попадает в расширяющую поверхностью кормовой части аппарата и областью свободного истечения газа из канала двигателя.

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

x 3.0, 3.8, 4.5. Первый участок образующей (до точки первого излома) задается выражением 1 0.057 x, остальные участки – кусочно-линейной функцией. Угол наклона второго участка образующей (между точками первого и второго излома) составляет 15 49', третьего - 16, четвертого - 24.

Вход в канал воздухозаборника на верхней образующей начинается в точке x 5.0. Длина канала воздухозаборника варьируется. Изменение длины канала достигается за счет увеличения длины нижней стенки канала:

выдвижения вперед и продления назад. Для различных расчетов длина канала изменяется от 1.0 до 2.3 единицы. Подобная геометрия летательного аппарата использовалась в работе [84]. В качестве характерного размера выбрана величина экспериментов, представленные в этой главе, выполнены для двух конфигураций вычислительной области: в первом случае расчетная область ограничена линией AA – длина канала составляет 1.0 единиц, нижняя кромка воздухозаборника находится на одной линии с верхней кромкой. Во втором случае расчет осуществляется во всей вычислительной области, содержащей область истечения из канала с верхней твердой границей.

Нижняя стенка канала воздухозаборника продлена вперед до точки x 4.7 и назад до точки x 7.0. Продление нижней стенки канала приводит к тому, что на выходе из канала образуется сверхзвуковое сопло.

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

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

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

теплопроводность задавались как функции температуры по степенному закону с показателем 0.76, вторая вязкость принималась равной нулю, показатель адиабаты 1.4, число Прандтля Pr 0.72. С учетом выбранных безразмерных параметров течение определялось числами Маха и Рейнольдса набегающего потока и углом атаки. В качестве начального распределения задавались значения невозмущенных скорости, плотности и вычислительной области.

В качестве критерия установления в расчетах, представленных в данной главе, выбрано условие:

Все выше сказанное в Главе III о характере поведения условия выхода на стационарный режим (2.6) остается верным и для данного исследования.

Сложная структура течения, которая образуется перед каналом при различных значениях определяющих параметров невозмущенного потока газа, влияет на скорость установления решения, а также на устойчивость схемы в целом. Для выбранного начального распределения абсолютно несогласованного распределения параметров газа, итерационный шаг выбирался из условия наискорейшего установления, варьировался в зависимости от значений параметров газа набегающего потока и составил (0.5 0.8) min(hx, hy ), параметр всюду полагался равным 1.0.

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

4.2 Результаты численных расчетов.

По предложенной в Главе I разностной схеме с оптимальным расщеплением (1.40) были проведены расчеты обтекания элементов в центральном продольном сечении летательного аппарата (в плоскости симметрии ЛА). После серии методических расчетов на последовательности сгущающихся сеток было определено оптимальное количество узлов расчетной области - большая часть расчетов выполнена на разностной сетке, содержащей 326 161 узлов, расчеты выполнены при различных числах Маха и углах атаки набегающего потока, изменении геометрии летательного аппарата и изменении краевых условий на несущей поверхности.

4.2.1 Расчеты течений при варьировании чисел Маха Обсуждение результатов расчетов и их детальный анализ начнем с рассмотрения результатов численного моделирования течений вблизи элементов ЛА при Re 104, 6 и изменением числа Маха набегающего потока M 0 4,6,8.

Рис. 49. Распределение плотности газа при различных числах Маха: а – На рис. 49 представлено поле плотности газа около исследуемого тела при различных числах Маха. Как хорошо видно из рисунка, при увеличении скорости набегающего потока угол наклона головного скачка уплотнения, возникающего при сверхзвуковом обтекании, с ростом числа Маха уменьшается. Образующаяся головная волна уплотнения около носка тела при числе M 0 4 не попадает в канал (рис. 49,а), при M 0 6 волна уплотнения попадает на нижнюю кромку воздухозаборника (рис. 49,б), а при M 0 8 попадает внутрь канала воздухозаборника (рис. 49,в). В точках излома верхней образующей возникают скачки уплотнения (см., например, рис 49,а), взаимодействующие с головным скачком уплотнения. Вблизи взаимодействующий с падающей головной волной. Плотность потока и давление в начальном сечении воздухозаборника резко возрастают. Ударная волна, распространяющаяся от нижней кромки канала, падая на верхнюю стенку воздухозаборника, вызывает отрыв потока, интенсивность которого зависит от числа Маха набегающего потока. С увеличением числа Маха увеличиваются размеры зоны возвратного течения, которая образуется на верхней стенке канала. О размерах и пространственном положении зон возвратного течения можно судить по рис. 50.

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

Для числа Маха набегающего потока M 0 4 высота зоны возвратного течения составляет 10%, для числа Маха M 0 6 уже 30%, а для M 0 8 около 50% высоты всего канала.



Pages:     || 2 |


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

«ДОМНИНА Алиса Павловна ЭНДОМЕТРИАЛЬНЫЕ СТВОЛОВЫЕ КЛЕТКИ: ПОЛУЧЕНИЕ, ХАРАКТЕРИСТИКА И ПРИМЕНЕНИЕ ДЛЯ СТИМУЛЯЦИИ РАЗВИТИЯ ЭНДОМЕТРИЯ КРЫС 03.03.04. – Клеточная биология, цитология, гистология ДИССЕРТАЦИЯ на соискание ученой степени кандидата биологических наук Научный руководитель : д.б.н., академик, Никольский Николай Николаевич Санкт – Петербург ОГЛАВЛЕНИЕ 1. СПИСОК СОКРАЩЕНИЙ....»

«БИКСОЛТ АЛЕКСАНДРА МОИСЕЕВНА ОРГАНИЗАЦИОННО-ПРАВОВЫЕ МЕХАНИЗМЫ ОГРАНИЧЕНИЯ КУРЕНИЯ ТАБАКА В РОССИЙСКОЙ ФЕДЕРАЦИИ 14.02.03. Общественное здоровье и здравоохранение Диссертация на соискание ученой степени кандидата медицинских наук Научный руководитель : академик...»

«Еременко Сергей Леонидович ЭКОНОМИЧЕСКОЕ ПОВЕДЕНИЕ РОССИЯН В ГЛОБАЛЬНОЙ КОМПЬЮТЕРНОЙ СЕТИ ИНТЕРНЕТ: СОЦИОЛОГИЧЕСКИЙ АНАЛИЗ 22.00.04 – социальная структура, социальные институты и процессы ДИССЕРТАЦИЯ на соискание ученой степени кандидата социологических наук Научный руководитель – доктор социологических наук Е.О. Кубякин Краснодар – Содержание Введение.. 1. Экономическое поведение россиян...»

«Ватолина Юлия Владимировна Гостеприимство как антропологический феномен Специальность: 09.00.13 – Философия культуры, философская антропология Диссертация на соискание ученой степени доктора философских наук Научный консультант доктор философских наук, профессор Б. В. Марков Санкт-Петербург 2014 2 СОДЕРЖАНИЕ ВВЕДЕНИЕ ГЛАВА 1. ГОСТЕПРИИМСТВО В ПРЕДЕЛЬНОЙ ИНТЕРПРЕТАЦИИ: СЦЕНАРИЙ ДЕКОНСТРУКЦИИ ГЛАВА 2. Я-ЧУЖОЙ: ПРЕДЕЛЫ...»

«Зайцев Владислав Вячеславович РАЗРАБОТКА И ИССЛЕДОВАНИЕ МЕТОДИКИ ПРОЕКТИРОВАНИЯ БАЗЫ МЕТАДАННЫХ ХРАНИЛИЩА ГЕОДАННЫХ Специальность 25.00.35 – Геоинформатика ДИССЕРТАЦИЯ на соискание ученой степени кандидата технических наук Научный руководитель д-р техн. наук, проф. А.А. Майоров Москва ОГЛАВЛЕНИЕ...»

«ТРУСОВА ВАЛЕНТИНА ВАЛЕРЬЕВНА ОЧИСТКА ОБОРОТНЫХ И СТОЧНЫХ ВОД ПРЕДПРИЯТИЙ ОТ НЕФТЕПРОДУКТОВ СОРБЕНТОМ НА ОСНОВЕ БУРЫХ УГЛЕЙ Специальность 05.23.04 – Водоснабжение, канализация, строительные системы охраны водных ресурсов ДИССЕРТАЦИЯ на соискание ученой степени кандидата технических наук Научный руководитель : доктор технических наук В.А. Домрачева ИРКУТСК ОГЛАВЛЕНИЕ...»

«ВЯТКИНА ВАЛЕРИЯ ВАЛЕРЬЕВНА ФОРМИРОВАНИЕ ЦЕЛЕВЫХ КОМПЛЕКСНЫХ ПРОГРАММ СОЦИАЛЬНОЙ ЗАЩИТЫ ДЕТЕЙ-СИРОТ В РЕГИОНЕ Специальность 08.00.05 – Экономика и управление народным хозяйством: региональная экономика; экономика, организация и управление предприятиями, отраслями, комплексами: сфера услуг...»

«UNIVERSITE JOSEPH FOURIER UNIVERSITE D’ETAT DE ROSTOV-SUR-LE DON DOCTORAT Physique de la Matire Condense et du Rayonnement Elena Nazarenko tel-00136821, version 1 - 15 Mar 2007 Structures locales de la magntite et de zirconates de type perovskite par diffraction rsonante et absorption X Thse dirige par Yves Joly et Rostislav Vedrinskii Date de la soutenance: le 25 janvier JURY L. Bugaev V. Dmitriev Rapporteur Y. Gufan Rapporteur Y. Joly K. Protassov Prsident R. Vedrinskii МИНИСТЕРСТВО...»

«БУШУЕВ Юрий Гениевич СТРУКТУРНЫЕ СВОЙСТВА ЖИДКОСТЕЙ С РАЗЛ ИЧНЫМИ ТИПАМИ МЕЖМОЛ ЕКУЛЯРНЫХ ВЗАИМОДЕЙСТВ ИЙ ПО ДАННЫМ КОМПЬЮТЕРНОГО МОДЕЛ ИРОВ АНИЯ 02.00.04 – физическая химия Диссертация на соискание ученой степени доктора химических наук Иваново 2001 ОГЛАВЛЕНИЕ стр. ВВЕДЕНИЕ 7 1. ПРИМ ЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ МОДЕЛИРОВАНИЯ СТРУКТУРЫ ЖИДКОСТЕЙ 1.1. Общие теоретические положения 1.2. Алгоритм Метрополиса 1.3....»

«ЗУЙКОВА АННА АЛЕКСАНДРОВНА ОСОБЕННОСТИ ЭТИОПАТОГЕНЕЗА ДЕЗАДАПТИВНЫХ ИЗМЕНЕНИЙ УЧАСТНИКОВ БОЕВЫХ ДЕЙСТВИЙ 14.03.03 – патологическая физиология диссертация на соискание ученой степени доктора медицинских наук Н.Новгород 2014 1 Оглавление. Введение. Глава 1. Современное представление об адаптации участников локальных вооруженных конфликтов после воздействия боевого стресса и травм. Глава 2. Исследование дезадаптивных изменений после воздействия боевого стресса. 2.1...»

«из ФОНДОВ РОССИЙСКОЙ ГОСУДАРСТВЕННОЙ БИБЛИОТЕКИ Сергеева, Галина Георгиевна 1. Прецедентные имена и понимание ик в молодежной среде 1.1. Российская государственная Библиотека diss.rsl.ru 2005 Сергеева, Галина Георгиевна Прецедентные имена и понимание ик в молодежной среде [Электронный ресурс]: Школьники 10-11 класса : Дис.. канд. филол. наук : 10.02.19.-М.: РГБ, 2005 (Из фондов Российской Государственной Библиотеки) Теория языка Полный текст: http://diss.rsl.ru/diss/05/0377/050377020.pdf...»

«ДМИТРИЕВА Диана Михайловна ПОВЫШЕНИЕ ЭФФЕКТИВНОСТИ СТРАТЕГИЧЕСКОГО УПРАВЛЕНИЯ ИНТЕГРИРОВАННОЙ КАЛИЙНОЙ КОМПАНИЕЙ Специальность 08.00.05 - Экономика и управление народным хозяйством (экономика, организация и управление предприятиями, отраслями, комплексами) ДИССЕРТАЦИЯ...»

«КАБИРОВ Валентин Рамильевич ОЦЕНКА ЭКОНОМИЧЕСКОЙ ЭФФЕКТИВНОСТИ РАЗРАБОТКИ ГРУППЫ ТЕРРИТОРИАЛЬНО-СБЛИЖЕННЫХ РУДНЫХ (МЕТАЛЛИЧЕСКИХ) МЕСТОРОЖДЕНИЙ Специальность 08.00.05 – Экономика и управление народным хозяйством (экономика, организация и управление предприятиями,...»

«Анкудинова Полина Михайловна ЭВОЛЮЦИОННОЕ СТАНОВЛЕНИЕ ЧЕЛОВЕКА В ФИЛОСОФСКО-АНТРОПОЛОГИЧЕСКИХ КОНЦЕПЦИЯХ ХХ ВЕКА: МИРОВОЗЗРЕНЧЕСКИЙ АСПЕКТ Специальность 09.00.13 – философская антропология, философия культуры Диссертация на соискание ученой степени кандидата философских наук Научный руководитель : доктор философских наук,...»

«Достовалов Дмитрий Николаевич СПЕЦИФИКАЦИЯ И ИНТЕРПРЕТАЦИЯ МОДЕЛЕЙ ПЕРЕХОДНЫХ ПРОЦЕССОВ В СИСТЕМАХ ЭЛЕКТРОЭНЕРГЕТИКИ 05.13.11 – Математическое и программное...»

«Кириллов Андрей Геннадьевич ПОЛИТИЧЕСКИЙ НАРРАТИВ: СТРУКТУРА И ПРАГМАТИКА (на материале современной англоязычной прессы) Специальность 10.02.04 – германские языки ДИССЕРТАЦИЯ на соискание ученой степени кандидата филологических наук Научный руководитель : доктор культурологии, профессор М....»

«Баклыков Герман Евгеньевич ПОВЫШЕНИЕ КАЧЕСТВА ФУНКЦИОНИРОВАНИЯ ПРОИЗВОДСТВЕННЫХ ПРЕДПРИЯТИЙ НА ОСНОВЕ СОВЕРШЕНСТВОВАНИЯ СИСТЕМЫ ТОВАРОДВИЖЕНИЯ И УПРАВЛЕНИЯ МАТЕРИАЛОПОТОКАМИ 08.00.05 – Экономика и управление народным хозяйством (стандартизация и управление качеством продукции) Диссертация на соискание ученой степени кандидата экономических наук...»

«Костюкевич Юрий Иродионович Компенсационные ионные ловушки с динамической гармонизацией для масс-спектрометра ионного циклотронного резонанса 01.04.17 – химическая физика, горение и взрыв, физика экстремальных состояний вещества диссертация на соискание ученой степени кандидата физико-математических наук Научный руководитель :...»

«ДЮЛИЧЕВА Юлия Юрьевна УДК 519.68 МОДЕЛИ КОРРЕКЦИИ РЕДУЦИРОВАННЫХ БИНАРНЫХ РЕШАЮЩИХ ДЕРЕВЬЕВ 01.05.01 – Теоретические основы информатики и кибернетики Диссертация на соискание ученой степени кандидата физико-математических наук Научный руководитель : доктор физико-математических наук, профессор ДОНСКОЙ Владимир Иосифович Симферополь – ОГЛАВЛЕНИЕ ВВЕДЕНИЕ.. Раздел 1. Методы синтеза и редукции решающих деревьев: обзор и...»

«РАДЬКО Сергей Иванович РАЗРАБОТКА И ИССЛЕДОВАНИЕ ЭЛЕКТРОТЕХНОЛОГИЧЕСКОГО ОБОРУДОВАНИЯ ДЛЯ ПЕРЕРАБОТКИ ТЕХНОГЕННЫХ ОТХОДОВ С ИСПОЛЬЗОВАНИЕМ ПАРОВОДЯНОГО ПЛАЗМОТРОНА Специальность: 05.09.10 – Электротехнология Диссертация на соискание ученой степени кандидата технических наук Научный руководитель – доктор технических наук,...»






 
2014 www.av.disus.ru - «Бесплатная электронная библиотека - Авторефераты, Диссертации, Монографии, Программы»

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