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

УДК 519.612:632.4

Б. И. Березин1, Д. В. Благодатских2 ДОСТАТОЧНЫЕ УСЛОВИЯ УСТОЙЧИВОСТИ

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

Используя свойство положительной определенности квадратичных сеточных форм, получены достаточные условия устойчивости схемы КАБАРЕ, в том числе и для практически важного случая доминирующей конвекции.

Ключевые слова: конвекция-диффузия, доминирующая конвекция, численные методы, схема КАБАРЕ, условие устойчивости.

1. Введение. Одной из наиболее важных и актуальных задач современной вычислительной математики является моделирование переноса пассивной примеси в газодинамических и гидродинамических потоках. В роли примесей могут выступать аэрозоли [1], газовые смеси [2], загрязнения в океане [3]. Перенос пассивной смеси является серьезной вычислительной проблемой, если в нем преобладает процесс конвекции (адвекции), как в приведенных выше случаях. В то же время алгоритмы численного моделирования атмосферных и океанических течений предполагают большой объем вычислительных ресурсов. Таким образом, от искомых алгоритмов требуется максимально возможная скорость исполнения на суперкомпьютерах с параллельной архитектурой.

Разностная схема CABARET (Compact Accurately Boundary-Adjusting high Resolution Technique) была предложена в работах А.А. Самарского и В.М. Головизнина [4, 5] изначально для случая одномерного линейного уравнения переноса. Первоначальный вариант этой схемы был получен в трехслойном виде. CABARET аппроксимирует исходную задачу со вторым порядком точности по пространству и по времени. Как было показано в статьях [4, 5], схема CABARET не вносит в результаты расчета какой-либо дополнительной схемной вязкости и обладает улучшенными дисперсионными свойствами. Также было установлено, что схема CABARET устойчива при числах Куранта от нуля до единицы. Схема обладает фиксированным компактным шаблоном и допускает эффективную реализацию на многопроцессорных вычислительных комплексах.

Позже в работе [6] схема CABARET была представлена в эквивалентной двухслойной форме для линейного одномерного уравнения переноса, допускающей обобщение на содержательный случай газовой динамики [7-9]. Эту двухслойную форму и будем рассматривать в качестве основы для обобщения схемы CABARET на многомерный случай.

В дальнейшем возник вопрос обобщения схемы CABARET на многомерные случаи [10, 11]. В статье [11] двумерная модификация схемы CABARET была исследована на устойчивость методом спектральных гармоник (метод Неймана). В этой работе было получено только необходимое условие устойчивости, к тому же модификация схемы кабаре, предложенная в [11], где потоковые переменные располагаются в узлах, пока не получила практического применения.

Кроме того, необходимо добавить, что даже в работах [4, 5] был проведен анализ устойчивости только для чисто гиперболического случая, без учета диффузии. Таким образом, существует потребность в получении достаточных условий устойчивости для схемы CABARET для многомерного уравнения конвекции-диффузии.

Несмотря на то что схема CABARET в двухслойном виде формально эквивалентна некоторой трехслойной схеме, она тем не менее не сводится к операторному виду, приведенному в работах А.А. Самарского и А.В. Гулина [12, 13]. Поэтому использовать напрямую упомянутые в данных работах теоремы об устойчивости трехслойных схем не представляется возможным. Далее будет иметь место подход, подобный тому, что применялся в статье [4], т.е. получение априорных оценок норм сеточных функций с помощью непосредственных алгебраических преобразований, с использованием операторного подхода в качестве вспомогательного.

1 Факультет ВМК МГУ, доц., к.ф.-м.н., e-mail: berezinQcs.msu.su

2 ИБРАЭ РАН, мл. науч. сотр., e-mail: blagodatQibrae.ac.ru

2. Реализация схемы CABARET для линейного уравнения переноса. Уравнение линейного переноса имеет вид

