Содержание к диссертации
Введение
Глава 1 Задача стефана для моделирования процессов с фазовым переходом 22
1.1 Классическая одномерная двухфазная задача Стефана . 22
1.1.1 Математическая формулировка задачи 23
1.1.2 Граничные и начальные условия 24
1.1.3 Начальная асимптотика 26
1.1.4 Метод численного решения 26
1.1.5 Аналитическое решение 30
1.2 Термодиффузионная задача Стефана 31
1.2.1 Математическая формулировка термодиффузионной задачи Стефана 31
1.2.2 Граничные и начальные условия 33
1.2.3 Метод численного решения 34
1.2.4 Гезультаты тестовых расчетов 38
Глава 2 Численное моделирование динамики роста снежно-ледов ого покрова водоемов 43
2.1 Моделирование динамики роста ледового покрова пресного водоема 45
2.1.1 Математическая постановка задачи 45
2.1.2 Начальные и граничные условия 46
2.1.3 Начальная асимптотика 47
2.1.4 Метод численного решения 48
2.2 Моделирование динамики роста ледового покрова минера лизированного водоема 50
2.2.1 Постановка задачи 50
2.2.2 Граничные и начальные условия 52
2.2.3 Метод численного решения 53
Глава 3 Результаты расчетов процесса льдообразования на водоемах Западной Сибири 55
3.1 Методика О. Дэвика 56
3.2 Расчет толщины ледяного покрова по методике В.В. Пио-тровича
3.3 Ледотермический режим Новосибирского водохранилища. Результаты расчетов 59
3.4 Ледотермический режим озер Чановской системы. Результаты расчетов 64
Заключение 70
Литература 71
- Математическая формулировка задачи
- Математическая формулировка термодиффузионной задачи Стефана
- Моделирование динамики роста ледового покрова минера лизированного водоема
- Расчет толщины ледяного покрова по методике В.В. Пио-тровича
Введение к работе
Актуальность темы. В связи с проблемой наблюдающегося глобального изменения климата все более актуальными становятся вопросы, связанные с изучением последствий воздействия изменения метеорологических условий на гидрологические и экологические процессы в водоемах и водотоках при их круглогодичной эксплуатации. Большое значение приобретает рациональное использование водных ресурсов в зимнее время, когда очень важным становится обеспечение достоверности прогноза процесса формирования ледяного покрова и его прочности, позволяющее определить начало и сроки эксплуатации ледовых трасс, проложенных по поверхности водных объектов.
Разработкой методов расчета толщины ледяного покрова занимались многие исследователи: Я.Л. Готлиб (1983), Л.Г. Шуляковский (I960,1969, 1972), В.В. Пиотрович (1968, 1969, 1970), Р.В. Донченко (1971, 1983), Т.В. Одрова (1979), Г. Эштон (1978,1986), А.И. Пехович (1983), Г.А. Трегуб (1994, 1997), И.Н. Шаталина (1990), В.М. Мишен (1997) и многие другие. К настоящему времени предложено большое количество формул и расчетных приемов для определения толщины ледяного покрова пресноводных водоемов и водотоков, опирающихся на основные методы расчета: эмпирический и теоретический методы.
Эмпирический метод основан на отыскании эмпирических связей толщины льда и отдельных факторов, определяющих изменение толщины ледяного покрова. В этом случае расчетные эмпирические соотношения получены на основании относительно тесной корреляции между некоторыми температурными характеристиками и толщиной льда. Как прави-
ло, эти характеристики привязаны к конкретному региону.
Теоретический метод основан на интегрировании исходных дифференциальных уравнений, описывающих физическую сущность нарастания толщины льда.
Основой для получения эмпирических формул послужило уравнение для подсчета теплового баланса, так называемое условие Стефана. В 1889 г. И. Стефан в работе о промерзании грунта поставил и решил следующую задачу (Рубинштейн Л.И., 1967).
Среда, находящаяся в двух фазовых состояниях (жидком и твердом) и проводящая теплоту исключительно посредством теплопроводности, заполняет полупространство х > 0. В начальный момент времени вся она находится при постоянной температуре Ті > 0С. На плоскости х = 0 поддерживается постоянная температура Ті < 0С, под воздействием которой происходит кристаллизация, протекающая изотермически при температуре Т = 0С без переохлаждения и с пренебрежимо малым объемным эффектом.
Требуется определить температуры Ui(x,t) и Ui(x,t) твердой и жидкой фазы и положение х = y(t) границы раздела фаз.
Подсчет теплового баланса приводит, как показал Стефан, к условию
Л ^ = (kl± - к2Щ . (1)
dt К1 дх 2 дх J x=y(t)
Здесь Л - скрытая теплота кристаллизации, отнесенная к единице массы; р - плотность образующейся фазы; к\ и hi - коэффициенты теплопроводности соответственно твердой и жидкой фазы.
Эмпирические формулы для определения толщины ледового покрова, полученные на основе решения задачи Стефана, имеют общий вид :
л.-се = <ч/:с(-02), (2)
«
где hice - толщина льда, см; (—#г) - сумма средних суточных температур воздуха (С) на высоте 2м над уровнем водоема за расчетный период времени; а - эмпирический коэффициент, определяемый по данным непосредственных наблюдений и отражающий в среднем те условия, которые имели место в период наблюдений (температуру воды, высоту и плотность снежного покрова, глубину водоема, скорость подледного течения воды). Однако ввиду различия этих факторов даже для отдельных участков водоемов и недостаточной продолжительности наблюдений указанный параметр существенно меняется. Поэтому все подобные формулы носят локальный характер. Подробный обзор работ, описывающих получение эмпирических формул различными авторами для конкретных водных объектов представлен в монографии Д.В. Козлова (2000).
Очевидно, что в подобных эмпирических формулах не находят в полной мере своего отражения такие важные факторы, как мощность и физические свойства снежного покрова, интенсивность теплопотока из водной массы, метеорологические условия. Тем не менее, снег за счет своих теплоизоляционных свойств оказывает особое влияние на рост льда, уменьшая влияние отрицательных атмосферных температур.
В своих исследованиях норвежский ученый О. Дэвик (Винников С.Д., Проскуряков Б.В., 1988) к процессу льдообразования подошел с позиций физики, также основываясь на уравнении Стефана. В предложенной им формуле учитывается теплообмен на верхней поверхности снежного покрова, теплофизические свойства снега. Его подход стал классическим и получил широкое развитие. В дальнейших разработках многие ученые уделяли больше внимания процессу теплообмена снежно-ледовой поверхности с атмосферой (В.В. Пиотрович), теплообмену с дном водоема, водным потокам, начальной стадии формирования льда, образованию шуги
(Л.Г. Шуляковский, Р.В. Донченко, С.Д. Винников, Б.В. Проскурякрв, Г.А. Трегуб, И.Н. Шаталина).
На основе обощения этих разработок сформированы методические рекомендации, выпущенные научноисследовательскими и проектными институтами ВНИИГ, ГГИ, Ленгидропроект.
С развитием вычислительных методов математического моделирования и компьютерных технологий последнее время стали чаще обращаться к теоретическим методам, более точно описывающим физические явления процесса льдообразования. Появилась возможность численно решать сложные системы дифференциальных уравнений (Цибульский В.Р. и др., 1993: Васильев О.Ф., Бочаров О.Б., Зиновьев А.Т., 1999: Васильев В.И. и др., 1997: Дебольская Б.И., 2003: Прокофьев В.А., 2004), описывающих процессы льдообразования. Не остались без внимания важные задачи прогнозирования ледотермических режимов водоемов-охладителей тепловых и атомных электростанций, влияние их эксплуатации на температурный и ледовый режим рек (Квон В.И., Филатова Т.Н.). Над задачами моделирования ледового режима бьефов Красноярской ГЭС успешно трудились В.М. Белолипецкий, С.Н. Генова, В.Б. Туговиков, Ю.И. Шо-кин (1991, 1993).
На юге Западной Сибири имеется значительное число минерализированных озер. Как показали проведенные исследования, даже слабая минерализация оказывает существенное влияние на процесс замерзания. А именно, при формировании ледового покрова в минерализированном водоеме в результате образования практически пресного льда в воде перед фронтом кристаллизации формируется слой с повышенным содержанием химических веществ, что существенно влияет на температуру фазового перехода, а тем самым и на весь процесс формирования ледо-
вого покрова (Власов В.П., 1968, 1976), что необходимо учитывать при разработке математической модели и выполнении соответствующих расчетов.
Процесс замерзания морских вод достаточно подробно изучен, описан в учебниках по океанологии и статьях (Смирнов Г.Н., 1974: Федоров К.Н., 1976: Саркисян А.С, Демин Ю.Л., Бреховских А.Л., Шаха-нова Т.В., 1986: Багно А.В., Гаращук Р.В., Залесный В.Б., 1996: Яковлев Н.Г., 2003). Гидродинамические и гидрохимические режимы морских вод отличается от режимов соленых озер, которые по солевому составу чаще являются солоноватыми. Поэтому способы описания процесса образования морских льдов не применимы к озерным.
Во всех этих работах, основанных как на эмпирических, так и на теоретических методах, снежный покров в процессе льдообразования принято учитывать через эквивалентный слой:
к-
i^ecv == Щсе ~г i^snowT 5 V"J
Ksnow
где hecv - эквивалентный слой снежно-ледового покрова, hmmv - высота снежного покрова на льду, kice и 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)} и Qs(i) = {f(t) < z < Н}, занятые соответственно жидкой и твердой фазами данного вещества с гладкой границей раздела z = /(). Распределение температуры в областях описывается уравнениями:
от, 2д2т,
Ж = а'9Ї5- (4)
-ж=а« -* (5)
Здесь Ti(t, z), Ts{t, z) - температура жидкой и твердой фаз соответственно; f(t) - неизвестная граница раздела фаз; of, a2s - соответствующие жидкой и твердой фазам коэффициенты температуропроводности; t -переменная по времени; z - переменная по пространству; Н - размер заданной области.
Для замыкания системы задаются граничные условия. На подвижной границе z = f(t) выполняются условия сопряжения:
1) условие Стефана, которое описывает тепловой баланс (скачком плотности при фазовом переходе пренебрегаем):
XpsVf = ki
-к
>*-% ^
Z=fi
z=h
2) равенство температуры среды температуре фазового перехода Т/ данного вещества, которая считается известной постоянной величиной. Для воды Tf = T* = 0С.
Здесь Vf(t) - скорость нарастания льда; Т* - температура замерзания (кристаллизации) чистого вещества; ki, ks - коэффициенты теплопроводности, соответственно для жидкой и твердой среды; Л - скрытая теплота плавления; pi, ps- соответствующая средам плотность; ls - толщина твердой фазы.
На внешних границах, как правило, считается известной либо температура, либо тепловой поток. На границе z = Н должна поддерживаться температура, при которой происходит кристаллизация. Для воды Ts(t,H)<0C.
Кроме граничных условий и условий сопряжения формулируются начальные условия. Начальное распределение температуры по толщине можно определить функцией или константой. Для задания положения подвижной границы при малых временных параметрах t ~ At, когда скорость продвижения границы теоретически равна бесконечности, зна-
чение толщины зарождающейся твердой фазы определяется из начальной асимптотики скорости движения граница раздела, предложенной А.Н. Тихоновым и А.А. Самарским (1972).
Метод численного решения основан на отображении областей с криволинейными границами в регулярные - метод спрямления фронта (Бу-дак Б.М. и др. 1966). При аппроксимации уравнений полученной системы используется неявная схема с направленными разностями для конвективных слагаемых. Алгебраическая система уравнений решается методом прогонки. Путём вычислительных экспериментов на последовательности сеток показана устойчивость предложенного алгоритма.
Пункт 1.2. посвящен описанию математической постановки и методу численного решения термодиффузионной задачи Стефана. В такой постановке задача давно успешно используется при моделировании процессов фазовых переходов в бинарных смесях (производство полупроводниковых материалов, очистка методом направленной кристаллизации, металлургическое производство), подробно описанных Б.Я. Любо-вым (1969, 1975), Н.А. Авдониным (1980), A.M. Мейрмановым (1986).
Здесь рассматривается двухфазная термодиффузионная задача Стефана о замерзании раствора. При образовании пресного льда перед фронтом кристаллизации образуется слой с повышенным содержанием примеси, что влияет на температуру замерзания.
Область Qi(t) — {0 < z < fi(t)} содержит раствор соли (солей). Распределение температуры и концентрации примеси описывается уравнениями:
основные обозначения сохранены, C(t,z) - концентрация примеси, d - ко-
эффициент диффузии соли в воде. В образующейся твердой фазе Qs(t) — {fi(t) < z < /2(^)} предполагаем отсутствие диффузии примеси, а распределение температуры аналогичное:
В области, содержащей две фазы, имеются две подвижные границы, положение которых описывается следующими формулами:
fi(t) = її = Н — kpls - граница между жидкой и твердой фазой, /2 СО — /l (t) + h ~ внешняя граница, мало перемещаемая за счет разности плотностей твердой и жидкой фаз, то есть плавучести льда (Од-рова Т.В., 1979). Здесь li, ls - толщина слоя жидкой и твердой фаз, соответственно; кр = &.
На внешней границе z — /2(^) происходит охлаждение раствора. На подвижной границе z — fi(t) выполняются условия сопряжения:
условие Стефана - уравнение (6);
равенство температур жидкой и твердой фаз неизвестной температуре замерзания Т/, определяемой уравнением:
Tf=T*- 7С/; (10)
3) баланс массы растворенного в воде вещества (Васильев В.И. и др.,
1997):
CfTt=-dTzz=h (11)
На границе z = 0 для температуры и примеси определены: условие изоляции (отсутствие притока тепла и примеси), или закон теплообмена со средой, или постоянное значение.
Tf(t) - в данном случае неизвестная температура фазового перехода; Сj{t) - искомое значение примеси на границе раздела фаз; у - равновесный коэффициент распределения примеси.
Кроме граничных условий и условий сопряжения формулируются начальные условия: распределение температуры по толщине в двух фазах, значение концентрации примеси (оно, как правило постоянно С = Со).
Начальные значения V/ и ls определялись путем численного эксперимента и согласованы с асимптотикой для однородной жидкости (п. 1.1).
Для численного решения задачи использовались метод спрямления фронта. Аппроксимация уравнений проводилась как и в п. 1.1. В каждой из областей была построена равномерная сетка. Так как температура замерзания Т/ и солености С/ на границе z = fi(t) заранее неизвестны, то здесь классический метод прогонки использовать невозможно. Кроме того, в связи с тем, что начальная стадия процесса приводит к возникновению больших градиентов температуры и солености, то применение методов ловли фронта в узел сетки или сквозного счета вызывает большие трудности. Поэтому для численной реализации задачи была разработана модификация метода (Воеводин А.Ф., Леонтьев Н.А., Петрова А.Г., 1982), основанного на методе встречной прогонки, позволяющего эффективно использовать безитерационные (точные) методы решения разностной задачи, несмотря на ее нелинейность.
Идея этой модификации состоит в том, что значения Ті, Ts, С на границе z = f(t)i записываются через прогоночные коэффициенты. После подстановки полученных выражений в уравнения (6), (10), (11) имеем квадратное уравнение относительно значения концентрации примеси на границе раздела фаз Cf.
1АТС) + (Ас - BT)CS + Вс = 0. (12)
Результатом решения этого уравнения является один удовлетворяющий физическим условиям корень:
c = Вт -Ac + у/(Ac - Вт)2 - 4jATBc
где выражения At, Вт, Ac, Be находятся однозначно через прогоночные коэффициенты.
Получив значение примеси на границе раздела z = fi(t), из уравнения (10) находим температуру фазового перехода, соответствующую температуре жидкой и твердой фаз на этой границе. Зная граничные условия, восстанавливаем значения температуры и примеси на новом шаге по времени.
Вторая глава содержит численные модели, описывающие динамику роста ледового покрова водоемов, построенные на основе классической и термодиффузионной задачах Стефана.
Для моделирования ледового покрова естественных водоемов недостаточно описать процесс фазового перехода воды в лед. Необходимо точно учитывать наиболее важные физические факторы, оказывающие влияние на рост льда. Такими факторами являются, например, снежный покров и минерализировашюсть водоема.
Снег за счет своих теплоизоляционных свойств оказывает особое влияние на рост льда, уменьшая влияние отрицательных атмосферных температур. Снежный покров по высоте имеет неоднородную структуру. Соответственно, изменяются его физические свойства, такие как плотность и теплопроводность. Плотность снега меняется со временем в зависимости от усадки, ветрового воздействия, плотности свежевыпавше-го снега. Это приводит к нелинейному распределению температуры в слое снега, что и было отмечено при исследовании водного и ледового режимов на Новосибирском водохранилище в 1982 г. экспедицией
Н И С И им.В.В. Куйбышева (Козлов Д.В., 2000). Следовательно, появилась необходимость модифицировать математическую модель и ввести новую область, описывающую процессы в снежном покрове.
В моделях используется описанная в первой главе базовая система уравнений двухфазной задачи Стефана. С появлением третьего слоя h{t) < z < /з(0> содержащего снег, образуется еще одна подвижная неизвестная граница z = /з().
О* water т & * water (лл\
Pwatercp ТГ, — &water ~7Г~2 ) \^)
дс ^С п *
-Ж = dM' (15)
О lice 2 & J-ісе ґт\
~дГ = а-^Г' (16)
O-Lsnow (и 0-Lsnow\ f17^
CpPsnow q! — -qZ I Ksnow ~7T I 5 Vі 4
Ш) = W) + hnow, hnaw{t) = ln(l + bl*snow{t)),
Twater, Tice, Tsnaw - температура в слое воды, льда, снега.
psnow = poeb^3^~z^ - формула Абе, b = 1.255,
ksnow = 2.9 10_6/9gnou; + 0.043 - зависимость, предложенная В.В. Пиатро-
вичем (1968) для нахождения теплопроводности снега, psn0w ~ плотность
снега, pQ - плотность свежевыпавшего снега, lSnaw ~ толщина слоя снега,
l*snaw - толщина слоя свежевыпавшего снега.
На дне водоема при z = 0 принимаем для температуры воды и примеси постоянные значения или условия изоляции. На границе z = f\(t) - условие равенства температур и условия сопряжения (6), (10), (11). На границе z = /2(^) задается условие равенства тепловых потоков со стороны слоя льда и слоя снега и равенство температур. Значение температуры между льдом и снегом рассчитывалось в ходе решения. Темпе-
ратура снега на внешней границе z = fz(t) задавалась равной атмосферной температуре, измеренной на расстоянии 2 м над землей - условие, удовлетворяющее зимнему метеорологическому режиму Сибирского региона (Пиотрович В.В., 1963).
В связи с тем, что рассматриваемые водоемы являются слабопроточными, то конвективные перенос субстанций не учитывается.
Метод решения изложенной задачи аналогичен методу, описанному в первой главе.
Для пресного водоема задача упрощается. При С = 0 температура фазового перехода Т* = 0. Уравнения (10), (11) и (15) не рассматриваются.
Также для пресного водоема предложено решение проблемы начальной асимптотики для скорости продвижения фронта, за счет введения новой переменной S(t) = lfce(t) (Атавин А.А., Гранкина Т.Б., 2004). После такой замены начальное значение /г-се можно положить равным 0.
Третья глава содержит описание общепринятых методик расчета толщины ледового покрова водоемов, результаты расчетов по предложенным математическим моделям для водоемов Западной Сибири: Новосибирского водохранилища - пресный водоем, озера Чаны - система озер и плесов различной минерализации. Также проведен сравнительный анализ результатов расчетов по инженерным методикам и математическим моделям и исследование области применимости приближенных инженерных методов для решения задачи динамики ледового покрова.
В заключении описаны полученные результаты:
разработана математическая модель для расчета процесса формирования ледового покрова водоемов различной минерализации с учетом таких важных физических факторов, как наличие снежного
покрова и зависимость температуры фазового перехода "вода-лед" от концентрации солей в воде;
разработан и реализован на ЭВМ численный метод решения многофазной термодиффузионной задачи Стефана, с температурой фазового перехода, зависящей от концентрации примеси;
разработаны и реализованы на ЭВМ методы расчета процесса формирования ледового покрова водоемов различной минерализации;
создана численная модель и проведены методические расчеты для исследования ледотермических процессов на реальных объектах Западной Сибири (Новосибирское водохранилище, озера Чановской системы), проведен сравнительный анализ результатов расчетов с приближенными методами и сравнение с натурными натурными данными, что показало хорошую точность и эффективность расчетов по созданной модели.
Научная новизна. Разработан численный метод для решения задач типа Стефана. Создана и апробирована математическая модель, а также комплекс программ для исследования процесса динамики ледяного покрова водоемов с различной минерализацией. Впервые в математической постановке задачи учтены пространственное изменение плотности и температуры снежного покрова, а также зависимость температуры фазового перехода от минерализации водной среды. Разработаны численные модели, позволяющие исследовать ледотермические процессы реальных объектов Западной Сибири.
Практическая значимость. Полученные результаты могут быть использованы для прогноза динамики роста ледового покрова на водоемах Западной Сибири в условиях изменения климата, уровня воды в водо-
емах и степени минерализации водной среды. Предлагаемые математические модели, методы численного решения и программное обеспечение могут быть использованы при проведении экологического мониторинга водных объектов и организации ледовых переправ. Также, разработанные модели и метод могут быть использованы при исследовании физических, биологических и технологических процессов, приводящих к необходимости численного решения задач фазового перехода вещества, содержащего примесь.
Объект и методы исследований. Объектом исследований являются процессы формирования ледового покрова Новосибирского водохранилища и озер Чановской системы. Выполнение исследований базируется на математическом моделировании и численных методах решения термодиффузионной задачи Стефана.
На защиту выносится:
математическая модель для расчета процесса формирования ледового покрова водоемов различной минерализации с учетом физических факторов, таких как наличие снежного покрова и зависимость температуры фазового перехода "вода-лед" от концентрации солей в воде;
численный метод решения многофазной термодиффузионной задачи Стефана;
методы расчета процесса формирования ледового покрова водоемов различной минерализации;
результаты расчетов, выполненные для исследования ледотермиче-ских процессов на ряде реальных объектах Западной Сибири.
Апробация работы. Результаты докладывались на 1-ой межвузовской научно - практической конференции студентов и молодых ученых Сахалинской области (Южно-Сахалинск, 1997 г.), XXXII научно - методической конференции преподавателей ЮСГПИ (Южно-Сахалинск, 1998 г.), IV школе-семинаре "Математические проблемы механики сплошных сред"ИГиЛ (Новосибирск, 2000), XXXIX Международной научной студенческой конференции "Студент и научно - технический прогресс"(Новосибирск, 2001 г.), Международной конференции "Вычислительные технологии и математическое моделирование в науке, технике и образовании "(Алматы, 2002 г.), Всероссийской конференции "Математические методы в механике природных сред и экологии"(Барнаул, 2002 г.), 61-й научно - практической конференции профессорско - преподавательского состава НГАСУ (Новосибирск, 2004 г.), 17-м Международном симпозиуме по льду (Санкт-Петербург, 2004 г.), Международной конференции по измерениям, моделированию и информационным системам для изучения окружающей среды "ENVIROMIS-2004"(Томск, 2004 г.), VI конференции "Динамика и термика рек, водохранилищ и прибрежной зоны морей"(Москва, 2004 г.), 3-й и 4-й конференциях молодых ученых ИВЭП СО РАН (Барнаул, 2004, 2005 гг.), научной конференции "Фундаментальные проблемы изучения и использования воды и водных ресурсов" (Иркутск, 2005 г.), научных семинарах лаборатории прикладной и вычислительной гидродинамики Института гидродинамики им. М.А. Лаврентьева и Новосибирского филиала ИВЭП СО РАН.
Автор выражает глубокую благодарность научному руководителю проф. д.ф.-м.н. А.Ф. Воеводину за постановку задачи, постоянно внимание к работе, к.т.н. А.А. Атавину и д.ф.-.м.н. В.И. Квону за плодотворное сотрудничество, ценные советы и организационное содействие. Также
выражает признание сотрудникам Новосибирского филиала ИВЭП СО РАН и Лаборатории прикладной гидродинамики ИГиЛ СО РАН за помощь в работе, консультации, своевременную критику и дружескую поддержку.
Математическая формулировка задачи
Классическая постановка задачи, используемая при математическом описании многообразных процессов кристаллизации, была предложена Й. Стефаном еще в 1889 г. Исторически она формулируется следующим образом.
Область Q, занятая изучаемым веществом, в момент времени t О разбивается некоторой гладкой поверхностью z — f(t), на две подобласти Qt(t) = {О z /( )} и Qs(t) = {f(t) z #}, занятые соответственно жидкой и твердой фазами данного вещества. Поверхность z = f(t) неизвестна и подлежит определению. Характеристики жидкой и твердой фазы строго не определялись.
Здесь Vf(t) - скорость перемещения искомой границы в направлении нормали к ней; Л - скрытая удельная теплота плавления (кристаллизации), равная минимальному количеству энергии, которое необходимо сообщить единице массы рассматриваемого вещества для того, чтобы перевести его из твердого состояния в жидкое.
Известны различные численные методы решения задачи Стефана, такие как методы сквозного счета, используемый в работе (Бондарев Э.А., Васильев В. И., 1984: Васильев В.И. и др., 1997), методы с явным выделением фронта (Белолипецкий и др., 1993: Овчарова А.С, 1994,1995,1999). В работе (Бондарев Э.А., Попов Ф.С., 1989) приводятся оценки точности наиболее употребляемых приближенных методов решения задачи Стефана. Воспользуемся методом "спрямления фронта"(Будак Б.М. и др., 1966). Проведем преобразования, позволяющему свести две области с подвижными границами в фиксированные. N - количество узлов сетки,h - шаг по пространству, т - шаг по времени, п - номер шага по времени, a2 = [of,а2) - соответствующий для уравнений (1.1.6)-(1.1.7) коэффициент температуропроводности, и = (Ті, Ts) - соответствующая для уравнений (1.1.6)-(1.1.7) функция, I = (k, ls) - соответствующая длина жидкой и твердой области, v = (vi, vs) - соответствующая для уравнений (1.1.6)-(1.1.7) скорость. В каждой из областей строилась равномерная сетка. Шаг по времени брался постоянным. Для численной реализации задачи использовался метод прогонки. Решения для Ті и Ts будем искать в виде: иГ1 = «ІМІ1 + /Йь і = 1, -N, (1.1.19) где прогоночные коэффициенты а и / считаем по рекуррентным формулам: afo = -, / = Аф +Д і = 2..N. (1.1.20) Согласно граничным условиям: (1.1.3) al2 = 0 и / = Ті, для условия (1.1.4) а12 = 1 и / = 0. Остальные значения: и +1 = Т , а2 = 0, Pi = T ,unN+1 = T2. Найдя значения прогоночных коэффициентов а 3 и Д в по системе уравнений (1.1.18) и выражениям (1.1.20) находим значение u"+1 по формуле (1.1.19), т.е. соответствующее решению системы основных уравнений (1.1.1) и (1.1.2) для Ті и Т3 на новом временном шаге. Положение свободной границы /() получаем из (1.1.16) и (1.1.6). Таким образом получаем численное решение двухфазной задачи Стефана.
Рост льда является теплоэнергетическим процессом. Поэтому для выведения расчетных формул, оценивающих толщину ледяного покрова водоемов, за основу стали брать уравнение теплового баланса (условие Стефана). л dl идТ1 идТ z=h dt dz В водоемах с малыми скоростями течения теплоприток от воды к нижней поверхности льда достаточно мал. Поэтому первым слагаемым в правой части условия Стефана можно пренебречь. При установившемся режиме замерзания распределение температуры в слое льда можно принять линейным. С этими допущениями мы получаем следующее упрощенное уравнение: is = a(t)k dt Xps проинтегрировав которое получаем известную квадратичную форму: гее \ Ргсе о 2к Нсе\Р) — / (a(r))dr. (1.1.21) Уравнение (1.1.21) согласуется с формулами, предложенными в рекомендациях и руководствах для термических расчетов и гидрологических прогнозов водохранилищ. Проведенное сравнение тестовых расчеты по описанному выше методу решения двухфазной задачи Стефана для теплофизических свойств сред - воды и льда - с точным решением упрощенной задачи (1.1.21) показало качественное совпадение результатов.
Термодиффузионная задача Стефана давно успешно используется при моделировании процессов фазовых переходов в бинарных смесях (производство полупроводниковых материалов, очистка методом направленной кристаллизации, металлургическое производство), подробно описанных Б.Я. Любовым (1969, 1975), Н.А. Авдониным (1980), A.M. Мейрмановым (1986).
Здесь рассматривается двухфазная термодиффузионная задача Стефана о замерзании раствора. При образовании пресного льда перед фронтом кристаллизации образуется слой с повышенным содержанием концентрации примеси, что влияет на температуру замерзания.
Благодаря уменьшению плотности воды при замерзании лед плавает на воде, незначительно возвышаясь над ее поверхностью (Одрова Т.В., 1979). Приняв форму льдины в виде прямоугольного параллелепипеда с основанием, равным единице, согласно закону Архимеда, подсчитаем соотношение ее надводной и подводной частей.
Математическая формулировка термодиффузионной задачи Стефана
Термодиффузионная задача Стефана давно успешно используется при моделировании процессов фазовых переходов в бинарных смесях (производство полупроводниковых материалов, очистка методом направленной кристаллизации, металлургическое производство), подробно описанных Б.Я. Любовым (1969, 1975), Н.А. Авдониным (1980), A.M. Мейрмановым (1986). Здесь рассматривается двухфазная термодиффузионная задача Стефана о замерзании раствора. При образовании пресного льда перед фронтом кристаллизации образуется слой с повышенным содержанием концентрации примеси, что влияет на температуру замерзания. 1.2.1 Математическая формулировка термодиффузионной задачи Стефана Область Qi(t) = {0 z f\(t)} содержит раствор соли (солей). Распределение температуры и примеси описывается уравнениями: dTi 2d2Tt аГ = а {1-2Л) Ж = dM- (L2-2) В образующейся твердой фазе Qs(t) = {fi(t) z /2(і)} предполагаем отсутствие диффузии примеси, а распределение температуры аналогичное: дТв гд2Т, ж=а-м (1-2-3) Здесь Ti(t, z), Ts(t, z) - температура жидкой и твердой фаз соответственно; C(t,z) - концентрация примеси; a], a2s - соответствующие жидкой и твердой фазам коэффициенты температуропроводности; d - коэффициент диффузии в жидкой фазе; Н - глубина заданной области.
Благодаря уменьшению плотности воды при замерзании лед плавает на воде, незначительно возвышаясь над ее поверхностью (Одрова Т.В., 1979). Приняв форму льдины в виде прямоугольного параллелепипеда с основанием, равным единице, согласно закону Архимеда, подсчитаем соотношение ее надводной и подводной частей. Обозначим толщину надводной части льдины через I, а общую толщину льдины ls, условие равновесия запишем как: kps = (I - Qph где левая часть характеризует силу, стремящуюся утопить лед, а правая - выталкивающую силу. Отсюда можно определить и общую толщину льда и соотношение надводной и общей толщины L — Pl р h Pi Т.е. над водой будет возвышаться (1 — kp)ls общей высоты льдины, кр = ps/pi - коэффициент, характеризующий плавучесть льда. В области, содержащей две фазы, имеются две подвижные границы, положение которых описывается следующими выражениями: fi(t) = її = Н — kpls - граница между жидкой и твердой фазой, І2 (f) — к + h внешняя граница, мало перемещаемая за счет разности плотностей твердой и жидкой фаз. Здесь /;, 18 - толщина жидкой и твердой фаз, соответственно. 1.2.2 Граничные и начальные условия На границе z = 0 для температуры и примеси определены: условие изоляции (отсутствие притока тепла и примеси): дТх (1.2.4) = 0. п 9С dz z=0 OZ 2=0 Можно задать закон теплообмена со средой, или постоянное значение. На границе z = f\(t) выполняются условия сопряжения: 1) условие Стефана, которое описывает тепловой баланс: ЬріаУ/ = h дТі dz идТ z=h(t) OZ z=h(t) (1.2.5) dls 2) равенство температур на подвижной границе температуре замерза ния: Tl\z=fi(t) = Ts\z=h{t) = Т/, Tf = T - yCf; 3) баланс массы растворенного в воде вещества: ц/1 ,дС dfi z=h (1.2.6) (1.2.7) На внешней границе z — /2( ) происходит охлаждение раствора Ts\z=f2 = Ta(t) Tf. (1.2.8)
Здесь Ta(t) - внешняя температура, температура на внешней границе z = І2Іі)\ Tj(t) - неизвестная температура фазового перехода; Vf(i) скорость нарастания льда; Т - температура замерзания (кристаллизации) чистого вещества; Cf(t) - неизвестное значение примеси на границе раздела фаз; у - равновесный коэффициент распределения примеси; ki, ks - соответствующие коэффициенты теплопроводности жидкой и твердой фаз; Л - скрытая теплота плавления; pi, ps - соответствующая жидкой и твердой фазам плотность.
Кроме граничных условий и условий сопряжения формулируются начальные условия. Начальное распределение температуры по толщине в двух фазах задается по линейному закону или равное константе. Начальное значение концентрации примеси, как правило постоянно С — CQ.
В каждой из областей строилась равномерная сетка. Шаг по времени брался постоянным. Так как температура замерзания Т/ и солености Cf на границе z = f\(t) заранее неизвестны, то здесь классический метод прогонки использовать невозможно. Кроме того, в связи с тем, что начальная стадия процесса приводит к возникновению больших градиентов температуры и солености, то применение методов ловли фронта в узел сетки или сквозного счета вызывает большие трудности. Поэтому для численной реализации задачи была разработана модификация метода (Воеводин А.Ф., Леонтьев Н.А., Петрова А.Г., 1982), основанного на методе встречной прогонки, позволяющего эффективно использовать бе-зитерационные (точные) методы решения разностной задачи на верхнем временном слое, несмотря на ее нелинейность.
Моделирование динамики роста ледового покрова минера лизированного водоема
На обширной территории страны, в частности, на юге Западной Сибири имеется значительное число минерализированных озер. Как показали проведенные исследования, даже слабая минерализация оказывает существенное влияние на процесс льдообразования. А именно, при формировании ледового покрова в минерализированном водоеме в результате образования практически пресного льда в воде перед фронтом кристаллизации формируется слой с повышенным содержанием химических веществ, что существенно влияет на температуру фазового перехода, а тем самым и на весь процесс формирования ледового покрова (И.Н.Трегуб, В.П. Власов, В.И. Васильев), что необходимо учитывать при выполнении расчетов. Для описания ледового режима соленых озер на базе термо-диффузиошюй задачи Стефана была разработана следующая модель. Здесь T - температура замерзания чистой воды, Tf(t) - фазовая температура, Cf{t) - значение примеси на границе раздела фаз, j - равновесный коэффициент распределения примеси, зависящий от химических параметров растворенного в воде вещества (Гороновский И.Т. и др. Краткий справочник по химии, 1987), Л - скрытая теплота кристаллизации воды.
Классическим стал подход норвежского ученого О. Дэвика (Винни-ков С.Д., Проскуряков Б.В., 1988) к процессу льдообразования. В предложенной им формуле учитывается теплообмен на верхней поверхности снежного покрова, теплофизические свойства снега. В дальнейшем опи сание процесса теплообмена снежно-ледовой поверхности с атмосферой было усовершенствованно В.В. Пиотровичем (1968, 1969, 1970). Приведем описание методик О. Дэвика и В.В. Пиотровича.
Климат района Новосибирского водохранилища определяется особенностями циркуляции атмосферы над Западной Сибирью и обладает чертами резко континентального. Характерной чертой является холодная суровая зима. По календарным срокам зимний период в районе новосибирского водохранилища в среднем продолжается со второй декады ноября по первую декаду апреля. Начало зимы характеризуется ледоставом на акватории и устойчивым снежным покровом на окружающей суше, увеличением дней (до 25 за месяц) со среднесуточной температурой ниже 0С и резким увеличением числа дней с метелями.
Необходимым условием начала льдообразования является расходование тепла и понижение температуры воды в водоеме до О С. Такое состояние, прежде всего, наступает в поверхностном слое и на мелководных участках вблизи урезов воды. Первые ледяные образования в виде заберегов и сала на водохранилище появляются в конце октября - первой декаде ноября. Наиболее раннее появление ледяных образований отмечено в середине октября, наиболее позднее - в середине ноября.
Ледяной покров в течение всей зимы сплошной. Распределение толщины льда по акватории водохранилища и интенсивность его нарастания зависят от гидрометеорологических условий. Если замерзание происходит быстро, в течение 2-4 дней, и при тихой погоде и устойчивых моро зах - начальная толщина льда сравнительно одинакова по всему водохранилищу, а поверхность льда в этом случае гладкая. Если замерзание происходит при резких колебаниях температуры воздуха и сильном ветре, образуются торосы и толщина льда неравномерна. Льдообразование растягивается на длительный период - до 12-14 дней. Замерзает в этом случае прежде всего верхняя часть водохранилища - от села Соколово до села Ордынское. Позднее на 5-6 дней замерзает нижний озеровидный участок, а в конце первой декады ноября на всем водохранилище обычно устанавливается сплошной ледостав. Наиболее интенсивное нарастание льда наблюдается в первый период замерзания, при отсутствии снега на льду. Толщина льда к моменту установления ледостава составляет в среднем 15-20 см (по годам колеблется от 7 до 35 см). В течение зимы происходит увеличение толщины льда с интенсивностью 15-35 см в месяц с начальный период (с ноября по январь) и 6-15 см в месяц и остальную часть зимы. К концу марта - началу апреля нарастание обычно прекращается. Средняя продолжительность ледостава изменяется от 160 дней в речной части до 175 в приплотинной части.
Наблюдения с помощью электроледомерных реек показали, что нарастание льда происходит с нижней поверхности. В среднем за многолетний период толщина льда речной части водоема из-за больших скоростей течения на 5-10 см меньше, чем в озерной.
Расчет толщины ледяного покрова по методике В.В. Пио-тровича
Оз. Чаны - самый большой естественный водоем Западной Сибири, находится на юге Западно-Сибирской низменности, в центральной части Барабинской степи. Озеро Чаны бессточное, заполнено горько-соленой водой, имеющей значительную минерализацию (до 7 г/дм3). Правильнее считать его системой озер, поскольку оно включает в себя озера Большие и Малые Чаны, Яркуль. Одной из гидрохимических особенностей озЧаны является большая неоднородность минерализации воды по акватории. По химическому составу вода озера относится к высокоминерализованным хлоридным водам. По минерализации вода озера относится к типу солоноватых. Минерализация воды в различные сезоны года неоднородны. Наибольшая минерализация наблюдается к концу зимы (март) и колеблется от 1.76 до 18.7 г/дм3. Озеро издавна славится рыбой, оно дает около 70%, а в отдельные годы около 80% годовой рыбной продукции по Новосибирской области.
Изучение зимнего режима озера проводилось главным образом в связи с заморными явлениями, приводящими к массовой гибели рыбы, осо бенно в годы низкого стояния уровней. Эти исследования имеют в основном биогидрохимическую направленность. Вопросы ледового и термического режима, как правило, освещены в литературе очень схематично. Между тем, изучение особенностей термики зимнего периода оз. Чаны имеет важное значение для выяснения гидрофизических и биогидрофизических процессов, протекающих в солоноватых водоемах в то время, когда они покрыты льдом. Зимний режим озера связан со сложными физическими процессами, развивающимися под воздействием гидрометеорологических факторов. В условиях сравнительной мелководности большая лощадь зеркала озера способствует активному теплообмену водоема с окружающей средой. Вследствие этого процессы осеннего охлаждения и замерзания на озере протекают относительно быстро и одновременно. К началу ледостава температура на дне озера мало отличается от температуры у его поверхности. Ветровое перемешивание приводит к установлению почти полной гомотермии, с температурой близкой к температуре замерзания. Температура у дна может держаться постоянной. В условиях суровой зимы, когда ледовый покров может достигать толщины более 1 м, значительная часть водной массы переходит их жидкой фазы в твердую. Минерализация воды только вследствие льдообразования возрастает на 30 - 50% и зависит от морфометрии плесов. На особенно мелководных участках минерализация может увеличится в несколько раз. Увеличение минерализации подледной воды обуславливает понижение температуры ее кристаллизации и приводит в целом к еще большему охлаждению водных масс. Различная степень минерализации разных районов озера в зимний период вызывает также пространственную термическую неоднородность водных масс. Так как оз. Чаны является бессточным водоемом, в который впадают реки Чулым и Каргат, то по мере удаления от зоны опресняющего влияния этих рек возрастает степень минерализации воды, а следовательно уменьшается температура замерзания.
Первые ледяные образования в виде заберегов, сала и шуги появляются почти одновременно во всех районах озера, обычно в третьей декаде октября, при понижении температуры воздуха до отрицательных значений, а температуры воды до температуры замерзания.
При дальнейшем понижении температуры воздуха ледовый покров появляется на открытых мелководьях. При скорости ветра до 6 м/с и сравнительно небольших отрицательных температурах воздуха (до —2, —3С) озеро покрывается льдом в течении 1-3 дней. В ноябре - декабре толщина льда на постах озера 50-70 см, а к концу зимы достигает 100-110 см.
Было установлено, что ледяной покров формируется главным образом за счет намерзания с нижней поверхности. Наростание ледяного покрова за счет образования снегового льда происходит редко, лишь в начальный период ледостава при условии большой снеговой нагрузки и небольшой толщины льда. Интенсивнее процесс нарастания льда идет в начальный период замерзания, когда среднее нарастание толщины льда за декаду составляет 9-14 см. К середине зимы рост толщины льда замедляется, становится равным 3-7 см за декаду.
На процесс роста льда на оз. Чаны оказывает влияние снежный покров. В начале ледостава он распределяется достаточно равномерно по всему озеру. Далее происходит частичный перенос снега с открытой поверхности в более защищенные участки. На равномерность снежного по крова влияет большое количество островов и растительность. Устойчивый снежный покров, как правило, устанавливается в первой декаде ноября. Число дней со снежным покровом в году колеблется от 150 до 160. Его высота достигает 24 см.
Основные черты температурного режима глубоководных зон под ледяным покровом характеризуются тем, что отрицательные значения температуры (—0.3, — 0.4С) наблюдались до глубин 3.0 -3.5 м. Между глубинами 3.5 - 4.5 м отмечались отрицательные и положительные значения температуры, ниже температуры была положительной (до 2С). В некоторых случаях нулевая изотерма опускалась почти до дна.
Озеро Яркуль имеет наибольшую глубину среди всех плесов Чано-воской системы. Оно отличается более постоянной температурой у придонных слоев (большая тепловая инерция за счет глубины - 7-8 м, песчано-илистое дно). Отрицательная температура поверхностного слоя вызывается повышением степени минерализации вследствие образования значительной толщины льда и соответствует температуре дальнейшего замерзания воды. Наблюдается увеличение мощности слоя воды с отрицательной температурой до 4.5- б м от поверхности.
