Адаптивні апроксимації методу скінченних елементів для задач еластостатики (стр. 1 )

УДК 519.6:539.3

АДАПТИВНІ АПРОКСИМАЦІЇ МЕТОДУ СКІНЧЕННИХ ЕЛЕМЕНТІВ

ДЛЯ ЗАДАЧ ЕЛАСТОСТАТИКИ

Г. Квасниця, Г. Шинкаренко

Львівський національний університет імені Івана Франка

вул. Університетська, 1, м. Львів, 79000, e-mail: *****@***

Запропоновано адаптивну схему методу скінченних елементів для знаходження наближених розв’язків варіаційної задачі еластостатики в термінах вектора зміщень. За фіксованого вибору класу апроксимацій для зміщень ця схема дає змогу керувати локальним згущенням/розрідженням тріангуляції для досягнення наперед заданої точності розрахунків.

Теоретичний базис цієї схеми ґрунтується на аналізі функціонала джерел похибки та його зв’язку з варіаційною задачею про похибку дискретизації Гальоркіна. Наближений розв’язок цієї задачі з використанням просторів бабл-функцій дає змогу сконструювати недорогі зручні апостеріорні оцінювачі енергетичної норми похибки апроксимацій МСЕ.

Можливості запропонованої схеми доведені результатами числових розв’язків двовимірної задачі з кусково-квадратичною апроксимацією зміщень на трикутних сітках.

Ключові слова: метод скінченних елементів, функціонал похибок, апостеріорні оцінки похибки, h-адаптування тріангуляцій.

1. ВСТУП

За останні 15 років в теорії і практиці методу скінченних елементів (МСЕ) досягнуто значного прогресу щодо побудови та дослідження так званих h-адаптивних схем. Загальна концепція таких схем має на меті розв’язування задачі оптимального керування тріангуляцією, здатної створити можливості відтворення з мінімальними обчислювальними витратами структури шуканого розв’язку варіаційної задачі з наперед заданою точністю. Хоча оригінальні цілісні результати цього підходу є в праці [6], перші застосування до задач еластостатики зроблено майже через десять років [12]; стосовно огляду досягнень та літературних джерел з цього приводу див. [11].

Головні труднощі ефективної реалізації концепції h-адаптивних МСЕ пов’язані із конструюванням

§ алгоритмів гнучкого тріангулювання областей визначення шуканих розв’язків;

§ надійних і зручних засобів оцінки апостеріорних похибок дискретних розв’язків.

Зростання попиту на засоби тріангулювання зумовило появу потужних методів побудови неструктурованих сіток на основі тріангуляцій Делоне, (див. [8]). Стосовно другої, теоретичної частини концепції h-адаптивності, то, на нашу думку, перспективними є дослідження, пов’язані з аналізом лишків (нев’язок) вихідних рівнянь, локалізованих над кожним скінченним елементом тріангуляції. Нещодавні результати [7, 9, 10] пов’язані із побудовою належних ваг для локальних лишків на основі спряжених функцій та двоїстості заклали міцний фундамент. Поряд із цим числові результати, див. напр. [2] свідчать, що h-адаптивні МСЕ здатні успішно конкурувати зі стабілізувальними схемами МСЕ [3], забезпечуючи монотонне прямування похибок дискретизації до нуля з оптимальними порядками швидкості збіжності.

Ми мали на меті поширити підхід до побудови h-адаптивних МСЕ [4] на крайові задачі теорії пружності. За основу аналізу взято варіаційну задачу для визначення похибки апроксимації МСЕ, лінійний функціонал рівняння якої визначений кусково-визначеним лишком рівнянь рівноваги. Для відшукання наближених розв’язків цієї задачі ми використали простори бабл-функцій, базиси яких визначені структурою вжитої тріангуляції та простору апроксимацій МСЕ. Оцінки енергетичної норми знайдених розв’язків природним шляхом приводять нас до недорогих і зручних локальних апостеріорних оцінювачів похибок апроксимацій МСЕ для вихідної задачі еластостатики. Вони ідентифікують залежність рівня похибки МСЕ як від геометричних характеристик МСЕ, так і розподілу незбалансованих сил на кожному скінченному елементі. Як застосування запропонованої методики, ми наводимо результати числового аналізу модельних задач.

2. ФОРМУЛЮВАННЯ ЗАДАЧ ЕЛАСТОСТАТИКИ

Нехай пружне тіло займає обмежену область точок евклідового простору з неперервною за Ліпшицем межею . Позначимо через одиничний вектор зовнішньої нормалі до ,

