TVD-схема на подвижной адаптивной сетке Текст научной статьи по специальности «Математика»

Вычислительные технологии Том 5, № 1, 2000

TVD-СХЕМА НА ПОДВИЖНОЙ АДАПТИВНОЙ СЕТКЕ *

В. Б. БАРАХНИН, В. Б. КАРАМЫШЕВ Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected], [email protected]

Н.В. Бородкин Новосибирский государственный университет, Россия

A second order approximation scheme for dynamically adaptive grids which is the generalization of Harten scheme is presented. We expose the conditions which guarantee that the scheme is TVD. The conditions which are sufficient to ensure that this scheme is TVD have been studied. The numerical tests of the scheme and the comparison of its characteristics with that of the explicit predictor-corrector scheme with self-adjusting approximate viscosity have been conducted.

В последние годы при численном решении задач механики жидкостей и газов широко применяются подвижные сетки, адаптирующиеся к особенностям течения. Это позволяет получать результаты высокой точности при сравнительно небольших затратах машинных ресурсов. Однако для использования динамически адаптивных сеток аппроксимируемые уравнения необходимо записывать в подвижной системе координат, что неизбежно приводит к усложнению их вида.

Недостатком конечно-разностных схем, имеющих порядок аппроксимации выше первого, является осциллирующее поведение численного решения в областях резких изменений производных и локальных экстремумов искомых функций. Существует несколько подходов к монотонизации конечно-разностных схем. Один из них состоит в том, что в схему вводится искусственная вязкость (см., например, [5]). Преимуществом этого подхода является простота его реализации. Основной его недостаток заключается в том, что параметр, отвечающий за величину искусственной вязкости, подбирается опытным путем, при этом слишком большая вязкость приводит к уменьшению точности расчета и сглаживанию профилей решения. Ввиду сказанного искусственную вязкость желательно согласовывать со свойствами решения в тех или иных подобластях. Некоторые особенности применения на подвижных сетках явной схемы предиктор — корректор с автоматически настраиваемой аппроксимационной вязкостью описаны в работах [1, 6].

Другой подход к построению неосциллирующих алгоритмов заключается в использовании схем с шаблонами, перестраивающимися в зависимости от численного решения (TVD-,

* Работа выполнена при финансовой поддержке Программы интеграционных фундаментальных исследований СО РАН, проект №43, и Российского фонда фундаментальных исследований, гранты №96-0100136, 97-01-00819.

© В. Б. Барахнин, В. Б. Карамышев, Н. В. Бородкин, 2000.

TVB-, ENO-схемы и др.). Эти схемы часто используются в последние годы для расчетов как на равномерных, так и на неравномерных сетках, однако их аккуратное обоснование проводится, как правило, для случая скалярного закона сохранения на равномерной сетке.

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

В настоящей статье нами предлагается монотонная схема второго порядка аппроксимации, являющаяся обобщением известной схемы Хартена [7] для случая скалярного закона сохранения, записанного в подвижной системе координат. Дано ее обоснование и проведено тестирование путем сравнения ее свойств со свойствами явной схемы предиктор — корректор с автоматически настраиваемой аппроксимационной вязкостью.

1. Исходные уравнения

Рассмотрим скалярный закон сохранения, описываемый уравнением вида

ut + (f (u))x = 0, t> 0, (1.1)

с начальным условием

u(x, 0) = ф(х), — то <x< то. (1.2)

В недивергентной форме уравнение (1.1) можно записать как

ut + a(u)ux = 0, (1.3)

где a(u) = df/du.

Если полная вариация начальной функции ф(х) ограничена, то полная вариация решения рассматриваемой задачи Коши не возрастает со временем [8], то есть для любых t2 > ti > 0 имеет место неравенство

TV(u(x,t2)) < TV(u(x,t1)),

где

СО

TV(u(x,tj)) = J

— О

[u]i — скачок функции u(x,tj) на разрыве в точке x^.

