Верификация математической модели процесса промерзания с уравнением сопряжения на подвижной границе Текст научной статьи по специальности «Математика»

Journal of Siberian Federal University. Engineering & Technologies, 2016, 9(1), 24-31

УДК 624.14: 519.633

Verification of Frost Penetration Mathematical Model with Conjugate Equation in Moving Boundary

Valentin M. Juravlev, Evgeny P. Khagleev* and Pavel E. Khagleev

Siberian Federal University 79 Svobodny, Krasnoyarsk, 660041, Russia

Received 26.11.2015, received in revised form 07.12.2015, accepted 29.01.2016

The article devote to experimental verification of Stefan problem mathematical modelling with conjugate equation in moving boundary instead of the traditional forth type boundary condition. The physical model design and procedure of the experiment described. The comparison and analysis of numerical simulation and physical modelling make to confirm accordance the mathematical model with conjugate equation to the real frost penetration process.

Keywords: temperature, latent phase transition heat, moving boundary, conjugation equation, similitude parameter.

Citation: Juravlev V.M., Khagleev E.P., Khagleev P.E.Verification of frost penetration mathematical model with conjugate equation in moving boundary, J. Sib. Fed. Univ. Eng. technol., 2016, 9(1), 24-31. DOI: 10.17516/1999-494X-2016-9-1-24-31.

Верификация математической модели процесса промерзания с уравнением сопряжения на подвижной границе

В.М. Журавлев, Е.П. Хаглеев, П.Е. Хаглеев

Сибирский федеральный университет Россия, 660041, Красноярск, пр. Свободный, 79

Статья посвящена экспериментальной проверке математической модели задачи Стефана с использованием уравнения сопряжения на подвижной границе вместо традиционно применяемого граничного условия четвертого рода. Описаны конструкция физической модели, методика и ход эксперимента. Проведены сопоставление и анализ результатов численного эксперимента и физического моделирования, которые подтвердили адекватность математической модели с уравнением сопряжения реальному процессу промерзания.

© Siberian Federal University. All rights reserved Corresponding author E-mail address: sfu118@mail.ru

*

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

Введение

В классической постановке задачи Стефана [1] на подвижной границе промерзания влажного материала записывают граничное условие четвертого рода (ГУ IV) с источником теплоты или так называемое условие Стефана:

- - ч — -К — = п ^к п\

Ч - Ьи - ; К / дп -Ки дп - пРнах' (1)

где / - соответственно температура мерзлой, талой зон влажного грунта и фазового

превращения воды в лед, °С; X/, соответственно коэффициенты теплопроводности грунта в мерзлой и талой зонах, М//(тК); п - нормаль к поверхности раздела фаз; 4 - координата подвижной границы раздела фаз, м; т - время, с; Qph = qpИ р Ж - скрытая теплота фазового превращения влажного грунта, Дж/м3; qph, р, Ж - соответственно удельная скрытая теплота фазового превращения, Дж/кг, плотность грунта в сухом состоянии, кг/м3, и весовая влажность грунта в долях единицы.

Задача Стефана с ГУ IV рода (1) относится к классу сопряженных задач теплообмена в гетерогенных системах с подвижными и неподвижными границами раздела фаз. Однако существует иной подход 1С описанию задачи Стефана, суть которого заключается в одинаковом представлении процессов переноса теплоты в сопряженных телах и на границах их раздела. К подобному заключению пришли авторы работ [2-4], объединив уравнение теплопроводности с ГУ IV рода (1). Полученное таким образом обобщенное уравнение теплопроводности с эффективной сглаженной теплоемкостью становится пригодным для описания всей расчетной обла-стис. Это позволило в сопряженной задаче теплообмена с фазовыми превращениями вещества на подвижной границе применить метод сквозного счета.

Но к этому заключенню о едином математическом описании гетерогенной системы можно прийти прямым путем, непосредственно наблюдая природу процессов переноса теплоты и фазовых превращений вещества. Действительно, теплофизические свойства веществ на границах раздела и внутри взаимодействующих тел по природе своей одинаковы, отличаются друг от друга лишь количественно. Механизмы переноса теплоты в этих телах и на их границах раздела также происходят по одним и тем же законам природы: теплопроводность, конвекция и излучение. Подобное единство свойств веществ и механизмов переноса теплоты предполагает одинаковое формализованное описание процессов теплообмена в виде уравнений энергии, в частном случае теплопроводности [5-7]. При этом в уравнении теплопроводности на подвижной границе, назовем его уравнением сопряжения, для влажного вещества (Ж > 0) появится источник теплоты в виде фазовых превращений влаги.

