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

УДК 532.546

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

А. Р. Павлов, И, Г, Ларионова, М. В, Михайлова

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

1. Задача о промерзании влажной пористой среды сводится к системе уравнений тепломассопереноса с фазовыми переходами [1,2].

Температурная задача:

дТ4 дх ,

ди>9

На границе раздела фаз имеет место условие типа Стефана:

Л21Х - Х1дХ = ь[т* - ^Т*)]§, Т = Т. (х = т- (з)

Граничные условия для случая промерзания с левой границы: дТ

Л— = «(Т - Тс), х = 0, * > 0, (4)

дТ

Л— = о, х = 1,г>о. (5)

дх

дТ д (

С

дх

дТ д ( дТ

дх Л дх

(1) (2)

© 2006 Павлов А. Р., Ларионова И. Г., Михайлова М. В.

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

Т(х,0) = ф), х < I. (6)

Введены следующие обозначения: С^Х^ — объемная теплоемкость и коэффициент теплопроводности г-й фазы, Т* — температура фазового перехода, —* — значение влажности на фазовой границе со стороны талой зоны, —^Т*) — такая же величина со стороны мерзлой зоны, <ш\(Т) — известная функция, выражающая количество незамерзшей воды при данной температуре, —2 — влажность по льду, Ь — скрытая теплота, выделяющаяся при промерзании куб. м. воды, а — коэффициент теплообмена с окружающей средой, Тс — температура окружающей среды, х = £(£) — положение фронта фазового перехода.

Влажностная задача. Рассматривается движение влаги как в талой, так и мерзлой зонах:

Ж = Ш (ьдю )• Т>Т,. (Г,

Ж = ) Т<Т*. (В,

где — = —I + На границе фазового перехода скачком замерзает свободная влага —* — —^Т*) = —г(Т*). Уравнение баланса массы на фазовой границе можно записать в виде

дтг дт

к— к1дХ = МТ*] <. (9)

л; л; ы/о

Граничные условия для влажности записываются в зависимости от фазового состояния границы: д—

= а*(т — ф), х = О, Т > Т*, (10)

дх д—

к2—— = а*(—1 — ф), х = 0, Т < Т*, (11)

дх

о

= 0, х = 1, Т>Т*, (12)

дх

-^ = 0, x = l,T < Т*. (13) dx

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

w(x, 0) = , 0 < x < l, (14)

w — влажность в талой зоне, ф — равновесная с окружающей средой влажность материала, а* — коэффициент влагообмена.

Приведенная математическая модель, в частности, описывает процессы тепломассопереноса при промерзании и оттаивании грунтов, интерес исследователей к которым не ослабевает [3,4].

Уравнения (1)—(3) обычно принято записывать в виде одного уравнения, заданного во всей области. Аналогично поступаем и с уравнениями для влажности. Заметим, что функция T(x,t) в первой задаче является непрерывной функцией. А решение второй задачи, если рассматривать влажность по жидкой фазе, терпит разрыв первого рода на фазовом фронте. Для того чтобы написать влажностную задачу относительно непрерывной функции, введем новую функцию:

w, T > Т*

w, т <т* .

Функция Ф(Т) при T ^ Т* — 0 принимает значение Ф(Т) = w = w± (Т*) + w2(T*) = w*, а при Т ^ Т* + 0 — значение, равное w*, т. е. она непрерывна в точке Т = Т*.

Имеет место следующее утверждение: уравнения (7)—(9) можно заменить одним уравнением

дФ д / дФ N д / 3w2 \ д.. .

-d^ = dx[kdx) — ex [kldx) + rn{Sw)' (15)

определенным во всей области 0 < x < l, где

k= i h, т>т*, ( 1, т = т* (x = m),

\k2, Т < Т*, \о, ТфТ* Справедливость утверждения при Т < Т* и Т > Т* очевидна. Для доказательства того, что уравнение (15) включает и условие (9), воспользуемся методикой работ [5,6]. Выберем область па плоскости XOT,

*

Ф =

Рис. 1.

представлющую достаточно малую окрестность границы фазового перехода с контуром ABC DA, внутри которой находится фазовый фронт EF (рис. 1).

Пусть слева от линии EF (х = £(t)) находится мерзлая область, справа — талая. Проинтегрируем уравнение (15) по данной области и воспользуемся формулой Грина

<9Ф 3w2 \ Í к— - k2-jX )dt+ (Ф - w2{T*)) dx = О,

