Донской государственный технический университет1, Научно-исследовательский институт многопроцессорных вычислительных систем им. академика А. В. Каляева Южного федерального университета2

Abstract

Работа посвящена разработке методов решения модельной задачи эвтрофикации вод мелководного водоема, учитывающей движение водного потока, микротурбулентную диффузию, гравитационное оседание, пространственно-неравномерное распределение температуры и солености, а также загрязняющих биогенных веществ, кислорода, фито- и зоопланктона и др. В качестве объекта моделирования были выбраны мелководные водоемы - Азовское море и Таганрогский залив. Численное решение задачи основывается на градиентном методе вариационного типа - методе минимальных поправок. Ускорение и эффективность расчетов достигается при использовании многопроцессорной вычислительной системы Южного федерального университета. Решение поставленной задачи водной экологии позволит прогнозировать изменения качества вод мелководных водоемов, а также изучать механизмы формирования в них зон с пониженным содержанием кислорода. Ключевые слова: математическая модель, эвтрофикация, метод минимальных поправок, параллельные вычисления, Азовское море.

-

Математическое моделирование процессов эвтрофикации в мелководных водоемах на многопроцессорной

вычислительной системе* 0 , дном H  H(x, y) и цилиндрической поверхностью  для 0  t  T0.   0  H  – кусочно-гладкая граница области G .

y 0 (1) где Si – концентрация i -й примеси, i  1,17 , индекс i указывает на вид субстанции, i  1,17 : 1 – сероводород H2S ; 2 – элементная сера S  ; 3 – сульфаты SO4  ; 4 – тиосульфаты (и сульфиты); 5 – общий органический азот N  ; 6 – аммоний NH4  (аммонийный азот); 7 – нитриты NO2 ; 8 – нитраты NO3  ; 9 – растворенный марганец DOMn ; 10 – взвешенный марганец POMn ; 11 – растворенный кислород O2 ; 12 – силикаты (SiO3 – метасиликат; SiO4 – ортосиликат); 13 – фосфаты PO4  ; 14 – железо Fe2  15 – кремнекислота ( H2SiO3 – метакремневая; H2SiO4 – ортокремневая); 16 – фитопланктон; 17 – зоопланктон; u  u,v,w  – скорость водного потока; wgi – скорость гравитационного осаждения i-й компоненты, если она находится во взвешенном состоянии; i,i – коэффициенты турбулентного обмена соответственно по горизонтальному и вертикальному направлениям;  i – химико-биологический источник (сток) или член, описывающий агрегирование (слипание-разлипание), если соответствующая компонента является взвесью.

К системе (1) добавим граничные условия: Si  0 на ,если Un  0; Si  0 на ,если Un  0;Si,z  0 на 0;Si,z  iSi на H, 2 n где i – коэффициент поглощения i -й компоненты донными отложениями.

Необходимо также добавить начальные условия:

Si t0  Si0(x, y,z), i  1,17.

С помощью модели ЭВ (1) – ( 3 ) могут быть описаны процессы аммонификации, нитрификации, нитратредукции (денитрификации), ассимиляции NH4 , окисления H2S , сульфатредукции, окисления и восстановления марганца, а также можно изучать механизм условий формирования заморов, прогнозировать изменения кислородного и биогенного режимов в мелководном водоеме. При построении модели были параметризованы процессы биогеохимических циклов химических элементов, ответственных за трансформацию аэробных условий в анаэробные [2].

Проведено исследование трехмерной модели ЭВ Азовского моря вида (1) – ( 3 ), получены достаточные условия существования и единственности ее решения, сформулированные в виде теорем.

Теорема 1. Пусть Si(x,y,z,t),  i  C2(Цt)  C(Цt), где Цt  G  (0  t  T0); i  const  0 ; U  u,v,w  wgi  , i(z)  C1(G); Sio  C(G), i  1,17 . Тогда при выполнении неравенств: max i,i  10 maGx  pi   0 для всех i  1,17 , где  i  pi Sj  Si   i G ( 3 ) ,i  j ; i  LSi  DSi  piSi ,

