Имитационная модель загрязнения взвешенными веществами мелководных водных объектов

ИМИТАЦИОННАЯ МОДЕЛЬ ЗАГРЯЗНЕНИЯ ВЗВЕШЕННЫМИ ВЕЩЕСТВАМИ МЕЛКОВОДНЫХ ВОДНЫХ ОБЪЕКТОВ

Полупанов Владимир Николаевич,

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

РФ, г. Керчь

 

АННОТАЦИЯ

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

 

Ключевые слова: водный объект, мелководная зона, адвекция, диффузия, взвешенные вещества, седиментация, концентрация, осадок, клеточные автоматы, двумерная модель, web-программа.

 

Введение

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

Здесь рассматривается двумерная имитационная модель, применимая для расчётов загрязнения на небольших участках водных объектов (водоёма, водотока) с учётом влияния мелководья и берегов. В модели используется кинетический механизм дискретных клеточных автоматов, локальность которых позволяет учитывать особенности мелкомасштабной турбулентности мелководья “простыми средствами без утраты сложности” [30]. Правила автомата основаны на обобщениях из источников, перечисленных в списке литературы.

Кратко описывается программный комплекс [22], реализующий предложенную модель.

1 Постановка задачи

Распространение взвешенных веществ в мелководных зонах в большинстве случаев имитируется с помощью двумерных сеточных моделей, в которых используются уравнения “мелкой воды” - упрощения уравнений Навье-Стокса и Эйлера [5, 14], применяющиеся для участков, горизонтальные размеры которых во много раз превышают глубину. Основные проблемы использования таких моделей возникают из-за недостаточности существующих теорий турбулентности, применимых для непротиворечивого замыкания исходных уравнений [14]. Для традиционного решения этих проблем с использованием полуэмпирических параметризаций необходимы натурные исследования, которые, как правило, отсутствуют. Кратко возникающие проблемы можно проиллюстрировать с помощью известной модели горизонтальной турбулентной диффузии (ГТД-модели) [19], описывающей перемещение частиц плоского облака инертной взвеси в турбулентном потоке, векторное поле скорости u которого имеет вид:

где ū - вектор адвективной составляющей, обеспечивающий перемещение (адвекцию) частиц облака в направлении основного потока;

u' - вектор диффузионной составляющей, вызывающая растекание (диффузию) облака за счёт мелкомасштабных разнонаправленных колебаний.

При постоянной глубине поле скорости u соответствует нормальному закону. Предполагается, что в этом поле действует мгновенный точечный источник с производительностью P по взвеси, порождающий облако, концентрация взвеси С в котором представима в виде гауссовой функции. Положение центра x₀ этой функции прослеживается во времени t с помощью уравнения dx₀/dt = ū(x₀, t). ГТД-модель позволяет воспроизвести усреднённую по ансамблю всех реализаций во времени концентрацию взвеси С в двумерном однородном и изотропном потоке жидкости с постоянной глубиной h, движущемся со средней скоростью ū в горизонтальной декартовой системе координат (x, y), с помощью дифференциального уравнения переноса и диффузии (адвекции-диффузии). Последнее для краткости можно записать, ограничившись направлением вдоль оси x и для случая, когда мгновенный точечный источник мощностью P находится в точке (x = 0; y = 0):

где w - доля осаждающейся взвеси; K_ℓ(t) — коэффициент горизонтальной турбулентной диффузии, часто задаваемый в упрощённом виде:

Решение для такой упрощённой задачи при t > 0 можно получить аналитически в виде гауссовой функции распределения усреднённой концентрации взвеси:

где σ² = 2K_ℓ·t – дисперсия; σ определяет полуширину облака загрязнения.

Выражение (4) позволяет наблюдать за состоянием плоского радиально-симметричного диффундирующего “облака” усреднённой по глубине концентрации взвеси в заданный момент времени t при соблюдении условия изотропии случайной составляющей u' в однородном поле скорости ū.

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

где K_ℓ(ū_i,j) — тензор коэффициентов горизонтальной турбулентной диффузии; M - масса взвеси.

Аналитические решения для квазиплоского потока отсутствуют, поэтому необходимо использовать численные аппроксимации. Однако, использование сеточных моделей в мелководных прибрежных зонах, как упоминалось, связано с трудностями. Проблемы сеточных моделей можно попытаться “обойти” с помощью альтернативных бессеточных методов, разработанных на основе клеточных автоматов (cellular automata, CA или КА) [13, 31].

