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

математическое моделирование mathematical modeling

© 0

Я Check for updates

Научная статья

УДК 519.6

https://doi.org/10.23947/2587-8999-2023-7-3-39-48

Математическая модель транспорта трехкомпонентной взвеси А.И. Сухинов1 , И.Ю. Кузнецова1,2 H

'Донской государственный технический университет, Российская Федерация, г. Ростов-на-Дону, пл. Гагарина, 1 2НИЦ супер-ЭВМ и нейрокомпьютеров, Российская Федерация, Ростовская область, г. Таганрог, пер. Итальянский, 106

Н [email protected]

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

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

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

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

Ключевые слова: транспорт взвеси, многокомпонентная взвесь, трехмерная модель гидродинамики, схемы расщепления, численные методы

Финансирование. Исследование выполнено за счет гранта Российского научного фонда № 22-11-00295. https:// rscf.ru/project/22-11-00295/

Для цитирования. Сухинов А.И., Кузнецова И.Ю. Математическая модель транспорта трехкомпонентной взвеси. Computational Mathematics and Information Technologies. 2023;7(3):39-48. https://doi.org/10.23947/2587-8999-2023-7-3-39-48

Аннотация

© А.И. Сухинов, И.Ю. Кузнецова, 2023

Original article

Mathematical Model of Three-Component Suspension Transport Alexander I Sukhinov1 , Inna Yu Kuznetsova1,2 H

1Don State Technical University, 1, Gagarin Square, Rostov-on-Don, Russian Federation Supercomputers and Neurocomputers Research Center, 106, Italian Lane, Taganrog, Russian Federation

H [email protected] Abstract

Introduction. Prediction of suspension deposition zones is required to assess and minimize the negative impact on the ecosystem of the reservoir during dredging within the framework of large-scale engineering projects, prediction of suspension deposition zones is required to assess and minimize the negative impact on the ecosystem of the reservoir. It is necessary to build a mathematical model that takes into account many factors that have a significant impact on the accuracy of forecasts. The aim of the work is to construct a mathematical model of transport of multicomponent suspension, taking into account the composition of the soil (different diameter of the suspension particles), the flow velocity of the water flow, the complex geometry of the coastline and bottom, wind stresses and friction on the bottom, turbulent exchange, etc.

Materials and Methods. A mathematical model for the transport of a multicomponent suspension and an approximation of the proposed continuous model with the second order of accuracy with respect to the steps of the spatial grid are described, considering the boundary conditions of the Neumann and Robin type. The approximation of the hydrodynamics model is obtained based on splitting schemes by physical processes, which guarantee fulfillment mass conservation for discrete model.

Results. The proposed mathematical model formed the basis of the developed software package that allows to simulate the process of sedimentation of a multicomponent suspension. The results of the work of the software package on the model problem of sedimentation of a three-component suspension in the process of soil dumping during dredging are presented.

Discussions and Conclusions. The mathematical model of transport of three-component suspension is described and software implemented. The developed software allows to simulate the process of deposition of suspended particles of various diameters on the bottom, and to evaluate its effect on the bottom topography and changes in the bottom composition. The developed software package also allows to analyze the process of sediment movement in the case of resuspension of multicomponent bottom sediments of the reservoir, which causes secondary pollution of the reservoir.

Keywords: suspension transport, multicomponent suspension, three-dimensional hydrodynamics model, splitting schemes, numerical methods

Funding information. The research was carried out at the expense of the grant of the Russian Science Foundation no. 22-11-00295. https://rscf.ru/project/22-11-00295/

For citation. Sukhinov AI, Kuznetsova IYu. Mathematical Model of Three-Component Suspension Transport.

Computational Mathematics and Information Technologies. 2023;7(3):39-48. https://doi.org/10.23947/2587-8999-2023-7-3-39-48

Введение. Реализация масштабных инженерных проектов, таких как строительство мостов, расширение акватории, доступной для судоходства, требует проведения работ, оказывающих существенное влияние как на рельеф дна, так и на экосистему водоема в целом. Например, при проведении дноуглубительных работ в воду поступает значительное количество взвеси, которая в процессе оседания на дно или вторичного взмучивания может негативно повлиять на продукционно-деструкционные процессы водной экосистемы [1-2]. Для оценки возможного ущерба, наносимого экосистеме в процессе дампинга грунта при проведении дноуглубительных работ, необходимо предварительно оценить области акватории, в которых будет происходить оседание взвеси и в которых возможно ее взмучивание, что ведет к вторичному загрязнению водного объекта. Для прогнозирования зон осаждения взвешенных частиц предложена математическая модель транспорта взвеси на основе системы начально-краевых задач, включающей расчет гидродинамических характеристик акватории и транспорта взвеси.

