CC BY
79
Вычислительные технологии
Том 18, № 3, 2013
Метод адаптивных сеток для одномерных уравнении мелкой воды*
Г. С. Хлкимзянов, Н.Ю. ШокинА Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]
На примерах для нелинейного скалярного уравнения обсуждаются вопросы построения разностных схем, сохраняющих монотонность численного решения, на равномерных и адаптивных сетках. Показаны важные свойства сохранения постоянного решения, стационарного и движущегося скачков. С помощью метода дифференциального приближения дано новое объяснение механизма возникновения нефизичных численных решений и проведена энтропийная коррекция. Представлен пример TVD-схемы, увеличивающей количество экстремумов. Предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках. Схема предиктор-корректор, построенная для нелинейного скалярного уравнения, применена для решения одномерных нестационарных уравнений мелкой воды. Исследованы свойства схемы и проведено её численное тестирование.
Ключевые слова: численное моделирование, конечно-разностная схема, нелинейное скалярное уравнение, предиктор-корректор, монотонная схема, дивергентная схема, энтропийная коррекция, адаптивная сетка, метод эквираспределения, нелинейные уравнения мелкой воды.
Введение
Достоверная оценка зон затопления при накате волн цунами на берег является актуальной проблемой математического моделирования. Для получения эффективного прогноза необходимы, кроме качественных, исследованных алгоритмов, большие массивы подробной информации о рельефе дна и суши в прибрежной зоне, о параметрах источника волны, данные наблюдений для верификации и т. д. Как правило, вся эта входная информация весьма неточна. Кроме того, модель мелкой воды, используемая для моделирования процесса наката волн цунами на берег, выводится при достаточно жёстких ограничениях. Поэтому для получения удовлетворительных результатов в указанных задачах использование схем повышенного порядка точности не имеет смысла и достаточны схемы первого или второго порядка. Это определяет необходимость построения качественных дивергентных схем первого или второго порядка аппроксимации, точных на возможно широком классе решений исходной задачи, в частности, сохраняющих состояние покоя, монотонность либо решения, либо инвариантов Рима-на и дающих численное решение, имеющее наиболее широкий набор свойств решения исходной задачи. Нелинейное скалярное уравнение является модельным для системы уравнений мелкой воды, поэтому оно использовано в первой части работы для исследования свойств конструируемых разностных схем, применяемых далее для решения задач о течениях в рамках модели мелкой воды.
* Работа выполнена при финансовой поддержке РФФИ (грант № 12-01-00721а), а также в рамках Программы интеграционных исследований СО РАН (проект № 42).
При решении задач с разрывными решениями нарушение условия неубывания энтропии [1 -3] приводит к нефизичному численному решению в виде ударной волны разрежения в задачах о трансзвуковых течениях идеального газа, а также в окрестности критической точки, расположенной в транскритической волне разрежения в задачах о течениях идеальной жидкости в рамках модели мелкой воды [4, 5]. При численном решении этих задач необходимо, чтобы выбранный метод сходился к физически верному (энтропийному) решению, что гарантируется выполнением дискретного аналога условия неубывания энтропии. Для выполнения данного условия можно использовать известный приём энтропийной коррекции, который был предложен в работах [6, 7] и применялся в [4, 8-10] для метода Роу ив [11-14] для других численных методов. В настоящей работе показаны важные свойства сохранения схемой предиктор-корректор постоянного решения и стационарного или, в случае подвижной сетки, движущегося скачка. С помощью метода дифференциального приближения дано новое объяснение механизма возникновения нефизичных численных решений и проведена энтропийная коррекция схемы.
Использование адаптивных подвижных сеток для схем первого или второго порядка аппроксимации повышает точность численного решения, но главная причина их применения в задачах моделирования процесса наката волн цунами на берег состоит в том, что использование разностных схем в методах с выделением подвижной линии уреза в принципе невозможно без подвижных сеток. Кроме того, в двумерном случае необходимо учитывать криволинейность границы расчётной области. Поэтому в настоящей работе большое внимание уделяется описанию особенностей моделирования на подвижных сетках. В частности, предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках (аналогичный подход для линейного уравнения переноса рассмотрен в [15]).
Если при использовании адаптивных сеток управляющая функция, регулирующая расстановку узлов сетки, зависит от разностных производных, которые осциллируют, то сетка будет иметь чередование длинных и коротких ячеек, что приведёт к потере точности численного решения. Одна из причин появления осцилляций разностных производных — увеличение схемой количества экстремумов решения при переходе с одного шага по времени на другой. Известно [7], что если начальная функция имеет ограниченную вариацию, то количество локальных экстремумов слабого решения не возрастает. Вместе с тем некоторые ТУВ-схемы, сохраняющие монотонность численного решения, могут увеличивать количество экстремумов решения [15-19]. В настоящей работе этому вопросу уделено внимание на примере схем Лакса, противопоточной и предиктор-корректор для невязкого уравнения Бюргерса.
Во второй части работы схема предиктор-корректор, построенная для нелинейного скалярного уравнения, применена для решения одномерных уравнений мелкой воды. Исследованы свойства схемы и проведено её численное тестирование.
1. Разностные схемы для нелинейного уравнения, сохраняющие монотонность численного решения
Данный раздел посвящён разностным схемам для однородного
и + [/(«)]* = 0
(1)
и неоднородного
М + [/(м)]х = (2)
нелинейных скалярных уравнений, которые являются модельными для системы уравнений мелкой воды в случае ровного горизонтального и неровного дна соответственно.
1.1. Схема предиктор-корректор на равномерной сетке
Сначала рассмотрим случай равномерной неподвижной сетки с узлами ж. = и ша-
лп
гом Л > 0. В явной схеме предиктор-корректор [20, 21] искомая сеточная функция и™ определённая в момент времени ¿п в узлах ж., находится из разностного уравнения
- + /^+1/2 /^-1/2 = 0 (3) т Л '
аппроксимирующего закон сохранения (1). Необходимые для выполнения этого шага потоки /*, отнесённые к полуцелым узлам ж^+1/2 = ж. + Л/2, определяются на шаге предиктор
/+1/2 - 2 (/3++1 + 7) /™+1 - /™
7+1/2 2 7+ 7 + -0+1 -О = 0 (4)
т * + а^-+1/2 Л =0, (4)
где /. = / (ип), т*+1/2 = 0.5т(1 + $п+1/2), т — шаг по времени, ^^+1/2 — схемный параметр. Разностное уравнение (4) есть результат аппроксимации дифференциального уравнения для потоков
/ + а(")/ = 0, (5)
которое получается после умножения (1) на функцию а (и) = /и(и).
Наличие некоторых важных свойств схемы (3), (4) напрямую связано с видом сеточной функции ап+1/2, аппроксимирующей функцию а(и). В настоящей работе способ вычисления функции ап+1/2 определяется требованием сохранения для сеточных функций правила дифференцирования /х = а(м)мх функций непрерывного аргумента, т. е. выполнением равенства
/П^+1/2 = аП+1/2Иа;,^+1/2 (6)
для разностных производных
/п _ /п ип _ ип
<-п = ^ мп = ".7+1
/ж,^+1/2 = л , "х,.+1/2 = Л .
Условию (6) удовлетворяет, например, сеточная функция [7]
п
/п _ / п
ап I 71+1 Мп пРи "п+1 = " а(ип) при мп+1 = мп,
применение которой позволяет переписать уравнение (4) с использованием разностных производных от искомого решения
/ п I
/7+1/2 = ^ _ 10 +
и получить следующий одношаговый вариант схемы (3), (4):
ип+1 ип /п /п
т 2Н 2Н
((1 + 0)аЧ)П 2 - ((1 + 0)а2ил)
¿-1/2
0. (9)
Схему (9) можно назвать канонической формой явных дивергентных двухслойных схем для уравнения (1), поскольку любая такая схема может быть записана в виде (9) при соответствующем выборе параметра 0. Например, при ^*П+1/2 — 0 схема (9) совпадает со схемой Лакса — Вендроффа
ип+1 /п /п
+
т
2Н
2Н
2
а ^ ¿+1 /2 1а ^¿-1/2
2
0,
при
¿+1/2 - ^0,^+1/2 = СП 1 > 0,
п
7+1/2
где Сп+1/2 — число Куранта,
Сп+1/2 = К+1/2|ж < 1, Ж =
10)
получается противопоточная схема первого порядка аппроксимации, которая может быть записана как в дивергентной, так и в недивергентной форме [22]
м™+1 - ил + (а + |а|
т
а - | а| п
о иж I + —о— иж = 0,
2 /¿-1/2 V 2 /¿+1/2
:12)
при ¿+1/2 - ^,¿+1/2
(Сп )2 (С7+1/2)
— 1 > 0 — схема Лакса
ип+1 - 2(ип+1 + ип-1) + /¿+1 - /¿-1
2Н
13)
при 0°+1/2 — -1 — абсолютно неустойчивая схема второго порядка аппроксимации по Н
ип+1 _ ип /п _ /п и_и! + /^'+1 /^-1 = 0
т
2Н
Явные схемы с несимметричными шаблонами также получаются из канонической формы (9). Например, при а > 0 и
пп ^¿+1/2
пп
^,¿+1/2 I
л - ¿1И\
где
¿+1/2 = (Iа|(1 - С)иж)"
:14)
¿+1/2
схема (9) переходит в противопоточную схему второго порядка аппроксимации на несимметричном шаблоне
и™+1 - ип 3 (аих)°—1/2 - (аиж)°-3/2 _ Т (аиж)°-1/2 - (аих)п-3/2
+
Т
Н
п
1
2
2
Если 9 = 0(h), то схема (9) аппроксимирует уравнение (1) со вторым порядком относительно т и h. При 9 = const = 0 порядок аппроксимации на гладких решениях — первый. Необходимое условие устойчивости для случая a(u) = const [15] требует, в частности, чтобы параметр 9 принимал неотрицательные значения. Поэтому далее предполагается выполненным ограничение
9+1/2 >
В общем случае (для произвольного параметра 9 = 0(h)) схема предиктор-корректор (3), (4) не сохраняет монотонность численного решения. Однако этот параметр можно подобрать так, что схема будет обладать указанным свойством. В работе [23] для выбора параметра 9 использовалось первое дифференциальное приближение (п.д.п.) схемы предиктор-корректор
т h2
ut + fx = + y [a(a2œ2 - 1)«x]xx (15)
и параметр 9 подбирался так, чтобы диссипативный член п.д.п. приводил к изменению дисперсионного члена в тех подобластях, где возможно появление осцилляций в численном решении. Например, если положить
(а(1 — lai œ)ux)
9 = h^-121 ; , (16)
œa Ux
то из (15) получается уравнение с дисперсионным членом
ah2
— (1 - |a|œ) (2 - |a|œ) Uxxx, 6
который в силу условия (11) отличается по знаку, например, от дисперсионного члена
ah2 2 2
6
(a œ — 1)ux
в п.д.п. схемы Лакса — Вендроффа, не сохраняющей монотонность численного решения. В [23] при некоторых дополнительных ограничениях на функцию (16) получена формула для сеточной функции 9n+1/2
0 пРи |^П+1/2 |<|gn+1/2-s| и g+1/2 ■ gn+1/2-s > °>
9П+1/2 Н (9°(Х - ^))n+1 /2 при |g+1/21 > |g+1/2-s I и gn+1/2 ' g+1/2-s > 0 (17) 90j+1/2 пРи g+1/2 ■ gj+1/2-s < °
4 = Sgn П™ A™ = gj+1/2-S
5 = Sgn aj+1/2' 4j+1/2 = g ,
и доказано, что при условии (11) схема предиктор-корректор (3), (4), (17) сохраняет монотонность численного решения.
Отметим два свойства схемы (3), (4), связанные с её точными решениями.
Лемма 1. Схема предиктор-корректор (3), (4) сохраняет постоянное решение: если u™ = U0 = const, то и u°+1 = U0.
где
Второе свойство связано с сохранением разностной схемой разрывных решений задачи Римана для уравнения (1) с начальной функцией
и°(х) = I Ц, Х < х°: Ц = Ц. (18)
Пусть функция потока является выпуклой, т.е. /"(и) > 0 (для вогнутой функции потока все дальнейшие рассуждения справедливы при противоположных соотношениях между Ц и Ц). Тогда для Ц > Ц. задача (1), (18) при £ > 0 имеет следующее решение в виде ударной волны:
«(х.«>={иг при х < хо+(19)
где Ш — постоянная скорость движения скачка, определяемая из условия Ренкина — Гюгонио:
ш = лад/). (20)
Цг - Ц
Если Ш = 0, то скачок является стационарным, и в этом случае
/ (Ц) = / (Ц). (21)
Решение (19) удовлетворяет условию неубывания энтропии [4, 5] и является физически корректным.
Известно [4, 5], что точное решение задачи Римана в виде стационарного скачка является и решением некоторых разностных схем. Покажем, что этим свойством обладают все схемы вида (9), причём не только второго порядка аппроксимации, но и первого, например, противопоточная (12).
Лемма 2. Схема предиктор-корректор (3), (4) сохраняет стационарный скачок (18), (21).
Доказательство. Пусть на п-м слое по времени (п > 0) скачок располагается между узлами х^0 и х^+1^, т.е.
ип = { Ц при 3 - 3°' (22)
'3 \ Цг при 3 > 3° + 1.
В силу определения (22) и равенства (21) получаем, что для всех 3 выполняется равенство (аих)™+1/2 = 0. Тогда из формулы (9) с учётом равенства (21) следует, что ип+1 — ип, т. е. решение не меняется при переходе с одного временного слоя на другой.
При Ц < Цг разрывное решение (19) в виде ударной волны разрежения является математически возможным для задачи (1), (18), но физически некоректным (неустойчивым) и нарушающим условие неубывания энтропии. Устойчивым же является непрерывное решение в виде центрированной волны разрежения
и(х,£) = Ц при х — х° + а(Ц)£, х — х°
а(и) = —-— при х° + а(Ц)£ — х — х° + а(Цг)£, £ > 0,
и(х,£) = Ц при х > х° + а(Цг)£,
которое не нарушает энтропийное условие [4, 5]. Например, для невязкого уравнения Бюргерса (/(и) = и2/2) это решение имеет вид
Ц при х < х0 + Ц¿, и(х, г) = ^ Х ^Х0 при хо + < х < хо + Цг, г > 0. (23)
при х > х0 + иг.
Поскольку схема предиктор-корректор при Ц < Ц дает нефизичное решение невязкого уравнения Бюргерса, то, как указано выше, необходима её корректировка для возможности расчётов решений со сменой знака величины а. Энтропийная коррекция сводится в итоге к принудительному вводу в схему дополнительной вязкости в тех подобластях, где значение |а| становится меньше некоторого порогового и исчезает численная диссипация [24]. Покажем, что такой же приём может быть получен на основе метода дифференциального приближения, причём с явной формулой для порогового значения.
Раскрывая скобки в правой части п.д.п. (15), получим, что схемная вязкость ц (коэффициент при второй производной ихх) выражается формулой
ц = 29а2 + у (3а2ж2 - 1)ил. (24)
Следовательно, при малых значениях а и при > 0 схемная вязкость становится отрицательной, что и является причиной появления нефизичного численного решения при расчёте волн разрежения (23) с критической точкой а = 0. Отметим, что для невязкого уравнения Бюргерса в волнах сжатия (Ц > Ц) выполняется условие < 0, поэтому отрицательная вязкость в окрестности точки с нулевой критической скоростью а = 0 не возникает, а значит, не возникает нефизичных численных решений.
Из формулы (24) следует, что для обеспечения неотрицательности вязкости ц можно заменить функцию ($а2)™+1/2 на следующую:
¿+1/2 пРи
(9а2)П+1/2 < ¿+1/2,
¿1/2 = { ¿+1/2 """ I СП+1/2 < 1/л/3, иП,^.+1/2 > 0, (25)
(9а2)П+1/2 в других случаях, где функция 9 вычисляется по формуле (17),
¿+1/2 = ^ ((1 - 3С2К)П+1/2 . (26)
С учётом коррекции схемной вязкости (25) шаг предиктор (8) запишется в виде
/п + /п т , л п
¿^=- 2 О2+(27)
1.2. Схема предиктор-корректор на подвижной сетке
Рассмотрим некоторую начально-краевую задачу для уравнения (1) на отрезке [0, /]. Предполагается, что узлы хп подвижной неравномерной сетки, покрывающей в каждый момент времени этот отрезок, связаны взаимно-однозначно соответствием х = х(^,£)
с узлами qj равномерной неподвижной сетки с шагом Н = 1/Ж, покрывающей единичный отрезок [0, 1]. На подвижной сетке схема предиктор-корректор для скалярного уравнения (1) имеет следующий вид:
'+1/2 2 7 +7 + (Л - =0, (28)
Т.+1/2 V 7 / .+1/2
/7+1/2 2 /+1 + . + (а/ - Г = 0 (29)
где
си)
2 ^ .7+1 1 Ь
Т.+1/2 V 7 / .+1/2
(/Ц)и+1 - (7ц)и + (/* - хии-%+1/2 - (/* - хии*).-1/2 = 0, (30)
т Н '
и _ и и / и ___
и _ -.'+1 и (-и _ /.+1 / _
„ и _ Ии /и _ /и хи+1 _ хи
Ии = и+1 и Ъи = Ъ .+1 хи х
ид.+1/2 = 1 , /д.+1/2 = , ,
Н ' +1/2 Н ' т
хи I хи хи хи 7и | 7и
Х*,.7 + ти _ Х.+1 _ и ти _ л+ 1/2 + 1/2
и _ 1 + 1 ти _ .7 + 1 .7 _ и т; _
х4,.+1/2 = 2 , ".7+1/2 = н = хд,.+1/2, 7 = 2 '
при этом сеточная функция аи+1/2 вычисляется по формуле (7) и выполняется равенство вида (6)
/п.:/+1/2 = (аид)и+1/2 • (31)
Дивергентное разностное уравнение (30) шага корректор аппроксимирует уравнение (1), записанное в новых координатах в дивергентной форме
(Ти)* + (/ - = 0.
Недивергентная форма последнего уравнения
1
7
И + Т (/д - ) = 0
используется при получении шага предиктор (28). Еще одно уравнение шага предиктор (29) получается в результате аппроксимации уравнения для потоков (5) в новых координатах ¿):
а
Л + 7 (/д - ) =
Выписанная схема (28)-(30) имеет при 9 = О(Н) второй порядок аппроксимациии на подвижной сетке. При этом на равномерной сетке (х* = 0, 7 = /) разностные уравнения (29), (30) полностью совпадают с приведёнными выше уравнениями (4), (3). Отличие от случая равномерной сетки состоит, в частности, в том, что на подвижной сетке необходимо вычислять не одну, а две предикторные величины и* и /*. Однако в уравнении шага корректор эти предикторные величины используются в комбинации /* - ^и*, для которой, согласно уравнениям (28), (29), справедливо равенство
(/* - хии*).+1/2 = (/ - М") где ( и
Р _ т * I (а - \
/7+1/2 = /7+1/2 - Т.+1/^1 7 )
Ид /
7+1/2
/п I /п /п = + ^
¿¿+1/2 = ■
и
2 ' ^'+1/2 2 "
В результате разностная схема предиктор-корректор (28)-(30) на подвижной сетке состоит из шага предиктор
и п+1 I и п
/ - /п
т"
а2
+ "ги ¿+1/2 /¿+1/2
(32)
где а = а - х*, и шага корректор
(7и)п+1 - (7и)п
+
/ - (х*и)п) - (/ - (х4и)п) ¿+1/2 V /
¿-1/2
т
н
0,
(33)
т. е. на шаге предиктор достаточно вычислять только одну вспомогательную величину ./¿+1/2.
Рассмотренный в предыдущем разделе способ выбора параметра 9 можно использовать и для монотонизации схемы предиктор-корректор (32), (33) на подвижной сетке. Для этого в п.д.п. данной схемы
(7и)* + (/ - х*и)
9 2
а2 97и9
н2 + (Г
/аж\ 2 а | - 1 ) щ
+
99
т2 +6
аа 2
7
7
н2 + 12
(34)
2 7 + ) и9 <Л:и99
следует взять параметр 9 таким, чтобы изменился знак коэффициента при третьей производной и999, например, в виде (16):
9=н
|а|ж
а ( 1--— I щ
жа2
7
и9
Тогда для вычисления сеточной функции 9¿+1/2 можно использовать ту же формулу (17), что и в случае равномерной сетки, с тем лишь отличием, что на подвижной сетке в = sgnап+1/2, а формулы (14), (11) модифицируются в следующие:
¿+1/2 = (|а|(1 - СК)п+1/2, Сп+1/2 = * (М)
ж I I < 1, ¿+1/2
(35)
при этом схема предиктор-корректор (17), (32), (33) будет сохранять монотонность численного решения и совпадать с ТУВ-схемой [25], являющейся обобщением схемы Хар-тена [7] на случай подвижных сеток.
Для схемы на подвижной сетке вязкость ц выражается более сложной, чем (24), формулой. Так, для невязкого уравнения Бюргерса из (34) получаем
т „ а2 Н2 / (аж \2 \ т2 (а \2/т а т . .„ т Ц = т9 — + — (3( —) - 1) и - — (-) Ц + -Л) + —^
2 V/
7
Н!
4
п
0
9
9
9
Численные эксперименты показали, что для исключения возможности появления нефи-зичных численных решений достаточно потребовать, чтобы в окрестности критической точки (а = 0) волны разрежения выполнялось неравенство
т^а2 Н2 / (азд\2 \
2 77 + ¥ 3(т) - 0 * 0
Отсюда получаем следующий аналог формулы (25):
ф
£+1/2
и
°7+1/2
2и
,а2 \
при
2и
9")
< Ли
7+1/2
Си+1/2 < 1/л/3, ии
д,7+1/2
> 0,
9— | в других случаях,
7)
7+1/2
(36)
где функция ^и+1/2 вычисляется по формуле (26) (с заменой на и использованием числа Куранта (35) на подвижной сетке).
Вследствие коррекции схемной вязкости шаг предиктор (32) претерпит изменения и запишется в виде формулы, аналогичной (27):
т / а2 \и
Л-+1/2 = /и+1/2 - 2((7 + ^Ч^
Одношаговый вариант схемы (32), (33) имеет дивергентный вид
(7и)и+1 - (7и)и (/ - х*и)и+1/2 - (/ - х*и)и_ 1/2
+
Н
2Н
(1 + 9)
7+1/2
- ( (1 + 9) 7Ид
7-1/2.
0,
(37)
и при подходящем выборе параметра 9 из схемы (37) могут быть получены любые явные двухслойные дивергентные схемы на подвижной сетке. Например, при вычислении параметра 9 по формуле (10), в которой используется число Куранта для подвижной сетки, получается противопоточная схема в дивергентной форме на подвижной сетке:
(7и)и+1 - (7и)и 1
т + Н
/ - - Н|а|иЛ - (/ - - Н|а|ид
2 /7+1/2 V 2 /7-1/2,
0. (38)
Путём эквивалентных преобразований схема (38) может быть представлена в виде (12) с односторонними разностями. В самом деле, используя равенство (31), запишем уравнение (38) в недивергентной форме
(7и)и+1 - (7и)и (а + |а|
+
т
2 И
аа + ( —^ И
7—1/2 V 2 V 7+1/2 7
ии и Х*,7+1/2 Х*,7—1/2
Н
0.
Геометрический закон сохранения [15]
Т'г+1 = 7/7 + т )
и д )7
(39)
и
и
2
2
и
и
позволяет использовать в производной по времени значение якобиана только с одного слоя по времени
хп хп
(7и)п = 7п+1ип - тип '¿+1/2 - ''¿-1/2 , что приводит к эквивалентной форме записи противопоточной схемы (38)
ип+1 - ип +
т ¿+1
а + |а| \п (а - |а|
-о-и 9 + -о-и9
2 7,-1/2 V 2 7,-+1/2
0.
Отметим, что утверждение леммы 1 остается справедливым и для подвижных сеток. Что касается леммы 2, то для случая подвижных сеток она может быть переформулирована в виде утверждения о том, что точное решение (19) задачи Коши (1), (18) является и точным решением схемы (32), (33), если левый и правый соседние со скачком узлы с номерами и + 1 движутся так, что выполняется равенство
п
xí'¿0+1/2
Ж (40)
Лемма 3. При условии (40) схема предиктор-корректор (32), (33) сохраняет движущийся скачок (18), (19).
Доказательство. Пусть на п-м слое по времени (п ^ 0) скачок располагается между узлами х™ и хп0+1, т.е. выполнено определение (22). В силу (20) и условия (40) получаем, что для всех ] выполняется равенство
(аи9 ¿+1/2 = (41)
Поэтому из (32) следует, что /¿+1/2 = /¿+1/2, а из (33) получаем
п+1 -( Т„лп т ( гп , гп ) , ¿х*и)п+1/2 (х*и)п_1/2 ¿ = ^^ - 2 ^ 9,^+1/2 + ^'¿^^ + т н '
Используя формулу дифференцирования произведения двух функций
(х*и)п+1/2 - (х4и)п-1/^ 1 1 -Н- = (х*9и)п + 2 Кх*и9)п + 1/2 + (х*и9)п'-1/^ ,
а также равенства (31) и (41), получаем формулу (7и)п+1 = + т(х*9)п) ип, которая в силу тождества (39) означает, что ип+1 = ип. □
Приведём некоторые результаты экспериментального исследования свойств схемы предиктор-корректор. На рис. 1,1, а показаны профили 1, 2 численного решения на отрезке [0, /] начально-краевой задачи для невязкого уравнения Бюргерса (/ = и2/2) с непрерывной начальной функцией
Ц при х < хг
и0(х) = ^ -- Ц I--Ц при хг < х < хг
^^г ,х г ,х г хг
Ц при х > хг.
1
Для лучшей детализации изображены только фрагменты профилей, расположенные вблизи точки разрыва решения. Графики соответствуют моменту времени £ =10 и получены при Ц = 1, Ц = -1, хг = 10, хг = 20, длине области I = 30, числе узлов N = 60 и постоянном шаге по времени т = к3ап^, где кзап — коэффициент запаса, значение которого во всех расчётах полагалось равным 0.2.
Точным (слабым) решением задачи является центрированная волна сжатия
{Ц при х < хг +
'X 'X гЦ
- при хг + Ц £ < х < жг + Цг £<£ *,
£ £ *
Ц при ж > хг + Ц
которая в момент градиентной катастрофы
'X г _ _ _ _
£ * г, х * хг + Цг£ * Хг + £ * Ц — иг
генерирует устойчивый скачок
«(х,£)Н Цг при х<х* + Ш(£ - £*), +
Ц при х > х* + Ж(£ - £ *), 2
При указанных выше значениях входных данных градиентная катастрофа происходит при £ * = 5 в точке х * = 15 и возникающий скачок является стационарным (Ш = 0).
Адаптивная сетка (см. рис. 1,1, б) строилась методом эквираспределения [15] с использованием управляющей функции
i _ i
<+1/2 = 1 + а: 'Г _ '' (42)
с коэффициентом = 15, параметром неявного сглаживания а = 60 и параметром в = 80, обеспечивающим гладкость траекторий узлов [15].
Из рис. 1,1, а видно, что схема предиктор-корректор (17), (32), (33) сохраняет монотонность численного решения как на равномерной, так и на подвижной сетке, при этом имеет место высокая разрешимость скачка, но решение (2) на адаптивной сетке аппроксимирует точное решение задачи (тонкая сплошная линия) лучше, чем решение (1) на равномерной сетке. Так, на равномерной сетке скачок размывается на две ячейки и имеет ширину d =1. На подвижной же сетке, адаптирующейся вначале к волне сжатия, а затем к стационарному скачку, последний не размывается вовсе и располагается всегда между одними и теми же двумя узлами подвижной сетки (положение стационарного скачка в точном решении показано на рис. 1,1, б штриховой линией). Отметим, что в расчётах по схеме предиктор-корректор без процедуры монотонизации, а также по схеме Лакса — Вендроффа осцилляции в численном решении после наступления градиентной катастрофы возникали всегда.
На рис. 1,11, а приведены профили 1, 2 численного решения рассмотренной выше задачи при использовании всех прежних значений параметров за исключением значения Ц = 0. В этом случае волна сжатия генерирует при £* = 10, х * = 20 скачок, движущийся с постоянной скоростью Ш = 0.5. Видно, что на равномерной сетке скачок размывается на четыре ячейки ^ = 2). Из профиля на адаптивной сетке следует, что размывание здесь также имеет место, но с шириной зоны ^ = 0.08) существенно
I
а б
II
а б
Рис. 1. Профили численного решения задачи о генерации стационарного (I) и движущегося (II) скачка, полученные на равномерной (1) и адаптивной (2) сетках в моменты времени £ = 10 (I) и £ = 20 (II) (а); траектории узлов адаптивной сетки (б)
меньшей, чем на равномерной сетке. Отметим, что для рассматриваемого численного эксперимента условие (40) не выполнялось, поскольку узлы сетки при своём движении медленно пересекали фронт скачка. Возможно, что при тщательном подборе параметров а1, в, 0, влияющих на степень адаптации сетки к решению, удалось бы построить сетку с выполненным на скачке условием (40) и, как следствие, с сохранением на адаптивной сетке скачка без размывания (как в лемме 3), но такая цель в данной работе не ставилась.
На рис. 2, а представлены профили численного решения задачи Римана о формировании волны разрежения из начального скачка (18) при хо = 15, Ц = -1, Ц = 1 и тех же значениях других параметров, как в предыдущих двух тестовых задачах. Гра-
а б
Рис. 2. Профили численного решения задачи о генерации центрированной волны разрежения, полученные в момент времени Ь = 10 на равномерной (1) и адаптивной (2) сетках с коррекцией вязкости, на равномерной сетке без коррекции (3) и с неполной коррекцией вязкости (4) (а); траектории узлов адаптивной сетки (б)
фик точного решения (23) изображён тонкой сплошной линией, траектории движения границ волны разрежения (23) показаны на рис. 2, б штриховыми линиями. Из рис. 2, а следует, что при использовании схемы предиктор-корректор с коррекцией вязкости решение 2 на адаптивной сетке аппроксимирует точное решение задачи лучше, чем решение 1 на равномерной сетке. Однако в последнем примере преимущество использования адаптивной сетки не столь безусловное, как в предыдущих тестах. Это связано с тем, что к моменту времени £ = 10 волна разрежения стала протяжённой и очень пологой, поэтому согласно выбранной управляющей функции (42) сетка к данному моменту времени стала почти равномерной (см. рис. 2, б). На рис. 2, а приведены также графики, полученные без использования коррекции вязкости 3 и с неполной коррекцией 4, когда значение 8 в формуле (25) заменялось на уменьшенное в десять раз. Видно, что в первом случае получается нефизичное решение в виде стационарного скачка (19), а во втором — нефизичное решение в виде двух волн разрежения, сопряжённых через стационарный скачок. Этот эксперимент показывает, что по крайней мере порядок величины (26), полученной на основе анализа п.д.п. рассматриваемой схемы, определён верно.
1.3. Схема для уравнения с источниковым членом
Для неоднородного уравнения (2) можно написать аналог разностной схемы (32), (33), основываясь на аппроксимации уравнений
1 а
иг + - (/ - ХгЩ) = д, ¡г + - (/ - Хиа) = ад (43)
для шага предиктор и уравнения
(-и)г + (/ - Хги)а = -д (44)
для шага корректор. Вначале выписываются неоднородные аналоги разностных уравнений (28)-(30)
¿1/2 - ип+1/2 ( /9 - х4и9
+ ^^ = ¿+1/2, (45)
¿+1/2 V 7 / ¿+1/2
/п+1/2 - ¿+1/2 . а(/9 - х*и9)
+ "9 т ' И = (ау)п+1/2, (46
где
т * 'V 7 I ^^)¿+1/2,
' ¿+1/2 V 7 /¿+1/2
(/и)п+1 - (7и)п + (/* - хпи^+1/2 - (/* - ХnU*)п-l/2 = (7^)*,
т Н ¿'
пп
х + х +1
хп
^п+1/2 = #(хп+1/2 ,ип+1/2), хп+1 / 2 = 2 , = ^(х*,^*,и*),
т т^1 - 7
¿* = ¿п + тт* = -(1 + г) х* = хп + ¿х™, = 7п + т*
2 ' ¿П ^ 1 ^ ¿¿Р ^ " ¿ ■ ■ п т
Видно, что теперь, в отличие от однородного уравнения, для выполнения шага корректор кроме предикторных величин /*+1/2, и*+1/2 в полуцелых узлах требуются также величины и* в целых узлах. Для их вычисления будем использовать следующий аналог разностного уравнения (45):
и* _ ип /г \
^ U'¿ , / /9 - х*и9 \__/ п .п ~п
тп
+ ^/9 ^^ = ^ (хп,^п, Йп) ,
48)
где
ип. 1 + ип- /п.-- /п
~п = ¿-1 ^ "¿+1 /п = •/¿+1 -/¿-1 ип = "¿+1 "¿-1 (49)
и = 2 , /9'п = 2Н , U9'¿ = 2Н , (49)
а функции вп определим по формуле вида (17), скорректированной с учётом подвижности сетки (35), при этом величины дп вычисляются в целых узлах
0п = (|а|(1 - СК)п, Сп = *(у)п , в^ = ё^п - 1.
Кроме того, как и в случае однородного уравнения, вместо вычисления двух величин / * и и * в полуцелых узлах можно вычислять только одну /, используя формулу вида (32)
(^) ¿+1/2 + (7и')п+1/2 = (а»>п+1/2
и переходя на шаге корректор от (47) к уравнению вида (33)
( 7и)п+1 _ ( 7и)п ( / - (х*и)п ) , - ( / - (х*и)п ) ,
-^ + "-¿+1/2 ^-= (^)*. (50)
тн
Если в = О(Н), то схема предиктор-корректор (48)-(50) для неоднородного уравнения (2) имеет второй порядок аппроксимации.
На рис. 3, а изображены графики численного решения начально-краевой задачи для неоднородного невязкого уравнения Бюргерса с правой частью
1 - 2и _2 (х - ¿/2 - х0
= сь- ---
п
а б
Рис. 3. Профили численного решения задачи о движении сглаженного скачка (51), полученные в момент времени Ь = 20 по схеме предиктор-корректор на равномерной (1) и адаптивной (2) сетках, по противопоточной схеме на равномерной сетке (3) (а); траектории узлов адаптивной сетки (б)
Точное решение задачи
- 2 * () (51)
представляет собой сглаженную ступеньку, движущуюся без изменения своей формы вправо со скоростью 1/2. Начальная функция задана формулой (51) при £ = 0 с точкой перегиба в точке х0 = 10 и значением V = 10-2. Все остальные входные данные взяты такими же, как в предыдущих тестах, за исключением двух параметров, влияющих на степень сгущения узлов и уменьшенных по причине гладкости решения до значений а\ = 5, а = 5. Из рис. 3, а видно, что численное решение на адаптивной сетке (линия 2) имеет неоспоримое преимущество по точности перед решением (линия 1), полученным по схеме предиктор-корректор на равномерной сетке с тем же числом узлов, которое в свою очередь точнее полученного с использованием противопоточной схемы (линия 3). В рассматриваемом примере высокая точность метода адаптивных сеток объясняется сильной концентрацией узлов сетки в окрестности подвижной точки перегиба, траектория движения которой показана на рис. 3, б штриховой линией.
1.4. О росте числа экстремумов при использовании ХУЮ-схем
Как было указано выше, некоторые ТУВ-схемы, сохраняющие монотонность численного решения, могут увеличивать количество экстремумов численного решения при переходе с одного временного шага к другому. Последнее является одной из причин появления осцилляций разностных производных, и при построении адаптивной сетки использование управляющей функции, зависящей от этих производных (например, управляющей функции (42)), приведёт к потере точности численного решения. Таким образом, для успешного применения схемы предиктор-корректор, которая является
ТУВ-схемой, необходимо знать, сохраняет ли она число экстремумов при переходе с одного временного шага на другой.
Сначала на примере схемы Лакса (13) для невязкого уравнения Бюргерса, при условии (11) являющейся ТУО-схемой, покажем, что количество экстремумов действительно может увеличиться на следующем шаге по времени. Для этого возьмём сеточную функцию с одним экстремумом — максимумом в узле с номером 3 = 0:
иП =
если если
3 = 0, 3 = 0,
(52)
и вычислим значения решения на (п + 1)-м временном слое:
0, если 3 < -2,
1/2 — ж/4, если 3 = —1,
ит " = < 0, если 3 = 0,
1/2 + ж/4, если 3 = 1,
0, если 3 > 2.
п+1
(53)
Легко проверить, что при выполнении неравенства (11) функция (53) имеет большее, чем ип, количество экстремумов: два локальных максимума и один локальный минимум. При этом, несмотря на рост количества экстремумов, полная вариация сеточного решения не может увеличиться, поскольку рассматриваемая схема обладает ТУВ-свой-ством:
ТУ (ига+1) =
Очевидно, что кроме схемы Лакса существуют и другие ТУВ-схемы, увеличивающие количество экстремумов (см., например, [12, 19]).
При том же условии (11) противопоточная схема (12) для невязкого уравнения Бюргерса, также относящаяся к классу ТУВ-схем, дает на (п + 1)-м слое по времени функцию с одним экстремумом
е1
— ип+1| = 2 = ТУ (ип).
и
п+1
0, если 3 <—1,
1 — ж/2, если 3 = 0,
ж/2, если 3 = 1,
0, если 3 > 2,
т.е. количество экстремумов не увеличивается, при этом полная вариация в этой схеме уменьшается: ТУ(ига+1) = 2 — ж < ТУ(ип). Такой же результат получается при использовании ТУВ-схемы предиктор-корректор (3), (4), (17), если коррекция вязкости (25) не используется. В противном случае также получаем решение с одним экстремумом
и?+1 = <
0,
ж(—3ж2 + ж + 2)/8, 1 + ж(3ж2 — ж — 6)/8, ж/2, 0,
если если если если если
3 < —2, 3 = —1, 3 = 0, 3 = 1, 3 > 2,
но с другим значением полной вариации
ж
ТУ(ига+1) = 2 + ^(3ж2 — ж — 6) < 2 — ж < ТУ(ип).
Отметим, что в работе [15] указан признак, позволяющий определить для заданной функции (52) сохраняет ли или увеличивает количество экстремумов ТУВ-схема, аппроксимирующая линейное уравнение переноса с постоянным коэффициентом. Достаточно полное исследование эволюции числа экстремумов при использовании схем, аппроксимирующих нелинейное скалярное уравнение (1), выполнено в [19].
2. Схема предиктор-корректор для уравнений мелкой воды
Система уравнений мелкой воды, описывающая течение несжимаемой жидкости со свободной границей, имеет следующий вид:
и + ^ = О, (54)
где
и =( Ни), * (и)=( Ни 2 +ия2/0 , О(Х, и)=( я!) ,
£ — время, и(х,£) — скорость, п(х, £) — отклонение свободной поверхности от невозмущённого уровня, Н(х) — глубина дна, отсчитываемая от невозмущённой свободной границы, Я = п + Н — полная глубина. Уравнения (54) записаны в безразмерных переменных, в которых ускорение свободного падения равно единице.
В новых координатах (д, £) (см. раздел 1.2) это уравнение можно записать как в виде дивергентного (44)
(7 и)4 + ^ — я4и)д = О, (55)
так и в виде недивергентных (43)
и + 1 (^ — ж4ид) = 10, {t + 1а (^ — ж4ид) = 1 аО (56)
уравнений, где через вектор О обозначена правая часть уравнения (54), представленная в новых координатах и умноженная на якобиан 7, а через а — матрица Якоби
и)=( Я°0 , а(и) = £(и)=( —и2 0+ Я 2и).
Отметим, что в координатах д, £ дно является, вообще говоря, подвижным, хотя в исходных переменных зависимость функции Н от времени не предполагалась. Этот факт следует из равенства Н(х) = Н(х(д,£)) = Н(д, £).
На шаге предиктор используются аппроксимации уравнений (56), предварительно переписанных в характеристической форме. Она получается следующим образом. Если использовать матрицы
l = 1( -Л2 m r = щ -1 m л=с Л 0 l Е\ -Al 1 ) ' r 2 V -Al Л2 J ' л ^ 0 Л2
где Л1 = u — л/И, Л2 = u — собственные значения матрицы a, то справедливо
представление a = ral, при этом rl = lr = e — единичная матрица. Умножим уравнения (56) слева на матрицу l
lut + 1l (fq — xtUq) = 1lG, lft + 1л£ (f, — xtu,) = -1^lG
и аппроксимируем полученные уравнения аналогично (45), (46):
^^^ ^^ + (I£ (,, - х,и„))П+1/2 = ( 1 (р-.£)«+1/2 £+1 + (1Л£ (Г. - х. и,) У = (1л£С)
(57)
-Р* _ -СП / л \ и / 1
11+1 /2 11+1 /2 ч\ /1
- х. ид)
'1+1/2 /1+1/2
где
ига + ига Р + 1П / 0
„и _ и1 + и1+1 пи _ 11 + 11+1 рга п / 1га) /-.га / 0
/о --I , О --I ' - 1 (и I ' /о - I /ТТ7 \и
1+1/2 = 2 ' 11+1/2 = 2 ' 11 = 1 Vй?/ ' О1+1/2 = I (ЯЛ,д)
^1+1/2
Ги = _4_ / -Ли,1+1/2 1
¿1+1/2 (ли Ли )2 \ _Ли 1
(Л2,1+1/2 - Л1,1+1/2) V Л1,1+1/2 1
Ли 1+1/2 (к = 1, 2) — собственные значения матрицы аи+1/2, приближающей а, Ли = ( Ли,1+1/2 0 ) ри = ( 1 + ^)+1/2 0
1/2 = 1 0 лп^; ' 1/2 Ч 0 1+^
91+1/2 (к = 1, 2) — схемные параметры.
Для дальнейших преобразований уравнений (57), (58) важным является способ аппроксимации матрицы а. Аналогично скалярному случаю (31) будем требовать сохранения равенства 1 я = аид на разностном уровне, т. е. выполнения равенства
1 П,1+1/2 = (аи,)П+1 /2 . (59)
Этому условию удовлетворяет, например, следующая матрица:
01
1 "1+1 + Н1+1/2 2и1+1/2
с собственными значениями
АП+1/2 = ( -"П"П+1 + НП+1/2 2м"+1 /2 ) = (^Л£)П+1/2 (60)
ЛП,1+1/2 = "П+1/2 («и+1/2)2 - «и"и+1 + НП+1/2' к =1' 2' (61)
при этом
НП+1/2 = (НП + НП+1)/2' "П+1/2 = ("и + "П+1)/2'
ии
,и = (р и )-1 = Л2,1+1/2 Л1,1+1/2 / -1 1
1+1/2 = 1^1+1/^ = 4 I _ли , ли
4 \ Л1 1+1 /2 Л2 1+
Используя выражение (60) и новые обозначения
ди = лп _ хП р рП = (ги )П
Л1+1/2 = Л1+1/2 Х.,1+1/2Е ' Р1+1/2 = )1+1/2
получаем формулы для предикторных величин и* и 1*
и*+1/2
Т ( 1
и - 2(1 ^(ЛР -¿О)
1+1/2
*
11+1/2
Т ( 1
1 - 2 (1 ЯрЛ(ЛР - ¿О)
1+1/2
и
и
которые применяются на шаге корректор в разностном уравнении вида (47)
(7и)П+1 — (7и)" + (Г — Х«и*) -+1/2 — (Г — Х'и*)1/2 = 0* (62)
Т + Н = 0" (62)
аппроксимирующем дифференциальное уравнение (55), при этом
0* = ( 0
В уравнении (62) величины и * и ^* используются в комбинации ^* — Ж'и *. Аналогично случаю скалярного нелинейного уравнения легко проверить, что
* — ж'и * )0+1/2 = {0+1/2 — (ж'и)"+1 /2,
где
^ 0+1/2
Т / 1 - -f — ^1 ^РЛ(Л Р — £0)
. (63)
0+1/2
Итак, схема предиктор-корректор состоит из шага предиктор (63) и шага корректор
(7и)Т+' — (7и)" + — <Х'и>^,+,/2 — - <Х'и>");-,/2 _ _*
Т + Н =с*'.
Для замыкания этой схемы необходимо указать способ вычисления величин Я*, Н^ и функций 0*+1/2. Функции ^ будем определять по формуле, аналогичной (17), но теперь входящие в эту формулу величины должны учитывать наличие двух собственных значений (61) и подвижность сетки аналогично тому, как это сделано в формуле (35) для схемы, аппроксимирующей нелинейное скалярное уравнение:
Г 0 пРи 1^+1/21 < кк+1/2-^к1 и ¿+1/2 ■ ¿+1/2-*к > 0,
^0+1/2 = < 90к0+1/2(1 — ек+1/2) пРи |^.?+1/2| > кк+1/2-^к1 и ¿+1/2 ' ¿+1/2-*к > 0, (64)
1 <0+1/2 пРи ¿+1/2 ■ ¿+1/2-^ < 0,
где
к = 1, 2, £?+1/2 = +1/2-^ /#7+1/2,
„А; _ л \(Л г< )" __I |Ак ^ - 1 /Й
I < 1, , с 0+1/2 Ск,0+1/2
р" 0+1/2 — компоненты вектора Р"+1/2, А" 0+1/2 — диагональные элементы матрицы ЛЛ"+1/2,
= А",0+1/2.
Для вычисления предикторной величины Я* применяется [26] аппроксимация в целых узлах ^ первого уравнения (56):
1 + 01
¿+1/2 = (|Ак 1(1 — СА)рк)"+1 /2, С",0+1/2 = « ( ЦЧ < 1, ^к,0+1/2 = С"--1
V 7 / 7+1/2 Ск,0+1/2
тт *
Я — Н 3 +
т/2
7 (Л2 — А1)
(А^Я, — А^Яи), + Я Н,) —
1 +
7(Л2 — А^ А1А2Я А2(Яи)9 '
(а^Я, — А2 (Яи)^ + Я Н,)
"
= 0, (65)
0
где, как и при получении уравнения (48), для определения всех разностных производных используются центральные разности и обозначения (49), функции вычисляются по формуле (64), но по компонентам вектора ри в узлах с целыми индексами. Для вычисления производной ^ применяется формула
h* — h* 1
hV = J+1/2h J-1/2, h**±1/2 = 4 [h(x™) + h(xn±1) + h(xn+1) + hj1)] .
Схема предиктор-корректор (63)-(65) на подвижной сетке является аналогом схемы (48)-(50) для неоднородного скалярного уравнения (2). В отличие от скалярного случая для нелинейных уравнений мелкой воды не удаётся строго обосновать свойство сохранения монотонности численного решения. Поэтому особую важность приобретают исследование свойств схемы и её численное тестирование.
Уравнения (54) сохраняют состояние покоя жидкости, т. е. если в начальный момент времени п(Х' 0) = 0, и(Х' 0) = 0, то и при £ > 0 система уравнений (54) имеет решение Н(Х'£) = Л(х), "(Х'£) = 0, соответствующее покою жидкости с невозмущённой свободной границей п(х'£) = 0. Это свойство желательно и для разностной схемы, поскольку в случае неровного дна уравнения мелкой воды пригодны для описания волн лишь небольшой амплитуды, в силу чего при наложении на такие волны значительных возмущений, вызванных неподходящей аппроксимацией, картина течения может заметно исказиться. Справедливы следующие леммы [26], касающиеся свойств рассматриваемой схемы на равномерной неподвижной сетке (х. = 0, 3 = /).
Лемма 4. Схема предиктор-корректор (63)-(65) сохраняет состояние покоя жид-n = 0, Hn = h(x), un
кости: если = 0, Hn = h(x), un = 0, то и ПП+1 = 0, Hn+1 = h(x), un+1 = 0.
Лемма 5. Для плоского дна h(x) = h0 = const схема предиктор-корректор (63)-
j = H0 = const, un
(65) сохраняет постоянное течение жидкости: если H™ = H0 = const, и™ = u0 =
const, то и Щ+1 = H0, «"+1 = u0.
В случае плоского горизонтального дна нелинейные уравнения мелкой воды (54) имеют следующее разрывное решение, описывающее движение устойчивого гидравлического скачка, обращённого вправо:
H(xt)=J H пРи x ^ xo + Wt u(xt) i U пРи x < xo + Wt, (66)
H (X,t) \ Hr при x>x0 + Wt, U(X,t) \ Ur при x>x0 + Wt, ( ) где Hi > Hr, W — скорость движения скачка,
и = иг + (Н - Нг VН+Нг ■ = иг + /Н .
На скачке выполняется соотношение
Ж(Нг - Нг) = НгЦ - Нгиг.
Скачок является стационарным, если Ж = 0. В этом случае НгЦ = Н^.
Лемма 6. Схема предиктор-корректор (63)-(65) сохраняет стационарный скачок (66).
Указанное в лемме 6 свойство схемы остаётся справедливым и при других аппроксимациях матрицы а, удовлетворяющих условию (59), например, когда в качестве ап+1/2 берется матрица Роу [27].
Применим схему предиктор-корректор для решения задачи о разрушении плотины. Предполагается, что дно является плоским и горизонтальным (Н(х) = 1), в начальный момент времени жидкость покоится, т. е. и(х, 0) = 0, а полная глубина имеет различные постоянные значения слева и справа от некоторой точки х0, 0 < х0 < 1:
Н (х, 0)
Н1, если х < х0, Н2, если х > х0.
Далее предполагается, что Н1 > Н2. В этом случае при £ > 0 возникает течение с двумя волнами: волной понижения, движущейся влево, и бором, движущимся вправо. Между волной понижения и бором возникает течение с постоянными скоростью и0 и глубиной Н0, которые подлежат определению. Точное решение поставленной задачи имеет вид
Н (х,£)
Н1, Н0,
Н2,
х — х0 £
если х < х1(£),
если х1(£) < х < х2(£),
если х2 (£) < х < хь(£), если х > хь(£),
и(х, £)
где
0,
2
л/НГ+
х — х0
£
Ц), 0,
если х < х1(£),
если х1(£) < х < х2(£),
если х2(£) < х < хь(£), если х > хь(£),
х1(£) = х0 — ¿\/Нъ х2(£) = х0 + £(2у/НТ — 3\/Н0), хь(£) = х0 + Ц/ ^
и0 = — VН0), Н0 является корнем уравнения
Н0 Н0 + Н2
Н2
2
(Н — Н2)у Н^у+Н + = 2УН1.
Ниже представлены результаты численного решения задачи, предложенной в [28] и характеризующейся большим начальным перепадом полной глубины: Н1 = 15, Н2 = 1. Расчёт проводился в области (0, 2) до момента времени £ = 0.15, при х0 = 1, N = 100. Шаг по времени определялся по формуле
кз
Н
тах I |А
0<7<М-1, к=1,2 *
п
й,0+1/2|
К Л +
0+1/2
где кзап — коэффициент запаса, кз
0.8.
2
Рис. 4. Профили полной глубины в момент времени Ь = 0.15 (а); траектории узлов адаптивной сетки (б)
На рис. 4, а представлены графики полной глубины в конечный момент времени при использовании схемы предиктор-корректор. Видно: как на адаптивной (штриховая линия), так и на равномерной (сплошная линия) сетке осцилляции в численном решении не возникают, хорошо передаются положение фронта бора и профиль волны понижения. Анализ увеличенных изображений показывает, что график полной глубины, полученный на адаптивной сетке, лучше аппроксимирует график точного решения (тонкая линия на рис. 4, а) и метод адаптивных сеток даёт лучшую разрешимость решения в окрестности разрыва, чем схема на равномерной сетке, при использовании которой получаются более размытые профили.
Отметим, что адаптивная сетка строилась при £ > 0 методом эквираспределения [15] с параметром в = 20 и управляющей функцией вида (42), в которой регулировка сгущения узлов выполнялась по градиенту функции п с использованием параметров а\ = 0.1 и а =15. Из рис. 4, б видно, что сетка имеет подвижные сгущения узлов в окрестности фронта бора, траектория движения которого изображена штриховой линией, и непосредственно в волне разрежения, подвижные границы которой также показаны штриховыми линиями. В то же время на участках постоянства решения сетка получается разрежённой.
На рис. 4, а сплошными линиями приведены два графика полной глубины, полученные на равномерной сетке. Один из них имеет нефизичный стационарный скачок, расположенный в точке первоначального разрыва х0 = 1, в которой собственное значение А1 обращается в нуль. Появление таких нефизичных скачков в волне разрежения при использовании различных разностных схем является известным фактом [4]. Другой график не имеет такого скачка. В последнем случае применялась коррекция вязкости по формуле (36) с заменой Л на Л1.
Заключение
Для однородного нелинейного скалярного уравнения выписана схема предиктор-корректор на равномерной неподвижной и неравномерной подвижной сетках. Выбор схемного
параметра, при котором схема сохраняют монотонность численного решения, и энтропийная коррекция проведены на основе исследования дифференциального приближения схемы. Рассмотрены свойства сохранения схемой постоянного решения и разрывных решений задачи Римана в виде стационарного скачка на неподвижной сетке и движущегося скачка на подвижной сетке. Предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках. Численно решена начально-краевая задача с непрерывной начальной функцией для невязкого уравнения Бюргерса, где точным решением является центрированная волна сжатия, в момент градиентной катастрофы генерирующая в зависимости от значений входных параметров стационарный или движущийся скачок. На примере численного решения задача Римана с точным решением в виде волны разрежения показано преимущество использования энтропийной коррекции для схемы предиктор-корректор. Для неоднородного нелинейного скалярного уравнения выписана схема предиктор-корректор на неравномерной подвижной сетке и представлены результаты численного решения начально-краевой задачи с известным точным решением для неоднородного невязкого уравнения Бюргерса. Для всех решённых задач продемонстрировано преимущество использования адаптивных сеток. На примере схемы Лакса для невязкого уравнения Бюргерса показано, что TVD-схемы могут увеличивать количество экстремумов при переходе с одного шага по времени на другой. Противопоточная схема и схема предиктор-корректор количество экстремумов не увеличивают.
Для одномерных уравнений мелкой воды выписана схема предиктор-корректор на неравномерной подвижной сетке. Так как строго обосновать свойство сохранения монотонности численного решения здесь не удаётся, то проведены исследование свойств схемы и её численное тестирование. Отмечены такие свойства схемы на равномерной неподвижной сетке как сохранение состояния покоя и в случае плоского горизонтального дна — сохранение постоянного течения жидкости и разрывного решения в виде стационарного скачка. На примере численного решения задачи о разрушении плотины на плоском горизонтальном дне продемонстрированы преимущества использования подвижных адаптивных сеток и энтропийной коррекции схемы.
Список литературы
[1] Oleinik O. Discontinuous solutions of nonlinear differential equations // Amer. Math. Soc. Transl. 1957. Ser. 2. Vol. 26. P. 95-172.
[2] Lax P. Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves. Society for Industrial and Appl. Math. Philadelphia, 1973.
[3] Smoller J. Shock Waves and Reaction-Diffusion Equations. Series of Comprehensive Studies in Mathematics. Vol. 258. New-York: Springer-Verlag, 1983.
[4] LeVeque R.J. Numerical Methods for Conservation Laws. Basel, Boston, Berlin: Birkhäuser Verlag, 2008.
[5] Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Berlin, New-York: Springer-Verlag, 2009.
[6] Harten A., Hyman J.M. Self-adjusting grid methods for one-dimensional hyperbolic conservation laws // J. of Comput. Phys. 1983. Vol. 50. P. 235-269.
[7] Harten A. High resolution schemes for hyperbolic conservation laws // Ibid. 1983. Vol. 49. P. 357-393.
[8] Roe P.L., Pike J. Efficient Construction and Utilisation of Approximate Riemann Solutions. Els. Sci. Publ. B. V. (North-Holland). INRIA, 1984.
[9] Roe P. L. Self-adjusting grid methods for one-dimensional hyperbolic conservation laws // SIAM J. Sci. Stat. Comput. 1992. Vol. 13, No. 2. P. 611-630.
[10] Dubois F., Mehlman G. A non-parameterized entropy correction for Roe's approximate Riemann solver // Numer. Math. 1996. Vol. 73. P. 169-208.
[11] Osher T.S. Riemann solvers, the entropy condition, and difference approximation // SIAM J. Numer. Anal. 1984. Vol. 21, No. 2. P. 217-235.
[12] Tadmor E. Numerical viscosity and the entropy condition for conservative difference schemes // Math. Comput. 1984. Vol. 43, No. 168. P. 369-381.
[13] Tadmor E. The numerical viscosity of entropy stable schemes for systems of conservation laws // Ibid. 1987. Vol. 49, No. 179. P. 91-103.
[14] Tadmor E. Entropy stability theory for difference approximations of nonlinear conservation laws and related time dependent problems // Acta Numerica. 2003. Vol. 12. P. 451-512.
[15] Хлкимзянов Г.С., Шокинл Н.Ю. Некоторые замечания о схемах, сохраняющих монотонность численного решения // Вычисл. технологии. 2012. Т. 17, № 2. С. 79-98.
[16] Breuss M. An analysis of the influence of data extrema on some first and second order central approximations of hyperbolic conservation laws // ESAIM: Math. Model. Numer. Anal. 2005. Vol. 39, No. 5. P. 965-993.
[17] Breuss M. The correct use of the Lax-Friedrichs method // Ibid. 2004. Vol. 38, No. 3. P. 519-540.
[18] Tang H., Warnecke G. A note on (2k + 1)-point conservative monotone schemes // Ibid. 2004. Vol. 38. P. 345-358.
[19] LeFloch P.G., Liu J.-G. Generalized monotone schemes, discrete paths of extrema, and discrete entropy conditions // Math. Comput. 1998. Vol. 68, No. 168. P. 1025-1055.
[20] Яушев И.К. О численном расчёте нестационарных течений газа в одномерном приближении в каналах со скачком площади сечения // Изв. СО АН СССР. Техн. науки. 1967. Т. 8, № 2. С. 39-48.
[21] Хлкимзянов Г.С., Шокин Ю.И., Блрлхнин В.Б., Шокинл Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.
[22] Murman E.M. Analysis of embedded shock waves calculated by relaxation methods // AIAA J. 1974. Vol. 12. P. 626-633.
[23] Шокин Ю.И., Сергеевл Ю.В., Хлкимзянов Г.С. О монотонизации явной схемы предиктор-корректор // Вестник КазНУ. Математика, механика, информатика. 2005. Т. 2. Спец. выпуск. С. 103-114.
[24] Tadmor E. The large-time behavior of the scalar, genuinely nonlinear Lax-Friedrichs scheme // Math. Comput. 1984. Vol. 43, No. 168. P. 353-368.
[25] Блрлхнин В.Б., Клрлмышев В.Б., Бородкин Н.В. TVD-схема на подвижной адаптивной сетке // Вычисл. технологии. 2000. Т. 5, № 1. С. 19-30.
[26] Shokin Yu.I., Sergeeva Yu.V., Khakimzyanov G.S. Predictor-corrector scheme for the solution of shallow water equations // Rus. J. Numer. Anal. Math. Model. 2006. Vol. 21, No. 5. P. 459-479.
[27] Roe P.L. Approximate Riemann problem solvers, parameter vectors, and difference schemes // J. Comput. Phys. 1981. Vol. 43, No. 2. P. 357-372.
[28] Fjordholm U., Mishra S., Tadmor E. Energy preserving and energy stable schemes for the shallow water equations // Foundations of Computational Mathematics. Proc. of FoCM held in Hong Kong 2008 / Eds. F. Cucker, A. Pinkus, M. Todd. London Math. Soc. Lecture Notes Ser. 2009. Vol. 363. P. 93-139.
Поступила в 'редакцию 4 апреля 2013 г.