На первом этапе проверки адекватности математической модели в новой постановке реальному процессу нами был проведен вычислительный эксперимент по изучению процесса промерзания влажного грунта [8]. Для сравнения данная задача была реализована в традиционной постановке с ГУ IV (1) с использованием известных численных методов решения: метода ловли фронта в узел сетки (МЛФУС) и метода выпрямления фронтов (МВФ) [4, 9-13]. Срав-

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

Необходимо отметить, что численное решение задачи с уравнением сопряжения на подвижной границе [8] обладает рядом достоинств перед известными методами. Во-первых, модель с уравнением сопряжения позволяет реализовать метод сквозного счета, как в работах [2-4], но в отличие от них в авторской модели представляется возможность явно выделять местоположение подвижной границы. Во-вторых, местоположение подвижной границы устанавливается без перестройки пространственно-временной сетки на каждом новом шаге по времени в отличие от методов МЛФУС и МВФ [9-13]. И, в-третьих, фазовое превращение влаги на подвижной границе рассматривается как процесс с непрерывным отслеживанием агрегатного состава вода-лед от начального до его конечного значения W0 > W > 0 [8], где W0, W - соответственно начальная и текущая весовая влажность материала в долях единицы. Тогда как в известных методах фазовое превращение влаги на подвижной границе «проглатывается» за один шаг по времени с фиксированием только крайних значений влажности: начальной W = W0 и конечной W = 0 [2-4, 9-13].

Для дальнейшей проверки адекватности математической модели в новой постановке наряду с проведенным вычислительным экспериментом [8] мы сравнили результаты численного решения задачи с результатами наблюдений за процессом промерзания влажного материала на физической модели.

Физическая модель процесса промерзания слоя влажного материала

В качестве испытуемого материала был использован шлак жидкого шлакоудаления с частицами 2 - 4 мм, полученный при сжигании бурого угля Ирша-Бородинского месторождения со следующими теплофизическими свойствами: W0 = 0,065; р, = 1350 кг/м3; = 0,25 Вт/(м-К); Cth = 1090 кДж/(м3К); р/ = 1348 кг/м3; 1/ = 0,27 Вт/(м-К); СЛ = 988 кДж/(м3-К).

Физическое моделирование процесса промерзания влажного шлака проведено на экспериментальной установке (рис. 1), воспроизводящей климатические условия г. Красноярска:

• продолжительность среднемноголетнего зимнего сезона составляет пять месяцев: с ноября по март включительно;

• среднезимняя температура наружного воздуха = минус 11 °С;

• сезонное промерзание грунтов в г. Красноярске составляет ^ = 2 м.

Геометрический масштаб моделирования в экспериментальной установке принят равным

Г = 1/10. При этом толщина слоя шлака в установке с учетом величины сезонного промерзания грунтов составила И3 = Г = 20 [2].

Для обеспечения временного подобия пр оцессов те плообмена в модели и реальном шлакоотвале использован критерий! тепловой гомохронности Фурье:

Ро=# ,2)

где а - коэффициент температуропроводности шлака, м2/с.

6

Рис. 1. Экспериментальная установка: 1 - корпус установки с перфорированным дном; 2 - слой шлака; 3 - датчики температуры; 4 - емкость с антифризом; 5 - крышка емкости с антифризом; 6 - кабель датчиков температуры; 7 - внешняя теплоизоляция; 8 - водонасыщенный дренаж; 9 - емкость с крупнообломочным материалом, залитым водой

Согласно этому критерию (2), учитывая, что а' = а, соотношение между временем развития процесса в модели и реальном шлакоотвале составит

г \2

' к л

т =

s.f

т = 0,01 • т (3)

Про тяженность зимнего периода в г. Красноярске, как отмечалось выше, составляет пять месяцев, или 3625 ч (725 ч в месяце). Согласно соотношению (3) протяженность «зимнего периода» в эксперименте на физической модели составит 36,25 ч (7,25 ч в месяце).

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

В сосуд плотно укладывался влажный шлак высотой И' = 0,20 м, а под перфорированным днищем сосуда слой дренажа толщиной И = 0,01 м. В процессе укладки шлака устанавливали десять цифровых термометров DS18В20 с шагом 0,02 м по высоте, а в приповерхностном сло е - два термометра с шагом 0,01 м (рис. 1). Контактные выводы термометров витой парой присоединялись к интерфейсной плате персонального компьютера.

Сосуд после засыпки дренажа и шлака был установлен в емкость диаметром 270 мм с гравием на дне, залитым водой до контакта шлак-дренаж. Установка в собранном виде помещалась в холодильную камеру.

