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

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

Том 21, № 6, 2016

Численный анализ весового метода конечных элементов для задачи теории упругости с сингулярностью

В. А. Руклвишников*, С. Г. Николаев Вычислительный центр ДВО РАН, Хабаровск, Россия *Контактный e-mail: vark0102@mail.ru

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

Ключевые слова: краевая задача с сингулярностью, весовой метод конечных элементов, Rv-обобщенное решение.

Введение

Обобщенное решение краевой задачи для эллиптического уравнения в двумерной области П с границей, содержащей входящие углы (i = 1,..., N), принадлежит классу

W2,+l3~£ (П), fa = min {faj}, где faj = ж/ji для задачи Дирихле или Неймана и faj = i=l,...,N

для смешанной краевой задачи; е — любое положительное число (см., например, [1]). В соответствии с принципом согласованных оценок [2, 3] приближенное решение по методу конечных элементов (МКЭ) или конечных разностей сходится к обобщенному решению задачи со скоростью 0(№) (fa < 1) в норме пространства (П).

Для решения задач теории упругости с сингулярностью разработаны различные численные методы, обладающие теми или иными преимуществами по сравнению с классическим методом конечных элементов, выделим некоторые из них:

— сглаженные методы конечных элементов (smoothed finite element methods) [4-6];

— разрывные схемы Бубнова — Галёркина (discontinuous Petrov — Galerkin framework) [7];

— метод граничных элементов [8, 9];

— методы с сингулярными суперэлементами (singular-super-elements) [10, 11];

— четвертьточечные конечные элементы (quarter-point finite elements) [12-14];

— методы со сгущением сеток и повышением степени аппроксимирующих полиномов [15-18];

— бессеточные (meshless) методы [19-22];

— смешанные методы, основанные на расширенном (extended finite element method) и сглаженном методах конечных элементов [23].

© ИВТ СО РАН, 2016

В [24] для создания численных методов без потери точности для краевых задач с сингулярностью решения предложено определять решение как -обобщенное. В статьях [25-29] для краевых задач с сильной сингулярностью решения, у которых интеграл Дирихле от решения расходится, установлены оценки скорости сходимости первого и второго порядков точности по шагу сетки приближенного решения по МКЭ к Д^-обобщенному решению в нормах весовых пространств С.Л. Соболева и Лебега соответственно. В работах [30 - 32] построен и исследован весовой метод конечных элементов для задач теории упругости в областях с входящим углом, скорость сходимости которого не зависела от размеров этого угла и была равна О (К).

В настоящей работе для системы уравнений Ламе с решениями двух типов (содержащим сингулярную и регулярную составляющие и только сингулярную составляющую) проведен сравнительный численный анализ весового метода конечных элементов, классического МКЭ и МКЭ с геометрическим сгущением сеток к точке особенности. Для первых двух методов нашли подтверждение теоретические оценки скорости сходимости. Кроме того, установлено, что при прочих равных вычислительных условиях на сетках большой размерности для МКЭ со сгущением сеток происходит "сбой", в то время как весовой МКЭ устойчиво позволяет находить приближенное решение с теоретической точностью. В итоге весовой метод конечных элементов при одинаковой теоретической скорости сходимости О (К) с МКЭ со сгущением сеток позволяет находить решение при размерности сеток выше на один-два порядка.

1. Основные обозначения

Пусть П = (-1,1) х (-1,1)\[0,1] х [— 1, 0] С В2 — ¿-образная область с границей дП, содержащая входящий угол величиной 3^/2 с вершиной в точке 0(0, 0), П = П и дП.

Обозначим через П' = {х Е П : (х^ + х"2)1/2 < ^ < 1} часть ^-окрестности точки (0,0), лежащую в П. Введем весовую функцию р(х), совпадающую с расстоянием до начала координат в П' и равную 8 для х Е П\П'. Через а(П, 8) обозначим весовое множество функций, удовлетворяющих условиям:

1) < с1(8 / р(х))а+к для х Е П', где к = 0,1 и с1 —положительная постоянная, не зависящая от к;