LSi  Si  div SiUi  , t

DSi  iSi  z i Szi  , 0   2 1 / L2x  1 / L2y  1 / L2z  ; Lx,Ly,Lz – пространственные максимальные размеры расчетной области; модельная задача эвтрофикации мелководного водоема ЭВ вида (1) – ( 3 ) имеет решение.

Теорема 2. Пусть Si(x,y,z,t),  i  C2(Цt)  C(Цt), i  const  0; U , i(z)  C1(G), i  1,17 . Тогда при выполнении неравенств: 2i 1 / L2x  1 / L2y   2vi / L2z  i для всех i  1,17 (где функции i определяются источниками загрязняющего вещества (ЗВ)) модельная задача ЭВ мелководного водоема вида (1) – ( 3 ) имеет единственное решение. Для дискретизации модели (1) – ( 3 ) построим в области решения задачи связную сетку h . hx,hy,hz – векторные параметры, характеризующие плотность расположения узлов: h  xj  jhx, yk  khy,zl  lhz; j  0,Nx;k  0,Ny;l  0,Nz;

Nxhx  Lx,Nyhy  Ly,Nzhz  Lz ,h  h   ,  tn  n,n  0,Nt , где i, j,k – индексы по направлениям x,y,z ; Nx,N y,Nz – количество узлов по координатным направлениям.

Расчетная область по пространственным направлениям x,y,z представляет собой объединение параллелепипедов. Дискретизируем модель (1) – ( 3 ) с помощью разностной схемы с весами:  2i  Snij1,k1,l  (1   )Snij,k1,l  2 Snij1,k,l  2(1   )Snij,k,l   i  Snij1,k1,l  (1   )Snij,k1,l   hy hy2  h1z  nij,k,l12  h1z  nij,k,l12  Snij1,k,l1  (1   )Snij,k,l1   Snij1,k,l  (1   )Snij,k,l  

hz   Snij1,k,l  (1   )Snij,k,l   Snij1,k,l1  (1   )Snij,k,l1    hz  i  0, hx  2i / un Ch ,hy  2i / vn Ch ,hz  2i / wn  wgni  Ch ,i  1,17. ( 6 ) Исследование погрешности аппроксимации разработанной схемы с весами вида ( 4 ) показало, что при   0.5 (схема Кранка-Николсон) она имеет наибольший порядок точности по временной и пространственным переменным: O  2  h 2  , где h  hx2  hy2  hz2 .

Входными данными предложенной выше модели эвтрофикации вод мелководного водоема вида (1) – ( 3 ) является поле вектора скорости водного потока, рассчитанное по моделям Сухинова А.И., Чистякова А.Е. [3–6]. 3. Метод решения сеточных уравнений Представим в виде операторного уравнения дискретную модель вида ( 4 ):

Au  f ( 7 ) с невырожденным оператором A , заданным в вещественном гильбертовом пространстве H . Рассмотрим неявную двухслойную итерационную схему вида

B yk1  yk  Ayk  f,k  0,1, (8)

 k с произвольным начальным приближением y0  H и невырожденным оператором B [7]. Любой двухслойный итерационный метод, построенный на основе схемы ( 7 ), характеризуется операторами A и B , энергетическим пространством HD , в котором доказывается сходимость метода, и набором итерационных параметров  k . Основным вопросом теории итерационных методов является вопрос об оптимальном выборе параметра  k [8–10].