Модель должна описывать перенос взвеси квазиплоскими течениями (2^') в расчётной КА-области с ячейками (клетками) квадратной формы (в терминах КА - окрестность Мура [31]) с набором скоростей D2Q9 (два измерения, девять соседей) с использованием некоторых идей метода дискретных решёточных уравнений Больцмана (lattice Boltzmann method, LBM) [15,34].

Для наименования модели в дальнейшем будет использоваться аббревиатура CAD2Q9 (CAD2Q9-модель), а для её программной реализации - CAD2Q9-SIST (Cellular Automata D2Q9-Simulation Suspension Transfer) [22]. Вычислительная эффективность CAD2Q9-модели обеспечивается за счёт идеи, похожей на идею «мгновенной максвелизации» скоростей турбулизованного потока [15], при которой компоненты скорости в клетках расчётной области напрямую доступны из распределений плотности.

В модели достаточно учитывать ветровые течения и сток, с учётом влияния рельефа дна и берегов, так как в большинстве случаев на эти факторы приходится основная доля течений бесприливного мелководного объекта на участках с ограниченными возможностями для развития волн [1, 6, 27].

Предполагаются выполнимыми следующие условия:

– ветер и расход воды считаются заданными;

– данных о батиметрии достаточно для восстановления глубин;

– численное решение, по аналогии с LBM, может состоять из нескольких этапов: на первых – предварительное решение подобное (4), а затем - окончательное, путём адаптации к граничным условиям и законам сохранения:

– на открытых границах - свободное протекание;

– для береговой линии - непроницаемость и отражение методом отскока;

– у дна – скольжение с торможением, зависящими от глубины и её градиентов;

– в каждой клетке - приближение к условия сохранения баланса импульсов;

– в целом по расчётной области - приближение к сохранению вещества.

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

где U_i,j - вектор турбулентного течения;

Ū_i,j - вектор основного потока;

U'_ij – отклонения от вектора основного потока.

2 Обзор источников

2.1 Течения мелководья. Сразу же следует отметить чрезвычайно небольшое количество сведений о вероятностных свойствах течений непосредственно мелководной прибрежной зоны (с глубинами до 5 м и до 1-2 км от берега), которое удалось обнаружить, хотя именно в этой зоне наиболее часто выполняются гидротехнические работы.

Ветровые и стоковые течения обычно рассматриваются отдельно, хотя во многом проблемы моделирования в прибрежных зонах морей и в водотоках совпадают, поскольку в обоих случаях используются приближения уравнений “мелкой воды”, только в первом случае движение воды обязано ветру, а во втором - уклону русла. Наиболее часто используемые упрощения приводят к разделению исходного течения на две составляющие: усреднённый по глубине поток и крупномасштабный вихрь или крупномасштабное квазипериодическое возмущение [5, 14]. Первая составляющая применительно к нашей задаче (установившееся течение) описывает векторы скорости основного потока Ū^W_i,j и Ū^Q_i,j или сумму этих векторов.

Можно выделить две основных теоретических проблемы вычисления скоростей основных потоков при использовании уравнений “мелкой воды”. Первая состоит в том, что для вычислении поля давления (или уровня свободной поверхности) используется остаточный член из уравнения неразрывности [14, 20], который по сути является подгоночным (не вполне физичным) параметром при нарушениях гладкости рельефа дна, свойственных прибрежным мелководным участкам.

Вторая проблема связана с определением закона сопротивления осреднённому по глубине турбулентному потоку. Для практических приложений используется либо постоянный коэффициент сопротивления, либо коэффициент сопротивления, зависящий от “шероховатости”, тогда как анализ многочисленных прямых инструментальных измерений даёт основание утверждать, что коэффициент сопротивления изменяется в различных водных объектах по ещё не установленным законам, а “шероховатость” - это не устоявшееся понятие, включающее в себя много факторов (“зернистой шероховатости”, скорости потока, геометрии и гидравлического сопротивления русла и др.), роль которых к настоящему времени не установлена [2, 3]. Традиционные методы расчётов с использованием коэффициентов шероховатости (по Маннингу, Павловскому, Гангилье-Куттера или табличных по Срибному, Чоу и т.п.) приводят к большим погрешностям (для скоростных коэффициентов Шези - от 30% до 100%) и считаются устаревшими [9].

В уравнениях “мелкой воды” типа квазилинейных (бездивергентных) уравнений Сен-Венана, традиционно используемых в речной гидравлике, также используются, упомянутые коэффициенты шероховатости, за счёт чего недопустимо упрощается физическая сторона процесса [12]. Кроме этого, огрубление геометрии потока за счёт чрезмерного сглаживания приводит к существенному искажению пространственной картины векторов скоростей течения [4, 12]. Такое огрубление исходной геометрии русла позволяет получать приемлемые результаты для достаточно протяжённых участков, при досточно подробных исходных данных об уклонах дна и водной поверхности, глубин, скоростей течений и расходах на нескольких гидростворах. Но при моделировании на ограниченном достаточно широком участке реки те допуски, что заложены в уравнениях Сен-Венана, являются некорректными [12]. Первоисточник этих некорректностей лежит в самой формулировке гипотезы Шези и последующей подгонки решений путём подбора валидирующего коэффициента (“эффективной шероховатости”) [2]. Ошибки и неоднозначность расчётов возрастают из-за больших погрешностей измерений уклонов и уровней при уменьшении длины участка реки [9].

Отсутствие исходных данных, требующихся для использования и верификации моделей на основе уравнений “мелкой воды” и, в первую очередь, достаточно точных измерений уклонов русла (свободной поверхности), связано с трудоёмкостью традиционных наземных и дистанционных гидрологических наблюдений, тогда как сейчас массово доступны новые, спутниковые и беспилотниковые измерения, с помощью которых возможно построение цифровых моделей водного объекта. Известны применения для этой цели малых беспилотных летательных и водных аппаратов (БПЛА [28] и БПВА [26]), оснащённых недорогими измерительными устройствами, прежде всего такими как: лазерный сканер, эхолот, профилограф скорости течения.

Следует отметить отсутствие сведений о вероятностных характеристиках турбулентности, применимых для моделирования на небольших участках рек и и узкой прибрежной полосы моря, водохранилища. Эту проблему также можно решать новыми методами, в том числе с использованием недорогих дрифтеров, оснащённых GPS-приёмником и GSM-связью [24].

Кроме сведений, полученных при использованиb уравнений “мелкой воды”, для определения скорости основного потока за счёт ветровой составляющей Ū^W_i,j можно использовать эмпирические обобщения, полученные по данным натурных экспериментов, проводимым для различных целей в прибрежных зонах морей и водных объектах типа озёр и водохранилищ [6, 18, 19, 23, 27, 29]. Эти обобщения позволяют сделать несколько наиболее важных для нашего случая выводов, применительно к установившимся течениям в мелководной прибрежной полосе, возбуждаемых ветром с относительно небольшими скоростями (до 8 м/сек), при которых разрешается выполнение гидротехнических работ.

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

2. Для ветровых течений характерна закономерность увеличения скоростей при уменьшении глубин и с приближением к берегу.

3. Эпюры скоростей по вертикали в большинстве случаев имеют вид, близкий к пологой сигмоиде и течения в большинстве случаев можно считать однонаправленными, поэтому возможно усреднение по вертикали.

Стоковая составляющая основного потока Ū^Q_i,j для нашего случая должна быть представлена с учётом распределения скоростей по площади (плановые течения). Использование для этой цели линеаризованных двумерных уравнений Сен-Венана приводит к чрезмерно большим ошибкам [3, 4, 12]. Поэтому на практике применяется несколько приближённых эмпирических подходов. При одном из наиболее популярных: имея в качестве заданных расход и глубины для поперечного створа (разреза, живого сечения), вычисляется средняя скорость потока и затем эта средняя скорость (Ū) и глубины (H_i) используются для расчета усредненной по глубине скорости потока на i-той вертикали (Ū_i) створа по известному методу Н.М. Бернадского. Однако использование традиционных формулы Маннинга и коэффициентов шероховатости обесценивают этот метод и в гидрологических расчётах допускается введение упрощённых фиктивных величин, выполняющих роль коэффициентов шероховатости [2]. Используются также упрощённые методы расчёта плановых течений применительно к конкретным практическим задачам русловой гидравлики [4,25]. Представляется, что применительно к нашей задаче целесообразно использовать упрощённый расчёт распределения скорости потока по живому сечению водотока с использованием новых возможностей построения цифровых моделей русла и эмпирических соотношений между гидроморфометрическими характеристиками [9].

Усреднённые значения турбулентных колебаний течений U'_i,j, как следует из обобщений, касающихся мелководных водных объектов [1, 5, 6, 10, 17, 18, 19, 23, 27, 29], можно объяснить влиянием крупномасштабных инерционных колебаний водных масс, постепенно становящихся трудно различимыми при приближении к берегу из-за усиления веса морфологических возмущений. На основании этих обобщений можно сделать следующие выводы.

1. В мелководных зонах прослеживается зависимость коэффициента горизонтального турбулентного обмена K_ℓ(U) от средней скорости основного потока [23, 29]. Разброс наблюдаемых значений вокруг зависимости можно объяснить тем, что при её построении не учитывалось влияние глубины [17,18].

2. В мелководных зонах не подтверждаются универсальные зависимости для открытого моря типа степенных «закона 4/3» Ричардсона-Обухова или «закона 5/3» Колмогорова [19], хотя и отмечается связь инерционных колебаний с характерными масштабами Ʌ, изменяющимися от сотен до нескольких тысяч метров в зависимости от размеров и морфологии водного объекта [1, 19, 29]. Для учёта такого влияния в прибрежной зоне можно использовать линейные или близкие к линейным зависмости от масштаба [19], где под масштабом понимается шаг расчётной сетки. Отмечается также возможность нелинейной зависимости между коэффициентами горизонтальной и вертикальной турбулентности [17, 18, 19].

3. Характерную для прибрежной зоны анизотропию горизонтальной турбулентности [6, 19, 27] можно объяснить, учитывая изменение с глубиной интенсивности вертикальной турбулентности при нелинейной зависимости между коэффициентами горизонтальной и вертикальной турбулентности [17, 18, 19].

В связи с последним выводом следует обратить внимание на предположение, вытекающее из осесимметричной теории С. Чандрасекхара (упоминается в [19]), применимое для квазиплоского потока. Cогласно этой теории связь между интенсивностями горизонтальной и вертикальной турбулентностями безсдвигового осесимметричного течения эллипсоидального вида можно получить из выражения для тензора второго ранга. Тогда качественно картину турбулентности можно представить в виде эллипсоида вращения, “сплющенность” которого зависит от глубины. При увеличении глубины увеличивается вес K_z и, соответственно, уменьшается вес K_ℓ. Если глубины уменьшаются, то это, наоборот, приводит к уплощению (распластованию) потока с учётом обратной связи между K_ℓ и K_z.

Эффект влияния обратной связи между горизонтальной и вертикальной диффузией учитывают введением некоторого эффективного коэффициента горизонтальной диффузии K_l^эф, который трактуется как добавочный к горизонтальному переносу взвеси и возникает при усреднении по глубине трёхмерного поля скорости с вертикальным градиентом (сдвигом) скорости течения [10,19]. В связи с этим в мелководной зоне повышается значимость точности восстановления глубин.

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

Изменчивость глубин и геометрия русла способствуют неустойчивости течения и возникновению когерентных структур [3] (сверхбольших вихрей или вихрей первого порядка по А.Н. Колмогорову), обладающих большой пульсационной энергией и играющих главную роль в процессах транспорта наносов. Несмотря на важность существования таких структур следует признать неизученность влияния этого явления на турбулентность, поэтому на практике вводится эмпирическое понятие “эффективная турбулентность” [3, 4]. Особенно важно, что когерентные структуры и характеристики макромасштабной турбулентности существенным образом связаны с геометрическими характеристиками русла и плановым распределением скоростей потока, а продольные масштабы существенно возрастают и могут на порядок превышать масштабы вертикальной и поперечной турбулентности [3].

Турбулентность в узкой прибрежной зоне существенно отличается от наблюдаемой в потоке при удалении от береговой линии. Перепады глубин вблизи бровки берегов способствуют увеличению неоднородности скоростей и нарушению двумерной анизотропии, возникают дополнительные (геоморфологические, гидрологоморфологические и гидрологоморфо-метрические [9]) источники турбулентности. Другими словами - резко возрастает дисперсия течений по направлению и скорости. Дисперсия течений всё более сказывается при приближении к берегу, когда большой вес в кинематике потока начинают оказывать перепады глубин [6] и всё большее значение для горизонтальной диффузии приобретает сила тяжести, когда транспортирующая способность потока вверх по уклону меньше, чем по уклону вниз. Наконец поток, достигая бровки берега, кардинально меняет направление в зависимости от угла между направлением потока и линией берега. В крайних случаях, при приближении угла к нулю поток направлен вдоль берега, а при направлении по нормали - поток отражается от берега. В промежуточных случаях направление потока меняется по более сложным законам. Представляется, что в общем случае поведение потока у береговой линии подчиняется законам, похожим на известные законы отражения по Френелю, имитируемым в имитационных моделях с использоваием условия отскока.

В связи с выделенными особенностями турбулентности возрастает роль учёта геометрии водного объекта и его локальных морфометрических характеристик. В то же время следует ожидать, что учёт таких особенностей может вызвать проблемы с соблюдением законов сохранения при построении модели, как это происходит при использовании линеаризованных двумерных уравнений [14, 20].

2.2 Распространение и седиментация взвеси. Следует отметить, что в нашем случае задача упрощается из-за небольших скоростей потока (менее 0,5 м/с), обычных при небольших скоростях ветра и немноговодных периодов, когда чаще всего выполняются гидротехнические работы. При таких скоростях обычно не учитываются влияние вертикального турбулентного обмена на скорость осаждения частиц взвеси и взмучивание у дна [32].

Особое внимание уделяется необходимости учёта полидисперсности при существенном различии размеров частиц фракций во взвеси, так как при использовании средневзвешенной скорости осаждения получаемые результаты могут существенно отличаться от результатов с учётом полидисперсности. В последнем случае обычно рассчитывается распространение нескольких невзаимодействующих монодисперсных шлейфов (облаков). Для получения итоговой концентрации независимых облаков используется принцип суперпозиции. В работе [10] метод суперпозиции облаков заменяется использованием так называемой «эффективной гидравлической крупности», зависящей от времени распространения взвеси. Однако манипулирование отдельными облаками оставляет более широкие возможности в случае учёта взаимозависимости облаков, например при флокуляции [32]. Флокуляция, наиболее знáчимо влияющая на результаты моделирования, может проявляться при наличии мелких пылевидных и глинистых вязких частиц полидисперсной взвеси, состоящей из нескольких фракций (fr =1,..,N).

Выделяют два основных типа осаждения полидисперсных частиц в зависимости от их концентрации: свободное осаждение и стеснённое осаждение. Свободное осаждение наблюдается при низкой концентрации частиц Ο(10⁰ мг/л), когда можно считать, что частицы оседают независимо друг от друга, не сталкиваясь и не взамодействуя друг с другом [32]. В этом случае скорость осаждения частиц размером 0,05 мм и менее рассчитывается непосредственно по формуле Стокса [32].

Для определения плотности воды используется эмпирическая завимость от температуры (T,°C) и солёности (S,‰) при стандартном атмосферном давлении, принятая для Международного уравнения состояния (УС-80) [7].

Стеснённое осаждение начинает проявляться при повышенных (Ο(10¹мг/л)) концентрациях мелких частиц [32]. В этом случае начинают сказываться две противонаправленные силы. Одна из них замедляет оседание из-за увеличения плотности образующейся смеси и сталкивания частиц. Другая, наоборот, ускоряет оседание при наличии агрегативно неустойчивых частиц из-за флокуляции, когда образуются более тяжёлые агрегаты (флокулы). Интенсивность образования флокул зависит в первую очередь от концентрации и в меньшей мере от солёности и других факторов.

Следуя работе Ван Рийна [32] можно использовать подход для расчёта осаждения частиц мелких (менее 0,05 мм) вязких фракций, при котором к формуле Стокса вводятся поправочные коэффициенты на затруднённое и ускоренное осаждение из-за флокуляции.

При определении толщины осадка следует учитывать недоуплотнение (разрыхлённость) осадка по сравнению с природной плотностью, вводя коэффициент консолидации [8].

3 Моделирование течений

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

Обозначим клеточную (расчётную) область в виде квадрата размером L с числом клеток n×n и размером клетки ∆x = L/n. Положение центральной клетки окрестности в пределах клеточной области можно определить с помощью индексов i (номер строки) и j (номер столбца), а в пределах окрестности - с использованием следующего механизма индексации относительно центральной клетки: i + l, j + m; l = - 1, 0, 1; m = - 1, 0, 1. Тогда компоненты D2Q9-набора скорости для имитации установившегося турбулентного течения в окрестности Мура с центральной клеткой с координатами i = 1, j = 1 можно представить в виде схем, показанных на рисунке 1.

а)	б)	в)
Рисунок 1. Становление D2Q9-набора компонент скорости течения: а) диффузионные компоненты, сгенерированные вектором основного потока; б) суммарные компоненты адвективно-диффузионного течения до адаптации; в) суммарные компоненты адвективно-диффузионного течения после адаптации (векторы основного потока: - на начальном этапе; - после адаптации)

 

