МОДЕЛИРОВАНИЕ СКОРОСТЕЙ ТЕЧЕНИЙ МЕЛКОВОДНЫХ ВОДОЕМОВ В ПАКЕТЕ FREEFEM++

MODELING THE CURRENT VELOCITY OF SHALLOW WATER IN FREEFEM++

Artem Sapozhnikov

student, Department of Mechanical Engineering, Marine Engineering and Transport, Far Eastern Federal University,

Russia, Vladivostok

Elena Amosova

professor, Far Eastern Federal University,

Russia, Vladivostok

АННОТАЦИЯ

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

ABSTRACT

The problem of finding the velocities of the currents of shallow water bodies, taking into account the determining physical factors, has been solved. A program written in the FreeFem ++ language allows you to simulate a two-dimensional field of water flow velocities on various geometries.

Ключевые слова: мелкая вода, метод конечных элементов, FreeFem++, математическое моделирование.

Keywords: shallow water, finite element method, FreeFem++, math modeling.

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

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

Все более активно модель мелкой воды находит свое применение для моделирования астрофизических систем [4, с. 32], [5, c. 32], метеорологических и климатических процессов в системе атмосфера-океан.

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

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

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

Для системы уравнений (1) устанавливаются следующие граничные условия:

1. Условие на свободной поверхности на входящей и выходящей границах:

Условие Дирихле для скорости и заданный уровень поверхности воды на входящей границе:

Условие «непроницаемой» стенки на выходящей границе:

2. Условие на твердых границах.

Условие «проскальзывания»:

Условие «прилипания»:

Также устанавливаются начальное условие для функции скорости и уровня поверхности воды в начальный момент времени:

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

Система уравнений (1) должна быть обезразмерена, поэтому после проведения такой процедуры получаем:

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

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

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

При добавлении краевого условия (2) и начального условия (6) задача для функции поверхности воды примет вид:

где  – граница входа, .

Согласно методу характеристик аппроксимация полной производной по времени происходит следующим образом:

Рассмотрим геометрию задачи. На входящей границе (рис. 1) задан вектор скоростей, который удовлетворяет краевому условию (2).

Рисунок 1. Заполнение области  в любой момент времени

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

Характеристики  определяются из задачи Коши:

где  – точка в области,  – вектор скорости потока.

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

Обозначим через  – время появления частицы с координатами  в область течения . Заметим, что  для тех частиц, которые при  заполняли область  для тех частиц, которые втекли в  через боковую поверхность  со скоростью .

Для точки  время появления частицы равно , в то время как для точки  время появление частицы равно .

Учитывая появление частичек в области не только в начальный момент времени, но и в другой момент времени с боковой поверхности, перепишем аппроксимацию (10):

где оператор конвекции:

В операторе конвекции  функция  удовлетворяет граничному условию (2).

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

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

где  – параметр регуляризации.

При добавлении  в первое уравнение системы (9) изменяется тип самого уравнения на параболический. Для корректности постановки задачи добавим начальное условие (6) и граничное условие для функции поверхности воды, тогда:

В задаче (14)  – решение, полученное методом характеристик.

Проведем дискретизацию по времени уравнения (14):

Запишем слабую формулировку поставленной задачи для функции свободной поверхности воды.

Пусть  – тестовая функция, такая что  Домножим уравнение (15) на тестовую функцию и проинтегрируем по области :

Используя формулу Остроградского-Гаусса к последнему интегралу и учитывая граничное условия для тестовой функции, а также начальное и граничное условия функции поверхности воды, получим слабую формулировку задачи:

где

Перейдем к рассмотрению системы уравнений Навье-Стокса. Проведем дискретизацию обоих уравнений:

где

Заметим, что в уравнениях присутствуют нелинейные слагаемые. Чтобы избавиться от нелинейности применим метод Ньютона:

Внеся  под знак корня и упростив получившееся выражение, получим:

Аналогично для второй компоненты:

Приступим к вариационной формулировке. Введем тестовые функции  и , такие что  Умножим уравнение Навье-Стокса для первой компоненты вектора скорости на тестовую функцию , а уравнение для второй компоненты – на , и, просуммировав получившиеся уравнения, проинтегрируем по области , учитывая начальные и граничные условия для функции скорости.

где

В слабой постановке (17) переменная  является решением, полученным при решении задачи (14).

Помимо полученных вариационных формулировок (14), (17) необходимо добавить задачи для нахождения функции дна и вектора скорости в начальный момент времени.

Начнем с функции дна. Ее нахождение обусловлено тем, что для лучшей сходимости решения функция дна должна быть гладкой, поэтому  должна удовлетворять уравнению Пуассона:

где  – произвольная функция множества, проекция дна.

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

На границе входа/выхода устанавливается условие «непроницаемой» стенки:

На твердых границах устанавливается условие Дирихле:

где  – произвольная функция или число.

Также после получения решения необходимо провести его нормировку, чтобы максимальная высота была равна 1:

Перейдем к задаче для вектора скорости в начальный момент времени , постановка которой обусловлена тем, что крайне сложно подобрать функцию скорости, которая бы удовлетворяла краевым условиям. Поэтому будем решать систему уравнений Навье-Стокса в стационарном случае с заданными краевыми условиями. Такая система уже рассматривалась ранее только для нестационарного случая, поэтому, чтобы не проводить вывод уравнений заново внесем некоторые изменения в вариационную постановку (17), а именно исключим интегралы, в которых присутствовали производные по времени, и функцию  заменим на функцию .

Результаты. Начнем рассмотрение получившихся графиков на прямоугольной геометрии. На рис. 3 и рис 4 представлены графики скоростей течений при коэффициенте придонного трения  и скорости ветра  м/с.

Рисунок 2. Функция дна акватории на прямоугольной геометрии

Рисунок 3. Скорость течения в начальный момент времени

Рисунок 4. Скорость течения в конечный момент времени

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

Теперь изменим функцию дна (рис. 7).

Рисунок 5. Функция дна акватории на прямоугольной геометрии

Задав вектор скорости ( м/c) и коэффициент придонного трения (), скорость течения получает вид:

Рисунок 6. Скорость течения в начальный момент времени

Рисунок 7. Скорость течения в конечный момент времени

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

Теперь изменим не только вектор скорости ветра, трение и функцию дна (рис. 14), а также геометрию (рис. 13).

Рисунок 8. Расширяющаяся геометрия

Рисунок 9. Функция дна на расширяющейся геометрии

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

На рис. 15 и рис. 16 представлен результат расчета, в решении которого скорость ветра равна 0 м/с, и трение на дне равно .

Проведя расчет, были получены следующие графики:

Рисунок 10. Скорость течения в начальный момент времени

Рисунок 11. Скорость течения в конечный момент времени

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

Перейдем к рассмотрению сужающейся геометрии (рис. 20).

Рисунок 12. Сужающаяся геометрия

Для данной геометрии представим следующую функцию дна (рис. 21):

Рисунок 13. Функция дна на сужающейся геометрии

Рассчитанное течение при скорости ветра равной 26.6 и придонном трении равном  представлено ниже (рис. 22, рис. 23):

Рисунок 14. Скорость течения в начальный момент времени

Рисунок 15. Скорость течения в конечный момент времени

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

Таким образом, все полученные результаты показывают работоспособность написанной программы.

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

2. На основе полученного алгоритма решения задачи, написана работоспособная программа в пакете FreeFem++, которая моделирует функцию дна, функцию свободной поверхности и скорость течения воды с учетом таких параметров как трение дна, трение поверхности воды, ветровые течения, турбулентный обмен.

3. Проведена оценка получившихся результатов.

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

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

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

1. Frederic Hecht in collaboration with Syl an Auliac, Oliver Pironneau, Jacques Morice, Antoine Le Hyaric, Kohji Ohtsuka, Pierre Jolivet – FreeFem++ Third Edition, Version 3.58-1
2. Georges Sadaka. FreeFem++, a tool to solve PDEs numerically. 2012.
3. Климович В.И., Воронков О.К., Давиденко В.М., Иванов С.Н. Численное моделирование приливов и распространения волн цунами // Инженерно-строительный журнал. 2018. № 6(82). С. 228–242.
4. Истомина М. А., Злотник А. А., Елизарова Т. Г. - статья «Гидродинамические аспекты формирования спирально-вихревых структур во вращающихся газовых дисках, "Астрономический журнал"», г. Москва, 2018.
5. «Shallow Water Magnetohydrodynamics in Plasma Astrophysics. Waves, Turbulence, and Zonal Flows» Arakel Petrosyan, Dmitry Klimachkov, Maria Fedotova and Timofey Zinyakov. 2020 by the authors. Licensee MDPI, Basel, Switzerland.
6. Удалов А. А. «Численное решение уравнения мелкой воды для ограниченного бассейна в пакете FreeFem++» // свидетельство о государственной регистрации программы для ЭВМ. 2021, Российская Федерация.
7. Сегерлинд Л. «Применение метода конечных элементов» — М.: Мир, 1979. — 392 С.
8. Адрин Гилл. Динамика атмосферы и океана: В 2-х т. Т. 1. Пер. с англ.: - М.: Мир, 1986. – 396 с., ил.
9. Адрин Гилл. Динамика атмосферы и океана: В 2-х т. Т. 2. Пер. с англ.: - М.: Мир, 1986. – 415 с., ил.
10. Н. Е. Вольцингер, К. А. Клеванный, Е. Н. Пелиновский. Длинноволновая динамика прибрежной зоны. Ленинград: Гидрометеоиздат, 1989.