МОДИФИКАЦИИ МЕТОДА РАСЩЕПЛЕНИЯ РЕШЕНИЯ ЗАДАЧ

ГАЗОВОЙ ДИНАМИКИ

Слюняев А.Ю.

Институт вычислительных технологий СО РАН, г. Новосибирск

 

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

При построении неявных разностных схем наиболее часто используются методы факторизации и расщепления (см., [1-7]), которые позволяют свести решение исходных многомерных задач к решению их одномерных аналогов или более простых задач. Решение одномерных задач по неявным разностным схемам методов приводит к необходимости обращения матриц, размерность которых зависит от числа уравнений в системе, а также от размерности задачи по пространству. Понятно, что решение задач по таким схемам также не является экономичным, поэтому представляется целесообразным построение таких численных алгоритмов, реализация которых даже в одномерном случае является экономичной и сводится, например, к скалярным прогонкам. Известно [2], что расщепление задачи или факторизация операторов в исходной многомерной задаче приводит к появлению дополнительных членов в разностной схеме – диссипативных членов и членов более высокого порядка и, как следствие, к ухудшению свойств численного алгоритма.

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

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

                                                                                                                                 (1)

где

здесь  – плотность,  – скорость,  – давление,  – температура газа,  - полная энергия газа,  – внутренняя удельная энергия газа,  – удельная теплоемкость газа при постоянном объеме. Система уравнений Эйлера замыкается уравнением состояния совершенного газа:

Система уравнений Эйлера также может быть записана в недивергентном векторном виде:

                                                                                                          (2)

здесь

.

Покажем связь между дивергентной и недивергентной формами записи уравнений Эйлера. Из (1), учитывая, что вектор функций  зависит от вектора функций , запишем:

обозначим  и , тогда получим:

отсюда получим связь между дивергентной и недивергентной формами записи уравнений Эйлера:

                                                                                                                           

Разностная схема. Вводя в расчетной области  разностную сетку и аппроксимацию дифференциальных операторов разностными операторами, перейдем от непрерывной модели (1) к конечно-разностному аналогу. Разностная схема примет вид:

                                                                                                        (3)

Она аппроксимирует систему уравнений (2) с порядком , а в случае установления – с порядком . В линейном приближении схема безусловно устойчива при  [2].  Матричный оператор , аппроксимирующий оператор , перепишется в следующем виде:

                                                                                                           

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

                                                                                                                              (4)

                                                                                                                       (5)

                                                                                                                             (6)

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

                                                                                 (7)

Вид операторов  определяется расщеплением исходного оператора по физическим процессам:

где матричные операторы запишутся в виде:

                                                                   

Разностная схема (7) аппроксимирует исходную систему уравнений (2) с порядком , а при установлении – с порядком , в линейном приближении при  является абсолютно устойчивой (см. [2]). Для реализации разностной схемы (7) рассмотрим эквивалентную ей схему в дробных шагах:

                                                                                                                              (8)

                                                                                                                 (9)

                                                                                                                           (10)

Разностная схема (7) в дробных шагах (8) – (10) реализуется скалярными прогонками [9].

Аппроксимация производных в разностных схемах проводилась разностным оператором  ( - порядок аппроксимации,  - пространственное направление) по следующим схемам: на дробных шагах аппроксимация первых производных выполнена по противопотоковой схеме с первым порядком точности [4], на нулевом дробном шаге аппроксимация производных выполнена по схеме ():

                     

 

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

Коротко рассмотрим результаты тестовых расчетов для задачи о распаде произвольного разрыва. Труба постоянного сечения заполнена покоящимся газом и перегорожена заслонкой. По обе стороны для газа определены параметры – давление, скорость и плотность. В момент времени  заслонка мгновенно убирается. Требуется описать изменение во времени параметров газа. Подробное описание задачи дано в [5]. Параметры для расчетов задавались следующие: число узлов на расчетной сетке составило , решение приведено на момент времени . На рисунках 1 и 2 представлено численное распределение плотности, полученное по предложенной схеме. Численное решение (дано треугольниками) достаточно хорошо повторяет структуру точного решения (сплошная красная линия). Как видно из рисунков, ширина переходной зоны на ударной волне составляет всего 4 узла. На рисунке 3 представлены распределения скорости газа, полученные с первым и вторым порядком точности. Очевидно, что решение полученное со вторым порядком (дано кружками) точности значительно лучше повторяет структуру течения и имеет меньшую зону «размазывания» на ударной волне, чем решение, полученное с первым порядком (дано квадратиками).

 

Рис. 1 Распределение плотности

Рис. 2 Распределение скорости

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

 

 

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

 

Рис. 4 Сопло Лаваля.

 

Если в задаче на входе задавать дозвуковое течение, а на выходе произвольное течение, то решение уравнений квазиодномерной задачи так же существует, но оно будет разрывным [5]. Для решения этой задачи в рамках одномерной системы уравнений (1) вводится замена переменных, стационарная система уравнений Эйлера изменит свой вид на следующий [5]:

