СЕМЕЙСТВО МЕТОДОВ ОБРАТНЫХ ХАРАКТЕРИСТИК Текст научной статьи по специальности «Физика»

УДК 519.63 Научная статья

https://doi.org/10.23947/2587-8999-2023-6-1-53-69 Семейство методов обратных характеристик И. Б. Петров .И, Д. И. Петров

Московский физико-технический институт (национальный исследовательский университет), Российская Федерация, Московская область, г. Долгопрудный, Институтский переулок, 9 H [email protected]

Аннотация

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

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

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

Ключевые слова: сеточно-характеристический метод, уравнения гиперболического типа, уравнение переноса, схема Бима-Уорминга, схема Лакса-Вендроффа, метод неопределенных коэффициентов, условия совместимости.

Финансирование. Работа выполнена в рамках проекта Российского научного фонда № 21-71-10015.

Для цитирования. Петров, И. Б. Семейство методов обратных характеристик / И. Б. Петров, Д. И. Петров // Computational Mathematics and Information Technologies. — 2023. — Т. 6, № 1. — С. 53-69. https://doi.org/10.23947/2587-8999-2023-6-1-53-69

Original article

A family of inverse characteristics methods

I. B. Petrov H, D. I. Petrov

Moscow Institute of Physics and Technology (National Research University), 9, Institutsky Lane, Dolgoprudny, Moscow Region, Russian Federation H [email protected]

Abstract

Introduction. The main idea of the grid-characteristic method is to take into account the characteristic properties, systems of hyperbolic equations, and the finite velocity of propagation of perturbations in the simulated media. Materials and methods. The simplest hyperbolic equation is a one-dimensional linear transfer equation. To increase the order of approximation of the grid-characteristic scheme to the second, you can use the Bim-Warming scheme. If we use a four-point pattern, we get a central Lax-Vendroff scheme. Difference schemes for the linear transfer equation can be obtained using the method of indefinite coefficients.

Results. The grid-characteristic scheme admits a conservative variant, which is relevant if there are discontinuities (shock waves, shock waves) inside the integration domain, while the original system of equations for a matrix with constant coefficients, in partial derivatives, should be presented in a divergent form.

© Петров И. Б., Петров Д. И., 2023

Discussion and conclusions. When numerically solving a three-dimensional problem, the construction is performed similarly, in the case of upper and lower bounds, after scalar multiplication of the scheme by eigenvectors, relations approximating the compatibility conditions with the first order of accuracy are obtained.

Keywords: grid-characteristic method, hyperbolic type equations, transfer equation, Beam-Warming scheme, Lax-Wendroff scheme, method of indefinite coefficients, compatibility conditions.

Funding information. The work was carried out within the framework of the Russian Science Foundation project No. 21-71-10015.

For citation. Petrov, I. B. Family ofmethods ofinverse characteristics / I. B. Petrov, D. I. Petrov // Computational Mathematics and Information Technologies. — 2023. — Vol. 6, no. 1. — P. 53-69. https://doi.org/10.23947/2587-8999-2023-6-1-53-69

Введение. Основной идеей сеточно-характеристического метода является учет характеристических свойств, систем уравнений гиперболического типа, конечной скорости распространения возмущений в моделируемой средах. Обзоры соответствующих работ приведены в [1-6].

Простейшим уравнением гиперболического типа является одномерное линейное уравнение переноса:

du du CD --+ a— = 0, a = const, a > 0, ^

dt dx

а простейшей разностной схемой, учитывающей характеристические свойства (1) — схема «уголок», или схема Куранта-Изаксона-Рисса:

un+l =(1 -о)и" + оы" о = —. (2)

m \ J m m-1" j

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

, \ы" +, - u" , a < 0 ,

п+1 n _ m+1 m' '

u = u — О <

m m

К - <—1, a < a (разность выбирается с учетом наклона характеристики):

un+; =un -T

\a + (u"-u" ) + a-(un -u" ;)!,

I \ m+1 m ! \ m m-1 ! I'

где « += 0,5 (a + |a|), a ~ = 0,5 (a - \a \ ); ="+1 = u"m(u"m+l-и"^ M (u"m+l-2u\ + u"m_i ) (схема с явным выделением диссипативного слагаемого, обеспечивающего ее устойчивость):

,n+1

u

= un -о(Ф" +1/2-Ф" 1/2),

m \ m +1/2 m -1/2

где

ф m+ш = 2 [«(<+i- < )-H(<+i- < )]

ф

n

(потоковая форма). т-1/2 2

В случае нелинейного уравнения переноса:

—+дС=о с=Ц- (3)

Й дх 2

схема «уголок», учитывающая направления характеристик, может быть представлена в следующем виде:

п+1 „ Т I /1+1 ~ /1 , и"т < 0 ' и = и--<

т т Н\г - /■" и" > 0.

Приведенная разностная схема имеет первый порядок аппроксимации по времени и координате. Для повышения порядка аппроксимации сеточно-характеристической схемы до второго, можно использовать схему Бима Уорминга, которая может быть получена на четырехточечном шаблоне * т - 2' -т -1' лт / путем квадратической интерполяции решения на п-ом временном слое по узлам х"т_2, х1,_ 1 > х^1:

1 [4: - <-1 )-H Um - um-1 )]

= < - о (С - <(1 - о)(< - 2u"m+ u"m-2).

(4)

Если использовать четырехточечный шаблон { , х"т_, хт, х"т+1 }, тополучим центральную схему Лакса-Вендроффа:

П О / П / П r\ n . n \

= U m —Z\Um + 1 - Um-1 )+~r\Um + 1 - 2Um + Um-1).

(5)

Эта схема, как и схема Бима-Уорминга, имеет второй порядок аппроксимации по времени и пространству О (т2 + к2), что проверяется разложением сеточных функций ыпт+1 в ряд Тейлора, и является устойчивой при выполнении условия Куранта т < к/а, что можно получить, используя спектральный признак устойчивости Неймана. При добавлении точки х к шаблону, соответствующему схеме Лакса-Вендроффа, получим схему 3-го порядка аппроксимации на четырехточечном шаблоне {хт+1 , х"т _, х"т, х"т+1 }:

= и" - — (2ип , +3ы" - бы" , +u" ,) + — (и" , -2ы" + un ,)-

