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

АПВПМ-2019

СРАВНЕНИЕ РАЗНОСТНЫХ СХЕМ РАЗЛИЧНЫХ ПОРЯДКОВ ТОЧНОСТИ НА АДАПТИВНЫХ СЕТКАХ ДЛЯ РАСЧЕТА ПОГРАНИЧНЫХ И ВНУТРЕННИХ СЛОЕВ

В, Д, Лиеейкин1,2, В, И, Паасонен1,2

1 Институт, вычислительных технологий СО РАН, 630090, Новосибирск 2 Новосибирский государственный университет, 630090, Новосибирск

УДК 519.6

Б01: 10.24411/9999-016А-2019-10045

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

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

Введение

Рассматривается первая краевая задача для обыкновенного дифференциального уравнения

-ем" + а(х)и' + ^(х,и)= 0, ж € (0, 1), ^(и,х) > 0, м(0) = и0, и( 1) = и1; (1)

где £ > 0 малый параметр.

Для численного решения задачи (1), являющейся модельной для многих задач механики, связанных с исследованием ударных переходов, пограничных и внутренних слоев и других узких зон чрезвычайно быстрого изменения искомых функций, рассматриваются три различные схемы на трехточечном шаблоне, аппроксимирующие уравнение соответственно с первым, вторым и третьим порядком относительно локального значения шага неравномерной сетки 0 = х0 < х\ < ... < х^ = 1, являющейся образом равномерной сетки при гладком отображении.

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

Обобщение алгоритма задания сеток для случая высокоточных схем построено по аналогии на основе априорных оценок старших производных решения и конструирования для генерации сеток преобразований необходимого класса гладкости [4]. При этом поскольку схемы порядков аппроксимации выше первого не

Работа выполнена при частичной финансовой поддержке Российского научного фонда (код проекта 17-72-30006).

!ЯВ.\ 978-5-901548-42-4

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

1 Компактные разностные схемы

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

-ели + аа±и + ^ = 0. (2)

первого порядка точности со стандартной аппроксимацией второй производной

д+и — Д_м

+ 7 7

ли = -—-, s = h+ + h_

и односторонней разделенной разностью а±и, выбираемой в зависимости от знака а(хг). Символами к+ и к- обозначены местные значения шагов сетки справа и слева от данного узла.

Во второй схеме первая производная аппроксимируется по трем узлам со вторым порядком, а главный член погрешности разностного аналога второй производной ли компенсируется дополнительными слагаемыми так, чтобы в целом получилась схема второго порядка аппроксимации

-ели + ааи + ^пи + (Е + ^А)^ = 0, й = к+ - к-, (3)

где

к- А+ + к+а- ^ а+а+ - а- А_ А =-+-+-, П = + ,-, в = к+ + к-,

я в/2

Е — ^^^^^^^^етпый оператор, а± — местные значения функции а(х) в полуцелых узлах справа и слева.

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

ю ю &

(-£ +-аа)ли + (а + — Ла)Ам +(--аг)пи + = 0, р = к_к+, (4)

6 12 3

где

V и , Р + 3р + ¿2

Х = Е + (--аг)а+--Л, г =-.

^3 ; 12 ' 36е

Построенные схемы являются модельными по отношению к схемам для решения нелинейных двумерных задач [5], [6], описывающих стационарные течения вязкой несжимаемой жидкости.

Следует ожидать что в данном случае порядок аппроксимации схемы (4) будет выше на единицу, а именно четвертым, так как для квазиравномероной сетки (т. е. являющейся образом равномерной при гладком преобразовании х(£,е) € С2(0, 1)) разность соседних шагов й = к+ - к- имеет второй порядок малости. Для обеспечения равномерной по е сходимости разностного решения к точному требуется равномерная по £ оценка производпьк решения относительно новой сеточной координаты £

dpu(x(£, е), е)

dp

< М, р < п (5)

с константой М, те зависящей от е.

В тестовых задачах под «точным» решением на приводимых графиках подразумевается численное решение, полученное на самой подробной сетке в процессе ее детализации. Для оценки С-нормы ошибки на сгущающихся сетках при неизвестном точном решении использовалась разность приближенных решений на двух последовательных сетках:

Sh = max |u2f — uf |, h =1/N.

i=0

где uN и u2N — приближенные решения на сетках с N и 2N шагами. По оценкам ошибки вычислялась апостериорная оценка реально наблюдаемого порядка точности ph = log2(Sh/Sh/2). Во всех таблицах оценки

