Бикомпактные схемы для многомерных гиперболических уравнений и их эффективная реализация Чикиткин Александр Викторович

Содержание к диссертации

Введение

Глава 1. Построение бикомпактных схем произвольного порядка аппроксимации 17

1.1 Вывод формул для коэффициентов схем 18

1.2 Тестовые расчёты

1.2.1 Скалярное уравнение переноса 22

1.2.2 Система уравнений акустики 25

Глава 2. Бикомпактные схемы для неоднородных уравнений переноса с источником 28

2.1 Вывод уравнений схемы 29

2.2 Гибридная схема

2.2.1 Построение гибридной схемы на основе схемы низкого порядка 34

2.2.2 Монотонная схема с наименьшей схемной вязкостью 36

2.3 Тестовые расчёты 42

Глава 3. Бикомпактные схемы для многомерных уравнений газовой динамики 50

3.1 Полудискретная бикомпактная схема 50

3.2 Полностью дискретные бикомпактные схемы

3.2.1 Схема первого порядка по времени 52

3.2.2 Схема третьего порядка по времени 53

3.2.3 Схема переменного порядка по времени 53

3.3 Тестовые расчёты 58

4 Глава 4. Эффективная реализация бикомпактных схем 65

4.1 Итерационная факторизация 65

4.1.1 Приближенная факторизация разностных уравнений 65

4.1.2 Сходимость итерационного метода 68

4.1.3 Тестовые расчёты 69

4.2 Параллельная реализация вычислительного алгоритма 71

4.2.1 Параллельный алгоритм бегущего счёта 71

4.2.2 Параллельный алгоритм факторизованной схемы 73

Заключение 78

Список основных обозначений 79

Список литературы

• Тестовые расчёты
• Построение гибридной схемы на основе схемы низкого порядка
• Схема первого порядка по времени
• Сходимость итерационного метода

Введение к работе

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

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

Также набирают популярность методы спектрального типа. В этих методах повышение порядка аппроксимации достигается за счёт применения интерполяции высокого порядка внутри ячейки и за счёт использования неравномерных сеток внутри ячеек (Чебышёвских сеток, узлов квадратур Гаусса-Лобатто и т.п.). Главное преимущество таких методов в том, что они дают очень высокую точность на гладких решениях даже при небольшом порядке точности. Это свойство является ключевым в задачах, описывающих распространение сигналов и волновые процессы.

Диссертация посвящена исследованию свойств, обобщению и эффективной реализации бикомпактных разностных схем, предложенных в работе (Рогов: 2010). Эти компактные схемы совмещают в себе свойства конечно-

разностных, конечно-объёмных и спектральных методов и обладают рядом полезных свойств и преимуществ:

Консервативность;

Простота с точки зрения формулировки и реализации;

Возможность повышения порядка точности без увеличения шаблона;

Безусловная устойчивость, пригодность для решения жёстких задач;

Минимальность шаблона, удобство параллельной реализации;

Высокая точность на гладких решениях, выше, чем у методов TVD, ENO, WENO высших порядков при решении некоторых задач;

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

Цели работы

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

2. Исследование свойств построенных схем.

3. Эффективная реализация алгоритмов на современных многоядерных процессорах.

4. Исследование эффективности схем на тестовых прикладных задачах.

Задачи работы

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

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

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

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

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

Основные положения, выносимые на защиту

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

10. Построены гибридные схемы на основе бикомпактных схем с интегральным средним по ячейке для многомерных неоднородных уравнений переноса.

11. Построены бикомпактные схемы и гибридные схемы на их основе для двумерных уравнений газовой динамики.

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

13. Доказана сходимость итерационного метода на основе приближенной факторизации.

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

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

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

используется наиболее простой подход геометрического разбиения расчётной области.

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

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

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

