Чисельні методи в задачах теплопровідності

Припустимо, що розв'язується краєва задача теплопровідності для рівняння Лапласа в області з границею . Границя складається з двох частин . На ділянці границі відомі значення температури, на ділянці границі відомі теплові потоки.

Рис. 2.1

Якщо на границі області відомі температура і тепловий потік , то температура в будь-якій внутрішній точці визначається по формулі [3]

(2.1)

Температура в будь-якій точці , яка належить границі області , визначається по формулі [3]

(2.2)

В цих формулах - фундаментальний розв'язок рівняння Лапласа, - похідна від фундаментального розв'язку в напрямку зовнішньої нормалі до границі області [4,5].

Коефіцієнт , якщо вузол належить гладкій границі. Для решти випадків значення коефіцієнту визначається в процесі розв'язування задачі [3].

Для двовимірного випадку [4,5]

, (2.3)

тобто

; (2.4)

; (2.5)

Припустимо, що границя розбита на елементів, тобто , тоді граничне інтегральне рівняння (2.2) можна записати у вигляді [3]

. (2.6)

Якщо тепер припустити, що значення функцій і є сталими в межах граничного елемента , то їх можна винести за знак інтеграла і ми одержимо вираз

. (2.7)

Інтеграли задають зв'язок між i-им вузлом і граничним елементом , по якому обчислюється інтеграл. Будемо позначати ці інтеграли , аналогічно інтеграли виду будемо позначати , тоді рівняння (2.7) запишеться у вигляді

(2.8)

Введемо позначення

(2.9)

Тоді рівняння (2.8) можна записати у вигляді

. (2.10)

Припустимо, що на границі S зафіксовано N точок (вузлів) і для кожної точки (вузла) записано інтегральне рівняння (2.8), тоді наступна система співвідношень, записана у матричному вигляді через рівняння (2.10), визначає зв'язок між значеннями функцій U і q на границі області

HU = GQ. (2.11)

Згідно з умовою задачі, на границі S відомі N₁ значень функції U і N₂ значень функції q, тому система співвідношень (2.11) дає можливість сформувати систему лінійних алгебраїчних рівнянь для знаходження невідомих значень функцій U і q на границі області. Для цього систему співвідношень (2.11) треба перетворити в СЛАР виду

, (2.12)

виконавши перенос коефіцієнтів при всіх невідомих в матрицю A. Вектор правої частини СЛАР одержимо, виконавши множення відповідних коефіцієнтів на відомі значення функцій U і q на границі області.

Розв'язавши СЛАР (2.12), ми одержимо всі значення функцій U і q на границі області і матимемо можливість знаходити значення температур і теплових потоків в будь-якій внутрішній точці по формулі

(2.13)

Формула (2.13) задає в інтегральному вигляді зв'язок між внутрішньою точкою області і значеннями функцій U і q на границі області. Формулу (2.13) доцільно записати у вигляді

(2.14)

Значення теплових потоків можна обчислити виконавши диференціювання формули (2.14), а саме

(2.15)

в якій x_k - координати досліджуваної точки, k = 1,2 для двовимірної і k = 1, 2,3 для тривимірної області.

Інтеграли і можна обчислити, використовуючи для всіх граничних елементів (крім елемента, на якому розташований вузол) прості квадратурні формули Гауса

(2.16)

(2.17)

в яких l _j - довжина елемента S _j, w_k - коефіцієнти квадратурної формули.

У випадку, коли вузол x _i розташований в межах граничного елемента S_i, присутність в розрахункових формулах виразу , зумовлює появу логарифмічної особливості підінтегральної функції і для обчислення інтегралів треба використовувати спеціальні квадратурні формули. Якщо значення функцій U і q є сталими в межах граничного елемента S_i, то інтеграли і легко обчислити аналітично. Діагональний елемент матриці ,наприклад, дорівнює нулю тому, що в даному випадку нормаль і поверхня граничного елемента ортогональні

(2.18)

Для діагонального елемента одержимо формулу

(2.19)

Використовуючи місцеву координату , пов'язану з елементом S_i так, як це показано на рис. 2.2, будемо мати