г г

ABCDA BC DA

FE

, dw\ , f / 8w\ 1 f , 3w2 , ki— dt + [k2—\dt - k2 —— dt dx J J \ dx J J dx

EF FE FE

j w* dx + j w dx — j w2(T*) dx = 0.

EF FE FE

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

/( дт ди!1 \ [

--~д—) — т(Т*) ах = О,

ЕЕ ЕЕ

откуда следует (9). Совершенно аналогично доказывается, что уравнения (1)^(3) заменяются одним уравнением

дТ д / дТ ) тдт2

ст=д-х{ХдХ)+1-дГ> (16)

определенным во всей области, где

с= { с, Т>Т., г аь Т>Т, \с, Т<Т*, { \2, Т<Т*.

Граничные и начальное условия (10)—(14) формулируются относительно функции Ф(х, £):

= - ф), х = 0,Т>Т*, (17)

к2^(Ф - Щ2) = а*(Щ - ф), х = 0, Т < Т*, (18)

дх

д

= х = I, Т > Т*, (19)

А(ф - щ) = 0, х = 1,Т < Т*, (20)

дх

Ф(х,0) = <*, 0< х < I. (21)

2. Разностная задача. Вводится равномерная по каждой переменной разностная сетка с шагами Н и т: х^ = гН, г = 0,1,..., х^ = N1 = I; ^ = зт, 3 = ОД,.. ,,зо.

На границе фазового перехода коэффициенты уравнений (15), (16) имеют разрывы первого рода и не определены на ней. Для построения разностных схем сквозного счета разрывные коэффициенты заменяются сглаженными (непрерывными) функциями. Сглаживание проводится в окрестности точки фазового перехода (Т* - А,Т* + А), например, линейной интерполяцией.

Прежде чем строить разностную схему для температурной задачи, исключим производную щ по £ с помощью уравнения (8):

Для этого уравнения интегро-интерполяционным методом строится итерационная разностная схема

т т+1 т т+1 т т

С г Тц = (А—.5 Тх)х,г + Щк2 )г_0 .5 (т )х)х,г, г < N - 1, (23)

где т — номер итерации.

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

Интегрируя уравнение (22) по х от х = 0 до х = 0.5Н, получим

дш\

Н—дТ

АдТ дТ _ , дт

- А— дх

0 дх . дх 0

- Ько

0.5

дх

Используя условия (4), (11), получаем разностное граничное условие

следующего вида:

т т т т

тА0,5 Тх,о + тЬ(к2)0.5(т)х,о

т+1 V

т

= ^С0( Т - Т)о + та( То - Тс)+ тЬа*("т - (24) Аналогично строится правое граничное условие:

т т+1 т т Н" т+1 V

тАN_0.5 Тх,м + тЬ(к2)м-о.5(т= С^ Т - Т)N■ (25)

Начальное условие выглядит так:

Т(хг,0) = ^(хг), г < N. (26)

Аналогичным образом строится итерационная разностная схема для влажностной задачи:

т тт т т

Ф- = (к(Ф)х)х - {к2{т2)

(27)

т т+1 Н т+1 т+1

(У о^)х^(Ф)-,^ «*(($) о-'М, Т0>П, (28)

т т Н т+1 т Н V

(к2)0.5(т)х,о=2Ф-,о + а*(т - ')о - 2 -,о, Т0 < Т*, (29)

m m+1 h m+1

- Al) n-o.6(*)x,n = ¿^-.n, Tn >T*, (30)

m m h m+1 v

- (k2) N -o .5 (w) x,N = -( $í,N - ^W) -,n), Tn < T* , (31) Ф(х< ,0) = ^*( x¿), ¿ = 0,1,...,N. (32)

T

mm

ф — решения разностных задач на текущем моменте времени, T , Ф — решения задач (23)-(26) и (27)-(32) соответственно. Введем две функции

m m m m

= T - T, my2 = Ф -Ф. Тогда для них получаем следующие задачи:

^-m+l m m+1 m m m

C ^^^^ ^)x)x + t((A - A)Tx)x + rL(k2{w1 - w)x)x

