Численное решение одномерной плоской задачи Стефана Текст научной статьи по специальности «Математика»

УДК 519.87

Н.Н. ЛАПИНА, В.Н. ПУШКИН

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОМЕРНОЙ ПЛОСКОЙ ЗАДАЧИ СТЕФАНА

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

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

Введение. Рассматривается одномерная модель распространения плоского фронта фазового превращения в твердой среде. Данная задача имеет множество приложений. В частности, она может быть рассмотрена в связи с поддержанием в надлежащем тепловом состоянии приустьевой зоны скважин, эксплуатирующихся в условиях вечной мерзлоты [1]. Модель основана на следующих представлениях. В каждый момент времени физическое пространство разделено на две зоны, различающиеся тем, что заполняющее их одно и то же твердое вещество находится в разном агрегатном состоянии и, соответственно, обладает разными теплофизическими свойствами. Последние полагаются однородными по пространству и неизменными во времени для каждой из зон. Область фазового превращения вещества из одного состояния в другое имеет малую толщину, что позволяет считать её поверхностью. При переходе через эту поверхность свойства среды меняются «скачком».

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

С

где С:

м 1, п

н С2Р 2 о с1р1

Уф

пРи У< Уф пРи У> Уф

ІЇ,

д0_

дт

м1,

- п 3 ■ н 3

О 1

= л

д 20 д у 2 '

при

—, при

У < Уф У > Уф

У:

(1)

г

И '

ИЪ 1р.

Здесь z - пространственная координата; t - время; И - продольный размер физической области; Т - температура в точке с координатой г в момент времени £ Т1, Т2 - масштабные значения температуры, первое из которых соответствует ее значению при г=0, а второе - при г=Ь; с , I , р - соответственно удельная теплоемкость, коэффициент теплопроводности и плотность вещества; индекс ф относится к величинам на границе раздела между зонами; 1 - к левой зоне, 2 -к правой; гф - координата фронта фазового превращения.

Кроме того, модель включает граничные условия первого рода:

при у = 0: 0 = 0, (2)

при у = 1: 0 = 1, (3)

что соответствует тому, что на одной из границ на протяжении всего периода наблюдения поддерживается температура Т1, а на другой - Т2.

Начальное распределение температуры полагается однородным:

0

л

1

т

при т = 0, У е (0;1] : 0 = 1. (4)

Дополнительно на "межфазной" границе требуется выполнение условий непрерывности температуры и теплового баланса, то есть при т > 0, у = Уф :

0 (Уф - 0) =0 (Уф + 0) =0 ф ;

э©

где V

ёуф

ёт

0

эУ

Т - Т І± її

Т - Т '

-л

Ч:

э©

эУ

Чф

С |Т2 - Т|'

Т

ТФ

(5)

(6)

температура фазового перехода;

qф - теплота фазового перехода в расчете на единицу массы вещества в состоянии 1;

V - безразмерная скорость перемещения фронта фазового превращения.

Заметим, что задача (1) - (6) соответствует процессу фазового превращения вещества, первоначально находившегося в однородном состоянии с температурой Т2, когда на левой границе, начиная с начального момента времени и в дальнейшем, устанавливается температура Т1. Подразумевается, что Т0 < Тф < Т1.

Из (4) в силу очевидного неравенства 0 <0 Ф < 1 следует, что начальная координата фронта Уф(0) = 0 .

Для расчета начальной стадии процесса, а также для тестирования численного алгоритма применим аналитическое решение задачи для полубесконечного пространства.

Для малых значений переменной времени Т влияние правого граничного условия (3) несущественно для развития процесса. Поэтому заменим его условием на бесконечном расстоянии:

при У = ¥ : 0 = 1. (7)

Будем искать решение задачи (1),(2),(4)-(7) в виде 0 (т; у) = I(п) , где П =

Уравнение (1) удовлетворяет функция:

м V

.2 ц

п Ь у ехр з--ч du + Ьфри

п о и 4 ш

Ж Си2 ц

(8)

В к

и 4Л ш

к

где Лф = к, соответственно, Уф = к4Т; V = 2т . Чтобы выполнялись граничные условия и