Процес стаціонарного деформування пружного тіла під дією масових сил та поверхневих напружень описує така крайова задача еластостатики:

(2.1)

Тут коефіцієнти пружності задовольняють звичайні умови симетрії та еліптичності

(2.2)

скрізь тут і далі за індексами, що повторюються, передбачена звичайна домовленість про підсумовування в межах від 1 до ;

Увівши простір кінематично допустимих векторів зміщень

(2.3)

сформулюємо відповідну (2.1) варіаційну задачу еластостатики в термінах зміщень:

(2.4)

де

(2.5)

та

За припущень, що

(2.6)

та умов (2.2) варіаційна задача (2.4) має єдиний розв’язок , і

(2.7)

де

(2.8)

є природною (енергетичною) нормою в просторі і  – відповідною нормою в спряженому просторі .

3. ДИСКРЕТИЗАЦІЯ ВАРІАЦІЙНОЇ ЗАДАЧІ

Перенесення розв’язування (2.4) в скінченновимірні підпростори із приводить до дискретизованих задач для відшукання апроксимацій Гальоркіна вигляду

(3.1)

Добре відомо, що в задачах еластостатики апроксимації Гальоркіна є ортогональними (відносно енергетичного скалярного добутку ) проекціями розв’язку на простори апроксимацій і в цьому разі їхні похибки

(3.2)

мають властивість найкращих наближень

(3.3)

Якщо простори апроксимацій побудовані за технологією МСЕ, то в більшості випадків апроксимації забезпечують бажану якість наближення до точного розв’язку задачі еластостатики (2.4). Однак якщо навантаження сильно нерегулярне, матеріал пружного тіла яскраво гетерогенний чи простежуються інші сингулярності конструкції, то такі апроксимації можуть досягатись надмірно великою ціною, і важко сказати, чи можна їх гарантовано знаходити із наперед заданою точністю. Ця обставина спонукає до шукання надійних критеріїв апостеріорної оцінки якості знайдених апроксимацій.

4. АНАЛІЗ ФУНКЦІОНАЛА ПОХИБОК

Поряд з енергетичною нормою із (2.7) наділимо простір допустимих зміщень стандартною нормою

спряжений простір  – нормою

Пропозиція 4.1 про функціонал похибок.

Нехай  – розв’язок задачі (2.4).

Для будь-якого фіксованого визначимо лінійний функціонал (джерел похибки) згідно з правилом

(4.1)

Тоді будуть правильними такі твердження

(!) (4.2)

де така, що

(4.3)

(!!) (4.4)

Доведення. Із визначення (4.1) та задачі (2.4) випливає, що

(4.5)

Звідси безпосередньо обчислюємо

і

(4.6)

Далі, приймаючи , із рівності (4.5) одержуємо, що

(4.7)

що разом із (4.6) приводить до двосторонньої оцінки норми функціонала в в термінах стандартної норми простору допустимих зміщень.

Нарешті, замінивши в наших оцінках (4.6) та (4.7) норму на енергетичну норму , отримаємо бажану рівність (4.4).

Для потреб адаптування в МСЕ важливим є такий наслідок.

Наслідок 4.1 про декомпозицію функціонала похибок.

Нехай – будь-яка регулярна тріангуляція області і , , за допомогою якої побудовано простір апроксимацій схеми МСЕ.

Тоді функціонал похибок (4.1) допускає таку декомпозицію

(4.8)

Розпочинаючи із наведення (4.5) за допомогою інтегрування частинами на кожному скінченному елементі, знайдемо, що

(4.9)

Тепер врахуємо рівняння рівноваги крайової задачі (2.1) та неперервність функції у разі переходу через спільну грань суміжних скінченних елементів. У результаті відповідних спрощень отримаємо задекларовану декомпозицію (4.8) функціонала похибок.

На цьому етапі ми здатні схарактеризувати зв’язок функціоналу (4.1) із похибкою дискретизації Гальоркіна.

Наслідок 4.2 про структуру функціонала похибок дискретизації Гальоркіна.

Нехай виконано умови пропозиції 4.1 та наслідку 4.1 і  – апроксимація Гальоркіна, знайдена як розв’язок задачі (3.1).

Тоді будуть правильними такі твердження

(4.10)

(4.11)

(4.12)

Доводять цей наслідок безпосередньо за допомогою обчислень на основі (4.5) при .

