TVD-МОДИФИКАЦИЯ МЕТОДА ГОДУНОВА 3-ГО ПОРЯДКА Текст научной статьи по специальности «Математика»

МАТЕМАТИКА И МЕХАНИКА

www.volsu.ru

DOI: https://doi.org/10.15688/mpcm.jvolsu.2021.4.2

УДК 519.6:533.7 Дата поступления статьи: 17.09.2021

ББК 22.19 Дата принятия статьи: 18.10.2021

TVD-МОДИФИКАЦИЯ МЕТОДА ГОДУНОВА 3-ГО ПОРЯДКА

Евгений Иванович Васильев

Доктор физико-математических наук, профессор кафедры фундаментальной информатики и оптимального управления, Волгоградский государственный университет vasilev@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Татьяна Анатольевна Васильева

Кандидат физико-математических наук, доцент кафедры фундаментальной информатики и оптимального управления, Волгоградский государственный университет vasilevaTA@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Аннотация. В этой работе представлено исследование свойств монотонности новой модификации метода Годунова, имеющей 3-й порядок аппроксимации по пространству и времени. Разностная схема метода является полностью дискретной, то есть основана на совместной дискретизации уравнений по пространству и времени без использования стадий Рунге — Кутта. Доказано TVD-свойство разностной схемы применительно к линейному скалярному уравнению переноса при использовании адаптивного лимитера центральных разностей. Предложены новые модификации лимитеров, существенно улучшающие точность решения в окрестности разрывов и локальных экстремумов. Представлена экономичная версия метода для Ш-уравнений газовой динамики с квадратичной реконструкцией параметров в трациционных переменных. Продемонстрирована эффективная работа метода на некоторых стандартных тестах.

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

см о

см <

Н

га

m

&

■я Ч

к о га PQ

S Ы

m

■я Ч

к о га

PQ

©

Введение

В последние два десятилетия прошлого века широкое распространение получили нелинейные разностные методы 2-го порядка аппроксимации по пространству и времени, обладающие повышенной устойчивостью на разрывных решениях. С начала века активно развиваются нелинейные методы более высокого порядка аппроксимации по пространству и времени. В своем большинстве такие методы являются полудискретными [6; 8; 11; 12], реализующими два этапа аппроксимации дифференциальных уравнений в частных производных. На первом этапе с помощью кусочно-полиномиальной реконструкции строится дискретный разностный оператор для производных по пространству, в результате чего получается система обыкновенных дифференциальных по времени уравнений размером по количеству ячеек сетки. На втором этапе к этой системе применяются методы Рунге — Кутта 3-го или 4-го порядка по времени.

Отдельная дискретизация по пространству и времени упрощает построение метода, однако она существенно расширяет шаблон разностной схемы и требует дополнительных мер для обеспечения TVD-свойства, что приводит к увеличению количества стадий в методе Рунге — Кутта и порождает заметную диссипацию численного решения [6; 12]. Использование метода Рунге — Кутта для аппроксимации по времени впервые предложено в [14]. Позже эти же авторы в [15] пришли к выводу, что совместная аппроксимация является более эффективной.

В работе [4] предложен полностью дискретный метод 3-го порядка, использующий совместную дискретизацию уравнений по пространству и времени без использования стадий Рунге — Кутта. Этот метод является расширением метода Годунова 2-го порядка [1] за счет подключения двух отдельных блоков коррекции потоков. Первый блок повышает аппроксимацию для линейных систем уравнений, а второй устраняет погрешность, возникающую из-за нелинейности уравнений.

В первом разделе представлено компактное описание метода [4] без подробных обоснований. Во втором разделе сформулирована и доказана теорема о монотонности разностной схемы метода применительно к модельному уравнению переноса при использовании адаптивного TVD-лимитера центральных разностей [9]. В заключение раздела приведен способ построения новых модифицированных лимитеров, удовлетворяющих условию теоремы о монотонности.

На тестовых примерах показано, что новые лимитеры улучшают точность решения как в окрестности сильных и слабых разрывов, так и в окрестности локальных экстремумов решения.

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

