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

УДК 532.546 + 536.425

А.А. Папин

Математические вопросы динамики снежного покрова*

A.A. Papin

The Mathematical Theory of Snow Cover Dynamics

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

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

Melting snow is considered as a three-phase medium consisting of water, air and ice. A mathematical model taking into account the phase transitions and capillary jump, is based on the equations of conservation of mass, energy, Darcy’s law and the Laplace formula. The problems of model validation are considered.

Key words: two-phase filtration, capillary jump, heat and mass transfer.

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

[3, 7, 8].

Следуя [5-8], будем рассматривать тающий снег как сплошную среду, состоящую из подвижной (ш) и связанной (Ь) воды, льда (*) и воздуха (а). Подвижная вода фильтруется в снеге, а связанная - неподвижна относительно порового скелета, движением водяного пара пренебрегают.

Каждая компонента V (V = т, Ь, г, а) смеси занимает объемную долю фи (в единице объема смеси). По определению имеем

0 < фи < 1, фф + ф™ + фь + фа = 1.

Пористость (доля пор в единице объема) может быть выражена из равенства

фр = ф™ + фь + фа = 1 - фф.

Также вводятся общая насыщенность воды = в™ + вь, остаточная насыщенность воды 3™% = вь и эффективная насыщенность воды

Se =

Sw - Sw і — S .

wi

ви

І- въ■

По определению приведенные (рф, р™, рь, ра) и истинные (р1, р^, рА) плотности для каждой компоненты связаны равенствами:

рф = ффр1, р™ = ф™ рш, рь = фьрш, ра = фарА.

Уравнение баланса массы для компоненты V:

дри д

г \ ( V У\ V

-дп + Ш(р ” > = т '

где т = ^ т->, т- интенсивность перехо-

да массы из ]-й в v-ю составляющую в единице объема в единицу времени; у1' - скорость фазы.

В медленном потоке, который преобладает в тающем снеге, ускорениями пренебрегают. Уравнения сохранения импульса для каждой компоненты V представляют в виде:

дР . V . *2 V 1 \ иШ ( V (Л \

--д^ + рд + рв =^ 2^ т (у - у ),

Объемные доли подвижной воды, связанной воды и воздуха (на единицу объема пор) вводятся следующим образом

в™ = ф™/фр, вь = фь/фр, ва = фа/фр.

*Работа выполнена при финансовой поддержке государственного задания Министерства образования и науки Российской Федерации №1.3820.2011 и гранта РФФИ (проект 13-08-01097).

где д - ускорение силы тяжести; р = ^^ р - плотность смеси; ри - давление; Ви - тормозящая сила, возникающая из-за межфазного взаимодействия.

Предполагается, что температура Т является общей для всех фаз и уравнение сохранения энергии имеет вид

Е

VCЯ DvT

p Cp Dt

p r

Здесь —. = dt + vv dZ, Cpf - теплоемкость при постоянном давлении смеси; rv - внешний поток тепла; KN - теплопроводность смеси; Liw - удельная теплота плавления льда.

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

pBw = pWd^w — (фw)2(vw — vi),

dz kw J

r?a A дф p j ^a\2/ a i\

pB = p ~3^ — ka(ф ) (v — v ), где pw, pa - вязкости воды и воздуха; kw ,ka -проницаемости, связанные с эффективными проницаемостями krw,kra равенствами krw = kw/к, kra = ka/k; k - внутренняя проницаемость (k = a exp ЬфР; a,b,p - параметры). Предполагается, что krw = krw (Se) (в частности, krw = (Se)n).

Согласно формуле Лапласа разность давлений воды и воздуха связана с капиллярным давлением pC:

pC = pA — pW > 0.

Капиллярное давление также зависит от Se, т.е. pC = pC(Se).

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

% + dzwV)} =т<; pW {ж + д ф vw >} = т'ь—m

-t дz

’ — (

-і + дz1

pW\ д(фъ^)\ = -mwb - (1 - n)mi;

дд

—ї (rpA) + —~z (фpaAva) = °°;

^ + pWg - фw (vw - vi) = -1-(nmi--z 1 kw v ’ 2фш v

-pA i,A

-mwb)(vi - vw) - —L- + pAg - 7-фa(va - vi) = О.

-z

ka

После некоторых преобразований полученной системы и привлечения соотношения

-7j -1

Т = Тм + 8 х 10-7рс - 5 х \о-‘в:

приходим к следующему уравнению для эффективной насыщенности

д

pW g -

-t (^Se) + 7W/k-z

(krw) +

1 д (,rw -p

7W/k -z \ -z

—Q

+

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

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

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

+ Му(рфщ) = ^21зФ, г = 1, 2, 3, з=1

3

1]г = -1ф], Ь] = 0; (1)