Для предложенной трехмерной модели транспорта многокомпонентной взвеси описан подход к аппроксимации непрерывной модели со вторым порядком точности относительно шагов пространственной сетки с учетом граничных условий второго и третьего рода. Предложенная математическая модель транспорта взвешенных частиц дополнена трехмерной моделью гидродинамики, позволяющей рассчитывать поля вектора скорости водного потока [3-4]. Предлагаемая математическая модель легла в основу разработанного программного комплекса, позволяющего моделировать процесс осаждения многокомпонентной взвеси. Приведены результаты работы программного комплекса на модельной задаче осаждения трехкомпонентной взвеси в процессе дам-пинга грунта при проведении дноуглубительных работ.

Материалы и методы

1. Постановка задачи. Для построения математической модели транспорта многокомпонентной взвеси воспользуемся уравнением диффузии-конвекции, которое может быть записано в следующем виде [3]:

К )'t+(ucr) x + К); +((w + ) Cr)z ' = (n (cr) x ') ' + (cr) ' +{y(cr) ,') ' + Fr , (1)

где сг — концентрация г-ой фракции взвеси, мг/л; У={м,у,м>} — компоненты вектора скорости течения водного потока, м/с; V ,г — скорость осаждения г-ой фракции взвеси, м/с; д, V — горизонтальная и вертикальная составляющие коэффициента турбулентного обмена, соответственно, м2/с; Ег — функция, описывающая интенсивность распределения источников г-ой фракции взвеси, мг/(лс).

Выражение (1) рассматривается при следующих начальных сг(х, у, 2, 0) = с0г(х, у, 2) и граничных условиях:

- на свободной поверхности: (с п)2 = 0

- вблизи поверхности дна: у(сг )г' = -ws,гсг;

- на боковой поверхности: (сг )п' = 0, если (V, п) > 0 и д(сг )п' = (V, п)сг, если (V, п) < 0 , где (V, п) — нормальная составляющая вектора скорости, п вектор нормали, направленный внутрь расчетной области.

Уравнение диффузии-конвекции (1) дополним трехмерной моделью гидродинамики мелководных водоемов [5] для расчета вектора скорости течения водного потока:

- уравнение движения (Навье-Стокса):

и;+иих +ш.'у +м>и'г = - р;/р+(ци;и +{]ш'у)'у + (Ю,

u;+uuv +vuW +wuZ =-py )р + (К )X + (LL) uy) У + (vuZ ) ',

<U;wX ww+yuww; =-p-/pu(i+wLУ u (j+w'Wu (vwVVZ u g,

(2)

у z г/ г х/х У'У

- уравнение неразрывности в случае переменной плотности:

р;+(р« Ух Уу = о. (3)

где Р — давление, Па; р — плотность, кг/м3; g — ускорение свободного падения, м/с2. Система уравнений (2)-(3) рассматривается при следующих граничных условиях:

- на входе: и = и0, V = v0, Рш' = 0, Уп' = 0;

- боковая граница (берег и дно): рдиП = -тх, рдVП = -ту, Уп = 0, Рп' = 0;

- верхняя граница: РДК =-т*, РД^ =-Ту, w =-ю - р '/(ря), Р„' = 0,

где а — интенсивность испарения жидкости, тх, ту — составляющие тангенциального напряжения.

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

0+М;+(ус);+= {д о^у+(д су'У+{ю/ )', (4)

где компонента скорости V неявным образом учитывает скорость осаждения рассматриваемой фракции взвеси ws.

Введем равномерную сетку по времени ют = и" = п т; п = 0, И,; И,т = Тгде т — шаг по времени, Ы{ — количество временных слоев, Т — длительность интервала моделирования.

Предположим, что расчетная область вписана в параллелепипед О = {0 < х < Ьх,0 < у < Ьу,0< г < Ь2}. Замыкание области О получим путем присоединения к О граней параллелепипеда, то есть определим как О = {0<х < Ьх,0 < у < £у,0 < = < }.

Таким образом приходим к цепочке начально-краевых задач:

(с")('+<1у(у-с" ) = <Цу(к • гас! с"), (5)

где k = { ц, ц, V} — коэффициент турбулентного обмена, (х, у, 2, 0 е О х [0 < / < Т], < t < сп (х, у, 2, tn-l) = с"-1 (х, у, 2, t"-1), (х, у, 2, (х, у, 2) е G.

При этом начальные и граничные условия запишутся в виде:

- начальное условие: с (х, у, 2, 0) = с0 (х, у, 2), (х,у,г) е О;

- граничное условие на свободной поверхности: (сг)г' = 0;, (х;у,х)е Е0;

- граничное условие надн^: =(с")' = п, (х, 2) е Ея ;

- граничное условие на боковой поверхности: (с")П = 0, если (Л^, п) > 0 и ц(сг )п' = (V, п)сг, если (V, п) < 0, где Е0 = {0 < х < Ь, 0 <у < Ь, 2 = 0} — верхняя грань параллелепипеда О, Ея = {0 < х < Ь, 0 <у < Ь , 2 = Ь} — нижняя грань параллелепипеда О..

Запишем член, описывающий конвективный перенос взвешенных веществ в симметричной форме. Такой подход при дискретизации задачи позволит построить разностный оператор, обладающий свойством кососимметричности [6]:

(V • сп) = -

дс" дсп дсп д( п) д( п) д( п)

-+ у-+ ы-+ — \исп 1+ — |усп 1+ — Ыс")

дх ду ¿Зг дхУ ду ' дхХ '

(6)

Для аппроксимации задачи (5) введем равномерную расчетную сетку <= = <ох х <вух юг

ю.

= {= : х = гИх; I = 0,Мх; кхМх=Ьх},

щ =\у.:у. = ]Ъ; ] = 0,N ; И N = I },

У Т J у J ' У' У У У У'

СО.

= = к=0,Nг■ к2м2= Ь2},

где к, ку, к^ — шаги по пространственным координатным направлениям Ох, Оу и О, соответственно, Ых, Ыу, Ых — количество узлов расчетной сетки по каждому из пространственных направлений, Ь, Ьу, Ь2 — длина расчетных интервалов по каждому из пространственных направлений.

Множество внутренних узлов расчетной сетки обозначим как юх, ю юх.

Аппроксимацию задачи (5) на пространственно-временной сетке сот х со выполним с заданием скоростей (коэффициентов конвективного переноса) в узлах сеток, сдвинутых на полшага по координатным направлениям О, и О .

Для оператора конвективного переноса, заданного в симметричной форме (6) в уравнении (5), имеем:

Ссп (и" (х + 0,5кх, у, г)^" {х+ кх, у, г)-и"(х-0,5кх, у, г )• с" (х - кх, у, г))+

2кх

+ — {у"(х,у + °,5ку ,г) с"(х,у+ ку,г)-у"(х,у- 0,5},,г)• с"(х,у - ку,г))+

2 к

(7)

+ — (х, у, г + 0,5 к2 )• сп (х, у,х + к2)-™п(х, у, г -0,5кг )• сп (х, у, 2 - к2)), 2к

где сп обозначает сеточную функцию с" = с(х,у,х,tn), (х= у, х) е ю, а через с" — достаточно гладкую функцию непрерывных переменных (х, у, х, /).

Для оператора диффузионного переноса в уравнении (5) имеем:

Б-п = —

ц (х + 0,5кх, у, г)-

п (х + кх, У, 2)---(x, У, 2)

-ц(х - 0,5кх, у, г ) -ц(х, у ~0,5ку, 2 )

"(x,У,г)-с"(х - кх,у,г)

1

+ — к

(

ц(х, у + 0,5ку, г )С

" (х, у + ку, г)-с" (х, у, г )

с" (х, у, 2)-с" (х, у - ку, 2 1 (

-г(х, у, х - 0,5кх )•

+ — к

у(х, у, 2 + 0,5к2)

с" (х, у, 2 + к2)-с" (х, у, 2)

(8)

2(х,У, х)-о"(х,у, х - кх)

к

С учетом записанных аппроксимаций (7)-(8) получим следующий вид аппроксимации уравнения (5) во внутренних узлах сетки:

—п —п-1

с - с 1

т + [и"(х+°>5Ь* >у,2 }с"(х+кх, у, 2)~и"(х - , у,2 У(х - , у, г)) +

+ 2к~ [[ V" (х У + °'5лу ' 2 у (х у + ьу , 2)- V" (х, у - °,5Лу, г )с" (х, у - к у, 2))) +

+ -1-[wn (х,У, 2 + °,5к2 )с (х,У, 2 + к2)- wn (х,у, 2 - °,5к2 )с" (х,у, 2 - к2=

2к2 v )

/ , \с (х + к , х, 2) - с (х, У, 2) / , \с (х, У, 2)- с (х - к , У, 2) ц(х + °,5кх,У,2)-^-- ц(х-°,5кх,У,2) ^ >у> > К-

кх кх

»(х,У + °,5ку, 2)(х,У,2)-с" (х,У,2) - ^х,У - °,5ку, 2)с" (х,у,2)-с" (х,У - ку,2)

1

+ — к

У

—п / \ —п / \ —п / \ —п / \

( \ с (х,У,2 + к2)-с (х,У,2) / ^с (х,У,2)-с (х,у,2 - к2) у(х, У, 2 + °,5к2)---—-4 ' - у(х, У, 2 - °,5к2)—4 '--—

и

к

к

х

х

к

к

Л

У

к

к

у

2

1

к

Л

+

к

У

2

2

Дополним полученную аппроксимацию (9) начальными и граничными условиями. Для задания граничных условий на дне, свободной и боковой поверхностях рассматриваемой области О, удобно ввести расширенную сетку [7] ш* = ш* х ш* х ш*, где:

Ч = & : x = ihx ; i =-1, Nx +1; hxNx = Lx], = {y. : y. = A ; j = -1, +1; hyNy= T }

У V J '7 J ■'У •> 'У ' У У У )

Ш* = {^ : 2к=кИ2; к = -1,Лг +1; к2Л2 = Ь2}. В дальнейшем будем предполагать, что:

еп{х,у, г) = 0, (х, у, z) е ш * \ ш , (10)

где ш \ ш — граничные узлы сетки ш .

Также считаем известными значения компонент вектора скоростей водной среды в узлах сетки ш \ ш с дробными значениями индексов. Для узлов сетки ш \ ш , которые находятся вне водоема, значения компонент вектора скорости задаются равными нулю.

В случае потоков на боковых гранях области О, совпадающих по направлению с внешними нормалями к граням (случай (V, п) > 0), имеют место граничные условия Неймана. Данный случай можно записать в виде:

ип (0,5кх, у, 2) + ми(- 0,5йх, у, 2) <0, и" (Ьх - 0,5кх ,у,2) + и"(Ьх + 0,5/, у, т) > 0, V"(х,0,5ку,г) + V"(х,-0,5ку,7) < 0, V"(х,0,5ку,7) + у"(х,-0,5ку,г) < 0, Запишем аппроксимацию граничных условий второго рода для оператора конвективного переноса. Рассмотрим случай х = 0, 0< у < L , 0 < 2 < Ь2. В этом случае в качестве разностной аппроксимации конвективного члена можно рассматривать выражение:

^(и"(0,5кх,y,z)ën(hx,y,z)-un(- 0,5hx, y,z)C"{-kx,y,z)). (12)

Выражение (12) аппроксимирует конвективный член с погрешностью 0(к2х). Помимо формы (12) аппроксимацию конвективного члена с погрешностью 0(к2х) можно записать в виде:

С (кх , 2)-с"(- кх, у,2) =о 2к '

откуда получаем:

€"(-кх ,у,2) = С" (кх,у,2). (13)

Из аппроксимаций (12) и (13) получаем:

С (о- ) О- (кх,у,2)(ип(0,5к=, у, N)- ип(- 0,5кх,у, 2)). (14)

'*=о 2кх

Аналогично записываются случаи х = Ьх, у = 0, у = Ьу, 2 = 0 (граничное условие на свободной поверхности). Получим:

с*с)^ "2Т'С"(ь*-К'У,+°>5К'У'*)-""(ь*-У,*)), (15)

Су {е") = • Сп{х,ку,2)\уп(х,0,5ку, 2)-уп(х,-0,5ку, 2)), (16)

|у-° 2пу

1 _

С,(оп)^ - ^ ■ с"(ХLy-ку,2)\vix,Ly + 0 ,5kу,2)- v"(x,Ly -0,5ky,z)), (17)

Cz (с" ) z=0 = ^Г ' ^ X У,h) "И*, y>°>5hz)-W У,-0,5Аг )). (18)

1 _,

/|г=0 2к

Запишем аппроксимацию граничных условий второго рода для оператора диффузионного переноса. Рассмотрим случай х = 0, 0 < у < Ь , 0 < 2 < Ь2. В этом случае в качестве разностной аппроксимации диффузи-

онного члена на расширенной сетке можно рассматривать выражение:

( -п

/ \| I

Б

0 h

- \

- м(-0,5hx, y, z )

cn (h-, y, z )-cn (0, y, z ) h-

P(0, y, 2)-cn{-hx,y, 2 )T

( ») 1 ) c (h-,У,z)-c"(0,y,z)

(c)-=0 " Г ^{°,5h-,y,z)-h----(19)

При выполнении первого условия из (11) получаем, что с"(-Нх,у,г) = с"(Их,у,г). Тогда с учетом последнего равенства и выражения (19), получаем следующую аппроксимацию оператора диффузионного переноса в случае

х = 0, 0 < у < Ь , 0 < г <Х :

У г

С ) х=0 = ,у, х )+ц(-0,5}, уо))(с"(кх, у,х}-с" (0, у, г ))). (20)

Аналогично записываются случаи х = Ь, у = 0, у = Ь. Например, в случае х = Ь, 0 < у < Ьу, 0 < г < Ь, при выполнении второго условия из (11) и с учетом равенства сп(Ьх +ИХ,у,г) = сп(Ьх- Их,у,г), получаем:

С ) " 12 + °>5К> + - У>2))Н1*- К, у, г)-с (4, У, г})). (21)

Аналогично для у = 0 и у = Ьу:

°уС) =0 "ЧлЫхОЫу, у)+цу*,-0,5йуо])(с"(сс,})-Сп(с,0,(22)

°у С ) , = АМх, ЬУ + 0,5ку, т)+ц(х, 1у- 0,5}, г)) (сп(х, 1у- ку, г)-(х, 1у, г))). (23)

