Об одной разностной схеме для уравнений газовой динамики Текст научной статьи по специальности «Математика»

12. Феллер В. Введение в теорию вероятностей и ее приложения. Т. 2. М.: Мир, 1984.

13. Esary J., Prochan F., Walkup D. Association of random variables with applications // Ann. Math. Stat. 1967. 38, N 5. 1466-1474.

14. Булинский А.В., Шашкин А.П. Предельные теоремы для ассоциированных случайных полей и родственных систем. М.: Физматлит, 2008.

15. Piterbarg V.I. Discrete and continuous time extremes of Gaussian processes // Extremes. 2004. 7, N 2. 161-177.

Поступила в редакцию 17.05.2017

УДК 519.6

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

А. В. Звягин1, Г. М. Кобельков2, М. А. Ложников3

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

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

A conservative difference scheme with linear dependence of the pressure on the density of gas is proposed for gas dynamics equations. The scheme allows us to simulate 1-D flows inside a cylindrical domain with time-variable cross-sections and guarantee the positive sign of the density function.

Key words: gas dynamics, conservative difference schemes, variable boundary, positive density.

1. Введение. Одной из актуальных прикладных задач гидро- и газовой динамики является расчет течения при закрытии клапанов и заслонок в трубопроводах [1, 2]. Для одномерного моделирования таких течений обычно используется метод характеристик. Клапан моделируется как место, где параметры течения (давление и скорость) разрывны, т.е. имеют разные величины с двух сторон от клапана. При этом система уравнений механики сплошной среды становится незамкнутой. Для того чтобы ее замкнуть, добавляются два эмпирических соотношения, связывающие скорость и давление перед клапаном со скоростью и давлением за ним. Эти соотношения зависят от технических характеристик конкретного клапана. В настоящей работе предлагается рассматривать клапан как часть трубы с переменным во времени сечением. Данная схема позволяет моделировать закрытие клапана в одномерном течении.

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

2. Постановка задачи. Рассмотрим систему уравнений, описывающую движение вязкого ба-ротропного сжимаемого газа (жидкости) в прямой трубе переменного радиуса:

д(Ар) + д(Ари) = Q

dt dx

2, я_ (1)

dt ox ox

1 Звягин Александр Васильевич — доктор физ.-мат. наук, проф. каф. газовой и волновой динамики мех.-мат. ф-та МГУ, e-mail: [email protected].

2 Кобельков Георгий Михайлович — доктор физ.-мат. наук, проф., зав. каф. вычислительной математики мех.-мат. ф-та МГУ, вед. науч. сотр. ИВМ РАН, e-mail: [email protected].

3 Ложников Михаил Андреевич — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: [email protected].

где р, и и р — плотность, скорость и давление жидкости соответственно, А = А(х, ¿) — площадь поперечного сечения трубы, а Ь = Ь(х, ¿) — периметр поперечного сечения трубы. Связь между давлением и плотностью положим линейной р = кр, где к — некоторая положительная константа. Сила трения во втором уравнении системы (1) моделируется экспериментальной формулой Аи|и|, где Л — положительная константа, зависящая только от числа Рейнольдса. Уравнения будем рассматривать в цилиндре Qт = [ХЬХ2] х [0,Т].

3. Априорные оценки. Дополним систему (1) краевыми

и(Х1) = и(Хг) = 0 (2)

и начальными условиями

р(х, 0) = po(x), 0) = u0(x). (3)

Не уменьшая общности, упростим запись, считая, что к = 1 и А = 1. Умножив первое уравнение

системы (1) скалярно в L2 на — \и2, а второе на и и сложив, получим

+ ((Ари)ь,и) - ^(Ари)х,и2^ + ((Ари2)х,и)+(Арх,и) + (и\и\Ь,и) =0. (4) Элементарные преобразования дадут