2) N1 и,а(П\П') ^ С2 > 0 с н°рмой

1/2

N^а(П) = I ^ I Р^ 2 ^

л д|Л|

где И = —т-, Л = (Л1, Л2), |А| = Л1 + Л2; Л1, Л2 — целые неотрицательные числа,

дх11 дх22

а — некоторое вещественное неотрицательное число.

Через ¿2 а(П,8) обозначим множество функций, которые подчиняются требованиям 1 и 2 с нормой

1/2

N1Ь,,а(П) = (У Р2аи2(1х \п

Множество IV!а(П,8) С а(П,8)

определим как замыкание по норме (1) множества Со(П,8) — бесконечно дифференцируемых финитных в П функций, для которых выполняются оба условия.

Будем говорить, что <р Е Ш2^(дП, 5), если существует такая функция Ф из Ш2 а(П, 5),

чтО Ф(^)|9П = <Р(Х) и \МЦ/>п,6) = Ф|вП=¥> ||Ф|Ца(п,6) .

Для соответствующих множеств вектор-функций будем использовать обозначения

W1,в(n,í), Г^П,*), Л¥ 1,в(П,£). 2. -обобщенное решение

Пусть и = (и\,и2) — вектор-функция перемещений, область П представляет собой однородное изотропное тело, а деформации малы. В П рассмотрим краевую задачу для системы уравнений Ламе с постоянными коэффициентами А и р относительно поля перемещений и:

-(2а^е(и)) + У(Аи)) = {, X Е П, (2)

щ = дг, г = 1,2, х Е дП. (3)

1

(дщ диз \ дху дхг I

Здесь е(и) — тензор деформаций с компонентами е^(и) = — I —--+

2 у

Предположим, что правые части (2), (3) удовлетворяют условиям

{ Е Ь2,/з(П,й), Цг Е (дП,5), г = 1, 2, р> 0 Введем билинейные и линейные формы:

а\(и, V) =

м . о -.дщд(р2VI) дщд(р2VI) дщ д(р2щ) дщ д(р2щ)

(А + 2Р)^----+ ----+ Л------+ р —----

ОХ\ ОХ\ ОХ2 ОХ2 ОХ2 ОХ\ ОХ\ ОХ2

«2(и, V):

дщд (р2и . дщд(р2 V2) . ,дщ д (р2и . дщ д (р2и

----+ р^----+ (А + ----+ Р^---

дх\ дх2 дх2 дх\ дх2 дх2 дх\ дх\

к(V) = у р "¡\Vxdx, 12(V) = у р "¡2Ь2<1х, п п

а(и, V) = (аг(и, v),a2(u, V)), /(V) = (к

Определение 1. Функцию и^ из множества W\v(П,$) будем называть -обобщенным решением задачи (2), (3), если почти всюду на дП она удовлетворяет краевому условию (3) и для любых V из (П,$) интегральное тождество

а (и^, V) = I (V) (5)

выполнено для любого фиксированного и, удовлетворяющего неравенству и > [5.

В работах [31, 32] для краевой задачи (2), (3) с однородными граничными условиями установлены существование и единственность Ки-обобщенного решения.

Теорема 1. Пусть выполнено условие (4). Тогда для любого и > ¡3 всегда найдется такой параметр 5, что Н»-обобщенное решение и^ задачи (2), (3) с однородными граничными условиями существует и единственно во множестве и(П,$). При этом справедлива оценка

(п) < \\ъ2ф(п), где с3 — положительная постоянная, не зависящая от ¥.

К настоящему времени построена законченная теория классических решений краевых задач с гладкими исходными данными (коэффициентами уравнения, правыми частями уравнения и граничных условий) и достаточно гладкой границей области [33-37].

На основе введения обобщенного решения проводились обширные исследования краевых задач с разрывными исходными данными и негладкой границей области в соболевских и различных весовых пространствах (см., например, [1, 38-41]). На базе метода Галёркина для нахождения приближенного обобщенного решения разработаны теории разностных схем, конечных объемов и метода конечных элементов [2, 3 42 и др.]

Будем называть краевую задачу задачей с сильной сингулярностью, если для нее нельзя определить обобщенное решение, оно не принадлежит пространству Соболева (Н1) или, другими словами, интеграл Дирихле от решения расходится. Нами в [24] для краевых задач с сильной сингулярностью предложено определять решение как Ки-обобщенное в весовом пространстве Соболева. Суть его состоит в том, что в интегральное тождество обобщенного решения вводится весовая функция в некоторой степени, равная расстоянию до точек сингулярности в окрестностях этих точек. Роль этой функции заключается в том, чтобы подавить сингулярность, возникающую в решении из-за особенностей задачи, и добиться сходимости интегралов в обеих частях интегрального тождества. Исходя из локального характера сингулярности весовая функция задается как расстояние до каждой точки сингулярности внутри окружности радиуса 8 с центром в этой точке и равной 8 вне ее. Показатель степени весовой функции в определении -обобщенного решения и весовое пространство, которому оно принадлежит, зависят от пространств, которым принадлежат исходные данные задачи, геометрических особенностей границы области (входящих углов) и смены типа граничных условий.

В [27, 28] для преобразованной системы уравнений Максвелла в области с входящим углом, у которой решение не принадлежит пространству Ш^, разработан весовой реберный метод конечных элементов на основе введения -обобщенного решения, у которого в отличие от других методов [43, 44] скорость сходимости была О (К) и не зависела от размеров сингулярности.

Предлагаемая методика введения -обобщенного решения позволяет не только эффективно находить решение для краевых задач с расходящимся интегралом Дирихле, но и для задач со слабой сингулярностью, когда решение задачи принадлежит и не принадлежит пространству Ж22.

3. Весовой метод конечных элементов

Для задачи (2), (3) построим схему весового метода конечных элементов на основе Ки-обобщенного решения. Для этого произведем квазиравномерную триангуляцию Тн области П и введем специальные весовые базисные функции.

Разбиваем П на конечное число треугольников К, называемых конечными элементами, с вершинами Рк (к = 1,... , N) — узлами триангуляции. Обозначим через Пн = и кетк ^ объединение всех элементов (К — максимальная из длин их сторон). Потребуем, чтобы разбиение удовлетворяло традиционным ограничениям, накладываемым на триангуляцию [26]. Через Р = {Рк^^ обозначим множество всех внутренних узлов триангуляции, а через Р = {Рк— множество узлов, принадлежащих дП.

Каждому узлу Рк Е Р поставим в соответствие функцию гфк вида

Фк (ж) = ри* (х)фк (х), к =1,...,п,

где фк (х) линейна на каждом конечном элементе и ^ (Р]) = , к,] = 1,... ,п, — символ Кронекера, V* — вещественное число.

Определим множество Ун как линейную оболочку, натянутую на систему базисных функций [фк}к=1, и обозначим через Vh = [Ун]2 соответствующее множество вектор-функций. Выделим в Vh подмножество

у^ = {V е Ун,уг(рк= о, г =1, 2} .

Конечно-элементная аппроксимация компонент вектора перемещений, связанная с построенной триангуляцией, имеет вид

п п

= ^ ^2к-\фк, и^ = ^ ^2к Фк, й] = р-и * !±1 ]) Сэ , ] = 1,..., 2п. к=1 к=1

Определение 2. Приближенным Н»-обобщенным решением краевой задачи (2), (3) по весовому методу конечных элементов будем называть такую функцию и^ е Ун, которая в узлах границы дО удовлетворяет краевому условию (3), для произвольных vh(x) е "Vн и V > ¡3 выполнено интегральное тождество а(и^,vh) = /(vfe), где ^ = (^ 1,^22).

В [30, 32] показано, что при использовании весового метода конечных элементов для нахождения -обобщенного решения в задачах теории упругости скорость сходимости приближенного решения к точному не зависит от размеров входящего угла на границе области и всегда остается равной О (К). В следующем разделе приведены результаты сравнительного численного анализа весового МКЭ с классическим методом конечных элементов и с МКЭ с геометрическим сгущением сеток двух типов для модельных задач (2), (3).

4. Результаты численных экспериментов

В области О рассмотрена задача Дирихле для системы уравнений Ламе (2), (3) с постоянными коэффициентами А = 3 и ^ = 5. В качестве решения этой задачи выбирались вектор-функции и = (и1,и2) двух видов.

Задача А. Компоненты решения и модельной задачи (2), (3) содержат только сингулярную составляющую

и1 = сов^) сов2(ж2)(ж2 + х2)а3051,

и2 = ео82(ж1) сов(ж2)(ж2 + х22)а3051.

Порядок сингулярности щ, и2 соответствует размерам входящего угла на границе области 7 = 3^/2 [45].

Задача Б. Решение и модельной задачи (2), (3) содержит сингулярную и регулярную составляющие, принадлежащие W2(О):

и1 = сов(ж1) сое2(х2)(х2 + ж2)а3051 + (х2 + х2),

и2 = сов2(ж1) сов(ж2)(ж2 + ж2)а3051 + (х! + ж2).

4.1. Сравнительный анализ для обобщенного и -обобщенного решений

Результаты численных экспериментов, приведенных в данном подразделе, получены с помощью программы "Проба-ГУ" [46] на регулярных сетках, построенных по следующей схеме:

1) область П разделялась на квадраты линиями, параллельными осям координат, отстоящими друг от друга на расстояние 1/Ы, где N — половина числа отрезков разбиения вдоль большей стороны;

2) каждый квадрат делился диагональю на два треугольника.

В этом случае величина шага К может быть определена как К = . Пример регу-

лярной сетки для N = 4 представлен на рис. 1, а.

Вычисления проводились для различных значений N. Оптимальные параметры 8, V, V* для Ки-обобщенного решения определялись с помощью комплекса [47]. Обобщенное решение устанавливалось интегральным тождеством (5) при и = 0.

Для найденных приближенных и2 = (и22,и22) — обобщенного и и^ = (и2г,и22) — Ки-обобщенного решений вычислялись соответственно ошибки е = (е1,е2) = (и\ — иI, и2 — и2) и еи = (еи>еи>2) = (щ — 1,и2 — 2). В табл. 1 и 2 приведены значения от-

Рис. 1. Пример регулярной (а) и сгущающихся сеток I (б) и II (в) (Ж = 4, к = 0.4)

Таблица 1. Зависимость относительной погрешности приближенного обобщенного ^ и Ки-обобщенного (^) (5 = 0.0029, и = 1.2, и * = 0.16) решений задачи А от шага сетки

2И 128 256 512 1024 2048 4096

к 1.105е-2 5.524е-3 2.762е-3 1.381е-3 6.905е-4 3.453е-4

^ 6.963е-2 1.52 4.579е-2 1.52 3.007е-2 1.52 1.972е-2 1.53 1.293е-2 1.53 8.476е-3

^ 7.011е-2 1.55 4.522е-2 1.64 2.756е-2 2.17 1.272е-2 2.21 5.745е-3 1.98 2.902е-3

Таблица 2. Зависимость относительной погрешности приближенного обобщенного ^ и Ки-обобщенного ^ (5 = 0.0029, и = 1.2, и * = 0.16) решений задачи Б от шага сетки

2И 128 256 512 1024 2048 4096

к 1.105е-2 5.524е-3 2.762е-3 1.381е-3 6.905е-4 3.453е-4

^ 2.849е-2 1.54 1.850е-2 1.53 1.205е-2 1.53 7.870е-3 1.53 5.146е-3 1.53 3.367е-3

^ 2.868е-2 1.57 1.827е-2 1.65 1.107е-2 2.16 5.117е-3 2.21 2.319е-3 1.98 1.171е-3

носительных погрешностей для обобщенного решения в норме пространства Соболева и для Я»-обобщенного решения в норме весового пространства Собо-

W1 и = 2

лева W1 и =

при различных значениях К для задач А и Б соответствен-

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

Рис. 2. Графики ^ для обобщенного (линия с квадратами) и щ для -обобщенного (линия с кругами) (5 = 0.0029, у = 1.2, у * = 0.16) решений задачи Л(а)иБ(б)в зависимости от числа разбиений 2М

Таблица 3. Количество, процентное отношение к общему числу и распределение узлов, в которых абсолютные погрешности | (г = 1, 2) нахождения компонент приближенного обобщенного решения задачи Л не меньше заданных предельных значений

2Ж

|в11

|е21

Распределение

Предельные е11 е21

значения % Кол-во % Кол-во

• > 5е—6 48.077 6 045 579 48.077 6 045 579

О > 1е—6 29.387 3 695 290 29.387 3695 290

О > 5е—7 6.724 845468 6.724 845468

О > 1е—7 9.624 1210192 9.624 1210192

О > 5е—8 2.564 322449 2.564 322 449

О > 0 3.624 455 743 3.624 455 743

4096

Таблица 4. Количество, процентное отношение к общему числу и распределение узлов, в которых абсолютные погрешности | (г = 1, 2) нахождения компонент приближенного ^-обобщенного решения задачи Л (5 = 0.0029, у = 1.2, = 0.16) не меньше заданных предельных значений

2Ж

1

21

Распределение

Предельные | К, 21

значения % Кол-во % Кол-во

• > 5е—6 0.033 4102 0.033 4102

О > 1е—6 0.764 96 075 0.764 96 075

О > 5е—7 2.457 308985 2.457 308985

О > 1е-7 21.704 2 729 186 21.704 2 729 186

О > 5е—8 12.589 1582 976 12.589 1 582 974

О > 0 62.454 7853397 62.454 7853 399

4096

Таблица 5. Количество, процентное отношение к общему числу и распределение узлов, в которых абсолютные погрешности |е^| (г = 1, 2) нахождения компонент приближенного обобщенного решения задачи Б не меньше заданных предельных значений

2N

N

|е21

Распределение

Предельные е1| С 21

значения % Кол-во % Кол-во

• > 5е—6 48.078 6 045 622 48.078 6 045 622

О > 1е—6 29.387 3 695 278 29.387 3 695 278

О > 5е—7 6.724 845 466 6.724 845 466

О > 1е—7 9.624 1210158 9.624 1210159

О > 5е—8 2.564 322 439 2.564 322 438

о > 0 3.624 455 758 3.624 455 758

4096

Таблица 6. Количество, процентное отношение к общему числу и распределение узлов, в которых абсолютные погрешности |е^| (г = 1,2) нахождения компонент приближенного Ки-обобщенного решения задачи Б (5 = 0.0029, и = 1.2, и* = 0.16) не меньше заданных предельных значений

к

Распределение

Предельные |е^д | |е^,21

значения % Кол-во % Кол-во

• > 5е—6 0.033 4108 0.033 4108

О > 1е—6 0.771 96 899 0.771 96 899

О > 5е—7 2.481 311996 2.481 311996

О > 1е—7 21.789 2 739 862 21.789 2 739 863

О > 5е—8 12.588 1 582 876 12.588 1 582 876

о > 0 62.339 7838 980 62.339 7838 979

4096

графики скорости сходимости обобщенного и Ки-обобщенного решений соответствующих задач, изображенные в логарифмической шкале. Пунктирная линия на рисунках соответствует сходимости со скоростью О (К).

В табл. 3-6 для заданных предельных погрешностей приведено количество узлов, в которых |е1|, |е2|, |е^д|, 1е„,21 принадлежат заданному диапазону, процентное отношение этого количества к общему числу узлов, а также представлены изображения распределения абсолютных погрешностей в области П.

4.2. Метод конечных элементов со сгущением сетки. Сравнительный анализ

Приведем результаты анализа погрешности нахождения обобщенного решения задач А и Б с помощью МКЭ на сгущающихся сетках двух видов (подробнее о сетках со сгущением см., например, в [16, 48, 49]).

Сетка I. Данное разбиение построено по следующей схеме:

1. В области П для заданного N строилась регулярная сетка, описанная в под-разд. 4.1.

2. Для каждого узла определялся уровень I = шах(|Ж — [(хг + 1)Ж]|), которому он

1=1,2

принадлежит. Здесь Хг (г = 1, 2) — исходные координаты узла на регулярной сетке, [•] — целая часть числа.

3. Новые координаты узлов сетки со сгущением вычислялись по формуле ([(хг + ад — N)1-1(1/Ы)1/к (г = 1, 2).

Сетка II. Построение этой сетки отличалось от процесса, описанного выше, спосо-

2

бом определения уровня. Здесь I = ^ — [(^ + ]|. При этом новые координаты

г=1

определялись только для узлов с I < N, а для остальных оставались прежними. Пример сеток I и II приведен на рис. 1, б и в соответственно.

Решение по методу конечных элементов с указанными сгущениями сеток при значениях параметра к меньше порядка сингулярности решения сходится с первым порядком по шагу сетки [16, 48].

Вычисления проводились для различных значений N и к. Для найденных приближенных обобщенных решений и2, и^, полученных соответственно на сетках I и II, в каж-

Таблица 7. Зависимость относительной погрешности приближенного обобщенного решения задачи А на сетке I от степени ее измельчения при различных к

2И 128 256 512 1024 2048 4096

К = 0.3

т 2.659е-2 2.00 1.332е-2 2.00 6.675е-3 1.91 3.501е-3 0.75 4.650е-3 0.27 1.741е-2

к 0.062263 1.979 0.031459 1.99 0.015812 1.995 0.007926 1.997 0.003968 1.999 0.001985

К = 0.4

т 2.111е-2 2.00 1.057е-2 1.99 5.302е-3 1.78 2.971е-3 0.53 5.559е-3 0.26 2.154е-2

к 0.044928 1.986 0.02262 1.993 0.011349 1.997 0.005684 1.998 0.002845 1.999 0.001423

К = 0.5

т 1.990е-2 1.99 1.001е-2 1.99 5.038е-3 1.71 2.940е-3 0.46 6.401е-3 0.25 2.513е-2

к 0.034611 1.99 0.017387 1.995 0.008714 1.998 0.004362 1.999 0.0021823 1.999 0.001092

К = 0.6

т 2.315е-2 1.92 1.204е-2 1.93 6.254е-3 1.70 3.678е-3 0.50 7.292е-3 0.26 2.818е-2

к 0.030169 1.993 0.015135 1.997 0.007580 1.998 0.003793 1.999 0.0018973 1.9996 0.0009489

Таблица 8. Зависимость относительной погрешности приближенного обобщенного решения задачи А на сетке II от степени ее измельчения при различных к

2Ы 128 256 512 1024 2048 4096

К = 0.3

Пи 2.392е-2 2.00 1.196е-2 2.00 5.982е-3 1.99 3.012е-3 1.46 2.059е-3 0.36 5.687е-3

к 0.05114 1.982 0.025805 1.99 0.012962 1.995 0.006496 1.998 0.003252 1.999 0.001627

К = 0.4

Пи 1.974е-2 2.00 9.879е-3 2.00 4.942е-3 1.97 2.511е-3 1.16 2.167е-3 0.30 7.154е-3

к 0.038606 1.988 0.019417 1.994 0.009737 1.997 0.004876 1.999 0.00244 1.999 0.001220

К = 0.5

Пи 1.954е-2 1.98 9.857е-3 1.99 4.963е-3 1.93 2.565е-3 0.94 2.726е-3 0.28 9.725е-3

к 0.031006 1.99 0.015564 1.996 0.007797 1.998 0.003902 1.999 0.001952 1.9995 0.000976

К = 0.6

Пи 2.339е-2 1.91 1.225е-2 1.92 6.386е-3 1.90 3.368е-3 1.14 2.966е-3 0.31 9.712е-3

к 0.025906 1.995 0.012987 1.997 0.006502 1.999 0.003253 1.999 0.001627 1.9997 0.000814

дом узле вычислялись ошибки ег = и — и^ и ец = и — и^. В табл. 7, 9 приведены значения относительных погрешностей для обобщенного решения задач А и Б в нор-

1 I ||е/ ||w1

ме пространства W!(П) при различных К и к, полученных на сетке I = ——-—2

V ||иМ

а в табл. 8, 10 — на сетке II щ = —^—2 . Кроме того, в таблицах представлены

1и|^1

2

отношения норм погрешностей и шагов при увеличении числа узлов в четыре раза. На рис. 3 представлены графики скорости сходимости обобщенного решения соответствующих задач для сеток I и II, изображенные в логарифмической шкале. Пунктирная линия на рисунках соответствует сходимости со скоростью О (К).

Таблица 9. Зависимость относительной погрешности приближенного обобщенного решения задачи Б на сетке I от степени ее измельчения при различных к

2И 128 256 512 1024 2048 4096

к = 0.3

-щ 9.851е-3 1.99 4.955е-3 1.97 2.510е-3 1.36 1.845е-3 0.33 5.639е-3 0.25 2.247е-2 к 0.062263 1.979 0.031459 1.99 0.015812 1.995 0.007926 1.997 0.003968 1.999 0.001985

к = 0.4

-щ 7.712е-3 1.99 3.870е-3 1.95 1.988е-3 0.98 2.034е-3 0.28 7.218е-3 0.25 2.866е-2 к 0.044928 1.986 0.02262 1.993 0.011349 1.997 0.005684 1.998 0.002845 1.999 0.001423

к = 0.5

-щ 7.625е-3 1.99 3.839е-3 1.92 1.995е-3 0.87 2.301е-3 0.27 8.676е-3 0.25 3.454е-2 к 0.034611 1.99 0.017387 1.995 0.008714 1.998 0.004362 1.999 0.002182 1.999 0.001091

к = 0.6

•щ 9.330е-3 1.92 4.849е-3 1.88 2.584е-3 0.91 2.847е-3 0.28 1.016е-2 0.25 4.001е-2 к 0.034611 1.99 0.017387 1.995 0.008714 1.998 0.004362 1.999 0.002182 1.999 0.001091

Таблица 10. Зависимость относительной погрешности приближенного обобщенного решения задачи Б на сетке II от степени ее измельчения при различных к

2И 128 256 512 1024 2048 4096

к = 0.3

гиг 5.963е-3 2.00 2.982е-3 2.00 1.492е-3 1.91 7.819е-4 0.77 1.013е-3 0.27 3.757е-3 к 0.05114 1.982 0.025805 1.99 0.012962 1.995 0.006496 1.998 0.003252 1.999 0.0016267

к = 0.4

^гг 6.349е-3 2.00 3.178е-3 2.00 1.591е-3 1.87 8.490е-4 0.67 1.263е-3 0.26 4.805е-3 к 0.038606 1.988 0.019417 1.994 0.009737 1.997 0.004876 1.999 0.00244 1.999 0.0012203

к = 0.5

^гг 7.441е-3 1.98 3.756е-3 1.98 1.894е-3 1.83 1.037е-3 0.60 1.717е-3 0.26 6.606е-3 к 0.031006 1.99 0.015564 1.996 0.007797 1.998 0.003902 1.999 0.001952 1.9995 0.0009763

к = 0.6

^гг 9.574е-3 1.91 5.000е-3 1.92 2.602е-3 1.85 1.409е-3 0.78 1.804е-3 0.27 6.660е-3 к 0.025906 1.995 0.012987 1.997 0.006502 1.999 0.003253 1.999 0.001627 1.9997 0.00081366

Рис. 3. Графики щ на сетке I (линия с квадратами) и ^ц на сетке II (линия с кругами) для задачи А (а) иБ (б) в зависимости от числа разбиений 2 Ж; к = 0.3

Выводы

1. Приближенное Я»-обобщенное решение задачи (2), (3) сходится к точному со скоростью О (К) в норме W1 у (О,, 5) в отличие от обобщенного, скорость которого имеет порядок О (К0 61) для классического МКЭ.

2. Величина абсолютной погрешности приближенного Я„-обобщенного решения в подавляющем числе узлов сетки на один-два порядка меньше, чем для обобщенного решения.

3. При одинаковых вычислительных условиях из-за мелкости сетки около точек сингулярности невозможно использовать МКЭ со сгущением сеток при больших значениях N, в то время как весовой МКЭ позволяет корректно находить решение с большой степенью точности.

Список литературы / References

[1] Кондратьев В.А. Краевые задачи для эллиптических уравнений в областях с коническими или угловыми точками // Тр. Московского матем. общества. 1967. Т. 16. С. 209-292. Kondrat'ev, V.A. Boundary value problems for elliptic equations in domains with conical or angular points // Transactions of the Moscow Mathem. Soc. 1967. Vol. 16. P. 227-313.

[2] Сьярле Ф. Метод конечных элементов для эллиптических задач. М.: Мир, 1980. 512 с. Ciarlet, P.G. The finite element method for elliptic problems. Amsterdam: North-Holland Publishing Company, 1978. 529 p.

[3] Самарский А.А., Лазаров Р.Д., Макаров В.Л. Разностные схемы для дифференциальных уравнений с обобщенными решениями. М.: Высш. шк., 1987. 296 c. Samarskiy, A.A., Lazarov, R.D., Makarov, V.L. Finite-difference schemes for differential equations with generalized solutions. Moscow: Vysshaya Shkola, 1987. 296 p. (In Russ.)

[4] Nguyen-Thoi, T., Vu-Do, H., Rabczuk, T., Nguyen-Xuan, H. A node-based smoothed finite element method (NS-FEM) for upper bound solution to visco-elastoplastic analyses of solids using triangular and tetrahedral meshes // Comput. Methods Appl. Mech. Engrg. 2010. Vol. 199. P. 3005-3027.

[5] Liu, G.R., Nguyen-Thoi, T. Smoothed finite elements methods. Boca Raton, USA: CRC Press, 2010. 688 p.

[6] Nguyen-Xuan, H., Liu, G.R., Bordas, S., Natarajan, S., Rabczuk, T. An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order // Comput. Methods Appl. Mech. and Engineering. 2013. Vol. 253. P. 252-273.

[7] Bramwell, J., Demkowicz, L., Gopalakrishnan, J., Qiu, W. A locking-free hp DPG method for linear elasticity with symmetric stresses // Numerische Mathematik. 2012. Vol. 122, No. 4 P. 671-707.

[8] Yan, A.M., Nguyen-Dang, H. Multiple-cracked fatigue crack growth by BEM // Computational Mechanics. 1995. Vol. 16. P. 273-280.

[9] Gospodinov, G., Drakaliev, P., Kerelezova, I. A singular boundary element method for the general corner case // Anniversary Sci. Conf. 50 Years Faculty of Hydrotechics of the Univ. of Architecture, Civil Engineering and Geodesy, 6-8 Oct. 1999, Sofia, Bulgaria. 1999. P. 113-119.

[10] Yosibash, Z., Schiff, B. A superelement for two-dimensional singular boundary value problems in linear elasticity // Intern. J. of Fracture. 1993. Vol. 62. P. 325-340.

[11] Szabo, B.A., Yosibash, Z. Numerical analysis of singularities in two dimensions. Pt 2: Computation of generalized flux/stress intensity factors // Intern. J. for Numerical Methods in Engineering. 1996. Vol. 39. P. 409-434.

[12] Henshell, R.D., Shaw, K.G. Crack tip finite elements are unnecessary // Intern. J. for Numerical Methods in Engineering. 1975. Vol. 9. P. 495-507.

[13] Barsoum, R.S. Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements // Intern. J. for Numerical Methods in Engineering. 1977. Vol. 11. P. 85-98.

[14] Sayantan, P., Rao, B.N. A modified quarter point element for fracture analysis of cracks // Indian J. of Engineering & Material Sciences. 2007. Vol. 14. P. 31-38.

[15] Szabo B., Babuska I. Finite element analysis. New York: John Wiley & Sons, 1991. 368 p.

[16] Apel T., Sandig A.-M., Whiteman J.R. Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains // Mathematical Methods in the Applied Sciences. 1996. Vol. 19, No. 1. P. 63-85.

[17] Morin, P., Nochetto, R.H., Siebert, K.G. Convergence of adaptive finite element methods // SIAM Review. 2002. Vol. 44, No. 4. P. 631-658.

[18] Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z. The finite element method: its basis and fundamentals. Sixth edition. Elsevier, 2005. 733 p.

[19] Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M. Meshless methods: A review and computer implementation aspects // Mathematics and Computers in Simulation. 2008. Vol. 79. P. 763-813.

[20] Belytschko, T., Gracie, R., Ventura, G. A review of extended/generalized finite element methods for material modeling // Modelling and Simulation in Materials Science and Engineering. 2009. Vol. 17(4):043 001 (24 p).

[21] Ivannikov, V., Tiago, C., Moitinho de Almeida, J.P., Diez, P. Meshless methods in dual analysis: theoretical and implementation issues // V Intern. Conf. on Adaptive Modeling and Simulation (ADMOS 2011) / Eds. by B. Tie, D. Aubry, P. Diez. 2011. P. 291-308. Available at: http://www.lacan.upc.edu/sites/default/files/ADM0S_2011.pdf

[22] Byfut, A., Schroder, A. hp-Adaptive extended finite element method // Intern. J. for Numerical Methods in Engineering. 2012. Vol. 89. P. 1392-1418.

[23] Bordas, S., Duflot, M., Le, P. A simple error estimator for extended finite elements // Communications in Numerical Methods in Engineering. 2008. Vol. 24, No. 11 P. 961-971.

[24] Рукавишников В.А. О весовой оценке сходимости разностных схем // Докл. АН СССР. 1986. Т. 288, № 5. С. 1058-1062.

Rukavishnikov, V.A. The weight estimation for the convergence rate of finite-difference schemes // Dokl. AN SSSR. 1986. Vol. 288, No. 5. P. 1058-1062. (In Russ.)

[25] Рукавишников В.А. Задача Дирихле с несогласованным вырождением исходных данных // Докл. РАН. 1994. Т. 337, № 4. С. 447-449.

Rukavishnikov, V.A. The Dirichlet problem with the noncoordinated degeneration of the initial data // Doklady Akademii Nauk. 1994. Vol. 337, No. 4. P. 447-449. (In Russ.)

[26] Rukavishnikov, V.A., Rukavishnikova, H.I. The finite element method for boundary value problem with strong singularity // J. of Computational and Applied Mathematics. 2010. Vol. 234, No. 9. P. 2870-2882.

[27] Rukavishnikov, V.A., Mosolapov, A.O. New numerical method for solving time-harmonic Maxwell equations with strong singularity // J. of Computational Physics. 2012. Vol. 231, No. 2. P. 2438-2448.

[28] Рукавишников В.А., Мосолапов А.О. Весовой векторный метод конечных элементов для одной задачи электромагнетизма с сильной сингулярностью // Докл. Академии наук. 2013. Т. 449, № 2. С. 144-148.

Rukavishnikov, V.A., Mosolapov, A.O. Weighted edge finite element method for Maxwell's equations with strong singularity // Doklady Mathematics. 2013. Vol. 87, No. 2. P. 156-159.

[29] Rukavishnikov, V.A., Rukavishnikova, H.I. On the error estimation of the finite element method for the boundary value problems with singularity in the Lebesgue weighted space // Numerical Functional Analysis and Optimization. 2013. Vol. 34, No. 12. P. 1328-1347.

[30] Рукавишников В.А., Николаев С.Г. Весовой метод конечных элементов для задачи теории упругости с сингулярностью // Докл. Академии наук. 2013. Т. 453, № 4. С. 378-382.

Rukavishnikov, V.A., Nikolaev, S.G. Weighted finite element method for an elasticity problem with singularity // Doklady Mathematics. 2013. Vol. 88, No. 3. P. 705-709.

[31] Рукавишников В.А., Николаев С.Г. Об Rv-обобщенном решении системы Ламе с угловой сингулярностью // Докл. Академии наук. 2015. Т. 263, № 2. С. 137-139. Rukavishnikov, V.A., Nikolaev, S.G. On the Rv-generalized solution of the Lame system with corner singularity // Doklady Mathematics. 2015. Vol. 92, No. 1. P. 421-423.

[32] Rukavishnikov, V.A. Weighted FEM for two-dimensional elasticity problem with corner singularity // Numerical Mathematics and Advanced Applications ENUMATH 2015. Lecture Notes in Computational Science and Engineering. 2016. Vol. 112. P. 411-419.

[33] Douglis, A., Nirenberg, L. Interior estimates for elliptic systems of partial differential equations // Comm. Pure Appl. Math. 1955. Vol. 8. P. 503-538.

[34] Agmon, S., Douglis, A., Nirenberg, L. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions, I // Comm. Pure Appl. Math. 1959. Vol. 12. P. 623-727.

[35] Agmon, S., Douglis, A., Nirenberg, L. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions, II // Comm. Pure Appl. Math. 1964. Vol. 17. P. 35-92.

[36] Хермандер Л. Линейные дифференциальные операторы с частными производными. М.: Мир, 1965. 380 c.

Hormander, L. Linear partial differential operators / Eds B. Eckmann, B.L. van der Waerden. Grundl. d. Math. Wiss., V. 116. Berlin: Springer-Verlag, 1963. 279 p.

[37] Лопатинский Я.Б. Об одном способе приведения граничных задач для системы дифференциальных уравнений эллиптического типа к регулярным интегральным уравнениям // Украинский матем. журн. 1953. Т. 5. С. 123-151.

Lopatinskiy, Y.B. A method of reducing boundary problems for systems of differential equations of the elliptic type to regular integral equations // Ukrains'kyi Matem. Zh. 1953. Vol. 5. P. 123-151. (In Russ.)

[38] Ладыженская О.А. Краевые задачи математической физики. М.: Наука, 1973. 407 с. Ladyzhenskaya, O.A. The boundary value problems of mathematical physics // Applied Mathematical Sciences. Vol 49. New York: Springer-Verlag, 1985. 322 p.

[39] Мазья В.Г., Пламеневский Б.А. Lp-оценки решений эллиптических краевых задач в областях с ребрами. // Тр. Моск. матем. общества. 1979. Т. 37. С. 49-93.

Maz'ya, V.G., Plamenevskii, B.A. Lp-estimates of solutions of elliptic boundary value problems in domains with ribs // Trans. Moscow Math. Soc. 1980. Vol. 1. P. 49-97.

[40] Grisvard, P. Elliptic boundary problems in nonsmooth domains. Monographs and Studies in Mathematics, 24. Boston, MA: Pitman (Advanced Publishing Program), 1985. 410 p. (ISBN: 0-273-08647-2).

[41] Dauge, M. Elliptic boundary value problems in corner domains. Smoothness and Asymptotics of Solutions // Lecture Notes in Mathematics. Vol. 1341. Berlin; Heidelberg: Springer, 1988. 267 p.

[42] Самарский А.А., Андреев В.Б. Разностные методы для эллиптических уравнений. М.: Наука, 1976. 352 c.

Samarskiy, A.A., Andreev, V.B. Finite difference methods for elliptic equations. Moscow: Nauka, 1976. 352 p. (In Russ.)

[43] Assous, F., Ciarlet, P., Segre, J. Numerical solution to the time-dependent Maxwell equations in two-dimensional singular domain: The singular complement method //J. Comput. Phys. 2000. Vol. 161. P. 218-249.

[44] Costabel, M., Dauge, M. Weighted regularization of Maxwell equations in polyhedral domains // Numer. Math. 2002. Vol. 93. P. 239-277.

[45] Rossle A. Corner singularities and regularity of weak solutions for the two-dimensional Lame equations on domains with angular corners // J. of Elasticity. 2000. Vol. 60, No. 1. P. 57-75.

[46] Николаев С.Г., Рукавишников В.А. Проба IV — программа для численного решения двумерных задач теории упругости с сингулярностью: Св. 2013616248 РФ // Программы для ЭВМ. Базы данных. Топологии интегральных микросхем. 2013. Бюл. № 3(84). Nikolaev, S.G., Rukavishnikov, V.A. Proba IV — program for the numerical solution of the two-dimensional elasticity problems with singularity: Sv. 2013616248 RF // Computer programs. Data bases. IC chips topology. 2013. Byul. No. 3(84). (In Russ.)

[47] Рукавишников В.А., Маслов О.В., Мосолапов А.О., Николаев С.Г. Автоматизированный программный комплекс для поиска оптимального набора параметров весового метода конечных элементов на вычислительном кластере // Computational Nanotechnology. 2015. № 1. C. 9-19.

Rukavishnikov, V.A., Maslov, O.V., Mosolapov, A.O., Nikolaev, S.G. Automated software complex for determination of the optimal parameters set for the weighted finite element method on computer clusters // Computational Nanotechnology. 2015. No. 1. P. 9-19. (In Russ.)

[48] Оганесян Л.А., Руховец Л.А. Вариационно-разностные методы решения эллиптических уравнений. Ереван: Изд-во АН Армянской ССР, 1979. 235 c.

Oganesyan, L.A., Rukhovets, L.A. Variational-difference methods for solving elliptic equations. Erevan: Izd-vo Akad. Nauk. Arm. SSR, 1979. 235 p. (In Russ.)

[49] Raugel, G. Resolution numerique par une methode d'elements finis du probleme de Dirichlet pour le laplacien dans un polygone // Comptes Rendus de l'Academie des Sciences. Series A — Mathematics. 1978. Vol. 286. P. 791-794.

Поступила в редакцию 12 января 2016 г., с доработки —16 мая 2016 г.

Numerical analysis of the weighted finite element method for the elasticity problem with singularity

Rukavishnikov, Viktor A.*, Nikolaev, Sergey G.

Computing Center FEB RAS, Khabarovsk, 680000, Russia * Corresponding author: Rukavishnikov, Viktor A. e-mail: vark0102@mail.ru

In this paper we consider the first boundary value problem for the Lame system in the two-dimensional nonconvex domain with boundary containing one reentrant corner. For this problem we define the solution as Rv-generalized one in the weighted set. On the basis of this solution we constructed a scheme of the weighted finite element method with accuracy independent of the corner size and with O(h) convergence rate for the approximate solution to the exact one.

We carried out a comparative numerical analysis of the weighted finite element method, the classical FEM, and the FEM with meshes geometrically refined to the point of singularity. For this purpose we considered model problems for L-shaped domain. Singularity order of solutions of these problems was set as consistent with the magnitude 3^/2 of the reentrant corner.

Performed numerical analysis has demonstrated that

— the approximate Rv-generalized solution of the problem converges to the exact one at O(h) rate in the weighted Sobolev norm, while the generalized solution (for the classical FEM) has O(h0'61) rate of convergence in the Sobolev norm;

— in the overwhelming majority of grid nodes, the absolute error of the approximate Rv-generalized solution is one or two orders of magnitude lower than for the generalized solution;

— the FEM with graded meshes fails on high-dimensional grids, but the weighted FEM stably allows finding approximate solution with the high accuracy under the same computational conditions.

Keywords: boundary value problem with singularity, weighted finite element method, Rv-generalized solution.

Received 12 January 2016 Received in revised form 16 May 2016

© ICT SB RAS, 2016