5. ВАРІАЦІЙНА ЗАДАЧА ДЛЯ ПОХИБКИ ДИСКРЕТИЗАЦІЇ ГАЛЬОРКІНА

Наведені вище результати достатньо повно характеризують похибку дискретизації Гальоркіна з позицій функціонала , визначеного в (4.1). Зокрема, з огляду на (4.11) ми тепер здатні сформулювати варіаційну задачу про відшукання похибки дискретизації Гальоркіна:

(5.1)

Знайти розв’язок (5.1) так само важко, як і вихідної задачі еластостатики (2.4). Тому тут ми вдамося до її дискретизації за схемою Гальоркіна в деякому скінченновимірному підпросторі із простору похибок . Щоб у майбутньому цю процедуру можна було виконати на практиці легко, ми пов’яжемо її реалізацію із тріангуляцією , за допомогою якої знайдено розв’язок дискретної задачі еластостатики (3.1), і сформулюємо таку задачу:

(5.2)

Тепер, повертаючись до наслідку (4.2), ми здатні визначити важливий результат.

Теорема 5.1 про коректність задачі про похибку дискретизації.

Нехай варіаційна задача еластостатики (2.4) має єдиний розв’язок і  – його апроксимація Гальоркіна, знайдена із (3.1) з довільним .

Тоді варіаційна задача (5.2) для відшукання апроксимації похибки дискретизації є коректно поставленою і

, (5.3)

. (5.4)

Отже, енергетична норма апроксимації похибки дискретизації Гальоркіна не перевищує норми функціонала джерел похибки, причому, як свідчить (5.4), різницю між цими нормами можна оцінити, якщо ми в стані обчислити норму . З цією метою можна скористатись інтерполяційними оцінками, типовими для апріорних похибок МСЕ (див. напр. [5]).

Нижче ми пропонуємо процедуру обчислення енергетичних норм

(5.5)

з використанням бабл-функцій для побудови просторів . Притаманна цим функціям ортогональність значно здешевлює знаходження (5.5) на кожному скінченному елементі насправді вона потребує на кожному з них розв’язування системи лінійних алгебричних рівнянь порядку .

6. БАБЛ-АПРОКСИМАЦІЯ ПОХИБКИ ДИСКРЕТИЗАЦІЇ МСЕ

На заданій тріангуляції розглянемо кусково-визначену функцію такого вигляду:

(6.1)

де  – точка, наприклад, центр ваги, скінченного елемента

Виберемо тепер в задачі (5.2) систему функцій за базис простору . Це дає нам змогу записати її розв’язок у вигляді лінійної комбінації

(6.2)

з невідомими векторами . Фізичний зміст векторів стає зрозумілим після підстановки в (6.2) і врахування (6.1):

(6.3)

Для визначення цих векторів ми підставляємо (6.2) у рівняння задачі (5.2) і послідовно приймаємо

(6.4)

У результаті цього отримаємо послідовність систем лінійних алгебричних рівнянь вигляду

(6.5)

коефіцієнти якої обчислюють згідно з правилами

(6.6)

Отже, використання бабл-функцій дає змогу задачу відшукання системи векторів декомпонувати в незалежне послідовне розв’язування систем лінійних алгебричних рівнянь -го порядку на кожному скінченному елементі тріангуляції . Цілком очевидно, що матриці цих систем рівнянь додатно визначені, що дозволяє знайти єдиним чином похибку (6.2).

Теорема 6.1 про оцінювачі апостеріорної похибки дискретизації.

Нехай  – апроксимація Гальоркіна для розв’язку задачі еластостатики (2.4), знайдена на тріангуляції . Нехай є послідовністю функцій з властивостями (6.1).

Тоді апроксимація похибки однозначно визначена лінійною комбінацією (6.2) шляхом розв’язування систем лінійних алгебричних рівнянь (6.5) і в цьому разі величини

(6.7)

та

(6.8)

є локальними апостеріорними оцінювачами похибки наближення розв’язку на кожному скінченному елементі та тріангуляції в цілому, відповідно.

7. Результати обчислювального експерименту.

Пропонований підхід застосовано до розв’язування модельної задачі з [10] про пружну деформацію конструкції з тріщиною. Квадратна пластинка з тріщиною рівномірно навантажена на частині верхньої межі та жорстко закріплена вздовж нижньої та правої сторін (див. рис. 1). Ця задача цікава тим, що її розв’язок має сингулярність на вершині тріщини, де простежується концентрація напружень з асимптотичною поведінкою вигляду , ( – відстань від кінця тріщини). В цьому прикладі пружні сталі матеріалу пластинки Н/м2 (алюміній), а інтенсивність поверхневого тиску Н/м2.

