Исследование процесса теплопроводности в пластине с внутренними источниками теплоты постоянной мощности Текст научной статьи по специальности «Математика»

Исследование процесса теплопроводности в пластине с внутренними источниками теплоты постоянной мощности

А.В. Еремин, К.В. Губарева, А.С. Шульга Самарский государственный технический университет

Аннотация: На основе совместного использования дополнительных граничных условий (характеристик) и интегрального метода теплового баланса получено приближенное решение краевой задачи теплопроводности с внутренними источниками теплоты постоянными во времени при граничных условиях третьего рода. Рассматривая тепловой поток на поверхности пластины, как новую искомую функцию, получено простое по форме аналитическое решение указанной задачи. Использование предлагаемого подхода возможно для решения дифференциальных уравнений в частных производных, не допускающих разделение переменных. Показано, что скорости распределения устремляются в бесконечность и имеют начальную скорость в определенный момент времени.

Ключевые слова: аналитическое решение, внутренние источники теплоты, нестационарная теплопроводность, краевая задача, граничные условия третьего рода, дополнительная функция, интеграл теплового баланса.

Возникая вследствие различных факторов (экзотермические химические реакции, протекание электрического тока по проводнику и др.) [1 - 4] внутренние источники теплоты могут оказывать существенное влияние на тепловые процессы в технике. В зависимости от соотношения мощности источников теплоты и интенсивности теплоотвода с поверхности, температура тела с течением времени может либо устанавливаться на некотором уровне, либо неограниченно возрастать (тепловое воспламенение) [4], либо снижаться (охлаждение).

Серьезные трудности представляет математическое описание процессов переноса теплоты в телах с внутренними источниками теплоты. Для получения решений краевых задач теплопроводности, описывающих данные процессы, используются точные (интегральных преобразований, метод Фурье и др. [5, 6]), приближенные аналитические методы (метод Галеркина, коллокаций [7, 8], интегральные методы теплового баланса [9 -10]), а также численные методы. Применяемые ранее методы имеют

1

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

В настоящей работе на основе совместного применения метода дополнительных граничных условий в интегральном методе теплового баланса и дополнительной искомой функции (теплового потока) отображены результаты разработки метода решения краевых задач теплопроводности для бесконечной пластины.

В данной работе рассмотрен метод на конкретном примере, для бесконечной пластины при граничных условиях III рода при нестационарном режиме и в результате воздействия источника теплоты постоянного во времени.

Математическая постановка задачи в данном случае имеет вид

дТ(х, т) д 2Т(х, т)

—= а-+ — (т> 0; 0 < х <8); (1)

дт дх су

Т (х, 0) = Т„; (2)

дТ Т) = 0, (3)

дх

+ а(Тср - Т(8, т)) = 0, (4)

где Тср - температура среды, К; ю - удельная мощность внутреннего

источника теплоты; х - координата; а - коэффициент теплоотдачи; у -плотность; т - время; X - коэффициент теплопроводности; а = Х/(ср) -коэффициент температуропроводности; Т - температура; р - плотность; с -удельная изохорная теплоемкость; Т0 - начальная температура, К; 8 -половина толщины пластины, м; X - коэффициент теплопроводности.

Задача (1) - (4) может быть представлена в безразмерном виде [5] (см. рис. 1)

д0(^,Ро) д 20(^,Ро)

д Ро 0&0) = 1; д©(0, Ро)

д^2

+ Ро (Ро > 0; 0 < £ < 1)

д^

д0(1, Ро)

д^

= 0;

+ Bi0(1, Ро) = 0:

(5)

(6)

(7)

(8)

