CC BY
25
УДК 519.6
Мотайло А. П.1, Хомченко А. Н.2, Тулученко Г. Я.3
Старший викладач кафедри вищоТ математики i математичного моделювання Херсонського национального mexHi4Hoao
унверситету, Херсон, УкраТна
2Д-р ф.-м. наук, професор, завдувач кафедри прикладноТ та вищоТ математики Чорноморського державного
ушверситету iменi Петра Могили, МиколаТв, УкраТна 3Д-р техн. наук, професор, професор кафедри вищоТ математики i математичного моделювання Херсонського
нацонального технчного унверситету, Херсон, УкраТна
_ПОБУДОВА БАЗИСУ Б1П1РАМ1ДИ_
У статт бтрамща вперше розглядаеться як 6-вузловий скшченний елемент (СЕ). Для побудови й бжвадратичного базису використовуються два рiзних шдходи: матричний cnoci6 та метод внутршньо! конденсаци базису бтрамщи як 7-вузлового СЕ. Перший пiдхiд дозволяе дослщити принципово можливу кiлькiсть базиав, а другий тако! можливостi не надае, але е бшьш економiчним. Показано, що пiсля задоволення традицшних вимог до базисних функцiй у МСЕ у бжвадратичних базисних функщях бiпiрамiди як 6-вузлового СЕ, як будуються за допомогою названих рашше пiдходiв, залишаеться рiзна кiлькiсть невизначених коефщенпв. Цi коефiцiенти надалi використовуються для надання базисним функщям спещальних властивостей, якi адаптують !х до розв'язання граничних задач iз рiвнянням Лапласа. У якостi критерiю прогностичного оцшювання апроксимацiйних властивостей СЕ у формi бiпiрамiди обрана величина слщу матрицi жорсткостi. Мiнiмiзацiя слiду матриц жорсткостi приводить до побудови одного i того ж бiквадратичного базису при обох тдходах.
На основi отриманого базису аналiзуються межi припустимих деформацш геометрично! форми бiпiрамiди. Вперше теоретично доведено, що юнуе СЕ, при використанш якого як тшрки скiнченно-елементноi сiтки, найкраща точнiсть досягаеться при вщхиленш геометрично! форми СЕ вщ правильного багатогранника, у даному випадку вщ октаедра. Знайдено критичне значення коефщента стиснення, яке забезпечуе мтмум слiду матрицi жорсткостi для бтрамщи з геометричною формою, що дослщжуеться.
Проведено обчислювальний експеримент, результати якого пщтверджують теоретичний прогноз властивостей бтрамщи як СЕ. Виявленi залежност дозволяють припустити доцiльнiсть застосування базиав бiльш високого порядку для СЕ у формi бiпiрамiди. Ключовi слова: метод скшченних елементiв, бiпiрамiда, слщ матрицi жорсткостi, тетраедрально-октаедральна решiтка.
НОМЕНКЛАТУРА
МСЕ - метод сюнченних елеменлв;
СЕ - сюнченний елемент;
K. - вузли СЕ у форм1 бМрамщи, i = 06;
a, b, h - параметри, що визначають геометричш роз-м1ри тша;
р - коефщент стиснення бМрамщи;
(x; y; z) - декартов! координати точки у тривишрно-му простора
Ni - загальна базисна функщя у форм1 повного пол-1ному другого порядку вщ трьох змшних;
к() - коефщенти при мономах у виразах базисних
функцш, i = 1;6, j = 1;10;
NS - базис бМрамщи як 7-вузлового СЕ;
NC - базис бМрамщи як 6-вузлового СЕ, отриманий в результата застосування операцп внутршньо1 конденсаци до базису NS;
NM - базис бМрамщи як 6-вузлового СЕ, побудова-ний матричним способом;
{а; ß; у} - вагов1 коефщенти для операцп внутршньо! конденсаци;
CoordKnots - матриця координат вузл1в бМрамщи;
Monom - вектор моном1в пол1ному другого порядку;
Coeff - матриця коефщ1ент1в базисних функцш бшрамщи, яю побудоваш матричним способом;
С - матриця жорсткостц
Trace - слщ матриц жорсткостц
V - компактний тополопчний простар;
¿2 (V) - середньоквадратична метрика;
C(V) - метрика простору неперервних на компакп V функцш;
s1 - похибка чисельного розв'язку у метриц L2 (V); s 2 - похибка чисельного розв'язку у метрищ C(V); T = T (x, y, z) - температура у довшьнш точщ бруса; T = T (x, y,const) - точний розв'язок задачц T = T (x, y, z) - чисельний розв'язок задачц V - область окремого скшченого елемента. ВСТУП
Головною перевагою застосування мток з включен-ням ком1рок у форм октаедр1в пор1вняно 1з виключно тет-раедральними мтками е зменшення об'ем1в розрахунюв i, як наслщок, зростання швидкост реал1заци МСЕ. Один 7-вузловий октаедр зам1нюе 8 чотиривузлов1 тетраедри.
Звичайно вимога дотримання р1вном1рност1 мтки е обтяжливою при реал1зацп методу СЕ [1]. Тому актуальною е задача дослщження впливу в1дхилень форми СЕ вщ правильного геометричного тша (в даному випадку вщ октаедра) на обчислювальш характеристики такого елемента. Особливо актуальною ця задача е для СЕ, яю використовуються у приграничному шар1, коли один (або кшька) вузл1в переносять на границю област для кращо! апроксимацп ii геометрп.
Метою дослщження е встановлення умов застосування бМрамщ у ск1нченно-елементних реш1тках. 1 ПОСТАНОВКА ЗАДАЧ1
Для досягнення мети дослщження необхщно розв'я-зати низку послщовних задач, що наводяться нижче. На основ1 заданих координат вершин бМрамщи (рис. 1):
К1,3(±а; 0; 0); K 2, 4(0; ± a; 0); K5 = (0; 0; pa);
K6 = (0; 0;-a); де a > 0, p > 0, a e R, p e R.
© Мотайло А. П., Хомченко А. Н., Тулученко Г. Я., 2016 DOI 10.15588/1607-3274-2016-4-4
Необхiдно знайти значення невiдомих коефщенпв у виразах базисних функцiй бiпiрамiди:
N = к+ к() х + y + к() z + к_4j) х 2 + kjj) y2 + + к£°z2 + k7j)xy + k8j)xz + к9°yz , j = 16, (1)
як 6-вузлового СЕ другого порядку, за двома тдходами: матричним способом (Nj = NMj) i з застосуванням опе-рацп внутршньо1 конденсацп (Nj = NCj).
Для цього обгрунтовано скласти перелж вимог до базисних функцiй бМрамщи iз вiдомих у МСЕ вимог до базисних функцш, який приводить до однозначного виз-начення коефщенлв базисних функцiй бiпiрамiди.
Показати, що функцiя слiду матрицi жорсткостi бМра-
мiди залежить вiд двох аргумента Trace = Trace(a; p) , яю визначають геометричну форму СЕ. Знайти мшмум функцп Trace(a; p) -— min для довiльного значення аргументу а та координати точки мшмуму.
Перевiрити на конкретних прикладах загальне теоре-тичне положення МСЕ про те, що кращими апроксима-цiйними властивостями володiе той базис, матриця жор-сткостi якого мае менший слад Оцiнку точностi набли-женого розв'язку, отриманого МСЕ iз використанням решггок з комiрками у формi бшрам^ та тетраедрiв, вiдносно точного розв'язку, отриманого методом Фур'е, здiйснити використовуючи середньоквадратичну метрику ¿2 (V) та метрику простору неперервних на компактi V функцш C (V).
2 ОГЛЯД Л1ТЕРАТУРИ
Вперше про застосування октаедра як СЕ стало ввдомо завдяки роботам [2, 3]. Кусково-лшшш [2] та квадратичнi базиснi функци [3] побудованi для вершин та точки перетину даагоналей багатогранника (тобто октаедр розпшдаеть-ся як 7-вузловий СЕ). При цьому автором роботи [3] вщзна-
чаеться, що октаедр i3 квадратичним базисом е неузгод-женим СЕ i3 4-вузловим (лiнiйним) тетраедром.
Як вщомо, СЕ без внутршшх вузлiв е бшьш економ-ними у планi використання обчислювальних ресурсiв при реалiзацп МСЕ. У робота [4] виконано операщю внутрш-ньо1 конденсацп та здшснено перехiд до нового базису октаедра як 6-вузлового СЕ.
Переваги застосування октаедрiв у сюнченно-елемен-тних мтках обмежуються частинами цих мток, де мають мiсце рiвномiрнi розбиття. Тому дощльною е побудова СЕ, який збер^ае переваги октаедра, але може застосо-вуватися в нерiвномiрних мтках.
Дослщжень щодо впливу деформацiй октаедра (при його трансформацп у бМрам^) на обчислювальнi характеристики такого СЕ у доступнш лiтературi не вияв-лено. Таким чином, юнуе необхiднiсть вивчення питан-ня наявност у нового математичного об'екту традицш-них для МСЕ властивостей.
Для побудови базисних функцiй СЕ нового виду зас-тосуемо матричний спомб, описаний, зокрема, у роботах [5, 7, 8], та метод внутршньо! конденсацп [6]. Перший дозволяе досл^жувати принципово можливу кшьюсть базисiв СЕ, що задовольняють певним вимогам. Другий е бiльш ращональним щодо економп обчислювальних ресурсiв порiвняно iз попереднiм способом, але обме-жуе рiзноманiггя базисiв СЕ серендипового типу лшшни-ми комбiнацiями базисних функцш вщповщного СЕ лаг-ранжевого типу.
Для прогнозування апроксимацiйних властивостей СЕ скористаемося критерiем на основi величини ^ду мат-рицi жорсткостi [8].
Критерп щодо величини припустимих вiдхилень гео-метрп СЕ вiд правильно1 форми багатокутниюв або бага-тогранникiв описаш у роботах [9-10]. 3 МАТЕР1АЛИ ТА МЕТОДИ
Побудова базису «шмрамми як 7-вузлового СЕ. Виве-демо рiвняння базисних функцiй бiпiрамiди як для СЕ з 7 вузлами. Будемо розглядати 8-гранну бМрам^, в якш вiдрiзок КоK5 (рис. 1) може бути довшьно! довжини:
Ко К 5 = p ■ a = b , де p > 0 i p e R.
Побудова перших шести базисних функцш е очевидною iз традицшно! для МСЕ вимоги рiвностi базисно1 функцп в своему вузлу одинищ i нулю - в ушх iнших вузлах [5]:
NSi, 3 = x(x ± a); NS2, 4 = y(y ± a);
NS5 =•
2a 2
1
2a2
-z(z + a); NS6 =-
1
-z(z -b). (2)
b(a + b) a(a + b)
Вираз останньо1 базисно1 функцп бМрамщи у вузлу kq знайдемо iз умови рiвностi одиницi суми вих базисних функцiй [5]:
nsq = 1 -£ nsj, де i = 1;6.
(3)
Рисунок 1 - Бтрамща як СЕ
i=1
1з (3), з урахуванням (2), маемо
1--V(p 2 + y2)+ z2 -a(p -1)z)
NSq = 1 - 2
pa
або
NSq =
(P +1)2 4 p
• F,
(4)
де F = 1 --
(x2 + y2
) (z - zo)2
1 a 2 ( p + 1)
— • a--
p4
2
2 (p + 1)2
4
1з виразу (4) видно, що поверхт рiвня базисно1 функцп NSq е елшсощами зi змiщеними по ом Oz центрами, якщо початок системи координат пов'язувати iз вузлом Kq.
Побудова узагальненого базису бiпiрамiди шляхом внутршньоТ конденсаци. Розподiлимо внесок базисно! функцп NSq (4), яка асоцшована з внутрiшнiм вузлом K о, мiж базисними функцiями, якi асоцшоваш iз зовшш-нiми вузлами, тобто здшснимо операцiю внутршньо! конденсацй [6]:
NCi = NSi + а • NS0 , i = 1..4; NC5 = NS5 +ß'NS0 , NC6 = NS6 +Y• NS0,
(5)
де 4а + ß + у = 1; 0 <а< 0,25.
Вимоги рiвностi базисно! функцп в своему вузлу оди-нищ i нулю - в умх iнших вузлах i вимога рiвностi суми базисних функцiй одиницi для функцш (5) виконуються автоматично [6]. Перевiрка виконання вимоги повноти базису (5) [7]:
6 6
х=Zxi • NCi(х;y;z); y=Zyi • NCi(x;y;z); i=1 i=1
z = Z zj • NCj (x; y; z); j=1
(6)
де (xj; yj; zj) - координати вузлiв бiпiрамiди: NCj (x; y; z) -базисш функцп бшрам^и пiсля операцп внутршньо! конденсацй, - показуе, що для виконання останшх рiвно-стей необхщно ввести додатковi спiввiдношення мiж зна-ченнями вагових коефiцiентiв:
ß=
1 - 4а p +1
Y =
p +1
•(1 - 4а).
(7)
Враховуючи iнтервал змiни вагового коефiцiента а: 0 < а < 0,25, маемо таю обмеження для коефщента в: 0 <в< 1/( p +1) та Y: 0 < y< p /(p +1).
Отже, базиснi функцп конденсованого базису визна-чаються за формулами (5) i (7).
Очевидно, що тсля задоволення усiх традицшних ви-мог до базисних функцш у МСЕ залишаеться невизначе-ним один параметр а. Ним можна скористатися для на-дання базису спещальних властивостей.
Побудова узагальненого базису бШрам1ди матрич-ним способом. Осюльки метод внутршньо1 конденсацп працюе виключно iз базисними функщями, якi поперед-ньо визначенi для елемента з внутршшм вузлом, тому доцiльно перевiрити iснування iнших базисiв бiпiрамiди з квадратичними функцiями.
Матриця координат вузлiв бiпiрамiди (рис. 1) мае виг-ляд:
f
CoordKnots =
a 0 0 1
0 a 0
- a 0 0
0 -a 0
0 0 pa
0 0 - a,
Базисш функцп будемо шукати у виглядi полiномiв другого степеню (1).
Вектор мономiв такого жшному мае вигляд:
2 2 2
Monom = (1; x; y; z; x ; y ; z ; xy; xz; yz).
1з коефщенлв всiх базисних функцiй утворимо мат-рицю коефiцiентiв:
f k (1) ko k1(1) k2) . . k9)
Coeff = k (2) k0 k(2) k (2) k2 k (2)
k (6) ^ k0 k16) k (6) k2 k (6) 9
(8)
За введених позначень вектор базисних функцiй зна-ходиться за формулою:
Т
NM = Coeff • Monom .
Врахуемо наявнi симетрп у бiпiрамiдi. Бiпiрамiда е симетричною вiдносно площини y = 0, тому базисш функцп NM\ i NMз повиннi бути парними за змшною y, тобто не мютити доданкiв iз змiнною y у непарних степенях. Це означае, що дорiвнюють нулю такi коефщенти:
k« = = kW = 0 i k23) = k<3) = k93) = 0.
Аналопчно через симетричнiсть бiпiрамiди в^нос-но площини х = 0 у виразах базисних функцш NM2 i NM4 дорiвнюють нулю таю коефщенти:
k{2)=k72) = kf=0 i k(4) = k74) = k84)=0.
Симетричнiсть бiпiрамiди вщносно тих же площин у рiвняннях базисних функцiй NM5 i NM6 приводить до рiвностi нулю коефшденпв:
k[5)=k25) = k75) = k85) = kf=0 i k[6)=k26)=k76)=kf=kf=0.
Таким чином матриця коефшденпв базисних функцiй (8) набувае вигляду:
Coeff =
Г k (1) k0 k1(1) 0 k3a) k (1) k4 k« k (1) k6 0 k8) 0
k (2) 0 0 k2(2) k3(2) k(2) k4 k5(2) k(2) k6 0 0 k9( 2)
к (3) k0 k1(3) 0 k3(3) k (3) k4 k53) k (3) k6 0 k8(3) 0
k (4) 0 0 k2(4) k3(4) k(4) k4 k5(4) k(4) k6 0 0 k(4) K9
к (5) k0 0 0 k3(5) k(5) k4 k55) k(5) k6 0 0 0
k(6) ^ k0 0 0 k3(6) k (6) k4 k5(6) k(6) k6 0 0 0
Вимога piBHOCTi одиницi базисно! функци у своему вузлу i нулю в умх iнших вузлах приводить до 6 систем piBrarn, якi в матричному виглядi записуються так:
Coeff (i, 1..10) • Monom _ in _ KnotsТ = E (i, 1..6), (9)
де i = 1;6; Е - одинична матриця pозмipом 6 x 6 ■
Матриця значень мономiв у вузлових точках мае виг-ляд:
a 0 0 2 a 0 0 0 0 0
0 a 0 0 a2 0 0 0 0
- a 0 0 2 a 0 0 0 0 0
0 -a 0 0 a2 0 0 0 0
0 0 pa 0 0 22 pa 0 0 0
0 0 -a 0 0 a2 0 0 0
Monom in Knots =
Вм шють систем (9) е сумiсними, але невизначеними, тому з них можемо однозначно знайти тшьки частину значень коефщенлв:
Coeff =
- a 2 pk6(1) 1 2a 0 a(1 - p)k® pk6X) + -Or 2a pk6(1)
- a2pk62) 0 1 2a a(1 - p)k62) pk6(2) pk^ + ^
- a 2 pk6(3) -1 2a 0 a(1 - p)k63) pk63)+ 2a pk6(3)
- a2pk64) 0 -1 2a a(1 - p)k64) pk6(4) pk^ + ^
a2 pk^5) +-- 0 0 a(1 - p)k65) + 7——-г pk65) - , (5) 1 pk6)
p +1
- a pk6 +
p +1
a(1 - P)k66) -
a(p +1)
Pkf -
a2 (p +1)
Pk66) -
a 2 (p +1) p
a2 (p +1)
kf k(3 k64 k((5 k(6
0 k(1) 0
0 0 k92)
(3)
0k
0 0 k9(4)
0 0 0 0 0 0
Врахуемо вимогу повноти базису (6). Напишемо тепер ll у матричному вигляди
Г 0
CoordKnots • Coeff =
10 0 0 00100 0 0 0 1 0
0 ^
0
0
(10)
Добуток матриць CoordKnots1 • Coeff доpiвнюе:
CoordKnots • Coeff =
a3p(k63) - k«) 1 0 a2 (1- p4« - kf) apP$ - k?) ap$ - kgj) ^ - &) 0 a(k« - kf)) a3p(k64) - kf) 0 1 a2 (1 - p)k62) - k64)) ap(k62) - k64) )= 0 ap^ - kf )= 0 ap^? - k64) )=0 0 0
1 0 0 a(k66) - pk65)) 0 0
00
U
a(k9 ) - ki9 ))
(11)
Тотожнiсть (10) iз врахуванням виразу ll право! частини (11) приводить до трьох узгоджених систем piвнянь, iз яких знаходимо:
k (1) = k (3) 1 k (2) = k (4) 1 .
k6 = k6 I . k6 = k6 I. k(5) = 1 k(6)
k(1) = k(3) J ; k (2) = k(4) J; k6 = ~k6 .
k8 8 J k9 9 J p
0
1
0
0
0
Шсля врахування розв'язюв (12) систем (10-11) матриця Coeff набувае вигляду:
Coef =
~2 Pk (3)
a pk> ' —
г 6
- a 2 pk64)
-a 2 pk63) -a 2 pk64)
- a 2 pk66) + —
F 6 p +1
- a2pk66) +-S-
6 p +1
1
2a 0
1
2a 0
0 0
0 1
2a 0
a(1 - P)k63) a(1 - P)k64) a(1 - P)k63)
1
2a
— a(1 - p)k6'
(4)
0 a(1 - p)k6(6) +
(
0 a(1 - p)k6(6) -
pk6(3) + 1 6 2a2
( +1) _J_
( + 1)
pk6(4)
k (3)-
pk6
1
2a2
k(4)
pk6
pk66) -pk66) -
pk (4)-
1
6 2a2
pk6(3)
p'-™ + ¿г
a 2 (p + 1)
2 (p + 1)
p'66) -
1
p'66) - 2
a
a 2 (p + 1) p
(p + 1)
k6(3) 0 '83) 0
k6(4) 0 0 k9(4)
k6(3) 0 k8(3) 0
k6(4) 0 0 k9(4)
1 '66) p6 0 0 0
k6(6) 0 0 0
(13)
Знайдемо частину коефщенлв i3 умови, що сума вмх базисних функцш дор1внюе 1. Для цього знайдемо суми елеменлв матрицi (13) по стовпцям:
k(3) = -k(4) - Р+1 k(6)
2 p
£ NM, = 1 ^ i=1
-83) = 0 -94) = 0
(14)
Матриця (13) набувае нового вигляду з додатковими вимогами (14):
Coeff=
' a2 '64) + ^k66) J 1 2a 0
- a2p'64) 0 1 2a
a2{ p'64) + -66) ] -1 2a 0
- a2p'64) 0 -1 2a
_ a2'66) + J-I 6 p +1 0 0
- a2p'66) + Л-p + 1 0 0
a -ppi'[ p'64) + ^ '66)
.(4)
a(1 - p)-6
a p'64) + pT '66)
a(1-p)k6(4) a(1- p)'66) + 1
a(1 - p)k6
(6)
a(p + 1) 1
a(p + 1)
k(4) + £11 k(6) 1 + -L -I k(4) + ^k(6) 1 --64) -^-66) 0 0 0
45 2 *6 J ' 2a2
pk6(4)
p'64) + £
-I '(4) + ^ '66)J + ¿2 -['64) + ^ '66)
pk6(4)
k(6) -k6 2
1
p-4+I?
1
k(.6)-
pk6(6)
a" (( + 1)
p
a2 (( + 1)
a
(( + 1)
p'66) --
2p
k
(4)
000
- k64) - ^ k66) 0 0 0 2p
2((+1)
k(4) k6
k6
k(6)
000 000 000
Отже, з урахуванням вмх вимог, базиснi функцп бшра-мщи знаходяться за формулою:
NM = Coeff • Monom . (15)
Таким чином, невiдомими залишаються значення двох коефщенлв k((4) i k|(6), якi можна визначити, вико-
ристовуючи специфiчнi вимоги до базисних функцш.
Побудова базисш бммрашми i3 мШмальним слiдом мат риц! жорсткостi. У загальному випадку для набору базисних функцш {Ni} слiд матрицi жорсткост С з оди-ничною матрицею коефщенпв пружностi [5] обчис-люеться за формулою:
Trace=ЕШ
i=1 к
3N, дх
2
I ду
dN, ~3z
dV. (16)
М!н!м!зац!я слду матриц! жорсткосп для базису 6inipa-мщи, який отримано матричним способом. Шдставимо ба-зиснi функци (15) у формулу (16). Отримана функцш сл!ду
матриц жорсткосп е функцию двох незалежних змiнних k|(4) i
'66): Trace = Trace^^l'66)). За своею структурою вона
може мати тшьки мшмум, тому за необхщною умовою юнуван-ня екстремуму функци двох змшних знаходимо, що [8]:
(4) =__1_ 10p3 - p2 + 20p - 5 ;
6 4a2 15p4 + p3 +18p2 -3p + 5'
'Г = 4 62
10p3 - p2 + 20p -5
a2 15 p5 +16 p4 +19 p3 +15 p2 + 2 p + 5
. (17)
А мшмально можливий слщ матрицi жорсткост за-лежить вщ параметра видовження p за формулою:
a 230p6 + 462p5 + 653p4 + 620p3 + 372p2 + 214p185 оч
Trace = 15--(p + 1( 2 p + +1-. (18)
a
a
2
2
+
lYIiiiiMhaiiiii сладу матриц! жорсткост для базису бШрамвди, який отримамо шляхом вмутршмьоТ комдем-сащь Пiдставимо базиснi функцп (4, 6) у формулу (16). Отримана функцiя слiду матрицi жорсткостi е функщею
одше! незалежно1 змшно! a: Trace = Trace (a). Також i3 необхщно! умови iснування екстремуму функцп одше! змшно! знаходимо значення a та шших вагових ко-ефщенпв за формулами (7):
Таблиця 1 - Похибка обчислювального експерименту
Похибка по облает! бруса V Структура скшченно-елементно! реш1тки
Тетраедрально-октае^ральна Тетраедральна
Базисш функцп б1п1рам!ди (4), (6) Кусково-лшшш 6азисы1 функцИ октаедра [12]
p = 1 p = 0,7584
е1 0,73 0,69 0,73 0,48
е2 3,27 3,42 3,27 3,75
Р
a = —
10p 3 - p 2 + 20p - 5
dTrace da
= 0 ^ Р =
Y =
4 (5p2 + 2p + 5)2 - p +1) 5p4 + 2p3 - 2p2 + 2p + 5
(p +1)(5p2 + 2p + 5)p2 - p + 1J 5 p 4 + 2 p 3 - 2 p 2 + 2 p + 5 (5p2 + 2p + 5)p2 - p +1)^
p+1
. (19)
Обчислення ^ду матрицi жорсткост за формулою (16) iз базисними функцiями (4) i ваговими коефщента-ми (19) приводять до виразу, який ствпадае iз виразом (18) ^ду матрицi жорсткостi для базису, що побудова-ний матричним способом (15, 17).
Дослщження на екстремум функцп ^ду матрицi жорсткостi (теплопровщноста) для бiпiрамiди (18) пока-зують, що мiнiмум досягаеться при р « 0,7584, i складае приблизно 2,4776. Перевiримо ефективнiсть застосуван-ня отриманого значення параметру р стиснення бшра-мiди при розв'язаннi методом сюнчених елементiв за-дачi теорп поля.
4ЕКСПЕРИМЕНТИ Розглянемо прямокутний брус
V = {(х,у, г) :0 < х < а,0 < у < Ь, 0 < г < к}, де а, Ь, к > 0, який виготовлено з iзотропного матерiалу. Одна з його граней тдтримуеться при температурi Т = / (у), решта
граней - при температурi Т = 0.
Стащонарний розподiл температури задовольняе рiвняння Лапласа:
ДТ = 0, (20)
де Д = д2/дх2 + д2/ду2 + д2/дг2 - оператор Лапласа,
Т = Т (х, у, г) - температура у довшьнш точщ бруса V, з граничними умовами Дiрiхле:
Т\ х=0 = Т у=0 = Т у=Ь = Т|г=0 = Т\г=к = 0,
T| x=a = f (y) .
(21)
Аналiтичний розв'язок задачi (20-21), отриманий у робот [11], порiвняемо iз чисельним, знайденим методом сюнчених елементiв iз використанням решiток тетра-едрально! та тетраедрально-октаедрально! структур [12]. 5 РЕЗУЛЬТАТЫ
У табл. 1 представлен розрахунки похибки чисельно-го розв'язку по областi бруса V лшшних розмiрiв
a = b = 1, h = 2 при заданштемпературi f (y) = 10y(b - y), де t0 = 20 у середньоквадратичному:
1
X (t (x, y, const) - T (x, y, z))
i=1
та у метрищ простору неперервних на компакп V функцiй C (V ):
е2 = /max max (T(x, y,const)-T(x, y,z))
2 ^1<i<n (x,y, z
де T = T (x, y, const) i T = T (x, y, z) - точний та чисельний розв'язки задачу V - область сюнченого елемента
( n \
n - кiлькiсть сюнчених елеменлв.
n
V = U V i=1
При цьому у ансамблi iз лiнiйним тетраедром [7] тес-тувалися базиснi функцп октаедра (бшрамщи), заданi формулами (5), (7) для p = 1 та p = 0,7584, а також кус-ково-лiнiйнi функцп [12]. 6 ОБГОВОРЕННЯ
Аналiзуючи результати чисельного експерименту, слщ вiдмiтити, що використання бiпiрамiди iз коефщен-том видовження p = 0,7584 у ансамблi iз лшшним тетраедром е ефективним за точшстю обчислень порiвняно з аналопчною за структурою решiткою iз комiрками у формi октаедрiв при p = 1. При цьому неузгоджешсть [7] сюнчених елементiв рiзноl топологл виявляеться бiльш слабким критерiем вибору базисних функцш, нiж мiнiмальний слщ матрицi теплопровiдностi (жорсткостi). Найменша (у даному експерименл) похибка у середньоквадратичному отримана при використанш тетраед-рально1 решiтки за рахунок бшьш дрiбних комiрок та по-требуе бiльшого часу для реалiзацil програмного модуля, створеного безпосередньо авторами статл. В свою
чергу, розрахунок похибки у метрицi простору C (V) вказуе на перевагу використання альтернативно! реши-ки у МСЕ iз використанням октаедрiв при p = 1. Також вщзначимо, що значення
е 2 = max max
1<i<n (x, y, z )v-
(t (x, y, const) - T(x, y, z))
досяга-
ються на гранях скiнчених елементiв незалежно вщ струк-тури решiтки. Останне означае, що може бути доцшьним побудова базису бiпiрамiди бiльш високого порядку.
Таким чином, нами вперше знайдено СЕ, для якого найкраща (у середньо-квадратичному сенс!) точшсть прогнозуеться при вдаутност правильно1 геометрично! форми.
У програмних продуктах, що реал!зують МСЕ (зокре-ма, у ANSYS) для характеристики ком!рок мтки викорис-товують показник асиметрп Skewness, який розрахо-вуеться !з м!ркувань обрання за еталон СЕ правильно! геометрично! форми [9, 10]. Для обчислень з високими вимогами до точност отримуваних розв'язюв традицш-но вважають припустимими в!дхилення у межах ±0,1 вщ об'ему правильного геометричного тла, при невисоких вимогах до точност розв'язюв вважають припустимими в!дхилення у межах ±0,25 вщ того ж об'ему
Тод! по аналоги !з традицшними рекомендащями на практищ слщ використовувати бМрамщи з коефщентом видовження 0,51 < p < 1,01 при невисоких вимогах до точност отримуваних розв'язюв i 0,66 < p < 0,86 - при високих вимогах до точност отримуваних розв'язюв. ВИСНОВКИ
У робот! вперше узагальнено СЕ у форм! октаедра до СЕ бшьш довшьно! геометрично! форми, а саме, бМра-мщи, що дозволяе використовувати новий СЕ у частинах мтки, де вщбуваеться згущення ком!рок, або у пригра-ничнш зон!. Це дозволяе розширити можливост практичного застосування СЕ у форм! бшрамщи у сюнчен-но-елементних мтках пор!вняно !з обмеженнями застосування СЕ у форм! октаедра виключно в р!вном!рних сгтках.
Для нового СЕ у форм! бМрамщи вперше обгрунто-вано перелж вимог до базисних функцш, який забезпе-чуе !х однозначне визначення. Показано, що залучення традицшних вимог до базисних функцш, при побудов! базимв бМрамщи за двома тдходами: матричним способом i методом внутршньо! конденсацп, не дозволяе остаточно визначити значення вмх коефщенлв, як! вхо-дять до !х складу. Це дозволяе надавати базисам бМрам-щи додаткових доцшьних у МСЕ властивостей.
Вперше показано, що при залучент вимоги мшшзацл сл!ду матриц! жорсткоста, базиси бМрамщи, як! побудо-ваш за двома вказаними вище тдходами, ствпадають. Це дозволяе рекомендувати метод внутршньо! конденсацп як бшьш доцшьний для практичного використання !з м!рку-вань економп обчислювальних ресурмв.
Авторами вперше у теорп МСЕ знайдено елемент у форм! бшрам^и, якому за величиною сл!ду матриц! жорсткост прогнозуються кращ! апроксимацшш влас-тивост не при правильнш геометричнш форм! (октаедра), а при вдаиленш вщ не!. У робот! вперше теоретично розраховаш припустим! в!дхилення бМрамщи вщ форми октаедра. Вони виявилися неспод!ваними i такими, що не вщповщають традицшним уявленням про умови досягнення найкращо! точност на р!зних мтках СЕ. Знай-деш меж! змши параметра стиснення бМрамщи дозво-ляють будувати сюнченно-елементш мтки !з включен-
ням ком1рок у форм1 б1п1рам1д i3 задовшьними обчис-лювальними властивостями.
Результати обчислювального експерименту i3 вико-ристанням реш^ок на основi дослiджуваного елемента тдтверджують ефективнiсть застосування бiпiрамiд i3 критичним значенням коефщенту стиснення. Очевидно, що отриманий СЕ заслуговуе бiльш широкого застосування у практищ МСЕ. ПОДЯКИ
Дослiдження виконанш у вiдповiдностi до прюритет-ного напрямку розвитку науки i техшки в Украш на пер-юд до 2020 року: «1нформацшш та комушкацшш технологи» за темами: «Геометричне моделювання на диск-ретних елементах правильно! просторово! форми» (номер державно! реестрацп 0111U007990) та «Розробка методiв комп'ютерного дослiдження математичних моделей фiзичних полiв» (номер державно! реестрацп 0111U007951). Фшансування здiйснюеться за рахунок особистих коштiв виконавцiв.
СПИСОК ЛГГЕРАТУРИ
1. Алгоритм построения трехмерной адаптированной сетки для задач аэродинамики, решаемых методом конечных элементов / [Ю. А. Крашаница, А. В. Бахир, В. А. Тараненко, Ю. С. Мащенко] // Открытые информационные и компьютерные интегрированные технологии. - 2014. - № 66. - С. 105-110.
2. Greiner G. Hierarchical tetrahedral-octahedral subdivision for volume visualization / G. Greiner, R. Grosso // The Visual Computer. -2000. - I. 16. - Р. 357-369.
3. de Bruijn H. Numerical Method for 3D Ideal Flow [Electronic resource] / Han de Bruijn - Access mode: http:// hdebruijn.soo.dto.tudelft.nl/jaar2010/octaeder.pdf.
4. Мотайло А. П. Базисы шестиузлового октаэдра [Электронный ресурс] / А. П. Мотайло. - Материалы международной научно-практической конференции «Перспективные научные исследования - 2011». Серия: Математика: Прикладная математика (17-25 февраля 2011 г.). - София, Болгария. - Режим доступа: http://www.rusnauka.com/Page_ru.htm.
5. Сегерлинд Л. Применение метода конечных элементов / Л. Сегерлинд. - М. : Мир, 1979. - 392 с.
6. Шопов П. Й. Метод конденсации для задач механики несжимаемых флюидов / П.Й. Шопов // Сердика : Българско математи-ческо списание. - 1984. - Т. 10. - С. 198-205.
7. Зенкевич О. Метод конечных элементов в технике / О. Зенкевич. - М. : Мир, 1975. - 541 с.
8. Секулович М. Метод конечных элементов / М. Секулович. -М. : Стройиздат, 1993. - 664 с.
9. Checking the Skewness // ANSYS Icepak 12.1: User's Guide [Electronic resource]. - Access data: http://orange.engr.ucdavis.edu/ Documentation12.1/121/ICEPAK/iceug.pdf
10. ANSYS Fluent. [Electronic resource]. - Access data: https:// www.sharcnet.ca/Software/Fluent6/html/udfnode1.htm
11. Несис Е. И. Методы математической физики / Е. И. Несис. -М. : Просвещение, 1977. - 199 с.
12. Мотайло А. П. О численном решении стационарной задачи теплопроводности методом конечных элементов на решетке тетраэдрально-октаэдральной структуры / А. П. Мотайло // Научные ведомости БелГУ Математика. Физика. - 2014. -№25(196), Вып. 37. - Белгород: «НИУ БелГУ», 2014. -С. 119-127.
Стаття надшшлп до редакци 29.06.2016.
Шсля доробки 04.07.2016.
Мотайло А. П.1, Хомченко А. Н.2, Тулученко Г. Я.3
'Старший преподаватель кафедры высшей математики и математического моделирования Херсонского национального технического университета, Херсон, Украина
2Д-р ф.-м. наук, профессор, заведующий кафедры прикладной и высшей математики Черноморского государственного университета имени Петра Могилы, Николаев, Украина
3Д-р техн. наук, профессор, профессор кафедры высшей математики и математического моделирования Херсонского национального технического университета, Херсон, Украина
ПОСТРОЕНИЕ БАЗИСА БИПИРАМИДЫ
В статье бипирамида впервые рассматривается как 6-узловой конечный элемент (КЭ). Для построения ее биквадратичного базиса используются два разных подхода: матричный способ и метод внутренней конденсации базиса бипирамиды как 7-узлового лагранже-вого КЭ. Первый подход позволяет исследовать принципиально возможное количество базисов КЭ, а второй такой возможности не предоставляет, но является более экономичным. Показано, что после удовлетворения традиционных требований к базисным функциям в МКЭ в биквадратичных базисных функциях бипирамиды как 6-узлового КЭ, которые строятся с помощью названых ранее подходов, остается разное количество неопределенных коэффициентов. Эти коэффициенты далее используются для придания базисным функциям специальных свойств, адаптирующих их к решению граничных задач с уравнением Лапласа. В качестве критерия прогностического оценивания аппроксимационных свойств КЭ в форме бипирамиды выбрана величина следа матрицы жесткости. Минимизация следа матрицы жесткости приводит к построению одного и того же биквадратичного базиса бипирамиды при обоих подходах.
На основе полученного базиса анализируются границы допустимых деформаций геометрической формы бипирамиды. Впервые теоретически доказано, что существует КЭ, при использовании которого в качестве ячейки конечно-элементной сетки, наилучшая точность достигается при отклонении геометрической формы КЭ от правильного многогранника, в данном случае от октаэдра. Найдено критическое значение коэффициента сжатия, которое обеспечивает достижение минимума следа матрицы жесткости для бипирамиды исследуемой геометрической формы.
Проведен вычислительный эксперимент, результаты которого подтверждают теоретический прогноз свойств бипирамиды как КЭ. Установленные закономерности позволяют предполагать целесообразность применения базисов более высокого порядка для КЭ в форме бипирамиды.
Ключевые слова: метод конечных элементов, бипирамида, след матрицы жесткости, тетраэдрально-октаэдральная решетка.
Motailo A. P.1, Khomchenko A. N.2, Tuluchenko G. Ya.3
1Senior Lecturer in Chair of Higher Mathematics and Mathematical Modeling, Kherson National Technical University, Kherson, Ukraine
2Dr. Sc., Professor, Head of Chair, Applied and Higher Mathematics, Petro Mohyla Black Sea State University, Mykolayiv, Ukraine
3Dr. Sc., Professor, Professor of Chair of Higher Mathematics and Mathematical Modeling, Kherson National Technical University, Kherson, Ukraine
THE CONSTRUCTING OF BIPYRAMID'S BASIS
The bipyramid for the first time is considered as a 6-knots finite element (FE) in the article. For the construction of her biquadratic base two different approaches are used: matrix method and internal condensation method for the bipyramid's base as a 7-knots Lagrange FE. The first approach allows to investigate the fundamentally possible amount of bases for FE, and second approach does not give such possibility, but it is more economical. It is shown that after satisfaction of the traditional requirements to the basic functions in FEM in the biquadratic basic functions of bipyramid as a 6-knots FE, which are built by means of the named before approaches, there are a different amount of the indefinite coefficients. These coefficients are used future for the giving of the special properties to the basic functions, which adapt them to the solving of a boundary problems with Laplace's equation. The value of trace of stiffness matrix is chosen as a criterion for the prognostic evaluation of approximation properties of FE in bipyramid's form. The minimization of the trace value of the stiffness matrix results in the construction of the same biquadratic base of the bipyramid at the both approaches.
On the basis of the got base the borders of the possible deformations of the bipyramid's geometrical form are analysed. It is first well-proven in theory, that FE exists, at using of that as a cell of finite-element mesh, the best accuracy is arrived at the deviation of FE geometrical form from a regular polyhedron, in this case from an octahedron. The critical value of the aspect ratio which provides the achievement of a minimum of trace of stiffness matrix for the bipyramid of the studied geometrical form is found.
A calculable experiment the results of that confirm the theoretical prognosis for properties of bipyramid as FE is conducted. The found regularities allow to suppose the expediency of application of bases with higher-order for FE in form bipyramid.
Keywords: finite element method, bipyramid, trace of stiffness matrix, tetrahedral-octahedral mesh.
REFERENCES
1. Krashanitsa Yu. A., Bahyr A. V., Taranenko V. A., Mashchenko Yu. S. Construction Algorithm of the Three-Dimensional Adapted Mesh for the Tasks of Aerodynamics by Finite Elements Method, Open Information and Computer Integrated Technologies, 2014, No. 66, pp. 105-110.
2. Greiner G., Grosso R. Hierarchical tetrahedral-octahedral subdivision for volume visualization, The Visual Computer, 2000, I. 16, pp. 357-369.
3. de Bruijn H. Numerical Method for 3D Ideal Flow [Electronic resource]. Han de Bruijn Access mode: http:// hdebruijn.soo.dto.tudelft.nl/jaar2010/octaeder.pdf.
4. Motailo A. P. Bases of a 6-knots octahedron [Electronic resource], Proceedings of the International Scientific-Practical Conference «Future Scientific Researches — 2011». Series: Mathematics Applied Mathematics (2011 February, 17-25). Sofia, Bulgaria. Access mode: http://www.rusnauka.com/Page_ru.htm.
5. Segerlind L. Application of Finite Elements Method. Moscow, World, 1979, 392 p.
6. Shopov P. I. Condensation Method for Mechanics Problems of Compressible Fluids, Serdica: Bulgarian Mathematical Journal, 1984, Vol. 10, pp. 198-205.
7. Zienkiewicz O. Finite Elements Method in Technique. Moscow, World, 1975, 541 p.
8. Sekulovich M. Finite Elements Method. Moscow, Stroyizdat, 1993, 664 p.
9. Checking the Skewness ANSYS Icepak 12.1: User's Guide [Electronic resource]. Access data: http://orange.engr.ucdavis.edu/ Documentation 12.1/121/ICEPAK/iceug. pdf
10. ANSYS Fluent. [Electronic resource]. Access data: https:// www. sharcnet.ca/Software/Fluent6/html/udf/node1.htm
11.Nesis E. I. Methods of Mathematical Physics. Moscow, Enlightening, 1977, 199 p.
12. Motailo A. P. About the Numeral Solution of Stationary Task of Heat Conductivity by the Finite Element Method on the Mesh with Tetrahedral-Octahedral Structure, Scientific lists of BelSU. Mathematics. Physics, 2014, No. 25 (196), I. 37. Belgorod: «NRU BelSU», 2014, pp. 119-127.