|у=ьу ку

Для свободной поверхности (случай у = 0) с учетом у(х, у, - 0,5кг) = 0 аппроксимация диффузионного оператора запишется в виде:

В2 (с") п - у(х,у,0,5к2)(с"(х,у,к2)-с"(х,у,0)). (24)

и-0 к2

Рассмотрим аппроксимацию оператора диффузионного переноса на донной поверхности (г = Ь) на расширенной сетке. Формально аппроксимация диффузионного члена может быть записана в виде:

Л (п ) 1 {( Т+ с, \ и*' У'^+^ЬТ (х>У'к ))

Л (С 1 ^ " (Х' У' ^ + У ^-}-

(сп (х,у,12) -сп(х, у,12-к2))

к

(25)

-у(х,у,Ь2 - 0,5})-

Воспользуемся аппроксимацией второго порядка точности 0(к2г) граничного условия третьего рода у(с")г = , (х, у, z) е Ея соотношением:

( ч(с"(х,у,4 + к2)-сп(х,у,Ь2 -к2У) _п( . (26) 2к,

Воспользуемся соотношением:

у(х, "у,

Подставляя (27) в (26), получим:

у(х, У, ь2) = 1 (у (х, У, Ь2 + 0,5}) + V (х, У, Ь2 - 0,5к2)) + О (к] ). (27)

(( , ч ( т А с. ))(с"У'+с"(ХУ'-)) -п( т\ (28)

(у(х,у,Ь2 + 0,5к2)+у(х,у,4 -0,5Аг))----= -ыс (х,у,Ь2).

4 к

Из (28) имеем:

сп(х,у,Ь2 + к2) = -^--^сп(х,у,Ь2) + сп(х,у,Ь2 -к2). (29)

1 2/ Нх,У,Ь2 + 0,5к2)+у(х,у,12 -0,5к2)) 2> 2 2> ^

Тогда аппроксимация диффузионного члена в уравнении диффузии-конвекции взвеси имеет вид:

(с" ) ^ - ^^-^ П'Т ^-, С" (х, у, 4 ) +

\z-Lz Ь VI х 1, Т _1_П Ы -I- 1/1V т; Т — П I

у(х,у, L2 + °,5Az) + у(х,у, Lz - 0,5})

+ -1 у(х,у, 12 + 0,5})(сп(х,у, 12 - к2)- сп(х,у, 12))- (30)

-у(х,у,ь2 - 0,5к2 Г 2, 2 '

к

Аналогичным способом можно получить аппроксимацию конвективного оператора при г = Ь:

С (с") ^Г 2 + )_ с" (х, у, )+ (31)

12—0 р! 1/1 V Л1 Т 4-П Ы 14- 1/1V л) Т — П Я/? I

^ '12—ь2 2к.

у(х,у,Ь2 + 0,5})+у(х,у, 12 - 0,5})

к

2

+ (х, у, Ь2 + 0,5йг) - ^ (х, у, Ь2 - 0,5Аг))с" (х, у, Ь2 - к2)).

Полученные аппроксимации операторов диффузионного (30) и конвективного (31) переносов в граничных узлах (при 2 = Ь) подходят для дна с различными морфологическими характеристиками («жидкое» дно, непроницаемое дно и т. п.) при задании соответствующим образом коэффициента турбулентного обмена V.

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

3. Аппроксимация трехмерной модели гидродинамики. Аппроксимацию модели (2)-(3) будем проводить на расчетной сетке ю = ют х юк. Для аппроксимации модели (2)-(3) воспользуемся схемами расщепления по физическим процессам [8]. Согласно данному методу, исходная модель гидродинамики (2)-(3) разобьется на три подзадачи [6; 9].

Первая подзадача представлена уравнением диффузии-конвекции, на основе которого вычисляются компоненты поля вектора скорости водного потока на промежуточном вре менном слое:

и и- + ийХ + + м^щ ' = ( + (к' )у + (УП2 x ,

т

v - v + uvX + vvy + wV2' = + (pT')v + (vv2)z,

W W + UWx '+ VWy '+ WWz '=(VWx)x + (vWX + ^

где и, V, w — компоненты вектора скорости на предыду 1цем временном слое; и,V, ^ — компоненты вектора скорости на промежуточном временном слое; и = си +(1 - с)м, с е [0,1] — весовой коэффициент или вес схемы.

На основе второй подзадачи вычисляется поле давления:

Р ' + Р" + Р ' = Р^Р , (р~)'* , (р~) у , (р"0 ''

1 хх ^ 1 УУ ^ 1 гг 2

Т Т Т Т

На основе третьей подзадачи по явным формулам вычисляются компоненты поля вектора скорости водного потока на следующем временном слое:

и — и 1 , V — V 1 , ^ — ^ 1 ,

-=--г , -=--г , -=--г ,

х' у' г'

т р т р т р где и, V, V — компоненты вектора скорости на текущем временном слое.

Аппроксимация задачи расчета поля скорости движения водной среды по пространственным переменным выполнена на основе метода баланса.

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

- скорость течения водного потока на основе системы уравнений (2)-(3);

- процесс транспорта взвешенных частиц в толще воды с учетом полученной скорости течения водного потока;

- процесс оседания взвеси на дне на основе модели (1)-(3).

В качестве примера работы программного комплекса приведем результаты численного моделирования задачи транспорта трехкомпонентной взвеси при моделировании процесса дампинга грунта в ходе проведения дноуглубительных работ.

Параметры расчетной области: длина — 1 км; ширина — 720 м; глубина — 10 м.

Параметры расчетной сетки: шаги по горизонтальной и вертикальной пространственным координатам составили 10 и 1 м, соответственно; расчетный интервал — 2 ч, шаг по времени — 5 с. Ось О направлена вдоль расчетной области, ось Оу — по ширине расчетной области, ось О2 — по глубине расчетной области (от 0 до -10 м, где отметка 0 соответствует поверхности воды, -10 - дну водоема).

Входные параметры модели: среднее расстояние от точки выгрузки грунта до дна водоема в районе проведения дноуглубительных работ составляет 5,5 м; область выгрузки грунта вдоль оси О (по длине водоема) размещена в пределах от 200 до 250 м; область выгрузки грунта вдоль оси Оу (по ширине водоема) размещена в пределах от 300 до 400 м; скорость течения на глубинах от 4 до 10 м составляла 0,075 м/с (течения направлены слева направо); плотность пресной воды при нормальных условиях — 1000 кг/м3; плотность взвесей — 2700 кг/м3; коэффициент формы частицы для всех трех взвесей — 0,222 (шарообразная форма); начальная вязкость воды — 1,002 мПа-с (при температуре 20 °С); диаметр частицы фракции А — 0,05 мм; скорость осаждения фрак-

ции А — 2,31 мм/с; процентное содержание фракции А — 20 %; диаметр частицы фракции Б — 0,04 мм; скорость осаждения фракции Б — 1,48 мм/с, процентное содержание фракции Б — 30 %; диаметр частицы фракции В — 0,03 мм; скорость осаждения фракции В — 0,83 мм/с, процентное содержание фракции В — 50 %.

На рисунке 1 представлены результаты моделирования процесса транспорта трехкомпонентной взвеси в толщи воды. Горизонтальная ось направлена вдоль течения, срез представлен посередине расчетной области, где наблюдается максимальная концентрация взвешенных частиц (в плоскости у = 360 м).

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

а

0,0

-2,5

-5,0

-7,5

мг/л

в

(0,22]

(22, 44] м

(44, 66] и, и

(66, 88] р о

(88, 110] отав

(110, 132] кв а

(132, 154] а н

(154, 176] и 1С

(176, 198] ^ л

(198, 220] и

0,0

-2,5

-5,0

-7,5

мг/л

В

(0,0; 4,4] (4,4; 8,8]

: (8,8; 13,2]

(13,2; 17,6]

■ (17,6; 22,0]

(22,0; 26,4]

(26,4; 30,8]

(30,8; 35,2]

(35,2; 39,6]

(39,6; 44,0]

-10,0

250 500 750 1000 Длина акватории, м

а)

-10,0

0 250 500 750 1000 Длина акватории, м б)

и

и р

н и

1С

а

0,0

-2,5

-5,0

-7,5

мг/л

(0,0; 6,6]

(6,6; 13,2]

(13,2; 19,8]

(19,8; 26,4]

■ (26,4; 33,0]

(33,0; 39,6]

(39,6; 46,2]

(46,2; 52,8]

(52,8; 59,4]

(59,4; 66,0]

а

0,0

-2,5

-5,0

-7,5

мг/л

■ (0,11]

(11,22]

■ (22,33]

(33,44]

(44,55]

□ (55,66]

■ (66,77]

(77,88]

(88,99]

(99,110]

0

-10,0

-10,0

250 500 750 1000

250

500 750 1000

Длина акватории, м Длина акватории, м

в) г)