1. Общий вид модификации схемы Годунова 3-го порядка

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

I + ^ = 0. (1)

от ОХ

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

по времени для г-й ячейки будет иметь вид

ипг+1 - иП I К+1) - / К-2) А

—-- +--2-— = 0, (2)

т п

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

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

д/

иг = -3(и)их , где 3(и) = — . (3)

ои

В методе 3-го порядка [4] выполняется квадратичная реконструкция, которая реализуется в виде двух этапов линейной реконструкции по схеме Горнера.

Для нелинейной /(и) квадратичная реконструкция потоков отличается от квадратичной реконструкции параметров и на величину 2-го порядка малости, которая в [4] учитывается введением дополнительных потоков N на последнем этапе коррекции

N = ^3^ - П23хих). (4)

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

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

1) Первый уровень левых и правых поправок во всех ячейках (индексы опущены)

й~ = и - Пих - т3(и)их, и+ = и + Пих - т3(и)их. (5)

6 3 6 3

2) Второй уровень левых и правых поправок во всех ячейках

= и - п^1- - Т^('и-)'^-, = и + - т^. (6)

Вычисление производных их на этих этапах осуществляется по трем точкам с использованием гладких ТУЭ-ограничителей [2; 10].

3) Предиктор, дающий решение и 2-го порядка

Щ - иП /(иг+1) - /(иг-1) ч ^

—^ +-2-~п-- = 0 , где иг+2 = Щь+^-ц) . (7)

4) Корректор с нелинейными потоками, дающий решение ип+1 3-го порядка

<+1 - U Ni+1 - Мг-i

+ —-Ч-1 = 0 , (8)

т h

где нелинейные потоки (4), например, по формулам

N = (J (иг) - J (ип))(иг - ип) ( J (иП+i) - J (иП))(иг+1 - ип) (9)

г+1 24 24 ' ()

Нелинейные потоки (4) можно вычислять с 1-ым порядком точности.

Формулы (5)-(9) описывают этапы одного шага для схемы Годунова 3-го порядка: предварительный шаг (7) с поправками первого (5) и второго уровня (6) плюс шаг коррекции (8) с нелинейными поправками третьего уровня (9).

2. TVD-свойство схемы 3-го порядка

В этом разделе проведем обоснование монотонности вышеописанной схемы 3-го порядка применительно к модельному уравнению переноса, свойства решений которого хорошо известны

ди ди „

— + а— = 0 , а = const > 0. (10)

dt дх

Обозначим через иг и иг приближенные значения и(х) в ячейке г в моменты t и t + т соответственно. В этих обозначениях схема Годунова 1-го порядка для (10) запишется в виде двухточечной схемы

т

иг = иг - у(иг - иг-i) , где v = а-, (11)

h

которая устойчива и монотонна при выполнении условий

0 < v < I. (12)

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

C(a, b) = (^М^И-»(М,1Ы>, »ТО, if ab> 0 , где ,(в, ь) = £±» (13)

\ 0, if ab < 0 ' V ' ; 2 v '

с адаптивным ограничивающим коэффициентом 0, зависящим от v

0 =--1-- ^ 1 < 0 < 2. (14)

max(v, 1 - v)

Заметим, что при 0 = 1 адаптивный лимитер (13) совпадает с обычным лимитером центральных разностей [9]. Функция C(a, b) удовлетворяет очевидным свойствам, которые будут использованы далее

|£(a, &)| < min(20|a|, 20|Ь|); C(a, Ъ) = £(Ь, a); (15)

С(а + |е|, b) ^ С(а, b); С( ka,kb) = кС(а, Ь). (16)

Для каждой ячейки г введем обозначение Ащ разностей против потока

Ащ = щ — щ_ь (17)

Для модельного уравнения переноса (10) отсутствуют нелинейные потоки (8) и решение задачи Римана %(v+, v~+l) = щ+. В результате схема (5)-(7) будет состоять из трех этапов 1

1) Un = и + -(1 - 2v)£(Au,, Аи+1), (18)

6

2) Vi = иг + -(1 - v)£(AU,, AUm), (19)

3) щ = n — vA Vi. (20)

