Розрахунок двовимірних стаціонарних температурних полів за наявності внутрішніх джерел тепла методом скінченних елементів Текст научной статьи по специальности «Математика»

ISSN 1994-7836 (print) ISSN 2519-2477 (online)

@ EE3 Correspondence author V. P. Karashetskyy volodymyr10@gmail.com

М. I. Дмитрусь, В. П. Карашецький

Нацюнальний лкотехшчний утверситет Украши, м. Львiв, Украта

РОЗРАХУНОК ДВОВИМ1РНИХ СТАЦ1ОНАРНИХ ТЕМПЕРАТУРНИХ ПОЛ1В ЗА НАЯВНОСТ1 ВНУТР1ШН1Х ДЖЕРЕЛ ТЕПЛА МЕТОДОМ СК1НЧЕННИХ ЕЛЕМЕНТ1В

Сформульовано краеву задачу розрахунку стацiонарного температурного поля за наявност внутрiшнiх джерел тепла, що описуеться диференцiальними рiвняннями. Для побудови скiнченно-елеменшоl моделi розрахунку розподшу температури всерединi двовимiрноl обласп, заповнено1 нелiнiйними безгiстерезисними анiзотропними середовищами, врахування зада-них вздовж ll границi граничних умов Дирихле та одш^дних граничних умов Неймана, використано лагранжевi трикутни-ки и-го порядку. Застосовано кубатурну формулу чисельного штегрування по площi трикутника на основi iнтерполяцiйного повного полшому для лагражевого трикутника другого порядку. З умови мшмуму функцюнала отримано рiвносильну нель нiйну систему алгебра1чних рiвнянь, яка розв'язуеться, як правило, ггерацшним методом Ньютона. Подано алгоритм визна-чення внеску кожного внутрiшнього скiнченного елемента (СЕ) у вектор нев'язок i матрицю Якобi системи рiвнянь з вико-ристанням елеменпв тензора диференщально! теплопровiдностi середовища. Наведено також алгоритм визначення внеску кожного граничного СЕ у вектор нев'язок i матрицю Якоб^ якщо вiн мае ильки один чи кiлька граничних вузлiв iз гранич-ними умовами Дирихле або мае два граничних вузли з однорщними граничними умовами Неймана i один граничний вузол iз граничними умовами Дирихле.

Krnuoei слова: стацюнарне температурне поле; внутрiшнi джерела тепла; теплопровдають; лагранжевий трикутник; метод скшченних елементiв; кубатурнi формули; метод Ньютона; тензор диференщально! теплопровiдностi; граничш умо-ви.

Вступ. З iнтенсивним розвитком комп'ютерних тех-нологiй особливого значення набувае в наш час матема-тичне моделювання рiзних фiзичних процеав. Для розв'язання задачi стацюнарно! теплопроввдносп у дво-вимiрнiй областi використовують рiзнi чисельнi методы, серед яких найбiльший розвиток отримав метод скшченних елеменпв (МСЕ).

Не зважаючи на велику к1льк1сть публiкацiй з МСЕ, у них ввдсутнш опис алгоритму формування системи алгебрачних рiвнянь з умови мiнiмуму функцiонала для розрахунку розпод^ температури всерединi дво-вимiрноl областi, заповнено! нелiнiйними безпстерезис-ними анiзотропними середовищами при наявносп внут-рiшнiх джерел тепла, та урахування заданих вздовж ll границ граничних умов Дирихле i однорщних граничних умов Неймана.

Актуальним завданням залишаеться також викорис-тання в якостi ск1нченних елементiв лагранжевих три-кутник1в 2-го порядку з метою тдвищення точностi розрахунк1в.

Мета роботи полягае в розробленш методики розрахунку двовимiрних стацiонарних температурних полiв за наявностi внутрiшнiх джерел тепла методом сшнчен-них елеменпв, яка б дала змогу враховувати розподш

1нформащя про aBTopiB:

Дмитрусь Мирослав 1ванович, acnipaHT кафедри шформацшних технологiй. Email: myroslavdm@gmail.com Карашецький Володимир Петрович, канд. техн. наук, доцент. Email: volodymyr10@gmail.com

Цитування за ДСТУ: Дмитрусь М. I., Карашецький В. П. Розрахунок двови/^рних стацiонарних температурних полiв за наявност

внутрiшнiх джерел тепла методом скшченних елеменпв. Науковий вкник НЛТУ УкраУни. 2017. Вип. 27(5). С. 134-138. Citation APA: Dmytrus, M. I., & Karashetskyy, V. P. (2017). Calculation of two-dimensional stationary thermal fields with internal heat sources of the finite element method. Scientific Bulletin of UNFU, 27(5), 134-138. https://doi.org/10.15421/40270527

НАТЫ

УКРЛ1НИ

Hl/IUB

Науковии BicHMK НЛТУУкраТни Scientific Bulletin of UNFU

http://nv.nltu.edu.ua https://doi.org/10.15421/40270527 Article received 22.06.2017 р. Article accepted 29.06.2017 р.

УДК 536.21

температури всередиш двовимiрно! обласп, заповнено! нелшшними безгiстерезисними анiзотропними середовищами.

Виклад основного матерiалу. Для краево! задачi розрахунку стацiонарного температурного поля за наяв-ностi внутрiшнiх джерел тепла, що описуеться дифе-ренцiальними рiвняннями (Yushko, Borshch, & Yushko, 2011):

divB(H) = J , (1)

Н = -gradU , (2)

у плоскш обласп D функцiонал F набувае вигляду

