Моделирование загрязнения мелководных зон

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

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

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

Метод. Использование численного дискретного метода на основе клеточных автоматов (cellular automata - CA или КА) [12], являющегося альтернативой традиционным сеточным методам интерпретации уравнений гидродинамики. Глобальности сеточных методов противопоставляется локальность КА, что более адекватно мелководной прибрежной зоне, где преобладает вес мелкомасштабной (морфометрической) турбулентности.

Кроме этого, можно использовать идеи метода дискретных решёточных уравнений Больцмана (lattice Boltzmann method, LBM) [6, 13]. Поскольку LBM можно считать разновидностью клеточных автоматов, то для будущей модели можно использовать логику, принятую для КА, а именно: модель должна имитировать перенос взвеси квазиплоскими течениями в расчётной КА-области с ячейками (клетками) квадратной формы (в терминах КА - окрестность Мура [12]) и с использованием набора скоростей D2Q9 (в терминах LBM - два измерения, девять соседей) [13], где D — означает пространственную размерность клеточной области, а Q — число связей в окрестности.

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

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

Перспективность подхода к реализации обусловлена тем, что КА соответствуют условиям так называемого «внутреннего параллелизма» [6], а современные браузеры поддерживают распараллеливание за счёт использования многоядерности (многопоточности) процессорами современных компьютеров [4].

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

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

Решение в каждой клетке расчётной области, по аналогии с LBM, может состоять из нескольких этапов: на первом получаем компоненты скоростей для гауссового радиально симметричного плоского облака [8], а затем - компоненты скорости, формирующие квазиплоское асимметричное облако, путём адаптации к граничным условиям и законам сохранения:

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

Задачу можно пояснить следующим образом.

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

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

При постоянной глубине поле скорости u соответствует нормальному закону. Если в этом поле действует мгновенный точечный источник с производительностью P по взвеси, то возникнет облако, концентрация взвеси С в котором представима в виде гауссовой функции. Положение центра x₀ этой функции прослеживается во времени t с помощью уравнения dx₀/dt=ū(x₀,t).

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

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

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

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

При постоянной глубине h и при перемещении взвеси на небольшое расстояние такое решение вполне можно было бы использовать на практике. Для определения адвективной составляющей имеются эмпирические уравнения типа ū=f(W,h) [10] и ū=f(Q,h) [2], где W - скорость ветра, Q - расход воды. Коэффициенты K_ℓ на мелководье можно определять, используя имеющиеся эмпирические зависимости типа K_ℓ=f(ū) [9, 11]. Величину σ можно определить из (1.3). Отсюда перемещение облака концентрации взвеси, порождаемого точечным источником, и осаждаемую долю получаем, задавая P и w. Многократный ввод P даёт возможность имитировать непрерывный источник.

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

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

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

Аналитического решения для (1.2') не существует, поэтому необходимо использовать численные методы.

Распространение взвешенных веществ в мелководных зонах в большинстве случаев имитируется с помощью двумерных сеточных моделей на основе уравнений “мелкой воды” - упрощения уравнений Навье-Стокса и Эйлера [1, 5], применяющиеся для водных объектов, горизонтальные размеры которых во много раз превышают глубину.

Однако использование таких моделей для прибрежных участков связано с трудностями из-за нерешённых проблем турбулентности [5] и чрезмерного сглаживания получаемых результатов [2,5]. Для традиционного решения этих проблем с использованием полуэмпирических параметризаций необходимы данные натурных исследований турбулентности в прибрежных зонах, которые, как правило, отсутствуют.

Проблемы сеточных моделей в нашем случае для ограниченного прибрежного участка можно попытаться “обойти” с помощью CAD2Q9-модели.

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

Вольцингер Н.E., Пясковский Р.В. Теория мелкой воды. Океанологические задачи и численные методы. Л.: Гидрометеоизат, 1977. 207 с.
Копалиани З.Д, Жук М.М. О перспективах создания методов оценки гидрологических и гидравлических характеристик неизученных рек на основе гидроморфологических зависимостей // - СПб.: РГГМУ, № 5. 2007, С. 86-97.
Котеров В.Н. Моделирование переноса взвешенных веществ на океаническом шельфе. Горизонтальное рассеяние. / В.Н. Котеров, Ю.С. Юрезанская // Журн. вычисл. матем. и матем. физ. ‑ 2010. ‑ Т. 50, № 2. ‑ С. 375-387.
Лабберс П. HTML5 для профессионалов / П. Лабберс, Б. Олберс., Ф. Салим // - М.: ООО «И.Д. Вильямс», 2011. - 272 с.
Лятхер В.М., Прудовский А.М. Гидравлическое моделирование // - М.: Энергоатомиздат, 1984. - 392 с.
Мачин Д. А., Четверушкин Б. Н. Кинетические и lattice Boltzmann схемы // Математическое моделирование. — 2004. — Т. 16, No 3. — C. 87–94.
Методика исчисления размера вреда, причинённого водным биологическим ресурсам. Утв. приказом Минсельхоза РФ от 31.03.2020 №167.
Озмидов Р.В. Диффузия примесей в океане. ‑ Л.: Гидрометеоиздат, 1991. - 279 с.
Пухтяр Л.Д. Турбулентные характеристики прибрежной зоны моря. / Л.Д. Пухтяр, Ю.С. Осипов. // Труды ГОИН. — 1981. — Вып.158. — С. 36-41.
Судольский Р.В. Динамические явления в водоемах // Л.: Гидрометеоиздат, 1991. — 263 с.
Толмазин Д.М. Проблемы динамики вод северо-западной части Черного моря / Д.М. Толмазин, В.А. Шнайдман, Ж.М. Ациховская. — К.: Наукова думка, 1969. — 130 с.
Тоффоли Т. Машины клеточных автоматов: пер. с англ. / Т. Тоффоли, Н. Марголус. — М.: Мир, 1991. — 280 с.
Wolf-Gladrow D. A. A lattice Boltzmann equation for diffusion // Journal of Statistical Physics. — 1995. — Vol. 79, No 5–6. — P. 1023–1032.