Перед моделированием движения жидкости ставится задача воспроизведения циркуляции жидкости в некой замкнутой или незамкнутой области, то есть нахождение модуля и направления скорости движения жидкости в каждой точке моделируемого объекта или системы. Существуют одномерные, двумерные и трехмерные модели, которые характеризуются количеством независимых координат, описывающих моделируемый объект или систему [2]. Наибольший практический интерес для воспроизведения циркуляции водных объектов, таких как озеро, водохранилище, река и другие, представляют трехмерные модели, которые могут быть реализованы в трехмерной прямоугольной (декартовой) системе координат.
Как правило, для задач, имеющих практическую ценность, при моделировании водных объектов или систем вода принимается как несжимаемаемая Ньютоновская жидкость. Такое допущение существенно упрощает уравнения Навье-Стокса, которые используют при моделировании движения жидкостей или газов.
Основная цель исследований – численно разрешить уравнения Навье-Сткоса в их нелинейном виде относительно простых физических величин с целью моделирования циркуляции жидкости для различных замкнутых и незамкнутых водных объектов, к которым можно отнести различные географические водные объекты.
Уравнения количества движения несжимаемой жидкости в консервативной форме в случае трехмерного пространства без участия внешних сил в безразмерных величинах можно записать в виде:
, (1)
, (2)
, (3)
где 
,
, 
– проекции неизвестной функции скорости 
на оси 
в декартовой прямоугольной трехмерной системе координат; 
– неизвестная функция давления; 
– число Рейнольдса.
Данные уравнения имеют в своем составе конвективные члены, представленные частными производными первого порядка по пространственным переменным, содержащимися в левой части. Эти члены уравнения учитывают перенос вихря, связанного с конвекцией и члены привносят в уравнение нелинейность. Также уравнение содержит диффузионные члены, представленные частными производными второго порядка по пространственным переменным в правой части уравнения, которые соответственно учитывают диффузионный перенос вихря.
Уравнение неразрывности для несжимаемой жидкости имеет вид:
. (4)
При моделировании жидкости систему уравнений (1)-(3) часто заменяют линеаризованными аналогами, выполнив замену мгновенной скорости на некую среднюю скорость потока в конвективных членах уравнения, тем самым упрощая данную систему. Среднюю скорость в данном случае можно считать постоянным коэффициентом, а, следовательно, уравнения становятся линейными. Такие линеаризованные аналоги имеют более простое численное решение. Такой упрощение было использовано при разработке математической модели Ладожского озера в работе [1].
Также в некоторых случаях используют упрощения, при которых нелинейные (конвективные) члены уравнения отбрасываются. В этом случае остаются только диффузионные члены уравнения и дальнейшее решение уравнений ведется без конвективных членов. Такое допущение также ведет к существенному упрощению численного решения этой системы уравнений. Однако модели, полученные таким образом, не всегда достаточно адекватно воспроизводят движение жидкости особенно для практических целей.
Произведем некоторые математические операции с исходными дифференциальными уравнениями (1)-(3), описывающими циркуляцию жидкости в трехмерном пространстве для получения формулы расчета неизвестной функции давления.
Продифференцируем систему уравнений по пространственным переменным: уравнение (1) – по 
, уравнение (2) – по 
, а уравнение (3) – по 
, и сложим полученные уравнения. Получим уравнение Пуассона [3], которое применяется для расчета функции давления в трехмерном пространстве:
. (5)
Приведем уравнение (5) к следующему виду, тождественно приравняв правую часть уравнения для простоты к члену 
, который в дальнейшем будет использован как «источниковый» член при нахождении сеточного аналога функции 
:
, (6)
где 
– дивергенция скорости; 
– оператор Лапласа.
Дивергенция скорости в принятых обозначениях имеет формулу:
. (7)
Согласно условию неразрывности дивергенция скорости равна нулю. В связи с этим можно было бы исключить члены, содержащие , из формулы (6) для расчета давления. Однако, поскольку для нахождения функции давления, уравнение (6) относительно будет решаться итерационными методами по указанным ниже причинам, то из-за недостаточной степени точности итерационного решения конечно-разностный аналог может быть не равен 0, то есть в узле сеточной области с номерами 
по осям трехмерной системы координат 
. Поэтому исключение из уравнения членов, содержащих , приведет к накоплению ошибки в уравнения количества движения, а также возможна и неустойчивость разностных аналогов данных уравнений в связи с накоплением этой ошибки, как указывается в работе [3]. Поэтому расчет членов, содержащих , несколько устраняет проблему. Выполняя расчеты с , имеется возможность приравнивать его сеточный аналог на следующем рассчитываемом временном шаге к нулю, тем самым не накапливая ошибку. В этом случае в конечно-разностном аналоге, аппроксимирующем производную дивергенции по времени необходимо принимать значение дивергенции на следующем временном слое равное нулю 
.
Для аппроксимации формулы (6) воспользуемся методом маркеров и ячеек, который был разработан авторами Харлоу и Уэлчем в 1965 году для решения дифференциальных уравнений численными методами. Описание метода представлено в работе [5].
Данный метод включают специфичную разностную сетку и специфичную структуру ячейки. Этот метод применяется для уравнений в простейших физических переменных, поэтому может быть применим для решения полученных уравнений. Трехмерная прямоугольная разностная сетка с позиционированием сеточных функций, используемая для аппроксимации непрерывных функций давления и компонентов скоростей, изображена на рисунке 1.
Рис. 1. Фрагмент разностная сетка для аппроксимации нелинейных членов уравнений
Особенность данной разностной сетки заключается в том, что сеточный аналог функции давления позиционируется в узлах трехмерной прямоугольной сетки (на рисунке такие узлы изображены крупными), а сеточные аналоги функций компонент скорости располагаются в дополнительных узлах (на рисунке они изображены мелкими), которые находятся посередине ребер трехмерной ячейки, образованной основными узлами сетки. Причем сеточный аналог проекции скорости на ось (
) располагается на ребре, параллельном оси , между узлами 
и 
. В свою очередь 
позиционируется на середине ребра, параллельного оси между узлами 
и 
, разделяя это ребро на две равных части, поэтому имеет номер узла 
по оси . Аналогичное позиционирование имеет и сеточный аналог проекции скорости на ось . Граничными узлами являются узлы, в котором позиционируется функция давления, значения скоростей задаются в этих узлах особенным образом, как будет показано ниже.
С учетом предложенной сетки конечно-разностный аналог дифференциального оператора 
, исходя из правил аппроксимирования «вперед» [4], имеет вид:
. (8)
Однако с учетом того, что , как было указано выше, для устранения невязки при итерационном решении получаемой системы алгебраических уравнений, конечно-разностный аналог в уравнении (6) примет вид:
. (9)
Конечно-разностные аналоги дифференциальных операторов 
,
, 
, 
, определяются следующими схемами, исходя из введенной сетки:
. (10)
. (11)
. (12)
. (13)
Значение 
можно определить осреднением по формуле:
. (14)
В связи с этим произведение величин 
можно представить как произведение средних величин:
. (15)
Заменим дифференциальный оператор правой части в уравнении (6) конечно-разностной схемой для расчета члена с учетом формул (10)-(13), (15):
. (16)
где конечно-разностный аналог дивергенции скорости имеет вид:
. (17)
После расчета члена 
в каждой ячейки сеточной области, необходимо произвести расчет сеточного аналога функции давления 
. Конечно-разностный аналог левой части уравнения (6) для расчета сеточного аналога функции давления имеет вид в результате центральной аппроксимации частных производных второго порядка:
. (18)
После решения уравнения (6) становится известным значения сеточной функции во внутренних узлах конечно-разностной сетки, которые необходимы для нахождения сеточных аналогов компонент скорости.
Аппроксимируем уравнения (1)-(3) конечно-разностными схемами с учетом введенной конечно-разностной сетки и выразим значение сеточного аналога функции компоненты скорости, получаемой на следующем временном шаге. На основании вышесказанного для уравнения (1) имеем:
. (19)
Отсюда выражая 
, получим явную одношаговую схему для расчета сеточного аналога функции 
на новом временном шаге:
. (20)
Аналогично получены формулы для расчета остальных сеточных аналогов функций компонентов скорости течения .
С помощью формул (16), (20) и аналогичных ей можно находить значения простых физических переменных в узлах используемой сеточной области. Полученные данные можно использовать при моделировании циркуляции жидкости в различных водных объектах и системах.
Вывод
Получены формулы, необходимые для моделирования водных объектов и систем, путем численного решения уравнений Навье-Стокса методом маркеров и ячеек. Решение произведено в простых физических переменных: давление и три проекции скоростей на трехмерной системе координат. При использовании соответствующих граничных и начальных условий данное решение позволит с высокой долей адекватности описывать различные водные объекты и системы.
Также следует отметить емкость полученных формул, что, безусловно, отразится на времени проведения расчетов, однако имеющиеся на настоящее время мощности вычислительных машин с многоядерными процессорами (кластеры и суперкомпьютеры) значительно ускорят процесс счета.
Работа выполнена при финансовой поддержке грантов РНФ 14-17-00740 «Озера России - диагноз и прогноз состояния экосистем при климатических и антропогенных воздействиях» и ПСР.1.11.1115А «Совершенствование деятельности в области преподавания географии в ВУЗе и в школе на основе внедрения новых информационных технологий».
Рецензенты:
Филатов Н.Н., д.г.н., профессор кафедры Географии Петрозаводского государственного университета, г. Петрозаводск;
Карпечко Ю.В., д.г.н., ведущий научный сотрудник лаборатории географии и гидрологии ИВПС КарНЦ РАН, г. Петрозаводск.