Теорема 1. Разностная схема (18)-(20) с функцией-ограничителем (13)-(14) и с числом Куранта (12) является TVD-схемой.

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

С(а, Ь) — С(с, а) = 2ц0а, где |ц| < 1. (21)

При объединении (18)-(20) в одну формулу применим это равенство два раза

£(Ащ, Ащ+i) — С(Аи—1, Аиц) = 2п0Ащ, (22)

£(Ащ, AUi+i) — £(AU¿_i, Ащ) = 2ц0Ащ. (23)

Получим одно соотношение для схемы с двумя коэффициентами |ц| ^ 1 и |п| ^ 1

щ = щ — v( 1 + (1 — v)H0(1 + 3(1 — 2v)n0)) Ащ. (24)

Таким образом схема (18)-(20) формально приведена к двухточечной схеме (11), но с переменным коэффициентом v¿, зависящим от решения

Vi = v( 1 + (1 — v)H0(1 + 3(1 — 2v)n0)), |ц| < 1, |п| < 1. (25)

Двухточечная схема с переменным коэффициентом v¿ будет TVD-схемой, если выполнены условия Vi ^ 0 и 1 — Vi ^ 0 [7]. С учетом (12) проверка этих условий сводится к проверке неравенств для двух функций

Р(v, ц,П) = 1 + ц0(1 — V)(1 + 1 П0(1 — 2v)) > 0 , (26)

3

Q(v, ц,п) = 1 — H0v(1 + 1 n0(1 — 2v)) > 0 . (27)

3

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

11 - 2v|

Р=-' 1 \ ^ 0 1. (28)

max(v, 1 — v)

Далее оценим величины с крышкой в (23).

1 V

Айг = Ащ + -(1 - 2у)2цвАиг = (1 + кг)Аиг, где ^. (29)

63

Обозначим для краткости А = С(Ащ, Ащ+1), А = £.(Ащ, Айг+1). Пусть Аи^ > 0, тогда из условий (16)

А = С((1 + кг)Аиг, (1 + кг+1)Аиг+1) < С((1 + кА)Аиг, (1 + кА)Аиг+1) = (1 + кА)А,

А = С((1 + кг)Аиг, (1 + кг+1)Аиг+1) > С((1 - кА)Аиг, (1 - кА)Аиг+1) = (1 - кА)А,

где кА = тах(1кг1, 1) < р/3. Эти два неравенства объединим в одно |А- A| < кА|А|. Случай Аиг < 0 приводит к точно такой же оценке. Аналогичные оценки выполенны и для величин В = С(Аи1-1, Ащ) и В = ^Ащ^, Ащ). В результате с учетом первого свойства функции-ограничителя (15) будем иметь

И - A| ^ 3и ^ ^ АъЪ Ф - B| ^ ||Б| ^ М Аи^. (30)

В обозначениях А и В система (22)-(23) приводится к виду

