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

УДК 536.21

РОЗРАХУНОК ДВОВИМ1РНИХ СТАТИЧНИХ БЕЗВИХРОВИХ ТЕПЛОВИХ ПОЛ1В МЕТОДОМ СК1НЧЕННИХ ЕЛЕМЕНТ1В Н.1. Дмитрусь1, В.П. Карашецький2

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

Ключовi слова: безвихрове теплове поле, теплопровщшсть, лагранжевий трикутник, метод скшченних елементiв, кубатурна формула, граничш умови.

Для краево! задачi розрахунку статичного безвихрового теплового поля, що описуеться рiвняннями:

ЖУЩ ) = 0, (1)

Н = -ътаи , (2)

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

F = | WdS , (3)

я

Б _ _

де: Ш = |ШБ; (4)

о

и - температура; Н, Б - розмщеш у площинi В вектори напруженосп теплового поля та густини теплового потоку; Я - площа областi В.

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

^ = 0. (5)

dU

Для побудови синченно-елементно!' моделi заповнимо область розрахунку В сукупнктю лагранжевих трикутникiв п -го порядку [2].

Нехай внаслвдок трiангуляцií двовимiрноí обласп розрахунку В отри-муемо М лагранжевих скiнченних елементш (СЕ). Кожному з них присво'мо порядковий номер т (т = 1,М) i локальну нумерацда вузлiв, згiдно з якою г -му вузлу т -го СЕ вiдповiдае номер тг. Для вск!' областi розрахунку встановимо откову (наскрiзну) нумерацiю К внутрiшнiх вузлiв i О граничних вузлш. По-точнi значения порядкових номерiв внутрiшнiх вузлiв позначимо г.

Для скiиченно-елементноí областi функщонал F з урахуванням (3), (4) на-буде вигляду:

1 аспiр. Н.1. Дмитрусь - НЛТУ Украши, м. Львгв;

2 доц. В.П. Карашецький, канд. техн. наук - НЛТУ Украши, м. Льв1в

де

м м

F = ^ ^т = ^ Шт , т=1 т=1

Шт = | WdS;

(6) (7)

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

1 У1

1

т 2

(8)

1 *2 У2 ■ 1 уз

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

и = (иь...,ий); и = (иь...,иЙ>. (9)

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

- - dF

Фи*] = ^г = 0.

dU*

(10)

Застосуемо для (7) кубатурш формули чисельного жегрування за площею лагранжевого трикутника [2] ■ Наприклад, у випадку використання лагранжевих трикутникiв другого порядку (п =2), кшьккть вузлш у яких р=6, отримаемо

1 6

F = — $ У Ш

1 т ^т>у т 3 >=4

(11)

де

Шт> = | HdB ■

(12)

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

и = иткт1* = кк^Лт* ,

де ит = (итЪ---,итб)

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

к = (1, х, у, ху, х2, у2) координатний вектор-рядок поточно! точки з координатами х, у;

хт1 ут1 хт1ут1 хт1

ут1

хт2 ут2 хт2Ут2 хт2 ут2

2 2

хт3 ут3 хт3ут3 хт3 ут3

2 2

хт4 ут4 хт4ут4 хт4 хт4

2 2

хт5 ут5 хт5ут5 хт5 ут5 22

хт6 ут6 хтбУт6 хт6 ут6

(13)

(14)

(15)

(16)

я

0

кт* =

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

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

H = -gradU = -^7-^7. (17)

ax ay

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

H ■ = -au I x y ■ = -U K(x) = -K(xU »■ (18)

11 xmi ~ У mi ~ ^ mmi* ~ Jvmiwm*; V10/

ax

H ymi = ^ |xmb ymi = -UmK)(ii* = —KmíUm*, (19)

ay

де: Kmx = kmx)km1; Km) = km*; (20)

KKm* = km1^^*; K m* = kmky*; (21)

kmx = \xmbymi = (0,1,0,ymi,2xmi,0) ; (22)

ax

kimi^ = \xmU ymi = (0,0,1, xmU 0,2ymi) ; (23)

ay

kmx*, km* - стовпщ, отриманi транспонуванням рядкiв (22), (23).

Диференцiюючи вираз (12) по вектору Um i враховуючи (18), (19), отри-маемо

^ dF 1 6 d Hml _ — 1 6 dH _

Фп* = —Um = 3 Sm£ —U— | (dH )B = - Sm£ ~UjmLBmi =

dU m 3 i=4dU m 0 3 i=4 dU m

1 6 dH dH 16 ^ ^

