РОССИЙСКАЯ АКАДЕМИЯ НАУК
Ордена Ленина Институт прикладной математики
им. М. В. Келдыша
Нелинейная монотонизация схемы К. И. Бабенко для численного решения квазилинейного уравнения переноса
Т. А. Александрикова
М. П. Галанин
Москва - 2003
Аннотация
Работа посвящена тестированию предложенного авторами варианта нелинейной монотонизации схемы К.И. Бабенко для численного решения квазилинейного уравнения переноса, а также его сравнению с другими известными конечно-разностными схемами. Для сравнения использовались явная и неявная схемы с левой разностью первого порядка аппроксимации, схема Lax-Wendroff'a второго порядка аппроксимации, а также монотонизованная схема “Кабаре”.
Приведенная информация об ошибках численного решения позволяет сравнить качество представленных схем. Сравнительный анализ ошибок численного решения показал, что предложенная схема дает более высокую точность в широком диапазоне изменения чисел Куранта, а также меньшее “размазывание” разрывов решений.
T.A. Alexandrikova, M.P. Galanin
Nonlinear monotonization of K. I. Babenko scheme for the numerical solution of the quasi - linear advection equation
Abstract
The paper is dedicated to testing the authors offered variant of nonlinear monotonised К.I. Babenko scheme for numerical solution of quasi - linear advection equation, as well as its comparison with other well - known finite - difference schemes. There were used for the comparison the implicit and explicit schemes with the left difference of first - order approximations, the Lax - Wendroff scheme of second order approximations, as well as monotonised "Cabaret" scheme.
The presented information about the errors of numerical solutions allows to compare a quality of the schemes. Benchmark analysis of the errors has shown that the offered scheme gives higher accuracy in the broad range of Courant numbers, as well as smaller "smudging" decision of shocks.
Введение
Уравнение переноса является представителем фундаментальных уравнений математической физики [1], широко используемым для описания движения сплошной среды. Многие задачи газовой динамики, магнитной гидродинамики, астрофизики связаны с решением уравнений такого типа. Учет движения элемента сплошной среды приводит к уравнению переноса квазилинейного типа. Основным следствием нелинейности этих уравнений является “опрокидывание” начального профиля за конечное время и возникновение ударных волн. Изучение закономерностей распространения ударных волн представляет особый интерес в исследовании течений сплошной среды. Наличие разрывов в решении представляет определенные трудности, при этом точность численного решения зависит от качества воспроизведения этих разрывов. На настоящий момент существует большое количество разностных схем для решения задач такого типа. Однако большинство из них дают “расплывчатое” решение, в том числе “размазанный” фронт ударных волн, а схемы повышенного порядка аппроксимации могут приводить к искажениям численного решения, в том числе осцилляциям, т.е. немонотонности решения [2]. Различными способами можно повысить качество численного решения. Один из них связан с использованием неоднородных конечно-разностных схем, в которых реализуется контроль за перемещением разрыва и изменение алгоритма вычислений в его окрестности. Однако основная проблема в использовании таких схем заключается в том, что для нелинейных уравнений гиперболического типа характерно появление разрывных решений даже при гладких начальных условиях, так что положение особых точек заранее неизвестно. К тому же число разрывов может меняться со временем и к каждому из них нужно приспосабливать алгоритм расчета отдельно. Поэтому особый интерес вызывает исследование однородных схем повышенного порядка аппроксимации, а также построение новых квазимонотонных схем на их основе. При этом существуют разнообразные способы монотонизации. В данной работе монотонизация схемы К. И. Бабенко проведена при помощи введения в разностное уравнение искусственной вязкости. Также в данной работе рассмотрена предложенная в [3] схема, монотонизованная с использованием алгоритма, основанного на знании области зависимости точного решения.
В настоящей работе представлен вывод предложенного авторами варианта нелинейной монотонизации схемы К.И. Бабенко для квазилинейного уравнения переноса. Также приведены результаты сравнительного анализа новой и некоторых других известных конечно-разностных схем (явной и неявной схем с левой разностью [2, 4], схемы Lax-Wendroff'a [4] и монотонизованной схемы “Кабаре” [3, 5, 6]). Сравнение численных результатов, полученных по рассмотренным схемам, проведено на системе тестов, аналогичной использованной в работах [4, 5, 7].
Решения получены для финитных начальных условий и различных чисел Куранта на пространственно-временном прямоугольнике достаточно больших размеров. Для определения ошибок численных решений использованы конечномерные анал оги норм пространств C, L1, L2. В работе представлены рисунки с точным и численными решениями для некоторых моментов времени, а также таблицы ошибок для тех же моментов времени и для всего временного промежутка. Такой способ представления результатов позволяет качественно и количественно оценивать характеристики схем.
Авторы выражают благодарность Елениной Т.Г. за внимание, советы и полезные обсуждения.
Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект РФФИ № 03-01-00461).
1. Постановка задачи
Будем рассматривать задачу Коши для квазилинейного уравнения переноса
, , (1.1)
с финитными начальными условиями u(x,0)= u0(x) вида:
1. “треугольник”,
2. “прямоугольник”,
3. “левый треугольник”,
4. “правый треугольник”,
5. “ступенька вниз”,
6. “ступенька вверх”.
Точное решение уравнения (1.1) имеет вид функционального соотношения
u(x,t)= f(x-ut),
где функция f определяется начальным возмущением при :
Поясним решение (1.1), полученное методом характеристик [8]. Заметим, что выражение является полной производной функции u(x,t):
вдоль характеристики G, которая задается уравнением
.
Из выражения (1.1) получаем . Это означает, что функция u(x,t) принимает постоянное значение
u(x,t) = u0()
на кривой G, уравнение которой имеет вид
x = + u0()t.
Приведем точные решения для указанных выше начальных условий:
“треугольник”:
“прямоугольник”:
“левый треугольник”:
“правый треугольник”:
“ступенька вниз”:
“ступенька вверх”:
Сравнение численного решения с точным будем проводить с использованием конечномерных аналогов норм в пространствах C, L1, L2 на =(-,+)[0,T]:
, , ;
и на R=(-,+) при t=ti :
, , .
2. Нелинейная монотонизация схемы К. И. Бабенко
Как известно [5], cхема К.И. Бабенко [9] имеет второй порядок аппроксимации по x и t, но является немонотонной. Cхема является неявной, однако это не мешает вычислять искомое решение на новом временном слое “бегущим счетом”. При числе Куранта, равном единице, для линейного уравнения переноса схема дает точное решение.
Для численного решения введем равномерную (для простоты) пространственно-временную сетку
где h, - шаги разностной сетки по x и t соответственно.
Здесь и далее будем использовать стандартные обозначения для сеточных величин: , где верхний индекс - номер временного слоя, нижний - номер узла по х. Решение на слое n считаем известным.
Перепишем уравнение (1.1) в виде
.
Схему Бабенко (“квадрат” [9]) для рассматриваемого уравнения можно записать следующим образом:
(2.1)
где , или
(2.2)
Пусть , тогда (2.2) перепишем в виде:
Или
(2.3)
где через =a/h обозначено число Куранта.
Заметим, что первые два слагаемых (2.3) дают обычную схему с левой разностью.
Проведем монотонизацию схемы, добавив в (2.3) слагаемые , , где - искусственная диффузия. Поскольку временные производные , на точном решении квазилинейного уравнения переноса аппроксимируют производные , соответственно, то введенные члены дают дополнительные диффузионные слагаемые. Таким образом, модифицированная схема запишется в виде
(2.4)
Заметим, что при ? 0 получим исходную схему (2.1), (2.3), а при ? 1 - схему с левой разностью.
Запишем (2.4) в виде
Где
. (2.5)
Предположим, что y 0 во всей пространственно-временной области, тогда схема (2.5) удовлетворяет принципу максимума [10], если выполнены неравенства
(2.6)
Отметим, что буквальное условие применимости принципа максимума отличается от второго неравенства (2.6): в нем фигурирует . Однако из второго неравенства (2.6) в предположении следует нужное условие.
Введем параметр и положим =(R,),
-1=(R-1,-1). Заметим, что (2.7)
Неравенства (2.6) выполнены, если функция =(R,) имеет вид [3]:
(2.8)
Введенная таким образом функция обеспечивает монотонность рассматриваемой схемы. Отметим, что третьей строке (2.8) соответствует отрезок, на котором схема сохраняет второй порядок аппроксимации для линейных и близким к ним профилей. Четвертая строка обеспечивает переход от нулевого к отрицательным значениям , что соответствует появлению в схеме антидиффузионных слагаемых, приводящих к уменьшению “размазывания” разрывов решения.
Параметр R* в данной работе подбирается экспериментально. Для численных расчетов, как и в [5], использовалось R* =1.2.
Уравнение (2.4) является нелинейным. Из (2.7) следует, что искомую величину можно вычислить по известному R:
(2.9)
Используя полученную зависимость =(R,), найдем соотношения для R. Распишем уравнение (2.4):
(2.10)
которое с учетом (2.7) приобретает вид
(2.11)
Обозначим правую часть (2.11) через b:
. (2.12)
Отметим, что (2.11) представляет собой нелинейное уравнение относительно R с известной правой частью. При этом левая часть (2.11) - монотонная функция R. Учитывая вид функции =(R, ), получим следующие расчетные формулы:
(2.13)
Замечания:
при значение b не определено, однако из (2.7) следует, что в этом случае R = 0.
из (2.13) видно, что R не определено при b = 0, но тогда из (2.12) и (2.10) непосредственно вычисляется искомое значение .
Таким образом, используя формулы (2.9), (2.12), (2.13), можно определить искомую величину . Отметим, что зависит от , поэтому численно величина определяется итерационным способом. В данной работе в качестве первого приближения выбирается равной , для которой определяется соответствующее значение г и по известному b находится значение R. Далее по формуле (2.9) определяется неизвестная величина , тем самым осуществляется переход на следующую итерацию при условии невыполнения критерия окончания итерационного процесса.
3. Другие разностные схемы
Для сравнения предложенной схемы в данной работе были рассмотрены некоторые широко известные конечно-разностные схемы, а именно: явная и неявная схемы с левой разностью, схема Lax-Wendroff'a, а также предложенный в [3] метод прыжкового переноса. В большинстве из них число Куранта выражается через решение не так, как в § 2. Однако вид его достаточно очевиден, так что мы сохраняем обозначение. Ниже представлена информация об основных свойствах этих схем. Более подробно о них можно узнать в указанной литературе.
Схема с левой разностью [2, 4]
Явная схема, имеющая первый порядок аппроксимации по пространству и по времени. Схема устойчива при числе Куранта 1, при этом в случае линейного уравнения переноса для = 1 дает точное решение.
Неявная схема с левой разностью [2, 4]
Абсолютно устойчивая схема, также имеющая первый порядок аппроксимации по пространству и по времени.
Схема Lax-Wendroff'a [4]
Явная схема, аппроксимирующая исходное уравнение со вторым порядком по времени и по пространству. Является устойчивой при 1, при =1 для линейного уравнения переноса дает точное решение.
Метод прыжкового переноса [3, 5, 6, 11]
В работе [3] предложена схема, построенная с помощью монотонизации схемы “Кабаре” [3], и на ее основе разработан алгоритм метода прыжкового переноса.
Трехслойная схема “Кабаре”, аппроксимирующая уравнение конвективного переноса со вторым порядком точности, является бездиссипативной и обладает улучшенными по сравнению с классическими линейными схемами дисперсионными свойствами [6]. Она устойчива при 1, дает точное решение в случае линейного уравнения переноса для = 1 и = 0.5, однако является немонотонной.
Монотонизация проводится при помощи “обрезания” искомого решения из требования выполнения принципа максимума в следующей форме: вычисляемое на новом временном слое значение должно быть ограничено максимальным и минимальным значениями из области зависимости решения. Для контроля возникающего при этом дисбаланса переносимой величины вводится функция “запаса”, что обеспечивает консервативность схемы.
4. Организация вычислительного эксперимента
Аналогично работам [4, 5, 7] численное решение поставленной задачи (1.1) получено для всех указанных форм начальных профилей при следующих значениях параметров: l=520, l1=0, l2=20, T=1000, h=1. Число Куранта () ограничено сверху величиной m, которая принимала значения: 0.1, 0.25, 0.5, 0.9. Для неявной схемы с левой разностью также были проведены расчеты при m = 3. Таким образом, временной шаг принимал те же значения, что и число Куранта. Для определения точности полученного при этом решения вычислялись относительные ошибки с использованием норм
, ,
интегрально по временному промежутку 0 t T, а также
, ,
локально для момента времени tj.
В рамках данной работы для проведения численного эксперимента была создана программа для персонального компьютера, реализующая алгоритмы поиска решения по указанным разностным схемам. На экран монитора выводится графическое изображение точного и численного решений для каждого временного слоя, а также относительные ошибки вычислений для нескольких промежуточных и конечного моментов времени. Такое представление результатов позволяет визуально и численно оценить качество полученного решения и тем самым сравнить используемые схемы.
5. Результаты численного исследования
Полученные в численном эксперименте результаты представлены следующим образом. Для каждой схемы, используемой для получения численного решения, приведены таблицы интегральных и локальных по времени ошибок для различных чисел Куранта. В каждой таблице указаны относительные ошибки вычислений для всех начальных условий и используемых норм, соответствующее число Куранта приведено в левом верхнем углу таблицы.