А - В = 2цдАщ , А - В = 2ц0(1 + к^Ащ. (31)

Вычитаем в (31) первое уравнение из второго, переходим к модулям и учитываем ранее выведенные оценки (29) и (30). Получим

2|ц - л|0|А^| = |1 - А - В + В - 2ц0^Аи^ <

< - A| + |В - В| + 2|цQkiАui| < 2pв|Аui|, откуда непосредственно следует связь между коэффициентами ц и п:

|ц - п| < р. (32)

Таким образом, обоснование ТУЭ-свойства схемы сводится к необходимости доказать, что неравенства (26) и (27) для функций Р(у, ц,п) и Q(v, ц,п) выполнены в диапазоне аргументов 0 < V < 1, |ц| < 1, |п| < 1, |ц - п| ^ р.

Для фиксированного V диапазон представляет собой шестиугольник в плоскости (ц, п), изображенный на рисунке 1.

Заметим, что функции Р и Q связаны центральной симметрией

д(1 - V,-ц, -п) = Р(V, ц,п), Р(1 - V, -ц, -п) = Q(V, ц,п),

поэтому доказательство достаточно провести для 0 ^ V ^ 2.

По аргументам (ц,п) функции Р и Q являются гармоническими, то есть Рцц+Рпп = = 0 и Qцц + Qпп = 0, следовательно, их максимум в шестиугольнике достигается на границе шестиугольника. Задача сводится к проверке экстремумов функций одного переменного на шести сторонах шестиугольника. На сторонах АВ, АР, СИ и ИЕ (рис. 1) функции линейны, на сторонах ВС и РЕ функции квадратичны с точкой экстремума за пределами диапазона. То есть на всех сторонах шестиугольника функции либо монотонно возрастают либо монотонно убывают. В этом случае для проверки неравенств (26) и (27) достаточно вычислить и проверить значения функций в вершинах шестиугольника

А(ц = -1,п = -1), В(ц = -1,п = -1+ р), С(ц =1 - р,п = 1), Р(ц =1,п = 1), Е(ц =1,п = 1 - р), Р(ц = -1 + р,п = 1).

Рис. 1. Область возможного изменения параметров ц и ц функции устойчивости в (25)

Опуская несложные выкладки, выпишем значения функций Р и Q в вершинах, используя для краткости обозначение q = р/3:

Р (А) = q,

Р(В) = (1 - p)q,

Р(С) = 1 + (1 - Р)(1 + q),

Р (D) = 2 + q,

Р (Е) = 2 + (1 - p)q,

Р(F) = (4 - p)q,

Q(A) = 1 + (1 - р)(1 - q),

Q(B) = 1 + (1 - р)(1 - q(1 - р)),

Q(C) = (5 - р - p2)q,

Q(D) = (2+ p)q,

Q(E) = (2 + 2p - p2)q,

Q(F) = 1 + (1 - P)2(1 - q).

С учетом диапазона изменения 0 ^ р ^ 1 видим, что значения функций во всех вершинах неотрицательны. Следовательно, схема (18)—(20) при условиях (12)-(14) обладает ТУЭ-свойством. Теорема доказана.

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

ös ds

s(a,b) = s(b,a), — > 0, s(ka,kb) = ks(a,b), s(a,a) = a, — да da

(33)

a=b

Здесь последнее свойство требуется для обеспечения 3-го порядка аппроксимации разностной схемы. Такими свойствами обладает известная функция среднего гармонического и ряд других функций [2; 5]. Для построения новых аналогичных лимитеров можно

использовать следующее утверждение: если в(а,Ь) при положительных аргументах удовлетворяет условиям (33), то этим же условиям удовлетворяет и функция

ё(а,Ъ) = 8(а,Ъ)(1 + г(1 -д)2д), где & = ^ Ь\,, 0 < г < 3^3. (34)

шах(а, Ь)

Добавка в (34) приближает лимитер к известному лимитеру "зирегЬее" [13], сохраняя гладкость и положительность в указанном диапазоне коэффициента .

Далее в тестовых расчетах использовались лимитеры К,(а, Ь), С (а, Ь), М(а, Ь), которые при аргументах одного знака (аЬ > 0) имеют вид

К(а, Ь) = з!§п(а)шт(2|а|, в^, Щ), 2Щ), (35)

С(а, Ь) = з1§п(а)шт(20|а|, в^, Щ), 2Щ), (36)

М(а, Ь) = з1§п(а)шт(20|а|,§(|а|, |&|), 20|&|), (37)

где з(а, Ь) и 0 из формул (13)—(14), а ё(а, Ь) из формулы (34) с коэффициентом г = 3у/3.

Известным недостатком ТУЭ-схем является эффект срезания локальных экстремумов численного решения. На рисунке 2 представлены тестовые расчеты влияния вышеприведенных лимитеров на численное решение модельного уравнения переноса (10) с а = 1 с треугольным профилем точного решения после длительного перемещения.

Начальный треугольный профиль состоит из левого линейного участка шириной на 0,2 и высотой 5 и правого разрыва. Расчеты проводились на отрезке [0,1] с периодическими граничными условиями на сетке со 160-ю ячейками и числом Куранта СЕЬ= 0,6. Результаты расчетов для трех лимитеров изображены на момент времени I = 5, то есть на момент, когда начальный профиль проехал расстояние 5 длин расчетной области.

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