Граничное условие I рода на верхней поверхности слоя шлака воспроизводилось путем размещения на этой поверхности тонкостенного алюминиевого цилиндрического сосуда с залитым в него антифризом «IceStream - 45°», охлажденным до минус 12 °С (рис. 1). Когда антифриз в ходе эксперимента нагревался до минус 7-9 °С, его замещали новой порцией охлажденного до минус 12 °С антифриза и т. д. В результате на верхней поверхности шлака воспроизводилось ГУ I рода в виде «пилообразной функции» (рис. 2, кривая 1).

-и

JJ It It >0 J6

Рис. 2. Температура на внешних поверхностях слоя шлака (ГУ I рода): 1 - температура на верхней поверхности слоя шлака t (0, т); 2 - температура на донной поверхности слоя t (20, т) (шлак-дренаж); 3 -ередняя по интервалам времени темпераарре наружной поверхности

Математическая модель иромерзания влажного шлака с использованием уравнения сопряжения

Математическую модель вышеописанной физической постановки задачи можно представить в виде следующих дифференциальных уравнений:

• для внутренних точек мерзлого и талого слоев влажного шлака уравнения теплопроводности вида

С да <з2? _ a d2t

С/ : / л. 2 ' Cth - - Pth „ 2 '

я дт J ду 22т (9y2

• уравнение сопряжения на подвижной границе

(4)

0 -

( дН ^

(т)+0

( д

- qphP

д (т)-о

dW

дт

(5)

уравнение движения для подвижнри границы! промерзания

«О =

W0 > W > 0 ;

кЧтА

I в (т) ± ¿В.- W - 0-

(6)

1В уравнении сопряжения (5), представляющем собой уравнеаве оеплопроводнести на подвижной границе, отсутсттует локальная производная, характеризующая изменение температуры во времени ^/дт). Вода во влажном шлаке находится в свободном, а не в связанном состоянии, иак, например, в глиниитых гру нтах. Поэтому тммпература фазового превращения воды в лед на подвижной границе - величина пастаяаная, отсюда ^/дт) = 0, и уравнение сопряжения на подвижной границе принимает вид (5).

Для замыкания системы дифференциальных уравнений (4) - (6) были приняты следующие краевые условия:

• начальное ус ловие

(7)

t(0, т) = fi(x);

(8)

t (hs,x) = /2(т).

(9)

Пилообразная фу нкция температуры на верхней (8) и донной (9) поверхностях слоя шлака (рис. 2, кривые 1, 2) в алгоритме задачи были представлены кусочно-линейными функциями вида

где t0n - температура на поверхности слоя в начале n-го интервала времени, °С; кп - угловой коэффициент на n-м интервеле; xn+i - конечное значение времени n-го интервала, с; т,- текущее оремя, с.

Таким образом, процесс теплообмена во всей гетерогенной системе с подвижной границей описывается с помощью уравнений единой физической природы - уравнений теплопроводности вида (4, 5). Единообразие математической модели (4)-(9) позволило применить метод сквозного счета без перестройки пространственно-временной сетки (Ду = const, Дт = const) с явным выделением подвижной границы промерзания [!8].

В вычислительном эксперименте была испо льзоваеа достаточно подробная пространственно-временная сетку. Шаг по пространственной координате принят равным Ду = 0,002 м, чта для слоя шлока толщиней 0,220 м (рис. 1) соответствовало ста узловым точкам. Валичина шигт по времени выбрана наименьшейиз двух возможныв значений для мерзлой и талой зон слоя:

где Дт < Ду2/2а,, , = 1, 2 - условие устойчивос ти для мерзлого г = 1 и талого , = 2 слоев шлака; ак = X,- / С,- - коэффициент температур опроводност и мерзлой или талой зон слоя шлака, м2/с. Таким образом, шаг временной сетки согласно (10) составил Дт = 6,19 с.

Анализируя результаты математического и физического моделирования (рис. 3), можно отметить, что матуматическая модель с уравнением сопряжения на подвижной границе достаточно точно, с расхождением 6 - 9 %, описывают общий ход процпс са промерзания.

Пенрвоначально в период от 0 до 8,75 ч наблюдалось быстрое промерзание влажного шлака (рис. а). Это объясняется высокими градиентами температуры, когда на верхней поверхности слоя шлака установилась температура в среднем на уровне минус 9,5 °С (рис. 2, кривая 3) при относительно тонком промерзающем слое шлака.

Затем с 9 до 21 ч температура на поверхности слоя шлака стала несколько выше, чем в первоначальный период, и установилась на уровне минус 9,0 °С (рис. 2, кривая 3) при продолжающемся нарастании £(т). При этом градиенты температуры в промерзающем слое снизились.

ДТ) = ЛИ + kn(Tп+1 - Т i).

Дт = тт{Дт;},

(10)

Обсуждение результатов эксперимента