Применение адаптивных сеток делает целесообразным переход к подвижной системе координат. Такой переход осуществляется посредством достаточно гладкого невырожденного преобразования координат

du(x, tj) dx

dx + ^2 |[u]i|

(1.4)

t = t, x = x(q, т),

(1.5)

причем при фиксированном т функция х = х(д, т) монотонно возрастает по переменной д. Это преобразование переводит на каждом временном слое Ьп равномерную сетку с длиной ячейки Ад, построенную в области пространственной переменной д, в неравномерную сетку в области пространственной переменной х. Совокупность узлов этих сеток обозначим соответственно {дг} и |х™}.

Далее мы будем обозначать переменную времени в новых координатах по-прежнему через *. Тогда уравнение (1.1) примет вид

(и1)4 + (/(и) - х4и)д = 0, (1.6)

а его недивергентная форма

1 Г 1 (/(и))„ - хщ =0, (1-7)

где 1 = хд — якобиан отображения (1.5). Очевидно, решение задачи Коши (1.6), (1.2) будет удовлетворять условию невозрастания полной вариации (1.4).

Рассмотрим в плоскости переменных (*,ч) прямоугольник Б с вершинами (¿га,дЗ-), (*п, ЧЗ+О, (*п+1, ЧЗ+О и (¿п+1, Ч]). Будем полагать, что искомая функция и задана в полуце-лых узлах пространственной сетки, и введем обозначение и(^т,дг+1/2) = и™^. Запишем разностный аналог закона сохранения в расчетной ячейке Б в виде

(и1 )Я+1/2Ач - (и1 )П+1/2Ач + (Л'+1 - (х*)п+1щ+^ А^а _ (Л' _ (х*)Пщ) А^П = °- (1-8)

Здесь (х*)П = (х™+1 — х™)/А*га, А*п = *п+1 — Ьп (далее для краткости записи индекс п у величин (х4)" и А*п будем опускать), щ, /г — аппроксимации функций и и /(и) в узле сетке с номером г, от выбора которых зависит конкретный вид конечно-разностной схемы. Следуя [7], будем называть (1.8) ТУБ-схемой, если

ТУ^(ига+1) < ТУ^(ига),

где

ГО

ТУ'‘(и”‘)= £ |А<ит|. А,и™ = и”+1/2 — и"-1/2.

2. Построение конечно-разностной схемы

Перейдем к получению конкретных конечно-разностных схем. Вначале построим схему с направленными разностями, имеющую первый порядок аппроксимации. Пусть в (1.8) используются аппроксимации вида

хт х™ 1

1+1/2 = ‘+1Ач ‘ . 4 = 2(и’-,/2 + <+1/2), (2.1)

и

/г = 2 (/(и?-1/2) + /М+1/2) _ Л ЫА,и”) , (2.2)

где Л = А*/Ад, V = Л(« _ (х^)г), А.ип = и™+1/2 _ и™_1/2,

Г (/(и?+1/2) _ /(и?_1/2))/А.ига. если А.иП = °.

«г = \

у а(и™_1/2), если Агип = 0.

Заметим, что использование аппроксимации (2.1) влечет выполнение равенства

п+1 ^п \ | /п+1 п\

_иЗ+1(хП+1 _ хП+1) + иЗ(хП _ хП) =

1

_2 1(хП+1 _ хП+1 )А3+1щП + (хП+1 _ хП)АЗщП) _ иП+1/2АЗ+1/2хП+1 + иП+1/2АЗ+1/2хП. (2-3)

С учетом этого мы можем преобразовать (1.8) к виду

п+1 п 1 п 1