Заметим, что повышение качества решения не требут значительных затрат, так как лимитер М получается из лимитера С простыми арифметическими операциями.

3. Реализация схемы 3-го порядка для уравнений газовой динамики

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

дЬ дх

где

р р

™ = ру , {= ру2 + р

(е + р)у _

(38)

(39)

Здесь р(х, Ь), р(х, Ь), у(х, Ь) — плотность, давление и скорость газа зависят от координаты х вдоль трубы и времени ¿. Полная энергия единицы объема газа е = р(у - 1)-1 + + ру2/2 (для воздуха у = 1,4). Вектор массовых переменных ™ (плотность, импульс, энергия) используется для записи законов сохранения. Вектор традиционных переменных и = (р, р, у) обычно используется для ввода/вывода и анализа результатов.

Ца,Ь)

о.з

0.4

0.5

0.Б

0.7 0.3

0.7 0.3

Рис. 2. Сравнение численных решений с разными лимитерами (35)-(37) для модельного уравнения переноса (10) с а = 1 на решении с треугольным профилем. Сетка из 160 ячеек на отрезке [0,1] с периодическими граничными условиями, СРЬ = 0, 6. Результаты показаны на момент времен £ = 5

Систему уравнений (38)-(39) можно записать с производными по времени от и

£ + 4й =0- (40)

от ох

где

" р" V 0 р

u = р , а = 0 V YP

V 0 р-1 V

(41)

Пусть выбрана равномерная по х сетка с шагом к, причем направление роста номеров ячеек г совпадает с направлением координаты х. Схема Годунова для системы (38) записывается в следующем виде:

w?+1 - w?

+

f (ui+1 ) - f (Ui- 1 ) h

0,

(42)

где — вектор параметров в ячейке г в момент Ьп + т. Величины с полуцелыми ин-

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

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

и? + и

и I 1

г+1

+ 1 )-

*г+1

i+ 2 2 '+ 2' 2

(43)

Коэффициенты матрицы А с полуцелыми индексами вычисляются по средним значениям параметров в соседних ячейках. Функция б1§п(А) вычисляется через разложение

А = КАК-1 ^ з^п(Л) = Кэ^п(Л) ВТ1,

(44)

где К и А — матрица собственных векторов и диагональная матрица собственных чисел матрицы А, Б^п(Л) = ё1ад(в1дп(Л1),..., 8^п(Лт)). Разложение (44) для (41) выглядит следующим образом:

А

1

ср

1

11 0 с2 0 ср

1

V - с 0 0

0 V 0

0 0 V + с

1 с-2

2 Ь

2 рС-1

0

2

1 -с-2

0 к-2 2рс-1

(45)

где квадрат скорости звука с2 = ур/р.

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

т = V шт(

к

Ы + С,

),

(46)

где коэффициент V, называемый числом Куранта — Фридрихса — Леви (СЕЬ), должен быть меньше единицы (V < 1). В таком виде схема (42) имеет 1-й порядок аппроксимации.

Для записи схемы 3-го порядка введем обозначения для разностей

Ди = < - 1

Этап 1.

1 1 т

и? = < + Д( т1Е - 3 к Л) С(К-1Диг, Д-1Дит)

(47)

(48)

где матрицы К, А и К разложения (45) вычисляются на и™.

Этап 2.

1 1 т -

V- = < + д( -1Е -- к А) £(д-1ди-,н~1 ди-+1)

у+ = и? + Д( + 1Е - — А) £(Д-1Ди+, Д-1Ди++1) ,

1 т

2 к

(49)

(50)

здесь матрицы К, А и Я 1 вычисляются соответственно на и и на и Этап 3.

л - w •

+

f (и+1) - f(и_ 1)

0- где и+1 = v¿+l)

(51)

т к

2

с

Этап 4.

w™+1 - w т

. + 2 - N*-1

h

N„

i+2

24Й - Ши - u) - 24(Jm - Ji)(ui+i - u)•

Здесь матрица производных 7 от потоков системы (38) по переменным и имеет вид

V 0 р

J

1

2pw

v2

