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

УДК 519.635.8:532.616

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ РОСТА ЛЕДЯНОГО ПОКРОВА В ПРЕСНЫХ И СОЛОНОВАТЫХ ВОДАХ*)

А. Ф, Воеводин, Т, Б, Гранкина

Введение.

В настоящее время много усилий направляется на изучение и освоение водных ресурсов северных территорий (водоемы, болота, реки) [1]. Данная работа нацелена на изучение режима льдообразования в слабопроточных водоемах [2].

В математическом отношении решение проблемы сводится к интегрированию уравнения теплопроводности в трех областях (рис. 1) с неизвестными подвижными границами (/, /з) и условиями сопряжения на этих границах, учитывая тепловой баланс и переменную температуру фазового перехода. Для решения используется метод «спрямления фронта», позволяющий решать уравнения в регулярных областях, соответствующих жидкой, твердой фазам вещества и слою снежного покрова. В результате этого преобразования в уравнениях теплопроводности появляется конвективное слагаемое, что наряду с образующимся погранслоем существенно влияет на выбор численного метода. Численно задача реализовывалась с помощью метода встречной прогонки.

Работа выполнена при проекта № 4.8 Программы фундаментальных исследований Президиума РАН.

© 2012 Воеводин А. Ф., Гранкина Т. В.

Рис. 1. Область решения задачи: 1 — вода, 2 — лед, 3 — снег 2. Постановка задачи

2.1. Основные уравнения. Рассмотрим слабопроточный водоем глубиной H (см. рис. 1). Распределение температуры и примеси в слое воды 0 ^ z ^ f (t), в образующихся ледовом f(t) ^ z ^ f (t) и снежном f2(t) ^ z ^ f (t) покровах описывается уравнениями теплопроводности:

U-Lw _ , О lw , .

Р'Ш Ср 0 , - /С-ш 0 Г, 5 (, J- J

сИ - т дz дС дС

дТгс 2 д Тгс /'о\

дТ8П д / дТ8П \

(4)

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

Тт, С) = 1000,08 + 0,0588Тад + 0, 797С - 0, 008Т^ - 0, 00325Т^С + 0,00013С2 + 0,0000477Т^ + 0,0000389Т?,С + 0,00000288С2ТШ - 0,00000006С3.

Плотность снега зависит от толщины снежного покрова и плотности выпадающего снега, определяется следующей зависимостью [4]:

p - p eKh(t)-z)

Hsu — yo^

Теплопроводность снега зависит от его плотности ksn = 2, 9 • 1О-pSn + 0,043.

Положение подвижных границ раздела сред вода — лед f (t), лед — снег f (t) и снег — атмосфера f (t) находим из соотношений

fl(t) — Iwa — Н lickp, кр — —,

pw

/2(t) = /l(t) + lic; f3(t)=h(t)+ls n,

где lw — толщина слоя воды, lic — искомая толщина слоя льда, lsn — толщина слоя снега.

Высота снежного покрова может быть известна из натурных измерений. Также ее можно определять из метеорологических данных, через количество свежевыпавших осадков l*sn (м) или восстанавливать через водный эквивалент W (мм):

isn=ln{1+hblL\ rm = o,mw^.

ь Po

Здесь Tw(t, z), Tic(t, z), Tsn(t, z) — температура воды, льда и снега соответственно, o С; C(t, z) — соленость, г/дм3; Ь = 1, 255; pw, pic, psn —

po

снега, кг/м3; а2 — коэффициент температуропроводности, м2/с; kw, kic, ksn — соответствующие среде коэффициенты теплопроводности, Вт/м°C; d — коэффициент диффузии соли в воде, м2/с.

2.2. Граничные и начальные условия. Для замыкания задачи ставятся краевые условия, определяющие внешнее воздействие на гра-

z

тепла и примеси:

AT... dC

= 0

dz

= »• £

z=0 dz

z

или

Ты = соп

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

На подвижной границе между водой и льдом г = /(Ь) должно выполняться условие сопряжения:

1) классическое условие Стефана, описывающее тепловой баланс:

дТгс

VdliC _ , дТи}

гс ТГ — 'ни о

оЬ дг

— кг,

дг

2) баланс массы растворенного в воде вещества:

л оь

дг

(5)

(6)

3) равенство температур сред и условие для температуры замерзания [5]:

Тгсе = Тш = Т/, Т/ = Т * - 7 С Л. (7)

Здесь Т* — температура замерзания чистой воды, Т* = 0°С; Т/(Ь) — фазовая температура (температура замерзания водоема); С/(Ь) — значение примеси па границе раздела фаз; 7 — равновесный коэффициент распределения примеси, зависящий от химических параметров растворенного в воде вещества [6]; А — скрытая теплота кристаллизации воды, Дж/кг.

На подвижной границе между льдом и снегом г = /(Ь) задаем равенство тепловых потоков и температур:

дТгс

дг

— ко

дТ

дг

Т

Тг

гси=/2

= Т.

ей |г=/2 •

(8)

2= Л