иЯ-1/2 = иП+1/2 _ 2 тп+1 (К' 1 + ^7) АЗ щП + 2 тп+1 (|^З+1| _ ^З+1) А3+1щП

21З+1/2 21З+1/2

или

иП+1/2 = иП+1/2 _ С_,ЗАЗщП + С+,З+1 АJ+1un, (2-4)

где 1

С±,з = 21+г- (|^1 т V?).

21ЗТ1/2

Как показано в [7], условием, достаточным для того, чтобы схема вида (2.4) была ТУБ-схемой, является выполнение неравенств

С_,З > 0, С+, > 0, С_,З + С+З < 1. (2.5)

Введя обозначение

С = тт(1_1/2, 1г+1/2),

получим, что условия (2.5) выполняются, если

IV,| < ¿п+‘- (2.6)

Таким образом, при выполнении неравенства (2.6) схема (1.8), (2.1), (2.2) является ТУБ-

схемой.

Рассмотрим построение для уравнения (1.6) разностной схемы второго порядка типа одношаговой схемы Лакса — Вендроффа (см., например, [9]). Разложим функцию и1 в ряд Тейлора по * на п + 1-м слое относительно п-го слоя:

А*2

(и1)ПН-1/2 = (и1 )п+1/2 + А* ((и1 ).)П+1/2 + -у- ((и1 )«)П+1/2 + 0(Д*3)- (2-7)

Используя дифференциальные следствия уравнений (1.6) и (1.7), получим, что

(и1 )* = _ (/(и) _ х*и)0 = _ (/(и))0 + (х*и)д

и

(и1 )Й = _ ((/(и) _ х*и)4)д = _ ((/(и))и и _ х*и* _ х**и)д

т

С учетом равенства

= (а(и) _ х*) - (/ _ х*щ) ) + (х**и)д.

А* \п хп+1 _ хп

х* + -2-х*^ . = г А* г + 0(А*2)

ч

имеем

Лі2 Г „Х-1/2 + <

^Ліх*и + ^и) = «+1 - хП)Иі-1/2 + Иі+1/2 + 0(Лі3 + ЛіЛд2).

2

Тогда

Лі2

Лі (М)4)”+1/2 + — (М)**);+1/2 =

Д./(иі+3/2) /(иі-1/2) + Лі Л ) иі+3/2 + иі+1/2 ( ) ИІ+1/2 + Иі-1/2 . +