где 0 = (Т - Тср /(Т0 - Тср) - безразмерная температура; £, = х/5 - безразмерная

координата; Ро = (ат)/52- критерий Фурье (безразмерное время); В1 = а5/^ -число Био; Ро = (ю52)/[Л,(Тср - Т0)] - критерий Померанцева.

16,0 е 12,0

8,0

4,0

--___« ...

0.2

_

0,25

0,5

0,75 I ! ,0

Рис. 1. - Схема теплообмена В соответствии с предлагаемым методом, введем в рассмотрение новую искомую функцию времени

(9)

_ ч д0(1, Ро) ф(Ро) = —1-1 = а,

д^

где а - угол между касательной к графику функции 0(£,, Ро) в точке = 1 и

координатной осью.

Возвращаясь к размерным величинам, выражение (9) можно представить в виде

Ф(т) = -V ™ (10)

T0 - Тср дх

Плотность теплового потока по закону Фурье на поверхности пластины определяется выражением

*т) = -хМ. (11)

дх

Подставляя (11) в (10) получаем g

ф(т) =тт^—q(T) = kФ), (12)

Л(T0 - Тср )

где k = const - некоторый коэффициент, определяемый масштабом системы.

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

Решение задачи (5) - (8) представим в виде полинома n-ной степени:

©&Fo) = Х b (Fo) -1 , (13)

i=1

где п е N - натуральное число, соответствующее количеству членов ряда (13); Ь^о) - неизвестные коэффициенты, зависящие от безразмерного времени.

Для получения решения исходной задачи (5) - (8) в первом приближении ограничимся тремя слагаемыми (п=3) в выражении (13). Выполнив подстановку граничных условий (7) - (9) в полином п-ной степени (13), получим систему трех простейших алгебраических уравнений с коэффициентами Ь1 ^о) . Система имеет следующий вид

Ь2 = 0;

Ы(Ь1 + Ь2 + Ь3) + Ь2 + 2Ь3 = 0; Ь2 + 2Ь3 - ф^о) = 0,

Решая систему, находим коэффициенты

Ь1(Р0) =_ф(Р°)(2 + В1). ь2(Ро) = 0; ьз(Ео) =

ф(Бо)

2В1 ' 24 7 ' 34 7 2 С учетом полученных коэффициентов полином (13) запишется в таком

виде