Условие останова для итерационного процесса полагалось следующее:  Расчеты проводились на сетке в 401 узел. В таблице 1 приведена оценка относительной погрешности вычислений в процентах для плотности, скорости и давления газа, число итераций, необходимых схеме для установления по заданному критерию останова, а также максимальное число Куранта (, где  – скорость звука), с которым схема позволяет проводить расчет.

Таблица 1. Относительная погрешность вычислений.

Количество итераций до установления

Число Куранта

0,6

0,4

0,3

330

7,21

 

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

 

Рис. 5 Распределение плотности в канале

Рис. 6 Распределение скорости в канале

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

 

 

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

Исходная система уравнений. Двумерный случай. Полная система уравнений Навье-Стокса для вязкого сжимаемого теплопроводного газа в дивергентном векторном виде в декартовой системе координат для двумерного случая может быть записана в следующем виде [6, 8, 9]:

                                                                                                      

где

                                                           

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

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

                                                                                                            

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

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

                                                              ,                                                             

где

                   

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

Дадим обобщение схемы (7) на уравнение Навье-Стокса для двумерного случая. В переменных  существует единственное расщепление матричных операторов:

                                                                                                                               (11)

где

                 

элементы  имеют вид:

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

                                                                                 

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

                                                                         (12)

Она аппроксимирует исходную систему уравнений с порядком , а при установлении – с порядком , в линейном приближении при  является абсолютно устойчивой для  [3]. Для реализации разностной схемы (12) рассмотрим эквивалентную ей схему в дробных шагах:

                                                                                                                            (13)

                                                                                                             (14)

                                                                                                                           (15)

Разностная схема (13) – (15) реализуется скалярными прогонками.

Течение газа в канале типа «воздухозаборник». Рассмотрим результаты расчетов, которые были проведены для задачи о течении вязкого сверхзвукового теплопроводного газа в канале типа «воздухозаборник». Канал представляет собой плоский канал, нижняя стенка которого начинается в точке , а вторая снесена по оси  относительно него на . Высота канала составляет , длина канала составляет  единиц. На стенках канала поставлены условия прилипания газа к стенкам, а так же условия теплоизоляции стенок канала. На выходе из канала поставлены мягкие условия [5]. На вход канала подается невозмущенный сверхзвуковой поток с параметрами  под нулевым углом атаки. Расчет проведен на сетке  узел (), итерационный шаг , параметр схемы , установление получено за 2500 итераций, время, затраченное на расчет, составило 122 минуты.

 

Рис. 8 Распределение модуля вектора скорости в канале

Рис. 9 Распределение температуры в канале

 

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

 

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

Рис. 11 Изолинии давления в канале

 

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

 

Рис. 12 Изолинии плотности в канале. Подобласть канала

Рис. 13 Изолинии давления в канале. Подобласть канала

 

Более подробно структуру течения канала и количественные характеристики можно представить по рисункам 14 и 15, на которых представлены распределения числа Маха в сечениях , а также плотность в тех же сечениях и дополнительно на стенках канала.

 

Рис. 14 Распределение числа Маха в продольных сечениях канала

Рис. 15 Распределение плотности газа в продольных сечениях канала

 

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

На рисунках 16 и 17 представлены поперечные распределения числа Маха и плотности газа. Из рисунков видно, что при продвижении газа по каналу толщина дозвуковой зоны газа увеличивается – от  в сечении  (на рисунках решение дано красным цветом), до в сечении .

 

Рис. 16 Распределения чисел Маха в поперечных сечениях канала

Рис. 17 Распределение плотности в поперечных сечениях канала

 

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

 

 

 

 

 

Рис. 18 Изолинии давления в канале для различных чисел Маха

Рис. 19 Распределение компонент вектора скорости в поперечном сечении

 

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

 Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант 05-01-00146) и Интеграционного проекта № 116.

 

 

Литература:

  1. Гильманов А.Н. Методы адаптивных сеток в задачах газовой динамики. – М.: Физматлит, 2000. – 247 с.
  2. Ковеня В.М. Разностные методы решения многомерных задач: Курс лекций. – Новосибирск: НГУ, 2004. – 146 с.
  3. Ковеня В.М., Лебедев А.С. Модификации метода расщепления для построения экономичных разностных схем //Ж. вычислительной математики и математической физики. 1994. Т. 34. – №6. – с. 886 – 897.
  4. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. – Новосибирск: Наука. Сиб. отд-ние, 1981. – 304 с.
  5. Лебедев А.С., Черный С.Г. Практикум по численному решению уравнений в частных производных. - Новосибирск: НГУ, 2000. – 136 с.
  6. Лойцянский Л.Г. Механика жидкости и газа. – М.: Наука, 1970. – 904 с.
  7. Марчук Г.И. Методы расщепления. – М.: Наука, 1988. – 264 с.
  8.  Роуч П. Вычислительная гидродинамика. – М.: Мир, 1980. – 616 с.
  9. Слюняев А.Ю. Об одном численном алгоритме решения уравнений Навье-Стокса вязкого сжимаемого газа. «Исследования и перспективные разработки в авиационной промышленности», статьи и материалы конференции. – М.: ОАО «ОКБ Сухого», 2005 – с. 105 – 111.