условия на границе раздела "фаз" (при V ф = к ), константы в (8) должны определяться из следующих выражений:

Ь2 = 0;

1 +

кч

2ЛГ

|ехР

к

С (к2 - и2) 4Л

ёи

Л

ехр

к2_ 4

= — ехр

1 Л

Ск2 4Л

ехр

ехр

С (к2 - и2)

4Л

к

ёи + | ехр

о

и

4

ёи

к^_

4

»,-

-о

+о

ф

1

к ( h2 = Ь11 ехр

и

4

2 \

du

Как видно, приведенные соотношения зависят от безразмерной координаты границы раздела Лф = к . Последняя, в свою очередь, определяется из трансцендентного уравнения:

^ =0 ф. (9)

В качестве модельных примем следующие значения входящих в задачу констант: 12 = 1,65Вт/(м • К) ; А1 = 1,32Вт/(м • К) ; с2 =1220 Дж/(кг • К) ; с1 =1640 Дж/(кг • К) ;

р 2 =1780 кг/м3 ; р1 =1800 кг/м3 ; qф =66740 Дж/кг ; h =10 м ; 7}=+6°С; Т = -2°С; Тф =0°С.

Решением уравнения (9) для базового набора параметров является к = 0,5003 . Соответственно, зависимость от времени координаты фронта фазового превращения в задаче с полубес-

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

Уф = 0,5003лУТ . (10)

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

а,в(к +1) = ОРС11) + dIe«k;1) + Ь,

где

= 1,2,

с =

1

У и 1 - У г

А

У ,+1- У,

(11)

при г < ^

/

при г > ^

1

У г - У г -1

_А_______

У г - У г ’

п с. + йпри1—г-^з;

П А т

а г = Н -

п с ■ + d. + СМри Гг1 , з О Ат

0 (к) Г- 1

Ат

п 0 (к) , з

пО А т

при г < з

С(гф У(к + 1))ф+ У'

,(к +1)

зф

Ат

Ь = 0 (к) Ч

С (г* У(к+ 1))ф+ У(

-1

Ат

+ qV

(к) .

Уф,1 = У<к) , V»)ЧАТ ; Ук) = Уф{Тк) ; V« = V{тк) ; т, = тк + Ат ; з - индекс контрольного объема, в котором происходит фазовое превращение.

<

>

>

г

з

з

Как видно, кроме температур система (11) включает неизвестное значение V(к) скорости фронта фазового превращения. Для решения полученной системы применялся метод сквозной «прогонки» в сочетании с методом «пристрелки» в отношении величины V(к), что позволя-

лл-ъ. п (к +1)

ет находить попутно очередную координату поверхности раздела фаз уф .

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

0 ф = °,75. Последнее из представленных на рис.1 распределений температуры близко к стационарному распределению, когда скорость распространения поверхности раздела "фаз" равна нулю. Для того чтобы рассчитать положение уфс) "стационарного" фронта достаточно в уравнении (1) производные по времени положить равными нулю и воспользоваться условиями теплового балан-

са (5), (б). В результате получаем соотношение Уф -

е

е ф іЛ(1 -е ф) •

Ґ

-1 г2 / 3 — 4

-Г”! ■ <— —1—1—1— —1—1—г —1—1—1— —.—1—1— 1 г —і—.—,— —і—1—1— —і—і—і— Г—1

0,8

0.7

® 0.6 га

Cl

?

І 0.5

С

s

Р

о.з

0,2

0,4 0,5 0,6

координата, у

а

y

Рис.1. Распределение температуры 0 в различные моменты времени т :

1 - т = 0,0002 ; 2 - т = 0,0777 ; 3 - т = 0,6557 ; 4 - т = 11,71

с 12

I я базовых параметров это дает Уфс) = — » 0,706. Численный расчет, что в некоторой

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

рис.2 представлены зависимости от безразмерного времени т безразмерных скорости V фрон.а фазового превращения и координаты самого фронта Уф , полученные в результате численного решения задачи (1)-(6). Здесь же в виде набора точек показано определяемое из соотношения (10) изменение с тече—... "р£00£динатау—'— границы раздела между зонами для

