Математическое моделирование процесса теплообмена с фазовым переходом

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ

УДК 536.421.4

Р. Г. ЗАЙНУЛЛИН

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ТЕПЛООБМЕНА С ФАЗОВЫМ ПЕРЕХОДОМ

Излагается применение метода вырожденных гипергеометрических преобразований (ВГГП) к решению одной нестационарной задачи теплообмена с фазовым переходом на примере процесса промерзания некоторой сплошной среды. Получено аналитическое решение задачи при специальных краевых условиях. В ходе решения устанавливается закон движения границы раздела двух фаз. Рассматривается случай, когда па уравнение для ядра интегрального преобразования, используемого при решении задачи, накладывается требование самосопряженности, что дает возможность применения метода разложения по собственным функциям. Фазовый переход; свободные границы; задача Штурма-Лиувилля

Краевые задачи теплопроводности в областях с подвижными границами за последние годы приобретают все большее значение как в теоретических, так и в прикладных разделах физики и математики. Количество работ, посвященных решению подобных задач, за последние годы стало заметно увеличиваться [1-9]. Аналитический подход при решении краевых задач теплообмена в системах со свободными границами относится к числу труднейших проблем в современной аналитической теории математической физики. Вследствие зависимости положения характеристического раздела области от времени к этому классу задач неприменимы классические методы дифференциальных уравнений в частных производных, так как оставаясь в рамках этих методов, не удается согласовать решение уравнения теплопроводности с движением границы фазового перехода.

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

Контактная информация: (347) 273-77-35

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

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

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

ПОСТАНОВКА ЗАДАЧИ

Рассмотрим применение метода вырожденных гипергеометрических преобразований (ВГГП), разработанного в [10, 11], к процессу теплообмена с фазовым переходом при затвердевании некоторой сплошной среды под действием плоского источника холода. В начальный момент времени среда обладает постоянной температурой t0 > 0 при х > а£,0. Внешняя поверхность среды обладает температурой te < 0. Образуются зоны промерзания (к = 1) и охлаждения (к = 2), и граница фазового перехода с течением времени продвигается внутрь среды. Математическую модель этого процесса для одномерной схемы можно представить в виде:

дtk

Эх

= ак

х є Бк (т): Д (т) =Хо < х < ^ (т)},

В2 (т) = (4і(т) < X < 52(т)};

т > о, £і(+0) = ^о > о, ^2(т) = о^і(т), а >1;

*1 (х,0) =

V

*2 (Х,0) =

^1

(а -1)5 0 *1(0,т)=*е;

*к(51(т),т)=0;

*2(52(т),т) = *0;

Э*1

= ОУ-

Эх Эх ёт

Здесь ак и \к - коэффициенты температуропроводности и теплопроводности в Дк(т), О -скрытая теплота кристаллизации, а V - плотность образующейся фазы, а - безразмерный параметр теплового влияния [12].

Предполагаем, что їк(х, т) и 51(т) являются решениями задачи (1-7) для всех 0 < т < Т < да, если:

1)

Э 2*к

Эх2

при 0 < т < Т; 2) 4(х,т) и

и

Эт

непрерывны в Дк

Э*к

Эх

непрерывны в

Дк(т) при

0<т<Т;

3) tk(x, т) непрерывны также при т = 0;

4) £,1(т) непрерывно дифференцируема при 0 < т < Т;

5) они удовлетворяют уравнениям (1-7).

НАХОЖДЕНИЕ СОБСТВЕННЫХ ФУНКЦИИ И ЗАКОНА ДВИЖЕНИЯ СВОБОДНОЙ ГРАНИЦЫ

Используя подстановку

х = >-51

и вводя вспомогательные функции

(1 - у)ґе при к = 1;

Тк(^ 5) = (к(У^Сй Ф - *

У -1

а -1

*0 при к = 2,

сведем задачу (1-6) к решению уравнения

52

= а,,

ЭТк

Эт

+ 51(т)5і (т)у

+

ЭУ

*е51(т)5^ (т)у при к =1;

'51(т)5(т )у при к = 2

1с однородными краевыми условиями

^,0) = 0;

Т,

Решение задачи (10-12) отыскиваем, используя метод конечных интегральных преобразований по у:

ик Сь т) = \Тк

У1

где у! = 0, у2 = 1 при к = 1 и у! = 1, у2 = а при

к = 2, с априори неизвестным ядром, полагая,

что постулируемые ниже свойства преобразования (13) имеют место равномерно относительно т.

Ядра Кк(у, у) преобразования (13) являются решениями уравнений

х

Э 2 [р к(У) Кк(У, т)]

к 3 2

ЭУ

Э[ур к( у ) Кк(У, т)]

- 51(т)51 (т) х

ЭУ

+ Ц 2,ур к(У) Кк(У, У) = 0

при однородных граничных условиях

К1

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

В дальнейшем, будем пользоваться формальными равенствами

Кк(У, У) = Кк ,у ( у); ц2( у) = ц 2,у; ик Сьт) = ик ,у (т),

*

а

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

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

акрк (У) = 51(т)5 (т)УРк (У).

Полагая, что

2^(т)£ (т) = Л2, получаем весовую функцию

Рк

V 4ак )

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

т

42(т) = | Л2^ = Л2т + 5 2.

0

Параметр Л подлежит определению через условие (7). Он определяется характером положения и движения свободной границы и является функцией тепло-физических характеристик и краевых условий. После наложения требования самосопряженности уравнение (14) можно представить в виде

_Э_

ЭУ