m б т + ^ m~ ' 2 m m m~ '

3

-—(ы" , + 3ы" , - 3ы" - ы" ).

6 V т+^ тт ш—2 у

(6)

Для повышения порядка точности схемы добавим к последнему шаблону точку х , при этом разностная схема 4-го порядка аппроксимации будет иметь вид:

2

,п+1 . п О /о о . п . . П \ . О

= К, -^^l - ^-1 - и:+2 + (16М-1 - 30< -3 16«:-, - u

2 )-

(7)

6u"m- ^u- 4 с-, )u:+2 )u"m-2).

Если рассматривается линейная акустическая система уравнений вида:

ди + 1 др о

аГ + ~Р "дх ~ ,

(8)

др 2 ди — + рс--= 0,

9/ дх

где и, р — локальная скорость среды и давление, р — плотность, с — скорость звука, то, как отмечалось во второй главе, ее можно привести к виду:

дг дг — + а— = 0

дt дх

--а— = 0,

^дt дх

где г = и + -—-, 5 = и - —--инварианты Римана.

ра ра

В этом случае сеточно-характеристический метод может быть представлен в следующем виде:

(9)

\с1=(1 - ^ )• rm + Ic1=(i - о )• < + о <+1

Соответствующий шаблон имеет вид:

(10)

n+1, m

n, m-1 O-

Jj n+1, m

cj-

-O n, m+1

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

n, m

n, m

Приведенные разностные схемы можно получить следующим образом.

/

и+1

> X

-И

—Хт—х

Рис. 2. Двухточечный шаблон (хп = О, хт= -Л)

На рис. 2 представлен двухточечный шаблон (х: = 0, хт= -И), на котором, в случае схемы 1-го порядка, строится полином 1-го порядка Р1(х)=а1х + а , коэффициенты которого находятся из условий:

Тогда получается:

а1=(<тгС-1)/ К а0=и"т.

п п

Р ( X )а.

-т-1 X + и"

Поскольку ит+1 = и ( х) = и (- Хт), то получается схема «левый уголок» (или КИР):

и"+ = и" -с (и" - и" ,).

т т \ т т-1 '

Если к этому шаблону добавить точку хпт+1 то получается полином второго порядка: Р2 (х) а а2х2 + а1 х + а0, коэффициенты которого находятся аналогичным образом (из условий прохождения полинома Р(х) через значения функции в узлах шаблона):

Ы2 ~ ^т+\ _ и и «т-1" " +П "

а 1 = (и т+а-иС,) /2И, тоа ит.

В этом случае получается схема Лакса-Вендроффа. При добавлении к данному шаблону точки хт-2, получается полином третьего порядка: Р3 ( х ) а а3 х3 + а 2 х2 + а1 х + а 0 с коэффициентами, которые находятся из тех же условий: , п ^ з

а3 = (11 т+1 + 3+з - 3-3-- - ит-2 ) 2 ^ ,

