РЕШЕНИЕ ДВУХФАЗНОЙ ЗАДАЧИ СТЕФАНА В ЭНТАЛЬПИЙНОЙ ПОСТАНОВКЕ СО СГЛАЖИВАНИЕМ КОЭФФИЦИЕНТОВ Текст научной статьи по специальности «Физика»

УДК 519.633

DOI: 10.18698/1812-3368-2021-4-4-23

РЕШЕНИЕ ДВУХФАЗНОЙ ЗАДАЧИ СТЕФАНА В ЭНТАЛЬПИЙНОЙ ПОСТАНОВКЕ СО СГЛАЖИВАНИЕМ КОЭФФИЦИЕНТОВ

B.И. Васильев1 М.В. Васильева1, 2

C.П. Степанов1 Н.И. Сидняев3 О.И. Матвеева4 А.Н. Цеева4

vasvasil@mail.ru vasilyevadotmdotv@gmail.com cepe2a@inbox.ru sidnyaev@bmstu.ru

yapniis@mail.ru

antseeva@mail.ru

1 СВФУ, Якутск, Российская Федерация

2 Техасский университет A&M, Колледж-Стейшен, Техас, США

3 МГТУ им. Н.Э. Баумана, Москва, Российская Федерация

4 ЯкутПНИИС, Якутск, Российская Федерация

Аннотация

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

Ключевые слова

Математическое моделирование, криолитозона, энтальпия, задача Стефана, сглаживание коэффициентов, скрытая теплота, численные расчеты

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

с использованием цементно-песчаных растворов Поступила 25.02.2021 с температурой 25 °С. Получены распределения тем- Принята 02.04.2021 пературы грунта для различных моментов времени © Автор(ы), 2021

Работа выполнена при финансовой поддержке РНФ

(проект РНФ № 19-11-00230)

Введение. Тепловые расчеты имеют важное значение в строительстве инженерных геотехнических сооружений и зданий в криолитозоне. Температурный режим (совокупность последовательных температурных полей в грунтовом массиве, соответствующих любым заданным моментам времени от начала расчета) рассчитывается как результат задаваемых на весь период расчета прогноза тепловых воздействий на верхней и нижней границах основания сооружения. Наиболее характерной особенностью этих процессов, вследствие которых их математические модели нелинейны и трудны для анализа, являются неизвестные заранее («свободные») границы между различными (талым и мерзлым) состояниями грунта. Результаты расчетов климатических явлений, имеющие важное значение для геотехнических сооружений на многолетнемерзлых основаниях (в криолитозоне), представляют собой непрерывное повторение одного и того же неизменного сезонного цикла. Этим обусловливаются периодические гармонические колебания всех величин, определяющих состояние основания геотехнического сооружения, в том числе и температуры. Однако и при других процессах также происходят периодические гармонические изменения состояния, вызванные постоянно повторяющимися воздействиями. Изменение температуры за время одного периода может иметь различный характер, например изменяться скачкообразно, непрерывно возрастать или убывать и т. д.

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

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

В научной литературе имеется несколько названий задачи Стефана: задача с подвижными границами (moving boundary problem), задача со свободными границами (free boundary problem), задача о фазовом переходе (phase change problem). Примеры физических процессов с фазовыми переходами — задача о таянии льда с подвижной границей между водой и льдом, задача о плавлении твердого вещества с неизвестной границей между твердой и жидкой фазами.

Принципиально новый подход, дающий возможность эффективно численно решать многомерные, многофронтовые задачи типа Стефана, в феврале 1964 г. на Втором всесоюзном съезде по теоретической и прикладной механике независимо предложили две группы отечественных математиков под руководством А.А. Самарского и Б.М. Будака, впоследствии названный энтальпийным методом [1, 2]. Основная идея энтальпийного метода заключается во введении понятия «эффективной» теплоемкости, содержащей скрытую теплоту фазового перехода, которая выделяется (поглощается) на поверхности раздела фаз. Такая постановка с использованием дельта-функции Дирака позволяет записать единое квазилинейное параболическое уравнение сохранения энергии сразу во всей области, при этом условие Стефана автоматически выполняется на всех поверхностях раздела фаз. Для численного решения полученной задачи ими предложены различные способы размазывания функции источника и полиномиаль-

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