В двухслойных итерационных методах вариационного типа для вычисления параметров  k не требуется никакой априорной информации об операторах схемы (8) (кроме условий общего вида A  A*  0, DB1A*  DB1A и т.д.), и построение этих методов основано на следующем принципе: если задано приближение yk , а yk1 находится по схеме (8), то итерационный параметр  k1 выбирается из условия минимума в HD нормы погрешности zk1  yk1  u , где u – решение уравнения ( 7 ). Последовательность yk , построенная по формуле (8), в которой параметры  k выбираются из указанного выше условия, является минимизирующей для квадратичного функционала вида I y  D y  u ,y  u . Этот функционал (в силу положительной определенности оператора D ) ограничен снизу, достигает минимума, равного нулю, при решении уравнения ( 7 ), т.е. при y  u . Выбор параметра  k1 из указанного условия обеспечивает локальную минимизацию функционала I y  при переходе от yk к yk1 , т.е. за один итерационный шаг. В случае явной схемы B  E  переход от yk к yk1 осуществляется по формуле yk1  yk   k1rk, rk  Ayk  f .

Для самосопряженного положительно определенного оператора A переход от yk к yk1 происходит по направлению rk , которое совпадает с направлением антиградиента для функционала A y  u ,y  u в точке yk . По направлению антиградиента происходит наибольшее убывание значения функционала. Параметр  k1 будем находить из условия минимума в HD нормы погрешности zk1  yk1  u [11].

Формула для вычисления итерационного параметра  k1 находилась в предположении, что оператор A невырожден. Выпишем сначала уравнение для погрешности: zk  yk  u , k  0,1, . Подставляя yk   k  u в схему (8), получим zk1  E   kB1A  zk , k  0,1, , z0  y0  u . Замена zk  D12xk позволяет перейти к уравнению, содержащему только один оператор: Поставленную выше задачу о выборе параметра  k1 можно сформулировать, используя равенство zk D  xk   D  норма в HD ,   норма в H  . Следующим образом: выбрать параметр  k1 из условия минимума нормы xk1 в пространстве H . Вычислим норму: xk1 2  E   k1C  xk, E   k1C  xk   xk 2  2 k1 Cxk,xk    k21 Cxk,Cxk   имеем Cxk,Cxk   0 и минимум нормы xk1 достигается при k21  1  Cxk, xk 2 / Cxk,Cxk  xk, xk  . Формула ( 11 ) определяет оптимальное значение итерационного параметра  k1 . Подставляя 1 xk  D 2zk в ( 11 ) получим:  k1  DB1Azk,zkxk  / DB1Azk,B1Azk  ,k  0,1, , ( 14 ) Учитывая, что Azk  Ayk  Au  Ayk  f  rk – невязка, а B1rk  k – поправка, формулу для параметра  k1 можно записать в следующем виде:  k1  Dk,zk  / Dk,k  ,k  0,1, , ( 15 ) а итерационная схема (8) – в виде явной формулы для вычисления yk1 :

yk 1  yk   k 1k, k  0,1,.... ( 16 ) Алгоритм, реализующий построенный метод, можно кратко описать следующим образом: 1) по заданному yk вычисляется невязка rk  Ayk  f; 2) решается уравнение для поправки Bk  rk; 3) по формуле ( 15 ) вычисляется параметр  k 1. 4) по формуле ( 16 ) находится новое приближение yk1.

Рассмотрим частные случаи двухслойных градиентных методов, которые мы будем использовать для решения модельной задачи эвтрофикации вод мелководного водоема. Каждый конкретный метод определяется выбором оператора D и имеет свою область применимости. Оператор D будем выбирать так, чтобы в формулу ( 15 ) для итерационного параметра  k1 входили только известные в процессе итераций величины.

Если оператор A самосопряжен и положительно определен в H , то для решения (8) можно использовать метод скорейшего спуска (МСС). Если оператор A несамосопряженный и невырожденный, а оператор B*A положительно определен, то можно использовать метод минимальных невязок (ММН) [12].

