Modeling of freezing processes by an one-dimensional thermal conductivity equation with fractional differentiation operators

Abstract


We have studied the Stefan problem with Caputo fractional order time derivatives. The difference scheme is built. The algorithm and the program for a numerical solution of the Stefan problem with fractional differentiation operator are created. For the given entry conditions and freezing ground parameters we have obtained the space-time temperature dependences for different values of parameter α. The functional dependences of the interface motion for the generalized Stefan conditions depending on the value of α are estimated. Finally we have found that the freezing process is slowed down during the transition to fractional derivatives.

Full Text

Введение. Задача Стефана описывает явления тепломассопереноса в средах с фазовым переходом, сопровождающимся выделением или поглощением тепла. Большой прикладной интерес представляет обобщение задачи Стефана для сред, в которых не выполняется принцип локального равновесия, что приводит к необходимости учета особенностей теплопереноса на межфазной границе с учетом нелокальных эффектов по времени (эффект памяти) [1-3] и по пространству (эффект пространственных корреляций) [3, 4]. Одно из направлений обобщения неравновесной термодинамики связано с развитием концепции фрактала. Процессы переноса тепла в этом случае могут быть описаны на основе дифференциальных уравнений в производных дробного порядка [5, 6]. Ранее [3] нами было показано, что учет эффектов памяти приводит к зависимости координаты межфазной границы от времени. В настоящей работе нами была поставлена задача сравнить пространственно-временные температурные поля, получающиеся при варьировании дробного параметра α для производной по времени (производная Капуто). В дальнейшем сопоставление с экспериментальными результатами позволит определить данный коэффициент для различных грунтов. 1. Математическая постановка задачи. Влажная фрактальная структура находится с некоторой постоянной температурой T0 в талом состоянии. Внезапно в начальный момент времени на поверхности устанавливается температура Tc , которая ниже температуры замерзания Tз . При этом с некоторой переменной толщиной ξ = f (t) образуется промерзший слой. Нижняя подвижная граница имеет всегда температуру замерзания Tз , а на границе происходит фазовый переход, на что требуется теплота фазового перехода Qf . При этом верхняя граница талой зоны имеет постоянную температуру замерзания, а нижняя - температуру грунта на большой глубине. В качестве математической модели процесса промерзания во фрактальных структурах рассмотрим следующую задачу: ∂ 2 T1 (x, t) , 0 < x < ξ(t), t > 0, ∂x2 ∂ 2 T2 (x, t) α ∂0t T2 (x, t) = D2 , ξ(t) < x < L, t > 0, ∂x2 T (x, 0) = T0 , ∂T (L, t) T (0, t) = Tc , = 0, t > 0, ∂x  T1 = T2 = Tз , x = ξ(t) : ∂T ∂T α ξ(t), λ1 1 - λ2 2 = Qf ∂0t ∂x ∂x α ∂0t T1 (x, t) = D1 (1) (2) 377 Б е й б а л а е в В. Д., А л и в е р д и е в А. А., М а г о м е д о в Р. A., М е й л а н о в Р. Р., А х м е д о в Э. Н. где 0 < α 1, ρ - плотность грунта, t = τ /τ0 , x = ξ/ξ0 - безразмерные время и координата, τ0 , ξ0 - характерные время и масштаб, D1 = a1 τ0 /ξ02 , D2 = a2 τ0 /ξ02 - безразмерные коэффициенты температуропроводности, Q - количество тепла, выделяемое или поглощаемое в процессе таяния льда или замерзания воды, постоянные температуры Tc < Tз < T0 , α ∂0t T (x, t) = 1 Γ(1 - α) t 0 Tt (x, s) ds (t - s)α - частная дробная производная Капуто. Как известно, при поиске аналитических решений дробных дифференциальных уравнений возникают большие трудности. Поэтому в настоящее время при решении дробных дифференциальных уравнений применяют как аналитические, так и численные методы. Численным методам решения краевых задач для дифференциальных уравнений с производными дробного порядка посвящены работы [6-13]. Унифицированным методом приближенного решения дифференциальных уравнений, применимым для широкого класса уравнений математической физики, является метод конечных разностей (или метод сеток) [14]. 2. Численная математическая модель. Для решения задачи (1), (2) методом сеток вводим равномерную сетку по пространственной переменной xm = mh, m = 0, 1, . . . , M, h = L/N, и неравномерную сетку по времени tn+1 = tn + τn+1 , n = 0, 1, . . . , N - 1, t0 = 0, tN = tкон , τn+1 > 0. Для дробной производной Капуто имеет место аппроксимация [15] α ∂0t T (xm , tn ) = 1 Γ(2 - α)τn n 1-α 1-α k+1 k - T1,m )(tn-k+1 - tn-k ) + O(τ ), (T1,m (3) k=0 а для производной ∂ 2 T (x, t)/∂x2 имеет место аппроксимация n+1 n+1 + T n+1 Tm+1 - 2Tm ∂ 2 T (xm , tn ) m-1 = + O(h2 ). ∂x2 h2 (4) Шаг по времени τn+1 , n = 0, 1, . . . , N - 1, нужно выбрать таким образом, чтобы за этот временной промежуток (от tn до tn+1 ) граница фазового перехода сдвинулась ровно на один шаг пространственной сетки. Тогда, воспользовавшись разностной аппроксимацией дробной производной (3), можно записать α ∂0t ξ(t) 1 ≈ Γ(2 - α)τn+1 n 1-α (ξk+1 - ξk )(t1-α n-k+1 - tn-k ) ≈ k=0 h ≈ Γ(2 - α)τn+1 378 n 1-α (t1-α n-k+1 - tn-k ) = k=0 htn1-α . Γ(2 - α)τn+1 Моделирование процессов промерзания одномерным уравнением . . . Воспользовавшись (3), (4), получим следующую разностную схему: 1 Γ(2 - α)τn n k+1 1-α k (T1,m - T1,m )(t1-α n-k+1 - tn-k ) = k=0 m = 1, 2, . . . , m∗ - 1, T1 D1 n+1 n+1 n+1 (T - 2T1,m + T1,m-1 ), h2 1,m+1 m=1 = Tc , T1 m=m∗ = Tз , (5) где m = m∗ - граница фазового перехода. Тогда разностная схема второй части уравнения будет иметь вид 1 Γ(2 - α)τn n 1-α k+1 k - T2,m )(t1-α (T2,m n-k+1 - tn-k ) = k=0 m = m∗ + 1, . . . , M - 1, T2 D2 n+1 n+1 n+1 ), + T2,m-1 - 2T2,m (T h2 2,m+1 m=m∗ = Tз , ∂T2 ∂x m=M = 0. (6) После преобразований уравнения (5) и (6) соответственно запишутся следующим образом: n+1 n+1 n+1 = Fm , + Bm T1,m+1 - Cm T1,m Am T1,m-1 ∗ m = 1, 2, . . . , m - 1, T1 m=1 = Tc , T1 m=m∗ = Tз , (7) где Am = B m = Fm и D1 , h2 Cm = t1-α t1-α n 1 1 =- T1,m + Γ(2 - α)τn Γ(2 - α)τn 2D1 t1-α 1 + , h2 Γ(2 - α)τn (8) n-1 1-α k+1 k - T1,m )(t1-α (T1,m n-k+1 - tn-k ) (9) k=0 n+1 n+1 n+1 = Fm , Am T1,m-1 - Cm T2,m + Bm T2,m+1 m = m∗ + 1, . . . , M - 1, T2 m=m∗ где Am = B m = Fm = - D2 , h2 Cm = t1-α t1-α n 1 1 T2,m + Γ(2 - α)τn Γ(2 - α)τn = Tз , ∂T2 ∂x m=M = 0, 2D2 t1-α 1 + , h2 Γ(2 - α)τn (10) (11) n-1 k+1 1-α k (T2,m - T2,m )(t1-α n-k+1 - tn-k ). (12) k=0 Проведем дискретизацию граничного условия в случае x = ξ(t): λ1 T1,m∗ - T1,m∗ -1 T2,m∗ +1 - T2,m∗ htn1-α - λ2 = Qf , h h Γ(2 - α)τn+1 то есть λ1 Tз - T1,m∗ -1 T2,m∗ +1 - Tз ht1-α n - λ2 = Qf . h h Γ(2 - α)τn+1 379 Б е й б а л а е в В. Д., А л и в е р д и е в А. А., М а г о м е д о в Р. A., М е й л а н о в Р. Р., А х м е д о в Э. Н. Таким образом, τn+1 = Qf h2 t1-α n . λ1 (Tз - Tm∗ -1 ) - λ2 (T2,m∗ +1 - Tз ) Γ(2 - α) (13) Шаг по времени зависит от температуры. Поэтому поле температуры можно определить методом простой итерации. Решение систем (7) и (10), где коэффициенты определяются согласно равенствам (8), (9) и (11), (12), будем на каждом временном слое находить методом прогонки в виде Tm = Qm+1 Tm+1 + Dm+1 , m = M - 1, M - 2, . . . , 0, (14) где Qm+1 = Bm , Cm - Qm Am Dm+1 = Am C m - F m , C m - Q m Am m = 1, 2, . . . , M - 1. (15) Алгоритм численного решения задачи 1) 2) 3) 4) 5) 6) Вводим входные данные. Находим расчетный шаг сетки по координате h = L/M . 0 , m = 1, 2, . . . , M . Вводим начальное поле температуры Tm Задаем начальное положение границы фазового перехода k = 1. Если время t tкон , то переходим к пункту 17, иначе - к пункту 6. Сохраняем температурное поле с предыдущего временного слоя: n Tm = Tm , m = 1, 2, . . . , M. 7) Задаем положение границы фазового перехода k = k + 1. 8) Сохраняем в новом векторе поле температуры с предыдущей итерации: s n Tm = Tm , m = 1, 2, . . . , M. 9) Определяем шаг по времени по формуле (13). 10) Задавая прогоночные коэффициенты Q1 и D1 на основе левого граничного условия, согласно (15) находим прогоночные коэффициенты. 11) Задаем температурное поле на границе фазового перехода Tks+1 = Tз . s+1 , m = k - 1, . . ., 12) Находим температурное поле до фазовой границы Tm 2, 1, по формуле (14). 13) Задавая прогоночные коэффициенты Q1 и D1 на основе левого граничного условия, согласно (15) находим прогоночные коэффициенты. s+1 14) Находим температуру на правой границе TM на основе правого граничного условия. s+1 , m = M - 1, 15) Находим температурное поле после фазовой границы Tm . . . , k по формуле (14). s+1 - T s | 16) Если maxm |Tm ε, то увеличиваем переменную времени на m шаг τ и переходим к пункту 6, иначе переходим к пункту 9. 17) Выводим результаты. 380 Моделирование процессов промерзания одномерным уравнением . . . 3. Результаты и их обсуждение. Для анализа полученного решения нами были выбраны начальные условия, типичные для зимнего промерзания, и взяты табличные теплофизические характеристики [15, стр. 133]: глубина грунта L = 0.3 м, теплофизические характеристики промерзшей зоны грунта - λ1 = 2.7 Вт/(м · К), ρ1 = 917 кг/м3 , a1 = 2090 Дж/(кг · К); теплофизические характеристики талой зоны грунта: λ2 = 0.6 Вт/(м · К), ρ2 = 1000 кг/м3 , a2 = 4220 Дж/(кг · К); характерные температуры T0 = 293 K, Tз = 273 K, Tc = 268 K, теплота фазового перехода Qф = 3.32 · 105 Дж/кг, влажность грунта w = 1 кг/кг. На рис. 1 приведены графики численного решения задачи при различных значениях параметра α и в произвольно выбранный момент времени t = 130 000 сек. Как видно из рисунка, в точке замерзания (0 ℃)1 температурная зависимость претерпевает излом, причем с уменьшением показателя дробной производной точка замерзания лежит ближе к поверхности. Рис. 1. Графики численного решения задачи (1), (2) при различных значениях параметра α в момент времени t = 130 000 сек [Figure 1. Graphs of the numerical solution of the problem (1), (2) for different values of the parameter α at time t = 130 000 sec] На рис. 2 и 3 приведены графики зависимости температуры от времени на глубинах 0.0015 м и 0.075 м при различных значениях параметра дробной производной. Как видно из рисунков, промерзание на глубине 0.0015 м происходит практически сразу, в то же время на глубине 0.075 м при любом значении α температура достаточно долго держится выше точки замерзания. В обоих случаях даже незначительное уменьшение α приводит к существенному замедлению скорости охлаждения, что характерно для сред с фрактальной структурой [2]. В таблице приведены моменты времени, соответствующие фазовым переходам, для различных координат при значениях параметра α = 1 и α = 0.9. Данные для α = 1 и α = 0.9 соответствуют функциональной зависимости ξ(t) ≈ 0.00015t1/2 и ξ(t) ≈ 0.0002t0.45 соответственно. Таким образом, фазовую границу можно задать функциональной зависимостью ξ(t) ≈ σ(α)tα/2 , где 0 < α 1. 1 Здесь и далее мы для удобства приводим температуру в градусах Цельсия. 381 Б е й б а л а е в В. Д., А л и в е р д и е в А. А., М а г о м е д о в Р. A., М е й л а н о в Р. Р., А х м е д о в Э. Н. Рис. 2. Графики зависимости температуры от времени на глубине 0.0015 м в различные моменты времени и при различных значениях параметра α и Tc = -5 ℃ [Figure 2. Graphs of temperature versus time at a depth of 0.0015 m at different times and for different values of the parameter α and Tc = -5℃] Рис. 3. Графики зависимости температуры от времени на глубине 0.075 м при различных значениях параметра α и Tc = -5◦ C [Fig. 3. Graphs of temperature versus time at a depth of 0.075 m for different values of the parameter α and Tc = 5◦ C] Time corresponding to phase transitions Depth x, m 0.0015 0.0030 0.0045 Time t, sec (α = 1) 93 296 606 Time t, sec (α = 0.9) 154 545 1190 0.0060 1023 2105 Для скорости движения межфазной границы, согласно обобщенному условию Стефана, имеем следующее выражение: α VΓ = ∂0t ξ(t) = 382 σ Γ(1 - α) t 0 ξ (s) σ(α)Γ(α/2 + 1) -α/2 ds = t . α (t - s) Γ(1 - α/2) Моделирование процессов промерзания одномерным уравнением . . . Таким образом, скорость движения фазовой границы является функцией, зависящей от времени и параметра производной дробного порядка. Фазовая скорость VΓ → 0 при t → ∞. Заключение. В работе на основе классической модели Стефана построена математическая модель процессов промерзания с учетом особенностей теплопереноса на межфазовой границе, учитывающая эффекты памяти и фрактальность среды. Разработан алгоритм и создана программа численного решения задачи Стефана с оператором дробного дифференцирования. Оценены функциональные зависимости движения межфазной границы для обобщенного условия Стефана в зависимости от значения дробного параметра α. Установлено, что переход к дробным производным позволяет описать замедление процесса промерзания грунта относительно классического решения. Конкурирующие интересы. Мы не имеем конкурирующих интересов. Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.