Рисунок 1 демонстрирует особенности подхода, состоящего из нескольких этапов. На начальном этапе ветер и (или) сток возбуждают первоначальный вектор основного потока во влияющих клетках, расположенных с наветренной стороны (выше по течению) от центральной клетки. Начальный вектор генерирует в клетках окрестности девять компонент одинаковой величины Ũ'_i+n,j+m, определяющих диффузию (рисунок 1-а).

Разложение вектора основного потока по проекциям на два соседних направления и суммирование со сгенерированными диффузионными компонентами даст компоненты скорости адвективно-диффузтонного потока в виде D2Q9-набора Ũ_i+n,j+m (рисунок 1-б). Горизонтальные компоненты Ũ_i+n,j+m при n,m ≠ 0 обеспечивают адвекцию и горизонтальную диффузию (адвективную диффузию), а вертикальная Ũ_i,j (n,m = 0) отвечает за тормозящий эффект, обусловленный гидродинамическим сопротивлением потоку.

Затем полученные девять компонент трансформируются в D2Q9-набор компонент U_i+n,j+m (рисунок 1-в) путём адаптации к перепадам глубин окрестности при сохранении первоначально заданного импульса, которым обладал D2Q9-набор Ũ_i+n,j+m. После адаптации горизонтальные компоненты U_i+n,j+m (n,m ≠ 0) каждого D2Q9-набора обеспечивают адвективную диффузию с учётом влияния рельефа дна, а вертикальная компонента U_i,j (n,m = 0) отвечает за тормозящий эффект, теперь уже с учётом сохранения импульса.

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

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

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