Рк(У )-

ЭКкл (у)'

^к ,у

эу

+

Рк

при однородных граничных условиях

К1у

Уравнение (16), при необходимости, можно привести к виду

э2 Кк,у

Эу2

+ 2 Рку~

эу

+

+

ц к,-

Кк,у (у) = 0.

Задача (16-17) есть задача Штурма-Лиу-вилля. Подстановками

Кк,у (у) =

22

Л

ик, у(г)ехР(-г =

Л2у

4а,

=Рку2

уравнения (14) приводятся к вырожденным ги-пергеометрическим:

Э 2ик,у

Эг2 " +12 ^ У Эг

- ьк и ,у(г) =0;

Ьк ,у = 1 -

V Л У

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

«1, У

подстановкой

Кк,у (У) =ик,у (У) х

х ехр(-в к I ^ = ик ,у(У) ехР^ -1 в кУ2) Л2

где вк =----уравнение (18) приводится к виду

4а

2

ц к,-

и!,у

}к\^к ,у

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

и1д

Решая задачу (19-20), находим собственные функции задачи (16-17) ортонормированные с весом рк(у), которые имеют вид

3

К1,у(у) = с- ур (, -2;г У ехр(- г); К 2,у (у ) = С2",1уАр | ь2,у, -2;г У ехр(- г),

где Ск'у - нормирующие множители, определяемые формулами

Ли + , 5 „ 1 („ I „ РГ

ф.,, + 1,5.р, Iехр(-р,)Е—-

с2 =

с1,у =

6Р1

п= (1,5 )„

X Е

п=1

где

*

с2 =

^ = 2В2 V Л2

(0,5)п

п (г 2п-1

а Р0 - Р

К

л

п!

п-1 1 п-1

\п = ^ = Е

т=0*к ,у + т *

у - 0,5 + т

В = -3*.яр| Ъ2л +1-;в2 | -

5

'ґ

- а| *2,, - 1 ) Рр| *2,1 + 5'^; в 21;

а

к

а

к

а

а

к

1

ак2

q = FI b2,у - ^2,:2; в2a2

= F|'2,y - ^,2; в 2 1, Po = FI '2,y ф в 2 |.

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

ОПРЕДЕЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ

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

33 Е

где

3 2 1 1 ,

АЕ (*2,т ,-; в 2 у2) = Е (42,у- -,-; в 2а2) х

3 2 3 2

х УЕ(Ъ2,1,2. в 2 У ) - аЕ(42,у ,2. в 2а ) х

1 1 2

х Е (42,т- --; в 2 у2).

Используя частные случаи вырожденной гипергеометрической функции [15]:

Е (4 у)=(пу)-0,524 - 2,5 ^ у х х Г(* - |)Щ-2Ь (^Т^у) - Д-2* (д/Зу )];

Е(4 -1,1; у) = п_0’52ь-1,5е0-5у х 2 2

х Г(4)[Д1-2* (^л/2у) + Д-24 (^2у )],

где Г(Ь) есть гамма-функция, в уравнениях (24) можно перейти к функциям параболического цилиндра

[ А-26^ ( л/2р1) -

- А-2^ ^л/2в7)]Г(Ь1>т - -2) = 0;

[^1-262 у

- А-262 „ (-йд/2Р 2 )А-262 „ (>/2Р 2 )] '

]x

г(*2,у - 2)Г(*2,у ) = 0.

При малых значениях безразмерного параметра вк способ нахождения собственных чисел спектральной задачи (16-17) представлен в статье [16]. Будем считать, что левая часть уравнения (21) порождается возмущенным дифференциальным оператором. Собственные значения задачи (21-22) с невозмущенным оператором хорошо известны и имеют вид

2

Мт,у

Для нахождения собственных значений Pk,y воспользуемся известным асимптомати-

ческим разложением

=(пг)2 + Е

2m

пу

m=0 у

2¥

+Е

d

2m

2m

a2 Ia 1У m=0 У

В работе [13] приведена система уравнений m Гm - 0,5 ^

Е

г=1

*

"2a,

=m-r I s=1

= ------1---- |2A2m Гm - ^ У V

У 2 - У, У, I 2 У

причем У, = 0, У2 = 1, ^ = C2as при k = 1 и у, =

= 1, У2 = a, c*a^ = d2a, при k = 2

, г ,

V У

m! Г

----------------, о = Е a „,

(m - г)!r! s=1

Ai(^ у )= Е

l

m=0

l+m

Bl, m (У ),

величины Bl? m(y) находятся из рекуррентных формул

Bi,m (у )=[(РкУ )2 + Pk \ Bi-2,m (y ) +

+ 2iBi-i,m-l (У )- Bi-2,m (У ),

Bi ,m (У ) = 0, если I < 0 или m < 0,

Bo,o(y)=1, i = v-!.

Заметим, что m-е уравнение системы (27) содержит коэффициенты с номерами не больше 2m-2, причем C*m-2 в первой степени. Это дает возможность последовательно найти все коэффициенты асимптотического разложения (26). Выпишем вычислительные формулы для нескольких первых коэффициентов:

1 - !

Co = jh (у)dУ, с2 = 4jg-2(у)dУ,

0 4 0

- 1 - 1

C4 = - j g-3 (У )dy+—j [g!(У)]2 ^,

= -—с2 +—У g,4 (y )dy +—x 4 2 64Jn 1V ' 32

a

0

C

6

1 -I

х j g! (у ) [g! (у )] dy+—j Lgtv )] dy

64

+

- a - a