задачи с бесконечной протяженностью физического пространства. Как показывают расчеты, влияние правого граничного условия в задаче (1)-(6) несущественно, пока фронт фазового перехода Уф не достигает значений близких к 0,2 (соответственно, пока переменная времени Т не больше 0,16). Об этом свидетельствуют, в частности, почти идентичные распределения температур на каждом шаге по времени в совместной расчетной области для каждой из задач.

Время, х

Рис. 2. Изменение поверхности раздела между зонами с различным агрегатным состоянием вещества: 1, 2 - соответственно графики скорости перемещения фронта и координаты фронта; 3 - график изменения с течением времени положения границы раздела "фаз"

Эти распределения определяются по формулам (8). Кроме того, при численном решении задачи (1)-(6) отслеживались изменения величин уф/и 2V•^/X , которые в задаче с бесконечной толщиной грунта представляют собой константу k . Оценивать относительную величину разницы между значениями Уф , полученными для каждой из двух моделей, не имеет смысла, потому что, во-первых, для малых значений времени малы и значения Уф , и, в силу того, что ЭВМ работает с ограниченным числом разрядов, при малых X относительная разность может быть довольно большой, а, во-вторых, начальная скорость перемещения фронта в силу модельного представления о разрывности начального распределения температур должна быть бесконечной, чего также невозможно добиться при численном решении с помощью ЭВМ. С учетом указанных обстоятельств, для проверки близости решений разумным представляется соотнесение величин ку = уф14Х и ку = 2V , получаемых в ходе численного решения задачи с конечным продольным размером физического пространства, со значением k , определяемым из (9). Заключение. Расчеты показывают, что на начальном временном промежутке относительные разницы между к у и к , а также ^ и к , не превышают 10% от к , причем чем больше X отличается от нуля, тем меньшей становится эта разница, и лишь тогда, когда начинает существенно сказываться условие на правой границе задачи (1)-(6) (соответствующие X намного

больше значения, которое соответствует концу начального промежутка), она приобретает явную тенденцию к возрастанию. Заметим, что при неограниченном возрастании переменной Т кривая

2 (см.рис.2) выходит на горизонтальную асимптоту у = уф (для базовых параметров у = 0,706 ). Последнее означает, что положение фронта стабилизируется. При этом кривая 3 асимптоты не имеет, то есть координата поверхности раздела зон в случае бесконечной расчетной области не имеет верхнего предела.

Библиографический список

1. Быков И.Ю. Математическая модель охлаждения мерзлых пород приустьевой зоны скважин в условиях естественной воздушной конвекции / И.Ю. Быков, В.Н. Пушкин, В.Н. Емельянов // Строительство нефтяных и газовых скважин на суше и на море. - 2008. - №9. - С.10-14.

2. Патанкар С.В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах / С.В. Патанкар. - М.: Изд-во МЭИ, 2003.

Материал поступил в редакцию 10.11.09.

N.N. LAPINA, V.N. PUSHKIN THE NUMERICAL SOLUTION OF ONE-DIMENSIONAL FLAT STEFAN'S PROBLEM

The approach to modeling and numerical solution of the problem on distribution of the phase transformation front in the frozen breeds which is quite actual for drilling and operation of chinks under the permafrost conditions is considered. The exact solution of the problem in case of semi-infinite physical space is presented. The area of its applicability is defined by means of a computing experiment.

ЛАПИНА Наталья Николаевна, старший преподаватель кафедры АИС Ухтинского государственного технического университета. Окончила УГТУ (2005).

Область научных интересов: математическое моделирование и численные методы решения задач тепломассообмена, возникающих в нефтегазовом деле Автор 13 научных работ.

Natalya-Lapina@mail.ru; nlapina@ugtu.net

ПУШКИН Виктор Наркистович (р. 1954), заведующий кафедрой АИС Ухтинского государственного технического университета, доцент (2002), кандидат физ.-мат.наук (2000). Окончил УГТУ (1977).

Научные интересы: математическое моделирование и численные методы решения задач гидродинамики и тепломассообмена.

Автор 50 опубликованных работ.

vpushkin@ugtu.net