1V3 Y-1v Y-1р + 3 Р^2

(52)

(53)

(Р, P, v)

(54)

Тестирование изложенного метода в работе [4] подтвердило 3-й порядок точности по пространству и времени.

На рисунке 3 представлены результаты расчетов изложенным методом для стандартного теста [16] с лимитерами (35) и (37). Начальные условия на отрезке [0,1] представляют собой разрыв в точке х = 0,41 с исходными данными

(Р l,pl,vl) = (8, 10, 0), (рд ,pr,vr) = (1,1, 0), у = 1, 4.

0

Рис. 3. Сравнение численных решений с разными лимитерами (35) и (37) для теста Сода. Сетка из 100 ячеек на отрезке [0,1] , СРЬ=0,6. Результаты для плотности р(ж) показаны на момент времени £ = 0, 27

Расчеты проводились на сетке из 100 ячеек с числом Куранта в (46) V = 0,6. Результаты для плотности р(ж) показаны на момент времени £ = 0, 27. Для правого варианта на рисунке 3 лимитер М.(а,Ь) использовался только для первой и второй компоненты вектора характеристических разностей Д-1Ди в формулах (48)-(50). Тем

самым его действие проявлялось в волне разрежения и контактном разрыве. Видно, что при использовании М.(а,Ь) точность локализации контактного разрыва и границ волны разрежения существенно выше, причем признаки каких-либо численных флуктуаций отсутствуют. Разные лимитеры для разных характеристических полей успешно использовались в [5] для двухфазных течений. На целесообразность такого подхода впервые указано в [7].

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

Заключение

Предложено семейство гладких нелинейных адаптивных лимитеров, основанных на лимитере центральных разностей. Показано, что их использование существенно улучшает точность решения в окрестности разрывов первых производных и локальных экстремумов. Доказана теорема о том, что модификация метода Годунова 3-го порядка [4] с этими лимитерами обладает TVD-свойством. Представлена экономичная версия метода для уравнений газовой динамики с квадратичной реконструкцией параметров в трациционных переменных. Метод имеет 3-й порядок аппроксимации по пространству и времени и является развитием W-метода 2-го порядка [1] за счет подключения двух отдельных блоков коррекции. Предложенный подход построения разностных схем 3-го порядка может быть использован для неоднородных и двумерных гиперболических систем нелинейных уравнений.

СПИСОК ЛИТЕРАТУРЫ

1. Васильев, Е. И. W-модификация метода С.К. Годунова и ее применение для двумерных нестационарных течений запыленного газа / Е. И. Васильев // Журн. вычисл. матем. и матем. физики. — 1996. — Т. 36, № 1. — C. 122-135.

2. Васильев, Е. И. W-модификация метода Годунова и ее приложения в моделировании газодинамических течений с ударными волнами: дис. ... д-ра физ.-мат. наук / Васильев Евгений Иванович. — Волгоград, 1999. — 213 с.

3. Годунов, С. К. Разностный метод численного расчета разрывных решений уравнений газовой динамики / С. К. Годунов // Матем. сб. — 1959. — Т. 47 (89), № 3. — C. 271-306.

4. Метод Годунова 3-го порядка аппроксимации для уравнений газовой динамики / E. И. Васильев, Т. А. Васильева, Д. И. Колыбелкин, Б. Красовитов // Математ. физика и компьютер. моделирование. — 2019. — Т. 22, № 1. — C. 71-83. — DOI: https://doi.Org/10.15688/mpcm.jvolsu.2019.1.6.

5. Численное моделирование и экспериментальное исследование влияния синерезиса на распространение ударных волн в газожидкостной пене / E. И. Васильев, С. Ю. Митич-кин, В. Г. Тестов, Х. Хайбо // Журнал технической физики. — 1997. — Т. 67, № 11. — C. 1-9. — DOI: https://doi.Org/10.1134/1.1258855.

6. Bianco, F. High-order central schemes for hyperbolic systems of conservation laws / F. Bianco, G. Puppo, G. Russo // SIAM J. Sci. Comput. — 1999. — Vol. 21, № 1. — P. 294-322. — DOI: https://doi.org/10.1137/S1064827597324998.

