CC BY
15
ЧЕБЫШЕВСКИЙ СБОРНИК Том 18 Выпуск 3
УДК 519.6 Б01 10.22405/2226-8383-2017-18-3-304-315
ОБ ОДНОЙ НЕЯВНОЙ РАЗНОСТНОЙ
СХЕМЕ ДЛЯ УРАВНЕНИЙ БАРОТРОПНОГО ГАЗА
Г. М. Кобельков1, А. Г. Соколов2
Аннотация
Для системы уравнений, описывающей течение идеального (вязкого) баротропного газа, предложена разностная схема, обеспечивающая положительность плотности. Доказано существование решения получающейся системы нелинейных уравнений при любых шагах сетки по времени и пространству. Предложен итерационный процесс для решения полученной системы на временном шаге.
Ключевые слова: Разностная схема, итерационный процесс, течение идеального (вязкого) баротропного газа.
Библиография: 7 названий.
IMPLICIT FINITE DIFFERENCE SCHEME FOR BAROTROPIC GAS EQUATIONS
G. M. Kobel'kov, A. G. Sokolov (Moscow)
Abstract
An implicit finite difference scheme approximated barotropic gas equations is proposed. This scheme ensures positivitv of density compared to previous methods. Existence of a solution to this scheme is proved for any time and space mesh-steps, an iterative method for solving the system of nonlinear equations on each time step is proposed.
Keywords: Difference scheme, iterative process, flow of an ideal (viscous) barotropic gas.
Bibliography: 7 titles.
1 Кобельков Георгий Михайлович, заведующий кафедрой вычислительная математика механико-математического факультета МГУ им. М. В. Ломоносова, ИВМ РАН, kobelkov@dodo.inm.ras.ru
2 Соколов Александр Германович, старший научный сотрудник кафедры вычислительная математика механико-математического факультета МГУ им. М. В. Ломоносова, shurunya@ mtu-net. ru
1. Введение
Построению разностных схем для уравнений газовой динамики и использованию этих схем для численного моделирования течений сжимаемого идеального (вязкого) газа посвящено огромное количество публикаций (см., например, [1, 2] и цитированную литературу). При этом какие-либо теоретические исследования практически отсутствуют. Так, например, естественный вопрос — а обеспечивает ли разностная схема положительность плотности (что следует из физики), остается открытым. Исключение составляет монография [3], где для уравнений вязкого сжимаемого газа в лагранжевых переменных в одномерном случае доказана неотрицательность плотности, и основанные на этом работы по построению разностных схем (см., например, [4] и последующие публикации), где аналогичный результат доказан для предложенной там разностной схемы,
В настоящей работе предложена неявная разностная схема, аппроксимирующая систему уравнений динамики баротропного газа, которая обеспечивает положительность плотности и сохранение баланса массы. Кроме этого, доказано существование решения разностной задачи при любых соотношениях шагов сетки по времени и пространству, а также выполнение энергетического неравенства. Предложен итерационный процесс для решения системы нелинейных уравнений, возникающих на временном шаге.
2. Постановка задачи и априорные оценки
Система уравнений, описывающая течений идеального баротропного сжимаемого газа имеет вид (см., например, [5])
др + д(ри) = о дх '
+ + £ = 0,Р = *Р< = еош ^ 1.
(1)
ЭЬдх дх
Для простоты в дальнейшем полагаем а =1, Уравнения (1) рассматриваются в цилиндре (^т = [0,1] х [0,Т] и дополняются краевыми и начальными условиями
и(0,Ь) = и(1,£) = 0, и(х, 0) = и0(х), р(х, 0) = ро(х) > 0. (2)
Приведем вывод энергетического тождества [6] для задачи (1), (2), который нам понадобится в дальнейшем при исследовании разностной схемы. Умножим первое уравнение (1) (уравнение неразрывности) екалярно Ь2 на — 2и2, а второе уравнение (1) — екалярно на и. Складывая результаты, по-
1„,2 2
лучим:
— 2(рг,и2) — 2((ри)х,и2) + ((ри), и) + ((ри2)х,и) + (рх,и) = 0.
Преобразуем левую часть следующим образом:
((ри)г,и) — 2((рг,и2) = 2[(рг,и2) + (р, (и2)^} = 11(р,и2).
Интегрируя последнее соотношение по ¿от 0 до ¿, имеем г
Г 1.1
дЛ,
Далее,
((ри)1,и) — )
2 [р(х,1)и2 (х,Ь) — ро(х)ио(х)2] йх.
{(ри2 )х,и) — (ри)х ,и2) = 1 ((ри)х ,и2) + (ри2,их)
(3)
1 (ри3)х йх = 0.
(4)
Наконец,
(р„и) = ^ (ри, (р-У-%) = (рг,р1-1) = ^ ((р^)<, 1) Интегрируя последнее соотношение по ¿, получаем
(рх,и)&
7 — 1
(р(х, £) — р2(х))в,х
Из соотношений (3)-(5) следует "энергетическое тождество'
2 р(1)и2(х,1)д>х +--р(х,£))д,х
Ч 1 — 1]
(5)
(6)
1 Рои20(1х +
1 — 1
p2(x)dx.
Примечание, В случае наличия вязких членов в уравнении движения (второе уравнение (1)), равенство (6) становится неравенством.
Наша цель — построить разностную схему, аппроксимирующую задачу (1), (2) и удовлетворяющую следующим условиям:
1, Если в начальный момент времени плотность р в узлах сетки положи-
тельна, то во все последующие моменты времени это свойство будет
выполняться.
1
4
1
1
1
1
1
1
1
2, Соотношение (6) при условии положительности плотности означает, что нормы
\Ш\\ыо,1], \\р1/2(^и(1)\\Ь2т равномерно ограничены по времени. Поэтому естественно потребовать от решения сеточной задачи, чтобы оно также было равномерно ограничено по времени,
3, Должен выполняться сеточный аналог закона сохранения массы,
4, Решение разностной схемы должно существовать.
3. Аппроксимация уравнения неразрывности
Введем на прямой равномерную сетку с шагом К = 1/М, хг = = 0,..., М] пусть т = Т/М — шаг по времени. Определим сеточные функции р™, 0 ^М — 1; и™, 0 ^ г ^ М, иП = и= 0, Будем считать, что за границами указанных индексов сеточные функции доопределены нулем. Как обычно, будем использовать обозначения
Ух = (Уг+1 — Уг)/к, Ьх = (Ьг — х)/к, VI = (ьп+1 — ьп)/т.
Аппроксимируем уравнение неразрывности (первое уравнение (1)) в г-и узле, используя следующие соображения:
а) Если щ ъ иг+1 имеют одинаковые знаки, то используем аппроксимацию "против потока", А именно, положим
(ри)х ~ ^ Рг+1Щ+1-Ргиг
РгУ-г + 1 Рг 1Щ 0 щ, и. + 1 > 0,
н , если иг,иг+1 < 0.
Нетрудно видеть, что в этом случае мы аппроксимируем член (ри)х с порядком О (К).
б) Если иг < 0, иг+1 > 0 (разрежение), то положим
/ ч , Рг(и+ — иг)
(ри) х = РхП + рих ~ ---.
К
Таким образом, аппроксимируется только второе слагаемое, а первое слагаемое полагаем равным нулю. Так как на концах отрезка длины К функция и принимает значения разных знаков, то на гладких решениях величина рхи на этом отрезке имеет порядок О (к), откуда следует, что аппроксимация (ри) х также имеет порядок О (К).
в) Если щ > 0, иг+1 < 0 (сжатие), то положим
/ ч , Рг+1 Пг+1 — РгПг (рг — Рг-^Ц (ри) х = Рхи + рих ~ ----+
КК
Первое слагаемое аппроксимирует (ри) х, в то время второе — аппроксимирует рхи. Однако, как уже отмечалось выше, на гладких решениях функция
скорости в этом случае на отрезке имеет порядок О (К). Таким образом, и в этом случае аппроксимация имеет первый порядок. Подытоживая рассмотренные выше случаи, имеем
PiU-i + l —О; ли-i ^ Г\
' h 1 1, если ui,ui+i > 0,
Pi+lUj + l-pjUj
( N л( N _ I н > если Ui, щ+\ < ° (ри)х ~ А{и)р = ^ {щ+н-щ) / п^ (7)
1 рг( г+н-г_1 , 6СЛИ Ч^ < 0, Щ+\ > 0,
Р1+1Щ+1~Р1- 1Щ , если щ > 0, щ+1 < 0
Рассмотрим двухслойную по времени неявную разностную схему, аппроксимирующую уравнение неразрывности. Будем использовать стандартные обозначения: рп = р, рп+1 = р и рг = где т — шаг по времени. Аппроксимируем первое уравнение (1) следующим образом
рг + А(и)р = 0. (8)
Запишем (8) в более удобном виде:
р^+1-Р? +
- pni++1:(-'Ui+1 + l'Ui+l\)+pni+1('Ui+1 + l'Ui+ll-'Ui + l'Ui\)-pni+1l('Ui + l'Ui\) _ 0 (9)
2h V, к '
0 ^ i ^ N - 1, П ^ 0.
Напомним, что сеточная функция р™+1 доопределена нулем при i _ — 1, i_ N.
Уравнение (9) неявно по р. Выпишем матрицу оператора А. Имеем (с точностью до множителя
А
( Ui + |«l| Ui — luil
—Ui — lull U2 + |«2|— Ui + luil 0 —U2 — lu2l
V
0
0
0 ...
U2 — lu2l 0
U3 + lu3l—U2 + lu2l U3 —|«з|
—u3 — 1чз1
0
AN-i,N—2
0
AN-2,N-i AN-i,N-i)
Матрица А обладает следующими свойствами:
- элементы на, диагонали неотрицательны,;
- элементы, вне диагонали неположительны,;
- сумма элементов по столбцу равна нулю.
Из (8) следует (I + тА(и))р = тр. Тогда из рассмотренного выше следует, что матрица I + тА(и) является М-матрицей, поэтому у обратной к ней матрицы все элементы будут положительны. Отсюда следует
0
0
Теорема 1, Пусть р0 > 0, ] = 0,..., М — I, тогда при любом п > 0 имеет место неравенство рп > 0, т.е. аппроксимация уравнения неразрывности (8) обеспечивает положительность плотности.
Из вида матрицы А следует
Теорема 2, Аппроксимация (8) уравнения неразрывности обеспечивает выполнение закона сохранения массы в сеточном, случае. А именно, имеет место равенство
М-1 М-1
^ Щ = ^ кР%п =1,...,^. (10)
3=0 3=0
Доказательство, Так как сумма элементов по столбцу матрицы А равна нулю, то отсюда следует утверждение теоремы.
4. Аппроксимация уравнения движения
Введем другую форму записи аппроксимации уравнения неразрывности, А именно, положим
{р} = [рг] — I Рх ^ёп(и), где [рг] = ■
заметим, что операторы {.}, [. ] зависят от сеточной функции и. Тогда по аналогии с (8) полностью неявная разностная схема для уравнения неразрывности примет вид
Рг + ({р}и)х = 0, ж = Иг, г = 0,...,М — 1. (11)
При аппроксимации уравнения движения будем использовать тот же принцип направленных разностей против потока, что и при аппроксимации уравнения неразрывности. Аппроксимируем исходную задачу полностью неявной разностной схемой:
(ри)г + ({ри}и)х + {р}(р1-1)х = 0, х = гк,г =1,...,М — 1. (12) 1 — 1
Неизвестными в уравнениях (11), (12) являются функции р и й, т.е. система уравнений (11), (12) содержит 2М — 1 уравнений и такое же количество неизвестных. Докажем равномерную ограниченность норм ||(рп)1/2ип||, Введем скалярное произведение (напомним, что функции рп и ип продолжены нулем на всю сеточную область) по формуле
М
(и, и) = ^^ киг
V
По аналогии е дифференциальным случаем умножим (11) екалярно на 2й2, а (12) — екалярно на тй. Вычитая из второго соотношения первое, получим
(13)
(14)
т((ри)г,и) - 2(рг,й2) + т(({рй}й)х,й) - 2(({р}й)х,й2)
+Ц ({ р}(р(-%,й) = 0.
Оценим выражения в (13), Имеем
т((ри\,и) - 2(Рг,й2) = (р,й2) - (р,ии) - 2(р,й2) + \(р,и2) = 1 [(р,й2) - (р,и2)] + 1 (р, (и - и)2).
Далее,
^({р}(р1-г)„й) = Ц((р1-г)„ {р}й)
= -^(Р-, ({Р}й)х) = Ц(р1-г,Рг) (15)
= ^(Р1, 1) - ¿л (Р1-, Р).
Оценим последний член в (15) с помощью неравенства Юнга с показателями
и 7:
1-1 1
Р1-г • Р ^ Р1 + 1Р1. Тогда из этого соотношения и (15) следует неравенство
^ ({Р}(Р1-)*,Ь) > [(Р1, 1) - (Р1, 1)] = ^ [(Р, 1) - (Р, 1)]. (16) Наконец, оценим последние два слагаемых в (13), А именно, докажем что
3 = (({рй}й)х,и) - 2 (({р}й)х,и2) ^ 0.
При оценке, для упрощения формул вместо р и и будем писать ржи. Воспользуемся формулой (9) для представления аппроксимаций ({р}и)х и ({ри}и)х. Имеем
-0.5(({р}и)х,и2) = 4н ^ [рг+г(\щ+1 \ -Щ+г)и2
- Рг (\иг+1 \ + Щ+г )и2 + рг(-\щ\ + Пг)и2 + р- (щ + \щ\)и2 (({ри}и)х,и) = 2н ^ \piniЫ+г - \щ+1 \)щ+1 щ
+Рг(иг+г + К+г\ + \щ\ - и^и2 - Рг-г(и-Щ + Щ-г\щ\)щ
(17)
Изменяя индексы в суммировании в формулах (17) и (18), получаем
3 = {({ри}и)х ,и) — 0.5(({р}и)х,и2) = 1 (1)рг, (19)
где I = ±(11 + 12 + 13 + Ь + Л + Л) и
11 = (1иг+11 + Щ+1)и2, 12 = (1Щ I— Щ)и.2, 1з = (Ы — Щ)и2-1,
14 = (1Щ+11 + Щ+1)и2+1, (20)
31 = 2(щ — 1иг1)игиг-1), 32 = — 2(1щ+11 + щ+1)щщ+1).
Покажем, что I > 0. Прежде всего заметим, что > 0, к = 1, • • • , 4. Рассмотрим восемь различных случаем. Имеем:
1) Пусть иг, иг+1, иг-1 > 0. В этом случае 31 = 0, 32 = —4и2+1иг,
11 = 2чц+1и2, ¡4 = 2и3+1. Тогда
I > 32 + 11 + Ь = 2щ+1(щ+1 — щ)2 > 0.
2) Пусть иг, иг+1, иг-1 ^ 0. В этом случае 32 = 0, 31 = 4иг-1и2,
12 = 2и21игI, 13 = 21иг1и2-1. Тогда
I >31 + 12 + 1а = 21иг1(иг-1 + Ы)2 > 0.
3) Пусть иг > 0, иг+1, иг-1 ^ 0, В этом случае 32 = 0, 31 = 0, т.е. I > 0,
4) Пусть иг ^ 0, иг+1, иг-1 > 0, В этом случае 32 > 0, 31 = 4иг-1и2, 12 = 2и21щI, 13 = = 21иг1и2-1. Тогда
31 + Ь + 13 = 21иг1(иг-1 + 1иг1)2 > 0.
5) Пусть иг-1 ^ 0, иг+1, иг > 0. В этом случае 31 = 0, 32 = 4и2+1иг, Д = 2щ+1и2, 14 = 2и3+1. Тогда
32 + Ь + 14 = 21т+11(т+1 + иг)2 > 0.
6) Пусть иг-1 > 0, иг, иг+1 ^ 0, В этом случае 32 = 0, 31 = 4иг-1и2 > 0, т.е. I > 0.
7) Пусть иг+1 ^ 0, иг, иг-1 > 0, В этом случае 32 = 0, 31 = 0, т.е. I > 0.
8) Пусть иг+1 > 0, иг, иг-1 ^ 0, В этом случае 32 = —4иги2-1 > 0, 31 = 4и—и2, 12 = 2и21иг1, 1з = 21щ1и2—1. Тогда
1>31 + 12 + 1з = 21иг1(иг-1 + 1иг1)2 > 0.
Итак, мы установили, что 1(г) > 0. Тогда из (19) и того факта, что р > 0, следует неравенство 3 > 0.
Из соотношения (13) и неравенств 3 > 0, (14), (16), таким образом, следует
Теорема 3, Для разностной схемы справедливо равенство
п— 1
1 („п (ип\2^ и+1 г+1 _ „1
{рп, (ип)2) + 1 £ {рг+1, (иг+1 — иг)2) + -—(рп, 1)
г=0 7 1 (20)
^ 2 (р0, (П0)2) + ф-1 (р0, 1).
5. Существование решения разностной схемы
Докажем, что решение (11), (12) при любых шагах сетки ги^ всегда существует, Действительно, введем новую переменную V = ри. Тогда уравнение (12) примет вид
ьг + (Щи/р)х + {р}{р-—1)-х = 0. (21)
Для доказательства разрешимости задачи (11), (21) воспользуемся теоремой Лере-Шаудера, Запишем (11), (21) в виде
Р = Р — т({^}й)х = ^ (p,v),
V = V — г [({Ъ}Ъ/р)х + {р}{р1—1)х] = Ш$).
(22)
Покажем, что все возможные решения задачи
р = \Р1(р,у),
ь = \Р2(р,ь), А е (0,1),
(23)
лежат в некотором шаре. Действительно, элементарными преобразованиями (23) сводится к виду
р + рг + ({й} р)х = 0,
(24)
^ ри + (ри^ + ({р'й}й)х + р(р> = 0
Умножая первое уравнение (24) екалярно на — 0.5й2, а второе — на и, после сложения получим
(р,й2) + 1 (р, (й)2) + -—-1 (р-, 1)
2у I -—1I
^ 1 {р,и2) + (Р-, 1) ^ 1 {р0, (и0)2) + -—1 (р0, 1). Кроме этого, в силу (10) и теоремы 1 имеем
М 1 М 1
Р > 0, кРг = Пр0
пр 0
=0 =0
Далее,
\axpi ^ У^ Рг.
М -1
тах '
г=о
Тогда используя (25) и приведенные выше соотношения, получаем
1 Р\\ 17Л + 2 И2 = IIР11I, + Щ рй\\2
ф-1 IIР11 + 2к IIР0 КЛ р,^2)
1 тах{1, Цр0Ць^} \^ Цр|| + 2 (р,й2)
(26)
^ 1 тах{1, Цр0Ць^} [1 {р0, (й0)2) + (р0,1)
Неравенство (26) означает равномерную по Л € (0,1) ограниченность норм решений операторного уравнения (23), Тогда из теоремы Лере-Шаудера [7] следует
Теорема 4, Решение задачи (11), (21) при любых тик существует.
Отметим, что по доказанному ранее, р > 0, поэтому й = и/р также существует.
Для нахождения решения рп+1, йп+1 задачи (11), (12) может быть использован следующий итерационный процесс
рк + 1 -рп,
+ ({рк+1 }йк )х = 0,
рк+1ик+г1 -р"и" + ({рк+1 йк+1 }йк+1 )х (27)
+ ^ {рк+1}((рк+1)1-1)-х = 0, йк I к=0 = йп;
здесь к = 0,1,... — помер итерации. Для решения первого уравнения (27) в одномерном случае можно использовать метод прогонки, а в многомерном — метод простой итерации. Для решения второго уравнения (27) используем метод простой итерации. Можно доказать, что при достаточно малом шаге по времени эти процессы сходятся.
6. Заключение
Таким образом для системы уравнений, описывающей течение идеального (вязкого) баротропного газа, предложена разностная схема, обеспечивающая положительность плотности. Доказано существование решения получающейся системы нелинейных уравнений при любых шагах сетки по времени и пространству. Предложен итерационный процесс для решения полученной системы на временном шаге.
Работа выполнена при частичной поддержке РФФИ (проект 17-01-838) и Программы ОМН РАН "Современные проблемы теоретической математики" ("Оптимальные методы решения задач математической физики"),
СПИСОК ЦИТИРОВАННОЙ ЛИТЕРАТУРЫ
1, Белоцерковский О.М.. Васильев М.О.. Ведерников А,Б,, Дымников В,П., Замышляев Б,В., Борис Юрьевич Крысапов Б.Ю., Николай Вадимович Ковшов 11.В.. Вячеслав Евгеньевич Куницын В.К.. Евгений Александрович Молоков К.Л.. Репин Л.К).. Сидоренкова 11.Л.. Ступнцкпй Е.Л., Холодов Я.А., Холодов А,С, О численном моделировании некоторых задач взаимодействия литосферы, гидросферы и атмосферы Земли, В кн.: Фрагменты истории и достижения ПАП РАН, 1986-2011, ПАП РАН, ООО I III "Полет Джонатана" , 2011, с, 14 - 71,
2, Попов И,В., Фрязинов И,В, Метод адаптивной искусственной вязкости численного решения уравнений газовой динамики, «М,: КРАС АНД, 2014», 2015.
3, Антонцев ('.II.. Кажихов А,В., Монахов В.Н, Краевые задачи механики неоднородных жидкостей, Новосибирск, Наука, 1983,
4, Амосов А,А,, А, А, Злотник А,А,, Разностные схемы второго порядка точности для уравнений одномерного движения вязкого газа. Ж, вычисл, матем, и матем, физ,, 1987, т, 27, №7, с,1032-1049
5, Куликовский А,Г., Погорелов 11.B.. Семенов А.К). Математические вопросы численного решения гиперболических систем уравнений, М,: Наука, 2001.
6, Lions P.-L, Mathematical Topics in Fluid Mechanics, V.2. Compressible Models, Oxford Lecture Series in Mathematics and its Applications 10, Oxford: Clarendon Press, 1998,
7, Ладыженская О, А, Математические вопросы динамики вязкой несжимаемой жидкости, М,: Наука, 1970,
REFERENCES
1, Belotserkovskii О,M,, VasiPev M,О,, Vedernikov A,В., Dymnikov V.P., Zamyshlvaev B.V., Krvsanov B.Yu.. Kovshov N.V., Kunitsvn Y.K.. Molokov E.A., Eepin A.Yu,, Sidorenkova N.A., Stupitskii E.L., Kholodov Ya.A,, Kholodov A,S, Numerical simulation of some problems of interaction of Earth's lithosphere, hydrosphere and atmosphere. In book: [Hystory fragments and achievements of the Design Automatization Institute of Russian Academy of Sciences, 1986-2011], Publishing House " Johnatan
flight" , Институт автоматизации проектирования РАН, ООО ИЦ " Полет Джонатана" , 2011, с, 14-71,
2. Popov I.V., Frvazinov I.V, Метод адаптивной искусственной вязкости численного решения уравнений газовой динамики Metod adaptivnoy iskusstvennoy vyazko-sti [Adaptive artificial viscosity 'method], Moscow: Krasand, 2014,
3. Antontsev S.N., Kazhikhov A.V., Monakhov V.N, Краевые задачи механики неоднородных жидкостей Kraevye zadachi mekhaniki neodnorodnykh zhidkostey [ Boundary value problems of inhomogenious flows]. Новосибирск, Наука, 1983.
4. Amosov A.A., Zlotnik А.А., Разностные схемы второго порядка точности для уравнений одномерного движения вязкого газа. Raznostnye skhemy vtorogo poryadka tochnosti dlya uravneniy odnomernogo dvizheniya gaza [Second order accuracy finite difference schemes for ID viscous gas flow] Ж. вычисл. матем, и матем. физ,, 1987, т. 27, №7, с. 1032-1049
5. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. Математические вопросы, численного решения гиперболических систем, уравнений. Matematicheskie voprosy chislennogo resheniya giperbolicheskikh sistem [Mathematical problems of numerical solution of systems of hyperbolic equations] M,: Наука, 2001.
6. Lions P.-L. Mathematical Topics in Fluid Mechanics. V.2. Compressible Models. Oxford Lecture Series in Mathematics and its Applications 10. Oxford: Clarendon Press, 1998.
7. Ladvzhenskava O.A. Математические вопросы, динамики вязкой несжимаемой жидкости Matematicheskie voprosy dinamiki vyazkoy neszhimaemoy zhidkosti [The Mathematical Theory of Viscous Incompressible Flow], M.: Наука, 1970.
получено 22.05.2017 принято в печать 14.09.2017
