Размещено на http://www.Allbest.ru/

Размещено на http://www.Allbest.ru/

Мiнiстерство освiти і науки України

Національний технічний університет України

Київський полiтехнiчний iнститут

Програма підготовки 6.050103 «Програмна інженерія»

Спеціальність 6.05010301 «Програмне забезпечення систем»

Напрямок підготовки 6.050101 «Комп'ютерні науки»

Спеціальність 6.05010101 «Інформаційні технології проектування»

Навчальний посібник

для студентів вищих навчальних закладів

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

О.К. Молодід

Київ - 2016



Рецензенти:

д.т.н., с.н.с. О.С.Богініч,

д.т.н., проф. Л.М.Девін,

Відповідальний редактор - к.т.н., доц. В.М.Медведєва

Рекомендовано Вченою радою Національного технічного університету України “Київський політехнічний інститут” (протокол №5 від 11 квітня 2016р).

Молодід О.К.

Л84 Чисельні методи в задачах теплопровідності: Навч. посіб. - К.: Знання України, 2015. - 84 с.: іл. - Бібліогр.: с. 80.

Посібник написано на основі курсу лекцій, що читався автором для студентів Київського політехнічного інституту спеціальностей "Програмне забезпечення систем", "Інженерія програмного забезпечення".

Викладено два основні методи побудови розв'язків краєвих задач теплопровідності - метод скінченних елементів (МСЕ) і метод граничних елементів (МГЕ). Наведені приклади та завдання для самостійної роботи студентів.

Призначено для студентів напрямків підготовки бакалаврів ”Комп'ютерні науки”, ”Програмна інженерія”, а також може бути використано студентами інших технічних спеціальностей.



Зміст

Вступ

1. Метод скінченних елементів

1.1 Основна концепція методу скінченних елементів

1.2 Дискретизація області

1.3 Інтерполяційний поліном 1-ї степені

для апроксимації функції U = U (x,y) на простому трикутнику

1.4 Інтерполяційний поліном 2-го степеня

для апроксимації функції U = U (x,y) на трикутнику 2-го порядку

1.5Природна система координат для трикутника

1.6 Побудова дискретного аналогу неперервної функції U = U (x,y)

1.7 Зведення краєвої задачі теплопровідності до варіаційної задачі

1.8 Матрична форма запису функціонала

1.9 Мінімізація функціонала

1.10 Розрахункові формули методу скінченних елементів

1.10.1 Формула для обчислення інтеграла І

1.10.2 Формули для обчислення інтеграла II

1.10.3 Формула для обчислення інтеграла III

1.10.4 Формули для обчислення інтеграла IV

1.10.5 Формули для обчислення інтеграла V

1.11 Формування СЛАР методу скінченних елементів

1.12 Особливості СЛАР методу скінченних елементів і її перетворення

1.13 Приклад розв'язування краєвої задачі теплопровідності методом скінченних елементів

1.14 Приклад обчислення матриць елементів

2. Метод граничних елементів

2.1 Двовимірні задачі для рівняння Лапласа

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

2.3 Лінійна апроксимація на елементі

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

2.5 Приклад розв'язування краєвої задачі теплопровідності методом граничних елементів

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

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

Відповіді

Література

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



Вступ

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

В методі скінченних різниць (МСР) розв'язування краєвих задач зводиться до розв'язування систем лінійних алгебраїчних рівнянь (СЛАР) високих порядків шляхом заміни похідних скінченно-різницевими виразами. Метод МСР є загальновідомим і добре описаним в літературі [11].

Метод скінченних елементів (МСЕ) [1,2,7] сформувався в 60-их роках минулого століття і на сьогодні є ефективним засобом розв'язування різноманітних задач фізики і техніки. Метод виник із структурного аналізу, розвинутого в техніці для статичного розрахунку конструкцій і будівель [2,10]. Звідти ідея “дискретизації“ була успішно перенесена на неперервні системи. Спочатку область застосування методу скінченних елементів обмежувалась розрахунком полів статичної деформації пружних систем: пластин, оболонок, арок, гідротехнічних споруд. Потім метод почали застосовувати в задачах динаміки пружних систем, нестаціонарних задачах теплопровідності, нелінійної дифузії, гідродинаміки [1,2].

Швидкий розвиток методу пов'язаний насамперед з його наглядністю. Переваги методу найкраще проявляються при розв'язуванні краєвих задач в обмежених областях неправильної конфігурації з складними умовами на границях. За допомогою типових елементів різної форми і степеня апроксимації можна досить точно змоделювати границю області [7]. При цьому граничні умови не накладають спеціальних умов на елементи, які використовуються.

Підхід до методу скінченних елементів, як до способу побудови моделі для чисельного аналізу системи, дозволяє розглядати його як одну з гілок діакоптики - загального методу дослідження систем шляхом їх дискретизації [2].

На практиці дослідження теплових процесів найчастіше зводиться до розв'язування краєвих задач для рівнянь Лапласа і Пуассона [4,12]. Згідно з теорією формули Гріна [4,5], гармонічну в області D функцію u можна побудувати, якщо відомі значення функції u і її нормальної похідної на границі S області D. Ці значення не можуть бути довільними. Зв'язок між ними задається через граничне інтегральне рівняння, на використанні якого базується метод граничних елементів (МГЕ) [3].