аа2=(1ит+1-:-и^н-1т^2 2И+

1п а и .

В результате получается схема Русанова третьего порядка точности. Если добавить к шаблону точку хт+2, то получается полином четвертой степени: Ри ( х) а аз х4 + а3 х3 + а2 х2 + а1 с коэффициентами:

а5 = (+ит-4-:+1-4u+m-1+u:-2)/2UИ 4,

1з = (2 и^-2 ита+ии^+ит^/Ш3,

а2 а (u6umu1 -зо-т+шт^-и-- -ит-)^2,

а,а(8 ит+и-В u-1-1-u-:uu+u-:-и^/12и-,

и, соответственно, схему четвертого порядка аппроксимации.

п

0

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

Например, на шаблоне {(__ _ х—_— —, ( t_, хт —, (t_ , хт _ — — — (___ х— _ — — , (t__ хт _ _ —} можно представить все линейные схемы с неопределенными коэффициентами {а°; у = -2 _ - 1, 0, 1} в следующем виде:

«Г = а°2 2 u"m_ i + а °3u"m + а 0 u"m+i.

После разложения сеточных функций в ряд Те йлора:

u"m+1 = Unm +т(и',)п + — «)" + — L >;¥ + 0'(т)(

m m \ t / m 2 t / m у V t / m V /

2 6

u"m+l = um+hu )m+у (u "X (<)m + °ih 4)

lm+\ ~ m

m-2 m

-h(u'ml+ у«1+ O'H ,

щих+2 hwx -44- (<!+o'(h4)

принимая во внимание следствия линейного уравнения переноса:

; = <, п"„и г2и-, п^и-гх,

<=иУ)и+, и1 = 12и1, ииих= и приравнивая коэффициенты в левой и правой частях уравнения (10), можно получить условия аппроксимации:

- для схемы первого порядка:

|А°—2 + А°—1 + А°° + А° = 1, 12A°—2 + А,— — А° = а,

А°_2 + А°_1 + А°° + А° — 1, 2А°2 + А°—1 — А° — а, (А°—2 + А°—1 + А° — а2,

(11)

(12)

для схемы второго порядка:

А°_2 + + А°° + А° — 1, 2А°—2 + А°—1 — А° — а, (А°—2 + А°—1 + А° — а2, 8Х°_2 + 1 — — а3,

(13)

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

Если найти все неопределенные коэффициенты из (13) и подставить их в (10), то можно найти первое дифференциальное приближение схемы:

У + X u' = X ( XV + 2 Х2т h + Vh2 - 2h3), * 24 V '

(14)

т. е. на рассмотренном шаблоне можно построить разностные схемы не выше третьего порядка аппроксимации (полагая О (т3 + "3).

В случае шаблона общего вида множество разностных схем может быть представлено в следующем виде:

ы"+1 = Ух'. (с ) ы"+'.;

т / 1 ' "у V / т+] 5

I = 1,0,-1,...; у = 0, ±1, ±2,... Условия аппроксимации первого порядка в этом случае можно записать в виде:

= 1, £(у - (а) = -а.

(15)

Условия аппроксимации второго и более высокого порядков будут иметь вид:

У - ¿я)Р а =(-а)'; Р = 2,3,... Так, для четырехточечного шаблона и двухслойной явной разностной схемы:

, X ; t , X 1 ; t ; X 1 ; ^ ; X ; ^ ; X , 1 }

' т ' ' т—2' ' т—1' ~ т~ ' ог+1 I

можно представить все семейство линейных разностных схем с неопределенными коэффициентами:

Ь, (] = —2, —1,0,1)

в следующем виде:

„,п+1 7-0 „л 7-0 „л 7-0 и 7-0 „л (16)

ит = Ь—2 • —2 + Ь— 1 • —1 + Ь0 • + Ь1 • "ш+1. (16)

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

| а°_2 + + + + = 1, |2а°2 + ° — = о, 4а- а а- а а = =с2, 8а+ аа° - а - =с3.

Эти условия можно легко получить, если разложить все значения проекции точного решения задачи при подстановке в ряд Тейлора относительно любой точки нашего сеточного шаблона, например, (Г, хт):

1+1 = и" + т(и1 )" + т(и„)" + Т(иш)" + О (т4),

п т V От 2 ^ Ш'т V /

«>т+1=«а | ит+и и )т+1 и )т+£ к, )т+о н)

а°и" = а°и",

0 т 0 т'

-И (и)" + — (и)" + — (и )" + О (И4)

V х'т Л V х / т /Т V / т \ )

- 2И (и, )"т + 2И2 (и, )"т - ^ ^ )"т + О (И4 )

Также нужно воспользоваться следствиями уравнения:

и, = -Ч , и„ = . ии, г= = Нхс,

и„, = ---и м^ = -2и и =-2и .

111 ххх' /гх ххх' /хх ххх

13и . и ^ = —и

ххх' Их

С использованием эти выражений можно избавиться от частных производных по времени, выражая через пространственные производные:

»2 2 „3 3

п+1 л / у . к Т / к Т / \» . ^ /_4\

«„ = " М— )и--~ (их )т + О (Т )■

Приравнивая значения ип , (и )п , (и )п , (и )п , в левой и правой частях получаем условия аппроксимации. При а -2 = 0 выводится разностная схема Лакса-Вендроффа, при а -2 = -у ( с -1) — схема Бима-Уорминга,

при а-2 = — (с2 -1) — единственная на рассматриваемом шаблоне схема третьего порядка аппроксимации (Ру-

6

санова).

Результаты исследования.

Предалагается рассмотреть случай линейной системы уравнений переноса вида:

^+А д~и = о, (17)

5/ дх

где А(пхт) — матрица с постоянными элементами и действительными собственными числами X, i = 1,..., п. Тогда эта матрица представляется в виде:

А =

где Л — диагональная матрица, состоящая из собственных чисел: Л = diag {Х1Х2, ..., Хп} = diag {X}; I =1, которые определяются из уравнения: det (А - ХЕ) = 0,

где Е — единичная матрица; ^ — матрица, строками которой являются левые собственные векторы матрицы А с точностью до их длины из системы линейных однородных уравнений

ю,. (А -ХЁ) = 0; I = \,...,п.

В этом случае (17) записывается в виде:

—+о-лп ^ = о,

дЬ дх

затем эту систему умножим на О:

„5и . „ ди „ (18)

О— + ЛП — = 0.

дх

При введении переменных Римана м = Ои последнее уравнение представляется в виде:

дм .дм _ (19)

— + Л— = 0,

дг дх

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

^ = (х /); 1 = «•

Функции называются инвариантами Римана системы (17); X, являются скоростями распространения возмущений в среде. Общее решение рассматриваемой системы дифференциальных уравнений в частных производных есть сумма п бегущих волн, каждая из которых распространяется со своей характеристической скоростью X.:

и=Ё6? ^ (* /) •

1=1

Следует отметить, что если векторы Ь, нормализованы, т. е. | Ь. | = 1, то величины м . можно интерпретировать как амплитуды соответствующих бегущих волн. Можно представить (18) в скалярном виде:

ди1 ч ди, „ . , (20)

а,—'- + Х1щ—'- = 0, I = 1,.., п. у '

д/ дх

Далее в области интегрирования { > 0 , 0 < х < X} вводится разностная сетка

{е = пт;п = 0,1,...;х,. =(т _ 1) А;т = 1,..,М;X = (М _ 1)h} и проводится аппроксимация (20) с помощью конечно-разностных соотношений с учетом наклона характеристик:

dt

(или с учетом знака X ):

5х =

n+1 n n _ n

щ Um _ Um тХю иш+1 + & f = 0 i = 1 n. (21)

i ill iJ ' т л

Последнее соотношение может быть представлено в матричном виде, если положить:

|Л| = {|} |(Л +-},Л-|(Л-|Л|) : (22)

Q(u:+1 - к) - сЛ+Q (ul| - ы"а) -сЛ-£2(ul + -«:) + xQ/ = 0. Поскольку матрица Q не вырождена, то полученную разностную схему можно представить в следующем виде:

„Т = К - oA (u"m+1/2 -u+1/2)-т/+|£2- —(-2u"m + (23)

1 / n I n \

= 2 (U™ ±1 + M »

и

m ±1/2

2

В этой схеме явно выделено слагаемое:

- о-1 Л О (<-- 2«;+<+1

которое обеспечивает ее устойчивость при выполнении условия КФЛ (Куранта-Фридрихса-Леви):

1

о <-т-т •

Можно сказать, что в исходную систему вводится дополнительный «диссипативный» член, обеспечивающий монотонность полученной схемы первого порядка аппроксимации, но являющийся также нефизичным источником, имеющим разностное происхождение (аппроксимационная вязкость).

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

и"+1(Е - стО—ЦО)и" + стО— 1Л+ Ои" , -стО—1м" . -т/.

т у \\}т т—1 т—1 ^

Рассматриваемый численный метод можно также представить в следующем виде:

1+1 1 _ и„ = и_, - о

(□ -1ЛО)Дм + ор (О-1 |Л|р О) Д2и,

где

п+1

м„ = и„ -

и^-с^ог'ло^Дм+а13 (о-1 |л|рд)д2«,

Эта схема имеет первый порядок аппроксимации при в = 1 и второй — при в = 2.

Можно представить данную схему, как схему с весовым коэффициентом а (0 < а < 1):

и"т+1 = и"т -с ( О-ЛО )д и + 2 с (а-11 Л|р О ) Д2и +(1-а )ст2 (О-1АО )д2и.

При а = 1 получается схема первого порядка аппроксимации, при а = 0 — второго. При 0 < а < 1 получается схема, имеющая формально первый порядок, но, при надлежащем выборе значения а, позволяющая улучшить описание численных решений с большими градиентами за счет уменьшения аппроксимации вязкости [7].

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

ди „ дФ ди ди , ч ди ди Л

— + Л-= — + —( А и )= — + А — = 0, (24)

дГ дх д Г дхУ ' д I дх

где А = ——; Ф = 24 и; Ф = Аи; поскольку рассматриваемая система является гиперболической, то: А =

ди

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

"т = а (ф:+1/2 - фт1/2 ) + а [Бт+1/2 (ы"т+1 - ы"т) / к - Бт_1/2 (< - ы"т - ) / к ], (25) где ф : ±1/2=2 (Ф т ±1 + Ф :).

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

Ви'Хх, где В = 2 О"1 |Л| О.

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

Нужно отметить, что среди схем первого порядка аппроксимации данная схема имеет минимальную аппрок-симационную вязкость. В случае трех независимых переменных {/, х, у} консервативная и неконсервативная схемы будут иметь, соответственно, следующий вид:

п+1 _ п _ Т Т ф п _ф п \ _ Т Т ^ п _ ^ п \ _ Т

- — ~ и ш1 — —"»(=1Фт-1 12,1) ) уТГ ^ т,1-1/2- —

' \ I-1/-,! (и Пт+- ~и ^ ) / V1-у(и "и1~ и 1-г,1 ) / " ]

+

+ — \о2 2П(ип м-и п\1Ъп - £»2 , -п(и" ,-и" , Л /й!

7 I т,1+1/2 у т,7+1 ш1 ) 2 т,7—1/22 \ т,7 т,7-1 - '

2

п 11-и--т,)-^^^^ лЛ,- о +1)]+

где Бк = ^ С!-1 |Лк| Ск; к = 1,2;

матрицы Лк, Л±к, О как и в одномерном случае, состоят из собственных значений и собственных векторов матриц А; в области интегрирования ^ > 0; х < х < X, 0 < у < У] введена разностная сетка

{п = пт, п = 0,1,...; хи = (т _ 1) к1; т = 1,..., М; уе = (1 _ 1) V =1,..Д,; X = (М _ 1) к„ У = (Ь _ 1) к2 ].

Аналогично строятся и разностные схемы для трех переменных.

Следует вернуться к соотношению (18), которое можно представить в скалярном виде:

где каждая компонента вектора:

ди , ди (сu . ди | „ . ,

ю,--+ X.ю.— = ю. I--+ X.— | = 0, г = 1,..,n,

д/ дх \ dt дх

ди . ди

--+ к, —

dt дх

является уравнением переноса для соответствующей компоненты вектора и:

ди. ди.

+ X. = 0; ] = 1,.., п. дt дх

Причем для каждого из этих уравнений можно построить уже известную схему:

«п =2>; (т, а , К++ ) п •

(27)

(28)

(29)

Далее нужно заменить вектор (28) его разностной аппроксимацией:

со ,м"+1 -Усх," (о )ю и

1 т л I! 1 .

I ди „ ди |

i = 1,.., И,

при этом очевидно, что каждое уравнение в (30) есть ,-ая компонента вектора:

„5 и . „ ом Q — +ЛО — = 0;

dt дх

т. е.

„ д и . „ д и Q-+ЛО- :

dt дх

Юi)юi '

Ю X+IXK )ю „-и

Y

m+Ц

(3°)

(31)

При умножении полученного соотношения на матрицу О ', получается общий вид разностной схемы для численного решения системы уравнений гиперболического типа (17):

ип+=Уп-1Б'< Оип++ (32)

т / , Ц т+ц>

где В ^ = diag {аЦ o¡ )} — диагональная матрица с компонентами, являющимися коэффициентами а1^, определяющими конкретный вид схемы.