г,3 = 1

k

Щ = -Ко —-(Уpi + p0g), i = 1,2,

2

p2- pi = pc(si,в), ^2s- = 1; (2)

i=i

дв

(E p0ic-ai) —і + (E p0civi)Ув

i=i

i=i

+ — (фeSev-) = mw>} - nm-.

= Злу(\сУв)+ Vдрa.. (3)

Здесь щ - скорость г-й фазы; рф - приведенная плотность, связанная с истинной плотностью р0 и объемной концентрацией аф соотношением рф = аФр0 (условие ^2®=1 а = 1 является следствием определения рф); 1]Ф - интенсивность перехода массы из ]-й в г-ю составляющую в единице объема в единицу времени; уФ = тзфщ -скорости фильтрации воды и воздуха; т - пористость снега; в1, в2 - насыщенности воды и воздуха («1 = тв1, «2 = тв2, аз = 1 - т); Ко -тензор фильтрации; к0ф - фазовые проницаемости (коф = коф(вф) > 0, коф|я,=0= 0); - динамиче-

ская вязкость; рФ - давление фаз; рс - капиллярное давление, д - вектор ускорения силы тяжести; в - температура среды (вФ = в, г = 1,2, 3,

ci = const > 0 - теплоемкость i-й фазы при постоянном объеме; v = const > 0 - удельная теплота плавления льда; Ас - теплопроводность снега (Ха = a,с + ЬаРс, Pc = Е3= 1 P0ai, ac = const > 0, bc = const > 0).

Система дополняется гипотезами [5]

U3 = 0, I13 = 113 (6), 112 = 0,

I23 = 0, p0 = const. (4)

После этих гипотез из уравнения неразрывности для льда следует дрз ^ m) = I31(9). В частности, можно считать, что пористость - функция температуры.

Классические задачи фильтрации о движении двух несмешивающихся несжимаемых жидкостей в пористой среде основаны на модели Маскета-Леверетта [9]. В большинстве задач пористость считается постоянной либо заданной функцией точки [10]. В [10, 11] построена теория для системы (1)-(3) в случае m = m(x). Важным моментом этой теории является доказательство классического принципа максимума для насыщенности

0 < s < 1.

В общем случае задача (1)-(4) является очень сложной и изучена лишь в автомодельной постановке [12-14]. Целью работы является анализ разрешимости системы (1)-(4) при заданной температуре 6(x,t), т.е. в случае зависимости пористости от x, t.

Задача двухфазной фильтрации в тающем снеге

Преобразование уравнений. Считая температуру заданной, систему (1)-(2) с учетом дополняющих гипотез (4) можно привести к виду:

д

— (ms1p1) + div(p01v1) = I31, (5)

д

— (ms2P0) + div(p2'V2)=0, (6)

-(p3(1 - m)) = I

13

I

31,

k

= -Ко^(Vp- + pig), i =1,2,

P2 - Pi = Pc(si,0).

(Т)

(В)

Сложив уравнения (5) и (7), получим преобразованное уравнение неразрывности первой фазы

д

(та1Р°1 + Рз(! — т)) + Лтрщ) = 0.

Складывая уравнения неразрывности (5) и (6), предварительно поделенные на р° и р2 соответственно, и учитывая равенство в1 + S2 = 1, приходим к соотношению

— (т + р3 (1 — т)\ + йтV = 0, (9)

дь V Р° )

в котором V = VI +«2 - вектор суммарной скорости фильтрации.

Введем новую искомую функцию [2]

P = Pi -

dpc —о2

1

ds и2 + -02

П2 Ml M2

de, (10)

где s = S1, к = ко1 + к°2, к°® = . Чтобы объ-

яснить такой выбор искомой функции, выразим предварительно вектор V из уравнений (8) через Ур1 и Vs