Обоснование корректности и численные методы решения различных актуальных прикладных задач Стефана обобщены в [3-6]. Для численного решения задачи Стефана в [7] предложен метод контрольного объема, применяемый для численного решения трехмерных задач теплопроводности с учетом фазовых переходов. Математическая модель пониженного порядка процессов теплопередачи грунтового теплообменника в многолетнемерзлые основания инженерных сооружений, значительно уменьшающая вычислительную работу, построена и численно реализована в [8]. Там же представлены численные результаты решения модельной задачи, иллюстрирующие эффективность предложенной математической модели и численного алгоритма процессов теплопередачи в грунте через систему охлаждающих труб. С помощью обобщенного многомасштабного разрывного метода Галеркина в [9] численно решена задача теплообмена с фазовым переходом в гетерогенных областях. Рассмотрены численные результаты для различных геометрий, демонстрирующие хорошую точность метода. Для определения напряженно-деформированного состояния мерзлых грунтов в [10] использовано численное моделирование методом конечных элементов задач тер-мопороупругопластичности с фазовыми переходами. Особое внимание уделено оттаиванию мерзлых грунтов, которое может вызвать дополнительные деформации основания инженерного сооружения и привести к потере его устойчивости. С использованием обобщенного метода конечных элементов в [11, 12] построена и численно реализована математическая модель процесса искусственного замораживания оснований зданий в криолитозоне, представлены численные результаты для двух модельных задач в двух-и трехмерной постановках. Разработке и численной реализации методом конечных элементов математической модели теплообмена отапливаемого сооружения в криолитозоне со слоистым основанием посвящена [13]. Три метода численного решения двухфазной задачи Стефана в энтальпий-ной постановке, основанные на различных способах сглаживания разрывных коэффициентов и размазывания функции источника, предложены в [14]. Приведены результаты сравнения численных расчетов модельной задачи с точным решением для одномерной постановки, а также расширение метода для решения двухмерной задачи.

Энтальпийная формулировка двухфазной задачи Стефана. Для

простоты изложения рассмотрим двухфазную задачу Стефана для одномерного квазилинейного параболического уравнения. Пусть определению подлежит функция иt), х еО. = [0, /], £ е(0, Т], Т > 0, удовлетворяющая квазилинейному дифференциальному уравнению с частными производными второго порядка параболического типа с разрывными коэффициентами:

(с (и) + £5 (и - и")= 4- {к (и)Iй| + /(и, х, £),

ох V дх) (1)

хеО, £е(0,Т].

Коэффициенты уравнения (1) зависят от искомой функции и £), терпят разрыв первого рода при и £) = и*, следовательно, задаются формулами

Г С1(и), и < и*, Г кх(и), и < и*,

с(и)= N ^ * к(и) = 1 N ^ *

[ с2(и), и > и, [ к2(и), и > и,

где 5 (•) — дельта-функция Дирака; и* — заданная постоянная величина, представляющая собой значение искомой функции и £) на границе раздела фаз (температура фазового перехода). Такая постановка задачи работает и для многомерных задач и в случае наличия нескольких границ раздела фаз, что часто встречается при решении актуальных задач криолито-зоны. Поскольку коэффициент левой части уравнения (1) при и £ ) = и*

обращается в бесконечность, численное решение данного уравнения возможно после размазывания функции точечного источника и сглаживания его коэффициентов по искомой функции и £). Отметим, что в [1, 2] предложены различные способы такого сглаживания и размазывания.

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

В качестве коэффициентов с использованием специальных функций зададим гладкие функции сд(и), кд(и), 5д(и) с одним параметром А :

С2 - С1