(2.20)

Рис. 2.2

Виконавши заміну змінних , при умові що , одержимо

(2.21)

Через те, що останній інтеграл дорівнює 1, одержимо

(2.22)

Якщо на границі області задано граничні умови 3-го роду, то це означає що в кожному вузлі задано лінійний вираз значень температури і теплового потоку у вигляді

eU + f q = d. (2.23)

Якщо в умові (2.23) e = 0, то матимемо задачу Неймана, якщо f = 0, то матимемо задачу Діріхле. Якщо і , то будемо мати граничну умову конвективного типу. Якщо умову (2.23) записати для всіх вузлів, то одержимо вираз

Q = D - E U, (2.24)

в якому вектор D і діагональна матриця E містять в кожному вузлі величини відповідно d / f і e / f. Підставляючи умову (2.24) в рівняння (2.11), одержимо систему рівнянь

(H + G E) = GD, (2.25)

яку можна подати у вигляді

AY = F. (2.26)

Розв'язавши систему (2.26), з допомогою умов (2.23) можна знайти значення функцій U і q в вузлах на границі області.

2.2 Двовимірні задачі для рівняння Пуассона

Розглянемо випадок, коли функції U і q змінюються на граничному елементі по лінійному закону і вузли розташовані на стику двох прямолінійних граничних елементів.

В будь-якій точці границі області виконується рівняння

(2.33)

яке можна записати в дискретному вигляді так

(2.34)

В формулі (2.33) коефіцієнт , якщо вузол належить гладкій границі. Для решти випадків значення коефіцієнту с_i визначається в процесі розв'язування задачі.

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

Рис. 2.4

(2.35)

Безрозмірна координата дорівнює = x / (l / 2), а функції ₁ і ₂ задаються формулами

, . (2.36)

Тепер інтеграли по довжині граничного елемента S _j в формулі (2.34) можна подати у вигляді наступної формули

(2.37)

в якій .

Для інтегралів, які стоять в правій частині формули (2.34) можна записати наступний вираз

(2.38)

в якому

Для того щоб записати в дискретному вигляді рівняння, яке відповідає i_му вузлу на границі, треба сформувати ті коефіцієнти рівняння, які враховують вклад двох сусідніх граничних елементів S _{j - 1} і S _j в цьому вузлі. В результаті одержимо вираз

, (2.39)

в якому кожен елемент дорівнює сумі елемента , обчисленого на S_j-1 і елемента , обчисленого на S _j, якщо використовувати нумерацію в напрямку проти стрілки годинника. Це справедливо і для обчислення елементів . Тому граничне інтегральне рівняння для i-го вузла в дискретному вигляді можна записати так:

, (2.40)

або простіше, у вигляді виразу

, (2.41)

в якому

Якщо формулу (2.41) записати для всіх вузлів, то одержимо систему NN рівнянь, яка визначає зв'язок між значеннями функцій U і q на границі S області D, і в матричній формі має вигляд

H U = G Q. (2.42)

Якщо в вузлі з номером i границя області не є гладкою, то рівність не виконується, і діагональні елементи матриці H треба обчислювати, використовуючи ту обставину, що якщо на всій границі задано сталу температуру, то теплових потоків через границю не буде. При цьому умова (2.42) для замкнутих областей матиме вигляд

H I = 0. (2.43)

Із рівняння (2.43) слідує, що сума всіх елементів рядка матриці H має дорівнювати нулю і це дає змогу обчислити діагональний елемент, якщо відомі всі позадіагональні елементи, а саме

(2.44)

Діагональні елементи G_{i j} матриці G для лінійних елементів можна обчислити аналітично так

(2.45)

Тут l _i - ₁ і l _i довжини двох граничних елементів, які контактують в i-му вузлі.

Використовуючи задані в задачі граничні умови, систему співвідношень (2.42) можна переписати так, що в результаті одержимо систему СЛАР виду

AY = F, (2.46)

в якій невідомими у векторі Y будуть всі невідомі значення краєвих умов задачі в вузлах на границі області.

2.4 Квадратична апроксимація на елементі

