CC BY
35
Вычислительные технологии
Том 13, Специальный выпуск 1, 2008
Разработка нестационарной математической модели вторичного загрязнения водного объекта с учетом материального обмена с донными отложениями*
В, В, КОЗЛОВ
Институт динамики систем, и теории управления СО РАН, Иркутск, Россия
e-mail: [email protected]
This paper describes a non-stationary mathematical model related to the processes of formation of the diffusion flow from the bottom sediments and the flow exchange inside the "bottom — water" system. The model accounts for the molecular diffusion of the decay products in the porous medium including macro-kinetics of such processes. Results of computations of several test problems are discussed.
Введение
Почти повсеместное обострение экологических проблем в водных объектах делает актуальным исследование возможности количественной оценки величины диффузионного потока растворенных органических загрязняющих веществ (ОЗВ) из донных отложений [1-3] и его влияния на качество воды водного объекта. Трудности организации полноценных натурных наблюдений за такими процессами заставляют прибегать к их математическому моделированию [3-7].
Органические загрязняющие вещества поглощаются биотой водного объекта, сорбируются на взвесях, поступают в донные отложения. Это приводит к уменьшению концентрации в воде ОЗВ. С другой стороны, донные отложения становятся длительным источником вторичного загрязнения. При этом надо учитывать, что ОЗВ подвержены биохимической трансформации в донных отложениях.
Предполагается, что значительные изменения величины и направления диффузионного потока растворенного ОЗВ через поверхность раздела вода — дно происходят лишь во время протекания переходных процессов (например, барьерная роль окисленного слоя донных отложений, препятствующая свободному диффузионному обмену растворенными веществами на границе вода — дно). Эти переходные процессы обусловлены сменой условий во внешней по отношению к системе придонная вода — донные отложения среде и приводят к значительным потокам ОЗВ в поверхностном слое осадков, В настоящей работе рассматривается формулировка математической модели и приводится численный алгоритм решения нестационарной краевой задачи молекулярной диффузии ОЗВ в пористой среде с подвижными границами. Приводятся результаты расчетов на модельном примере,
* Работа выполнена при финансовой поддержке Российского гуманитарного научного фонда (проект № 06-02-00055а), Российского фонда фундаментальных исследований (грант № 07-06-12023), Междисциплинарного интеграционного проекта СО РАН № 40.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
1. Постановка задачи
Донные отложения представляют собой естественную насыщенную водным раствором различных ОЗВ пористую среду, формирующуюся в водоеме путем последовательного сложения большого количества осаждающихся на дно частиц различных размеров, генезиса и состава. Эта среда ограничена снизу коренным ложем водоема, а сверху — поверхностью раздела с водной толщей.
Поскольку точное количественное описание структуры реальной пористой среды и диффузионного переноса в ней сопряжено со значительными трудностями, при макроскопическом подходе она обычно рассматривается как некоторая фиктивная гомогенизированная среда — континуум. Эффективные параметры данного континуума выбираются таким образом, чтобы обеспечить равенство диффузионных потоков в нем и в реальной пористой среде при равных прочих условиях.
Вводится локальная подвижная система координат, жестко связанная с границей раздела вода — дно. Начало координат совмещено с этой границей, положительная ось Ог направлена вниз, в толщу осадков. Поскольку изменение концентрации на единицу длины по вертикали в донных отложениях значительно больше изменения величины концентрации по другим направлениям, диффузионный процесс будем рассматривать как одномерный, В ходе осадконакопления (его скорость У,) система координат перемещается вверх параллельно самой себе со скоростью Уь.
1.1. Уравнение молекулярной диффузии ОЗВ в пористой среде
В случае переменных коэффициентов диффузии, пористости и скорости осадконакопления с учетом обратимости процессов сорбции и десорбции основное уравнение макрокинетики процессов в донных отложениях для ОЗВ (оно широко используется для описания кинетики химических реакций, осложненных массопереносом за счет молекулярной диффузии [1, 8, 9]) имеет вид
d[G(z,t)Cb(z,t)] _ д_ dt dz
д
- — [G(z, t)Vb(z, t)Cb(z, t)] + fb(z, t), (1)
где D34) = D0pm, 1, 3 < m < 3 — эффективный коэффициент молекулярной диффузии, p — пористость, Do — коэффициент молекулярной диффузии растворенного вещества в воде при бесконечном разбавлении [8], В общем случае fb(z,t) описывает совокупное действие источников — химических реакций, физико-химических и биологических процессов, влекущих за собой изменение во времени концентрации рассматриваемого ОЗВ [1, 7, 9], а также процесс распада исходного органического вещества. Здесь под Cb(z,t) понимаем концентрацию растворенного ОЗВ в пористом осадке, выраженную в единицах химического потребления кислорода (ХПК), При значительной скорости сорбционного процесса принимается, что равновесие между жидкой и твердой фазами устанавливается мгновенно и описывается линейной изотермой. Процесс равновесной сорбции описывается функцией G(z, t) = 1 + Ka:
K = , Kai, 0 < z < z*(t),
Ka ^ Ka2, z*(t) < z < L.
Здесь г* (£) — толщина (мощность) окисленного слоя донных отложений.
Плотность полного потока вещества в донных отложениях вычисляем по формуле
дСь (М) ■
Fb(z, Ь)полн
дх
+ [G(z,í)Vь(z,í)Cь (М)]
Выражение для полного потока через поверхность донных отложений имеет вид
ВД) = ^ (х,*)пол н|г=о-
Начальное вертикальное распределение концентрации растворенного органического вещества в осадке Сь(х,Ь)1=0 = Ф(х),
Функция Ф(х) находится идентификацией параметров [10] или получается методом установления при численном решении.
На нижней границе антропогенного слоя донных отложений с координатой х = Ь градиент концентрации равен нулю:
дСь (х,Ь)
дх
0.
х=Ь
Граничное условие на границе вода — дно Сь(х,Ь)|^=0 = С(ж*,£). Здесь ж = ж* — положение (координата) в водном объекте точки расчета потока из донных отложений, концентрация С (ж*, Ь) находится из решения уравнений блока качества воды (в расчетах принята одна из модификаций модели Стритера — Фелпса [4, 5, 11] с учетом процессов сорбционного обмена).
1.2. Поглощение кислорода донными отложениями
В отличие от ОЗВ кислород не имеет положительных источников в толще донных отложений, Он поступает в осадки через поверхность вода — дно. Изменение содержания кислорода в водах придонного слоя неглубокого водоема вызывает сопряженные колебания концентрации кислорода в поровых растворах верхнего слоя осадков. Это влечет за собой изменение во времени глубины проникновения кислорода в донные отложения и окислительно-восстановительной обстановки в них, что не может не отразиться на течении всего комплекса физико-химических процессов в верхнем горизонте осадков.
Диффузионное проникновение кислорода Ск(х, Ь) в донные отложения неглубокого водоема, сопровождающееся его поглощением, скорость которого пропорциональна его текущей концентрации, запишем в следующем виде [1, 3]:
дСк ОМ) = д дЬ дх
д
- —[н(м)аом)] + <мм)- (2)
В начальный момент времени в колонке донных отложений задано вертикальное распределение концентрации растворенного кислорода Ск(х, ¿)^=0 = ^(х), функция ^(х) находится идентификацией параметров или получается методом установления при численном решении.
На нижней границе антропогенного слоя донных отложений с координатой х = Ь градиент концентрации равен нулю:
дСк (М)
дх
х=Ь
На границе вода — дно Ок(г,£)|х=0 = С02(х*,£), х = х* — положение (координата) в водном объекте точки расчета потока из донных отложений. Концентрация растворенного кислорода в придонной воде С02 (х*, £) находится из решения уравнений блока качества воды.
Функция Фк(г, £) описывает кинетику поглощения кислорода органическим веществом донных отложений в соответствии с законом действующих масс [9, 11],
1.3. Изменение мощности поверхностного окисленного слоя осадков
На основании экспериментальных исследований [1] показано, что под границей окисленного слоя содержание растворенного кислорода невелико и составляет около Б* = 0.0005 мг/см-3. Это обстоятельство позволяет считать концентрацию 02 на нижней границе окисленного слоя при прочих равных условиях практически постоянной. Тогда изменение во времени положения в осадке линии уровня, отвечающей данной концентрации, описывает динамику во времени толщины поверхностного окисленного слоя, обусловленную колебаниями содержания кислорода в придонной воде: Н(г, £) = Ск(г, £) — Б* = 0
Выражение, описывающее скорость движения нижней границы окисленного слоя, приобретает вид
(г* Н
(И Н'
X
*
(3)
Начальное положение г* ^=0 = г* определяется по начальному вертикальному распределению растворенного кислорода <^(г).
:
Х=Х
1.4. Распад исходного валового органического вещества
Исходное органическое вещество Сорг(г,£) (выраженное в единицах ХПК), захороненное в осадках водных объектов, постепенно распадается, что приводит к уменьшению содержания органического углерода в толще донных отложений от верхних горизонтов к нижним горизонтам. Не отрицая всей сложности и многообразия совокупности реакций распада органического вещества в донных отложениях, полагаем, что в целом этот процесс в первом приближении может быть удовлетворительно описан кинетическим уравнением первого порядка. Концентрация органического углерода в донных отложениях на глубине г0 связана с начальной его концентрацией на этой же глубине уравнением первого порядка
=-кСщ^оЛ (4)
Сор = Со (го) •
Здесь к — константа скорости распада; £0 — время, в течение которого отложился слой осадков мощностью г0.
При постоянной скорости накопления осадков происходит равномерное увеличение толщины слоя донных отложений г = УъЬ. Время пребывания захороненного органического вещества определяется глубиной его залегания в грунтовой колонке. Предельный вариант распределения исходного органического вещества в осадках по вертикали имеет вид Сорг(г) = С° ехр(—к/Уь£).
2. Численный метод
В общем случае получим краевую задачу (1), (2) для уравнений параболического типа с граничными условиями на подвижных границах (3), (4) (граница раздела вода — дно и граница окисленного слоя донных отложений). Для дискретизации по пространственной переменной используются схема с разностями против потока (называемая также схемой с донорными ячейками или разностными уравнениями с положительными коэффициентами) для конвективных членов и центральные разности для диффузионной составляющей [12, 13]. Принцип замороженных коэффициентов [14] позволяет раздельно решать задачу для распределения концентрации ОЗВ и затем определять положение подвижной границы. На каждом интервале времени т = £п+1 — Ьп коэффициенты, входящие в систему дифференциальных уравнений, вычисляются по значениям искомой функции, найденной на предыдущем интервале времени. Для этой задачи рассматривается один из вариантов итерационных методов с набором параметров Чебышева (метод с ЧНП, известный также под названием метода Ричардсона) [12, 14].
Считаем, что отрезок времени [Т0, Т1] разбит та интервалы (¿п-1,£п), в пределах которых коэффициенты системы удовлетворяют условиям устойчивости и сходимости. За основу численной реализации задачи на интервале времени т = Ьп — £п-1 берется явная схема локальных итераций .111-М [15, 16]. Результаты расчетов для тестовых примеров [1] показали, что эта схема для параболического уравнения с постоянными коэффициентами является устойчивой, имеет первый порядок аппроксимации по времени, порядок аппроксимации по пространству совпадает с порядком аппроксимации оператора задачи.
3. Результаты расчетов для донных отложений
В настоящее время весьма актуальным становится расчет и прогноз изменения содержания растворенных веществ в водных объектах. Эту задачу рассмотрим на примере изменения в водном объекте содержания ОЗВ, которое поступает в водный объект из донных отложений [11]. Рассмотрим упрощенную гидрологическую модель, представляя водный объект одной камерой (однокамерная модель) при заданных постоянных объеме и расходе. Задано изменение концентрации кислорода в придонной воде. В начальный момент ОЗВ в водном объекте отсутствует. Процессы распада ОЗВ в водном объекте не учитываются. Пористость донных отложений изменяется по линейному закону. Константа сорбции в верхнем окисленном слое в два раза больше величины сорбции в нижнем анаэробном слое. В расчетах принято, что в анаэробных и аэробных условиях скорость распада исходного органического вещества различается в два раза.
Результаты расчетов, представленные на рисунке, показывают, что во время протекания переходных процессов, связанных со сменой окислительно-восстановительной обстановки (кривая изменения концентрации кислорода в придонной воде), значительно возрастает величина суммарного потока ОЗВ из донных отложений. Первый пик (значительное возрастание потока) связан с наступлением анаэробных условий (отсутствие кислорода в придонной воде). Второй (отрицательный) пик — поступление кислорода в придонный слой воды. При этом суммарная величина вторичного загрязнения водного объекта монотонно возрастает (кривая прогноза загрязнения — концентрация условного ОЗВ).
0.05
10
« g -0.01 -
g о \
I- CL. \
£ = -0.02 --
Время протекания процесса в сутках
Прогноз загрязнения водоема органическим загрязняющим веществом при заданном ходе кислорода в придонном слое воды (однокамерная модель): ■ ■ ■ — поток ОЗВ, г/м2 в сутки; --прогноз концентрации ОЗВ, мг/л;--концентрация кислорода, мг/л
Приведенные выводы соответствуют результатам, полученным в работе [1] при реализации двухслойной модели фосфатного фосфора.
Список литературы
[1] Мизандронцев И.Б. Химические процессы в донных отложениях водоемов. Новосибирск: Наука, 1990. 176 с.
[2] Бреховских В.Ф., Вишневская Г.Н., Ганшина H.A. и др. О сезонной смене приоритетных факторов, определяющих интенсивность потребления кислорода грунтами водохранилища долинного типа // Вод. ресурсы. 2003. Т. 30, Jfa 1. С. 61-66.
[3] козлов в.в. Разработка и идентификация нестационарной математической модели распространения загрязняющих веществ в водном объекте // Тр. Междунар. конф. "Вычисл. и информ. технологии в науке, технике и образовании". Павлодар, 2006. Т. I. С. 641-649.
[4] Железняк M.II. Математические модели миграции радионуклидов в каскаде водохранилищ // Системный анализ и методы математического моделирования в экологии. Киев, 1990. С. 48-58.
[5] Васильев О.Ф., Воеводин А.Ф. Математическое моделирование качества воды в системах открытых русел // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики. 1975. Вып. 22. С. 73-88.
[6] Бреховских В.Ф. Проблемы качества поверхностных вод в бассейне Северной Двины.
М.: Наука, 2003. 233 с.
[7] Семчуков А.Н., Квон В.И., Квон Д.В., Зонов С.Д. Математическое моделирование турбулентных течений, процессов переноса тепла и вещества и их биохимической трансформации в равнинных водохранилищах // Матер. Междунар. конф. "Вычислительные и информационные технологии в науке, технике и образовании". Ч. 1. Алматы-Новосибирск, 2004. С. 18-23.
[8] Макрокинетика процессов в пористых средах (топливные элементы) / Ю.И. Чизмаджев, B.C. Маркин, М.Р. Тарасевич, Ю.Г. Чирков. М.: Наука, 1971. 363 с.
[9] МАРРИ Дж. Нелинейные дифференциальные уравнения в биологии. Лекции о моделях. М.: Мир, 1983. 397 с.
[10] Идентификация моделей гидравлики / Г.Д. Бабе, Э.А. Бондарев, А.Ф. Воеводин, М.А. Каниболотский. Новосибирск: Наука, 1980. 160 с.
[11] КОЗЛОВ В.В. Разработка и идентификация нестационарной математической модели распространения загрязняющих веществ в водном объекте. Иркутск, 2006. (Препр. ИДСТУ. № 1).
[12] Самарский A.A. Теория разностных схем. М.: Наука, 1989. 616 с.
[13] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.
[14] Годунов С.К., Рябенький B.C. Разностные схемы. М.: Наука, 1973. 400 с.
[15] Жуков В.Т. Явно итерационные схемы для параболических уравнений // Вопр. атомной науки и техники. Сер.: Мат. моделирование физических процессов. 1993. Вып. 4. С. 40-46.
[16] Зайцев Ф.С., Костомаров Д.П., Курбет И.И. Применение явных итерационных схем для решения кинетических задач // Мат. моделирование. 2004. Т. 16, № 3. С. 13-21.
Поступила в редакцию 25 января 2008 г.
