Ю Ящук - Адаптивний алгоритм чисельного дослідження задачі теорії пружності - страница 1

Страницы:
1  2 

ВІСНИК ЛЬВІВ. УН-ТУ Серія прикл. матем. інформ. 2010. Вип. 16. C. 96-105

VISNYK LVIV UNIV. Ser. Appl. Math. Inform. 2010. Is. 16. P. 96-105

УДК 517.958:519.6

АДАПТИВНИЙ АЛГОРИТМ ЧИСЕЛЬНОГО ДОСЛІДЖЕННЯ ЗАДАЧІ ТЕОРІЇ ПРУЖНОСТІ

Ю. Ящук

Львівський національний університет імені Івана Франка, вул. Університетська, 1, Львів, 79000, e-mail: yuriy.yashchuk@gmail.com

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

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

1. ВСТУП

Метод скінченних елементів (МСЕ) широко розповсюджений потужний чисельний засіб розв'язання задач математичної фізики [5, 12]. Адаптивні схеми цього методу допомагають оптимізувати обчислювальні ресурси для досягнення надійності та ефективності [9, 12]. У цьому разі виникає проблема розв'язування задачі на несумісному розбитті, ефективним вирішенням якої є використання мортарних функцій [6, 7, 10]. Важливим етапом адаптивного процесу є оцінка похибки результату. Найпоширеніші способи проведення такої оцінки або використовують відновлені значення невідомих функцій, або в явному чи неявному вигляді використовують нев'язки скінченноелементної апроксимації. Різні способи оцінок наведено у працях [3, 9, 11]. Оцінки апроксимацій з використанням мортарних функцій опрацьовано у [8] разом із процесом адаптації, який ґрунтується на такій оцінці. Наша мета - дослідити можливості побудови одного нового способу оцінки похибки та його застосовності до розробки адаптивних схем.

Поряд із МСЕ для знаходження чисельних розв'язків задач теорії пружності активно використовують метод граничних елементів, зокрема його прямий варіант (ПМГЕ) [0]. У цьому випадку дискретизується не весь об'єм тіла, а лише його границя, що дає підстави значно зменшити порядок результуючої системи рівнянь. Важливим є також той факт, що в задачах пружності ПМГЕ дає змогу обчислити переміщення та напруження з однаковим порядком точності, тоді як у МСЕ обчислені напруження мають точність на порядок нижчу, ніж переміщення. Тому важливо дослідити залежність між локальними похибками результатів обох методів (переміщень і напружень). Це допоможе використовувати порівняння результатів МСЕ та прямого методу граничних елементів (ПМГЕ) для оцінювання похибки кожного з методів, та побудови алгоритму адаптації сітки.

2. ФОРМУЛЮВАННЯ ЗАДАЧІ

Розглянемо двовимірну область ££ , обмежену границею Г. Позначимо через n зовнішню нормаль до границі. У згаданій області статичну двовимірну задачу

© Ящук Ю., 2010

АДАПТИВНИЙ АЛГОРИТМ ЧИСЕЛЬНОГО ДОСЛІДЖЕННЯ^

пружності записуємо в переміщеннях u = (u1,u2) так:

juhu + +ju)grad divu = 0, x є ££ (1)

0 u (2)

p = Po> xє Гр■

Тут Л , - сталі Ляме; ru - частина границі, на якій відомо переміщення u0; Г  - частина границі, на якій відомо зусилля p0, причому ru U Гp = Г , Ги П Гp = 0

Зусилля p = (p1,p2)T виражають через напруження a = (a11,a22,a12)т та нормаль

n = (n1, n2 )T

pt = І]аЛ; i = 1,2. (3)

j=1

Компоненти  тензора  напружень   виражають   через  компоненти вектора переміщень

^ Эх j dxt

И i   -л -л       i r

ЭХ1 Эх­

i, j = 1,2, (4)

де Stj - символ Кронекера.

Співвідношення Коші набувають вигляду 1

( 3uL +<Эм1 ^

i, j = 1,2 (5)

Застосуємо до чисельного розв'язування задачі (1-5) метод скінченних елементів з використанням мортарних функцій і прямий метод граничних елементів.

3. ВИКОРИСТАННЯ МОРТАРНИХ ФУНКЦІЙ У МСЕ

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

Нехай область ££ розбита на підобласті Q.i, i = 1,K так, що Ц П ££j =0 ,

K _

i Ф j та U Ц = ££ . Спільну границю двох підобластей Гі} = ЭЦ П d££j називатимемо

