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

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

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

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2016. № 08. С. 165-181.

DOI: 10.7463/0816.0843482

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

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

08.07.2016 22.07.2016

УДК 517

Граничные интегральные уравнения многомерной задачи Стефана и их численное решение

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

Московский технический университет связи и информатики,

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

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

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

1. Введение и постановка задачи

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

Используемые обозначения полностью соответствуют [1]. Согласно [1] задача Стефана приводится к виду с т.н. сосредоточенной теплоемкостью, учитывающей теплоту фазового перехода [8-12]:

Р [ с + rnnS (и - uf) ] = d iv (X g r ad u) + f (M, t) ,M EV,t> 0 , (1a)

du

X— + Ru = q, MeS,t>0, (1б)

дп

< const < , (1 в)

где 5(x) - дельта-функция Дирака, V - область m-мерного евклидового пространства Em, не имеющая разрезов (V\V = дV), с кусочно-гладкой границей S = дV, дV = д V; u - функ-

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

Поверхность 8/в момент времени I к ограничивает область

Уг: ¥/= {М е V| и(М, £) > иг }.

Функции р (М ) , с (М) , А (М) ,/ (М, 0 , гпя(М),(3 (М ) считаются таковыми, что существует классическое решение задачи (1) и(М/) и классическое решение задачи (1) при гпл = 0, которое обозначим И(М/). Согласно [1-9] это будет выполнено, например, если р (М ) , с (М) , А (М) ,/ (М, 0 , гпл(М),// (М ) будут положительными константами, / е ¿200,ч е I 2(5) ,и 0( М ) е С1 (V) .

Если нормальная к 8/ скорость перемещения фронта фазового перехода гп ( М , 0 > 0 , то выполнится условие:

V/ ( £х ( £ 2) ,0<^<Г2 < 0 , (2)

где 0 < £ < 0 - интервал времени, на котором происходит расширение области плавления.

Задача (1) приводится к интегральному уравнению минимальной размерности [1]. Пусть G(M,N,t) - функция Грина линейной краевой задачи, получающейся из (1) при гпл = 0. Введем в рассмотрение функцию /М), являющуюся наименьшим корнем уравнения и(М,/М)+0) = и/+0, М е V . В точках М е V , где такого решения не существует, положим значение /М) бесконечным. Функция /М) имеет смысл времени фазового превращения вещества в точке М (в рассматриваемом случае из твердого в жидкое).

В силу определения /М) и при выполнении (2):

V/( £) = {М е V | £г(М) < £ }, £ < 0.

Откуда

( £) = {М е V | £г (М ) = £ } , £ < 0 ,

так как

V £ < 0 , М е К : 1 (и (М, £) - и7) = 1 (£ - ^ (М) ) ,

Учитывая, что при

г = г^М): IV(М, £) = иг - Я (М, £г (М ) ) , получаем из (1) нелинейное интегральное уравнение Фредгольма 2 рода относительно функции /М):

| р(Л) гпл(Х)С (м, л, ^(м) - 1 (гДм) - гДл)) сЯ^ = у

= Я (м, £г (м) ) - и/; М е К. (3)

2. Исследование интегрального уравнения (3)

Существование кусочно-непрерывного решения уравнения (3) следует из его вывода. Для доказательства единственности любого кусочно-непрерывного решения (3) в области Vy ( в ) покажем, что справедливо

Утверждение 1. Пусть z(M) - кусочно-непрерывное решение уравнения (3), удовлетворяющее условиям:

1) V M6V0C V,z (М) < в ;

-> grad z(M)

2) вектор U о = 1--—— является вектором внешней единичной нормали к гра-

I Cj Т*СХ> СХ> Z (У*/ ) I

нице области

П = {Me V | z (М) < t }, t < в ,

тогда

z (М) = ty (М) , V М e Vf ( в ) , V0 = V.

Доказательство:

Введем в рассмотрение функцию:

R(M,T) = Н(М,Т)~ J p(N)rnjI(N)G(M,N,T ~ z(N))l(T - z(N))dVN>M Е V,

v

Поскольку z(M) есть решение (3), то R( М , z ( М ) ) = и^, М e V0.

