Метод адаптивных сеток для одномерных уравнений мелкой воды Текст научной статьи по специальности «Математика»

Аннотация научной статьи по математике, автор научной работы — Хакимзянов Гаяз Салимович, Шокина Нина Юрьевна

На примерах для нелинейного скалярного уравнения обсуждаются вопросы построения разностных схем, сохраняющих монотонность численного решения, на равномерных и адаптивных сетках . Показаны важные свойства сохранения постоянного решения, стационарного и движущегося скачков. С помощью метода дифференциального приближения дано новое объяснение механизма возникновения нефизичных численных решений и проведена энтропийная коррекция . Представлен пример TVD-схемы, увеличивающей количество экстремумов. Предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках. Схема предиктор-корректор , построенная для нелинейного скалярного уравнения , применена для решения одномерных нестационарных уравнений мелкой воды. Исследованы свойства схемы и проведено её численное тестирование.

Похожие темы научных работ по математике , автор научной работы — Хакимзянов Гаяз Салимович, Шокина Нина Юрьевна

Adaptive grid method for one-dimensional shallow water equations

The questions of construction of difference schemes, which preserve monotonicity of a numerical solution, on uniform and adaptive grids are discussed using the examples for nonlinear scalar equation. Important properties of preservation of constant solution, stationary and moving shocks are shown. A new explanation of generation of nonphysical numerical solutions and the entropy fix are done using the method of differential approximation. An example of a TVD scheme, which increases a number of extremums, is given. A new approach to construction of any explicit two-layer divergence schemes on moving grids is suggested. The predictor-corrector scheme is used for numerical solution of one-dimensional non-stationary shallow water equations. Investigation of scheme properties has been done combined with its numerical testing.

Текст научной работы на тему «Метод адаптивных сеток для одномерных уравнений мелкой воды»

Метод адаптивных сеток для одномерных уравнении мелкой воды*

Г. С. Хлкимзянов, Н.Ю. ШокинА Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: khak@ict.nsc.ru, nina.shokina@ict.nsc.ru

На примерах для нелинейного скалярного уравнения обсуждаются вопросы построения разностных схем, сохраняющих монотонность численного решения, на равномерных и адаптивных сетках. Показаны важные свойства сохранения постоянного решения, стационарного и движущегося скачков. С помощью метода дифференциального приближения дано новое объяснение механизма возникновения нефизичных численных решений и проведена энтропийная коррекция. Представлен пример TVD-схемы, увеличивающей количество экстремумов. Предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках. Схема предиктор-корректор, построенная для нелинейного скалярного уравнения, применена для решения одномерных нестационарных уравнений мелкой воды. Исследованы свойства схемы и проведено её численное тестирование.

Ключевые слова: численное моделирование, конечно-разностная схема, нелинейное скалярное уравнение, предиктор-корректор, монотонная схема, дивергентная схема, энтропийная коррекция, адаптивная сетка, метод эквираспределения, нелинейные уравнения мелкой воды.

Достоверная оценка зон затопления при накате волн цунами на берег является актуальной проблемой математического моделирования. Для получения эффективного прогноза необходимы, кроме качественных, исследованных алгоритмов, большие массивы подробной информации о рельефе дна и суши в прибрежной зоне, о параметрах источника волны, данные наблюдений для верификации и т. д. Как правило, вся эта входная информация весьма неточна. Кроме того, модель мелкой воды, используемая для моделирования процесса наката волн цунами на берег, выводится при достаточно жёстких ограничениях. Поэтому для получения удовлетворительных результатов в указанных задачах использование схем повышенного порядка точности не имеет смысла и достаточны схемы первого или второго порядка. Это определяет необходимость построения качественных дивергентных схем первого или второго порядка аппроксимации, точных на возможно широком классе решений исходной задачи, в частности, сохраняющих состояние покоя, монотонность либо решения, либо инвариантов Рима-на и дающих численное решение, имеющее наиболее широкий набор свойств решения исходной задачи. Нелинейное скалярное уравнение является модельным для системы уравнений мелкой воды, поэтому оно использовано в первой части работы для исследования свойств конструируемых разностных схем, применяемых далее для решения задач о течениях в рамках модели мелкой воды.

* Работа выполнена при финансовой поддержке РФФИ (грант № 12-01-00721а), а также в рамках Программы интеграционных исследований СО РАН (проект № 42).

При решении задач с разрывными решениями нарушение условия неубывания энтропии [1 -3] приводит к нефизичному численному решению в виде ударной волны разрежения в задачах о трансзвуковых течениях идеального газа, а также в окрестности критической точки, расположенной в транскритической волне разрежения в задачах о течениях идеальной жидкости в рамках модели мелкой воды [4, 5]. При численном решении этих задач необходимо, чтобы выбранный метод сходился к физически верному (энтропийному) решению, что гарантируется выполнением дискретного аналога условия неубывания энтропии. Для выполнения данного условия можно использовать известный приём энтропийной коррекции, который был предложен в работах [6, 7] и применялся в [4, 8-10] для метода Роу ив [11-14] для других численных методов. В настоящей работе показаны важные свойства сохранения схемой предиктор-корректор постоянного решения и стационарного или, в случае подвижной сетки, движущегося скачка. С помощью метода дифференциального приближения дано новое объяснение механизма возникновения нефизичных численных решений и проведена энтропийная коррекция схемы.

