CC BY
73
Вычислительные технологии
Том 9, № 6, 2004
ПОСТРОЕНИЕ МОНОТОННЫХ СХЕМ НА ОСНОВЕ МЕТОДА ДИФФЕРЕНЦИАЛЬНОГО
ПРИБЛИЖЕНИЯ*
Ю.И. Шокин, Ю.В. Сергеева, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: shokin@ict.nsc.ru, yulya@gorodok.net, khak@ict.nsc.ru
The new approach based on analysis of a differential approximation of a scheme for constructing of monotonic difference schemes of the second order approximation is considered. Monotonic scheme for the nonlinear transport equation is given. The generalization for the shallow water equations is carried out.
Введение
В настоящее время ТУБ-схемы и их многочисленные модификации применяются при решении многих задач с разрывными решениями. Причина столь большой популярности этих методов заключается в том, что они дают неосциллирующие профили решения, высокую разрешимость в области разрывов и сохраняют высокую точность в областях гладкости решения.
Современные ТУБ-схемы высокого порядка аппроксимации основаны на тех или иных способах восстановления (реконструкции) значений функций на гранях ячеек по их значениям в центрах соседних ячеек. При этом шаблон схемы является переменным и зависит от поведения численного решения. Алгоритмы реконструкции основываются на использовании специальных ограничителей потоков [1, 2], которые строятся так, чтобы схема с ограничителями обладала ТУБ-свойством [3] (свойством невозрастания полной вариации численного решения) и, как следствие, сохраняла монотонность численного решения.
В настоящей работе монотонизацию схем второго порядка аппроксимации предлагается выполнять не на основе построения ограничителей потоков, а путем анализа дифференциального приближения схемы. Этот подход к построению монотонных схем демонстрируется на явной схеме предиктор-корректор [4] для нелинейного скалярного уравнения и для системы уравнений мелкой воды с одной пространственной переменной.
* Работа выполнена при поддержке Программы интеграционных исследований СО РАН (проект 2003-5) и Программы поддержки ведущих научных школ РФ (грант НШ-2314.2003.1).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
1. Монотонная схема предиктор-корректор для скалярного уравнения
Рассмотрим явную схему предиктор-корректор для скалярного уравнения
и + [/(«)]* = 0.
На шаге предиктор
f* _ _( fn I fn
Jj 2(Jj-1/2 + J j+1/2
T
) f
_ + an Jj+1/2 + aj
nn - fj-
j-l/2
h
0
(1)
(2)
вычисляются потоки /*, отнесенные к целым узлам х^ (границам ячеек) равномерной сетки с шагом Н. В уравнении (2) т* = 0.5т(1 + в"), т — шаг по времени, в" — параметр схемы, вообще говоря, меняющийся от узла к узлу и от одного временного слоя к другому,
/Т+1/2 = / (^+1/2),
nn - f j-
j+1/2
u
j+1/2 Uj—1/2
j 1/2 при иП+1/2 = иП-
j—1/2'
а(иП+1/2) пРи и"+1/2 = и"-1/2'
а(и) = /и(и). Уравнение (2) есть результат аппроксимации уравнения
/г + а(и)/х = 0, (3)
которое получается после умножения (1) на функцию а (и). На шаге корректор
un+1 _ un Uj+1/2 U j+1/2 + Jj+1
** J j+1 - J j
T
h
0
(4)
находятся искомые величины un+1/2, определенные в полуцелых узлах xj+1/2 = Xj + h/2 (центрах ячеек).
При в = 0 выписанная схема совпадает со схемой Лакса — Вендроффа. Если в = O(h), то схема (2), (4) аппроксимирует уравнение (1) со вторым порядком относительно т и h. Необходимое условие устойчивости, полученное для случая а(и) = const, требует, в частности, чтобы параметр в принимал неотрицательные значения. Поэтому далее предполагается выполненным следующее ограничение:
j > 0.
Схему (2), (4) можно записать в виде
un+1 _ un U j+1/2 Uj+1/2
T
+ C- ,jU>n,j _ C+,j+1U:
lX,j + 1 = 0,
где
и
Un+1/2 Un-1/2
h
C
±,j
ж(1 + в^а1^)2 т oOX, , ж = т/h = const.
n
а
j
1
2
Выполнение неравенств
C-j > 0, C+j > 0, C-j + C+j < 1 (6)
ж
достаточно [3] для того, чтобы рассматриваемая схема была TVD-схемой. Поэтому при выполнении условий
—< |an|ж < , 1 (7)
1 + в™ - 1 j 1 - y^TTef v ;
и в = 0(h) рассматриваемая схема будет схемой второго порядка аппроксимации, сохраняющей монотонность численного решения.
Отметим, что для произвольной неотрицательной функции e(x,t) такой, что в = 0(h), добиться выполнения условий (7) практически невозможно. Даже если a(u) = const и условия (7) выполнены и, следовательно, рассмотренная схема сохраняет монотонность, она и в этом случае имеет существенный недостаток, заключающийся в том, что на мелких сетках для сохранения второго порядка аппроксимации и сохранения монотонности расчет приходится вести с числами Куранта, близкими к единице, причем чем мельче сетка, тем ближе к единице должно быть число Куранта.
Далее мы будем предполагать, что сеточная функция в™ зависит от численного решения на n-м временном слое. Оказывается, что для такого переменного параметра схема предиктор-корректор может сохранять монотонность решения уже при любых числах Куранта, подчиняющихся условию
|an|ж < 1. (8)
Для выбора параметра в будем использовать метод дифференциального приближения [5, 6]. Предполагая, что функция в(х) и все ее производные имеют порядок 0(h), получим следующее выражение для первого дифференциального приближения (ПДП) схемы (2), (4):
т h2
ut + fx = ^(вa2Ux)x + — [а(а2ж2 - 1)uJxx . (9)
В подобластях, где возникает угроза появления осцилляций численного решения, необходимо изменить дисперсию разностной схемы. При в = 0(h) основной вклад в дисперсию разностной схемы будет вносить второе слагаемое правой части (9). Подберем функцию в так, чтобы первое слагаемое правой части ПДП частично или полностью компенсировало дисперсионный член
ahh2
— (а2ж2 - 1) Uxxx (10)
либо давало такой вклад в ПДП, который приводит к смене знака коэффициента при третьей производной uxxx. Пусть для определенности а > 0 и аж < 1. Тогда в (10) коэффициент при uxxx будет иметь отрицательный знак. Если, например, положить
в(х) = h[а(1 - ^^x
3а2ж^
x
то дисперсионный член (10) полностью компенсируется первым слагаемым правой части ПДП (9). Если функцию 0(ж) выбрать, например, в виде
*(*) = Н К1 - ажК]*, (11)
а
то с учетом первого слагаемого правой части (9) дисперсионный член ПДП
ah2
— (1 - аж) (2 - аж) Uxxx 6
будет иметь уже положительный коэффициент. Таким образом, выбирая ту или иную функцию в(х), можно управлять дисперсией разностной схемы.
Для определенности далее мы будем рассматривать функцию в(х) в виде (11) и по-прежнему предполагать, что a > 0. Ясно, что менять дисперсию схемы во всей области не имеет смысла. Поэтому формула (11) требует дальнейшего уточнения. Учтем теперь условие неотрицательности функции в(х). Тогда в тех подобластях, где производные [а(1 — аж)их]х и ux имеют разные знаки, можно принять в = 0 (схема Лакса — Вендроффа) либо положить в = Ch, C = const > 0 (схема предиктор-корректор с малой аппроксима-ционной вязкостью). Рассмотрим далее первый случай. Тогда вместо (11) получим
a< \ Л [а(1 — аж)их]х Л в(х) = max h——---——, 0 .
\ а2жих у
Кроме того, необходимо позаботиться об ограниченности функции в(х) при h ^ 0. Можно поступить, например, так. Если
[а(1 - аж)их]х > 1
a2 ^Ux
h'
то в качестве значения функции 0(х) использовать число, которое для данных ж и а > 0 удовлетворяет неравенствам (7). К примеру, положим
в = -L -1.
аж
Тогда указанный способ выбора функции приводит к формуле
0(я) = min
1 - аж / [а(1 - аж)мх]х п max ( h——2-— , 0
аж \ а жи
Вид сеточной функции зависит от выбора аппроксимационной формулы для вычисления производной [а(1 — аж)их]х. При ап > 0 заменим ее, например, следующим конечно-разностным выражением (противопоточная аппроксимация):
а™ (1 — а™ж) их. — а™—1 (1 — а™-1Ж «х.—1
Н '
Поступая аналогичным образом при а < 0, приходим к одной из многих возможных формул для вычисления сеточной функции в случае произвольного знака функции а(м):
0 при а. = 0 или |¿7| < |<.-—8| , Мх.мх.—8 > 0,
а.1 — жа2) Их. — (|а.—^ — жа2-^ Мх,.—8 -—-^--—-- при а. = 0
жа. (12)
и |<.1 >
> 0,
в;
|aj | - жа2
при aj = 0 и uxjuX)j--s < 0,
жа2
где в = sgn а.,
& = ^ГТ(1 - |а. И
2
и опущен верхний индекс п у сеточных функций ап и «П..
Покажем, что при условии (8) схема предиктор-корректор (2), (4) с сеточной функцией в, заданной формулой (12), будет сохранять монотонность численного решения. Для этого преобразуем выражение для потока (индекс п опущен):
с * _ 1 ( г I г __2/
Л* = 2 (/.7+1/2 + Л -1/2 - та2(1 + в. К,,") =
1
2 (У7+1/2 + /1/2 - |а. |Них. + 2Н& - а2тв.. (13)
Числа и имеют одинаковые знаки (ввиду условия (8)), поэтому выражение для потока можно представить в виде
/* = 2 (Л+1/2 + У7 —1/2 - |а. |Них,. + 2И-5/2) ,
где
# = Г )тт(|&|, |(&.+1|) при &&.+1 > О,
Л+1/2 \ 0 при &< 0.
Поскольку числа 1/2 и ^.+1/2 не могут иметь разные знаки,
Ь+1/2 - 1/2| < тах(|#._ 1/2 |, |^.+1/2|) < ^ |
и
Н|^| < Ж = ^ (1 -|а,|в) <1|1, (14)
где
^.+1/2 - #7 —1/2 ,
-г- при Их,, = 0,
7. = ^ Них,.
0 при = 0.
Из оценки (14) следует, что числа а. и ^М = а. + Н^. одинакового знака, поэтому имеет место равенство
1
Л* =0 У.+1/2 + Л —1/2 + Н (&'_1/2 + #.+1/2 - |
Используя это выражение для потока, схему (4) можем записать в виде (5)
.+1/2 И.+1/2 + | + У? ип - |^.+1| - У?+1 ип
Т +2 2 Их,.+1
при этом С_. > 0, С+. > 0, поэтому если выполнено условие
«К | < 1, (15)
то она будет ТУБ-схемой и будет сохранять монотонность решения.
Воспользовавшись при мх. = 0 оценкой (14), получим
ж|.|< ж|а.1 + 2(1 — |а.|ж)) '
Следовательно, выполнение неравенства
ж|а.1 + 2(1 — |а.|ж)) < 1 (16)
достаточно для выполнения условия (15). Легко проверить, что при условии (8) неравенство (16) выполняется.
Отметим, что если функция в задается по формуле (12), то ТУО-схема (2), (4) совпадает с известной схемой Хартена [3] с ограничителем типа тттс^
2. Схема предиктор-корректор для уравнений мелкой воды с одной пространственной переменной
Применим рассмотренный выше способ управления дисперсией к схеме предиктор-корректор, аппроксимирующей на равномерной сетке с шагом Н уравнения мелкой воды:
V* + ¥х = О. (17)
Здесь
у=(Н„); Е = (); О= (нН
м(ж, ¿) — скорость жидкости; Н(ж, ¿) = п(ж, ¿) + Н(ж) — ее полная глубина; п(ж, ¿) — возвышение свободной поверхности над невозмущенным уровнем г = 0 и г = —Н(ж) — функция, задающая дно бассейна.
На шаге корректор находятся величины Уп++11/2 в полуцелых узлах на (п + 1)-м слое по времени. Для этого используется аппроксимация дивергентной системы (17):
уп+1 — V™ т?* т?*
* .+1/2 У.+1/2 — _ п* Пйч
- + -Н- = О.+1/2' (18)
Т
Потоки Е* вычисляются на шаге предиктор. Для этого аппроксимируется уравнение для вектора потоков (аналог уравнения (3))
Е* + АЕх = АО, (19)
полученное в результате умножения уравнения (17) на матрицу Якоби А = дЕ/дV. Пусть Л&, (к = 1, 2) — собственные значения матрицы А, Л — диагональная матрица с элементами Л1, А2 на диагонали, И(V) — правые собственные векторы, соответствующие собственным значениям Л&, Я^) — матрица, столбцами которой являются эти векторы, Ь = Я-1. Тогда А2 = ЯЛ2Ь и уравнение (19) можно переписать так:
ЬЕ4 + Л2ЬУх = ЬАО'
Аппроксимируем его следующим разностным уравнением:
р* — 1 (Р"- , + Р" , ) (Д-1Ь)га 7 2 ( з+1/2 + з-1/2) + (Л2Р)га = .
з т / 2 3 3
Здесь Р" = V" з; — диагональная матрица с элементами 1+0Пз-, 1+$2 з на диагонали;
з > 0, = 0(Л,). Отсюда получаем окончательную формулу для вычисления вектора потоков:
F = 2
FJ+i/2 + F"-i/2 - т(ЯЛ2ВР)" + т (AG)n' . (20)
Последняя формула аналогична формуле (13) для потока f *, поэтому представляется возможным и при вычислении функций применить аналоги формулы (12), использованной при решении скалярного уравнения:
^fc j
0 при Afc j = 0 или |<7fcj | < |gfc j_s| , Pfc jPfc j-s > 0
- fffc j-'fc , j P , j A , j1 - «Afc j
2j j при Afcj = 0 и | > |gk,j-s| , PkjPkj-s > 0,
^Afc,i pfc'i (21)
при Afc'j = 0 и Pfc'jPfc'j-s < 0
где 5 = Л^з; к =1, 2, Р1 и Р2 — компоненты вектора Р; индекс п опущен и
9к,з = ^ (1 - 1Л^,з И •
Расчеты тестовых задач с разрывными решениями подтвердили, что схема (20), (18), (21) сохраняет монотонность профилей сеточных функций.
Заключение
На примере явной схемы предиктор-корректор описан новый подход для конструирования монотонных нелинейных разностных схем второго порядка аппроксимации, основанный на исследовании дифференциального приближения схемы. Указана одна из возможных формул задания аппроксимационной вязкости, при использовании которой построенная схема совпадает с известной ТУБ-схемой Хартена. Известные и широко используемые на практике ТУБ-схемы с другими ограничителями потоков также могут быть получены на основе предлагаемого подхода.
Список литературы
[1] Sweby P.K. High resolution schemes using flux limiters for hyperbolic conservation laws // SIAM J. Numer. Anal. 1984. Vol. 21, N 5. P. 995-1011.
[2] Куликовский А.Г., ПогорЕлов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001.
[3] Harten A. High resolution schemes for hyperbolic conservation laws //J. Comput. Phys. 1983. Vol. 49, N 3. P. 357-393.
[4] ЧИСЛЕННОЕ моделирование течений жидкости с поверхностными волнами / Г.С. Хакимзянов, Ю.И. Шокин, В.Б. Барахнин, Н.Ю. Шокина. Новосибирск: Изд-во СО РАН, 2001.
[5] Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, Сиб. отд-ние, 1985.
[6] Шокин Ю.И., ХАКИМЗЯНОВ Г.С. Введение в метод дифференциального приближения: Учеб. пособие. Новосибирск: НГУ, 1997.
Поступила в редакцию 21 октября 2004 г.