Рис. 1. Концентрация взвеси в толще воды спустя 2 часа после выгрузки: а) — общая концентрация взвеси; б) — фракция А; в) — фракция Б; г) — фракция В

Обсуждение и заключения. В работе представлена трехмерная математическая модель транспорта многокомпонентной взвеси, дополненная трехмерной моделью гидродинамики мелководного водоема. Представленная модель учитывает состав грунта (различный диаметр частиц взвеси), скорость течения водного потока, сложную геометрию береговой линии и дна, сгонно-нагонные явления, ветровые течения и трение о дно, турбулентный обмен, что позволяет повысить точность моделирования.

Аппроксимация предложенной модели транспорта многокомпонентной взвеси на основе трехмерного уравнения диффузии-конвекции выполнена со вторым порядком точности относительно шагов пространственной сетки с учетом граничных условий второго и третьего рода. Аппроксимация трехмерной математической модели гидродинамики выполнена на равномерной прямоугольной расчетной сетке с использованием схем расщепления по физическим процессам.

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

0

0

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

Список литературы

1. Матишов Г.Г., Ильичев В.Г. Об оптимальной эксплуатации водных ресурсов. Концепция внутренних цен. Доклады Академии наук. 2006;406(2):249-251.

2. Ковтун И.И., Проценко Е.А., Сухинов А.И. и др. Расчет воздействия на водные биоресурсы дноуглубительных работ в Белом море. Фундаментальная и прикладная гидрофизика. 2016;9(2):27-38.