Работа над диссертацией была поддержана грантом Правительства РФ по постановлению N 220 по договору N 11.G34.31.0072, заключенного между Министерством образования и науки РФ, ведущим ученым и МФТИ, грантом РФФИ в рамках научного проекта N 16-31-00296 мол_а, а также стипендией Президента молодым учёным и аспирантам на 2016-2018 годы.

Степень достоверности и апробация работы Доказательства свойств исследуемых методов напрямую следует из основных результатов теории чис-

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

16. Международная конференция «Разностные схемы и их приложения», по-свящённая 90-летию профессора В.С. Рябенького, Москва, ИПМ, 27-31 мая 2013 года

17. 57-я научная конференция МФТИ, г. Долгопрудный, МФТИ, 24.11.2014-29.11.2014

18. European Workshop on High Order Nonlinear Numerical Methods for Evolutionare PDEs: Theory and Applications, University of Trento, Италия, 16.03.2015-20.03.2015

Публикации по теме диссертации Основные результаты диссертации опубликованы в 7 работах, 4 из которых – в изданиях из списка ВАК.

Личный вклад Все результаты работы, кроме отдельно оговоренных случаев, получены лично диссертантом под научным руководством д.ф.-м.н. Б. В. Рогова. В работах из списка ВАК с соавторами вклад соискателя состоит в следующем: в работе [1] – обобщение бикомпактных схем на двумерные нелинейные уравнения газовой динамики, сравнение различных формул для весового коэффициента в гибридной схеме и выбор наилучшего варианта, реализация алгоритма, тестовые расчёты двумерных задач; в работе [2] – построение гибридной схемы для стационарного неоднородного уравнения переноса на основе бикомпактных схем с интегральным средним, сравнительный анализ разных видов гибридизации, проведение тестовых расчётов; в работе [3] – построение гибридной схемы для нестационарного многомерного неоднородного уравнения переноса на основе бикомпактной схемы с интегральными средними, построение метода коротких характеристик для многомерного случая, исследование диссипативных свойств методов низкого порядка, тестовые расчёты; в работе [4] – разработка параллельного алгоритма, программная реализация параллельного алгоритма, проведение тестовых расчётов на современных многоядерных процессорах, исследование масштабируемости алгоритма.

Структура и объем работы Диссертация состоит из введения, четырех глав, заключения, списка основных обозначений, списка иллюстраций, списка таблиц и приложения. Основной текст занимает 89 страниц.

Тестовые расчёты

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

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

В спектральном методе = 1 и матрица линейной системы имеет (-1) ненулевых элементов. Если теперь предположить, что вычислительные затраты пропорциональны квадрату числа ненулевых элементов, что справедливо, если система решается итерационным методом, становится очевидным, что следует выбирать тот или иной метод, учитывая точность и затраты одновременно.

Например, чтобы получить точность 10-12 в данной задаче, необходимо решать задачу методом 4-го порядка с 1024 ячейками. На один шаг по времени уходит 1.2 10-2 c. Для достижения той же точности можно решить эту задачу на одной ячейке методом 28-го порядка с = 15, один шаг по времени – 3.310-4 с. Таким образом, для достижения данной точности в данной задаче, предпочтителен спектральный метод, он даёт выигрыш во времени в 100 раз.

Для достижения точности 10-3 метод 4-го порядка на сетке из 8 ячеек тратит уже 510-4 с на шаг. Спектральный метод c = 7 – 210-4 с на шаг. В этом случае удобнее пользоваться методом 4-го порядка, как более простым. Для других задач метод другого порядка может оказаться оптимальным с точки зрения вычислительных затрат.

Зависимость ошибки от числа узлов в логарифмической шкале для метода 4-го порядка и коллокационного метода на одной ячейке

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

На рисунках 2,3 показано сравнение точного решения с численным решением в момент времени = 1. Красными маркерами изображены граничные (целые) узлы расчётных ячеек. В первом случае используется метод с 6-ю узлами в ячейке (формально 10-й порядок), всего 6 ячеек, во втором случае – метод с 28-ю узлами в ячейке (формально 54-й порядок), с одной ячейкой. В обоих случаях ошибка мала. Видно, что в задачах с гладким решением методы спектрального типа с большим числом степеней свободы внутри расчётной ячейке позволяют получить очень высокую точность на маленькой сетке.