В методі МГЕ спочатку визначаються невідомі величини на границі області, шляхом дискретизації цієї границі і побудови СЛАР відносно невідомих величин з квадратною повністю заповненою матрицею [3]. На відміну від методу МСЕ, розвязок всередині області в методі МГЕ є повністю неперервним.

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

Посібник написано на основі курсу лекцій, що читався автором для студентів Київського політехнічного інституту спеціальностей “Програмне забезпечення систем”, “Інформаційні технології проектування”.

Метою навчального посібника є поглиблення розуміння студентами лекційного матеріалу та закріплення у них знань по математичному моделюванню теплових процесів і чисельному вирішенню краєвих задач теплопровідності перспективними чисельними методами МСЕ і МГЕ.

Структуру посібника і його зміст визначено з урахуванням того, що його будуть використовувати студенти для самостійної роботи, а також під час курсового та дипломного проектування.



1. Метод скінченних елементів

1.1 Основна концепція методу скінченних елементів

Основна ідея методу скінченних елементів полягає в тому, що будь-яку неперервну величину, таку, як температура, тиск і переміщення, можна апроксимувати дискретною моделлю, побудованою з використанням кусково-неперервних функцій, визначених на скінченній кількості підобластей [2].

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

1. В області визначення неперервної величини фіксується скінченна кількість точок. Ці точки називаються вузловими точками або просто вузлами.

2. Значення неперервної величини в кожній вузловій точці вважається змінною, яку треба знайти.

3. Область визначення неперервної величини ділять на скінченну кількість підобластей, які називають елементами. Ці елементи мають спільні вузлові точки і всі разом апроксимують форму області.

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

На рис. 1.1 проілюстровано побудову дискретної моделі у випадку, коли неперервна величина є функцією двох змінних. В двовимірному випадку найчастіше використовують елементи у формі трикутника. Використання елементів у формі трикутників дозволяє будувати дискретні моделі неперервних величин, які задані в областях складної форми [2,7].

Рис. 1.1

Зрозуміло, що для апроксимації всієї поверхні Z = F(x,y) треба в площині XOY помістити ряд трикутників (елементів) і для кожного елемента виконати побудову апроксимації так, як це представлено на рис. 1.1. В результаті неперервна величина Z = F(x,y) буде апроксимована кусково-неперервною дискретною моделлю.

1.2 Дискретизація області

При розв'язуванні краєвих задач методом скінченних елементів область постановки краєвої задачі ділять на елементи. Поділ області на елементи називають дискретизацією області. В двовимірному випадку найчастіше використовують елементи у формі трикутника. Використання елементів у формі трикутників дозволяє будувати дискретні моделі неперервних величин, які задані в областях складної форми [2,7].

Прикладом області складної форми може бути область, зображена на рис. 1.2. Зрозуміло, що найпростіше діляться на елементи області простої геометричної форми. На рис. 1.3 зображено поділ на елементи прямокутної області.

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

Рис. 1.2

Рис. 1.3

Рис. 1.4

Розглянутий спосіб дискретизації області складної форми можна запрограмувати. Рисунки рис. 1.2, 1.3, 1.4 ілюструють реалізацію цього способу дискретизації.

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

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

1.3 Інтерполяційний поліном 1-го степеня для апроксимації функції U = U (x,y) на простому трикутнику

Нехай задана неперервна функція U = U (x,y).

Фіксуємо трикутник з вершинами i, j, k.

Треба побудувати апроксимацію функції U = U (x,y) на трикутнику з вершинами i, j, k поліномом першого степеня у вигляді:

U^(e) (x,y) = [N^(e) (x,y)][U^(e)]; (3.1)

В розгорнутому вигляді апроксимація (3.1) запишеться

U^e) (x,y) = [N_i^(e) (x,y); N_j^(e) (x,y); N_k^(e) (x,y);] . (3.2)

Значок ^(e) означає, що апроксимація дійсна всередині елемента.

N_i^(e) (x,y), N_j^(e) (x,y), N_k^(e) (x,y) - функції форми елемента, які треба знайти;

U_i, U_j, U_k - значення функції U = U (x,y) у вершинах i, j, k елемента;

Графічна постановка задачі апроксимації представлена на рис 1.5.

Інтерполяційний поліном 1-го степеня запишемо у вигляді:

(x,y) = x + ₃ y. (3.3)



Рис. 1.5

Поставимо умову, щоб значення полінома (x,y) у вузлах елемента i, j, k точно співпали з значеннями функції U (x,y) в цих вузлах. Одержимо систему СЛАР для знаходження _3.

. (3.4)

Якщо ввести в розгляд обернену матрицю, то розв'язок системи (3.4) запишеться у вигляді:

(3.5)

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

(x,y) = [a_iU_i + a_jU_j + a_kU_k] + [b_iU_i + b_jU_j + b_kU_k]x +

+ [c_iU_i + c_jU_j + c_kU_k]y (3.6)

Якшо в (3.6) перегрупувати елементи, то одержимо:

(x,y) = [a_i + b_i x + c_i y] U_i + [a_j + b_j x + c_j y] U_j +

+ [a_k + b_k x + c_k y] U_k; (3.7)

Якщо вираз (3.7) порівняти з виразом (3.2) то можна зробити висновок, що задача розв'язана. При цьому функції форми мають вигляд:

N_i^(e) (x,y) = a_i + b_i x + c_i y,

N_j^(e) (x,y) = a_j + b_j x + c_j y, (3.8)

N_k^(e) (x,y) = a_k + b_k x + c_k y