Головна мета нашого обчислювального експерименту – дослідження характеристик запропонованої вище процедури -адаптування з використанням оцінювачів похибки (6.7) та (6.8). Ітераційний процес адаптування розпочинався з дуже грубої тріангуляції , зображеної на рис. 2. Цю та всі наступні тріангуляції , , перебудовували за допомогою пакета triangle, розробленого Дж. Шевчуком (див. http://www. cs. cmu. edu/~quake/triangle. html).

Якщо хоча б на одному зі скінченних елементів біжучої тріангуляції виконувалась умова

(7.1)

то цей та решта подібних йому скінченних елементів із підлягали подрібненню в такий спосіб, щоб елементи нової тріангуляції не містили кутів, менших за 25о. Деякі зі знайдених тріангуляцій показано на рис. 3 та 4.

Оскільки умова (7.1) повинна врешті-решт приводити до тріангуляції з рівномірно розподіленою між елементами похибкою, то рис. 3, 4 надають якісне уявлення про динаміку реагування цієї адаптивної схеми МСЕ на різні види сингулярностей розв’язку розглядуваної задачі еластостатики.

Рис. 2. Початкова тріангуляція Т0

Рис.1. Геометрія тестового прикладу

Рис. 4. Структура тріангуляції Т15.

Рис. 3. Структура тріангуляції Т5.

На початку (див. рис. 3) простежується інтенсивне згущення сітки в околі правого верхнього кута, у досить великому околі якого накопичена значна частина розподіленої в пластині енергії. На цій стадії сингулярність на кінці тріщини виявляє себе менше, оскільки ця особливість локалізована в дуже малому околі вершини тріщини. На рис. 4 виявляються знайдені адаптуванням інші особливості точного розв’язку: в околі середини верхньої сторони (зумовлена стрибком навантаження у цій точці) та лівої нижньої вершини квадрата (зумовлена зміною виду однорідних крайових умов задачі).

Деякі кількісні результати нашого експерименту наведено в табл. 7.1 та 7.2. У цілому вони відображають, що після виконання десяти кроків адаптування обчислені характеристики виходять на асимптотичні показники зі стабільними значеннями порядків швидкості збіжності.

Тут  – номер ітерації;  – загальна кількість трикутників у тріангуляції ;  – загальна кількість вузлів тріангуляції;  – діаметр біжучої тріангуляції;  – порядок швидкості збіжності схеми МСЕ в нормі () та в енергетичній нормі () відповідно.

Наприклад, показники швидкості збіжності МСЕ в -нормі та енергетичній нормі близькі, відповідно, до другого та першого порядків відносно параметра дискретизації

, (7.2)

де  – загальна кількість вузлів біжучої тріангуляції . У цьому разі порядки швидкості збіжності обчислювали згідно з правилом

, (7.3)

де  – одна із норм або .

Таблиця 7.1.

Збіжність обчислених характеристик -адаптування з кусково-квадратичними апроксимаціями МСЕ

0	16	45	0,55902	252,8130		1,4417
1	26	67	0,50000	176,6930	1,80	0,7333	3,40
2	51	124	0,50000	141,7940	1,14	0,5275	1,98
3	105	242	0,50000	115,4270	0,93	0,3686	1,62
4	219	486	0,50000	83,4260	0,93	0,1706	1,79
5	353	766	0,35355	56,3991	1,06	0,1016	1,87
6	612	1297	0,29702	38,8304	1,11	0,0739	1,77
7	861	1806	0,29702	32,0846	1,12	0,0535	1,78
8	1066	2233	0,27951	28,8475	1,11	0,0329	1,94
9	1252	2611	0,27951	28,0356	1,08	0,0315	1,88
10	1315	2740	0,27951	27,9104	1,07	0,0315	1,86
11	1352	2819	0,27951	27,6135	1,07	0,0300	1,87
12	1376	2869	0,27951	26,9855	1,08	0,0266	1,92
13	1395	2908	0,27951	26,9024	1,07	0,0263	1,92
14	1410	2939	0,27951	26,8129	1,07	0,0262	1,92

Крім того, в табл. 7.2 наведені деякі додаткові величини, контрольовані в ході обчислювального експерименту, а саме – потенціальна енергія деформування;  – значення -інтеграла;  – максимальне із значень модуля зміщень у вузлах тріангуляції;  – максимальне із значень похибки зміщень у центрах ваг трикутників.