О численном решении стационарной задачи теплопроводности методом конечных элементов на решетке тетраэдрально-октаэдральной структуры Текст научной статьи по специальности «Физика»

О ЧИСЛЕННОМ РЕШЕНИИ СТАЦИОНАРНОЙ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ НА РЕШЕТКЕ ТЕТРАЭДРАЛЬНО-ОКТАЭДРАЛЬНОЙ СТРУКТУРЫ

А.П. Мотайло

Херсонский национальный технический университет, Бериславское шоссе, 24, Херсон, 73008, Украина, e-mail: akilehzria@ukr.riet

Аннотация. В работе установлена сходимость метода конечных элементов на решетке тетраэдрально-октаэдральной структуры с применением кусочно-линейных координатных функций шестиузловсн'о октаэдра. На примере стационарной задачи теплопроводности показана эффективность использования октаэдра в качестве ячейки пространственной решетки.

Ключевые слова: октаэдр, кусочно-линейные узловые функции, тетраэдрально-октаэдральная решетка.

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

Анализ предшествующих публикаций, цели статьи. Авторы |1| используют модель октаэдра с семью узлами интерполяции (рис.1) и кусочно-линейные узловые функции:

= 1 — |x| — |y| - |z|; ^1,3 = (x ± |x|)/2; ^2,4 = (y ± |y|)/2; ^5,6 = (z ± |z|)/2. Автором статьи построены несколько различных систем узловых координатных функций для шеетиузловой модели октаэдра (без центрального узла), изучены интерполяционные качества каждого набора функций |2|, В работе |3| установлено, что октаэдр с кусочно-линейными узловыми координатными функциями интерполирует точно линейные ноля и может быть использован для решения краевых задач математической физики с дифференциальным оператором второго порядка.

Цель статьи - на примере задачи о распределении температуры в сплошной среде установить сходимость метода конечных элементов на решетке тетраэдрально-октаэдральной структуры (рис. 2) с использованием кусочно-линейных узловых координатных функций октаэдра, результаты проверить численным экспериментом.

Основная часть. Рассмотрим прямоугольный брус V = {(x,y,z) : 0 < x < а, 0 < У < b, 0 < z < h}, оде a,b,h > 0, изготовленный из изотропного материала, одна из

граней которого поддерживается при заданной температуре Т = f (у), остальные — при температуре Т = 0 (рис.3).

У

Рис. 1. Октаэдр |х| + |у| + \z\ < 1.

Рис. 2. Фрагмент решетки тетраэдрально-октаэдральной структуры.

Стационарное распределение температуры в пространственной области V удовлетворяет уравнению:

ДТ = 0, (1)

где Д = д2/дх2 + д2/ду2 + д2/дг2 — трехмерный оператор Лапласа, Т = Т(х,у,г) -температура в произвольной точке (х, у, г) брус а V, с граничными условиями Дирихле:

T = T =T =T = T

x=0 y=0 y=b z=0

z=h

0; T

= f (y)

(2)

Аналитическое решение задачи (1), (2) найдено в сечении Vxy = {(x,y) : 0 < x < а, 0 < y < b} брус a V плоскостью z = const [4]:

T(х,У) = Y^f

sh(nnx/b) = J'1 sh(nna/b)

sin(nny/b) ,

где

= ^ ] ${у)ы&(жпу/Ъ)(1у. 0

В эквивалентной трехмерной вариационной постановке МКЭ задача (1), (2) сводится к нахождению функции Т(х,у,г), минимизирующей функционал

X =

(дТ/дх)2 + (df/dy)2 + (dff/dz)2) dxdxydz

где Т = Т(х, у, г) - функция множества допустимых в области V пробных функций. При этом, согласно [5], функция Т(х, у, г) является допустимой в области V, если Т(х, у,г) непрерывна и имеет кусочно-непрерывные частные производные в этой области.

Рис. 3. Брус V = {(х, у, г) :0 < х < а, 0 < у < Ь, 0 < г < К}.

Разобьем область V на / конечных элементов Пг в форме октаэдров и тетраэдров.

I

Тогда элементный вклад х°г в функционал х = ^^ Х- определяется равенством:

Г=1

Х- =11 I ((дТПг /дх)2 + (дТПг/ду)2 + (дТ- /дг)2) дьхдьхудьх ,

