HP-вариант метода коллокации и наименьших квадратов для эллиптических задач со старшими производными второго порядка на треугольных сгущающихся сетках

Численный метод коллокации и наименьших квадратов (МКНК) развивается уже более тридцати лет. За это время работы в основном сосредоточены на трех областях исследований. Первая охватывает задачи гидродинамики и тепловых процессов, включая моделирование лазерной сварки [1] и течения полимерных жидкостей [2]. Второе направление связано с анализом напряженно–деформированного состояния (НДС) пластин при изгибе [3–6]. Третья задача посвящена разработке и постоянному совершенствованию новых вариантов МКНК, их верификации на классических тестовых задачах для 2D стационарных уравнений и сравнению с другими приближенными методами. К числу таких задач обычно относятся краевые задачи для уравнения Пуассона [5, 7–9], бигармонического уравнения [3, 4, 9], уравнений На-вье– Стокса [10, 11], системы уравнений изгиба пластин в рамках теории Рейсснера – Миндлина (ТРМ) [5, 6, 9]. По результатам многих работ [4–6, 9] МКНК неоднократно сравнивался с классическими численными методами, включая метод конечных элементов [12], метод конечных разностей (МКР) [13, 14] и коллокационные бессеточные [15, 16] и сеточные методы [17], и в ряде тестовых задач демонстрировал более высокую точность или другие преимущества, в том числе более простую компьютерную реализацию.

На практике реализация МКНК весьма гибкая. В каждой ячейке задается полином с неизвестными коэффициентами, после чего последовательно формируются:

• 1)    уравнения коллокации во внутренних точках, обеспечивающие выполнение исходного дифференциального уравнения;

• 2)    условия согласования в точках сопряжения на границах соседних ячеек, необходимые для склейки полиномиального решения;

• 3)    краевые условия, записанные в точках, лежащих точно на границе области ∂Ω.

Совокупность всех уравнений формирует, как правило, переопределенную глобальную систему линейных алгебраических уравнений (СЛАУ), решение которой определяет коэффициенты аппроксимирующего полинома. Характеристики метода – точность, порядок сходимости, обусловленность и вид СЛАУ – существенным образом зависят от степени p и типа полинома [4,9], формы ячеек, значений весовых коэффициентов в линейной задаче наименьших квадратов [3, 6, 10], количества уравнений приближенной задачи [10], а также от расположения точек коллокации, согласования и записи краевых условий [9, 10]. Такая вариативность позволяет исследовать различные версии метода и сравнивать их как между собой, так и с другими численными схемами.

Настоящая статья является продолжением работы [5], в которой был разработан коллокационный метод с квадратной СЛАУ при p = 4 (h-KM ₄ ) на треугольных сгущающихся сетках для решения эллиптических задач со старшими производными второго порядка. Построение этого варианта основано на алгоритме аппроксимации, вытекающем из МКНК и использующем уравнения коллокации, согласования и краевые условия. В дальнейшем на основе аналогичного подхода был создан вариант коллокационного метода при p = 7, который также приводит к квадратной СЛАУ; соответствующая работа принята к публикации, но пока не вышла в печать. Для продолжения анализа МКНК на треугольных сгущающихся сетках применительно для данного класса задач необходимо исследовать другие случаи, приводящие к переопределенным СЛАУ, что позволит получить более точное представление о поведении метода при различных p и детальнее охарактеризовать его численные свойства.

1.    Постановка задачи

Рассмотрим сначала краевую задачу Дирихле для уравнения Пуассона в двумерной области Q С R2 с границей dQ d2u + ^u = f. (x,y) e Q.(1)

∂x2

u = g. (x.y) e dQ.(2)

где u(x.y) - искомое решение, f (x.y ) и g(x.y ) - заданные функции.

Перейдем теперь к модели изгиба однослойной пластины в рамках ТРМ. В безразмерных переменных система уравнений имеет вид [5, 9, 18, 19]:

-

1,

