CC BY
52
УДК 551.578:536.421:517.95
Численное решение двумерной задачи движения воды и воздуха в тающем снеге*
В.А. Гоман, А.А. Папин, К.А. Шишмарев
Алтайский государственный университет (Барнаул, Россия)
A Numerical Solution of a Two-Dimensional Problem of Water and Air Movement in Melting Snow
V.A. Goman, A.A. Papin, K.A. Shishmarev Altai State University (Barnaul, Russia)
Показано численное решение двумерной задачи движения воды и воздуха в тающем снеге. В качестве математической модели используются уравнения сохранения массы для каждой фазы, уравнения двухфазной фильтрации Маскета-Леверетта для воды и воздуха, уравнение сохранения энергии для тающего снега и уравнение движения льда. На поверхности снежного покрова задаются насыщенность воды, температура воздуха (выше температуры плавления льда), атмосферное давление воздуха, скорость воды и воздуха, на границе контакта с промерзшим грунтом вода отсутствует, заданы температура воздуха (ниже температуры плавления льда) и давление. Приведен обзор моделей снежного покрова разного уровня детализации физических процессов. Рассмотрены балансовые модели, обсуждаются проблемы получения натурных данных. Описан алгоритм численного решения двумерной задачи движения воды и воздуха в тающем снеге. Задача сводится к системе из трех уравнений относительно температуры, «приведенного» давления и насыщенности водной фазы. Для полученной системы уравнений рассмотрена начально-краевая задача и построена конечно-разностная схема на основе метода переменных направлений. Проведены расчеты тестовой задачи, определены насыщенность и температура при заданных начальных приближениях и проведен графический анализ результатов. Для насыщенности воды установлено свойство конечной скорости распространения возмущений.
Ключевые слова: тающий снег, пористая трехфазная среда, фазовые переходы, численный расчет, метод переменных направлений.
БОТ 10.14258/^уа8и(2014)1.2-01
* Работа выполнена при финансовой поддержке государственного задания Минобрнауки РФ № 2014/2 и гранта РФФИ № 13-01-98016.
In this paper, a numerical solution of a two-dimensional problem of water and air movement in melting snow is presented. The equation of mass conservation for each phase, the equations of two-phase filtration of Musket-Leverette for water and air, equation of energy conservation for melting snow and equation of motion of ice were used as a mathematical model. Formulation of the problem is given with the following parameters: water saturation, air temperature (above the melting temperature of ice), atmospheric air pressure, water and air velocities on the snow surface, air temperature (below the melting temperature of ice), and pressure with no water at the frozen ground surface. Among other things, an overview of snow cover models for different detail levels of physical processes is presented. Balance models of snow cover are considered, and the problem of obtaining field data is discussed. Presented numerical solution of the two-dimensional problem of water and air movement in melting snow. The problem is reduced to a system of three equations for temperature, "reduced"pressure and saturation of a water phase. The obtained system of equations is considered to be an initial-boundary value problem, and the finite-difference scheme based on the alternating directions method is elaborated. Test calculations and validation are performed with saturation and temperature defined for given initial approximations. Graphical analysis of the obtained results is described. Also, finite perturbation velocity for water saturation is estimated.
Key words: melting snow, three-phase porous media, phase transitions, numerical solution, method of alternating directions.
1. Постановка задачи. Снег представляет собой пористую среду, твердый каркас которой составляют частички льда. В процессе таяния в пористой среде происходит совместное движение воды (г = 1), воздуха (г = 2) и льда (г = 3). В снеге происходят постоянные фазовые превращения, ведущие к перераспределению фазовых масс. Для описания процесса используются уравнения сохранения массы для каждой фазы, уравнения двухфазной фильтрации Маскета-Леверетта для воды и воздуха, уравнение сохранения энергии для тающего снега (в пренебрежении сублимацией и обменом массами между водой и воздухом) и уравнение движения льда [1, 2]:
Ot(Piai) + div(p°aiUi) = J] Iji, ¿ = 1,3,
j=i
3
iji = -iij, iij = o, i,j=1
(1)
1, 2,
P2 - Pi = Pc(si, 0),
2
E si
i=i
(2)
i=i
(J^PiCiai)— + (J^PiCiViWe
i=i
= div(AcV0) + v
dp3a3 dt '
(3)
a3p°3(^+u3^) = -gp03a3-^ +
+M3
a2 u3
dz2
- -T
j
du3 _ dz
(u3 - uj),
(4)
A
1 ~ аз
аз
~Рз,
где ai — концентрация г-й фазы; р0 — плотность г-й фазы; ui — скорость г-й фазы; iij — интенсивность перехода массы из j-й в г-ю составляющую в единице объема в единицу времени; pi - приведенная плотность, связанная с р0 и ai соотношением pi = aiр0; щ = msiUi — скорости фильтрации воды и воздуха; si, S2 — насыщенности воды и воздуха (ai = msi, а2 = ms2, аз = 1 — m); v3 = (1 — m)u3 — расход льда; K0 — тензор фильтрации; k01, k02 — относительные фазовые проницаемости воды и воздуха (koi = koi(si) > 0, koi|s =0= 0); pi — динамическая вязкость; pi — давление г-й фазы; pc — капиллярное давление; g = (0,0, —g) — вектор ускорения силы тяжести; 0 — температура среды (0i = 0, г = 1, 2, 3); ci = const > 0 — теплоемкость г-й фазы при постоянном объеме; v = const > 0 — удельная теплота плавления льда; Ac — теплопроводность снега
3
(Ac = ac + hcp"2, pc = P0ai, ac = const > 0,
i=i
hc = const > 0); A* — коэффициент пропорциональности; k — безразмерный коэффициент, характеризующий геометрическую форму фазы и меняющийся от 2 до 3.
В общем случае величины Iij, р0, в-ij должны быть заданы (истинные плотности, как правило, функции температуры и давления фаз). В частности, в работе [1] для замыкания системы (1)-(4) использовались гипотезы: ii2 = 0, I23 = 0, р0 = const, г = 1, 2, 3. Причем 1з1 считается функцией температуры I31 = I31 (0).
Близкие к приведенной выше постановке задачи рассматривались в [2, 3].
Описанию фазовых переходов посвящены работы [2-7].
В настоящее время все большее применение находят модели, описывающие взаимодействие снежного и ледового покрова. Единой модели ледового покрова, пригодной для прикладных задач, не существует, что является существенным препятствием для надежного прогнозирования поведения ледового покрова. Наиболее изученным является поведение ледового покрова под действием поверхностных или изгибно-гравитационных волн, когда ледовый покров моделируется упругой плавающей пластиной [8, 9].
Балансовые модели строятся на основе системы вертикально осредненных уравнений снежных процессов [10]. Как правило, эта система включает в себя описание изменения высоты снежного покрова в зависимомсти от времени, содержание льда, воды, плотность снега, учитывает таяние снега, сублимацию, замерзание талых вод и метаморфизм снежного покрова [11].
Большое количество работ посвящено сбору и анализу натурных данных [12-28].
2. Численное решение. В данном разделе описан численный алгоритм решения двумерной задачи движения воды и воздуха в тающем снеге. За основу взята система уравнений (1)-(3), описывающая снег как многофазную пористую среду, но без учета скорости движения льда (v3 = 0). Следуя [29], введем приведенное давление p
1
dpc k02
P = Pi —
d£
M2
( hm. 1
V Ml
Рассматриваемая система уравнений сводится к системе из трех уравнений относительно приведенного давления р, насыщенности водной фазы и температуры в
д_
т
(р?т5) = § (р0т) +
div(K Vp + f) =
т I -ту - 1 Pi
(5)
JL
dt
+ K0aVs + /0)),
3
1
3
3
(Eli Р°сгаг) |f + (£ti РМ) V0
= div(AcW) +
(7)
где
/ = к J + *2vPc + g кгРи
К = kK0, /г = fcoi + k02, k0i = —;
Mi
7dpc k02 ds k
Z- (pirns) = % (p°3m) +
at
(8)
+div(p?(KQaVs - bv + F)),
KiK-1 = ]tQ1KQ(kKQ)-1 = /CQlk-1 = b(s),
F = fo - bf.
С учетом равенства (5) уравнение (8) запишем в виде
ds fp$M „ , А дв
+div(KQaVs + f) - b's(Vs)v.
mrn =
p(x, H,t)= P2(x,H,t) - Pc(0,в+)-1
(11)
_ f dpc fc02 J (fcoi+fco2)
de,
s(x, 0,t)=0, в(х, 0,t) = в-, p(x, 0,t) = p2(x, 0,t) - pc(0,в-)-1
_ f dpc k02 It.
J (fcoi+fco2)
(12)
ds = 0, дв = 0, dp
dx dx dx
ds = 0, дв = 0, dp
dx dx dx
= 0 (x = L); = 0 (x = 0).
(13)
(14)
Кг = Коког, г = 1, 2.
Система уравнений (5)-(7) получена при следующих гипотезах: истинные плотности р0 постоянны, фазовые переходы воды и льда в воздух отсутствуют, пористость т считается заданной функцией температуры т = т(9).
Правая часть уравнения (6), с учетом представления для скорости V = -(К ^р+/) = VI +«2, может быть переписана относительно в, V [29]
Для численного решения системы (5), (7), (9) используется метод переменных направлений. В области И = {0 < х < Ь, 0 < у < Н}, зададим начальные условия для в и в
в(х, у, 0) = во(х, у), в(х, у, 0) = во(х, у), (10)
и краевые условия
в(х,н,г) = 0, в(х,н,г) = в+,
Начальное состояние снега задается распределениями насыщенности воды и температуры. На верхней и нижней границах насыщенность равна 0, а температура имеет значения в+ (больше плавления льда) и в- (меньше замерзания воды) соответственно. Область, занимаемая снегом, определяется по значению пористости (если пористость 0 < m < 1, то в данной точке присутствует снег).
Численное решение проводится в предположении Kq = const, а относительные фазовые проницаемости определяются по формулам
k01 = sni, k02 = (1 - s)n2.
Пористость задается следующим образом [30]:
1, в > в+, m = { m- + т1(в - в-), в- < в < в+, т-, в < в-,
где П1, П2, т(в-) G (0,1), m1 = (1 - т-)/(в+ --в-) — заданные параметры.
Капиллярное давление является заданной функцией насыщенности и температуры pc(s, в) = = Y (Фо(в).
Численные значения параметров задачи перечислены в таблице 1.
Рассмотрим вычислительную сетку на области Q. Разобьем область на N1 частей по x и на N2 частей по у. Тогда
xi+1 - xi = Дx = L/N1 = l,
У3+1 - Уз = Ду = H/N2 = h.
Запишем систему (5), (7), (9) в следующем виде:
дв
W(s, в)— = div(XcVe) - V(v)Ve, (15)
о \ дв div(KVp + f)= (16)
ds дв
т— = B(s) — +div(K0aVs+F) + U(s,v)Vs, (17)
dt
где
W(s, в) = pQciai + vp03m^ ,
V (V) = (¿/
/pQcifVi
B(s) = ( g (1 + b) - b - s ) me, U(s, v) = -b'sv.
1
Рассмотрим разностную схему переменных направлений для определения температуры
W3(C* - *S)/(r/2) = А- , - е':р)/1
п+ -
Г.П+ -
i+Ai
\п (Qn+2 Дп+2\//2_|_ \п /пп Qn\/h2
т
+ ((Fx )n+ij — (Fx)n-ij )/2/+((Fy j — (Fy )nj-i )/2h+
I / 77- ^H" 2" I OTl 2 I Tl 2 \ /7
+ (nx2 si+1j + Px2sij + 7x2si-1j )/1 + + «2 sj+1 + вп2 sj + 7y2sj-1
переход со слоя n + т> на слой п + 1 описывается следующим образом:
~К.. 1 ~ &ij-i +
13—
+ Wi 0ij+i + eni°ij- + 7^10j-i )/h
переход со слоя n + ^ на слой n + 1 описывается следующим образом:
-А?.. ЛС -С- )/^
+ (nyi0ij+i + + 7yni ji )/h,
ПП1 = (|(Vx )n | — (Vx )n )/2, eni = —|(Vx)?,-1,
7xni = (l(Vx )n | + (Vx )n )/2.
' ij I 1 V x 7 i j ,
Разностная схема для определения p имеет вид
+ ((/x )n+ij — (/x)n-ij )/21 + ((/y ) j+i — (/y ) j-i )/2h—
(р0/р1 — 1) m,(j — 0Jj-)/t,
где Т1 - итерационный параметр. Переход со слоя п + 7} на слой п + 1 описывается следующим образом: 1
=
= (РЙ -Р?5 (С3 )Л2+
+(РЙ1! -ОМ2 ■-■(р^1 -РЯ )м2+
+ ((/x )n+ij — (/x)n-ij )/21 + ((/y ) j+i — (/y ) j-i )/2h—
— (р3/р0 — 1) m,(j — 0g-)/t. Разностная схема для определения s имеет вид
- s?.)/(t/2) = щ^1 ~
n+
n+7
П+7
'i+ij^i+lj
j ij
'i-ij'
Msp1 - sl+h/(r/2) = -
n+i
n+
(p n+1
n+i
+ ((Fx )?+ij — (Fx)?-ij )/21+((Fy )ij+i — (Fy )ij-i )/2h+
I / Tl iT-H" 2 I OTl n+ 2 I Tl 712 \ / 7 I
+ (nx2 si+1j + Px2sij + 7x2si-1j )/f +
+ + ву2 + )/Л,
ПП2 = (1(их| - (их)/2, вх2 = -|(их
7хп2 = (|(их 1 + (их )П )/2.
Алгоритм вычисления следующий: сначала находятся значения температуры, затем соответствующие давление и насыщенность. Значения на слое находятся методом прогонки по ж, а значения на п + 1 слое — прогонкой по у.
Рис. 1. График изменения температуры
Рис. 2. График изменения насыщенности воды
2
При расчетах использовались следующие значения: р0 = 103 кг м-3, р0 = 1.2928 кг м-3, р3 = 0.925 х 103 кг м-3 - истинные плотности воды, воздуха и льда, К0 = 10-9 м2 - тензор фильтрации, р1 = 0.17865 кг м-1 с-1, р2 = 1720 х 10-8 кг м-1 с-1 - коэффициенты динамической вязкости воды и воздуха, с1 = 4.18 х 103 кг м2 с-2,
с2 = 29.2 х 103 кг м2 с-2, с3 = 2.04 х 103 кг м2
-2
с-2 - теплоемкости воды, воздуха и льда при постоянном давлении, V = 6009/18 х 10-3 Дж кг -удельная теплота плавления льда, в- = 268.15 К, в+ = 273.15 К.
На рисунках 1-2 представлены графики примеров расчетов изменения температуры и насыщенности для тестовой задачи с заданным полем
скоростей. Вычисления проводились на равномерной сетке с 25 узлами как по х, так и по у. Длина и высота снежного покрова равны 1 м, шаг по пространственным переменным составил 0.0416. Задавались начальные распределения температуры и насыщенности во и во. На верхней и нижней границе задавались краевые условия первого рода.
Численный расчет насыщенности с временным шагом т = 0.1 представлен на рисунке 2. За конечное время значение насыщенности воды становится равным нулю (вода замерзает).
Авторы статьи признательны С.С. Кузикову за обсуждение задачи и конструктивные замечания.
Библиографический список
1. Папин А.А., Коробкин А.А., Гоман В.А. Движение воды и воздуха в тающем снеге / / Известия Алт. гос. ун-та. — 2012. —№1/1(73).
2. Трофимова Е.Б. Математическая модель снежного покрова как многофазной среды // Тр. IV Всесоюзн. гидролог. съезда. — 1976. — Т. 6.
3. Денисов Ю.М. Перенос тепла и влаги в почве (неподвижной пористой среде) // Тр. САНИГМИ. — 1968. — Вып. 39(54).
4. Lou D., Hammond W.D. Heat and Mass Transfer for Ice Particle Ingestion Inside AeroEngine // Journal of Turbomachinery. — ASME. — 2011. — July. — Vol. 133.
5. Tantserev E., Cristophe Y. Galerne, Podladchikov Y. Multiphase flow in multi-component porous visco-elastic media // The Fourth Biot Conference on Poromechanics. — 2009.
6. Sellers S. Theory of water transport in melting snow with a moving surface // Cold Regions Science and Technology. — 2000. — №31.
7. Gray J.M.N.T. Water movement in wet snow // Phil. Trans. R. Soc. Lond. A. — 1996.
8. Коробкин А.А., Папин А.А., Шишма-рев К.А. Поведение ледового покрова канала под действием поверхностных волн // Известия Алт. гос. ун-та. — 2012. — №1/1(73).
9. Коробкин А.А., Папин А.А., Шишма-рев К.А. Аналитическое и численное исследование квазиизотермической задачи взаимодействия ледового покрова канала и поверхностных волн // Известия Алт. гос. ун-та. — 2012. — №1/2(73).
10. Gelfan A.N., Pomeroy J.W., Kuchment L.S. Modeling Forest Cover Influences on Snow Accumulation, Sublimation, and Melt // Journal of Hydrometeorology. — 2004. — October. — Vol. 5.
11. Kuchment L.S., Romanov P. Gelfan A.N., Demidov V.N. Use of satellite-derived data for characterization of snow cover and simulation of snowmelt runoff through a distributed physically based model of runoff generation // Hydrol. Earth Syst. Sci. - 2010. - October. - Vol. 10.
12. Koskinen J.T., Pulliainen J.T., Luojus K.P., Takala M. Monitoring of Snow-Cover Properties During the Spring Melting Period in Forested Areas // IEEE Transactions on Geoscience and Remote Sensing — 2010. — January. — Vol. 48.
13. Marechal N., Guerin E., Galin E., Merillou N., Merillou S. Heat Transfer Simulation for Modeling Realistic Winter Sceneries // Eurographics. — 2010. — №2. — Vol. 29.
14. Farinotti D., Magnusson J., Huss M., Bauder A. Snow accumulation distribution inferred from time-lapse photography and simple modelling // Hydrological Processes. — 2010. — March. — Vol. 24.
15. Reusser D.E., Zehe E. Low-cost monitoring of snow height and thermal properties with inexpensive temperature sensors // Hydrological Processes. — 2011. — December. — Vol. 25.
16. Woods R.A. Low-cost monitoring of snow height and thermal properties with inexpensive temperature sensors // Advances in Water Resources — 2009. — Vol. 32.
17. Sturm M., Holmgren J., Liston G. A seasonal snow cover classification-system for local to global applications //J Clim. — 1995. — Vol. 8(5).
18. Floyd W., Weiler M. Measuring snow accumulation and ablation dynamics during rain-on-snow events: innovative measurement techniques // Hydrological Processes. — 2008. — October. — Vol. 22.
19. Green R.O., Painter T.H., Roberts D.A., Dozier J. Measuring the expressed abundance of the three phases of water with an imaging spectrometer over melting snow // Water Resources Research. — 2006. — Vol. 42.
20. Luce C.H., Tarboton D.G., Cooley K.R. Sub-grid parameterization of snow distribution for an energy and mass balance snow cover model // Hydrological Processes. — 1999. — Vol. 13.
21. Hock R., Holmgren B. A distributed surface energy-balance model for complex topography and its application to Storglaciaren, Sweden // Journal of Glaciology. — 2005. — Vol. 51. — № 172.
22. Prevost M., Barry R., Stein J., Plamon-don A.P. Snowmelt modeling in a balsam fir forest: comprasion between an energy balance model and other simplified models // Can. J. For. Res. — 1991. — Vol. 21.
23. Munro D.S. A surface energy exchange model of glacier melt and net mass balance // International Journal of Climatology. — 1991. — Vol. 11.
24. Gallee H. Air-snow interactions and the surface energy and mass balance over the melting zone of west Greenland during the Greenland Ice Margin Experiment // International Journal of Climatology — 1997. — June. — Vol. 102. — № D12.
25. Mundy C.J., Barber D.G., Michel C. Variability of snow and ice thermal, physical and optical properties pertinent to sea ice algae biomass during spring // Journal of Marine Systems. — 2005. — November. — Vol. 58.
26. Theriault J.M., Stewart R.E., Milbrandt J.A., Yau M.K. On the simulation of winter precipitation types // Journal of Geophysical Research. — 2006. — September. — Vol. 111.
27. Finger D., Pellicciotti F., Konz M., Rimkus S., Burlando P. The value of glacier mass balance, satellite snow cover images, and hourly discharge for improving the performance of a physically based distributed hydrological model // Water Resources Research — 2011. — July. — Vol. 47.
28. Takala M., Pulliainen J., Huttunen M., Hallikainen M. Detecting the onset of snow-melt using SSM/I data and the self-organizing map // International Journal of Remote Sensing. — 2008. — February. — Vol. 29. — №3.
29. Коробкин А.А., Папин А.А., Хабахпаше-ва Т.И. Математические модели снежно-ледового покрова. — Барнаул, 2013.
30. Папин А.А. Краевые задачи двухфазной фильтрации. — Барнаул, 2009.