где Пг — октаэдр ПО или тетраэдр ПГ

г

г-,

Тпг I ЕГ=1 ^гТг , (х, у,г) е Пг;

\о , (х, у, г) е Пг

- пробная функция, N = А"г(х, у, г) — узловые координатные функции Пг, Тг Т(хг,уг,гг) — значения температуры Т в узлах (хг, уг,гг) КЭ Пг, пг — число узлов П причем пг = 6 если Пг = П^, и пг = 4, если Пг = ПГ-Положим, что

лт, ч /л^(х,у,л), (х,у,л) е г =Т7б; А'¿(х,у,л)= < __(.3)

[Мх,у,г) , (х,у,г) е ПГ = 1> 4,

где ЖЬг(х,у,г) — кусочно-линейные функции [2], соответствующие узлам октаэдра, а Ьг(х,у,г) — линейные функции, объёмные координаты тетраэдра [4]. При этом (х,у,г) получаются из функций

(1 + 2|£|± 3£ - |п| - К |)/6 , г = 1, 3; ^г(С, П, С) ={ (1 + 2|п| ± 3п - |С| - К|)/6 , г = 2, 4; (4)

(1 + 2|<|± 3( — | — |п|)/6 , г = 5, 6

заменой (x,y,z)T = [J] ■ (£,ц,()T + (x0,r, y0,r,z0,r)T, где [J] — матрица преобразования локальной системы координат O^nC в глобальную Oxyz, (x0,r , y0,r, z0,r) - координаты начала O£n( для октаэдра ПО в системе Oxyz.

Утверждение. Функция T(x,y,z) = NjTj, где N = Ni(x,y,z) определяются равенствами (3), (4), n — число узлов разбиения V = (Jr=1 Пг, является допустимой в области V.

□ Покажем, что функция T(x, y, z) = ^™=1 NT непрерывна в области V = (Jr=1 Пг. Действительно, из определения функции TПг (x, y, z) в области V следует, что

Т(х, у, z) = ТПг(х, у, z), (х, у, z) EÜr, i = 1,1.

Из равенств (3)-(4) следует, что для любого г = 1,1 функция ТПг(х, у, z) непрерывна в области Пг:

1) если Пг = ПГ, то TПг (x, y, z) = a1r) + a2r)x + a^y + a4r)z непрерывна как линейная в области конечного элемента для любого г = 1,1\,

2) если Пг = ПО, то

Tnr(x,y, z) = ßr) + ß2r)x + ß^y + ß4r)z + ß5r)|x - xo,r| + ß6r)|y - yo,r| + ß7r)|z - zo,r|

непрерывна как кусочно-линейная в области для любого г = 1,1-2 (а\Г\ ßj^ ~ действительные числа, которые однозначно определяются координатами узлов КЭ Пг [6, 7],

/i + /2 = 0. ' _

В межэлеменшой зоне, т.е. на гранях смежных элементов, для любых г, s = 1,1

равенство TQr(x,y,z) = TQs (x, y, z) справедливо на непустых множествах ПГ П П [6]

i

и ПО П ПО [3]. Таким образом, функция TПг(x, y, z) непрерывна в области V = Пг.

г=1

Последнее означает, что частные производные функции TПг (x, y, z) в области V существуют и определяются равенствами:

<]Г {x,y,z)eü\,i--