Например, в случае рассмотрения разностной сеточно-характеристической схемы, получится (22), (23), (24):

или, в другом виде, С = ЕП)и"а+2П)(и"а_1 -и"т)-о (П1 Л^)(<+1 -и"а); для полученной ранее для скалярного уравнения схемы Лакса-Вендроффа (5):

(33)

+ -2-(ф-1 - Ф:+1 ) + —(п- |Л| п) (и I-1 - 2и1 + и1+1). (34)

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

,+1=u-+f (ф^ - ф- )+f (q-1л2 q) (« (! -2и :+и:+,)+

(35)

+ % (Q-1GAQ)(M(, — 2Q _( +

и(>, - и(

(q-1c?|A|Q)(m(,^ ^^^ Q1 + 6 и 1+1+и ^+2),

где О = diag {я1(о1)], , = 1, ..., п — диагональная параметрическая матрица.

Данная схема также может быть представлена как схема «предиктор-корректор»:

Так, при

=и К-1-Ф —)+у^Но)^-2 и"т+и "т+1),

«й'^+п-'^л2-уНЦ«й-1-2«й+«й11)+ +а1аа- {й"т-2 и т + й^)- -С--2 ;;1( чъгиц+1)].

8,

получается схема 3-го порядка точности (аналог схемы Русанова).

?

/К

т-1/2 т+1/2

Т—Г

п+1

т-1

к т т+1

Рис. 2. Схема предиктор-корректор

Сеточно-характеристический метод может быть представлен также как метод потоков, или как интегро-интер-поляционный метод, если уравнение (23) будет записано в интегральной форме:

цГди + _ ах = ийх _фа() = о, (36)

и далее будет аппроксимироваться последний интеграл:

/ +]ффа- / «-л-/фф = о.

х-1/2 хт-112 гп

Из последнего равенства следует, если проведется аппроксимация интегралов методом средних, получится:

^т+1/2 т

(интегро-интерполяционный метод), или же:

ип++ = ип - о (фи+!/2 -Ф"+!/2)

т т ^ у т+!/2 т-!/2}

(37)

(38)

(метод потоков).

Стоит отметить, что в зарубежной литературе эти методы получили название методов конечных объемов. Понятно, что (38) представляет собой семейство численных методов, свойства которых зависит от способа вы-

л+1/2

числения потоков фв+1/2 ; например, линейного уравнения переноса (1) потоки могут быть вычислены следующим образом, но с полуцелым верхним индексом

Ф"

Ф"

= 2\\["(и ™++ - 4 ш )" П(и * ++ и* )] + 2 [ \

П+1 п+1 \ I \(

и т+\/2 " И» )-И(

п+1 _ п+1

и т+1

)]

= 2 {{ [\ и"ш ~ ит-1)-\а\(и"т~ и1 —1 [ а(ип - ип\)-\а\(ип1 - И^)^.

С2 — 1

6

X

п

Аналогично в потоковой форме могут быть представлены методы (33), (34) и др.

В двумерном случае сеточно-хараткеристический метод первого порядка аппроксимации может быть представлен в следующем виде [2]:

п-1-1 п , 1 г.

Uml = Uml + bl

где

+ TLl

(40)

him, = Gl [(О," Ц " J " (Ц Ц (" J [((ОХ, (Um" " U m, J " foXOL (Um+ " U m, )"

b ~ml = о,

л+ = 1 kk Л" = 2(Лkk Лk = (Ml- лk = (xi

2 1 г чг-' 2'

1 = 1,..., N — диагональные матрицы; №. — ихсобственные числа; £1к = {со*.} — неособенные матрицы, строками которых являются линейно-независимые собственные левые векторы сок. матриц Ак

Для матриц Ак = {ак.}, в случое системы уравнений механики деформируемого твердого тела:

1к. = V, + /., 1 = 1, ..., 7, к = 1,2.

1 к ^ V ' ' ' '

Используя обычные кинематические соотношения для компонент симметричного тензора скоростей деформаций [8] в фиксированной в пространстве криволинейной ортогональной системе координат х1 и х2, х3:

1 1 d U + _Н + ±

2 Н„ dx н Н dx= И„

m m

£ UL^

s=1 HS dxs

1

2И„

И, еи„

V„-- + ■

2, и = 1,2,3.

Выбирая определяющие соотношения в виде:

d у* _ d^jtdi±d^L. иd °t _ v о

- z^ %

dx dx„

i,j = 1,2,3,

й/ Л Н1 йх1 Н2 йх2 Н3 йх3 т ,„=1 нужно записать замкнутую систему двумерных нестационарных уравнений в виде:

ди + а ди + А ди _ .

дг 1 дх1 2 дх2 ^

Строго говоря, вместо субстанциональной производной дв. / dt следует использовать, так называемую, производную Яуманна для компонент девиатора тензора напряжений типа:

3

dS dS / _

1

СО,.,. = — J 2

dt

г

-Sjk )>

—'---- I, S.. = ст.. -5..-> с

' - - j -j

dx.

dx.

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

Здесь символ и = {у1, V , а11, о22, с33, Т} — вектор искомых переменных, включающий компоненты у1, у2 вектора скорости V (по осям — х1, х2 соответственно), отличные от нуля компоненты в.. симметричного тензора напряжений и температуру Т;х1, х2, и) = /2,/п,/12,/22,/ъъ,/Т} — вектор правых частей со следующими компонентами в криволинейной ортогональной системе координат х1, х2, х3:

А = = + -Л = = +-

1 22)

P HH 2

1 C - 13 3-)

P H 2H3

дхд

дх 0

H, H з

н Hj

H

H

2 dHj + 1 5H3H V2

dx2 H3 dx2 ) HXH 2

i 2 dH 2 + 1 5H3H"

L H2 dx1 H3 dx1 J H 2H

8H,

дх,

V 5H2 2 -

сЗх,

А = g-,11^2 дн2 (q( 2-gj- ^^^ H=JH2 дх2 2HxH2

dH,

dx.

■ + V2

<3H2

c3xj

+■ +

HlH2 dxJ H3

Vi 5H3

V2 5H3

H"j <2xj H2 dx2

dH 2 *

^XJ у

dH

дх2

P c

P С

11 v2 fJH-L

Hj H2 dx2

HH 2

dJH+i v dH2

' T" 2

cX0

5x,

, V22 v2 <'H2 + °33

HXH2 dxx

H,

V 1 ^3. + . V2

dH.

VHi 5xi

H2 5x2

Здесь Fl, F2 — компоненты вектора массовых сил, с Т — внутренняя энергия для термоупругой среды, с — удельная теплоемкость материала, Т — температура, Q — объемная плотность источников тепла, Н,, I =1, 2, 3— коэффициенты Ляме, характеризующие выбираемую ортогональную криволинейную систему координат, р — плотность, определяемая уравнением состояния твердого тела, например:

1п Р = -( 3К )-1 £ а,,- 3а Т, Р0 1=1

р0 — плотность материала в недеформируемом состоянии, К — коэффициент объемного сжатия.

В соответствии с предполагаемой здесь двумерностью деформируемого состояния, перемещения материальных точек в направлении х3 отсутствуют и<;13=(!23 = у! = 0, д / дх2 = 0. Матрицы Ак = {а,к }, к = 1, 2, , = 1, ..., 7 при этом имеют вид:

А = ^

1 Н

А = А_

Н 2

0 - 1/Р 0 0 0 0

0 0 — 1/Р 0 0 0

-Чш1 —(?1112 + 91121)/2 П 0 0 0 0

?1211 — ( ?1212 + ?1221) /2 0 ^ 0 0 0 5

?2211 — ( ?2211 + ?2221) /2 0 0 0 0

?3311 —(?3312 + ?3321) /2 0 0 0 V! 0

- СТ117Р1 - СТ117 Р1 0 0 0 0 ^

0 — 1/Р 0 0 0 0

0 ^ 0 —1/Р 0 0 0

-(91112 + 91121)/2 —?1111 ^ 0 0 0 0

— (?1212 + ?1221 ) / 2 — ?1211 0 ^ 0 0 0

— (?2211 + 92221 ) / 2 —?2211 0 0 V2 0 0

?3311 — (?3312 + ?3321) /2 0 0 0 0

- СТ12 / Р2 — СТ22 / Р 2 0 0 0 0

Для принятой в настоящей работе модели Прандтля — Рейсса компонеты тензора четвертого ранга д,ш упру-гопластического материала имеют вид [8]:

Чцы = ^8у8И -

§ у

§Л §у - §ук Ъа )-

I ЦУ

у к1

рс К к

где X, д — параметры Ляме, k — предел текучести на сдвиг, 5т — символы Кронекера, девиатор тензора напряжений

: с -5 С

тп тп ~ / , ~ ^

3 X=1

у = (3Х + 2д) а (а — коэффициент линейного расширения материала при нагреве). I определяется из условия пластичности Мизеса

|0 при 5 = sll + 5222 + 5222 + £2 < 2к2, 11 при 5 > 21с2.

I =

Определяющие соотношения при I = 0 являются обычным «законом Гука» для упругих изотропных материалов. При у = 0, Q = 0 в первые 6 уравнений не входит температура и их можно решать независимо от уравнения энергии:

к к .У к-У 7

ак + (ак - 4<4У'

к к У 2 = У 6 =

ак- (<*1-4<*кк'

У к=-Ук =Ук = 0.

> к к (( к к кг кг \\ . к к !( к к к к \\ . к к (( к к к к \\

кк = «13«24 ((«31«42 - «32«41 )) + «1 4«25 Х ((«41«32 " «42«31 )) + «13«25 ((«31«32 - «32«31 )) >

к = 1,2,

Здесь:

Ц=

q2=

raj 1 ra1l га1з ra 14 ra15 0 0

Ц ra2i 1 га2з ra 24 ra25 0 0

0 0 ra 32 ra 34 1 0 0

= 0 0 ra43 ra 44 0 1 0

Ю1 5 0 0 ra 53 ra 54 0 0 1

ral1 1 -ra 23 -ra24 -ra25 0 0

1 ra1l -ra13 -ra14 -ra15 0 0

rali 1 га1з ra ?4 2 ra 15 0 0

1 2 ra 22 2 ra 23 21 ra 24 ra2 25 0 0

0 0 1 ra2 34 ra2 35 0 0

0 0 0 ra 44 44 ra 45 1 0

0 0 0 ra 54 ra 555 0 1

1 ra 2i 2 ra 23 2 ra 24 2 -ra25 0 0

2 ral1 1 2 "ran 2 -ra14 "ra?5 0 0

2 - a!

a,,

(у2 )2-a?

ю,., =

(yi)

(y; )-a?

2

к _ к к . к к 1 _ ] г\

^ ^ — 1 ^ 11 I ^лг a с 1, k —

al3® ,'l ra

k ;

У/ k a22 ra

k

®,-5 = k

4 a42 P/'l a4lPi2

®/з = -4^4--Ты-

a4iPl2 a4lP

4l 2

^ 2 _ a52 Pil - a5lPi2

a22p4l

k k k k a, ,ю,., + a,.ran

= l,2,

^k _ a32Pll allPl2

a3lP42 — al2 P 4l

a42 P/l a4lP 2

5l H 42

a52lP22 - a52P4/ ''5 a^ - allP

Pff = al j ю15 + al j ^6 + alj ^ в2 = a2jЮ]з + a62.,-^2 + a7 j® 27 ' j = I2, ' = 3,4,5

о={ ^}.

В этом случае для обратных матриц:

где:

Рп = Р 22 =

P'I Pl' 0 0 0 Plk2 Plkl

k Prn к P22 0 0 0 A- P22 к P2l

P3I Рз2 P3k3 0 0 P32 - РЗ1

ph P42 р4З 0 0 P42 - P4I

P2kl P52 Р2З 0 0 P52 P2I

Рб^ Рб2 Рб2з l 0 Рб2 - Рб^

P7I P772 Р7З 0 l P772 - p7I

( a2) 2 2i 2 = ю 12 2А2 ' Jp-12i = ю 2 i 2 A22 '

a- =1-ю; 2

„ 1 „1 „1 1

Рз1

2А1

„1 _ ^СлСк

132

2А1

1 _ 2515^24 Рзз ~

2А1

1 ш23 Ра1 = —

11

ю ,, - со „ со „

2А1

2

С0„, =

2l

l2

a

2

CO., =

со,, =

ll

l2

2

k

= ю15ю3з ~ „1 = 2-)_®2З

„42 2А; Аз А; :

и р1 _ 1)1 .з-24 с®ю3з

Р51 = -ГТ1-' =52 "

2А; ^52 2А;

, _ 1214 ® 2з ®1з 11 24

.г 53 _

И, аналогично:

А1

Рр ^^^^^(,з (с-®-®)5®1з))(2А1) >

Рп ^ГО^Нз0®^®0 )-®Р2,з(®114И 2Л1) ^^

[3 = [-1-2,3 2—3-2о - И^-®) --2-2,4 (Р1зС02о И 1--2С )) )А1 ) ' ) А- =С0115 (|^^4,Р1Сз ■— СО -ЗзСО 31:.-)) + ю) (со Ъ-3- - Ю-.Ю 4 ) + (р)Ю^ -О-зЮ®).

2 _ ~ ^з! „ 1 _ <°2з<°1з

Р43~ Д2 ' Л4" 2 Д? ■

2 _ -2 _ «пЗ^" ^З ЮМ

Рб2 , ¡2 ' Ръ3

2 2 1 53 А2

Р2х = \_42,- («24 - ) - °22 2,4 (Ю25 - -23®35 )] ) 1 2

Р4 = [ю,2-2,4 (о>15 - 0-222 ) - «Б22,5 ^ - и-22 )] (2А2 ) ,

„1 Г 2 2 2 2 2 2 \ 2 2 2 2 2 „2 \ "1 ( \ 2 \ Р,3 = [<3)^4 [А " С0153--123 - " Ю<- 2,5 («М®2 -^и™2 )_|(Д ) ,

I = 6,7,

А2 = Ш-, (юг2 -3 3Й2^СЙ44 ) -4- 2 (о^-О2 ) + (Юм- ~ 031250-, ) •

При построении расчетных формул на границах прямоугольной (в координатах t, х1, х2 ) области интегрирования ограничимся рассмотрением только верхней (х2 = 0 ) и нижней (х2 = 1) границ, имея в виду, что остальные границы (х1 = 0, х1.) часто являются плоскостью (или осью) симметрии или периодичности решения либо выбираются таким образом, чтобы за рассматриваемое время t < t возмущения не достигали этих границ. Обобщение на случай более сложных условий на границах не представляет принципиальных трудностей и аналогично рассматриваемому ниже.

При умножении скалярно на собственные векторы (Ю2)" получатся соотношения:

к21х+=в2 = (о ?£ис+т/:)±

+

7- I, (ю 2 )т («.т + «О )1 = Ъ -..,7,

к2

аппроксимирующие с первым порядком точности условия совместности вдоль линий пересечения характеристических поверхностей системы и координатной плоскости х1 = х (с уравнениями дх = А2, dt):

ю2и, + А2.ю2и . = ю2. /-Аи Л г = 1,.., 7.

1 t 1 х2 1 v 1 х27' ' '

Как известно, число граничных условий для гиперболической системы уравнений типа определяется числом отрицательных (положительных) собственных значений матрицы Аи на верхней (соответственно, нижней) границе области интегрирования. В рассматриваемых ниже задачах на верхней границе х1 = 0 имеет место А27< А26< 0, на нижней границе х2 = 0 — соответственно, Х21 > № > 0 и; следовательно, на каждой из этих границ требуется постановка двух граничных условий, которые записываются в виде:

х и1,..., и7) = 0, 1 = 1, 2 при х2 = 0, Фl.(t, х1, и1,..., и7) = 0, 1 = 167 при х2 = 1,

причем необходимо, чтобы det Ф 0, det Ф 0,

где, соответственно, 0_=||ю1 Ю2Ю3...ю7|| , □+=||ю1...ю5Юбю?| .

Здесь ю , = {дф, / ди, , . . . , дф, / ди7},, = 1,2,6,7, а , = 1,...7 — собственные векторы матрицыЛ2. Для рассматриваемых ниже задач граничные условия выбирались полулинейными и после их аппроксимации имели вид:

ф, =ю,^"+1,уа+' -g¡(г,)= 0,

i = 1,2 при х2 = 0, i = 6,7 при х2 = 1. Возможная схема расщепления для случая двух переменных:

В .chi „ди — + A — + B— = 0: dt дх dy

п (ди , ди + _ ( ди в1„ ди | >

Pi + Р2 = ii(3sep= >o i

= в i " 1 "ш+вв" ■

+i-1 2"2m/>

где иа1 и и2а1 находятся из численного решения двух одномерных систем уравнений:

ди „ 1 ,ди „ ди „ 1 ,ди — + Р-1Л— = 0и— + вгЛ— = 0. дt дх дt ду

т-> """ "" - ди ^ ииk В случае пространственной динамическои задачи:--+ У Ak-= 0

Pit ' ' Р)л>*

ди.

dt

дх.

матрицы имеют вид: riJkl =

=1 (q,yk/+j)1 * k;

A -

A -

A„ =

13 ,CT22 , CT23 °33>U} - вектор-сто лбец искомых величин.

V 1 0 0 —1/P 0 0 00 0 0

0 V1 0 0 — 1/P 0 00 0 0

0 0 V1 0 0 — 1/P 00 0 0

- 01111 — r 1112 0 V1 0 0 0 0 0 0

- 01211 — r 11221122 0 0 V 1 0 0 0 0 0

- 01311 — r 2221122 0 0 0 V 1 0 0 0 0

— 02211 — r ' 1312 0 0 0 0 V1 0 0 0

- 02311 — r 3212 0 0 0 0 0 V1 0 0

- 03311 — r 3331122 0 0 0 0 00 V 1 0

-°11 / P — °12/ P 0 0 0 0 00 0 V 1

V 2 0 0 —1/P 0 0 00 0 0

0 V2 0 0 — 1/P 0 00 0 0

0 0 V2 0 0 — 1/P 00 0 0

—r 1112 _ 01122 0 V2 0 0 0 0 0 0

—' r 11221122 — 01232 0 0 V2 0 0 0 0 0

— r 11331122 ~ 01322 0 0 0 V2 0 0 0 0 5

— r 212122 02222 0 0 0 0 V2 0 0 0

— r 22312 _ 02322 0 0 0 0 0 V2 0 0

— r 33312 _ 03322 0 0 0 0 00 V2 0

" °22 7 P 0 0 0 0 00 0 V2

V3 0 0— 1/P 0 0 00 0 0

0 V3 0 0 — 1/P 0 00 0 0

0 0 V3 0 0 — 1/ P 00 0 0

— r 1113 — r 11123 — 01133 V3 0 0 0 0 0 0

— r 11221133 — r 11222233 — 01233 0 V3 0 0 0 0 0

— r 1313 — r 1323 — 01333 0 0 V3 0 0 0 0

— r 2213 — r 2222223 -0 '2233 0 0 0 V3 0 0 0

— r 2313 — r 223323 — ^2333 0 0 0 0 V3 0 0

— r 33313 — r 3323 -0 3333 0 0 0 00 V3 0

"°13 /P "' °23/ P - °33 1 P 0 0 0 00 0 V3

и

V, , v , v

2> '3

12

Одна из возможных сеточно-характеристических схем для численного решения приведенной системы может быть представлена в виде:

и::,: = иЩ - о1О-1ЛАти + о [ О |Л!|т ОА2ти - о2О-1!Аги +

Iп+1 т1р

+ от2От |Лт|Т ОтА2и - озО-Л 1Ари + о^Оз|Л3|г ОтА2

_о ^ ( п _ п \

ит ~ т+1,1, р ит-1,1,р) ,

и = 0,5 (и"т,1,р+1 - и"т,1,р-1 )' и т = 0,5 (и"т+и, р _ 2 и"т1р + и"т-1,1, р ) ' и I = 0, 5 (и т,1+1,р - 2ит1р и т,1-\,р ) ,

= 0,5 (и

п О " _1_ '

т,1, р+1 _ 2 ^ т1р т,1, р-

1),

кк (ткк—шаги по времени итрем 1соординатам; к=1, 2, 3), тли

..«+1 ..« I _

и . = и , + о

т1р т1р

• ^-1 >+А- 1р А" и -(О-111О1 )тр А

+ о„

т1р

["(О^А ) А"и -(О- 11+О2) А+и 1" + о3Г(О- 11+А ) А"и -(О- 11+О3) А

|_\ 22 т)т1р т V 2 2 2 !т1р т _ 3 |_\ 3 3 т)т1р т V 3 3 3 ) т1р '

+

А± и т,1±1,р ит1р ,

А± = и , ,, - и , .

р т,1,р±1 - тр

Схема имеет первый порядок точности при у = 1 (что и было реализовано в работе), второй — при у = 2; здесь явно выделен «демпфирующий» (вязкий) член, придающий схеме устойчивость и положительную (по Фридрих-су [5]) определенность (или монотонность — в одномерном случае); при реализации нет необходимости находить матрицы Л±к; схему несложно гибридизировать (т. е., например, полагать у = 1 в областях с большими градиентами параметров течения и у = 2 — в «гладких» областях). Преимущество записи сеточно-характеристической схемы состоит в уменьшении количества арифметических действий.

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

>п)и~2'л(а2) А иТ}± с(V)) ))) *т»> ' = и.,/ -

'т1р ± V2 т!р ' ^ I 3 V ' 'т1р

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

га /М + к;3га3и' = -гаи (4К + Пм' ) ■

В случае одномерной системы уравнений газодинамики:

ди . ди _ — +А — = 0

д1 дх

вектор искомых функций и матрица А имеют следующий вид [2]:

и р 0

и =)

^ А = 1

р" 1 дР

дs

0

и

" 1 др Р д

Р / Р

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

Х1 = и + с; %2 = и; = и - с, где с — скорость звука в газе; матрицы из собственных векторов, собственно, будут:

Q =

ю,

ш.

ю,

dp dp -t- p c —

д p дг

p 0 - P2 dp dp

- p e

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

1. Рождественский, Б. Л. Системы квазилинейных уравнений и их приложения к газовой динамики / Б. Л. Рождественский, Н. Н. Яненко. — Москва : Наука, 1968. — 546 с.

2. Магомедов, К. М. Сеточно-характеристические методы / К. М. Магомедов, А. С. Холодов / 2-е изд. — Москва : Юрайт, 2017 г., — 312 с.

3. Куликовский, А. Г. Математические вопросы численного решения гиперболических систем уравнений / А. Г. Куликовский, Н. А. Погорелов, А. Ю. Семенов. — Москва : Физматлит, 2012 г. — 656 с.

4. Favorskaya, A. V. Innovation in Wavt Processes Modeling and Decision Making Grid-Charakteristic Method and Applications / A. V. Favorskaya, I. B. Petrov // Springer. — 2018. — 270 p.

5. Fridrichs, K. O. Symmetrics hyperbolic linear differential equation / K. O. Fridrichs // Communications on Pure and Applied Mathematics — 1954. — Vol. 7, no. 2. — P. 345—392.

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

7. Петров, И. Б. Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики деформируемого твердого тела / И. Б. Петров, А. Г. Тарасов, А. С. Холодов // Журнал вычислительной математики и математической физики. —1990. — Т. 30, № 8. — С. 1237-1244.

8. Петров, И. Б. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом / И. Б. Петров, А. С. Холодов // Журнал вычислительной математики и математической физики. 1990. — Т. 30, № 8. — С. 1237-1244.

Поступила в редакцию 02.02.2023. Поступила после рецензирования 01.03.2023. Принята к публикации 02.03.2023.

Об авторах:

Петров Игорь Борисович, член-корреспондент РАН, доктор физико-математических наук, профессор, Московский физико-технический институт (национальный исследовательский университет) (РФ, 141701, Московская

область, г. Долгопрудный, Институтский переулок, 9), Math-Net.ru, MathsciNet, eLibrary.ru, ORCID, ResearchGate, [email protected]

Петров Дмитрий Игоревич, кандидат физико-математических наук, лаборатория прикладной вычислительной геофизики МФТИ, Московский физико-технический институт (национальный исследовательский университет) (РФ, 141701, Московская область, г. Долгопрудный, Институтский переулок, 9), Math-Net.ru, eLibrary.ru

Конфликт итересов:

Авторы заявляют об отсутствии конфликта интересов.

Все авторы прочитали и одобрили окончательный вариант рукописи.