- Лі--------2Л----------+ Л І0^1---------2-------(х^'-----------------2-1 +

Лі2 + 2А^

аі+1 - (х*)і+1 / /(иП+3/2) - /(иП+1/2) - ( ) и?+3/2 - иП+1/2

Д*+1 V Лд (Х*)^'+1 Лд

- аі - (х*)і / /(Ип+1/2) - /(иП-1/2) - (х ) <+1/2 - <-1/2 ^ Лд ( Лд

1Ш І 1Ш

<Л_1 !О +

+ 0(Лі3 + Лд3)

(здесь и далее полагается Д™ = -—-^—^—¿+1/2, Лі = О(Лд)). Отбрасывая члены порядка

0(Лі3 + Лд3), из разложения (2.7) получим одношаговую схему Лакса — Вендроффа для подвижных адаптивных сеток в форме (1.8), где Д+1/2 и иг аппроксимируются согласно (2.1), а

/г = 1 /(<-1/2) + /(<+1/2) - Д М^иЛ . (2.8)

Очевидно, данная схема не удовлетворяет условию ТУБ.

Проведем монотонизацию полученной схемы. Как и выше, запишем иг согласно (2.1), а / аппроксимируем следующим образом:

■г = 2 (у(иі—1/2) + /(иП+1/2) + д (^г-1/2 + #г+1/2 - N + 7г|ЛгиП)^ • (2-9)

Здесь

2 і */ V г— 1/2/ 1 о \~г+1/2/ 1 д

«г+1 тіп(|^і|, |#г+1І), если #¿#¿+1 > 0,

#г+1/2 = \

0, если #г#г+1 < 0,

1 Л | 1 / ^ д п ч #г+1/2 - 1/2

#г = 2 - ДП Ы ^ Лги , «г = Sgn(#г), ^ = -------Л"ЦП------•

Докажем второй порядок аппроксимации полученной схемы. Для этого достаточно установить (см. [7]), что разность правых частей выражений (2.8) и (2.9) есть величина порядка 0(Лд2), т.е.

#¿-1/2 + #¿+1/2 - + 7г1 - дп(^г)2^ ЛгиП = 0(Лд2).

Последнее равенство можно переписать как

#г-1/2 + #г+1/2 + (|^г1 - |^г + 7г|) ЛгиП = ^N1 - дп(^г)2^ ЛгиП + 0(Лд2). (2^10)

Используя методику работы [7], из определения #¿+1/2 легко получить, что #¿-1/2 = #г + О^?2^ #г+1/2 = #г + °(Л?2).

Отсюда сразу следует (2.10), поскольку

|(Ы - IV + 7г|)ЛгМп| < І7гЛгип| = |#г+1/2 - #г-1/2І = 0(Лд2).

Итак, мы показали, что схема (1.8), (2.1), (2.9) имеет второй порядок аппроксимации. Кроме того, так как числа #г-1/2 и #г+1/2 не могут иметь разные знаки, то

|#г+1/2 - #г—1/21 < тах(|#г-1/21, |#г+1/2|) < |#г1

и

|7г| < |#г|/|ЛгиП| < ^

N - ДП(^г)

2

1----------IV |

Т п'

Используя тождество (2.3), запишем рассматриваемую схему в виде

(2.11)

ип+1 = ип____________________

и+1/2 ИІ+1/2 2 тп+1

2Т,+1/2

1 V, + 7, | + V, + №+1/2.~ #,—11/2 I Л,- и”+

Л, ип

п+1

2Т,+1/2

(і«,+!+7,+.! - *,+, - л,+1Ип,

то есть

и

"+1/2 = иП+1/2 - 2 тп+1 С + ^М) Л,иП + 2 тп+1 (|-“+1| - ,+1) Л,+1иП

2Т,+1/2 2Т,+1/2

где Vм = V + 7г. Таким образом, рассматриваемая схема приведена к форме (2.4), причем в определении коэффициентов С±,, величины V, заменяются величинами V“. Ввиду этого схема (1.8), (2.1), (2.9) является ТУБ-схемой, если

|, < ^п+1.

Используя неравенство (2.11), получим оценку

(2.12)

I, < IV,I + І7,I < IV,И 1 + ^

1 - ^ IV, |

Т П I */1

Пусть

IV,| < ¿П

(на неподвижной сетке это условие эквивалентно условию (2.6) для базовой схемы (1.8), (2.1), (2.2)). Тогда имеет место неравенство

3,

то есть (2.12) выполняется, если

| V, I < -шіп(^п,^п+1)

3

(2.13)

Неравенство (2.13) является достаточным условием того, чтобы схема (1.8), (2.1), (2.9) была ТУБ-схемой.

1

Замечание 1. Условие (2.13) можно переписать в виде

^Гн75!“•> - (x<)jI £ 2• (2.14)

(Ах)?

(Ах)П+1/2 _ ш:л(хП ____ хП хп _______ хп хп+1 ____ хп+1 хп+1 __ хп+1)

где (Ах)? — 111111(Х? х?___1, х?+1 X? , х? х?_1, х?+1 X? ).

Замечание 2. Для того, чтобы построенная схема была е-схемой (т. е. удовлетворяла энтропийному неравенству), возможна ее регуляризация заменой в формуле (2.9) |^ + 7^ на ф(^ + 7г); а в определении заменой |^| на ф(^), где

z2 + є2

если |zI < є,

_ { 2е

|г|, если |г| > е,

е — малый положительный параметр (см. [7]). При этом условие (2.13) остается в силе.

3. Тестирование алгоритма

Опишем численные эксперименты, демонстрирующие преимущества построенной схемы. В этих экспериментах для построения подвижных адаптивных сеток нами использовался модифицированный метод эквираспределения (см., например, [1, 6]). Сетка строилась путем решения конечно-разностного аналога уравнения

d хП+1 _ хП

—(w(xn,ira)xn+1) = e (—---------un )• (3-1)

Здесь

w =1 + a0|u| + a1 |ux| — (3.2)

управляющая функция, a0, a1, в — неотрицательные параметры.

Несколько слов о практическом выборе величины At. Проблема состоит в том, что в оценку (2.14) входят значения xn+1 _ Х—1 и я^1 _ xn+1, определяемые после выбора At. Поэтому, если первоначальный выбор At не привел к выполнению условия (2.14), то указанную величину следует уменьшить (например, в два раза), построить новую сетку К+1}, вновь проверить выполнение условия (2.14) и т.д. При достаточно малых значениях At условие (2.14) будет выполнено, так как в расчетной области, являющейся компактом в пространстве (q,t), функции a(u) и xt ограничены сверху, а значения (Ax)n+1/2 при заданном числе узлов сетки N отделены от нуля снизу в силу невырожденности преобразования координат.

В описываемых ниже расчетах мы использовали более простой прием выбора At. Поскольку применяемый нами метод построения сеток не приводит к резкому изменению сетки на очередном слое по времени, в формуле (2.14) мы полагали (Ax)™+1/2 = = min(xn _ xn_1, xn+1 _ xn), а функция (xt)j аппроксимировалась разностью назад, при этом At бралось с коэффициентом запаса, примерно равным 0.5. На практике это обеспечивало выполнение условия (2.14).

В качестве тестовых уравнений рассматривались уравнение переноса, для которого f (u) = au, a = const, и уравнение Хопфа, иногда называемое уравнением Бюргерса. В этом уравнении f (u) = u2/2.

u(x, 0)

Для расчетов по уравнению переноса с параметром a = 1/2 были взята начальная функция вида

0, если 0 < x < xi,

A1 sin2(w1(x — x1)), если x1 < x < x2,

A1 cos(w2(x — x2)), если x2 < x < x3,

A2(x — x3)/(x4 — x3), если x3 < x < x4,

A2(x5 — x)/(x5 — x4), если x4 < x < x5,

0, если x5 < x < 21,

где w2 = A2/(A1(x4 — x3)), x2 = x3 — 3n/(2w2), x1 = x2 — n/(2w1). Выполнение этих условий обеспечивает существование производной начальной функции в точках x1, x2, x3. В расчетах полагалось A1 = 3/4, A2 = 3/2, w1 = 1/2, x3 = 9, x4 = 11, x5 = 13.

На рис. 1, а приведены графики точного и численного решения рассматриваемой задачи в момент времени t = 11. Видно, что решение, полученное по схеме (1.8), (2.1),

(2.9), не осциллирует. На рис. 1, б изображены траектории узлов использованной в расчете подвижной адаптивной сетки с числом узлов N = 33, построенной при а0 = 2.1, а = 0.15, в = 25.

Для уравнения Хопфа проводилось сравнение свойств построенной схемы со свойствами явной схемой предиктор — корректор с автоматически настраиваемой аппроксимацион-ной вязкостью [1, 6]. Были рассмотрены гладкое и разрывное решения уравнения Хопфа, взятые из [3]. Начальные данные для гладкого решения следующие:

u(x, 0) = a1 + (a2 — a1) exp[—(32/x)2], 0 < x < 256.

Точное решение этой задачи непрерывно и имеет вид

u(x, t) = u(z, 0),

где z есть решение нелинейного уравнения

x = z + u(z, 0)t = z + (a1 + (a2 — a1) exp[—(32/z)2])t.

В расчетах полагалось a1 = 0.5, a2 = 1.5.

На рис. 2 приведены графики точного и численных решений рассматриваемой задачи в момент времени t = 120. Численное решение на рис. 2, а получено по схеме (1.8), (2.1), (2.9); на рис. 2, б — по схеме с автоматически настраиваемой аппроксимационной вязкостью.

На рис. 2, в изображены траектории узлов использованной в расчете по схеме (1.8), (2.1),

(2.9) подвижной адаптивной сетки с числом узлов N = 33, построенной при а0 = 0, а1 = 0.2, в = 5.

Для разрывного решения начальные данные брались в виде

1, если 0 < x < 2,

u(x, 0) = < (8 — x)/6, если 2 < x < 8,

У 0, если 8 < x < 16.

Точное решение при 0 < t < 6 имеет вид

1, если 0 < x < 2 + t, u(x,t)=^ (8 — x)/(6 — t), если 2 + t < x < 8,

0, если 8 < x < 16.

Рис. 1. Результаты расчетов для уравнения переноса.

В момент времени Ь = 6 происходит градиентная катастрофа. Если Ь > 6, то

1, если 0 < х< 5 + Ь/2,

и(х, і) =

0, если 5 + і/2 <х < 16,

т. е. при і > 6 решение разрывно.

На рис. 3 приведены графики точного и численных решений рассматриваемой задачи в момент времени і = 12. Численное решение на рис. 3, а получено по схеме (1.8), (2.1), (2.9); на рис. 3, б — по схеме с автоматически настраиваемой аппроксимационной вязкостью. На рис. 3, в изображены траектории узлов использованной в расчете по схеме (1.8), (2.1),

(2.9) подвижной адаптивной сетки с числом узлов N = 33, построенной при а0 = 0, аі = 0.1, в = 5.

а

4. Заключение

Таким образом, нами построена консервативная схема второго порядка аппроксимации для расчетов на подвижных адаптивных сетках. При выполнении неравенства (2.13) схема удовлетворяет условию ТУБ. Численные эксперименты показали, что построенная схема обеспечивает отсутствие осцилляций решения на разрыве и достаточно точно передает движущийся фронт волны.

Авторы выражают глубокую признательность Г. С. Хакимзянову за ценные советы и

замечания, позволившие улучшить статью.

Список литературы

[1] Барахнин В. Б. Конечно-разностные схемы для численного решения задач теории мелкой воды с использованием адаптивных сеток. В “Вычислительные технологии”. ИВТ СО РАН, Новосибирск, 4, №11, 1995, 38-50.

[2] ГильмАнов А. Н., КулАчковА Н. А. Численное исследование двумерных течений газа со скачками методом TVD на физически адаптивных сетках. Матем. моделирование, 7, №3, 1995, 97-107.

[3] Гужев Д. С., КАлиткин Н. Н. Уравнение Бюргерса — тест для численных методов. Там же, 7, №4, 1995, 99-127.

[4] Карамышев В. Б. Монотонные схемы и их приложения в газовой динамике. Новосибирский университет, Новосибирск, 1994.

[5] Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений и их приложения к газовой динамике. Наука, М., 1978.

[6] BARAKHNIN V. B., KHAKIMzYANov G. S. On the application of adaptive grids to the numerical solution of one-dimensional problems in the shallow-water theory. Russ. J. of Numer. Anal. and Math. Modelling. 10, No. 5, 1995, 373-391.

[7] Harten A. High Resolution Schemes for Hyperbolic Conservation Laws. J. of Comput. Phys. 49, No. 3, 1983, 357-393.

[8] Lax P. D. Hyperbolic system of conservation laws and the mathematical theory of shock waves. SIAM Regional Series on Appl. Math., No. 11, 1973.

[9] Lax P. D., Wendroff B. Difference schemes for hyperbolic equations with high order of accuracy. Commun. Pure Appl. Math. 17, No. 3, 1964, 381-398.

Поступила в редакцию 23 декабря 1998 г., в переработанном виде 14 января 1999 г.