інтерфейсом. Кожна підобласть Q.l (незалежно від інших підобластей) буде розбита на скінченні елементи Q.ie з визначеною апроксимацією шуканих функцій. Позначимо через Vl простір, натягнутий на базисні функції скінченноелементного розбиття Q.i . Зауважимо, що ніяких додаткових умов на сумісність скінченноелементного розбиття на межі підобластей не накладається. За такої апроксимації на інтерфейсі Гу може порушуватися умова неперервності переміщень

u , тому вимагається виконання умови слабкої неперервності [7]

J -Uj)(pkdTv = 0; к = 1,...,M , (6)

де рк = рк (x), x єГ у , k = 1,...,M - набір базисних мортарних функцій на деякому розбитті інтерфейсу Гу. Зокрема, у [10] пропонують брати за вузли розбиття інтерфейсу Г у вузли скінченноелементного розбиття Ц або Ц j, які припадають на Г у , та побудувати на них класичний мортарний базис, схематично зображений на рис. 1.

Рис. 1. Класичний мортарний базис

Позначимо через Л простір, натягнутий на базис рк, а також задамо простір

K

V = Y\V . Після згаданих апроксимацій можна записати таку варіаційну задачу, що

відповідає задачі (1-2).

Знайти такі (uh,Лк) є (X,Л) , що

a (uh,vh) + b,vh ) = f (vh),  vh є X,

b(Vh,uh ) = 0,  /и„ єЛ, де a(;■) та f (•) - відповідна білінійна та лінійна форми задачі (1)-(2); білінійна форма b (■, ■) враховує умову слабкої неперервності (6)

b (jU, u )= J (u - иу )//й?Гу . (8)

Застосувавши метод Бубнова-Гальоркіна [5], знаходимо наближений розв'язок задачі (7).

4. ДОСЛІДЖЕННЯ ПОХИБОК

Далі за норму приймаємо L2 - норму деякого вектора a по області S

a   = /І aTadS , (9)

або у випадку скалярної величини a : ||a||S =    a2dS .

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

Of =^П +CT222 +CT323 + 2СТ122. (10)

АДАПТИВНИЙ АЛГОРИТМ ЧИСЕЛЬНОГО ДОСЛІДЖЕННЯ^

Величина о33 у випадку плоского деформованого стану залежить від о11 та

Л

Г(°11 + o22

(11)

33   2 + ju)

Позначимо через of (x) та ob (x) тензори напружень у деякій точці x,

обчислені відповідно за методом скінченних і граничних елементів. Тоді для кожної точки x визначимо ефективне значення різниці обчислених напружень так:

-о,

F 22

+ 2 (о,

+

Норма такої різниці напружень на елементі набуде вигляду ІІАо„

■і

Ao„Bd Qe

(12)

(13)

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

\AOFB

AoWB[1]d Q

J 1d Qe

(14)

Тобто під Qe  треба розуміти площу елемента Qe.

Аналогічно можна обчислити середньоквадратичне ефективне напруження обчислене за МСЕ та ПМГЕ, відповідно

\°f\

\\°b\

JOBeff 2d Q

J 1d Qe

(15)

(16)

Оскільки oB (x) обчислюється з ліпшим порядком точності, ніж oF (x), то

певне наближення відносної похибки обчислення напружень на елементі можна визначати за такою формулою:

Г)=­

AOFB

(17)

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

= A~°FBq, _IIa°fb|„71 N1

nFB - II      II    /11^ II      . (18)

\°А J\H\

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