і кожна функція форми виражається через компоненти стовпчиків матриці, оберненої до матриці СЛАР (3.4).

Запис інтерполяційного полінома у вигляді (3.1)

де [N^(e)] - матриця функцій форми елемента є універсальним для елемента будь-якого типу. Такий запис дозволяє скоротити громіздкі аналітичні викладки при розгляді теорії методу скінченних елементів [1,2,7].

Крім значень функції U = U (x,y) на елементі, при розв'язуванні задач треба знати значення частинних похідних і на елементі.

Для того, щоб одержати апроксимацію частинних похідних функції U = U(x,y) на елементі треба продиференціювати формулу (3.2).

(3.9)

, (3.10)

де [B^(e)] -матриця частинних похідних функцій форми елемента.

Запис виразу для апроксимації частинних похідних функції U = U (x,y) у вигляді

[g^(e)] = [B^(e)][U^(e)] (3.11)

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

При використанні простого трикутника апроксимація градієнта функції U = U(x,y) має особливо простий вигляд. У цьому випадку згідно (3.8)

(3.12)

тобто при апроксимації частинних похідних і на простому трикутнику матриця градієнтів [B^(e)] є константою і має вигляд

[B^(e)] = (3.13)

1.4 Інтерполяційний поліном 2-го степеня для апроксимації функції U = U (x,y) на трикутнику 2-го порядку

Рисунок рис.1.1 показує, що апроксимація функції U = U(x,y) поліномом 1-го степеня є досить грубою. Підвищити точність апроксимації на трикутному елементі можна двома способами. Можна просто зменшити розміри елемента. При цьому збільшується кількість елементів і зростає об'єм обчислювальної роботи.

Можна збільшити степінь інтерполяційного полінома. Якщо виконується апроксимація функції U = U(x,y) поліномом 2-го степеня, то інтерполяційний поліном беруть у вигляді:

(4.1)

Для знаходження 6-ти коефіцієнтів поставимо умову, щоб значення полінома (4.1) точно співпали зі значеннями функції U = U(x,y) в 6-ти вузлах трикутного елемента. Вузлами елемента 2-го порядку будуть вершини трикутника i, j, k і середини його сторін, тобто точки l, m, n. В методі скінченних елементів прийнято здійснювати обхід вузлових точок елемента проти годинникової стрілки в послідовності i, l, j, m, k, n. Графічна постановка задачі побудови інтерполяційного полінома 2-го степеня на трикутному елементі 2-го порядку представлена на рис. 1.6.

Рис. 1.6

Для побудови інтерполяційного полінома 2-го степеня у вигляді

, (4.2)



треба для полінома (4.1) дослівно повторити міркування, наведені в пункті 1.3 при побудові полінома 1-го степеня.

1.5 Природна система координат для трикутника

В математичній літературі використовують також терміни L - кординати, барицентричні координати. В L - координатах розташування точки М на трикутнику визначається трьома числами L_1, L₂, L₃.

Якщо 1, 2, 3 - вершини трикутника, то координата L₁ це відношення

Рис. 1.7 Рис. 1.8.

MM' - відстань від точки М до сторони 1.

h₁ - висота, опущена з вершини 1 на сторону 1.

Координати L₂ і L₃ визначаються аналогічно. Якщо з точкою М пов'язати трикутник, основою якого є сторона 1, як на рис. 1.8, то легко помітити, що

Три координати в двовимірному випадку не можуть бути незалежними.

Зв'язок між координатами L₁, L₂, L₃ дається рівністю

L₁ + L₂ + L₃ = 1;