Пусть оператор A самосопряжен и положительно определен в H . Для метода скорейшего спуска (МСС) D  A . Оператор B должен быть положительно определен в H . Учитывая соотношения Azk  Axk  f  rk и A  A* из ( 15 ) получим формулу для итерационного параметра  k1 в неявном методе скорейшего спуска:  k1  rk,k  / Ak,Ak  ,k  0,1, . ( 17 ) Для случая явной двухслойной схемы ( 11 ) B  E получим k  B1rk  rk и формула для  k1 примет вид:

 k1  rk,rk  / Ark,Ark  ,k  0,1, . (18) При решении модельной задачи (1) – ( 3 ) методом скорейшего спуска по невязке параметр  k1 рассчитывается по формуле (18), а метод скорейшего спуска по поправке реализуется с помощью расчетной формулы ( 17 ). Метод минимальных поправок (ММП) можно применять для решения модельной задачи (1) – ( 3 ) в случае любого несамосопряженного невырожденного оператора A , кроме того требуется положительная определенность оператора B*A . Для метода минимальных невязок D  A*A . Формула для итерационного параметра  k1 в методе минимальных поправок имеет вид:  k1  (Ak,rk)/(Ak, Ak),k  0,1, . (19) В случае явной схемы (8) B  E  требуется положительная определенность оператора A, а формула для  k1 примет вид:

 k1  Ark,rk  / Ark,Ark  ,k  0,1, . (20) Был проведен анализ эффективности описанных выше градиентных методов вариационного типа на основе численного решения модельной задачи эвтрофикации вод мелководного водоема вида (1) – ( 3 ). В таблице 1 приведены результаты сравнения скоростей сходимости двухслойных градиентных методов вариационного типа, используемых для решения задачи эвтрофикации вод мелководного водоема вида (1) – ( 3 ): минимальных невязок (ММН), минимальных поправок (ММП), скорейшего спуска по невязке (МССН), скорейшего спуска по поправке (МССП). Таблица 1. Сравнение скоростей сходимости методов вариационного типа

3 IS  47,IX  43 IS  53,IX  46 IS  22,IX  19

4 IS  48, IX  42 IS  53,IX  45

IS  21,IX  19 IS  76,IX  63 IS  86,IX  72 IS  49,IX  43

IS  76,IX  62 IS  85,IX  71

IS  48,IX  40 IS  67,IX  60 IS  74,IX  67 IS  35,IX  29

IS  67,IX  59 IS  74,IX  66

IS  34,IX  28 IS  78,IX  66 IS  86,IX  73 IS  51,IX  42

IS  78,IX  65 IS  86,IX  72 IS  49,IX  41 p

n 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

1 IS  48,IX  44 IS  53,IX  48 IS  21,IX  21 IS  17,IX  15 IS  83,IX  65 IS  87,IX  73 IS  51,IX  51 IS  40,IX  38 IS  69,IX  62 IS  71,IX  68 IS  35,IX  32 IS  24,IX  23 IS  68,IX  64 IS  86,IX  74 IS  53,IS  49 IS  40,IX  38 ММП ММН 2 IS  47,IX  43 IS  52,IX  47 IS  22,IX  20 IS  14,IX  13 IS  77,IX  64 IS  87,IX  73 IS  50,IX  46 IS  38,IX  35 IS  67,IX  61 IS  74,IX  67 IS  35,IX  30 IX  21,IX  20 IS  79,IX  67 IS  87,IX  73 IS  52,IX  45 IS  36,IX  35 МССП МССН В таблице 1 используются следующие обозначения: p – номер итерации; n – номер временного слоя; I – количество итераций, необходимое для сходимости метода решения СЛАУ для уравнения, описывающего изменение концентрации  ,   S,X , где S  S5,X  S16 .