3. Сухинов А.И., Чистяков А.Е., Атаян А.М. Математическая модель процесса осаждения на дно многокомпонентной взвеси и изменения состава донных материалов. Известия Института математики и информатики Удмуртского государственного университета. 2022;60:73-89. https://doi.org/10.35634/2226-3594-2022-60-05

4. Сухинов А.И., Кузнецова И.Ю., Чистяков А.Е. и др. Исследование точности и применимости разностной схемы для решения задачи диффузии-конвекции при больших сеточных числах Пекле. Вычислительная механика сплошных сред. 2020;13(4):437-448. https://doi.Org/10.7242/1999-6691/2020.13.4.34

5. Кузнецова И.Ю., Сухинов А.И., Чистяков А.Е. и др. Математическая модель гидродинамики устьевых районов. Сборник трудов международной научной конференции «Актуальные проблемы прикладной математики, информатики и механики». 2021:960-965.

6. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. Москва: URSS. 2009:248.

7. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. Москва: Наука. 1978:592

8. Белоцерковский О.М., Гущин В.А., Щенников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости. Журнал вычислительной математики и математической физики. 1975;15(1):197-207. https://doi.org/10.1016/0041-5553(75)90146-9

9. Тихонов А.Н., Самарский А.А. Уравнения математической физики. 7-е изд. Москва: Наука: Издательство Московского университета, 2004:798.