где с > 0. Здесь (р — переносимая вдоль оси х скалярная величина и с — скорость переноса. Введем сетку по пространству и по времени

Шп = {Хг = ih, ^оо < i < +00} U {xi+i/2 = ih + /i/2, ^00 < i < +00},

Шт = {tQ = 0 < i1'2 < t1 < ... < tn < fn+l/2 < tn+1 < ... < tM = T},

где h — размер ячейки по пространству, г = tn+1 — tn, t1/2 = tn + т/2, г • M '/".

Рассматриваются два типа переменных — консервативные ф, относящиеся к узлам с дробными индексами, и потоковые (р для узлов с целыми индексами. Уравнение (1) по схеме CABARET в ее традиционной двуслойной форме записывается так:

,п+1/2 / п —1/2

^»+1/2 ~ ^i+1/2 + сУ?+1 ~ У? = 0 (2)

т h '

=2<++1/22-^- (3)

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

3. Схема CABARET для многомерного линейного уравнения переноса. Рассмотрим N- мерное уравнение переноса

д(р А д(р

тг- + / са —— = 0, (4)

dt ^ дха w

а=1

где са > 0.

Введем в iV-мерном пространстве равномерную по каждому из направлений расчетную сетку с координатами узлов

= {xta = iaha, ^ОО <ia< +00} U +1/2 = iaha + ha/2, ^OO < %a < +00},

где а € [l,...,iV], ha — размер ячейки в направлении а, и сетку по времени шт, такую же, как для одномерного уравнения переноса.

Рассмотрим множество расчетных ячеек fi^. Каждая расчетная ячейка представляет собой N-мерный параллелепипед с центром в хс = • • • > ж"а+1/2' • • •' ж^+1/г)Т и с центрами "гра-

ней" в xia = ... ,х?а, ■ ■ ■ 1/г)Т- Введем сокращение для суммирования по индексу %а —

Xia+ia = {x\1+i/2-, ■ ■ ■ х?а+п ■ ■ ■ > ж^+1/2)Т- Аналогично одномерному случаю введем два типа переменных — консервативные фс = ф(хс), относящиеся к центрам расчетных ячеек, и потоковые сpia = cp(xia), tpia+1 = (p(xia-1-1), относящиеся к центрам "граней".

Дискретизируем уравнение (4) по схеме CABARET в ее двухслойной форме:

,п+1/2 ; п —1/2 N п п

а=1

а

На примере двумерного случая разностная схема (5), (6) выглядит так:

/ 71+1/2 /71—1/2 71 71 71 71

Уг+1/2, з + 1/2 ~ Уг+1/2, з + 1/2 ^г+1, 3 + 1/2 ~ ^г, 3 + 1/2 ^г+1/2,з + 1 ~ ^г+1/2,з _ »

| Сх . 'Су — и,

т Ах А у

п+1 _ 9 /П+1/2 п

П+1,3 + 1/2 — ' г+1/2,3 + 1/2 ^г,3 + 1/2!

п+1 _ 9 /П+1/2 _ п

Pi+l/2,j+l ~ ' г+1/2, j-\-l/2 ^¿+1/2,^

где ¿i = г, ¿2 = i; ''i = сж, С2 = су, h\ =Д ж, /12 =Д у. Расчетная ячейка в данном случае предствляет собой прямоугольник. Консервативная переменная относится к центру прямоугольника, потоковые переменные — к центрам его ребер.

4. Анализ устойчивости схемы CABARET для чисто гиперболического случая. Исследуем разностную схему (5), (6) на устойчивость по начальным данным, при этом будем полагать сеточные функции финитными, т. е. обращающимися в ноль вне некоторого ограниченного множества. Преобразуем (6) в (7):

= Ж+1/2 - Г С+1/2 = «V1! + Ю/2, tf = 2ф'Г1/2 - 4>V ^ I = К + 1 + ^а"1)/2-

(7)

Перепишем (5):

JV

фп +1/2 _ C-l/2 + £ Гв(¥,« +1<| _ ) = 0, (8)

а=1

где га = сат/ка.

Умножив (8) на фс+1^2 + ф™ получим

N

(с+1/2)2 - (Фс~1/2)2 + Егж+1/2+Фп~1/2тх+1в - о- (9)

а=1

Распишем фс+1^2 + фс_1^2 = {ф^+х + + (^+1 + У^"1)/2' пользуясь соотношениями (7),

и получим тождество

ф'2+1/2 +ф™-1'2 = 9 = 9*, (10)

где 9 = 0.5[(^+1 + у»?"1) + + V?«)]. = О-бК^+Л + ^+1) + К + уГ1)]-

Преобразуем —

2у>?«+1«~2у>?° = 0.5[(^+!1а + &а+1а) - + у,?-1) + ^в+1в - <р1 - ч>£Ъа + ^Г1] =

= о.5[(^+!1а + &а+1а) - + ^с1) + (^в+1в + у,?-1) - +

Воспользовавшись тождеством (10) и предыдущим равенством, получим

= 0.5 { + <рЪ + 1а) - (^ + ^ + Ыа + 1а + <Х) - « + 1а + О] =

= 0.25 {(¥>£+1а + &а+1а)2 ~ + <_1)2 + «+1а + <_1)2 " «++1в + О2} •

Распишем с помощью (7) (■г/>с+1,/2)2 — {фс~1^2)2'

п

(С+1/2)2 - (Фс~1/2)2 = Е [К + 1а + О + К + 1а + ^а"1)] /4ЛГ

и представим (8) в виде ^ N ^ N

а= 1

а=1 а=1

Е {К++1а + < + 1а)2 " К + + K + la + О' " K + la + О'} = 0>

а= 1

ИЛИ

п / . \ п

S ÍV - Га) + + Е r«(<PiSla + < + 1а)2 =

а=1 ^ ' а=1

п / 1 \ п £ ÍV " Га) + + Е + vV)2- (11)

п— 1 N / п— 1

а=1 а=1

Просуммировав равенство (11) по всем пространственным ячейкам, с учетом финитности сеточных функций получим

^п+1/2 _ ^п —1/2 (12)

где

п

П / . N П

En+l'2 = Е Е - ч («1-+ю2 + Е Е '«(«i- +

Slft а=1 ^ ' flh а=1

^-1/2 = £ ¿ (1 _ Гв) (V«+le + + £ ¿ Гв(^«+1в + -¿j:

о, л—1 \ / о, л—1

„п-1 \2

' alfíc + lc "I" V

fi/j а=1 7 ílti а=1

— положительно определенные квадратичные формы при условии

< 1/JV. (13)

Если снова использовать свойство (7), то (12) можно переписать

п / 1 \ п

ЕП+1/2 =ЕЕ4(^-Г") (^+1/2)2 + Е Е »-«K + la + < +

а=1 ^ ' fi/j а=1

п / 1 \ п д«-1/2 = ^ ^ 4 П (С-1/2)2 + ^ ^ Гв(< + 1в + ^ )2.

а=1 ^ ' П^ а=1

Выражения вида £'п+1/2, Еп~1!2 можно рассматривать в качестве нормы сеточных функций при выполнении условия (13) и, следовательно, равенство (12), согласно [12, 13], является достаточным условием устойчивости разностной схемы (5), (6).

5. Схема CABARET для многомерного уравнения конвекции^ диффузии. Перейдем теперь к многомерному уравнению конвекции-диффузии

гч N „ N

sr LL - V^ LL (ЛА\

dt+ 2^С«дха- 2' ^

а=1 а=1 "

где /X ^ 0.

Дискретизируем его по схеме CABARET:

/П+1/2 i п —1/2 JV / n _ п Ч n п

^—^— + Ес« ^ =^Е+1/2 + (i-Е«г1/2, (15)

Т а=1 а а=1 а=1

= 2С"1/2 - ^Г1

(16)

'ia + l — "Ve ' _

где 8 € [0,1], = (V'c+u - + ^c-ia)/h2a, причем с + 1а = н + 1/2,... ,ia + +3/2,... ,iN + 1/2,

с - 1а = ¿i + 1/2, ...,ia~ 1/2,..., iN + 1/2. Перепишем (15):

N , п п %

^„+1/2 _ фП — 1/2 + £ Га(ср1 + 1а ~<р?а)= /trlíjAa^ + 1/2 + (1 " í) Е 1/2 ) • i")

а=1 а=1 а=1

N

Введем, согласно [12, 13], обозначение Афс = — ^ Лсфс. Поступим с (17) аналогично тому, как

а=1

поступили с (8), умножим на фс+1^2 + фс 1^ и просуммируем по всем пространственным ячейкам П^: Еп+1/2 ^ Еп-1/2 = ¿^п+1/2^+1/2 +(!_£) ^п-1/2 ^п+1/2 + ^п-1/2^ ? (18)

N

где А = — Ка — как известно из [12, 13], положительно определенный, самосопряженный оператор, а=1

обладает свойством ограниченности по норме

Распишем (Афс ' , ф,

А <Д0 Е,

$1-1/2 „£»+1/2

N г

а=1 "

(19)

(а^"1/2, ф^1'2 + Фс~1/2) =1(а (С+1/2 + Фс'1'2) + А [ф'Г1'2 ~ Фс+1/2) , С+1/2 + Фс~1/2) =

= | (А (с+1/2 + С-1/2) ,С+1/2 + ФГ1'2) +1 (А$"1/2,$"1/2) -1 (А^1/2,^1/2

Аналогично

(а$+1/2,с+1/2 + фпс~1'2) = \{а {фпс+1/2 + фг1'2) + а (с+1/2 - Фс'1'2),С+1/2 + С"1/2) = = \ (а (с+1/2 + Фс~1/2) ,Фс+1/2 + ФГ1'2) + \ - \ (афг1/2,ФГ1/2) ■

Получим, что в правой части (18) выполнено

5 (Аф^1'2^1'2 + ф^-1'2) + (1 - (а^"1/2^С+1/2 + Фс~1/2) =

[а (фпс+1'2 + ФГ1'2), С+1/2 + С"1/2) + (1 - 25) (АФГ1/2,ФГ1/2) +

= 0.5

+(25-1) (А$+1/2,$+1/2

В итоге (18) можно переписать в виде

Еп+1/2 _ £«-1/2 = ^

п / 1 \ п

Еп+1/2 = £ £ 4 (1 - гв) (^с+1/2)2 + Е Е Г^га + <+1в)2 - 2(1 - 25) С+1/2) ,

П^ а=1 ^ ' П^ а=1

п / 1 \ п ^„-1/2 = ЕЕ4(^^) №-1/2)2 + Е Е + сЛ „)2 - 2(1 - 25) (агс-1/2,гс~1/2) ,

П^ а=1 ^ ' П^ а=1

= 2цт (а (^+1/2 + ФГ1/2), С+1/2 + С"1/2) ^ о,

или в операторной форме

Еп+т = ^¿4 (1 (С+1/2,С+1/2) +

Пь а=1 ^ '

п

+ Е Е »•«К+1« + <.+1 «)2 - 2(1 - 25) ,

П^ а=1

Еп~1'2 = ЕЕ4 (<Фс~1/2жг1/2)+

пк а=1 ^ '

п

+ Е Е + <Г+1а)2 - 2(! - 25) {АФГ1'2ж~1'2) ■

а= 1

Условие (20) есть достаточное условие устойчивости, если фигурирующие в ней квадратичные формы положительно определены, и их соответственно можно рассматривать как нормы сеточных функций.

Выражения £'п+1/2; Еп~1!2 — положительно определенные квадратичные формы при условии выполнения операторного неравенства

* ( 1 \

2 Е ( ы ~Га)Е ~ "^ (21)

а= 1 ^ '

Неравенство (21) будет заведомо выполнено при любом /х ^ 0, если га ^ 1/Ы и 8 € [0.5,1]. Для справедливости же неравенства (21) тогда, когда 8 < 0.5, необходимо выполнение неравенств

1 - 0, А (22)

ДГ ЛГ 4 ;

2 S (jV Г")

а= 1

Кроме того, используя свойство (19), получим

26) Л, EZA^A. (23)

2 S (JV — Г")

а=1

Неравенства (23) совместны при условии

/хг( 1 - 2<5) о ^ ^

2£ (*-r«)

Дп1^

а=1

Следовательно, неравенства (23) верны тогда, когда

1-гв^г£(1-2«5), (24)

где г^ = 0.5/хг Д^1 — параболическое число Куранта, при Д0= 4^ 1/Л,2. Условие (24) есть достаточное условие устойчивости схемы (15), (16), если 8 < 0.5.

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

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

1. Чубарова Н.Е., Горбаренко Е.В., Незваль Е.И. и др. Аэрозольные и радиационные характеристики атмосферы во время лесных и торфяных пожаров в 1972, 2002 и 2010 гг. в Подмосковье // Изв. РАН. Физика атмосферы и океана. 2011. 47. № 6. С. 790-800.

2. Еланский Н.Ф., Мохов И.И., Беликов И.Б. и др. Газовые примеси в атмосфере над Москвой летом 2010 г. // Изв. РАН. Физика атмосферы и океана. 2011. 47. № 6. С. 729-738.

3. Дианский H.A., Гусев A.B., Фомин В.В. Особенности распространения загрязнений в северозападной части Тихого океана // Изв. РАН. Физика атмосферы и океана. 2012. 48. № 2. С. 247-266.

4. Головизнин В.М., Самарский A.A. Разностная аппроксимация конвективного переноса с пространственным расщеплением временной производной // Матем. моделирование. 1998. 10. № 1. С. 86-100.

5. Головизнин В.М., Самарский A.A. Некоторые свойства разностной схемы Кабаре // Матем. моделирование. 1998. 10. № 1. С. 101-116.

6. Головизнин В.М., Карабасов С. А., Кобринекий И.М. Балансно-характеристические схемы с разделенными консервативными и потоковыми переменными // Матем. моделирование. 2003. 15. № 9. С. 29-48.

7. Головизнин В.М. Балансно-характеристический метод численного решения одномерных уравнений газовой динамики в эйлеровых переменных // Матем. моделирование. 2006. 18. № 11. С. 14-30.

8. Karabasov S.A., Goloviznin V. М. Compact accurately boundary adjusting high-resolution technique for fluid dynamics // J. Comput. Phys. 2009. 228. P. 7426-7451.

9. Goloviznin V.M., Hynes T.P., Karabasov S. A. CABARET finite-difference schemes for the one-dimensional Euler equations // Mathematical Modelling and Analysis. 2001. 6. P. 210-220.

10. Karabasov S. A., BerloffP.S., Goloviznin V. M. CABARET in the ocean gyres//Ocean Modelling. 2009. 30. P. 155-168.

11. Кострыкин С.В. Об одном варианте многомерного обобщения схемы "кабаре" // Матем. моделирование. 2010. 22. № 2. С. 69-82.

12. Самарский А. А., Гулин А. В. Устойчивость разностных схем. М.: Наука, 1973.

13. Самарский А. А. Введение в численные методы. М.: Наука, 1987.

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

50

ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2013. № 4

SUFFICIENT CONDITIONS FOR STABILITY OF THE KABARE-APPROXIMATION OF THE MULTIDIMENSIONAL CONVECTION-DIFFUSION EQUATIONS ONTO ORTHOGONAL CALCULATION GRIDS

Berezin B. I. Blagodatskikh D. V.

The sufficient conditions for stability of KABARE scheme including the case of dominating convection are obtained by utilizing the positive definite feature of mesh quadratic forms.

Keywords: convection-diffusion, dominating convection, numerical methods, KABARE scheme, stability condition.