0(^о) = ¿(О ф(Бо), (14)

с2 В1 + 2

где /1(С) = ----- функция, зависящая от координаты. Найденная

2 2В1

функция удовлетворяет граничным условиям (7) - (9) при любых значениях

ф(Р°).

Для приближенного удовлетворения исходного дифференциального уравнение (5), проинтегрируем его в пределах изменения пространственной координаты, т.е. составим интеграл теплового баланса [9 - 10]

д 20(О,Ро)

| ае(^о) =1

0 дБо * 0

+ Ро

. (15)

дО2

Вычисляя интеграл, получаем обыкновенное дифференциальное уравнение вида

йф(Бо) ( 1 1Л

+

V В1 3 у

+ ф(Бо) + Ро = 0, (16)

из решения, которого находим

Бо

1 1

-+-

ф(Бо) = С1е В1 3 _Ро, (17)

где С1 - константа интегрирования.

Подставляя (17) в (14), получаем

©(^о) = /1(С) (С^ _ Ро) (18)

^ 1 1

где К =--+ —.

В1 3

Для выполнения начального условия (6) потребуем ортогональности невязки к координатной функции и составим его невязку

1

| [0(£, 0) -1]/1(£) д^ = (Ы + 3)С1 - 3(Ро - Ы + РоЫ) = 0. (19)

0

Из решения уравнения (19) определим константу интегрирования

_ 3(Ро-Б1 + РоБ1) В й й С

С1 =-б—3-. Выражение (18) с учетом найденной константы

представляет собой решение задачи (5) - (8) в первом приближении и может быть записано в следующем виде

3(Ро -Б1 + РоБ1)

0(£, Бо)

Б1 + 2Л

ч 2 2Б1 J

Б1 + 3

-е Т - Ро

(20)

Графики изменения безразмерной температуры по формуле (20) представлены на рис. 2, 3.

Рис. 2. - График распределения температуры в пластине.

- - приближенное решение,------численное решение; 1,3 -

номер приближения; Б1 = 0,5; Ро = 15

14,С 0 10,5

7,0 3,5

0 0,25 0,5 0,75 Ро 1

Рис. 3. - График распределения температуры в пластине.

--приближенное решение,----численное решение; 1,3-

номер приближения; В1 = 0,5; Ро = 15 Невысокая точность полученного в первом приближении решения объясняется недостаточным количеством числом членов полинома (12). Для нахождения решения задачи (1) - (4) во втором приближении будем использовать шесть слагаемых в выражении (13). Кроме граничных условий (7), (8), (9) при нахождении неизвестных коэффициентов Ь (Бо) в выражении (13) будут использоваться граничные условия, физический смысл которых заключается в выполнении исходного дифференциального уравнения (1) и выражений, полученных после его дифференцирования в точках О = 0 и С = 1. Отметим, что в работах [9] показано, что если исходное уравнение выполняется лишь в граничных точках, то это приводит к его выполнению и внутри области.

Для получения решения задачи (5) - (8) во втором приближении будем использовать шесть членов ряда (13) Ь(Бо) (п = 6). В этом случае, для определения неизвестных коэффициентов Ьi (Бо) помимо условий (7), (8), (9) будем использовать три дополнительных условия, одно из которых будет получено для координаты С = 0 и два для С = 1.

Записывая уравнение (5) в точке = 1 с учетом граничных условий (8) и (9) найдем первое дополнительное граничное условие для второго приближения

1 дф(£,Бо) д 20(^Бо)

Б1 ^ ^ = (23)

2 Ро = 0. (21)

Б1 дБо д^2

Для того, чтоб найти второе граничное дополнительное условие необходимо продифференцировать по пространственной переменной исходное уравнение (5)

д 20(^о) = д 30(^о) (22)

дБод^ д^3 ( )

Соотношение (22) в точке = 1 с учетом граничных условий (8) и (9)

представляет собой второе дополнительное условие и записывается в виде

дф(^о) = д 30(^,Бо) дБо = д^3

Далее найдем третье дополнительное граничное условие. Для этого следует продифференцировать исходное уравнения (5) по пространственной переменной в = 0 . С учетом граничного условия (7), уравнение принимает следующий вид

д 30(03Бо) = 0. (24)

д^3

Для нахождения неизвестных коэффициентов Ь (Бо) поставим полином (13) в граничные условия первого приближения (7) - (9) и дополнительные условия (21), (23), (24). Неизвестные коэффициенты Ь (Бо) (п = 6) находятся из решения системы шести алгебраических уравнений

Ь,(Ро) = Б'2 - 8 дФ<Ро> -Ф<Ро><40 + 28Б'> - Ро; Ь2(Ро) = 0;

1 40Б1 дБо 40 5 2

Ь3(Ро) = ф(Бо) - + Ро; Ь4(Бо) = 0;

^^ 12Б1 дБо 2

1

,,„. В12 - 4 ^ф(Бо) Ро

Ь5<Ро) =^ВТ -ф<Ро) "У

В12 - 3 ¿ф(Го) ф(Бо) Ро Ь6(Бо) =---——^ + ^—1 + —

15В1 й^о 5 5 С учетом выше найденных коэффициентов и подставляя уравнение (13) в интеграл теплового баланса (15) получим дифференциальное уравнение второго порядка

Л2ф(Ро)(В12 - 9) - ^ф(Ро)(39В1 + 90) - 90В1ф(Ро) - 90В1Ро = 0. (25)

Решение, которого имеет вид ф(Бо) = С1 exp(K1Fo) + С2 exp(K2Fo) - Ро, (26)

3(13В1 + 740В137169В1Т7420В1+900 + 30)

где К1

К 2 =

2(В12 - 9)

3(13В1 - V 40В13 + 169В12 + 420В1 + 900 + 30)

2(В12 - 9)

Выражение (14) после подстановки в него (26) может быть представлено в виде

0(О, Бо) = /&) С1 exp(KlFo) + /&) С2 exp(K2Fo), (27)

где

/1(0)=

К1 В1К1 1

V 5В1 15 5 у В1К1 1 К1 7

С5 +

В1К1 К1 1

8

2В1 2

С4 +

К1 В1К1

2В1

40 В1 5В1 10'

/2(С) =

