CC BY
36
Вычислительные технологии
Том 14, № 4, 2009
Формулы замыкания для компактных схем в неоднородных областях*
В. И. ПААСОНЕН Учреждение Российской академии наук Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: paas@ict.nsc.ru
Данная работа является продолжением цикла статей, посвященных применению для решения краевых задач с разрывными коэффициентами компактных схем, в которых внешние и внутренние граничные условия аппроксимируются с повышенным порядком. В работе формулируется замыкание этого алгоритма в виде специальных высокоточных формул вычисления решения на пересечениях границ раздела сред, в угловых точках внешней границы и в других особых узлах сетки, покрывающей неоднородную область клетчатой структуры.
Ключевые слова: односторонняя аппроксимация потока, многоточечные граничные условия, многоточечная аппроксимация потока, компактная схема, неоднородная область, схема высокого порядка точности, аппроксимация граничных условий.
Введение
В [1, 2] разработаны теоретические основы обычной и параллельной реализации разностных методов, опирающихся на компактные схемы внутри однородных подобластей и на многоточечные одномерные аппроксимации для потоков во внешних и внутренних граничных условиях. В одномерном случае, а также для наиболее простых двумерных задач, например, для задачи Дирихле в областях, полученных разрезанием прямоугольника на полосы, занятые различными материалами, технология [1, 2] самодостаточна и никаких дополнительных ухищрений не требуется. Если же область представляет собой произвольное конечное объединение однородных прямоугольников, то в узлах, где сопрягается более двух материалов или где граница раздела двух материалов имеет излом, возникает необходимость специальных вычислительных процедур. Такая же ситуация имеет место в тех узлах на внешних границах, где встречаются два материала или где граничный узел соответствует излому границы области, если, конечно, в этих узлах поставлены не условия Дирихле. В данной работе предлагается замыкание алгоритма для вычисления решения в таких проблемных узлах, исключающее потерю точности.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-00264). © ИВТ СО РАН, 2009.
1. Постановка задачи
Если расчетная область для задачи Дирихле составлена из конечного числа одинаково ориентированных прямоугольников, а границы раздела сред параллельны друг другу, т.е. неоднородность имеет не клетчатую, а полосчатую структуру, то сочетание компактных схем повышенной точности внутри однородных подобластей и многоточечных односторонних аппроксимаций той же точности для аппроксимации потоков на границах раздела сред [3] составляет замкнутый алгоритм, не требующий дополнительных вычислительных процедур. Схемы в однородных подобластях расщепляются на трехточечные одномерные, на первом шаге (вдоль полос) по границам раздела сред прогонки не проводятся. На втором дробном шаге прогонки поперек полос осложняются тем, что границам раздела сред вместо трехточечных соотношений соответствуют изолированные многоточечные. Проблема решается применением либо метода исключений Гаусса для пошагового превращения многоточечного условия в трехточечное [1], либо параллельного алгоритма [2]. Численные эксперименты [3] подтверждают совпадение теоретического порядка сходимости с практически наблюдаемым. Если же неоднородность имеет структуру более общего вида, подобно изображенной на рисунке сплошными линиями, то внутри области появляются проблемные узлы, где сопрягается более двух материалов (точки А). Описанный выше алгоритм не дает решения в этих узлах, и его необходимо доопределять как-то иначе, но без потери точности. Однако экстрапо-
Е
С
С
Е
В
А
В
В
А
В
4 ¡1-х /з+у 1 Ь+х
3 Ь-у 2
С
Е
Схема неоднородной области с особыми точками и шаблон разностной схемы
ляция, например, — не самая хорошая идея, поскольку она оказывается источником локальных ошибок, несмотря на формально повышенный порядок экстраполяционных формул. Если на внешней границе заданы не условия Дирихле, а условия Неймана или условия третьего рода, то проблемные точки возникают также и на границе. Это точки сопряжения материалов О, С и угловые точки В и Е (см. рисунок).
Для достижения упрощения логики будущего алгоритма продолжим границы раздела сред на всю область, на рисунке продолжения обозначены штриховыми линиями. Заметим, что внешние границы также полезно продолжить, дополнив таким образом любое объединение одинаково ориентированных прямоугольников до одного большого прямоугольника, рассеченного вдоль и поперек на отдельные клетки. Откажемся от аппроксимации дифференциального уравнения на фиктивных (штриховых) границах, заменив ее аппроксимациями равенства потоков, аналогично соотношениям на реальных границах раздела сред. Эта возможность, теоретически обоснованная в [4], в действительности продуцирует естественные "мягкие" условия гладкости в виде равенства односторонних производных, ввиду того, что по разные стороны штриховой линии теплофизические характеристики одинаковы. Но этот важный шаг приводит к равноправности всех границ клеток и позволяет формально отождествить между собой узлы вида А, и не только их, что сильно упрощает логику алгоритма. Обратим внимание на различие всех трех узлов А на рисунке — в разных точках этого типа граничат два, три и четыре материала. Из рисунка видно, что и другие типы узлов обладают этим свойством.
Расщепим двумерную факторизованную задачу на две одномерные. Сначала выполним прогонки по горизонтальным линиям сетки, пропуская все границы раздела, реальные и фиктивные, затем — по всем вертикальным, также кроме всех границ. После этого остается найти решение лишь на вертикальных границах. Это можно сделать по явным формулам из многоточечных одномерных граничных условий на этих границах, но, к сожалению, не во всех узлах — на пересечениях границ, т. е. в узлах А, В, С, О и Е это сделать невозможно. В самом деле, рассмотрим, например, узел О на правой границе. Вычисление решения в нем из граничного условия требует задания коэффициента теплопроводности, а он терпит разрыв на границе раздела. Использование усредненного значения допустимо для схемы второго порядка, а для схемы четвертого порядка необходимы более точные вычисления. Многоточечная экстраполяция по границе в точку О тоже не обеспечивает нужную точность, к тому же в этом способе нет однозначности — экстраполировать можно по границе раздела или по внешней границе сверху либо снизу. Те же проблемы имеют место во всех угловых узлах клеток. Таким образом, для замыкания изложенного выше алгоритма необходимо организовать вычисление решения в проблемных узлах на верхнем слое, предположив, что решение на нижнем слое известно всюду, а решение на верхнем слое — всюду, кроме этих изолированных узлов в углах всех клеток. Заметим, что на верхнем слое всюду в окрестности этих узлов решение уже известно, поэтому, если удастся сформулировать высокоточное неявное разностное соотношение в проблемных узлах, воспользоваться им можно будет как явной формулой.
2. Схема сквозного счета
Формулы расчета решения в проблемных узлах получим, обобщив формулы простой схемы сквозного счета [5]. Эта схема имеет вид системы четырех простых разностных
уравнений, записанных отдельно для каждой из четырех ячеек, образующих шаблон, изображенный в правом нижнем углу рисунка. Обозначим предельные значения потоков в окрестности центрального узла через
^ = А ^ = А
д- = АкдХ' ^ = Акду'
где индекс к = 1, 2, 3,4 означает номер ячейки сетки в данном шаблоне, с нумерацией по часовой стрелке и началом в первом квадранте, а коэффициент теплопроводности А& предполагается постоянным в пределах каждой ячейки. Координатное направление х или у и принадлежность к одной из четырех ячеек определяют восемь различных значений потока в центральной точке шаблона. Шаги сетки будем считать переменными и обозначим через , и Т_у, Т+У их местные значения.
Запишем разностную аппроксимацию уравнения теплопроводности для каждой из четырех ячеек:
с1 дт = (А1Л+-и - 9-) + (А1Л+уи - + Л' (1)
С2ди = (А2А+Жи - д-) + ^ («2 - А2А_уи) + /2' (2)
Сз ди = тЬ (д- - АзЛ_-и + («3 - Азд_у и + /з' (3)
— у
с4= ^ (д- - А4Л_-И) + т+У (А4Л+УИ - + Л. (4)
Здесь операторы Л+ и Л_ означают обычные двухточечные разделенные односторонние разности вперед и назад по соответствующей переменной, а / — плотность объемного источника тепла в центральном узле со стороны ячейки с номером к. Каждое из уравнений (1)-(4) аппроксимирует с первым порядком уравнение теплопроводности
ди (д2и д2и\
Ск т = \ дХ2 + дУ^; + /к
в соответствующей ячейке. Систему (1)-(4) следует дополнить условиями непрерывности потоков
д- = д- = = ду, ду = дз (5)
через границы смежных ячеек в горизонтальном и вертикальном направлениях.
В [5] такого рода система рассматривалась для уравнения с переменными теплофи-зическими характеристиками в криволинейной ортогональной системе координат. Для наших же целей достаточно декартовой системы и кусочно-постоянных коэффициентов. Разностная схема из системы уравнений (1)-(5) получается интегрированием по всем участвующим ячейкам. Пусть, например, центральный узел шаблона является внутренним. Тогда все четыре уравнения имеют смысл, поскольку все ячейки заняты материалами, вообще говоря, различными. Умножив уравнения (1)-(4) на площади соответствующих ячеек и сложив полученные результаты, а затем разделив полученное соотношение на суммарную площадь четырех ячеек, в силу непрерывности потоков (5) получим
ди
с— = Л-и + Лу и + /' (6)
где
Лх
Лу
2
к+х + к—X 2
к+у + к—у
(Л+хА +х Л —X А —X) ) Л+х
(Л+уА+у — Л—уу) , Л+у :
Л^+у + Л2^ —у
к+у + к—у
Л1Л.+Х + Л4к—х
к+х + к—X
Л—х
Л—у
Л4к+у + Лзк—у к+у + к—у
Л2к+х + Лзк—х
к+х + к—х
, (7)
. (8)
Здесь с и f — осредненные по четырем ячейкам теплоемкость и плотность источников. Выражение для теплоемкости имеет вид
с =
к+хк+у с1 + к+хк— у С2 + к—хк—у Сз + к —хк+у с4
(к+х + к—х)(к+у + к—у) ;
и плотность источников вычисляется точно так же.
После этого остается ввести дискретизацию времени с шагом т, и тогда симметричная аппроксимация по двум временным слоям превращает (6) в аналог схемы Кранка— Николсон второго порядка по времени.
Аналогично получаются разностные соотношения во всех граничных точках, когда на границе поставлены не условия Дирихле, иначе случай становится тривиальным. Например, предположим, что центральный узел шаблона граничный. Пусть это будет, для определенности, точка В на рисунке. Тогда первая, третья и четвертая ячейки в шаблоне граничные, а второй ячейки нет. Следовательно, из системы выпадает уравнение (2). Остальные интегрируем, как прежде, только в дополнение к условиям (5) для
2з
исключения потоков дх и ду следует воспользоваться пределами граничных условий на граничных ребрах, когда точка стремится по ним к центральному узлу шаблона.
Если в окрестности граничного узла отсутствуют три ячейки из четырех (пусть это будет вторая, третья и четвертая ячейки), тогда в системе (1)-(4) остается только первое уравнение, а в системе (5) — первое и третье, все остальные фиктивны. И в этом случае все потоки благополучно исключаются с привлечением граничных условий для потоков ^х к граничным ребрам, которые получаются при стремлении точки по границе в центр шаблона, который в данном случае является угловой точкой Е области.
В случае горизонтального или вертикального участка границы (узлы С и Д) в системе присутствуют ровно две ячейки из четырех (верхние, нижние, левые или правые), из системы (1)-(4) выпадает два уравнения, а из (5) — одно, связывающее потоки в выпавших ячейках. Внешние потоки, нормальные к границе, опять исключаются из граничных условий. Во всех перечисленных случаях процедура одинакова — интегрирование уравнения по присутствующим ячейкам шаблона, нормировка результата к площади шаблона, дискретизация времени, аппроксимация уравнения по двум слоям со вторым порядком по времени.
3. Высокоточные формулы замыкания алгоритма
Описанная схема однородна, поэтому реализуется сквозным образом, но имеет всего первый порядок точности. Однако позитивным свойством является ее универсальность в том смысле, что характер расположения внешних границ и внутренних границ раздела сред по ребрам шаблона может быть совершенно произвольным. Поэтому представляется весьма заманчивым сохранить это полезное свойство универсальности и в то же время повысить точность формул (1)-(4), а результирующие схемы, получаемые из них
в точности описанным выше способом, использовать в качестве формул для вычисления с высокой точностью решений в проблемных узлах. Тем самым мы получим замкнутый высокоточный алгоритм для широкого класса краевых задач в неоднородных областях рассматриваемого типа. Именно это и делается ниже.
Заметим, что прямая замена операторов Д±х, Д±у многоточечными аппроксимациями первых производных с повышенным порядком совершенно не годится, поскольку в этом случае вообще теряется аппроксимация уравнения теплопроводности. В самом деле, в оригинальных формулах (1)-(4) она достигается именно за счет главной части остаточного члена разностей вперед и назад. Например, разложение первого слагаемого правой части (1) в окрестности центра шаблона дает
2 /, „ -| \ 2 (^ I ди + Н+х д2и
(л'д+- - «о=£ (Чш+¥Ш+Н - «О=А.+
и+х V 1 Н+Л V дх 2 дх2
и так аналогично для каждого слагаемого уравнений (1)-(4).
Отсюда видно, что для достижения цели необходимо в уравнения (1)-(4) вместо односторонних разностей Д±хи, Д±уи подставить всюду высокоточные аппроксимации выражений вида
ди Н±х д2и ди Н±у д2и , ,
— ±--— ± —--. (9)
дх 2 дх2' ду 2 ду2'
А это можно сделать с любой точностью.
Запишем, например, для правосторонней разности по х в общем виде одномерное разностное выражение вида
Д+хи(хиУз) = ^ вк и(х1+к, уj). (10)
к=0
Выбрав достаточно большое натуральное в (для компактных схем четвертого порядка достаточно взять в = 4) и определив подходящим образом набор коэффициентов вк, получим нужный порядок аппроксимации. А именно, если коэффициенты удовлетворяют системе линейных алгебраических уравнений
8
вк гк = ¿у + Zl52j, х? = х*+к - х*, ( = 0,(11)
к=0
где 8ту — символ Кронекера, а г. = — х^ = Н+х, то
ди Н+х д2и . .
Д+хи = дх + "Т дх + 0(Н+х)-
Очевидно, все г к попарно различны, следовательно, определитель Вандермонда системы (11) отличен от нуля, поэтому она имеет единственное решение. Его легко вычислить в общем виде по правилу Крамера. В важном частном случае, когда сетка равномерна в пределах прямоугольных однородных подобластей (клеток), решение имеет вид
вк = (—Н^С? (1 + 1 £ (-1)У (к - () сА , к = 0, во = - £ вк. (12) кН+х У к у = 1 ( У к=1
Очевидно, что для левосторонней производной шаг Н+х следует заменить на Н-х и поменять общий знак выражения.
Подстановка операторов такого вида в (1)-(4) даст уравнения, аппроксимирующие уравнение теплопроводности с порядком s в соответствующей элементарной ячейке. Тогда построение искомых высокоточных формул замыкания общего алгоритма совершенно совпадает с описанной выше процедурой построения двухслойной разностной схемы сквозного счета.
Список литературы
[1] PAASONEN V.I. Compact difference schemes for inhomogeneous boundary value problems // Russian J. Numer. Analys. and Math. Modelling. 2004. Vol. 19, N 1. P. 65-81.
[2] ПААСОНЕН В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. T. 8, № 3. C. 98-106.
[3] Ичетовкин Д.А., Паасонен В.И. Численное исследование независимой аппроксимации граничных условий на решениях с разрывами производных // Вычисл. технологии. 2009. В печати.
[4] Паасонен В.И. Параллельные алгоритмы на основе мягких внутренних граничных условий // Вычисл. технологии. 2006. T. 11. Спецвыпуск, посвященный 85-летию со дня рождения Н.Н. Яненко, ч. 2. C. 21-27.
[5] Паасонен В.И. Моделирование тепловых процессов в неоднородных конструкциях с источниками тепла // Изв. СО АН СССР. Серия технических наук. 1980. T. 1. C. 108-112.
Поступила в редакцию 12 января 2009 г.
