CC BY
28
Вычислительные технологии
Том 15, № 5, 2010
О применении компактных схем для уравнения колебаний в кусочно-однородных средах*
В. И. Паасонен
Институт вычислительных технологий СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия e-mail: paas@ict.nsc.ru
Работа является продолжением цикла статей по применению компактных схем с универсальной высокоточной аппроксимацией граничных условий на линиях раздела сред для решения краевых задач с разрывными коэффициентами и заключается в распространении технологии, разработанной ранее для параболических и эллиптических уравнений, на краевые задачи для уравнения колебаний в составных областях. Предлагаются формулы повышенной точности для вычисления решения в особых узлах сетки (в углах области, в точках сопряжения нескольких материалов и т.д.), замыкающие высокоточный алгоритм в целом.
Ключевые слова: односторонняя аппроксимация потока, многоточечные граничные условия, компактная схема, кусочно-однородная среда, неоднородная область, схема высокого порядка точности, аппроксимация граничных условий.
В [1] были разработаны обычные и параллельные алгоритмы реализации разностных методов для уравнения теплопроводности в составных областях, опирающиеся на компактные схемы внутри однородных подобластей и на многоточечные одномерные аппроксимации потоков, входящих во внешние и внутренние граничные условия. Данная технология самодостаточна лишь для слоистых областей, а для областей более сложной структуры, имеющей особые точки (точки сопряжения трех или четырех материалов, изломы границ, уступы и т.д.), на каждом шаге счета необходимы дополнительные локальные вычисления в этих точках. Соответствующие формулы для параболических задач, замыкающие вычислительный алгоритм для широкого класса областей и обеспечивающие произвольный порядок точности, предложены в [2]. В настоящей работе на базе компактной схемы [3] технология [1] с идеей замыкания [2] распространена на случай краевых задач для уравнения колебаний в кусочно-однородных средах.
1. Постановка задачи
Для решения уравнения колебаний
* Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00264) и Интеграционного проекта СО РАН № 103.
д2п д2 u
дх2 ^ ду2 ^ ^ ^
в однородной области самой экономичной является, вероятно, компактная схема [3], имеющая четвертый порядок точности относительно шагов пространственной и временной сетки, В двумерном случае схема имеет вид
Р АгА^1 2^ + иП~1 = л ^ + Л2 + ЩА Л1 л2) + (1)
где
1 / ч ( Ъ2 Ъ2 \ т2 Jn+1 — 2 ^ п + /•
1 Л _2 / „ г,2А л т?п _ I тр , Ъ1 л , Ъ2 л \ т | т ^ + J
Аг = Е--(\г2/р-^)Аг, Рп=(е + -±А1 + -1А2)р +
п— 1
12^ " V " у 12 12 у 12 т2
Для ее применения в неоднородных областях, составленных из однородных подобластей специального вида с постоянной плотностью р и постоянным модулем Юнга Л, представляется целесообразным распространить идею [1] использования внутри однородных подобластей высокоточной схемы, а на границах раздела — непосредственной аппроксимации граничных условий с тем же порядком аппроксимации с применением одномерных односторонних аналогов первых производных на шаблоне с достаточным числом узлов.
Будем предполагать, что неоднородная область составлена из однородных подобластей с постоянными характеристиками, а однородные подобласти имеют вид объединения конечного числа прямоугольников, одинаково ориентированных вдоль осей. На внешней границе заданы граничные условия, согласованные между собой и с начальными данными для решения и и его производной по времени. На границах раздела сред перемещения и натяжения в нормальном к границе направлении предполагаются непрерывными:
= и+, 1=1. (2)
Продолжим все вертикальные и горизонтальные границы до пересечения со всеми границами ортогонального направления, заключив таким образом реальную область в прямоугольник, разделенный в обоих направлениях на клетки, часть которых может лежать вне области. Введем прямоугольную сетку, равномерную в пределах каждой клетки и согласованную с границами. Для удобства добавим по одному ряду фиктивных ячеек со всех четырех сторон прямоугольника, что позволит рассматривать любой узел сетки единообразно как центральную точку шаблона, окруженную четырьмя соседними ячейками, вне зависимости от того, является ли она внутренней или граничной, В построенной области клетчатой структуры наряду с действительными внутренними границами раздела сред имеются фиктивные, разделяющие одинаковые среды, С целью упрощения логики алгоритма эти границы удобно рассматривать как равноценные действительным, задавая на них также условия (2), которые в силу совпадения значений модуля Юнга по разные стороны границы раздела переходят на фиктивных границах в условия равенства односторонних нормальных производных, т, е, в условия гладкости решения.
При формальном совпадении условий (2) с условиями непрерывности температуры и потока в тепловых задачах представляется возможным аппроксимировать их совершенно так же, как это сделано в [1] при применении высокоточных схем для уравнения теплопроводности, т, е, односторонними разностями на многоточечных одномерных
шаблонах. При этом в каждой клетке сетка должна быть не слишком крупной, чтобы шаблон односторонней разностной аппроксимации первой производной не выходил за ее пределы, В частности, для схемы (1), чтобы сохранить четвертый порядок аппроксимации натяжения, необходимо аппроксимировать по пяти точкам. Поэтому в каждой клетке должно быть не менее четырех шагов сетки в каждом направлении.
Из самой схемы расщепления факторизованной схемы (1) при естественном определении дробного шага
„."+1 2ип + к"—1 мга+1/2 _ ^ ^_~Г »
т 2
следует, что для промежуточной величины во всех узлах вертикальных границ, кроме углов клеток, возникает такое же многоточечное одномерное разностное граничное условие.
Схема (1) расщепляется на две совокупности одномерных задач обычным образом, причем для обоих координатных направлений матрица одномерной системы является почти трехдиагональной, имеющей лишь фиксированное, по числу границ раздела сред, число изолированных строк с девятью ненулевыми элементами. Такие системы легко сводятся элементарными преобразованиями "длинных" строк к трехдиагональ-ным с последующим применением прогонок либо реализуются параллельно по клеткам [1] без предварительного приведения к трехдиагональному виду. Заметим, что по границам раздела сред (по границам клеток) прогонки осуществить невозможно, так как нет способа аппроксимации граничных условий в тех углах клеток, где сопрягается до четырех различных материалов. По этой причине после прогонок в горизонтальном и вертикальном направлениях значения решения на верхнем слое определяются всюду, кроме вертикальных внешних и внутренних границ.
Затем решение на верхнем слое определяется явно всюду на этих линиях, кроме углов клеток, непосредственно из многоточечных разностных граничных условий, записанных на верхнем слое и разрешенных относительно значения решения на границе. Для завершения процедуры вычисления решения на временном шаге требуется определить его значения в углах клеток. Если данный угол принадлежит внешней границе и на ней поставлено условие Дирихле, вычисления не требуются. Если угол фиктивный, т, е, лежит вне области, вычисления также не нужны, В противном случае расчеты необходимы, причем с адекватной точностью, С этой целью построим специальные формулы замыкания, аналогичные предложенным в [2].
2. Формулы замыкания повышенной точности
Формулы расчета решения в проблемных узлах получим путем обобщения обычного способа получения схемы сквозного счета, имеющей вид системы четырех разностных уравнений, записанных отдельно для каждой из четырех ячеек, образующих шаблон. Рассмотрим стандартный шаблон типа "ящик" с центром в фиксированном углу клетки. Обозначим предельные значения натяжений в центральном узле шаблона
пк-\ди пк-\ди
где индекс к = 1, 2, 3, 4 означает номер ячейки сетки в данном шаблоне с нумерацией по часовой стрелке и началом в первом квадранте, а А& — модуль Юнга для каждой из
четырех соседних ячеек. Координатное направление х или у и принадлежность к одной из четырех ячеек определяет восемь различных предельных значений натяжений, когда точка стремится по ребру ячейки к центру шаблона. Местные значения шагов сетки обозначим через Н—х, Н+х и Н—у, к+у.
Запишем разностную аппроксимацию волнового уравнения для каждой из четырех ячеек:
= - <£) + (^Д+уИ - + /ь (3)
Р2— = — (А2А+ЖМ -(£)+ — (д2у - А2А_уи) + /2, (4)
дг2 Н+х
д2и 2
дг2 ~ Н+х
д2и 2
дг2 ~ Н — х
д2и 2
2
Рз— = — Й " АзА_жм) + — - А3А_ум) + /3, (5)
Н — у 2
= — - А4А_ЖМ) + — (А4А+УМ - + /4. (6)
Если бы Д+ и Д- были здесь простыми двухточечными операторами односторонних разделенных разностей, то каждое из уравнений (3)-(6) аппроксимировало бы волновое уравнение в своей ячейке с первым порядком. Обобщение же состоит в том, что эти операторы означают специальные многоточечные разделенные односторонние разности вперед и назад по соответствующей переменной. Например, разность вперед по переменной х задается в виде
1 5
к иг+к^,
Н+
ь+х к=0
где коэффициенты в+ однозначно определяются из условия аппроксимации выражения
ди 1г+х д2и дх 2 дх2
с порядком Для их вычисления используется система уравнений
]гв+к = + , з = о, 1,...,в, к=0
где — символ Кронекера, Решение системы уравнений имеет вид
# = -Цг^ ^ + "Л. МО, # = (7)
Аналогично определяется левосторонняя разность
10-Д-ж^г,;/ = Т Е
Н—х к=-в
от которой требуется аппроксимация с порядком 8 выражения
ди Н—х д2и
дх 2 дх2
Ее коэффициенты отличаются от (7) только знаком:
в— = -в+, к = 0,1,...,в.
При описанном определении операторов Д+ и Д_ погрешность аппроксимации волнового уравнения в каждой ячейке составляет величину О (к8). Для нашего случая достаточно взять в = 4, при этом каждое из уравнений системы (3)-(6), как и основная схема (1), имеет четвертый порядок аппроксимации.
Систему (3)-(6) следует дополнить условиями непрерывности натяжений
оХ
Ях,
оХ
ЯХ ,
о1
Я2,
Яу
ЯХ
(8)
на границах смежных ячеек в горизонтальном и вертикальном направлениях. Разностная схема из уравнений (3)—(6), (8) получается их интегрированием по всем непустым ячейкам шаблона с последующей аппроксимацией по дискретной временной переменной, Пусть, например, центральный узел шаблона является внутренним. Тогда пустых ячеек нет и все четыре уравнения имеют смысл. Умножив уравнения (3)-(6) на площади соответствующих ячеек и сложив полученные результаты, а затем разделив полученное соотношение на суммарную площадь четырех ячеек, в силу условий (8) получим аппроксимацию волнового уравнения с погрешностью О (к8):
д2и „ „ „
(9)
где
2 (Х-\-х^-\-х — Х-хА_х)
к+Х + к — Х
Л
+Х
Л\к+у + Л2 к_у к+у+ к—у
Л_х
+ Аз Н-ь к+у+ к—у
а
2 (Л+уА+?/ — \-yA_y) _ Х\Ь+х + Л4Н-а
' _ кх + н.х
Лу
Аг^+ж + Аз ¡1-2
к+Х + к —Х
р и f — усредненные по четырем ячейкам плотность и правая часть (9), Выражение для плотности имеет вид
р
к+хк+ур1 + к+Хк— уР2 + к—хк—уРх + к—хк+уР4
(к+х + к—х)(к+у + к—у) '
аналогично вычисляется и правая часть.
Теперь необходимо преобразовать (9) в трехслойную схему, вводя дискретизацию временной переменной с шагом т и усредняя слагаемые в правой части (9) по трем слоям, Заметим, что данная задача существенно отличается от параболического случая, где достаточно усреднения полусуммой, как в схеме Кранка—Ннколсон, поскольку там
т
О(т4 + к4), в силу чего при усреднении в правой части (9) важно добиться также четвертого порядка по времени, С этой целью воспользуемся произволом в выборе параметра
симметричного усреднения по трем слоям (с весами в (1 — 2,5) и в), Ввиду симметрии второй порядок имеет место при любом в Нетрудно проверить, что при специальном значении веса , = 1/12 достигается четвертый порядок:
ип+1 _ 2ип + ип-1 п+1 + 10ц"+ Ц"-1 /га+1 + 10/" + /"-1
Р -у-2 ^ 1 ^ • ' ^
Уравнение (10) относится только к углам клетки и занимает шаблон типа "большой крест" по трем временным слоям. Оно и используется для вычисления решения в углу. Для этого соотношение (10) необходимо разрешить относительно значения решения в центре шаблона на верхнем слое, В результате получится именно явная формула, так как решение на верхнем слое всюду вне углов клеток к этому моменту уже найдено на предыдущих этапах. Можно было поступить иначе — построить трехслойные раз-
1/12
ячейкам, и в этом случае приходим точно к соотношению (10),
Аналогично получаются разностные соотношения в любых внешних углах клеток, если там поставлены не условия Дирихле, а значения нормальных производных к граням ячеек. Пусть, например, узел является граничным и пусть в шаблоне только вторая ячейка пустая. Тогда из системы выпадает уравнение (4), и при интегрировании трех других уравнений (или их дискретных аналогов) для исключения величин д2 и дУ потребуется воспользоваться условиями на граничных ребрах, когда точка границы стремится к центральному узлу шаблона.
Если в окрестности граничного узла отсутствуют три ячейки из четырех (пусть это будет вторая, третья и четвертая ячейки), тогда в системе (3)-(6) остается только первое уравнение, а в системе (8) — первое и третье, все остальные фиктивны, И в этом случае все дополнительные величины при интегрировании исключаются путем привлечения граничных условий для ?2" При сопряжении материалов на горизонтальном или вертикальном участке границы в системе присутствуют всего две ячейки из четырех, тогда из системы (3)-(6) выпадает два уравнения, а из (8) — одно, связывающее величины в фиктивных ячейках, и система снова замыкается. Во всех случаях процедура построения формул для вычисления решения в узле, где сопрягается до четырех различных материалов, одинакова — интегрирование уравнений (3)-(6) по присутствующим ячейкам шаблона и усреднение результирующего уравнения по трем слоям с 1/12
Построенные формулы заполняют пробел, связанный с неопределенностью в вычислении решения в узлах сопряжения нескольких материалов и в других особых точках, и существенно расширяют круг краевых задач, которые могут быть решены с повышенным порядком точности с помощью компактной схемы для уравнения колебаний. Очевидно, что совершенно аналогично проблема решается и при большем числе пространственных переменных, В трехмерном пространстве, например, уравнений вида (3)-(6) для отдельных ячеек пространственного шаблона будет восемь, значений натяжения в центре шаблона по три для каждой из восьми ячеек, и для них двенадцать естественных связей на границах раздела сред. Из данной системы уравнений предложенным способом получается формула повышенной точности для вычислений во внутренних углах клеток. Во внешних особых точках это также удается сделать с привлечением граничных условий.
Список литературы
[1] Паасонен В.И. Параллельный алгоритм для компактных схем в неоднородных областях // Вычисл. технологии. 2003. Т. 8, № 3. С. 98-106.
[2] Паасонен В.И. Формулы замыкания для компактных схем в неоднородных областях // Там же. 2009. Т. 14, № 4. С. 93-99.
[3] Валнуллнн А.Н., Паасонен В.И. Экономичные схемы повышенного порядка точности для многомерного уравнения колебаний // Численные методы механики сплошной среды. 1970. Т. 1, № 1. С. 17-30.
Поступила в редакцию 16 ноября 2009 г.