Метод минимальных поправок (ММП) был выбран в качестве основного метода в виду его наибольшей скорости сходимости, согласно данным, приведенным в таблице 1. 4. Параллельный вариант метода решения сеточных уравнений Для геометрического разбиения расчетной области с целью равномерной загрузки вычислителей (процессоров) МВС использовался метод k-means, основанный на минимизации функционала суммарной выборочной дисперсии разброса элементов (узлов расчетной сетки) относительно центра тяжести подобластей: Q  Q( 3 ). Пусть Xi – множество расчетных узлов сетки, входящих в i  ю подобласть, i  1,..., m , m  заданное количество подобластей. Q( 3 )   1  d2(x,ci)  min , где ci  1  x – центр подобласти Xi , а d(x,ci)– расстоi Xi xXi Xi xXi яние между расчетным узлом сетки x и центром подобластиci в Евклидовой метрике. Метод kmeans является сходящимся только тогда, когда все подобласти будут примерно равны.

Результат работы метода k-means для модельных двумерной и трехмерной областей представлен на рис. 2.

Рис. 2. Результат работы алгоритма k-means для разбиения модельной двумерной области на 9, 38, 150 подобластей; для трехмерной на 6 и 10 Опишем работу максиминного алгоритма. В качестве центров подобластей выбирает расчетные узлы сетки следующим образом: 1) первый центр – первый расчетный узел области; 2) второй центр находится в расчетном узле сетки, расположенном на максимальном расстоянии от первого центра; 3) если количество подобластей больше 3-х, то каждый следующий центр находится на максимальном удалении от ближайшего центра.

Опишем алгоритм работы k-means. 1) Выбираются начальные центры подобластей при помощи максиминного алгоритма. 2) Все расчетные узлы разбиваются на m клеток Вороного по методу ближайшего соседа, т. е. текущий расчетный узел сетки x  Xc , где Xc – подобласть выбирается из условия x  sc  min x  si , где sc – центр области Xc .

1im 3) Рассчитываются новые центры по формуле s(ck1)  X1(ik) xX(ik) x . 4) Проверяется условие остановки s(ck1)  s(ck) для всех k  1,..., m . Если условие остановки не выполняется, то осуществляется переход на пункт 2 алгоритма.

Для организации обмена данными в вычислительном процессе требуется найти все точки, лежащие на границе каждой подобласти. Для этой цели использовался алгоритм Джарвиса (задача построения выпуклой оболочки). Был сформирован список соседних подобластей для каждой подобласти и разработан алгоритм пересылки данных между подобластями [13,14]. При решении СЛАУ методом минимальных поправок для расчета итерационного параметра  использовался метод сдваивания.

Синхронизация алгоритма решения задачи (1) – ( 3 ) требуется только в ММП при переходе на следующую итерацию. 5. Описание программного комплекса

Для решения модельной задачи эвтрофикации мелководного водоема (1) – ( 3 ) был создан комплекс прикладных программ, позволяющий производить расчеты концентраций ЗВ, фито- и зоопланктона в областях сложной формы (Азовское море и Таганрогский залив) [15].

С помощью созданного комплекса программ можно осуществлять:  совершенствование и внедрение системы комплексного рыбохозяйственного мониторинга в водоемах (наблюдение, оценка и прогноз состояния режима экосистем, кормовой базы и запасов промысловых объектов);  разработку, согласование предложений и мероприятий по обеспечению оптимального режима, сохранения биоразнообразия промысловых ресурсов, экосистем мелководных водоемов;  совершенствование методологии природоохранных исследований, разработка новых, апробация и внедрение перспективных методов изучения состояния водных экосистем и отдельных компонентов;  разработку и совершенствование методов диагностики токсического воздействия ЗВ на гидробионты, в том числе ранней и дифференциальной диагностики токсикоза, а также поиск средств антидотной защиты водных экосистем;  организацию и проведение исследований по выявлению тенденций и закономерностей изменения состояния водных экосистем под воздействием антропогенных факторов, разработка предложений и мероприятий по снижению и предупреждению таких воздействий;  оценку ущербов рыбному хозяйству, наносимых разными видами хозяйственной деятельности, разработку предложений по предотвращению, уменьшению и адекватной компенсации ущербов.

