CC BY
75
ПРИКЛАДНАЯ МАТЕМАТИКА
УДК 519.6:517:958:533.7
М.Д. Михайлов, В.В. Трапп
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА САМООЧИЩЕНИЯ ЗАГРЯЗНЕННОГО УЧАСТКА РЕКИ
Главной задачей моделирования водных экосистем является получение адекватного количественного представления о биотическом круговороте и взаимодействии компонент экосистемы в соответствии с общими законами сохранения и биологическими закономерностями.
Неблагоприятная экологическая ситуация, сложившаяся в Томской и ряде других областей, ведет к актуализации проблемы загрязнения водных объектов в результате сброса в них сточных вод от промышленных предприятий и населенных пунктов.
Сегодня необходимо создать экологические нормативы и определить пути сохранения и усиления самоочищающихся механизмов и функций природных систем. Для успешного планирования, оперативного управления и контроля над качеством водной среды, необходимо комплексное описание гидродинамических, гидрохимических и гидробиологических процессов в водоемах. Такое описание проводится с использованием методов математического моделирования, которые позволяют количественно описывать круговорот и взаимодействие компонент, входящих в экосистему загрязненного участка реки.
Целью данной работы является изучение математических моделей, описывающих процесс самоочищения загрязненного участка реки и применение численных методов к решению такого рода задач.
Предлагаются модификации моделей Кемпа и более сложной модели, предложенной в работе [1].
Модификации моделей
Начнем с модификации модели Кемпа [2]. В модели описывается взаимодействие между РК (растворенный кислород),который является одним из основных показателей жизнедеятельности водной среды, и БПК(биохимическое потребление кислорода), за счет которого, в основном, уменьшается содержание РК в реке и, следовательно, ухудшается качество воды. Данная модель учитывает постоянную скорость ухода БПК за счет осаждения или адсорбции, скорость поступления РК за счет фотосинтеза, осуществляемого водорослями, а также скорость поступления БПК в проточ-
ную систему за счет вымывания из придонных иловых отложений.
Модель Кемпа представляет собой задачу Коши для ОДУ первого порядка:
— = кЬ — к2В — Р\
Лі
сЬ
Лі
— -к, Ь — к з Ь + р,
(1.1)
В(0) — В0 , Ь(0) — Ь0 . Модификация модели (1.1) на двумерный случай осуществлялась вводом в нее операторов кон-
д 5 д2 д2
вективного----1— и диффузионного------------1----
дх ду дх2 ду2
членов.
В результате получается система двух дифференциальных уравнений в частных производных:
дВ (дБ дВ\ , т , ^ „ ^ (д2В дБ
—к • Ь—к2 • В—Р+ВВ ■
д
дЬ
и
(
+К •
--------1------
дх ду
\
дЬ дЬ
дх дУ„
——к • ь—кз • ь+р+в
дХ ду2 (дь дЬ
дХ
с соответствующими начальными
Ь(х, у,0) — Ь0 , В(х, у,0) — В0
и граничными условиями
дВ дВ
дх х— 0 дх
дЬ — дЬ
дх х—0 дх
— 0
— 0
дВ
ду
дЬ
ду
у—0
у—0
дВ
ду
дЬ
ду
у—1у
ду2 (1.2)
(1.3)
— 0, — 0.
у—1у
(1.4)
Здесь коэффициенты конвекции де-
фицита кислорода и концентрации органического вещества соответственно, а ВО, В^ - коэффициенты диффузии дефицита кислорода и концентрации органического вещества. Будем рассматривать случай, когда = Ж , ВО = Вь .
Опишем вторую, более сложную, модель. В ней учитывается круговорот органических веществ, концентрация растворенного кислорода, интенсивность солнечной радиации, круговорот
х
биогенных элементов, концентрация фитопланктона и моллюсков [1].
Математически модель описывается следующей системой ОДУ:
сИ -
—^ = А + _с - А1 - swd - к • Ь • (1 - е гч1) ,
С, °ге К ’
= -к •Ь_,
■ = -к • Ь, • (1 - е ^),
ч
гч1 '
с?
Сч
С1
ч
гч1 '
= _5 + аг • к • Ьоге • (1 - е гч )-а-g р Б/ ,
= а1 • (б - ч) + а2 • Р/ - 8Чп -
- ч - ч -а2 • к • I • (1 - е ^) - С/• Бг • (1 - е ч ) -
- с • Б • (1 - е гчт) ,
т т У ’ ?
В
Ж,
= Рг - _С - sw/ - Сг • Бг • (1 - е г/)
СБт Я* Т Б Я - Я* —т =-------Т - _т • Бт
л и
(1.5)
(1.6)
Я
с начальными условиями Ь (0) = Ь , Ь_ (0) = Ь_ , 1, (0) = Ь
ог^у ' Orgo ’ 1Э V ._0 М V _
g(0) = g0 , ч(0) = чо ,
Б/ (0) = / Бт (0) = Бт0.
Здесь Ь - концентрация органического вещест-
ва, Ь_ - концентрация исходного органического вещества, Ь, - концентрация стойкой фракции,
ч- концентрация растворенного в воде кислорода, g- концентрация биогенных элементов, Бт -концентрация биомассы моллюсков, Б/ - концентрация биомассы фитопланктона, А- скорость поступления органического вещества со стоком, 80- скорость отмирания фитопланктона, А1- скорость осаждения на дно органического вещества,
__Ч_
к • Ьогг • (1 - е щ1) - скорость разложения органического вещества, ^С-скорость пополнения извне запаса биогенных элементов,
- _я_
ag • к • Ь • (1 - е щ1) - скорость поступления биогенных элементов из разлагающегося вещества, а • g • р • Б/ -скорость потребления биогенных
элементов фитопланктоном, а1 • (б - ч) -скорость обмена кислородом между водой и атмосферой. Для модификации модели введем операторы
диффузионного
8 8
конвективного --------1-- ]
8х 8у
82 82
-----1----членов, получая двумерную краевую
8х2 8у2
задачу, описываемую системой из семи дифференциальных уравнений в частных производных:
8Х (8Ьгг 8Ьгг Л (82 X 82Ь Л
- + Ж
8х
8у
= ВГ
8х
8у - _±_
+ А + _С - А1 - swd - к • Ь • (1 - е гч1)
ог? 4 '
- + Ж
+ Ж
( 8Ь_ + 8Ь^ Л 8х 8у
Г 81^ + 81^ Л
8х 8у
= ВГ
( 2Ь._
82Ь_ Л
Л. 2
- к • Ь
= Ог
8 (15) Ь Л
2 Ь_, 2 Ь_
х 2
у 2
- к • Ь_, • (1 - е гч1)
8g (
+ Ж
8,
8? + 8% 8х 8у
Л
= Ог
(8;£ + 8^Л
8х2 8у2
+ _д +
+ ag • к • Ьог? ■ (1 -е гч1)-а • g рБ
—+Ж 8
+-
дх ду
= Ог
( 82 ч 82 ч Л
—1 +----------1
ч 8х2 8у2 у
+ а1 •(б - ч) +
+а2 • Р/ - щп-а2 • к • ^ог? * (1 - е ч ) -
- С/ • Бг • (1-е Г/) -Ст • Бт • (1-е Гчт)
-/ "/
СБ/
С,
СБт
т
С,
■ + Ж
( 8Б г 8Б г Л
гчт >
• Б •(
тт (я2 Р Я2 73 Л
+
8х 8у
= Ог
V
82 Б/ 82Б/
8х
2
у
2
+
■ + Ж
+ Р/ - _С - sw/ - С/ • Б/ • (1 - е гч/ )
(82Бт 82Бт Л
(
8Бт 8Бт
+
8х 8у
= Ог
х 2
+
8у2
+
+ Я Т Б Я-Я
+-----Т - _т • Бт-----
и т Я
с соответствующими начальными
Ьо?? (^ у,0) = Ьог? 0 , Ц*(^ у,0) = Ь_ 0
Ь_,(x, у,0) = Ь_,о, г(x, у,0) = г0 ,
ч(х, у,0) = ч0 , Б/ (х, у,0) = Б/0 , Бт (x, у,0) = Бт0
и граничными условиями
(1.7)
(1.8)
8Ьо
8г
8Ьо
8г
*.=1
+
+
2
8
+
ч
ч
ч ч
г=0
дЬ,
ді
дЬ.
ді
— 0, дЬ*
ді
дЬ„
д^
ді
дВ
ді
ді
—0
ді
дВт
1—1
дВ/
ді
ді
дВт
ді
і—іі
і—іі
дд ’ ді
— 0, — 0,
і—0 дд ді
ді —0
— 0,
і — (х, у).
(1.9)
Численная реализация двумерных моделей
Для решения задач (1.2)-(1.4) , (1.7)-(1.9) используется метод дробных шагов и неявная разностная схема переменных направлений [3]. Область, в которой ищется решение, покрывается равномерной сеткой с шагами кх, ку. Для реализации схемы переменных направлений наряду со значениями искомой сеточной вектор-функции
и (х, і) на слоях п и п+1, которые обозначим и
—п+1
и и соответственно, вводится промежуточное
значение, т.е. значения и при і=іп+ т/2. Переход от слоя п к слою п+1 совершается в два этапа с шагом т/2. Схема для расчетов:
и
— и
0,5т
П + у» п
+ Ж(А1и /2 +Л2и ) —
(2.1)
— ВВ •Л1и /2 + ВВ -Л2и +ф ,
п+1 п+
и — и /2 — п+—п+1
• + Ж(А1и /2 + А 2 и ) —
0,5т
п+ У- п+1 п
— ВВ •Л1и /2 + ВВ -Л2и + р ,
(2.2)
где
/' п г п г п г п
А / — 3,к 3—1,к а / п _______________________ • 3,к ^ 3,к—1
х
/■ п гу г п . г п
3'+1,к — 213 ,к + 13—1,к
и2
/ п г) /* п I £ п
— 3 ,к+1 13 ,к 13 ,к—1
Л 21 — И1
(2.3)
р - значения правых частей уравнений на слое п в точке (х],ук ).
Для системы (1.1) - и = (Г, Ь)Т, для системы (1.7) - и = (Ь , Ь_ , Ь_, , ?, q, Бг , Бт )Т .
Схема (2.1) неявна по направлению х и явна по у, схема (2.2) явна по х и неявна по у.
Запишем разностную аппроксимация для краевых условий для задачи (1.2) - (1.4)
В(х у,0) — В0(х у),
Ь(х,у,0) — Ь0(х,у), (х,у) (2.4)
п+1 п п+1 п
В3,0 _ В3,1 , Ь3,0 _ Ь3,1 и для задачи (1.7) - (1.9)
Ь01 (Х, у,0) — Ь010(Х, у) ,
К(x,у,0) — Ь,0(x,у),
Ьі(x, у,0) — Ьі 0( x, y), і(x, у,0) — 10( x, у), д( x, у,0) — д0( x, у),
В1(х у,0) — В/ 0(х у),
Вт (^ у,0) — Вт0 (^ У), (Х, у) (25)
-г п + 1 _ -г п -г п + 1 _ -г п
ЬогЯ3,0 — ЬогЯ3,1 , і$ і ,0 — і ,1 ,
п+1 п
Ьїі3 ,0 — Ьїі3 ,1
3 — іпіл, 3 — 3
п+1 п п+1 п
В г — В г , В Л — В .,.
/3,0 / 3,1 ’ т3,0 т3,1
Реализация схемы переменных направлений основана на последовательном определении зна-
п+ X п+1
чений и и и . Решение находится методом прогонки сначала по одному направлению (по переменной х), потом по другому направлению (по переменной у).
Для исследования вопросов аппроксимации, устойчивости и сходимости численного решения задач (1.2)-(1.4) , (1.7)-(1.9) использовалась теория, предложенная в [4].Согласно этой теории, рассматриваемую схему можно записать в опера-торно- разностном виде:
Е + 2А(1)Іи /2 —І Е — ^А(2) и + т • /хп
Е + тА(2) |ип+1 —І Е — -А(1) Ь"^2 + т/2п,
2
<(“) _
Т . (1) и ~п+ V-
2 )
А =Ла+Ла, а = 1,2.
(2.6)
Сформулируем условия устойчивости схемы переменных направлений в виде теоремы.
Теорема. Пусть в схеме (2.6) постоянные операторы А(а) > 0 ,а = 1,2. Тогда для разностного решения имеет место следующая оценка устойчивости по начальным данным и правой части:
■х: |/11 Ч /Ц)
к=0
(2.7)
Доказательство теоремы рассмотрено в [4].
(Е +Т А(2) !ип+1 < (Е +Т А(2) ^и °
1 2 ) 1 2 )
і—0
і—і
і—0
і—0
і—і
і—0
і—0
у
п
Рис. 1. Изменение концентрации органического вещества (точечная модель)
сут
Рис. 2. Изменение концентрации биомассы фитопланктона, биомассы моллюсков, растворенного
в воде кислорода (точечная модель)
Схема переменных направлений имеет второй порядок аппроксимации по времени и по пространству, однако, использование мягких горничных условий дает суммарную аппроксимацию первого порядка по времени и по пространству. Тогда, согласно теореме Лакса, имеет место сходимость [5].
Обсуждение результатов
Для численной реализации точечных моделей
(1.1), (1.5)-(1.6) применялся неявный метод Эйлера. В качестве органического вещества рассматривался фенол, т.к. он является основным из загрязнителей р. Томь. Известно, что предельно допустимая концентрация (ПДК) фенола составляет
0,001 мг/л [6]. Начальная концентрация Ь0 = 0,28мг / л, Г0 = 0 мг / л,
Ь = 0,28мг / л,
огг 0 ’ ’
Ьт — 1мг / л, Ьо00 — 1мг / л, — 1мг / л,
д — 8мг / л, В/ — 1мг / л, Вт0 — 1мг / л.
(3.1)
Расчеты проводились для 14 суток. Согласно полученным результатам (рис.1), по модели Кемпа
(1.1) ПДК фенола достигается на 6 сутки, а полное очищение участка реки от загрязнителя достигается к концу 6 суток. По второй, более сложной модели (1.7), ПДК достигается к концу 3 суток, полное очищение происходит на 4 сутки.
Для оценки работы моделей также проводились расчеты при загрязнении 30 мг/л. В этом случае очищение произошло быстрее по модели Кемпа (на 8 сутки), а по модели (1.5)-(1.6) на 8 сутки концентрация органического вещества составляла 27,2 мг/л . Дефицит кислорода (модель
(1.1)) достигает наибольшего значения - 0,027 мг/л к концу 1 суток, а затем падает до нуля.
На рис.2 показано изменение концентрации для ч (растворенного в воде кислорода), Бп
(биомассы моллюсков), Б/ (биомассы фитопланктона).
На рис.3 - изменение концентрации для Ь._ (исходного органического вещества), Ь, (стойкой фракции), г (биогенных элементов).
^ сут
Рис. 3.Изменение концентраций исходного органического вещества, стойкой фракции и биогенных элементов (точечная модель)
Рис. 4. Концентрация органического вещества (Ъог& мг/л)
(одномерная модель) (х, км; 1, сут)
Из рис. 2 видим, что биомасса фитопланктона и биомасса моллюсков растут;
на 14 сутки биомасса фитопланктона достигает 8,76 мг/л, моллюсков - 2,03 мг/л. Концентрация растворенного в воде кислорода достигает 1,2 мг/л.
Концентрация исходного органического веществ падает на 14 сутки до 0,856 мг/л, концентрация стойкой фракции на 14 сутки равна 0,8578 мг/л. Изменение концентрации биогенных элементов почти не происходит.
Наряду с точечными моделями и их модификациями на двумерных случай, описанными в п.1, п.2, проводились расчеты для модифицированных моделей (1.1), (1.5)-(1.6) на одномерный случай. Задача самоочищения решалась с помощью неявной схемы первого порядка по И (шаг по пространству) и первого порядка по т (шаг по времени).
Рассматривался участок реки длиной 100 км. Скорость течения реки бралась равной 120 км/сут, а коэффициент диффузии - 2,5 км2/сут. Начальное
■ Начальная концентрация
концентрация на первые сутки
- концентрация на 2 сутки
- концентрация на 3 сутки
концентрация на 4 сутки
Рис. 5. Изменение концентрации органического вещества (одномерная модель)
загрязнение Ь0=Ь,
ог%0
0.28 мг/л, начальное зна-
чение стойкой фракции Ь_,0=1 мг/л на первых 4 км реки, на остальном участке они полагались равными нулю. Остальные параметры задавались равномерно на всем рассматриваемом участке.
На рис. 8-11 показано поведение Ьогг , Ь_, , Ь,
Рис. 6. Изменение концентрации стойкой фракции (Ь%1,мг/л)(одномерная модель) (х, км; I, сут)
ч, О для одномерной пространственной модели. Рис. 4-5 и рис. 6-7 показывают изменение концентрации органического вещества и стойкой фракции в течении 14 суток. Из рисунков видно, что обе величины убы-
вают с течением времени. Значение ЬоГг падает до нуля на 4 су-
1,1
0,9
0,7
£
1_
0,5
ІЛ
-1
0,3
0,1
-0,1
х, км
Рис. 7. Изменение концентрации стойкой фракции (одномерная модель)
♦ начальная концентрация
—концентрация на 1 сутки
—Д— концентрация на 2 сутки
—концентрация на 3 сутки
—концентрация на 4 сутки
—Ж— концентрация на 8 сутки
х, км
Рис. 9. Изменение концентрации органического вещества (одномерная модель)
100
14 0
Рис. 11. Изменение растворенного в воде кислорода (д, мг/л) (одномерная модель) (х, км; 1 сут)
50 5
Рис.12. Распределение органического вещества (Ьогі, мг/л). Начальная концентрация (двумерная модель) (х, у, км)
тки, Ь_, с течением времени убывает медленно и к 10 суткам уменьшается до 0,88 мг/л.
На рис. 8-9 и 10 приведено изменение концентрации органического вещества (Ь) и дефицита кислорода (Г) для одномерной пространственной модификации модели (1.1) в течении 30 суток. Качественно картина поведения Ь совпадает с поведением ЬоГг для более сложной модели. Заметим только, что в модели Кемпа процесс самоочищения идет медленнее. Дефицит кислорода в первые сутки, когда концентрации органического вещества в реке достаточно велика, увеличивается, а затем начинает уменьшаться до нуля.
На рис. 12-14 приведены значения Ьогг в различные моменты времени для двумерной модели. Рассматривался прямоугольный участок
реки длиной 100 км и шириной 5 км. Скорость течения реки бралась равной 120 км/сут, а коэффициент диффузии 2,5 км2/сут.
Усложненная модель создавалась для описания процесса самоочищения, когда концентрация органического вещества достаточно высока.
В данной работе рассматривалось Ьогг=30 мг/л и Ьог?=0.28 мг/л. Анализ полученных результатов показал, что усложненная модель одинаково хорошо работает как для больших, так и для малых концентраций загрязнения. При малых концентрациях органического вещества расчеты по моделям Кемпа и усложненной дают близкие результаты. При больших концентрациях органического вещества модель Кемпа дает результаты, не соответствующие действительности.
СПИСОК ЛИТЕРАТУРЫ
1. Биологические процессы и самоочищение на загрязненном участке реки (на примере Верхнего Днепра). - Минск.: Изд-во БГУ, 1987.
2. Белолипетцкий В.М., Шокин Ю.И. Математические модели в задачах охраны окружающей среды.
- Новосибирск: Издательство «ИНФОЛИО - пресс», 1997.
3. Самарский А .А. Введение в теорию разностных схем. - М.: Наука,1971.
4. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. - М.: Эдиторал УРСС, 1999.
5. Вержбицкий В.М. Основы численных методов: Учебник для вузов. -М.:Высш. шк., 2002.
6. Экологический мониторинг. Состояние окружающей природной среды Томской области в 1998 году. - Томск, 1999.
□ Авторы статьи:
Михайлов Михаил Дмитриевич
- ст. преп. каф. вычислительной математики и компьютерного моделирования Томского Государственного университета
Трапп
Виктория Викторовна
- дипломант каф. вычислительной математики и компьютерного моделирования Томского Государственного университета
УДК 622.1:528 (043.2)
В. И. Снетков
КРИТЕРИЙ ОПТИМАЛЬНОСТИ ПРИ АППРОКСИМАЦИИ ГЕОЛОГИЧЕСКИХ ПОКАЗАТЕЛЕЙ С ПОМОЩЬЮ РЯДОВ ФУРЬЕ
Для оценки тенденций в изменении геологического показателя по линии или в пространстве обычно применяют фильтрующие или аппроксимирующие процедуры, позволяющие уменьшить влияние случайной изменчивости. К ним относятся различного рода скользящие средние, полиномиальные функции, сплайны, гармонический анализ Фурье и другие [1-3, 5-7, 9]. Наибольшее распространение в геологии и горном деле получил метод Фурье, однако при его практическом применении существуют трудности, связанные с выбором числа гармоник, необходимых для правильного описания закономерного изменения изучаемого показателя.
Как известно, анализ Фурье применяется при исследовании временных процессов и представляет собой метод разложения сигнала или случайной функции на отдельные гармоники. В геологии и горном деле принято считать сигнал состоящим из трех частей: линейного тренда, различных периодических или циклических компонент и случайной компоненты. Тренд данных обнаруживается подбором соответствующей прямой регрессии. Считается, что тренд существует, если угол наклона прямой является значимым. Отметим, что дисперсионный анализ не всегда способен правильно решить эту задачу, и поэтому практика показала, что тренд необходимо всегда устранять вне зависимости от значимости коэффициентов регрессии. Устранение тренда обязательно для применения метода Фурье для того, чтобы осуществить преоб-
разование исходного временного ряда в последовательность отклонений от подобранной прямой линии. Оставшаяся часть временного ряда состоит из сигнала (периодической компоненты) и шума (случайной компоненты).
Итак, одной из основных задач анализа Фурье является выделение главной периодической компоненты, что эквивалентно выделению закономерного изменения показателя из композиции в случайном процессе.
Основные положения анализа Фурье и принципы его применения следующие. Любая периодическая функцию может быть разложена в так называемый ряд Фурье [2,3], который для простых последовательностей имеет вид:
^ ( 2пжх- Уп’лхс-^
Ъ =^\апСо*-ПХ + ^п-Я ,
п=Л Я Я '
где п - гармоническое число; ОСп, /Зп - коэффициенты гармоник; Я - длина волны основной гармоники.
В действительности, обычно величина Я неизвестна и в качестве основной длины волны выбирают произвольное значение Ь. Если выбрать Ь большим или равным длины изучаемого ряда данных и вычислить достаточное число гармоник, то получим оценку периодичности, присутствующей в наших данных. Произвольный выбор значения Ь скорее всего будет неправильной