Проведя преобразования, обратные осуществленным при выводе (3), получим, что функция /?(М, т) необходимым образом является решением краевой задачи:

рeR + ргпл<5 (т - z (М) ) = di v(Я flfrad R) + f (М, т) , М e V, т > 0, (4а)

Яр

+ £ R = q ,М e 5,т > 0 , (4б)

R (М,0 ) = 0, | R ( М , т) | < const < оо , М e V, (4в)

R(M,Z(M)) = uflME

откуда

R(M,t) = Н(М, т) — J | p(N)ran(N)G(M,N,T - ~ z(N))dVNd$,M Е V,t>0.

0 к

Учитывая известное тождество:

£4(т) f (iV,т) d Vv = 4 О,(т) f (iV,т) Vn (Vт) d5VV + I(T) £f (iV,т) d vv

где - нормальная составляющая вектора скорости движения границы области

ш (т) в точке # в момент времени т (нормаль внешняя к ш ( т) ), находим, что

R (М,т) = Я ( М ,т) - J0T/5z (т)р (VV) r л(V) G( М, VV,т - О vn (JV,<f)dS^f^ e 7, (5)

Где

Sz(t) = дй)(т\й)(т) ={M E V\z(M) < т},

так как на поверхности Б* (г) :

2 (М) = г, то Я(М, г) = ирМ Е Б2(г) С V0,г < 0. (6)

Учитывая известные свойства теплового потенциала простого слоя (каковым является второе слагаемое в (5)) и непрерывную дифференцируемость Н(М,г), получаем, что при переходе через поверхность Б* (г) нормальная производная от Я (М , г) терпит скачок первого рода:

п дя+ л дя~

= Р ^ МЕБ* , (7)

где п2- вектор единичной нормали к Б*(г), внешней по отношению к ш(г) .

Равенства (4), (6), (7) означают, что для г < 0 функция Я (М , г) удовлетворяет задаче Стефана (1) и в силу единственности ее решения

д(М,т) = и{М, т),М е У,т < в.

Поскольку наименьший корень уравнения и(М,г + 0)= и+0, М Е V, единственный, то 2 (М) = г г (М) , V М Е V г ( 0 ).

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

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

3.1 Метод ячеек

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