Через те, що на границі області задані теплові потоки, у вираз для точного розв'язку входить довільна стала, значення якої визначається з додаткових умов. У нас це ілюстративна задача, покладемо для зручності Const = 55;

Нехай:

Тоді постановка задачі має вигляд:

Рівняння:

Область:

Краєві умови:

Розв'язок:

Для застосування методу граничних елементів треба знати значення функції і її нормальну похідну на границі області постановки крайової задачі. В даному випадку ми маємо задачу Неймана. Відомі теплові потоки на границі області, треба знайти значення температур на границі області. Сітка методу граничних елементів і значення нормальної похідної q(x,y) в вузлах сітки представлені на рис. 2.8.

Рис. 2.8

В даному випадку:

- границя області розбита на 28 елементів;

- використовується елемент 2-го порядку, тому маємо 56 вузлів;

На рис. 2.8, щоб не затемнювати рисунок, жирним шрифтом проставлено номери вузлів тільки в кутових точках.

При такій нумерації вузлів маємо СЛАР 56 рівнянь, 56 невідомих.

Матриця СЛАР є квадратною повністю заповненою;

- через те, що це задача Неймана, то для визначеності СЛАР будемо вважати відомим значення функції в вузлі №29.

Результати виконаних обчислень представлені на наступних рисунках.

Рис. 2.9 Рис. 2.10

Результат, представлений на Рис. 2.9, одержано при використанні сітки із 28 граничних елементів, 56 вузлів, довжина елемента дорівнює 0,5. На рисунку жирним шрифтом проставлено номери вузлів. Числа напроти вузлів це різниця між точним значення функції в вузлі і значенням , одержаним після розв'язування СЛАР.

Результат, представлений на рис. 2.10, одержано при використанні сітки із 112 граничних елементів, 224 вузлів, довжина елемента дорівнює 0,125. На рисунку жирним шрифтом проставлено номери вузлів, які відповідають вузлам, представленим на рис. 2.9. Аналогічно, числа напроти вузлів це різниця між точним значення функції в вузлі і значенням , одержаним після розв'язування СЛАР. Видно, що точність визначення значень функції, дуже помітно виросла. Видно також, що найгірша точність спостерігається в вузлах, які належать граничним елементам, які контактують з кутовою точкою.

В вузлі з номером 17 (див. Рис. 2.9) має місце стрибкоподібна зміна значення нормальної похідної і цим зумовлено зниження точності безпосередньо в околі кутової точки. Аналогічна картина спостерігається і в кутових точках з номерами 1, 17, 45. В кутовій точці 29 стрибка значення нормальної похідної немає і немає погіршення точності.

2.6 Розрахунок температурного поля в області з ребрами

В сучасних технічних пристроях, з метою інтенсифікації теплопередачі та зменшення габаритів теплообмінної апаратури, широко використовуються ребристі (розвинуті) поверхні. На сьогодні сфера їх використання значно розширилася [12,13,14].

Рис. 2.11 Стінка, оребрена з однієї сторони

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

В інженерній практиці розрахунок ребристих поверхонь виконується по одномірній або двомірній моделі, в залежності від геометричних розмірів стінки та ребра, а також необхідної точності розв'язання задачі. При цьому найчастіше використовують систему рівнянь балансу теплових потоків, що проходять через ребристу стінку [13]. Розглядається стаціонарна задача у лінійній постановці. На практиці виконується розрахунок таких ребер, які використовують в конструкціях (прямокутні, трапецієвидні та ін.)

Математично задача зводиться до побудови розв'язку краєвої задачі для рівняння Лапласа в області з ребром. Саме наявність ребра ускладнює побудову розв'язку методом скінченних різниць (МСР) і методом скінченних елементів (МСЕ).

У випадку тонкого ребра для досить точного розв'язування задачі методом (МСР) треба використовувати дуже дрібну сітку і це дуже збільшує об'єм обчислень.

У випадку, коли треба дослідити вид температурного поля в тонких ребрах, метод скінченних елементів (МСЕ) має ряд дуже істотних переваг. У методі (МСЕ) при дискретизації області можна застосовувати елементи різного розміру. Можна на ребрі використовувати дрібні елементи, а поза ребром використовувати елементи більшого розміру.

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

