CC BY
5
Вычислительные технологии
Том 23, № 2, 2018
О монотонности и точности схемы КАБАРЕ при расчете обобщенных решений с ударными волнами
О. А. КовыркинА1'*, В. В. ОСТАПЕНКО1,2
1Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия 2 Новосибирский государственный университет, Россия *Контактный e-mail: olyana@ngs.ru
Изучены монотонность и точность модифицированной схемы КАБАРЕ, аппроксимирующей квазилинейную гиперболическую систему законов сохранения. Получены условия, при которых эта схема сохраняет монотонность разностного решения относительно инвариантов линейного приближения аппроксимируемой системы. В качестве конкретного примера рассмотрена аппроксимация системы законов сохранения теории мелкой воды. На примере этой системы показано, что подобно TVD-схемам повышенного порядка аппроксимации на гладких решениях схема КАБАРЕ, несмотря на высокую точность при локализации ударных волн, снижает свой порядок сходимости в областях их влияния.
Ключевые слова: гиперболическая система законов сохранения, монотонность и точность схемы КАБАРЕ, теория мелкой воды, прерывные волны.
Введение
Для численного решения гиперболических уравнений предложена трехслойная по времени и двухточечная по пространству схема Upwind Leapfrog [1], которая имеет второй порядок аппроксимации на гладких решениях, является явной и условно устойчивой при числах Куранта г Е (0,1]. Детальный анализ этой схемы проведен в работах [2, 3], где с учетом кососимметричности своего пространственного шаблона она названа схемой КАБАРЕ. Основные достоинства этой схемы связаны с тем, что она задана на компактном пространственном шаблоне, является обратимой по времени и точной при двух различных числах Куранта г = 0.5,1, что наделяет ее уникальными диссипатив-ными и дисперсионными свойствами [3].
Для численного решения уравнений одномерной газовой динамики [4, 5] разработан балансно-характеристический вариант схемы КАБАРЕ [6], который с учетом коррекции потоковых переменных (необходимой для монотонизации разностного решения на ударных волнах) показал высокую точность при расчете классического теста Blast Wave [7]. В настоящее время для численного моделирования пространственно многомерных газодинамических [8-10] и гидравлических [11] течений широко применяется двухслойная по времени форма записи схемы КАБАРЕ [12]. Монотонность этой схемы при аппроксимации линейного уравнения переноса в одномерном случае изучалась в [13], а в многомерном — в [14]. Условия монотонности двухслойной по времени схемы КАБАРЕ, аппроксимирующей скалярный закон сохранения с выпуклым потоком, исследовались в [15].
© ИВТ СО РАН, 2018
В настоящей работе результаты, полученные в [13-15], использованы для построения схемы КАБАРЕ, аппроксимирующей квазилинейную гиперболическую систему законов сохранения и обеспечивающую сохранение монотонности разностного решения относительно инвариантов линейного приближения аппроксимируемой системы. В качестве конкретного примера рассмотрена монотонная аппроксимация этой схемой системы законов сохранения первого приближения теории мелкой воды. На примере этой системы показано, что подобно ТУВ-схемам повышенного порядка аппроксимации на гладких решениях [16- 18] схема КАБАРЕ, несмотря на высокую точность при локализации прерывных волн, снижает свой порядок сходимости в областях их влияния. Приведены тестовые расчеты, иллюстрирующие эти свойства построенной схемы.
1. Схема КАБАРЕ, аппроксимирующая
гиперболическую систему законов сохранения
Рассмотрим задачу Коши для квазилинейной гиперболической системы законов сохранения
и + £(и)л = 0, и(ж, 0)= (1.1)
где и(ж, ¿) — искомая, а £(и) и ч(х) — заданные гладкие вектор-функции, содержащие т компонент. Недивергентная форма записи системы (1.1) имеет вид
и4 + А(и)иж = 0, А(и) = £ и (и). (1.2)
Гиперболичность системы (1.2) означает, что все собственные значения АДи) матрицы А(и) действительны и соответствующие им левые собственные векторы Р(и) образуют базис в пространстве Кт. Будем предполагать, что все функции АДи) являются знакопостоянными, т. е. принимают только положительные или только отрицательные значения. Будем также считать, что каждая дифференциальная форма Р(и)^и интегрируема, в силу чего система (1.2) допускает запись в форме инвариантов
(и)
гМ + = 0' < =1'™' (Ь3)
где вектор инвариантов w = (/ш\,/ш2, ...,'шт) системы (1.3) связан с вектором базисных переменных и невырожденным преобразованием w = Ш(и).
Применяя результаты, полученные в [13, 15], построим модификацию двухслойной по времени схемы КАБАРЕ [12], которая при аппроксимации задачи Коши для гиперболической системы законов сохранения (1.1) обеспечивает сохранение монотонности разностного решения относительно инвариантов линейного приближения этой системы. Схему КАБАРЕ будем строить на прямоугольной разностной сетке
{Ху, Ьп] : Ху = ]А, ¿„+1 = Ьп + тп, ¿о = 0, (1.4)
в которой А — постоянный шаг сетки по пространству, а тп — переменный шаг сетки по времени. В этой схеме используются потоковые и™ = u(xj, Ьп) и консервативные Ц/+1/2 = и(ж^+1/2, Ьп) векторные переменные, заданные соответственно в целых х^ и полуцелых х¿+1/2 = х^ + А/2 пространственных узлах разностной сетки.
Пусть и™, и+1/2 — известное численное решение задачи (1.1) на п-м временном слое ¿п, при п = 0 — сеточная аппроксимация начальной функции у(ж). Применяя
подход, предложенный в [15], численное решение и™+1, и™_+11/2 на (п + 1)-м временном слое ¿га+1 будем находить по схеме КАБАРЕ в четыре этапа. На первом этапе по разностным уравнениям
и«.+ 1/2 ига сга сга
/9 Т1 + 1/2 — 1о
'+1/2 '+1/2 + 1 ^ =0, Г = I(ига), (1.5)
тп/2 А ' ' 4 ^
вычисляются значения консервативных переменных и™++11//^ = и(ж^+1/2, ¿га+1/2).
тт ^га+1/2 п( га+1/2ч га+1/2
На втором этапе находятся значения потоков I^ = I (и^- ), где и = , ¿га+1/2), которые применяются на третьем этапе для определения из разностных уравнений
ттга+1 _ ига ^га+1/2 /.га+1/2
ц'+1/2 ц'+1/2 + 1 " -1 -= 0 (1.6)
Тга А
значений консервативных переменных ига+11/2
На четвертом этапе находятся значения потоковых переменных и*^1.
2. Вычисление потоков на втором этапе
и потоковых переменных на четвертом этапе
В начале второго этапа по формулам
wгa=ж (ига), wгa+l/2=ж (ц^), wгa++í//22=^ (ига++!/2) (2.1)
определяются численные значения вектора инвариантов, которые используются для вычисления компонент этого вектора в целых узлах на временных слоях (п + 1/2) и (п+1). Опишем процедуру этого вычисления для некоторого инварианта 'ш^, переносимого со скоростью \г вдоль характеристик г-го семейства. Для краткости в приводимых ниже формулах индекс г будем опускать.
Предварительное значение инварианта ^га+1 определяется в зависимости от знака функции А(и). При выполнении неравенства
А(и) > 0 (2.2)
предварительное значение инварианта +1 вычисляется по формуле
г^1 = 2<12 - (2.3)
а при выполнении противоположного неравенства
Л(и) < 0 (2.4)
— по формуле
г^1 = 2<++1/2 - ^га+1. (2.5)
После этого находится предварительное значение инварианта
г«га+1/2 = (< + г«га+1) /2. (2.6)
При помощи стандартного двухстороннего ограничителя
u < т,
F (u, т, М)=< u, т <u < М,
М, u > М,
значения w™+1/2 с учетом результатов работы [13] корректируются по формуле
w™+1/2 = F (<+1/2, т?, М, (2.7)
в которой
т™ = min {w]-1/2, max(w™, ^]-1/2)) , М™ = max (w"J-1/2, min(w™, ^™-1/2)) при выполнении неравенства (2.2) и
т™ = min (Ц+1/2, max(w], ^+1/2)) , М™ = max (w"+V2, min(w™, ^+1/2)) при выполнении неравенства (2.4), где
га = wj—1/2- (i- ^;_+11//22) w-1 n = w"+1/2- (i- ^р w"+1
^i-1/2 = n+1/2 , ^i+1/2 = , n+1/2.
га+1/2 ' ' и IV" | га+1/2|
^-1/2 ^¿+1/2 1
га+1/2 = Тп | л га+1/2 | га+1/2 = тп , л га+1/2 |
^-1/2 = А |Л^-1/2 |, ^+1/2 = А |Л?+1/2 |.
После вычисления вектора инвариантов w™+1/2 во всех узлах ] полуцелого временного слоя по формулам
fn+1/2 = f /™+1/2\ , un+1/2 = w-1 , (2.8)
где Ш-1 - оператор, обратный к оператору Ш, находятся численные значения потоков £»г+1/2, которые на третьем этапе используются в разностном уравнении (1.6). На четвертом этапе по формуле
W™+1 = 2w™+1/2 - w? (2.9)
находится второе предварительное значение вектора инвариантов. Компоненты этого вектора корректируются по формуле
w™+1 = F (w^1, т™+1, М7+1) , (2.10)
в которой
rn? = min (w^, w^) , М™ = max (w^, w^) , где w;+/2 и w™+11/2 - соответствующие компоненты векторов
W™_+i/2 = W (U-i/2) , wJ+11/2 = W (U^) .
После определения вектора инвариантов w™+1 по формуле
u™+1 = W-1 (w^1) (2.11)
находится вектор потоковых переменных на (п + 1)-м временном слое. На этом завершается вычисление разностного решения ига^1, ига+11/2 на (п + 1)-м временном слое. Шаг разностной схемы по времени га выбирается из условия устойчивости
гА
(2.12)
га ^ |( Л^)га+1/2| '
где г е (0,1), ( Лг)га+1/2 = Мига+1/2).
Если в результате применения формул коррекции (2.7) и (2.10) предварительные
га+1/2 —га+1/2 га+1 ~га+1
значения инвариантов не изменяются, т. е. = и^ и = к^- , то построен-
ная модификация схемы КАБАРЕ имеет второй порядок классической аппроксимации в смысле тейлоровского разложения на гладких решениях. Численные расчеты показывают, что формулы (2.7) и (2.10) приводят к коррекции инвариантов в нескольких узлах, расположенных в окрестности сильного разрыва, где отсутствует локальная (поточечная) сходимость разностного решения к точному. Если в результате такой коррек-
га+1/2 / -га+1/2 га+1 / —га+1
ции и^ = и^ или = ць , то порядок классической аппроксимации схемы
КАБАРЕ снижается до первого. Как показано в следующем разделе, коррекция инвариантов (2.7) и (2.10) обеспечивает монотонность разностного решения относительно инвариантов линейного приближения аппроксимируемой гиперболической системы (1.1).
3. Монотонность схемы КАБАРЕ
Линейным приближением системы (1.2) относительно постоянного решения u0 = const называется линейная гиперболическая система
5ui . . <9ui
-U + Ж*,) ^ = 0, (3-l)
получаемая в первом приближении по малому параметру е ^ 1 при подстановке в систему (1.2) решения u(;r, t) = u0 + eu1(;r, t). Линейная система (3.1) имеет полный набор инвариантов w^(u1) = lt(u0)u1, что позволяет записать ее в виде m линейных уравнений переноса
—+ Aj(uo) — =0, 1=^. (3.2)
Если схема КАБАРЕ (1.4)-(1.6), (2.1)-(2.11) аппроксимирует линейную гиперболическую систему (3.1), то входящий в нее численный поток fra является линейной вектор-функцией fra = A(u0)(u1)™. Умножая с учетом этого векторные разностные уравнения (1.5) и (1.6) слева на левый собственный вектор 1, соответствующий собственному значению А матрицы A(u0), получаем линейные разностные уравнения
„ га+1/2 _ „ п „га _ „га
Wj+1/2 Wj+1/2 + a„J+1 W3 = 0
Гга/2 А ,
„,,"+1 га+1/2 га+1/2
^■+1/2 - ^'+1/2 + ' = 0,
тп А
аппроксимирующие соответствующее уравнение (3.2) и описывающие изменение численных значений инварианта = 1и1 на первом и третьем этапах реализации схемы КАБАРЕ. Отсюда с учетом формул (2.1)-(2.11) и результатов работы [13] следует сохранение монотонности численных значений инварианта при переходе с п-го на (п + 1)-й временной слой. Это означает, что построенная схема КАБАРЕ сохраняет монотонность разностного решения относительно инвариантов линейного приближения аппроксимируемой системы.
4. О точности разностных схем сквозного счета
В классической работе [19], широко известной в связи со схемой распада разрыва, введено понятие монотонности разностной схемы и показано, что среди линейных разностных схем нет монотонных схем повышенного порядка аппроксимации. Дальнейшее развитие теории разностных схем сквозного счета для гиперболических систем законов сохранения в значительной степени было направлено на преодоление этого "запрета Годунова". В результате разработаны различные классы разностных схем, в которых повышенный порядок аппроксимации на гладких решениях и монотонность (при аппроксимации линейной системы и скалярного закона сохранения) достигались за счет нелинейной коррекции потоков, приводящей к нелинейности этих схем даже при аппроксимации линейного уравнения переноса. Перечислим основные классы таких схем: МиБСЬ-схемы [20], ТУБ-схемы [21], КЕБ-схемы [22], ШЕШ-схемы [23], САБАЯЕТ-схемы [12]. Основное достоинство этих схем заключается в том, что они с высокой точностью локализуют ударные волны при отсутствии существенных нефизических осцилляций.
При построении указанных выше схем повышенный порядок аппроксимации понимается в смысле тейлоровского разложения на гладких решениях, что не гарантирует аналогичного повышения точности при расчете разрывных решений. В работах [16- 18] было показано, что различные типы ТУБ-схем имеют не более чем первый порядок локальной (поточечной) сходимости в областях влияния ударных волн и тем самым, по существу, схемами повышенной точности не являются. Такое снижение порядков сходимости свидетельствует о том, что в этих схемах происходит потеря точности при передаче условий Гюгонио через размазанные фронты ударных волн. В работах [24, 25] для ТУБ-схемы формально второго порядка, а в [26] для трех разных типов ШЕКО-схем формально пятого порядка точность передачи схемой условий Гюгонио через фронт ударной волны оценивается путем определения порядка сходимости интеграла от разностного решения (а не от его модуля, как в норме Ь1), что соответствует сходимости в соответствующей негативной норме [27]. В [24-26] показано, что в ТУБ- и ШЕКО-схемах такой порядок интегральной сходимости снижается до первого на интервалах интегрирования, одна из границ которых находится в области влияния ударной волны.
Одна из основных причин такого снижения точности заключается в том, что минимаксная коррекция потоков, характерная для этих схем, вызывает снижение гладкости разностных потоков, что в свою очередь приводит к снижению порядка аппроксимации е-условий Гюгонио на фронтах ударных волн [28]. В то же время, как показано в [24, 25], классические немонотонные схемы повышенной точности, имеющие аналитические функции численных потоков и, как следствие, с повышенной точностью
аппроксимирующие е-условия Гюгонио, сохраняют повышенный порядок сходимости в негативной норме при интегрировании по областям, содержащим сильные разрывы. В результате эти немонотонные схемы, в отличие от схем типа ТУВ, сохраняют повышенный порядок сходимости в областях влияния ударных волн, несмотря на заметные схемные осцилляции на их фронтах.
Далее в данной работе указанными выше методами исследуется точность модифицированной схемы КАБАРЕ.
5. Методы оценки точности схемы КАБАРЕ
Предположим, что задача Коши (1.1) имеет единственное обобщенное решение и (ж, ¿), ограниченное при £ > 0. Для приближенного определения порядка сходимости к этому точному решению разностного решения, получаемого по схеме КАБАРЕ на последовательности сжимающихся сеток, будем использовать равномерные разностные сетки, в которых шаг по времени является постоянным и с учетом условия устойчивости (2.12) определяется по формуле
гД . .
т = т 1П тп =-——-г. (5.1)
п тах |( Л, )П+1/2|
Поскольку в схеме КАБАРЕ основными переменными, определяющими точность разностного решения, являются консервативные переменные Ип+1/2, заданные в полуцелых пространственных узлах ж^-+1/2 = ж^- + Д/2 = (] + 1/2)Д разностной сетки, то для корректного применения метода Рунге необходимо, чтобы полуцелые узлы основной базисной сетки совпадали с некоторыми полуцелыми узлами более мелких сеток. Это означает, что пространственный шаг более мелкой сетки должен быть как минимум в 3 раза меньше, чем пространственный шаг более крупной сетки. Отметим, что при оценке точности ТУВ- и ШЕКО-схем [24 - 26], в которых основные переменные заданы в целых узлах разностной сетки ж^- = ^'Д, минимальное отношение пространственных шагов более крупной и более мелкой сеток равно 2.
Зафиксируем пространственный шаг Д равномерной базисной сетки (1.4), (5.1) и рассмотрим последовательность разностных решений (Ик)п+1/2 = и(жк+1/2, ¿П) задачи Коши (1.1), получаемых по схеме КАБАРЕ на последовательности сжимающихся сеток
Д
{ж£, О : ж£ = зДк, ¿П = птк, п > 0; Дк = -т, тк = -т, к > 0.
Выберем момент времени Т = тт = 3ктТк > 0, для которого путем линейной интерполяции доопределим эти разностные решения до непрерывных по ж функций
к
ик(ж, Т) = (ик)3-Т/2 + ж-Д71/2 ((ик)3+?/2 - (ик,
где ж € [жк_1/2, жк+1/2]. Зафиксируем интервал [а, Ь] и зададим интегралы
ь ъ
и(Т,а, Ь)= и(ж,Т)^ж, Ик(Т, а, Ь)= Ик(ж, Т)¿ж.
(5.2)
Дадим следующее определение. На интервале [а, Ь] последовательность разностных решений {(и к)™+1/2} с порядком р < 2 сходится в интегральной норме к точному решению и (ж, ¿), если с точностью до о(А^) выполнено условие
5Ик(Т, а, Ь) = И к(Т, а, Ь) - и(Т, а, Ь) = СА£, (5.3)
где вектор-функция С не зависит от Ак.
Ограничение р < 2 связано с тем, что с учетом линейной интерполяции (5.2) интегралы И к(Т, а, Ь) вычисляются по формуле трапеций, которая имеет второй порядок точности на гладких функциях. Если рассчитываемое точное решение и (ж, ¿) не имеет особенностей, то порядок интегральной сходимости для разностной схемы повышенной точности на гладких решениях будет равен двум. В ином случае порядок такой сходимости на интервале [ , ], содержащем особенности точного решения или пересекающемся с их областями влияния, может снижаться ниже двух, что зависит от точности передачи схемой условий Гюгонио через размазанные фронты ударных волн [24, 25].
Для приближенного определения порядка интегральной сходимости в случае, когда точное разрывное решение и (ж, ¿) заранее неизвестно, необходимо провести три расчета по схеме КАБАРЕ с достаточно малыми шагами А0 = А, А1 = А/3, А2 = А/9 и воспользоваться правилом Рунге. Вычитая из формулы (5.3) эту же формулу, в которой индекс к заменен на к + 1, получаем
Шк,к+1(Т, а, Ь) = Ик(Т, а, Ь) - И+1(Т, а, Ь) = С(А? - Ак+1). (5.4)
Взяв отношение модулей равенств (5.4) при к = 1 и к = 0, находим
|Ш1.2(Т,а,6)| = А? - А? = (А/3)" - (А/9)" _ (1' |<Шод(Т, а, &)| А? - А? А? - (А/3)?
Отсюда получаем искомую формулу
(3)'
(Г м 1 1 ^UlЛT,a, b)l p = p{T,a, 6) = log1/^———-—. (5.5)
1 |oUo;i(T, a, b)l
Если предположить, что последовательность разностных решений U k (х, T) с порядком а сходится в некоторых точках х к точному решению u(x,T), т.е. в этих точках с точностью до о(Д£) выполнено условие
SUk(х, T) = Uk(х, T) - u(x, T) = еД£, (5.6)
где вектор-функция c не зависит от Дk, то порядки такой локальной (поточечной) сходимости а = a(x,T), определяемые по формуле, аналогичной (5.5), сильно осциллируют в областях влияния ударных волн, что не позволяет эффективно оценить их реальные значения. Как показано в [16, 29], подобная ситуация имеет место также и для других разностных схем повышенной точности, в том числе TVD-схем [21]. Это означает, что в разностных схемах сквозного счета повышенной точности классическая локальная сходимость (5.6) отсутствует не только в окрестностях ударных волн, но также и в значительной части их областей влияния. Поэтому, следуя [29], локальную точность схемы КАБАРЕ при расчете разрывных решений будем определять путем
вычисления порядков сходимости осредненного разностного решения Ик(ж,Т, е), получаемого по рекуррентной формуле
X+£
Ик(ж,Т,= ^ У Ик-1(£,Т,И°°(ж,Т,е) = Ик(ж,Т), г = М, (5.7)
где ¡и£ — заданные параметры осреднения.
Предположим, что в некоторых точках ж осредненное разностное решение Ик (ж, Т, е) с порядком Р сходится к осредненному точному решению иг(ж,Т, е), определяемому по рекуррентной формуле, аналогичной (5.7):
Ж+£
и (ж, Т, = и,-1(С,Т, и°(ж, Т, е) = и(ж,Т), ¿=1,7.
Это означает, что в таких точках ж с точностью до о(Д^) выполнено условие
Шгк (ж, Т, е) = Ик (ж, Т, е) - иг(ж,Т, е) = еД£, (5.8)
где вектор-функция е не зависит от Дк. Для приближенного определения порядка локальной сходимости Р, по аналогии с (5.5), получаем формулу
| ¿И12 (ж, Т, е)| , ч
р = р(х,T,е) = к^/з , <5-9)
в которой
Шк,к+1 (ж, Т, е) = Ик (ж, Т, е) - Ик+ (ж,Т, е), к = 0,1. (5.10) Из формулы (5.10) при к = 0 с учетом формулы (5.8) при к = 0,1 имеем
Ш°д(ж,Т, в) = е (Д? - Д?) = еД^ (1 - (1/3)") .
Отсюда с учетом (5.9) находим
1
_ Ш°,1 (ж,Т,е) Ш°,1 (ж,Т,в)/ |Ш1,2(ж " " е = . „ '-,„ =-^- 1 - ,
/ |Ш1,2(ж,Т,е)|\ _ ^ | Ш°д(ж,Т, е)^ .
Д/» (1 - (1/3)/») Д3 ^ | ш°,1(
Подставляя это значение вектор-функции е в формулу (5.8), в которой к = 0, получаем приближенное выражение
(|^Иг (ж Т е)| \ -1
1 -Щ^Т^) (5.11)
для вектора локальных дисбалансов осредненного разностного решения на базисной сетке с пространственным шагом Д = Д°.
В ряде случаев, наряду с оценкой точности вычисления базисных переменных (5.11), полезно иметь информацию о точности, с которой находятся инварианты решения.
По аналогии с формулой (5.11) для определения дисбалансов осредненного численного значения вектора инвариантов w = Ш(и) на базисной сетке с пространственным шагом А = А0 имеем приближенную формулу
1
| (К, (Х,Т, £)1
6wl(х, Т, е) = 6wl0(x,Т, е) = 6wl0Л(x,Т, е) I 1 - —
/ I ^1,2 (Х,Т, е)| \ ^ | ^0,1 (Х,Т, е)1)
I ^0,1 (х,т, е)1
^к,к+1(x, Г1, £) = ^(x, Г1, £) - ^+1(x,Т, ¿0, к = 0, -1,
где осредненные значения wk(х,Т, е) вычисляются по рекуррентной формуле, аналогичной (5.7),
wk(x,Т, £) = 1 I wk-1(£,Т, е)(%, w0k (Х,Т, £) = Wк (х,т), г = ТЛ,
Ж —£
в которой w к(х,Т) = Ш(Ик(х,Т)).
6. Система уравнений теории мелкой воды
В качестве конкретной гиперболической системы (1.1) выберем систему уравнений первого приближения теории мелкой воды [30], для которой
и = ( * ) , <(и) = ( ф/к + 9Л2/2 ) . <6-Ч
где Ъ,(Ь,х) и ,х) — глубина и расход жидкости; д — ускорение свободного падения. При записи этой системы в форме инвариантов (1.3) получаем
А1 = V - с, Л2 = V + с, т1 = V - 2с, и>2 = V + 2с,
где V = д/к — горизонтальная скорость жидкости; с = у/дК — скорость распространения малых возмущений. Будем рассматривать только такие решения системы (1.1), (6.1), которые описывают докритическое течение жидкости
Н < с ^ А1 < 0, А2 > 0.
Рассмотрим для системы (1.1), (6.1) задачу Коши с периодическими начальными данными
(2и Х 1 ( (2жХ ж\ V
у(х, 0) = а81п( — + -! , К(х, 0) = —(а 81п ( — + - ) + Ь ) , (6.2)
которым соответствуют следующие начальные значения инвариантов:
/2жх ж \
'ш1(х, 0) = -Ь, и>2(х, 0) = 2у(х, 0) + Ь = 2а вт ( + — ) + Ь,
где а = 2, Ь = 10 и X = 9.99. Точное решение этой задачи моделируется численным расчетом по схеме КАБАРЕ на мелкой сетке с пространственным шагом А = 0.005
и безразмерным параметром = 10. Профили глубины, получаемые в этих расчетах на отрезке [0, X] длины периода, показаны сплошными линиями на рис. 1, а в моменты времени Т = 0.5 и Т = 1. Кружками изображены результаты расчета по схеме КАБАРЕ на сетке с пространственным шагом Д = 0.2.
Из этих расчетов следует, что при Т ^ 0.5 в результате градиентной катастрофы в точке ж 6 отрезка [0,Х] формируется ударная волна, распространяющаяся в положительном направлении оси ж. При Т =1 эта ударная волна достигает точки ж ~ 9 и за ее фронтом образуется область влияния (след волны), лежащая внутри интервала (4,9). Аналогичные прерывные волны образуются на всех отрезках [/X, (/ + 1)Х], / € Z. К моменту времени Т = 1 . 5 эти прерывные волны проходят расстояние, превышающее длину периода, и их области влияния покрывают всю ось ж. Из рис. 1, а (справа) следует, что схема КАБАРЕ с высокой точностью локализует ударную волну, размазывая ее на два пространственных интервала разностной сетки при отсутствии нефизических осцилляций на фронте волны.
Точками на рис. 1, а показаны порядки интегральной сходимости р(Т, ж^ ,Х), определяемые по формуле (5.5), а на рис. 1, б — порядки локальной сходимости Р(ж^, Т, /, е), где I = 5, е = 3Д, определяемые по формуле (5.9). На рис. 2 в логарифмической шкале приведены относительные локальные дисбалансы
5 У (ж,- ,Т,е,/)
| ¿И (ж, Т, е)| |И (ж, Т, е)|
(6.3)
х о
х о
а
Рис. 1. Результаты численных расчетов на моменты времени Т = 0.5 (слева) и Т = 1 (справа): а — профиль глубины Н, получаемый при расчете на сетке с пространственным шагом А = 0.005 (сплошная линия, моделирующая точное решение); значения глубины Н, получаемые при расчете на сетке с шагом А = 0.2 (кружки); порядки интегральной сходимости р (точки), определяемые по формуле (5.5); б — порядки локальной сходимости Р (точки), определяемые по формуле (5.9)
вычисления модуля вектора иг (х, Т, е) осредненных базисных функций Н и д, а на рис. 3 в логарифмической шкале показаны относительные локальные дисбалансы
вычисления модулей компонент осредненного вектора инвариантов ж1 (х,Т,е). Значе-
ния параметров I и е выбирались путем проведения серии тестовых расчетов с тем, чтобы максимально уменьшить амплитуду осцилляций порядков локальной (поточечной) сходимости в области влияния ударной волны (рис. 1, б, справа).
На рис. 2 значения дисбалансов (6.3) в момент времени Т = 0.5 изображены кружками, а в момент времени Т = 1 — точками. На рис. 3 значения дисбалансов (6.4) инварианта Шх показаны треугольниками, а инварианта и^ — кружками. Расчеты интегральных и локальных порядков сходимости, а также относительных локальных дисбалансов проводились на базисной сетке с пространственным шагом А = 0.0135, что соответствует 740 пространственным ячейкам сетки на отрезке [0,Х] длины периода.
На рис. 1-3 результаты этих расчетов показаны для каждого 15-го пространственного узла ] = 15г разностной сетки. Выбор такого достаточно малого значения для пространственного шага А базисной сетки обусловлен тем, что при больших значениях этого шага, например при А = 0(0.1), получаемые по формулам Рунге порядки интегральной и локальной сходимости заметно зависят от величины А, особенно в области влияния фронта ударной волны. Расчеты показали, что начиная со значения А = 0.0135 такая зависимость при дальнейшем уменьшении пространственного шага базисной сетки становится достаточно слабой.
Из рис. 1 (слева) можно заключить, что в момент времени Т = 0.5, когда ударная волна еще не сформировалась, схема КАБАРЕ сохраняет второй порядок как интегральной, так и локальной сходимости почти во всей расчетной области. Осцилляции локальной сходимости при х ^^ 5.2 и X ~ 7.4 связаны с тем, что схема КАБАРЕ, подобно ТУО-схемам [21], может снижать порядок сходимости в окрестностях экстремумов точного решения. Для преодоления этого недостатка схемы КАБАРЕ можно использовать методы, предложенные в [31-33].
\g\6V\
-7
0
2
4
6
8
10 *
Рис. 2. Относительные локальные дисбалансы (6.3) на моменты времени Т = 0.5 (кружки), Т = 1 (точки)
Из рис. 1, а (справа) следует, что схема КАБАРЕ, подобно ТУБ- и ШЕКО-схемам [24-26], имеет первый порядок интегральной сходимости на отрезках [xj ,Х], правая граница которых расположена перед ударной волной, а левая граница — за ударной волной в области ее влияния. В результате порядок локальной сходимости в области влияния ударной волны также снижается приблизительно до первого порядка (см. рис. 1, б, справа). Снижение точности схемы КАБАРЕ в области влияния ударной волны приводит к тому, что в этой области возрастают значения локальных дисбалансов (6.3) осредненного разностного решения (рис. 2 и 3). Причем максимальное возрастание происходит не только на правой границе области влияния, которая примыкает к фронту ударной волны, но и в окрестности ее левой границы, где точное решение является достаточно гладким.
Как показано в [24, 25, 29], аналогичное снижение точности на границе области влияния ударной волны, противоположной ее фронту, демонстрируют и другие разностные схемы повышенной точности, как типа ТУБ [21], так и немонотонные схемы с гладкими функциями численных потоков такие, как схема Лакса — Вендроффа [34], схема Русанова [35], схема МакКормака [36] и компактная схема из работы [29]. Причина этого заключается в том, что в точном решении задачи Коши (1.1), (6.1), (6.2) из малой окрестности точки градиентной катастрофы выходят расходящиеся характеристики первого семейства, распространяющиеся в расширяющейся окрестности левой границы области влияния ударной волны и переносящие в эту окрестность значения инварианта и^. Поэтому такое снижение точности разностных схем сквозного счета аналогично снижению их точности при расчете центрированных волн разрежения [37].
Из рис. 3 следует, что инварианты и т2 вычисляются с различной точностью. Причем эта точность существенно снижается, когда характеристики, вдоль которых переносится инвариант, проходят через область больших градиентов (рис. 3, а) или пересекают линию фронта ударной волны (рис. 3, б). В области влияния ударной волны инвариант вычисляется с заметно более низкой точностью (рис. 3, б), поскольку именно он переносит в эту область информацию о точности, с которой разностная схема аппроксимирует условия Гюгонио на фронте волны. Поэтому снижение в данной
Рис. 3. Относительные локальные дисбалансы (6.4) инвариантов и> 1 (треугольники) и -Ш2 (кружки) на моменты времени Т = 0.5 (а), Т = 1 (б)
области точности вычисления вектора базисных переменных u = (h,q) прежде всего связано со снижением точности вычисления инварианта W\.
Таким образом, схема КАБАРЕ подобно TVD- и WENO-схемам [24-26], но в отличие от немонотонных схем Русанова и МакКормака, а также в отличие от компактной схемы из [29] имеет не более чем первый порядок точности при расчете разрывных решений с ударными волнами.
Заключение
Из результатов данной работы, а также работ [24-26, 29] следует, что в настоящее время в теории и практике разностных схем сквозного счета сложилась следующая альтернатива: невозможно одновременно с высокой точностью монотонно локализовать сильные разрывы и сохранить повышенный порядок сходимости в областях их влияния. Поэтому в случае численного моделирования сложных течений с большим числом ударных волн различной амплитуды следует использовать монотонные схемы формально повышенной точности, обладающие хорошей разрешимостью на фронтах этих волн. В то же время при необходимости с высокой точностью определять параметры течения в областях влияния изолированных ударных волн можно рекомендовать применение на достаточно мелких сетках немонотонных схем повышенной точности с гладкими функциями численных потоков (типа схемы Русанова), обеспечивающих в этих областях скорость сходимости на несколько порядков выше, чем схемы типа TVD.
Благодарности. Работа выполнена при финансовой поддержке РНФ (грант № 16-1110033).
Список литературы / References
[1] Iserles, A. Generalized leapfrog methods // IMA J. of Numer. Anal. 1986. Vol. 6, No 3. P. 381-392.
[2] Головизнин В.М., Самарский А.А. Разностная аппроксимация конвективного переноса с пространственным расщеплением временной производной // Матем. моделирование. 1998. Т. 10, № 1. С. 86-100.
Goloviznin, V.M., Samarskii, A.A. Finite difference approximation of convective transport equation with space splitting time derivative // Matematicheskoe Modelirovanie. 1998. Vol. 10, No 1. P. 86-100. (In Russ.)
[3] Головизнин В.М., Самарский А.А. Некоторые свойства разностной схемы "Кабаре" // Матем. моделирование. 1998. Т. 10, № 1. С. 101-116.
Goloviznin, V.M., Samarskii, A.A. Some characteristics of finite difference scheme "cabaret" // Matematicheskoe Modelirovanie. 1998. Vol. 10, No 1. P. 101-116. (In Russ.)
[4] Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений. М.: Наука, 1978. 687 c.
Rozhdestvenskii, B.L., Yanenko, N.N. Systems of quasilinear equations and their applications to gas dynamics. Moscow: Nauka, 1978. 687 p. (In Russ.)
[5] Куликовский А.Г., Погорелов Н.В., Семенов Ф.Ю. Математические вопросы численного решения гиперболических систем уравнений. М: Физматлит, 2001. 607 c. Kulikovskii, A.G., Pogorelov, N.V., Semenov, A.Yu. Mathematical Aspects of Numerical Solution of Hyperbolic Systems. Moscow: Fizmatlit, 2001. 607 p. (In Russ.)
[6] Головизнин В.М. Балансно-характеристический метод численного решения уравнений газовой динамики // Докл. АН. 2005. Т. 403, № 4. С. 459-464.
Goloviznin, V.M. Balanced characteristic method for systems of hyperbolic conservation laws // Dokl. Math. 2005. Vol. 72, No 1. P. 619-623.
[7] Woodward, P., Colella, P. The numerical simulation of two-dimensional fluid flow with strong shocks // J. Comput. Phys. 1984. Vol. 54, No. 1. P. 115-173.
[8] Karabasov, S.A., Goloviznin, V.M. New efficient high-resolution method for nonlinear problems in aeroacoustics // AIAA J. 2007. Vol. 45, No. 12. P. 2861-2871.
[9] Semiletov, V., Karabasov, S.A. CABARET scheme with conservation-flux asynchronous time-stepping for nonlinear aeroacoustics problems //J. Comput. Phys. 2013. Vol. 253. P. 157-165.
[10] Faranosov, G., Goloviznin, V., Karabasov, S., Kondakov, V., Kopiev, V., Zaitsev, M. CABARET method on unstructured hexahedral grids for jet noise computation // Comput. Fluids. 2013. Vol. 88. P. 165-179.
[11] Karabasov, S.A., Berloff, P.S., Goloviznin, V.M. Cabaret in the ocean gyres // Ocean Modelling. 2009. Vol. 30. No. 2. P. 155-168.
[12] Новые алгоритмы вычислительной гидродинамики для многопроцессорных вычислительных комплексов / В.М. Головизнин, М.А. Зайцев, С.А. Карабасов, И.А. Короткин. М.: Изд-во Моск. ун-та, 2013. 467 c.
New CFD algorithms for multiprocessor computer systems / V.M. Goloviznin, M.A. Zaitsev, S.A. Karabasov, I.A. Korotkin. Moscow: Izd-vo Mosk. Univ., 2013. 467 p. (In Russ.)
[13] Ковыркина О.А., Остапенко В.В. О монотонности двухслойной по времени схемы Кабаре // Матем. моделирование. 2012. Т. 24, № 9. С. 97-112.
Kovyrkina, O.A., Ostapenko, V.V. On monotony of two layer in time cabaret scheme // Math. Model. Comput. Simul. 2013. Vol. 5, No. 2. P. 180-189.
[14] Ковыркина О.А., Остапенко В.В. О монотонности схемы КАБАРЕ в многомерном случае // Докл. АН. 2015. Т. 462, № 4. С. 385-390.
Kovyrkina, O.A., Ostapenko, V.V. On the monotonicity of the CABARET scheme in the multidimensional case // Dokl. Math. 2015. Vol. 91, No. 3. P. 323-328.
[15] Зюзина Н.А., Остапенко В.В. О монотонности схемы КАБАРЕ, аппроксимирующей скалярный закон сохранения с выпуклым потоком // Докл. АН. 2016. Т. 466, № 5. С. 513-517.
Zuzina, N.A., Ostapenko, V.V. On the monotonicity of the CABARET scheme approximating a scalar conservation law with a convex flux // Dokl. Math. 2016. Vol. 93, No. 1. P. 69-73.
[16] Остапенко В.В. О сходимости разностных схем за фронтом нестационарной ударной волны // Журн. вычисл. математики и матем. физики. 1997. Т. 37, № 10. С. 1201-1212. Ostapenko, V.V. Convergence of finite-difference schemes behind a shock front // Comput. Math. and Math. Phys. 1997. Vol. 37, No. 10. P. 1161-1172.
[17] Casper, J., Carpenter, M.H. Computational consideration for the simulation of shock-induced sound // SIAM J. Sci. Comput. 1998. Vol. 19, No. 1. P. 813-828.
[18] Engquist, B., Sjogreen, B. The convergence rate of finite difference schemes in the presence of shocks // SIAM J. Numer. Anal. 1998. Vol. 35. P. 2464-2485.
[19] Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Матем. сборник. 1959. Т. 47. № 3. С. 271-306.
Godunov, S.K. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics // Matematicheskii Sbornik. 1959. Vol. 47, No. 3. P. 271-306. (In Russ.)
[20] Van Leer, B. Toward the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method // J. Comput. Phys. 1979. Vol. 32, No. 1. P. 101-136.
[21] Harten, A. High resolution schemes for hyperbolic conservation laws //J. Comput. Phys. 1983. Vol. 49. P. 357-393.
[22] Nessyahu, H., Tadmor, E. Non-oscillatory central differencing for hyperbolic conservation laws //J. Comput. Phys. 1990. Vol. 87, No. 2. P. 408-463.
[23] Jiang, G.S., Shu, C.W. Efficient implementation of weighted ENO schemes //J. Comput. Phys. 1996. Vol. 126. P. 202-228.
[24] Ковыркина О.А., Остапенко В.В. О сходимости разностных схем сквозного счета // Докл. АН. 2010. Т. 433, № 5. С. 599-603.
Kovyrkina, O.A., Ostapenko, V.V. On the convergence of shock-capturing difference schemes // Dokl. Math. 2010. Vol. 82, No. 1. P. 599-603.
[25] Ковыркина О.А., Остапенко В.В. О реальной точности разностных схем сквозного счета // Матем. моделирование. 2013. Т. 25, № 9. С. 63-74.
Kovyrkina, O.A., Ostapenko, V.V. On the practical accuracy of shock-capturing schemes // Math. Models and Comput. Simulations. 2014. Vol. 6, No. 2. P. 183-191.
[26] Михайлов Н.А. О порядке сходимости разностных схем WENO за фронтом ударной волны // Матем. моделирование. 2015. Т. 27, № 2. С. 129-138.
Mikhailov, N.A. On convergence rate of WENO schemes behind a shock front // Math. Models and Comput. Simulations. 2015. Vol. 7, No. 5. P. 467-474.
[27] Самарский А.А. Теория разностных схем. М.: Наука, 1983. 616 c.
Samarsky, A.A. The theory of difference schemes. Moscow: Nauka, 1983. 616 p. (In Russ.)
[28] Остапенко В.В. О конечно-разностной аппроксимации условий Гюгонио на фронте ударной волны, распространяющейся с переменной скоростью // Журн. вычисл. математики и матем. физики. 1998. Т. 38, № 8. С. 1355-1367.
Ostapenko, V.V. Finite-difference approximation of the Hugoniot conditions on a shock front propagating with variable velocity // Comput. Math. and Math. Physics. 1998. Vol. 38, No. 8. P. 1299-1311.
[29] Остапенко В.В. О построении разностных схем повышенной точности для сквозного расчета нестационарных ударных волн // Журн. вычисл. математики и матем. физики. 2000. Т. 40, № 12. С. 1857-1874.
Ostapenko, V.V. Construction of high-order accurate shock-capturing finite difference schemes for unsteady shock waves // Comput. Math. and Math. Physics. 2000. Vol. 40, No. 12. P. 1784-1800.
[30] Стокер Дж.Дж. Волны на воде. Математическая теория и приложения. М.: Изд-во иностр. лит., 1959. 617 c.
Stoker, J.J. Water waves: The mathematical theory with applications. N.Y.: Wiley, 1957. 600 p.
[31] Карабасов С.А. О возможностях методов второго порядка аппроксимации на примере модельных задач газо- и гидродинамики // Матем. моделирование. 2010. Т. 22, № 7. С. 93-120.
Karabasov, S.A. On the powerof second-order accurate numerical methods for model problems of gas- and hydrodynamics // Math. Models and Comput. Simulations. 2011. Vol. 3, No. 1. P. 92-112.
[32] Зюзина Н.А., Остапенко В.В. Модификация схемы кабаре, обеспечивающая ее сильную монотонность и повышенную точность на локальных экстремумах // Докл. АН. 2014. Т. 457, № 3. С. 268-273.
Zuzina, N.A., Ostapenko, V.V. Modification of the CABARET scheme ensuring its strong monotonicity and high accuracy on local extrema // Dokl. Math. 2014. Vol. 90, No. 1. P. 453-457.
[33] Зюзина Н.А., Остапенко В.В. Модификация схемы КАБАРЕ, обеспечивающая ее повышенную точность на локальных экстремумах // Матем. моделирование. 2015. Т. 27, № 10. С. 21-31.
Zuzina, N.A., Ostapenko, V.V. Modification of the Cabaret scheme ensuring its high accuracy on local extrema // Math. Models and Comput. Simulations. 2016. Vol. 8, No. 3. P. 231-237.
[34] Lax, P., Wendroff, B. Systems of conservation laws // Commun. on Pure and Appl. Math. 1960. Vol. 13. P. 217-237.
[35] Русанов В.В. Разностные схемы третьего порядка точности для сквозного счета разрывных решений // Докл. АН СССР. 1968. Т. 180, № 6. С. 1303-1305.
Rusanov, V.V. Difference schemes of third order accuracy for the through calculation of discontinuous solutions // Dokl. AN SSSR. 1968. Vol. 180, No. 6. P. 1303-1305. (In Russ.)
[36] MacCormack, R.W. The effect of viscosity in hypervelosity impact cratering // AIAA Paper. 1969. P. 69-354.
[37] Остапенко В.В. О сильной монотонности нелинейных разностных схем // Журн. вы-числ. математики и матем. физики. 1998. Т. 38, № 7. С. 1170-1185.
Ostapenko, V.V. On the strong monotonicity of nonlinear difference schemes // Comput. Math. and Math. Physics. 1998. Vol. 38, No. 7. P. 1119-1133.
Поступила в редакцию 17 сентября 2017 г., с доработки — 9 октября 2017 г.
On monotonicity and accuracy of CABARET scheme for calculation of weak solutions with shocks
KOVYRKINA, OLYANA A.1'*, OSTAPENKO, VLADIMIR V.1'2
1Lavrentyev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia 2 Novosibirsk State University, Novosibirsk, 630090, Russia * Corresponding author: Kovyrkina, Olyana A., e-mail: olyana@ngs.ru
We studied the monotonicity and the accuracy of the modified CABARET scheme approximating the quasilinear hyperbolic system of conservation laws. We obtained conditions under which this scheme preserves the monotonicity of the finite-difference solution with respect to invariants of the linear approximation of the considering system. As a specific example, we considered the approximation for the system of conservation laws for shallow water theory.
To determine the accuracy of the scheme in the regions influenced by shocks that propagate with variable velocity and behind the fronts of which a non-constant solution is formed, series of test calculations on the sequence of contracting grids were carried out. These calculations enabled us to apply the Runge rule to define the order of convergence of the finite-difference solution. To estimate the accuracy of the modified CABARET scheme, the orders of its integral and local convergence are calculated, as well as local disbalances in calculating of the absolute value of averaged basis functions
© ICT SB RAS, 2018
vector and local disbalances in calculating of the absolute value of the components of the averaged vector of invariants. The order of the integral convergence makes it possible to estimate the accuracy of finite-difference schemes in translation of the Rankine-Hugoniot conditions across a non-stationary shock wave front. In this case, the order of integral convergence is calculated for the finite-difference solution itself (rather than for its absolute value, as in the L1 norm) on spatial intervals crossing the shock wave.
It is shown that, like TVD-schemes of high order approximation on smooth solutions, the CABARET scheme being formally of the second order, in spite of high accuracy in the localization of shocks, reduces the order of integral convergence to the first order on the integration intervals, one of the boundaries of which is in the region influenced by a shock. The main reason of this decrease in accuracy is that the flux correction used in the CABARET scheme to its monotonization leads to the decrease in the smoothness of the finite-difference fluxes, which in turn leads to the decrease in the order for approximation of Rankine-Hugoniot conditions on the shocks. As a result, the local accuracy of the CABARET scheme in the regions of shocks influence also decreases to about the first order.
Keywords: hyperbolic system of conservation laws, monotonicity and accuracy of CABARET scheme, shallow water flows, shocks.
Acknowledgements. This research was supported by the Russian Science Foundation (grant No. 16-11-10033).
Received 17 September 2017 Received in revised form 9 October 2017
