CC BY
156
Вычислительные технологии
Том 17, № 2, 2012
Некоторые замечания о схемах,
>к
сохраняющих монотонность численного решения*
Г. С. Хлкимзянов, Н.Ю. ШокинА, Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]
На простых примерах схем для уравнения переноса с постоянным коэффициентом обсуждаются общие проблемы, возникающие при построении адаптивных сеток и конструировании дивергентных разностных схем, сохраняющих монотонность численного решения. В частности, приведены пример схемы с положительными коэффициентами, которая не сохраняет монотонность численного решения, и пример ХУБ-схемы, увеличивающей количество экстремумов.
Ключевые слова: численное моделирование, конечно-разностная схема, уравнение переноса, предиктор-корректор, монотонная схема, дивергентная схема, адаптивная сетка, метод эквираспределения, результаты расчётов.
Введение
В статье С.К. Годунова [1] показано, что для уравнения переноса
не существует схем с постоянными коэффициентами, сохраняющих монотонность численного решения и имеющих второй порядок аппроксимации. Например, при использовании схемы Лакса — Вендроффа заданная при £ = 0 "ступенька" искажается на следующих слоях по времени нефизичными осцилляциями. Но схему Лакса — Вендроффа можно модифицировать так, чтобы она стала обладать ТУВ-свойством [2], а значит, сохраняла монотонность численного решения. При этом коэффициенты монотонизиро-ванной схемы уже не будут постоянными, они могут зависеть от численного решения и его разностных производных, т. е. модифицированная схема станет нелинейной. Стандартный способ монотонизации схемы Лакса — Вендроффа заключается в переключении на противопоточную схему в тех узлах сетки, где возникает угроза появления ос-цилляций численного решения. В ТУВ-схемах переключение производится с помощью ограничителей [3]. Современные ТУВ-схемы высокого порядка аппроксимации основаны на тех или иных способах реконструкции значений функций на границах ячеек по их значениям в центрах соседних ячеек [4]. Отметим, что идея локального переключения с одной схемы на другую с целью сохранения свойства монотонности численного решения уравнения (1) использовалась намного раньше (см., например, [5]) Хартена [2].
* Работа выполнена при финансовой поддержке РФФИ (гранты № 10-05-91052-НЦНИ и 12-01-00721) и Совета по грантам Президента РФ по государственной поддержке ведущих научных школ РФ (грант № НШ-6293.2012.9).
ди ди
——+ а— = 0, а = const, dt дх
(1)
В настоящей работе рассматривается другой подход к построению монотонных разностных схем, основанный на исследовании их дифференциальных приближений. При этом основное внимание уделено тем вопросам, которые не были затронуты в публикациях [6-8], в частности, проблемам, связанным со свойствами монотонности и дивер-гентности схем на подвижных неравномерных сетках, с построением сеток, адаптирующихся к разрывным решениям, с разрешимостью уравнений для сетки.
Многие схемы, сохраняющие монотонность численного решения, дают на разрывных решениях осциллирующие профили разностных производных [9], что может негативно сказаться при использовании метода адаптивных сеток: если управляющая функция, регулирующая расстановку узлов, зависит от разностных производных, а они осциллируют, то адаптивная сетка будет иметь чередования длинных и коротких ячеек, что приведет к потере точности численного решения. Осциллирование разностных производных численного решения может быть вызвано, в частности, и "нефизичным" ростом количества экстремумов решения (даже при использовании ТУВ-схем) при переходе с одного шага по времени на другой. Для устранения проблем, связанных с резкими изменениями значений управляющей функции, предлагается использовать процедуру неявного сглаживания данной функции, после применения которой в окрестность разрыва попадает большее количество узлов и происходит плавное изменение длин соседних ячеек, что способствует лучшему воспроизведению решений с разрывами.
1. Монотонизация явных схем на равномерных сетках
Продемонстрируем идею нового подхода к монотонизации разностных схем на явной схеме предиктор-корректор [10, 11] для уравнения (1):
и* — 7/П 1,п+1 ип и* — и*
3+1/2 и3+1/2 , ^ п п и 3 ~ и 3 , 3+1/2 3 —1/2 _п --- + аих 3+1/2 — о, - + а--- — и,
Т3+1/2
т
К
где т — шаг по времени, К — шаг равномерной сетки с узлами ж3 — ]К, вспомогательные величины и*+1/2, определенные в полуцелых узлах £3+1/2 — ж3 + К/2, относятся к моменту времени £ — £п + т*+1/2, £п — пт,
и
3+1/2
иП+1 + ип п
2 , их,3+1/2
и3+1 и3 т
К
т3*+1/2 — 2(1 + ^п+1/2),
9 — параметр, меняющийся, вообще говоря, от узла к узлу и от одного временного слоя к другому. На шаге "корректор" вычисляются искомые величины ип+1, определенные в целых узлах ж3.
1.1. Каноническая форма двухслойных явных схем
Одношаговый вариант схемы (2) имеет следующий вид:
ип+1 ип ип _ ип
из - из + аи3_+1/2 и3 —1/2
т
К
а2т ~2Л
(1 + 9)их) - (1 + 9)их
х +1/2 х —1/2
0. (3)
Схема (3) является в определенном смысле канонической формой явных двухслойных схем для уравнения (1), поскольку любая такая схема может быть записана в виде (3)
п
п
при соответствующем значении в. Например, при в = 0 схема (3) совпадает со схемой Лакса — Вендроффа
п+1 п ип ип 2
— " з "3+1/2 — "^—1 /2 а2т / п п ч
Т + а Н ^3+1/2 — = 0' (4)
при
в = во = С — 1 > 0, (5)
где С — число Куранта,
Т
С = |а|ж < 1, ж = —,
Н
получается противопоточная схема первого порядка аппроксимации
— — "п + а + |а|„п + а — |а|„п
Т +2 "х,3-1/2 + 2 "х,3+1/2
1
при в = в^ = —- — 1 > 0 — схема Лакса С 2
1
«п+1 — о ("п+1 + "п-1) "п , -ип ,
3 2 3+1 3 + из+1/2 "3-1/2 = 0
Т Н
при в = -1 — абсолютно неустойчивая схема с центральной разностью
"п+1 "п "п "п
из — из , и3+1/2 3-1/2
+ а ^17 ,—^^ = 0.
ТН
Заметим, что в канонической форме (3) записываются не только трехточечные схемы с симметричным шаблоном. Например, при а > 0 и
вп = в 1 х,з-1/2
вз+1/2 = во I 1 — "п-
их,3+1/2.
(3) переходит в противопоточную схему второго порядка аппроксимации на несимметричном шаблоне
п+1 п 3ип _ ип 2 „,п 0„,п I „,п
из — и , 3 х,3—1/2 их,з-3/2 а2Т " — 2и,--1 + и,--2
—~ + а-2-= "2--Н2-. (10)
Выбирая то или иное значение параметра в, можно получить явные схемы с различными свойствами: первого или второго порядка аппроксимации, абсолютно или условно аппроксимирующие, условно устойчивые или абсолютно неустойчивые, сохраняющие или не сохраняющие монотонность численного решения. Остановимся здесь только на вопросе о влиянии параметра в на свойство монотонности двухслойных явных схем.
1.2. Постоянный параметр в
Сначала рассмотрим случай постоянного параметра в = const. Тогда схема (3) может быть записана в виде
un+1 = b-iun-1 + boun + biun+1 (11)
с постоянными коэффициентами 6а (а = -1, 0,1)
6-1 = а2^2(1+ в) + , bo = 1 - aW(1 + в), bi = а2^2(1+ в) - . (12) Необходимыми условиями устойчивости схемы являются ограничения
в > 0, С < , 1 , (13)
> ' < VTT^
которые при условии (6) могут быть записаны в виде
0 < в < eL. (14)
Для сохранения монотонности численного решения необходима и достаточна [1] неотрицательность коэффициентов 6а, что эквивалентно выполнению неравенств
во < в < eL. (15)
При в = const = 0 все схемы (3) имеют первый порядок аппроксимации на гладких решениях, а если параметр в принадлежит интервалу монотонности (15), то они сохраняют монотонность численного решения. Например, в схеме Лакса (в = в^) коэффициенты
, 1 + аж , 1 — аж
6-1 =-, 61 =-
1 2 ' 1 2
при условии (6) положительны, поэтому она сохраняет монотонность численного решения.
На рис. 1, a показаны границы интервала монотонности схемы (3) и точки, соответствующие числу Куранта C = 0.8 и значениям в = 0 (1), в = в0 (2), в = в^ (3). На рис. 1, б приведены профили численного решения в момент времени t =10 для а = 1, длины области l = 30, h = l/N, N =150, и кусочно-постоянной начальной функции
у \ Г 1 при ж < ж*, _ , ...
Uo(x) = S п " 0 < ж* < l, (16)
w [0 при x > ж*, v 7
ж* = 10. Поскольку схема (3) не позволяет вычислить приближенное решение uNv+1 в самом правом узле сетки жv = l, необходимо для этого узла использовать дополнительное разностное уравнение, например схему второго порядка (10).
Из рис. 1, б видно, что схема Лакса — Вендроффа, соответствующая параметру в = 0, не принадлежащему интервалу монотонности (15), дает осциллирующий профиль 1. Противопоточная схема и схема Лакса, которым отвечают точки, лежащие на нижней в = в0 и верхней в = в^ границах интервала монотонности, дают соответственно монотонные профили 2 и 3. При этом численное решение, полученное по схеме Лакса, имеет осциллирующие профили разностных производных, поскольку данная схема не является "сильно" монотонной [9].
а б
Рис. 1. Область монотонности схемы (3) (а); профили численного решения в момент времени £ = 10 при разных значениях параметра в (б)
Заметим, что если схема не сохраняет монотонность численного решения, то это не означает, что она все монотонные функции переводит в немонотонные. Возьмем, например, схему Лакса — Вендроффа (4), которая не является схемой, сохраняющей монотонность численного решения, и рассмотрим множество монотонно возрастающих на n-м слое по времени функций ,
<J+1/2 > 0, j
:i7)
удовлетворяющих при всех j условию
j '' 0,
i8)
где
S Sgn a,
Представим схему (4) в следующем виде:
U"j+1/2 U"j-1/2 h '
ат / \ а2т2
ип = ип 2 \ип,з'+1/2 + иж,з-1/^ + ипж-.
Записав аналогичное равенство для узла х-+1 и вычтя первое равенство из второго, получим выражение
|а|т
Un,j+1/2 = Un,j+1/2 + 2 (sC ^ S%:c,j+1 (sC + 1) S%:c,j
из которого в силу предположения (18) и ограничения на число Куранта (6) следует, что свойство монотонного возрастания (17) сохранится и на (n + 1)-м слое по времени для всех функций, удовлетворяющих условиям (18). Нетрудно проверить, что при a > 0 схема Лакса — Вендроффа сохраняет на (n + 1)-м слое по времени монотонность всех
монотонно возрастающих функций, выпуклых вверх (мПх < 0), и монотонно убывающих функций, выпуклых вниз (мПх > 0), а при а < 0 — наоборот: монотонно возрастающих функций — выпуклых вниз, монотонно убывающих — выпуклых вверх.
Следующее замечание касается связи дисперсии схемы со свойством сохранения монотонности численного решения. Поскольку причиной появления осцилляций численного решения является дисперсия разностной схемы, то имеет смысл рассмотреть вопрос о выборе такого параметра в, при котором дисперсия схемы будет по возможности минимальной. С этой целью рассмотрим второе дифференциальное приближение схемы (3)
дм дм 52м 53м
"тгт + а— = с^т-^ + сз~ з, (19)
от дх дх2 дх3
где
С2 = ^в, сз = ^ (а2ж2(3в + 1) - 1) . (20)
Отсутствие дисперсии в решении второго дифференциального приближения возможно только при нулевом значении коэффициента с3 [12], т. е. при
в
в = л* = у. (21)
Назовем схему (3) с параметром в = в* "бездисперсионной". Она имеет следующий вид:
м?+1 - 1 (м?+1 + 4мп + мп ,) мп _ мп о
' 6 ' + ^». = 0. (22)
т К 3
На рис. 1, б профиль решения, полученного по схеме (22), соответствующий точке 4 на рис. 1, а, изображен линией 4. Видно, что осцилляции отсутствуют. Тем не менее это пример схемы, которая, несмотря на отсутствие дисперсии во втором дифференциальном приближении, не сохраняет монотонность численного решения. Нетрудно проверить, что при а > 0 коэффициенты схемы (22) выражаются через число Куранта по формулам
2С2 + 3С + 1 , 2С2 - 3С + 1 , 2(1 - С2)
о_ 1 =-, о1 =-, о0 =-.
1 6 ' 1 6 ' 0 3
Предположим, что для числа Куранта выполнены неравенства
0.5 < С < 1. (23)
Тогда 61 > 0, Ь0 > 0, но 6-1 < 0, откуда следует [1], что "бездисперсионная" схема не обладает свойством сохранения монотонности численного решения.
Отсутствие осцилляций на графике численного решения 4 (см. рис. 1, б) объясняется тем, что схема (22), как любая схема первого порядка аппроксимации, обладает значительной диссипацией (однако меньшей, чем противопоточная схема), поэтому возникающие в начале расчета немонотонности численного решения быстро подавляются благодаря действию диссипативных механизмов схемы.
Приведем простейший пример возникновения осцилляций в численном решении при использовании "бездисперсионной" схемы. Возьмем монотонно возрастающую функцию
мп Г 0 если з < ° (24)
м =1 1, если з > 0, ( )
и вычислим значения решения на (п + 1)-м временном слое:
{0, если з < -2,
6Ь если 3 = -1, (25)
00 + о1, если з = 0, 47
1 , если з > 1 .
Предположим, что а > 0 и для числа Куранта выполнены неравенства (23). Тогда
=(С - 1/2)(С - Р < 0, 0 <60 + =(1 - ^ + 5/2) < 1
33
Отсюда следует, что на (п + 1)-м слое решение стало немонотонным с одним строгим локальным минимумом в узле с номером з = - 1 .
Таким образом, даже при отсутствии дисперсии в решении второго дифференциального приближения схема может не сохранять монотонность численного решения. Это связано с тем, что дисперсионные члены могут содержаться в дифференциальных приближениях более высокого порядка.
1.3. "Квазипостоянный" параметр в
Если параметр в является "квазипостоянным", т. е. не зависит от х и но зависит от шагов К и т так, что в = О(К) (или в = О(т)), то коэффициенты схемы (3) будут зависеть только от шагов сетки. Тогда при выполнении условий (15) и в = О(К) > 0 порядок аппроксимации схемы (3) будет второй, и она будет сохранять монотонность численного решения [6]. Тем не менее в данном случае схема имеет существенный недостаток, заключающийся в том, что при измельчении сетки для сохранения второго порядка аппроксимации (чтобы при К ^ 0 в ^ 0) и монотонности расчет приходится вести с числами С, близкими к единице, при этом чем мельче сетка, тем ближе к единице должно быть число Куранта, поскольку из неравенства (15) следуют ограничения
1 <С< 1
1 + в - - v/ГT0'
Ясно, что при решении практических задач такое жесткое условие выполнить не удастся, поскольку в данном случае расчет обычно производится с такими числами Куранта, которые с некоторым запасом удовлетворяют условию устойчивости.
1.4. Переменный параметр в
В работах [6, 7] показано, что при специальном выборе переменного неотрицательного параметра в = О (К) схема второго порядка аппроксимации (3) может сохранять монотонность решения при любых числах Куранта, подчиняющихся условию (6). Для выбора в использовалось первое дифференциальное приближение (п.д.п.) схемы (3)
5м дм а2т д Лдм\ аК2 , 2 2 ^ д3м д£ дх 2 дх V дх / 6 ^ 7 дх3'
и параметр в подбирался так, чтобы диссипативный член п.д.п. либо частично или полностью компенсировал дисперсионный член, либо давал такой вклад в п.д.п., который
приводит к смене знака коэффициента при третьей производной. В результате была получена [7] следующая формула:
0 пРи K,j+1/2l ^ Kj+1/2-s| и ^+1/2^+1/2-* > 0
j+1/2 Н ~9 ^ - £+1/2) при |un,j+1/2| > |un,j+1/2-s| и «^+1/2^+1/2- > 0 (27)
0 при Ua;)j+1/2Un,j+1/2-s < 0
ГДе 0 = COnst > °> S = Sgn a> j+1/2 = Ua;,j+1/2—s/ua;,j+1/2.
Сделаем некоторые замечания об использовании переменного параметра 0. При использовании формулы (27) схема (3) становится схемой с переменными коэффициентами. На следующем примере покажем, что для схем с переменными коэффициентами условие неотрицательности коэффициентов не является достаточным для сохранения монотонности численного решения: коэффициенты могут быть положительными, но схема не будет сохранять монотонность численного решения. Пусть решается задача Коши для уравнения
ut + a(x)ux = 0,
где a(x) — строго возрастающая положительная ограниченная функция: 0 < a(x) < 1 и a' > 0. Рассмотрим для решения этой задачи схему с переменными коэффициентами
(«+1 + un-1)
«Г1 -- j + Л un+1 - un-1
т + = 0' (28)
где а. = а (ж.), ж. — узел равномерной сетки. Выписанная схема является аналогом схемы Лакса (8), которая сохраняет монотонность численного решения. Будем считать, что для любого ] выполнено условие
жа. < 1, (29)
гарантирующее устойчивость схемы (28) в равномерной норме по начальным данным. Запишем схему (28) в виде
«п+1 = + Ьо. + &1. «п+1, (30)
1 + жа. 1 — жа.
= —^ > о, Ьо,. = о, = —^ > 0.
Здесь коэффициенты снабжены дополнительным индексом ], поскольку они являются переменными и изменяются при переходе от одного узла к другому.
В силу условия (29) коэффициенты положительны, однако схема (28) не сохраняет монотонность численного решения. Так, используя монотонно возрастающую функцию (24), получим, что на (п + 1)-м слое по времени
(0, если ] < -2,
Ь1,_1, если ] = -1,
61)0, если ] = 0,
1, если ] > 1.
где
Но 61,-1 > 61,0, поэтому сеточная функция не является монотонно возрастающей.
Приведенный пример показывает, что для схем с переменными коэффициентами необходимы другие признаки монотонности, поскольку признак монотонности 6а > 0, используемый для схем с постоянными коэффициентами, в данном случае непригоден. Для схем с переменными коэффициентами справедлива следующая теорема [7]. Теорема 1. Пусть коэффициенты разностной схемы (30) удовлетворяют в каждом узле х равенству
6-1,3 + 60,3 + 61,3 = 1.
Тогда выполнение при всех ] условий
6±1,з > 0, 6-1,3 + 61,3-1 < 1
необходимо и достаточно для того, чтобы схема (30) с переменными коэффициентами сохраняла монотонность численного решения.
Ранее на примере "бездисперсионной" схемы было показано, что анализ дифференциального приближения не дает полной информации о свойствах схемы и при отсутствии дисперсии в решении второго дифференциального приближения возможно нарушение монотонности численного решения. Поскольку выбор формулы (27) для переменного параметра в также был основан на анализе дифференциального приближения, свойство сохранения монотонности требует строгого доказательства. Для этого можно использовать утверждение [8], следующее из теоремы 1. Теорема 2. Выполнение условий (6) и
во < в < 3 вь
достаточно для того, чтобы схема (3) с переменным параметром (27) сохраняла монотонность численного решения.
Отметим, что при использовании схемного параметра (27) с постоянной в = в0 схема предиктор-корректор (3) переходит в известную ТУВ-схему Хартена [2]. Все известные и широко используемые на практике ТУВ-схемы также могут быть получены на основе анализа п.д.п. схемы (3) как на равномерных, так и на подвижных неравномерных сетках.
2. О некоторых свойствах схем на подвижных сетках
Точность численного решения задачи может заметно возрасти при использовании адаптивных сеток, однако при конструировании схем на адаптивных сетках возникает много проблем. Некоторые из последних будут рассмотрены в настоящем разделе.
2.1. Схема предиктор-корректор на неравномерной подвижной сетке
Схема предиктор-корректор, аппроксимирующая уравнение (1) на неравномерной подвижной сетке, имеет вид
31/2 - 3+1/2 + / а у Ч» =0 (^)Г1 - , (^У*)з+1/2 - (^У*3-1/2 =0
т*+1/2 V/ V 3+1/2 т к
где
,.n + v.n xn I xn xn+1 _
j + 1 T Uj __
„.n _ t>j + 1 „n _ j j
2 , = T
Vj+1/2 = 2 , a = a Xi, Xt,j+1/2 = 2 , Xtj = ,
тп = жп = Х/+1 тп = "/'+1/2 + "'_1/2 = „ = Х/+1 Х/_1 (33)
"/'+1/2 = хд,/+1/2 = т , = 2 == 2Т . (33)
Здесь Т = 1/Ж — шаг равномерной сетки <<покрывающей единичный отрезок <5 = [0, 1], жп — узлы неравномерной подвижной сетки на отрезке Й = [0, /], при этом узлы жп являются образами узлов д. Е <^ при некотором гладком невырожденном преобразовании координат
ж = ж(д,г), (34)
ж(0,*) = 0, ж(1,*) = 1, (35)
в каждый момент времени * взаимно-однозначно отображающим < на Й, " = > 0 — якобиан преобразования координат, ж — скорость движения узлов, ^(д, *) = и(ж(д, *),*).
Аналогом формулы (3) является следующая запись схемы (31) в одношаговом варианте:
("^)П+1 - (^)п , (а«)П+1/2 - (а«)П_1/2 Т + Т
й2 \п 1
= 0. (36)
T 2h
a2 Nn / a2
(1 + 0) - (1 + 0)-Vq
J / j+1/2 V J J j-1/2
Исследование устойчивости схемы (36) с помощью метода замораживания коэффициентов [13] приводит к следующим ограничениям:
0n+1/2 > 0, max (Vl + 0 C)n < 1, (37)
j+V j=0,...,N—A /j+1/2" V 7
где Cn+1/2 — число Куранта на подвижной сетке, при этом предполагается выполненным неравенство вида (6):
CW = «(-)n+1/2 < 1 (38)
В случае равномерной неподвижной сетки схема (36) совпадает с (3), а условия (37) переходят при 0 = const в указанные ранее условия (14).
Для того чтобы схема (36) сохраняла монотонность численного решения, достаточно [6] использовать, например, следующий параметр:
0 при |gjn+1/21 < |gn+1/2-S| и j+1/2 ■ j+1/2-s > 0,
j+1/2 H (00(1 - ^))n+1/2 при |gn+1/2| > |gn+1/2-s| и gn+1/2 ' ^+1/2- > 0, (39)
0n,j+1 /2 при <7j+1/2 ' gn+1/2-s < 0,
где
s = sgnan+1/2, °n,j+1 /2 = T^n 1,
C
n
j+1/2
gjn+1/2 = (I a | (1 - C)vq )n+1 /2, j+1/2 - gn+1/2-s
q / I —1—| / 2 ~ .J I -•-/ ~ /1 n
gj+1/2
(40)
Нетрудно проверить, что для выбранного параметра (39) условия (37) выполняются для всех чисел Куранта (38), а на равномерной сетке формулы (39) и (27) совпадают.
2.2. Геометрический закон сохранения
Покажем, что формулы вычисления якобиана и скорости узлов должны быть строго согласованными. Пусть Кп+1/2 — шаг сетки между узлами жп и хп+1, — такой же
шаг на (п+1)-м слое, = ж™+1 — х™+1. Если скорости движения 3-го и (3 + 1)-го узлов
на отрезке времени [¿п, ¿га+1] обозначить через с3- и с3+1 соответственно, то положение узлов на новом временном слое будет определяться формулами
Хп ++ ТС3, Хп+1 Хп+1 ++ ТС3+1,
вследствие чего с необходимостью должно выполняться равенство
3+1/2 = 3+1/2 + Т (сз+1 — С,-) , (41)
которое означает, что ячейка на новом временном слое полностью определяется ее положением на текущем слое и скоростями движения границ. Поэтому равенство (41) называется геометрическим законом сохранения [14, 15].
Учитывая равенство Кп+1/2 = К/п+1/2, легко проверить, что при использовании формул (32), (33) для вычисления якобиана и скоростей узлов геометрический закон (41) выполняется и может быть записан в виде разностного уравнения для якобиана отображения (34):
п
тга+1 _ т
",'+1/2 "3+1/2 Х4,3 + 1 — х4,3 = о (42)
К
Справедливо также тождество
тга+1 тга хп _
"3 "3 Х*,3'+1/2 Х ¿,3—1/2 = о (43)
т К
Выполнение геометрического закона сохранения является необходимым условием согласованности разностных формул, связанных с вычислением характеристик подвижной неравномерной сетки, и означает не только согласованность разностных формул для якобиана и скоростей движения узлов, но и гарантию того, что разностная схема (31), как и дифференциальное уравнение (1), имеет в качестве своего точного решения постоянную функцию, что отвечает общему требованию о сохранении разностной схемой как можно боольшего числа свойств аппроксимируемого дифференциального уравнения. Отметим, что получение согласованных разностных формул в многомерных задачах [11] уже не так тривиально, как в рассмотренном одномерном случае.
2.3. О дивергентных схемах на подвижной сетке
Известно [16], что на разрывных решениях недивергентные схемы могут расходиться, поэтому свойство дивергентности является одним из общих требований, предъявляемых к разностным схемам. Хотя на шаге предиктор схемы (31) аппроксимируется недивергентная форма уравнения переноса (1), записанного в новых переменных
V + а^д = 0, (44)
дивергентность схемы достигается на шаге корректор в результате аппроксимации уравнения переноса в дивергентной форме
(+ (аг)„ = 0. (45)
Покажем простое решение вопроса о построении любых явных двухслойных дивергентных схем на подвижной сетке. Рассмотрим, например, аналог противопоточной схемы (7) на подвижной сетке. Если для построения этой одношаговой схемы использовать аппроксимацию уравнения (44) с односторонними разностями, направление которых определяется в каждом узле сетки знаком величины а, то вряд ли стоит ожидать, что полученная схема будет дивергентной. Поэтому поступим так же, как в случае равномерных сеток, а именно, дивергентные схемы будем получать на основе дивергентной канонической формы (36), выбирая подходящее значение параметра 9. Например, согласно (5) противопоточную схему на подвижной сетке можно получить из (36) при
пп 93+1/2 = 0О,3+1/2.
(46)
Использование формулы (46) приводит к следующей дивергентной форме противопо-
точной схемы:
("г)П+1 — ("*)? + 1
т + К
аг — К |а|гд] — (аг — К |а|гд
2 / 3+1/2 V 2 /3-1/2
0.
(47)
Путем эквивалентных преобразований схема (47) может быть записана в более "привычном" виде (7). Так, преобразование выражения в квадратных скобках приводит к записи схемы через односторонние разности
(Уг)"+1 — ("г)п /а + |а|
+
2
аа
+ I —^ V,
3-1/2
— V
п Х4,3+1/2 Х4,3-1/2
2 "Я I "3
2 ' 3+1/2
К
а геометрический закон сохранения (43) позволяет использовать в производной по времени значение якобиана только с одного слоя по времени
(Уг)п = "п+1гп — тг
П Х*,3 + 1/2 Х*,3-1/2
К
(48)
что приводит к следующей форме записи противопоточной схемы (47):
гП+1 — г"
+
"П+1
а + |а| 2
+
3-1/2
а - | а|
3+1/2
0.
(49)
Противопоточная схема (49) аппроксимирует уравнение (44) с первым порядком, а на равномерной сетке она переходит в схему (7).
Аналогично можно построить дивергентную противопоточную схему второго порядка. Ради простоты изложения будем предполагать, что а > 0, и в соответствии с формулой (9) выберем параметр 9 в случае подвижной сетки следующим образом:
пп
93+1/2
пп
90,3+1/2
1
3-1/Л.
3+1/2 )
Тогда
( а2 \п 1 /а2^ \
(^(1 + 9)уг^ = - (агд)П+1/2 — (аг,3-1/2 + ^
3-1/2
п
п
0
п
п
1
Я
Я
п
и схема (36) становится дивергентной противопоточной схемой второго порядка аппроксимации
(7у)П+1 - 7*)? , (ду)?+1/2 - М?-1/2 к
т + л 2
3+1/2 - 3-1/2
к
3-1/2 - (а^3-3/2 к
т 2к
а2
7^ V 7\
7 ' 3-1/2 / 3-3/2
а2
Недивергентную форму этой схемы можно записать в виде
(51)
73 - (^)П ( аП+1/2 - аП-1/2
т
3
+ ' к ' + 3-1/2 " 3-3/2 -
1 2'
т 2к
7^
п /а2 4 п
3-1/2 V 7 / ¿-3/2_
или, с учетом формулы (48), как
_
т
1
3
1
+ 7+1 1 2(а^)П-1/2 - 2(^)П-3/2-
т 2к
а2
7^ V 7\
7 ^ 3-1/2 V7 /3-3/2
а2
0.
(52)
На равномерной сетке схема (52) переходит в (10).
При использовании схемы предиктор-корректор (31) формула вида (52) применяется для вычисления недостающего граничного значения функции ■ип+1. Например, при а > 0 краевое условие м(0,£) = ^0(£) для уравнения (1) задается в точке х = 0 и "0+ полагается равным (¿п+1). Во внутренних узлах хп+1 = 1,... , N - 1) значения определяются на шаге корректор (31), а в граничном узле х*^1 = I — по формуле (52).
п
п
0
2
0
п
п
3. Особенности построения сеток,
адаптирующихся к разрывным решениям
При расчете разрывных решений на подвижных сетках возникают проблемы, связанные с разрешимостью уравнений для сетки и с адаптацией сетки к разрывному решению. Часто метод построения сетки не позволяет сосредоточить достаточное число узлов сетки в зонах скачкообразного изменения решения. Кроме того, возникающие осцилляции численного решения приводят к возникновению осцилляций значений функции, управляющей плотностью размещения узлов, особенно если в этой функции используются разностные производные высокого порядка. Наличие осцилляций в управляющей функции приводит к немонотонному изменению длин соседних ячеек: происходит чередование длинных и коротких ячеек сетки, в результате чего понижается точность расчета. Некоторые из указанных особенностей конструирования адаптивных сеток для разрывных решений будут продемонстрированы ниже на примере метода эквираспре-деления.
3.1. Метод эквираспределения для построения неравномерных сеток
В методе эквираспределения отображение (34) при £ = 0 определяется как решение уравнения
Их 0)х)д = 0, (53)
удовлетворяющее краевым условиям (35). Здесь £) — заданная управляющая функция. Для определенности предположим, что она задана в виде
£) = 1 + «1|их(х,£)|, (54)
и на п-м временном слое используем следующий сеточный аналог этой функции:
wn+i/2 = 1 +
h j+1/2
j = 0,...,N - 1. (55)
Для вычисления координат x0 узлов неравномерной сетки на начальном слое по времени используется какой-либо итерационный метод (в настоящей работе — метод последовательных приближений) поиска решения нелинейной разностной схемы
1 / __x^ _ x^ \
- (^w°+i/2 j+1h j - W0-1/2 j h j-1j =0, j = 1,...,N - 1, X0 = 0, xN = l, (56)
аппроксимирующей задачу (53), (35). Отметим, что решение разностной задачи (56) удовлетворяет принципу эквираспределения
W+m^V—- = Ch = const, j = 0,..., N - 1, (57)
-j+1/2 -
где
N-1
С = £ И)0+1/2 . (58)
3=о
Предположим, что сетка на п-м слое по времени построена и на ней вычислено решение ■и™. Для нахождения координат х™+1 на новом временном слое используется конечно-разностный аналог уравнения
Их,г)хд )д = вж*, (59)
где в — положительный параметр, отвечающий за уменьшение осцилляций траекторий узлов.
Далее на примере построения начальной сетки будут отмечены трудности, возникающие при использовании метода эквираспределения в случае разрывной начальной функции и0(х).
3.2. Пример с резким изменением длин ячеек в окрестности разрыва
При использовании начальной функции (16) управляющая функция (55) принимает значение 1 во всех ячейках, кроме ячейки , ], в которую попадает точка разрыва ж*, точнее, для которой выполняется неравенство жко < х* < жко+1:
1 при 0 < 3 < к0 — 1,
I п а ,
^•+1/2 ={ 1 + „ ^ „ пРи 3 = ко, (60)
•^ко+1 ко
1 при к0 <3 < N — 1.
Построенная при £ = 0 неравномерная сетка удовлетворяет принципу эквираспреде-ления (57) с постоянной (58)
С = I + «1, (61)
т. е. для шагов сетки имеют место следующие выражения:
кС при 3 = ко,
кз+1/2 = Х3+1 - х3 = ( кС. - а1 при 3 = к°! (62)
Видно, что все ячейки сетки, кроме одной, имеют одинаковую длину кС. = Ах + ка1, несколько большую длины Ах = равномерной сетки с тем же количеством узлов. Ячейка же [хко, хко+1] имеет длину меньшую Ах, а именно, кко+1/2 = Ах - (1 - к)а1. Увеличивая значение а1, длину ячейки [х^0, х^0+1 ] можно сделать сколь угодно малой, но тем не менее получающаяся сетка станет неудовлетворительной в силу того, что соседние с рассматриваемой ячейки будут значительно превышать ее по длине. Таким образом, использование управляющей функции (55) ведет к нарушению условия плавного изменения длин соседних ячеек сетки, следствием чего является потеря точности численного решения при расчете на адаптивной сетке.
Отметим, что из требования х&0 < х^0+1 вытекает ограничение на параметр а1:
I , \
а < N-1. (63)
При нарушении условия (63), т.е. при обращении неравенства в равенство, ячейка [хк0, хк0+1] вырождается в точку.
3.3. Пример с отсутствием решения нелинейной задачи (56)
Покажем, что решение нелинейной задачи (56) с управляющей функцией (60) может не существовать даже при выполнении условия (63). Если решение задачи (56), (60) существует, то координаты узлов неравномерной сетки, согласно полученной формуле (62), вычисляются так:
3(I + «1) при 3 < ко,
хз = ^ N
N (I + «1) - «1 при 3 > ко + 1. В силу этого условие х&0 < х* < х^0+1 можно переписать в следующем виде:
N (I + «1) < х* < ко^+ 1 (I + а1) - аь
Отсюда получаем, что существует целое число ко, 0 < ко < N - 1, удовлетворяющее неравенствам
(х* + а^ x*N
- 1 < ко <
I + а1 I + а1
Пусть I = 1 и х* = 0.5. Тогда число ко удовлетворяет неравенствам
(1 + 2а^ N
- 1 < ко <
2(1 + а) 2(1 + а)
Наконец, в силу условия (63)
1
0 < а <
N - 1 получаем, что
N N
у - 1 < ко < у. (64)
Отсюда видно, что при четных N не существует целого числа к0, удовлетворяющего условию (64), следовательно, и не существует решения задачи (56), (60).
3.4. Осциллирование управляющей функции
Некоторые схемы, сохраняющие монотонность численного решения, могут давать осциллирующие профили разностных производных [9], что может негативно сказаться при использовании метода адаптивных сеток: если управляющая функция, например вида (55), зависит от разностных производных, а они осциллируют, то адаптивная сетка будет иметь чередования длинных и коротких ячеек, что в итоге приведет к потере точности численного решения.
Осциллирование разностных производных численного решения может быть вызвано, в частности, и ростом количества экстремумов при переходе с одного шага по времени на другой. Часто встречается мнение, что ТУВ-схемы не могут увеличивать количество экстремумов. На примере схемы предиктор-корректор (2) с постоянными коэффициентами (12) покажем, что несмотря на то что она сохраняет монотонность численного решения в области монотонности 9 Е [90, 9^] и поэтому является ТУВ-схемой первого порядка аппроксимации, количество экстремумов может увеличиться на следующем шаге, если 9 Е (95, 9^], где
93 = 9о + 2 (9ь - 9о). (65)
Для этого рассмотрим сеточную функцию с одним экстремумом — максимумом в узле с номером j = 0:
п , 0, если j = 0, , ,
и? = < „ „ (66)
j 11, если j = 0 и вычислим значения решения на (n + 1)-м временном слое:
иГ1
0, если j < -2,
b1, если j = -1,
b0, если j = 0,
b-1, если j = 1,
0, если j > 2.
(67)
Легко проверить, что при 9 Е (9s, 9l] выполняется неравенство
b0 < min {b_i, bi} ,
означающее, что функция (67) имеет большее, чем un, количество экстремумов: два локальных максимума и один локальный минимум. Отметим, что несмотря на рост количества экстремумов, полная вариация сеточного решения не может увеличиться, поскольку при 9 Е (9s, 9l] все рассматриваемые схемы обладают TVD-свойством:
оо
TV(ига+1) = £ u+i - ип+1| = 2 (b-1 + bi - bo) = 2 - 4bo < 2 = TV(un).
Отсюда, в частности, следует, что при использовании схемы Лакса (6^ > "ступенька" (24) останется монотонной на (п + 1)-м шаге, а производная
и
х,3+1/2
0,
1/Ь,
если если
3 = -1, 3 = -1,
с одним экстремумом перейдет в функцию
0,
и
п+1 Х,/+1/2
0, (1 0,
если 3 < — 3,
если 3 =- ■2,
если 3 =- 1,
если 3 = 0,
если 3 > 1
(68)
(69)
с тремя экстремумами (двумя локальными максимумами и одним локальным минимумом).
Итак, параметр 68, которому соответствует точка 5 на рис. 1, а, делит схемы с постоянными коэффициентами, сохраняющие монотонность численного решения, на два подмножества: количество строгих экстремумов функции (66) остается прежним при 6 Е [6о, 68) и возрастает при 6 Е (68, 6^]. Значение 6 = является граничным: экстремум остается единственным, но нестрогим, размытым на три узла (и^1 = и
п+1
и
га+Ъ
3.5. Сглаживание управляющей функции
Для устранения указанных выше проблем может быть использовано сглаживание управляющей функции. При сглаживании функции т в область разрыва попадает большее количество узлов и происходит плавное изменение длин соседних ячеек, что позволяет лучше воспроизводить численное решение. Кроме того, расширяется диапазон допустимых значений параметра а1 из формулы (55).
Как показали многочисленные эксперименты, хорошие результаты дает неявная процедура сглаживания [17]
^■+1/2 = ^-+1/2 - ^+1/2 + 2(^-1/2 + тЗ+3/2), 3 = - 2, (70)
где а — произвольный неотрицательный параметр, при 3 = 0 и 3 = N — 1 полагается т = т.
3.6. Результаты расчетов
На рис. 2, а приведены профили 1, 2 численного решения задачи с разрывной начальной функцией (16), полученные по схеме предиктор-корректор с применением описанной процедуры монотонизации. Графики соответствуют моменту времени £ =10 и получены при а =1, длине области I = 30 и числе узлов N = 150. Расчет на адаптивной сетке (рис. 2, б) выполнялся с использованием управляющей функции (55) с коэффициентом а1 = 10, при в = 150, а = 100. Видно, что осцилляции в профилях решения отсутствуют и решение 2 на адаптивной сетке аппроксимирует точное решение задачи
лучше решения 1 на равномерной неподвижной сетке. Кроме того, процедура монотонизации приводит к лучшему описанию разрыва даже на равномерной сетке по сравнению с результатом 3, полученным по самой лучшей из рассмотренных в разделе 1.2 "бездисперсионной" схеме (22).
На рис. 3, а приведены профили численного решения в момент времени £ = 3 начально-краевой задачи для уравнения (1) с точным решением
и(х,£) = е-25(ж-жо-«*02. (71)
а б
Рис. 2. Профили численного решения (а) и траектории узлов адаптивной сетки (б) в задаче с разрывной начальной функцией
аб
Рис. 3. Профили численного решения (а) и траектории узлов адаптивной сетки (б) в задаче с гладкой начальной функцией
В начальный момент времени вершина волны (71) располагалась в точке х0 = 1. Значения других параметров были следующими: а =1, I = 5, N = 150. Адаптивная сетка на рис. 3, б строилась с помощью управляющей функции
<+1/2 = 1 + ао К+1/2| , 3 = 0,..., N - 1, (72)
при этом полагалось а0 = 20, в = 20, а = 10. В расчетах использовалась схема предиктор-корректор с переменным параметром 6. Видно, что решение 2 на адаптивной сетке визуально неотличимо от точного решения и имеет несомненное преимущество по точности перед решением 1, полученным на равномерной неподвижной сетке.
Заключение
В работе рассмотрена процедура монотонизации явных двухслойных схем с помощью специального схемного параметра, выбор которого основывается на исследовании дифференциальных приближений схем. Подробно изучен вопрос о влиянии постоянного, "квазипостоянного" и переменного схемного параметра на монотонность двухслойных разностных схем. Для постоянного схемного параметра построен пример схемы с отсутствием дисперсии в решении второго дифференциального приближения, но не сохраняющей монотонность численного решения. Приведен пример схемного параметра, при котором явная двухслойная схема на подвижной неравномерной сетке является монотонной. Показана связь согласованной аппроксимации якобиана и скоростей движения узлов сетки с геометрическим законом сохранения. Предложен новый подход к построению дивергентных схем на подвижных сетках.
Кроме того, авторами исследованы особенности построения сеток, адаптирующихся к разрывным решениям. Некоторые из этих особенностей продемонстрированы на примере метода эквираспределения. Приведены примеры, связанные с разрешимостью уравнений для сетки и качеством адаптации сетки к решению. На примере схемы предиктор-корректор с постоянными коэффициентами показано, что ТУВ-схемы могут увеличивать количество экстремумов. Для устранения указанных проблем построения адаптивных сеток предложено использовать сглаживание управляющей функции, при котором в область разрыва попадает боольшее количество узлов и происходит плавное изменение длин соседних ячеек, что позволяет лучше воспроизводить разрывное решение задачи. Приведены примеры, показывающие преимущества метода адаптивных сеток.
Подходы к решению проблем построения монотонных разностных схем и адаптивных сеток продемонстрированы на простых примерах для уравнения переноса с постоянным коэффициентом, однако они обобщаются и на другие задачи. Так, предложенный подход монотонизации был применен для построения монотонных разностных схем для нелинейного уравнения переноса и одномерных уравнений мелкой воды [18], в том числе для одномерных задач наката-отката как на равномерных, так и на подвижных неравномерных сетках [19]. Кроме того предложенная технология монотонизации численного решения обобщена для задач о распространении и трансформации поверхностных волн в рамках плановой модели мелкой воды [20].
Список литературы
[1] Годунов С.К. Разностный метод численного расчёта разрывных решений уравнений гидродинамики // Математический сборник. 1959. Т. 47. С. 271-306.
[2] Harten A. High resolution schemes for hyperbolic conservation laws // J. of Comput. Phys. 1983. Vol. 49. P. 357-393.
[3] Sweby A. High resolution schemes for hyperbolic conservation laws // SIAM J. of Numerical Analysis. 1984. Vol. 21, No. 5. P. 995-1011.
[4] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001.
[5] Гольдин В.Я., Кллиткин Н.Н., ШишовА Т.В. Нелинейные разностные схемы для гиперболических уравнений // Журн. вычисл. матем. и математической физики. 1965. Т. 5, № 5. С. 938-944.
[6] Сергеева Ю.В., ХАкимзянов Г.С. Об использовании дифференциального приближения при построении монотонных схем // Вычисл. технологии. 2004. Т. 9. Спец. выпуск: Тр. Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. С. 139-149.
[7] Шокин Ю.И., Сергеева Ю.В., ХАкимзянов Г.С. О монотонизации явной схемы предиктор-корректор // Вестник КазНУ. Математика, механика, информатика. 2005. Т. 2. Спец. выпуск. С. 103-114.
[8] Shokin Yu.I., Sergeeva Yu.V., Khakimzyanov G.S. Construction of monotonic schemes by the differential approximation method // Russian J. of Numerical Analysis and Math. Modelling. 2005. Vol. 20, No. 5. P. 463-481.
[9] ХАкимзянов Г.С., Шокин Ю.И., Барахнин В.Б., ШокинА Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.
[10] Остапенко В.В. О монотонности разностных схем // Сибирский матем. журн. 1998. Т. 39, № 5. С. 1111-1126.
[11] Яушев И.К. О численном расчёте нестационарных течений газа в одномерном приближении в каналах со скачком площади сечения // Изв. СО АН СССР. Техн. науки. 1967. Т. 8, № 2. С. 39-48.
[12] Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985.
[13] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1978.
[14] Trulio J.G., Tigger K.R. Numerical solution of the one-dimensional hydrodynamic equations in an arbitrary time-dependent coordinate system // Tech. Rep. UCLR-6522. Univ. of California, Lawrence Radiation Laboratory. 1961.
[15] Thomas P.D., Lombart C.K. Geometric conservation law and its application to flow computations on moving grid // AIAA J. 1979. Vol. 17, No. 10. P. 1030-1037.
[16] Самарский А.А. Теория разностных схем. М.: Наука, 1977.
[17] Похилко В.И., Тишкин В.Ф. Однородный алгоритм расчёта разрывных решений на адаптивных сетках // Математическое моделирование. 1994. Т. 6, № 11. С. 25-40.
[18] ХАкимзянов Г.С., ШокинА Н.Ю. Численное моделирование поверхностных волн, возникающих при движении подводного оползня по неровному дну // Вычисл. технологии. 2010. Т. 15, № 1. С. 105-119.
[19] Bautin S.P., Deryabin S.L., Sommer A.F. et al. Use of analytic solutions in the statement of difference boundary conditions on a movable shoreline // Russian J. of Numerical Analysis and Math. Modelling. 2011. Vol. 26, No. 4. P. 353-377.
[20] ХАкимзянов Г.С., ШокинА Н.Ю. Определение зоны затопления при накате длинных волн на берег // Труды шестого совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. Алматы: Казак университета, 2009. С. 305-314.
Поступила в 'редакцию 27 февраля 2012 г.
