Введение
Актуальность темы. В связи с проблемой наблюдающегося глобального изменения климата все более актуальными становятся вопросы, связанные с изучением последствий воздействия изменения метеорологических условий на гидрологические и экологические процессы в водоемах и водотоках при их круглогодичной эксплуатации. Большое значение приобретает рациональное использование водных ресурсов в зимнее время, когда очень важным становится обеспечение достоверности прогноза процесса формирования ледяного покрова и его прочности, позволяющее определить начало и сроки эксплуатации ледовых трасс, проложенных по поверхности водных объектов.
Разработкой методов расчета толщины ледяного покрова занимались многие исследователи: Я.Л. Готлиб (1983), Л.Г. Шуляковский (I960,1969, 1972), В.В. Пиотрович (1968, 1969, 1970), Р.В. Донченко (1971, 1983), Т.В. Одрова (1979), Г. Эштон (1978,1986), А.И. Пехович (1983), Г.А. Трегуб (1994, 1997), И.Н. Шаталина (1990), В.М. Мишен (1997) и многие другие. К настоящему времени предложено большое количество формул и расчетных приемов для определения толщины ледяного покрова пресноводных водоемов и водотоков, опирающихся на основные методы расчета: эмпирический и теоретический методы.
Эмпирический метод основан на отыскании эмпирических связей толщины льда и отдельных факторов, определяющих изменение толщины ледяного покрова. В этом случае расчетные эмпирические соотношения получены на основании относительно тесной корреляции между некоторыми температурными характеристиками и толщиной льда. Как прави-
ло, эти характеристики привязаны к конкретному региону.
Теоретический метод основан на интегрировании исходных дифференциальных уравнений, описывающих физическую сущность нарастания толщины льда.
Основой для получения эмпирических формул послужило уравнение для подсчета теплового баланса, так называемое условие Стефана. В 1889 г. Й. Стефан в работе о промерзании грунта поставил и решил следующую задачу (Рубинштейн Л.И., 1967).
Среда, находящаяся в двух фазовых состояниях (жидком и твердом) и проводящая теплоту исключительно посредством теплопроводности, заполняет полупространство х > 0. В начальный момент времени вся она находится при постоянной температуре Ti > 0°С. На плоскости х = 0 поддерживается постоянная температура Т\ < 0°С, под воздействием которой происходит кристаллизация, протекающая изотермически при температуре Т = 0°С без переохлаждения и с пренебрежимо малым объемным эффектом.
Требуется определить температуры U\{x,t) и U2(x,t) твердой и жидкой фазы и положение х = y(t) границы раздела фаз.
Подсчет теплового баланса приводит, как показал Стефан, к условию
dt V дх дх) x=y{t)
Здесь Л - скрытая теплота кристаллизации, отнесенная к единице массы;
р - плотность образующейся фазы; к\ и кч - коэффициенты теплопроводности соответственно твердой и жидкой фазы.
Эмпирические формулы для определения толщины ледового покрова, полученные на основе решения задачи Стефана, имеют общий вид :
где h^ - толщина льда, см; '?(—62) - сумма средних суточных температур воздуха (°С) на высоте 2м над уровнем водоема за расчетный период времени; а - эмпирический коэффициент, определяемый по данным непосредственных наблюдений и отражающий в среднем те условия, которые имели место в период наблюдений (температуру воды, высоту и плотность снежного покрова, глубину водоема, скорость подледного течения воды). Однако ввиду различия этих факторов даже для отдельных участков водоемов и недостаточной продолжительности наблюдений указанный параметр существенно меняется. Поэтому все подобные формулы носят локальный характер. Подробный обзор работ, описывающих получение эмпирических формул различными авторами для конкретных водных объектов представлен в монографии Д.В. Козлова (2000).
Очевидно, что в подобных эмпирических формулах не находят в полной мере своего отражения такие важные факторы, как мощность и физические свойства снежного покрова, интенсивность теплопотока из водной массы, метеорологические условия. Тем не менее, снег за счет своих теплоизоляционных свойств оказывает особое влияние на рост льда, уменьшая влияние отрицательных атмосферных температур.
В своих исследованиях норвежский ученый О. Дэвик (Винников С.Д., Проскуряков Б.В., 1988) к процессу льдообразования подошел с позиций физики, также основываясь на уравнении Стефана. В предложенной им формуле учитывается теплообмен на верхней поверхности снежного покрова, теплофизические свойства снега. Его подход стал классическим и получил широкое развитие. В дальнейших разработках многие ученые уделяли больше внимания процессу теплообмена снежно-ледовой поверхности с атмосферой (В.В. Пиотрович), теплообмену с дном водоема, водным потокам, начальной стадии формирования льда, образованию шуги
(Л.Г. Шуляковский, Р.В. Донченко, С.Д. Винников, Б.В. Проскурякрв, Г.А. Трегуб, И.Н. Шаталина).
На основе обощения этих разработок сформированы методические рекомендации, выпущенные научноисследовательскими и проектными институтами ВНИИГ, ГГИ, Ленгидропроект.
С развитием вычислительных методов математического моделирования и компьютерных технологий последнее время стали чаще обращаться к теоретическим методам, более точно описывающим физические явления процесса льдообразования. Появилась возможность численно решать сложные системы дифференциальных уравнений (Цибульский В.Р. и др., 1993: Васильев О.Ф., Бочаров О.Б., Зиновьев А.Т., 1999: Васильев В.И. и др., 1997: Дебольская Б.И., 2003: Прокофьев В.А., 2004), описывающих процессы льдообразования. Не остались без внимания важные задачи прогнозирования ледотермических режимов водоемов-охладителей тепловых и атомных электростанций, влияние их эксплуатации на температурный и ледовый режим рек (Квон В.И., Филатова Т.Н.). Над задачами моделирования ледового режима бьефов Красноярской ГЭС успешно трудились В.М. Белолипецкий, С.Н. Генова, В.Б. Туговиков, Ю.И. Шо-кин (1991, 1993).
На юге Западной Сибири имеется значительное число минерализированных озер. Как показали проведенные исследования, даже слабая минерализация оказывает существенное влияние на процесс замерзания. А именно, при формировании ледового покрова в минерализированном водоеме в результате образования практически пресного льда в воде перед фронтом кристаллизации формируется слой с повышенным содержанием химических веществ, что существенно влияет на температуру фазового перехода, а тем самым и на весь процесс формирования ледо-
вого покрова (Власов В.П., 1968, 1976), что необходимо учитывать при разработке математической модели и выполнении соответствующих расчетов.
Процесс замерзания морских вод достаточно подробно изучен, описан в учебниках по океанологии и статьях (Смирнов Г.Н., 1974: Федоров К.Н., 1976: Саркисян А.С, Демин Ю.Л., Бреховских А.Л., Шаха-нова Т.В., 1986: Багно А.В., Гаращук Р.В., Залесный В.Б., 1996: Яковлев Н.Г., 2003). Гидродинамические и гидрохимические режимы морских вод отличается от режимов соленых озер, которые по солевому составу чаще являются солоноватыми. Поэтому способы описания процесса образования морских льдов не применимы к озерным.
Во всех этих работах, основанных как на эмпирических, так и на теоретических методах, снежный покров в процессе льдообразования принято учитывать через эквивалентный слой:
и. h — h _i_ h гсе
"¦ecu — "¦fee ~r IT'snow т, '
где hecv - эквивалентный слой снежно-ледового покрова, hsnow - высота снежного покрова на льду, к{се и ksnow - коэффициенты теплопроводности льда и снега, соответственно. При этом предполагается линейное распределение температуры в слое льда и снега, что далеко не соответствует действительности.
Снежный покров по высоте имеет неоднородную структуру. Соответственно, изменяются его физические свойства, такие как плотность и теплопроводность. Плотность снега меняется со временем в зависимости от усадки, ветрового воздействия, плотности свежевыпавшего снега. Это приводит к нелинейному распределению температуры в слое снега.
В работе Л.И. Рубинштейна (1967) подробно изложена история развития "проблемы Стефана". Задачи Стефановского типа, помимо изучения
ледового режима водоемов, получили широкое распространение в других областях науки. Они успешно используются при моделировании процессов фазовых переходов в бинарных смесях (производство полупроводниковых материалов, очистка методом направленной кристаллизации, металлургическое производство), подробно описанных Н.А. Авдониным (1980), Рубинштейном Л.И., 1967 A.M. Мейрмановым (1986), при описание роста кристаллов (Карслоу Г., Егер Д., 1964: Любов Б.Я., 1969, 1975: Флеминге М., 1977), при описании диффузионного переноса вещества в зоне реакции и.т.д.
Известны различные численные методы решения задачи Стефана, такие как метод сквозного счета, используемый в работе (Бондарев Э.А., Васильев В.И., 1984: Васильев В.И. и др., 1997), методы с явным выделением фронта (Будак Б.М. и др., 1966: Воеводин А.Ф., Леонтьев Н.А., Петрова А.Г., 1982: Петрова А.Г., 1983: Белолипецкий В.М. и др., 1993; Овчарова А.С., 1994, 1995, 1999: Журавлева Е.Н., 1998). В работе (Бондарев Э.А., Попов Ф.С., 1989) приводятся оценки точности наиболее употребляемых приближенных методов решения задачи Стефана.
Целью работы является разработка математических моделей и эффективных численных методов для исследования процесса формирования ледового покрова водоемов различной минерализации в условиях Западной Сибири.
Основные задачи формулируются следующим образом:
• разработка математических моделей и их численная реализация для исследования процесса формирования ледового покрова в водохранилищах при наличии снежного покрова с учетом минерализации воды;
• разработка численного метода решения двухфазной термодиффузи-
онной задачи Стефана с температурой кристаллизации, зависящей от концентрации примеси;
• проведение расчетов ледотермического режима реальных объектов: Новосибирское водохранилище, озера Чановской системы;
• сравнительный анализ результатов расчетов по разработанным моделям и известным инженерным методикам и сравнение с натурными измерениями.
В работе основное внимание уделено разработке математических моделей для исследования процесса динамики ледяного покрова водных объектов с учетом важных физических факторов: наличие снежного покрова и зависимость температуры фазового перехода "вода-лед"от концентрации солей в воде, построению адекватных разработанным моделям численных методов и созданию реализующих эти методы программ для ПЭВМ, созданию численных моделей для исследования ледотерми-ческих процессов на реальных объектах Сибири (Новосибирское водохранилище, озера Чановской системы).
В первой главе пункт 1.1. посвящен описанию и разработке алгоритма численного решения одномерной двухфазной классической задачи Стефана, моделирующей кристаллизацию чистого вещества (Флеминге М., 1977; Рубинштейн Л.И., 1967). Рассматриваются области Qi(t) = {О < z < f(t)} и ?ls(t) = {f(t) < z < Н}, занятые соответственно жидкой и твердой фазами данного вещества с гладкой границей раздела z = f{t). Распределение температуры в областях описывается уравнениями:
10
Здесь Ti(t, z), Ts(t, z) - температура жидкой и твердой фаз соответственно; f(t) - неизвестная граница раздела фаз; af, a2s - соответствующие жидкой и твердой фазам коэффициенты температуропроводности; t -переменная по времени; z - переменная по пространству; Н - размер заданной области.
Для замыкания системы задаются граничные условия. На подвижной границе z = f(t) выполняются условия сопряжения:
1) условие Стефана, которое описывает тепловой баланс (скачком плотности при фазовом переходе пренебрегаем):
XpsVf = h—
-kf
z=f
i
dz
z=h
2) равенство температуры среды температуре фазового перехода Т/ данного вещества, которая считается известной постоянной величиной. Для воды Tf = T* = 0°C.
Здесь Vf(t) - скорость нарастания льда; Т* - температура замерзания (кристаллизации) чистого вещества; hi, ks - коэффициенты теплопроводности, соответственно для жидкой и твердой среды; Л - скрытая теплота плавления; pi, ps- соответствующая средам плотность; ls - толщина твердой фазы.
На внешних границах, как правило, считается известной либо температура, либо тепловой поток. На границе z = Н должна поддерживаться температура, при которой происходит кристаллизация. Для воды Ts(t,H)<0°C.
Кроме граничных условий и условий сопряжения формулируются начальные условия. Начальное распределение температуры по толщине можно определить функцией или константой. Для задания положения подвижной границы при малых временных параметрах t ~ At, когда скорость продвижения границы теоретически равна бесконечности, зна-
11
чение толщины зарождающейся твердой фазы определяется из начальной асимптотики скорости движения граница раздела, предложенной А.Н. Тихоновым и А.А. Самарским (1972).
Метод численного решения основан на отображении областей с криволинейными границами в регулярные - метод спрямления фронта (Бу-дак Б.М. и др. 1966). При аппроксимации уравнений полученной системы используется неявная схема с направленными разностями для конвективных слагаемых. Алгебраическая система уравнений решается методом прогонки. Путём вычислительных экспериментов на последовательности сеток показана устойчивость предложенного алгоритма.
Пункт 1.2. посвящен описанию математической постановки и методу численного решения термодиффузионной задачи Стефана. В такой постановке задача давно успешно используется при моделировании процессов фазовых переходов в бинарных смесях (производство полупроводниковых материалов, очистка методом направленной кристаллизации, металлургическое производство), подробно описанных Б.Я. Любо-вым (1969, 1975), Н.А. Авдониным (1980), A.M. Мейрмановым (1986).
Здесь рассматривается двухфазная термодиффузионная задача Стефана о замерзании раствора. При образовании пресного льда перед фронтом кристаллизации образуется слой с повышенным содержанием примеси, что влияет на температуру замерзания.
Область ?li(t) — {0 < z < fi{t)} содержит раствор соли (солей). Распределение температуры и концентрации примеси описывается уравнениями:
основные обозначения сохранены, C(t,z) - концентрация примеси, d - ко-
12
эффициент диффузии соли в воде. В образующейся твердой фазе ?ls(t) — {h(t) < z < h{t)} предполагаем отсутствие диффузии примеси, а распределение температуры аналогичное:
а (9)
В области, содержащей две фазы, имеются две подвижные границы, положение которых описывается следующими формулами:
fi(t) = li = Н — kpls - граница между жидкой и твердой фазой, /2CO — hit) + h ~ внешняя граница, мало перемещаемая за счет разности плотностей твердой и жидкой фаз, то есть плавучести льда (Од-рова Т.В., 1979). Здесь li, ls - толщина слоя жидкой и твердой фаз, соответственно; кр = &•.
На внешней границе z = /г(^) происходит охлаждение раствора. На
подвижной границе z — fi(t) выполняются условия сопряжения:
1) условие Стефана - уравнение (6);
2) равенство температур жидкой и твердой фаз неизвестной температуре замерзания Т/, определяемой уравнением:
Tf=T*- 7С/; (10)
3) баланс массы растворенного в воде вещества (Васильев В.И. и др., 1997):
= ~d
дС
dt dz
(И)
На границе z = 0 для температуры и примеси определены: условие изоляции (отсутствие притока тепла и примеси), или закон теплообмена со средой, или постоянное значение.
Tf(t) - в данном случае неизвестная температура фазового перехода; Cf(t) - искомое значение примеси на границе раздела фаз; 7 - равновесный коэффициент распределения примеси.
13
Кроме граничных условий и условий сопряжения формулируются начальные условия: распределение температуры по толщине в двух фазах, значение концентрации примеси (оно, как правило постоянно С = Со).
Начальные значения V/ и ls определялись путем численного эксперимента и согласованы с асимптотикой для однородной жидкости (п. 1.1).
Для численного решения задачи использовались метод спрямления фронта. Аппроксимация уравнений проводилась как и в п. 1.1. В каждой из областей была построена равномерная сетка. Так как температура замерзания Т/ и солености С/ на границе z = fi(t) заранее неизвестны, то здесь классический метод прогонки использовать невозможно. Кроме того, в связи с тем, что начальная стадия процесса приводит к возникновению больших градиентов температуры и солености, то применение методов ловли фронта в узел сетки или сквозного счета вызывает большие трудности. Поэтому для численной реализации задачи была разработана модификация метода (Воеводин А.Ф., Леонтьев Н.А., Петрова А.Г., 1982), основанного на методе встречной прогонки, позволяющего эффективно использовать безитерационные (точные) методы решения разностной задачи, несмотря на ее нелинейность.
Идея этой модификации состоит в том, что значения Tj, Ts, С на границе z = f(t)i записываются через прогоночные коэффициенты. После подстановки полученных выражений в уравнения (6), (10), (11) имеем квадратное уравнение относительно значения концентрации примеси на границе раздела фаз
jATC2f + (Ас - &r)Cf + Вс = 0. (12)
Результатом решения этого уравнения является один удовлетворяющий физическим условиям корень:
14
Вт - Ac + J(Ac — Вт)2 - 4тАтBe
где выражения Ат, Вт, Ас, Вс находятся однозначно через прогоночные коэффициенты.
Получив значение примеси на границе раздела z = fi(t), из уравнения (10) находим температуру фазового перехода, соответствующую температуре жидкой и твердой фаз на этой границе. Зная граничные условия, восстанавливаем значения температуры и примеси на новом шаге по времени.
Вторая глава содержит численные модели, описывающие динамику роста ледового покрова водоемов, построенные на основе классической и термодиффузионной задачах Стефана.
Для моделирования ледового покрова естественных водоемов недостаточно описать процесс фазового перехода воды в лед. Необходимо точно учитывать наиболее важные физические факторы, оказывающие влияние на рост льда. Такими факторами являются, например, снежный покров и минерализировашюсть водоема.
Снег за счет своих теплоизоляционных свойств оказывает особое влияние на рост льда, уменьшая влияние отрицательных атмосферных температур. Снежный покров по высоте имеет неоднородную структуру. Соответственно, изменяются его физические свойства, такие как плотность и теплопроводность. Плотность снега меняется со временем в зависимости от усадки, ветрового воздействия, плотности свежевыпавше-го снега. Это приводит к нелинейному распределению температуры в слое снега, что и было отмечено при исследовании водного и ледового режимов на Новосибирском водохранилище в 1982 г. экспедицией
15
Н И С И им.В.В. Куйбышева (Козлов Д.В., 2000). Следовательно, появилась необходимость модифицировать математическую модель и ввести новую область, описывающую процессы в снежном покрове.
В моделях используется описанная в первой главе базовая система уравнений двухфазной задачи Стефана. С появлением третьего слоя h{t) < z < /з(0> содержащего снег, образуется еще одна подвижная неизвестная граница z = /з(?).
ОJ-water __ т & ¦'¦water
PCp 777 — ^
&-lice 2 ® ¦'¦ice
п
ВТ f) I f)T
UJ-snow __ (и UJ-snow
CpPsnow 7ГТ ~" I fc 7^
ot oz \ oz )
Mt) = /2(t) + hnew, hnaw(t) = i ln(l + &CnoJt)),
Tsrww — температура в слое воды, льда, снега. psnow = poeb^3^~z^ - формула Абе, b = 1.255,
= 2.9 • 10~6/9gnou; + 0.043 - зависимость, предложенная В.В. Пиатро-
вичем (1968) для нахождения теплопроводности снега, psnow - плотность снега, ро - плотность свежевыпавшего снега, lsnaw - толщина слоя снега, l*now - толщина слоя свежевыпавшего снега.
На дне водоема при z = 0 принимаем для температуры воды и примеси постоянные значения или условия изоляции. На границе z = fi(t) - условие равенства температур и условия сопряжения (6), (10), (11). На границе z = /г(^) задается условие равенства тепловых потоков со стороны слоя льда и слоя снега и равенство температур. Значение температуры между льдом и снегом рассчитывалось в ходе решения. Темпе-
16
ратура снега на внешней границе z = /з(?) задавалась равной атмосферной температуре, измеренной на расстоянии 2 м над землей - условие, удовлетворяющее зимнему метеорологическому режиму Сибирского региона (Пиотрович В.В., 1963).
В связи с тем, что рассматриваемые водоемы являются слабопроточными, то конвективные перенос субстанций не учитывается.
Метод решения изложенной задачи аналогичен методу, описанному в первой главе.
Для пресного водоема задача упрощается. При С = 0 температура фазового перехода Т* = 0. Уравнения (10), (11) и (15) не рассматриваются.
Также для пресного водоема предложено решение проблемы начальной асимптотики для скорости продвижения фронта, за счет введения новой переменной S(t) = lfce(t) (Атавин А.А., Гранкина Т.Б., 2004). После такой замены начальное значение 1{се можно положить равным 0.
Третья глава содержит описание общепринятых методик расчета толщины ледового покрова водоемов, результаты расчетов по предложенным математическим моделям для водоемов Западной Сибири: Новосибирского водохранилища - пресный водоем, озера Чаны - система озер и плесов различной минерализации. Также проведен сравнительный анализ результатов расчетов по инженерным методикам и математическим моделям и исследование области применимости приближенных инженерных методов для решения задачи динамики ледового покрова.
В заключении описаны полученные результаты:
• разработана математическая модель для расчета процесса формирования ледового покрова водоемов различной минерализации с учетом таких важных физических факторов, как наличие снежного
17 |