dx \ ß{f + ßf ■ sgn(x - x0<j), (x, у, z) e Щ, j

dT _\ Q'f, (x, y,z) eQl,i =

дУ \ß? + ß(63) ■ sgn(у - y0j), (x, y,z)e Щ, j =

dT_(af, {;x,y,z) eü\,i =

öJ " | ßf + ß? ■ Sgl!(z - Zo,j) , (X, y, z) еЩ,з =

h;

h;

I1;

h

l2;

При этом dT/dx, dT/dy, dT/dz непрерывны в области V \ П°! и кусочно-непрерывны

r=1

12

в области П^ т.е, dT/dx, dT/dy, dT/dz кусочно-непрерывны в области V. U

Заметим, что условие допустимости функции T(x, y, z) в области V является условием сходимости МКЭ к точному решению в рассматриваемой задаче |6, стр.3211. При этом сходимость МКЭ означает, что при неограниченном уменьшении размеров КЭ величина х стремится к нулю.

Минимизация функционала х = Г= Xn предполагает дифференцированне по T (г = 1,/г), как по неизвестной величине, функции х(Т\, Т2,..., Тп), связанное с вычислением интегралов по каждому элементу Пг:

дХПг 1 [[[ (dNp dNq dNp dNq dNp dNA nr

где p, q — глобальные номера узловых координатных функций КЭ Пг,

Процедура вычисления интеграла (5) по объему тетраэдра ПГ описана в [6, 8]. Интегрирование по октаэдру ПО выполним в локальных координатах:

Г [Г (dNp dNq dNp dNq dNp dNq\

/ / / ---+ ---+ ---dxdydz =

JJJ \ dx dx dy dy dz dz J n°

dNp dNq + dNp dNq + dNp dNq dx dx dy dy dz dz (no)*

где (ПО) * _ октаэдр ПО в системе координат O£nC и

/щ>\

dx dNp ду dNp V dz )

/ dx dx 5x\

д£ dn 5C

ду <9у <9y

dn <9(

dz dz Sz

5C

fdNp\ <9? dNp

dNp

\9( J

= ДО]

| J ,

При этом, учитывая кусооновнепрерывный характер частных производных функций А^, область октаэдра (ПО) * может быть представлена в виде (ПО)* = и8=1 П3, где П — равновеликие по объему прямоугольные тетраэдры, внутри которых переменные С, п, С сохраняют знак. Например, интегрирование по области П1 = {(С,П,С) : 0 < С < 1, 0 < п < 1 - С, 0 < С < 1 - С - п} гДе С, П, С > 0 а функции дАр/дС, дАр/дп, дАр/д( непрерывны, принимает вид:

•/У,/ V дх дх ду ду дг дг )

Пх

1-?

где

1

= у d£JdпJ т1 [З]1 ЛБф № 0 0 0

5 -1 -1 -1 -1 -Г

В\ = - 1—1 5 -1-1 -1 -1 6 1 -1 -1 -1 -1 5 -1,

- матрица частных производных узловых координатных функций октаэдра ПО, определенных в области П1.

Очевидно, что подынтегральное выражение не зависит от п, С ПРИ этом объем тетраэдра равен

1 1-? 1-?-ч

I I ^ I ¿С = 1. 0 0 0

Аналогично можно получить объемные интегралы но ПЛ,, 5 = 2,8. Тогда

где

Бо

дх дх ду ду дг дг )

6

1

6

8=1

1 1 -5 1 1 1 -1 5 -1 -1 -1 -1 -1 -1 -1 -1 5 -1

Бз

1

6

11 11 -1 -1

-5 1 1 1 -5 1

-1 -1 5 -1

Б4

Бб

1

6

1

6

5 1 -1

-1 -1 -1 -1 -1 1 1 -5 1 1 -1 -1 -1 5 -1

1

1

1 -5 1 1 1 15 -1 -1 -1 -1 1 1 1 1 -5

Б..

б7

1

6

1

6

5 -1 -1 5 11

-1 -1 -1 -1 11

1 1 -5 1 1 1 1 1 -5 1 111 11 -5,

-1 -1 -1 -1 1 -5

1 1

Б«

1

6

5 -1 -1 -1 -1 -1 1 1 1 -5 1 1 1 1 1 1 1 -5

Полученная формула (6) является формулой численного интегрирования но октаэдру ПО-

Условия минимума дх/дТ = 0 функционача х приводят к решению системы линейных алгебраических относительно [Т] = (Т1,Т2,...,Тп) уравнений [к][Т] = 0 где [к] -глобальная матрица теплопроводности, элементы которой определяются равенствами Нц = (дхПг/дТг)1 ЬЗ = 1)«- Удовлетворяя главным граничным условиям (2) путем корректирования отдельных строк матрицы [к] и нулевого вектора в правой части системы, получаем численное решение в виде Т(х, у, г) = ^™=1 ^Т^,

Д с

Рис. 6. Температура бруса V вдоль осей сечений К = г/2, К = г/4, К = г/8: ■ ■ ■ — аналитическое решение (с точностью до 10-6); х — решение МКЭ на решетке тетраэдрально-октаэдральной структуры; о — решение МКЭ на решетке тетраэдральной структуры.

Дня нахождения численного решения МКЭ задачи (1), (2) с использованием тетраэдральной и тетраэдрально-октаэдральной решеток автором статьи составлены программы в среде Maple. На рис. 4 значения температуры бруса найдены в направлениях осей y = b/2 и x = а/2 в сечениях z = h/2 (4a, 46); z = h/4 (4в, 4r); z = h/8 (4д, 4e). Результаты численного эксперимента соответствуют размерам а = 1, b = 2, h =10 и условию на границе: f (y) = t0y(b — y) для t0 = 20.

Анализируя результаты численного эксперимента, заметим, что решение задачи, иолучешюе на решетке тетраэдрально-октаэдральной структуры (8864 тетраэдра, 344 октаэдра, 2121 узел интерполяции), достаточно близко к точному и решению МКЭ с использованием тетраэдральной решетки (24576 тетраэдров, 4913 узлов интерполяции). При этом время работы программы (в секундах) в нервом случае 2621.05, во втором - 6402.11. Последнее означает, что, используя октаэдр с кусочно-линейными узловыми координатными функциями, в рассматриваемой задаче можно сократить время и объем вычислений.

Выводы. В работе установлена сходимость МКЭ па решетке тетраэдралыю-окта-эдралыюй структуры с применением кусочно-линейных координатных функций шести-узлового октаэдра; получены формулы численного интегрирования но октаэдру и решение МКЭ стационарной задачи теплопроводности дня прямоугольного изотропного бруса. Результаты численного эксперимента свидетельствуют о том, что использование октаэдра в качестве ячейки пространственной решетки в данной задаче оптимизирует время и объем вычислений, сохраняя достаточно высокую их точность.

Литература

1. Grosso R., Grcincr G. Hierarchical Meshes for Volume Data Computer Graphics International 1998 (CGI'98). 1998. P.761-771.

http://dblp.uni-tricr.dc/db/conf/cgi/cgil998.html // GrossoG98.

2. Мотайло А.П., Хомченко А.И. Порпзняльний ана;ш базисгв октаедра /7 Матер1али. IX м1жнародно1 науково-практичноТ конференцп «Найновшп науков1 доеягнення - 2013». Серш: Математика: Фкадка. Сучаеш шформащйш технологи' (7-15 березня 2013 р.). Соф1я, Болгар1я: «Вял ГРАД-БГ» ООД, 2013. Т.21 С.28 33.

3. Мотайло А.П., Хомченко А.Н. Интерполяция кусочно-линейными функциями на решетках тетраэдрально-октаэдральной структуры /7 Математичне та комп'ютерне моделю-вання. Серш: Фкадко-математичш науки: зб. наук. пр. №8. Кам'янець-Подгл. нац.ун-т ¿м. I.Опенка; 1н-т кибернетики iM. В.М.Глушкова, Нац.акд.наук Укра'ши. Кам'янсць-Подьл.: К-ПДУ ¿м. 1вана Опенка, 2013. С.139-150.

4. Несис Е.И. Методы математической физики / М.: Просвещение, 1977. 199 с.

5. Норри Д., де Фриз Ж. Введение в метод конечных элементов / Пер. с англ. / М.: Мир, 1981. 304 с.

6. Зенкевич О. Метод конечных элементов в технике / М.: Мир, 1975. 541 с.

7. Хомченко А.Н., Мотайло А.П. Две модели кусочно-линейной интерполяции на октаэдре // Проблем« шформацшних технологи!: зб. наук. пр. №1. Херсон: ХНТУ, 2011.

С.' '47 50.

8. Сегерлинд Л. Применение метода конечных элементов / М.: Мир, 1979. 392 с.

ABOUT NUMERICAL CALCULATION OF STATIONARY HEAT CONDUCTION PROBLEM BY THE FINITE ELEMENT METHOD ON THE LATTICE TETRAHEDRAL-OCTAHEDRAL STRUCTURE

A.R Motailo

Khersori National Technical University, Beryslav Highway, 24, Kherson, 73008, Ukraine, e-mail: akilehzna@ukr.net

Abstract. Convergence of numerical approximations on the basis of the finite element method is established when the lattice of tetrahedral-octahedral structure is used by six nodal octahedron with pieeewise linear functions. The efficiency of octahedron as the cell of space lattice is shown on the basis of the example of stationary heat conduction problem.

Key words: octahedron, pieeewise linear nodal function, tetrahedral-octahedral lattice structure.