Исследования вихревых течений в ограниченных водоемах: результаты и перспективы

УДК 519.63, 532.5.031, 551.465

Шавлюгин Александр Иванович

Владивостокский государственный университет экономики и сервиса Владивосток, Россия

Исследования вихревых течений в ограниченных водоемах: результаты и перспективы

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

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

На протяжении нескольких лет на кафедре электроники ВГУЭС и в лаборатории математического моделирования управления научных исследований и опытно-конструкторских работ проводятся исследования вихревых течений в океане. Методическую основу этих работ составляет метод контурной динамики (МКД), который на протяжении около тридцати лет считался одним из наиболее популярных и продуктивных инструментов исследований в теоретической гидродинамике и океанологии. МКД опирается на предположение о кусочно-постоянном распределении завихренности (потенциальной завихренности (ПЗ)) в потоке, что позволяет определять поле скорости посредством контурного интегрирования по границам (фронтам завихренности (ФЗ)) областей постоянной ПЗ (вихревых пятен (ВП)). Это фактически позволяет понизить размерность решаемой задачи и конструировать эффективные численные алгоритмы для изучения эволюции течений.

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

A.И. Шавлюгин Исследования вихревых течений в ограниченных водоемах ...

информацию о первом десятилетии развития МКД в основном для задач классической гидродинамики можно найти в [2].

Океанологические приложения МКД начали развиваться в середине 80-х годов прошлого века в работах группы под руководством

B.Ф. Козлова в Тихоокеанском океанологическом институте ДВНЦ (ДВО) АН СССР (РАН). Начиная с первой работы, в которой МКД был применен к задаче о взаимодействии топографической и свободной завихренности, В. Ф. Козлов с сотрудниками демонстрировали многообразие возможностей для применения метода в задачах геофизической гидродинамики и развивали оригинальные алгоритмические приемы для повышения его эффективности. Обзор исследований, выполненных этой группой с помощью МКД, содержится в работе В. Ф. Козлова [3].

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

МКД-модель круглого баротропного океана

Рассматривается бесконечно глубокий, однородный по плотности океан, поле течения в котором одинаково во всех горизонтальных плоскостях z = const . Предположим, что в начальный момент времени все указанные поверхности оказались разбитыми на множество областей Dt (ВП), в каждой из которых ПЗ П принимает постоянное значение П При сделанном предположении функция тока рассматриваемого течения может быть найдена по формуле:

Y

Di

где (x, у) и (£ n) - декартовы прямоугольные координаты точек наблюдения и интегрирования соответственно,

G - зависящая от граничных условий функция Грина для оператора Лапласа.

Эволюционное уравнение (закон сохранения) для ПЗ:

dp = n t + иП x + vn у = 0

at

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

При помощи формулы Грина можно свести задачу определения поля течения к контурным интегралам по границам Ci ВП, а известные формулы связи функции тока с геострофическими скоростями у ,

V = x¥x

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

~Г = и ~Т< = V

dt , dt .

В случае, когда океан ограничен непроницаемой цилиндрической поверхностью, сечение которой перпендикулярной к оси плоскостью имеет форму круга радиуса Я0, функция Грина принимает вид

*

Яр

(X У;Х,П) =

2р

1п Я - 1п

Яо

где использованы обозначения

Я = [(Х-X)2 + (Л-у)2 ]2

Я

*

(X - X )2 + (Л - у )2 р р = (X2+п2)2

1 1

Р

Детали преобразования интегралов по области в (1) в контурные, основанного на введении вспомогательной функции, можно найти в [4]. Отметим, что впервые баротропная МКД-модель для круглой области была рассмотрена в [5].

Квазигеострофическая МКД-модель круглого двухслойного океана

Детальное описание модели двухслойного океана в квазигеострофи-ческом приближении можно найти в [6]. В этом случае, как показано в [7], давления в слоях (линейно связанные с функцией тока) определяются в виде линейной комбинации интегральных слагаемых вида (1) с ядрами, представляющими функции Грина для операторов Лапласа и Гельмгольца. Задача о построении функции Грина для оператора Гельмгольца в круге была решена в работе [8], в которой было получено выражение

а„ (г,р) = -

Р

¥ К (кЯ ) 1

I 1п

п=орп1п (кЯо) 2 ,

где (г, ф), (р, в) - полярные координаты точек наблюдения и интегрирования соответственно,

к - бароклинный параметр (отношение характерного горизонтального масштаба движения к внутреннему радиусу деформации Россби), то = 2, Иг = 1,г = 1, 2,...,

Кп(х), 1п(х) - функция Макдональда и модифицированная функция Бесселя соответствующего порядка.

Детали преобразования интегралов по области в контурные в формулах геострофических скоростей и некоторые другие подробности численной реализации модели можно найти в работах Шавлюгина А. И. [8].

Численный алгоритм МКД

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

1. Дискретизация ФЗ при помощи множества опорных точек

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

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

3. Вычисление интегралов в выражениях скоростей по формуле трапеций. Поскольку подынтегральные функции содержат логарифмические и степенные особенности, для построения работоспособной численной версии МКД-моделей требуется предварительное аналитическое устранение особенностей, что достигается при помощи интегрирования по частям.

4. Интегрирование уравнений движения (3) по схеме Рунге-Кутта второго порядка точности.

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

6. Перераспределение опорных точек вдоль границ ВП с целью обеспечения их равномерного распределения по длине контуров. Данная процедура обеспечивает повышение точности работы алгоритма в процессе эволюции, которая сопровождается, как правило, локальными растяжениями и сжатиями отдельных участков ФЗ.

Некоторые подробности применяемого численного алгоритма можно найти в работах [4, 8], а также в более ранних публикациях группы под руководством В.Ф. Козлова [9 - 11].

Численные результаты

1. Исследование линейной устойчивости осесимметричных вихревых течений и моделирование нелинейной стадии развития неустойчивости.

Изучение гидродинамической устойчивости различных течений является одной из основных областей исследований в классической и геофизической гидродинамике. Вследствие этого существует множество методов и подходов как для аналитического, так и для численного решения указанных задач. Описываемые в настоящей статье результаты были получены на основе подхода, который для моделей баротропного и двухслойного неограниченного по горизонтали океана подробно изложен в работах [7, 12].

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

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

Баротропный компенсированный вихрь Рэнкина

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

П г2 + П (г2 - г2) = о

пенсации 1 1 2 2 1 течение в невозмущенном начальном со-

стоянии отсутствует всюду за пределами кольцевой оболочки. Подобные модельные структуры гораздо лучше воспроизводят наблюдаемые в реальных океанских вихрях радиальные профили скорости, чем просто вихрь Рэнкина, поле скорости которого достаточно медленно спадает по мере удаления от его центра (обратно пропорционально расстоянию). Исследования нелинейной стадии развития неустойчивости компенсированного вихря Рэнкина в модели неограниченного по горизонтали океана были впервые выполнены в [12], где было установлено, что исходная вихревая структура трансформируется в центральное ядро часто квазим-ногоугольной формы, окруженное сформировавшимися из компенсирующей оболочки вихрями-сателлитами, количество которых зависит от доминирующей в процессе моды возмущения. Исключение составляет случай второй моды, когда в зависимости от геометрических параметров задачи возможны два сценария развития процесса: первый из них укладывается в описанную выше схему и приводит к формированию структуры, которая называется триполем; во втором случае центральное ядро постепенно растягивается и значительно позже оболочки делится на две части, в результате чего образуются два диполя, которые расходятся в противоположные от центра структуры стороны благодаря присущему диполям механизму самодвижения.

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

На рисунке 1 в указанные безразмерные моменты времени представлена эволюция неустойчивого компенсированного вихря Рэнкина

(^ = °Д5, Г2 = °^5, П1 = ^°, П 2 =-о,92о(45)), подтверждающая сформулированное выше наблюдение.

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

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

157,5 163,5

Рис. 1. Эволюция неустойчивого компенсированного вихря Рэнкина

Двухслойный осесимметричный вихрь

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

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

Исследование неустойчивости двухслойных осесимметричных вихрей в круглом океане было выполнено в работе [8], в которой автор описывает некоторые эффекты влияния береговых границ. В частности, было обнаружено, что в процессе распада исходных партнеров на несколько двухслойных вихревых пар наблюдается отражение этих аналогичных баротропным диполям структур от границ бассейна и их возвращение в центр области. Еще один пример эволюции неустойчивых хетонов для случая доминирования третьей моды возмущения приведен на рис. 2. Этот пример отличается от описанного в [8] эксперимента тем, что сформировавшиеся мелкомасштабные двухслойные вихри после достижения береговых границ не отражаются сразу же от них, а совершают вдольбе-реговое движение, обмениваясь партнерами по вихревой паре, и только затем возвращаются в центр бассейна.

35 70 85 95

103 111 117,5 123

Рис. 2. Эволюция неустойчивого двухслойного вихря

2. Исследование эволюции эллиптического ВП. Некоторые модельные объекты, не обладающие такой степенью симметрии, как рассмотренные выше, привлекают повышенное внимание исследователей в силу

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

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

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

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

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

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

3. Исследование нелинейного взаимодействия ансамблей ВП. Задачу о взаимодействии тождественных ВП можно отнести к классическим и достаточно подробно изученным задачам гидродинамики. Исследования, проводившиеся в данном направлении начиная с 70-х годов прошлого века с помощью разнообразных лабораторных методик и численных моделей, обнаружили тенденцию вихревых образований к увеличению размеров посредством парного взаимодействия (pairing или merging). Аналитические и численные исследования, выполненные преимущественно с помощью МКД, позволили установить различные критерии инициации процесса слияния ВП. В частности, в [17, 18] на основе найденных численно стационарных форм взаимодействующих ВП и исследования их устойчивости был сделан вывод, что решения, для которых отношение радиуса эквивалентного круглого ВП к расстоянию между центрами тяжести партнеров превышает 0,32, являются неустойчивыми и их эволюция приводит последовательно к деформации границ ВП, филаментации (образованию тонких нитей завихренности) и последующему слиянию.

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

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

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

Рис. 3. Трехвихревое взаимодействие томсоновского ансамбля ВП, локализованных в верхнем слое двухслойного океана.

4. Построение стационарных состояний ВП. Под стационарным состоянием (V-state в английской терминологии) ВП в гидродинамике понимается такая его форма, которая остается неизменной в некоторой системе координат. В частности, возможны ротационные стационарные состояния, когда вихревая структура вращается твердотельно с постоянной угловой скоростью, например эллипс Кирхгофа, и трансляционные стационарные состояния, когда объект перемещается прямолинейно и равномерно, например вихревой диполь. Исследованию стационарных состояний за последние 30 лет было посвящено множество публикаций, начиная с пионерской работы Г.С. Дима [19]. Для отыскания решений было разработано множество алгоритмов, которые с успехом применялись в различных модельных задачах.

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

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

vX - uy - W

где точками обозначены производные по параметру v. Форма ВП ищется в полярных координатах, причем зависимость полярного угла от пара-

в=p

метра предполагается известной (линейной) N (N - число точек на

контуре), а полярный радиус определяется в виде усеченного ряда Фурье

^, 2mv

Р = Z bn cos--

с неизвестными коэффициентами n N . Условие (4) записыва-

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

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

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