7. Harten, A. High resolution schemes for hyperbolic conservation laws / A. Harten // J. Comput. Phys. — 1983. — Vol. 49. — P. 357-393.

8. Kurganov, A. A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations / A. Kurganov, D. Levy // SIAM J. Sci. Comput. — 2000. — Vol. 22, № 4. — P. 1461-1488. — DOI: https://doi.org/10.1137/S1064827599360236.

9. van Leer, B. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection / B. van Leer // J. Comput. Phys. — 1977. — Vol. 23. — P. 276-299.

10. van Leer, B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method / B. van Leer // J. Comput. Phys. — 1979. — Vol. 32. — P. 101-136.

11. Bryson, S. On the total variation of high-order semi-discrete central schemes for conservation laws / S. Bryson, D. Levy // J. of Scientific Computing. — 2006. — Vol. 27. — P. 163-175. — DOI: https://doi.org/10.1007/s10915-005-9046-8.

12. McCorquodale, P. A high-order finite-volume method for conservation laws on locally refined grids / P. McCorquodale, P. Colella // Comm. App. Math. and Comp. Sci. — 2011. — Vol. 6, № 1. — P. 1-25.

13. Roe, P. L. Characteristic-based schemes for the Euler equations / P. L. Roe // Annu. Rev. Fluid Mech. — 1986. — Vol. 18. — P. 337-365. — DOI: 10.1146/annurev.fl.18.010186.002005.

14. Shu, C.-W. Efficient implementation of essentially non-oscillatory shock-capturing schemes / C.-W. Shu, S. Osher // J. Comput. Phys. — 1988. — Vol. 77. — P. 439-471. — DOI: 10.1016/0021-9991(88)90177-5.

15. Qiu, J. Finite difference WENO schemes with Lax-Wendroff-type time discretizations / J. Qiu, C.-W. Shu // SIAM J. Sci. Comput. — 2003. — Vol. 24, № 6. — P. 2185-2198. — DOI: 10.1137/S1064827502412504.

16. Sod, G. A. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws / G. A. Sod // J. Comput. Phys. — 1978. — Vol. 27, № 1. — P. 1-31. — DOI: https://doi.org/10.1016/0021-9991(78)90023-2.

REFERENCES

1. Vasil'sev E.I. W-modification of Godunov's method and its application to two-dimensional non-stationary flows of a dusty gas. Comp. Maths Math. Phys., 1996, vol. 36, no. 1, pp. 122-135.