У випадку ізотропних областей, для побудови температурного поля в області з ребрами, зручно використовувати метод граничних елементів (МГЕ).

Розв'язування краєвої задачі методом (МГЕ) дозволяє зменшити на одиницю розмірність задачі, виконуючи всі обчислення для границі досліджуваної області, а також одержати розв'язок неперервний в області постановки краєвої задачі.

Метод (МГЕ) дає змогу одержувати значення температури і теплового потоку в будь-якій точці досліджуваної області і це, в свою чергу, дозволяє одержувати значення ряду параметрів теплообмінних пристроїв, потрібних при їх конструюванні, наприклад таких як:

§ коефіцієнт теплопередачі для стінки з ребрами;

§ температурне поле в міжреберному каналі;

§ температурний напір при основі ребра;

§ коефіцієнт впливу ребра

і вибрати оптимальні характеристики оребрення [12,13].

Результати, представлені на рисунках рис. 2.12 - 2.17, ілюструють можливості методу граничних елементів (МГЕ).

Рис. 2.12

Рис. 2.13

Рис. 2.14

Рис. 2.15

Рис. 2.16

Рис. 2.17

Завдання для лабораторних робіт

В задачах з номерами 1-10, область постановки краєвої задачі прямокутник. Розміри області і сітка елементів представлені на рисунку. Краєві умови задано на певних елементах.

Стрілки показують напрямок теплового потоку, може бути нагрівання, або охолодження області. Число біля стрілок це інтенсивність теплового потоку. В окремих вузлах може бути задана температура. Задані також коефіцієнти теплопровідності у напрямку координатних осей K_x і K_y.

Ламана лінія зі стрілкою означає конвекцію. Число біля стрілки - коефіцієнт конвекції. При наявності конвекції, вказана температура оточуючого середовища T_cp.

Вважається, що на границі області поза вказаними елементами задано теплоізоляцію.

В задачах з номерами 11 - 30, представлену на рисунку сітку елементів пропонується відобразити на заданий чотирикутник з вершинами ABCD і розв'язати відповідну задачу.

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

Вважається, що нумерація вузлів виконується вздовж короткої сторони прямокутника. Всі задачі були попередньо розв'язані і перевірені при виконанні лабораторних робіт студентами. Для контролю одержаного розв'язку в розділі “Відповіді” наведено значення температур в перших п'яти вузлах дискретної сітки скінченних елементів.

Відповіді

Література

Основна література

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

2. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975 г.

3. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов в технике. - М.: Мир, 1987 г.

4. Тихонов А.Н., Самарский А.А. Уравнения математической физики.- М.: Наука. Гл. ред. физ.- мат. лит., 1977 г.

5. Владимиров В.С. Уравнения математической физики, изд. 2-е, Наука,1971г.

Додаткова література

6. Деклу Ж. Метод конечных элементов. - М.: Мир,1976г.

7. Корнеев В.Г. Схемы метода конечных элементов высоких порядков точности. - Л.: Из-во ЛГУ, 1977г.

8. Галлагер Р. Метод конечных элементов. Основы. - М.: Мир,1984 г.

9. Норри Д., де Фриз Ж. Введение в метод конечных элементов. - М.: Мир,1981 г.

10. Безухов Н.И, Лужин О.В. Применение методов теории упругости и пластичности к решению инженерных задач. - М.: Высшая школа,1974 г.

11. Самарский А.А. Теория разностных схем. Учеб. пособие.- М.: Наука 1977.

12. Теория тепломассообмена / Под ред. А.И. Леонтьева. - М.: Изд-во МГТУ им. Н.Э. Баумана, 1997 г.

13. Ройзен Л.И., Дулькин Н.И. Тепловой расчет оребренных поверхностей. - М.: Энергия, 1977 г.

14. Керн Д., Краус А. Развитые поверхности теплообмена.- М.: Энергия, 1977

15. Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение. М., Мир, 1998.

Предметний покажчик

А