(0fu -+ (о„ - от.

F22 + 2 (о.

T22

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

Ао,„

Ы L/INI

(19)

(20)

Аналогічні до  (19) та (20)  формули можна записати для напружень, обчислених за ПМГЕ.

+ (о,

- OT

+ в

B33       T33 і AObtII

+2 к

Ы \J\H\

(21) (22)

(23) (24) (25)

(26) (27) (28)

Для визначення того, наскільки адекватно AoFB чи AuFB відображає похибку, відповідно oF чи uF, обчислимо так звані індекси ефективності [11]

За цим самим принципом обчислимо віднесені похибки переміщень

AuFB (x) =V(uF1 -uB1 )2 +(uF2 -uB2 )2 Au,T (x) ^(u,1 - un )2 +(u, 2 - ut 2 )2 AuBT (x) = V(uB1 - uT1 )   +(uB2 - uT2 )

Au

\\ub\\at

/\\a\\

II Au FTII II    FT1 la

/!» II

W/i II

/\\a\\

II t||a/

II II

IIAubt| L,

/INI

Ы la/

IIM

АДАПТИВНИЙ АЛГОРИТМ ЧИСЕЛЬНОГО ДОСЛІДЖЕННЯ^

1Аа*Л а

в = \-р- (29)

\АаА L,

в = h~^(30)

Чим ближчим до одиниці є такий індекс, тим ліпше величина, що є в чисельнику, наближає реальну похибку, яка є в знаменнику. Отож, близькість індексу до одиниці означає, що віднесена різниця r/rB чи fjrB адекватно відображає реальну віднесену похибку напружень або переміщень, відповідно.

5. ЧИСЕЛЬНІ РЕЗУЛЬТАТИ

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

, Ро

1 wvW^y

Рис. 2. Схема тестової задачі

Позаяк метою проведених досліджень є побудова адаптивних сіток, то згущення розбиття проводили відповідно до відомого розташування зон зміни граничних умов і геометрії задачі. У нашому випадку зоною з високими градієнтами є ліва частина об'єкта. Вузлами дискретизації границі у ПМГЕ вибирали вузли скінченноелементної (СЕ) дискретизації, які потрапили на границю. Для кожного розбиття обчислювали локальні індекси ефективності (29) чи (30). Порівнювали переміщення та напруження, отримані МСЕ та ПМГЕ на рівномірних і нерівномірних розбиттях.

Рівномірне розбиття об'єкта задачі проводили, відповідно, на 8, 32 та 128 прямокутних скінченних елементи. Схему одного з цих розбиттів показано на рис. 3. Результати підтверджують, що на рідкій сітці похибка обчислених переміщень ПМГЕ набагато менша, ніж похибка обчислених переміщень МСЕ. Після згущення СЕ сітки результати МСЕ та ПМГЕ стають ближчими, а на сітці зі 128 елементів отримані переміщення практично збігаються. Як наслідок, різниця в переміщеннях між методами погано характеризує реальну похибку переміщень (рис. 3).

Варто зауважити, що деякий зв'язок між віднесеною різницею переміщень і реальною віднесеною похибкою переміщень все ж простежується. fjrB та fjrT пропорційні між собою; коефіцієнт пропорційності змінюється залежно від розбиття. Оскільки результати ПМГЕ та МСЕ теж пропорційні між собою, то це ж можна сказати і про fjrB та fjBT . Подібний зв'язок можна простежити і на нерівномірних сітках.

а

0 2

Рис. 3. Віднесені різниці fjFB (а) та індекси ефективності в (б)

Розглянемо результати аналогічних досліджень, проведених для похибок напружень. Їх обчислюємо за ПМГЕ з таким самим порядком точності, як і переміщення - O (h2), тоді як за МСЕ цей порядок є O (h) (завдяки чисельному

диференціюванню отриманих переміщень). Оскільки напруження за МСЕ та ПМГЕ отримуються різними за природою способами, між ними не спостерігають лінійної залежності, як це було у переміщеннях. Зважаючи на значну різницю в похибках, стає очевидним, що віднесена різниця nFB мало відрізняється від похибки nFT , тож індекс ефективності близький до одиниці (рис. 4). Навіть на розбитті з восьми елементів в лежить в діапазоні від 0.95 до 1.05, що свідчить про можливість використання різниці між методами для оцінки похибки напружень МСЕ. На розбитті із 192 елементів можна чітко простежити розподіл f]FB та в: віднесені похибки зростають з наближенням до лівого верхнього та лівого нижнього кутів, де відбувається стрибкова зміна граничних умов. За цим самим принципом розподілено індекси ефективності.

Наведемо деякі результати отримані на нерівномірних сітках. Нагадаємо, що згущення проводили в околі особливих точок, а саме кутів об'єкта, в яких відбувається стрибкова зміна граничних умов. Схематично одне з цих розбиттів зображено на рис. 5. Область об'єкта була поділена на 3 підобласті, у кожній з яких використовували рівномірну сітку. Обчислення на таких несумісних розбиттях проводили з використанням мортарних функцій. Це приводить до зміни розподілу похибок по об'єкту: біля інтерфейсу точність обчислених напружень гірша, оскільки тут апроксимовані переміщення мають розрив.

Страницы:
1  2 


Похожие статьи

Ю Ящук - Адаптивний алгоритм чисельного дослідження задачі теорії пружності