На верхней границе между снегом и атмосферой г = / (Ь) задаем внешнюю температуру Та(Ь), равную температуре воздуха па высоте 2 м над поверхностью снежно-ледового покрова.

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

Начальное значение концентрации примеси постоянно: C = Cq- Температура выпадающего снега постоянна и равна атмосферной.

Одной из особенностей подобных задач является наличие в начальный момент времени только одной фазы — жидкой. Будем считать, что при малых временных параметрах (t « At) движение фронта замерзания подчиняется следующему закону:

hc(t) = aVt, Vf = ^j-,

Vf — скорость движения фронта, a — некоторая константа, вычисляемая по формулам, предложенным А. Н. Тихоновым и А. А. Самарским [7], зависящая от теплофизических свойств жидкости, начальной температуры жидкости и внешней температуры па границе z = H (атмосферной температуры).

3. Метод численного решения

Для решения воспользуемся методом «спрямления фронта» [8]. Перейдем к новым координатам:

0<6<li 6 = AW-/'-!(«) ("=1'2'3):

/о (t) соответствует нижней границе z = 0.

Запишем систему уравнений (1)—(4) в новых координатах: dT k d2T dT

2 ^-t-w __u J-w UJ-w /l/\

w o, — oT2 Vw ' V /

dt pw cp d£i

/2 8C _ nd2C 8C /ч

j2 ^^ic _ 2 ^Tic QTic f

ic~dT~a,ic~dlТ~ЩсЖ' 1 j

dT k d2T dT

2 и ± sn _ гь8П u ± sn и ± sn f

где

vw = vic = (kp - (,2)hc—r-,

dt dt

dL„

dls

= (kp - 1%п-тг -Ь—тг

1 dls

dt dt ' psncp d£3 '

Здесь vw, vic, vsn — скорость изменения температуры и примеси соответственно.

Условия сопряжения на подвижных границах примут вид

dlic lit

1

A pi

k dT

ww

kpCr

kic dTic

lw d£i

dlic

kic dTic

«1=1 lic db

d dC

«2 = 0

licdb

dt

«2 = 1

lw d^i

k dT

™sn sn

«1=1

lsnd£3

(5')

(60 (80

«

4. Разностные уравнения

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

un+1 - un

= К-

n^1 - 2ИП+1 + пП+}

i

i+i

т

(vi - |vi

n

1

2|vi|nun+1

-v

"i-1

2 hr

которые перепишем следующим образом:

Aiun+l - Ciun+1 + BU^1 = -Fi,

i = 2,...,Nr

(9)

где Лг, Бг, Сг, определяются однозначно в зависимости от искомых функций; Мг — количество узлов сетки; г = (1, 2, 3) — номер области, ВОда — лед — снег, соответственно; Нг — шаг по пространству; т — шаг по времени; п — номер шага по времени; К — соответствующий для уравнений (10^(40 коэффициент; и — соответствующая для уравнений (10^(40 искомая функция; V — соответствующая для уравнений

(1')-(4') скорость. Система уравнений имеет первый порядок аппроксимации по пространству и времени.

В каждой из областей строилась равномерная сетка. Шаг по времени брался постоянным. Для численной реализации задачи использовался метод встречной прогонки [9].

В области 1 решения для Тш и С будем искать в виде

<+1 = , г = 1,..., N, (10)

где прогоночные коэффициенты считаем по рекуррентным формулам

из,с _ _^г_ о11>,с _ -^-гРг ^

аг+1 ~ „ ии,с л ■> Уг+1 ~ „ ии,с л 1

С - аг' Аг т С - аг' Аг согласно граничным условиям (6), (7): а^'с = 1, = 0В областях 2, 3 осуществляем сквозную прогонку, решение Тгс, Тзп ищем в виде

= , * = 1, (П)

здесь коэффициенты

л. Ев ' + Р

пгс.вп ~ 1 г

вг ' =

^ гС'Зп р 1

хг+1 Ег Сг — аг+1 Ег

согласно граничным условиям: а|Пг+1 = 0, вгПг+г = Та- Коэффициенты аг]Сг+1, в/Сг+1 Для расчета температуры в слое льда определяем через прогоночные коэффициенты аЩ в{п ■ Индексы ш, ге, зп соответствуют слою воды, льда, снега.

Соотношения (10) при г = N и соотношение (11) при г = 1 далее будут использоваться в условиях сопряжения для вычисления потоков тепла па границе я = /1(2), которые аппроксимируются в первой области

дТ„:

дЦ дС

_ Т1 ~ ТМГ _ Г/(1 - аЪг+1) ~ Рмг+1

~ Ъ ~ Л ''

1 ш ,1гш

^ С/ - Сдгг _ С( 1 - а%г+1) - /Здгг+1

~ Ъ ~ Л '

1 1 ш 1 ш

во второй области

дТ>

д£

Т2 - _ Т.г(а2 - !) + 02

«=0

Ъ-

Ъгг

Ъгг

Эти выражения для потока тепла и примеси на границе подставим в условия сопряжения (5'), (6'). Тогда получим систему из трех уравнений

V) = АтТ/ + Вт, V) С/ = Ас С/ + Вс, Т/ = Т» - 7С/,

из которых следует квадратное уравнение относительно С/

7Ат С) + (Ас - Вт) С/ + Вс = 0.

В результате решения этого уравнения получаем один удовлетворяющий физическим условиям корень [2]

С)

Вт -Ас + \](Ас - Вт)2 - 47АТВС

27 Ат

где

Ат

1

Л р.

Вт -

/ Ъ

<"Ш "и 1

(1 - а " к

'

1)

/Ъ

кг

-(агг -1)

Ас

Лр

¿(1 - + 1 Ь 1 Ь

Рдгг+1 - 7 7 12

Вс

¿в

'Дгг+1

Ь 1

' ' р

Получив значение примеси на границе, из уравнения (7) находим температуру фазового перехода, соответствующую температуре жидкой и твердой фаз на границе. Зная граничные условия, по (10) и (11) восстанавливаем значения температуры и примеси на новом слое. Тем самым заканчивается расчет на одном шаге по времени.

Рис. 2. Расчет толщины снежно-ледового покрова для Новосибирского водохранилища, 1976-1977 гг.

Z, м

Рис. 3. Расчет толщины снежно-ледового покрова для оз. Яркуль,1999-2000 г.

5. Результаты расчетов

Требуемые входные данные: шаг по времени т; Nr (r = 1,2,3) — количество узлов сетки в каждой области; H — глубина водоема; Ta(t) — массив метеоданных или функция, определяющая атмосфер-

Рис. 4. Расчет толщины снежно-ледового покрова для оз. Яркуль, 2002 2003 гг.

ную температуру; 1зп или 1*п — высота снежного покрова или суточные измерения осадков в водном эквиваленте, или измеренная толщина све-жевыпавшего снега; расчетное время; требуемые начальные условия — распределение температуры и примеси в средах.

Для расчета были выбраны следующие объекты: пресное Новосибирское водохранилище и минерализованное озеро Яркуль Чановской системы озер. Для определения начальной даты расчета анализировались температура воздуха и скорость ветра за несколько суток, так как общепринятое условие начала ледостава — среднесуточная температура воздуха Та < —50С и скорость ветра меньше 5 м/с.

На рис. 2 приведены результаты расчетов, проведенных по данным метеопоста Ордынское, и сопоставление с натурными измерениями. Средняя глубина водохранилища в районе метеостанции 10 м, средняя скорость течения 0,2 м/с. Расчет толщины ледяного производился в период с ноября по март 1976-1977 гг.

0.9

Чс

0.4

0.6

0.1

0.2

03

0.5

0.8

0.7

о У.......................,,,.,,,...............

7-Х1. 22-Х1. 7-ХИ. 22-ХП. 6-1. 21-1. 5-И. 20-П. 7-Ш. 22-Ш.

Рис. 5. Расчеты толщины ледового покрова для Новосибирского водохранилища, 1988-1989 гг.

Приведены результаты расчетов процесса формирования ледового покрова оз. Яркуль и сопоставление с натурными измерениями (рис. 3, 4). Входными данными служили метеоданные по метеостанции н.п. Купино зимы 1999-2000 гг. (рис. 3), 2002-2003 гг. (рис. 4). Начальная минерализация оз. Яркуль задавалась равной С = 5 г/дм3

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

[10].

ЛИТЕРАТУРА

1. Ямал — проблемы развития. Институт проблем освоения Севера СО РАН. Тюмень, 1993.

2. Воеводин А. Ф., Гранкина Т. В. Численное моделирование роста ледяного покрова в водоеме // Сиб. журн. индустр. математики. 2006. Т. 9, № 1. С. 47-54.

3. Tenth report of the joint panel on oceanographic tables and standards. UNESCO Technical Papers in Marine Sei., 1981. No. 36, UNESCO, Paris.

3. Винников С. Д., Проскуряков В. В. Гидрофизика. Л.: Гидрометеоиздат, 1988.

4. Васильев В. П., Максимов А. М., Петров Е. Е., Ципкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука. Физматлит, 1997.

5. Гороновский И. Т. и др. Краткий справочник по химии, 1987.

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

7. Вудак В. М., Гольдман П. Л., Успенский А. В. Разностные схемы с выпрямлением фронтов для решения многофронтовых задач типа Стефана // Докл. АН СССР, 1966. Т. 167, № 4. С. 735-738.

9. Воеводин А. Ф., Леонтьев П. А., Петрова А. Г. Термодиффузная задача о кристаллизации шара // Динамика сплошной среды, 1982. №5. С. 118-123.

10. Савкин В. М., Двуреченская С. Я., Сапрыкина Я. В., Марусин К. В. Основные гидролого-морфометрические и гидрохимические характеристики озера Чаны //Сиб. экологический журн. 2005. Вып. 2. С. 183-192.

г. Новосибирск

1 февраля 2012 г.