Поля скоростей водного потока, рассчитанные в работах [16], относятся к входным данным для моделей гидробиологических процессов, описанных в первой главе. Для математического моделирования гидробиологических и гидродинамических процессов в трехмерной области сложной формы – Азовское море и Таганрогский залив использовались последовательно сгущающиеся прямоугольные сетки размерностями: 251  351  15, 502  702  30 , 1004  1404  60, .

Начальное распределение загрязняющих веществ и планктона было учтено в форме, соответствующей пространственно-временным масштабам моделируемых процессов. Реализованный алгоритм численного решения позволяет свободно варьировать граничные условия, вид управляющих функций и значения соответствующих параметров [17]. Понимание механизма функционирования системы и знание еe основных характеристик позволили для преодоления трудностей при настройке программы использовать феноменологический подход. Эффективность такого подхода достаточно высока ещe и потому, что поведение системы часто определяется точностью не отдельных параметров, а их соотношений.

Калибровка и верификация разработанной модели эвтрофикации вод мелководного водоема проводились на основе экологических данных по Азовскому морю, полученных в ходе научноисследовательских экспедиций, проводимых учеными ЮФУ, начиная с 2000 года. В ходе исследований экватории Азовского моря изучались: виды и концентрации основных загрязняющих воды Азовского моря веществ; пространственные распределения солености и температуры; кислородный режим; видовой состав фито- и зоопланктона; механизмы возникновения и развития заморов в центрально-восточной части водоема.

Обработка экспедиционных данных заключалась в оцифровке, пересчете в стандартные единицы, классификации для использования в разных модельных задачах гидробиологии моря. На рис. 4 приведена схема разработанного программного комплекса. Рис. 4. Схема программного комплекса Программный комплекс включает: блок управления, базы океанологических и метеорологических данных, системы интерфейсов, системы ввода – вывода и визуализации. Программный продукт обладает удобным интерфейсом, что позволяет вводить нужную информацию в диалоговом режиме. Построенный программный комплекс имеет универсальный характер и привязка его к условиям конкретных объектов и районов осуществляется, как правило, на уровне входной информации. Это значит, что для практического использования модулей требуется создание специальной информационной базы, содержащей сведения о физико-географических и климатических условиях исследуемых объектов, о параметрах, определяющих источники примесей.

При разработке программного комплекса использовался язык высокого уровня С++. Были разработаны две версии параллельного алгоритма решения модельной задачи эвтрофикации вод мелководного водоема. Первая предназначена для ЭВМ с операционной системой Windows, основана на технологии создания в этой системе дополнительных потоков. Вторая может использоваться для кластерных систем и основана на технологии передачи сообщений (MPI). 6. Результаты численных экспериментов

При моделировании пространственно-неоднородных процессов эвтрофикация вод Азовского моря (1) – ( 3 ) учитывалась внешняя периодичность, приводящая к усложнению системы. При этом колебания плотности планктона могут быть столь велики, что не могут быть объяснены случайными флуктуациями, и визуальная картина такова, что сравнительно небольшие площади высокой плотности ("пятна", "облака") разделены пространствами с низкими плотностями, иногда не фиксируемыми стандартными методами наблюдений. Особенно ярко это явление выражено в тех местах водоема, для которых характерна потребность в биогенных элементах. При моделировании процессов эвтрофикации учитывался вегетационный период фитопланктона.

Диффузные процессы действуют в направлении сглаживания пространственного распределения и рассеивания "пятен". Одна из попыток объяснения парадокса стабильности "пятен" с помощью численных экспериментов заключается в предположении об активном передвижении гетеротрофных организмов (зоопланктона и рыб) в направлении градиента "пищи", что обеспечивает закрепление пространственной неоднородности биогенных веществ в водной среде. Устойчивая неоднородность пространственного распределения может быть, например, обусловлена диффузионными процессами и наличием у фитопланктона механизма эктокринного регулирования, т.е. регулирования скорости роста посредством выделения в среду биологически активных метаболитов.