Якщо в системі координат xOy вершини трикутника 1, 2, 3 мають координати відповідно (x₁, y₁₎, (x₂, y₂), (x₃, y₃), то зв'язок між координатами xOy і L - координатами деякої точки M(x,y) задається формулами:

(5.1)

Якщо систему (5.1) розв'язати відносно L₁, L₂, L₃, то одержимо

(5.2)

S^(e) - площа елемента.

Формули (5.2) є виразами для функцій форми простого трикутника, які перед цим були одержані у вигляді формул (3.8). Таким чином, координатні змінні L₁, L₂, L₃ є функціями форми для простого трикутного елемента.

N_i = L₁; N_j = L₂; N_k = L₃; (5.3)

При використанні методу скінченних елементів доводиться обчислювати інтеграли по площі елемента і вздовж його сторін. Перевагою використання L _ координат є наявність інтегральних формул, які спрощують обчислення таких інтегралів [2].

(5.5)

Z - довжина сторони елемента

А - площа елемента

1.6 Побудова дискретного аналогу неперервної функції U = U (x,y)

В пункті 1.3 показано, як будується апроксимація неперервної функції U = U (x,y) на трикутному елементі 1-го порядку. Для побудови дискретного аналога функції U = U (x,y) у всій області її визначення треба:

- поділити область на елементи;

- перенумерувати елементи і вузли;

- побудувати апроксимації функції U = U (x,y) на кожному елементі, тобто для кожного елемента треба побудувати функції форми, враховуючи при цьому нумерацію елементів і вузлів;

- записати дискретний аналог функції U = U (x,y).

Виконаємо побудову дискретного аналога неперервної функції U = U (x,y) в області, представленій на рис 1.9. Номери елементів зображено в кружках.

Рис. 1.9

На елементі ⁽¹⁾ замінюємо функцію U = U (x,y) поліномом першого степеня, побудованим таким чином, що значення полінома співпадають з значеннями функції в вузлах 1, 2, 4. Тому для елемента ⁽¹⁾ апроксимаційний поліном матиме вигляд:

(6.1)

Аналогічні формули можна записати для решти елементів. Для функції U = U (x,y), визначеної на всій сукупності елементів, дискретний аналог запишеться у вигляді:

(6.2)

Вираз (6.2) можна записати коротше:

^(e) (x,y) = [N₁^(e) N₂^(e) N₃^(e) N₄^(e) N₅^(e) N₆^(e)][U]; (6.3)

e = 1, 2, 3, 4 - номери елементів;

Вираз (6.2) для всієї області можна записати ще коротше:

(x,y) = [N][U]; (6.4)

(x,y) - дискретний аналог неперервної функції U = U (x,y);

[N] = [N(x,y)] -матриця функцій форми елементів.

1.7 Зведення краєвої задачі теплопровідності до варіаційної задачі

Нехай є деяка просторова область D з границею S.

Нехай границя S складається з трьох частин S = S₀ + S₁ +S_2;

Для області D з границею S поставлена наступна краєва задача:

Знайти розв'язок рівняння

(7.1)

який задовольняє наступним граничним умовам.

На ділянці границі S₀ задана температура:

. (7.2)

На частині границі S₁ + S₂ задано граничні умови у вигляді:

(7.3)

[n_x; n_y; n_z] - вектор нормалі в точці границі;

h - коефіцієнт конвекції;

q - густина теплового потоку;

U_cp - температура оточуючого середовища;

В квадратних дужках (7.3) записано звичайний вираз для теплового потоку, який тече через границю області. Цей потік тепла урівноважується:

потоком тепла q, який може бути поданий в тіло;

втратами тепла за рахунок обміну з оточуючим середовищем - конвекції;

Конвекція - процес переносу теплоти при переміщенні об'ємів рідини або газу із зони з однією температурою в зону з іншою, при цьому перенос теплоти нерозривно пов'язаний з переносом самого середовища [12]. Цей тепловий потік пропорціональний різниці температур поверхні тіла і оточуючого середовища. Коефіцієнт теплообміну h визначають експериментально. Його називають ще коефіцієнтом зовнішньої теплопровідності [12].

В варіаційному численні доводиться, що функція U = U(x,y,z), яка є розв'язком поставленої краєвої задачі, надає мінімальне значення наступному функціоналу:

(7.4)

Інтегральний вираз (7.4) кожній функції (P) ставить у відповідність деяке число. Якщо функція (P) співпадає з розв'язком поставленої краєвої задачі U(P), то інтегральний вираз набуває мінімального значення. Таким чином, побудова розв'язку поставленої краєвої задачі теплопровідності зводиться до знаходження функції, яка надає мінімум функціоналу (7.4).

1.8 Матрична форма запису функціонала

Вводимо в розгляд матриці:

(8.1)

(8.2)

Тепер функціонал (7.4) можна записати у вигляді:

(8.3)

Припустимо тепер, що в області D введена сітка елементів і на кожному елементі побудована апроксимація ^(e) функції U(P), яка є розв'язком краєвої задачі. Якщо це виконано, то на кожному елементі функцію U(P) і її частинні похідні можна замінити їх апроксимаціями по формулах:



(8.4)

де [Ф] - вектор вузлових значень функції U(x,y), по яких побудована апроксимація.

Якщо вирази (8.4) підставити у (8.3), то функціонал (7.4) буде виражено через вузлові значення [Ф]. При цьому для складової функціонала (7.4) на елементі одержимо наступний вираз (8.5).

(8.5)

Вираз для функціонала (7.4) одержимо просумувавши (8.5) по всіх елементах у вигляді

(8.6)

де E - кількість елементів.

1.9 Мінімізація функціонала

Значення функціонала (8.6) залежить від вектора вузлових значень [Ф]. Для того щоб знайти мінімальне значення функціонала (8.6) на сукупності вузлових значень [Ф] треба виконати наступні дії:

- компоненти ф_і вектора [Ф] вважають незалежними змінними;

- знаходять вирази для частинних похідних

і прирівнюють їх нулю

(9.1)

Сукупність виразів (9.1) дає систему лінійних алгебраїчних рівнянь для знаходження точки мінімуму функціонала (8.6). Через те, що інтегральний функціонал (8.6) можна подати у вигляді суми складових по елементах, на які розбита область, вирази для частинних похідних в системі (9.1) можна формувати поступово в міру обрахунку кожного елемента, використовуючи наступні формули:

(9.2)

(9.3)

(9.4)

Інтеграли в (9.3) визначають матрицю теплопровідності елемента, а інтеграли в (9.4) - вектор навантаження на елемент. Формули (9.3) і (9.4) одержують при диференціюванні виразу (8.5) на окремому елементі. Остаточну СЛАР методу скінченних елементів одержують після урахування вкладів всіх елементів, на які розбита область у вигляді

(9.5)

Алгоритм формування СЛАР (9.5) описано в п. 1.11.

1.10 Розрахункові формули методу скінченних елементів

Припустимо, що треба розв'язати краєву задачу теплопровідності в обмеженій області D з граничними умовами, зображеними на рис. 1.10.

Рис. 1.10

Вводимо наступні позначення:

1) S₀ - частина границі області на якій задано температуру.

2) S₁ - частина границі області на якій задано тепловий потік.

Тепловий потік задано, якщо задана густина теплового потоку q (кількість тепла, що протікає через одиницю довжини (або площі) границі за одиницю часу [12]).