—V = К1 ^р1 + р° д) + К2 (Vp2 + Р° д), где К® = К°к°®, и с учетом р2 = р1 + рс получим

- v = Ki^pi + p.g) + К2(У(п + Pc) + p2g)

= k- Iу I pi - f % k2dei +

ds k

+ /Уis -k2de+itypc\ +1 >^-p^g

(e K-p-їь) .

(11)

Заметим, что соотношение (11) получено с использованием следующих равенств

др

Vpc(х,.?) = Vpc + Vpc,

дs

1 1

Г dpc kо2 dpc kо2 ж dpc kо2 „

У ~ds~d = У-es-de - ss-ys

где в Vpc символ V применяется только по переменной х, входящей явно (например, при п = 3,

VPc(X,s) = (y^Pc, дХ-Pc, дх;р^)). Таким обра-

c

зом

'= -(KVp + f) = v(s,p),

(12)

2

где f = KJ У ds. Щ2 de + K2 Vpc + К^о g),

К = кКо, с помощью подстановки (10) в уравнение (11) вектор V представляется через Чр и s и не зависит от Vs.

Аналогично с учетом (10) выразим вектор VI из уравнений (8) через Чр и Vs

- v. = Ki(ypi + pig) 1

= к. |У |p + 1 ^^-ї^1 + pill =

V ^ , I V7dPc -о2 JC dPc -о2^ , о^ l

Kl IVp + У - ds—+ pill ’

1

1

v

откуда, полагая а = — дт к01£02 и

/о = К1 / V дс к0- + К1р° д, получим

в

—«1 = К^р + К° aVs + /з = — «1 ^,р). (13)

Пользуясь (12), найдем

К^р = —К1К-1(« + V),

КХК-1 = к°1Ко(кК°)-1 = к°1 к-1 = Ь^) е (0,1),

следовательно, представлению (13) можно придать эквивалентную форму

—V1 = KоaVs — Ь« + Р, Р = /° — Ь/. (14)

Подстановкой (13) в уравнение неразрывности для первой фазы и (12) в соотношение (9) приходим к системе уравнений относительно .?,р

Зл«(К 'Чр + /)

д ( ° \ д ( о )+

— (рт) = ^гЛрз,т)+

дт{ Р0 — 1 дЬ

(15)

дЬ 1 дЬ

+ ^«(р^К^р + KоaVs + /°)), (16)

а при подстановке (14) в (6) - к эквивалентной системе относительно ■?,р,«

Зл«(К Vр + V)

дт( Р0 — О

дЬ

(17)

дЬ

+ Зл«(р\ (К° aV s — Ь« + Р)). (18)

Капиллярное давление и относительные фазовые проницаемости обладают свойствами

др — — —

—c■ < 0, к = к°1 + к°2 > 0, /го^^од = 0,

дs

и поэтому а(х^) > 0 при s е (0,1) и а(х, 0) = а(х, 1) = 0.

Таким образом, (15)-(16) (или (17)-(18)) представляет собой систему, состоящую из эллиптического уравнения для р(х, Ь) и вырождающегося при .? = 0,1 параболического уравнения для ,?(х, Ь).

Пусть снег занимает конечную область О с кусочно-гладкой границей Г = дО, Г = Г^У Г2, где Г1 - граница с непроницаемой поверхностью, Г2 - граница с воздухом.

Положим От = О х [0, Т], Г1т = Г1 х [0, Т], Г2т = Г2 х [0, Т], Гт = Г1Т и Г2Т. Граничные условия на Г 1т берутся в виде

V = (KVp + /') = «1 = (К1Vр + К° aV s + /о) = 0,

(х, Ь) е Г 1т, (19)

На Г2Т заданы приведенное давление и насыщенность

р = р°(х,Ь), s = s°(х,Ь), (х,Ь) е Г2т. (20)

Начальное условие для насыщенности имеет

вид

,?(х, 0) = .?°(х, 0), х е О.

(21)

Определение решения. Предположим, что все заданные функции к°®(.?), т(х,Ь), р^х^) и К (х), через которые выражаются коэффициенты уравнений (15)—(16), определены при всех (х,.?) и удовлетворяют условиям

(г) ||к° 1^), к°2^), ||с[о 1 < М^р^Цж^ < М;

(гг) М-1 < [т, к, (К°£, £)] < М;

(ггг) тг > 0;

(гу) 0 < (а,к° 1,к°2)^ е (0, 1);

а^од = к°1(0) = к°2(1) = 0.

Функции р°(х,Ь) и ,?°(х,Ь), входящие в граничные условия, предполагаются заданными при (х, Ь) е О х [0, Т].

(у)

дsn

дs

1,Пг

; 11^112,пт; ||pо, ^°Ц2

< М;

(ьг) 0 < sо(x, Ь) < 1;

Отметим, что в силу определения а, /, /° и Р имеем

1п

/ /

Р

кк01 кк02 кк01 кк02 кк02 кк01 кк02

< М° (М).

Ограниченные измеримые в О х [0, Т] функции .?(х,Ь),р(х,Ь) назовем обобщенным решением задачи (15)-(16), (19)-(21) если выполнены следующие условия:

а) 0 < .? < 1 почти всюду в От;