На рисунке 5. показаны результаты расчета концентрации загрязняющего биогенного вещества для модели динамики вредоносной водоросли (начальное распределение полей течений при северном ветре). Приведенные ниже рисунки отражают влияние структур течений водного потока в Азовском море на распределение загрязняющего биогенного вещества и фитопланктона.

Рис. 5. Распределение концентрации загрязняющих веществ в различные моменты времени Результаты моделирования динамики фитопланктона в Азовском море представлены на рис. 6. ( N – номер итерации, начальное распределение полей течений водного потока при северном направлении ветра).

Рис. 6. Распределение концентрации фитопланктона в различные моменты времени Белым цветом выделены максимальные значения концентраций биогенного вещества (азота) и фитопланктона черным – минимальные.

Критерием проверки адекватности построенных моделей гидробиологии мелководного водоема служила оценка погрешности моделирования с одновременным учетом натурных данных по имеющимся n замерам, которая вычислялась по формуле:   n Sk nat  Sk 2 n Sk2 nat , где Sk nat – значение концентрации загрязняющей приk1 k1 меси, полученное с помощью натурных экспедиционных измерений; Sk – значение, рассчитанное с помощью модели (1) – ( 3 ).

Рассчитанные при различных ветровых ситуациях концентрации загрязняющих веществ и планктона принимались к рассмотрению, если относительная погрешность не превышала 30%. 7. Заключение

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

Отличительными особенностями разрабатываемых алгоритмов, реализующих поставленные гидробиологические модельные задачи, являются: высокая производительность, достоверность и точность получаемых результатов. Высокая производительность достигается за счет использования эффективных численных методов решения сеточных уравнений, ориентированных для применения на параллельных вычислительных системах в реальном и ускоренном масштабах времени. Достоверность достигается за счет учета определяющих физических факторов, таких как: сила Кориолиса, турбулентный обмен, сложная геометрия дна и береговой линии, испарение, стоки рек, динамическое перестроение расчетной области, ветровые напряжения и трение о дно, а также за счет учета отклонения значения поля давления от гидростатического приближения. Точность достигается применением подробных расчетных сеток, учитывающих степень «заполненности» расчетных ячеек, а также отсутствием неконсервативных диссипативных слагаемых и нефизичных источников (стоков), возникающих в результате конечно-разностных аппроксимаций.

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

Анализ показал, что в результате разработки программного комплекса удалось повысить точность прогнозов изменения концентраций загрязняющих веществ, фито- и зоопланктона в мелководном водоеме на 10 – 15% в зависимости от решаемой модельной задачи водной экологии. Литература 1. Якушев Е.В., Сухинов А.И. и др. Комплексные океанологические исследования Азовского моря в 28-м рейсе научно-исследовательского судна «Акванавт» // Океанология. 2003. Т. 43, №1. С. 44–53. Mathematical modeling of eutrophication processes in shallow waters on multiprocessor computer system A.I. Sukhinov1, A.V. Nikitina2, A.E. Chistyakov2, I.S. Semenov2, A.A. Semenyakina2,

D.S. Khachunts2 The paper covered the development of solution methods for model of eutrophication of shallow water in view the movement of the water flow, microturbulent diffusion, gravitational settling, spatially uneven distribution of temperature and salinity, and pollutants of nutrients, oxygen, phyto- and zooplankton, etc. Shallow waters of the Azov Sea and Taganrog Bay were selected as the object of simulation. Numerical solution of the problem is based on gradient method of variation type – the method of minimal corrections. The acceleration and efficiency of calculations is achieved with using multiprocessor computer system of Southern Federal University. The solution of the problem of water ecology will allow to predict changes of water quality of shallow reservoirs, and to study the mechanisms of formation of zones with low oxygen content.

References

Yakushev E.V. , Sukhinov A.I. Complex Oceanographic research of the Sea of Azov in the 28-th flight of the research vessel "Aquanaut" // Oceanology. 2003 . Vol. 43 , No. 1. P. 44 - 53 .