Далее излагаются основные алгоритмы процедуры, составляющие CAD2Q9-модель с учётом обобщений обзора в п. 2.1.

3.2 Основной поток. Векторы основного потока_i,j за счёт ветровой составляющей определяются с использованием эмпирической формулы [27]:

(5)

где _i,j - скорость потока (м/мин), усреднённая по глубине H_i,j (м);

W - скорость ветра (м/с) на высоте 10 м;

b₀ = 0.70, b₁ = 0.33, b₂ = 1.5.

Направление в модели задаётся в диапазоне от 0° до 360°, для ветра - “в компáс”, для основного потока - “из компáса”.

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

Вычисляем плановую скорость потока. Для этого вводим относительную глубину (ћ_i) на вертикалях, расположенных равномерно на поперечном створе:

где – средняя глубина по живому сечению.

Затем пределяем скорость на i-той вертикали:

(7)

где q_i – коэффициент:

λ₁ = ~ 1,8÷2,0; λ₂ ~ 2÷3.

Значения λ₁ и λ₂ на каждом из створов уточняются исходя из наилучшего приближения расчётного расхода воды к исходному. Такую оптимизацию параметров можно осуществить, например, с помощью подпрограммы (Javascript-функции) [21]. С помощью λ₁ имитируется торможение из-за гидродинамического сопротивления потоку: уменьшение (увеличение) λ₁ соответствует увеличению (уменьшению) торможения, при λ₁ = 2 торможение отсутствует. На вид эпюр скоростей оказывают также влияние глубины, градиенты глубин и величина скорости потока, учитываемые с помощью λ₂. При возможности следует подбирать начальные значения λ₁ и λ₂, добиваясь совпадения с эпюрами, построенными по натурным данным. При отсутствии наблюдений можно использовать значения «по умолчанию» ( λ₁=1,85; λ₂=2).

Вид эпюр скоростей, рассчитываемых по (6)-(8), в пределах точности расчётов соответствует эпюрам по данным натурных наблюдений и рассчитываемым по методу Г.П. Скребкова, которые, "по форме и существу [25]" соответствует уравнению А.В. Караушева [8].

Далее выполняются расчёты векторов скоростей основного потока в клетках по всей длине русла. Не вдаваясь в подробности опишем лишь основные этапы алгоритма:

– определяются характерные точки линии максимальных глубин (тальвега) ;

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

– рассчитываются расстояния характерных точек тальвега от начала участка водотока и упорядочивается массив этих точек по возрастанию расстояния;

– последовательно определяются очередные две точки тальвега, удаляясь от начала водотока в сторону наиболее удалённой точки, строятся в этих точках поперечные створы и затем - полигоны, ограниченные створами и берегами водотока (в начале и конце водотока - створами, берегами и границами клеточной области), и по заданным расходу и восстановленным глубинам определяются в каждой клетке внутри очередного полигона скорости основного потока с помощью выражений (6)-(8).

Описанный выше алгоритм приводит к приемлемым результатам на примерах расчётов для Кубани и Дона и их притоков в нижнем течении. Удовлетворительные результаты показывает также сопоставление расчётных и исходных значений на участке реки Доммель (участок длиной около 250 м), расположенный вблизи границы Бельгии и Нидерландов (подробности - на сайте pvn.org.ru).

Полученные векторы_i,j и_i,j за счёт ветровой и стоковой составляющей складываются для получения результирующего вектора основного потока_i,j.

3.3 Крупномасштабная турбулентность. На основании обобщений из п. 2.1 можно получить формулу

(9)

где Ũ' - приближение для компоненты (м/мин) модуля скорости Ư_k;

n – число клеток в стороне квадрата клеточной области;

Ʌ – масштаб крупномасштабных колебаний (по умолчанию Ʌ = L/2);