m ^

+ TL((A2 - ^(w)x)x - r(C - C)T-t, (33)

m m m

тaq.5 yxx,o + r(A - A)o.5Tx,0

m m m h ■mm+1

+rL(k2) 0.5 (w - W )x,0 + rL(k2 - k2) 0.5 (W )x,0 = -jCü У1,0 +Tomyl,0 + rLa*(- w)o + Th{C0 - CQ)T-fi, x = 0;

mm

tan-o.5 yix,N Vм/

m m m

+t(A - A) N-0.5 Tx,N + TL{k2) n -0 .5 (W - W )x,N

т m / ч h m m+1

+TL(k2 - k2 )n -0.5 (W )x,N = C N У1 ,N

h m _

-t — {cN - Cn)T-,n, X = l,

m+1 ,m,m+\ . m ,чж

y2=T(k( y2 )x)x + T(k - k)9 x

-T{k2(W2 - W2) x)x - T{{k2 - k2){W2 )x]

m m m

T k . y x, T k - k . x,

= НСш)о + ™*("Йо, х = О, Т0 >Т*,

т т т

т{к2)0.5(т - т)х,о + т(к2 - к2)0.5(т)х,о

= НГ$о + та*(" - ^о, х = 0, Т < Т*,

т т т

-т{к\)N _о .б( У2 )x,N - т(к\ - к\ )N _о

= Н("Й)N, х = I, TN >Т*,

т т т

-т{к2 ) N _0 .5 (ы-ш - т )x,N - т(к2 - к2) N _0 .5 (т )x,N

= НСш)N, х = I, TN < Т*. Далее используются скалярные произведения сеточных функций [7]

N _1 N

г

N _1

[Уг, 'Ш] = УгПг+ N ^ '

:УоПо.

г

тУ

формулу Грина [7]

,""+1 т+1 ,","+1 ,т+1ч ,

(С У\, У1) + т(А{ У1 )хл У1 )х]

т т т т

— т У1 ^ [А N _0 .51 У1 + (А - А N _0 .5 Тх^

" т " т+1 " т+1 + N_о.5(т -т+ -к2)N_о.5(т)х^] -т У\,о[А0А У\)х,о

т т т т

+ (А - А)0.5Тх,о + Цк2)0.5(т - т) + - к2)о.5(т)х,о]

т т+1 т т т+1

- ^(А - А)ТхД ^- тЬ(к2(т! - т)хЛ У1)х]

- тЬ{(к2 - к2)Юх, СУЬх] - ^(С - С)ТЪ тУЬ.

Отсюда, пользуясь граничными условиями для уравнения (33), придем к соотношению

т 7 т

,7^+1 т+1 /т+1 ч2

(С У\, У\) + т:СМ У1N

т:Со( У\ ,о) +ПМ У1 )хЛ У1)*

Ь т+1 /т^ тч \гг т+!

= 2Т У1NСN - Ск)- т У1 >°

а у1 -о + Ьа*( Ш - ^1)0

Ь — — '

— {Со - С0)Тё,о

т т+1 т т т+1