Построение гибридной схемы на основе схемы низкого порядка

Полудискретные схемы для уравнения переноса (33) выводятся методом прямых путем дискретизации пространственных производных на минимальном шаблоне. В 1D случае пространственный шаблон состоит из двух узлов Xj,Xj+\, в 2D случае - из четырех узлов (xj,yk),(xj+i,yk),(xj,yk+i),(xj+\,yk+i), а в 3D случае - из восьми узлов (xa,y ,Zj), а = j,j + 1, /3 = к,к -\- I, 7 = , + 1. Благодаря тому, что в каждом пространственном направлении шаблон является двухточечным, схемы названы бикомпактными. Вывод схемы рассмотрим для 3D случая. Поскольку он будет проводиться на минимальном пространственном шаблоне, то он справедлив в общем случае неравномерных декартовых сеток с шагами hx = hx,j + 1/2 = xj+i - xj, hy = hytk+i/2 = yk+i - Ук, hz = hz i+1/2 = zi+1 - zh первой и второй производных. Далее для простоты Для получения полудискретной схемы повышенного порядка аппроксимации по простран-ственным переменным наряду со значениями искомой функции и в узлах шаблона вводятся вспомогательные величины: 1. интегральные средние по 12 ребрам 3D ячейки G = [XJ,XJ+I] х [ук,Ук+і] х [ г,- г+і] Mj+l/2,/3,7 и имеют четвертый порядок аппроксимации по пространственным переменным. Уравнения (40a)-(40h) записываются для полуцелых узлов, центров пространственных ячеек (xj+i/2,yk+i/2, Zi+\/2). Операторы А есть операторы осреднения по ячейке, операторы Лі,Л2 есть разностные аналоги изложения будут представлены численные эксперименты с использованием 2D бикомпактных схем. Полудискретная форма этих схем для двумерного аналога уравнения (33) с с = 0 состоит из четырех обыкновенных дифференциальных уравнений Ах Аущ + аАуАхи + ЬА\Ахи + Ау Ах ои = AyAxQ (43a) АуАхщ + аАу Ах2и + ЪА\А\и + АуАхаи = AyAxQ (43b) А\Ахщ + aA\Axlu + ЬАуАхи + А\Ахаи = A\AXQ (43c) А\Ахщ + аА\А2и + ЬАуАхи + А\Ахаи = A\AXQ (43d) Уравнения бикомпактной схемы с интергральным средним по ячейке не могут быть получены из условий коллокации, в отличие от бикомпактной схемы со значением в полуцелом узле. Но в одномерном случае бикомпактную схемы с интегральным средним по ячейке можно также рассматривать как метод Рунге-Кутты решения ОДУ.

Система ОДУ (43a) - (43d) не содержит пространственных производных и определяет бикомпактную схему метода прямых для линейного неоднородного уравнения переноса в случае двух пространственных направлений. Как и для случая одной пространственной переменной, полученная динамическая система четырех уравнений вида Vt + Av = F, (51) может быть решена маршевым по дискретным переменным методом без пространственного расщепления. При выбранных знаках коэффициентов а 0 и b 0 направление маршевого счета –– в сторону увеличения независимых переменных. Возможность бегущего счета подтверждается структурой матрицы перехода А полученной системы, собственные значения которой определяются следующей формулой:

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

Опишем алгоритм построения гибридного решения для нестационарного случая. На данном шаге по времени tn сначала вычисляется решение во всей области по бикомпактной схеме с интегральными средними, которое мы обозначим ив. В каждой пространственной ячейке содержится 9 значений: 4 значения в вершинах ячейки, 4 значения интегральных средних по рёбрам и интегральное среднее по всей ячейке. Те же величины вычисляются по монотонной схеме низкого порядка точности, обозначим их и . Вопрос об оптимальном выборе схемы A будет рассмотрен ниже. После этого для каждой из девяти пар значений (и , ив) гибридное решение вычисляется по формуле:

Здесь С1 0 - константа, вообще говоря, зависящая от исследуемой задачи, но изменяющаяся в довольно узком диапазоне, m rri0 — 1, где rri0 - порядок аппроксимации по времени, є Ю-10 - вычислительный параметр, позволяющий избежать деления на ноль. Заметим, что формула (53) определяет оператор послойного перехода, который является локальным: в нем используются значения и , ив только в рассчитываемой пространственно-временной точке. Отметим, что используется формула для коэффициента w, отличная от формул, ис-пользовавшихся ранее. Выражение (54) более универсально, т.к. не зависит от шага по времени и содержит локальную нормировку на характерную величину решения \и \. Можно использовать также глобальную нормировку: С1\и — U \ W = : jj (55) maxj;fc \и \ + є здесь максимум берётся по значениям во всех узлах на текущем слое. В приведённых ниже расчетах обе нормировки дают графически неразличимые результаты, по умолчанию используется локальная нормировка.

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

Согласно формулам (53), (54) в областях гладкости решения гибридная схема будет сохра-нять порядок схемы В. Последнее было проверено на ряде тестовых задач. Вблизи разрыва решения и величина ин будет близка к величине и , определяемой по монотонной схеме А низкого порядка аппроксимации в согласии с теоремой Годунова. Поэтому качество численного решения вблизи разрыва (протяженность зоны размазывания разрыва) будет опреде-ляться величиной аппроксимационной вязкости монотонной схемы А. Ниже рассматривается задача, связанная с поиском схемы А с наименьшей диссипацией.

Схема первого порядка по времени

Метод построения компактной схемы опишем на примере численного решения задачи для двумерного квазилинейного уравнения: du du (90) dQ(u) dF{u) dG{u) L2U = 0, L211 = 1 h dt дх dy dQ{u) dp (и) dQ{u) dLr{u) ;—;— 0) Ъ(и) =:—:— 0 du du du du dQ{u) dp \u) du du u(x,y,0) = и (x, y), x 0, у u(0,y,t) = fi1(y,t),u(x,0,t) = fi2(x,t),t 0