References

1. Matishov GG, Ilyichev VG. On optimal exploitation of water resources. The concept of internal prices. Reports of the Academy of Sciences. 2006;406(2):249-251. (In Russ.).

2. Kovtun II, Protsenko EA, Sukhinov AI, et al. Calculating the Impact on Aquatic Resources Dredging in the White Sea. Fundamental and AppliedHydrophysics. 2016;9(2):27-38. (In Russ.).

3. Sukhinov AI, Chistyakov AE, Atayan AM, et al. Mathematical model of process of sedimentation of multicomponent suspension on the bottom and changesin the composition of bottom materials, Izvestiya Instituta Matematiki i Informatiki Udmurtskogo Gosudarstvennogo Universiteta, 2022;60:73-89. (In Russ.). https://doi.org/10.35634/2226-3594-2022-60-05

4. Sukhinov AI, Kuznetsova IYu, Chistyakov AE. Study of the accuracy and applicability of the difference scheme for solving the diffusion-convection problem at large grid peclet numbers. Computational Continuum Mechanics. (In Russ.). https://doi.org/10.7242/1999-6691/2020.13.4.34

5. Kuznetsova IYu, Sukhinov AI, Chistyakov AE, et al. Mathematical model of hydrodynamics of estuarine areas. Proceedings of the International scientific conference "Actual problems of applied Mathematics, computer science and Mechanics". 2021:960-965. (In Russ.).

