МІНІСТЕРСТВО ВИЩОЇ ОСВІТИ РОСІЇ
АРХАНГЕЛЬСЬКИЙ лісотехнічний інститут
К а ф е д р а т е п л о т е х н и к і
РОЗРОБКА ПРОГРАМИ
ДЛЯ РІШЕННЯ неодновимірному СТАЦІОНАРНИХ ЗАВДАНЬ
Теплопровідність Чисельні методи С
ВИКОРИСТАННЯМ Консервативно-Різницеві схеми
А Р Х А Н Г Е Л Ь С К
1 9 9 3
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
О Г Л А В Л Е Н Н Я
Вступ ................................ ... ... ... ... ... ... ... ... ... ... ... ... ... .......
1.Основні положення методики побудови консервативно-різницевої схеми при вирішенні неодновимірному задач стаціонарної теплопровідності .......... . ... ... ... ... ... ... ... ...........
2. Методика підготовки і рішення задачі на ЕОМ .... ... ... ... ... ...
2.1. Постановка задачі, розробка математичної моделі ................................... ... ... ... ... ... ... ... ... ... ... ... ... ... .....
2.2. Вибір методу чисельного рішення ....... ... ... ... ... ... ... ... ......
2.3. Розробка алгоритму та структури ......... ... ... ... ... ... ... ... ......
2.4. Написання програми і підготовка її до введення в ЕОМ ..................... ... ... ... ... ... ... ... ... ... ... ... ... ........ .......
2.5. Тестування, налагодження програми і рішення на ЕОМ
Література ....................... ... ... ... ... ... ... ... ... ... ... ... ... ... ................
В В Е Д Е Н Н Я
Базовий рівень підготовки інженера-енергетика в галузі інформатикита обчислювальної техніки визначається необхідним набором знань, уміньі навичок використання ЕОМ для вирішення різних технічних завдань.
Фахівці цієї категорії, крім уміння використовувати прикладнепрограмне забезпечення, повинні бути програмується користувачем, тому щоїхня професійна діяльність пов'язана з виконанням великої кількостітеплотехнічних розрахунків.
Для дотримання принципу фундаментальності вищої освіти роботапобудована на базі розгляду питань застосування ЕОМ для вирішенняосновних задач теорії теплообміну. До однієї з таких завдань ставиться завдання,пов'язана з визначенням температурного поля не одновимірних тел чисельними методами.
Розглянемо методику підготовки та рішення зазначеної задачі наперсональному комп'ютері.
1. О С Н О В Н И Е П О Л О Ж Е Н Н Я М Е Т О Д И К И
П О С Т Р О Е Н Н Я К О Н С Е Р В А Т И В Н О-Р А З Н О С Т Н О Й С Х Е
М И ПРИ Р І Ш Е Н І І Н Е О Д Н О М Е Р Н И Х З А Д А Ч С Т А Ц І О Н
А Р Н О Й Т Е П Л О П Р О В О Д Н О С Т І
Визначення температурного поля в будь-який момент часу єосновним завданням теорії теплопровідності. Для ізотропного тіла (зпостійним по різних напрямках коефіцієнтом теплопровідності () вонаможе бути описана диференціальних рівнянням теплопровідності
? T + Qv/(= 1/a * (dT/d (()),
(1)
де Т - температура; а - коефіцієнт температуропроводності, а = (/ ((* c);
(- Щільність матеріалу, с - питома теплоємність при постійномутиску,? -позначення оператора Лапласа (? = d/dx + d/dy + d/dz - вдекартових координатах x, y, z); (- час, Qv - об'ємна щільністьтеплового потоку.
Рівняння теплопровідності є математичним виразом законузбереження енергії в твердому тілі.
При вирішенні завдання до диференціальних рівнянь теплопровідностінеобхідно додати крайові умови. В опис крайових умов входять:поле температур для якого-небудь попереднього моменту часу
(початкові умови), геометрія тіла (геометричні умови),теплофізичні характеристики тіла (фізичні умови) і законтеплообміну між поверхнею тіла і навколишнім середовищем (граничніумови).
Якщо процес теплопровідності не тільки стаціонарний
(dT/d (tay) = 0), але і відбувається без тепловиділення усередині матеріалу (Qv =
0), то рівняння приймає вигляд
? (Т) = 0.
(2)
Зважаючи на складність і трудомісткість рішення неодновимірному завданьтеплопровідності аналітичними методами в інженерній практиці найбільшчасто використовують наближені. Один з них - метод кінцевих різниць,безпосередньо базується на диференціальному рівняннітеплопровідності і граничних умовах, становить найбільший інтерес.
В даний час значного поширення набули кінцево -різницеві методи, побудовані з використанням відомих законівзбереження. У цьому випадку різницеві схеми отримали назвуконсервативні. Такий підхід до побудови схеми, який зберігає фізичнусутність завдання, переважно суто аналітичного підходу,що полягає в безпосередній запису диференціальних рівняньзвичайно-різницеві аналогами.
Слід зауважити, що теорія кінцево-різницевих чисельних методівє самостійним розділом обчислювальної математики і широкопредставлена в спеціальній літературі [1,2,]. З основними методамипобудови кінцево-різницевих схем, алгоритмами розрахунку, програмнимзабезпеченням стосовно задач теплообміну можна ознайомитися внавчальній літературі [3,4,5].
При викладенні вказаного методу особливу увагу приділено фізичнійзмістом побудови консервативної різницевої схеми та її реалізації на
ПЕОМ в задачах теплопровідності.
При використанні чисельного методу з консервативною різницевоїсхемою тверде тіло розбивають на елементарні обсяги. Передбачається,що маса такого елементарного об'єму зосереджується в його центрі,званому вузлом. Для кожного вузла на основі закону збереження енергіїскладається рівняння теплового балансу, яке включає значення всіх теплових потоків на межах обсягів (осередків). Якщо осередок прилягає доповерхні тіла, то вирази для визначення теплових потоків маютьописувати теплообмін між тілом і навколишнім середовищем, тобто враховуватиграничні умови. Після виконання перетворень з рівняннямитеплового балансу отримують алгебраїчні рівняння для температури вкожному вузлі. Оскільки число вузлів і число комірок співпадають, тоутворена система алгебраїчних рівнянь є звичайно -різницевим аналогом диференціального рівняння теплопровідності ізамінює його з відповідними граничними умовами. Такий підхід доскладання кінцево-різницевого аналога, пов'язаний з тепловим балансом,дозволяє отримувати правдоподібні рішення навіть при грубому виборівідстані між вузлами (розміру комірки сітки).
Розглянемо деякі конкретні приклади складання кінцево -різницевих схем для вузлів двовимірної задачі теплопровідності. У цьомувипадку рівняння (2) набуває вигляду dT/dx + dT/dy = 0.
(3)
Внутрішня область типового двовимірного тіла показана на рис.1. < p> Рис.1. Розташування вузла всередині двовимірного тіла товщиною б.
Кожен елементарний прямокутник (комірка сітки) має довжину-х івисоту-у в напрямках осей х та у. Внутрішній вузол, позначенийсимволом 0, оточений чотирма сусідніми вузлами: 1,2,3,4. Кондуктивнаперенесення теплоти, що насправді відбувається в твердому тілічерез поверхні y * б и x * б (б-товщина тіла) будемо вважати як перенесеннятеплоти від відповідних вузлів до центрального. У сталихумовах рівняння балансу теплових потоків для вузла 0 за відсутностівнутрішнього тепловиділення буде мати вигляд
Q (1-0) + Q (2-0) + Q (3-0) + Q (4-0)
= 0, (4)
де Q (I-0) - тепловий потік; індекс (I-0) вказує напрямок перенесення ввузлах.
Для визначення кондуктивної теплового потоку може бути застосованийзакон Фур'є
Q = - lamda * F *dT/dn, (5)
де Т - температура, n - напрямок перенесення теплового потоку, F --поверхню, через яку переноситься тепловий потік.
для побудови розрахункової схеми градієнт температури у виразі (5)замінимо різниці температур в сусідніх вузлах. У цьому випадку перший членвираження (4) набуде вигляду
Q (1-0) = y * б * (T [1] -
T [0])/x. (6)
Тут градієнт температури визначається на кордоні двох вузлів 1 та 0,мають температури відповідно Т [1] і Т [0].
Аналогічні рівняння можуть бути отримані і для інших трьохчленів рівняння (1):
Q (2-0) = x * б * (T [2] - T [0])/y,
(7)
Q (3-0) = y * б * (T [3] - T [0])/x,
(8)
Q (4 -0) = x * б * (T [4] - T [0])/y.
(9)
Точність апроксимації градієнта залежить від розміру осередку. Якщоосередок має квадратну форму, то рівняння теплового потоку стаєнезалежним від форми тіла.
Підставляючи залежності (6 )...( 9) у вираз (4), можна побачити,що при постійному коефіцієнт теплопровідності для квадратної сітки
(x = y) воно зводиться до співвідношення між температурами в розглянутомувузлі та прилеглих:
T [1] + T [2] + T [3] + T [4] - 4 * T [0] = 0.
(10)
Вираз (10) застосовується до всіх внутрішніх вузлів.
Розглянемо вузол, розташований на поверхні твердого тіла, товщиноюб в двомірної задачі (рис.2).
Ріс.2.Расположеніе вузлів на поверхні двовимірного тіла, омиваного рідиною
Нехай вузол 0, розташований на кордоні твердого тіла,контактує з навколишнім середовищем, що має температуру ТC. Інтенсивністьтеплообміну з навколишнім середовищем характеризується коефіцієнтом тепловіддачіalfa. Вузол 0 може також обмінюватися кондуктивна потоком теплоти зтрьома сусідніми вузлами: 1,2,3. У цьому випадку тепловий баланс для вузла 0 запишеться наступним чином:
Q (1-0) + Q (2-0) + Q (3-0) + Q (c-0) = 0,
(11)
де Q (c 0)-тепловий потік, що передається від середовища вузлу 0 конвекцією.
За законом Ньютона - Ріхман < p> Q (c-0) = alfa * F * (T [c] - T [0]).
(12)
В результаті перетворень виразу (11), за аналогією з ранішевиконаними, для внутрішнього вузла, отримаємо
y * б * (T [1]-T [0])/x + (x/2) * б * (T [2]-T [0] )/y
+ (X/2) *
* б * (T [3]-T [0])/y + alfa *y * б * (Tc-T [0]) = 0. (13)
Співвідношення (13) значно спрощується при виборі квадратноїсітки. У цьому випадку при постійному коефіцієнт теплопровідності воноприводиться до виду
T [1] + 0,5 * (T [2] + T [3]) + Bi * Tc - (2 + Bi) * T [0] =
0, (14)
де Bi = alfa * x/lamda - число Біо.
Нижче наведені рівняння теплового балансу при інших граничнихумовах для двовимірних тел (x = y):
Вузол Схема
Розрахункове
рівняння
.....|/ Т
. 2 */Е
. |/П
Плоска поверх--+--.---- + - |/Л ність з тепло-|. |/Про ізольованою x. * -? - * |/І кордоном |. 1 0 |/З
0,5 (T [2] + T [3]) +
-+--.---- + - -?/Про
+ T [1] -2 * T [0] = 0
. |/Л
. -> + X ++---? x +