k_Ʌ – коэффициент для учёта изменения интенсивности эффективной турбулентности при имеющейся площади S (м²) водного зеркала (k_Ʌ=L/S^1/2) и, дополнительно при наличии стока, - за счёт дисперсии основного потока при изменении ширины водотока (k_Ʌ=k_Ʌ·l_v/(S/l_v);

l_v – длина участка водотока по тальвегу (м);

K_ℓ(Ū) - коэффициент горизонтальной турбулентной диффузии (м²/мин), вычисляемый по формуле с использованием выводов работы [23]:

(10)

где τ - лагранжев временной масштаб (по умолчанию τ = 1.7 мин);

При выводе формулы (9) использован дискретный вариант для (3)

(3΄)

где Δσ² - приращение горизонтальной дисперсии за время Δt.

В формуле (9) учтены выводы из п. 2.1 о введении понятия “эффективной турбулентности” путём дополнительного учёта влияния дисперсии основного потока за счёт морфологических особенностей водного объекта введением величины, которую можно именовать “эффективным коэффициентом турбулентности”:

(3΄΄)

Векторы D2Q9-набора компонент скорости квазиплоского потока с учётом (9) и разложения по двум из восьми горизонтальных направлений ранее полученной скорости основного потока вычисляются с использованием векторного суммирования:

где w_k – весовые коэффициенты, перераспределяющие интенсивности горизонтальной и вертикальной турбулентности при изменении глубины:

h_Ekm – толщина квазиоднородного слоя (по умолчанию: h_Ekm=

Δh_ekm= h_Ekm/9; w_l,m=0 = 8, при h_i,j/Δh_ekm > 8) .

Горизонтальные компоненты (l,m ≠ 0) обеспечивают адвекцию и горизонтальную диффузию, а вертикальная (l,m = 0) - тормозящий эффект, обусловленный гидродинамическим сопротивлением при данной глубине.

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

Рассчитываются новые значения скоростей Û_k с учетом перепадов глубин:

(13)

где q_i+l,j+m — рассчитывается с помощью функции логистического типа:

(14)

где λ₁ ~ 1,8÷2,0; λ₂ ~ 2÷3;

ħ_i+l,j+m — относительный перепад глубин:

(15)

Коэффициент λ₁ имитирует торможение горизонтальных компонент. При λ₁ = 2 торможение отсутствует. Коэффициент λ₂ влияет на условные скорости в зависимости от относительных перепадов глубин вокруг центральной клетки. У береговой линии (при H_i+l,j+m = 0) происходит обнуление компонент скорости. Коэффициенты λ₁ и λ₂ подбираются эмпирически (по умолчанию λ₁ = 1,85; λ₂ = 2) или с использованием оптимизационной процедуры [21].

Условие сохранения баланса импульса проверяется с помощью величины:

(16)

где S - обозначает сумму импульсов (горизонтальных скоростей) соответственно для Û_i+l,j+m и Ũ_i+l,j+m, при l,m ≠ 0:

		(17)

При |ΔŜ| > 0 выполняется упрощённая процедура восстановления баланса применением правила:

(18)

где ω_i+l,j+m - весовые коэффициенты:

(19)

Правила (16)-(18) повторяются с подстановкой (18) в (16) до достижения значения ΔŜ = 0.

Обычно расчёты для ОВОС выполняются в условиях меженных расходов и слабого ветра, вызывающих невысокие скорости течения. Но и в этом случае появляются составляющие течений типа компенсационных, возникающие в сеточных двумерных моделях из-за избыточного давления при достаточно больших колебаниях уровня свободной поверхности [5, 14, 20]. Иногда скорости таких составляющих могут приводить к ситуации, в которой недостаточно применения (18) для восстановления баланса. Возникновение такой ситуации проверяется с помощью величины:

(20)

где V - обозначает сумму импульсов с учётом вертикальной составляющей скорости соответственно для Ǔ_i,j и Ũ_i,j:

(21)

где S - вычисляется по формуле (17).

При |ΔV| > 0 выполняется процедура восстановления баланса применением правила: U_i,j увеличивается или уменьшается до достижения значения ΔV = 0.

С применением последнего правила заканчивается получение набора компонент U_i+l,j+m, определяющего турбулентное течение в каждой клетке. Векторная сумма этих компонент даст скорость адвекции Ū. Вычитание из U_i+l,j+m проекций Ū_i+l,j+m полученного вектора на ближайшие горизонтальные направления даст значения компонент U'_i+l,j+m, ответственных за диффузию. Таким образом, можно считать решённой задачу получения разложения (1').

3.5 Об учёте граничных условий. В расчётной области могут присутствовать два типа границ: береговая линия и открытая граница. У береговой линии должно происходить отражение компонент течения, близкое к отражению по Френелю, а на открытой границе - свободное протекание. Для выполнения этих условий в обоих случаях достаточно формул (13) - (21).

3.6 Определение времени. Учёт времени требует дополнительных вычислительных затрат, поскольку в КА-модели используются скорости течения, размеры клеток и номера циклов, без явного упоминания времени. Ход времени в каждой клетке может различаться, то есть время является локальной характеристикой. Такое естественное воспроизведение асинхронности процессов в различных зонах клеточной области возникает при использовании вещественных переменных для описания состояния клеток в отличие от классических КА-моделей с ограниченным набором состояний [30].

Возможны различные варианты исчисления времени, единого для всей области. В настоящем варианте модели используется упрощённый алгоритм определения времени с усреднением по всей расчётной области, исходя из преобладания слабых течений в расчётах для ОВОС, когда основной поток можно считать однородным. Тогда можно определить интервал времени в клетке Δt(I)_i,j, для каждого I-го цикла

(22)

и единое время для всех клеток области

(23)

где I = 0 cоответствует первому выбросу взвеси из клетки-источника.

4 Моделирование загрязнения

4.1 Адвективная диффузия взвеси. Следует отметить, что в CAD2Q9-модели предусмотрена и программно реализована возможность манипуляций со взвесью, состоящей из частиц нескольких фракций (полидисперсная взвесь), когда перемещение частиц каждой фракции имитируется в виде отдельного монодисперсного шлейфа (облака), а для получения суммарного количества вещества используется принцип суперпозиции. Такой подход равнозначен расчётам с использованием так называемой «эффективной гидравлической крупности» [10], однако манипулирование отдельными облаками оставляет более широкие возможности для развития модели в случае необходимости учёта взаимозависимости облаков, как это происходит при флокуляции [32]. В дальнейшем для краткости изложения полидисперсность учитывается лишь при необходимости.

Далее рассматриваются правила (алгоритмы) клеточного механизма, с помощью которого имитириуется распространение взвеси при известном D2Q9-наборе компонент турбулентного потока U_i,j.

На рисунке 2 показаны схемы расположения основных величин, используемых в CAD2Q9-модели с индексами клеток D2Q9-окрестности для случая, когда центральная клетка находится в клетке (1,1) клеточной области.

а)	б)	в)
Рисунок 2. Схемы расположения величин в CAD2Q9-модели: а) глубины, компоненты скорости течения и скорости осаждения взвеси; б) весовые коэффициенты перемещения и осаждения взвеси; в) масса взвеси и масса образовавшегося осадка

 

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

На рисунке 2 представлены следующие величины: U_i+l,j+m – компоненты скорости турбулентного течения (см. рисунок 1-в); Usub_i+l,j+m - скорости осаждения взвеси; H_i+l,j+m – глубина; W_i+l,j+m - весовые коэффициенты для учёта влияния компонент скорости на горизонтальное перемещение массы взвеси; Wsub_i+l,j+m - весовые коэффициенты для учёта влияния скорости осаждения на массы взвеси и образующегося осадка; Msup_i+l,j+m(I) - масса взвеси на I-том цикле; Msub_i+l,j+m(I) - масса образовавшегося осадка на I-том цикле.

Клеточный автомат представляет собой дискретную модель процесса распространения взвеси, развивающегося в каждой клетке двумерной клеточной области с окрестностью Мура [30] на очередном цикле I. Используется КА первого порядка, когда значения в клетках очередного цикла (I+1) зависят от значений лишь одного предыдущего I-го цикла. Количество взвешенного вещества М_i,j в клетке с координатами (i, j) рассчитывается по одному и тому же алгоритму (КА-правилу), зависящяму от компонент скорости турбулентного течения U_i+l,j+m, глубин H_i+l,j+m и скоростей осаждения взвеси Usub_i+l,j+m в клетках окрестности.

Пусть схема на рисунке 2 соответствует ситуации на цикле I. Тогда КА-правила для расчета масс взвешенного и осаждённого вещества в центральных клетках окрестности можно представить следующим образом:

где i = 0, 1, 2, ..., n-1; j = 0, 1, 2, ..., n-1; l = -1, 0, 1; m = -1, 0, 1;

Used(I)_i+l,j+m - скорость осаждения взвеси (формула (42), п.3.4.2);

Msrc_isrc,jsrc ‑ взвесь, поступающая из источника в клетке с индексами i = isrs, j = jsrc в течение некоторого числа циклов I = 1,2,… ;

