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

Наука к Образование

МГТУ им. Н.Э. Баумана

Сетевое научное издание

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 10. С. 419-437.

Б01: 10.7463/1015.0814769

Представлена в редакцию: Исправлена:

© МГТУ им. Н.Э. Баумана

15.09.2015 29.09.2015

УДК 517

Интегральные уравнения задачи Стефана и их приложение при моделировании оттаивания грунта

АруТЮНЯН Р. В.1'2' ' гоЪ57@тайл1

:МГТУ им. Н.Э. Баумана, Москва, Россия 2Московский технический университет связи и информатики,

Москва, Россия

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

Ключевые слова: задача Стефана, интегральные уравнения, численный метод, оттаивание грунтов

Введение и постановка задачи Стефана

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

Наиболее известны конечно-разностные методы решения задачи Стефана.

В [9] моделировалась динамика температурного поля грунтов основания здания, установленного полами по грунту в условиях сплошной мерзлоты в районах с сезонными промерзаниями и оттаиваниями почвы. Результаты моделирования позволяют анализировать опасные воздействия здания на грунт.

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

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

Соответствующая краевая задача для уравнения теплопроводности называется задачей Стефана и имеет следующий вид:

рс^ = й/ V(А дгай и) + /(М, О , М Е V\БГ , t > 0, (1а)

~ 1^ = РГЛ, и+=и~ =иг, М Е Бг , t > 0, (1б)

+ ри = ц,М ЕБ; и(М,0) = и0 (М),\ и(М,г) \ < <х>,М ЕV,t > 0, (1в)

где V - область пространства Ет, не имеющая разрезов, с кусочно-гладкой границей 8; и -функция температурного поля при наличии фазовых превращений на поверхности 8/, определяемой условиями (1), и/ - температура фазового перехода, индекс "-" соответствует области жидкой, а индекс "+" - области твердой фазы, ио(М) - начальное

ди—

распределение температур в области V, = д г а й и (М ±, ^ п^ (М , ^ , где п^ (М , £) - вектор

единичной нормали в точке М поверхности 8/, в момент времени I, внешней по отношению к области V/: ¥/= { М Е V \ и (М, 0 > и ^ } ,

д г а й и (М 0 = И т и-м, и (М , 0 , д г а й и (М+, 0 = И т N -м, и (М , 0 , М Е Б^, ^ > 0,

иеуш меу\у„л

ди—

— = д г а й и (М ±, 0 пг (М , 0 , vn (М, 0 = V (М , О пг (М, 0 ,

где V (М , £) - скорость движения поверхности в точке М Е Б ^ при I > 0.

Функции р (М),с (М),А (М)(М,£), гпл(М),р (М) будем считать таковыми, что существует классическое решение задачи (1) и(М/) и классическое решение задачи (1) при гпл = 0, которое обозначим И(М/). Это будет выполнено, например, если р (М),с(М),А (М),[ (М,г),гш (М),(3 (М) - положительные константы, /ЕЬ2^),цЕ I2 (Б) ,и0(М)Е С1 (V) [1-8].

Уравнение (1) приводится к эквивалентному виду с т.н. сосредоточенной теплоемкостью, учитывающей теплоту фазового перехода [7,8]:

р [с + гш S (и — Uf) ] = d i v (Л д г a d и) + f (М , t) ,М EV,t > 0, (2а)

Л-^ + (]и = ц, М Е S,t > 0, (2б)

и(М , 0) = и0(М), \и(М , t)\ < const < да, М EV, (2в)

где 5(х) - дельта-функция Дирака.

Представим u(M,t) как сумму u = w + H, где H - решение краевой задачи без учета стефановской нелинейности. Можно показать, что w удовлетворяет условиям:

pew + ргш [ 1 (w + Н — и^ — 1 (и0 — u^] = div(Л grad W),М EV,t > 0, (3а)

Л-^ + Рw = 0,М ES; w(М,0) = 0, \w (M,t) \ < сonst , М Е V , (36)

где 1(x) - единичная функция Хевисайда, W (М , t) = J^w (М, т) dx.

Пусть G(M,N,t) - функция Грина краевой задачи, получающейся из (3) при гпл = 0. Тогда решение (3) удовлетворяет уравнению:

t

W(M,t) = 11 p(N)rm(N)G(M,N,t -T)[l(ii0(JV) - uf) - l(ii(JV,r) - uf)]dVNdx,

0 v

ME V,t> 0.

Продифференцировав последнее равенство по времени, получим (с учетом известных свойств функции Грина):

w(M,t) = prf]JI[l(u0 — Uf) - l(u-uf)]+f'fv P(iV)^(iV)|G(M,iV,t-T)[l(uo(iV) -

u/) — 1 (и(N ,т) — и/)]dVNdx,М E V,t > 0. (4)

Введем в рассмотрение функцию tfM), являющуюся наименьшим корнем уравнения u(M,tf(M)+0) = Uf +0, М Е V . В точках М Е V , где такого решения не существует, положим значение tfM) бесконечным. Функция tfM) имеет смысл времени фазового превращения вещества в точке M (в рассматриваемом случае из твердого в жидкое). В силу определения tf(M): V f (t) = {М EV \ t f (М) < t } , t < в.

Откуда Sf (t) = {М E V \ tf (М) = t } , t <в , а так как V t <в, М Е V : 1 (и (М, t) — и/) = 1 (t — t f (М) ) , то (4) преобразуется к виду:

w(МЛ) = —Jv р(N)rm(N)G(M,N,t — tf (N)) 1 (t — tf (N))dVN, М Е V,t < в. (5) Учитывая, что при t = tfM): w(М, t) = и/ — Н(М , tf (М)), получаем из (5) нелинейное интегральное уравнение Фредгольма 2 рода относительно функции tfM):

fv р(Л0 rfl/l(N)G (м, N, tfW - tf(N)) 1 (t^M) - tf(N)) dVN = = Н (м , tf (М) ) —щ, М Е V. (6)

1. Методы численного решения нелинейного интегрального

уравнения (6)

1.1. Метод ячеек

Осуществим разбиение области V на элементы V], у=1,...,У (У может быть бесконечным), .

В пределах у : tf (М) - ty=const, р(М)гил(М) - р ¡гпл],Н(МЛ) - Н)(0. Тогда получаем следующую систему нелинейных уравнений относительно неизвестных (¿ь...,

г^хц( и - ^) 1 - ^) = Щ(^ -^Л = 1.....М; (7)

ху(!) = I pyC/■G(M¿,iV,t)dK^v,M¿ ЕУ0

где - узел коллокации, как правило, находящийся в центре ячейки .

1.2. Алгоритмы решения системы нелинейных уравнений (7)

Если размерность системы (7) бесконечна, число непересекающихся связных подмножеств конечно и не возрастает во времени, то эффективным может быть

следующий метод решения. Исходным является следующее свойство:

V/ е у/,у е дУр: ь < ц, - ty) = о. Обозначим для краткости через , тогда алгоритм примет вид:

1. а) задаем конечное число .

б) определяем начальное множество узлов , для которых

= Нц^ч*) = иТ> = = =

в) задаем число . Если компоненты вектора будут различаться менее, чем на , то они будут считаться равными.

2. Формируем множество узлов до , граничных к о.

3. Определяем размерность до: п = с1 1 т до, до ={г,...,]п}.

4. Присваиваем р значение, равное 1.

5. Полагаем хт 1П = х

6. Вычисляем

тах> Р — 1< тт\£тах> 21> ■■■ > гр-1) >Р > 17. Находим корень уравнения:

- иг - X Гш]!с) хЬЛ2Р - - ь) = °<

}ео

лежащий в интервале , .

8. Увеличиваем значение р на единицу.

9. Проверяем не исчерпаны ли элементы д<, что равносильно неравенству р>п. Если да, то переходим к п.12, иначе к п.6.

10. Полагаем г* = гш Пур(Е д а1р.

11. Проверяем выполнение неравенства . Если оно справедливо, то вычисления завершены, иначе переходим к п.12.

12. Формируем множество узлов

ш = {/р е д<| | гр - г*| < г}.

13. Расширяем множество <, присоединяя к нему множество ш.

14. Переходим к п.3.

Из рассмотрения описанного метода следует, что в нем строится неубывающая последовательность значений компонент вектора решения системы (12)

Ьй:Ь£ < ■ ■ ■ < Ь£ < ■ ■ ■

11 --- 1П -

1.3. Вычисление сингулярного слагаемого в сумме (7)

Слагаемое х £ у ( Ь — Ьу) при / = _/' является сингулярным, так как функция С( М£, М£, Ь) неограниченна в окрестности нуля . При этом вклад данного слагаемого относительно значителен. По этим причинам существует целесообразность вычисления сингулярного слагаемого аналитическими методами. Для этого в соответствующей ячейке К£ функция

Ьу (М) аппроксимируется линейной частью ряда Тейлора: Ьу(М ) « Ьу (М£) + к • ММ1.

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

центром в начале координат. Рассмотрим вычисление сингулярного слагаемого в одномерном случае (для краткости изложения считаем размер ячеек постоянным).

Хи = I р&в (х£,х, t/(x) - 1 (t/(x) - йх,

XI 2

к к X е [х£ --,х£ + -],

В области сингулярности рсС(М,М, Ь) ~( 2 л/тгаЬ) Ш ехр (—( 4 а Ь) ). Для одномерного случая: Ьу (х) « Ьу ( х £) + к(х — х £),

* (х - х£)2 "

Х1+7

ехр

1 2, пак(х - хЛ } 2л/7га/сх

XI 2

|8а/с

и

2л/пак(х — х£) ^ 2л1пакх

Н (1 + о(1)), Л^О.

2пак

При вычислениях можно использовать для оценки значения к конечно-разностную

С/ОсО- С/О^-!)

аппроксимацию производных: к —-.

2. Примеры практического применения описанного метода

2.1. Задача Коши в Ет

При постоянных :

Н (М , 0 = ¡Ет и0(М) С 0( М, М , 0 йУ„ + £ ¡Ет ^ С о (М , М, t - т) йУ„йт. (8)

* (9)

(10) (11)

С 0(М,М,г)=( 2 л1паг) е хр ( - / (4 а 0 ) ,а = Л / (р с) . Для метода ячеек (7) компоненты матрицы системы равны:

и

■") = '/2

\ 2 л/а? \ 2л/аГ

2.2. Моделирование оттаивания мерзлых грунтов

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

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

При математическом моделировании рассматриваются два этапа процесса оттаивания:

1) Нагрев грунта до температуры плавления (твердая фаза);

2) Нагрев оттаявшего и плавление мерзлого грунтов (жидкая фаза).

На 1 -ом этапе температура на поверхности изменяется от отрицательных значений до 0° С при помощи СВЧ энергии.

Рис. 1. Схема процесса

1-й этап продолжается до появления межфазной границы, то есть до достижения в верхнем слоя грунта температуры плавления льда (0° С). На этом этапе решается уравнение теплопроводности вида:

дщ л д \ „ , ч „ „

р1С1~Ш = Л1к1+ 1 > 0 < х < Н' (12)

с начальным условием:

Щ = инач(х) пРи г = (13)

и граничными условиями:

¿Щ ,, ч г, ди1

—1 = к (и - и ) при х = 0 и —1 = 0 при х = п

¿х °р дх

(14)

где с - удельная теплоемкость, р - плотность, Л - теплопроводность, щ - температура грунта (все параметры относятся к мерзлому грунту); г - время, х - координата, < -мощность внешнего источника тепла.

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

Рис. 2. Схема процесса оттаивания грунта

На втором этапе решается система уравнений теплопроводности с учетом фазового перехода (то есть задача Стефана) вида:

дщ д Щ Р- Сг^2 = Л2—2 + 2- (х, г),° < х < ^),

ддг „ д дх' . ч . ч

РЛ —1 = Л —Г + <1 (х, г), ) < х < п

дг

дх

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

и граничными условиями:

Щ2 = Щ/, Щ1 = Щ1кон

¿и

Л—1 -Л2~Г = Р1Гш—Г- , Щ1 = и2 = иГ пРи х = ^(0 ах ах аг

дщ п ,

—1 = 0 при х = п .

дх

- = к(щ - и ) при х = 0,

ах

аф)

(15)

(15а)

(156 ) (15в) (15г )

<

где Q = А • е 2ох, априорно неизвестная функция х = 77^) определяет положение

