Вычислительные технологии, 2020, том 25, № 5, с. 66^79. © ФИЦ ИВТ, 2020 Computational Technologies, 2020, vol. 25, no. 5, pp. 66-79. © FRC ICT, 2020
ISSN 1560-7534 elSSN 2313-691X
ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ
D01:10.25743/ICT.2020.25.5.006
Характер сходимости схем при расчете на адаптивных сетках задач со слоями
В. Д. ЛИСЕЙКИН1'*, В. И. ПААСОНЕН1'2
1 Федеральный исследовательский центр информационных и вычислительных технологий СО РАН, Новосибирск, Россия
2 *
Поступила 28 февраля 2020 г., доработана 2 августа 2020 г., принята в печать 28 августа 2020 г.
Проведено сравнение качества решений модельного уравнения второго порядка с малым параметром, полученных по трем различным разностным схемам на специальных адаптивных сетках, явно задаваемых координатным преобразованием, а также на равномерных сетках в новых переменных, соответствующих этому преобразованию. Исследуются схемы второго порядка точности с диагональным преобладанием и без него и простейшая противопотоковая схема. На основе оценок погрешностей сделаны прогнозы относительно свойств решений, подтвержденные анализом и численными экспериментами. Показано, что схема второго порядка аппроксимации с диагональным преобладанием сходится равномерно по малому параметру со вторым порядком лишь в частном случае, когда коэффициент при старшей производной мал только в слое; если же он мал также и вне слоя, порядок сходимости первый. Установлено также, что схема без диагонального преобладания имеет существенно более качественные решения без осцилляций в новых переменных на равномерной сетке, чем в соответствующих им исходных физических координатах. В противоположность ей схемы с диагональным преобладанием не чувствительны к выбору системы координат.
Ключевые слова: равномерная сходимость, адаптивные сетки, пограничный слой, диагональное преобладание, малый параметр.
Цитирование: Лисейкин В.Д., Паасонен В.И. Характер сходимости схем при расчете на адаптивных сетках задач со слоями. Вычислительные технологии. 2020; 25(5):66-79. Б01:10.25743/1СТ.2020.25.5.006.
Введение
Работа посвящена исследованию равномерной по малому параметру сходимости разностных схем для модельной задачи
-ф,е)ихх + а(х)их + ^(х,и) = 0, ж е (0,1) и(0) = и0, и(1) = иъ в линейном случае имеющей вид
-ф,е)ихх + а(х)их + с(х)и = /(х), х е (0,1) и(0) = и0, и(1) =
Здесь и — решение дифференциальной задачи, ее производная ^ > с0 > 0 (или, в линейном случае, функция с(х) > с0 > 0) отделена от нуля снизу, а функция у,(х,е) > 0 стремится к нулю хотя бы для одного фиксированного значения х € [0,1], когда стремится к нулю малый параметр е > 0, Для таких задач характерно наличие пограничных и внутренних слоев различных типов, для успешного расчета которых целесообразно использовать специальные неравномерные сетки, сгущающиеся в слоях по особым законам, Способам построения специальных сеток и их применению для решения задач с малым параметром посвящено множество исследований (см., например, работы [1-4] и помещенные в них обзоры).
Существует множество различных разностных схем на равномерных и неравномерных сетках для решения поставленной задачи, однако входящая в оценку ошибки константа обычно зависит от нижней границы изменения ^(х, е), стремясь к бесконечности, если хотя бы для одного значения х € [0,1] значение ^(х,е) ^ 0 при е ^ 0, и в итоге сходимость схемы может оказаться неравномерной по малому параметру е. Для расчета слоев более предпочтительными принято считать схемы, равномерно сходящиеся при стремлении к нулю малого параметра, когда константа в оценке ошибки остается ограниченной сверху при любых е ^ 0, Для большинства схем на равномерных и не только на равномерных сетках достичь этого не представляется возможным. Однако для отдельных схем (например, для простейшей противопотоковой схемы первого порядка аппроксимации) на специальных адаптивных сетках с "удачным" масштабированием в окрестности слоев это возможно [1-3]. На неравномерной сетке противопото-ковая схема имеет вид
, . . а(х) + |а(ж)| . а(х) — |а(ж)| . , . „. .
х, е)Ли + 1 ; 2 и + 2 А+и + С(Х)и = ^(Х), ^
где и — решение разностной задачи; х — произвольный внутренний узел сетки; Л = 2(Д+ — Д_)/в — разностный аналог оператора двойного дифференцирования; Д+ и Д- — операторы правой и левой разделенных разностей; в = + Н_ — сумма местных значений шага сетки справа Н+ и слева Н- от данного фиксированного узла.
На специальных сетках, учитывающих специфику поведения решения в слоях, в работе [4] при постоянном коэффициенте ^(х,е) = £ теоретически доказана равномерная по малому параметру сходимость схемы (1) с первым порядком. Она следует из оценок погрешности ||Ф|| < сопst/(М^,)) в слое и ||Ф|| < сопвне стоя, где N — число шагов сетки, и оценок устойчивости в смысле выполнения неравенства ||и|| < eonst ^||Ф|| в слое и неравенства ||и|| < сопэ^ФЦ вне слоя. При этом в качестве главного момента использовался тот факт, что при индексной записи схемы в виде
— АгЧг+1 + ВгЩ — Сш_1 = г = 1,. . . ,М — 1,
коэффициенты схемы А^ ВС г положительны и выполняется условие строгого диагонального преобладания В^ > А^ + С^
Для большинства схем второго порядка, построенных на основе комбинаций левой и правой разделенных разностей, например для схемы
. /\А /\ А Д— + , .
—^(х,е)Ли + а(х)Ди + с(х)и = / (х), Д =-, (2)
а также для трехточечных схем третьего порядка точности, это свойство не выполняется при произвольных шагах, хотя численное исследование на последовательности
специальных сгущающихся сеток, проведенное в работе [5], подтвердило не только их практическую применимость для расчета слоев, но даже определенное преимущество перед схемой с направленными разностями. Тем не менее отсутствие диагонального преобладания дает основание сомневаться в их равномерной сходимости. На практике это проявляется в том, что на грубых сетках решения схемы (2) осциллируют и по реальной точности они хуже решений противопотоковой схемы (1), но при увеличении числа узлов сетки постепенно достигается преимущество схемы (2), При этом порог детализации сетки, необходимой для достижения превосходства схемы (2) над противопотоковой схемой (1), отодвигается все более по мере уменьшения вязкости.
Однако среди схем второго порядка точности есть уникальные примеры схем с диагональным преобладанием. Это схема Ильина [6] и схема Булеева [7], которые без преувеличения являются гениальными модификациями двух простейших схем, формально реализованными введением специальных множителей при разностной аппроксимации второй производной. Первая схема предназначена исключительно для равномерных сеток, по построению гарантирует близость разностных решений к типичным точным решениям в слое экспоненциального типа и одновременно превращает схему с центральной разностью (т. е, схему (2) при Ь+ = Ь- = Ь) в схему с диагональным преобладанием, сохраняя при этом второй порядок аппроксимации. Вторая схема модифицирует про-тивопотоковую схему (1), изящно компенсируя главный член разложения погрешности односторонней аппроксимации конвективного слагаемого и одновременно сохраняя структуру и свойства исходной схемы (1), Схема Булеева на равномерной сетке имеет вид
1 . о>(%) + л а(х) — |а(ж)| . . . „. . , .
-Ф,е)^ „, ,-^и + ( ) 1 ( )| А-и + ( ) ' ( )1 Д+и + с(х)и = Кх), 3
1 + Н/1л(х, е) 2 2
где Н = |а(ж)|Ь/2, В случае неравномерной сетки схема Булеева сохраняет вид (3), если условиться под символом Ь подразумевать шаг слева Ь- при положительном а(х) и шаг справа Ь+ — при отрицательном а(х). По существу эта схема также противопото-ковая, но с подправленным выражением для коэффициента вязкости, В представлении схемы (3) в индексном виде ввиду очевидного неравенства |а|±а > 0 легко устанавливаются положительность А^ и С г и тождеет во В г = А^ + С + Сг, откуда в силу положительности функции с(х) следует строгое диагональное преобладание. При этом благодаря удачному выбору искусственного коэффициента при второй производной реализуется точная компенсация линейного по шагу Ь± члена погрешности аппроксимации потока, чем и достигается второй порядок точности на равномерных и квазиравномерных [8] сетках,
В данной работе на основе анализа погрешностей схем прогнозируется поведение
к нулю и анализируются результаты численных экспериментов на последовательности сгущающихся сеток по расчету экспоненциальных слоев.
Расчеты проводились как на адаптивных сетках в физических переменных, так и в новых переменных на равномерных сетках — прообразах адаптивных сеток при координатных преобразованиях. Это позволило параллельно исследовать вопрос о различии в качестве расчетов в старых и новых координатах. Задача в новых переменных имеет вид, аналогичный исходному:
- (л(х(£, е), е)иц + а^, е)Щ + С1(£, е)и = , е), £е (0,1) и(0) = и<>, и(1) = иъ
где новые коэффициенты уравнения выражаются через старые в виде
х''
аг = а{х{С,£))^^ + ф{£,£),£), а = [4]2 c{x{C,£)), ¡х = [4]2 / {х{£,£)).
& £
Принципиальное отличие данной постановки состоит в том, что сетка по новой переменной £ равномерна, а производные решения по пей до определенного порядка равномерно ограничены. При построении разностных схем в новых переменных сложная производная второго порядка в преобразованном дифференциальном уравнении раскрывалась и объединялись подобные слагаемые с первыми производными. Например, схема (2) на равномерной сетке в новых переменных выглядела так:
-^{х{£,£),£)Лои + аг{^,£)Аои + сл{^,£)и = ¡г{£,е),
где До — оператор центральной разноети, а Л — трехточечный аналог оператора двойного дифференцирования по £, Аналогично в новых переменных строились две другие схемы. При расчетах в новых переменных производные х'^ъ х'^, входящие в выражения коэффициентов уравнения, аппроксимировались стандартными разностными отношениями со вторым порядком.
1. Координатное преобразование, генерирующее сетку
В основу алгоритма построения сетки кладется модификация [5] метода [9], гарантирующая ограниченность на ней производных до необходимого порядка (равного порядку производной в главном члене погрешности аппроксимации схемы). Для этого вблизи £ = 0 определяется бесконечно дифференцируемая функция х = ф{£, е) в виде
х
ф{С,£) = £к ({1 - В0-1/Л - 1) , 0 < £ < £о
(4)
где к — масштаб экспоненциального слоя; И = {1-£кВ)/£0 > 1+тг > 1, В = А/{1+пА)] при этом А — положительная константа, отделенная от нуля.
Для значений £0 < х < 1 функция (4) доопределяется с соблюдением гладкости класса Сг[0,1], С этой целью вычисляется необходимое число производных функции (4) в точке £0 и на их основе составляется полиномиальное разложение по формуле Тейлора, а затем функция, склеенная из (4) и полинома, растягивается или сжимается так, чтобы значение £ = 1 отображалось точно в ж = 1. Полученное отображение сгущает сетку вблизи х = 0 зависит от ряда пара метров х = Ф{^,£, А,к,1,п,^0) и имеет вид
Г Сгф{^,£), 0 < £ < ¿0,
С\
X = Ф{£,е,... ) = <
Ф{^0,£) + ф'{^0,е){^ - £0)+
+2 Ф' '{£,0,е){Ш - Ы2 + ... + +1 Ф{1){£0,е){Ш - + Ы£ - Ы^1
(5)
£0 < £ < 1,
где I < щ с0 > 0, например с0 = 1; константа сг > 0 выбирается, как сказано выше, из условия Ф(1,е,...) = 1. Число £0 в (5) означает долю длины отрезка, отображаемую в предполагаемую зону погранслоя, границей которого условно считается точка, где модуль производной решения по физической координате не больше некоторой кон-
етанты, не зависящей от малого параметра, а в точках внутри погранелоя производная стремится к бесконечности при малом параметре, стремящемся к нулю, В частности, при равномерной сетке по £ число £0 означает долю общего числа шагов сетки, попадающих в погранелой.
Данное преобразование (5) учитывает оценки старших производных и является обобщением на случай п > 2 отображения из работы [2], предназначенного исключительно для схем первого порядка точности. Параметр п означает порядок производной в главном члене погрешности аппроксимации схемы, так что базовое преобразование (5) является универсальным, т, е, обеспечивает любой класс гладкости и годится для схем произвольных порядков точности. Его универсальность проявляется также в пригодности для любых экспоненциальных слоев, для этого даже нет необходимости менять
А
при расчете степенных слоев при некотором ограничении в виде неравенства на пара-А
Отметим, что для численного решения задач с малым параметром часто применяются также разностные сетки Н.С. Бахвалова [1] и Г.И. Шишкина [3], однако они предназначены исключительно для экспоненциальных слоев, при этом зависят от значения функции а(х) в точке х = 0, что вынуждает корректировать параметры сеточного преобразования для каждой отдельной задачи. Это особенно затруднительно в двумерных задачах с переменной шириной пограничного слоя, где требуется предварительно пробными расчетами экспериментально выстраивать зависимость свободного параметра преобразования от продольной координаты переменного по ширине слоя. Предлагаемая нами формула (5) свободна от этого недостатка,
2. Оценка погрешности схем
На неравномерной сетке для аналога второй производной имеет место разложение
2 с1 V + сР
Л и = -(Д+ - Д-)и = ихх + -иххх + ^—ихххх + 0(Ь3), в 3 12
и слева:
5 = Ь+ + Ь-, d = Ь+ - Ь-, р= Ь+Ь-.
Следовательно, погрешность схемы (3) на достаточно гладких решениях исходного уравнения имеет вид
Ф(ж) = —^ (ихх(х) + ^иххх(д1)\ + а(х)их(х)-^ + Н \ 3 )
-нихх(х) + <ЩЬ±ихххШ + с(х)и - ¡(х), 6
где Ь± означает л ибо Ь- (при положительном а(х)).; либ о Ь+ (в противном случае), а ^ и — некоторые значения из про межутка (х - Ь-,х + Ь+), Используя исходное уравнение, после приведения подобных при второй производной получим
Н 2 йи2 п(т)Ь2
ф(я) =--~^ихх{х) - , ^ + -Ч^и^Ш. (6)
^ + Н 3(и + Н) 6
Из выражения (6) видно, что при любом фиксированном не малом параметре ^{х, е), когда производные решения имеют порядок константы, погрешность составляет величину О (К) па произвольной неравномерной сетке и вели чину порядка О (К2) па квазиравномерной сетке. Квазиравномерная сетка является, по определению, образом равномерной сетки при гладком строго монотонном преобразовании с ограниченной второй производной [8]; для таких сеток разность соседних шагов <1 является, очевидно, величиной второго порядка малости.
Однако при исследовании порядка аппроксимации следует учитывать и характер сетки, и возможное стремление ^ к нулю. Рассмотрим, например, поведение погрешности в экспоненциальном слое масштаба к = 1 на специальной адаптивной сетке, генерируемой координатным отображением (5), Очевидно, производные в слое имеют порядки
их = 0{1//л), ихх = 0{1/^2), иххх = 0{1//л3).
Шаги сетки в слое составляют величины К± = 0{^/М), а разность соседних шагов <1 = 0{^/М2), где N — число шагов сетки. Анализ порядков величин в выражении погрешности схемы (3) показывает, что в этом случае в слое Ф ~ (М2^)-1. Аналогичные рассуждения в отношении схем (1), (2) приводят к выражениям погрешностей
$ $ \о>{х) \'Р
Фф = -»-иххх{<21) + нихх{д2), ФИ = -»-иххх{<21) + у иххХ{Я2)
3 3 6
и к оценкам внутри слоя Ф ~ (М^,)-1 и Ф ~ (М2^,)-1 соответственно. Таким образом, погрешности схем (2), (3) обратно пропорциональны ^ (как и погрешность противопо-токовой схемы (1)), но они на порядок выше относительно шага сетки, А поскольку схема (3), в отличие от схемы (2), имеет еще и диагональное преобладание, это позволяет надеяться на ее равномерную сходимость со вторым порядком (по аналогии с результатом [2] для противопотоковой схемы).
Эта надежда рассеивается по причине, вовсе не связанной с поведением решения в слое, которая заключается в отрицательном влиянии компенсирующего множителя на оценку погрешности схемы (3) именно вне слоя, в зоне умеренных градиентов. Действительно, вне слоя шаг сетки К = 0(1/М), разность соседних шагов <1 = 0(1/М2), а все производные имеют порядок константы, поэтому схемы (1) и (2) естественным образом имеют погрешности соответственно 0(1/М) и 0(1/М2), В противоположность им из-за наличия искусственного множителя при второй производной погрешность схемы Булеева (3) имеет вне слоя оценку 0(1/М), что снижает порядок точности схемы до первого.
На основании проведенного анализа следует ожидать, что в задачах, где вязкость мала, но все изменения решения происходят исключительно в слое, как при простом ламинарном обтекании тела, когда вне слоя все производные близки к нулю, отмеченное выше нарушение оценки может и не иметь существенного значения. Если же вариация решения вне слоя значительна, а ^вместе с е остается малым вне погранелоя (например, при зависимости ^(х,е) = £т), то от схемы (3) не следует ожидать точности второго порядка, и источником ошибок при этом является первый порядок точности вне слоя.
Однако при некоторых специальных зависимостях ^(х, е), когда ^ мало только в слое, а вне слоя не стремится к нулю вместе с £, погрешность вне слоя оценивается как 0(1/М2), в таком случае следует ожидать более точных расчетов по схеме (3), Примером таких ситуаций может служить зависимость вида ^(х,е) = (х + е)т, при которой вязкость мала только в окрестности х = 0 а при х > 0 не стремится к пулю.
Оценим погрешность схемы (3) в новых переменных на равномерной сетке, В этом повторяя проведенные выше рассуждения, получим
ф(л- Н2 ТТ ЬУ п (П ) + «^Ш-ТТ (П ) Н = |а^)|Ь
= -%(0 - 12(^ + Н)^(П1) + ~(П2), Н = —Г~.
Отсюда видно, что погрешность составляет величину 0(Ь2/¡) при всюду малых значениях Если же у мало только в слое, то вне слоя Ф(£) = 0(Ь2),
3. Результаты численных экспериментов
Для моделирования описанных выше ситуаций в качестве тестов взяты краевые задачи вида
-уихх -их + и = 5^(4^), у = (е + вх)2, и(0) = 0, и(1) = 1,
где параметры вид принимают значения 0 и 1 в различных сочетаниях. Во всех вариантах экспоненциальный слой (масштаба к = 1 относнтельно у или, что то же, масштаба к = 2 относительно е) располагается в окрестности х = 0, Вариант 5=1 моделирует ситуацию, при которой вариация решения вне слоя значительна, а при д = 0, напротив, решение вне слоя близко к линейной функции. При 5 = 0 вязкость у = е2 мала во всей области (0 < х < 1), а при в = 1 вязкость мала только в окрестности нуля, в этих вариантах разнятся порядки точности схемы (3),
Параметры алгоритма построения сетки во всех расчетах были следующими: значе-о = 0. 5
метр п задаемся на единицу больше, т. е. п = I + 1, а свободный параметр А задавался равным двойке.
На всех рисунках за точные условно приняты решения, полученные по наиболее точной схеме (2) на сетке, самой детальной из использованных (с числом шагов N = 768 = 3х 28), При этом контролировалась сходимость решений, вычисленных по разным схемам, к общему пределу при сгущении сетки,
В таблицах N — число шагов сетки, а и р^ при ] = 1, 2, 3 суть апостериорная С
В задаче 1 значения параметров в = д = 0, Результаты расчетов по трем схемам приведены на рис, 1, Видно, что на грубой сетке схема (2), не обладая диагональным преобладанием, осциллирует. При детализации сетки осцилляции затухают и схема дает результаты, превосходящие по точности расчеты по схемам (1) и (3), В табл. 1 приведены результаты расчета задачи 1 на последовательности сгущающихся сеток. Очевидно, малость параметра вязкости во всей области действительно не позволяет по схеме Бу-леева (3) реально получить второй порядок точности; результаты расчета по ней даже несколько хуже, чем по схеме первого порядка точности, В то же время схема (2) демонстрирует второй порядок точности и существенно более хорошие результаты,
В задаче 2 значения параметров э = 0 1 = 1 т-е- вязкость по-прежнему мала всюду в области, а периодическая правая часть генерирует вне слоя более сложную картину решения с существенной вариацией. Результаты расчета задачи 2 приведены на рис, 2 и в табл. 2, Соотношение результатов по различным схемам здесь примерно такое же, как
Рис. 1. Численные решения задачи 1 при количестве шагов сетки N = 24 (слева) и 48 (справа) Fig. 1. Numerical solutions of problem 1 with the number of steps N = 24 (left and N = 48 (right)
Таблица 1. Оценки ошибки и порядка точности, полученные на последовательности сеток при решении задачи 1 по трем схемам. Вязкость у = е2 = 10-4
Table 1. Estimates of the error and order of accuracy obtained on a sequence of grids by solving problem 1 using three schemes. Viscosity ^ = e2 = 10-4
N ¿1 Pi ¿2 P2 ¿3 'Рз
48 1.02E-02 6.62 1.05E+00 -0.07 1.32E-02 6.24
96 5.35E-03 0.93 6.18E-03 7.42 7.06E-03 0.90
192 2.74E-03 0.97 4.51E-04 3.78 3.65E-03 0.95
384 1.39E-03 0.98 1.13E-04 2.00 1.86E-03 0.98
768 6.97E-04 0.99 2.83E-05 2.00 9.36E-04 0.99
Рис. 2. Численные решения задачи 2 при количестве шагов сетки N = 48 (слева) и 96 (справа) Fig. 2. Numerical solutions of problem 2 with the number of steps N = 48 (left) and N = 96 (right)
Таблица 2. Оценки ошибки и порядка точности, полученные на последовательности сеток при решении задачи 2 по трем схемам. Вязкость ^ = е2 = 10-4
Table 2. Estimates of the error and order of accuracy obtained on a sequence of grids by solving problem 2 using three schemes. Viscosity ^ = e2 = 10-4
N ¿1 Pi ¿2 Р2 ¿3 Рз
48 2.52Е-01 1.99 2.05Е+00 -1.03 3.19Е-01 1.81
96 1.07Е-01 1.24 2.12Е-01 3.27 1.56Е-01 1.03
192 4.99Е-02 1.10 4.27Е-02 2.31 7.05Е-02 1.15
384 2.40Е-02 1.05 1.01Е-02 2.07 З.ЗОЕ—02 1.09
768 1.18Е-02 1.03 2.51Е-03 2.01 1.59Е-02 1.05
Таблица 3. Оценки ошибки и порядка точности, полученные на последовательности сеток при решении задачи 3 по трем схемам. Вязкость ^ = (е + ж)2, е = 0.01
Table 3. Estimates of the error and order of accuracy obtained on a sequence of grids by solving problem 3 using three schemes. Viscosity ^ = (e + x)2, e = 0.01
N ¿1 Pi ¿2 Р2 ¿3 Рз
48 1.59Е-02 5.97 1.10Е-01 3.18 2.58Е-03 8.60
96 8.15Е—03 0.97 1.80Е-03 5.94 8.04Е-04 1.68
192 4.11Е-03 0.99 4.55Е-04 1.99 2.38Е-04 1.75
384 2.06Е-03 0.99 1.14Е-04 2.00 6.76Е-05 1.82
768 1.03Е-03 1.00 2.85Е-05 2.00 1.84Е-05 1.88
Таблица 4. Оценки ошибки и порядка точности, полученные на последовательности сеток при решении задачи 4 по трем схемам. Вязкость у = (е + ж)2, е = 0.01
Table 4. Estimates of the error and order of accuracy obtained on a sequence of grids by solving problem 4 using three schemes. Viscosity ^ = (e + ж)2, e = 0.01
N ¿1 Pi ¿2 Р2 ¿3 Рз
48 1.19Е-01 3.15 2.54Е-01 2.09 7.91Е-02 3.66
96 6.63Е-02 0.84 1.80Е-02 3.82 2.54Е-02 1.64
192 3.53Е-02 0.91 4.42Е-03 2.03 7.32Е-03 1.79
384 1.82Е-02 0.95 1.10Е-03 2.01 1.97Е-03 1.89
768 9.26Е-03 0.98 2.75Е-04 2.00 5.12Е-04 1.95
в предыдущей задаче: противопотоковая схема (1) даже точнее схемы (3), а схема без диагонального преобладания (2) на совсем грубых сетках уступает обеим схемам в точности, но по мере детализации сетки восполняет отставание и демонстрирует второй практически наблюдаемый порядок точности и существенно более точные результаты.
В задаче 3 значения параметров э = 1, д = 0, т. е. вязкость мала лишь в окрестности нуля, а правая часть нулевая. Результаты расчета задачи 3 приведены в табл. 3. В данной задаче у схемы Булеева нет проблем с погрешностью вне слоя, апостериорная оценка порядка точности по ней, хотя и медленно, но стремится к двойке, и результаты ранжируются в соответствии с порядком точности — обе схемы второго порядка дают примерно одинаковые результаты, превосходящие по точности результаты расчета по схеме первого порядка.
Рис. 3. Численные решения задачи 4 при количестве шагов сетки N = 24 (слева) и 48 (справа) Fig. 3. Numerical solutions of problem 4 with the number of steps N = 24 (left) and N = 48 (right)
irl
и
-0.5
Расчет в новых переменных по схеме (2) M
Точное = £2 = 10-8 f
-—N= 24
N= 48
— N=96 S J ! W F i /
v\ / /
и
X
1 о
0.5
X 1
0.5
Рис. 4. Сравнение расчетов по схеме (2) в исходных физических (слева) и в новых (справа) переменных
Fig. 4. Comparison of calculations according by scheme (2) in original physical (left) and new (right) variables
В задаче 4 параметры q = э = 1, т. е .у, = 0(е2) является малой лишь в окрестности нуля, вне слоя ^ = 0(1), а вариация решения вне слоя значительна. Результаты расчетов приведены па рис. 3 и в табл. 4. Здесь ошибки по схемам второго порядка также примерно одинаковы. На грубой сетке схема (2), как и в прежних вариантах, проявляет осцилляции, исчезающие при детализации сетки.
Заметим, что во всех задачах осцилляции решения схемы (2) асимптотически затухают при увеличении числа шагов сетки, по порог детализации сетки, при котором результаты можно считать удовлетворительными, отодвигается все дальше по мере уменьшения вязкости. Этот результат свидетельствует о том, что в физических переменных схема (2) пе сходится равномерно по малому параметру.
Дня сравнения помимо расчетов па адаптивных сетках в физических переменных проводились также расчеты па соответствующих равномерных сетках в новых переменных. Результаты решения задачи 2 дня схемы (2) при очень малой вязкости
Рис. 5. Сравнение расчетов в исходных физических и новых неременных, полученных на эквивалентных сетках при ß(x, е) = (х + е)2 (вверху) и при ß(x, е) = £2 (внизу) Fig. 5. Comparison of calculations in the original physical and new variables obtained on equivalent grids for ß(x, e) = (x + e)2 (top) and ß(x, e) = e2 (bottom)
ß = e2 = 10-8 приведены на рис. 4. Они ясно свидетельствуют о предпочтительности использования повой переменной при расчетах по пей. Действительно, при расчете в физических координатах удовлетворительное совпадение с точным решением достигается лишь на сетке с числом шагов N = 192, хотя и здесь еще заметны слабые осцилляции, а па более грубых сетках решение совершенно подавляется осцилляциями, В то же время в новых переменных осцилляций пет даже па самых грубых сетках. Заметим, что преимущество расчетов в новых переменных наблюдалось только для схемы (2), не имеющей диагонального преобладания. Для двух других исследованных схем (имеющих диагональное преобладание) результаты расчетов в исходных переменных па адаптивной сетке и в новых переменных па равномерной сетке приблизительно эквивалентны по точности. Это иллюстрирует рис. 5, из которого видно, что лишь дня схемы (2) выбор переменных имеет существенное значение дня точности.
Заключение
Причиной неудовлетворительных результатов расчета в физических переменных по схеме (2) при малой вязкости является плохая обусловленность соответствующей си-
етемы уравнений, вызванная отсутствием диагонального преобладания на грубой сетке при больших градиентах решения в слое. При детализации сетки нарушение диагонального преобладания становится менее значительным, асимптотически исчезая при N ис улучшением обусловленности системы осцилляции затухают, В новых пе-
ременных, устраняющих слои, производные равномерно ограничены всюду, шаг сетки постоянный, поэтому для плохой обусловленности нет таких причин, как в физических переменных.
Схемы (1) и (3), имея диагональное преобладание, в любых координатах и при любых шагах порождают хорошо обусловленные системы, поэтому расчеты по ним одинаково удовлетворительны в исходных физических и новых переменных.
На основании проведенного исследования подтверждается достаточно надежная работа схем с диагональным преобладанием, при этом схема Булеева (3), имея теоретически второй порядок точности, в задачах с малой вязкостью всюду в области на практике обнаруживает только первый порядок и дает результаты, сравнимые по точности с расчетами по противопотоковой схеме (1) первого порядка.
Схема (2) на сгущающихся сетках всегда показывает второй практически наблюдаемый порядок точности, но сходимость по ней в физических координатах не является равномерной по малому параметру, ее преимущество в точности перед другими схемами проявляется только при достаточной детализации сетки, когда нарушение диагонального преобладания слабое или вовсе отсутствует.
Расчеты показывают, что схема (2), не имеющая диагонального преобладания, обеспечивает более стабильные результаты на равномерной сетке по новой переменной, относительно которой производные решения во всей области становятся умеренными, чем на ее образе — адаптивной сетке по физической переменной. Для двух других исследованных схем (с диагональным преобладанием), напротив, практически безразлично, в исходных физических или преобразованных переменных вести расчеты.
Благодарности. Исследование выполнено при финансовой поддержке РФФИ в рамках научного гранта № 20-01-00231,
Список литературы
[1] Бахвалов Н.С. К оптимизации методов решения краевых задач при наличии пограничного слоя. Журнал вычислительной математики и математической физики. 1969; 9(4):841-859.
[2] Лисейкин В.Д. О численном решении уравнений со степенным погранслоем. Журнал вычислительной математики и математической физики. 1986; 26(12):1813—1820.
[3] Miller J.J.K., O'Riordan Е., Shishkin G.I. Finited numerical methods for singular perturbation problems. Singapure, New Jersey, London, Hong Kong: World Scientific; 2012: 191.
[4] Liseikin V.D. Layer resolving grids and transformations for singular perturbation problems. Utrecht; Boston: VSP; 2001: 284.
[5] Лисейкин В.Д., Паасонен В.И. Компактные разностные схемы и адаптивные сетки для численного моделирования задач с пограничными и внутренними слоями. Сибирский журнал вычислительной и прикладной математики. 2019; 22(1):41-56.
[6] Ильин A.M. Разностная схема для дифференциального уравнения с малым параметром при старшей производной. Математические заметки. 1969; 6(2):237-248.
[7] Булеев Н.Н. О численном решении двумерных уравнений эллиптического типа. Численные методы механики сплошной среды. Новосибирск, СО АН СССР, 1975; 6(3):8-28.
[8] Калиткин Н.Н., Альшин А.В., Альшина Е.А., Рогов В.Б. Вычисления на квазиравномерных сетках. М.: Наука, Физматлит, 2005; 224.
[9] Liseikin V.D. Grid generation for problems with boundary and interior layers. Novosibirsk: NGU; 2018: 296. ISBN:978-5-4437-0822-5.
Вычислительные технологии, 2020, том 25, № 5, с. 66-79. © ФИЦ ИВТ, 2020 ISSN 1560-7534
Computational Technologies, 2020, vol. 25, no. 5, pp. 66-79. © FRC ICT, 2020 elSSN 2313-691X
COMPUTATIONAL TECHNOLOGIES
DOI: 10.25743/ICT.2020.25.5.006
Convergence behavior of popular schemes in case of calculating on adaptive grids problems with layers
Liseikin Vladimir D.1'*, Paasonen Vict or 1.1,2
1 Federal Research Center for Information and Computational Technologies, Novosibirsk, Russia
2
*
Received February 28, 2020, revised August 2, 2020, accepted August 28, 2020
Abstract
The paper compares solution quality to some model second-order equation with a small parameter obtained through three different schemes both on special adaptive grids specified explicitly by coordinate transformations eliminating layers and on uniform grids in a new coordinate related to the transformations. The schemes up to second order in physical and transformation variables both with a diagonal and not diagonal dominance and the simplest counter-flow scheme are analyzed. Predictions of a solution behavior based on estimates of solution errors are described, which are confirmed by numerical experiments and proofs. It is established, in particular, that the scheme of the second order with a diagonal dominance converges uniformly if the coefficient before the second derivative is small at the points of the boundary layer only. It was also demonstrated for the schemes without a diagonal dominance, mach better solutions without oscillations are obtained on uniform grids in new variables than on corresponding adaptive grids in the original physical coordinates.
Keywords: uniform convergence, adaptive grid, boundary layer, diagonal dominance, small parameter.
Citation: Liseikin V.D., Paasonen V.I. Convergence behavior of popular schemes in case of calculating on adaptive grids problems with layers. Computational Technologies. 2020; 25(5):66-79. DOI: 10.25743/ICT.2020.25.5.006. (In Russ.)
Acknowledgments. This research was supported by RFBR (grant No. 20-01-00231). References
1. Bakhvalov N.S. On the optimization of the methods for solving boundary value problems in the presence of a boundary layer. USSR Computational Mathematics and Mathematical Physics. 1969; 9(4):139-166.
2. Liseikin V.D. Numerical solution of equations with a power boundary layer. Computational Mathematics and Mathematical Physics. 1986; 26(6): 133-139.
3. Miller J.J.K., O'Riordan E., Shishkin G.I. Finited numerical methods for singular perturbation problems. Singapure, New Jersey, London, Hong Kong: World Scientific; 2012: 191.
4. Liseikin V.D. Layer resolving grids and transformations for singular perturbation problems. Utrecht; Boston: VSP; 2001: 284.
5. Liseikin V.D., Paasonen V.I. Compact difference schemes and layer-resolving grids for numerical modeling of problems with boundary and interior layers. Numerical Analysis and Applications. 2019; 12(l):37-50.
6. Il'in A.M. Differencing scheme for a differential equation with a small parameter affecting the highest derivative. Mathematical Notes. 1969; 6(2):596-602.
7. Buleev N.N. On numerical solution of two-dimensional equations of elliptic type. Chislennye Metody Mekhaniki Sploshnoy Sredy. 1975; 6(3):8-28. (In Russ.).
8. Kalitkin N.N., Alshin A.B., Alshina E.A., Rogov V.B. Computation on kvaziuniform grids. Moskow: Nauka, Physmathlit, 2005; 224. (In Russ.).
9. Liseikin V.D. Grid generation for problems with boundary and interior layers. Novosibirsk: NGU; 2018: 296. ISBN:978-5-4437-0822-5.