CC BY
34
Вычислительные технологии
Том 15, № 1, 2010
Численное исследование независимой аппроксимации граничных условий на решениях с разрывами производных*
Д. А. Ичетовкин Новосибирский государственный университет, Россия e-mail: neo@211.ru
В. И. Паасонен
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: paas@ict.nsc.ru
Работа является продолжением цикла статей, посвященных численному исследованию разностных схем для параболических и эллиптических уравнений с разрывными коэффициентами, в которых граничные условия аппроксимируются с произвольным порядком непосредственно, без привлечения дифференциального уравнения и его продолженной системы. Вопросы теоретического обоснования метода и его возможные приложения ранее обсуждались на страницах журнала "Вычислительные технологии", а предметом настоящей работы является подтверждение порядков сходимости в численных экспериментах на решениях с разрывами первых производных.
Ключевые слова: односторонняя аппроксимация потока, многоточечные граничные условия, многоточечная аппроксимация потока, компактная схема, неоднородная область, аппроксимации граничных условий, схема высокого порядка точности, высокоточная схема.
Введение
Известно, что многие дифференциальные уравнения в частных производных второго порядка, в том числе и с переменными коэффициентами, в ортогональных координатных системах могут быть аппроксимированы с четвертым порядком относительно шагов пространственной сетки на компактном шаблоне, не выступающем за пределы 3 х 3 х ... х 3-точечного шаблона (см., например, [1, 2]). Для применения высокоточных схем и в задачах с разрывными коэффициентами необходимо решить вопрос об аппроксимации с адекватной точностью внутренних граничных условий. Среди существующих подходов для решения этой проблемы отметим методы [3, 4]. Первый базируется на требовании гладкости потока в точках разрыва производной, и его обобщение на двумерный случай проблематично, второй основан на разностных аппроксимациях закона сохранения в балансных ячейках.
*Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00264) и интеграционного проекта СО РАН № 103.
© ИВТ СО РАН, 2010.
Авторы настоящей статьи считают, что с точки зрения универсальности подхода (относительно характера дифференциальных уравнений, порядка точности, типа граничных условий) для построения простого, гибкого, удобного для пользователя и хорошо структурированного алгоритма целесообразно в граничных условиях непосредственно использовать односторонние разностные аналоги первых производных. При этом необходимая точность в граничных условиях может быть достигнута выбором достаточного числа точек в односторонней аппроксимации потока, В общем виде такой подход сформулирован и исследован в [5], а в [6] предложена эквивалентная параллельная технология для его реализации. Метод устойчив и универсален, пригоден для условий Неймана, для условий третьего рода и условий равенства потоков на границах раздела различных сред.
Представленный здесь материал является продолжением данного цикла работ и посвящен численному исследованию алгоритмов такого типа на одномерных и двумерных задачах теплопроводности с разрывными коэффициентами. Необходимость численного эксперимента вызвана тем, что иногда наличие разрывов в коэффициентах не позволяет получить в расчетах ошибку ожидаемого порядка, вычисленного из теории, которая обычно предполагает достаточную гладкость решения. Предлагаемые здесь численные исследования были проведены с целью подтверждения теоретических результатов путем изучения реально наблюдаемых порядков сходимости.
1. Постановка одномерной разностной задачи
Пусть х0, х1;..., хг — возрастающая последовательность значений пространственной координаты х, представляющих собой координаты границ однородных слоев х^-1 <х<х^-, = 1, ...,г) с постоянными теплофизичеекими характеристиками. Внутри каждого слоя решение и(х, £) удовлетворяет уравнению теплопроводности на интервале времени 0 < £ < Т
ди д2и „ .
На внешних границах х = хо, х = хг поставлены граничные условия первого, второго или третьего рода (объединенные общей формулой)
ди
а на внутренних границах х = х1; ...,хг-1 заданы условия равенства потоков
ди ди
где индексы + и — означают пределы справа и слева.
Во внутренних узлах слоев уравнение теплопроводности аппроксимируем либо обычной чисто неявной схемой точности 0(т + к2), либо компактной схемой точности 0(т2 + к4):
„ и+1 — „.и к2
сА = АА ип- 5Г+1/2, А = Е + а Ь2 А, Б = Е + ^ А,
где Е — тождественный оператор, Л — обычная аппроксимация двойного дифференцирования, а вес с целью повышения порядка точности выражается следующим образом:
1 а Ат а =---, а = ——■
12 2' сЬ2
Таким образом, во внутренних узлах любого однородного слоя разностное уравнение приводится к симметричным трехточечным (временной индекс опущен) соотношениям
аи— + Ьпг + аПг+1 =
где ^ — результат операций над значениями функции и на нижнем слое и правой части f уравнения теплопроводности.
Так как теплофизичеекие характеристики и шаг сетки в пределах слоя постоянны, то с, А и Н, а также коэффициенты разностного уравнения следовало бы снабдить индексом, означающим номер слоя, которому принадлежит текущий узел сетки. Здесь и ниже будем его опускать, подразумевая зависимость этих величин от номера слоя.
Первые производные, входящие во внешние и внутренние граничные условия, аппроксимируем, применяя односторонние разделенные разности на (в + 1)-точечном шаблоне, Разность "вперед" имеет вид
^£ = 1 V- ^ к
к=0
актн.
Здесь Т/ — оператор едвига, аак- коэффициенты разностного оператора. Максимально возможный порядок аппроксимации равен при этом
Т
ак = -—--Ск8, к ф 0, а0 = - У
к=1
где С*к — число сочетайий из в по к.
Аналогично выглядит разность "назад" порядка в:
Д- 1 А 1 '
Н —Н^ к Н
к=0 к=0
Иначе говоря, коэффициенты противоположно ориентированных односторонних высокоточных разностей, как и в случае простейших разностей первого порядка, различаются только знаком,
В узлах, совпадающих с границами раздела сред, разностное уравнение запишем с помощью (в + 1)-точечных односторонних разностей:
. Дй . Д-й
А+—щ - А-—щ = О,
Н+
где индексы + и — означают принадлежность к правому или левому слою относительно данной границы раздела сред. На внешних границах также воспользуемся односторонними аналогами первых производных порядка 0(Нй):
щ—щ + ЦоПо = 00, —им - цгим = -фг. Н+
При сравнении методов различного порядка условимся использовать двух- и трехточечные граничные условия в сочетании с чисто неявной схемой, а в сочетании с компактной схемой — четырех- и пятиточечные граничные условия. Таким образом, при естественном предположении т = 0(к2) получим разностные схемы первого, второго,
к
односторонних разностей на различных шаблонах.
2. Многомерные разностные краевые задачи
В двумерном случае для простоты будем рассматривать только задачу Дирихле для уравнения теплопроводности в области, составленной из полос, занятых материалами с различными постоянными характеристиками. При кусочно-постоянных коэффициентах разностное уравнение, аппроксимирующее с погрешностью 0(т2 + к4) уравнение теплопроводности во внутренних узлах подобластей, получается факторизацией схемы с весами при специальных значениях весов осреднения
- _ ак 12 2 '
Лт и
ак = щ> к
1, 2.
На границах раздела сред использованы одномерные разностные условия равенства потоков порядка 5, описанные выше. При этом, как и в одномерном случае, рассмотрим варианты 5 = 1 и 5 = 2 для схемы второго порядка из = 3и 5 = 4 для компактной схемы четвертого порядка с фиксацией естественного предельного соотношения между тк
порядка 0(к5) для зпачений 5 < 4, Сравнение проводилось также с расчетами по схеме сквозного счета (СС) приближенной факторизации.
Для схемы в дробных шагах по линиям сетки с учетом разностных граничных условий на каждом временном шаге стандартным способом формируются одномерные разностные задачи такого же типа, что и в одномерном случае. Каждая из одномерных систем линейных алгебраических уравнений, к которым сводится общая задача, имеет почти трехдиагональную структуру. Ниже в качестве примера приведены два фрагмента расширенной матрицы системы для частного случая 5 = 4, Для компактности записи здесь использованы обозначения
в
^0 -1
к1
-ан
1,
А_
/Г
-а
31
11
и,
-а,-.
Первый фрагмент соответствует окрестности левой границы многослойного пакета:
^0 + во в1 в2
а+ а+
0 а+
0 0 а+
0 0 0
0 0
0 0
а+ а+ 0
—Фо
^4
0 0 0 0
0
Второй фрагмент иллюстрирует структуру матрицы в окрестности любой границы раздела сред:
а_ Ь- а- 0 0
0 а- Ь- а- 0
0 а- Ь- а-
0 0 а- Ь-
0 7- 7з" 7- 7-
0 0 0 0
0 0 0 0
0 0 0 0
7о
0 0 0 0 0 Я-з
0 0 0 0 0 Я-2
а- 0 0 0 0
+ 7+ 7+ 72+ 7+ 7+ 0 0
а+ Ь+ а+ 0 0
0 а+ Ь+ а+ 0
0 0 а+ Ь+ а+ 0 -^г+3
0 0 а+ Ь+ а+ Кг+4
К
1-4
Если разбивать каждый слой не грубее, чем на 5 шагов, то "длинные" условия на соседних границах не будут "переплетаться" между собой, что позволяет с помощью стандартных операций метода Гаусса понизить их длину, вычитая из длинных строк комбинации соседних трехточечных строк, В результате таких операций матрица системы превращается в трехдиагональную, после чего можно воспользоваться обычным методом прогонки. Именно так проводилась программная реализация, хотя существует и другой эквивалентный параллельный алгоритм [6],
3. Результаты численных экспериментов
С целью исследования практически наблюдаемых порядков сходимости схем на предмет сравнения их с теоретическими были проведены численные эксперименты. Сравнивались результаты, полученные по четырем различным схемам с односторонними многоточечными разностями в граничных условиях — с первого до четвертого порядка аппроксимации. Эти схемы обозначаются в описании тестов цифрами, соответствующими порядку точности, С целью достижения жестких условий тестирования шаги в слоях были взяты равными, а не подбирались из соображения равенства тепловых сопротивлений слоев, и использовалась двумерная квадратная сетка. Пространственные и поверхностные источники тепла отсутствовали.
При детализации сетки начиная с начального шага к0 = 0.1 наблюдалось стремление численных решений, полученных по различным схемам, к одному пределу, при этом, как и ожидалось, сходимость была быстрее для схем более высокого порядка. Поэтому при сравнении результатов на сетках с шагами к и к/2 за точное решение условно принималось численное решение, полученное по схеме четвертого порядка аппроксимации па сетке с шагом к/4. Всюду ниже термин "точное решение" имеет именно этот смысл. В одномерной задаче для контроля проводилось до восьми дроблений шага вдвое. В двумерной задаче такая детализация сетки была бы обременительна, особенно если учесть, что в этом случае для сохранения постоянства отношения т/к2 шаг по времени каждый раз уменьшался бы в четыре раза. Оценка практически наблюдаемого порядка точности проводилась на сгущающихся сетках по коэффициенту убывания ошибки
Гп= ъ = Чиь-и\\с ; который теоретически при малых к должен стремиться к степепям двойки 25.
В одномерном случае рассматривался двухслойный пакет единичной толщины с границей раздела сред посередине. Теплоемкость в обоих слоях принималась равной единице, а коэффициенты теплопроводности различались в сто раз (А_ = 2 и А+ = 0.02 соответственно для левого и правого слоев).
На рис. 1 показаны результаты расчетов на момент времени £ = 0.8, полученные на различных сетках (с шагом Н = 0.1 и 0.05) по схемам первого и четвертого порядка аппроксимации для смешанной задачи с начальным распределением температур
1Г(х,0) = ^ ехр(—ж2)(1 - 100ж) и внешними граничными условиями
1 д~и
и{ 0, г) = -(1 + 8ш(4тг<)), ^(1, г) = 49/е.
Из рис. 1 видно, что даже на очень грубой сетке расчеты по схеме 4 весьма близки к точному решению. Расчеты по схемам второго и третьего порядка точности занимали промежуточное положение между приведенными крайними результатами и, чтобы не загромождать рисунок, их графики не представлены.
В таблице приведены нормы ошибок на сгущающихся сетках для тестируемых схем и их отношения. Для всех порядков обнаруживается стремление отношений к 2й, а это свидетельствует о том, что несмотря на наличие существенных изломов в точном решении теоретический порядок точности на практике подтверждается.
Таблица и рис. 2 иллюстрируют сходимость численных решений на сгущающихся сетках для задачи Дирихле в двумерном случае. Область представляет собой единичный квадрат, разделенный по вертикали на два равных прямоугольника, занятых материалами с теми же теплофизическими характеристиками, что и в одномерном варианте. Условия на границе и начальные данные задавались из выражения
и{х,у,г) = \{ 1 + вш(4^)) (¡^т + ^т)-
Заметим, что точным решением эта функция не является, она используется лишь для генерирования согласованных начальных и граничных условий. Сплошной линией
Рис. 1. Численное решение одномерной задачи при пяти (а) и десяти (б) шагах в слое; сплошная
линия — точное решение, квадратики — результаты расчета по схеме 4, Ч--результаты расчета
по схеме 1
Экспериментальная оценка порядка точности в одномерной и двумерной задаче
Схема расчета Ошибка 6/г Ошибка £/¿/2 Коэффициент ть
Одномерная задача
1 /г = 1/20 0.04099 0.01881 2.17931
2 0.03648 0.00903 4.03675
3 0.01207 0.00212 5.68965
4 0.00742 0.00061 12.0609
1 /г = 1/40 0.02684 0.01338 2.00677
2 0.02132 0.00503 4.23852
3 0.00595 0.00079 7.53160
4 0.00215 0.00013 16.0440
Дщр черная, задача
СС 0.077821 0.030338 2.565067
1 /г = 1/20 0.142020 0.049197 2.886750
2 0.077260 0.021209 3.642719
3 0.063559 0.010129 6.274994
4 0.042807 0.003269 13.09257
СС 0.071975 0.024872 2.893750
1 /г = 1/40 0.148847 0.051145 2.910278
2 0.057868 0.013118 4.411326
3 0.027224 0.003122 8.720084
4 0.012520 0.000776 16.13031
Рис. 2. Численное решение двумерной задачи при пяти (а) и десяти (б) шагах в слое
изображены изолинии точного решения на момент времени £ = 1, а маркеры имеют тот же смысл, что и в одномерном случае. Здесь также решение по схеме 4 существенно ближе к точному, чем по схеме 1.
Результаты исследования порядка сходимости дня данного двумерного теста приведены в таблице. Здесь сравнение проводилось также со схемой сквозного счета. Видно стремление отношения погрешностей к целым степеням двойки с показателем, равным порядку аппроксимации, что согласуется с теоретическими результатами.
Полученные результаты свидетельствуют о том, что схемы с односторонней аппроксимацией потоков в граничных условиях в кусочно-однородных областях имеют порядки точности, совпадающие с теоретическими. Это объясняется тем, что предложенный алгоритм основан на непосредственной аппроксимации потоков односторонними разностными аналогами без привлечения самого дифференциального уравнения или его продолженной системы, в силу чего полная гладкость решения фактически не требуется, необходима лишь его достаточная гладкость в пределах каждого однородного включения отдельно. По этой же причине в отличие от методов, в которых разностные аппроксимации граничных условий зависят от дифференциального уравнения, в рассматриваемом случае не требуется гладкость потока на границах раздела сред.
Список литературы
[1] Паасонен В.И. Компактные схемы для уравнений второго порядка с конвективными членами // Вычисл. технологии. 1998. Т. 3, № 1. С. 55-66.
[2] Paasonen V.l. Compact schemes for system of second-order equations without mixed derivatives // Rus. J. Numer. Analvs. and Math. Model. 1998. Vol. 13, No. 4. P. 335-344.
[3] Коскин П.И. Схема повышенной точности для уравнения теплопроводности с разрывными коэффициентами // Численные методы механики сплошной среды. 1979. Т. 10, № 2. С. 85-96.
[4] Ильин В.П. Балансные аппроксимации повышенной точности для уравнения Пуассона // Сибирский мат. журн. 1996. Т. 37, № 1. С. 151-169.
[5] Paasonen V.l. Compact difference schemes for inhomogeneous boundary value problems // Rus. J. Numer. Analvs. and Math. Model. 2004. Vol. 19, No. 1. P. 65-81.
[6] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. Т. 8, № 3. С. 98-106.
Поступила в редакцию 22 декабря 2008 г., в переработанном виде — 8 августа 2009 г.
