Научная статья на тему 'Математическое моделирование процесса самоочищения загрязненного участка реки'

Математическое моделирование процесса самоочищения загрязненного участка реки Текст научной статьи по специальности «Математика»

CC BY
255
75
i Надоели баннеры? Вы всегда можете отключить рекламу.

Аннотация научной статьи по математике, автор научной работы — Михайлов Михаил Дмитриевич, Трапп Виктория Викторовна

Предлагаются модификации моделей Кемпа и более сложной модели путем добавления в системы обыкновенных дифференциальных уравнений, описывающих эти модели, операторов конвективного и диффузионного членов. Производится изучение и сравнение указанных моделей с помощью численных методов, основанных как на явных, так и на неявных разностных схемах.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Математическое моделирование процесса самоочищения загрязненного участка реки»

ПРИКЛАДНАЯ МАТЕМАТИКА

УДК 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.

Здесь Ь - концентрация органического вещест-

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

ва, Ь_ - концентрация исходного органического вещества, Ь, - концентрация стойкой фракции,

ч- концентрация растворенного в воде кислорода, 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у - _±_

+ А + _С - А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Б/

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

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

и граничными условиями

(1.7)

(1.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

Реализация схемы переменных направлений основана на последовательном определении зна-

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

п+ 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 суток. Качественно картина поведения Ь совпадает с поведением ЬоГг для более сложной модели. Заметим только, что в модели Кемпа процесс самоочищения идет медленнее. Дефицит кислорода в первые сутки, когда концентрации органического вещества в реке достаточно велика, увеличивается, а затем начинает уменьшаться до нуля.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

На рис. 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пжх- Уп’лхс-^

Ъ =^\апСо*-ПХ + ^п-Я ,

п=Л Я Я '

где п - гармоническое число; ОСп, /Зп - коэффициенты гармоник; Я - длина волны основной гармоники.

В действительности, обычно величина Я неизвестна и в качестве основной длины волны выбирают произвольное значение Ь. Если выбрать Ь большим или равным длины изучаемого ряда данных и вычислить достаточное число гармоник, то получим оценку периодичности, присутствующей в наших данных. Произвольный выбор значения Ь скорее всего будет неправильной

i Надоели баннеры? Вы всегда можете отключить рекламу.