границы оттаивания и находится в процессе решения системы дифференциальных уравнений.

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

2.2.2. Решение задачи методом интегральных уравнений 2.2.2.1. Характеристики грунтов

Для широкого класса грунтов, значения теплофизических характеристик приблизительно одинаковы для мерзлых и талых грунтов [9,10]. В качестве примера рассмотрим следующие данные / - влажность грунта). Характеристики талого грунта:

/= 0,02 и р=1800 кг/м3: Х=0,47 Вт/м-К, с=778 Дж/(кг-К); /= 0,12 и р =1800 кг/м3: Х=0,85; с=1278 Дж/ (кг-К); Характеристики мерзлого грунта:

/=0,03 и р =1700 кг/м3: Х=0,56 Вт/м-К, с=750 Дж/ (кг-К); /= 0,12 и р =1800 кг/м3: X =0,78 Вт/м-К; с = 1222 Дж/ (кг-К); Таким образом, погрешность усредненных значений пренебрежимо мала. Среднемассовую теплоту таяния грунта считаем равной гпл=335000/;

2.2.2.2. Оттаивание грунта при помощи поверхностного нагрева

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

Выражение для температурного поля согласно описанному в статье методу имеет

вид:

и

ехр

(х~хм)

-ехр

(х+хм)

I Jлa(t-tf(XN))