б) Vp е 1*2,&>(), aVs е Ь2(От), где при 0 < s < 1 (в этом случае а > 0) функция аV.? определяется формулой

а\^ =

дР^г, s)

д.9

в

Vu, и^) = J к°2ЬЗе, и е Ь2(От);

в) р = р°(х, Ь), s = s0(x, Ь) при (х,Ь) е Г2 х [0, Т];

г) для произвольных допустимых функций Ф,ф, таких что

ф(х,Ь) е Ш1(От), ф(х) е Wl(О), ф(х, Ь) = 0, ф(х, Ь) = 0 (х, Ь) е Г2т, ф(х, Т) = 0,

(

а

при почти всех і Є [0,Т] выполняются равенства

(22)

£2 — <<тз))і,ф)аі

(І тф а,

+ <щ, Уф)пь ■ (23)

Соотношения (22)—(23) получены путем домно-жения (15)—(16) на ф и ф и интегрированием по соответствующим областям. В дальнейшем будем называть задачу (15)-(16), (19)-(21) задачей I.

Принцип максимума. Продолжим каждую функцию в равенстве (23) вне промежутка [0, 1] согласно формуле

I<х, 0), в < 0,

I <х, в), 0 < в < 1,

I<х, 1), в > 1,

и, кроме того, заменим а<х, в) на а<х, в) =

аФ<х, в) + є, є > 0, а V и Ь на их усреднения и Ьь0 по х и в соответственно. Тогда (23) преобразуется к виду

£2 - <<тв)г,ф)Яь = ( ^

ї*<х,в)

(Рз \

* = Р0т,ф -

' Яь

Оценивая третье слагаемое снизу, а четвертое сверху, с учетом знаков остальных приходим к неравенству

г

Мтщ2,п <сI \\^тщ2,пЛт, к(0) = 0,

откуда следует, что я = 0 и, следовательно, s < 1.

Аналогично проводится доказательство .9 > 0. Рассмотрим вместо ф функцию в = шах{—s, 0}. Полагая у = —s, у = тах{у, 0} получим

1

1

2<т,У )п* + 2ІІУтУІІ2,п + <коаУу,у)п*+

Рз

+ (< РІ - Ь(Э - 0) тиУ)

