Численные модели для исследования двумерных ветровых течений в замкнутых водоемах
П.В. Белолипецкий
Красноярский государственный университет
E-mail:
belolip@icm.krasn.ruПостановка задачи Численный алгоритм
Введение
Сложность процессов и явлений гидродинамики требует проведения как натурных экспериментов, так и математического моделирования. Математическое моделирование позволяет получить дополнительную информацию, которая поясняет результаты экспериментов, а в некоторых случаях является единственно доступной. Использование математического моделирования позволяет реализовать идею виртуального испытательного стенда.
Данная работа посвящена разработке программного комплекса на основе методов вычислительной гидродинамики для исследования гидрофизических процессов в замкнутых водоемах. Математическая модель основана на двумерных уравнениях движения стратифицированной жидкости в приближениях Буссинеска и твердой крышки, уравнениях конвекции-диффузии для температуры и солености воды.
Для построения численного алгоритма применяются методы расщепления по физическим процессам и проекционный, в результате которых решение системы уравнений в частных производных сводится к последовательному решению отдельных уравнений. Получаемые уравнения сводятся методами конечных элементов и конечных разностей к системам линейных алгебраических уравнений. Матрицы этих систем являются разреженными и хранятся в компактной форме. Реализован ряд итерационных методов решения систем, в том числе и с помощью параллельных вычислений.Написаны следующие подпрограммы на языке С++: подпрограмма дающая изображение исследуемых величин в виде графиков, обеспечивающая сохранение, восстановление, управление параметрами и выполнение других полезных функций; подпрограмма генерации прямоугольной сетки со сгущением узлов в выделенных областях; подпрограмма решающая системы линейных алгебраических уравнений(в том числе и параллельная версия). Подпрограмма с помощью предыдущих подпрограмм реализующая численный алгоритм решения поставленной задачи. Проведены тестовые расчеты.
В дальнейшем планируется: добавить к модели уравнения, описывающие поведение некоторой биоты, распространить рассмотренные численные методы на решение трехмерных задач и на другие приложения вычислительной гидродинамики.
Обычно течения в замкнутых водоемах формируются под действием силы ветра. Течения оказывают влияния на распределения физических, химических и биологических компонент, которые в свою очередь влияют на течения. Рассматривается численный алгоритм для исследования двумерных в вертикальной плоскости ветровых течений и распределений компонент в замкнутых стратифицированных водоемах. Математическая модель основана на уравнениях турбулентных течений неоднородной жидкости в приближениях твердой крышки и Буссинеска с переменным коэффициентом вертикального турбулентного обмена. Основные уравнения [1]:
(1)
(2)
Здесь
– составляющие скорости течения воды в направлениях Оx и Оz,
– время,
- плотность воды,
- характерная плотность воды,
– давление,
– ускорение свободного падения,
– температура воды,
– соленость воды,
- коэффициенты турбулентного обмена.
Уравнения дополняются начальными и граничными условиями. Граничные условия на поверхности (z=0)
(3)
на дне (z=H)
(4)
Здесь
- внешняя нормаль к поверхности,
– напряжение трения ветра,
– полный поток тепла через свободную поверхность,
– теплообмен с ложем водоема,
– поток соли через свободную поверхность,
– массообмен c ложем водоема.
Для оценки коэффициента вертикального турбулентного обмена применяется формула Прандтля-Обухова [1].
Численный алгоритм основан на методе расщепления по физическим процессам и состоит из двух частей [2,7]. Поля скоростей и давления находятся из решения задачи (1), (3), (4) по известным значениям плотности. Затем из решения задачи (2), (3), (4) определяются температура, соленость и плотность воды. Пользуясь методами расщепления и проекционным методом задачу можно свести к последовательному решению уравнений эллиптического, параболического и гиперболического типов.
Эллиптические и параболические уравнения решаются с помощью методов Галеркина и конечных элементов [6,7].
В соответствии с методом Галеркина неизвестные функции представляются в виде разложения по первым элементам некоторого базиса. Коэффициенты при базисных функциях подлежат определению. Подставляя разложения вместо неизвестных функций в уравнения, получаем некоторую невязку, зависящую от неизвестных коэффициентов. Приравнивая нулю интегралы по области, в которой ищется решение задачи, от произведения невязки на базисные функции, получаем систему линейных алгебраических уравнений для определения коэффициентов. Таким образом минимизируется невязка.
Сечение водоема приближенно представляется в виде прямой свободной поверхности и кусочно-линейного дна, на котором отрезки соединяют известные значения глубины. На эту схематизацию накладывается прямоугольная сетка, из которой посредством выкидывания элементов, полностью не лежащих или большая часть которых не лежит внутри области, получается сетка для расчетов. Написана подпрограмма, реализующая генерацию сетки со сгущением в нужных областях, т.е. в областях с большими градиентами.
В соответствии с методом конечных элементов с каждым узлом сетки связывается одна базисная билинейная функция, равная единице в этом узле и нулю на всех элементах кроме тех, которые содержат данный узел. Поэтому на каждом четырехугольном элементе не равными нулю являются четыре базисные функции, которые в системе координат, связанной с данным элементом , имеют следующий вид:
Здесь
- размеры элемента,
- локальные координаты.
Затем применяем формулу Грина
к членам уравнений, содержащим вторые производные. Цель этого шага - уменьшить порядок производной в уравнении, что является стандартным приемом в методе конечных элементов. Это необходимо, т.к. базисные функции линейны по координатным направлениям.
Граничные условия учитываются следующим образом. Если на части границы поставлено условие Неймана, то в соответствующей части интеграла по границе в уравнениях производная функции по нормали заменяется на заданное значение. Если на части границы поставлено условие Дирихле, то соответствующая часть интеграла по границе в уравнениях не учитывается. Зато при формировании системы уравнений в строках, соответствующих неизвестным, лежащим на этой границе, все элементы кроме диагонального зануляются, диагональный берется равным единице, а в правую часть подставляется известное из условий Дирихле значение.
Так получаются системы линейных алгебраических уравнений относительно неизвестных коэффициентов. Матрицы полученных систем линейных уравнений являются сильно разреженными и поэтому хранятся в компактной форме [4]. Хранится только следующая информация - ненулевые элементы матрицы, номера их столбцов, указатели на первые элементы каждой строки, для каждого элемента указатель на следующий элемент в строке.
При таком виде хранения системы можно решать итерационными методами, такими как простой итерации, Якоби, Зейделя, верхней релаксации, наискорейшего спуска, минимальных невязок [3,5]. Все эти методы были реализованы двумя способами на основе массивов и на основе списков. Причем для всех кроме методов Зейделя и верхней релаксации написана также параллельная версия на основе массивов. Механизм распараллеливания такой: основной операцией распараллеливаемых методов является умножение матрицы на вектор. Оно-то и распараллеливается: вначале по всем процессам раскидываются матрица и вектор. Затем каждый процесс умножает свою часть строк матрицы на вектор и результаты собираются в главном процессе. Если требуется проведение следующей итерации, то всем процессам нужно передать обновленный вектор. Причем по процессам раскидывается не весь вектор, а только та его часть, на которую умножается соответствующая часть матрицы. Полноценного сравнения эффективности методов пока проведено не было и основные результаты были получены при использовании метода Зейделя для одного процессора на основе списков.Гиперболические уравнения возникают в связи с явлением конвективного переноса и имеют вид
Здесь
- переносимая величина,
- компоненты скорости.
Применение метода Галеркина с конечными элементами для уравнения конвекции-диффузии приводит к осциллирующим решениям [6]. Поэтому пользуются его модификацией - методом Петрова-Галеркина. В данной работе был выбран наиболее простой способ решения уравнений переноса - конечно-разностный метод разностей против потока:
здесь
при 
при 
при 
при 
Явная схема была выбрана потому, что ограничения на шаг по времени в решаемой задаче из условий устойчивости получаются несущественными по сравнению с ограничениями из соображений точности.
В нашей задаче для нахождения поля скоростей из задачи (1.1), (1.3), (1.4) применяется метод расщепления. Были запрограмированны следующие варианты метода. Во всех вариантах вначале отдельно рассчитывается конвективный перенос:
Затем для каждого варианта выполняются следующие этапы.
В первом варианте метода решаются уравнения:
Во втором:
В третьем:
В четвертом:
Расчеты для модельных задач показали, что наиболее эффективным является четвертый вариант.
Алгоритмы тестировались на модельных параболических
и эллиптических
уравнениях в прямоугольной области с точными решениями вида
и
соответственно. При тестировании в эти формулы подставлялись конкретные значения постоянных коэффициентов, а на границах ставились различные комбинации условий первого и второго рода.
Гиперболические уравнения тестировались на невязкой конвекции пассивной скалярной характеристики [6]. Решалось уравнение
с начальным условием
при 
при 
Начальная форма функции выглядит как конус. Ось конуса проходит через точку 
, а радиус основания равен 
.
Составляющие скорости предполагаются известными и задаются формулами
Точное решение соответствует вращению конуса относительно начала координат по круговым траекториям без искажения первоначальной формы. Искажение формы конуса может служить легко наблюдаемым качественным критерием точности конкретно выбранной вычислительной схемы.
Для тестирования схем расщепления использовались следующие два теста. Исходные уравнения (1) переписываются в переменных
- вихрь,
- функция тока.
Они связанны с исходными переменными соотношениями
Уравнения принимают вид:
,
.
Тест первый:
- выбирается так, чтобы точное решение в прямоугольной области 
имело следующий вид:
.
На границах задаются условия Неймана и Дирихле.
Тест второй:
установившиеся ветровые течения однородной жидкости в прямоугольной области 
без учета конвективных членов,
Здесь уравнения принимают вид:
,
,
Пусть 
тогда из второго уравнения 
Учитывая граничные условия (1.3), получаем при 
Откуда 
Следовательно 
Будем искать функцию 
в виде 
Здесь 
- решение задачи
,
при 
,
при 
- решение задачи
при 
при 
Получим
- коэффициенты разложения в ряд Фурье по синусам функции 
При 
.
В плоскости симметрии 
поэтому 
.
Численные решения в обоих случаях хорошо согласуются с построенными аналитическими решениями.
Разрабатываемый программный комплекс будет применяться для исследования гидрофизических и химико-биологических параметров замкнутых водоемов на примере озера Шира. Подпрограммы могут без изменений или с небольшими изменениями использоваться для решения различных задач.
Работа выполнена при финансовой поддержке ФЦП “Интеграция” проект № Э3-137. The research described in this publication was made possible in part by Aw
ard No. KY-002-X1 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF).