CC BY
29
Вычислительные технологии
Том 15, № 6, 2010
Численное исследование высокоточных схем в областях клетчатой структуры*
Д. А. Ичетовкин Новосибирский государственный университет, Россия e-mail: neo@211.ru
В. И. Паасонен
Институт вычислительных технологий СО РАН, Новосибирск, Россия Новосибирский государственный университет e-mail: paas@ict.nsc.ru
Работа является продолжением цикла статей, рассматривающих метод расчета краевых задач в кусочно-однородных областях с независимой высокоточной аппроксимацией потоков в граничных условиях. Ранее этот метод исследован авторами в случае задачи Дирихле для уравнения теплопроводности в слоистых областях. В данной работе представлено численное исследование метода для смешанных краевых задач в кусочно-однородных областях более общей, чем слоистые, структуры.
Ключевые слова: односторонняя аппроксимация потока, многоточечные граничные условия, компактная схема, неоднородная область, аппроксимации граничных условий, высокоточная схема, схема высокого порядка точности, многоточечная аппроксимация потока.
Введение
В [1] был предложен универсальный метод решения краевых задач в составных областях, основанный на аппроксимации граничных условий с любым порядком точности, охватывающий как традиционные схемы второго порядка, так и компактные высокоточные схемы. Теоретические оценки сходимости были подтверждены [2] численно на одномерных и двумерных задачах Дирихле для уравнения теплопроводности в областях слоистой структуры. Однако в случае неоднородных областей более общего вида, составленных из конечного числа однородных одинаково ориентированных прямоугольников, а также для более общих граничных условий при реализации разностной краевой задачи возникает проблема, связанная с наличием особых узлов. Это уступы подобласти, точки сопряжения трех или четырех материалов, а на внешней границе — точки сопряжения двух материалов и изломы границы. Для таких узлов в методике [1] непосредственно не предусмотрен способ вычисления решения, а доопределение решения в особых узлах применением различных экетраполяционных процедур может приводить к потере устойчивости и точности. Кроме того, такое доопределение решения в силу
*Работа выполнена при поддержке РФФИ (грант № 08-01-00264) и Интеграционного проекта СО РАН № 103.
неоднозначности способа и направления экстраполяции делает метод в целом несколько неопределенным,
В [3] была разработана концепция пересчета решения в особых узлах с любым порядком, достаточно перспективная ввиду того, что алгоритм пересчета не опирается на экстраполяцию, а основан на интегрировании по четырем соседним ячейкам разностных аппроксимаций дифференциального уравнения. Кроме того, процедура в этом случае является однозначной и универсальной, т, е, единой для любого типа особых точек. Настоящая работа посвящена численному исследованию разностных схем различных порядков аппроксимации для уравнения теплопроводности в составных областях описанного выше типа на основе алгоритма [1] с привлечением формул замыкания [3] для решения в особых точках области. Цель состоит в тестировании формул замыкания алгоритма и в экспериментальном подтверждении соответствия практически наблюдаемых и теоретических порядков сходимости в задачах с разрывными коэффициентами,
1. Концепция
Пусть двумерная неоднородная область представляет собой объединение однородных подобластей, составленных из конечного числа прямоугольников, одинаково ориентированных по осям, с постоянными в пределах подобласти теплофизическими характеристиками, Внутри каждой из однородных подобластей выполняется уравнение теплопроводности с постоянными коэффициентами, на различных участках внешних границ поставлены условия первого, второго или третьего рода, а на границах раздела сред — условия непрерывности температуры и нормальных к границе потоков. Чтобы смешанная краевая задача была корректна, граничные условия предполагаются согласованными.
Продолжив все границы раздела сред и все внешние границы в вертикальном и горизонтальном направлениях до пересечений со всеми границами ортогонального направления, переведем логически сложную неоднородную область в более простую область клетчатой структуры, где каждая клетка однородна, причем некоторые клетки могут быть и пустыми. Ясно, что при продолжении границ возникают фиктивные границы раздела сред между клетками, занятыми одним и тем же материалом, В этом случае с целью унификации алгоритма будем формально рассматривать их как различные материалы, т, е, будем задавать на фиктивных границах условия непрерывности потока, которые ввиду совпадения коэффициентов теплопроводности означают в данном случае условие гладкости в форме равенства левых и правых производных. Построим двумерную прямоугольную сетку, согласованную с внешними и внутренними границами так, чтобы в пределах каждой клетки сетка была равномерной в обоих направлениях, Добавим слева и справа, сверху и снизу по одной фиктивной линии сетки, т, е, окаймим область со всех сторон слоем шириной в одну ячейку. Это позволяет рассматривать внутренние и граничные узлы всех типов единообразно как центральные узлы 9-точечного шаблона типа "ящик". Последнее важно для унификации алгоритма вычисления решения в особых точках, которые в такой модифицированной постановке с клетчатой структурой области представляют собой точки пересечения границ клеток, В этом случае все типы особых узлов (точки сопряжения трех или четырех материалов, точки сопряжения двух материалов и внешней среды, уступы границы подобласти или внешней границы, углы области) находятся в углах клеток.
Разностные схемы, используемые для внутренних узлов клеток, общеизвестны и здесь не приводятся. Это две схемы приближенной факторизации — обычная чисто неявная схема точности 0(т + к2) на шаблоне "крест" и компактная схема точности 0(т2 + к4) на 9-точечном шаблоне со специальными значениями весов. Внутренние и внешние граничные условия аппроксимировались непосредственно путем замены нормальной производной па одностороннюю разделенную разность, записанную по (в + 1) точкам, с погрешностью аппроксимации 0(кя), максимально возможной на таком шаблоне, Для схемы второго порядка рассматривались варианты граничных условий с в =1 и ^ = 2, а для компактной схемы — более "длинн ые" условия с« = 3и в = 4, Таким образом, тестировались четыре различные разностные схемы, порядок которых менялся от единицы до четырех.
Внешние граничные условия записывались на одномерном шаблоне, состоящем из в + 1 узлов, а шаблон внутренних граничных условий занимал 2в + 1 узлов. Шаги сетки не должны быть чрезмерно крупными для того, чтобы в каждой клетке в обоих
в
потоков могут выступать за пределы границы однородности, что недопустимо, 2. Алгоритм реализации
В результате расщепления факторизованной схемы получается две совокупности одномерных линейных алгебраических систем с ленточной матрицей почти трехдиаго-нального вида. Отличие от трехдиагональной матрицы создают "длинные" строки, соответствующие многоточечным разностным граничным условиям. Например, для схемы четвертого порядка аппроксимации первая и последняя строки матрицы имеют по пять ненулевых элементов, а несколько изолированных строк, соответствующих границам раздела сред, содержат по девять элементов, отличных от нуля. Такая матрица за конечное число шагов элементарными преобразованиями приводится к трехдиагональной, для чего достаточно по отдельности преобразовать каждую "длинную" строку в трехточечную локальным применением процедуры исключения Гаусса, После этого задача решается обычной прогонкой.
На первом шаге прогонкой вычисляются промежуточные значения дробного шага по одному из направлений, например, по всем вертикальным линиям сетки, кроме границ клеток. На втором шаге проводим прогонку по всем горизонтальным линиям сетки также кроме граничных линий. Таким образом, после второй прогонки решение на верхнем слое остается не вычисленным на горизонтальных границах клеток. Множество узлов сетки этого множества состоит из двух подмножеств — подмножества не угловых узлов клеток (не являющихся точками пересечений с вертикальными границами клеток) и подмножества углов клеток (точек пересечения границ клеток). Для определения решения в первом подмножестве достаточно воспользоваться описанными выше разностными аппроксимациями граничных условий на верхнем временном слое, выразив из них значения искомой функции на границе. После этого решение на верхнем слое становится известным всюду, кроме углов клеток. Для определения решения в них использовались специальные формулы [3], аппроксимирующие исходное уравнение в узлах на стыке до четырех различных материалов. Не приводя здесь эти формулы, заметим, что они получены интегрированием разностных уравнений произвольного порядка аппроксимации по пространственным переменным по четырем смежным ячей-
(4в + 1)
3. Результаты численных экспериментов
С целью тестирования формул замыкания и численного исследования практически наблюдаемых порядков сходимости были проведены численные эксперименты. Расчеты проводились по четырем схемам рассмотренного типа — с первого до четвертого порядка аппроксимации. Эти схемы при описании результатов имеют номера, соответствующие порядку точности.
Тестовые задачи создавались следующим образом, В качестве начального распределения температур выбиралась гладкая функция двух переменных. На границах формировались граничные условия, согласованные между собой в угловых точках области и с начальными данными. Оценка практически наблюдаемого порядка точности проводилась на последовательности трех сгущающихся сеток по коэффициенту убывания ошибки
_ _ \\Цн - Цн/Лс
£н/2 \\иь,/2 — и^/4||с'
который теоретически при малых к должен стремиться к степеням двойки 2я, Чтобы число т/к2 оставалось постоянным, шаг по времени т при делении пространственного шага пополам уменьшался в четыре раза.
Точное решение в тестах не было известно. Однако во всех тестах при детализации сетки замечена сходимость результатов по всем схемам к одному пределу, поэтому за точное решение при расчете условно принималось численное решение, полученное на сетке с шагом к/4 по схеме четвертого порядка, которая сходилась быстрее других при уменьшении шагов сетки.
Первой задачей была задача Дирихле в единичном квадрате, в верхней правой четверти которого коэффициент теплопроводности А1 = 0.02, а в остальной части Л2 = 2, Теплоемкость принималась равной единице для обоих материалов. Начальные данные и граничные условия генерировались функцией
В этой задаче формула замыкания [3] действует в единственной точке — в центре квадрата, в котором расположен уступ одной из сред,
В табл. 1 приведены результаты численного исследования порядка сходимости по схемам различных порядков точности для данной задачи. Здесь показаны нормы ошибок на сгущающихся сетках для всех схем и коэффициенты убывания ошибок. Для всех порядков обнаруживается стремление коэффициента убывания ошибки к 2я, что соответствует теории.
Таблица 1. Экспериментальная оценка порядка точности в первой задаче при Н = 1/40
Схема расчета Ошибка е^ Ошибка 6^/2 Коэффициент г^
1 0.038298 0.019225 1.9921
2 0.030408 0.007102 4.2817
3 0.003418 0.000411 8.3181
4 0.003212 0.000201 15.9834
На рис, 1, а изображена поверхность, построенная по результатам расчета па момент времени £ = 1 по схеме четвертого порядка точности на сетке с шагом к = 0.025, На поверхности хорошо виден излом решения па границе раздана сред. Соответствующие изолинии приведены па рис, 1, б. Сплошные изолинии соответствуют "точному" решению, а маркерами отмечены результаты на более грубой сетке с шагом к = 0.1,
Вторая задача — смешанная краевая задача в единичном квадрате. Здесь линия раздана сред совпадала с границей нижней левой четверти квадрата, В пей коэффициент теплопроводности А1 = 0.02, в остальной области А2 = 2, В начальный момент времени температура во всей области была нулевая. На левой и нижней границах задавалась нормальная производная температуры
дх ду 10'
па остальной части внешней границы ставились условия Дирихле
Щх, и) = ^(993 - 158*2 - я + 60), Щ1,у,г) = ^(41 у3 - 62у2 + у + 20)
дня верхней и правой границ соответственно.
В данном эксперименте имеются четыре особые точки, в которых использовались формулы замыкания. Это углы нижней левой четверти квадрата.
Таблица 2. Экспериментальная оценка порядка точности во второй задаче при Н = 1/40
Схема расчета Ошибка 6/г Ошибка 6^/2 Коэффициент Г/г
1 0.030194 0.017058 1.7701
2 0.020274 0.005132 3.9507
3 0.004082 0.000511 7.9901
4 0.0020264 0.000128 15.83151
Рис. 2. Решение второй задачи на момент времени £ = 1 (а) и его изолинии (б)
В таб.и. 2 (см. с. 85) приведены результаты исследования порядка сходимости дня рассматриваемой задачи. Здесь также очевидна близость коэффициента убывания ошибки к 2я, что согласуется с теорией.
На рис. 2 показаны решение на момент времени £ = 1 и его изолинии.
Полученные результаты свидетельствуют о том, что несмотря па наличие разрывов производных решения дня краевых задач в кусочно-однородных областях клетчатой структуры применение форму.:: замыкания в особых точках по приводит к неустойчивости или локальным возмущениям, и разностные схемы с односторонней аппроксимацией потоков, входящих в граничные условия, па практике обнаруживают порядки точности, совпадающие с теоретическими.
Список литературы
[1| Паасопеп В.II. Параллельный алгоритм для компактных схем в неоднородных областях /7 Вычиел. технологии. 2003. Т. 8, № 3. С. 98 106.
[2] Ичетовкип Д.А., Паасопеп В.И. Численное исследование независимой аппроксимации граничных условий на решениях с разрывами производных /7 Там же. 2010. Т. 15, № 1. С. 77 84.
[3] Паасопеп В.И. Формулы замыкания для компактных схем в неоднородных областях /7 Там же. 2009. Т. 14, № 4. С. 93 99.
Поступила в редакцию 16 ноября 2009 г.