Для вывода уравнений бикомпактной схемы можно использовать общий подход, основан-ный на коллокационных многочленах, который описан в главе один, но здесь мы построим компактную схему для уравнения (90) так же, как и в одномерном случае, методом прямых совместно с интегро-интерполяционным методом. Это позволяет получить уравнения схемы из интегральных соотношений для исходного уравнения и его дифференциальных следствий. Для этого в области пространственных переменных П = {(х,у),х 0, у 0} вводится основная сетка шх х ujy, состоящая из целых узлов (xj,yk), где шх = {xj,j 0}, шу = {yk, к 0}. Дополнительно в этой области вводятся вспомогательные сетки шх х ш1, ш х х шу, ш х х ш1, состоящие из полуцелых узлов, где ш х = \Xj+1/2,j 0, ш = \Ук+1/2, к 0. Наряду с основным уравнением (90) рассматриваются также его дифференциальные следствия: dL2U п— = ох dL2U — = 0 (91) оу d2L2U п = дху На шаблоне, состоящем из девяти узлов (XJ, а +і/г, j+i) х (Ук, Ук+і/2, Ук+і), четырёх целых узлов и пяти полуцелых узлов, нетрудно получить четыре дифференциально-разностных уравнения, выполняющиеся с точностью до O hfnax) на точном решении уравнения (90), где hmax = mcLx( ii:J+b j/,fc+l), -icj+l = Xj+\ — Xj. Первое из этих уравнений имеет вид d v ч 1 „ 1 „ — (A0A0Q j-\-l/2,k+l/2 + Т AoAQFj+i/2tk+l/2 + T - -o o(jj+l/2,k+l/2 = 0 (92) u ilx,j-\-l vy}k-\-l и аппроксимирует на сетке основное уравнение, проинтегрированное по двумерной ячейке [XJ,Xj+\] х [ук,Ук+і]. Остальные три уравнения — ( AQAQQj+i/2tk+l/2 + 7 -o 2- J+l/2,fc+l/2 + 7 Ao okj+l/2,fc+l/2 = 0 (93) u ilx,j-\-l vy}k-\-l d .XAV \ 1 Л«ЛЖП 4 ЛІКЛЧ — (A0l\0Qj+i/2:k+l/2) + T 0 0- j +l/2,fc+l/2 + 7 - о 2 +1/2,А;+1/2 = 0 (94) OX vxJ-\-l vy}k-\-l d (ArAV \ 4 v 4 T « — AQAQQj+i/2tk+l/2 + 7 0 2- J+l/2,fc+l/2 + 7 0 2 +1/2 +1/2 = 0 (95) LIL xJ-\-l "y}k-\-l аппроксимируют уравнения (91), также проинтегрированные по этой ячейке. В уравнениях (92)- (95) сеточные операторы AQ, Af, AQ есть двумерные аналоги одномерных сеточных операторов первой производной, второй производной и интегрального среднего по ячейке, они определены на двумерных сеточных функциях, заданных на множестве полуцелых узлов: /\,-,7/ I 1 1 I 1 /О = 7/ I 1 1 I 1 /о 111 I 1 /о О 3 \ i-/ 2,к-\-1/ 2 3 т jk-\-V / 2 3 к г -/ 2 2 J+l/2,/c+l/2 Llj-\-l,k-\-l/2 Llj-\-l/2,k-\-l/2 T Llj,k-\-l/2 ( ) l/2,fc+l/2 V j +l,/c+l/2 r tLj-\-l/2,k-\-l/2 JJj,k-\-l/2 )/v где E - единичный оператор. Приведённые выше выкладки для скалярного уравнения можно полностью повторить и в случае системы уравнений dQ{u) dF{u) dG{u) п 1 n 1 n = 0 (97) at ox ay где u(x,t) - искомая вектор-функция с т компонентами, Q(u),F(u),G(u) - заданные вектор-функции размерности т. В результате можно получить системы обыкновенных диф-ференциальных уравнений (ОДУ), аналогичные (92)- (95).

Заменяя производные по времени в системе ОДУ (92)- (95) какими-либо разностными аппроксимациями, получим семейство компактных разностных схем для численного решения уравнения (90) и системы (97).

Если производные по времени на расчётном слое t = tn+\ аппроксимировать разностями назад, то получим базовую разностную схему, состоящую из четырёх уравнений (нижние индексы (j + 1/2, к + 1/2) опущены у всех величин): AyAxQn+1 + rxJ+1AyAxFn+1 + rVtk+1AxAyGn+1 = AyAxQn AyAxQn+1 + 4rxJ+lAyAxFn+1 + ry,k+lAxAyGn+1 = AyAxQn (98) AxAyQn+1 + rxJ+lAyAxFn+1 + 4ry,k+lAxAyGn+1 = AxAyQn AxAyQn+1 + irxJ+1AyAxFn+1 + Ary,k+1AxAyGn+1 = AxAyQn Здесь rxJ+i = r/hxJ+i, ry k+i = r/hytk+i, T = tn+i - tn - временной шаг. Разностная схема (98) есть неявная схема Эйлера, которая имеет первый порядок аппроксимации по времени и является L-устойчивой, а, следовательно, и абсолютно устойчивой. В дальнейшем верхний индекс п + 1 при величинах на рассчитываемом временном слое опустим.

Отметим, что дифференциально-разностные уравнения (92)- (95) получены интегро-ин-терполяционным методом из дифференциальных уравнений дивергентного типа. Такой способ получения этих уравнений обеспечивает консервативность базовой разностной схемы.

Схема (98) может быть решена с помощью алгоритма бегущего счёта следующим образом. На рассчитываемом временном слое рассмотрим угловую двумерную прямоугольную разностную ячейку с нижним левым узлом (хо,Уо). Из граничных условий известны значения сеточной функции и в пяти узлах (хо, уо), (х\/2, уо), (х\, уо), (хо, уі/г), (%о, Ух). Из четырёх разностных уравнений (98) при j = 0, к = 0 находим значения функции и в остальных четырёх узлах ячейки. Затем переходим к соседней ячейке либо в направлении оси х, либо в направлении оси у и проводим аналогичный расчёт. Фактически схема (98) является бикомпактной в каждом пространственном направлении по терминологии работы [6], полуцелые узлы сетки являются вспомогательными в схеме. Для реализации схемы значения искомой функции в полуцелых узлах нужны лишь на начальном временном слое = 0 и на границах расчётной области. Следствием бикомпактности схемы является сохранение порядка аппроксимации схемы на сетке, состоящей из целых узлов, при переходе от равномерной пространственной сетки к неравномерной сетке. Отметим, что система уравнений (98), не являясь факторизованной схемой, может решаться методом бегущего счета по любому из пространственных направлений, что характерно для факторизованных схем. Это следствие того, что схема имеет первый разностный порядок по каждой дискретной пространственной переменной. В предельных случаях, когда в уравнении (90) либо = 0, либо = 0 разностная схема (92)-(95) переходит в соответствующие одномерные разностные схемы.

Сходимость итерационного метода

В общем случае, способ расщепления нужно выбирать в зависимости от задачи. Тестовые расчёты задачи с периодической плотностью показали, что при расщеплении (+,+),(-,-) число шагов по времени для достижения максимальной точности на выбранной сетке примерно в 100 раз больше, чем при расщеплении (+,-),(-,+) . Очевидно, это связано с физикой задачи: в точном решении скорость направлена вдоль оси x, а скорость направлена против оси y. Поэтому, если в решаемой задаче есть выбранное направление, или заранее известны направления скорости, это нужно учитывать при выборе подходящего потокового расщепления.

Кроме того, в задачах, где присутствуют сильные ударные волны, необходимо увеличивать значение константы расщепления по направлению, нормальному к линии разрыва. Это позволяет получить неосциллирующее решение. Рис. 18. Численное решение задачи Case 17 по бикомпактной схеме. Давление показано цветом, плотность 30-ю изолиниями (0.53 - 1.98, шаг 0.05), скорость – векторами. 4

С увеличением размерности задачи d число уравнений в каждой ячейке растёт как 2d. При решении методом бегущего счёта, в каждой ячейке решается система размерности 2 с помощью прямого метода Гаусса, который имеет сложность С х (2d)3 = С х 2м. Уже при d = 3 это приводит к сильному росту вычислительных затрат по сравнению с одномерным случаем. Среди способов, которые позволяют избежать такого роста, можно упомянуть расщепление по направлениям (local one-dimensional scheme - LOD scheme), метод переменных направлений и методы приближенной факторизации линейных систем разностных уравнений. Большинство этих способов приводят к падению порядка точности до 2-го.

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

Рассмотрим факторизацию бикомпактной схемы для двумерного неоднородного уравнения переноса, для простоты выкладок будем считать коэффициент поглощения о постоян-ным: L2(u) = 0, L/2(u) = щ + аих + Ьиу + ои — Q = 0 (117) Уравнения бикомпактной схемы имеют вид (43a)-(43d): AxAyut + аАуАхи + ЬА\Ахи + Ау Ах ои = AyAxQ (118) АуАхщ + аАу Ах2и + ЪА\А\и + АуАхаи = AyAxQ А\Ахщ + аА\Ах1и + ЬАуАхи + А\Ахаи = A\AXQ A\Axut + аА\А2и + ЬАуАхи + А\Ахаи = A\AXQ с операторами: J\- ХЬ-іЛ-Л /О R -,/ —— U і /о Q A M/3,fc+l/2,7 — M/3,fc+l/2,7 AiMj+i/2,/3,7 = («i+1,/3,7 - Щрп)/1іх t\ -. 11 о і = { 11 а і \ л ) І П 1 ll /5, +1/2,7 V P, +J-,7— /3,fc,7 / / У Л2 (? +і/2,/з,7 = 6(tij__i,/5,7 — 2г/ +1/2й7 + jJf3J-y)/hx \У n ( ел—У і \ /12 A2W/3, +l/2,7 — 0( /3,/c+l,7 ZW/3,/c+l/2,7 "+" иРАу)/Пу (119) Рассмотрим базовую схему, когда система (119) интегрируется неявным методом Эйлера (каждую стадию DIRK метода можно рассматривать как шаг метода Эйлера): (120) (АХАУ(1 + аг) + таАуАх + тЪА\Ах)ип+ = АуАх(ип + rQn+ ) (АУАХ(1 + ат) + таАуАх2 + тЪА\Ах)ип+1 = АуАх(ип + rQn+l) (АУАХ(1 + ат) + таАуАх + гЬЛАж)мга+1 = Л ж(ига + rQn+1) (Л (1 + тт) + таЛ Л + т&ЛЛ?)ига+1 = Л Л (и11 + rQra+1) Введем следующие операторы: (121) B"(s) = (1 + 0.5ат)Аа + sA" B s) = (1 + 0.5 7т)Л" + sA2 C"(s) = 0.5(7 Аа + sA" C%(s) = 0.5(7Л" + sAg, a = x,y Операторы, стоящие в левых частях уравнений (120), допускают следующее представление: (122) АхАу{1 + ат) + таАуАх + тЪА\Ах = В\{Ът)Вх(ат) — т Cf(b)Ci(a) АуАх{1 + ат) + таАуА2 + тЪА\Ах = В\{Ът)В2(ат) — т СУ(Ь)С2 (а) А\АХ{1 + ат) + таА\Ах + тЬАуАх = В2{Ът)Вх(ат) — т С2{Ъ)Сх(а) A A l + ат) + таА\А2 + гЬА А = В2{Ът)В2(ат) — т С2{Ъ)С2(а) С учетом равенств (122) уравнения (120) можно записать таким образом: В\{Ът)Вх(ат)ип+ = АуАх(ип + TQU+ ) + т С\{Ъ)Сх(а) и га+ Ву(Ьт)Вх(ат)ип+1 = АуАх(ип + TQH+1) + т2С\{Ъ)Сх2{а)ип+1 Ву(Ьт)Вх(ат)ип+1 = АуАх(ип + rQra+1) + т2С ( (а)?/ (123) В2{Ът)В2(ат)ип+ = А\Ах(ип + TQU+ ) + г С2{Ъ)С2 (а) и га+1 Укажем на очевидные связи между операторами в левых и правых частях уравнений (123): " = Аа + тС"(а) 2 = А" + тС2(а), а = х,у Из структуры системы (123) видно, что её можно приближенно факторизовать при ма-лых т, отбросив последние члены в правых частях. Этот факт позволяет использовать для решения системы (123) следующий итерационный метод (далее вместо ип+1 используется й, вместо – ):

Уравнения (128) следуют из уравнений (126), если к ним применить оператор осреднения . В итоге, на каждом временном слое нужно решить 2( + ) одномерных задач. Каждая одномерная задача решается методом бегущего счета. Факторизованная схема имеет два преимущества перед исходной схемеой: во-первых, при решении одномерных задач обращаются матрицы размерности 22 вместо 2 2; во-вторых, одномерные задачи вдоль разных сеточных линий на каждой итерации можно решать независимо, поэтому алгоритм может быть эффективно распараллелен, в отличие от исходной схемы с более сложной зависимостью данных.