Psrc_isrc,jsrc ‑ мощность источника (масса взвеси в единицу времени);

γ ‑ коэффициент, исправляющий нарушения симметрии, свойственное окрестности Мура [13, 31] (γ = 2^-1/2, k = 0, 2, 6, 8);

Ū_i,j – модуль вектора скорости основного потока;

∆r_i,j – расстояние, проходимое в клетке в направлении вектора скорости основного потока.

Осаждённая часть вещества, накопленная в клетке за все циклы и используемая далее для расчета толщины образовавшегося осадка, вычисляется суммированием по всем циклам:

(29)

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

– вещество, поступившее из источника к текущему циклу

(30)

– вещество, содержащееся в облаке на текущем цикле

(31)

– вещество, выпавшее в осадок за все циклы

(32)

– вещество, вынесенное за пределы области на открытых границах за все циклы

(33)

Условие сохранения вещества можно записать в виде:

(34)

Необходимо отметить, что при больших перепадах глубин условие (34) может нарушаться, как это происходит и в двумерных сеточных моделях [14, 20]. При возникновении такой ситуации на очередном цикле потеря или излишек вещества восполняется или изымается пропорционально количеству вещества в клетках. С помощью таких поправок условие (34) соблюдается с заданной точностью (по умолчанию с точностью 5 %).

Концентрация взвеси C (мг/л) и толщина осадка Z (мм) рассчитываются с использованием упрощённых методов, соответствующих задачам, решаемым с помощью CAD2Q9-модели. После окончания очередного цикла:

	(35)
	(36)

где ρ_d - природная плотность осадка (г/см³ ).

4.2 Образование осадка. Здесь приводится пояснение расчёта скорости осаждения Used(I)_i+l,j+m, используемой выше по тексту в формуле (25). При этом следует вспомнить, что в настоящем варианте модели учитываются особенности полидисперсности взвеси, содержащей частицы нескольких фракций (fr =1,..,N) с различающимися размерами d_fr, когда возможна флокуляция частиц.

В CAD2Q9-модели при вычислении скорости осаждения Used_fr(I)_i,j учитываются два типа осаждения мелкодисперсных частиц в зависимости от их размера d_fr и весовой концентрации частиц фракции C_fr(I)_i,j: свободное осаждение и стеснённое осаждение. Свободное осаждение наблюдается при низкой концентрации частиц, когда можно считать, что частицы оседают независимо друг от друга, не сталкиваясь и не взамодействуя друг с другом. В этом случае скорость осаждения частиц размером 0,05 мм и менее рассчитывается непосредственно по известной формуле Стокса.

Для учёта стеснённости, проявляющейся при повышенных концентрациях мелких вязких частиц (менее 0,05 мм) используется подход [32], при котором к скорости осаждения по Стоксу вводятся поправки на затруднённое осаждение и флокуляцию.

Алгоритм расчёта скорости осаждения приводится в приложении А.

5 Реализация и апробация

Программная реализация КА-модели представляет собой Javascript-версию, имеющую наименование CAD2Q9-SIST (Cellular Automata D2Q9-Simulation Suspension Transfer, клеточно-автоматная D2Q9-имитация переноса взвеси) [22]. Имеется два варианта реализации: в среде браузера и десктопный вариант. Второй вариант использует платформу nw.js [36] и обладает преимуществом за счёт возможности работать с системными API операционных систем Linux, Windows и Mac OS.

Важным преимуществом по сравнению с программами, основанными на сеточных моделях, являются пониженные требования к машинным ресурсам, что, наряду со снижением трудоемкости за счет особенностей КА-модели, обеспечивается использованием веб-технологий в среде клиента (имеются ввиду современные Mozilla Firefox, браузеры на основе Chromium: Chrome, Яндекс.Браузер и MS Edge, а также десктопный вариант с применением nw.js). Для разработки использован язык Javascript в среде HTML4 с элементами HTML5 [11]. Построение и отображение интерактивных графиков осуществляется с помощью некоторых функций Javascript-библиотеки Highсharts JS [35]. Для сохранения и последующей интерпретации результатов используется клиентская база данныхна основе HTML5 API IndexedDB [34].

В CAD2Q9-SIST используется ещё одно преимущество, состоящее в использовании многоядерности (многопоточности) современных процессоров за счёт возможности запуска программы в нескольких окнах браузера (платформы nw.js). Таким образом можно, например, имитировать распараллеливание алгоритма, запуская распространение каждого монодисперсного облака подисперсной взвеси в одельном окне, реализуя принцип суперпозиции через базу данных посредством API IndexedDB.

Краткое описание работы программы приводится в приложении Б.

Сравнение с аналитическими решениями (4) в рамках ГТД-модели (2) проводилось при отсутствии и наличии адвекции. Распределение концентрации взвеси в КА-модели не трудно воспроизводить практически не отличающимся от получаемого по (4), изменяя параметры, принятые “по умолчанию”, для различных глубин и разрешений. С уменьшением глубины облако становится более плоским, так как сказывается усиление горизонтальной диффузионной составляющей скорости за счет уменьшения вертикальной. При уменьшении разрешения облако закономерно становится более плоским за счет усреднения. Впрочем, следует признать, что соответствие между моделями получить не трудно, поскольку в ГТД-модели допускается изменение параметров Ʌ и K_ℓ(t) в достаточно широких пределах [19].

На многочисленных примерах моделирования показано, что CAD2Q9-SIST воспроизводит отклонения горизонтальной турбулентной диффузии от изотропии, присущие прибрежным зонам [6, 19], что позволяет получать реалистичные характеристики распределения взвеси в таких зонах.

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

Заключение

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

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

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

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

В программе необходимо задавать константы, cвязанные с параметрами турбулентности. Удовлетворительные результаты можно получать, используя значения, заданные по умолчанию. Однако желательны уточнения по данным гидрологических наблюдений. Традиционные наземные и дистанционные измерения трудоёмки, недостаточно точны и подробны [9]. Необходимо использовать преимущества применения малых беспилотных водных и летательных аппаратов (БПВА [26] и БПЛА [28]), оснащённых соответствующими измерительными устройствами.

Представляется важным обратить внимание на слабую изученность вероятностных свойств турбулентности в мелководных прибрежных зонах. Проблему можно решать новыми методами с использованием недорогих дрифтеров, оснащённых GPS-приёмником и GSM-связью [24]. CAD2Q9-SIST можно использовать для планирования и нтерпретации наблюдений.

 

Приложение А

скорость осаждения взвеси

В CAD2Q9-модели учитывается два типа осаждения мелкодисперсных частиц в зависимости от их размера и концентрации: свободное осаждение и стеснённое осаждение. В первом случае скорость осаждения частиц размером 0,05 мм и менее рассчитывается непосредственно по формуле Стокса

(А.1)

где ω_fr,0 - скорость осаждения частиц фракции fr (м/с); ρ_fr,s - плотность частиц (кг/м³); ρ_w - плотность воды (кг/м³); g - ускорение силы тяжести (9,80665 м/с²); d_fr - диаметр частиц (м); ν – кинематическая вязкость воды (м²/с).

Для определения плотности воды используется эмпирическая завимость от температуры(T,°C) и солёности (S,‰) [7].

(А.2)

где b₀ = 8,24493·10^-1; b₁ = -4,0899·10^-3; b₂ = 7,6438·10^-5; b₃ = -8,2467·10^-7; b₄ = 5,3875·10^-9; c₀ = -5,72466·10^-3; c₁ = 1,0227·10^-4; c₂ = -1,6546·10^-6; d₀= 4,8314·10^-4;
ρ_T - плотность чистой воды:

(А.3)

где a₀ = 999,842594; a₁ = 6,793952·10^-2; a₂ = -9,095290·10^-3; a₃ = 1,001685·10^-4;
a₄ = -1,120083·10^-6; a₅ = 6,536332·10^-9.

Кинематическая вязкость воды (м²/с) равна

(А.4)

где μ - динамическая вязкость (Па·с=кг/(м·с) пресной воды, рассчитываемая по формуле Пуазейля:

(А.5)

Стеснённое осаждение проявляется при повышенных (Ο(10¹мг/л)) концентрациях мелких частиц. Следуя [32] используется следующий подход для расчёта осаждения частиц мелких (менее 0,05 мм) вязких фракций. Формула для расчёта скорости стеснённого осаждения имеет вид:

(А.6)

где w_fr,0 – скорость осаждения частицы мелкой фракции по Стоксу (А.1);

φ_hs, φ_floc – поправки на затруднённое опускание и флокуляцию.

Расчёт φ_hs выполняется с использованием эмпирической формулы:

(А.7)

где c_fr, c_gel - объёмные концентрации частиц фракции и ила.

Объёмная концентрация фракции вычисляется по известной массовой концентрации (35) и заданной плотности частиц ρ_fr (т/м³)

(А.8)

Осадок в виде ила (C_gel~130-1722 кг/м³) образуется из наиболее мелких фракций песка (0,05 мм) и вязких пылевидных и глинистых фракций (менее 0,05 мм) на заключительной стадии затруднённого осаждения, когда только что образовавшийся ил имеет максимально возможную пористость (η) в самом начале процесса уплотнения. В настоящем варианте программы по умолчанию c_gel ≈ 0,65.

Ускорение осаждения из-за флоккуляции рассчитывается с использованием эмпирической формулы

где d_sand – диаметр наиболее мелкой фракции песка (0,05 мм); d_fr – диаметр частиц фракций размером менее 0,05 мм; w_fr,0 – скорость осаждения отдельных частиц в чистой воде по Стоксу (А.1); C_fr – массовая концентрация частиц взвеси (C_fr=c_fr·ρ_fr,s); C_gel – массовая концентрация свежеобразовавшегося ила (C_gel=c_gel·ρ_gel); c_gel – объёмная концентрация ила.

Для расчёта концентрации ила в настоящем варианте модели используется упрощённая формула

где ρ_sand,d - массовая плотность сухих частиц (скелета) песка размером 0,05 мм;
k_cons – коэффициент консолидации, учитывающий разрыхлённость донного осадка (по умолчанию k_cons = 1,08).

Приложение Б

Описание работы с программой CAD2Q9-SIST

Описание работы с программой приводится с использованием результатов моделирования загрязнения при разведочном бурении на двух скважинах в Тузлинской промоине Керченского пролива. Применяются значения параметров “по умолчанию”, приведённые при описании КА-модели. Входные данные представлены в полях формы на рисунке Б.1 и близки к данным, упоминаемым на общественых слушаниях об исследованиях грунта перед предполагаемым строительством мостового перехода.

Для бурения используется буровая платформа типа “Ирбен”. В первые сутки проводятся работы на скважине №1, во вторые – на скважине №2. Взвесь вбрасывается в воду при извлечении обсадной трубы. Продолжительность вбросов, производительность по взвеси, характеристики фракций рассчитаны по исходным данным:

– глубина каждой скважины ‑ 70 м;

– скорость извлечения трубы ‑ 9 м/ч;

– диаметр трубы ‑0,146 м;

– толщина грунта, налипающего на трубу, ‑0,01 м;

– плотность грунта: ρ =1,82 т/м³; ρ_s =2,74 т/м³; ρ_d =1,32 т/м³;

– суммарная доля мелкодиспесных фракций (0,05 мм и менее) в грунте – 50%;

– фракции состоят из невязких частиц с размерами (мм): 0,05; 0,005; 0,001;

– соотношение долей (%) мелкодиспесных фракций: 40:40:20;

– гидравлические крупности частиц в программе рассчитывается по алгоритму, изложенному при описании модели, при: ρ_s=2,74 т/м³; t_w= 15°C; S = 12‰;

– плотность образующегося осадка принята равной 1,32 т/м³ при коэффициенте разрыхления 1,08 [8];

- ветер - северного направления со скоростью 4,5 м/с (по многолетним данным на ГМС “Опасное”), а сток направлен из Азовского моря в Чёрное с расчётным значением расхода воды через Тузлинскую промоину в размере 60 м³/с (при расчёте расхода использовались результаты работы [31]).

В верхней части формы на рисунке Б.1 предусмотрено поле для выбора наименования водного объекта с отображением карты-схемы объекта в поле слева. Местоположение квадрата (клеточной области) на карте-схеме можно выбирать перемещением мыши. Задавая различное количество клеток, можно изменять разрешение клеточной области. Координаты клетки с выбросами задаются в пределах перемещаемого квадрата.

На рисунке Б.1 показано также горизонтальное меню с пунктами:

- «Глубины» — загружается батиметрия, выполняется интерполяция и отображение поля глубин (предварительная оцифровка батиметрии района осуществляется с использованием вспомогательной Javascript-программы, входящей в состав комплекса);

- «Створы» — прокладываются продольный и поперечный створы;

- «Тальвег» — определяются критические точки тальвега;

- «Течения» — моделируются и отображаются скорости течений;

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

- «Глубины с илом» — иногда, при большом количестве осадка и малых глубинах, требуется перерассчитать глубины с учётом образовавшегося осадка для отображения и анализа и, возможно, перерассчитать течения и продолжить моделирование с изменившимися полями глубин и течений;

- «Открыть БД», «Добавить ил в БД», «Удалить ил из БД», «Сложить илы из БД», «Протекло концентрации - в БД», «Площ. с протёкш. к-циями -> в БД», «Площ. с протёкш. к-циями <- БД», «Сложить площади с к-циями из БД» —используются для получения характеристик образовавшегося осадка и площади под заданными величинами интегральной (“протёкшей” [16]) концентрации по результатам, сохраняемым в базе данных;

На рисунке Б.2 показан результат выполнения пункта меню «Глубины» с цветовым отображением карты глубин. На карте видны также линии продольных и поперечных створов, построенных с помощью кнопки «Створы».

Глубины интерполируются с использованием оптимизированного под поставленную задачу алгоритма, известного как метод обратных взвешенных расстояний (ОВР или IDW - inverse distance weighting). Входная батиметрия, представленная в GEOJSON-формате, готовится с использованием вспомогательной программы, входящей в состав комплекса CAD2Q9-SIST. Длина стороны квадрата составляет 1505 м, размер клетки - 3,76 м.

На рисунке Б.3 показан фрагмент интерактивной карты векторов адвекции, получаемой после нажатия кнопки «Течения».

 

Средняя скорость адвекции при усреднении по всей клеточной области составляет 0,076 м/с. Максимальная скорость наблюдается у берегов и достигает 0,2 м/с.

На рисунках Б.4-Б.13 приведены фрагменты карт и графиков, отображаемых в процессе моделирования, с помощью которых можно представить особенности пространственно-временной динамики распространения взвеси. Длительность одного расчётного цикла (шаг по времени) составляет 61 с.

 

 

 

 

 

 

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

1. Альтман Э.Н. Особенности горизонтальной турбулентности в мелководном проливе. /Э.Н. Альтман, Б.Л. Лагутин, Д.М. Толмазин // Басс. ГМО Черного и Азовского морей, вып. 4. Л.: Гидрометеоиздат, 1966. ‑С.75-87.
2. Барышников Н.Б. Гидравлическое сопротивление речных русел. - СПб, изд. РГГМУ, 2003. 147 с.
3. Боровков В.С. Русловые процессы и динамика речных потоков на урбанизированных территориях. Л.: Гидрометеоиздат. 1989. 286 c.
4. Волынов М.А. Влияние плановой геометрии речного русла на диффузию и дисперсию примеси // Фундаментальные исследования. – 2013. – № 6-3. – С. 535-540.
5. Вольцингер Н.E., Пясковский Р.В. Теория мелкой воды. Океанологические задачи и численные методы. Л.: Гидрометеоизат, 1977. 207 с.
6. Грузинов В.М. Прикладная океанография / В.М. Грузинов, Е.Б. Борисов, А.В. Григорьев // - Обнинск: изд-во «Артифекс», 2012. - 384 с.
7. ГСССД 76-84 Таблицы стандартных справочных данных. Морская вода // - М.: Изд-во стандартов, 1986.
1. 8. Караушев А.В. Теория и методы расчёта речных наносов // Л.: Гидрометеоиздат, 1977. -272 с.
8. Копалиани З.Д, Жук М.М. О перспективах создания методов оценки гидрологических и гидравлических характеристик неизученных рек на основе гидроморфологических зависимостей // - СПб.: РГГМУ, № 5. 2007, С. 86-97.
9. Котеров В.Н. Моделирование переноса взвешенных веществ на океаническом шельфе. Горизонтальное рассеяние. / В.Н. Котеров, Ю.С. Юрезанская // Журн. выч. матем. и мат. физ. ‑ Т. 50, № 2. 2010. - С. 375-387.
10. Лабберс П. HTML5 для профессионалов / П. Лабберс, Б. Олберс., Ф. Салим // ‑М.: ООО «И.Д. Вильямс», 2011. ‑ 272 с.
11. Липатов И.В. Оценка гидродинамических условий при ликвидации разливов нефти / И.В. Липатов, А.Е. Пластинин // Вестн. ГУ им. Макарова, вып. 5. - С 127-134
12. Лобанов А.И. Модели клеточных автоматов // Компьютерные исследования и моделирование. - 2010. - Т. 2, № 3. - С. 273-293.
13. Лятхер В.М., Прудовский А.М. Гидравлическое моделирование // - М.: Энергоатомиздат, 1984. - 392 с.
14. Мачин Д. А., Четверушкин Б. Н. Кинетические и lattice Boltzmann схемы // Математическое моделирование. — 2004. — Т. 16, No 3. — C. 87–94.
15. Методика исчисления размера вреда, причинённого водным биологическим ресурсам. Утв. приказом Минсельхоза РФ от 31.03.2020 №167.
16. Новиков Е.А. О турбулентной диффузии в потоке с поперечным градиентом скорости // Прикл. мат. и мех., 1958, т. 22, вып. 3, с. 576-579.
17. Озмидов Р.В. О зависимости коэффициента горизонтальной турбулентной диффузии в ветровых течениях от относительной глубины водоёма. - Изв. АН СССР, сер. геофиз., 1959, №8, с. 164-181.
18. Озмидов Р.В. Диффузия примесей в океане. ‑ Л.: Гидрометеоиздат, 1991. - 279 с.
19. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. Пер. с англ. - М.: Энергоатомиздат, 1984. - 152 с.
20. Полупанов В.Н. Алгоритм случайного поиска глобального минимума целевой функции // КОНеМ - 2012: матерiали II всеукраїнського наукового семiнару (м. Полтава, 7-8 вересня 2012 р.). - 84 с.
21. Полупанов В.Н. Программный комплекс «CAD2Q9-SIST». Свид-во о рег. программы для ЭВМ RU 2021611187. Патентное ведомство: Россия, 20
22. Пухтяр Л.Д. Турбулентные характеристики прибрежной зоны моря. / Л.Д. Пухтяр, Ю.С. Осипов. // Труды ГОИН. — 1981. — Вып.158. — С. 36-41.
23. Сильвестрова К.П. Возможности использования GPS-дрифтеров для исследования течений на шельфе Чёрного моря / К.П. Сильвестрова, С.А. Мысленков, А.Г. Зацепин, Е.В. Краюшкин, В.И. Баранов, Т.Е. Самсонов, С.Б. Куклев // Океанология, 2016, Т. 56, № 1. С. 159-166.
24. Скребков Г.П. Метод расчёта распределения скоростей по ширине слабодеформированного участка реки// Метеорология и гидрология, №3,1972. - С. 75-81.
25. Смуров А.Е Применение технологий и оборудования беспилотных водных аппаратов в картографировании и моделировании / А.Е. Смуров, С.А. Тесленок // [Электронный ресурс]. - URL: http://journal.mrsu.ru/earth/geodeziya-kartografiya-i-geoinformatika (дата обращения 26.05.2021).
26. Судольский Р.В. Динамические явления в водоемах // Л.: Гидрометеоиздат, 1991. — 263 с.
27. Федоровский А.С. Новые возможности исследования рек с помощью БПЛА / А.С. Федоровский, Дальневосточное отдел. РАН // Научное издание “Физика геосфер”, вып. 2, Владивосток, 2020. С. 177-191.
28. Толмазин Д.М. Проблемы динамики вод северо-западной части Черного моря / Д.М. Толмазин, В.А. Шнайдман, Ж.М. Ациховская. — К.: Наукова думка, 1969. — 130 с.
29. Тоффоли Т. Машины клеточных автоматов: пер. с англ. / Т. Тоффоли, Н. Марголус. — М.: Мир, 1991. — 280 с.
30. Шапиро Н.Б. К теории течений в Керченском проливе // Научн. основы компл. использов. прир. ресурсов шельфа. Севастополь, МГИ НАН Украины. - 2005. - С. 320-332.
31. Van Rijn, L.C. Principles of sediment transport in rivers, estuaries and coastal seas / Leo C. van Rijn – Amsterdam: Aqua Publication – I11. 1993, p.690.
32. Wolf-Gladrow D. A. A lattice Boltzmann equation for diffusion // Journal of Statistical Physics. — 1995. — Vol. 79, No 5–6. — P. 1023–1032.
33. Использование IndexedDB [Электронный ресурс]. ‑ URL: https://developer.mozilla.org/ru/docs/Web/API/IndexedDB_API/Using_IndexedDB (дата обращения 26.05.2021).
34. HIGHCHARTS [Электронный ресурс]. ‑ URL: http://highcharts.com (дата обращения 26.05.2021).
35. NW.js [Электронный ресурс]. - URL: https://nwjs.io/ (дата обращения 26.05.2021).