-

Н(х, 0 = (ит — и0)ег/с (-

X \

: +и0.

<2л/аГ/

Интегральное уравнение (6) для рассматриваемого случая принимает вид:

и,

ехр

(х

4а (гг(х) -

ехр

(х + У

4а (^(х) -

(¿Хдг,

2 ¡па (^(х) -

х > 0.

Полученное интегральное уравнение решается аналитически. Решение отыскивается в виде: ¿у(х) = кх2, где к — неизвестный постоянный параметр.

Обозначим р = ^Цр тогда интегральное уравнение приводится после подстановки ¿у (х) = х2/ ( 4 ар ) к нелинейному уравнению:

2гпл _ г77/4

= (% — ^0)ег/с(л/р) + и0--—д/р/я I (ехр(—р^2х) — ехр(—рЛд2х)) йх.

с ) о

Рассматриваемое нелинейное уравнение после отделения корня успешно решается методом хорд и касательных.

Результаты расчетов. При влажности 12% и плотности грунта (как мерзлого, так и

о

талого) р =1800 кг/м .

Я 1=0,78 Вт/м ■ К; с 1=1222 Дж/ (кг- К); Я2=0,85 Вт/м- К; с2=1278Дж/ (кг- К); Соответствующие средние значения: Я=0,815 Вт/м- К; с=1250 Дж/ (кг- К); (погрешность около 4%).

Значения температур в градусах Цельсия: и0=—10; ит=2000; иу=0. Вычисленные значения параметров:

_7

при невязке уравнения 10 : k = 249664; p = 2,764. Для оттаивания грунта до глубины в 1 м потребуется 69,35 час. Для сравнения, если не учитывать скрытую теплоту таяния (т.е. при гпл = 0):

_7

при невязке уравнения 4,5-10 : k = 174985; p = 3,944. Для оттаивания грунта до глубины в 1 м потребовалось бы 48,61 час.

2.2.2.3. Сравнение численного решения методом интегральных

уравнений и точного решений

щЛ = 1,

i

и 7

#¿(0 = (ит - и0)ег[с + щ.

= -

егП

егИ

2 л/а1/

2л/аГ / \ 2л/аГ / \ 2л/аГ

Выражение сингулярного слагаемого ( ):

В =

*и(!) = егП

к \ Вехр{-В2)

8 ак

(2/ -1) ^

+ В2ег/с(В),

, с учетом равенства к

■:В =

(24>

при .

П " V2а( Г- Сц)'

В таблице 1 представлены значения времен таяния на различной глубине. Шаг по координате равен 1 см. х£ — координаты фронта фазового перехода, — результаты расчета методом интегральных уравнений, ¿у (х£) — точные значения времен плавления.

Таблица 1. Значения времен таяния на различной глубине

х;, см час ^ (х) , час

4,5 0, 14 0,14

9,5 0,63 0,63

19,5 1,46 1,46

24,5 2,63 2,64

29,5 4,16 4,16

34,5 8,25 8,25

39,5 10,81 10,82

44,5 13,72 13,73

49,5 16,98 16,99

54,5 20,58 20,60

59,5 24,53 24,55

64,5 28,83 28,85

69,5 33,48 33,50

74,5 38,47 38,49

79,5 43,80 43,83

84,5 49,49 49,52

89,5 55,52 55,55

94,5 61,90 61,93

99,5 68,62 68,66

Данные таблицы 1 показывают высокую точность рассматриваемого численного метода. Затраты машинного времени составляют порядка 1 сек на 1 шаг по координате. Примечательная особенность рассматриваемых задач - даже значительное увеличение шага интегрирования не приводит к существенному росту погрешности.

2.2.2.4. Оттаивание грунта при помощи СВЧ нагрева

Данный метод является наиболее новым и эффективным. Сравнение СВЧ -оттаивания с другими способами подготовки мерзлых пород к экскавации позволило сделать вывод о перспективности этого способа в тех случаях, где требуется скоростное проведение работ: проходка шурфов, экстренное вскрытие при авариях, забивка свай и т.п. [9-12].

Плотность мощности тепловыделения на единицу объема при СВЧ нагреве диэлектрика: где - угловая частота ЭМП, - мнимая часть комплексной

относительной диэлектрической проницаемости диэлектрика, -диэлектрическая постоянная, Е - напряженность электрического поля. Если проводимость у материала мала, или частота настолько велика, что ( ), то диэлектрический нагрев

является доминирующим механизмом передачи энергии ЭМП в вещество.

Коэффициент, характеризующий глубину проникновения ЭМП в грунт, равен

мнимой части волнового числа: а = ^ 8 V ( I £г I — ^ е£г)/ 2■

2.2.2.5. Расчет ЭМП в грунте при оттаивании

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

жидкой фазе, - твердой фазе, - подстилающей поверхности, -

диэлектрические проницаемости соответствующего слоя.

На первом этапе осуществляется нагрев диэлектрика — «твердая фаза», после образования жидкой фазы начинается второй этап оттаивания.

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

Нормированные на единицу комплексные амплитуды электрического и

магнитного полей в -ом слое плоскослоистой структуры представим в виде:

Е&) = Ь^'^2 + це^.Н^г) = --*-,

Щ

где Ъо = 1,% = 0, / = 0, 1,...,М;к £ = ш^- комплексное волновое число, мнимая часть которого характеризует затухание электромагнитной волны в слое диэлектрика; Ъ ¿, ту -коэффициенты прохождения и отражения соответственно, а - волновое сопротивление. Из условия непрерывности электрического и магнитного полей на границах слоев получим систему уравнений для определения коэффициентов прохождения и отражения на первом этапе:

1 + г0 = Ъ2 + г2,-=-,

и/0 и/2

„ Ь2е~1к^ - г2е^2 ь е-]кЗп2

Ъ2е~] 2 2 + г2е] 2 2 = Ъъе~]к^2,—---= —-,

\м2 \л?3

и на втором этапе:

1 + г0 = Ь1 + гх,-=-,

Ъ1в-+ пе^1"1 = Ь2е~1к2П 1 + г2е^2П\ _ Ъ2е~1к2,11 — г2е^к2}11

ъ е-^12 + г еул2л12 = ь3е~1кзп^,—---= —-,

ц/2

где .

Решение этой системы позволяет определить комплексные амплитуды поля падающей и отраженной волны. Более детально вопросы расчета ЭМП СВЧ в многофазной среде при плавлении диэлектрика рассмотрены в [13].

2.2.2.6. Краевые условия I рода для уравнения теплопроводности

Характеристики грунта те же, что и в предыдущем пункте. Значения температур в градусах Цельсия: на поверхности - ит=10; на большой глубине - и0=-10; температура фазового перехода - иу=0.

В таблице 2 представлены значения времен таяния на различной глубине. Шаг по координате равен 1 см. - координаты фронта фазового перехода, - результаты расчета методом интегральных уравнений. Исходные данные для расчета:

Вещественная часть относительной диэлектрической проницаемости мерзлого грунта: ;

вещественная часть относительной диэлектрической проницаемости талого грунта:

Я е гг2 = Я е £г х( 1 — /) + 7 О/, мнимая часть относительной диэлектрической проницаемости мерзлого грунта:

/ш £г х = О, 4;

мнимая часть относительной диэлектрической проницаемости талого грунта:

/ ш £г 2 = / ш £г х( 1 — /) + 1 5 /,

циклическая частота:

суммарный коэффициент поглощения ЭМП: Ъ = а х + а2 ,

ч В

амплитуда напряженности электрического поля: Ет = 2 ■ 1 О л -

м

2.2.2.7. Температура при нулевой теплоте таяния

находится посредством решения краевой задачи:

рсН1 = АЯХХ + Ае~Ъх,х > 0; Я(0, 0 = ит, Я(х, 0) = и0, где свободный член описывает диэлектрический нагрев при помощи СВЧ. Решение краевой задачи находится аналитически и имеет вид:

(х \ Ас ^^

гЗм)+ щ+ ~яР~ ~ ^

На рис.3. показан график рассматриваемого распределения температур.

Рис.3. График функции Н (х, £) при t =10 час.

2.2.2.8. Результаты расчетов

В таблице 2 представлены значения времен таяния на различной глубине. Шаг по координате равен 1 см. - координаты фронта фазового перехода, - результаты расчета методом интегральных уравнений.

Таблица 2. Значения времен таяния на различной глубине

х;, см час

0.05 0.27

0.10 0.40

0.15 0.56

0.20 0.79

0.25 1.12

0.30 1.57

0.35 2.20

0.40 3.04

0.45 4.16

0.50 5.63

2.2.2.9. Краевые условия II рода

Характеристики грунта и исходные данные для расчета те же, что и в предыдущем пункте. Значения градиента температур на поверхности и на большой глубине

нулевые: ^ = 0. Начальное распределение: u (х, 0 ) = и0 , х > 0 . Решение соответствующей краевой задачи при нулевой теплоте фазового перехода находится аналитически и имеет вид:

На рис.4. показан график рассматриваемого распределения температур.

2.2.2.10. Результаты расчетов

Рис. 4. График функции Н (х, при t =1 час.

В таблице 3 представлены значения времен таяния на различной глубине. Шаг по координате равен 1 см. — координаты фронта фазового перехода, — результаты расчета методом интегральных уравнений.

х, см час

0.02 0.15

0.03 0.21

0.04 0.25

0.05 0.28

0.06 0.31

0.07 0.33

0.08 0.35

0.09 0.37

0.10 0.40

0.11 0.42

0.12 0.45

0.13 0.48

0.14 0.52

0.15 0.56

0.16 0.60

0.17 0.64

0.18 0.69

0.19 0.74

0.20 0.79

0.25 1.12

0.30 1.57

0.35 2.20

0.40 3.04

0.45 4.16

0.50 5.63

3. Заключение

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

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

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

Недостатком метода является трудоемкость построения функции Грина в случае переменных свойств среды.

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

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

1. Каменомостская С.Л. Проблема Стефана // Математический сборник. 1961. Т. 53, № 4. С. 489-514.

2. Мейрманов A.M. Задача Стефана. Новосибирск: Наука, Сиб. отд., 1986. 187 с.

3. Лыков А.В. Тепломассообмен. М.: Энергия, 1971. 560 с.

4. Данилюк И.И. О задаче Стефана // Успехи математических наук. 1985. Т. 40, вып. 5 (245). С. 133-185.

5. Рубинштейн Л.И. К вопросу о численном решении интегральных уравнений задачи Стефана // Известия вузов. Математика. 1958. № 4. С. 202-214.

6. Dauzhenka T.A., Gishkeluk I.A. Quasilinear Heat Equation in Three Dimensions and Stefan Problem in Permafrost Soils in the Frame of Alternating Directions Finite Difference Scheme // Proceedings of the World Congress on Engineering (WCE 2013), July 3-5, 2013, London, U.K. Vol. 1. Available at:

http://www.iaeng.org/publication/WCE2013/WCE2013pp1 -6.pdf , accessed 01.09.2015.

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

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

9. Гласко А.В., Федотов А.А., Сидняев Н.И., Храпов П.В., Мельникова Ю.С. Моделирование динамики температурного поля грунтов основания здания в криолитозоне // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2011. № 12. Режим доступа: http://technomag.bmstu.ru/doc/274059.html (дата обращения 01.09.2015).

10. Мельникова Ю.С. Математическое моделирование управления нестационарным температурным полем в двухфазных средах // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 2. Режим доступа: http://technomag.bmstu.ru/doc/330390.html (дата обращения 01.09.2015).

11. Рожин И.И. Численное моделирование переходных процессов в прикладных задачах теплопроводности с фазовыми превращениями: автореф. дис. ... канд. физ.-мат. наук. Якутск, 2005. 15 с.

12. Теплофизические свойства грунтов // MINKOR: сайт компании. Режим доступа: http://minkor.ru/upload/spravochnik/190609-2.pdf (дата обращения 01.09.2015).

13. Анфиногентов В.И. Математическое моделирование СВЧ нагрева диэлектриков: дис. ... докт. техн. наук. Казань, 2006. 340 с.

Science^Education

of the Bauman MSTU

Science and Education of the Bauman MSTU, 2015, no. 10, pp. 419-437.

DOI: 10.7463/1015.0814769

Received: 15.09.2015

Revised: 29.09.2015

ISS N 1994-0408 © Bauman Moscow State Technical Unversity

Integral Equations of the Stefan Problem and Their Application in Modeling of Thawing Soil

R.V. Arutjunjan1'2'* 'rob57@maiLru

1Bauman Moscow State Technical University, Moscow, Russia 2Moscow Technical University of Communications and Informatics,

Moscow, Russia

Keywords: Stefan problem, integral equations, numerical methods, thawing soils

The paper offers a new numerical method for solving the Stefan problem, based on the reduction of the initial boundary value problem for nonlinear integral equation equivalent to the minimum dimension. A feature of the method is that the time-varying function of phase transition is used as a major unknown function. This parameterization is an efficient computational technique enabling the considerably simplified method of solution. Computational qualities of this method are investigated by the example of solving the practical problems of thawing frozen soils under thermal action produced by various ways.

The method advantages are most evident in solving problems in an unbounded domain, as unlike traditional finite difference method there is no need to truncate the system of equations with a large possibly infinite dimension.

Despite the complexity to construct Green's function, in the case of variable properties of the environment there is a broad class of problems concerning the simulation of phase transitions. Their solution can be successfully found by described method of non-linear integral equations.

The method also provides additional convenience in Richardson extrapolation and in solution of Stephen inverse problems, in particular, problems of front control of phase transition.

References

1. Kamenomostskaya S.L. On Stefan's problem. Matematicheskii sbornik, 1961, vol. 53, no. 4, pp. 489-514. (in Russian).

2. Meirmanov A.M. Zadacha Stefana [Stefan problem]. Novosibirsk, Nauka Publ., 1986. 187 p. (in Russian).

3. Lykov A.V. Teplomassoobmen [Heat and mass transfer]. Moscow, Energiya Publ., 1971. 560 p. (in Russian).

4. Danilyuk I.I. On the Stefan problem. Uspekhi matematicheskikh nauk, 1985, vol. 40, no. 5, pp. 133-185. (English version of journal: Russian Mathematical Surveys, 1985, vol. 40, no. 5, pp. 157-224. DOI: 10.1070/RM1985v040n05ABEH003684 ).

5. Rubinshtein L.I. The numerical solution of the integral equations of the Stefan problem. Izvestiya vysshikh uchebnykh zavedenii. Matematika, 1958, no. 4, pp. 202-214. (in Russian).

6. Dauzhenka T.A., Gishkeluk I.A. Quasilinear Heat Equation in Three Dimensions and Stefan Problem in Permafrost Soils in the Frame of Alternating Directions Finite Difference Schem. Proceedings of the World Congress on Engineering (WCE 2013), July 3-5, 2013, London, U.K. Vol. 1. Available at:

http://www.iaeng.org/publication/WCE2013/WCE2013pp1 -6.pdf , accessed 01.09.2015.

7. Samarskii A.A., Moiseenko B.D. An economic continuous calculation scheme for the Stefan multidimensional problem. Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki, 1965, vol. 5, no. 5, pp. 816-827. (English version of journal: USSR Computational Mathematics and Mathematical Physics, 1965, vol. 5, iss. 5, pp. 43-58. DOI: ).

8. Budak B.M., Sobol'eva E.N., Uspenskii A.B. A difference method with coefficient smoothing for the solution of Stefan problems. Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki, 1965, vol. 5, no. 5, pp. 828-840. (English version of journal: USSR Computational Mathematics and Mathematical Physics, 1965, vol. 5, iss. 5, pp. 59-76. DOI: 10.1016/0041-5553(65)90005-4 ).

9. Glasko A.V., Fedotov A.A., Sidnyaev N. I., Hrapov P.V., Mel'nikova Yu.S. Modeling the dynamics of temperature field soil base of the building in permafrost. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2011, no. 12. Available at: http://technomag.bmstu.ru/doc/274059.html , accessed 01.09.2015. (in Russian).

10. Mel'nikova Yu.S. Mathematical simulation of time-dependent temperature field control in two-phase mediums. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2012, no. 2. Available at:

http://technomag.bmstu.ru/doc/330390.html , accessed 01.09.2015. (in Russian).

11. Rozhin I.I. Chislennoe modelirovanie perekhodnykh protsessov v prikladnykh zadachakh teploprovodnosti s fazovymi prevrashcheniyami. Avtoreferat kand. diss. [Numerical simulation of transient processes in applied heat transfer problems with phase transitions. Abstract of cand. diss.]. Yakutsk, 2005. 15 p. (in Russian).

12. Teplofizicheskie svoistva gruntov [Thermalphysic properties of soils]. MINKOR: company website. Available at: http://minkor.ru/upload/spravochnik/190609-2.pdf , accessed 01.09.2015. (in Russian).

13. Anfinogentov V.I. Matematicheskoe modelirovanie SVCh nagreva dielektrikov. Dokt. diss. [Mathematical modeling of microwave heating of dielectrics. Dr. diss.]. Kazan, 2006. 340 p. (in Russian).