About the authors

Vetlugin D Beybalaev

Dagestan State University; Institute of Geothermal Problems, Dagestan Scientific Center of RAS

Email: kaspij_03@mail.ru
43a, M. Gadzhiev st., Makhachkala, 367025, Russian Federation; 39a, Shamilya av., Makhachkala, 367030, Russian Federation
Cand. Phis.& Math. Sci., Associate Professor; Associate Professor, Dept. of Applied Mathematics; Senior Researcher, Lab. of Mathematical Modeling of Geothermal Objects

Abutrab A Aliverdiev

Dagestan State University; Institute of Geothermal Problems, Dagestan Scientific Center of RAS

Email: aliverdi@mail.ru
43a, M. Gadzhiev st., Makhachkala, 367025, Russian Federation; 39a, Shamilya av., Makhachkala, 367030, Russian Federation
Dr. Phis.& Math. Sci., Professor; Professor, Dept. of Applied Mathematics; Head of the Laboratory, Lab. of Mathematical Modeling of Geothermal Objects

Ramazan A Magomedov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

Email: ramazan_magomedov@rambler.ru
39a, Shamilya av., Makhachkala, 367030, Russian Federation
Senior Researcher; Lab. of Mathematical Modeling of Geothermal Objects

Rashid R Meilanov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