2. Sukhinov A.I. , Nikitina A.V. , Chistyakov A.E. Simulation scenario of biological rehabilitation of the Azov Sea / / Mathematical modeling. 2012 . Vol. 24 , No. 9. P. 3- 21 .

3. Sukhinov A.I. , Chistyakov A.E. , Alekseenko E.V. Numerical implementation three-dimensional model of hydrodynamics for shallow water basins on superficialities system // Mathematical modeling . 2011 . Vol. 23 , No. 3. P. 3- 21 .

4. Sukhinov A.I. , Chistyakov A.E. Parallel implementation of a three-dimensional model of hydrodynamics of shallow water bodies on superficialities system // Computational methods and programming: New computing technologies . 2012 . T. 13. P. 290 - 297 .

5. Belotserkovsky O.M. Turbulence: new approaches . M.: Science , 2003 . 286 p.

6. Samarskii A.A. Theory of difference schemes . M.: Science , 1989 . 616 p.

7. Samarskii A.A. , Nikolaev E.S. Methods of solving grid equations . M.: Science , 1978 . 592 p.

Konovalov A.N. The theory of the alternating-triangular iterative method // Siberian mathematical journal . 2002 . 43 :3. P. 552 - 572 .

9. Sukhinov A.I. , Chistyakov A.E. Adaptive modified alternating triangular iterative method for solving grid equations with non-selfadjoint operator / / Mathematical modeling. 2012 . Vol. 24 , No. 1. P. 3- 20 .

10. Chistyakov A.E. Theoretical estimation of speed up and efficiency of parallel implementation of the steepest descent PTM // Izvestiya SFEDU . Technical Sciences . 2010 . No. 6 ( 107 ). P. 237 - 249 .

11. Chistyakov A.E. , Hachunts D.S. , Nikitina A.V. , Protsenko E.A. , Kuznetsova I.Yu. A library of parallel iterative methods solvers of linear algebraic equation for the problem of convection-diffusion-based decomposition in one spatial direction // Modern problems of science and education . 2015. No. 1 . URL: http://www.science-education. ru /121-19510 (accessed: 4 . 06 . 2015 ).

12. Nikitina A.V. Numerical solution of problems of dynamics of toxic algae in the Taganrog Bay // Izvestiya SFedU . Engineering sciences. 2010 . No. 6 ( 107 ). P. 113 - 116 .

13. Nikitina A.V. modelling of biological kinetics, stabilizing the ecological system of the Gulf of Taganrog // Izvestiya SFedU . Engineering sciences. 2009 . No. 8 ( 97 ). P. 130 - 134 .

14. Sukhinov A.I. , Chistyakov A.E. , Bondarenko Yu .S. Error estimate of the solution of the diffusion equation on the basis of the schemes with weights // Izvestiya SFedU . Engineering sciences. 2011 . No. 8 ( 121 ). P. 6- 13 .

15. Sukhinov A.I. , Chistyakov A.E. , Semenikhina A.A. , Nikitina A.V. Parallel realization of the tasks of the transport of substances and the recovery of the bottom surface on the basis of the schemes of high order of accuracy // Computational methods and programming: new computing technologies . 2015 . T. 16. P. 256 - 267 .

16. Nikitina A.V. , Rudneva T.V. , Kamyshnikova T.V. , Bokareva T.A. , Duryagina V.V. To the formation of kill zones in the Eastern part of Azov Sea coast // Modern problems of science and education . 2015. No. 1 . URL: http://www.science-education. ru/121-19509 (accessed: 4 . 06 . 2015 ).

17. Nikitina A.V. , Puchkin M.V. , Semenov I.S. , Sukhinov A.I. , Ougolnitsky G.A. , Usov A.B. , Chistyakov A.E. Differential-game model of prevention of fish kill in shallow ponds // Managing large systems , Issue 55. M.: IPU Russian Academy of Sciences . 2015 . P. 343 - 361 .