F = |{(¥ - СУЗ , (3)

X

Н

де: W = |BdH ; (4)

0

С = Ш ; (5)

де: Н , В - розмщет у площиш Б вектори напруже-носп температурного поля та густини теплового потоку; и - температура в будь-якш точщ обласп Б ; J -нормальна до площини Б проекщя вектора густини теплового потоку внутрiшнiх джерел тепла; X - площа областi Б.

= 0.

(6)

Розподш температури U всередиш обласп D, що MiHiMi3ye функцiонал F, забезпечуе розв'язання краево! задача Умова мшмуму функцюнала (3) набувае такого вигляду

dF dU '

Для побудови сшнченно-елементно! моделi заповни-мо область розрахунку D сукупшстю лагранжевих три-кутнишв n -го порядку (Karashetskyi, 2007).

Нехай внаслвдок трiaнгyляцп двовимiрноl обласп розрахунку D отримуемо M лагранжевих сшнченних елементiв (СЕ). Кожному з них присво!мо порядковий номер m (m = 1,M ) i локальну нyмерaцiю вyзлiв, згiдно з якою i -му вузлу m -го СЕ вiдповiдaе номер mi. Для вае! облaстi розрахунку встановимо сiтковy (нас^зну) нyмерaцiю R внyтрiшнiх вyзлiв i G граничних вyзлiв. Поточнi значения порядкових номерiв внyтрiшнiх вуз-лiв позначимо r.

Для сшнченно-елементно! облaстi фyнкцiонaл F з урахуванням (3)-(5) набуде вигляду

F = Y Fm = Y (Wm - Cm ),

де:

W

Cm

-- J WdS;

Sm

J CdS ;

(7)

(8) (9)

де: Sm - площа m -го СЕ, що визначаеться через коор-динати його вершин у прямокутнш системi координат за формулою

Sm = idet m 2

1 x1 y1

1 x2 y2

1 x3 y3

(10)

Утворимо R ^рний вектор-рядок i вектор-стовпець температури U у внутршшх вузлах:

U = (U1,...,UR); U* = (U1,...,UR)*.

(11)

Умова мiнiмyмy фyнкцiонaлa F з урахуванням (6) рiвносильнa нелшшнш системi алгебра!чних рiвнянь

- - dF *[U.]=dF- = 0.

dU,

(12)

Застосуемо для (8), (9) кубатурну формулу чисель-ного штегрування за площею лагранжевого трикутника (Karashetskyi, 2007). Наприклад, у рaзi використання лагранжевих трикутник1в другого порядку (n =2), кшь-к1сть вyзлiв у яких p=6, отримаемо

1 6

Fm = 3 Sm XWmi — Cj ),

i=4

Hn,

де:

Wm, = J BdH ;

(13)

(14)

k * =

1 xm1 ym1 xm1 ym1 xm1 У2т1

1 xm2 УП2 xm2ym2 xm2 vJi

1 xm3 Уп3 xm3ym3 xm3 VJ3

1 xm4 ym4 xm4ym4 x2 xm4 x2 xm4

1 xm5 Ут5 xm5ym5 xm5 УП5

1 xm6 Ут6 xm6Vm6 xm6 Ут6

(19)

• координатна матриця, рядки яко1 е координатними векторами вигляду (18) у вузлах m -го СЕ; Um*, к,, km -вдаовдао вектори-стовпцi i матриця, трaнспоновaнi вiдносно Uт, к , km,.

У локальнш прямокyтнiй системi координат вектор H напруженосп температурного поля пов'язаний iз температурою U спiввiдношенням

я „, dU-г- dU

H = —gradU =--i--j .

dx dy

(20)

Проекцп Hx, Hy вектора H у mi -му вyзлi з ураху-ванням (16) i (20) набувають вигляду:

H

dU

ix v = —U k(x) = —K(x)U •

-v mi^Smi m mi* mi m* ">

dx

dU - - - -

H = X v = —U K(y) = — K(y>U '

ymi ~ \xmiiymi U mKmi* Kmi U m* >

dy

де:

K(x) = к(x,k-\; K(y' = к{У'к",;

mi mi m* > mi mi m* >

Kx* = k~lk(x*; K(У* = к~1к{У*;

mi* m mi* 1 mi* m mi* > - дк

kw =—\x .,y . = (0,1,0,y .,2x .,0); i i i i i dx

- дк

кw =—\x .,y . = (0,0,1,x . ,0,2y .);

i i i i i dy

(21) (22)

(23)

(24)

(25)

(26)

к(У* - стовпцi, отримаш транспонуванням рядк1в (25), (26).

Диференцшючи вираз (13) по вектору U m i врахову-ючи (21), (22), отримуемо

фт, = -^-J =1 Sm Y-^L J ((dH )B — UmiJmi) = Ут dUm 3 mi=4dUm J ^ ' mi

=1S Yf— J U

3 m Y f dUn mi j dU,

= 1 s y fd-H^B + HlB — J dU

3 m yi dUm ^ j ■

dUm dUm

16 - - 1 -

= — 3Sm Y (Kj1*bxj + KJ*)Bymi ) — 3SmJm 3 i=4 3

(27)

C. = U.J.. (15)

mi mi mi \ '

Представимо зaлежнiсть температури U в межах m -го СЕ повним полшомом другого степеня

U = UmKlk, = ^JUJ,, (16)

де: Um = (Um1,...,Uj6) (17)

• вектор-рядок значень температури U у вузлах m -го СЕ;

к = (1, x, y, xy, x2, y2) (18)

• координатний вектор-рядок поточно! точки з координатами x , y ;

де: Jm* = (Jm1,-, Jj3)* (28)

- вектор-стовпець заданих значень J у вузлах m -го СЕ; Bmi, Bmi, Bymi - вiдповiдно вектор густини теплового потоку i його склaдовi в mi -ому вузл^ як1 визнача-ються за значеннями проекцiй (21), (22) вектора напру-женостi i характеристикою теплопровiдностi нелшшно-го безгiстерезисного середовища, що виражаеться век-торним рiвнянням або двома скалярними рiвняннями:

в = B[H ]; (29)

Bx = Bx[Hx, Hy ]; By = By[Hx, Hy ]. (30)

Нелшшну систему рiвнянь (12) розв'язують, як правило, ггерацшним методом Ньютона.

Для визначення внеску m -го СЕ у систему рiвнянь (12) потрiбно:

m=1

n=1

S

• знайти вектор фт* на кожнш гтераци за формулою (27);

• за таблицею в1дпов1дност1 локально! i мтково! нумерацп встановити номери г вузл1в, яю зб1гаються з вузлами т1,...,тр ;

• кожний елемент вектора фт*, який в1дпов1дае г -му внутршньому вузлу, внести в1дпов1дно в г -те р1вняння системи (12).

Повну систему рiвнянь (12) отримаемо, виконавши цю процедуру для всiх М елеменпв. У цьому разi вик-ладена вище процедура використовуеться на етат фор-мування вектора нев'язок. Дещо трудомiсткiшою опера-цiею е складання для векторно! функцп ф* = (ф,...,фЯ)* матрицi Якобi р розмiрностi Я х Я . Виведемо загаль-нi вирази, що використовуються для ще! мети.

Диференщюючи вираз (27) за вектором ит*, отри-муемо матрицю розмiрностi 6*6

(y) dBymi )

_ Лфт, _ 1 S V|K(x) dBxmi + K

dUm, 3 f_í dUm, dUm,

З урахуванням (21), (22), (30) i (31) маемо

(31)

1 6 _ _ _ Vm _ -3Sm L^'m^xxmiK'mi' + ^xymi^i' ) + ,4

3 i_4 (32)

+K(y)(K K(x) + K K(y)))

mi*y yxmiKmi yymiKmi ))

dH

Де: k(x> _ dHm,K(y)

dH„

dUm* dUm*

тензора диференщально! теплопровiдностi середовища

dB 3B'

к_ dB _

dH

дИх dHy

dBy B dH 3H„

К Kx¡

К К

(33)

обчислюваш в i -ому вузлу m -го СЕ.

Для безгiстерезисного середовища на основi теоре-ми взаемносп (Silvester, Cabayan, & Browne, 1973)

т]'. Повну матрицю Якобi р отримаемо, виконавши цю процедуру для кожного з М синченних елеменпв об-ластi розрахунку Б.

У разi використання лагранжевих трикутник1в 1-го, 3-го i 4-го порядив, потрiбно застосувати вщповщщ ку-батурнi формули чисельного iнтегрування (Karas-hetskyi, 2007) i провести вивщ основних залежностей за викладеною вище методикою.

Уздовж гранищ областi Б повиннi бути задаш гра-ничнi умови Дирихле (значення потенцiалу и ) або од-норiднi граничнi умови Неймана

H _ о,

dn

(36)

де: Нп - нормальна складова вектора Н на одиничний вектор п зовшшньо! нормалi до гранищ обласп Б .

Для визначення внеску кожного т -го СЕ у вектор нев'язок ф* i матрицю Якобi р, якщо вш мае один чи калька граничних вузлiв iз граничними умовами Дирихле, потрiбно на кожнiй ггерацп для кожного вузла тБ iз граничними умовами Дирихле враховувати, що значення итБ задане i постiйне, тому:

dF

. т

dU„,

_ 0 ;

; Km (j, k _ x, y) - елементи

VmDp

VmpD

дищ ф

dUm

-_ о, p _ 1,6;

-_ 0, p _ 1,6.

(37)

(38)

(39)

Л = Л , тому Л . = Л ..

xy yx > J xymi yxm

У pa3Í Í30TponHoro нелiнiйного середовища тензор диференщально! теплопровщносп визначаеться за формулою (Dyshovyj, 1983).

x = \^pc0slnx + Лsin2Пх (Лр - A)cosncosПуJ (34)

[(Лр- Л )cosny cosn ЛР cos2 Пу + Лsin2 ПуJ'

де: Л = dBB, Лт = B - вiдповiдно pадiальна диференщ-

dH H

альна i тангенщальна теплопровщшсть середовища; П (l = x,y) - кути мiж вектором H або B[H] i вщповщ-

но ортами l, j локально! декартово! системи координат.

Для лшшного iзотpопного середовища маемо таку piвнiсть Лр=Лг= B/H = Лс, тому тензор диференщально! теплопроввдносп набувае такого вигляду

Л = {Л Л}. (35)

Для того, щоб визначити вклад m -го СЕ в матрицю Якобi (р, потpiбно обчислити матрицю q>m на кожнiй перацп за формулою (32) i тдсумувати всi !! елементи з ввдповвдними елементами матpицi р, враховуючи, що елемент pmj належить ns -й клгтиш матpицi р, де n , s - сiтковi номери вузлiв з локальними номерами mi i

Розглянемо визначення внеску кожного т -го СЕ у вектор нев'язок ф* i матрицю Якобi р, який межуе своею стороною iз границею обласп Б i мае тiльки два граничних вузли з граничними умовами Неймана, ос-ильки для СЕ п -го порядку ильисть Рн таких вузлiв

визначаеться за формулою

Ры = п . (40)

Умова Неймана (36) на гранищ обласп Б набувае вигляду

Н„ = Нп = Нп + НуПу = 0, (41)

де пх, пу - проекцп одиничного вектора п дотично! до

гранищ обласп Б.

З урахуванням (21), (22) запишемо (41) для кожного iз двох вузлiв iз граничними умовами Неймана т -го СЕ такий вираз

" (42)

(nxKmX), + nyK^Um- _ 0,

де

K(x) _

KmN*

K(x)

KmN,

K

(x)

К (y),

K (y)

KmN,

K

(y)

(43)

пpямокутнi матpицi pозмipностi 2x6, рядками яких е вектори-рядки, визначенi за виразами (23) для вузлiв з граничними умовами Неймана m -го СЕ; N1, N2 - вiд-повщно початкове i к1нцеве значення локальних номе-piв вузлiв з граничними умовами Неймана m -го СЕ.

Представимо (42) у виглядi

КРт, = 0, (44)

(45)

--n K(xl, + n K( у,

x mN* y mN*

де k

^ m¡

прямокутна матриця розмiрностi 2x6.

Утворимо вектор-стовпець значень потенщалу у вузлах m -го СЕ з умовами Неймана

ÜmN* _ (ÜmN1 ,ÜmN2 )*

(46)

i вектор-стовпець значень потенщалу у вах iнших вузлах, що не ввiйшли до складу ÜmN *, (внутрших вузлах i вузлах з умовами Дирихле)

Um*. = (Um*,-,UmL4)., (47)

де * , L4 - вiдповiдно початкове i шнцеве значения ло-кальних номерiв внутрiшнiх вузлiв i вузлiв з умовами Дирихле m -го СЕ.

У виразах (46), (47) потенщали вузлiв розташовуемо в порядку зростання 1'х локальних номерiв. З урахуван-ням (46) i (47) вираз (44) можна представити як

kmL^JmL. + ^„М* = 0, (48)

де: kmL - прямокутна матриця розмiрностi 2x4, утворе-на iз тих стовпщв матрицi kmn, номери яких збiгаються з локальними номерами елеменпв вектора ÜmL, ; kmN -

квадратна матриця розмiрностi 2x2, утворена iз тих стовпцiв матрицi kmn, номери яких збiгаються з локальними номерами елеменпв вектора ÜmN*.

З рiвияния (48) визначаемо

UmN* _ kmNkmLÜmL* _ GmÜmL* ■ Gm _ kmNkmL

(49)

де Gm = kmNkmL (50)

прямокутна матриця розмiрностi 2x4, що зв'язуе зна-чення потенцiалу у вузлах з граничними умовами Неймана з його значеннями в iнших вузлах m -го СЕ.

Виконавши транспонування у виразi (49), отри-маемо такий вираз

UmN = -UmLGm*. (51)

Вектор-стовпець Üm значень потенщалу у вузлах m -го СЕ визначаеться через вектори-стовпщ ÜmL*, ÜmN. у виглядi

Üm* = GmNÜmN* + GmLÜ mL* (52)

або у разi транспонованих векторiв

Üm = ÜmNGmN* + Ü„fi^ , (53)

де: GmN, GmL - прямокутш матрицi розмiрностi вщпо-ввдно 6x2 i 6x4, елементами яких е постшш числа 0 або 1; GmN», GmL* - матрищ, транспоноваиi вiдносно мат-

риць GmN , GmL .

Представляючи функщонал m -го СЕ у виглядi

складно1 функцiï Fm = Fm[Üm ] або Fm = Fm[ÜmN ,ÜmL ] i

ховуючи (51), (53), отримаемо

dFm _ dÜm dFm +

dÜmL dÜmLdÜm dÜnL dÜ„N dÜ„

_ Gm

dFm

m

' dÜ..

dFm dFm dFm Ü

Gm*GmN* „Ü _ (GmL* ~ Gm*GmN*) _ Qm ,f% _ 0.тфт

dÜm dÜm dÜm

де

Qm _ GmL* Gm*GmN *

(54)

(55)

прямокутна матриця розмiрностi 4x6, що забезпечуе ви-конання умов Неймана при переходi вщ вектора-стов-пця фп* до вектора-стовпця фт1* розмiрностi 4.

Аналогiчно, представляючи вектор-стовпець фп* у

виглядi фт* _ ^j*[Üm*] або Фт* _ 4*[ÜmN*,ÜmL*] i вр^ову-ючи (49), (52), (54), отримаемо квадратну матрицю роз-мiрностi 4x4 вигляду

PmL _

d܄

d Щт*)

dÜ mL*

Qm

d܄

._ Q ф dÜm* + Qm Ü Ü

dÜ m* dÜmL*

+Qm

dÜ„ dÜm

____ Q ТгП_ G — Q _

dÜ m* dÜmN* dÜmL* "" dÜm* ^ " dÜm* ^ "

Gntfim _ (56)

= Qm

~(Gml - GmNGm ) _ QmPmQm*

де Qm* _ GmL - GmNGm (57) - матриця, транспонована вщ-носно матрицi Qm.

Матрищ Qm, Qm* у виразi (56) забезпечують вико-нання граничних умов Неймана при переходi вщ матри-цi <рт до матрищ pmL за наявностi в п -му СЕ вузлiв iз граничними умовами Неймана.

Осшльки матриця pm симетрична, то, зпдно з (56), матриця < mL також симетрична.

Внесок кожного п -го СЕ, який мае вузли iз граничними умовами Неймана, у вектор нев'язок ф* i матрицю Якобi p на кожнш iтерацiï визначаеться вiдповiдно з Фт1* i PmL за правилами, встановленими рашше для внутрiшнiх СЕ.

Висновки

1. Сформульовано краеву задачу розрахунку стащ-онарного температурного поля при наявносп внутрш-нiх джерел тепла.

2. Побудовано ск1нченно-елементну модель розрахунку розпод^ температури всерединi двовимiрноï областi, заповненоï нелiнiйними безгiстерезисними аш-зотропними середовищами з урахуванням заданих вздовж ïï граиицi граничних умов Дирихле та однорщ-них граничних умов Неймана.

3. Застосовано кубатурну формулу чисельного штег-рування по площi трикутника на основi iигерполяцiйно-го повного полiному для лагражевого трикутника другого порядку.

4. Отримано з умови мiнiмуму функц1онала рiвно-сильну нелiнiйну систему алгебрачних рiвнянь, яка розв'язуеться, як правило, ггерацшним методом Ньютона.

5. Подано алгоритм визначення внеску кожного внутрiшнього СЕ у вектор нев'язок i матрицю Якобi системи рiвнянь з використанням елементiв тензора ди-ференцiальноï теплопровщносп середовища.

6. Наведено алгоритм визначення внеску кожного граничного СЕ у вектор нев'язок i матрицю Якоб^ якщо вiн мае тiльки один чи дешлька граничних вузлiв з гра-ничними умовами Дирихле або мае два граничних вуз-ла з однорiдними граничними умовами Неймана i один граничний вузол з граничними умовами Дирихле.

Перелш використаних джерел

Dyshovyj, R. V. (1983). Raschet staticheskogo magnitnogo polja v nejavnopoljusnyh jelektiicheskih mashinah difFerencialnym se-tochnym metodom. Abstract of candidate dissertation for technical sciences. Lvov, 18 p. [in Russian]. Karashetskyi, V. P. (2007). Kubaturni formuly chyselnoho intehru-vannia za ploshcheiu trykutnyka na osnovi interpoliatsiinykh povnykh polinomiv. [Cubature formulas for numerical integration through triangle area with full interpolation polynomials]. Scientific Bulletin of ÜNFÜ, 17(7), 275-280. [in Ukrainian]. Silvester, P., Cabayan, H. S., & Browne, B. T. (1973). Efficient techniques for finite element analysis of electric machines. JEEE Trans. PAS, 92(4), 1274-1281. Yushko, S. V., Borshch, O. Ye., & Yushko, M. A. (2011). Statsionar-na teploprovidnist: navch. posibn. Kharkiv: НТУ "ХП1", 80 p. [in Ukrainian].

m

М. И. Дмитрусь, В. П. Карашецкий

Национальный лесотехнический университет Украины, г. Львов, Украина

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

Сформулирована краевая задача расчета стационарного температурного поля при наличии внутренних источников тепла, которая описывается дифференциальными уравнениями. Для построения конечно-элементной модели расчета распределения температуры внутри двумерной области, заполненной нелинейными безгистерезисными анизотропными средами, учета заданных вдоль ее границы граничных условий Дирихле и однородных граничных условий Неймана, использованы лагранжевые треугольники n-го порядка. Применена кубатурная формула численного интегрирования по площади треугольника на основе интерполяционного полного полинома для лагранжевого треугольника второго порядка. Из условия минимума функционала получена равносильная нелинейная система уравнений, решаемая, как правило, итерационным методом Ньютона. Представлен алгоритм определения вклада каждого внутреннего конечного элемента (КЭ) в вектор невязок и матрицу Якоби системы уравнений с использованием элементов тензора дифференциальной теплопроводности среды. Приведен также алгоритм определения вклада каждого предельного КЭ в вектор невязок и матрицу Якоби, если он имеет только один или несколько граничных узлов с граничными условиями Дирихле или имеет два граничных узла с однородными граничными условиями Неймана и один граничный узел с граничными условиями Дирихле.

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

M. I. Dmytrus, V. P. Karashetskyy

Ukrainian National Forestry University, Lviv, Ukraine

CALCULATION OF TWO-DIMENSIONAL STATIONARY THERMAL FIELDS WITH INTERNAL

HEAT SOURCES OF THE FINITE ELEMENT METHOD

The boundary value problem of calculating the stationary temperature field with internal heat sources, which is described by differential equations, is formulated. To construct a finite-element model for calculating the temperature distribution within a two-dimensional region filled with nonlinear anisotropic environments without hysteresis, considering of the Dirichlet boundary conditions and the homogeneous Neumann boundary conditions given along its boundary, were used n-order Lagrangian triangles. The cubature formula of numerical integration over the area of a triangle is used based on the interpolation complete polynomial for the Lagrangian triangle of the second order. From the condition of the minimum of the functional, an equivalent non-linear system of equations is obtained, which is solved, usually, by the iteration method of Newton. An algorithm is presented for determining the contribution of each inner finite element (FE) to the vector of discrepancies and the Jacobi matrix of the system of equations using the elements of the differential thermal conductivity tensor of the environment. An algorithm is also given for determining the contribution of each boundary FE to vector of discrepancies and the Jacobi matrix if it has only one or more boundary nodes with the Dirichlet boundary conditions or has two boundary nodes with the homogeneous Neumann boundary conditions and one boundary node with Dirichlet boundary conditions.

Keywords: stationary temperature field, internal heat sources, thermal conductivity, Lagrangian triangle, finite element method, cubature formulas, Newton's method, differential thermal conductivity tensor, boundary conditions.