d0 =-7 j h2 (У ^ d 2 =—(-ГТ j g22 (У )dУ,

a -1- 4(a -1) -

2

1 a 1 a d4=8(0-1) jg 3 (У )dv+16(0-1) f k2 (У)] dy’

5 5 a

= — d2 +---т j g 2 (У) dy +

' 4 2 64(a-1)J*2^"

2

aa

32(0-1) jg 2 Су h 2 (У)] dv+6^ j[g 2(у )]2 dv•

Здесь

hk(У ) = (P кУ)2 + Pk,

1 y 2

gk(У ) = hk(У)--------------j hk(У ^

-/ 2 У I у1

причем у! = 0, у2 = 1 при к = 1 и у! = 1, у2 = а при к = 2,откуда

&

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

^04в+в,, ^ = -1 в4, с4 =-^в6+-1-в4,

45 945 12

2 q^i 1 о^| 1 о 4

с6 =--------р, +-----р, +------р,

6 2835 18 16

d0 = 3 (j2 + a + 1) в2 + p,

d2 = -1 (a -1)2 I°2 + - j + О Р

945

2

d4 =т^->-1)4| a2 + 31a+1 |p62 + ^pl,

d 6 =-(a -1)4 х

2835

,4 313 з 189 2 313 , ,08

х I a +-------a +--------a +---------a +1 I в 2 +

128 64 128 2

+ — (a -1)21 a2 + - a +1 | p2 + — p2, 18 V 'I 4 Г2 16

Числовые значения этих коэффициентов можно найти только после вычисления параметра Л. Используя значения полученных ко-

эффициентов, асимптотический ряд (26) можно записать в следующем виде:

2

^ = М2- + с„ + % + C- + % + о(4 |, y €К,

ТУТ

a

щ

2

,, + d„+d2+d4+d6+О

a-1У У2 У4 У6

У е N;

При умеренных значениях pk собственные числа могут быть вычислены с помощью уравнений (22) или (23) на ЭВМ в программной системе Maple. При больших значениях параметра рк для нахождения собственных чисел 2

Цку спектральная задача (24-25) решается путем преобразования (24) к виду

uk,у

Ук,у = (3 - 4bk,y)Р-1 .

В теории асимптотических методов точку Ук, Y называют простой точкой поворота уравнения (28) или точкой перехода. Полагая, что pk достаточно большой параметр для нахождения собственных функций vk, Y(y) воспользуемся известным асимптоматическим разложением [14]:

uk,y (У) =

эр( y, pk1

Эу

+A,,ri (v Pk) + o[-L

an I оm+1

n=1 pk v pk

Mk,TAi

+ pk 3

£ B,n (V P-n ) +

on I о m+1

n=0 pk V pk

2

Mk,TAi!

где

A,n (y, pk)=1

х

Эру

Эу

х

эру (y, p-1) 1

У Cn j Vk ,y "2 э 2Bn-1

эУ ъе dt

By ,n (y,p -1) =--

py

эУ

X

2

2

0

X

+

X

X

X

2

2

у

X I

Ук ,у

5,

ЭУ

2 э 2 л^в-1) Эг2

А,0( у, в -1) = с,

Э5У

ЭУ

Л,-1(У,в-1) = 0 = Бъ-1(у, в-1), Сп = 0, Ау,2п+1 (У,в-1) = 0, БЪ2п (У,в-1) = 0, п е N

5у (У, в-') =

3 У

-1 2 3

V ук ,Г

2

1 ] 3

22 у -у

88П(У - Ук,Д

Соответствующие вычисления приводят к следующим соотношениям

2

5, (У, Р-1)=-[ 3 ]3 X

(

X

V

2 у / 2 2 У, агссо8-----------------л у - у у

к ,у У V к ,у ^ ^

к,

у

0 < У < У1,уЛ < У < У2,Y,

5, (у, в-1)=[ 3 ]3 X

Г

X

у^у2-уI -у2к,у!п:

У + л У - У

к ,у

У,

У1 у < у < 1, У2,у < у < a,

Э5Х(у,в-^ _ Г 4 ] 3 I ^2 „2

ЭУ

— Ку -у X

3 11 к,т

2 2 2 Ур - у*,т - у к,у1п"

22

У + ЛУ - У

1

] 3

У

2

ми= Б1

в13 5у

2

^у=-Л1

в13 5у Ь-1)

М 2,у = Б1

У

Г 2

в23 5у

л

^2,у = -Л1

в23 5у

2 ^у\1'Ы-’2

V У

]

1

Тогда собственные значения ц2 в асим-

к,

птотическом приближении (29) с точностью до порядка о(вк) являются решениями следующего характеристического уравнения

Г 2 . ] Г 2 . ]

Б1

в к15у

2

кк У V

] Г 2

-Г

= 0,

- Л1 вк15у

V у V у

где у1 = 0, у2 = 1 при к = 1, у1 = 1, у2 = а при к = = 2. При больших значениях аргументов в уравнении (30), особенно при к = 2, следует воспользоваться известными асимптоматиче-скими формулами для функций Эйри [17]:

* •/ \ 1 -0,5 -0,25 -£ / 1\п >— п

Л1(г)~- п , г ’ е ^ £ (-1) спС ,

2

Б1(г)~ п-0,5г-°’2У £ (-1)^-п.

Здесь

Г(3п + 0,5) 2 2 , ,

----- —7——^, С = -г2, г >> 1,

п!54п Г(п + 0,5) 3 11

< 1 п.