= ~Sm^ (~jj=x Bxmi + 'Bymi)--Sm^, (Kmi*Bxmi + КППу7*ВУп7) , (24)

3 i = 4 dUm dUm 3 i =4

де Bmi, Bxmi, Bymi - вiдповiдно вектор густини теплового потоку i його складовi в mi -му вузлi, яю визначаються за значениями проекцiй (18), (19) вектора нап-руженостi i характеристикою теплопровiдностi нелiнiйного безгiстерезисного середовища, що виражаеться векторним рiвняииям або трьома скалярними рiв-няннями:

В = B[H ], (25)

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

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

Для визначення внеску п -го СЕ в систему рiвнянь (10) потрiбно: • знайти вектор fm* на кожнш iтерацií за формулою (24);

• за таблицею вщповщноси локально!' i cítkobo'í нумерацп встановити номери r вузлiв, якi збiгаються з вузлами mi,...,mp ;

• кожний елемент вектора fm*, який вщповщае r -му внутрiшньому вузлу, внести вщповщно в r -те рiвняння системи (10).

Повну систему ршнянь (10) отримаемо, виконавши цю процедуру для bcíx M елементш. У цьому випадку викладену вище процедуру використовують на еташ формування вектора нев'язок. Бшьш трудомкткою операщею е складання для векторно!' функцй' f = (fi,...,fR)* матриц Якоб1 р розм1рност1 RхR. Виведе-мо загальш вирази, що використовуються для ще!' мети.

Диференцдаючи вираз (24) за вектором Um*, отримуемо матрицю розм1р-носп 6х6

Pm = dfr = -\SmhKm>'§r + Km±dfT). (27)

dU m* 3 г=4 dU m* dUm*

З урахуванням (18), (19), (26) i (27) маемо

( =— 1 s V(K(x)(l dH xmi + 1 dHymi) + K(y)(l dH xmi + 1 dH ymiy. = Pm — (Kmi*(/lxxmi ,f} ~l~/lxymi - )_l~ Kmi*(/lyxmi ,f} + /lyymi . fT ))~

3 i=4 dUm* dUm* dUm* dUm*

6

Jm/

\dBx dBx

dB = dHx dHy

dH " dBy dBy

[H dHy

1 6 _ _ ^^ _ _

= ~ (Кт^-^ххт^пп + 1xymiK<nti) + Kmi^Ay&ni^mi + ^уутшК-тУ), (28) 3 i=4

де 1jkmi(j, k = x, y) - елементи тензора диференщально! теплопровiдностi середо-вища

' дБх дБх '

Я*. Key 1 , (29)

*yx *yy

Гх дИу

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

Для безгктерезисного середовища на ochobí теореми B3aeMHOCTÍ [3]

*су = *ух , тому *ymi = *yxmi.

У випадку iзотропного нелiнiйного середовища тензор диференщально!' теплопровщносп визначають за формулою [1]

\лр cos2 hx +1 sin2 hx (1p - 1) cos hx cos hy ] [(1p-1r)cOShy cos hx 1pCOs2hy +1 sin2 hy

dB B_

~dH' 1= И

лопроввднкть середовища; h( = x, y) - кути мiж вектором И або B[H] i вщпо-вiдно ортами I, j локально! декартово! системи координат.

Для тшйного iзотропного середовища 1p = 1K = B / И = 1s, тому тензор диференщально! теплопроввдносп набувае вигляду

_ЧП_-

1= J 'lx 1 '¡x Vvp 'KJ^w'ix^^'ly I (30)

s2hy +1 sin2 h

де: 1p= , 1 = B - вщповщно рaдiaльнa диференцiaльнa i тaнгенцiaльнa теп-

Для того, щоб визначити вклад т -го СЕ в матрицю Якобi р, потрiбно об-числити матрицю рт на кожнiй iтерацií за формулою (28) i пiдсумувати во 11 елементи з вiдповiдними елементами матриц р, враховуючи, що елемент ртр належить пя -й клггиш матрицi р, де п, я - сiтковi номери вузлiв з локальними номерами т1 i тр.

Повну матрицю Якобi р отримаемо, виконавши цю процедуру для кожного з М скiнченних елементiв областi розрахунку Б. У випадку використання лагранжевих трикутниюв 1-го, 3-го i 4-го порядкiв потрiбно застосувати ввдпо-вiднi кубатурнi формули чисельного iнтегрування [2] i провести вивщ основних залежностей за викладеною вище методикою.

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

Нп=-ди = 0,

дп

(32)

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

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

др

■ т = 0; (33)

дит

ртБр =

дит

= о, р = 1,6;

(34)

ртрБ =

_ ^ Утр _

= о, р = 1,6.

(35)

дитБ

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

Ры = п . (36)

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

Нп = Нп = Нп + Нупу = 0, (37)

де пх, пу - проекцп одиничного вектора п дотично! до границi областi Б.

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

(пхК$* + пуК^-рт* = 0,

де:

К(л?* =

К(х)

л тЫ-1

К(х)

К (у) =

КтЫ* -

К (у) К (у)

(38)

(39)

прямокутш матриц розмiрностi 2x6, рядками яких е вектори-рядки, визначенi за виразами (20) для вузлш з граничними умовами Неймана т -го СЕ; М1, М2 -вiдповiдно початкове i кiнцеве значения локальних номерiв вузлш з граничними умовами Неймана т -го СЕ. Подамо вираз (38) у виглядi

ктпОт* = 0, (40)

де ктп = Пх^тМ* + ПУКгтМ*

(41)

прямокутна матриця розмiрностi 2x6. Утворимо вектор-стовпець значень потенциалу у вузлах т -го СЕ з умовами Неймана

итМ* = (0 тМ!0 mN2)* (42)

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

0тЬ* = (ити,..,итЬ4>, (43)

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

У виразах (42), (43) потенщали вузлш розташовуемо в порядку зростання !х локальних номерiв. З урахуванням (42) i (43) вираз (40) можна подати як

кт1°тЬ* + кт№ тМ* = ^ (44)

де: ктЬ - прямокутна матриця розмiрностi 2x4, утворена iз тих стовпцiв матри-цi ктп, номери яких збiгаються з локальними номерами елемеитiв вектора 0тЬ*; kmN - квадратна матриця розмiрностi 2x2, утворена iз тих стовпщв матрицi ктп, номери яких збкаються з локальними номерами елеменлв вектора 0тМ*. З рiвняння (44) визначаемо

0тМ* =-к:тккт1ртЬ* = -0m0mL*, (45)

де От = ктмкть (46)

прямокутна матриця розмiрностi 2 х 4, що зв'язуе значення потенцiалу у вузлах з граничними умовами Неймана з його значеннями в iнших вузлах т -го СЕ. Виконавши транспонування у виразi (45), одержимо

ЛтМ =-ЛтЬОт* . (47)

Вектор-стовпець 0т значень потенцiалу у вузлах т -го СЕ визначаеться через вектори-стовпщ 0тЬ*, 0тN* у виглядi

°т* = От№ тМ* + ОтЬ° тЬ* (48)

або у випадку транспонованих векторш

Лт = °тМОтМ* + 0mL0mL*, (49)

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

Представляючи функцюнал m -го СЕ у виглядi складно! функцл Fm = Fm[Um] або Fm = Fm[ÜmN,Umb] i враховуючи (47), (49), отримаемо

2 dFm d Üm dFm dÜmN dÜm dFm ^ dFm ^ ^ dFm - + —T2-^T2--= GmL*—T7%--Gm*GmN*-

m ы1 m , ы ^ mN vi-j m m ^ ^ m r^ r^ m

fL* =-2-= —2--2--I--2--2--2— = GmL*-2—

dÜmL ^ÜmLdÜm dÜmL dÜ-N düm dÜ,

dFm dFm -= (GmL* — Gm*GmN*) j = Qm j = Qmfm*, (50)

dÜm dÜm

де Qm = GmL* — Gm<GmN* (51)

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

Аналогiчно, представляючи вектор-стовпець fm* у виглядi fm* = fm*[Üm*] або fm* = fm*[ÜmN*,ÜmL*] i враховуючи (45), (48), (50), одержимо квадратну мат-рицю розмiрностi 4x4 вигляду

j= dfmL* = d(Qmfm*) = q dfm* = q d fm* dÜ m* + q d fm* dÜ m* dÜmN* = dÜmL* dÜmL* dÜmL* dÜm* dÜmL* dÜm* dÜmN* dÜmL*

= Qm -, GmL — Qm , GmNGm = Qm .fj (GmL — GmNGm) = QmpmQm* , (52) dÜm* dÜm* dÜm*

де Qm* = GmL — GmNGm (53)

матриця, транспонована вiдносно матрицi Qm.

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

Внесок кожного m -го СЕ, який мае вузли з граничними умовами Неймана, у вектор нев'язок f* i матрицю Якобi р на кожнш ггерацп визначаеться вщпо-вiдно з fiL* i PirL за правилами, встановленими рашше для внутрiшнiх СЕ.

Лггература

1. Дышовый Р.В. Расчет статического магнитного поля в неявнополюсных электрических машинах дифференциальным сеточным методом : автореф. дисс. на соискание учен. степени канд. техн. наук / Р.В. Дышовый - Львов, 1983. - 18 с.

2. Карашецький В.П. Кубатурш формули чисельного штегрування за площею трикутника на осжга штерполяцшних повних полшомгв // Науковий вiсник НЛТУ Украши : зб. наук.-техн. праць. - Льв1в : РВВ НЛТУ Украши. - 2007. - Вип. 17.7. - С. 275-280.

3. Silvester P. Efficient techniques for finite element analysis of electric machines / P. Silvester, H.S. Cabayan, B.T. Browne // JEEE Trans. PAS. - 1973. - Vol. 92, № 4. - Pp. 1274-1281.

Надшшла до редакцп 30.06.2016р.

Дмитрусь М.И., Карашецкий В.П. Расчет двумерных статических безвихревых тепловых полей методом конечных элементов

Выведены основные формулы метода конечных элементов для краевой задачи расчета двумерных статических безвихревых тепловых полей в областях, заполненных нелинейными безгистерезисными анизотропными средами с использованием лагранжевых треугольников 1-4 порядков, кубатурных формул численного интегрирования и с уче-

том граничных условий Неймана и Дирихле. Рассмотрен алгоритм определения вклада каждого конечного элемента в вектор невязок и матрицу Якоби нелинейной системы уравнений, которая решается методом Ньютона.

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

Dmytrus M.I., Karashetskyy V.P. Calculation of Two-dimensional Static Non-eddy Thermal Fields of the Finite Element Method

The results of our research provide basic formulas of the finite element method for boundary value problem calculating of two-dimensional static non-eddy thermal fields filled with nonlinear without hysteresis anisotropic environments using Lagrangian triangles 14 order, cubature formulas of numerical integration and with considering of boundary conditions of Neumann. Dirichlet was obtained. The algorithm determining the contribution of each finite element in the vector of residuals and matrix Jacobi nonlinear system of equations solved by Newton's method was considered.

Keywords: non-eddy thermal field, thermal conductivity, Lagrangian triangle, finite element method, cubature formula, boundary conditions, Newton's method.

УДК 621.325:669.539.43

ТРИВИМ1РНЕ ПРЕДСТАВЛЕННЯ ЗОБРАЖЕНЬ З ВИКОРИСТАННЯМ IX ФРАКТАЛЬНИХ РОЗМ1РНОСТЕЙ

1.М. Журавель1, В.М. Максимович2

Оптичш мжроскопи вже давно набули застосування у рiзних галузях - вщ медицини i бюлогн до неруйшвного контролю на виробництв^ Доступшсть цифрових вщеокамер призвела до появи нового класу задач, пов'язаних з оптичними мшроскопами. Насампе-ред потрiбно видшити завдання автоматичного аналiзу мшрозображень. Розвиток комп'ютерно! техшки призвгв до можливост виршення одного з найважливших зав-дань автоматичного аналiзу - побудови тривишрних моделей об'екпв за !х зображен-ням. Запропоновано пiдхiд до побудови тривимiрних зображень шляхом використання фрактальних розмiрностей. Основною перевагою запропонованого методу е те, що для побудови рельефу поверхш використовують тiльки одне зображення.

Ключовi слова: оброблення зображень, 3Э моделювання, фрактальна розмiрнiсть, морфометрична карта.

Вступ. Рiзномаштш дослiдження в галузях матерiалознавства, медицини, архггектури, побудова систем "вiртуальноí реальности, наприклад тренажерiв транспортних засобiв, потребують об'емного представления тривишрних об'eктiв. Звичайно, що таке тривимiрне представлення е уявним, оскiльки реаль зуеться, здебшьшого, на двовимiрнiй площинi дисплея мошгора. Але, незважа-ючи на це, воно е бшьш зручним для сприйняття людиною та дае змогу виявити такi особливостi об'екта, яю невидимi при його двовимiрному представленш. Зважаючи на зазначене вище, можна стверджувати, що задача тривимiрного представлення об'екпв е актуальною.

Аналiз вщомих пiдходiв та постановка задачi. Шдход1в до вирiшения цiеí задачi е багато. Серед них видшимо найвiдомiшi - метод побудови триви-

1 ст. наук. с^роб. 1.М. Журавель, канд. техн. наук - НУ "Львгвська полггехнка";

2 проф. В.М. Максимович, д-р техн. наук - НУ "Львгвська полггехнка"