2. W-modifikatsiya metoda Godunova i ee prilozheniya v modelirovanii gazodinamicheskikh techeniy s udarnymi volnami: dis. ... d-ra fiz.-mat. nauk [W-Modification of Godunov's Method and Its Application to Modeling of Gas-Dynamic Flows with Shock Waves. Doctor Dissertation]. Volgograd, 1999. 213 p.

3. Godunov S.K. Raznostnyy metod chislennogo rascheta razryvnykh resheniy uravneniy gazovoy dinamiki [A Difference Method for the Numerical Calculation of Discontinuous Solutions of the Equations of Hydrodynamics]. Matem. sb. [Mat. Sbornik], 1959, vol. 47 (89), no. 3, pp. 271-306.

4. Vasilev E.I., Vasilyeva T.A., Kolybelkin D.I., Krasovitov B. Metod Godunova 3-go poryadka approksimatsii dlya uravneniy gazovoy dinamiki [The Third-Order Approximation Godunov Method for the Equations of Gas Dynamics]. Matemat. fizika i kompyuter. modelirovanie [Mathematical Physics and Computer Simulation], 2019, vol. 22, no. 1, pp. 71-83. DOI: https://doi.org/10.15688/mpcm.jvolsu.2019.L6.

5. Vasil'ev E.I., Mitichkin S.Yu., Testov V.G., Khaibo K. Numerical simulation and experimental research on the effect of syneresis on the propagation of shock waves in a gas-liquid foam. Technical Physics, 1997, vol. 42, no. 11, pp. 1-9. DOI: https://doi.org/10.1134/L1258855.

6. Bianco F., Puppo G., Russo G. High-Order Central Schemes for Hyperbolic Systems of Conservation Laws. SIAM J. Sci. Comput., 1999, vol. 21, no. 1, pp. 294-322. DOI: https://doi.org/10.1137/S1064827597324998.

7. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws. J. Comput. Phys., 1983, vol. 49, pp. 357-393.

8. Kurganov A., Levy D. A Third-Order Semidiscrete Central Scheme for Conservation Laws and Convection-Diffusion Equations. SIAM J. Sci. Comput., 2000, vol. 22, no. 4, pp. 1461-1488. DOI: https://doi.org/10.1137/S1064827599360236.

9. van Leer B. Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection. J. Comput. Phys., 1977, vol. 23, pp. 276-299.

10. van Leer B. Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov's Method. J. Comput. Phys., 1979, vol. 32, pp. 101-136.

11. Bryson S., Levy D. On the Total Variation of High-Order Semi-Discrete Central Schemes for Conservation Laws. J. of Scientific Computing, 2006, vol. 27, pp. 163-175. DOI: https://doi.org/10.1007/s10915-005-9046-8.

12. McCorquodale P., Colella P. A High-Order Finite-Volume Method for Conservation Laws on Locally Refined Grids. Comm. App. Math. and Comp. Sci, 2011, vol. 6, no. 1, pp. 1-25.

13. Roe P.L. Characteristic-Based Schemes for the Euler Equations. Annu. Rev. Fluid Mech., 1986, vol. 18, pp. 337-365. DOI: 10.1146/annurev.fl.18.010186.002005.

14. Shu C.-W., Osher S. Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes. J. Comput. Phys., 1988, vol. 77, pp. 439-471. DOI: 10.1016/0021-9991(88)90177-5.

15. Qiu J., Shu C.-W. Finite Difference WENO Schemes with Lax--Wendroff-Type Time Discretizations. SIAM J. Sci. Comput., 2003, vol. 24, no. 6, pp. 2185-2198. DOI: 10.1137/S1064827502412504.

16. Sod G.A. A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. J. Comput. Phys., 1978, vol. 27, no. 1, pp. 1-31. DOI: https://doi.org/10.1016/0021-9991(78)90023-2.

THE TVD MODIFICATION OF THIRD-ORDER GODUNOV METHOD

Evgeniy I. Vasilyev

Doctor of Physical and Mathematical Sciences, Professor, Department of Fundamental Informatics and Optimal Control, Volgograd State University vasilev@volsu.ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Tatyana A. Vasilyeva

Candidate of Physical and Mathematical Sciences, Associate Professor, Department of Fundamental Informatics and Optimal Control, Volgograd State University vasilevaTA@volsu.ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Abstract. A study of the monotonicity properties of a new modification of the Godunov method, which has a third order of approximation in space and time, is presented. The difference scheme of the method uses a simultaneous discretization of conservation laws in space and time without of Runge — Kutta stages. This method is a development of the second order Godunov method through the connection of two additional fluxes correction procedures. The first procedure increases the order of approximation for linear systems of equations. The second procedure uses the finite differences of the Jacobi matrix of the system of equations and eliminates the second-order error that occurs due to the nonlinearity of the equations. The TVD property of the difference scheme is

strictly proved in relation to the linear scalar transfer equation when using an generalized adaptive limiter of central differences. New modifications of limiters are proposed that significantly improve the accuracy of the solution in the vicinity of discontinuities and local extremes. New limiters are obtained from the known ones using a simple function of local smooth deformation. An economical version of the third-order Godunov method for gas dynamics equations with quadratic reconstruction of parameters in primary variables is presented. The use of primary variables for reconstruction significantly reduces the FPU time during calculations. The advantages of the third-order method in terms of the accuracy of the solution in the vicinity of contact discontinuities and the boundaries of the rarefaction wave are demonstrated on standard tests. The proposed approach to constructing third-order difference schemes can be used for inhomogeneous and two-dimensional hyperbolic systems of nonlinear equations.

Key words: Godunov's method, third order, TVD-properties, limiter functions, conservation laws.