л л 17 m т j 1» in.

Рис. 3. Вычислительный эксперимент по промерзанию слоя влажного шлака: 1 - кривая промерзания на физической модели; 2 - кривая промерзания в математической модели с уравнением сопряжения на подвижной границе

В результате процесс промерзания замедлился, кривые промерзания 1 и 2 стали более пологими (рис. 3).

На следующей стадии с 21 до 32 ч температура на поверхности шлака понизилась и установилась на уровне минус 10,5 °С. В результате темп промерзания слоя влажного шлака сохранился. Зотем с 32 ч до конца эксперимеера в 36,25 ч температ°ра на поверхности шлака снова повысилась до минус 9,5 °С (рис. 2, кривая 3), что привело к снижению интенсивности промерзания (рис. 3).

Заключение

Реал изация ороемотиче ской модраи с использеаанием усавнероо српряжерия нр подвижной границе дает достаточно хорошее совпадение с результатами физического эксперимента. Предложенный метод решения задачи Стефана, прошедший экспериментальную проверку методом физического моделирования и ранее выполненной проверкой известными методами численного дешения [8], может быть применен в практических расчетах при изучении процессов теплообмена в гетерогенных системах, сопровождаемых фазовыми превращениями вещества.

1 СНиП 2.02.01-83* «Основания зданий и сооружений».

2 Физические величины, принятые; соштрихами, относятся к модели.

Список литературы

[1] Основы мерзлотного прогноза при инженерно-геологических исследованиях / под ред. В. А. Кудрявцева. М.: Изд-во МГУ, 1974, 431 с. [The basics of prognosis of permafrost engineering geological studies. Ed. By V.A.Kudriavtsev. Moscow State University, 1974, 431 p. (in Russian)]

[2] Самарский А. А., Моисеенко Б. Д. Журн. вычисл. математики и мат. физики, 1965. 5(5), 816-827 [Samarskii A.A., Moiseenko B.D. Computational Mathematics and Mathematical Physics, 1965. 5(5), 816-827 (in Russian)]

[3] Будак Б. М., Соловьева Е. Н., Успенский А. Б. Журн. вычисл. математики и мат. физики, 1965, 5(5), 828-840 [Budak B.M., Solov'eva E.N., Uspenskii A.B. Computational Mathematics and Mathematical Physics, 1965, 5(5), 828-840 (in Russian)]

[4] Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с. [Samarskii A.A., Vabishchevuch P.N. Computational heat transfer, Moscow, Editorial URSS, 2003. 784 p. (in Russian)]

[5] Хаглеев, Е. П., Хаглеев Ф. П. Гидродинамика больших скоростей: межвуз. сб. КрПИ. Красноярск, 1992, 63-69 [Khagleev E.P., Khagleev P.F. Hydrodynamics of high speeds, Krasnoyarsk, 1992, 63-69 (in Russian)]

[6] Хаглеев Е. П., Хаглеев П. Е. Труды II Российской национальной конференции по теплообмену. Т. 6. Интенсификация теплообмена. Радиационный и сложный теплообмен. М.: Изд-во МЭИ, 1998, 216-219 [Khagleev E.P., Khagleev P.F. Proceedings of II Russian national conference on heat transfer, Vol. 6, Moscow, National Research University "MPEI", 1998, 216-219 (in Russian)]

[7] Khagleev E. P. J. Sib. Fed. Univ. Eng. technol, 2013, 5(6), 485-497.

[8] Khagleev E. P. J. Sib. Fed. Univ. Eng. technol, 2014, 7 (8), 894-910.

[9] Будак Б. М., Гольдман Н. Л., Успенский А. Б. Решения задач типа Стефана. М.: Изд-во МГУ, 1972, вып. 2, 3-23 [Budak B.M., Goldman N.L., Uspenskii A.B. Solving Stefan-type problems, Moscow State University, 1972, vol. 2, 3-23 (in Russian)]

[10] Crank J. Free and Moving Boundary Problems. Oxford, Clarendon Press, 1984, 425 p.

[11] Douglas J, Gallie G. M. On the numerical integration of a parabolic differential equation subject to a moving boundary condition. Duke Math. J., 1955, 22(4), 557-572.

[12] Javierre-Perez E. Literature Study: Numerical methods for solving Stefan problems Delft University of Technology, Report 03-16, 2003, 94 p.

[13] Красношлык Н. А., Богатырев А. О. Вкник Черкаського унiверситету. Серiя При-кладна математика. Ыформатика, 2011, 194, 11-31 [Krasnoshlyk N.A., Bogatyrev A.O. Herald of Cherkasy University. Applied mathematics, computer science, 2011, 194, 11-31 (in Ukrainian)]