В пределах У| : гг(М) « г^сопБ^ р(М)гил (М) « р ¡гпл],Н(М,г) « ^ (г). Тогда получаем следующую систему нелинейных уравнений относительно неизвестных (¿ь...,

1%г-^хц(и - ^ 1 (г - ^ = Н(и) -щЛ = 1.....N; (8)

(О = I е

V]

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

Построим конструктивное доказательство возможности построения для задачи (1)

корректного приближенного метода, основанного решении на системы вида:

/

Ч ' ц'

1 %!Гвл]7с,Уа(Ъ - 1 (и - ъ) = Н ( 1д -игл = 1.....N. (9)

Для этого введем в рассмотрение схему метода прямых, аппроксимирующую краевую задачу (4) в [1]:

рэд = (Л Н]Л/Н) 1+р^I[ 1 (ио г - щ) - 1 (щ - щ)] Л Е Vй, (10)

В*(м/*) = 0Л ЕБ* = дV11,

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

Щ = Щ + Н(Ми 0. "о« = щШд-Возможность построения такой схемы прямых следует из известных результатов [1]. Решения (10) удовлетворяют условиям:

I I ^ I I г2(^)<со = с 0 п 5 £,

при

К = тах1еуп ттуе^|М,Му| -> 0, ||шЛ - Р^Н,^ -> 0, где - оператор сужения, ,

I I I I Л)=^ /о'Е у е ^2( £) | V/ Iй £ I Ъ I = ^ I/,

Для равномерной сетки | V/ | = /гт. Соотношения (10) перепишем в виде:

V* = А^" + Р\Р? = гпл ¿/^ [ 1 (ио; - иг) - 1 (и г - иг) ] , / е V", (11)

где

Пусть у^ у ( £) - коэффициенты нормальной фундаментальной матрицы У[{] , удовлетворяющей матричной задаче Коши: У = А * У, У [ 0 ] = / , где - единичная матрица:

По формуле Коши:

После дифференцирования получаем:

= ^ [1(и0| - "/) - - "/)] + Еуе^Т*/о - т)[1(«у(т) - и,) -

1 (и 0 у - иу) ] йт, / е V*, (12)

Введем в рассмотрение вектор , соответствующий выше определенной функции

£/ (М ): £* = ( £х.....£„) ,

где - наименьший неотрицательный корень уравнения

+ 0) = иу + ОД > 0. Если для какого-то г такого корня не существует, то полагаем £ ^ = оо. Тогда (1 2) переписываются в виде:

"Уц

м

Wi(t) = -ZjeуНГл]7cjytj(t - tj) 1 (t - tj),i EVh,t< eh. (13)

где [0, ]- интервал времени, для которого выполнено свойство:

Vf (t±) Q Vf (t2), 0 <t±<t2< Oh. (14)

V/ (t) = {iEVh \ и (t)>uf}.

При (13) обращается в систему нелинейных уравнений относительно :

Z J ev 7 V n (ti - tj) 1 (ti - tj) = Hi (ti) -ur,ie Vh. (15)

3.2 Исследование численного метода

Существование и единственность решения системы уравнений (15) в области показывается аналогично тому, как это было сделано ранее для функции t f (M) . Сходимость th к t f (M) оценивается ниже.

Утверждение 2. Пусть существуют независящие от h < h0 числа , такие что

Zjeve1 { ftj 1\ Uj__Uf \ у + V(мj) \ u(M ._)_Uf \ У} \Vj \ <K

где подмножество vq Q Vh, определяется из условия:

max{tittf(Mi)} < 9lt

IEV01 '

Тогда

\im\\th - Phtf\\ = 0,p > 1.

iP(v01)

Доказательство: Так как

Si1 (r-ti) dr = (t-ti) 1 (t-ti) ,

следовательно, при

tf(Mi) -ti = /Ц1 [ 1 ( ui(r) - Uf) - 1 ( u(Mi,r) - uf)] dr,

откуда

\ tf (Mi) -ti \ P = {/o 1 \ 1 (u (r) - uf) - 1 (u (Mi,r) - uf) \ drf, p>1. в силу неравенства:

j^-|l(uf(T) - Uf) - l(u(Mi,T) - Uf)I drj 1 fdl

^ Q- l^ico - uf) - lo(m;,t) - Uf) | dr, t>l J о

получаем, что

По свойствам единичной функции:

^ I £л-Рл£у I I 1/Р „ р>1 .

У

|1(х + у)-1(х)| < - , г > О,{х,у} с Е1.

х

Поэтому

г*г(М{)

4

щ{т)-и{Мьт)

щ{т) - щ

ь рб 1

с1т + г

щ{т)-и{Мьт)

и(Мьт)-иг

йт

Следовательно, с учетом неравенства Коши-Буняковского, получаем неравенство:

" Т"1г(?в1)

й- ЧМ/'Т)Г атЩ}1/3 С\мМмрт) -

иг\ +

1/5

- и(М;,г)Г5 йЦЯ) {Е,-^ /^Ь« -

иу I -""(йт I I },«,(?> 0,1 +1 = 1 .

у 2

Выбираем Г ^ = у, Г = —у, тогда Я = 1 + -, Г5 = 2 ,

1 + 2 У

ЬЬ-РНг II < ^2/(2+^11/1 _ р/1 11 г\(ув1) ~ И

| Со/Р _ ^ К<УнУС° " 5'

а так как

и при

ип _ Рпи = ]Л,П _ Рп]л/

Л — 0: I I и* - Р'и I I , , — 0,

:

то утверждение доказано.

Условие Утверждения 2 упрощается, если не рассматривать однофазные задачи Стефана. Тогда постоянная К может быть выбрана из условия:

У [ Г7-Т-¡у№|</Г,0<у<1.

Оценка величины I I £* - Р*£у I I ^ ^ ^ по форме остается прежней.

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

Если размерность систем (8), (15) бесконечна, число непересекающихся связных подмножеств конечно и не возрастает во времени, то эффективным может быть следующий метод решения. Исходным является следующее свойство:

VI е V/1,) е дУр-. и < ц, - = 0.

Обозначим для краткости через , тогда алгоритм примет вид:

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

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

^1(^1) = = Нч(.гч) = иГ'г11 = = гч =

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

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

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

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

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

6. Вычисляем

тах> Р — 1<

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

%0р) - иГ - X ГпЛ7/су */руОр - - = 0,

}ЕО

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

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

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

