CC BY
10
Вычислительные технологии
Том 21, № 5, 2016
Варианты метода коллокации и наименьших невязок
__U f* U U
повышенной точности в области с криволинеинои границей
В. П. Шапеев1'2, В. А. Беляев1'2'*
1 Новосибирский национальный исследовательский государственный университет, Россия 2Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия e-mail: belyaevasily@mail.ru
Предложены и реализованы новые варианты метода коллокации и наименьших невязок повышенной точности для численного решения краевых задач для уравнений с частными производными в областях с криволинейной границей. Их реализация и численные эксперименты выполнены на примерах решения уравнений Пуассона и бигармонического уравнения. Решение второго уравнения использовано для моделирования напряженно-деформированного состояния изотропной пластины, находящейся под действием поперечной нагрузки.
Ключевые слова: метод коллокации и наименьших невязок, неканоническая область, повышенный порядок аппроксимации, уравнение Пуассона, бигармоничес-кое уравнение.
Введение
В вычислительной математике немало задач, для которых при ограниченных вычислительных средствах можно добиться такой малой погрешности решения методом повышенной точности (МПТ), которая недостижима методом низкой точности (МНТ). При этом зачастую, когда требуется найти решение задачи с некоторой заданной погрешностью, применением МПТ за счет выбора параметров метода это можно сделать за существенно меньшее время счета, чем с использованием МНТ [1,2]. Несомненно, для многих нужд в технике не требуется повышенная точность расчета. Однако в истории астрономических открытий и в настоящее время в задачах моделирования физических процессов можно назвать ряд примеров, когда требовалось изобретать вычислительные методы повышенной точности, чтобы решить конкретную задачу с требуемой точностью. Собственно в самой вычислительной математике методы повышенной точности позволяют строить решения задач, которые могут применяться в качестве тестов для выяснения возможностей вновь создаваемых и существующих методов [3-6].
Многие численные методы, дающие высокую точность решения краевых задач для уравнений с частными производными (PDE) в канонических областях, не имеют аналогов в областях с криволинейной границей. Во втором случае их создание бывает связано с серьезными затруднениями. Это можно видеть на примере спектрального метода. В разностных методах есть несколько приемов построения схем повышенной точности
© ИВТ СО РАН, 2016
решения РВЕ в областях с криволинейной границей. Один из приемов в разностных методах заключается в том, что для расчета значения решения в узлах разностной сетки около границы расчетной области используются "законтурные" (лежащие вне расчетной области) узлы разностной сетки [7,8]. Относительно простые подходы при построении вариантов повышенной точности метода коллокации и наименьших невязок (КНН) в областях с криволинейной границей использованы в [9,10].
В настоящей работе при построении варианта повышенной точности метода КНН в области с криволинейной границей реализованы более простые, несколько различающиеся между собой два подхода использования законтурных частей ячеек прямоугольной сетки, отсеченных криволинейной границей области решения задачи. Чтобы такой подход был применим, как и в случае других численных методов, предполагается, что решаемое уравнение справедливо не только в расчетной области, но и в малой окрестности ее границы. Рассматриваемые здесь варианты метода КНН и реализующие их компьютерные программы оказались проще, чем реализованные в [2,10]. Возможности предлагаемого здесь подхода сначала были проверены на решении достаточно простой и хорошо обусловленной задачи Дирихле для уравнения Пуассона, а затем на решении краевых задач для бигармонического уравнения, которые описывают напряженно-деформированное состояние (НДС) под различными нагрузками изотропной пластины, имеющей в плане форму фигуры с криволинейной границей. Задача о НДС прямоугольной пластины методом КНН с повышенной точностью решена в [11]. Можно отметить, что многие явления в реальном мире, которые моделируются с помощью численных методов, происходят в областях с криволинейными границами. Отсюда настоятельная потребность иметь возможность моделировать эти явления в областях с формой, соответствующей действительности, чтобы получать достоверные результаты.
Описание, характеристика, история возникновения, некоторые возможности метода КНН и библиография по нему приведены в [5,6, 12-15]. Здесь отметим, что в методе КНН исходной дифференциальной задаче путем проектирования в конечномерное линейное функциональное пространство сопоставляется дискретная задача, аппроксимирующая исходную. Решение дискретной задачи в каждой ячейке расчетной сетки, которой накрывается область решения задачи, ищется в виде линейной комбинации с неопределенными коэффициентами базисных элементов пространства. В качестве последнего в данной работе взяты пространства полиномов второй и четвертой степени. В методе КНН уравнениями дискретной задачи являются уравнения коллокации — требование удовлетворения решением дискретной задачи уравнениям дифференциальной задачи, условия согласования на общих сторонах, принадлежащих двум соседним ячейкам, и краевые условия на границе, если ячейка является граничной. Нелинейные уравнения дифференциальной задачи предварительно линеаризуются. В каждой ячейке выписывается переопределенная система уравнений, после подстановки в которую искомого приближенного решения из нее получается переопределенная система линейных алгебраических уравнений (СЛАУ) для определения коэффициентов представления в линейном функциональном пространстве решения дискретной задачи. Решение переопределенной СЛАУ ищется из требования минимизации функционала невязки, составленного из невязок всех его уравнений. В качестве функционала обычно берется сумма квадратов невязок. Его минимум может отыскиваться разными способами. Например, он достигается на построенном ортогональным методом псевдорешении переопределенной СЛАУ [14]. Детальное описание алгоритма метода КНН будет дано ниже.
1. Задача Дирихле для уравнения Пуассона
1.1. Постановка задачи и описание метода
Рассмотрим задачу Дирихле для уравнения Пуассона
I
Аь = /(хх, х2), (хх, х2) е О, V!<т = 9(х1, х2)
(1)
в двумерной области О с криволинейной границей 5О (рис. 1), где V(хх, х2) — искомая функция, /(хх, х2) и д(хх, х2) — известные функции. Для тестовых расчетов реализованного здесь подхода построения варианта метода КНН взята область, составленная из частей кругов разных радиусов при различном взаимном расположении (рис. 1, справа). Ее граница выделена жирной линией. Реализация рассматриваемого метода решения других краевых задач для уравнения Пуассона осуществляется аналогично описанному далее.
Впишем расчетную область О в прямоугольник, две стороны которого параллельны оси хх и проходят через принадлежащие 6О точки с минимальной и максимальной координатами по оси х2, а две другие стороны, параллельные оси х2, проходят через принадлежащие 6О точки с минимальной и максимальной координатами по оси хх. Полученный прямоугольник покроем регулярной сеткой с прямоугольными ячейками размера 2Кх х 2К2 (рис. 2 и 3). Прямоугольные ячейки, которые полностью расположены внутри области, в том числе и те, стороны которых имеют общие части с 5О, назовем внутренними. Куски, отсеченные границей области 5О от прямоугольных ячеек сетки и принадлежащие расчетной области, назовем граничными ячейками. Здесь имеют место различного вида нерегулярные граничные ячейки. Их стороны, принадлежащие расчетной области, назовем внутренними, а часть дуги границы области 6О, расположенную между отсеченными сторонами нерегулярной ячейки, назовем условно ее внешней стороной. Некоторые нерегулярные граничные ячейки могут быть треугольного вида с двумя прямыми сторонами, прямым углом между ними и противолежащей внешней стороной — дугой границы области. Другие нерегулярные ячейки могут быть с одной, тремя или четырьмя прямыми сторонами в расчетной области. Прямоугольную ячейку сетки, от которой границей области отсечена нерегулярная ячейка, назовем материнской.
Х2 А
Х2 А
1
1
О
Рис. 1. Область решения задачи
Для удобства реализации метода в каждой внутренней и материнской ячейке введем локальные координаты
(хг - Ху) (Х2 - X2j )
Уг =-;-—, У2 = '
кг
(2)
где (хг^, Хц) — центр у-й ячейки; у = 1,..., N, N — количество расчетных ячеек, в которых строится приближенное решение. Задача (1) после замены (2) в локальных переменных в ячейке Щ примет вид
1 д2 и 1 д2 и
Щ2 ~ду2 + Щ Ъу2 = ^ + Ш Ы1, Х2 * + У2]г2), + Ш Ыг' Х2 * + ^2) £ Щ, (3)
и = д (хц + уг кг, Х23 + У2 ^2), (хц + Уг кг, хц + У2 к2) £ $ Щ П Щ.
В каждой 2-й ячейке сетки приближенное решение и^ задачи (3) ищем в виде линейной комбинации базисных элементов пространства полиномов второй степени от двух переменных. Здесь в качестве базисных элементов ф^ (г = 1,..., 6) взяты мономы:
инз (уг, У2) = Ьцфг = Ъц + ЪцУг + Ь31 у2 + Ъцуг2 + Ь^У22 + b6jугУ2.
(4)
ъ=1
Коэффициенты Ъц определяются из уравнений коллокации, условий согласования на общих сторонах, принадлежащих двум соседним ячейкам и краевых условий на 6Щ, если ячейка является граничной.
Уравнения коллокации в каждой ячейке выписываются в точках , ^) и имеют
вид
1 а2 ) + Л я^ф) = ;
к!
дуг2
ду22
(5)
Рис. 2. Типы ячеек: а — внутренняя ячейка, б-д — четыре типа граничных ячеек; • — точки записи уравнений коллокации, штрихи — точки записи условий согласования, х — точки записи краевых условий
где Хк] = + hk^к], к = 1, 2. Если в нерегулярной ячейке выписываются уравнения коллокации (это уточняется далее), то используются точки, которые соответствуют точкам коллокации в материнских ячейках.
В качестве условия согласования решения в двух точках на общей стороне между соседними ячейками требуем непрерывность линейной комбинации с весами функции и и ее производной по нормали
ди+
кг—--+ к2и+
ОПп
к
ди дп1
+ к2и .
(6)
Здесь п^ — внешняя нормаль к границе ]-й ячейки; и+ и и~ — пределы функции при стремлении изнутри и извне к границе ^'-й ячейки соответственно; кх, к2 — весовые параметры.
Если сторона ячейки совпадает с некоторым куском границы области, в нескольких точках на ней выписываются краевые условия
к3и = кзд(хг, Х2),
(7)
где к3 — весовой параметр краевого условия в переопределенной системе уравнений.
Как показано в [16], значения весовых параметров влияют в некоторых пределах на обусловленность СЛАУ, из которой определяются коэффициенты приближенного решения, скорость сходимости итераций и величина погрешности приближенного решения задачи. Но погрешность приближенного решения определяется прежде всего порядком аппроксимации дифференциальной задачи дискретной задачей и зависит от значения старшей степени полиномов в базисе используемого для этого пространства [6]. Как показывает практика решения краевых задач методом КНН, во многих случаях можно положить кх = 1, к2 = 1, к3 = 1 и получить решения, сходящиеся на последовательности сеток, без необходимости поиска их оптимальных значений.
В каждой расчетной ячейке запишем четыре уравнения коллокации в четырех точках с локальными координатами (±д, ±5). В данной работе численные эксперименты проведены при д = 0.5. В граничных ячейках точки коллокации могут выходить за границу 8О.
Рис. 3. Фрагмент расчетной области: • — точки записи уравнений коллокации; штрихи — точки записи условий согласования; х — точки записи краевых условий; о — центры граничных ячеек, расположенные вне области О
Для записи условий согласования в любой ячейке или краевых условий для граничной ячейки ее каждая прямолинейная и внешняя стороны делятся на три отрезка, длины которых составляют 1/4, 1/2 и 1/4 части от длины стороны ячейки. На концах среднего (в указанном здесь порядке) отрезка записываются условия согласования или краевые условия соответственно, если взятая сторона оказалась внешней. Дуга внешней стороны ячейки может делиться на отрезки указанной длины приблизительно. В случае, если граничная ячейка не имеет три прямолинейные стороны, то для компенсации недостающих условий согласования на ее внешней границе дополнительно приблизительно равномерно расставляются точки для записи краевых условий в таком количестве, чтобы число уравнений, определяющих решение в ней, как и в других ячейках, было двенадцать. В случае, когда граница области такова, что при движении в пределах любой нерегулярной ячейки вдоль дуги ее внешней стороны касательная к дуге поворачивается на небольшой угол, можно точки в количестве записываемых краевых условий расставить на хорде дуги, а затем спроектировать их на нее ортогонально к хорде.
Здесь выделено четыре типа ячеек, каждая из которых имеет четыре (а, б), три (в), две (г) и одну (д) соседнюю ячейку соответственно (см. рис. 2). На рис. 3 изображен фрагмент расчетной области.
Объединяя уравнения коллокации, условия согласования и краевые условия (в случае граничных ячеек), в каждой ячейке относительно неизвестных коэффициентов Ь^ получим переопределенную СЛАУ вида
^ внь.
гз = I = 1,..., 12,
г=1
которую условно можно назвать "локальной". Для того чтобы определить, что понимается под решением этой системы, рассмотрим функционал
ф = Е (Е- я)
1=1 \1=1 )
ИЩ] —
Он представляет собой сумму квадратов невязок всех уравнений в ячейке на приближенном решении задачи. Решение системы (8) на каждой итерации находится из условия минимума функционала Ф по коэффициентам Ь^, г = 1,..., 6 [5]. Система уравнений, полученная объединением уравнений во всех ячейках расчетной области (глобальная СЛАУ), решается в итерационном процессе Гаусса — Зейделя. В ней одна "глобальная итерация" состоит из последовательного решения локальных СЛАУ (8) во всех ячейках области. При построении решения в каждой ячейке матрица системы из двенадцати уравнений приводится к верхнетреугольному виду ортогональным методом Гивен-са [14]. В численных экспериментах в качестве начального приближения здесь взяты Ь^ = 0.4. Итерационный процесс продолжается до тех пор, пока не выполнится условие
тах \Ьг]п+1 - Ъгзп1 <е, (10)
гз
где Ь^ — г-й коэффициент полинома (г = 1,..., 6), аппроксимирующего решение в ячейке с номером ] на п-й итерации. Величина е — заданная константа, называемая псевдопогрешностью решения. Ее значение выбирается таким, чтобы погрешность решения глобальной СЛАУ была существенно меньше погрешности аппроксимации исходной задачи.
Если в расчетной сетке имеется или появилась после измельчения шагов сетки нерегулярная ячейка, у которой хотя бы один из ее линейных размеров (высота или ширина) мал по сравнению с одним из линейных размеров соседней ячейки, то глобальная СЛАУ задачи может быть плохо обусловлена. Такие ячейки (в том числе вытянутые) для краткости условно назовем "малыми". Они могут быть причиной низкой точности приближенного решения задачи на содержащей их сетке. Чтобы избежать таких ситуаций, здесь используется идея присоединения этих ячеек к соседним. Если геометрический центр материнской ячейки лежит в расчетной области П, то такую ячейку назовем самостоятельной, иначе — несамостоятельной. Предлагается несамостоятельную ячейку всегда присоединять к той соседней самостоятельной ячейке, с которой она имеет наибольшей длины общую сторону при сравнении со сторонами, общими с другими соседними ячейками. Если несамостоятельная ячейка имеет одинаковую длину с несколькими общими сторонами самостоятельных соседних с ней ячеек, то ее присоединяем к ячейке с наименьшими индексами [г,]). Через эту сторону в присоединенную ячейку(ки) продолжается решение из самостоятельной ячейки, к которой ее (их) присоединили. Нерегулярная ячейка, имеющая только одну прямолинейную сторону, всегда присоединяется к соседней ячейке. При этом, согласно сформулированному правилу, к самостоятельной ячейке могут быть присоединены одна, две или три малые (несамостоятельные) ячейки. В этом случае внешняя сторона объединенной ячейки состоит из внешних сторон всех объединенных в ней ячеек. На этой дуге приблизительно равномерно расставляется столько точек для записи краевых условий, чтобы в сумме число точек записи уравнений коллокации, условий согласования и краевых условий в объединенной ячейке было равно 12. Эти краевые условия включаются в переопределенную систему уравнений этой ячейки. На рис. 3 изображена часть расчетной области, покрытая сеткой с прямоугольными ячейками. В получившихся расчетных ячейках по сформулированным правилам расставлены точки записи уравнений коллокации, условий согласования и краевых условий.
1.2. Численные эксперименты и результаты расчетов
В численных экспериментах на сходимость решения на последовательности сеток при измельчении шагов сетки вдвое рассматривалась задача (1) в областях, границы которых составлены из частей окружностей разного радиуса при различном взаимном расположении (см. рис. 1, справа). Такие области здесь взяты в качестве примера. Реализованные алгоритмы без изменений могут применяться и в других неканонических областях, в частности, границы которых заданы кусками сплайнов. Они также могут быть состыкованными с кривыми других видов.
В качестве тестовых решений брались аналитические функции, подстановкой которых в уравнения подбиралась их правая часть (в том числе в разд. 2) для того, чтобы эти функции удовлетворяли рассматриваемым дифференциальным уравнениям. Краевые условия брались из значений тестовых функций и их производных на 8П.
В табл. 1-4 приведены значения погрешности решения в равномерной норме
II и - ин ||с= тах тах 1и[хи,Х23) - иы[хи,х23)1, (11)
где N — количество расчетных ячеек в области; Ьг — количество равномерно распределенных точек [х\8,х2з) для подсчета погрешности в г-й ячейке; и — точное решение
задачи (1); и^ — приближенное решение в области П; и^ — приближенное решение в г-й ячейке. В численных экспериментах в одном случае взяты следующие значения параметров: Ь^ = 100 в регулярной ячейке, е = 10-10, к1 = 1, к2 = 5, к3 = 1.5. При этом погрешность решения незначительно отличалась от погрешности в случае, когда к1 = 1, к2 = 1, к3 = 1.
Таблица 1. Результаты численного эксперимента с тестовым решением и(х1, х2) = еХ1 + ех2 задачи (1) в расчетной области, составленной из кругов (ж1 — 0.5)2 + (ж2 — 0.5)2 = 0.52, (ж1 — 0.25)2 + (Ж2 — 0.5)2 = 0.552, (ж1 — 0.75)2 + (ж2 — 0.5)2 = 0.552
Размер сетки Погрешность Количество итераций
16 16 5.953 ■ 10-5 157
32 32 1.493 ■ 10-5 447
64 64 3.742 ■ 10-6 1778
128 > 128 9.370 ■ 10-7 6793
Таблица 2. Результаты численного эксперимента с тестовым решением и(х1, Ж2) = ео8(ж1)+ео8(ж2) задачи (1) в расчетной области, составленной из кругов (ж1—0.5)2+(ж2—0.5)2 = 0.52, (Ж1 — 0.25)2 + (Ж2 — 0.5)2 = 0.52, (ж1 — 0.75)2 + (ж2 — 0.5)2 = 0.52
Размер сетки Погрешность Количество итераций
32 > 32 7.466 ■ 10-6 438
64 > 64 1.866 ■ 10-6 1745
128 > 128 4.672 ■ 10-7 6663
256 > 256 1.179 ■ 10-7 25158
512 > 512 2.920 ■ 10-8 94 428
Таблица 3. Результаты численного эксперимента с тестовым решением и(х1, Ж2) = ео8(ж1)+еов(ж2) задачи (1) в расчетной области, составленной из кругов (ж1—0.5)2+(ж2—0.5)2 = 0.52, (Ж1 — 0.25)2 + (Ж2 — 0.5)2 = 0.52, (ж1 — 0.75)2 + (ж2 — 0.5)2 = 0.52
Размер сетки Погрешность Количество итераций
32 > 32 7.54 ■ 10-6 170
64 > 64 1.87 ■ 10-6 628
128 х 128 4.68 ■ 10-7 2609
256 > 256 1.18 ■ 10-7 7960
512 512 2.9^ 10-8 21788
Таблица 4. Результат численного эксперимента с тестовым решением и(ж1, Ж2) = еХ1 + еХ2 задачи (1) в расчетной области, составленной из кругов (ж1 — 0.5)2 + (ж2 — 0.5)2 = 0.52, (ж1 — 0.25)2 + (Ж2 — 0.5)2 = 0.52, (Ж1 — 0.75)2 + (ж2 — 0.5)2 = 0.52, е = 0.5 ■ 10-13
Размер сетки Погрешность Количество итераций
4 х 4 7.212 ■ 10-7 74
8 х 8 3.339 ■ 10-8 81
16 х 16 2.101 ■ 10-9 144
32 х 32 1.264 ■ 10-10 212
1.3. Ускорение итерационного процесса решения СЛАУ
Для уменьшения времени построения приближенного решения использован вариант известного, основанного на подпространствах Крылова [14,17,18] метода ускорения сходимости итерационного процесса решения СЛАУ. Результаты такого численного эксперимента представлены в табл. 3.
Приведенные в табл. 1-3 результаты показывают, что численное решение задачи (1) на последовательности сеток сходится и с высокой точностью совпадает с точным аналитическим решением. При этом скорость сходимости итерационного процесса с использованием подпространств Крылова в несколько раз выше, чем в случае, когда ускорение не применялось.
1.4. Вариант метода КНН с базисом в пространстве полиномов четвертой степени
Рассмотрим вариант метода КНН с базисом в пространстве полиномов четвертой степени. В каждой ]-й ячейке решение ищем в виде линейной комбинации базисных элементов — мономов от двух переменных:
15
и(у1, У2) = ^ & = ^ + ^У1 + У2 + Ц-У12 + Ь5эУ2 +
г=1
+b6j У1У2 + Ьтз У12У2 + Ь8з У1У22 + Ц' У12У22 + Ь^ У13 + Ьщ У23 + Ъ12э У\У2 +
+Ь1Ы У1У23 + Ьщ У14 + Ь^ У24. (12)
В данной работе в каждой ячейке 16 уравнений коллокации записывались в 16 равномерно распределенных точках внутри ячейки. На каждой внутренней прямолинейной стороне ячейки равномерно распределялись четыре точки записи условий согласования. На внешней стороне каждой расчетной граничной ячейки (в том числе и объединенной), как и выше, равномерно распределялись точки записи краевых условий. При этом они брались в таком количестве, чтобы число уравнений, определяющих решение в ячейке, было 32. Результаты численного эксперимента с использованием подпространств Крылова представлены в табл. 4. Здесь на сетке 32 х 32 погрешность решения равна 1.264 • 10-10, хотя граница области криволинейная.
2. Бигармоническое уравнение 2.1. Постановка задачи и описание метода
Пусть требуется определить напряженно-деформированное состояние изотропной пластины с криволинейной границей, закрепленной на краях и находящейся под действием поперечных нагрузок произвольного вида. Пластина предполагается тонкой, что позволяет свести трехмерную задачу теории упругости к двумерной задаче теории тонких пластин. Следуя Тимошенко [19], разрешающую систему уравнений запишем в перемещениях, определив которые, можно будет полностью восстановить НДС пластины. В этом случае величины прогибов пластины определяются из решения бигармониче-ского уравнения
—^Т" + 2 аад + = , (Х1,Х2) е а (13)
Здесь т(х1,х2) — прогиб срединной поверхности пластины; д(х1 ,х2) — поперечная нагрузка; И = ЕК3/12(1 — и2) — жесткость пластины при изгибе; К — толщина пластины, Е,и — модуль Юнга и коэффициент Пуассона изотропного материала пластины; П — проекция срединной поверхности пластины на плоскость (х1,х2). Введем обозначения а, Ь для длины и ширины прямоугольника, включающего в себя область П. На каждом куске границы (края) пластины может быть задано одно из следующих условий закрепления:
дт
т = 0, -7— = 0 — защемленный край; (14)
дп
т = 0, Мпт = 0 — шарнирно закрепленный край; (15)
Мпт = 0, Упт = 0 — свободный край. (16)
Здесь и далее п — внешняя нормаль к границе области дП; Мп — дифференциальный оператор второго порядка; Уп — дифференциальный оператор третьего порядка. Дифференциальные операторы Мп, Уп взяты из [19]. Условие (16) может быть задано только на некоторой части границы, а на оставшейся ее части должны быть заданы другие условия, иначе задача некорректна.
Область П покрывается расчетной сеткой с прямоугольными ячейками П (г=1,.., N, N — число расчетных ячеек сетки). Далее к поставленной задаче применим алгоритм, аналогичный описанному выше с использованием пространства полиномов четвертой степени. Отличия будут в записи уравнений дискретной задачи.
В каждой ячейке П делается переход к локальным переменным. Тогда на к-й итерации для определения решения ц в П выписываются:
• уравнения коллокации — требования удовлетворения приближенного решения уравнению (13) в заданных точках коллокации
§^ + §Щ = М, (,„Ы 6 П.; (17)
К ду4 ду2ду2 К2 ду4 И
• условия согласования — условия склейки решения ц в некоторых точках на общей стороне с соседней ячейкой с решением в ней та(у, где Кп = К1, если направление п совпадает с направлением оси у1, и Кп = К2, если направление п совпадает с направлением оси 2,
1 дп? 1 дт л-
т+й к. ж = та*+р. кп -таг ■ ^ е дП-\дП- (18)
Ш + - Ш = ^ + - £ ^ (- *> ^ДдП; (19)
• краевые условия в некоторых точках на внешних границах ячеек.
Расстановка точек для записи уравнений коллокации, условий согласования и краевых условий аналогична описанной в подразд. 1.4. При этом количество всех перечисленных точек в ячейке 32, а определяющих решение уравнений 48, так как в каждой точке для записи условий согласования или краевых условий записывается по два уравнения.
2.2. Численные эксперименты с пластиной под специальной нагрузкой
Рассмотрим шарнирно закрепленную на краях прямоугольную пластину размера а х Ь, находящуюся под действием нагрузки д (Па) вида
, = 105 8Ш (^) 8,п (^). (20)
В этом случае задача имеет точное аналитическое решение [19]
да4Ь4
т(х1,х2) =-к. (21)
1 2 (а2 + Ь2)2 К 1
Пусть проекция срединной поверхности пластины на плоскость (х1, х2) имеет форму области П, составленной из частей кругов (х1 — 0.5)2 + (х2 — 0.5)2 = 0.52, (х1 — 0.25)2 + (х2 — 0.5)2 = 0.52, (ж1 — 0.75)2 + (х2 — 0.5)2 = 0.52, как показано на рис. 1 справа. Для того чтобы воспользоваться приведенным аналитическим решением в качестве теста, в численном эксперименте на границе 8П зададим точные значения функции (21) и изгибающего момента Мга^. Ранее в [11] методом КНН была решена задача для прямоугольной пластины. В расчетах использовались следующие параметры: а =10 м, Ь = 9.984234м, к = 0.1 м, Е = 200 ГПа, и = 0.28, е = 10-10. Здесь значения параметров р1 и р2 в условиях согласования взяты из работы [11].
Для наблюдения сходимости приближенного решения проведены численные эксперименты на последовательности сеток с уменьшением шага сетки вдвое. Для расчета погрешности в каждой расчетной ячейке равномерно в каждом направлении распределялось 100 точек. Относительная погрешность численного решения во всей области вычислялась в этих точках по формуле
_ шахХ1,Х2 1т(х1,х2) — и)н(жь^) /00л ^Мг,М2 = -1—/--, (22)
шахЖ1,Ж2 1т(х1,х2)1
где М1 х М2 — размер сетки.
В табл. 5 приведены результаты численного эксперимента с использованием подпространств Крылова. Из нее видно, что сходимость решения начинается уже на достаточно грубых сетках. Эти результаты по точности незначительно отличаются от результатов работы [11]. На более мелкой сетке 128 х 128 в погрешности решения возрастает доля ошибок округления, что неизбежно сказывается на точности вычисления решения. На рис. 4 с использованием для наглядности различных масштабов численных значений величин прогиба, координат х1 и х2 показана форма прогиба рассматриваемой пластины.
Таблица 5. Результат численного эксперимента
Размер сетки Погрешность Количество итераций
8 х 8 2.79 • 10-2 181
16 х 16 4.30 • 10-3 448
32 х 32 3.96 • 10-4 620
64 х 64 3.27 • 10-5 2117
128 х 128 1.75 • 10-5 9278
О
Рис. 4. Величина прогиба и> деформированной пластины
В настоящей статье рассмотрены также два других подхода к построению вариантов метода КНН, которые алгоритмически проще в формулировке и реализации. В первом из них вытянутые четырехугольные ячейки остаются самостоятельными, даже если центр материнской ячейки лежит вне расчетной области П. На сетке 64 х 64 (в рамках аналогичного численного эксперимента, приведенного выше) значение погрешности равно 5.56 • 10-5. Во втором случае введем критерий недопустимой малости ячейки. Пусть г\ и г2 — длины сторон граничной треугольной ячейки. Если Тг < С(2^)2 (г = 1,2), то эта ячейка считается недопустимо малой и теряет самостоятельность. В ней не строится отдельный кусок решения. На ее внешней стороне приближенно равномерно расставляются восемь точек записи краевых условий, пронумерованные в определенной последовательности, например в направлении обхода . Далее четыре из записанных краевых условий с номерами 1, 3, 5 и 7 присоединяются к переопределенной системе уравнений той соседней ячейки, к которой точка с номером 1 ближе, а уравнения с номерами 2, 4, 6 и 8 присоединяются к системе уравнений второй соседней ячейки (к которой точка номер 8 ближе). Константа С зависит от конкретной задачи и подбирается в результате численных экспериментов. Нерегулярная ячейка, имеющая только одну прямолинейную сторону, присоединяется к соседней ячейке. На ее внешней стороне равномерно расставляются четыре точки записи краевых условий, которые включаются в переопределенную систему объединенной ячейки. На сетке 64 х 64 (в рамках аналогичного численного эксперимента) значение погрешности равно 5.73 • 10-5. Видно, что лучший результат получен (см. табл. 5) тогда, когда несамостоятельные ячейки (потерявшие центр материнской ячейки) присоединяются к соседним самостоятельным ячейкам.
Все приведенные расчеты проводились на одном вычислительном узле персонального компьютера, имеющего частоту 2.80 ГГц.
Замечания
1. В случае наличия на границе области точек ее излома желательно изначально построить сетку так, чтобы они оказались на сторонах ячеек сетки. При решении задачи на последовательности сеток с кратным мельчением ячеек точки излома границы продолжат оставаться на их сторонах. В ином случае можно ожидать большую погрешность в решении задачи.
2. Не составляет труда область решения задачи с криволинейной границей накрывать сеткой из нерегулярных треугольных ячеек и сформулировать вариант метода КНН с ее применением [9]. Однако, когда сетка прямоугольная в неканонической области, как в приведенном здесь исследовании, алгоритм поиска соседних ячеек и их перебора в области значительно проще, чем в случае применения треугольных сеток. В первом случае затрачивается меньше операций на компьютере при решении задачи и программа работает быстрее. По этой причине (и не только) преимущество использованной здесь сетки перед сетками с треугольными ячейками возрастает при решении задач с повышенной точностью в неканонических областях на адаптивных сетках [10].
3. При отработке алгоритма программы, реализующей предложенный подход, для рассматриваемых здесь уравнений были построены методом КНН решения краевых задач для PDE с повышенной точностью в различных трапециевидных областях. Это также использовано для расчета НДС трапециевидных пластин, что представляет практический интерес. Изложение этих расчетов из-за ограничений на объем статьи не приводим.
4. Здесь затруднительно указать ограничения на вид неканонических областей и весь класс решаемых в них задач, для которых будет успешно применим предложенный подход. Однако очевидно, что этот класс достаточно широк.
Таким образом, предложен и реализован подход к построению вариантов метода КНН решения краевых задач для PDE в областях с криволинейной границей с повышенным порядком точности. Приближенное решение задач искалось в пространстве полиномов второй и четвертой степеней. Для апробации метода проведен ряд численных экспериментов. Показано, что приближенные решения, полученные методом КНН, с высокой точностью совпадают с аналитическими решениями тестовых задач.
Благодарности. Работа выполнена при поддержке РФФИ (грант № 13-01-00227-a). Список литературы / References
[1] Липавский М.В., Толстых А.И. Об одной мультиоператорной схеме десятого порядка и ее применении в прямом численном моделировании // Журн. вычисл. математики и матем. физики. 2013. Т. 53, № 44. C. 600-614.
Lipavskii, M.V., Tolstykh, A.I. Tenth-order accurate multioperator scheme and its application in direct numerical simulation // Comput. Math. and Math. Phys. 2013. Vol. 53, № 4. P. 455-468.
[2] Ш^апеев А.В., ШШапеев В.П. Разностные схемы повышенной точности для решения эллиптических уравнений в области с криволинейной границей // Журн. вычисл. математики и матем. физики. 2000. Т. 40, № 2. С. 223-232.
Shapeev, A.V., Shapeev, V.P. High-order accurate difference schemes for elliptic equations in a domain with a curvilinear boundary // Comput. Math. and Math. Phys. 2000. Vol. 40, No. 2. P. 213-221.
[3] Botella, O., Peyret, R. Benchmark spectral results on the lid-driven cavity flow // Computers & Fluids. 1998. Vol. 27, No. 4. P. 421-433.
[4] Shapeev, A.V., Lin, P. An asymptotic fitting finite element method with exponential mesh refinement for accurate computation of corner eddies in viscous flows // SIAM J. Sci. Comput. 2009. Vol. 31, No. 3. P. 1874-1900.
[5] Слепцов А.Г. Коллокационно-сеточное решение эллиптических краевых задач / / Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1991. Т. 5(22), № 2. С. 101-126.
Sleptsov, A.G. Collocation-grid solution of elliptic boundary-value problems // Modelirovanie v Mekhanike. 1991. Vol. 5(22), No. 2. P. 101-126. (In Russ.)
[6] Исаев В.И., ШШапеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье — Стокса // Журн. вычисл. математики и матем. физики. 2010. Т. 50, № 10. 1758-1770.
Isaev, V.I., Shapeev, V.P. High-accuracy versions of the collocations and least squares method for the numerical solution of the Navier —Stokes equations // Comput. Math. and Math. Phys. 2010. Vol. 50, № 10. P. 1670-1681.
[7] Gibou, F., Fedkiw, R. A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem // J. of Chemical Physics. 2005. Vol. 202, No. 2. P. 577-601.
[8] Xu-Dong Liu, Fedkiw, R.P., Kang, M. A Boundary Condition Capturing Method for Poisson's Equation on Irregular Domains // J. of Chemical Physics. 2000. Vol. 160, No. 1. P. 151-178.
[9] Слепцов А.Г., Шокин Ю.И. Адаптивный проекционно-сеточный метод для эллиптических задач // Журн. вычисл. математики и матем. физики. 1997. Т. 37, № 5. С. 572-586. Sleptsov, A.G., Shokin, Yu.I. An adaptiv grid-projection method for elliptic problems // Comput. Math. and Math. Phys. 1997. Vol. 37, No. 5. P. 558-571.
[10] Беляев В.В., ШШапеев В.П. Метод коллокаций и наименьших квадратов на адаптивных сетках в области с криволинейной границей // Вычисл. технологии. 2000. Т. 5, № 4. С. 12-21.
Beljaev, V.V., Shapeev, V.P. The collocation and least squares method on adaptiv grids in a domain with a curvilinear boundary // Comput. Technologies. 2000. Vol. 5, No. 4. P. 12-21. (In Russ.)
[11] Голушко С.К., Идимешев С.В., ШШапеев В.П. Метод коллокаций и наименьших невязок в приложении к задачам механики изотропных пластин // Вычисл. технологии. 2013. Т. 18, № 6. С. 31-43.
Golushko S.K., Idimeshev S.V., Shapeev V.P. Application of collocations and least residuals method to problems of the isotropic plates theory // Comput. Technologies. 2013. Vol. 18, No. 6. P. 31-43. (In Russ.)
[12] Семин Л.Г., Слепцов А.Г., ШШапеев В.П. Метод коллокаций — наименьших квадратов для уравнений Стокса // Вычисл. технологии. 1996. Т. 1, № 2. С. 90-98.
Semin, L.G., Sleptsov, A.G., Shapeev, V.P. Method of collocations — least squares for Stokes equations // Comput. Technologies. 1996. Vol. 1, No. 2. P. 90-98. (In Russ.)
[13] Исаев В.И., ШШапеев В.П. Развитие метода коллокаций и наименьших квадратов // Тр. ИММ УрО РАН. 2008. Т. 14, № 1. С. 41-60.
Isaev, V.I., Shapeev, V.P. Development of the collocations and least squares method // Trudy Instituta Matematiki i Mekhaniki UrO RAN. 2008. Vol. 14, No. 1. P. 41-60. (In Russ.)
[14] Шапеев В.П., Ворожцов Е.В., Исаев В.И., Идимешев С.В. Метод коллокаций и наименьших невязок для трехмерных уравнений Навье — Стокса // Вычисл. методы и программирование. Разд. 1. Вычисл. методы и приложения. 2013. Т. 14, № 1. С. 306-322. Shapeev, V.P., Vorozhtsov, E.V., Isaev, V.I., Idimeshev, S.V. The method of collocations and least residuals for three-dimensional Navier —Stokes equations // Numerical Methods and Programming. Section 1. Numerical Methods and Applications. 2013. Vol. 14, No. 1. P. 306-322. (In Russ.)
[15] Shapeev, V. Collocation and least residuals method and its applications. EPJ Web of Conferences 108, 01009 (2016). DOI: 10.1051/epjconf/201610801009.
[16] Исаев В.И., Ш^апеев В.П., Еремин С.А. Исследование свойств метода коллокации и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье —Стокса // Вычисл. технологии. 2007. Т. 12, № 3. С. 53-70.
Isaev, V.I., Shapeev, V.P., Eremin, S.A. An investigation of the collocation and the least squares method for solution of boundary value problems for the Navier —Stokes and Poisson equations // Comput. Technologies. 2007. Vol. 12, No. 3. P. 53-70. (In Russ.)
[17] Слепцов А.Г. Об ускорении сходимости линейных итераций II // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1989. Т. 3(20), № 5. C. 118-125. Sleptsov, A.G. On convergence acceleration of linear iterations. II // Modelirovanie v Mekhanike. 1989. Vol. 3(20), No. 5. P. 118-125. (In Russ.)
[18] Saad, Y. Numerical methods for large eigenvalue problems. Manchester Univ. Press, 1991. 358 p.
[19] Тимошенко С.П., Войновский-Кригер С. Пластины и оболочки. М.: Физматгиз, 1963. 536 c.
Timoshenko, S., Woinowsky-Krieger, S. Theory of plates and shells. Moscow: Fizmatgiz, 1963. 536 p. (In Russ.)
Поступила в 'редакцию 5 мая 2016 г., с доработки —10 июня 2016 г.
Versions of high order accuracy collocation and least residuals method in the domain with a curvilinear boundary
Shapeev, Vasily P.1,2, Belyaev, Vasily A.1,2'* 1 Novosibirsk State University, Novosibirsk, 630090, Russia
2Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Belyaev, Vasily A., e-mail: belyaevasily@mail.ru
The paper proposed and implemented new versions of the method of collocation and least residuals (CLR) for the numerical solution of boundary value problems for PDE in domains with a curvilinear boundary. Their implementation and numerical experiments are performed on the examples of the biharmonic and Poisson equations. The solution of the biharmonic equation is used for simulation of the stress-strain state of an isotropic plate under the action of the transverse load. In the present study we apply the idea of using parts of the cells of a regular grid outside the domain, which are cut off by the boundary for the constructing the CLR methods. It is assumed that the solution has no singularities on the boundary and in a certain small neighborhood of it. The differential equation for the problem is true not only in the computational domain, but also in a small neighborhood of its boundary. The original domain with a curvilinear boundary is inscribed into a rectangle and then covered by a regular grid with rectangular cells, so that the points of fracture of boundary will be on the sides of the cells. We use the idea of joining the "small" irregular cells to the adjacent cells in order to reduce the condition number for the global system of linear algebraic equations.
© ICT SB RAS, 2016
110
В.П. WlaneeB, B.A. BexaeB
It is shown that the approximate solutions obtained by CLR converge with high order of accuracy thus accurately match the analytical solutions of test problems.
Keywords : collocations and least residuals method, non-canonical area, high order approximation, Poisson's equation, biharmonic equation.
Acknowledgements. The work has been supported by the RFBR (grant No. 13-0100227).
Received 5 May 2016, Received in revised form 10 June 2016