точности схем, по которым получены результаты.

Важной характеристикой алгоритма явного задания узлов адаптивной сетки является равномерность модуля приращения решения на шаге, т. е. оценка |u^j —u^ | < m/N, где m те зависит от е и N. Практически во всех расчетах это свойство имело место. Нарушалось оно лишь при грубой сетке, не достаточно детальной для данного значения вязкости, и обычно проявлялось в осциллирующем характере приближенных решений.

2 Базовое преобразование, генерирующее адаптивную сетку

Для построения разностных сеток, сгущающихся в экспоненциальных слоях возле точки х0 = 0, используем преобразование класса Сг[0, 1], устраняющее в слое экспоненциальные особенности масштаба к до порядка п. Базовое преобразование, отображающее взаимно-однозначно отрезок 0 < £ < 1 та отрезок 0 < х < 1, строится следующим образом. На части отрезка 0 < £ < £0 определяется предложенная в [7] функция

* = ф(0 = ^((1 -<%)-1/а—1), 0<е<Со, (в)

где <1 = (1 — )/£о > 1 + Ш1 > 1, Р = а/(1 + па), при этом а — произвольная положительная константа, отделенная от нуля (а > т2 > 0), а к — масштаб слоя. Функция (6) устраняет экспоненциальную особенность

0 < < 1

этого полинома достаточно знать значения функции (6) и ее производных до порядка I в точке £ = £0, которые выражаются явно

= ^ ^С^ + Л ...(— ^ £ка{п-г)/{1+па) г> 1. /?ч

а\а / \а )

1

растяжением (сжатием) обеспечить равенство единице значения склеенной функции на правом конце £ = 1. В результате получается составное базовое преобразование

С1Ф(0 , 0 < £ < £о ,

жь(£,е,а,к,1,п, £о) =\ Г ' „л п (8)

, Ь < £ < 1

где у = £ — £0, I < п, с0 > Константа с0 > 0 обеспечивает строгую монотонность составной функции, а с1 > 0 выбир^тся из условия хь(1, £, а, к,...) = 1.

Преобразование (8) учитывает оценки старших производных и является обобщением на случай п > 2

0

долю отрезка, относящуюся к ветви (6) составного преобразования отображаемую в погранслой. В случае равномерной сетки по £ число £0 означает долю числа шагов сетки, попадающих в погранслой. Числа п (по-

в зависимости от порядка аппроксимации схемы р (обычно I = р, п = р + 1).

Преобразование (8) является довольно универсальным для устранения экспоненциальных особенностей функции и(х, е) до порядка п, так как его константа а >т> 0 годится при любых значениях а(0). Кроме того, преобразование (8) при некоторых ограничениях на константу а пригодно также и для устранения степенных сингулярностей.

3 Результаты численных исследований

Задача 1. Рассмотрим задачу с экспоненциальным слоем масштаба к = 1 в окрестности левой границы х = 0, который реализуется щи всюду отрицательном а(х), отделенном от нуля константой (см. [1] и [3]).

Таблица 1: Результаты численного решения задачи 1 по схемам 0(кр) (р = 1, 2, е = 0.001

N Р1 ¿2 Р2 ¿3 Р3

80 2.40е-02 5.384 1.51е-02 6.053 1.62е-02 5.951

160 1.20е-02 0.991 3.44е-03 2.129 3.40е-03 2.250

320 б.Юе-03 0.982 8.28е-04 2.056 1.18е-03 1.531

640 3.06с-03 0.998 2.05е-04 2.015 2.09е-04 2.489

1280 1.53е-03 0.997 5.11е-05 2.004 1.58е-05 3.728

2560 7.66е-04 0.999 1.28е-05 2.001 9.84е-07 4.005

Сетка определяется с помощью функции (8) при параметрах к = 1, а = 2, Со = 1/2. Краевая задача имеет

—ей" — (1 + х)и' + хи + ес8(4этж) = 0, м(0) = 0, и(1) = 1.

В таблице 1 сравниваются расчеты по трем схемам на последовательности сеток при е = 1е — 03. В первой колонке приведено число шагов, затем пара колонок содержит апостериорные оценки С-нормы ошибки ^ и порядка точности р\, вычисленные по результатам, полученным на двух последовательных сетках по схеме первого порядка аппроксимации. Следующие две пары колонок имеют тот же смысл, но в отношении схем второго и соответственно третьего порядков аппроксимации. На рисунке 1 приведены результаты расчетов

Рис. 1: Сравнение результатов расчетов при N = 20 (слева) и N = 40 (справа). Вязкость е = 0.001.

задачи 1 на сетках с числом шагов N = 20 (слева) и N = 40 (справа) по схемам первого, второго и третьего порядков.

Как показывают результаты таблицы 1, на грубых сетках схема третьего порядка аппроксимации может давать не лучшую точность, но при детализации сетки расчеты по схеме третьего порядка становятся более точными. Это объясняется тем, что в главном члене погрешности порядок производной выше у схем более высокого порядка точности, а локальные значения производных высших порядков в погранслое могут быть существенно большими, чем значения младших производных, и различие компенсируется только детализацией сетки за счет более высокого показателя степени в главном члене погрешности 0(кр) схем повышенной точности. Именно это происходит на практике.

Заметим, что если погранслой расположен не на левой, а на правой границе (реализуется при отделенном от нуля положительном а(х)), то в качестве отображения, генерирующего сетку, можно взять функцию

у(£, £, а, к, I, п, Со) = 1 — хь(1 — £, а, к, I, п, Со) (9)

Задача 2. Рассмотрим случай, когда слои одинакового масштаба к = 1 расположены на концах отрезка (0,1). Такая конфигурация реализуется при а(0) < 0 а(1) > 0 и а!(ж0) > 0 если а(ж0) = 0. Например, это имеет место при возрастающей функции а(х), меняющей знак внутри отрезка (0,1). Для расчета была взята следующая задача

—ей" + (х — и' + и + 8т(2^ж) = 0, м(0) = ^ , «(1) = 0.

Таблица 2: Результаты численного решения задачи 2 по схемам 0(ЬР) (р = 1, 2, е = 1е — 03

N ¿1 Р1 ¿2 Р2 ¿3 Р3

40 8.33е-02 2.815 9.39е-02 2.789 2.60е-01 1.013

80 6.08е-02 0.454 2.05е-02 2.192 3.75е-02 2.797

160 3.77е-02 0.687 5.04е-03 2.028 4.11е-03 3.188

320 2.15е-02 0.808 1.25е-03 2.010 2.63е-04 3.966

640 1.14е-02 0.915 3.12е-04 2.001 1.63е-05 4.011

1280 5.88е-03 0.960 7.81е-05 2.001 1.01е-06 4.004

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

Ра0шеа т т аашоёайё паоёа

Рис. 2: Решение задачи 2 на сетке с N = 20 (слева) и Ж = 40 (справа) при е = 0.001 разными схемами

Это можно осуществить как суперпозицией левостороннего и правостороннего преобразований (8). (9). так и отображением, склеенным из двух отображений, построенных на основе базовых. Последнее имеет вид

*(0 = 2(1+^))> V = 2£ — 1, 0 <£< 1, _( ) [ — У(—V,■■■), —1 <Ч<0,

\ у(ч,.■■), 0 <Ч< 1,

а у ................ отображение (9). В нашем примере параметры отображения а = 2, Со = 0.5. На рисунке 2 представлены

= 0.001

собой, заметно некоторое размазывание решения монотонной схемой первого порядка точности.

Более подробное представление о результатах расчета задачи 2 на последовательности адаптивных сеток можно составить по приведенным в таблице 2 данным.

Задача 3. Для построения сеток для задач со степенными пограничными слоями можно использовать

а

задаче с экспоненциальным погранслоем.

Пусть в нашем уравнении Ри(х,и) > с > 0 а(0) = 0 а'(0) > 0 а(1) < 0 и а'(хо) > 0 для каждого хо, для которого а(х0) = 0, 0 < х0 < 1. В этом случае производные решения, согласно [1] и [3], оцениваются следующей формулой:

|и(р)(х, е)| < М [1 + еь/2(е1/2 + х)-Ь-р] , р <п, 0 <х < 1 , (10)

где 0 < т < Ь < с/а'(0), т. е. решение имеет единственный степенной погранслой масштаба 1/2 возле хо = 0. Поэтому для построения сетки, сгущающейся в слое, можно использовать преобразование хь(£,£,а, 1/2,...) из (8), в котором 0 < т < а < Ъ/п2.

В качестве примера рассматривается следующая задача

—ей" — 8х(х — 1/2)3и' + 16м — exp(ж) = 0, м(0) = 0 , и(1) = 0.2 .

Для этого примера с = 16 о!(0) = 1, поэтому можно положить

I = 2 , Ъ =16 , £0 = 1/2 , а = 16/п2 , р = 16/(п(16 + п)).

Результаты расчета степенного слоя (задача 3) на сгущающихся равномерных и адаптивных сетках приведены в таблице 3.

Таблица 3: Результаты расчета степенного слоя масштаба 1/^и е = 0.001.

равномерная сетка адаптивная сетка

N ¿1 ¿2 ¿3 ¿1 ¿2 ¿3

20 9.87е-04 2.24е-03 1.40е-02 2.00е-03 7.07е-04 7.03е-03

40 2.87е-03 2.03е-03 8.06е-03 8.61е-04 3.62е-04 5.61е-03

80 3.09е-03 2.46е-03 2.15е-03 8.72е-04 З.ЗЗе-04 1.07е-03

160 1.88е-03 1.54е-03 2.87е-04 5.58е-04 7.83е-05 1.55е-04

320 6.00е-04 4.70е-04 1.99е-05 3.21е-04 2.30е-05 1.36е-05

640 1.97е-04 1.27е-04 1.29е-06 1.75е-04 6.53е-06 9.06е-07

1280 6.73е-05 3.24е-05 8.15е-08 9.13е-05 1.65е-06 7.02е-08

Заключение

Формат данной работы не позволяет привести более полно результаты многочисленных тестов, проведенных авторами, однако характер зависимостей ошибок от порядка точности схем во всех расчетах примерно такой же, как в приведенных тестах. Помимо рассмотренных здесь задач были были проведены расчеты экспоненциальных слоев с различными масштабами, внутренних слоев, логарифмических и смешанных [4]. Численные эксперименты свидетельствуют об эффективности симбиоза компактных схем и специальных адаптивных сеток, явно задаваемых на основе априорных оценок производных решения. В большинстве расчетов реальная точность тем выше, чем выше порядок точности схем, аномалии наблюдаются лишь при весьма грубых сетках. Это объясняется тем, что используемые оценки производных носят асимптотический характер при к ^ 0, и рекомендации по распределению узлов сетки тем более точны, чем детальнее сетка. Кроме того, значения констант в оценках производных высших порядков неизвестны, в то время как они могут влиять на ошибку при крупном шаге сетки.

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

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

[1] Лисейкин В. Д. Оценки производных решения дифференциальных уравнений с пограничными и внутренними слоями // Сибирский математический журнал. 1992. Т. 33, № 6, С. 106-117.

[2] Liseikin V. D. Grid Generation Methods// Springer, third ed., Berlin, 2017.

[3] Liseikin V. D. Layer Resolving Grids and Transformations for Singular Perturbation Problems// VSP, Utrecht, 2001.

[4] Лисейкин В. Д., Паасонен В. И. Компактные разностые схемы и адаптивные сетки для численного моделирования задач с пограничными и внутренними слоями // Сиб. Журн. Вычисл. математики. 2019. Т. 22, № 1, С. 41-56.

[5] Паасонен В. И. Схема третьего порядка аппроксимации на неравномерной сетке для уравнений Навье-Стокса// Вычислительные технологии. 2000. Т. 5, № 5, С. 78-85.

[6] Глуховский А. С., Паасонен В. И. Компактные разностные схемы для уравнений Навье-Стокса на неравномерных сетках // Марчуковские научные чтения 2017. Труды международной конференции «Вычислительная и прикладная математика 2017». Изд. ИВМиМГ. Новосибирск. 25 — 30 июня 2017. С. 211-217.

[7] Лисейкин В. Д. О численном решении уравнений со степенным пограничным слоем// Журн. Вычисл. матем. и матем. физики. 1986. Т. 26, № 12, С. 1813-1820.

[8] Лисейкин В. Д. О численном решении сингулярно возмущенных уравнений с точками поворота // Журн. Вычисл. матем. и матем. физики. 1984. Т 24, № 12, С. 1812-1818.

Лисейкин Владимир Дмитриевич —д.ф.-м.н., вед. науч.сотр. Института

вычислительных технологий СО РАН;

e-mail: [email protected];

Паасонен Виктор Иванович — к.ф.-м.н., ст. науч. сотр. Института

вычислительных технологий СО РАН;

e-mail: paasQict.nsc.ru;

Дата поступления — 30 апреля 2019 г.