— (KоaVs + Р, Vф')Qt — (CVhVs, ф)дь —

(> (ж — 0 т"ф)я; (24)

где с = (Ь°)«, Ь0 = Ь^.

Ограниченные измеримые функции

s(x,Ь),p(x,Ь) назовем обобщенным решением вспомогательной задачи 1е, если они обладают свойствами (б-г) определения обобщенного решения, в котором (23) заменено на (24).

Лемма 1. Пусть выполнены условия (г) — (Vг) обобщенного решения; s,p - обобщенное решение вспомогательной задачи 1е. Тогда почти всюду в О х [0,Т] для s выполнено следующее неравенство

0 < s < 1.

Доказательство. В уравнении (24) положим ф = к = шах^ — 1, 0} . Тогда

рз

(mts, s)nt + (mst, s)nt = (-^тг, s)nt — (KоaVs+

р1

/ р3

+ Р, Vs)nt — (cvhVs, s)nt — (Ь^ (-* — 1)тг, s)Qt,

р1

Из последнего равенства выводим (Р = 0 при

s > 1)

1

Отсюда, проводя аналогичные доказательству s < 1 оценки слагаемых, получаем, что у = 0 и, следовательно, s > 0.

Разрешимость. Сначала будет установлена разрешимость регуляризованной задачи 1е, а затем предельным переходом по е и Н в тождестве (24) доказана и теорема существования решений исходной задачи I. Задача I исследуется в общем случае, когда в области течения допускаются застойные зоны ^ = 0 или s = 1) и тем самым уравнение (16) может вырождаться.

Теорема. При выполнении условий (г) — (Vг) определения решения существует по крайней мере одно обобщенное решение задачи I.

Доказательство теоремы следует [2] и проводится в несколько этапов.

1. Построение галеркинских приближений вспомогательной задачи 1£ в виде

N

(Р,Ь = ^2 а!Офк(х) + sо(X,t), к=1

N р! (х,Ь = ^2 ь! ОФк (х)+ р°(/,г), к=1

где функции фк ,фк нормированы следующим образом

(фк, фг )п = 5к, (Vфк, Vфi)n = 5к.

Для определения неизвестных а!, Ь! получаем следующую нелинейную систему уравнений

2<тьВ2)п+І^т4'22,п+<КоаУв, +<сіїЧв, в)п* +

Оо

+ X1 - ЬНоX1-----1)ть, в)п* = °-

РІ

3,а^

йі

N

^а^ а^к + вк, а1^' <0)

3 = 1

0

N

j=i

в которой

ajk — (mtфj, ф ^

- (cfhУфj,фk)» - (KоaУфj, Уфk)»,

А= (1 »„«)» — (F,yn) » -

(bh0 ^p3 - ^ mt^k^j

-(~Ж'т,ф^ - (mtSо,Фk)»,

Hjk = (-y.j, y.k)»,

vk = - (л, y.k)

+

+

'dm( pi - 1 dt

фk I - (КУп, y.k)».

Для определения Ь^! получаем линейную алгебраическую систему уравнений, а для а! приходим к задаче Коши для системы обыкновенных дифференциальных уравнений.

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

»

»

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

1. Бэр Я., Заславски Д., Ирмей С. Физикоматематические основы фильтрации воды. — М., 1971.

2. Антонцев С.Н., Кажихов А.В., Монахов В.Н. Краевые задачи механики неоднородных жидкостей. — Новосибирск, 1983.

3. Blagovechshenskiy V., Eglit M., Naaim M. The calibration of avalanche mathematical model using field data // Natural Hazards. — 2002. - N. 2.

4. Naaim M., Gurer I. Two-phase Numerical Model of Powder Avalanche Theory and Application. // Natural Hazards — 1998. - N. 117.

5. Кучмент Л.С., Демидов В.Н., Мотови-лов Ю.Г. Формирование речного стока. Физикоматематические модели. — М., 1983.

6. Anderson E.A. Hydro-17 - Snow Model. NWSRFS Users Manual. Part II.2 // National Weather Service. NOAA. DOC. Silver Spring. MD.

— 1996.

7. Anderson E.A. Development and testing of snow pack energy balance equations // Water Resources Research. — 1968. - V. 4, №1.

8. Fowler A.C. An introduction to

mathematical modeling. Mathematical Institute. -Oxford, 2002.

9. Gray J.M.N.T. Water movement in wet snow // Philosophical Transactions: Mathematical, Physical and Engineering Sciences. — 1996. — V. 354, №1707.

10. Sellers S. Theory of water transport in melting snow with a moving surface // Cold Regions Science and Technology. — 2000. — V. 2000, №31.

11. Жумагулов Б.Т., Зубов Н.В., Монахов В.Н., Смагулов Ш.С. Новые компьютерные технологии в нефтедобыче. Алмата, 1996.

12. Папин А.А. Разрешимость модельной задачи тепломассопереноса в тающем снеге // Прикладная механика и техническая физика. — 2008.

- Т. 49, №4.

13. Папин А.А. Краевые задачи двухфазной фильтрации. - Барнаул, 2009.

14. Папин А.А., , Коробкин А.А., Гоман В.А. Движение воды и воздуха в тающем снеге // Известия АлтГУ. — 2012. - №1/1(76).