Якщо потік тепла направлено всередину області, то маємо нагрів.

Якщо потік тепла направлено зовні області, то маємо охолодження.

Якщо густина потоку тепла q = 0 на ділянці границі, то маємо теплоізоляцію.

3) S₂ - частина границі області на якій відбувається обмін теплом з оточуючим середовищем, температура якого відома і дорівнює t_{середовища} = t_ср.

4) По всьому об'єму області D можуть діяти джерело тепла (F(x,y) - інтенсивність джерел тепла [4,12]).

На рис. 1.10. для простоти зображено прямокутну область D. Зрозуміло, що область D може мати будь-яку форму і наші міркування від цього не зміняться.

При розв'язуванні поставленої задачі методом скінченних елементів необхідно виконати наступні дії:

- в області D вводиться сітка, область D ділиться на елементи, в області D вводиться сукупність вузлових точок.

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

СЛАР треба сформувати. Для формування СЛАР необхідно виконати обрахунок кожного окремого елемента області D. Це означає, що для кожного елемента області D треба побудувати:

– матрицю теплопровідності елемента.

– вектор навантаження на елемент.

Відомо, що матриця теплопровідності елемента обчислюється по формулі:

Вектор навантаження на елемент обчислюється по формулі:

Ці формули записано в найзагальнішому вигляді і вони дійсні для кожного елемента будь-якого типу. Вигляд формули залежить від розв'язуваної задачі. Якщо в задачі, наприклад, відсутня конвекція, то коефіцієнт h = 0 і вигляд розрахункової формули спрощується. Алгоритм реалізації цих формул залежить від типу елементів, які використовуються при дискретизації області.

Розглянемо алгоритми обчислення інтегралів, які присутні у вищенаведених формулах у випадку, коли елементом є простий трикутник з трьома вузлами.

1.10.1 Формула для обчислення інтеграла І

Якщо елементом є простий трикутник, то матриця частинних похідних [B^(e)] буде постійною на елементі і тому добуток матриць [B^(e)]^Т[D^(e)][B^(e)] можна винести за знак інтеграла.

Через те, що матриця градієнтів має вигляд

[B^(e)] = ;

обчислення інтеграла I полягає в знаходженні добутку матриць [B^(e)]^Т[D^(e)][B^(e)] і множенні кожної компоненти цього добутку на площу елемента S.

1.10.2 Формули для обчислення інтеграла II

Якщо елемент примикає до частини границі області D, на якій відбувається конвекція, то для побудови матриці теплопровідності треба обчислювати інтеграл II:

В цьому інтегралі [N^(e)] - матриця функцій форми елемента.

[N^(e)] = [N_i^(e)(x,y); N_j^(e)(x,y); N_k^(e)(x,y)];

N_i^(e) (x,y) - функції форми елемента.

Тепер доцільно пригадати, що:

- якщо задано трикутний елемент з вершинами (вузлами) i, j, k;

- в вузлах елемента задані значення функції U_i, U_j,U_k;

- побудовані функції форми елемента N_i^(e), N_j^(e), N_k^(e);

то значення функції U(x,y) в будь-якій точці елемента обчислюється за формулою:

U(x,y) = [N_i^(e)(x,y); N_j^(e)(x,y); N_k^(e)(x,y)] ;

По визначенню функції форми елемента мають наступні властивості:

- для кожного вузла елемента з номером i

N_i^(e) (x,y) = 1 в цьому вузлі;

N_i^(e) (x,y) = 0 вздовж сторони трикутника протилежної вузлу i;

Аналітичні вирази для функцій форми елемента можна записати, скориставшись природною системою координат для трикутного елемента. Природна система координат для трикутного елемента задається формулами:

Якщо ці рівняння розв'язати відносно змінних L₁, L₂, L₃ то одержимо вирази для функцій форми елемента N_i^(e)(x,y), N_j^(e)(x,y), N_k^(e) (x,y), записані в основній системі координат.

Звідси маємо висновок, що в природній системі координат функції форми простого трикутного елемента мають вигляд:

Враховуючи це, можемо написати для інтеграла II наступну формулу:

Функції форми залежать від x і y, тому добутки виду N_i N_j не можна винести за знак інтеграла. Значення інтеграла залежить від того, на якій стороні елемента має місце конвективний теплообмін. Враховуючи всі ці міркування можна записати наступні розрахункові формули для обчислення інтеграла II. При цьому використовується формула (5.5) для обчислення інтегралів вздовж сторони елемента в природній системі координат.

вздовж сторони i j;

- довжина сторони між i і j;

вздовж сторони j k;

- довжина сторони між j і k;

вздовж сторони k i;

- довжина сторони між k і i;

1.10.3 Формула для обчислення інтеграла III