Использование адаптивных подвижных сеток для схем первого или второго порядка аппроксимации повышает точность численного решения, но главная причина их применения в задачах моделирования процесса наката волн цунами на берег состоит в том, что использование разностных схем в методах с выделением подвижной линии уреза в принципе невозможно без подвижных сеток. Кроме того, в двумерном случае необходимо учитывать криволинейность границы расчётной области. Поэтому в настоящей работе большое внимание уделяется описанию особенностей моделирования на подвижных сетках. В частности, предложен новый подход к построению явных двухслойных дивергентных схем на подвижных сетках (аналогичный подход для линейного уравнения переноса рассмотрен в [15]).

Если при использовании адаптивных сеток управляющая функция, регулирующая расстановку узлов сетки, зависит от разностных производных, которые осциллируют, то сетка будет иметь чередование длинных и коротких ячеек, что приведёт к потере точности численного решения. Одна из причин появления осцилляций разностных производных — увеличение схемой количества экстремумов решения при переходе с одного шага по времени на другой. Известно [7], что если начальная функция имеет ограниченную вариацию, то количество локальных экстремумов слабого решения не возрастает. Вместе с тем некоторые ТУВ-схемы, сохраняющие монотонность численного решения, могут увеличивать количество экстремумов решения [15-19]. В настоящей работе этому вопросу уделено внимание на примере схем Лакса, противопоточной и предиктор-корректор для невязкого уравнения Бюргерса.

Во второй части работы схема предиктор-корректор, построенная для нелинейного скалярного уравнения, применена для решения одномерных уравнений мелкой воды. Исследованы свойства схемы и проведено её численное тестирование.

1. Разностные схемы для нелинейного уравнения, сохраняющие монотонность численного решения

Данный раздел посвящён разностным схемам для однородного

нелинейных скалярных уравнений, которые являются модельными для системы уравнений мелкой воды в случае ровного горизонтального и неровного дна соответственно.

1.1. Схема предиктор-корректор на равномерной сетке

Сначала рассмотрим случай равномерной неподвижной сетки с узлами ж. = и ша-

гом Л > 0. В явной схеме предиктор-корректор [20, 21] искомая сеточная функция и™ определённая в момент времени ¿п в узлах ж., находится из разностного уравнения

