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

Вычислительные технологии

Том 15, № 1, 2010

Численное исследование независимой аппроксимации граничных условий на решениях с разрывами производных*

Д. А. Ичетовкин Новосибирский государственный университет, Россия e-mail: [email protected]

В. И. Паасонен

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

e-mail: [email protected]

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

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

Введение

Известно, что многие дифференциальные уравнения в частных производных второго порядка, в том числе и с переменными коэффициентами, в ортогональных координатных системах могут быть аппроксимированы с четвертым порядком относительно шагов пространственной сетки на компактном шаблоне, не выступающем за пределы 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 г.