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

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
Cand. Phis.& Math. Sci., Associate Professor; Associate Professor, Dept. of Applied Mathematics; Senior Researcher, Lab. of Mathematical Modeling of Geothermal Objects 43a, M. Gadzhiev st., Makhachkala, 367025, Russian Federation; 39a, Shamilya av., Makhachkala, 367030, Russian Federation

Abutrab A Aliverdiev

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

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

Ramazan A Magomedov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

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

Rashid R Meilanov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

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

Enver N Akhmedov

Institute of Geothermal Problems, Dagestan Scientific Center of RAS

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

References

Supplementary files