CC BY
69
УДК 519.6:532.5
А.В. Никитина, И.С. Семенов
МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ЭВТРОФИКАЦИИ МЕЛКОВОДНОГО
ВОДОЕМА
Работа посвящена разработке методов решения модельной задачи эвтрофикации мелководного водоема, учитывающей движение водного потока, микротурбулентную диффузию, пространственно-неравномерное распределение температуры и солености в мелководных водоемах - Азовское море и Таганрогский залив. Численное решение задачи основывается на методе минимальных поправок. Решение задачи эвтрофикации позволит прогнозировать качеством вод мелководных водоемов, а также изучить условий формирования зон с пониженным содержанием кислорода в воде (заморных явлений).
Математическая модель; эвтрофикация; метод минимальных поправок; Азовское море.
A.V. Nikitina, I.C. Semenov MODELING OF PROCESSES OF EUTROPHICATION OF SHALLOW POND
The work is devoted to the development of methods of solution of the model problem of eu-trophication of shallow water body, taking account of the movement of the water flow, микротурбулентную diffusion of spatially-nonuniform distribution of temperature and salinity in the shallow waters, the Azov sea and Taganrog Bay. Numerical solution of the problem based on the method of minimal amendments. Solution of the problem of eutrophication will allow to predict the quality of the water of shallow water bodiesв as well as to examine the conditions for the formation of zones with reduced oxygen content in the water (kills phenomens).
Mathematical model; eutrophication; the method of minimum amendments; the sea of Azov.
Введение. Азовское море - внутреннее море на востоке Европы. Это самое мелкое море в мире, его глубина не превышает 13,5 метров. В силу своей мелко-водности оно наиболее подвержено эвтрофикационным процессам. Известно, что эвтрофикация (др.греч. euxpo9ia - хорошее питание) - это обогащение рек, озёр и морей биогенами, сопровождающееся повышением продуктивности растительности в водоёмах. Эвтрофикация может быть как результатом естественного старения водоёма, так и антропогенных воздействий. К основным химическим элементам, способствующим эвтрофикации, можно отнести фосфор и азот. Рельеф и течения моря находятся в сильной взаимосвязи друг с другом и оказывают влияние на биогенный режим.
Постановка задачи. При построении модели эвтрофикации вод Азовского моря использовались работы [1], [2], посвященные моделированию гидрохимических гидробиологических процессов. Учитывался тот факт, что при штилях и близких к ним ветровых ситуациях возникают анаэробные условия в придонных слоях Азовского моря. Восстановление поверхностного водонасыщенного ила влечет за собой высвобождение в раствор (кроме сероводорода) сульфатов, двухвалентного марганца и железа, органических соединений, аммония, силикатов и фосфатов.
Предположим, что область решения задачи G представляет собой замкнутый бассейн, ограниченный невозмущенной поверхностью моря Z 0, дном ZH = ZH (x, y) и боковой поверхностью О (2 = 2H иZ0 ио).
Модель основана на системе уравнений [3]
dS,. dS,. dS,. , \dS, Л „ d ( dS) (1)
OS OS OS / \OS , „ О I OS I
—L + u—L + v—L + (w - w ■)—L = u.AS. H--1 v.—L I + w.,
Ot Ox Oy 1 Oz ' Oz { ' Oz 1
где Б{ - концентрация ¡-й компоненты; и,V,w - компоненты вектора скорости водного потока; и = (и,V, w); wg - скорость гравитационного осаждения ¡-й компоненты, если она находится во взвешенном состоянии; А - двумерный оператор Лапласа; - источник (сток) данной компоненты или член, описывающий агрегирование (слипание-разлипание), если соответствующая компонента является взвесью; ¡Л1, - коэффициенты турбулентного обмена соответственно по горизонтальному и вертикальному направлениям; у/ - химико-биологический источник, индекс ' указывает на вид субстанции, .е{1,...,15}:
1 - сероводород (И2 Б);
2 - элементная сера (5);
3 - сульфаты (Б04);
4 - тиосульфаты (и сульфиты);
5 - общий органический азот (N);
6 - аммоний (4) (аммонийный азот);
7 - нитриты (02);
8 - нитраты (И03);
9 - растворенный марганец (ВОЫп);
10 - взвешенный марганец (РОМп);
11 - растворенный кислород (02);
12 - силикаты (Б'03 - метасиликат; Б'04 - ортосиликат);
13 - фосфаты ( Р04);
14 - железо (И 2 Б);
15 - кремнекислота (И БЮ - метакремневая; И2Б'0А - ортокремневая).
Пусть п - вектор внешней нормали к поверхности X , и - нормальная по
отношению к X составляющая вектора скорости водного потока. Начальные условия для модели (1) задаются в виде
4=0 = Б'°(^ Ь ' е {^...Д5} , (2)
х = (х, у, г) ,х е О. Добавим граничные условия:
^л с
Б; = 0 на cr, если ип < = о на а, если ип > 0;
дп п
^ = 0 на ^0; ^ = на Ъ , (3)
дг дг
где £. - коэффициент поглощения ¡-й примеси донными отложениями.
Входными параметрами для модели (1) - (3) являются компоненты вектора скорости водной среды, которые описываются гидродинамической моделью [4]—[6].
Построение и исследование дискретной модели. Для численной реализации построенной модели вводится равномерная сетка:
={{П = m, X = ¡К , У1 = , г, = Щ; п = 0..N,' = 0..N, 1 = 0..Иу, Л = 0..N; N,7 = Г,ИХ = 1х,Ыуйу = 1у,ИД = /г} = йТхаКкЛ , где 7 - шаг по времени; кх, йу - шаги по пространству; И, - верхняя граница по времени; Их, Иу - границы по пространству; 1х, 1у, /г - характерные
размеры модельной области О .
Линеаризация задачи (1)-(3) выполнена на основе метода Ньютона. Для дискретизации модели эвтрофикации использовались схемы с центральными разностями следующего вида [7]:
5Д+1 + иБр +1 + vS р +1 + (w - w„ ) Б?р +1 = м (( + Бр) ) + (у,Б(р.)) + N • (4)
()) (г)х (г)у ^ *г> (¡)1 ^Л ()хх (г)уу/ \ ' (г'
Здесь Б1, 1 е {1,...,15} - значения соответствующих функций в узлах сетки на п + 1-м временном слое; Б() = сБ() +(1 -с)Б(), се [0,1] - вес схемы; р - номер итерации в итерационном процессе [5], [6].
Погрешность аппроксимации математической модели эвтрофикации мелководного водоема равна О ( +1\к\\2) в случае с = 1 / 2, где ||й|| = ^й2 + й2 + й2 [8], [9].
Достаточное условие устойчивости и монотонности разработанной модели (1)-(3) определяется на основе принципа максимума при ограничениях на шаг
временным координатам 7 < ш1п {71,72,..., 715} :
Ж)-1,((1 -с)||2М/+ 2М/йу2 + 2к/йг2-у)-1 | /е {1.....15},
где N1=ш N.
Методы решения модельной задачи эвтрофикации мелководного водоема. Представим в виде операторного уравнения дискретную модель (4):
Аи = / (5)
с невырожденным оператором А, заданным в вещественном гильбертовом пространстве И. Рассмотрим неявную двухслойную итерационную схему вида
В Ук+1 - Ук + Аук = /, к = 0,1,... (6)
7+1
с произвольным начальным приближением у е И и невырожденным оператором В.
Любой двухслойный итерационный метод, построенный на основе схемы (6), характеризуется операторами А и В, энергетическим пространством Ип, в котором доказывается сходимость метода, и набором итерационных параметров 7к . Основным вопросом теории итерационных методов является вопрос об оптимальном выборе параметра 7к [10]-[12].
В двухслойных итерационных методах вариационного типа для вычисления параметров 7, не требуется никакой априорной информации об операторах схемы (6) (кроме условий общего вида А = А* > 0,(ОВ1 А)* = ОВ1 А и т.д.), и построение этих методов основано на следующем принципе: если задано приближе-
7 = шип
ние ук, а ук+1 находится по схеме (6), то итерационный параметр Тк+1 выбирается из условия минимума в нв нормы погрешности гк+1 = у +1 — и , где и - решение уравнения (5).
Последовательность ук, построенная по формуле (6), в которой параметры Т выбираются из указанного выше условия, является минимизирующей для квадратичного функционала вида I(у) = (Б(у — и), у — и). Этот функционал
(в силу положительной определенности оператора Б) ограничен снизу, достигает минимума, равного нулю, при решении уравнения (5), т.е при у = и . Выбор параметра Тк+1 из указанного условия обеспечивает локальную минимизацию функционала I(у) при переходе от ук к ук+1, т.е. за один итерационный шаг. В случае явной схемы (В = Е) переход от ук к ук+1 осуществляется по формуле
ук+1 = ук Т к+1Гк ,
где Г = Аук — /.
Для самосопряженного положительно определенного оператора А переход от у к у 1 происходит по направлению — г , которое совпадает с направлением
антиградиента для функционала (А( у — и), у — и) в точке у . По направлению антиградиента происходит наибольшее убывание значения функционала. Параметр Т будем находить из условия минимума в Н нормы погрешности
г = у — и [13].
к+1 -'к+1 -1
Формула для вычисления итерационного параметра тк+1 находилась в предположении, что оператор А не вырожден. Выпишем сначала уравнение для погрешности: ^ = ук — и, к = 0,1, .... Подставляя ук =тк + и в схему (6), получим гк+1 = (Е — ткВ^А)гк, к = 0,1, ..., г0 = у0 — и. Замена гк = Б 2хк позволяет перейти к уравнению, содержащему только один оператор:
X = 5 X , (7)
к+1 к +1 к ' 4 '
где 5 = Е — тС, С = Б 2 (БВ 1 А)Б 2.
Поставленную выше задачу о выборе параметра Т + можно сформулировать, используя равенство =||х^||, следующим образом: выбрать параметр Т+ из
условия минимума нормы х в пространстве Н . Вычислим норму х :
к+г( = ((Е — Т+С) Хк,(Е — Т+С) Хк ) = (8)
= I к IГ— 2тк+1 (Схк, х)+тк2+1 (Схк, Схк) =
-]2
= (Схк, Схк)
(Cxk. xk) (Cxk, Cxk)
+i x 112 _ C' xk)2
(Cxk, Cxk)
Т+1 =-
Так как оператор А не вырожден, то не вырожден и оператор С. Поэтому для любого х имеем (Сх , Схк) > 0 и минимум нормы х +1 достигается при
(Схк, хк) (9)
( Схк , Схк ) Подставляя (9) в (8), получим:
||хк+Л = Рихк||, (10)
где
г = 1 (11)
'к+1 ' (Схк, Схк)(хк, хк)
Итак, формула (9) определяет оптимальное значение итерационного парамет-1
ра т . Подставляя х = Б2 г в (9) получим:
к +1 к к
Тк+1 = ^, к = 0,1,...
к+1 (БВ1Агк, В1Агк)
Учитывая, что Агк = Аук — Аи = Аук — / = гк - невязка, а В~1гк = а>к -поправка, формулу для параметра Тк +1 можно записать в следующем виде:
Т+ = , к = 0,1,..., (12)
(Ба>к ,а>к)
а итерационную схему (6) - в виде явной формулы для вычисления ук +1:
ук+1 = ук —Тк+1^к, к = 0,1, ... (13)
Алгоритм, реализующий построенный метод, можно кратко описать следующим образом:
1) по заданному у вычисляется невязка гк = Аук — /;
2) решается уравнение для поправки Вюк = гк;
3) по формуле (12) вычисляется параметр тк+1;
4) по формуле (13) находится новое приближение ук +.
Описание двухслойных градиентных методов. Рассмотрим теперь частные случаи двухслойных градиентных методов, которые мы будем использовать для решения модельной задачи вида (1)-(3). Каждый конкретный метод определяется выбором оператора Б и имеет свою область применимости. Оператор Б будет выбираться так, чтобы в формулу (12) для итерационного параметра тк+1 входили
только известные в процессе итераций величины [14], [15].
Пусть оператор А самосопряжен и положительно определен в Н . Метод скорейшего спуска (МСС) характеризуется следующим выбором оператора Б : Б = А. Оператор В должен быть положительно определен в Н. Учитывая
соотношения Аг = Аук — / = гк и А = А , из (12) получим формулу для итерационного параметра тк+1 в неявном методе скорейшего спуска:
Тк+1 , к = 0,1,... (14)
(Ао\ ,0\)
Для случая явной двухслойной схемы (6) (В = Е) получим щ = В~1гк =
(гк, гк)
формула для 7,+1 примет вид
7,+1 =^^, к = 0,1,... (15)
( Агк, Гк )
При решении задачи (1)-(3) методом скорейшего спуска по невязке параметр 7к рассчитывается по формуле (15), а метод скорейшего спуска по поправке реализуется с помощью расчетной формулы (14).
В методе скорейшего спуска будем минимизировать норму погрешности
1
+1 в энергетическом пространстве иА :Ц= (Агк,гк)2.
Опишем использование метода минимальных невязок (ММН) для решения модельной задачи эвтрофикации мелководного водоема вида (1)-(3). Этот метод можно применять в случае любого несамосопряженного невырожденного оператора А. Положительная определенность каждого в отдельности операторов А и В не предполагается, требуется лишь положительная определенность оператора
В А. Метод минимальных невязок определяется следующим выбором оператора
О : О = А* А.
Формула (6) для итерационного параметра 7 в методе минимальных невязок имеет вид
_ _ (Ащк, г,)
к = 0,1, ...
(Аюк, Аюк)
В случае явной схемы (6) (В = Е) требуется положительная определенность оператора А, а формула для 7к имеет вид
7 = (Агк, гк) , = 01 7к+1 = (Агк,АГкГ к = 0,!'...
Минимизируем норму невязки оператора О , получим:
\\гк IО=(Огк, гк)=(А* Агк, гк)=11 Агк IГ=I\гк 112.
Итак, для рассматриваемого метода норма погрешности в И О равна норме
невязки, которая вычислялась в итерационном процессе и использовалась для контроля его окончания.
Опишем использование метода минимальных поправок (ММП) для модельной задачи эвтрофикации мелководного водоема вида (1)-(3). Этот метод можно применять для решения уравнения (5) с несамосопряженным, но положительно определенным оператором А. Требуется, чтобы оператор В был самосопряженным, положительно определенным и ограниченным. Метод минимальных поправок определяется следующим выбором оператора О : О = А* В 1А.
Формула (12) для итерационного параметра 7 в методе минимальных поправок имеет вид
7 = (АЩ,Щ) , , = 0,1,... (16)
(В Ащ, Ащ )
В случае явной схемы (6) (В = Е) методы минимальных поправок и минимальных невязок совпадают.
Минимизируем норму поправки в H для выбранного оператора D, получим:
B
II zk IID = D,Zk) = (A*B lAzk,Zk) = O,rk) = (Щ) =11O IIB .
Норма поправки в H может вычисляться в итерационном процессе и ис-
B
пользоваться для контроля его окончания.
Заключение. С помощью экспедиционных исследований проведена первичная верификация модели эвтрофикации Азовского моря. Создан исследовательско-прогнозский комплекс, объединяющий базу данных по Азовскому морю и пакет прикладных программ [16]—[19].
С помощью модели (1)-(3) могут быть описаны такие процессы, происходящие в мелководном водоеме как: аммонификации, нитрификации, нитратредукции (денитрификации), ассимиляции NH4, окисления H2S, сульфатредукции, окисления и восстановления марганца.
Подобные математические модели могут быть использованы для разработки возможных сценариев реабилитации мелководных водоемов с целью восстановления их экосистем до естественного уровня [17].
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Сухинов А.И, Никитина А.В. Математическое моделирование и экспедиционные исследования качества вод в Азовском море // Известия ЮФУ. Технические науки. - 2011.
- № 8 (121). - С. 62-73.
2. Матишов Г.Г., Фуштей Т.В. К проблеме вредоносных «цветений воды» в Азовском море // Электронный журнал «Исследовано в России». - С. 213-225.
3. Деболъская Е.И., Якушев Е.В., Сухинов А.И. Формирование заморов и анаэробных условий в Азовском море // Водные ресурсы. - 2005. - Т. 32, № 2. - С. 171-183.
4. Сухинов А.И., Чистяков А.Е., Алексеенко Е.В. Численная реализация трехмерной модели гидродинамики для мелководных водоемов на супервычислительной системе // Математическое моделирование. - 2011. - Т. 23, № 3. - C. 3-21.
5. Никитина А.В. Модели биологической кинетики, стабилизирующие экологическую систему Таганрогского залива // Известия ЮФУ. Технические науки. - 2009. - № 8 (97).
- С. 130-134.
6. Якушев Е.В., Сухинов А.И., Лукашев Ю.Ф., Сапожников Ф.В., Сергеев Н.Е., Скирта А.Ю., Сорокин П.Ю., Солдатова Е.В., Фомин С.Ю., Якубенко В.Г. Комплексные океанологические исследования азовского моря в 28-м рейсе научно-исследовательского судна "Акванавт" (июль-август 2001 г.) // Океанология. - 2003. - Т. 43, № 1. - С. 44-53.
7. Никитина А.В. Численное решение задачи динамики токсичных водорослей в Таганрогском заливе // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 113-117.
8. Sukhinov A.I., Sukhinov A.A. Reconstruction 0f 2001 Ecological Disaster in the Azov Sea on the Basis of Precise Hydrophysics Models. Parallel Computational Fluid Dynamics, Mutidisciplinary Applications, Prcoceedings of Parallel CFD 2004 Conference, Las Palmas de Gran Canaria, Spain, ELSEVIER, Amsterdam-Berlin-London-New York-Tokyo, 2005.
- P. 231-238.
9. Никитина А.В., Долгой В.Е. Построение пространственно-неоднородных математических моделей в Таганрогском заливе // Известия ЮФУ. Технические науки. - 2008. - № 1 (78).
- С. 177-178.
10. Никитина А.В., Третьякова М.В. Моделирование процесса альголизации мелководного водоема путем вселения в него штамма зеленой водоросли Chlorella vulgaris bin // Известия ЮФУ. Технические науки. - 2012. -№ 1 (126). - C. 128-133.
11. Чистяков А.Е. Теоретические оценки ускорения и эффективности параллельной реализации ПТМ скорейшего спуска // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 237-249.
12. Сухинов А.И., Никитина А.В., Чистяков А.Е. Моделирование сценария биологической реабилитации Азовского моря // Математическое моделирование. - 2012. - Т. 24, № 9.
- С. 3-21.
13. Сухинов А.И., Никитина А.В., Чистяков А.Е., Семенов И.С. Математическое моделирование условий формирования заморов в мелководных водоемах на многопроцессорной вычислительной системе // Вычислительные методы и программирование. - 2013. - Т. 14.- С. 103-112.
14. Сухинов А.И., Никитина А.В., Чистяков А.Е. Восстановление качества вод Азовского моря с помощью численного моделирования // Труды международной науч.-практ. конференции «Преобразование Таганрога - ключ к возрождению России», 29-30 января 2013 г. - С. 135-137.
15. Никитина А.В., Чистяков А.Е., Фоменко Н.А. Применение адаптивного модифицированного попеременно-треугольного итерационного метода для численной реализации двумерной математической модели движения водной среды // Электронный научно-инновационный журнал "Инженерный вестник Дона". - 2012. - C. 4.
16. Никитина А.В., Чистяков А.Е., Семенов И.С. Расчет распространения токсичной водоросли в Азовском море на вычислительной системе с использованием многопоточности в среде Windows. Зарегистрирован в Реестре программ для ЭВМ: N 2012614681 от 25.05.2012.
17. Сухинов А.И., Никитина А.В., Семенов И.С. Расчет для модели взаимодействующих фитопланктонных популяций в Азовском море на вычислительной системе с использованием базы экспедиционных данных. Зарегистрирован в Реестре программ для ЭВМ: N 2012614678 от 25.05.2012.
18. Сухинов А.И., Никитина А.В., Семенов И.С. Расчет для модели взаимодействия планктона и промысловых рыб в мелководном водоеме на вычислительной системе с использованием библиотеки программ эффективного решения сеточных уравнений. Зарегистрирован в Реестре программ для ЭВМ: N 2012614677 от 25.05.2012.
19. Сухинов А.И., Никитина А.В., Чистяков А.Е., Царевский В.В., Фоменко Н.А. Программный комплекс решения сеточных уравнений для трехмерных задач диффузии-конвекции-реакции итерационными методами. Зарегистрирован в Реестре программ для ЭВМ: N 2012614680 от 25.05.2012.
Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.
Никитина Алла Валерьевна - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: nikitina.vm@gmail.com; 347928, г. Таганрог, пер. Некрасовский, 44; тел.; 89515168538; кафедра высшей математики; к. ф.-м. н.; доцент.
Семенов Илья Сергеевич - e-mail: flanker555@yandex.ru; тел.: 89085029807; кафедра высшей математики; аспирант.
Nikitina Alla Valervevna - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: nikitina.vm@gmail.com; 44, Nekrasovsky, Taganrog, 347928, Russia; phone: +79515168538; the department of higher mathematics; cand. of pthis.-math. sc. associate professor.
Semenov Ilya Sergeevich - e-mail: flanker555@yandex.ru; phone: +79085029807; the department of higher mathematics; postgraduate student.
УДК 532.5.031
Н.А. Фоменко
МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ МОРСКИХ ВОЛН НА МЕЛКОВОДЬЕ И ВОЗМОЖНОСТЬ ИХ ФОКУСИРОВКИ
Цель данной работы заключается в разработке двумерной математической модели волновой гидродинамики, учитывающей геометрию рельефа дна и функцию возвышения уровня. Дискретизация задачи по временному координатному направлению выполнена на основе метода поправки к давлению. Для аппроксимации систем уравнений по пространст-