(

cA(u) = С\ 1 + erf + D5a("); (2)

kA (u) = k + ^ {1 + erf , (3)

V

где

§д (и)= е"(и - и*)2/(2 А2). (4)

у/2кА

Следовательно, при таком подходе фазовые превращения не проходят мгновенно, а происходят в малом интервале температуры (и* -А, и* + А). Разрывные коэффициенты уравнения (1) заменяются ограниченными гладкими функциями температуры, в результате получается квазилинейное параболическое уравнение:

(сд (и) + Обд (и - и*) )4Ы = кА (и) + f(u, х, г), х еО, г е (0, Т]. (5)

дг ох V дх у

Для того чтобы определить единственное решение уравнения (2), зададим граничные условия, например, граничное условие третьего рода

Ои

-а,1(и) -Р1&Д (и) —= gl(u, г), х = 0, г е (0, Т ]; (6)

ох

Ои

а2(и) + Р2^д(и)— = £2(и,г), х = I, г е (0,Т], (7)

ох

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

и(х ,0) = и0(х), х еО. (8)

Здесь аг(и), gi(х, г), г = 1,2 — заданные функции; рг, г =1,2 — постоянная величина, принимающая значения 0 или 1. Таким образом, граничные условия (6), (7) в зависимости от аг-(и), представляют собой взвешенную комбинацию граничных условий первого-третьего родов и условия Стефана — Больцмана. Если задано условие первого рода, то рг =0, аг (и) = 1, gi (и, г) = gi (г), г =1,2, в остальных случаях рг =1, г = 1,2.

Отметим, что (2) является квазилинейным параболическим уравнением с положительными ограниченными гладкими коэффициентами.

Неявная разностная схема в случае гладких коэффициентов. Аппроксимация. Для численного решения одномерной задачи Стефана (5)-(8) используем метод конечных разностей [14]. В области определения задачи введем равномерную пространственно-временную сетку

®Ит = ®И х®т,

®И ={хг = гИ, г = 0,1,..,п; И = 1 /п}; шт ={г^ = ;т, ] =0,1,..,/; т = Т/]}.

Начально-краевой задаче (5)-(8) с гладкими коэффициентами (2)-(4) на построенной сетке ®нт поставим в соответствие линеаризованную неявную разностную схему, построенную интегро-интерполяционным методом [14]. На каждом временном слое во внутренних узлах пространственной сетки имеем конечно-разностный аналог дифференциального уравнения (5) в виде системы трехточечных нелинейных уравнений:

3 - 1 f . - . - Л

ei yi -yi _ 1

т h

ai yui - yi jy{ - yj-i aM— a —Г"

+fi3,

i _1,2,.. ,n -1, 3 _1,2,.., /. (9)

Здесь

ej _ cA(y{ -1), a/+i_ кд ^ 3 + 3 j,

d3 _ f (yi_1,Xi, t/ ), i _ 1,2,...,n -1, / _1,..., J.

Выпишем конечно-разностные аналоги обобщенных граничных условий (6):

eiyLzÄL_2a/ybyLJ^O+ä+g/ приР_1,/_ 1,2,—,j;

т h2 h h

y0 _g1(t/) при ß _ 0, / _1,2,..., J.

(10)

Аналогично строим дискретный аналог граничного условия на правой границе (7):

/yn - yn-1 2a//LZÄ а2 (уП ) . г/ . g2 (yn, ti)

n _ 2W n _ „ _ l/n "T ,

[yi _ g2(t/) при ß _0, / _ 1,2,..., J.

e^Z^^ _2anin=L^n--^ + fj + ~ v -, приß _1,/ _1,2,... ,J;

x h2 h h

(11)

Полученную систему замыкаем дискретным аналогом начального условия (8):

уI = щ(х1), I = 0,1,..., п. (12)

Погрешность аппроксимации разностной схемы (9)—(12) при решении дифференциальной задачи (5)-(8) имеет порядок 0(х + Н2). Таким

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

Здесь выбор значения параметра А на каждом временном слое предлагается автоматизировать с использованием следующей процедуры.

1. Определение пространственного интервала, содержащего границу

раздела фаз или, другими словами, нахождение i* такого, что yj*yj*+1 < 0.

2. Принятие А =| yj*+1 - yj*_ 11.

Численные расчеты. Приведем результаты сравнения предложенной разностной схемы на примере решения двухфазной автомодельной задачи Стефана с точным решением [15] при входных данных: l = 8 м; T = 107 с;

Pill

u(0, t) = -5 °С; —(l,t) = 0, u0(x) = u0 = 5 °С; u =0; k1=2,21 Вт/(м-К);

dx

k2=0,59 Вт/(м-К); d = 1,89-106 Дж/(м3-К); c2=4,12-106 Дж/(м3-К);

D = 3,33-108 Дж/(м3-К); E, = yVt. Здесь у — корень трансцендентного уравнения:

k1 gg-T2/(4«1)2 + k2UQg-^2/(4a2)2 + d к =0

оцФ (у/2а1) оц(1 -Ф (у/2а0) 2 ,

где a1 = 4k 1 / С1; й2 = 4h / С2; Ф(-) — интеграл ошибок. Вычисленное с использованием итерационного метода секущих значение равно у = = 0,00023897230346. Для получения точности |у"5+1 -уs |<10"12 понадобилось провести девять итераций, т. е. итерации сходились достаточно быстро.

Графики автомодельного (точного) решения в финальный момент времени T = 115,741 сут и границы раздела фаз при t е [0,115,741] сут приведены на рис. 1. На левом графике четко заметно, что первая производная искомой функции по пространственной переменной x на границе раздела фаз u (^(T), T) = u претерпевает разрыв первого рода в соответствии с условием Стефана.

Проведем численную реализацию линеаризованной разностной схемы (9)-(12), построенной для численного решения задачи с гладкими коэффициентами (5)-(8) с автоматическим выбором интервала сглаживания А =| yj*+\ - yj*_11 по температуре, где y^y^ < 0.

и(х, Т) \

4

2

О -2 -4 -6

Рис. 1. Графики автомодельного решения (а) и границы раздела фаз (б)

Результаты расчетов с использованием линеаризованной неявной разностной схемы (9)-(12) приведены на рис. 2: на рис. 2, а — графики численного (кривая синего цвета) и точного (кривая черного цвета) распределений температуры в финальный момент времени, на рис. 2, б — графики численно найденной методом линейного интерполирования в том интервале, в концах которого искомая функция имеет разные знаки у\ 1 у/^1 < 0 (кривая синего цвета), и точной границы раздела фаз (кривая черного цвета). Расчеты проведены на достаточно грубой пространственно-временной сетке п = 100, ] = 100. Непосредственно за границей раздела фаз наблюдается несущественный «перегрев» в малой области, занятой водой. Граница раздела фаз определяется с небольшой погрешностью, за исключением малых значений времени.

Г %

4

2

0 -2 -4

-6

Рис. 2. Графики решения при использовании чисто неявной разностной схемы и точного решения (а), вычисленной и точной границы раздела фаз при п = 100, ] = 100 (б)

Аналогичные результаты численных расчетов, проведенных с использованием линеаризованной разностной схемы (9)-(12) при п = 200, J = 100, приведены на рис. 3. Отклонение численного решения от точного чуть меньше, чем в предыдущем случае более грубой пространственной сетки. Погрешность определения границы раздела фаз достаточно малая даже при малом времени. Таким образом, для серийных расчетов достаточно ограничиться использованием линеаризованной разностной схемы на достаточно подробной пространственной сетке.

1 2 3 4 5 6 7 X О 20 40 60 80 100 £ а б

Рис. 3. Графики решения при использовании линеаризованной разностной схемы и точного решения (а), вычисленной и точной границы раздела фаз

при п = 200, ] = 100 (б)

Численное исследование температурного режима основания сооружения после установки новой сваи. Постановка задачи. Рассмотрим математическую модель, описывающую динамику распределения температуры с учетом фазовых переходов поровой влаги в лед и обратно при некоторой заданной температуре фазового перехода Т *в области О = . Здесь (£) = {х | х еО, Т(х, t)> Т*} — область, занятая талым грунтом, где температура превышает температуру фазового перехода; О" ^) = = {х | х еО, Т(х, t)< Т*} — область, занятая мерзлым грунтом. Фазовый переход происходит на границе раздела талого и мерзлого грунтов 5 = ).

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

дТ

(Ф (Т) + тР1Ь5 (Т - Т*)) — - аММТ)Т) = /, х еО, I е (0, ии], (13)