39a, Shamilya av., Makhachkala, 367030, Russian Federation
Junior Researcher; Lab. of Mathematical Modeling of Geothermal Objects

Enver N Akhmedov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

Email: aen-code@yandex.ru
39a, Shamilya av., Makhachkala, 367030, Russian Federation
Junior Researcher; Lab. of Mathematical Modeling of Geothermal Objects

References

  1. Liu Junyi, Xu Mingyu Some exact solutions to Stefan problems with fractional differential equations // J. Math. Anal. Appl., 2009. vol. 351, no. 2. pp. 536-542. doi: 10.1016/j.jmaa.2008.10.042.
  2. Мейланов Р. П., Бейбалаев В. Д., Шахбанова М. Р. Прикладные аспекты дробного исчисления. Saarbrücken: Palmarium Academic Publishing, 2012. 135 с.
  3. Алхасов А. Б., Мейланов Р. П., Шабанова М. Р. Уравнение теплопроводности в производных дробного порядка // ИФЖ, 2011. Т. 84, № 2. С. 309-317.
  4. Самко С. Г., Килбас А. А., Маричев О. И. Интегралы и производные дробного порядка и некоторые их приложения. Минск: Наука и техника, 1987. 688 с.
  5. Нахушев А. М. Дробное исчисление и его применение. М.: Физматлит, 2003. 272 с.
  6. Tadjeran C., Meerschaert M. M., Scheeffler H.-P. A second-order accurate numerical approximation for the fractional diffusion equation // Journal of Computational Physics, 2006. vol. 213, no. 1. pp. 205-213. doi: 10.1016/j.jcp.2005.08.008.
  7. Meerschaert M. M., Tadjeran C. Finite difference approximations for two-sided spacefractional partial differential equations // Applied Numerical Mathematics, 2006. vol. 56, no. 1. pp. 80-90. doi: 10.1016/j.apnum.2005.02.008.
  8. Лафишева M. M., Шхануков-Лафишев М. Х. Локально-одномерная разностная схема для уравнения диффузии дробного порядка // Ж. вычисл. матем. и матем. физ., 2008. Т. 48, № 10. С. 1878-1887.
  9. Алиханов А. А. Разностные методы решения краевых задач для волнового уравнения с дробной производной по времени // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2008. № 2(17). С. 13-20. doi: 10.14498/vsgtu606.
  10. Бейбалаев В. Д. Численный метод решения задачи переноса с двусторонней производной дробного порядка // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2009. № 1(18). С. 267-270. doi: 10.14498/vsgtu643.
  11. Бейбалаев В. Д. Математическая модель теплопереноса в средах с фрактальной структурой // Матем. моделирование, 2009. Т. 21, № 5. С. 55-62.
  12. Таукенова Ф. И., Шхануков-Лафишев М.Х. Разностные методы решения краевых задач для дифференциальных уравнений дробного порядка // Ж. вычисл. матем. и матем. физ., 2006. Т. 46, № 10. С. 1871-1881.
  13. Головизнин В. М., Короткин И. А. Методы численного решения некоторых одномерных уравнений с дробными производными // Дифференц. уравнения, 2006. Т. 42, № 7. С. 907-913.
  14. Бейбалаев В. Д., Абдуллаев И. А., Наврузова К. А., Гаджиева Т. Ю. О разностных методах решения задачи Коши для ОДУ с оператором дробного дифференцирования // Вестник дагестанского государственного университета. Сер. 1. Естественные науки, 2014. № 6. С. 53-61.
  15. Кузнецов Г. В., Шеремет М. А. Разностные методы решения задач теплопроводности. Томск: Томск. политехн. ун-т, 2007. 172 с.

Statistics

Views

Abstract - 28

PDF (Russian) - 5

Cited-By


Article Metrics

Metrics Loading ...

PlumX

Dimensions

Refbacks

  • There are currently no refbacks.

Copyright (c) 2017 Samara State Technical University

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies