CC BY
10УДК 519.63:536.2
Математическое моделирование переходных термических процессов с учётом явлений тепловой релаксации и термического демпфирования
Л. С. Петрова, Е. В. Заец
В статье представлена математическая модель с двухфазным запаздыванием на основе уравнения теплопроводности гиперболического типа, учитывающего явления тепловой релаксации и термического демпфирования. Получено численное решение гиперболической задачи теплопроводности с учётом конечной скорости распространения тепла и демпфирования температуры. Рассмотрена реализация метода сеток с использованием трёхслойной неявной разностной схемы при решении задачи нестационарной теплопроводности для переходных термических процессов. Описан алгоритм расчёта температурного поля с использованием уравнения с двухфазным запаздыванием. Представлены результаты расчётов температурных полей в биологических тканях по уравнению теплопроводности гиперболического типа для расчётной области в пределах двухфазной зоны.
Ключевые слова: математическая модель, уравнение теплопроводности гиперболического типа, модель с двухфазным запаздыванием, численные методы, метод прогонки.
1. Введение
В настоящее время исследование быстропротекающих (высокоскоростных) переходных термических процессов является перспективным направлением в области интенсификации технологических процессов с достижением сверхвысоких и сверхнизких температур, использования лазерного облучения и материалов со сложной структурой.
В современных научных исследованиях (Ю. А. Кирсанов, А. Ю Кирсанов, А. Е. Юда-хин [1], Л. А. Фильштинский, Т. А. Киричек [2], Э. М. Карташов [3], H. Z. Poor, H. Moosavi, A. Moradi [4] и др.) отмечается ограниченность применения для быстропротекающих переходных процессов классической теории теплопроводности, основанной на гипотезе Фурье, содержащей допущение о бесконечной скорости распространения тепла, и приводится более сложная зависимость между тепловым потоком и градиентом температуры, которая учитывает тепловую релаксацию и термическое демпфирование [1, с. 70]:
q^т = , (1)
^ = -ЛУ Т-тТ ЛдУТ дт дт
где Тц, тт - два несвязанных времени запаздывания (т^ - время тепловой релаксации, тт -время термического демпфирования), д - плотность теплового потока, УТ - градиент температуры, т- время, Л- теплопроводность.
дТ
Сочетание формулы (1) и уравнения теплового баланса ер— = -ёгуд позволило иссле-
дт
дователям получить гиперболическое уравнение теплопроводности, учитывающее скорость изменения плотности потока тепла и скорость изменения градиента температуры [5, с. 550]:
дТ д 2Т дТл
Я
я = aAl T + , (2)
От 4 От2 v От J
где а = — - коэффициент температуропроводности, ДТ = У2Т - div(grad Т) - оператор еР
Лапласа.
В исследованиях российских авторов для уравнения (1) используется название - уравнение Максвелла-Каттанео-Лыкова [6, с. 18] или уравнение с двухфазным запаздыванием [2, с. 164], иностранными исследователями это уравнение называется уравнением двухфазного запаздывания (ёиа1-рЬа8е-1а§) [5, с. 550].
В настоящее время исследователями рассматриваются математические модели с запаздыванием в различных прикладных областях науки и техники. Использование дифференциальных уравнений в частных производных с эффектом запаздывания при описании математических моделей динамических процессов в электродинамике, механике деформированного твердого тела, медицине, экологии и экономике отражено в работах Ю. Ф. Долгия, П. Г. Суркова [7], В. Г. Пименова, А. Б. Ложникова [8], А. И. Цаплина, И. Л. Никулина [9], Н. А. Быкова, О. К. Овчинниковой [10] и других.
Использование аналитических методов решения задач нестационарной теплопроводности с двухфазным запаздыванием на основе интегральных преобразований представлено в работах [1-5, 11-14]. Применение численных методов для решения гиперболической задачи теплопроводности с учётом конечной скорости распространения тепла и демпфирования температуры остается мало исследованным [9, 10]. В работе Н. А. Быкова и О. К. Овчинниковой [10, с. 144-146] реализуется метод ловли фронта в узел сетки с переменным шагом по времени (метод с выделением границы раздела фаз), при этом математическая постановка задачи сводится к системе двух уравнений теплопроводности параболического типа. Исследователями А. И. Цаплиным, И. Л. Никулиным [9, с. 117-121] используется метод сквозного счета, но с введением эффективной теплоемкости уравнение теплопроводности для расчетной области в пределах двухфазной зоны приводится к стандартному виду. Исследований по тематике применения численных методов с реализацией неявной разностной схемы при решении задач гиперболической теплопроводности с учётом явлений тепловой релаксации и термического демпфирования нами не обнаружено.
2. Задача нестационарной теплопроводности для переходных термических процессов в биологических тканях с учётом конечной скорости распространения тепла и демпфирования температуры
Значительный интерес у исследователей вызывает математическое моделирование быст-ропротекающих переходных термических процессов в биологических тканях. Реализация аналитических методов решения задач теплопроводности в условиях быстропротекающего термического процесса в материалах с низкой теплопроводностью представлена в исследованиях F. Filipoiu, A. I. Bogdan, I. M. Cârstea [11], H. Z. Poor, H. Moosavi, A. Moradi [4], Y. Zhang [12], K. C. Liu, P. J. Cheng, Y. N. Wang [1з], E. A. N. Al-Lehaibi [14]. В работах [4] и [13] описано применение метода разделения переменных и операционного метода для расчёта температурного поля с использованием уравнения с двухфазным запаздыванием. Матема-
тическое моделирование на основе уравнения теплопроводности параболического типа для расчётной области в пределах двухфазной зоны представлено в исследованиях [11] и [12]. Использование операционного метода при решении задачи теплопроводности с двухфазным запаздыванием на основе перехода к системе двух уравнений теплопроводности (параболического и гиперболического типов) описано в работе [14]. Использование операционного метода при решении краевых задач с уравнением теплопроводности гиперболического типа, учитывающим явления тепловой релаксации и демпфирования температуры, в случае постоянного, периодического и импульсного тепловых потоков на границе сопряжено с возрастанием трудоёмкости выполняемых преобразований и актуализирует проблему использования численных методов решения.
Математическая модель быстропротекающего процесса теплопередачи в тканях кожи с постоянным тепловым потоком на поверхности включает в себя уравнение теплопроводности с двухфазным запаздыванием [4, с. 1460]:
д2Т (х,т) , чдГ (х,т) / / \ \
—^^ + (рс + тчрьсьшь ) —+ рьсьшь (т (х,т)-та) =
TqPC- 2
дт
= я
д2Т (х, т) д3Т (х, т)
--Ь Тт-
дх2
+ Q + Q Ьт дQmet ЬТ *°Qext Ь Qmet Ь Qext bTq Ь т
начальные условия:
дх дт 0 < х < L, 0 <т<тт;
Т (х, 0) = Та, дТ (х,т)
дт
дт
дт
т=0
= 0, х е[0, L];
условия на границах:
-1
дТ (х,т)
дх
х=0
= qo,
дТ (х,т)
дх
= 0, 0 <т<тт,
х=L
(3)
(4)
(5)
(6) (7)
где Т(х, т)- температура ткани в точке х в момент времени т , Я — коэффициенты теплопроводности ткани, c — теплоемкость ткани, р — плотность ткани, cъ — теплоемкость крови, P — плотность крови, шъ — скорость перфузии крови, Та — температура крови, Qmet — метаболическое тепловыделение в ткани, Q^ — мощность внешних тепловых источников, q0 — плотность теплового потока. В этом исследовании Qmet = const, Qe:ct = 0.
3. Метод сеток с использованием трёхслойной неявной разностной схемы при решении гиперболической задачи теплопроводности для расчётной области в пределах двухфазной зоны
Реализуя метод сеток для решения задачи (3) - (7), вводим прямоугольную пространственно-временную сетку Охт = |х =(/ -1) к, / = 1, N; = ]к, ] = 0,М|. Для получения трёхслойной неявной разностной схемы используем формулы соответствующих разностных ана-
логов частных производных первого и второго порядков по времени, производной второго порядка по пространству, представленные в работе [15, с. 33-36]. Применение теории разностных схем, описанной в работе [15, с. 71-80], позволило получить разностный аналог для
смешанной производной третьего порядка с погрешностью аппроксимации O (АЛ2):
д3Т (*,т), Т-Л1 - Т+1 - 2 Т+1 - Т/) + Т--+1 - Т-1
Iх,т) _
2 ~ 2 дх дт к к
(8)
Разностное уравнение в этом случае имеет вид:
т/+1 - 2Т/ + Т/-1 , ,Т/ +1 - Т/ ( /+1 \
т-Ре-2-+ (Ре + ТЧРЪеЪ™Ъ )-]-+ РЪеЪ™Ъ (Т - Та )
к к
= Я
гт{+11 - 2Т/+1+Т-I1 ТУ - Т{1 - 2 (Т/+1 - Т/)+Т-+1 - Т-1
-9-+ ТТ -9-
V Л Л к
+ О^швг, (9)
где / = 2, N -1, ] = 1, М -1.
Модифицируя метод прогонки для семиточечной разностной схемы, систему (9) приводим к системе линейных алгебраических уравнений, имеющей трёхдиагональную структуру:
А—- — - - Т— +1+ С - — = ^, (10)
где А- = с—- =я( к2 + ктТ), В— = к2 (т- + к)(ре + крЪеЪшЪ) + 2Як(к + тТ), ^ = к2тцреТ/ 1 +
+ЯтТк (т+ Т/+1) - т/ 2ЯтТк + 2к1тчРе + кк2 (ре + т-РъеЪЩ ) - к2к2 (РЪеЪшЪТа + Qmet) .
Используя основное соотношение прогонки Т-[+1 =аг-_^ -Т/+1 + Д-_ 1 для уравнения (10), получаем формулы для расчёта прогоночных коэффициентов:
Щ =-А-, р. = С -1 - . (11)
В— - с— -щ-1 В - с -Щ-1
Определяя коэффициенты Г— на первом временном слое, используем соотношение для
фиктивных узлов Т 1 = Т0, полученное из аппроксимации второго начального условия:
^ =Яттк(т-1 + Т°-1)-Т? 2Яттк + к1тчРе + кк2 (Ре + тчрЪеЪЩ )]-к2к2 (рЪеЪ™ЪТа + Qmet) (12)
Для получения начальных прогоночных коэффициентов щ, Д используем разностное уравнение (9) и соотношение для значений сеточной функции в фиктивных узлах, полученное из аппроксимации левого граничного условия с погрешностью О (к2):
Щ+1 = ^-2ЯкТ-, (13)
к (т- + к)(ре + АръСъШъ) + 2Як(к + т)
,+1 2кк2-0 + 2ЯтТк(т/ -Т{ ) + к2ре т/ (2тц + к)- т/ 1т-
+ =-ж-
к (т- + к)(ре + кръеъШъ) + 2Як(к + Тт) кк2рЪеЪшЪ (Т\ т- + кТа) + йи^2
+ ~2-------. (14)
к (т- + к)(ре + кръеъШъ) + 2Як(к + т)
Начальные прогоночные коэффициенты на первом временном слое определяем, используя соотношение для фиктивного узла 7} 1 = Т\0 :
2Лттк (Т0 - Т20 ) + к^реТ} (тд + к) + кк2рьоьшь (т\% + кТа ) + кк2 (2д0 + (¿тк)
Р\ =-9-. (15)
к (Гц + к^(ре + кръеътъ) + 2Лк(к + Тт)
Формулы для вычисления значений температуры на правой границе получаем, используя аппроксимацию правого граничного условия с погрешностью О (к2) и значение сеточной
функции в фиктивном узле Т— = т] на первом временном слое:
+1 _ кРе-2тдтм -кЧ)-кк2РъеЬшЬ + кТа) ] 2Лк(тт + к)(а^-\ -\)-к2(тд + к)(ре + крьеьШь)
2Лттк (ц-\ - Ц)- (тк2к2 - 2Лкры-\ (тт + к)
+-———2-; (16)
2Лк(тт + к)(&]-\ -\)-к (тц + к^(ре + крьеьШь)
' 2 „ гтО ( , , \ ,, 2 ( гтО
т\ _ -к рет° (тц + к) - кк рьеьть (тцтЫ + кта
2Лк(тт + к-\ -\)-к2 (тд + к)(ре + крьеь^ь)
| 2Лттк(т°-т0)-Отк2к2 -2Лкрн(тт + к) 2Лк(тт + к){а -\)-к2 (тд + к)(ре + крьеьШь)
При реализации модификации метода прогонки, предусматривающей использование процедуры расчёта значений сеточной функции на промежуточном и первом временных слоях, на этапе прямой прогонки для определения прогоночных коэффициентов используются формулы (11), (13) - (15). Расчёт температурного поля на этапе обратной прогонки осуществляется с применением формул (16) и (17) для вычисления значений температуры на правой границе, а также с применением основного соотношения прогонки для остальных узлов на текущем временном слое.
Отличие знаменателей прогоночных коэффициентов от нуля и выполнение условия
| < \, I = \, N - \, исключающего быстрый рост погрешности округления, обеспечивает успешное применение метода прогонки с учётом безусловной устойчивости предложенной неявной разностной схемы с погрешностью аппроксимации О (к + А2 ) .
4. Результаты расчётов температурных полей с использованием уравнения с двухфазным запаздыванием
Реализация программы численного решения задачи (3) - (7) осуществлялась в системе MathCAD, в среде программирования Visual Studio Express Edition на языке C# и в среде Dev-C++.
Сравнение и анализ результатов расчётов температурных полей на основе предложенного численного решения проводились при следующих значениях теплофизических характеристик и геометрических условиях [4, с. 1465, 13, с. 65]: Л = 0,235 Вт/(м• К), р = 1190 кг/м3 ,
ръ = 1060 кг/м3 , с - 3600 Дж/(кг • К), сь = 3770 Дж/(кг • К), Wb = ръшъ = 0,5 кг/(м3с),
Та = 37 °С, 0?иег - 368,1 Вт/м3 , ^ = 5000 Вт/м2 , = 2000 Вт/м2 , £ = 6 мм, тт = 10 с,
тп = 16 с.
ч
Результаты расчётов температурных полей с постоянным тепловым потоком на поверхности 41 = 5000 Вт/м2 (Т1), 42 = 2000 Вт/м2 (Т2) представлены в виде графиков распределения температуры по толщине слоя при т = 15 с, т = 45 с и функции изменения температуры при х = 0,1 мм, х = 1,6 мм на рис. 1 - 2.
70 65 (Т10)160
еч
55
(Т20)1
50
■ * • •
(Т21)1 ... 45
40
35
- 3 - 3 - 3 - 3 - 3 - 3
1x10 3 2x10 3 3x10 3 4x10 3 5x10 3 6x10 3
Рис. 1. Распределение температуры по толщине слоя: при т = 15с (Т10, Т20) и при т=45с (Т1ц, Т2^)
115
105
(ТЧ 95
(ТЧ 85
Ш » » 9 (ТЧ 75
(Т2,). 65
55
45
35
0 20 40 60 80 100 120 140 160 180 200
Рис. 2. Изменение температуры во времени: при х=0,1мм (Т10, Т20 ) и при х=1,6мм (Т1ц, Т2^)
0
х
т
Сопоставимость представленных результатов с результатами расчётов, полученными в работе [4, с. 1469-1471] с применением аналитических методов на основе метода разделения переменных, обосновывает достоверность полученных результатов.
5. Заключение
Исследование явлений теплопередачи в условиях быстропротекающего термического процесса в материалах с низкой теплопроводностью (например, в биологических тканях) имеет большое значение и может найти разнообразное применение, в том числе для разработки медицинских технологий, конструирования оборудования и многих других целей.
Разработанная математическая модель с двухфазным запаздыванием может использоваться при анализе технологических процессов, основанных на теплопередаче быстропроте-кающих переходных термических процессов в тканях кожи с постоянным тепловым потоком на поверхности, для расчётов температурных полей в биологических тканях с использованием уравнения теплопроводности гиперболического типа, учитывающего явления тепловой релаксации и термического демпфирования.
Литература
1. Кирсанов Ю. А., Кирсанов А. Ю., Юдахин А. Е. Решение краевой гиперболической задачи теплопроводности для переходного термического процесса // Материалы X школы-семинара молодых учёных и специалистов академика РАН В. Е. Алемасова «Проблемы тепломассообмена и гидродинамики в энергомашиностроении». 2016. С. 70-73.
2. Фильштинский Л. А., Киричек Т. А. Применение неклассических моделей теплопроводности для расчётов тепловых полей в элементах конструкций // Авиационно-космическая техника и технология. 2005. № 7 (23). С. 162-170.
3. Карташов Э. М. Математические модели теплопроводности с двухфазным запаздыванием // Инженерно-физический журнал. 2016. № 2 (89). С. 338-349.
4. Poor H. Z., Moosavi H., Moradi A. Analysis of the DPL bio-heat transfer equation with constant and time-dependent heat flux conditions on skin surface // Thermal Science. 2016. V. 20, № 5. P.1457-1472.
5. Кирсанов Ю. А. Влияние тепловой релаксации и термического демпфирования на переходные процессы при циклических граничных условиях // Теплофизика высоких температур. 2017. № 4 (55). С. 549-555.
6. Юдахин А. Е. Исследование моделей теплопроводности в условиях быстропротекающего термического процесса в низкотеплопроводном твердом теле: дис. ... канд. техн. наук. Казань. 2017. 141 с.
7. Долгий Ю. Ф., Сурков П. Г. Математические модели динамических систем с запаздыванием: учебное пособие. Екатеринбург: Изд-во Урал. ун-та, 2012. 122 с.
8. Пименов В. Г., Ложников А. Б. Численный метод моделирования управляемого уравнения теплопроводности с запаздыванием // Вестник Тамбовского университета. Серия. Естественные и технические науки. 2013. № 5-2 (18). С. 2635-2636.
9. Цаплин А. И., Никулин И. Л. Моделирование теплофизических процессов и объектов в металлургии: учебное пособие. Пермь: Изд-во Перм. гос. техн. ун-та, 2011. 299 с.
10. Брыков Н. А., Овчинникова О. К. Численное решение задачи плавления твердого вещества // Успехи современной науки. 2017. № 3 (6). С.144-147.
11. Filipoiu F., Bogdan A. I., Cárstea I. M. Computer-aided analysis of the heat transfer in skin tissue // Proceedings of the 3rd WSEAS Int. Conference on finite differences - finite elements -finite volumes - boundary elements. C. 53-59.
12. Zhang Y. Generalized dual-phase lag bioheat equations based on nonequilibrium heat transfer in living biological tissues // International Journal of Heat and Mass Transfer. 2009. V. 52. P. 4829-4834.
13. Liu K. C., Cheng P. J., Wang, Y. N. Analysis of non-Fourier thermal behavior for multi-layer skin model // Thermal Science. 2011. V. 15 (1). P. 61-67.
14. Al-Lehaibi E. A. N. The skin tissue of the human head subjected to thermal diffusion // Mathematical Problems in Engineering. 2018. [Электронный ресурс]
URL: https://doi.org/10.1155/2018/8781950 (дата обращения 08.02.2019).
15. Самарский А. А. Теория разностных схем: учебное пособие для вузов. М.: Наука, 1977. 656 с.
Статья поступила в редакцию 19.02.2019.
Петрова Лилия Сергеевна
к.пед.н., доцент по специальности «Математическое моделирование, численные методы и комплексы программ», доцент кафедры «Высшая математика», Омский государственный университет путей сообщения (ОмГУПС) (644046, Омск, пр. Маркса, 35), e-mail: petrov.306@mail.ru.
Заец Евгений Валерьевич
магистрант кафедры «Теплоэнергетика», Омский государственный университет путей сообщения, e-mail: zaetsevgenijmagistr@gmail.com.
Mathematical modeling of transition thermal processes taking into account effect of thermal relaxation and thermal damping.
L. S. Petrova, E. V. Zaets
The article presents a mathematical model with dual phase lag based on the heat conduction equation of hyperbolic type, taking into account effect of thermal relaxation and thermal damping. The numerical solution of the hyperbolic heat conduction problem is obtained taking into account the final speed of heat spread and temperature damping. The implementation of the grid method using a three-layer implicit difference scheme in solving the problem of non-stationary heat conductivity for a transition thermal processes were treated. An algorithm for calculating the temperature field using a dual phase lag equation is described. The results of calculations of temperature fields in biological tissues according to the hyperbolic thermal conductivity equation for the calculation area within the dual phase zone are presented.
Keywords: mathematical model, heat conduction equation of hyperbolic type, model with dual phase lag, numerical methods, sweep method.