где Ь — удельная теплота фазового перехода; 5(Т- Т*) — дельта-функция Дирака.

Для коэффициентов уравнения имеем следующие соотношения:

cp(T ) =

c р , T < T c+Р+, T > T \

UT ) =

Х~, T < T ; , T > T

Здесь с+, р+, с , р — удельная теплоемкость и плотность талого и мерзлого грунтов.

Рассматривается процесс распространения теплоты в насыщенной пористой среде, поэтому для теплофизических характеристик запишем

с "р" =(1 - т)с5ср5с + тсф1, с+р+ =(1 - ш)с5ср5с + тсн,рн,,

где т — пористость; ее, V, г — индексы, обозначающие скелет пористой

среды, воду и лед соответственно. Для коэффициентов теплопроводности талого и мерзлого грунтов запишем аналогичные соотношения:

Х~ = (1- т)Х¡с + тк[, Х+ =(1- т)Хс + тк„.

На практике фазовые превращения не проходят мгновенно и могут происходить в малом интервале температуры [Т* - А, Т* + А]. Разрывные коэффициенты уравнения (13) заменим достаточно гладкими функциями температуры:

(cp)a(T) = c р + m(c+р+- c р )

(

1 + erf

T - T

42,

42

_(T - T *)2/(2ct2).

KO^

XG (T ) = + m(X+-X~ )1

2

(

1 + erf

( t _ T * ^

42Ö2 j

Таким образом, получили уравнение для температуры:

3T

(ср)ст(T) — - div^a(T)T) = 0, X eQ, t e (0, tmax].

of

(14)

Уравнение (14) является многомерным квазилинейным параболическим уравнением с гладкими коэффициентами.

Дополним уравнение (14) начальным условием Т(х,0) = То(х), х еО. Граничные условия формулируются следующим образом.

1

e

1. На верхней (дневной) поверхности грунта задается граничное условие, учитывающее температуру воздуха региона (граничное условие третьего рода):

= а(Т - Таг), x еГ(,

СП

где а — коэффициент конвективного теплообмена; Тагг — температура наружного воздуха.

2. На боковых и нижних границах задаются условия отсутствия теплообмена:

1 дТ н г

-к— = 0, x еГ5.

дн

Аппроксимация методом конечных элементов. Квазилинейное параболическое уравнение (14) с соответствующими граничными и начальными условиями аппроксимируется с использованием метода конечных элементов в сочетании с чисто неявной линеаризованной конечно-разностной аппроксимацией по времени [14]. Это означает, что при дискретизации коэффициенты, зависящие от искомой функции, берутся с предыдущего временного слоя. Запишем вариационную постановку задачи для каждого временного слоя: найти Т е Н1 такую, что

а (Тн+1, у) = Ц(у), уеН1,

где

а(Тн+1, у) = - | (ф)а (Тп )Tn+1wdx + { V (Тн )Tn+1wdx + { аТн+1 ds;

1 а а т{

Ц(у) = -1 (ср)аTnwdx +1 аТагт ds. ТП г.

Для численного решения задачи необходимо перейти от непрерывной вариационной задачи к дискретной. Для этого введем конечномерные пространства У/ е Н1, У/ е НО и определим на них следующую задачу: найти Т/ е Ун такую, что а (ТП+1, у) = Ц(у), V е Ун, где

аОГ1, у) = -1 (ф)а (Тп )T¡n+1wнdx (Тп)Tнn+1wнdx +{ аTn+1ds;

1 а а т{

Цу) = - I (ф)ст Т^/^ + | аТаг^. ТП г.

Здесь используем линейные базисные функции.

Численные результаты. Проведем расчет влияния нового строящегося объекта на существующие строительные конструкции — сваи существующего сооружения, попадающие в зону строительства нового объекта. Таким образом численно смоделируем растепляющий эффект строящегося объекта на существующее свайное поле.

Условия воздействия нового строительства: цементно-песчаный раствор для заполнения пазух между грунтом и сваей с температурой 25 °С; сваи с поперечным сечением 40x40 см, заглубление свай на 10 м от поверхности грунта.

Рассмотрим задачу установки новых свай на расстоянии 3 м от свай существующего объекта (рис. 4). Расчетная область состоит из нескольких слоев грунта, в ней имеются три сваи существующего здания и устанавливаются две сваи строящегося объекта. Расчетные значения температуры грунта приведены в табл. 1 (результаты измерений температуры грунта перед началом укладки новых свай). Теплофизические характеристики взяты из СП 25.13330.2012 и представлены в табл. 2. Расчеты проводились на длительность 50 сут с шагом по времени т = 1 сут (24 ч).

Рис. 4. Расчетная геометрия:

1,2 — сваи строящегося и существующего объектов

Таблица 1

Расчетные значения температуры грунта

Глубина, м Температура грунта, °С Глубина, м Температура грунта, °С

1,0 -0,7 7,0 -2,9

2,0 -1,9 8,0 -3,0

Окончание табл. 1

Глубина, м Температура грунта, °С Глубина, м Температура грунта, °С

3,0 -2,4 9,0 -2,9

4,0 -2,8 10,0 -2,7

5,0 -3,0 11,0 -2,5

6,0 -3,1 12,0-15,0 -2,5

Таблица 2

Расчетные значения теплофизических характеристик грунтов

Объемная

Мощ- теплоемкость Теплопроводность Удельная

Элемент ность грунта (ср) -10~6, грунта k, Вт/(м-К) теплоемкость

грунта, м Дж/(м3-К) l -10-3, Дж/м3

талого мерзлого талого мерзлого

Щебень 0,1 2,22 2,09 1,86 2,16 22 600

Песок:

средней

крупности 1,0 3,09 2,11 1,92 2,15 101 600

то же,

мерзлый 13,9 2,97 2,23 2,28 2,52 71 600

Температура на дневной поверхности основания сооружения задавалась с учетом амплитуды колебания температуры воздуха, принятой по данным метеостанции Якутска. Зависимость температуры воздуха от времени (в течение 2 лет) показана на рис. 5.

Рис. 5. Зависимость температуры воздуха Tair от времени

Приведем результаты численных расчетов при установке фундамента в летнее время с использованием цементно-песчаных растворов при температуре 25 °С. Распределения температуры грунта для различных моментов времени приведены на рис. 6. Установка дополнительных свай при строительстве нового объекта не оказывает негативного влияния на температурный режим грунта свайного основания здания. Остывание раствора в пазухах между стволом сваи и стенками буровой скважины нового объекта и восстановление температурного режима грунтов в основании нового объекта происходит через 10 сут (см. рис. 6). Кроме того, через 15 сут температура вокруг новой сваи становится такой же, как и на глубине 3 м, происходит полное восстановление температурного режима грунта.

Рис. 6. Распределения температуры грунтов в расчетной области в начальный момент времени (а), через 10 (б), 15 (в) и 50 сут (г) (представлены распределения температуры на глубине 5 м, линия синего цвета — температура по середине сваи, красного — между двух свай)

Заключение. Для численной реализации энтальпийной постановки задачи Стефана, заключающейся во введении «эффективной» теплоемкости, содержащей скрытую теплоту фазового перехода, которая выделяется (поглощается) на поверхности раздела фаз, предложено размазывание функции точечного источника и разрывных коэффициентов гладкими функциями.

В одномерном случае построена однородная неявная разностная схема сквозного счета с автоматическим выбором параметра сглаживания. Сравнение с автомодельным решением показало хорошую точность предложенного численного метода. Метод также реализован для решения трехмерной задачи Стефана методом конечных элементов. Сравнение результатов расчетов, выполненных с использованием предложенного численного метода решения трехмерной задачи, с данными натурного эксперимента подтвердили достаточно высокую точность численного моделирования.

ЛИТЕРАТУРА

[1] Самарский А.А., Моисеенко Б.Д. Экономичная схема сквозного счета для многомерной задачи Стефана. Журн. вычисл. матем. и матем. физ., 1965, т. 5, № 5, с. 816-827.

[2] Будак Б.М., Соболева Е.Н., Успенский А.Б. Разностный метод со сглаживанием коэффициентов для решения задач Стефана. Журн. вычисл. матем. и матем. физ., 1965, т. 5, № 5, с. 828-840.

[3] Danilyuk I.I. On the Stefan problem. Russian Math. Surveys, 1985, vol. 40, no. 5, pp. 157-223. DOI: https://doi.org/10.1070/RM1985v040n05ABEH003684

[4] Васильев В.И. Численное интегрирование дифференциальных уравнений с нелокальными граничными условиями. Якутск, ЯФ СО АН СССР, 1985.

[5] Васильев В.И., Максимов А.М., Петров Е.Е. и др. Тепломассоперенос в промерзающих и протаивающих грунтах. M., Наука, 1996.

[6] Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. M., Либроком, 2009.

[7] Крылов Д.А., Сидняев Н.И., Федотов А.А. Математическое моделирование распределения температурных полей. Математическое моделирование, 2013, т. 25, № 7, с. 3-27.

[8] Vasilyeva M.V., Stepanov S.P., Sirditov I.K. Reduced dimension model for heat transfer of ground heat exchange in permafrost. J. Phys.: Conf. Ser., 2017, vol. 937, art. 012056. DOI: https://doi.org/10.1088/1742-6596/937/1/012056

[9] Stepanov S., Vasilyeva M., Vasil'ev V. Generalized multiscale discontinuous Galerkin method for solving the heat problem with phase change. J. Comput. Appl. Math., 2018, vol. 340, pp. 645-652. DOI: https://doi.org/10.1016/j.cam.2017.12.004

[10] Sivtsev P.V., Kolesov A.E., Sirditov I.K., et al. The numerical solution of thermo-poroelastoplasticity problems. AIP Conf. Proc., 2016, vol. 1773, art. 110010.

DOI: https://doi.org/10.1063/L4965014

[11] Stepanov S.P., Vasilyeva M.V., Vasilyev V.I. Numerical simulation of the convective heat transfer on high-performance computing systems. AIP Conf. Proc., 2016, vol. 1773, art. 110011. DOI: https://doi.org/10.1063/1.4965015

[12] Васильев В.И., Васильева М.В., Сирдитов И.К. и др. Математическое моделирование температурного режима грунтов оснований фундаментов в условиях многолетнемерзлых пород. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2017, № 1 (70), с. 142-159.

DOI: https://doi.org/10.18698/1812-3368-2017-1-142-159

[13] Vasilyeva M., Stepanov S., Spiridonov D., et al. Multiscale Finite Element Method for heat transfer problem during artificial ground freezing. J. Comput. Appl. Math., 2020, vol. 371, art. 112605. DOI: https://doi.org/10.1016/jxam.2019.112605

[14] Vasil'ev V., Vasilyeva M. An accurate approximation of the two-phase Stefan problem with coefficient smoothing. Mathematics, 2020, vol. 8, iss. 11, art. 1924.

DOI: https://doi.org/10.3390/math8111924

[15] Самарский А.А. Введение в теорию разностных схем. М., Наука, 1971.

Васильев Василий Иванович — д-р физ.-мат. наук, профессор, заведующий научно-исследовательской кафедрой «Вычислительные технологии» СВФУ (Российская Федерация, 677000, Якутск, ул. Белинского, д. 58).

Васильева Мария Васильевна — канд. физ.-мат. наук, доцент научно-исследовательской кафедры «Вычислительные технологии» СВФУ (Российская Федерация, 677000, Якутск, ул. Белинского, д. 58); сотрудник Техасского университета A&M (США, Техас, TX 77843-3368, Колледж-Стейшен, Байззел-стрит 400).

Степанов Сергей Павлович — канд. физ.-мат. наук, доцент научно-исследовательской кафедры «Вычислительные технологии» СВФУ (Российская Федерация, 677000, Якутск, ул. Белинского, д. 58).

Сидняев Николай Иванович — д-р техн. наук, профессор, заведующий кафедрой «Высшая математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, корп. 1).

Матвеева Ольга Иннокентьевна — канд. техн. наук, директор ЯкутПНИИС (Российская Федерация, 677000, Якутск, ул. Дзержинского, д. 20).

Цеева Анастасия Николаевна — канд. техн. наук, заведующая отделом оснований и фундаментов ЯкутПНИИС (Российская Федерация, 677000, Якутск, ул. Дзержинского, д. 20).

Просьба ссылаться на эту статью следующим образом:

Васильев В.И., Васильева М.В., Степанов С.П. и др. Решение двухфазной задачи Стефана в энтальпийной постановке со сглаживанием коэффициентов. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2021, № 4 (97), с. 4-23. DOI: https://doi.org/10.18698/1812-3368-2021-4-4-23

NUMERICAL SOLUTION OF THE TWO-PHASE STEFAN PROBLEM IN THE ENTHALPY FORMULATION WITH SMOOTHING THE COEFFICIENTS

V.I. Vasilyev1 vasvasil@mail.ru

M.V. Vasilyeva1' 2 vasilyevadotmdotv@gmail.com

S.P. Stepanov1 cepe2a@inbox.ru

N.I. Sidnyaev3 sidnyaev@bmstu.ru

O.I. Matveeva4 yapniis@mail.ru

A.N. Tseeva4 antseeva@mail.ru

1 North-Eastern Federal University, Yakutsk, Russian Federation

2 Institute for Scientific Computation, Texas A&M University, College Station, USA

3 Bauman Moscow State Technical University, Moscow, Russian Federation

4 Yakutsk State Design Research Institute of Construction, Yakutsk, Russian Federation

Abstract

To simulate heat transfer processes with phase transitions, the classical enthalpy model of Stefan is used, accompanied by phase transformations of the medium with absorption and release of latent heat of a change in the state of aggregation. The paper introduces a solution to the two-phase Stefan problem for a one-dimensional quasilinear second-order parabolic equation with discontinuous coefficients. A method for smearing the Dirac delta function using the smoothing of discontinuous coefficients by smooth functions is proposed. The method is based on the use of the integral of errors and the Gaussian normal distribution with an automated selection of the value of the interval of their smoothing by the desired function (temperature). The discontinuous coefficients are replaced by bounded smooth temperature functions. For the numerical solution, the finite difference method and the finite element method with an automated selection of the smearing and smoothing parameters for the coefficients at each time layer are used. The results of numerical calculations are compared with the solution of Stefan's two-phase self-similar problem — with a mathematical model of the formation of the ice cover of the reservoir. Numerical simulation of the thawing effect of installing additional piles on the existing pile field is carried out. The temperature on the day surface of the base of the structure is set with account for the

Keywords

Mathematical model, cryolithozone, enthalpy, Stefan's problems, coefficient smoothing, latent heat, numerical calculations

amplitude of air temperature fluctuations, taken from the

data of the Yakutsk meteorological station. The study

presents the results of numerical calculations for concrete

piles installed in the summer in large-diameter drilled

wells using cement-sand mortars with a temperature Received 25.02.2021

of 25 °C. The distributions of soil temperature are ob- Accepted 02.04.2021

tained for different points in time © Author(s), 2021

This work was supported by the Russian Science Foundation (project RSF no. 19-11-00230)

REFERENCES

[1] Samarskii A.A., Moiseyenko B.D. An economic continuous calculation scheme for the Stefan multidimensional problem. USSR Comput. Math. Math. Phys., 1965, vol. 5, iss. 5, pp. 43-58. DOI: https://doi.org/10.1016/0041-5553(65)90004-2

[2] Budak B.M., Sobol'eva E.N., Uspenskii A.B. A difference method with coefficient smoothing for the solution of Stefan problems. USSR Comput. Math. Math. Phys., 1965, vol. 5, iss. 5, pp. 59-76. DOI: https://doi.org/10.1016/0041-5553(65)90005-4

[3] Danilyuk I.I. On the Stefan problem. Russian Math. Surveys, 1985, vol. 40, no. 5, pp. 157-223. DOI: https://doi.org/10.1070/RM1985v040n05ABEH003684

[4] Vasilyev V.I. Chislennoe integrirovanie differentsial'nykh uravneniy s nelokal'nymi granichnymi usloviyami. [Numerical integration of differential equations with nonlocal boundary conditions] Yakutsk, YaF SO AN SSSR Publ., 1985.

[5] Vasilyev V.I., Maksimov A.M., Petrov E.E., et al. Teplomassoperenos v pro-merzayushchikh i protaivayushchikh gruntakh [Heat and mass transfer in frozen and melted soils]. Moscow, Nauka Publ., 1996.

[6] Samarskiy A.A., Vabishchevich P.N. Vychislitel'naya teploperedacha [Computational heat and mass transfer]. Moscow, Librokom Publ., 2009.

[7] Krylov D.A., Sidnyaev N.I., Fedotov A.A. Mathematical modelling of temperature distribution. Matematicheskoe modelirovanie, 2013, vol. 25, no. 7, pp. 3-27 (in Russ.).

[8] Vasilyeva M.V., Stepanov S.P., Sirditov I.K. Reduced dimension model for heat transfer of ground heat exchange in permafrost. J. Phys.: Conf. Ser., 2017, vol. 937, art. 012056. DOI: https://doi.org/10.1088/1742-6596/937/1Z012056

[9] Stepanov S., Vasilyeva M., Vasil'ev V. Generalized multiscale discontinuous Galerkin method for solving the heat problem with phase change. J. Comput. Appl. Math., 2018, vol. 340, pp. 645-652. DOI: https://doi.org/10.1016/jxam.2017.12.004

[10] Sivtsev P.V., Kolesov A.E., Sirditov I.K., et al. The numerical solution of thermo-poroelastoplasticity problems. AIP Conf. Proc., 2016, vol. 1773, art. 110010.

DOI: https://doi.org/10.1063/L4965014

[11] Stepanov S.P., Vasilyeva M.V., Vasilyev V.I. Numerical simulation of the con-vective heat transfer on high-performance computing systems. AIP Conf. Proc., 2016, vol. 1773, art. 110011. DOI: https://doi.org/10.1063/L4965015

[12] Vasilyev V.I., Vasilyeva M.V., Sirditov I.K., et al. Mathematical modeling of temperature regime of soils of foundation on permafrost. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2017, no. 1 (70), pp. 142-159 (in Russ.). DOI: https://doi.org/10.18698/1812-3368-2017-1-142-159

[13] Vasilyeva M., Stepanov S., Spiridonov D., et al. Multiscale Finite Element Method for heat transfer problem during artificial ground freezing. J. Comput. Appl. Math., 2020, vol. 371, art. 112605. DOI: https://doi.org/10.1016/j.cam.2019.112605

[14] Vasil'ev V., Vasilyeva M. An accurate approximation of the two-phase Stefan problem with coefficient smoothing. Mathematics, 2020, vol. 8, iss. 11, art. 1924.

DOI: https://doi.org/10.3390/math8111924

[15] Samarskiy A.A. Vvedenie v teoriyu raznostnykh skhem [Introduction into theory of difference schemes]. Moscow, Nauka Publ., 1971.

Vasilyev V.I. — Dr. Sc. (Phys.-Math.), Professor, Head of the Research Department of Computing Technologies, North-Eastern Federal University (Belinskogo ul. 58, Yakutsk, 677000 Russian Federation).

Vasilyeva M.V. — Cand. Sc. (Phys.-Math.), Assoc. Professor, Research Department of Computing Technologies, North-Eastern Federal University (Belinskogo ul. 58, Yakutsk, 677000 Russian Federation); employee, Institute for Scientific Computation, Texas A&M University (College Station, TX 77843-3368, 400 Bizzell St.).

Stepanov S.P. — Cand. Sc. (Phys.-Math.), Assoc. Professor, Research Department of Computing Technologies, North-Eastern Federal University (Belinskogo ul. 58, Yakutsk, 677000 Russian Federation).

Sidnyaev N.I. — Dr. Sc. (Eng.), Professor, Head of the Department of Higher Mathematics, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5/1, Moscow, 105005 Russian Federation).

Matveeva O.I. — Cand. Sc. (Eng.), Director of Yakutsk State Design Research Institute of Construction (Dzerzhinskogo ul. 20, Yakutsk, 677000 Russian Federation).

Tseeva A.N. — Cand. Sc. (Eng.), Head of the Department of Basements and Foundations, Yakutsk State Design Research Institute of Construction (Dzerzhinskogo ul. 20, Yakutsk, 677000 Russian Federation).

Please cite this article in English as:

Vasilyev V.I., Vasilyeva M.V., Stepanov S.P., et al. Numerical solution of the two-phase Stefan problem in the enthalpy formulation with smoothing the coefficients. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2021, no. 4 (97), pp. 4-23 (in Russ.). DOI: https://doi.org/10.18698/1812-3368-2021-4-4-23