При к = 2: 1 < у2,у < а, г = в235у (а,в21), Л1(в22 5у (а, в 2-1))~1 п-0,5 г-°^25е-5 (1 - 0(С-1));

Б^5у(а,в2-1))~ п-0,5г-°’25е5(1 - о(с-1)),

и уравнение (30) в более подробных обозначениях принимает вид

2

]] 3

2Л1

ГГ

3

4

V V

агссо8

-^3 - 4Ь2,у- в 2) в

в2

3 - 46.

2,

У

2

V

У

2

2

п=0

п=0

Сп =

2

2

2

3

к,

2

( (

ехр

Яд/(°2Р2 - 3 + 4Ъ2,У )р2 - (3 - 4Ъ2,У )Х

х 1п

°л/р2 + д/°2Р2 - 3 + 4Ъ.

2,У

V V

р - 4Ъ2, ,

• Ві

((

3

4

V 'V

агоео8

-^3 - 4Ъ2,у- в2) в

в2

3 - 4Ъ

2,

2

ЛЛ 3

//

/у

= 0.

Для сохранения точности порядка о(в2'): о(с-1 )< ор1) достаточно выполнение условия 42 > в2 , т. е.

ал1а2 - У2, у - У2,у 1п

2,_ а + л/°2 - У22,у

> 2,

У2,

где У 2,; =

3 - 4Ъ.

2,

в2

ПОСТРОЕНИЕ ТЕМПЕРАТУРНОГО ПОЛЯ

Преобразованием (13) задача (10-12) приводится к виду

дик,у

к ,у \ Эх

- +

где

икл (0) = 0,

яи = ^р^\г5;ві |;

Нк ^ (т)’

и,У

7 бси V1,1 2

я,, =

ґоЛ2

2,

х I

п=0

6 (і - а)С

х

2,

(а2п+3 -1) д ^4" +

+ 3а (а2п -1)ра

РП

п!

Л2 (- 0,5 )п Решение задачи (25-26) получается в виде

я

к,

42

ехр

1 »-3 Л , 1 2

т 14

Осуществив преобразования, обратные (13), и учитывая (8), (9), для распределения температурных полей получим следующие формулы:

Ік ,у (Х т ) = 1ик ,у (ъ)Кк ,у

( х Л

у=1

V 41(т) у

+

+

1 -

41

іе при к = 1;

V 4і

-1

а -1

при к = 2,

где функции Кі

(хЛ

к ,у

4і(т).

определяются форму-

^ и

лами (23) и выражаются через комбинации вырожденных гипергеометрических функций [13].

Из условия (7) получаем соотношение для нахождения параметра Л:

X

7=1

^

ЭУ 12

ЭУ

-^ -\ґе = -оуЛ2. а -1 1 е 2

При т = 0 получим, что

Л2 = | ^ 2^0

уо V1 - а

С учетом того, что параметр Л не зависит от времени т, решение задачи (31-32) принимает вид:

^(т) =

я,

к,

2

Ик ,у

(

4 0

Л2т + 42

И’ к ,у Л

0 у

Это решение соответствует случаю, когда т^-0 и т ^х>. При т ^ к> для нахождения параметра Л из условия (7) получается следующее соотношение:

X

7=1

Я2,уЭК1,у

А<і 2 ^

Иі,

у

ЭУ

1

22 И 2,У

ЭУ

-\ґе =1 оуЛ2. а-1 1 е 2

В данном случае, при т ^х>, решение (33) стремится к автомодельному, которое будет иметь следующий вид:

Ік (х, т) = -Х

я,

к,

2 Кк ,1 7=1 Ик ,у

(Л

Лл/т

+

+

1 -

Лл/Г

Іе при к = 1;

Лл/т

а -1

при к = 2.

х

І

х

0

2

2

1

п

т

х

І

X

0

1

Представлением найденных значений Л и

ц 2ув (33) завершается решение поставленной

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

ОБОСНОВАНИЕ ПОЛУЧЕННОГО РЕШЕНИЯ

При решении задачи (19-20), а значит и (14-15), собственные функции определяются из однородных граничных условий вида (20), причем в первом из этих условий используется требование ограниченности при г^0 решения задачи. При этом свойства ортогональности собственных функций и вещественности собственных значений основываются на свойстве самосопряженности дифференциального оператора Ь уравнения (16) для класса функции, имеющих непрерывную вторую производную в промежутке (уь у2):

У 2

| {М - )лу = °

У1

где ЬКк ,у( у) = ^-Эу

Рк(У)-

ЭКкЛ. у)

к,

ЭУ

Самосопряженные дифференциальные формы Ь/ и Lg удовлетворяют тождеству Лагранжа [19]:

М - gLf =\рЖ (g, /)],

где & (g, Г) =

,(g,Г)е С2[у^у2].

Тогда

У 2

| {М - gLf )лу = р&

У1

V 2

VI

Если функция / и g удовлетворяют однородным граничным условиям (17), то оператор Ь будет самосопряженным за счет того, что

&

У2

У1

= 0.

Из теории задачи Штурма-Лиувилля [2023] следует, что ряды (33) сходятся в среднем равномерно по т, если интеграл

\\тк{У, т)]2 Рк(У) лу

равномерно ограничен относительно т, а функция Тк(у, т) дважды непрерывно дифференцируема в [уь у2] и удовлетворяют граничным условиям задачи Штурма-Лиувилля (16-17). Сходимость в среднем означает, что

У 2

Пт (

У1

Тк

у=1

X Рк (У) ЛУ = а

X

то есть ортонормированная с весом рк(у) система собственных функций Кк у(у) является полной [24] в рассматриваемом классе функций, которому принадлежит Тк(у, т). Необходимым условием полноты системы функций Кк У(у) является свойство замкнутости [24] этой системы, для которой выполняется равенство Парсеваля:

£2 тк2 Ыт )рк {у )лу = Ё ик, y{т),

1 у=1

доказываемое с помощью условия (36).

Если функция Тк(у, т) удовлетворяет сформулированным требованиям, то формула (33) представляет обратное преобразование. Из единственности разложения (33) следует взаимно-однозначное соответствие между функциями Тк(у, т) и ик(у, т). Поэтому, решив преобразованную задачу, как следствие, получаем решение исходной задачи.

ПОЛУЧЕНИЕ ЧИСЛОВЫХ РЕЗУЛЬТАТОВ

Для построения температурного поля в зоне промерзания согласно формуле (35) получается соотношение

^(х, т) = (1 - у)*е + ГеРхуе(1 у )в1;

. ^;Р. Ж,,3,|5,у2

ё

у=1

(1 -ь1,у)ь1,уЩ ь1,у +1,-2;р11°1,у

где у =

Лл/т

а величина о1 у является резуль-

татом дифференцирования Щ\ Ь1 у,~;Р1

по

первому параметру

¥ (Ь ) Р." П-1

Еч 1,У/" вт

/! г-4 | Ё

1

1

=0т + Ь

=1 (1,5)" "! т=0Ь1у + т

= У{" + Ь1,у)-

т=0

1,

^{Ь1,у) есть пси-функция.

Пользуясь рекуррентным соотношением для ВГГФ [17] можно записать

3

Ь1,у — 1Щ\ Ь1,У^;Р

5

Ь1,УЩ\Ь1,У +1,—;Р1 I+ -т -;р1 I = 0.

2

"

х

"

"-1

С учетом первого из уравнений (24) полу-

чим

Щ

Тогда формула (37) заметно упрощается:

2

^(х, т) = (1 - у)?е - ГеР1 уе(1-у )в1 X

3

¥ Щ

X Ё

у=:

2

1^1 ^ Ь1,у )01,у

Для построения температурного поля в зоне охлаждения согласно формуле (35) запишем соотношение

'о ,л-2 е _р2 у2

X Ё

И,

а -1

X

а.-Ь С2-‘ др( Ь2-1-р2-у21.

у=11 - Ьк,у

Выражение для Н2,у, данное в уравнении (31), может быть преобразовано к виду

и л2

н = с~1 ______

Н2 х С2,у6{1 - а)

5

аЩ Ь2,у, 2 ;Р2а I Щ1Ь2,у, 2 ;Р2

3а

4Р2 {1,5 - Р2,у)

3 -1

2, 2"

Ь2,у , ^ ;Р2а I

Что означают р и я указано в формулах (23), они включаются в нормирующие множители С2-,1 . Формула (40) содержит две разности.

За счет применения рекуррентных соотношений для ВГГФ [17]

{1 + а - V )Щ {а, V; г)- аЩ {а +1, V; г)+

+ {V - 1)Щ {а, V -1; г ) = 0;

vF {а, V; г )-уЩ {а -1, V; г)-гЩ {а, V +1; г )= 0

ко второй разности с учетом второго из уравнений (24) выражение в правой части формулы (40) сворачивается в полные производные

н = с ~1 ________{0Л у.

2,у 2,т 8{а-1){1,5 -Ь2,у) Р2

х

а^'у\Ь2:(,3; ра2)-АЩ (Ь2,у,|; Р2

Выражения для нормирующих множителей

С

к, также сворачиваются в полные производные

Сс,=-£Щ (ь4;р.)эгF^Ь1л•!;р I;

е~Р ( 3

С2 =-—АЩ к -;р2

4Р 2

2

1

Э

2

3

2 2,

С учетом полученных соотношений формула (39) преобразуется к виду

.2!

'2 {х,т) = {у -1)-^ - ?0р1ехр\{1 у2) р2 1

АЙ Ь

X Ё

г

а -1 3

;Р2 у

2{а -1)

X

X

1{1- Ь2,у ){1,5 - Ь2,у) О 2у

АЩ'у ку^^а2) г---—--------— -1

ащ; ку,3^

X

где у =

Лл/т

а величина о2 у является резуль-

татом дифференцирования АЩ\ Ь2у,—;Р2 | по

первому параметру

О 2,у = Ё

"=1

Х'"" {а - я)

{1.5),

"-1 X Ё

1

{Ь2,у - 0,5)"

=0 т + Ь2,

^ 2" 1Р0 - Р)x "-1 1

X Ё---------------------

т=0 т + Ь2,у

0,5

{0,5)

Р

-X

X

"

2 . "!;

Ё +1, = У{" + Ь2,у )-¥{Ь2,у ),

т=0т + Ь

Ё

2,

1

-0,5

т=0т + Ь2у

У{" + Ь2,у - 0,5)- у{Ь2у - 0,5)1

величины р, я, р0, я0 указаны в формулах

3 2

АЩ^\ Ь2,у^;р2у2 | = 4Р2

4Р

2

- ( 2 - 42-т ) Я У2 Щ ( Ь,,?Р = У21 +

+ а{1- Ь2,у )руЩ (Ь2,у - у,-|;р2у21

Начиная с у = 4 в формуле (43) используется первое из рекуррентных соотношений (41). На его основе для производных по первому параметру от ВГГФ получаются новые рекуррентные формулы.

2

х

X

X

я

оо

а

"

Э -Щ(Ь2,у,-2;Р2а2

2 ^ I ЭЬ2л

-2Щ(Ь2,У + 2,-2;р2а2)-(| + Р2а2 |X

XЩ(Ь2,у + 1,|;р2а2 1 (2 -Ь2,у IX

без учета порядка малости О

- - Ь1, 4 1

-0,5 Л

x |Ь2, Г Ё к" + 2 + *2,Т ) - ^2 + Ь2, Г Я-{|"

ВГГФ. Характеристическое уравнение (24) при этом принимает вид

{3 - 4Ь1,у)р1 = 0,

решения которого при Ь

1,у '

2 + 2Ьгг + р 2 а 21 "ГУ ьЦ |р 2Я 2):

{1,5)" "!

2,У/-' {1,5)" "!

Ё Ц" + 1 + Ь2,у )- ¥{1 + Ь2,у )]

^ Щ ^\’\^2а 2 ) = 2 Щ (Ь+М;М2|-- ( 2 + р2а 2 V ( ^ !,!; р 2а 2 ) -

Г (

- Ь

2, У

Ь 2„ +1 |Ё

2 I" = 1

У\ " + | + Ь2,у

Л

- ^\ - + Ь2,у

X

X

(1,5 + Ь2,т) {р2а2У1 -(1

(0,5)„

"!

-\ 2 + 2Ь2,у + р2а | X

X

^ " + 2 + Ь2,у 1 - | + Ь2,у

,5 +

X

X

{0,5+ь2„), ы:

С0,5)„ "!

Десятикратное применение этих формул влечет потерю примерно одной значащей цифры [17], но согласно расчетам это не повлияет на сходимость ряда (42) в силу его достаточно быстрой сходимости. Воспользовавшись асимптотическими представлением ВГГФ можно дать оценку суммы рядов в формулах (38), (42). Когда рУ < 10 и Ьку ^ - да имеет место следующая асимптотическая формула [17]

3-.2 ) Г 1-2

0,5

Щ\Ь,2;р; )=ехр\^р; Ix

X-

1п{;У {3 - 4Ь) р) У у/ {3 - 4Ь) р

(

1 + О

Тогда выражение под знаком суммы ряда в формуле (38) преобразуется так

Щ(ЬмЛ; р1 у 2

2

{1 - Ь1,у){1,5 - Ь1,у) О1 ,у

ехр

\0,5 р1 {у2 -1)]

2 У

X

•—^ ) { 1,У----) в1п0{3 - 4Ь1 у) р1у)

1 - *1,у ){1,5 - *1,т) ^ 1^Р1//

{щ)2

р1

-3

у е Ы/у.

Подставляя их в формулу (38), имеем

V3 - 4*1,, 1

{1 - ьЛи - *1,,) V

С другой стороны

{1 - Ь1,,К1,5 - *1,)

у / °1,у

ехр

\0,5р1 {у2 -1)] р - 4Ь1,,

2У ,

откуда следует, что при Ь1у ^ - да (у^да) решение (38) сходится не хуже числового ряда

Ё “т. То же самое утверждение имеет место и

7=1 У

для решения (42) в зоне охлаждения. Таким образом, получены окончательные формулы (38), (42), на основании которых нужно представить графики изменения температуры по пространственной и временной координатам. Сначала определимся с исходными данными. Значение безразмерного параметра теплового влияния для фазового перехода вода-лед можно взять равным а = 4 [12]. Для краевых условий зададим следующий температурный диапазон: 'е = -20° С, '0 = +10° С.

Теплофизические характеристики будем задавать в интернациональной системе: о =

= 335 103 ^; а1 = 1,255-10-6 —; а2 = 0,151 х кг с

х10-6 —; X = 2,44-^; X 2 = 0,63-^;

с м • К м • К

кг

V = 0,917 10 —. Прежде всего, найдем значе-

м

ние параметра Л по формуле (34).

Л2 = 0,30405 •Ю-6 —

и получим значения безразмерных параметров рк: р1 = 0,06057; р2 = 0,5034.

в

"=1

с

Т аблица 1

У b1,Y 01,у Ь2,у Ь2,у Gy

1 -39,99063660 0,012651708 -0,551623071 -235,7559131 -82,51310488

2 -162,2002809 -0,003162886 -2,311227245 53,60204448 73,24499588

3 -365,8821689 0,001405588 -5,038454622 -22,71512493 -193,5626213

4 -651,0367001 -0,000781960 -8,848920381 12,59507523 187,5437346

5 -1017,663924 0,000495600 -13,74908758 -8,015189847 -147,7108239

6 -1465,763853 -0,000815000 -19,73895690 5,476158764 175,4499255

7 -1995,336490 -0,011084899 -26,81833274 -8,501136619 -187,1514118

8 -2606,381838 -0,037830000 -34,98709674 440,7455479 159,9870429

9 -3298,899987 1,120000000 -44,24518293 -2108,65151 -176,2078465

10 -4072,890668 58,3034344 -54,59255349 5,025477  105 182,5251447

11 -4928,354152 7176,476967 -66,02918587 -2,659997  107 -164,3589950

Формулы (26) нахождения собственных чи-

сел с точностью до порядка О ются к виду

1'

8

преобразу-

V а1 J

= (nY)2 + Pi +1P2 +

(

+-

1

v

11

---+--2 +---4

45 12y2 16y4

P4

++

1

(

21 +-

9y4V105 2y2

P6+

2835/

P8;

Л

+-

V a2 J

24

=9 (nY)2 + P2 + 7P2 +

У

1

+

1

3

+—

у

5 + 12y2 16y4

Pi +

99 4Л

70 + T2y

R6 + 2151 R8-P2 + 80y6 в2;

При данных значениях параметра pk собственные числа задачи (21-22) будем находить из уравнений (24) с использованием программной системы Maple. В зоне промерзания для вычисления собственных чисел можно пользоваться и уравнением (24), и формулой (44), причем результаты в пределах заданной точности будут совпадать даже в том случае, если в формуле (44) мы ограничимся первыми тремя слагаемыми. Значения величин ok,Y в формулах (38), (42) также вычисляем в программной системе Maple. Кроме того, воспользуемся обозначением

AF; | Р2а2

I b2,r2;P2

-1

и заполним табл. 1. Основываясь на этой вспомогательной таблице, составим табл. 2 изменения температуры вдоль безразмерной пространственной координаты в момент времени т0 = 10 ч (при этом значении у = 1 будет соответствовать расстоянию х = 0,10462м) и табл. 3 изменения температуры вдоль временной координаты на расстоянии х0 = 3,3084" 10 2 м от начала координат. В эту же таблицу внесем зависимость положения свободной границы от времени, а также соответствующие значения безразмерной пространственной координаты. Все значения температуры даются в градусах Цельсия (оС), расстояние - в метрах (м), а время - в часах (ч). При этом вычисления производятся в СИ. Заметим, что ряды в формулах (38), (42) достаточно быстро сходятся. Скажем, для получения значения температуры в оС с точностью до 10 2 в зоне промерзания достаточно четырех, а в зоне охлаждения - восьми первых слагаемых. По составленным таблицам строим графики распределения температурного поля вдоль пространственной и временной координат, а также график зависимости положения свободной границы от времени (рис. 1-3).

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

1

2

У

2

2

1

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

ВЫВОДЫ И ЗАКЛЮЧЕНИЕ

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

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

Т  ¥ .

2. Решению задачи соответствует параболический закон движения свободной границы:

ад=ТЛ т + £,2 , где постоянный параметр Л находится из условия Стефана. При данных краевых условиях и теплофизических характеристиках скорость движения свободной границы определяется соотношением 2,1(т) =

= 0,5514 • 10-3 \/т в Международной системе СИ единиц физических величин.

3. Собственные значения оператора Штур-ма-Лиувилля зависят от параметра Л. В зависимости от величины в к= (4ак )-1 Л2 собственные числа находятся тремя различными методами, два из которых асимптотические, а два из них сочетаются с численными методами.

4. На основе асимптотического представления по первому параметру ВГГФ удается установить, что ряды в полученных решениях на бесконечности сходятся не хуже числового ряда ^ п~3 . Численные эксперименты показыва-

П=1

ют, что в зоне промерзания для получения значения температуры в градусах Цельсия с точностью до 10-2 достаточно четырех слагаемых, а в зоне охлаждения - восьми слагаемых.

Т аблица 2

У *1( У, т0) У 0 ,т 2 У *2 (У, т0 ) У 0 ,т 2

0,0 -20,00000000 1,0 0,000000000 2,0 8,718479689 3,0 9,838480529

0,1 -17,59810588 1,1 1,703297768 2,1 8,740267435 3,1 9,917735615

0,2 -15,61189533 1,2 3,210989542 2,2 8,763608958 3,2 10,01492898

0,3 -13,63326830 1,3 4,507429355 2,3 8,846878326 3,3 10,11085213

0,4 -11,66085318 1,4 5,609424009 2,4 9,005563962 3,4 10,17569521

0,5 -9,697227174 1,5 6,543905060 2,5 9,212867772 3,5 10,18639913

0,6 -7,743456676 1,6 7,323484020 2,6 9,420658622 3,6 10,14038258

0,7 -5,795889319 1,7 7,938320581 2,7 9,587881460 3,7 10,05929434

0,8 -3,854869694 1,8 8,370183465 2,8 9,700597157 3,8 9,982033493

0,9 -1,925009667 1,9 8,618348478 2,9 9,774103999 3,9 9,951289471

1,0 0,000000000 2,0 8,718479689 3,0 9,838480529 4,0 10,00000000

Т аблица 3

т О т ^1(т)= Лл/г

1 0,000000000 0,033084000 1,00000000

2 -5,657672112 0,046787842 0,70710678

3 -8,185291291 0,057303169 0,57735031

4 -9,697227174 0,066168000 0,50000000

5 -10,73245396 0,073978073 0,44721362

6 -11,49850489 0,081038919 0,40824830

7 -12,09486237 0,087532036 0,37354369

8 -12,57610973 0,093575683 0,35355340

9 -12,97505034 0,099252000 0,33333333

10 -13,31274403 0,104620794 0,31622778

11 -13,60341043 0,109727215 0,30151135

12 -13,85704181 0,114606338 0,28867514

13 -14,08089164 0,119286058 0,27735010

14 -14,28036923 0,123788993 0,26726124

15 -14,45960152 0,128133781 0,25819889

16 -14,62179901 0,132336000 0,25000000

17 -14,76950180 0,136408827 0,24253563

18 -14,90474940 0,140363525 0,23570226

19 -15,02920080 0,144209813 0,22941573

20 -15,14422118 0,147956146 0,22360680

21 -15,25094551 0,151609934 0,21821789

22 -15,35032627 0,155177715 0,21320072

23 -15,44316948 0,158665290 0,20851441

24 -15,53016245 0,162077837 0,20412415

25 -15,61189533 0,165420000 0,20000000

5. Из графика распределения температурного поля по пространственной координате видно, то температурная кривая в зоне охлаждения выходит за пределы граничных условий, т. е. Ї (х, т )> /0 . Это наглядно демонстрирует тот факт, что процесс фазового перехода сопровождается выделением теплоты при кристаллизации, что невозможно продемонстрировать в случае классических постановок задачи Стефана.

Как видно из (38) и (42), нами получены компактные вычислительные формулы, и они содержат ряды, которые хорошо сходятся, что является очень удобным при построении температурных полей.

Рис. 1. Кривая распределения температуры г(у, т0) (в (°С)) по безразмерной пространственной координате у в момент времени т = 10 ч

г

Рис. 2. Кривая распределения температуры г(х0, т) (в( С)) по временной координате т (в (ч.)) на расстоянии х0 = 0,33084 м от начала координат

£1

т

Рис. 3. Кривая зависимости положения свободной границы £,1

СПИСОК ЛИТЕРАТУРЫ

1. Рубцов Н. А., Слепцов С. Д., Савинова Н. А. Численное моделирование однофазной задачи Стефана в слое с прозрачными и полупрозрач ными границами // Прикладная механика и техническая физика. 2006. Т. 47, № 3. С. 84-91.

2. Dancer E. N., Du Yihong. A uniqueness theorem for a free boundary problem // Proc. Amer. Math. Soc. 2006. V. 134, № 11. P. 3223-3230.

3. Borisovich A., Friedman A. Symmetry-breaking bifurcations for free boundary problems // Indiana Univ. Math. 2005. V. 54, № 3. P. 927-947.

4. Baconneau O., Lunardi A. Smooth solutions to a class of free boundary parabolic problems // Trans. Amer. Math. Soc. 2004. V. 356, № 3. P. 987-1007.

5. Карташов Э. М., Рубин А. Г., Ожерелко-ва Л. М. Проблема Стефана // Математические модели физических процессов: Сб. науч. тр. 10 Меж-дунар. конф. Таганрог, 2004. С. 88-92.

6. Калиев И. А., Вагапова Э. В. Задача Стефана как коэффициентная обратная задача // Современные проблемы физики и математики: Тр. Всероссийск. науч. конф. 2004. Т. 1. Уфа, 2004. С. 43-49.

7. Ушакова В. И., Клочков А. В. Поведение решений однофазной задачи Стефана при больших значениях времени // Известия вузов. Математика. 2006. № 11. С. 55-60.

8. Бижанова Г. И., Солонников В. А. О задачах со свободными границами для параболических уравнений второго порядка // Алгебра и анализ. 2000. Т. 12, № 6. С. 98-139.

9. Карташов Э. М. Аналитические методы решения краевых задач нестационарной теплопроводности в областях с движущимися границами // Изв. РАН (Сер. Энергетика). 1999. № 5. С. 3-34.

10. Шафеев М. Н. Решение одной плоской задачи Стефана методом ВГГП // Инженерно-физический журнал. 1978. Т. 34, № 4. С. 713-722.

11. Шафеев М. Н. Решение одной нелинейной задачи методом ВГГП // Известия вузов. Математика. 1980. № 12(233). С. 73-75.

12. Хакимов Х. Р. Замораживание грунтов в строительных целях. М.: Госстройиздат, 1962. 257 с.

13. Дикий Л. А. Формулы следов для дифференциальных операторов Штурма-Лиувилля // Успехи математических наук. 1958. Т. 13, № 3(81). С. 111-143.

14. Федорнюк М. В. Ассимптоматические методы для линейных обыкновенных дифференциальных уравнений. М.: Наука, 1983. 352 с.

15. Прудников А. П., Брычков Ю. А., Мари-чев О. И. Интегралы и ряды. Дополнительные главы. М.: Наука, 1986. 800 с.

16. Зайнуллин Р. Г. Об одном аналитическом подходе к решению одномерной задачи переноса тепла со свободными границами // Известия вузов. Математика. 2008. № 2. С. 24-31.

17. Абрамович М., Стиган И. Справочник по специальным функциям с формулами, графиками и таблицами. М.: Наука, 1979. 832 с.

18. Слейтер Л. Д. Вырожденные гипергеомет-рические функции. М.: ВЦ АН СССР, 1968. 178 с.

19. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1972. 736 с.

20. Коддингтон Е. А., Левинсон Н. Теория обыкновенных дифференциальных уравнений. М.: ИЛ, 1958. 474 с.

21. Левитан Б. М. Разложение по собственным функциям дифференциальных уравнений второго порядка. М.-Л.: Техтеоретиздат, 1950. 159 с.

22. Березанский Ю. М. Разложение по собственным функциям самосопряженных операторов. Киев: Наукова думка, 1965. 798 с.

23. Котляков Н. С., Глинер Э. Б., Смирнов М. М. Уравнения в частных производных математической физики. М.: Высшая школа, 1970. 712 с.

24. Никифоров А. Ф., Уваров В. Б. Специальные функции математической физики. М.: Наука,

ОБ АВТОРЕ

Зайнуллин Рифат Гильметди-нович, ст. преп. каф. математики. Дипл. физик (КГУ,1982). Иссл. в обл. разработки анали-тич. методов решения спектральных задач с фаз. переходом для уравнений парабол, типа, иссл. и матем. моделирование процессов диффузии углеводорода в металл.

1984. 344 с.