aTTf Cf (It ++ д~(дт +') =

2(1 + v) \Хх \Хх     J dy \dy     J J

(2 —	v)в 2 д 2 ф у + —    v 2 ) dxdy v )в 2 д 2 ф х + —    v 2 ) dxdy	β 2	∂ 2 φ x	+ 12(1 + +    " 12(1 +	∂ 2 φ x	K
12(1 (2 —		12(1 — β 2	v 2 ) dx2 ∂ 2 φ y		v ) dy2 ∂ 2 φ y	2(1 + v ) K
12(1		12(1 —	v 2 ) dy2		v ) dx2	2(1 + v )

∂w

( Я+ ^фх/ = 0.

∂x

∂w

(^ + фу) =°,

где K = 5/6 - сдвиговой коэффициент Тимошенко, v - коэффициент Пуассона, w(x.y) - прогиб пластины, ф _х (x.y) и ф _у (x.y ) - углы поворота нормали срединной поверхности, в = t/L, t - толщина пластины, L - характерный размер области Q.

Обезразмеривание выполнялось по формулам: x = x/L.y = y/L.w = w • Et ³ /(qL ⁴ ), ф х = ф х • Et ³ /(qL ³ ),ф y = ф у • Et³/(qL³ ), E - модуль упругости, q - поперечная нагрузка. В (3) и далее тильды над переменными опущены. В расчетах значения ν, t, L, q, E брались постоянными.

Краевые условия для защемленного контура имеют вид [5, 9, 18, 19]:

^w = °.   ^ф п = ^ф х ^п х + ^ф у n y = °.  ^ф s = ^ф у n x — ^ф х ^п у = °.             ⁽⁴⁾

Замечание 1. ТРМ обладает рядом особенностей по сравнению с классической теорией Кирхгофа – Лява (ТКЛ) [18]. Во-первых, в ТРМ перед старшими производными искомых функций углов поворота ф _х (x,y) и ф _у (x,y ) возникает малый параметр в . В то время как в разрешающей системе уравнений для изгиба пластин в рамках ТКЛ малый параметр отсутствует, однако в ней присутствуют производные четвертого порядка [18], что также создает дополнительные трудности при численном решении из-за ухудшения обусловленности задачи [4]. Кроме того, для системы (3), описывающей постановку ТРМ, характерен так называемый эффект сдвигового запирания: при t ^ 0 численное решение все хуже удовлетворяет требованиям, допускающим нулевые деформации поперечного сдвига [17, 20, 22]. Существует несколько подходов к устранению этого эффекта [22]. Один из них – использование смешанной постановки, в которой перерезывающие силы Q _x и Q _y вводятся как дополнительные неизвестные наряду с w , φ _x и φ _y . Такие постановки позволяют эффективно решать задачи для тонких пластин при использовании полиномиальных пространств невысокой степени [6, 17]. Другой подход заключается в применении полиномов относительно высоких степеней в исходной постановке, что также позволяет снизить проявления эффекта запирания [20]. Второе ключевое отличие ТРМ от ТКЛ состоит в том, что для сравнительно толстых пластин (L/t < 20) расчеты НДС по ТРМ дают более точные результаты [6, 18].

Целью данной статьи является разработка, реализация, верификация и анализ hp-варианта МКНК (hp-МКНК), основанного на алгоритме из [5], для численного решения краевых задач (1), (2) и (3), (4). Рассматриваются тестовые примеры с заранее известными особенностями решения, в окрестности которых строятся сгущающиеся треугольные сетки.

2.    Описание hp-МКНК

Построенный в данной работе hp-МКНК опирается на алгоритм, предложенный в статье [5], краткое описание которого дадим ниже на примере задачи Дирихле для уравнения Пуассона. Основные отличия реализованного здесь hp-МКНК от h-КМ 4 заключается в использовании полиномиальных пространств иных степеней (р = 2, 3, 5, 6, 8), в формировании переопределенных СЛАУ вместо системы с квадратными матрицами, а также в другом количестве точек записи аппроксимирующих уравнений.

Область Q покрывается треугольной сеткой, состоящей из N ceiis ячеек, сгенерированной в пакете Gmsh (рис. 1 а) и сгущающейся вблизи особенностей решения задачи (рис. 1 б). В каждой ячейке вводится своя локальная система координат:

x-x        y-y cj                       cj hj , ^2 =    hj , где (xcj ,ycj) - координаты центра j-й ячейки, j = 1,..., Ncells, а hj - радиус описанной вокруг нее окружности.

В каждой ячейке численное решение представляется в виде полинома с неизвестными коэффициентами c _{i 1 i 2 ,}j :

p p - i 1

U hj (&.&) = £g C i . i 2 ,j S* & ■                       (5)

Подставляя (5) в (1) и домножая обе части на р_с >  0, получаем уравнения коллокации:

которые записываются в N _col точках коллокации, ≪ равномерно ≫ распределенных внутри треугольной ячейки согласно алгоритму [5] (рис. 1 а).

б

Рис. 1 . Фрагмент области со сгущающейся к внутренней окружности треугольной сетки при p = 2, N col = 3, N match = N bound = 2, где о обозначает точки коллокации, х - точки согласования, □ - точки записи краевых условий (а); график решения u(x,y) = e ^{-10 (x} +y ) с большими градиентами в окрестности (0, 0) из примера 1 (б)

а

Условия согласования задаются в N _match равномерно расположенных точках согласования на каждой общей стороне соседних ячеек (рис. 1 а) и имеют вид:

P m o U hj + P m i

∂u hj

∂ n j

du h

= P m o U h + P m i d n j ,

где u h - приближенное решение в соседней ячейке, n j- - внешняя нормаль к границе j -й ячейки, P m o ,P m i >  0.

Краевые условия, домноженные на р ь >  0, записываются в граничных ячейках в N bound равномерно расположенных на части д Q, ограниченной вершинами ячейки (рис. 1 а) и имеют вид:

P b U hj = P b g.                                        (8)

Замечание 2. При решении системы (3), (4) уравнения коллокации вводятся аналогично (6), краевые условия – аналогично (8), а условия согласования (7) записываются для каждой функции w , φ x и φ y . Конкретные значения N col , N match , N bound приведены ниже.

Объединяя уравнения (6), (7) и (8) во всех ячейках области, получаем в общем случае переопределенную глобальную СЛАУ, решение которой определяет коэффициенты c i ₁ i ₂ ,j .

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

Порядок сходимости hp -МКНК определялся по следующей формуле:

E 2

^lOg 2 E

R [

^log 2 J •

,

^N cells,1

N cells,2

где Е 1 и Е ₂ - погрешности, подсчитанные на двух сетках с количеством ячеек N _ceiis, ₁ (подробная сетка) и N _cells, ₂ (грубая сетка).

Для примеров 1 и 3 рассматривалась относительная погрешность ||E _r^u||_ю , а для примера 2 - абсолютная ||EU||^, вычисленные по следующим формулам:

I E U » «

max |uex(xi,yi) - uh(xi,yi)\ i=1,...,M max |uex(xi,yi)| i=1,...,M

|Eu|U = max |uex(xi,yi) - uh(xi,yi)\, i=1,...,M где uex – точное решение задачи, uh – численное кусочно-полиномиальное решение, являющееся объединением всех решений (5) в каждой j-й ячейке, j = 1,..., Nceiis, (xi,yi), i = 1,..., M, — вершины треугольных ячеек.

Для решения разреженной глобальной переопределенной СЛАУ использовался ортогональный метод, реализованный в библиотеке SuiteSparse [23]. Параллельные вычисления выполнялись на видеокарте Nvidia GeForce RTX 3070 Ti Laptop 8 GB с помощью технологии CUDA. В SuiteSparse QR-разложение проводилось с использованием отражений Хаусхолдера, применяемых локально к фронтальным подматрицам мультифронтального алгоритма.

Показатель сгущения сетки η в примере № 1 равнялся 10, в примере № 2 – 2, 16, 64, а в примере № 3 – 8. Сгущающиеся сетки, на которых здесь было проведено сравнение с работой [5], брались из статьи опубликованной в 2024 году.

Пример 1. Задача Дирихле для уравнения Пуассона с большими градиентами. Рассмотрим (1), (2) в Q = ( — 0, 5, 0, 5) ² с тестовым решением u ^ex (x,y) = e ^-1 ooo(x +y ) , которое образует выраженный пик вблизи (0,0) и меняется слабо в остальной части Q (рис. 1 б). Подставляя u ^ex в (1), получаем правую часть f (x,y) = 4(10 ⁶ (x ² + y ² ) - 10 ³ )e ^-1ooo(x 2 + у ^{2 )} .

В дальнейшем в табл. 1 приводятся результаты, полученные с помощью разработанного здесь hp-МКНК и сравниваемые с h-KM ₄ [5]. Для каждой степени полинома p указываются размеры локальных матриц, из которых следует и размер глобальной СЛАУ. Например, для p = 2 размер переопределенной глобальной СЛАУ равен (N ceiis • 9) x (N ceiis • 6). В примере 1 и 2 были взяты следующие веса: p c = h j , p _m 0 = 1, p m i ^h j , ^p b ¹.

В hp -МКНК здесь и в примерах 2 и 3 использовались следующие параметры:

N col   ³, N match    N bound   ^{2 ;}

N col   ³, N match    N bound   ^{3 ;}

N coi ⁼ 10 ^, N match ⁼ N bound ⁼ 5;

N coi ⁼ 15 ^, N match ⁼ N bound ⁼ 5;

N col = ^{28 ,} N match = N bound = 7.

полученных данных позволяет сделать следующие выводы.

• -    p = 2:

• -    p = 3:

• -    p = 5:

• -    p = 6:

• -    p = 8:

Анализ

• 1.    Сходимость. Порядок сходимости R в среднем не ниже p -го при четных p и (p - 1)-го - для нечетных p. Подобное поведение отмечалось и ранее для МКНК [7,11] и других сеточных коллокационных методах [17]. При этом отметим, если использовать

• 2.    Эффективность высоких степеней. Использование полиномов высокой степени в задаче с гладким решением позволяет существенно повышать точность на грубых сетках. Так, при p = 8 точность порядка 10 ^-8 достигается уже на сетке N _cells = 1660

• 3.    Случаи с квадратными СЛАУ. Были рассмотрены конфигурации, при которых матрица СЛАУ становится квадратной:

  – p = 2: N _col = 3, N _match = N _bound = 1;

  – p = 3: N _col = 1, N _match = N _bound = 3;

  – p = 5: Ncol = 6, Nmatch = Nbound = 5 или Ncol = 15, Nmatch = Nbound =2;

  – p = 6: Ncol = 1, Nmatch = Nbound = 9 или Ncol = 10, Nmatch = Nbound =6;

  – p = 8: Ncol = 15, Nmatch = Nbound = 10 или Ncol = 21, Nmatch = Nbound

  Для части этих вариантов (p = 2, 6, 8 и p = 5 при N _col = 15, N _match = N _bound = 2 и N _cells = 88) число обусловленности µ в спектральной норме, вычисленное с помощью библиотеки Eigen, равнялось то . В других случаях обусловленность была существенно больше, чем при использовании переопределенных СЛАУ. Так, при p = 3 значения ^ = 2, 56 • 10 3 и ^ = 1, 69 • 10 4 на грубых сетках N cells = 88, 396, что на порядок больше по сравнению со случаем применения МКНК при параметрах из табл. 1: ^ = 1, 90 • 10 2 и ^ = 1, 03 • 10 3 . Для p = 5 ситуация еще чувствительнее: ^ = 1, 02 • 10 6 и ^ = 1, 22 • 10 6 (N coi = 6, N match = N bound = 5) по сравнению с ^ = 6,10 • 10 2 и ^ = 3, 63 • 10 3 (N co = 10, N _match = N _bound = 5) на аналогичных сетках. Погрешность также значительно возрастала: так, при p = 5 на сетке N cells = 6678 имеем Ц Е ^ ||_ю = 2, 64 • 10 ^-6 против 1, 07 • 10 ^-4 для квадратной СЛАУ.

• 4.    Прочие комбинации параметров. Для всех рассмотренных значений p изучались и другие варианты переопределенных СЛАУ. Из-за нехватки места подробности не приводятся; отметим лишь, что используемые здесь параметры являются ≪ опти-мальными ≫ в том смысле, что увеличение числа уравнений не давало заметного прироста точности, а уменьшение – приводило к ее ухудшению. Например, при N _col = 10, N match = N bound = 8 погрешность на сетке N^ is = 88 равнялась | E U | ^ = 18, 5 при ^ = 1, 30 • 10 7 , а на сетке N _cells = 396 — | E U | ^ = 1, 39 • 10 ^-1 при ^ = 8, 87 • 10 9 . В то время как для оптимального варианта МКНК погрешность была Ц Е ^ ||_ю = 6, 23 • 10 ^-2 при ^ = 1, 80 • 10 3 на сетке N _cells = 88 и | E U | ^ = 3, 75 • 10 ^-4 при ^ = 1, 01 • 10 4 - на сетке N cells = 396.

в качестве полиномиального решения прямое произведение полиномов степеней p , расставлять точки коллокации в корнях полиномов Лежандра, использовать сетки с квадратными ячейками, то R для p будет оптимальным и равняться (p + 1), начиная с определенной степени полинома, которая зависит от старшей степени производной в исходном уравнении [9].

при общем количестве степеней свободы DOF = 1660 • 45 = 74 700, тогда как для h- КМ 4 аналогичного уровня требуется N cells = 101714 и DOF = 101714 • 15 = 1 525 710. Время решения СЛАУ на компьютере 12th Gen Intel(R) Core(TM) i7-12700H 2.70 GHz, DIMM DDR5 2400 MHz 16 Gb с видеокартой NVIDIA GeForce RTX 3070 Ti Laptop, 8 GB в первом случае составило 20,141 с против 51,437 с для h-КМ 4 .

Таблица 1

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

N cells	| E U ■	R	| E U l ~	R	| E ,“ I	∞	R
	p = 2, 9 x 6		p = 3, 12 x 10		p=4	, 16 x 16 [5]
88	1,04E+0	–	5,72E-1	–	4,25E-1		–
396	3,71E-1	1,37	2,92E-2	3,96	7,35E-3		5.39
1660	7,92E-2	2,15	8,55E-3	1,71	1,58E-4		5,35
6678	2,96E-2	1,41	1,59E-3	2,42	5,85E-6		4,73
25712	8,01E-3	1,94	3,56E-4	2,22	3,54E-7		4,16
101714	2,08E-3	1,96	9,32E-5	1,95	1,82E-8		4,31
	p = 5, 25 x 21		p = 6, 30 x 24		p=	8, 49 x 45
88	2,24E-1	–	6,23E-2	–	1,24E-2		–
396	2,33E-3	6,07	3,75E-4	6,80	6,70E-5		6,94
1660	3,92E-5	5,70	2,55E-6	6,96	4,51E-8		10,2
6678	2,74E-6	3,82	1,86E-8	7,07	6,49E-11		9,40
25712	1,01E-7	4,90	2,70E-10	6,28	–		–
101714	5,81E-9	4,15	–	–	–		–

Отдельно исследовался случай p = 2 с модифицированным весовым множителем p c = 10 • h j вместо p c = h j ; заметного влияния на точность он не оказал: погрешность изменялась от 1, 01 до 1,11 • 10 - 3 на сетках с N cells = 88,..., 101 714.

Пример 2. Задача Дирихле для уравнения Пуассона с разрывом производных. Рассмотрим (1), (2) в Q = (0, 1) ² при f (x,y) = 1 и g(x,y) = 0. Известно, что в такой постановке решение u(x,y ) имеет разрывы вторых производных в угловых точках области. Точное решение представляется в виде равномерно и абсолютно сходящегося ряда [14]:

u ex& y) =

- 16 π ⁴

∞∞

ЕЕ m=1 n=1

sin((2m — 1)nx) sin((2n — 1)ny)

(2m - 1)(2n - 1)((2m - 1) ² + (2n - 1) ² )^.

В табл. 2 приведены результаты расчетов из [5] для h-КМ ₄ и реализованного hp-МКНК при параметрах сгущения сетки n = 2,16, 64. При этом для p = 8, п = 64 и N cells = 26482 вычисления не удалось выполнить из-за нехватки памяти используемого компьютера. Как видно из [14], при использовании равномерных сеток наибольшая погрешность наблюдается в окрестностях угловых точек, что наглядно проявляется на рис. 3 отмеченной работы. Поэтому применение сгущающихся сеток в данной задаче является особенно актуальным.

	Таблица 2 Результаты численных экспериментов в примере 2
N cells	IIE UIU     R     N cells     ^E U ^ ^     R     N cells     ^E U I» R
	p = 2, n = 2              p = 2, n = 16              p = 2, n = 64
296	2,59E-4   –     498     1,97E-4   –     678     1,84E-4   –
1116	6,17E-5   2,16 1418    6,40E-5   2,14 1716    6,68E-5   2,18
4176	1,83E-5    1,84 5558    1,79E-5    1,86 6730    1,82E-5    1,90
16468	5,60E-6    1,72 21849   5,10E-6    1,83 26482   5,18E-6    1,83
	p = 3, n = 2              p = 3, n = 16              p = 3, n = 64
296	5,84E-5   –     498     3,56E-6   –     678     3,64E-6   –
1116	1,26E-5   2,31   1418    4,74E-7   3,85 1716    6,86E-7   3,59
4176	3,08E-6   2,13 5558    5,21E-8   3,23 6730    3,53E-8   4,34
16468	7,64E-7   2,03 21849   1,23E-8   2,10 26482   2,77E-9   3,71
p	= 4, n = 2 [5]             p = 4, n = 16 [5]             p = 4, n = 64 [5]
296	–         –     498     5,88E-7   –     678     6,16E-7   –
1116	–         –     1418    7,72E-8   3,88 1716    7,28E-8   4,38
4176	–         –     5558    1,80E-8   2,13 6730    2,69E-9   4,68
16468	–         –     21849   4,41E-9   2,05 26482   2,73E-10 3,34
	p = 5, n = 2              p = 5, n = 16              p = 5, n = 64
296	1,27E-5   –     498     2,36E-7   –     678     2,30E-8   –
1116	2,78E-6   2,28 1418    5,24E-8   2,87 1716    3,67E-9   3,95
4176	6,75E-7   2,14 5558    1,14E-8   2,23 6730    7,28E-10 2,36
16468	1,66E-7   2,04 21849   2,72E-9   2,09 26482   1,69E-10 2,13
	p = 6, n = 2              p = 6, n = 16              p = 6, n = 64
296	8,28E-6   –     498     1,15E-7   –     678     9,43E-9   –
1116	1,82E-6   2,28 1418    3,39E-8   2,33 1716    2,17E-9   3,16
4176	4,47E-7   2,12 5558    7,52E-9   2,20 6730    4,77E-10 2,21
16468	1,10E-7   2,04 21849   1,74E-9   2,13 26482   1,18E-10 2,03
	p = 8, n = 2              p = 8, n = 16              p = 8, n = 64
296	3,60E-6   –     498     6,68E-8   –     678     4,15E-9   –
1116	7,79E-7   2,30 1418    1,47E-8   2,89 1716    9,47E-10 3,18
4176	1,89E-7   2,14 5558    3,20E-9   2,32 6730    2,04E-10 2,24
16468	4,65E-8   2,04 21849   7,59E-10 2,10 26482   –         –

Проведенное исследование позволяет сделать следующие выводы.

• 1.    Порядок сходимости. В целом порядок сходимости R ^ 2, который стабилизируется при увеличении N cells и уменьшении параметра сгущения η. Это согласуется с результатами МКР [14] и МКНК [8] на равномерных сетках, где порядок также был равен двум вне зависимости от порядка схемы или степени полинома p.

• 2.    Эффект сгущения сетки. Сгущение сетки в окрестностях особенностей значительно повышает точность решения. Так, при p = 8: ^E U ||_ю = 4,15 • 10 ^-9 для п = 64 и N cells = 678 и | E UIL =4, 65 • 10 — 8 для П = 2 и N cells = 16 468. Таким образом, при увеличении сгущения сетки точность возрастает на порядок, несмотря на уменьшение числа ячеек на два порядка.

• 3.    Сравнение полиномов различной степени. Поскольку глобальный порядок сходимости R ^ 2 определяется гладкостью решения, выигрыш от увеличения степени полинома p при фиксированном числе степеней свободы DOF менее выражен, чем в примере 1. Так, для p = 8 значение погрешности ||EU|| ^ = 2, 04 • 10 ^-10 при DOF = 6730 • 45 = 302 850, а для p = 4 - ^E u | _м = 2, 73 • 10 ^-10 при DOF = 26482 • 15 = 397 230. Различие в точности оказалось относительно небольшим, что соответствует теоретическим ожиданиям: при наличии особенностей преимущество высоких значений p ограничено.

Для сравнения, в [8] на равномерной квадратной сетке и с использованием полиномов Чебышева 16-й степени в hp-МКНК погрешность составляла 4, 28 • 10 ^-8 . Аналогично, в МКР шестого порядка [14] на равномерной сетке 40 x 40 достигнута точность IE U |к = 9, 33 • 10 ^-7 при DOF = 1600. В настоящей работе достигнута точность IE U |к = 9 , 43 • 10 ^-9 при DOF = 678 • 28 = 18984. Если учесть, что в МКР [14] порядок сходимости равен двум, то экстраполяция показывает: достижение точности порядка 10 ^-9 потребовало бы DOF > 10 ⁵ , что подчеркивает эффективность сгущенных сеток.

Пример 3. Изгиб защемленной кольцевой пластины в рамках ТРМ. Рассмотрим задачу (3), (4) в кольце с точным решением [5, 19]

log r   r ² log r    r ²       C 2 r ²                        qr ²       qr ⁴

w = C ¹ \GtK — 4D ⁺ 8dJ — ~ — C ³ logr ⁺ C ⁴ — 4GtK + 64D, (9)

• Ф„    = (C2 + ^С + — log Г -     ^ X, Фу = fC2 + ^C + — log Г -     ^ У, (10)

• ^{x    2}    r ²    2D        16D   ^{, y 2} r ²    2D        16D ^,

где r = x ² + y ² — полярный радиус, G = E/(2(1 + ν)) — модуль сдвига, константы C 1 -C 4 определяются из решения СЛАУ размера 4 x 4 с учетом краевых условий (4).

В вычислениях (см. табл. 3) использовались границы области x ² ₁ + x ² ₂ = 5 ² и X i + X 2 = 1 2 , а параметры материала принимались равными E = 2 • 10 5 , v = 0, 28. Весовые множители выбирались следующим образом: в уравнениях коллокации – h _j ² ; в условии согласования при w, φ _x , φ _y – 1, а перед производными по нормали – h j ; для краевых условий – 1. По умолчанию N _col , N _match и N _bound совпадали со значениями, использованными в примерах 1 и 2.

Основные выводы по результатам расчетов.

• 1.    Сходимость и точность. Порядок сходимости R в среднем не ниже p-го при четных p и (р — 1)-го - для нечетных р. Колебания R связаны с эффектом сдвигового запирания, малого параметра β ² при старших производных φ x и φ y , а также ростом числа обусловленности для высоких p и достаточно подробных сеток. hp-МКНК демонстрирует более высокую точность, чем h-KM ₄ при p >  6, что ожидаемо, поскольку решение в области Ω бесконечно гладкое. Однако с учетом поведения сходимости оптимальными и рекомендуемыми значениями являются p = 4, 5, 6.

• 2.    Эффект сдвигового запирания. Известно, что при t ^ 0 сходимость численных методов проявляется лишь на достаточно подробных сетках [17, 20, 22]. Это под-

• тверждается и в данном примере: при p = 2 при t = 1, β2 = 10-2 сходимости не наблюдается; t = 10-3, β = 10-8 сходимости нет ни для одного из рассматриваемых

• 3.    Влияние переопределения СЛАУ. Отдельно были протестированы варианты метода с квадратными матрицами при p = 3 и p = 5:

p. Дополнительным осложняющим фактором является рост числа обусловленности

СЛАУ с увеличением степени полинома p.

Таблица 3

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

N cells	IE W | ^    R    | Е ф х k R					IE W к    R   ЦЕ ф х к   R
	p = 3, t = 1, β 2 = 10 -2					p=4, t= 1, β 2 = 10 -2 [			5]
254 820 2972 11374	3,73E-2 1,36E-2 3,85E-3 1,10E-3	1,72 1,96 1,87	3,66E-2 1,14E-2 3,68E-3 1,01E-3	1,99 1,75 1,93		1,29E-3 6,41E-5 2,22E-6 1,13E-7	5,12 5,22 4,43	2,15E-3 1,67E-4 6,18E-6 3,12E-7	4,36 5,12 4,44
	p=4, t= 10 -1 , β 2 = 10 -4 [					p=4, t= 10 -2 , β 2 = 10 -6 [				5]
254 820 2972 11374	2,52E-2 1,65E-3 9,34E-5 5,47E-6	5,22 4,46 4,22	3,49E-2 1,11E-3 1,02E-4 6,69E-6	5,88 3,70 4,05		2,43E+0 1,85E+0 1,11E-2 3,14E-4	0,46 7,94 5,31	1,79E-4 1,54E+0 1,10E-2 2,71E-4	0,25 7,67 5,51
	p=5, t= 10 -1 , β 2 = 10 -4					p=5, t= 10 -2 , β 2 = 10 -6
254 820 2972 11374	1,70E-2 6,01E-4 8,77E-6 5,27E-7	5,70 6,56 4,19	1,32E-2 3,97E-4 9,15E-6 6,18E-7	5,97 5,85 4,02		1,21E+0 1,89E-1 9,33E-3 6,59E-5	3,16 4,67 7,38	1,21E+0 1,97E-1 9,00E-3 5,17E-5	3,09 4,79 7,69
	p=6, t= 10 -1 , β 2 = 10 -4					p=6, t= 10 -2 , β 2 = 10 -6
254 820 2972	1,21E-3 7,25E-6 1,92E-7	8,73 5,63	1,10E-3 7,91E-6 1,68E-7	8,42 5,98		1,16E-2 1,20E-3 6,78E-5	3,87 4,46	9,76E-3 1,35E-3 7,15E-5	3,37 4,56
	p=8, t= 10 -1 , β 2 = 10 -4					p=8, t= 10 -2 , β 2 = 10 -6
254 820 2972	3,36E-5 4,22E-8 7,23E-9	11,39 2,74	2,80E-5 6,33E-8 5,47E-9	10.39 3,80		2,08E-4 9,01E-5 4,41E-5	1,42 1,10	2,01E-4 8,51E-5 3,37E-5	1,46 1,43

• 1)    при p = 3, t =1 и N _cells = 11374 удалось получить ||EW||_ю = 5, 62 • 10 ^-2 , что хуже более чем на порядок по сравнению с применением переопределенной СЛАУ;

• 2)    при p = 5, t = 10 ^-1 (тонкая пластина) и N cells = 11374 значение погрешности | E W ||^ = 5, 65 • 10 ^-3 , уступая варианту с переопределением почти на четыре порядка (табличное значение - 5, 27 • 10 ^-7 ).

• 4.    Особенность задач типа (9) , (10) . Решения включают логарифмический член log r, который стремится к то при r ^ 0. Несмотря на то, что точка ветвления (0, 0) / Q, ее наличие препятствует сходимости некоторых бессеточных коллокацион-ных методов [15, 16], использующих глобальную полиномиальную аппроксимацию на прямоугольнике, включающем исходную области с криволинейной границей; влияет на структуру погрешности даже в сеточных методах [2]. Поэтому для подобных задач необходимо сгущать сетку к отверстию, что и было реализовано.

Таким образом, переопределение СЛАУ оказывает существенное влияние на точность, особенно для тонких пластин.

Заключение

В работе представлен hp-МКНК, основанный на локальных полиномиальных аппроксимациях на треугольных сетках. На трех тестовых примерах с выраженными особенностями решения – большими градиентами, разрывами производных и наличием точки ветвления - показано, что метод обеспечивает высокую точность и при p >  6

превосходит h-КМ ₄ . Путем анализа числа обусловленности, степени переопределения СЛАУ и сравнения результатов расчетов определены оптимальные комбинации параметров hp-МКНК. Для решения задач НДС наиболее практичными являются степени p = 4, 5, 6. В задачах с локальными особенностями продемонстрированы явные преимущества hp-МКНК по сравнению с h-КМ 4 и МКР [14], реализованного на равномерной сетке, благодаря возможности сочетать сгущение сетки и увеличение степени аппроксимации.

Метод характеризуется универсальностью компьютерной реализации: варьируя степень p и количество точек записи уравнений N col , N match , N bound , можно эффективно использовать возможности hp-подхода и выбирать подходящие варианты СЛАУ под конкретную задачу. Перспективным направлением дальнейших исследований является разработка и анализ hp-МКНК для краевых задач с производными высоких порядков, в частности для бигармонического уравнения, широко встречающегося в задачах механики сплошных сред [4, 13].

Работа выполнена в рамках государственного задания ИТПМ СО РАН (номер гос. регистрации 124021400039-8).