CC BY
9
Известия Института математики и информатики Удмуртского государственного университета
2022. Том 60. С. 73-89
УДК 519.6
© А. И. Сухинов, А. Е. Чистяков, А. М. Атаян, И. Ю. Кузнецова, В. Н. Литвинов, А. В. Никитина
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОЦЕССА ОСАЖДЕНИЯ НА ДНО МНОГОКОМПОНЕНТНОЙ ВЗВЕСИ И ИЗМЕНЕНИЯ СОСТАВА ДОННЫХ МАТЕРИАЛОВ
В работе рассмотрены 2Э и 3Э модели транспорта взвешенных частиц, учитывающие следующие факторы: движение водной среды; переменную плотность, зависящую от концентрации взвеси; мно-гокомпонентность взвеси; изменение геометрии дна в результате осаждения взвеси. Аппроксимация трехмерного уравнения диффузии-конвекции выполнена на основе схем расщепления на двумерную и одномерную задачи. В работе используются дискретные аналоги операторов конвективного и диффузионного переносов в случае частичной заполненности расчетных ячеек. На основе функции заполненности описывается геометрия расчетной области. Использована схема, представляющая собой линейную комбинацию разностных схем «крест» и «кабаре» с весовыми коэффициентами, полученными в результате минимизации погрешности аппроксимации. Данная схема предназначена для решения задачи переноса примеси при больших сеточных числах Пекле. Приведены результаты численных экспериментов, из которых сделаны выводы о преимуществе 3Э модели транспорта многокомпонентной взвеси по сравнению с 2Э моделью. Выполнены численные эксперименты по моделированию процесса осаждения многокомпонентной взвеси, изучено его влияние на рельеф дна и изменение его состава.
Ключевые слова: модель транспорта взвеси, переменная плотность, схема «кабаре», схема «крест», изменение рельефа дна, параллельные алгоритмы.
001: 10.35634/2226-3594-2022-60-05 Введение
В последние годы участились явления природного и антропогенного характеров, вызывающих негативные экономические и тяжелые экологические последствия. Катастрофический шторм в Керченском проливе в ноябре 2007 г. привел к крушению более чем 20 судов. Разливы нефти привели к загрязнению береговой линии и донных отложений нефтепродуктами и другими вредными и опасными веществами. Соединения нефтепродуктов в виде битумов и смол были обнаружены на побережье Черного и Азовского морей, протяженностью более 200 км в 2008-2011 гг., что явилось одним из последствий этой катастрофы. 24-25 сентября 2014 г. произошло затопление прибрежных районов Азовского моря вследствие штормового нагона, уровень воды в районе порта г. Таганрога поднялся более чем на 4 метра, что нанесло значительный урон экономике региона. В конце ноября 2019 года из-за сильного восточного ветра над Азовским морем (10-15 м/с с порывами до 20 м/с) произошел катастрофический угон вод Таганрогского залива, при этом уровень воды в реке Дон снизился до критической отметки (-110 см). В связи с обмелением реки, приустьевой зоны моря и судоходных каналов было полностью прекращено судоходство. Кроме того, сильный ветер поднимал песок с обмелевшего дна моря в воздух и переносил его на значительные расстояния. Пыльная буря накрыла приморские города и поселения, нанося вред здоровью населения, проживающего в прибрежной зоне.
Для предупреждения возникновения и уменьшения последствий опасных и катастрофических явлений необходимо уметь выявлять критические состояния природных систем,
при которых возможно появление чрезвычайных ситуаций, строить оперативные и научно оправданные прогнозы и своевременно передавать полученные результаты лицам, которые могут принимать решения [1,2]. Для построения прогнозов, обладающих предсказательной ценностью, необходимо создать относительно недорогой, точный и быстро работающий инструментарий, учитывающий информацию, поступающую от административно и технически различных систем мониторинга. Инструментарием, оптимальным по сравнению с натурными исследованиями, является математическое моделирование гидродинамических и гидробиологических процессов [3-6]. Необходимы современные средства, позволяющие повысить точность моделирования: взаимосвязанные нелинейные 3Б модели гидродинамики и гидробиологии, гранично-адаптивные призматические сетки, разностные схемы повышенного порядка точности, верификация и валидация численных алгоритмов их реализующих.
Процессы гидрофизики водоемов могут быть описаны в виде сложных систем начально-краевых задач с уравнениями в частных производных, включающих уравнения гидродинамики, транспорта тепла и солей, взвешенного вещества. В результате полученные математические модели являются трехмерными, нестационарными и существенно нелинейными. Для решения подобного класса задач зачастую используют методы реализации уравнения диффузии-конвекции. Кроме того, актуальной задачей является разработка разностных схем для решения задач гидродинамики в случае преобладания конвективного переноса над диффузионным, что характерно при возникновении некоторых опасных явлений природного характера, в том числе штормовых нагонов, вызывающих взмучивание и вторичное загрязнение вод изучаемого водоема. В связи с ограниченным временем построения прогнозов изменения экологического состояния водоема и требуемой точностью является актуальным проведение расчетов с использованием высокопроизводительных систем, что влечет за собой необходимость разработки параллельных алгоритмов, ориентированных на супер-ЭВМ [7-9].
§ 1. Постановка задачи
Для описания транспорта многокомпонентных взвесей воспользуемся системой уравнений диффузии-конвекции, которая может быть записана в следующем виде [10,11]:
где сг — концентрация г-ой фракции взвеси [мг/л], V = — составляющие поля
вектора скорости течения водного потока [м/с], ,ш3>г — скорость осаждения г-ой фракции взвеси [м/с], ^ и V — горизонтальный и вертикальный коэффициенты турбулентной диффузии, соотвественно [м2/с], ^ — функция, описывающая интенсивность распределения источников г-ой фракции взвеси [мг/л-с].
На свободной поверхности поток в вертикальном направлении равен нулю, то есть (дсг¡дг) = 0. Вблизи поверхности дна V (дсг¡дг) = — гш3уГсг. На боковой поверхности отсутствует поток (дсг¡дп) = 0, если (V, п) ^ 0 и взвесь выходит за границу расчетной области V (дсг¡дп) = (V, п) сг, если (V, п) < 0, где п — вектор нормали, направленный внутрь расчетной области.
Плотность водной среды рассчитывается согласно формуле:
где Уг — объемная доля г-ой фракции взвеси, р0 — плотность пресной воды при нормальных условиях, ру,г — плотность г-ой фракции взвеси, Я — количество фракций взвеси.
у
(1.2)
Уравнение транспорта наносов в случае многокомпонентной взвеси запишется в виде [12]:
(1-ё)Щ + &у ( УУЛть) = сИу ( V Я ) • V ' "Ч-,.. (1.3)
) 17=1 81П ^ ) 7=1 Р"
где Я — глубина водоема, е — усредненная по фракциям пористость донных отложений, ть — вектор касательного тангенциального напряжения на дне водоема, тЬс г — критическое значение тангенциального напряжения для г-ой фракции взвеси, тЬсг = агв1п <^0, аг — коэффициент для г-ой фракции взвеси, — угол естественного откоса грунта в водоеме, кг = кг (Я, х, у, г) — нелинейный коэффициент, определяемый соотношением:
_ Аи)йг
ГЬг ==
((Рг - Р0) д4)
в
Ч--grad Я
в1П <^0
в-1
, (1.4)
где рг и Аг — плотность и характерный размер частиц г-ой фракции взвеси, соответственно, ш — усредненная частота волн, А и в — безразмерные постоянные.
Уравнение (1.3) рассматривается при следующем начальном и граничном условиях:
Я (х, у, 0) = Я0 (х, у) и ЯП = 0.
§ 2. Аппроксимация задачи транспорта взвешенных частиц
Рассмотрим аппроксимацию трехмерного уравнения диффузии-конвекции
4 + (ис)х + МУ + М! = (^4)1 + (^4) у + )!. (2.1)
Построим равномерную сетку шт по времени с шагом г:
шт = = пт, п = 0, Л^, Л^т = Т} .
Для уравнения (2.1) воспользуемся схемами расщепления на двумерную и одномерную задачи:
„га+1/2 — сп
' и (<?)',+ ь(<?)'у=(1л(<гуя!у1в+(»(<?,у
+ и (сп)1 + V (сп)У = (сп)У 1 + (сп, (2.2)
„га+1/2 '
где сп, сп+1/2 и сп+1 — значения концентрации однокомпонентной взвеси на текущем временном слое , промежуточном временной слое ¿п+1/2 и на последующем временном слое ¿п+1, соответственно.
Для численной реализации дискретной математической модели поставленной задачи (2.2) вводится пространственная сетка [13]:
"7, = {= '//.,.. Уз = ;)ЬУ\ /' = 0, МХ1 ) = 0, Ку; Ыгт = Г. \,1> г = I, - ^уку = I,,} ,
где т — шаг по времени, Нх, Ну — шаги по пространству, N — верхняя граница по времени, Ыу — границы по пространству, 1х, 1У — характерные размеры расчетной области. Для аппроксимации однородного уравнения (2.2) будем использовать схемы расщепления по пространственным координатным направлениям:
„п+1/2 - „
- + ?/. (г'1) , ,
х У^ >х) х
т
„п+1 — сп+1/2
+ и (сп)х = (сп)^х , (2.4)
+ V (сп+1/2)У = (сп+1/2)^У , (2.5)
Введем коэффициенты q0, q1, , q4, описывающие заполненность областей, находящихся в окрестности ячейки. Значение q0 характеризует заполненность области
{х е (»¿-1/2, »¿+1/2), У е (У]-1/2, У]+1/2)}, ql - А : До П{х ^ }, q2 - Аз : А П{х ^ ^}, ф - £>3 : А п {у ^ у,}, - £>4 : А) П {у ^ у,-}.
Заполненные части областей Ат будем обозначать Пт, где т = 0,4. Согласно [13] коэффициенты qm можно вычислить по формулам:
/ ч _ , ч _ + 4+1,] + ^¿+1, ]+1 + 0¿,j+1 _ ^¿+1,] + ^¿+1, ]+1
- ^Г^' — -4-' - -9-'
ч _ Ог,з + Ч, з+1 , ч _ 4+1,,;+1 + 4,7+1 / ч _ 4,7 + 4+1, 7
~~ 9 ' — 9 ' — 9 '
где o¿,j• — степень заполненности ячейки (г,]).
Дискретные аналоги операторов конвективного ucX и диффузионного переносов случае частичной заполненности ячеек и граничных условий Сп = апс + вп могут быть записаны в следующем виде [15]:
! Ч ' ! Ч 4+1,] 4,] , / Ч 4,] C¿— 1,]
Ыг,,"ИСх ~ Шг,Щ+1/2,з-ТТ- + Ыг,,"И»-1/2
Ы*," )'г - (VI I 2,/ 'г+112 ^ - ('/•->>:..,•/': I 2
j 2Ь у ' ¿—1/2>j 2Ь
хх
4+1,]' 4,]' / ч C¿,j 4— 1,]
Чи ¿,]/^+1/2,] ,2 ^^ 1/2,] 12
' х ' х
/ \ / \ + Дх
кх
Для аппроксимации модели (2.1) будем использовать схему, полученную как результат линейной комбинации схем «кабаре» и «крест», при этом будем учитывать функцию заполненности ячеек [16]:
(а) разностная схема для уравнения (2.4), описывающего перенос вдоль направления Ох, запишется в виде:
о "+1/2 _ ^ ^ _ ^
Щ2л,з + чо,м 4,7 4,7 р. 4,7 4-1,3 ,
+ I 2. / /2,.;. /-^Т--Г
3 Т —ь-ч^^ 3,
^ _ ^
/ \ C¿+1 7 C¿ 7
+ «¿+1/2,7 111111 Щ1Л,з-,Ч2,г,з) -¡7Г-- +
х
7------V,! з + з
А" _ А" А" _ А"
'^¿+1,] Ч,] 0 Ч,]' c¿— 1,]
2^+1/2,^91,-Т^-:--^1^г-1/2,зЯ.2Л,3
кх -^-ч^™ кх
I I ахЧ,] + вх
- |<?1,М - 0_2Л,з\ -^-> «М ^ °>
х
где Ах= ]2 _ ^ ^ /т.
(Ь) разностная схема для уравнения (2.5), описывающего перенос вдоль направления Оу, запишется в виде:
2 + q ™+1 _ П+1/2 П+1/2 _ П+1/2
¿Ч4Л,3 + Чол,з 4,3 4,3 К Ч7_4,3-1
—з--;— + 1 —шу— +
cn+1/2 cn+1/2 2А га+1/2 + А cn+1/2q
4,3 + 1 ~ 4,1 , У Ч,7 — 1 <?4,м + ^У4,7 1о,'
3//,.
I • / Ч ¿,] + 1 ¿,] I
+ 4,^+1/2 П1111 {дзл^, д^)-5Т--ь
Чу
cn+1/2 _ «+1/2 cn+1/2 _ cn+1/2
c¿,j+1 c¿,j п,, „ C¿,j C¿,j —1
Г) ¿,] + 1 ¿,] Г)
*Иг,3+1 /21зЛ,з-^--1 /2<1АЛ,3 '
,2
уу
"+1/2 + в
I | аУC¿,j + вУ ^ А
- |<?3,м - -7-, ^ 0,
ку
где Лу<Г/2 = - <71/2) /т-
Для того, чтобы получить разностные схемы, аппроксимирующие уравнения (2.4), (2.5) при < 0 и Ьг^ < 0 из представленных аппроксимаций, необходимо соответствующие координатные оси Ох и Оу направить в противоположные стороны. Уравнение (2.3) решается методом прогонки.
§ 3. Расчет изменения состава донных материалов
Для работы алгоритма необходимо сформировать следующие массивы: Дг,3п. — массив, описывающий уровень заполненности ячейки (г^п^-); (Д)г 3 п, — массивы, описывающие концентрации /-ой фракции взвеси в ячейке (г,пг,3); пг,3 — массив, описывающий количество уровней заполненных ячеек, расположенных ниже данной.
Алгоритм расчета изменения состава донных материалов.
1. Рассчитывается толщина слоя взвеси, выпавшей в осадок, для каждой фракции согласно формуле:
)г,3 = Т (°1 )г,3,к ^,г/Р1,
где р1 — плотность /-ой фракции взвеси, узел (г,^, к) находится на дне расчетной области.
2. Рассчитывается толщина слоя взвеси, выпавшей в осадок, для всех фракций:
^ ,3 = X] )г ,3 • I
3. Если толщина слоя больше нуля > 0), то рассчитывается концентрация каждой
фракции в слое выпавшей взвеси:
^,3 = №)г,3 •
4. Выпавшая взвесь размещается между ячейками на основе следующего алгоритма. Пока толщина слоя больше нуля > 0), выполняются следующие действия:
4.1. Рассчитывается уровень слоя осадка в ячейке (г,^,пг)3) после осаждения взвеси:
4.2. Рассчитывается уровень слоя осадка, который была способна принять ячейка
гг,3 hz hz •
4.3. Если ячейка (г, пг,3-) после осаждения взвеси переполняется (^,3 > ^), то выполняется расчет концентраций фракций взвеси и переход к заполнению следующей ячейки.
4.3.1. Пересчитываются концентрации для каждой фракции в ячейке (г,^п^):
(D¿)г)J)ni,J = (^Мг,," Ч] + (В1)г,3,щ,] ^.п;,;^•
4.3.2. Ячейка (¿,^,пг,3) помечается как заполненная: (Д)г3-п = 1.
4.3.3. Осуществляется переход к заполнению следующей ячейки пг , 3 = пг , 3 + 1.
4.4. Если слой взвеси помещается в ячейку (sjj ^ то выполняется пересчет концентраций фракций взвеси описанным ниже способом:
4.4.1. Пересчитываются концентрации каждой фракции в ячейке (i, j, njj):
4.4.2. Пересчитывается уровень заполненности ячейки (i, j, n^j):
(D1)i)j,ni,J = /hz.
4.5. Уменьшается часть слоя взвеси, не размещенная между ячейками:
dhij — dhij rij *
§ 4. Параллельная реализация модели транспорта многокомпонентной взвеси
В рамках данной работы построен параллельный алгоритм, реализующий поставленную трехмерную задачу переноса веществ (1.1) на многопроцессорной вычислительной системы Научно-технологического университета «Сириус» под управлением системы MPI. Кластер включает в себя 24 процессора (вычислительных узла) семейства Intel(R) — Xeon(R) Gold 5118 CPU 2.30GHz. Тип накопителя каждого компьютера — SSD объемом 3.84 TB. Выполнено сравнение результатов работы построенного параллельного алгоритма для разного числа задействованных процессоров. Данные о расчетных сетках и времени работы программных компонент в секундах представлены в Таблице 1.
Таблица 1. Результаты работы параллельного алгоритма
Расчетная сетка 1 2 4 8 12 16 20 24
200x200x100 1.3167 0.7747 0.3401 0.2053 0.1392 0.1352 0.1118 0.0974
1000x1000x100 55.303 30.031 16.640 6.890 4.601 4.421 3.268 2.900
Для повышения эффективности параллельного алгоритма выполнена декомпозиция расчетной области по двум пространственным направлениям [17,18]. Использован способ разбиения на р1 прямоугольников вдоль одного направления и р2 — вдоль другого [19]. Данный способ декомпозиции позволяет уменьшить объем передаваемых данных. При использовании декомпозиции по одному направлению объем пересылок равен 2 х р х х N, где р — число задействованных процессоров [20]. В случае декомпозиции по двум направлениям объем пересылок равен 2 х (р1 х + р2 х N) х N, где и N — количество расчетных узлов вдоль направлений осей Ох, Оу и О, соответственно.
Параллельный алгоритм расчета транспорта многокомпонентной взвеси с декомпозицией расчетной области по двум пространственным направлениям.
1. Вызов функций, которые возвращают р — число процессов внутри коммуникатора и т — номер текущего процессора в частной группе коммуникатора.
2. Задаются значения р1 и р2 (условно р1 — количество процессоров вдоль направления Ох, р2 — вдоль направления Оу).
3. Рассчитываются значения номеров процессоров вдоль каждого из направлений т1 и т2:
т
т\ = т пхшр1, т\ = —
р1
4. Рассчитываются Ж11 и Ж22 — номера начальных элементов и N и N — размеры фрагмента расчетной области вдоль осей Ох и Оу:
N11 =
N
22
тг (Жг - 2)
Р1
т2 (Я, - 2) Р2
N1
N2
(Ш1 + 1) (Ях _ 2)
_ N11 + 2
_ N22 + 2.
5. Выполняется декомпозиция расчетной области.
6. Начало цикла по времени.
7. Расчет элементов массива (соответствует последовательной программе).
8. Обмен данных вдоль направления Ох.
8.1. Передача/прием данных в сторону уменьшения номера процессора:
• если т1 > 0, то осуществляется передача данных (передаются элементы с индексом г = 1) процессору с номером р _ 1;
• если т1 < р1 _ 1, то осуществляется прием данных (принимаются все элементы с индексом г = N _ 1) от процессора с номером р +1.
8.2. Передача/прием данных в сторону возрастания номера процессора:
• если т1 < р1 _ 1, то осуществляется передача данных (передаются элементы с индексом г = N _ 2) процессору с номером р +1;
• если т1 > 0, то осуществляется прием данных (принимаются все элементы с индексом г = 0) от процессора с номером р 1.
Объем пересылаемых данных равен (^ _ 2) х (Жг _ 2).
9. Обмен данных вдоль направления Оу (аналогичным образом).
9.1. Передача/прием данных в сторону уменьшения номера процессора:
• если т2 > 0, то осуществляется передача данных (передаются элементы с индексом ] = 1) процессору с номером р _ р1;
• если т2 < р2 _ 1, то осуществляется прием данных (принимаются все элементы с индексом ] = ^ _ 1) от процессора с номером р + р1.
9.2. Передача/прием данных в сторону возрастания номера процессора:
• если т2 < р2 _ 1, то осуществляется передача данных (передаются элементы с индексом ] = ^ _ 2) процессору с номером р + р1;
• если т2 > 0, то осуществляется прием данных (принимаются все элементы с индексом ] = 0) от процессора с номером р _ р1.
Объем пересылаемых данных равен (^ _ 2) х (N2 _ 2).
10. Конец цикла по времени.
11. Вывод данных.
В таблице 2 представлены результаты работы описанного параллельного алгоритма для различного числа процессоров и различной декомпозиции расчетной области по двум пространственным направлениям.
Таблица 2. Результаты работы параллельного алгоритма с декомпозицией расчетной области по двум пространственным направлениям
200 х 200 х 100 1000 х 1000 х 100
р Р\ Р'2 время, с ускорение эффективность время, с ускорение эффективность
16 16 1 0.1352 9.74 0.6757 4.421 12.51 0.7818
16 8 2 0.1232 10.69 0.6854 3.52 15.71 0.9821
16 4 4 0.1123 11.72 0.7881 3.499 15.8 0.9878
20 20 1 0.1118 11.78 0.5891 3.268 16.92 0.846
20 10 2 0.1104 11.93 0.5965 3.203 17.27 0.8633
20 5 4 0.0969 13.6 0.6798 3.153 17.54 0.877
24 24 1 0.0974 13.51 0.5631 2.9 19.07 0.7946
24 12 2 0.0913 14.42 0.6007 2.846 19.43 0.8096
24 8 3 0.0837 15.74 0.6559 2.831 19.53 0.8139
24 6 4 0.0741 17.78 0.7407 2.77 19.97 0.8319
На рис. 1 представлены графики ускорения разработанного параллельного алгоритма в зависимости от числа задействованных вычислителей и размера расчетной сетки. Максимальное число использованных вычислителей — 24, максимальный размер расчетной сетки составил 1000 х 1000 х 100 расчетных узлов. Дополнительно на рис. 1 построен график линейного ускорения, которое представляет собой эталонное ускорение (эффективность равна 100%). В точках пересечения линейного ускорения и графика ускорения параллельного алгоритма эффективность вычислений приближается к 100%, а в некоторых точках превышает данную отметку, что свидетельствует о наличии суперлинейного ускорения разработанного параллельного алгоритма. Из рис. 1 видно, что построенный алгоритм на малых сетках (размером до 200 х 200 х 100 расчетных узлов) показал более низкую эффективность по сравнению с результатами, полученными на больших расчет-
24
Сетка 200 х 200 х 100
Сетка 1000 х 1000 х 100
Линейное ai
18
12
о
s
и s в ¿г г / о \ * g * 1 ' о 0 < , °8° 1
X О о 1 Л 0 ° Г« 0 , О о ) о о о „ о > о ° 0 } о
j/o
о
12
18
24
Количество процессоров
Рис. 1. Зависимость ускорения от числа вычислителей
ных сетках (размером от 1000 х 1000 х 100 расчетных узлов). Разработанный на основе технологии MPI параллельный алгоритм, используемый для решения трехмерных задач
диффузии-конвекции, на больших расчетных сетках показал высокую эффективность, которая составила более 75 %.
§ 5. Численное решение задач транспорта взвеси
Для изучения процессов оседания многокомпонентной взвеси разработан программный комплекс, включающий 2Б и 3Б математические модели и их численную реализацию. На основе 3Б модели (1.1) проведен численный эксперимент для расчетной области с параметрами: длина — 1 км; ширина — 720 м; шаг по горизонтальной пространственной координате — 10 м; шаг по вертикальной пространственной координате —1м; расчетный интервал — 2 ч. На рис. 2 представлена концентрация взвешенных частиц (г/л). Приведены значения поля концентрация взвеси в сечении расчетной области плоскостью, проходящей через точку выгрузки грунта, образованную векторами направленными: вертикально и вдоль течения. Течения водного потока направлены слева направо. Средняя эффективная скорость седиментации взвесей в районе проведения дноуглубительных работ задавалась следующим образом: = 2.042 мм/с. Среднее расстояние от точки выгрузки до дна водоема в районе проведения дноуглубительных работ составило 5.5 м. Скорость течения на глубинах от 4 до 10 м составила 0.075 м/с.
0
10
207.224
155.418
103.612
51.806
200 400 600 800
Рис. 2. Концентрация взвешенных частиц
Среднее время нахождения взвеси во взвешенном состоянии (время осаждения) составляет примерно 2 693 с (примерно 45 мин). Исходя из входных данных модели, можно заключить, что в период осаждения за счет конвективного переноса взвесь распределится примерно на 202 м вдоль течения (что составляет около 20 шагов сетки). Свой вклад в распространение взвеси внесет и диффузионный перенос, но разброс взвеси за счет диффузионного переноса не превысит 200 м. Из рис. 2 видно, что взвесь осаждается на расстоянии 100-300 м от места выгрузки, что согласуется с ожидаемым результатом. На рис. 3 приведена концентрация взвешенных частиц (г/л), усредненная по глубине водоема.
При моделировании пространственно-временной структуры и изменчивости полей концентраций взвесей на основе 2Б модели не учитывается распределение концентрации взвеси по вертикальной компоненте. В зонах проведения дноуглубительных работ максимальная концентрация может составлять более 50г/л, а максимально допустимая концентрация не должна превышать 50 мг/л, следовательно применяемая модель должна обеспечивать
5
0
0
72.855
54.641
400
36.428
200
18.214
0
0
0
200
400
600
800
Рис. 3. Концентрация взвешенных частиц, усредненная по глубине водоема
точность расчетов более 0.1%. При использовании двумерных моделей взвесь осаждается значительно медленнее. Например, предположим, что взвесь равномерно распределена по глубине. Возьмем шаг по времени, равный половине времени осаждения. Согласно модели, за первый шаг осаждается 1/2 взвеси, за второй — 1/4 взвеси. После двух шагов по времени во взвешенном состоянии не должно оставаться частиц, но согласно 2Б модели 1/4 взвеси еще не осаждено. При малом значении шага по времени т и расчете на время осаждения модель будет показывать, что в водоеме во взвешенном состоянии остается 1/е или 37 % частиц. Для того чтобы концентрация взвеси упала от С1 до с2, где с2 значительно меньше с1, модель покажет, что для осаждения взвеси затрачивается больше времени при малых шагах т.
Сравнение результатов экспериментальных расчетов на основе 2Б и 3Б моделей проведем на примере описанной выше модельной задачи дампинга грунта. На рис. 4 представлены сравнительные графики распределения концентрации взвешенных частиц через 2 часа после выгрузки, полученные на основе 2Б и 3Б моделей для различных скоростей осаждения взвеси и различного состава фракций. По вертикальной оси показана концентрация взвешенных частиц (для наглядности используется логарифмическая шкала). Горизонтальная ось проходит через зону проведения дноуглубительных работ и направлена вдоль течения.
Полученные экспериментальные результаты хорошо согласуются с приведенными выше теоретическими расчетами. Из рис. 4 видно, что вблизи зоны проведения дноуглубительных работ 2Б модель показывает более высокую концентрацию взвешенных частиц. На расстоянии 50 метров от района проведения работ 2Б и 3Б модели показывают идентичные результаты. При расчетах распространения взвеси на расстояние 120 метров результаты расчетов на основе 2Б и 3Б моделей отличаются уже на порядок. На основе расчетов при помощи 3Б модели по сравнению расчетами, выполненными на основе 2Б модели, видно, что взвесь осаждается значительно быстрее.
На рис. 5 представлены графики изменения гранулометрического состава дна при различных начальных концентрациях взвешенных частиц.
При моделировании транспорта взвешенных частиц рассматривалось осаждение двух фракций. Скорость осаждения фракции А равна 2.4 мм/с, процентное содержание фракции А в пылеватых частицах равно 36%. Скорость осаждения фракции Б — 1.775 мм/с, процентное содержание фракции Б — 64 %.
с, g/1 10
1
0.1
0.01
0.001
10 20 30 40 50 60 70 80 90 100 110 120 130 140 ™
Рис. 4. Распределение концентрации взвешенных частиц через 2 часа после сброса грунта, полученное на основе 2D и 3D моделей
Горизонтальная ось на рис. 5 проходит через зону проведения дноуглубительных работ и направлена вдоль течения. На рис. 5, а и 5, б по вертикальной оси отложена глубина водоема (в метрах), ось направлена вертикально вниз. На рис. 5, в и 5, г по вертикальной оси отложен уровень донных отложений (в миллиметрах), ось направлена вертикально вверх. На рис. 5, а цветом представлена концентрация фракции А в воде. На рис. 5, б представлены значения концентрации фракции Б в воде. На рис. 5, в и 5, г цветом отражен долевой состав фракций А и Б в донных отложениях соответственно.
Из рис. 5 видно, что более тяжелая фракция А осаждается ближе к зоне проведения дноуглубительных работ и залегает глубже в осадочных породах, чем фракция Б.
§6. Заключение
Разработан комплекс объединенных математических моделей, описывающих пространственно-трехмерные процессы переноса многокомпонентных взвешенных частиц и двумерные процессы транспорта донных материалов. Предложенные модели гидрофизики учитывают следующие факторы: движение водной среды; переменную плотность, зависящую от концентрации взвеси; многокомпонентность взвеси; изменение геометрии дна в результате осаждения взвеси. Предложен алгоритм расчета изменения состава донных материалов в случае многокомпонентной взвеси. Разработка данного алгоритма позволит моделировать движение наносов в случае повторного взмучивания многокомпонентных донных отложений водоема.
Аппроксимация трехмерного уравнения диффузии-конвекции выполнена на основе расщепления на двумерную и одномерную задачи. В работе используются дискретные аналоги операторов конвективного и диффузионного переносов в случае частичной заполненности ячеек. Геометрия расчетной области описывается на основе функции заполненности. При этом использована схема, представляющая собой линейную комбинацию разностных схем «крест» и «кабаре» с весовыми коэффициентами, полученными в результате минимизации погрешности аппроксимации. Данная схема предназначена для решения задачи переноса примеси при больших сеточных числах Пекле.
В рамках данной работы построен параллельный алгоритм, реализующий поставленную трехмерную задачу переноса веществ на многопроцессорной вычислительной системе Научно-технологического университета «Сириус». Разработанный на основе технологии MPI параллельный алгоритм, используемый для решения трехмерных задач диффузии-конвекции, реализованный на больших расчетных сетках (размером от 1000 х 1000 х 100
73.38
2
4
55.04 36.69 6 18.35 0
133.1 99.82 66.55 33.72 0
0
200 400 600 800
0
200 400 600 800
(а) концентрация взвеси в воде (фракция А) (б) концентрация взвеси в воде (фракция Б)
0.56 0.42 0.28 0.14 0
0.75 0.56 0.38 0.19 0
0
0
200 400 600 800
(в) концентрация взвеси в донных отложени- (г) концентрация взвеси в донных отложениях (фракция А) ях (фракция Б)
Рис. 5. Концентрация взвешенных частиц и гранулометрический состав дна
0
0
расчетных узлов) показал высокую эффективность, составившую более 75 %.
Приведены результаты численных экспериментов, на основании которых сделаны выводы о преимуществе 3Б модели транспорта многокомпонентной взвеси по сравнению с 2Б моделью. Разработанный программно-аналитический инструментарий позволяет проанализировать процесс движения наносов в случае взмучивания многокомпонентных донных отложений водоема, вызывающий вторичное загрязнение водоема, ухудшение его экологического состояния. Выполнены численные эксперименты по моделированию процесса осаждения многокомпонентной взвеси, а также по изучению его влияния на рельеф дна и изменение его состава.
Финансирование. Исследования выполнены при финансовой поддержке Российского научного фонда в рамках научного проекта 21-71-20050.
СПИСОК ЛИТЕРАТУРЫ
1. Ковтун И. И., Проценко Е. А., Сухинов А. И., Чистяков А. Е. Расчет воздействия на водные биоресурсы дноуглубительных работ в Белом море // Фундаментальная и прикладная гидрофизика. 2016. Т. 9. № 2. С. 27-38.
2. Чернявский А. В. Трансформация донных зооценозов в районе Григоровской свалки грунта // Дноуглубительные работы и проблемы охраны рыбных запасов и окружающей среды рыбохо-зяйственных водоемов. Астрахань, 1984. С. 208-210.
3. Лаврентьева Г. М., Суслопарова О. Н. Итоги рыбохозяйственных мониторингов, проводимых в восточной части Финского залива с целью оценки воздействия гидротехнических работ на гид-робионтов // Фонды ГосНИОРХ. 2006. Т. 1. Вып. 331. С. 5-11.
4. Иванова В. В. Экспериментальное моделирование заваливания зообентоса при дампинге грунтов // Фонды ГосНИОРХ. 1988. Вып. 85. С. 107-113.
5. Корректировка «Проекта на разработку месторождения песков «Сестрорецкое», расположенного в Финском заливе Балтийского моря» в связи с реконструкцией карьера. СПб., 2012. 50 с.
6. Морозов А. Е. Донная фауна малых рек и влияние на нее взвешенных веществ дренажных вод // Фонды ГосНИОРХ. 1979. Вып. 2. С. 108-113.
7. Гайджуров П. П., Савельева Н. А., Дьяченков В. А. Конечно-элементное моделирование совместной работы оползня скольжения и защитного сооружения // Advanced Engineering Research. 2021. Т. 21. № 2. С. 133-142. https://doi.org/10.23947/2687-1653-2021-21-2-133-142
8. Соловьев А. Н., Бинь До Т., Лесняк О. Н. Поперечные колебания круглого биморфа с пьезоэлектрическим и пьезомагнитным слоями // Вестник Донского государственного технического университета. 2020. Т. 20. № 2. С. 118-124. https://doi.org/10.23947/1992-5980-2020-20-2-118-124
9. Соловьев А. Н., Тамаркин М. А., Тхо Н. В. Метод конечных элементов в моделировании центробежно-ротационной обработки // Вестник Донского государственного технического университета. 2019. Т. 19. № 3. С. 214-220. https://doi.org/10.23947/1992-5980-2019-19-2-214-220
10. Самарский А. А., Вабищевич П. Н. Численные методы решения задач конвекции-диффузии. М.: Едиториал УРСС, 1999.
11. Белоцерковский О. М., Гущин В. А., Щенников В. В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Журнал вычислительной математики и математической физики. 1975. Т. 15. № 1. С. 197-207. http://mi.mathnet.ru/zvmmf8315
12. Sukhinov A. I., Chistyakov A. E., Protsenko E.A., Sidoryakina V. V., Protsenko S.V. Accounting method of filling cells for the solution of hydrodynamics problems with a complex geometry of the computational domain // Mathematical Models and Computer Simulations. 2020. Vol. 12. No. 2. P. 232-245. https://doi.org/10.1134/S2070048220020155
13. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.
14. Sukhinov A. I., Chistyakov A. E., Kuznetsova I. Y., Protsenko E.A. Modelling of suspended particles motion in channel // Journal of Physics: Conference Series. 2020. Vol. 1479. Issue 1. 012082. https://doi.org/10.1088/1742-6596/1479/1/012082
15. Сухинов А. И., Чистяков А. Е., Проценко Е.А. Разностная схема для решения задач гидродинамики при больших сеточных числах Пекле // Компьютерные исследования и моделирование. 2019. Т. 11. № 5. С. 833-848. https://doi.org/10.20537/2076-7633-2019-11-5-833-848
16. Гущин В. А., Никитина А. В., Семенякина А. А., Сухинов А. И., Чистяков А. Е. Модель транспорта и трансформации биогенных элементов в прибрежной системе и ее численная реализация // Журнал вычислительной математики и математической физики. 2018. Т. 58. № 8. С. 110-137. https://doi.org/10.31857/S004446690002007-8
17. Воеводин В. В., Воеводин Вл. В. Параллельные вычисления. СПб: БХВ-Петербург, 2010.
18. Четверушкин Б. Н., Якобовский М. В. Вычислительные алгоритмы и архитектура систем высокой производительности // Препринты ИПМ им. М. В. Келдыша. 2018. № 52. 12 с. https://doi.org/10.20948/prepr-2018-52
19. Сухинов А. И., Чистяков А. Е., Шишеня А. В., Тимофеева Е. Ф. Предсказательное моделирование прибрежных гидрофизических процессов на многопроцессорной системе с использованием явных схем // Математическое моделирование. 2018. Т. 30. № 3. С. 83-100. http://mi.mathnet.ru/mm3950
20. Sukhinov A. I., Chistyakov A. E., Filina A. A., Nikitina A. V., Litvinov V. N. Supercomputer simulation of oil spills in the Azov Sea // Bulletin of the South Ural State University. Series «Mathematical Modelling, Programming and Computer Software». 2019. Vol. 12. No. 3. P. 115-129. https://doi.org/10.14529/mmp190310
Поступила в редакцию 12.08.2022
Принята в печать 18.10.2022
Сухинов Александр Иванович, д. ф.-м. н., профессор, чл.-корр. РАН, заведующий кафедрой математики и информатики, Донской государственный технический университет, 344002, Россия, г. Ростов-на-Дону, пл. Гагарина, 1. ORCID: https://orcid.org/0000-0002-5875-1523 E-mail: sukhinov@gmail.com
Чистяков Александр Евгеньевич, д. ф.-м. н., профессор, Донской государственный технический университет, 344002, Россия, г. Ростов-на-Дону, пл. Гагарина, 1. ORCID: http://orcid.org/0000-0002-8323-6005 E-mail: cheese_05@mail.ru
Атаян Ася Михайловна, ассистент, Донской государственный технический университет, 344002, Россия, г. Ростов-на-Дону, пл. Гагарина, 1. ORCID: https://orcid.org/0000-0003-4629-1002 E-mail: atayan24@mail.ru
Кузнецова Инна Юрьевна, старший преподаватель, Южный федеральный университет, 344006, Россия, г. Ростов-на-Дону, ул. Большая Садовая, 105/42. ORCID: https://orcid.org/0000-0003-1996-1605 E-mail: ikuznecova@sfedu.ru
Литвинов Владимир Николаевич, к. т. н., доцент, доцент, Донской государственный технический университет, 344002, Россия, г. Ростов-на-Дону, пл. Гагарина, 1. ORCID: https://orcid.org/0000-0001-8234-3194 E-mail: litvinovvn@rambler.ru
Никитина Алла Валерьевна, д. т. н., доцент, профессор, Донской государственный технический университет, 344002, Россия, г. Ростов-на-Дону, пл. Гагарина, 1;
профессор, Южный федеральный университет, 344006, Россия, г. Ростов-на-Дону, ул. Большая Садовая, 105/42.
ORCID: http://orcid.org/0000-0001-7257-962X E-mail: nikitina.vm@gmail.com
Цитирование: А. И. Сухинов, А. Е. Чистяков, А. М. Атаян, И. Ю. Кузнецова, В. Н. Литвинов, А. В. Никитина. Математическая модель процесса осаждения на дно многокомпонентной взвеси и изменения состава донных материалов // Известия Института математики и информатики Удмуртского государственного университета. 2022. Т. 60. С. 73-89.
A. I. Sukhinov, A. E. Chistyakov, A. M. Atayan, I. Yu. Kuznetsova, V. N. Litvinov, A. V. Nikitina Mathematical model of process of sedimentation of multicomponent suspension on the bottom and changes in the composition of bottom materials
Keywords: suspension transport model, variable density, Upwind Leapfrog difference scheme, Standard Leapfrog difference scheme, bottom topography change, parallel algorithms.
MSC2020: 65Q10,65Y05
DOI: 10.35634/2226-3594-2022-60-05
The paper considers 2D and 3D models of transport of suspended particles, taking into account the following factors: movement of aqueous medium; variable density depending on the suspension concentration; multicomponent character of suspension; changes in bottom geometry as a result of suspension sedimentation. The approximation of the three-dimensional diffusion-convection equation is based on splitting schemes into two-dimensional and one-dimensional problems. In this work, we use discrete analogues of convective and diffusion transfer operators in the case of partial cell occupancy. The geometry of the computational domain is described based on the occupancy function. The difference scheme used is a linear combination of the Upwind and Standard Leapfrog difference schemes with weight coefficients obtained by minimizing the approximation error. This scheme is designed to solve the problem of impurity transfer at large grid Peclet numbers. Based on the results of numerical experiments, conclusions are drawn about the advantage of the 3D model of multicomponent suspension transport in comparison with the 2D model. Computational experiments have been performed to simulate the process of sedimentation of a multicomponent suspension, as well as its effect on the bottom topography and changes in its composition.
Funding. The study was funded by Russian Science Foundation, project number 21-71-20050.
REFERENCES
1. Kovtun 1.1., Protsenko E.A., Sukhinov A. I., Chistyakov A. E. Calculating the impact on aquatic resources dredging in the White Sea, Fundamentalnaya i Prikladnaya Gidrofizika, 2016, vol. 9, no. 2, pp. 27-38 (in Russian).
2. Chernyavskii A. V. Transformation of bottom zoocenoses in the area of the Grigorovskaya landfill, Dnouglubitel'nye raboty i problemy okhrany rybnykh zapasov i okruzhayushchei sredy ry-bokhozyaistvennykh vodoemov (Dredging and problems of protection of fish stocks and the environment of fishery reservoirs), Astrakhan: 1984, pp. 208-210 (in Russian).
3. Lavrenteva G.M., Susloparova O.N. Results of fisheries monitoring conducted in the eastern part of the Gulf of Finland to assess the impact of hydrotechnical work on aquatic organisms, Funds of GosNIORKh, 2006, vol. 1, issue 331, pp. 5-11 (in Russian).
4. Ivanova V.V. Experimental modeling of zoobenthos collapse during soil dumping, Funds of GosNIORKh:, 1988, issue 85, pp. 107-113 (in Russian).
5. Adjustment of the "Project for the development of the sand deposit "Sestroretskoye", located in the Gulf of Finland of the Baltic Sea" in connection with the reconstruction of the quarry, Saint Petersburg: LLC Eco-Express-Service, 2012, 50 p. (In Russian).
6. Morozov A.E. Bottom fauna of small rivers and the influence of suspended solids of drainage waters on it, Funds of GosNIORKh, 1979, issue 2, pp. 108-113 (in Russian).
7. Gaidzhurov P.P., Saveleva N.A., Dyachenkov V. A. Finite element modeling of the joint action of flow slide and protective structure, Advanced Engineering Research, 2021, vol. 21, no. 2, pp. 133-142 (in Russian). https://doi.org/10.23947/2687-1653-2021-21-2-133-142
8. Solov'ev A.N., Binh Do T., Lesnyak O.N. Transverse vibrations of a circular bimorph with piezoelectric and piezomagnetic layers, Vestnik of Don State Technical University, 2020, vol. 20, no. 2, pp. 118-124 (in Russian). https://doi.org/10.23947/1992-5980-2020-20-2-118-124
9. Soloviev A.N., Tamarkin M. A., Tho N.V. Finite element modeling method of centrifugal rotary processing, Vestnik of Don State Technical University, 2019, vol. 19, no. 3, pp. 214-220 (in Russian). https://doi.org/10.23947/1992-5980-2019-19-2-214-220
10. Samarskii A.A., Vabishchevich P.N. Chislennye metody resheniya zadach konvektsii-diffuzii (Numerical methods for solving problems of convection-diffusion), Moscow: Editorial URSS, 1999.
11. Belotserkovskii O.M., Gushchin V. A., Shchennikov V. V. Use of the splitting method to solve problems of the dynamics of a viscous incompressible fluid, USSR Computational Mathematics and Mathematical Physics, 1975, vol. 15, issue 1, pp. 190-200. https://doi.org/10.1016/0041-5553(75)90146-9
12. Sukhinov A. I., Chistyakov A. E., Protsenko E.A., Sidoryakina V. V., Protsenko S.V. Accounting method of filling cells for the solution of hydrodynamics problems with a complex geometry of the computational domain, Mathematical Models and Computer Simulations, 2020, vol. 12, no. 2, pp. 232-245. https://doi.org/10.1134/S2070048220020155
13. Samarskii A. A., Nikolaev E. S. Metody resheniya setochnykh uravnenii (Methods for solving grid equations), Moscow: Nauka, 1978.
14. Sukhinov A. I., Chistyakov A. E., Kuznetsova I. Y., Protsenko E.A. Modelling of suspended particles motion in channel, Journal of Physics: Conference Series, 2020, vol. 1479, issue 1, 012082. https://doi.org/10.1088/1742-6596/1479/1/012082
15. Sukhinov A. I., Chistyakov A. E., Protsenko E.A. Difference scheme for solving problems of hydrodynamics for large grid Peclet numbers, Computer Research and Modeling, 2019, vol. 11, no. 5, pp. 833-848 (in Russian). https://doi.org/10.20537/2076-7633-2019-11-5-833-848
16. Gushchin V. A., Sukhinov A. I., Nikitina A. V., Chistyakov A. E., Semenyakina A. A. A model of transport and transformation of biogenic elements in the coastal system and its numerical implementation, Computational Mathematics and Mathematical Physics, 2018, vol. 58, issue 8, pp. 1316-1333. https://doi.org/10.1134/S0965542518080092
17. Voevodin V. V., Voevodin Vl. V. Parallel'nye vychisleniya (Parallel computing), Saint Petersburg: BKHV-Peterburg, 2010.
18. Chetverushkin B.N., Yakobovskiy M. V. Numerical algorithms and architecture of HPC systems, Keldysh Institute Preprints, 2018, issue 52, 12 p. (In Russian). https://doi.org/10.20948/prepr-2018-52
19. Sukhinov A. I., Chistyakov A. E., Shishenya A. V., Timofeeva E. F. Predictive modeling of coastal hy-drophysical processes in multiple-processor systems based on explicit schemes, Mathematical Models and Computer Simulations, 2018, vol. 10, issue 5, pp. 648-658. https://doi.org/10.1134/S2070048218050125
20. Sukhinov A.I., Chistyakov A.E., Filina A. A., Nikitina A. V., Litvinov V.N. Supercomputer simulation of oil spills in the Azov Sea, Bulletin of the South Ural State University. Series "Mathematical Modelling, Programming and Computer Software", 2019, vol. 12, issue 3, pp. 115-129. https://doi.org/10.14529/mmp190310
Received 12.08.2022 Accepted 18.10.2022
Aleksandr Ivanovich Sukhinov, Doctor of Physics and Mathematics, Full Professor, Corresponding Member, Russian Academy of Sciences, Head of the Department of Mathematics and Computer Science, Don State Technical University, Gagarin square, 1, Rostov-on-Don, 344002, Russia. ORCID: https://orcid.org/0000-0002-5875-1523 E-mail: sukhinov@gmail.com
Aleksandr Evgen'evich Chistyakov, Doctor of Physics and Mathematics, Professor, Don State Technical University, Gagarin square, 1, Rostov-on-Don, 344002, Russia. ORCID: http://orcid.org/0000-0002-8323-6005 E-mail: cheese_05@mail.ru
Asya Mikhailovna Atayan, Assistant Lecturer, Don State Technical University, Gagarin square, 1, Rostov-on-Don, 344002, Russia.
ORCID: https://orcid.org/0000-0003-4629-1002 E-mail: atayan24@mail.ru
Inna Yur'evna Kuznetsova, Senior Lecturer, Southern Federal University, ul. Bol'shaya Sadovaya, 105/42,
Rostov-on-Don, 344006, Russia.
ORCID: https://orcid.org/0000-0003-1996-1605
E-mail: ikuznecova@sfedu.ru
Vladimir Nikolaevich Litvinov, Candidate of Engineering, Associate Professor, Assistant Professor, Don State Technical University, Gagarin square, 1, Rostov-on-Don, 344002, Russia. ORCID: https://orcid.org/0000-0001-8234-3194 E-mail: litvinovvn@rambler.ru
Alla Valer'evna Nikitina, Doctor of Engineering, Associate Professor, Professor, Don State Technical University, Gagarin square, 1, Rostov-on-Don, 344002, Russia;
Professor, Southern Federal University, ul. Bol'shaya Sadovaya, 105/42, Rostov-on-Don, 344006, Russia; ORCID: http://orcid.org/0000-0001-7257-962X E-mail: nikitina.vm@gmail.com
Citation: A. I. Sukhinov, A. E. Chistyakov, A. M. Atayan, I. Yu. Kuznetsova, V. N. Litvinov, A.V. Nikitina. Mathematical model of process of sedimentation of multicomponent suspension on the bottom and changes in the composition of bottom materials, Izvestiya Instituta Matematiki i Informatiki Udmurtskogo Gosu-darstvennogo Universiteta, 2022, vol. 60, pp. 73-89.