10. Полагаем х* = тт д ахр.

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

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

ш = [)р Е до| \ хр - х*\ < е}.

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

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

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

ггг < ■ ■ ■ <гг < ■ ■ ■

11 - - 'п -

Таким образом, разыскиваемое решение автоматически удовлетворяет условию: г г < 0п, хотя значение 0 априорно неизвестно.

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

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

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

5. Вычисление сингулярного слагаемого в суммах (8) и (15)

Слагаемые хг у ( г - гу),угу( г - гу) при / = ] являются сингулярными, так как функция С (Мг,Мг,г) неограниченна в окрестности нуля г. При этом вклад данного слагаемого относительно значителен. По этим причинам существует целесообразность вычисления сингулярного слагаемого аналитическими методами. Для этого в соответствующей ячейке У функция г у ( М ) аппроксимируется линейной частью ряда Тейлора:

гг(М) « гг(М г) + к ■ МШ[.

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

Уи = I со (хих, гг{х) - 1 (гг{х) - йх,

хь 2

к к х е [х1--,х1 + -],

где (согласно асимптотике для области сингулярности)

С о(М,Ы,г) « (2 л/паТ)~техр(- г^ /(4аг)). Для одномерного случая: ,

Х:+П/2

ехр

Уи

(х - х{)

2

Л/2

Г 4ак(х-хЛ 1(к(х_ ^ ах = Г ^[¿А ^ = ^

3 . 2л лак(х-х;) ) 2л1пакх

ехР. I к

2т1пак(х - X/) ./ 2л!пакх \J8ak

Х(-й/2 М у и п \ \

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

2пак

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

£/(*!)-£/(>¿-1) к~-л-"

6. Двумерная задача

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

Решение соответствующей краевой задачи при нулевой теплоте фазового перехода находится аналитически и имеет вид:

Г Г Г PoS(xN,yN)

H{M,t)= u0(N)G0(M,N,t)dVN + G0(M,N,t-x)dVNdx,

E2 ° E2 PC

где

схр( ^М-^2НУМ-УЫ)2)

g o(M,iv,t) =—^---,a = Я / (p с ) ,

S (xw,yw) — Д вум e p н ая дельта-функция,

Po — мощность нагревательного элемента на единицу толщины слоя. Если начальное распределение температур постоянно:

ио(.Ю = ио = const, то выражение функции Я (х, у, t) равно:

2 _l,,2n

где £1(х) - интегральная экспоненциальная функция. Компоненты матрицы (см. [1]) равны:

= <(0<22(0,4(0 = V2

Расчет осуществляется при помощи решения системы (8), (15).

h

erf' 2Vat / еГД 2Vat

Вычисление сингулярного слагаемого в (8), (15)

В [1] данная задача была решена для одномерного случая. Рассмотрим с целью вычисления сингулярного члена для пространственной области следующий интеграл

ехр (--(Х ~ Х*)2 + ~ Уй

/ = [ —-—^М* + 1 (гх{х - х{) + гу{у - у4)) йхйу,

I 4па[гх(х - хд + гу(у - уд] ^ у >

в котором функция ^ (х,у) заменена на линейную часть ее ряда Тейлора с центром в точке (Х1,У1), гх = д-^(хьуд, гу = д-^(х1,у1), (х,у) е У^.

Предполагая, что ячейка является квадратом с центром в точке (х^у^), находим:

/ х2+у2 \

Д = / Н/2 / Н/2 ех? \ -4 а(txX+tyy)) г х + у

Д )—н/2 3-Н/2 4 тта«Хх+1уу) 1 ^ ^ Уу) у'

Разобьем интеграл на две части: / = /0 + Д. Интеграл /0 вычисляется по области круга, вписанного в ячейку (рис.1):

г2п ГН/2 еХМ 4а|У Сатиру 2 г—\ ( Н \1

Д0 = /0 /о 4па | V ^ | 5 ш<р 1 (5 Ы(Р ) аЫу=2^[ 1-^хР{ Г ( ,

I V гг \ = у[% + Ц,

Этот круг делится на две части прямой ^х + Ьуу = 0 (на рис.1 изображена голубым цветом). Вклад в интеграл дает только закрашенная половина круга, в которой ^х + гуу > 0.

Интеграл Д вычисляется по области ячейки вне круга, которая на рис.1 закрашена сиреневым цветом. Эта часть также дает вклад в общий интеграл, так как в этой области также I хх + Ьуу > 0. Форма соответствующей области сложна для точного вычисления. Используем тот факт, что участки данной области относительно близко прилегают к соответствующей полуокружности. Для вычисления интеграла Д используем приближенную квадратурную формулу:

(1^2 тсЪ?'\

У ~ ~8~)~ + ~ ^^

где {С) - среднее значение функции С на соответствующей полуокружности, а сомножитель перед выражением в угловых скобках есть площадь рассматриваемой области, закрашенной на рис.1. в сиреневый цвет. В результате выражение интеграла Д примет вид:

и

л ехр —

/1 )Г\ И г- \ ва\4 1{\зт<р] ^ ~ \2 8/ п2а\Ч ■'о зтср ^

Рис.1. Разбиение ячейки на области Полученные одномерные интегралы эффективно вычисляются по формуле трапеций.

Аналитическое решение задачи Стефана

С учетом осевой симметрии данных задачи интегральное уравнение (3) в рассматриваемом двумерном случае упрощается и имеет вид:

/о

ГГдг

и1

2а(гг(г) -

= Н (г, £/0")) — — Г —-—-—-

У ; С1 2а (гДг)-^

■ш)

■Гд г

где Н (г, £) = и0 + ' 'о 00 — модифицированная формула Бесселя 1-го ро-

да.

Решение интегрального уравнения разыскиваем в виде: , где - неиз-

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

иг = Я(1Д) -

№л

2 Лк

I

Г (2ак(1 - х2)) ехр ( 4ак(1 — х2))

1-х2

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

Найденное точное аналитическое решение целесообразно использовать в качестве контрольного варианта для проверки вычислительных качеств вычислительного метода (15).

Результаты вычислений

Вычисления осуществлялись при шаге 10 см, Р0 = 1 0 0 Вт/м, и 0 = — 5 0С. Результаты численного решения представлены на рис. 2 и в табл. 1. Погрешность численного решения

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

и ОО II О____

о о о о о о

О О О О О О О О

ооооооооо оХ оооооооооо о\ ооооооооооо оооооооооооо ^

оооооооооооооо\

сссссссссссссс] СССССССССССССС1

ссссссссссссссс ссссссссссссссс

СССССССССССССС! СССССССССССССС!

оооооооооооооо/ оооооооооооооо/ ооооооооооооо/ оооооооооооо у ооооооооооо о >6 ооооооооооо// ссссососсс^г оооооооо О^о

О О О О О О ■Л д д ------

Рис. 2. Численное и аналитическое решение задачи Стефана (188 и 540 точек)

Таблица 1

Количество точек Численное решение Точное решение Погрешность, %

188 23417,90 24266,80 3,5

240 29967,62 30903,88 3,0

540 6886078,64 6989668,75 1,5

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

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

Вывод

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

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

1. Арутюнян Р.В. Интегральные уравнения задачи Стефана и их приложение при моделировании оттаивания грунта // Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 10. С. 419-437. DOI: 10.7463/1015.0814769

2. Каменомостская С.Л. О задаче Стефана // Математический сборник. 1961. Т. 53. № 4. С. 489-514.

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

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

5. Данилюк И.И. О задаче Стефана. УМН. 1985. Т. 40. Вып. 5(245). С. 133-185.

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

7. Taras A. Dauzhenka and Igor A. Gishkeluk. 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 2013. Vol. 1, WCE 2013, July 3-5, 2013, London, U.K.

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

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

10. Nekrasov S.A. The integral equations of the Stefan problem. II // Differential Equations. 1996. Vol. 32. No. 9. Pp. 1254-1258.

11. Nekrasov S.A. The integral equations of the Stefan problem. I // Differential Equations. 1996. Vol. 32. No. 8. Pp. 1118-1125.

12. Nekrasov S.A. Modeling of phase transitions of the first kind by the method of integral equations in the case of a stationary moving surface source // Journal of Engineering Physics and Thermophysics. 1994. Vol. 66. No. 6. Pp. 674-677. DOI: 10.1007/BF00867972

Science ¿Education

of the Bauinan MSTU

Science and Education of the Bauman MSTU, 2016, no. 08, pp. 165-181.

DOI: 10.7463/0816.0843482

Received: 08.07.2016

Revised: 22.07.2016

© Bauman Moscow State Technical Unversity

Boundary Integral Equations for Multidimensional Stefan Problem and Their Numerical Solution

R.V. Arutjunjan1'* "robSTigmaiLru

:Moscow Technical University of Communications and Informatics,

Moscow, Russia

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

The article offers a new numerical method for solving the multidimensional Stefan problem based on the reduction of the initial boundary value problem for the equivalent nonlinear integral equation of the minimum dimension. A feature of the method is to use the time of phase transition as a major unknown function. This parameterization is an efficient computational technique to have a greatly simplified method of solution. A computing quality of the method is investigated in the context of solving the practically important problems of thawing the frozen soils with different methods of thermal influence. The method advantages are most evident, when solving the problems in an unbounded domain, as apart from the traditional finite-difference methods there is no need to truncate a system of equations that have a large, possibly infinite, dimension.

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

The method also provides additional comfort in Richardson extrapolation, when solving the inverse Stephen problems, in particular the problems of the phase transition front control.

The article presents a detailed study of the considered problem correctness and the proposed method of its solution. Gives a strict mathematical proof of the assertions that there is the only solution available for the formulated multivariate nonlinear integral equation of the Stefan problem.

There is a proved statement on the convergence of the proposed numerical method of cells for the approximate solution of integral equations of the Stefan problem.

The article focuses on the singularity of the integral equation of the boundary value problem under consideration. Provides the accurate and practically acceptable asymptotic estimates of singular members of nonlinear integral equation of the Stefan problem.

The test two-dimensional boundary Stefan problem describing the process of thawing soil has been solved. The simulation results prove high quality of the computational method.

References

1. Arutyunyan R.V. Integral Equations of the Stefan Problem and Their Application in Modeling of Thawing Soil. Nauka i obrazovanie. MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2015, no. 10, pp. 419-437. (in Russian).

DOI: 10.7463/1015.0814769

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

3. Meyrmanov A.M. Zadacha Stefana [Stefan's problem]. Novosibirsk, Nauka Publ. (Sibirskoe otdelenie), 1986. 187 p. (in Russian).

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

5. Danilyuk I.I. On the Stefan problem. UspekhiMatematicheskikh Nauk, 1985, vol. 40, no. 5(245), pp. 133-185. (in Russian). (English version of journal: Russian Mathematical Surveys, 1985, vol. 40, no. 5, pp. 157-223. (in Russian).

DOI: 10.1070/RM1985v040n05ABEH003684 )

6. Rubinshteyn 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).

7. Taras A. Dauzhenka and Igor A. Gishkeluk. 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 2013. Vol. 1, WCE 2013, July 3-5, 2013, London, U.K.

8. Samarskiy 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. (in Russian). (English version of journal: USSR Computational Mathematics and Mathematical Physics, 1965, vol. 5, no. 5, pp. 43-58.

DOI: 10.1016/0041 -5553(65)90004-2 )

9. Budak B.M., Solov'yeva E.N., Uspenskiy 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. (in Russian). (English version of journal: USSR Computational Mathematics and Mathematical Physics, 1965, vol. 5, no. 5, pp. 59-76. DOI: 10.1016/0041-5553(65)90005-4 )

10. Nekrasov S.A. The integral equations of the Stefan problem. II. Differential Equations. 1996. Vol. 32. No. 9. Pp. 1254-1258.

11. Nekrasov S.A. The integral equations of the Stefan problem. I. Differential Equations. 1996. Vol. 32. No. 8. Pp. 1118-1125.

12. Nekrasov S.A. Modeling of phase transitions of the first kind by the method of integral equations in the case of a stationary moving surface source. Journal of Engineering Physics and Thermophysics. 1994. Vol. 66. No. 6. Pp. 674-677. DOI: 10.1007/BF00867972