- ^ [{Ар\, и2^ + ((.Ари\,и) = ~ (Ар, и2) , -\(^Ари)х,и2^ +((Ари2)х,и) =0. Преобразуя пятый член в левой части (4) в предположении, что р > 0, будем иметь

(5)

(Арх,и) = (^px,Apu^j = ((Inр)х,Ари) = (Inp,(Ap)t) = ((Aplnp)t - Ap(lnp)t,l) =

.P" ' J — —— - - -rv-r,.,-, (6)

= ((Apln p)t - Apt, 1) = ((Apln p)t - (Ap)t + pAt, 1) = ((Ap(ln p - 1))t + pA(ln A)t, 1).

Таким образом, равенство (4) с учетом (5) и (6) примет вид 1 d

__ (Ар, и2) + ({Ар{In р - l))t, 1) + (рА{In A)t, 1) + [\u\L, и2) = 0. (7)

Если функция (ln A)t непрерывна на отрезке [Xi,X2] при любом фиксированном t ^ 0, то предпоследнее слагаемое в (7) можно оценить с помощью неравенства

(pA(ln A)t, 1) ^ (pA, 1) min (ln A)t = (poAo, 1) min (ln A)t. (8)

[Xi ,X2 ] [Xi ,X2]

Используя соотношение (8), формально проинтегрируем (7) по времени на отрезке [0,Т]. Получим "энергетическое" неравенство

T X2 X2 т

1

2 (^Ар, и2^ + (Ар(\пр - 1), 1) + J j u2\u\Ldxdt + j p0A0dx J ^min ^(ln A)tdt ^

o Xl Xl

2

T

< i (Aopo,u^ + (A0po(ln po - 1), 1). (9)

Оценим величину ш1п(1п А)^ в некотором частном случае. Пусть функция А удовлетворяет условию

А* ^ -в при х € [Х1,Хг], í € [0,Т],

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

[Xi,X2]

A(x,t) ^ A0,min — et на отрезке [X^X2] и

„ At в в

(In A)t = 4>->-

A A Ao,min — et' В этом случае последнее слагаемое левой части (9) можно оценить с помощью неравенства

X2 T X2 T X2

/ poAodx / min (In Ä)tdt~^ / poAodx / —--— dt = In 0,I"in- / poAodx. (10)

[Xi,X2] Ao ,min et A0,min J

Xi o Xi o Xi

Выражение (10) позволяет в случае линейно убывающей по времени функции A(x,t) оценить последний член левой части неравенства (9) через минимальную площадь поперечного сечения канала в момент времени T.

4. Разностная аппроксимация. Введем сетку по пространству с шагом h = (X2 — Xi)/N, а именно Qp = [xi = X1 + ih, i = 0,..., N — 1}, Qu = [xi = X1 + ih, i = 0,..., N}. Положим шаг по времени т = T/M. Пусть в момент времени пт сеточная функция pn определена на Qp, а un определена на Qu. Будем считать, что эти функции продолжены нулем на всю сеточную прямую. Разностную аппроксимацию системы (1)—(3) будем строить таким образом, чтобы она обладала следующими свойствами.

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

N—1 N — 1

Е

i=0 i=0

hAfpf = £ hAn+1Pn+1,

где АпП = А(гЬ, пт).

2. Плотность рп должна быть строго положительна.

3. Сеточные аппроксимации (Ари)х и (Ари2)х должны быть сопряжены, т.е. должен выполняться сеточный аналог равенства

((Ари2)х, и) - i ({Ари)х, u2^j = 0.

¡XI М) — I

4. Решение разностной схемы должно существовать.

5. Из соотношения (9) следует ограниченность нормы ||(Ар)1/2и||^2 на любом конечном временном интервале. Решение разностной схемы должно обладать аналогичным свойством.

Аппроксимация уравнения неразрывности. Будем аппроксимировать уравнение неразрывности аналогично [4]:

если щ и щ+1 имеют одинаковые знаки, то используем аппроксимацию "против потока":

{piAi+lUi+l - рг-1Ащ --- при щ, щ+1 > 0;

л ^ л

Pi+lAi+lUi+l - piAiUi --- при щ, щ+1 < 0;

если и,, > 0, а щ+1 < 0, то

/ , ч л . / л \ Р,+1А,+1Щ+1 - р,Ащ (р, - р—)Ащ

(рАи)х = рхАи + р{Аи)х &----Ь

hh

если ui < 0, а ui+1 > 0, то

срАч)х = РхАч + р{Аи)х =

Данная аппроксимация имеет порядок О (К). Как обычно, обозначим й = ип+1, и = ип, их = щ, Щ = у'п+1~'и'п ; в безындексной записи и = щ, = и рассмотрим полностью неявную раз-

ностную схему для задачи (1)-(3). Для удобства введем сеточный оператор {■} следующим образом:

)рг-1, если щ > 0;

{Р}г = \ / п

I рi, если щ < 0.

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

(Ар); + ({р}Аи)х = 0. (11)

Точно так же как в [4], доказывается, что аппроксимация (11) обеспечивает выполнение сеточного закона сохранения массы, т.е.

N—1 N—1

£ р0 = £ МПрП, (12)

i=0 i=0 аО \ П „ Лп

и при выполнении условий р" > 0 и А? > 0 для любых I = 0, N — 1 ип ) 0 имеет место неравенство р™ > 0 при любых г = 0Д-1ип>0.

Аппроксимация уравнения движения. Второе уравнение системы (1) будем аппроксимировать следующим образом:

(Ари^ + -\{р}Ай(щ_1)+й)) +и\А{р} 1п—^--Ь \u\u\L = 0. (13)

2\ /х н р(—1)

Докажем, что норма ||(Ар)1/2и||£2)^ равномерно ограничена по времени на любом конечном времен-

"2 '

ном интервале. В самом деле, умножим уравнение (11) на — ¿тй2, а уравнение (13) на тй и сложим.

В результате получим равенство

-т\ ^ , й2) + тк (^Мр) 1п -у—> + тЛ (ü\ü\L, ü) = 0. (14)

Нетрудно показать, что

-\r^(Ap)t,ü2^ +t((Apu)t,ü) = ^t[aP,u2^ + ^(aP,(ü-u)2^, (15)

а сеточные аппроксимации (Apu)x и (Apu2)x кососимметричны, т.е.

^^({pjiü^.!) "4(({РМй)ж>й2) =°- (16)

С помощью соотношения ln x ^ x — 1 упростим предпоследнее слагаемое в (14):

т (^Ä{p}lnj£—,üj =r((lnp)ä,i{p}ü) =T((Ap)t\np,l) = т ((Aplnp)t - Ap(lnp)t,l) =

= т ((Apln p)t,l) - ^Apln^lj > т ((Aplnp)t,l) - (A(p - p),l) = т ((Aplnp)t - (Ap)t + pAt,l) >

At At

> т ((Ap( ln p - l))t, 1) + т(Ар, 1) min -f = т ((Ap( ln p - l))t, 1) + т(А0р 0,1) min -f. (17)

Пр A пр A

Просуммировав соотношение (14) по всем шагам по времени и воспользовавшись соотношениями (15)—(17), получим энергетическое неравенство. Таким образом, доказана следующая теорема.

Теорема. Решение разностной задачи (11), (13) удовлетворяет неравенству

1 1 n—1 , 2\

- (Апрп, (ип)2) + к(Апрп(Ырп - 1), 1) + - ^ ( Акрк, (uk+1 +

U—п V '

n—1 лк+1 _ Ак n

+К(А0ро, 1) £ Afc+1 + r\J2(\uk\Lk, (ик)2) < к=0 р к=1

< \ (а°р°, {и0)2) + к(а°р°(lnp° - 1), 1).

Точно так же как и в дифференциальном случае, оценим величину ^fc=o min А Ак+^ ПРИ нек0~ торых дополнительных ограничениях. Обозначим minnp A0 через A0,min. Пусть выполнены следующие соотношения:

^4o,min А — —--

т < -—--- и - ^ — в при к = 0, п — 1, (18)

в(п + 1) т v 7

где в — некоторая положительная константа. Тогда, используя оценку Ak ^ A0,min — твк, получаем

h Ak+l " hAk+l" h - т9(-к+^'

Отметим, что при условиях (18) выполнено неравенство -j-<1, А; = 0, п — 1. С помощью

неравенства у ^ оценим правую часть (19). Имеем

У^_^тв_ > у^ _ 1д ^0,min - тв{к + 1) _ Л,min ~ тв{п + 1)

Л),min - TÖ{k + 1) " ^ П Ao.min " тв(к + 1) - тв ~ П A0,min " тв

Данное соотношение позволяет оценить предпоследнее слагаемое левой части энергетического неравенства в случае линейно убывающей по времени функции A(x, t) с помощью минимальной площади поперечного сечения канала на n-м шаге по времени.

5. Существование решения. Пусть выполнено условие

(A0P0,1) — (An, 1) > 0 при любом п > 0. (20)

Аналогично [4] докажем, что система уравнений (11), (13) с краевыми условиями (2) и начальными условиями (3) имеет решение при любых т и h. Сделаем замену r = Ap и v = Apu. Тогда система примет вид

rt + ({V/A}AV/V^ = 0,

1 \ 1 ^ V/A ~ (21)

vt + - [{r/A}Av (i(_i)/f(_i) +v/f) /г) + n-A{r/A) In-'—--h Xv/r\v/r\L = 0.

24 Jx h f(—1)/A(—1)

Таким образом, система (21) представима в виде операторного уравнения

V = F1(f,V), V = F2(r,V).

Для доказательства существования решения воспользуемся теоремой Лере-Шаудера, а именно докажем, что решение системы

г = eF1(r,v), r = eF2(f,v) (22)

равномерно ограничено по в при в € (0,1). Перепишем систему (22) в виде 1в

Ар + (Ap)t + ({р}^«) = 0,

—^-Äpü + (Apu)t + - ( {p}Äü(üi_i\ + и)) + n\-Ä{p} In -J?---h Aü|ü|L = 0.

тв 2V Jx h p(—1)

Умножим первое уравнение на — \тй2, а второе на тй и сложим. Так как функции А(х,1) и Ь(х,1)

2

положительны, а pk > 0 при условии p0 > 0, то

1 — в f Va а2\ 1 — в V A — A 1,- л2

(Ар, и2) + к——(Ар,In р) + к(А0ро, 1) min —---Ь - (Ар, й2] + к(Ар(\п р - 1), 1) ^

V/в Пр A 2 V /

< ^ (Ар, u2) + к(Ар(1п р - 1), 1) <

1 f \ n—1 Ak+1 — Ak

^ -(A0p0,ulj + к(Лоро(1про - 1), 1) - и(Аоро, 1) 2_,™in—~Шл—• (23)

к=0 р

Оценим второе слагаемое левой части соотношения (23), учитывая условие (20). Имеем

(Ар, In р) ^ [Ар = (Ар - А, l) = (Аоро, 1) - (А, l) > 0. (24)

Для последнего слагаемого левой части (23) справедливо соотношение

(ip(In Р - 1), l) ^ (Ар (l - 0 - U, l) = - (Л l) . (25)

Используя оценки (24) и (25), а также сеточный закон сохранения массы (12), получим \\v\\l2th = (-Apü,Äpü) ^ rnaxАр ■ (Ap,ü2^ ^ ^\\Aopo\\blth (Ap,ü2) ^

< ^\\Aopo\\bhh

I ~ n Ak+l - Ak'

(A0po, (uo)2) + 2n (A0po(lnp0 - 1), 1) + (А, 1) - (Аоро, 1) min

к=0 Пр

Пр Ak+1

Величина Г ограничена в норме за счет свойства консервативности (12) разностной схемы. Таким образом, выполнены условия теоремы Лере-Шаудера, значит, задача имеет решение при любых т и Н.

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

дом) = (1)'

Ь(х,г) =2пЯ(х,г), (26)

А(х,г) = п(п(х,г))2,

где параметры ко и Т^ соответствуют скорости закрытия клапана и полному времени закрытия клапана, а через I обозначена характерная длина клапана. Формулы (26) верны на отрезке х € [-1,1], вне этого отрезка положим Ь(х, Ь) = 2п и А(х, Ь) = п. Будем считать, что граничные узлы и их соседи лежат вне клапана. Тогда сеточные граничные условия зададим следующим образом:

(ри)х (ХьЬ) = (ри)х (Х2,Ь) = 0,

р (X ,Ь) = п + ^ (и (ХьЬ)) ,р (Х2,Ь) = г 2, ( )

где Г1 > Г2 — некоторые постоянные, а функцию ^(и), соответствующую паспортной характеристике насоса, определим так:

^(и) = Со - С1и2, где Со, С > 0. (28)

Начальные условия и величину Г2 можно определить, решив стационарную задачу Коши, т.е. поло-

лучае

^(рц) = 0

йх '

*&й + Л+2\и\и\=0, (29)

р(0) = Г1, и(0) = и0.

жив = 0 и = 0. В этом случае уравнения примут вид

Итак, будем решать систему (11), (13), (26) с граничными условиями, определяемыми формулами (27), (28), и начальными условиями, определяемыми решением задачи Коши (29). Для нахождения решения задачи (рп+1,ип+1) на каждом следующем шаге по времени можно использовать внутренний итерационный процесс по к:

Ап+1„к+1 _ Апр

Л Р Л Р , / г „АН-1.1 лп+1 к

+ ({рк+1}Ап+1ик) =0,

+ -2({рк+1}Ап+1ик+1 (г,к11} +

Ап+1 „к+1,,к+1 _ лппп„.п 1

Л Р а Л Р а I _ [ {рк+1\АП+1Ук+1 (цк ' »-к

1 рк+1

+ К-г ик |Ьп+1 = 0,

Н р(—1)

р0к+1 = Г1 + ^(ик), р^1 = Г2,

рк+1ик+1 _ .0+1.0+1 = 0 рк+1ик+1 _ рк+1 и к+1 = 0

р0 и0 р1 и1 = 0, ^ —1 =

На рисунке приведены результаты расчетов в области [-10,10] на сетке из 500 узлов, шаг по времени равен 10—3. Графики представляют собой зависимости плотности от времени перед клапаном (в 225-м узле, кривая 1) и сразу после клапана (в 275-м узле, кривая 2) при различных значениях параметров к0 и А. Во всех расчетах используются одинаковые значения То = 100, I = 1, П = 1, к = 1, Со = 0 и С1 = 0,1.

Зависимость плотности перед и за клапаном от времени при следующих значениях параметров ко и А: а - к0 = 0,05; А = 0, 05; Ь - 0, 2; 0, 05; с - 0,4; 0,1; й - 0, 4; 0, 2

Согласно рисунку, a, b при увеличении скорости закрытия клапана усиливаются ударные волны. Графики на рисунке, c, d соответствуют разным значениям коэффициента трения Л, видно, что при увеличении коэффициента трения ударные волны быстрее угасают. Таким образом, приведенные зависимости показывают, что расчетное поведение давления перед и за клапаном качественно хорошо согласуется с реальным.

Работа выполнена при поддержке РФФИ (проект № 17-01-00838).

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

1. Chaudhry M.H. Applied Hydraulic Transients. 3d ed. N.Y.: Springer-Verlag, 2008.

2. Фокс Д.А. Гидравлический анализ неустановившихся течений в трубопроводах. М.: Энергоиздат, 1981.

3. Кобельков Г.М., Соколов А.Г. Об одной неявной разностной схеме для уравнений баротропного газа // Чебышёвский сборник. 2017. 18, № 3. 271-279.

4. Имранов Ф.Б., Кобельков Г.М., Соколов А.Г. О разностной схеме для уравнений баротропного газа // Докл. РАН. 2018. 478, № 4. 388-391.

Поступила в редакцию 22.10.2017

УДК 517.938.5

ИНВАРИАНТЫ ФОМЕНКО-ЦИШАНГА ТОПОЛОГИЧЕСКИХ БИЛЬЯРДОВ, ОГРАНИЧЕННЫХ СОФОКУСНЫМИ ПАРАБОЛАМИ

В. В. Ведюшкина 1

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

Ключевые слова: интегрируемая система, бильярд, лиувиллева эквивалентность, инвариант Фоменко-Цишанга.

A topological (Liouville) classification of integrable billiards in locally flat compact domains bounded by arcs of confocal parabolas is obtained using methods of the Fomenko-Zieschang theory on invariants of integrable systems.

Key words: integrable system, billiard, Liouville equivalence, Fomenko-Zieschang invariant.

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

Интегрируемость бильярда в эллипсе была замечена Дж. Д. Биркгофом [1]. Оказалось, что для любой фиксированной траектории-ломаной ее звенья лежат на прямых, касательных к некоторой квадрике (эллипсу или гиперболе), софокусной с граничным эллипсом. Таким образом, помимо длины вектора скорости (который сохраняется при абсолютно упругих отражениях на границе) вдоль траекторий сохраняется параметр софокусной квадрики, т.е. система обладает двумя независимыми интегралами. В книге [2] В. В. Козлов, Д. В. Трещев заметили, что интегрируемость бильярда сохранится, если перейти к плоским областям, которые ограничены дугами эллипсов и гипербол одного софокусного семейства и на границе которых нет точек излома с углами Щ-. В этом случае все

1 Ведюшкина Виктория Викторовна — канд. физ.-мат. наук, ассист. каф. дифференциальной геометрии и приложений мех.-мат. ф-та МГУ, e-mail: [email protected].