6. Samarskiy AA, Vabishevich PN. Numerical methods for solving convection-diffusion problems. Moscow: URSS. 2009:248. (In Russ.).

7. Samarskiy AA, Nikolaev ES. Methods for solving grid equations. Moscow: Nauka. 1978:592. (In Russ.).

8. Belotserkovsky OM, Gushchin VA, Schennikov VV. Splitting method applied to solving problems of dynamics of viscous incompressible fluid. Journal of Computational Mathematics and Mathematical Physics. 1975;15(1):197-207. (In Russ.). https://doi.org/10.1016/0041-5553(75)90146-9

9. Tikhonov AN, Samarskiy AA. Equations of mathematical physics. 7th ed. Moscow: Nauka: Moscow University Press, 2004:798. (In Russ.).

Об авторах:

Сухинов Александр Иванович, член-корреспондент РАН, доктор физико-математических наук, профессор, директор НИИ Математического моделирования и прогнозирования сложных систем, Донской государственный технический университет (РФ, 344003, г. Ростов-на-Дону, пл. Гагарина, 1), ORCID, MathSciNet, eLibrary.ru, ResearcherID, ScopusID, [email protected]

Кузнецова Инна Юрьевна, старший преподаватель кафедры интеллектуальных и многопроцессорных систем, НИЦ супер-ЭВМ и нейрокомпьютеров (РФ, 347900, г. Таганрог, пер. Итальянский, 106), ORCID, Scopus ID, elibrary.ru, ResearcherID, [email protected]

Заявленный вклад соавторов:

Все авторы сделали эквивалентный вклад в подготовку публикации.

Поступила в редакцию 02.08.2023 Поступила после рецензирования 24.08.2023 Принята к публикации 25.08.2023

Конфликт интересов

Авторы заявляют об отсутствии конфликта интересов.

Все авторы прочитали и одобрили окончательный вариант рукописи.

About the Authors:

Alexander I Sukhinov, Corresponding Member of the Russian Academy of Sciences, Doctor of Physical and Mathematical Sciences, Professor, Director of the Research Institute for Mathematical Modeling and Forecasting of Complex Systems, Don State Technical University (1, Gagarin Sq., Rostov-on-Don, 344003, RF), ORCID, MathSciNet, eLibrary.ru, ResearcherID, ScopusID, [email protected]

Inna Yu Kuznetsova, Senior Lecturer of the Department of Intelligent and Multiprocessor Systems, SIC Supercomputers and Neurocomputers (106, Italian Lane, Taganrog, 347900, RF), ORCID, ScopusID, elibrary.ru, ResearcherID, [email protected]

Claimed contributorship:

all authors have made an equivalent contribution to the preparation of the publication.

Received 02.08.2023 Revised 24.08.2023 Accepted 25.08.2023

Conflict of interest statement

the authors do not have any conflict of interest.

All authors have read and approved the final manuscript.