Якщо по площі (об'єму) елемента діють джерела тепла з інтенсивністю F, то їх дія враховується при обчисленні вектора навантаження на елемент в інтегралі III. Інтеграл III обчислюється за формулою:

V - об'єм елемента;

Отже, тепло, згенероване в елементі, розподіляється порівну між його вузлами.

1.10.4 Формули для обчислення інтеграла IV

Інтеграл IV входить у формулу обчислення вектора навантаження на елемент у випадку, коли на границі області присутній тепловий потік з інтенсивністю q. Фізично це означає, що на ділянці границі області спостерігається нагрів, охолодження або теплоізоляція [4,12].

Розглянемо поле температури U. Візьмемо елемент поверхні (dS) з певним чином направленою нормаллю і підрахуємо кількість тепла dQ, що протікає через елемент dS в напрямку за проміжок часу dt. Тепло тече від більш нагрітих частин тіла до менш нагрітих і тим скоріше, чим швидше зменшується температура. Кількість тепла dQ пропорціональна dS, dt і градієнту температур в напрямку . Якщо k коефіцієнт теплопровідності, то справедлива формула [4,12]:



Кількість тепла dQ буде величиною додатньою у випадку, коли , тобто, коли в напрямку температура зменшується. Тому в теплотехніці прийнято вважати що [12]:

- якщо потік тепла тече в тіло, то приймається q < 0;

- якщо потік тепла витікає з тіла, то приймається q > 0;

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

- довжина сторони елемента між вузлами i і j;

Справедливе наступне правило.

Якщо між вузлами i і j елемента має місце тепловий потік, то для побудови вектора навантаження на вузли i і j елемента треба обчислити кількість тепла, що протікає через сторону ij елемента і поділити її пополам. Одержана величина і буде навантаженням в кожному з вузлів.

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

1.10.5 Формули для обчислення інтеграла V

Інтеграл V ідентичний інтегралу IV з урахуванням заміни q на добуток hU_cp.

1.11 Формування СЛАР методу скінченних елементів

Припустимо, що розв'язується задача теплопровідності, умова якої представлена на рис. 1.14. В даній задачі область розбита на 24 елементи і в області маємо 20 вузлів. Для розв'язування задачі треба сформувати СЛАР в якій невідомими будуть значення температур в вузлах. Для формування СЛАР в памяті комп'ютера слід передбачити двовимірний масив з 20-ти рядків і 21 стовпчика. Позначимо цей масиа A[1..20, 1..21] і в самому початку роботи заповнимо його нулями.

При формуванні СЛАР для кожного елемента треба побудувати матрицю теплопровідності і вектор навантаження на елемент згідно з формулами, наведеними в п. 1.10.

Припустимо, що в процесі формування СЛАР ми дійшли до елемента №7 і треба врахувати вклад цього елемента в масив A[1..20, 1..21] (див. рис. 1.11)

Рис. 1.11

Виконавши обрахунок елемента 5, 9, 10, одержимо матрицю теплопровідності цього елемента M = { m_ij }_{i,j = 1}³ і вектор навантаження на елемент f = [f₅ f₉ f₁₀]^T. Матриця теплопровідності і вектор навантаження на елемент пов'язані з нумерацією елемента так, як це представлено в (11.1)

5 9 10 f

5 m₁₁ m₁₂ m₁₃ f₅

9 m₂₁ m₂₂ m₂₃ f₉ (11.1)

10 m₃₁ m₃₂ m₃₃ f₁₀

Вклад елемента 5, 9, 10 в підсумкову СЛАР враховується так:

A[5, 5]: = A[5, 5] + m₁₁;

A[9, 5]: = A[9, 5] + m₂₁; (11.2)

A[10,10]: = A[10,10] + m₃₃;

A[5,21]: = A[5,21] - f₅;

A[9,21]: = A[9,21] - f₉;

A[10,21]: = A[10,21] - f₁₀;

1.12 Особливості СЛАР методу скінченних елементів і її перетворення

Можна довести, що матриця СЛАР, яка виникає при розв'язуванні задач теплопровідності методом скінченних елементів є симетричною і додатньо визначеною [4,12].

При використанні великої кількості елементів матриця СЛАР буде сильно розрідженою, тобто більшість її коефіцієнтів будуть нулями. Очевидно, що коефіцієнт A_m,j в m - му рядку матриці буде відмінним від нуля лише тоді, коли вузли з номерами m і j будуть вершинами якогось спільного для них елемента.

Якщо задача розв'язується для прямокутної області і нумерація вузлів виконується вздовж короткої сторони прямокутника так, як це представлено на рис. 1.3., то легко помітити, що між номерами вузлів, які належать одному і тому ж трикутнику, зберігатиметься певна закономірність. Існування такої закономірності зумовлює ту обставину, що матриця СЛАР буде стрічковою. В стрічковій матриці всі ненульові коєфіцієнти і деякі нульові знаходяться між двома лініями, паралельними головній діагоналі. Відстань між головною діагоналлю і цими лініями називається шириною полоси напівстрічки) матриці.

Рис. 1.12

Всі компоненти матриці за межами цієї полоси дорівнюють нулю і їх не треба зберігати в пам'яті ПК. Грамотно написана програма розв'язування СЛАР з симетричною стрічковою матрицею при побудові розв'язку СЛАР використовує лише числа з верхньої напівстрічки матриці.

Ширина напівстрічки обчислюється за формулою

B = (R + 1) S;

R - максимальна по номерах величина найбільшої різниці між

номерами вузлів в окремому елементі.

S - кількість невідомих в кожному вузлі.

Мінімізація величини B пов'язана з мінімізацією R і виконується при послідовній нумерації вузлів при пересуванні у напрямку найменшого розміру області.

Підсумкова система рівнянь методу скінченних елементів має вигляд



[K] [Ф] = [F];

Цю систему рівнянь треба відкорегувати, якщо деякі компоненти вектора [Ф] відомі.

В більшості задач теплопровідності деякі значення шуканих величин є заданими [12]. Якщо деяка компонента ф_i вектора [Ф] є відомою, то СЛАР треба відкорегувати так, щоб при розв'язуванні СЛАР одержати правильну відповідь, не змінюючи розмірностей матриці [K] і вектора [Ф]. Процедуру корегування СЛАР найпростіше зрозуміти на конкретному прикладі. Якщо, наприклад, відоме значення ф₅ то для перетворення СЛАР треба виконати наступні дії:

1. - всі компоненти 5-го рядка матриці СЛАР, крім діагонального, прирівнюються нулю. Діагональна компонента залишається незмінною;

K_5,j = 0; для j = 1, 2,.... N; j ? 5;

2. - компонента f₅ вектора [F] замінюється на добуток K_5,5ф₅;

3. - решта рівнянь СЛАР, крім рівняння 5, перетворюються відніманням добутка K_i,5ф₅ із f_i і підстановкою

K_i,5 = 0; для i = 1, 2,.... N; i ? 5.

1.13 Приклад розв'язування краєвої задачі теплопровідності методом скінченних елементів

Задача

Є циліндр з прямокутним перерізом

Через грань в циліндр подається кількість тепла Q, а через грань із циліндра витікає кількість тепла Q. На решті граней задана теплоізоляція. Знайти температури всередині циліндра, вважаючи, що теплові потоки рівномірно розподілені по поверхні граней.

Графічно умова задачі зображена на рис. 1.13.

Рис. 1.13

Математичний запис умови краєвої задачі:

Рівняння:

Область:

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

Відомий точний розв'язок:

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

Нехай:

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

Рівняння:

Область:

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

Розв'язок:

В термінах методу скінченних елементів умова краєвої задачі представлена на рис. 1. 14.

Рис. 1.14



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

- область розбита на 24 елементи;

- в області маємо 20 вузлів.

Нумерацію вузлів здійснюємо вздовж короткої сторони прямокутника області. На рис. 1.14 проставлено номери лише для частини вузлів щоб не затемнювати рисунок.

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

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

- стрілки показують напрямок теплового потоку на стороні елемента;

- число біля стрілок - це лінійна густина теплового потоку на стороні елемента, яка в розрахункових формулах позначена символом q;

- задана температура у вузлі №16;

Треба сформувати матрицю СЛАР, здійснивши процес ансамблювання системи, так як це описано в п. 1.11. Треба сформувати вектор правої частини СЛАР, скориставшись формулами наведеними в п. 1.10. Через те, що в даній задачі на ділянках границі області задано теплові потоки, то вектор правої частини СЛАР, матиме вигляд, представлений на рис. 1.15.

Числа на рис. 1.15 це компоненти вектора правої частини СЛАР, вирахувані згідно з правилом, наведеним в п.1.10.4. Після виконання ансамблювання, одержимо СЛАР, яка відповідає даній задачі. Цю СЛАР треба відкорегувати, врахувавши ту обставину, що нам відома температура в вузлі №16. Для корегування СЛАР треба скористатися алгоритмом наведеним, в п.1.12.

Рис. 1.15



Верхня напівстрічка СЛАР, яка дає розв'язок поставленої задачі, дається в таблиці 1. Перший стовпчик в таблиці - це діагональ СЛАР. Сьомий стовпчик в таблиці - це вектор правої частини СЛАР.

Таблиця 1

1.	2.50	-1.25	0.00	0.00	-1.25	0.00	-2.50
2.	5.00	-1.25	0.00	0.00	-2.50	0.00	-20.00
3.	5.00	-1.25	0.00	0.00	-2.50	0.00	-20.00
4.	2.50	0.00	0.00	0.00	-1.25	0.00	-10.00
5.	5.00	-2.50	0.00	0.00	-1.25	0.00	15.00
6.	10.00	-2.50	0.00	0.00	-2.50	0.00	0.00
7.	10.00	-2.50	0.00	0.00	-2.50	0.00	0.00
8.	5.00	0.00	0.00	0.00	-1.25	0.00	0.00
9.	5.00	-2.50	0.00	0.00	-1.25	0.00	15.00
10.	10.00	-2.50	0.00	0.00	-2.50	0.00	0.00
11.	10.00	-2.50	0.00	0.00	-2.50	0.00	0.00
12.	5.00	0.00	0.00	0.00	0.00	0.00	67.50
13.	5.00	-2.50	0.00	0.00	-1.25	0.00	15.00
14.	10.00	-2.50	0.00	0.00	-2.50	0.00	0.00
15.	10.00	0.00	0.00	0.00	-2.50	0.00	135.00
16.	5.00	0.00	0.00	0.00	0.00	0.00	270.00
17.	2.50	-1.25	0.00	0.00	0.00	0.00	7.50
18.	5.00	-1.25	0.00	0.00	0.00	0.00	0.00
19.	5.00	-1.25	0.00	0.00	0.00	0.00	0.00
20.	2.50	0.00	0.00	0.00	0.00	0.00	67.50

Розв'язавши цю СЛАР, ми одержимо поле температур в області.

Поле температур представлено на рис. 1.16.

Рис. 1.16

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

1.14 Приклад обчислення матриць елементів

Припустимо, що методом скінченних елементів розв'язується краєва задача теплопровідності, умова якої представлена на рис. 1.17.

Рис. 1.17

В задачі краєві умови задані для елементів з номерами 4, 7, 18, 21. На ділянках границі області поза цими елементами задано теплоізоляцію.

Для елемента 4 на стороні k - i задано конвективний теплообмін з оточуючим середовищем з коефіцієнтом теплообміну h = 24.0;

Для елемента 7 на стороні i - j задано конвективний теплообмін з оточуючим середовищем з коефіцієнтом теплообміну h = 32.0. Конвективний теплообмін відбувається з оточуючим середовищем, температура якого Tcp = 67.4⁰C;

На стороні елемента 18 між вузлами 12 і 16 задано нагрів, тобто на сторону елемента 18 подається тепловий потік лінійна густина якого q = 64;

На стороні елемента 21 між вузлами 18 і 19 задано охолодження, тобто через сторону елемента 21 між вузлами 18 і 19 з області витікає тепловий потік лінійна густина якого q = 16;

Елемент N 4. Конвекція h^(e) = 24.00; Вузли

i 2. (0.0000; 1.0000) Коеф. теплопр.

j 7. (1.0000; 2.0000) Kxx = 2.80;

k 3. (0.0000; 2.0000) Kyy = 1.40;

Записуємо СЛАР (3.4) для побудови апроксимації на елементі 4.

Одержимо розв'язок СЛАР (3.4) у вигляді (3.5)

Записуємо матрицю частинних похідних [B^(4)] у вигляді (3.13)

[B^(4)] =

Довжина сторони (k - i) на якій задана конвекція = 1.00; Tcp = 67.40;

Матрицю теплопровідності

Матрицю конвекції елемента одержуємо,обчисливши елемента одержуємо,обчисливши інтеграл I. інтеграл II.

2. 0.7000 0.0000 -0.7000 8.0000 0.0000 4.0000

7. 0.0000 1.4000 -1.4000 0.0000 0.0000 0.0000

3. -0.7000 -1.4000 2.1000 4.0000 0.0000 8.0000

Матриця теплопровідності Вектор навантаження на вузли елемента 2,7,3. 3 і 2 одержимо, обчисливши інтеграл V, у вигляді

2. 8.7000 0.0000 3.3000 f₂ = - 808.8;

7. 0.0000 1.4000 -1.4000 f₇ = 0.0;

3. 3.3000 -1.4000 10.1000 f₃ = - 808.8;

Елемент N 7. Конвекція h^(e) = 32.00; Вузли:

i 5. (1.0000; 0.0000) Коеф. теплопр.

j 9. (2.0000; 0.0000) Kxx = 2.80;

k 10. (2.0000; 1.0000) Kyy = 1.40;

Записуємо СЛАР (3.4) для побудови апроксимації на елементі 7.

Одержимо розв'язок СЛАР (3.4) у вигляді (3.5)

Записуємо матрицю частинних похідних [B^(7)] у вигляді (3.13)

[B^(7)] =

Довжина сторони (i - j) на якій задана конвекція = 1.00; Tcp = 67.40;

Матрицю теплопровідності Матрицю конвекції елемента одержуємо обчисливши елемента одержуємо,обчисливши інтеграл I. інтеграл II.

5. 1.4000 -1.4000 0.0000 10.6667 5.3333 0.0000

9. -1.4000 2.1000 -0.7000 5.3333 10.6667 0.0000

10. 0.0000 -0.7000 0.7000 0.0000 0.0000 0.0000

Вектор навантаження на вузли. Матриця теплопровідності 5 і 9 одержимо, обчисливши елемента 5,9,10. інтеграл V, у вигляді

5. 12.0667 3.9333 0.0000 f₅ = - 1078.4;

9. 3.9333 12.7667 -0.7000 f₉ = - 1078.4;

10. 0.0000 -0.7000 0.7000 f₁₀ = 0.0;

Елемент N 18. Нагрів q = - 64.00; Вузли:

i 11. (2.0000; 2.0000) Коеф. теплопр.

j 16. (3.0000; 3.0000) Kxx = 2.80;

k 12. (2.0000; 3.0000) Kyy = 1.40;

Записуємо СЛАР (3.4) для побудови апроксимації на елементі 18.

Одержимо розв'язок СЛАР (3.4) у вигляді (3.5)



Записуємо матрицю частинних похідних [B^(18)] у вигляді (3.13)

[B^(18)] =

Довжина сторони (j - k) на якій задано нагрівання = 1.00; q = - 64.0;

Матрицю теплопровідності Інтеграл II на елементі 18 елемента одержуємо обчисливши дорівнює нулю через те, що інтеграл I. конвекції на елементі немає.

Вектор навантаження на вузли

Матриця теплопровідності 16 і 12 одержимо, обчисливши елемента 11,16,12. інтеграл IV.

11. 0.7000 0.0000 -0.7000 f₁₁ = 0.0;

16. 0.0000 1.4000 -1.4000 f₁₆ = - 32.0;

12. -0.7000 -1.4000 2.1000 f₁₂ = - 32.0;

Елемент N 21. Охолодження q = + 16.00; Вузли:

i 14. (3.0000; 1.0000) Коеф. теплопр.

j 18. (4.0000; 1.0000) Kxx = 2.80;

k 19. (4.0000; 2.0000) Kyy = 1.40;

Записуємо СЛАР (3.4) для побудови апроксимації на елементі 21.

Одержимо розв'язок СЛАР (3.4) у вигляді (3.5)



Записуємо матрицю частинних похідних [B^(21)] у вигляді (3.13)

[B^(21)] =

Довжина сторони (j - k) на якій задано нагрівання = 1.00; q = + 16.0;

Матрицю теплопровідності Інтеграл II на елементі 21 елемента одержуємо обчисливши дорі...