( К

+

К1 В1К2 1

V

5В1 15 5

С5 +

вк К 1

8

2В1 2

С4 +

2В1 12

В1К2 ^

2 +1

,, 2 В1К

С2 +

1К

7

у

40 В1 5В1 10

+1 12 у

С2 + :

+

Составляя невязку начального условия и требуя ортогональности невязки к каждой координатной функции _/!(£) и /2(0, получаем систему двух алгебраических уравнений

\ (0(^,0) - 1)Ш d^ = 0,

0

|(0(^,0) -1)Л(0 d^ = 0.

из решения, которой найдем С1, С2.

Например, при Б1 = 0,5 и Ро = 15, константы интегрирования равны С = 14,0 и С2 = 0,072.

Выражение (27) с учетом найденных констант интегрирования представляет решение задачи (5) - (8) во втором приближении. Для дальнейшего повышения точности необходимо увеличивать число членов ряда (13). Так, в третьем приближении будем использовать девять членов ряда, в четвертом - двенадцать и так далее. Для нахождения неизвестных коэффициентов будем использовать дополнительные граничные условия. Общие формулы для их определения имеют вид д2к _10(0,Ео)

■2 к-1

= 0;

д^2

1 дк-1ф(Бо) дк+10(1,Бо)

Б1 дБок-1 д^к+1 '

- Б. дк-10(1,Ео) = д2к-10(1,Ео)

1 д^к-1 = д£2к-1 ,

где к = 1, 2, 3 ....- номер приближения.

На рис. 2, 3 представлены графики изменения безразмерной температуры. Анализируя график, можно сделать вывод, что расхождения полученных результатов не превышает 3% в диапазоне 0,1 < Бо < да .

На рис. 4 представлен график функции ф(Бо). Отмечается, что точность определения искомой функции ф(Бо) (теплового потока на поверхности) значительно выше точности определения температуры. Расхождение результатов расчетов теплового потока ф(Бо) по формуле (9) с численным решением в диапазоне безразмерного времени 0,1 < Бо < да не превышает 5% уже в первом приближении.

о ф

0,75

1,50

-2,25

,00

3 1

0 0,125 0,25 0,375 Ро 0,5 Рис. 4. - Графики изменения функции ф(Бо):---приближенное

решение;-----численное решение; 1,3 - номер приближения;

Ы = 0,5; Ро = 15

Выражая координату как функцию времени Бо и температуры 0(^,Бо), соотношение (14) можно записать в виде

^(Fo) =

У

2Bi0(^, Fo) + 9(Fo)(2 + Bi)

(28)

Bi9(Fo)

Выражение (27) позволит для любых конкретных 0(£,,Fo) = const построить графики изотерм.

Безразмерные скорости движения изотерм и = d£, / dFo по координате в зависимости от времени Fo представляют собой первые производные по времени от выражения (28), а вторые производные - ускорения a = d/ dFo2. Формула скорости будет иметь вид

и^о) =

ф(Fo)Р®^ 0(Fo)

дРо

дРо

ф2(Ро)

"V

2_ + 20(Ро) +1 Bi ф(Ро)

Графики изменения изотерм и скоростей, использованием формул (28), (29) дано на рис. 5, 6.

(29)

полученных с

1,0

£ 0,75

0,50

0,25

0

@=9/ 10 / П/12 / 1 Ъ/ 14/

0,5 0,75 1,0 1,25 Бо 1,5 Рис. 5. - Графики изменения изотерм в пластине при граничных условиях 3-

го рода (В = 0,5; Ро = 15).

0,5 0,75 1,0 1,25 Ро 1,5 Рис. 6. - Распределение скоростей изотерм в пластине при граничных условиях 3-го рода (В = 0,5; Ро = 15) Изотермы имеют вид наклонных к поверхности (£, = 0) пластины

кривых. Каждая из изотерм имеет начальную скорость и0 и появляется на поверхности тела в определенный момент времени. Далее, под воздействием

постоянного во времени источника Po = 15 скорости изотерм устремляются к бесконечным значениям.

Исследование выполнено при финансовой поддержке РНФ в рамках научного проекта № 18-79-00171 и Совета по грантам Президента РФ в рамках научного проекта МК-2614.2019.8.

Выводы

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

2. Использованный метод особенно эффективен при определении плотности теплового потока на поверхности тела, погрешность определения которого в диапазоне времени 0,1 < Fo < да составляет не более 5% уже в первом приближении.

3. Полученные результаты могут быть применены на практике инженерами - энергетиками для составления теплового баланса, нахождения приближённого значения теплового потока и др.

Литература

1. Вилюнов В.Н. Теория зажигания конденсированных веществ. Новосибирск: Наука, 1984. 190 с.

2. Clerk Maxwell J. A Treatise on Electricity and Magnetism, 3rd ed., vol. 2. Oxford: Clarendon, 1892, Pp.68-73.

3. Франк - Каменецкий Д. А. Диффузия и теплопередача в химической кинетике. М.: Наука, 1967. 492 с.

4. Кудинов В.А., Еремин А.В., Стефанюк Е.В. Критические условия теплового взрыва в пластине с нелинейным источником теплоты // Проблемы машиностроения и надежности машин. 2016. №1. С. 44-49.

5. Лыков А.В. Теория теплопроводности. М.: Высшая школа, 1967. 600 с.

6. Цирельман Н.М. Прямые и обратные задачи тепломассопереноса. М.: Энергоатомиздат, 2005. 392 с.

7. Rao T.D., Chakraverty S. Modeling radon diffusion equation in soil pore matrix by using uncertainty based orthogonal polynomials in Galerkin's method // Coupled Systems Mechanics. 2017, vol. 6 (4), Pp. 487-499.

8. Nourgaliev R., Luo H., Weston B., Anderson A., Schofield S., Dunn T., Delplanque J.R. Fully-implicit orthogonal reconstructed Discontinuous Galerkin method for fluid dynamics with phase change // Journal of Computational Physics, 2016, vol. 305, Pp. 964-996.

9. Кудинов В.А., Еремин А.В., Кудинов И.В. Получение аналитического решения задачи Стефана с учетом абляции на основе определения фронта температурного // Инженерно-физический журнал. 2012. №6. С. 1332-1342.

10. Mitchell S.L., Myers T.G. Application of Standard and Refined Heat Balance Integral Methods to One-Dimensional Stefan Problems // Siam Review. 2010, vol. 52 (1), Pp. 57-86.

References

1. Vilyunov V. N. Teoriya zazhiganiya kondensirovannykh veshchestv [Ignition theory condensed substances]. Novosibirsk, Nauka, 1984. 190 p.

2. Clerk Maxwell J. A Treatise on Electricity and Magnetism, 3rd ed., vol. 2. Oxford: Clarendon, 1892, Pp.68-73.

3. Frank-Kamenetskiy D.A. Diffuziya I teploperedacha v khimicheskoy kinetike [Diffusion and Heat Transfer in Chemical Kinetics]. M.: Nauka, 1967. 481 p.

4. Kudinov V.A., Eremin A.V., Stefanyuk E.V. Journal of Machinery Manufacture and Reliability. 2016. № 1. Pp. 44 - 49.

5. Lykov A.V. Teoriya teploprovodnosti [Heat conduction theory], M.: Vysshaya shkola, 1967. 600 p.

6. Tsirelman N.M. Pryamyye i obratnyye zadachi teplomassoperenosa. Moscow, Energoatomizdat, 2005. 392 p.

7. Rao T.D., Chakraverty S. Coupled Systems Mechanics. 2017, V. 6 (4), Pp. 487-499.

8. Nourgaliev R., Luo H., Weston B., Anderson A., Schofield S., Dunn T., Delplanque J.R. Journal of Computational Physics, 2016, V. 305, Pp. 964-996.

9. Kudinov V.A., Eremin A.V., Kudinov I.V. Journal of Engineering Physics and Thermophysics. 2012. № 6. Pp. 1441-1452.

10. Mitchell S.L., Myers T.G. Siam Review. 2010, V. 52 (1), Pp. 57-86.