- г((Л - Х)Тх,( Ух - тЬ(к2(шх - ш )х, ( У\)х]

- тщ12 - х, СУЬ*] - 4(С - ОД, "ш). (35)

Коэффициенты уравнения (22) являются функциями С = С(Т, шх,ш2), к\ = кх{ш), к2 = к2(шх,ш2), Л = Л{Т), и пусть они имеют ограниченные производные по указанным переменным

С С дС,т ,ЗС т дС т

С - С = -—(ш - Ш) + я—(Ш - ш2) + —(I -дш дт2 дТ

но

т ¿Шл т т т ¿Шл т

Ш2= Ф- Шх, Шх - Шх = ——Ух, Ш - Ш2 = У2--— Ух-

аТ аТ

Окончательно получаем

т ( дС дС \ 3,Шх т дС т дС т

С - С= (дШх - дШ)1тУ1 + дШ:У2+дтУ1-

Аналогично

к к — ( дк2 дк2 \ ¿шх т дк2 т

2 2 дшх дш2) ¿Т У дш2 У2'

дкх т т дЛ т

^ - к1=дШУ3' Л - Л=дТУ1-

Учитывая полученные выражения, оцениваем сверху правую часть (35) с помощью е-неравенства

1 ,1 9 Ь2 \а* Ь\ ^ еа Л--.

е

Для сокращения записи все постоянные обозначим одной буквой М, не указывая их структуры, как это обычно делается при выводе априорных оценок решения разностных задач [7].

Окончательно придем к соотношению

/ ■ "77 л гт+1 т+1п / • \ о \ г/т+1\ /т+1 п

(тш С - тМ)[ У1, ^ + -Зе)[( ^)х Д У\)х\

-х г < г^ ^ л //'"г \ ч л гт т лч

< тМ^1, У1] + ((У1 )хД^)х] + [У2, (36)

Аналогичные выкладки производим для задачи (34). Умножив ска-

У

имеем соотношение

,т+1 т+1 ,""} ,т+1 -,

( У2, Ы +ПМ У2) хЛ У2) х\

т+1 т т+1 т

= т У2,N[кN_0.б(( У2^)х) + (к - к)N_0.бФх^

т т т

- (к2)N_о.5(т - т)x,N - (к2 - к2)N_о.5(т)х^]

т+1 т т+1 т т т

- т У2,0[ко.б( У2,о)х + (к - к)о.бФх.о - {к2)0.5(т2 - т)х,о

т т т т+1

- {к2 - к2)0.5(т)х,о] + т(к(т2 - т)х, ( У2)х]

+ т((к2 - ^(т)х, Гй)х] - т((к - к)Фх, СЙ)х]■ (37)

При оценке сверху правой части (37) учитывается фазовое состояние границ и используются соответствующие граничные условия из (34). Таким путем получаем в общем случае следующее соотношение:

гт+1 т+1-, . . . пт+1\ 1т+\ п

[ У2, У2\ + (тш к2 -Зе)т(( У2)х Д У2)х\

-х г < Г™ т Л / \ \ Л / \ \ Л гт т -,ч

< тМ[й, У1] + «У1) х, Ы )х] + ( Ы) х Д^) х] + [У2, У2]). (38)

В правых частях (36) и (38) воспользуемся неравенством

4

(Ух^х] < ^"[у^]

и объединим их. При достаточно малых е и т, удовлетворяющих условию

шш^тМ > О,

получим

гт+1 т+1 гт+1 т+1 /гт т гт т 1ч

[ Уь Ы < Мт([Уь У1] +

Отсюда при достаточно малом т < М следует сходимость итерацион-

ной схемы, т. е. доказана

Т

Если функции С(Т,шиш2), Л(Т),к1(ш), к2(ш\,и)2) имеют ограниченные производные, то при достаточно малом т решения итерационных

тт

Т,

ЛИТЕРАТУРА

1. Лыков А. В. Явления переноса в капиллярно-пористых телах. М.: Изд-во техн.-теоретич. лит., 1954.

2. Иванов И. С. Тепло- и массоперенос в мерзлых горных породах. М.: Наука, 1969.

3. Гамаюнов И. И., Гамаюнов С. Н. Перенос тепла и влаги при промерзании грунтов // ИФЖ. 2004. Т. 77, № 5. С. 72-81.

4. Бровка Г. П., Иванов С. И. Расчет температурных полей в грунте с фазовыми переходами вода-лед в спектре температур // ИФЖ. 2004. Т. 77, № 6. С. 112119.

5. Годунов С. К., Рябенький В. С. Разностные схемы. М.: Наука, 1973.

6. Коновалов А. И. Задачи фильтрации многофазной несжимаемой жидкости. Новосибирск: Наука, 1988.

7. Самарский А. А. Теория разностных схем. М.: Наука, 1983.

г. Якутск

28 апреля 2006 г.