аппроксимирующего закон сохранения (1). Необходимые для выполнения этого шага потоки /*, отнесённые к полуцелым узлам ж^+1/2 = ж. + Л/2, определяются на шаге предиктор

где /. = / (ип), т*+1/2 = 0.5т(1 + $п+1/2), т — шаг по времени, ^^+1/2 — схемный параметр. Разностное уравнение (4) есть результат аппроксимации дифференциального уравнения для потоков

которое получается после умножения (1) на функцию а (и) = /и(и).

Наличие некоторых важных свойств схемы (3), (4) напрямую связано с видом сеточной функции ап+1/2, аппроксимирующей функцию а(и). В настоящей работе способ вычисления функции ап+1/2 определяется требованием сохранения для сеточных функций правила дифференцирования /х = а(м)мх функций непрерывного аргумента, т. е. выполнением равенства

для разностных производных

Условию (6) удовлетворяет, например, сеточная функция [7]

ап I 71+1 Мп пРи "п+1 = " а(ип) при мп+1 = мп,

применение которой позволяет переписать уравнение (4) с использованием разностных производных от искомого решения

и получить следующий одношаговый вариант схемы (3), (4):

((1 + 0)аЧ)П 2 - ((1 + 0)а2ил)

Схему (9) можно назвать канонической формой явных дивергентных двухслойных схем для уравнения (1), поскольку любая такая схема может быть записана в виде (9) при соответствующем выборе параметра 0. Например, при ^*П+1/2 — 0 схема (9) совпадает со схемой Лакса — Вендроффа

где Сп+1/2 — число Куранта,

получается противопоточная схема первого порядка аппроксимации, которая может быть записана как в дивергентной, так и в недивергентной форме [22]

о иж I + —о— иж = 0,

— 1 > 0 — схема Лакса

при 0°+1/2 — -1 — абсолютно неустойчивая схема второго порядка аппроксимации по Н

ип+1 _ ип /п _ /п и_и! + /^'+1 /^-1 = 0

Явные схемы с несимметричными шаблонами также получаются из канонической формы (9). Например, при а > 0 и

схема (9) переходит в противопоточную схему второго порядка аппроксимации на несимметричном шаблоне

и™+1 - ип 3 (аих)°—1/2 - (аиж)°-3/2 _ Т (аиж)°-1/2 - (аих)п-3/2

Если 9 = 0(h), то схема (9) аппроксимирует уравнение (1) со вторым порядком относительно т и h. При 9 = const = 0 порядок аппроксимации на гладких решениях — первый. Необходимое условие устойчивости для случая a(u) = const [15] требует, в частности, чтобы параметр 9 принимал неотрицательные значения. Поэтому далее предполагается выполненным ограничение

В общем случае (для произвольного параметра 9 = 0(h)) схема предиктор-корректор (3), (4) не сохраняет монотонность численного решения. Однако этот параметр можно подобрать так, что схема будет обладать указанным свойством. В работе [23] для выбора параметра 9 использовалось первое дифференциальное приближение (п.д.п.) схемы предиктор-корректор

ut + fx = + y [a(a2œ2 - 1)«x]xx (15)

и параметр 9 подбирался так, чтобы диссипативный член п.д.п. приводил к изменению дисперсионного члена в тех подобластях, где возможно появление осцилляций в численном решении. Например, если положить

то из (15) получается уравнение с дисперсионным членом

— (1 - |a|œ) (2 - |a|œ) Uxxx, 6

который в силу условия (11) отличается по знаку, например, от дисперсионного члена

в п.д.п. схемы Лакса — Вендроффа, не сохраняющей монотонность численного решения. В [23] при некоторых дополнительных ограничениях на функцию (16) получена формула для сеточной функции 9n+1/2

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

и доказано, что при условии (11) схема предиктор-корректор (3), (4), (17) сохраняет монотонность численного решения.

Отметим два свойства схемы (3), (4), связанные с её точными решениями.

Лемма 1. Схема предиктор-корректор (3), (4) сохраняет постоянное решение: если u™ = U0 = const, то и u°+1 = U0.

Второе свойство связано с сохранением разностной схемой разрывных решений задачи Римана для уравнения (1) с начальной функцией

и°(х) = I Ц, Х < х°: Ц = Ц. (18)

Пусть функция потока является выпуклой, т.е. /"(и) > 0 (для вогнутой функции потока все дальнейшие рассуждения справедливы при противоположных соотношениях между Ц и Ц). Тогда для Ц > Ц. задача (1), (18) при £ > 0 имеет следующее решение в виде ударной волны:

где Ш — постоянная скорость движения скачка, определяемая из условия Ренкина — Гюгонио:

Если Ш = 0, то скачок является стационарным, и в этом случае

Решение (19) удовлетворяет условию неубывания энтропии [4, 5] и является физически корректным.

Известно [4, 5], что точное решение задачи Римана в виде стационарного скачка является и решением некоторых разностных схем. Покажем, что этим свойством обладают все схемы вида (9), причём не только второго порядка аппроксимации, но и первого, например, противопоточная (12).

Лемма 2. Схема предиктор-корректор (3), (4) сохраняет стационарный скачок (18), (21).

Доказательство. Пусть на п-м слое по времени (п > 0) скачок располагается между узлами х^0 и х^+1^, т.е.

'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)

Поскольку схема предиктор-корректор при Ц < Ц дает нефизичное решение невязкого уравнения Бюргерса, то, как указано выше, необходима её корректировка для возможности расчётов решений со сменой знака величины а. Энтропийная коррекция сводится в итоге к принудительному вводу в схему дополнительной вязкости в тех подобластях, где значение |а| становится меньше некоторого порогового и исчезает численная диссипация [24]. Покажем, что такой же приём может быть получен на основе метода дифференциального приближения, причём с явной формулой для порогового значения.

Раскрывая скобки в правой части п.д.п. (15), получим, что схемная вязкость ц (коэффициент при второй производной ихх) выражается формулой

ц = 29а2 + у (3а2ж2 - 1)ил. (24)

Следовательно, при малых значениях а и при > 0 схемная вязкость становится отрицательной, что и является причиной появления нефизичного численного решения при расчёте волн разрежения (23) с критической точкой а = 0. Отметим, что для невязкого уравнения Бюргерса в волнах сжатия (Ц > Ц) выполняется условие < 0, поэтому отрицательная вязкость в окрестности точки с нулевой критической скоростью а = 0 не возникает, а значит, не возникает нефизичных численных решений.

Из формулы (24) следует, что для обеспечения неотрицательности вязкости ц можно заменить функцию ($а2)™+1/2 на следующую:

(9а2)П+1/2 в других случаях, где функция 9 вычисляется по формуле (17),

С учётом коррекции схемной вязкости (25) шаг предиктор (8) запишется в виде

1.2. Схема предиктор-корректор на подвижной сетке

Рассмотрим некоторую начально-краевую задачу для уравнения (1) на отрезке [0, /]. Предполагается, что узлы хп подвижной неравномерной сетки, покрывающей в каждый момент времени этот отрезок, связаны взаимно-однозначно соответствием х = х(^,£)

с узлами qj равномерной неподвижной сетки с шагом Н = 1/Ж, покрывающей единичный отрезок [0, 1]. На подвижной сетке схема предиктор-корректор для скалярного уравнения (1) имеет следующий вид: