CC BY
24
2019, Т. 161, кн. 1
С. 5-23
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
УДК 53.072+519.95 doi: 10.26907/2541-7746.2019.1.5-23
МОДИФИЦИРОВАННЫЙ МЕТОД КОНТУРНОЙ ДИНАМИКИ И МОДЕЛИРОВАНИЕ ВИХРЕВЫХ СТРУКТУР
В.Ю. Белашов1, О.А. Харшиладзе2
1 Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия 2 Тбилисский государственный университет им. И. Джавахишвили, г. Тбилиси, 0179, Грузия
Аннотация
В работе рассматривается один из наиболее эффективных методов моделирования вихревых структур, описываемых двумерным уравнением переноса завихренности и уравнением Пуассона для функции тока, - метод контурной динамики, базирующийся на представлении вихревого потока в виде вихревых областей конечной площади. Предлагается модификация метода контурной динамики, минимизирующая погрешности, возникающие при прямом его применении к исследованию вихревых структур. Приводятся примеры результатов численных экспериментов по исследованию динамики взаимодействия вихревых структур для различных конфигураций их взаимного расположения, знаков завихренности и расстояний между границами вихревых областей конечной площади.
Ключевые слова: вихри, вихревые области конечной площади, моделирование, гидродинамика, модифицированный метод контурной динамики, режимы взаимодействия, явление квазивозврата, фазовое перемешивание
Введение
При решении любой задачи неизбежно возникает вопрос о выборе оптимального метода ее решения. При этом выбор метода определяется рядом факторов, среди которых необходимо отметить его временные затраты, простоту (при соблюдении условия адекватности решаемой физической задаче), надежность и универсальность. Для задач, изучаемых в рамках теории вихревых движений, выбор методов численного решения зачастую ограничен набором известных схем, базирующихся на конечно-разностной аппроксимации дифференциальных уравнений в частных производных, описывающих движение жидкости, газа и плазмы. Однако такой подход не всегда приемлем в силу ограничений, которые накладываются на разностные сетки. Моделирование вихревых явлений с помощью некоторых специальных методов зачастую связано с использованием мощных компьютеров и суперкомпьютеров, поскольку вычисления в этом случае требуют хранения данных на 1010 - 1011 слоях памяти и быстродействия в несколько террагерц [1]. Поэтому, несмотря на широкое распространение вычислительной техники и сравнительную ее доступность, поиск наиболее оптимального метода моделирования для каждой конкретной задачи остается актуальным и в настоящее время. Если говорить об аналитических решениях уравнений гидрогазодинамики, то еще Ламб [2] отмечал, что «движение тела в завихренном потоке жидкости - интересная, но, к сожалению, трудноразрешимая задача», поскольку уравнения Навье-Стокса имеют явное решение лишь в нескольких частных случаях нетурбулентного движения жидкости, которые практического интереса не представляют. Тем не менее разработка
теории турбулентных явлений всегда оставалась одной из актуальных проблем. Примером может служить работа [3], в которой явление турбулентности изучается с позиций статистической механики, или работа [4], где рассматриваются теоретические основы проблемы изучения вихревых движений с использованием аналитических методов.
Существенный прогресс в области исследования вихревых явлений наметился в конце 60-х годов прошлого века, когда появились достаточно мощные компьютеры и родилась новая наука, называемая рядом исследователей вычислительной динамикой жидкости [5, 6]. С тех пор исследование вихрей значительно продвинулось вперед, но и по сей день многие важные проблемы остаются нерешенными, поскольку охватить все стороны вихревых движений вряд ли возможно в силу неисчерпаемости их многочисленных комбинаций. Поэтому, возвращаясь к вопросу о выборе универсальной, адекватной и простой (в смысле доступности использования с применением имеющихся вычислительных средств) методики исследования, необходимо еще раз отметить его важность. В настоящей статье рассматривается один из эффективных методов, построенный на представлении вихревых структур в виде вихревых областей конечной площади (ВОКП) [7], в связи с чем описание их динамики существенно упрощается. Для численного исследования предлагается использование специального алгоритма, называемого алгоритмом контурной динамики (КД) [8], и рассматривается его модификация, позволяющая моделировать не только отдельные вихревые структуры, но и исследовать динамику N-вихревых систем, состоящих из отдельных вихрей, в зависимости от их относительного положения, порядка симметрии, величины и знака завихренности. В работе рассматриваются также некоторые результаты моделирования вихревых структур различной конфигурации с помощью модифицированного алгоритма КД.
1. Основные уравнения и методы моделирования вихревых явлений
В простейшем случае несжимаемой вязкой среды уравнениями, описывающими движение жидкости, являются уравнения Навье-Стокса
+ (уУ)у = — Ур + Уу = 0,
от р
где V - кинематическая вязкость, р - плотность среды, р - давление. Поскольку рассматривается вихревое движение, от этих уравнений, исключая давление и вычисляя ротор от обеих частей, мы переходим к уравнению переноса для завихренности и уравнению Пуассона для функции тока. В результате получаем нестационарное уравнение для завихренности среды
+ (уУ)с = VV2C,
где £ - завихренность, определяемая выражением £ = [V, у] , и уравнение Пуассона
Аф = -С. (1)
Скорость в данном случае определяется как V = [V, ф], где ф - функция тока. Таким образом, получается система уравнений, описывающая движение несжимаемой среды. В двумерном случае (который представляет особый интерес, например, в случае жидкости) эта система имеет достаточно простую форму
ду С + (уУ)С = V У2С, Аф = -С, V = [V, фег ], где ех - единичный вектор внешней нормали.
Моделирование процессов и явлений, описываемых этими уравнениями, осуществляется с применением специальных методов, к которым можно отнести: конечно-разностные схемы, метод «частиц в ячейке» (Р1С-модель), метод «водяного мешка», дискретных вихрей, конечных элементов и ряд других [9, 10]. Отметим некоторые особенности этих методов для того, чтобы выявить их положительные и отрицательные стороны, а также целесообразность использования при решении конкретных задач.
В отношении конечно-разностных схем, являющихся в настоящее время наиболее распространенными, нужно отметить, что этот подход, хотя и универсален, но имеет ряд недостатков, связанных с появлением искусственной дисперсии и диссипации на сетке. Например, в предложенной в [9] схеме Лакса-Вендроффа для уменьшения численной диффузии мы вынуждены находить промежуточные значения в дробные моменты времени. Не говоря уже о том, что для расчета необходимо накладывать на исследуемую область сетку с очень малым шагом, что, в свою очередь, приводит к необходимости введения больших ограничений на шаг по времени, который должен удовлетворять условию устойчивости Куранта -Фридрихса-Леви. В результате значительно увеличивается время расчетов и возрастают требования к используемой вычислительной технике.
Другой метод моделирования - так называемый метод «частиц в ячейке» (Р1С-модель) - был первоначально использован для моделирования плазменных явлений, но нашел также применение и в гидродинамике. При использовании Р1С-модели двумерная несжимаемая среда представляется как система «вихревых частиц» [9]. При этом элементы завихренности заменяются частицами, и плотность таких частиц в пространстве определяет распределение завихренности. Таким образом, используя бесстолкновительную модель «частиц в ячейке», можно моделировать двумерную завихренную несжимаемую среду. В частности, с помощью РТС-модели может эффективно исследоваться неустойчивость Кельвина-Гельмгольца, при развитии которой образуются вихревые структуры.
Еще один из распространенных и активно используемых методов исследования вихревых явлений - метод дискретных вихрей, подробно описанный в [11], основан на представлении завихренной среды в виде отдельных точечных вихрей, между которыми происходит взаимодействие. Следует, однако, заметить, что этот метод имеет один существенный недостаток, заключающийся в необходимости просчитывать все вихри не только по их границе, но и внутри ВОКП, тогда как это делать не обязательно, поскольку завихренность внутри вихревой области конечной площади одинакова во всех точках.
Рассматриваемый в работе метод КД является обобщением метода «водяного мешка», традиционно использующегося при изучении динамики плазмы, описываемой уравнением Власова [12]. Основная идея метода КД заключается в моделировании исследуемой среды с помощью функции распределения некоторой физической величины, задающей ее конфигурацию. Например, для идеальной жидкости с локальным вихревым возмущением это распределение завихренности, представляющее собой контур, внутри которого жидкость вращается с постоянной угловой скоростью. В этом случае для моделирования эволюции жидкости исчезает необходимость рассматривать всю область, занимаемую средой (что приходится делать при использовании конечно-разностного метода) или область внутри вихря (как в методе дискретных вихрей), достаточно численно решать уравнения для границы вихревой области (или областей), что значительно экономит машинное время и память. К этому надо добавить отсутствие искусственных дисперсии и диссипации, характерных для разностных схем.
2. Алгоритм метода контурной динамики
Основные положения метода описаны в работе [8], однако при его прямом использовании возникает ряд трудностей, поэтому для корректного применения к моделированию вихревых структур метод КД требует некоторой модификации. С целью более полного выяснения сути предлагаемого подхода, а также надлежащего его обобщения приведем ниже основные положения алгоритма метода КД, ограничившись двумерным случаем. Исходными уравнениями в методе контурной динамики являются уравнение переноса вихря (для идеальной среды)
дС дС дС
Ж + идХ + % = 0 (2)
и уравнение Пуассона для функции тока (1).
Систему уравнений (1), (2) необходимо также дополнить выражениями для компонент скорости
и = дф/ду, V = —дф/дх.
Как уже было сказано, идея метода КД состоит в том, что рассматривается взаимодействие между границами областей, причем завихренность каждой области является постоянной (но необязательно одинаковой для всех вихревых областей), благодаря этому размерность задачи снижается на единицу. Области положительной завихренности (или циркуляции) соответствуют вращению конвективных жидких элементов против часовой стрелки, поэтому мы будем изучать правостороннюю координатную систему, когда вектор £ направлен вдоль оси е2 = ех х еу. Аналитическое решение уравнения Пуассона (1) для функции тока ф имеет вид
ф(х,у) = — ^П [1пг]С(С,п), (3)
где
г = [(х — С)2 + (у — П)2]1/2, (4)
1п г - функция Грина уравнения (1).
Значение скорости в любой точке потока можно получить дифференцированием интеграла (3): и = Ухе2ф = ехдуф— еудхф. Проводя замену переменных ду ^ —дп и дх ^ —д^, исходя из (4), получим
и = 2тт II ¿С^пС[ехдп[1пг] — еуд?[1пг]]. (5)
Интегрируя далее (5) по частям, найдем, что
и = — 11 ^¿п [1пг][ехдпС — еуд5С], (6)
1
2П
ех х
Ц ^П [1п г] У С.
Индекс С у оператора Гамильтона, означает следующее: у = ехд^ + еудп и Ух = ехдх + еуду . Интеграл в (6) содержит две части, описывающие: 1) область внутри контура, где уС = 0, которая не влияет на скорость, и 2) тонкую полосу на границе контура, где уС = 0, которая дает единственный конечный вклад в изменение скоростей вихрей. Далее мы вводим на контуре локализованную ортогональную координатную систему с осями в и q, как это показано на рис. 1. Здесь площадь ^¿в = ¿С ¿п. Это преобразование с функциональным определителем модуля
и
(якобианом) справедливо, если контур является кусочно-постоянным (поскольку тангенциальная и нормальная составляющие могут меняться дискретно). Тогда
п = ^С| и С = ^,
Z = [C]^(q)cos(n,e) = [Z]¿(q)sin 0,
(7)
dn Z = [C]^(q)cos(n,n) = — [Z]¿(q)cos 0,
где [Z] = Zo — Z/ (см. рис. 1, а) и ó(q) - дельта-функция, для которой аргумент q положителен в направлении en. Компоненты направляющих косинусов есть cos(n, £) = sin 0 и cos(n, n) = — cos 0. Если мы заменим компоненты VС, используя соотношения (7), и выполним интегрирование по q, то получим
1 N
" 7 f [m r] [ex cos + ey sin 07 ] usj,
7=1
(8)
u = — V[Z]j </ [Inr] [ex cos 0j + ey sin 0j] dsj, 2n 7
1 N /
u = 2П ^[Z7 / [ln r] [ex + ey dnj]'
7
7 = 1
Здесь [Z] j - значение [С] на контуре j. Таким образом, области постоянной завихренности заменяются распределением источников с логарифмическими силами вдоль контуров (обозначенных индексом j), окруженных конечным числом Nc вихревых областей в поле завихренности.
2.1. Пространственная дискретизация. Будем предполагать, что контур, ограничивающий область постоянной завихренности, представляет собой многоугольник с расположенными на нем N узлами и перейдем к исследованию их взаимодействия. Как показано на рис. 1, б, вдоль каждого сегмента длиной hn между узлами n и n +1, f аппроксимируется следующим образом:
r2 = (x — е„ — е')2 + (y — nn — е' tg 0„)2
или f2 = гП(1 — 2Z^n + Z2), где rn - расстояние от n-го узла до точки (x,y) (то есть другого узла), фп - соответствующий угол, и
Mn =cos(^n — 0n), Z = (е — Cn)/(rn cos 0n) = е'/(rn cos 0n), (9)
а tg 0n = (yn+i — yn)/(xn+i — xn). Для рассматриваемых нами вихревых областей конечной площади уравнения (8) могут быть дискретизированы следующим образом:
N
u cos 0n + ey sin 0n], (10)
n= 1
где
1
(Ди)„ = 4Л[СК | [1п(?)2]^ (11)
о
представляет собой вклад скорости от линейного сегмента, находящегося между узлами п и п + 1. После несложных алгебраических преобразований получим
(Ди)п = [1п(г„)2] + (1 - Мп^-1г„)[1п(д„)] - 2+
+ 2(1 - МП)1/2(г„М„)аг^(в„)}, п =1,Ж; (12)
а)
Рис. 1. Иллюстрация к алгоритму метода КД. Основные обозначения
где
qn = 1 — (Ьп/г„) + (Ьп/г„)2 (13)
вп = (Ьп/гп) (1 — ^)1/2 [1 — М^п/гп)]-1. (14)
Заметим, что Ие {Сп[1п(1 — ЬпСп/гп)]} = ^п[1п^/2)] — (1 — ^)1/2 агС^(вп), где (п и
Сп^ - корни уравнения 1 — + £2 =0. Если значение п = 1 сопоставить с точкой с координатами (х,у), то на границе будем иметь
(Ди)* = ^[С]^ {[1п(Ь*)] — 1}, п = N; (Ди)1 = -П [С ]Ь1{[1п(Ь1)] — 1}, п = 1.
Вклад от узлов, которые удовлетворяют условию Ьп/гп ^ 1, может быть получен из асимптотического разложения (12) или непосредственно из (11):
(Ди)п = -^[С]Ьп{[1п(гп)] — 1/2^п(Ьп/гп) + 1/6(1 — ^п)(Ьп/гп)2 + 4п
+ 1/4^(1 — 4/3^п)(Ьп/гп)3 + 0(Ьп/гп )4 }.
Подобная процедура может быть использована для нахождения взаимодействия не только между узлами одного контура, но и между узлами разных, пространственно разнесенных, контуров.
2.2. Временная дискретизацияю. Рассмотрим теперь точку (хт,ут) на границе ВОКП. Каждый узел этой вихревой области вращается с локальной скоростью жидкости
N
и(хШ] ут) = хт = Хтеж + утеу = ^(Ди)жп|е1 СОв + еу вШ (15)
п=1
где (Дм)т„ вычисляется по формуле (10) с заменой г„ на г„т = [(хт — ж„)2 + + (ут — уп)2]1/2 , и ф„ = (ут — уп)/(жт — !„). Поскольку изучается консервативная система, то можно использовать центрированную по времени схему с перешагиванием с погрешностью аппроксимации 0(т2)
N
т+1 — хт-1 =2т]Т (Д<( СО8 ^ вх + В1п ^ еу). (16)
х _
т ^т /
п= 1
Устойчивость схемы (16) гарантируется условием Куранта - Фридрихса - Леви т < шах{м,.
Таким образом, задавая начальную конфигурацию вихревых областей конечной площади, мы можем исследовать временную эволюцию системы ВОКП.
2.3. Схема расчета. Алгоритм моделирования эволюции ВОКП состоит из нескольких этапов.
1. Для каждого контура задается количество узлов N (которое должно выбираться таким образом, чтобы контур можно было приближенно считать непрерывной кривой), первоначальную форму, и завихренность.
2. В соответствии с выражением (15) вычисляется вклад всех контуров в изменение скорости каждого узла на текущем контуре, при этом для нахождения используется одна из схем, предложенных в работах [8, 7, 13]:
• в простейшем случае [8] (Дм)тп вычисляется как
(Ди)тп = 71[С]^т{[1п(^т)] — 1}, если т = п — 1 или т = п;
2п
(Ди)тп = ]^т{[1п(гт)] — 1}, если т = п — 1 или т = п;
2п
• (Дм)тп может быть также получено по формулам [7]:
(Ди)тп = —[С]^т{[1п(^т)] — 1}, если т = п — 1 или т = п;
2п
(Дм)тп = ^[С [1п(Гт)2 + (1 — Мт^т/Л-т) М^т) — 2+
+ 2(1 — мт)1/2(гтМт)аг^ вт], если т = п — 1 или т = п где мт, , вт вычисляются по формулам (9), (13) и (14);
• могут быть также использованы формулы, полученные в работе [13]:
£
(Дм)тп = 7Г~ М(1 + Ат„)1п(Гт„+1) — Атп 1п(Гтп-1) — 1 +
2п
+ Втп arctg 2Вт" + п|Втп|Н(—Стп)],
стп
где
л _ (хп Хт)(хп+1 хп) + (уп ут)(уп+1 уп)
Атп — 72 ,
7п
Т-, (хп Хт)(уп+ 1 - Уп) - (Уп - Ут)(Хп+1 - Ж„)
72 7п
г2 I г2
_ тп + т+1п -|
Стп — 7п 11
Н - функция Хевисайда
Н(г) — I1' * >
^ г <
Если п — т и Атп — Втп — °, то (2.3.) преобразуется в (Ди)тт —
— ^~[С]|[1п(Гтт+1)] - 1} , а в случае п — т - 1 , то есть Атт-1 — -1 , Втт-1 — °, 1
имеем (Ди)тт-1 — 7Т- [С]{[1п(Гтт-1)] - 1} . 2п
3. Процесс продолжается либо до выхода на стационарный режим, либо до достижения заданного в эксперименте значения времени.
3. Модификация метода КД
Для того чтобы описанный выше метод КД давал устойчивые решения, его необходимо модифицировать, поскольку, во-первых, схема с перешагиванием дает ощутимую погрешность на достаточно больших временных промежутках, а, во-вторых, как показывают эксперименты, эволюция ВОКП приводит к «разбеганию» узлов контура, что недопустимо, так как при этом нарушается главное условие -контур должен быть кусочно-непрерывным, то есть как можно более точно аппроксимировать непрерывную линию.
Первая проблема может быть решена с помощью методики, предложенной в работе [12]. Для устранения погрешности в методе с перешагиванием следует синхронизировать значения переменных на двух временных слоях. При таком подходе предполагается, что на некотором временном слое известны координаты и скорости всех узлов контура. Тогда асинхронная компонента может быть отфильтрована с помощью усредненных переменных для момента времени р +1/2, причем р должно быть четным. Эти переменные находятся тривиально:
хр+1/2 — 1 {хр+1 + хр }, ир+1/2 — 2 {ир+1 + ир}.
Далее, задавая приращение за время т/2:
хР+1 — хр+1/2 + тир+1/2 хр — хр+1/2 _ тир+1/2 2 ' 2 '
вычисляем новое положение системы контуров. Синхронизация производится через заданный промежуток времени.
Вторая трудность может быть преодолена введением дополнительных узлов в местах «разрыва» контура, когда 7п > 27тах для сохранения точности вычислений в области «разрыва» вводится новая точка С с координатами
хс — (ха + хв )/2, ус — (уа + Ув )/2,
Рис. 2. Три случая образования «двойных» кривых
расположенная между соседними точками А и В. При этом, если используется схема с перешагиванием, дополнительные узлы необходимо добавить и на двух предыдущих временных слоях. Следует также отметить, что возможна и обратная ситуация - при сильном сближении двух узлов один из них можно исключить из расчетов, в этом случае точность вычислений не изменится, но будет экономиться машинное время.
Еще одним способом оптимизации расчетов является так называемое «подправ-ление» контуров, описанное в [12]. С течением времени линии одного и того же контура или линии, принадлежащие разным контурам, могут сблизиться друг с другом настолько, что площадь между ними будет стремиться к нулю. Тогда эти линии становятся одной кривой, которую приходится проходить дважды. Если значение функции по обе стороны от «двойной» кривой одинаково, то, исключив этот участок, мы не изменяем состояния системы.
На рис. 2 показаны возможные варианты образования «двойных» кривых. Рассмотрим случай, представленный на рис. 2, а. Если кривая образует петлю, то ее изменяют следующим образом: точки С и О соединяют между собой, при этом образуется новая замкнутая кривая, начинающаяся в точке А и оканчивающаяся в точке В. Далее, в зависимости от рассматриваемой задачи, точки между С - А и В - О либо заменяют одной кривой, либо полностью исключают из расчетов. Ситуация, изображенная на рис. 2, б, отличается от предыдущей тем, что необходимо рассматривать два участка, при этом соединяются соответственно точки С и А и точки В и О. На рис. 2, в показан случай, когда двойная кривая разделяет две области с разными значениями функции /. Здесь можно соединить точки С и А, а также В и О и определить кривые, которые соединяют точки А и В двух других контуров, для которых А/ = /3 — /1.
С учетом сделанных замечаний схема расчета для модифицированного метода КД будет содержать два дополнительных этапа и примет следующий вид.
1. Для каждого контура задается количество узлов Nj, его первоначальная форма, и завихренность.
2. В соответствии с выражением (15) вычисляется вклад всех контуров в изменение скорости каждого узла на текущем контуре.
3. Выполняется корректировка на узлах - вводятся дополнительные и/или исключаются лишние.
4. Проводится синхронизация значений скорости на трех временных слоях для устранения погрешности схемы «с перешагиванием».
5. Повторяются пп. 2-4 либо до выхода на стационарный режим, либо до достижения заданного в эксперименте значения времени.
4. Диагностика метода
Для того чтобы проверить адекватность метода КД, в качестве тестового был выбран простейший случай - эллиптический вихрь Кирхгофа. Для него имеется аналитически решение [2], согласно которому период обращения вихря является функцией отношения большой а и малой Ь полуосей эллипса:
г(ь/„)—. (1-)
В соответствии с формулой (17) для эллипса с отношением полуосей 1/2 период колебаний границы будет равен Т — 28.4. Тестировалось еще одно свойство системы: движение частиц жидкости в вихре должно происходить в том же направлении, что и направление вращения, но с большей скоростью. Взяв в качестве начального возмущения эллиптический вихрь с а/Ь — 1/2, мы промоделировали с помощью модифицированного метода КД его эволюцию. Результаты численного эксперимента показали, что период обращения, наблюдаемый в экспериментах, совпадает с аналитическим решением (17), при этом отчетливо наблюдается движение частиц на границе: линия, соединяющая первый и средний узлы контура, совершает вращательные движения с вдвое большей скоростью.
В работах [13, 14] проводилась диагностика и проверка устойчивости метода КД в его «классической» формулировке, при этом было показано, что на границе контура появляются слабые колебания узлов, рост амплитуды которых незначителен на рассматриваемых временных масштабах. Тестирование описанной выше модификации алгоритма проводилось нами путем сравнения результатов моделирования для вихрей с порядками симметрии т — 1, 2, 3 с известными результатами [7], полученными стандартным методом КД. При этом было отмечено, что указанная выше неустойчивость отсутствует. В связи с этим можно утверждать, что модифицированный метод КД применим для исследования эволюции вихревых структур на значительных временных интервалах.
5. Результаты моделирования вихревых структур различных конфигураций
Для иллюстрации эффективности рассмотренного выше метода нами было проведено численное моделирование эволюции вихревых структур типа ВОКП с помощью модифицированного алгоритма КД. В численных экспериментах исследовалась зависимость характера взаимодействия вихревых структур от нескольких параметров: от расстояния между границами ВОКП, от числа и конфигурации их взаимного расположения, направления вращения и порядка симметрии вихрей. Сначала мы рассмотрели простейшую конфигурацию - два вихря с одинаковыми знаками завихренности и различными расстояниями г между границами. Численные эксперименты показали, что при расстоянии, большем критического -3й/4 (где й - диаметр вихря), вихри вращаются вокруг общего центра (рис. 3а), причем скорость этого вращения зависит от расстояния г. При этом происходит заметное изменение формы вихрей - от круглой до формы вытянутого эллипса, но с течением времени происходит возврат к первоначальной форме (так называемое
Рис. 3. Взаимодействие двух вихрей с (1 = £2 = —1 при начальной дистанции между ними: а) 5 = 2й; б) 5 = й/2
Рис. 4. Взаимодействие круговых ВОКП с различными знаками завихренности (£ = — 1 и С =1)
явление «квазивозврата»). При уменьшении расстояния между вихрями происходит качественное изменение в их поведении - они начинают взаимодействовать, и с течением времени образуется общая вихревая область, которая сама состоит из завихренностей более мелкого масштаба, что можно видеть на рис. 3, б).
Далее было исследовано взаимное влияние двух вихрей с разными знаками завихренности, но одинаковой интенсивности. Как видно из рис. 4, их взаимодействие принципиально отличается от описанного выше случая. Вихри сближаются и движутся в направлении, перпендикулярном оси, соединяющей их центры. Направление и скорость движения вихрей зависят от знаков завихренностей и взаимного расположения, причем чем больше завихренность, тем скорость движения выше.
В следующей серии экспериментов мы изучили более сложную конфигурацию, впервые рассмотренную в задачах моделирования плазменных вихревых систем [15] - четыре вихря, расположенных вдоль одной прямой. При одинаковом расстоянии между вихрями, равном й/2, наблюдался эффект слияния двух внутренних вихрей, которые затем разрывались двумя крайними. В дальнейшем, в центре образовывался небольшой вихрь, вероятно (в пределах времени счета), устойчивый (рис. 5, а). При увеличении дистанции между вихрями до расстояния, равного й, характер их взаимодействия менялся (рис. 5, б). Система разбивалась на пары, внутри которых происходило взаимодействие, описанное выше для случая г < < 3й/4. Образовавшиеся при этом две вихревые структуры вращались как целое вокруг общего центра. При увеличении расстояния до образовывалась устойчиво вращающаяся вокруг общего центра структура из разбитых на пары, но не взаимодействующих друг с другом вихрей.
На рис. 6 представлены результаты для случая, когда два внутренних вихря находятся на расстоянии й/2, а два крайних вихря отстоят от них на г — й. Можно
/ = о / = 10 % б)
^ %
/ =20 • / = 32 -*
• /
Рис. 5. Взаимодействие четырех линейно расположенных при £ = 0 вихрей с = —1 при: а) 8г = 6/2; б) й = 2 6
Рис. 6. Взаимодействие четырех ВОКП с одинаковыми знаками завихренности (£ = — 1) при г = 6/2 и г = 6
видеть, что в процессе эволюции такой системы имеет место взаимодействие внутренних вихрей, тогда как два крайних во взаимодействии не участвуют, вращаясь вокруг центральной вихревой структуры. Подобное поведение, по-видимому, объясняется тем, что при появлении по краям двух вихрей в результате смещения «центра масс» системы возникает дополнительное «притяжение», которое заставляет внутренние вихри сблизиться и вступить во взаимодействие, которое в случае двух вихрей, находящихся на таком же расстоянии друг от друга, не возникает.
Еще одна конфигурация, изученная нами, - это три вихря, расположенные в вершинах равностороннего треугольника [16]. При одинаковых знаках завихрен-ностей и расстоянии между границами ¿/2 три вихря, сливаясь, образуют один большой вихрь, изображенный на рис. 7. При увеличении расстояния до г = 6 они устойчиво вращаются вокруг общего центра. Для конфигурации из четырех вихрей с одинаковыми знаками завихренности, расположенных в вершинах квадрата, характер эволюции качественно совпадает.
Дополнительно были проведены эксперименты для вихрей с различными знаками завихренностей. На рис. 8 показано движение вихрей в зависимости от их взаимного расположения: для случая, изображенного на рис. 8, а, как и можно было ожидать, эволюция подобна поведению пары вихрей с противоположными знаками - они попарно разлетаются в противоположные стороны. Однако результаты численного эксперимента говорят о том, что взаимодействие вихрей с противоположными знаками завихренностей происходит интенсивнее, чем для «одинаковых» вихрей. Особенно это заметно на рис. 8, б, где пара противоположно вращающихся ВОКП «уходит» в правую сторону, а два других вихря пытаются ее «догнать».
Рис. 7. Взаимодействие трех вихрей с одинаковыми знаками завихренности (£ = —1) при г = (/2
Рис. 8. Взаимодействие четырех вихрей при 5г = (: а) = Сз > 0, (2 = С4 < 0; б) С1 = С2 > 0, Сз = С4 < 0
Следующие серии численных экспериментов были посвящены изучению эволюции вихревых областей, имеющих различные порядки симметрии - эллиптических и «треугольных» ВОКП. Примеры результатов представлены на рис. 9. Видно, что на ВОКП треугольной формы с течением времени появляется точка заострения, которая затем вытягивается в нить, что приводит в результате к разрушению вихря. С эллиптическим вихрем этого не происходит, несмотря на то что форма его тоже меняется. Это говорит о большей устойчивости вихрей с меньшим порядком симметрии.
Суммируя основные результаты, можно сделать следующие выводы: во-первых, основную роль во взаимодействии вихрей играет расстояние между ними - чем меньше расстояние, тем сильнее их взаимное влияние, переходящее во взаимодействие при достижении критического расстояния; во-вторых, поведение вихрей с разными знаками завихренностей принципиально отличается от поведения «одинаковых» вихрей, в том числе характеризуется большей интенсивностью; в-третьих, взаимное расположение вихрей влияет на поведение структур, причем влияние это приводит к качественным отличиям в характере их взаимодействия; и, в-четвертых, от формы вихря зависит его устойчивость к разрушению, которая особенно проявляется при взаимодействии вихрей с разными порядками симметрии.
Заключение
Мы показали, что при изучении вихревых явлений важным является выбор оптимального метода исследования для каждой конкретной задачи. Среди специальных методов моделирования эволюции и динамики взаимодействия ВОКП наибо-
Рис. 9. Взаимодействие вихрей с разными порядками симметрии (т = 3 и т = 2), но одинаковыми знаками завихренности (£ = —1)
лее предпочтительным с точки зрения минимизации временных затрат является метод контурной динамики, который, кроме прочего, обладает рядом дополнительных преимуществ, например отсутствием искусственных диффузии и диссипации, которые появляются при численном решении с помощью конечно-разностных схем, кроме того, представление среды в виде отдельных вихревых областей конечной площади, между которыми происходит взаимодействие, исключает необходимость рассматривать все пространство, в котором находятся вихри, а также вычислять значения скоростей и завихренностей внутри этих областей (как в методе дискретных вихрей), что существенно ускоряет процесс моделирования. Однако при прямом использовании алгоритма метода КД появляются погрешности, которые связаны с «разрывом» контуров, а также с погрешностью метода «с перешагиванием», с помощью которого вычисляется временная эволюция ВОКП. Поэтому метод КД требует некоторой модификации, которая была нами осуществлена. Результаты численных экспериментов, полученных с помощью модифицированного алгоритма КД, позволяют заключить, что этот метод является весьма эффективным при исследовании вихревых явлений, наблюдаемых в средах, где присутствуют локальные вихревые области, между которыми происходит взаимодействие. Выполненная нами модификация стандартного алгоритма метода КД позволяет также, наряду с изучением эволюции отдельных вихревых образований, эффективно изучать эволюцию и динамику взаимодействия N-вихревых систем различных пространственных конфигураций, состоящих из областей конечной площади, в зависимости от пространственного расположения, порядка симметрии, величины и знака завихренности отдельных вихревых образований.
Отметим, что в класс вихревых структур, которые можно исследовать методом КД, кроме ВОКП и точечных вихрей, можно включить также вихревые пелены [18-20]. Например, «нити» завихренности, которые наблюдаются после точек заострения (рис. 9), можно отделить от ВОКП и рассматривать как вихревые пелены конечной протяженности. Кроме того, вихревые пелены естественным образом возникают при изучении устойчиво стратифицированной (слоистой) среды со сдвигом скорости, описываемой уравнениями Буссинеска. Например, в [7] обнаружено, что если начать с ВОКП и областей с кусочно-постоянной плотностью, то затем образуются вихревые пелены различной интенсивности, которые движутся вместе с линиями разрывов плотности. Их интенсивность изменяется в соответствии с изменением угла между касательной к линии разрыва и нормалью к ним.
Получаемые в численных экспериментах с помощью модифицированного метода КД результаты, наряду с их очевидной значимостью для адекватной интерпретации эффектов, связанных с турбулентными процессами в газах и жидкостях (в частности, вихревых движений в атмосфере Земли с учетом кориолисовых сил), могут быть также полезны при описании турбулентных процессов в плазме (например, при описании плазмы непрерывной моделью или моделью кулоновски
взаимодействующих квазичастиц и заряженных «нитей», а также при изучении динамики альфвеновских вихрей в космической плазме [17, 21-24].
Благодарности. Работа выполнена за счет средств субсидии, выделенной в рамках государственной поддержки Казанского (Приволжского) федерального университета в целях повышения его конкурентоспособности среди ведущих мировых научно-образовательных центров, и при поддержке Национального научного фонда им. Шота Руставели (SRNF) (проект № RF17_252).
Литература
1. Белоцерковский О.М., Опарин А.М.. Численные эксперименты в турбулентности от порядка к хаосу. - М.: Наука, 2000. - 223 с.
2. Ламб Г. Гидродинамика. - М.: Гостехиздат, 1947. - 900 с.
3. Монин А.С., Яглом А.М. Статистическая гидромеханика: в 2 ч. - М.: Наука, 1967. -Ч. 1. - 720 с.; Ч. 2. - 638 с.
4. Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. 1. - М.: Наука, 1963. - 584 с.
5. Роуч П.Дж. Вычислительная гидродинамика. - М.: Мир, 1980. - 616 с.
6. Вычислительные методы в гидродинамике / Под. ред. Б. Олдер. - М.: Мир, 1967. -384 с.
7. Дим Г. Забуски Н. Стационарные V-состояния, их взаимодействие, возврат и разрушение // Солитоны в действии / Пер. c англ. под ред. А.В. Гапонова-Грехова. - М.: Мир, 1981. - С. 289-304.
8. Zabusky N.J., Hughes M.N., Roberts K.V. Contour dynamics for the Euler equations in two dimensions // J. Comput. Phys. - 1979. - V. 30, No 1. - P. 96-106. - doi: 10.1016/0021-9991(79)90089-5.
9. Поттер Д. Вычислительные методы в физике. - М.: Мир, 1975. - 392 с.
10. Березин Ю.А., Федорчук Н.П. Моделирование нестационарных плазменных процессов. - Новосибирск: Наука, 1993. - 357 с.
11. Белоцерковский С.М., Гиневский А.С. Моделирование турбулентных струй и следов на основе дискретных вихрей. - М.: Физ.-мат. лит., 1995. - 365 с.
12. Белашов В.Ю., Сингатулин Р.М. Алгоритм метода контурной динамики и моделирование вихревых структур. - Казань: КГЭУ, 2003. - 39 с. - Деп. ВИНИТИ 11.02.2003 г., № 272-В2003.
13. Overman E.A., Zabusky N.J. Coaxial scattering of Euler-equation translating V-states via contour dynamics // J. Fluid Mech. - 1982. - V. 125. - P. 187-202. - doi: 10.1017/S0022112082003309.
14. Baker G.R. A study of the numerical stability of the method of contour dynamics // Philos. Trans. R. Soc., A - 1990. - V. 333. - P. 391-400. - doi: 10.1098/rsta.1990.0167.
15. Belashov V.Yu., Singatulin R.M. Dynamics of vortex type wave structures in plasmas and fluids. Plasma Physics // AIP Conf. Proc. - 2003. - V. 669. - P. 609-612. - doi: 10.1063/1.1594004.
16. Белашов В.Ю., Харшиладзе О.А. Модифицированный метод КД и моделирование вихревых структур // Двенадцатая ежегод. конф. «Физика плазмы в Солнечной системе». - М.: ИКИ РАН, 2017. - С. 160.
17. Belashov V.Yu. Modeling of dynamics of vortex structures in continuous media // Astrophys. Aerosp. Technol. - 2016. - V. 4, No 3. - P. 28.
18. Сэффмэн Ф.Дж. Динамика вихрей. - М.: Науч. мир, 2000. - 376 с.
19. Belashov V. Yu., Belashova E.S., Kharshiladze O.A. Modified CD method and simulation of vortical structures in a plasma and fluids //J. Phys. Chem. Biophys. - 2018. - V. 8. -P. 61. - doi: 10.4172/2161-0398-C1-028.
20. Belashov V. Yu., Kharshiladze O.A. Numerical modeling of interaction of vortex structures in fluids and plasmas // VIII Annu. Meet. of the Georgian Mechanical Union. Book of Abstracts. - Tbilisi: Tbilisi Univ. Press, 2017. - P. 31-32.
21. Pokhotelov O.A., Stenflo L., Shukla P.K. Nonlinear structures in the Earth's magnetosphere and atmosphere // Plasma Phys. Rep. - 1996. - V. 22, No 10. - P. 852-863.
22. Belashov V. Yu. Interaction of N-vortex structures in a continuum, including atmosphere, hydrosphere and plasma // Adv. Space Res. - 2017. - V. 60, No 8. - P. 1878-1890. - doi: 10.1016/j.asr.2017.05.028.
23. Belashov V. Yu. Numerical study of interaction of vortex structures in plasmas and fluids // J. Phys. Chem. Biophys. - 2017. - V. 7, No 3, Suppl. - P. 49. - doi: 10.4172/2161-0398-C1-022.
24. Belashov V. Yu., Belashova E.S., Kharshiladze O.A. Nonlinear wave structures of the soliton and vortex types in complex continuous media: Theory, simulation, applications // Lecture Notes of TICMI. - Tbilisi: Tbilisi Univ. Press, 2018. - V. 18. - 90 p.
Поступила в редакцию 28.09.18
Белашов Василий Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник НИЛ исследований ближнего космоса Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: vybelashov@yahoo.com
Харшиладзе Олег Автандилович, доктор физико-математических наук, доцент, Центр космических исследований Института геофизики им. М. Нодия Тбилисский государственный университет им. И. Джавахишвили
ул. Чавчавадзе, д. 1, г. Тбилиси, 0179, Грузия E-mail: o.kharshiladze@mail.ru
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2019, vol. 161, no. 1, pp. 5-23
doi: 10.26907/2541-7746.2019.1.5-23
The Modified Method of Contour Dynamics and Modeling of Vortical Structures
V.Yu. Belashova*, O.A. Kharshiladzeb**
aKazan Federal University, Kazan, 420008 Russia bIvane Javakhishvili Tbilisi State University, Tbilisi, 0179 Georgia E-mail: * vybelashov@yahoo.com, ** o.kharshiladze@mail.ru
Received August 28, 2018 Abstract
This paper considers one of the most effective methods for modeling vortical structures, which are described by the 2-dimensional equation of carry of a vortex and by the Poisson equation for a flow function, namely, the contour dynamics method based on representation of a vortical stream by finite-area vortical regions. A modification of the contour dynamics method minimizing the errors arising at its direct application to description of vortical structures has been elaborated. The examples of the results of numerical experiments on the study of the dynamics of interaction of vortical structures for various configurations of their relative positioning, signs of vorticity, and distances between borders of the finite-area vortical regions have been presented.
Keywords: vortices, finite-area vortical regions, modeling, hydrodynamics, modified contour dynamics method, regimes of interaction, quasi-recurrence phenomenon, phase intermixing
Acknowledgments. The work is performed according to the Russian Government Program of Competitive Growth of Kazan Federal University and supported by the Shota Rustaveli National Science Foundation of Georgia (SRNF) (project no. RF17_252).
Figure Captions
Fig. 1. Illustration of the algorithm of the contour dynamics method. Basic designations.
Fig. 2. Three cases of formation of double curves.
Fig. 3. Interaction of two vortices as the initial distance between them equals: a) S = 2d; b) S = d/2.
Fig. 4. Interaction of the circular of finite-area vortical regions with different signs of the vorticity: Z = —1 and Z = 1.
Fig. 5. Interaction of four linearly disposed at t = 0 vortices with Zi = —1 for: a) Si = d/2; b) Si = 2d.
Fig. 6. Interaction of four finite area vortical regions with the same signs of the vorticity (Z = —1) for r = d/2 and r = d.
Fig. 7. Interaction of three vortices with the same signs of the vorticity (Z = —1) for r = d/2
Fig. 8. Interaction of four vortices for Si = d: a) Zi = Z3 > 0, Z2 = Z4 < 0; b) Zi = Z2 > 0, Z3 = Z4 < 0.
Fig. 9. Interaction of vortices with different orders of symmetry (m = 3 and m = 2), but with the same signs of vorticity (Z = — 1).
References
1. Belotserkovskii O.M., Oparin A.M. Chislennye eksperimenty v turbulentnosti ot poryadka k khaosu [Numerical Experiments in Turbulence from Order to Chaos]. Moscow, Nauka, 2000. 223 p. (In Russian)
2. Lamb H. Hydrodynamics. New York, Dover Publ., 1945. 735 p.
3. Monin A.S., Yaglom A.M. Statisticheskaya gidromekhanika [Statistical Hydromechanics]. Moscow, Nauka, 1967. Pt. 1, 720 p.; Pt. 2, 638 p. (In Russian)
4. Kochin N.E., Kibel' I.A., Roze N.V. Teoreticheskaya gidromekhanika [Theoretical Hydromechanics]. Pt.1. Moscow, Nauka, 1963. 584 p. (In Russian)
5. Roache P.J. Computational Fluid Dynamics. Hermosa, Albuquerque, NM, 1982. 612 p.
6. Fundamental Methods in Hydrodynamics. Alder B.J., Fernbach S., Rotenberg M. (Eds.). New York, Acad. Press, 1964. 386 p.
7. Deem G.S., Zabusky N.J. Stationary V-states, interactions, recurrence, and breaking. In: Solitons in Action. Lonngren K., Scott A. (Eds.). Academic, 1978, pp. 277-293.
8. Zabusky N.J., Hughes M.N., Roberts K.V. Contour dynamics for the Euler equations in two dimensions. J. Comput. Phys., 1979, vol. 30, no. 1, pp. 96-106. doi: 10.1016/0021-9991(79)90089-5. (In Russian)
9. Potter D. Computational Physics. London, John Wiley & Sons Ltd, 1973. 392 p.
10. Berezin Yu.A., Fedorchuk N.P. Modelirovanie nestatsionarnykh plazmennykh protsessov [Simulation of Nonstationary Plasma Processes]. Novosibirsk, Nauka, 1993. 357 p. (In Russian)
11. Belotserkovskii S.M., Ginevskii A.S. Modelirovanie turbulentnykh strui i sledov na osnove diskretnykh vikhrei [Modeling of Turbulent Jets and Traces on the Basis of Discrete Vortices]. Moscow, Fiz.-Mat. Lit., 1995. 365 p. (In Russian)
12. Belashov V.Yu., Singatulin R.M. Algoritm metoda konturnoi dinamiki i modelirovanie vikhrevykh struktur [Algorithm of the Contour Dynamics Method and Simulation of the Vortical Structures]. Kazan, KGEU, 2003, 39 p. Dep. VINITI on Feb. 11, 2003, no. 272-B2003. (In Russian)
13. Overman E.A., Zabusky N.J. Coaxial scattering of Euler-equation translating V-states via contour dynamics. J. Fluid Mech., 1982, vol. 125, pp. 187-202. doi: 10.1017/S0022112082003309.
14. Baker G.R. A study of the numerical stability of the method of contour dynamics. Philos. Trans. R. Soc., A, 1990, vol. 333, pp. 391-400. doi: 10.1098/rsta.1990.0167.
15. Belashov V.Yu., Singatulin R.M. Dynamics of vortex type wave structures in plasmas and fluids. Plasma Physics. AIP Conf. Proc., 2003, vol. 669, pp. 609-612. doi: 10.1063/1.1594004.
16. Belashov V.Yu., Kharshiladze O.A. A modified CD method and simulation of vortical structures. Dvenadtsataya ezhegod. konf. "Fizika plazmy v Solnechnoi sisteme" [Proc. 12th Annu. Conf. "Plasma Physics in the Solar System"]. Moscow, IKI Ross. Akad. Nauk, 2017. 160 p. (In Russian)
17. Belashov V.Yu. Modeling of dynamics of vortex structures in continuous media. Astro-phys. Aerosp. Technol., 2016, vol. 4, no. 3, p. 28.
18. Saffman P.G. Vortex Dynamics. Cambridge Univ. Press, 1993. 376 p.
19. Belashov V.Yu., Belashova E.S., Kharshiladze O.A. Modified CD method and simulation of vortical structures in a plasma and fluids. J. Phys. Chem. Biophys., 2018, vol. 8, p. 61. doi: 10.4172/2161-0398-C1-028.
20. Belashov V.Yu., Kharshiladze O.A. Numerical modeling of interaction of vortex structures in fluids and plasmas. Proc. VIII Annu. Meet, of the Georgian Mechanical Union. Book of Abstracts. Tbilisi, Tbilisi Univ. Press. 2017, pp. 31-32.
21. Pokhotelov O.A., Stenflo L., Shukla P.K. Nonlinear structures in the Earth's magnetosphere and atmosphere. Plasma Phys. Rep., 1996, vol. 22, no. 10, pp. 852-863.
22. Belashov V.Yu. Interaction of N-vortex structures in a continuum, including atmosphere, hydro-sphere and plasma. Adv. Space Res., 2017. V. 60. P. 1878-1890.
23. Belashov V.Yu. Numerical study of interaction of vortex structures in plasmas and fluids. J. Phys. Chem. Biophys, 2017, vol. 7, no. 3, suppl., p. 49. doi: 10.4172/2161-0398-C1-022.
24. Belashov V.Yu., Belashova E.S., Kharshiladze O.A. Nonlinear wave structures of the soliton and vortex types in complex continuous media: Theory, simulation, applications. In: Lecture Notes of TICMI. Tbilisi, Tbilisi Univ. Press, 2018, vol. 18. 90 p.
Для цитирования: Белашов В.Ю., Харшиладзе О.А. Модифицированный метод / контурной динамики и моделирование вихревых структур // Учен. зап. Казан. ун-та. \ Сер. Физ.-матем. науки. - 2019. - Т. 161, кн. 1. - С. 5-23. - doi: 10.26907/25417746.2019.1.5-23.
For citation: Belashov V.Yu., Kharshiladze O.A. The modified method of contour dynamics and modeling of vortical structures. Uchenye Zapiski Kazanskogo Universiteta. \ Seriya Fiziko-Matematicheskie Nauki, 2019, vol. 161, no. 1, pp. 5-23. doi: 10.26907/25417746.2019.1.5-23. (In Russian)
