CC BY
5Нелинейная динамика. 2011. Т. 7. № 1. С. 75-100. Полнотекстовая версия в свободном доступе http://nd.ics.org.ru
УДК: 551.465.7
М8С 2010: 37Б45, 76Ш5, 76В47, 76Е20
О лагранжевом переносе вблизи осциллирующего вихря в набегающем потоке
А. Е. Гледзер
Рассмотрено поведение жидкости в окрестности особой точки разрушенной нестационарными возмущениями сепаратрисы. Проводятся при различных частотах возмущений численные расчеты движения жидких частиц на первую половину оборота вокруг вихря. Предложено аналитическое представление движения отдельной частицы, откуда выводятся критерии захвата и выхода частиц. Анализируется соответствие полученных результатов с известным законом о ширине стохастического слоя. Показывается связь полученных аналитических выражений с мельниковским интегралом.
Ключевые слова: динамический хаос, вихревые структуры, стохастический слой
1. Введение
Значительный интерес представляют стационарные конфигурации взаимодействия однородного потока и точечных вихрей в присутствии твердого тела. Тогда могут возникать неподвижные особые точки, в которых производные от функции тока равны нулю. Через эти особые точки проходит сепаратриса, которая разделяет течение на области, друг от друга полностью отделенные в смысле обмена частицами жидкости. Когда появляются зависящие от времени возмущения, то эта картина полностью разрушается. Оказалось, что. используя лишь численный счет, разобраться в поведении частиц трудно, так как в нелинейной системе малые отклонения в начальных данных для частицы и малые изменения в частоте внешних колебаний могут заметно сказаться на траектории. Возникает, как указа-
Получено 2 декабря 2010 года После доработки 22 декабря 2010 года
Работа выполнена при финансовой поддержке программы Президиума РАН №17 «Математические методы нелинейной динамики» и РФФИ (грант 080500942 «Фундаментальные основы нелинейной динамики»).
Гледзер Алексей Евгеньевич aegledzer@gmail.com
Институт физики атмосферы им. А. М. Обухова PAH 119017, Россия, г. Москва, Пыжевский пер., д. 3
но в книге Пуанкаре [1], сложная структура двух кривых, которые заменяют разрушенную возмущениями сепаратрису, этакие «призраки» ее.
Ранние работы 90-х годов XX столетия активно используют интеграл Мельникова [2, 3] для установления частоты внешних возмущений, оптимальной для активации переноса лагранжевой примеси. Среди этих работ можно указать [4, 5], где исследуется течение, образующееся в результате натекания осциллирующим по величине скорости потоком на круговой цилиндр с присутствием в потоке точечного вихря, а также [6], которая исследует простейшее течение с отличной от нуля завихренностью, подвергнутое 5-образным периодическим возмущениям во времени. В ней уделяется в целом внимание толщине стохастического слоя, в котором и происходит обмен жидкостью между проточной и вихревой областями.
Где-то к концу 80-х или к началу 90-х годов XX столетия была предложена техника the lobe dynamics (см. книгу [7], статью [8]), которой в книге [9] присвоено название «динамика лепестков». Данная математическая теория рассматривает сегменты устойчивого и неустойчивого многообразий, которые расположены между двумя последовательными точками пересечений многообразий. Эти сегменты (похожие на мягкие лепестки цветка) ограничивают физические объемы жидкости. Далее производятся циклы сечения Пуанкаре, чтобы проследить эволюцию. Эта техника включает в себя различные теоремы, касающиеся вычислений обмена жидкими частицами на сечении Пуанкаре при помощи структуры устойчивого и неустойчивого многообразий, и результаты, касающиеся связи мельниковского интеграла и геометрических свойств «лепестка».
На основе данной техники в 90-х годах появились работы по моделированию более сложных систем, чем академические функции тока обтекания вихрей и диполей. Например, аналитическая кинематическая модель Bower (1991) о меандрирущем характере течения Гольфстрим и ее численное обоснование на основе уравнения эволюции баротропного ква-зигеострофического потенциального вихря на ^-плоскости с периодическими граничными условиями в зональном и меридиональном направлениях [10, 11]. Эти работы основываются на технике «лепестков». Статья [10] ограничивается фактически рассмотрением траекторий различных частиц и рисованием переходов «лепестков», помеченных для наглядности различными цветами, в течение двух периодов возмущения. Другой конфигурацией по исследованию хаотической динамики является известный в океанологии «двойной вихреворот» (wind-driven double gyre), индуцируемый ветром. Можно отметить [12] для трехслойной модели и [13] для уравнений «мелкой» воды в одном слое, в обоих случаях задается известное направление ветра, меняющееся по широте. В результате расчетов получаются циклоническое и антициклоническое направления движения, которые взаимодействуют и приводят в результате к меандрам. В зависимости от силы ветра можно получить тот или иной режим течения: периодическое течение, непериодическое и хаотическое. В работе [12] вычисляют устойчивое и неустойчивое многообразия, обсуждая попутно методы численного интегрирования исследуемой системы обыкновенных дифференциальных уравнений (Runge-Kutta-Fehlberg, Runge-Kutta четвертого порядка), рисуют «лепестки», показывают их динамику. В работе [13] дополнительно исследуются собственные числа для особой точки.
Особенно интересными для задач лабораторного моделирования атмосферных течений являются статьи [14, 15], которые посвящены исследованию баротропных возмущений зонального потока со сдвигом, в частности, Бикли-струи.
Из работ геофизического направления отметим [17-20] (см. обзор в книге [9]). В работе [16] исследуется течение в полукруглом бассейне, в котором на стационарную вихревую планетарно-топографическую составляющую накладывается нестационарный безвих-
ревой поток. Проводились численные расчеты трассеров и вычислялся показатель Ляпунова. В [17] изучался топографический вихрь над подводной возвышенностью при учете кориолисовой силы. Возвышение моделировалось осесимметричным гауссовым рельефом. Полная функция тока представляется в виде суммы колеблющегося во времени однородного потока и планетарно-топографической составляющей. Образуется единственная гиперболическая особая точка. Если внешнее течение подвергнуто малому колебанию, то будет происходить хаотическая адвекция. Авторы вычисляют функцию Мельникова (в [17] указывается на невозможность ее аналитического представления), детально разъясняется динамика маркеров, найдена оптимальная частота перемешивания. Далее рассмотрена с помощью сечения Пуанкаре хаотизация ядра вихревой области, проиллюстрировано применение теории КАМ-торов. В статьях [18, 19] рассмотрен вихрь над 5-образной горой, похожее течение рассматривается в работе [8] в качестве объекта демонстрации расчетов с помощью функции Мельникова и теории «лепестков». Кроме этого (см. [9]) исследуют еще резонансы, КАМ-торы и хаотическое «море», анализируется струя Бикли, т. е. сложный вариант меандрирующей струи, как в работе [11].
В указанных работах [18-20], в частности, используется метод, изложенный в статье [21], для оценки толщины стохастического слоя вблизи критической точки, для которого применение интеграла Мельникова затруднительно. При этом предложены поправки для уравнения, полученного в [21], с учетом применения метода не только для простейшего гомоклинического соединения, но и для наличия двух критических точек, т. е. для гетеро-клинического соединения.
В настоящей работе с учетом подхода работ [18, 19] проведена более развитая, чем в этих работах, модификация уравнений, выведенных в [21]. Некоторые результаты были получены в [22, 23]. Далее результаты сравниваются с высокоточным прямым численным расчетом уравнений лагранжевой динамики и с простейшим приближением из [21]. Установлена связь этого приближения с интегралом Мельникова.
2. Представление нестационарных вихрей
Рассмотрим лагранжеву динамику в двумерном течении несжимаемой жидкости, описываемой зависящей от времени функцией тока ф(p, в, t) (в полярных координатах x = = pcos в, y = p sin в):
dв дф dp і дф ,
= в? ш = и‘> = -т' {2Л)
где функция ф задана двумя формами:
ф = и(1)р 1--------т- sin0 + ii(t)ln
V p J W (2.2)
U (t) = U0 U (t), U (t) = 1 + euU (ut); K (t) = KqK (t), К (t) = 1 + et K (ut);
ф = A(t)e-p2 — psin в, A(t) = A0 + eaA(ut). (2.3)
Функция тока (2.2) соответствует обтеканию нестационарным однородным полем скорости со значением ux = —U(t) на бесконечности кругового цилиндра радиуса ао, причем вокруг цилиндра наложена также осциллирующая вблизи значения Ко циркуляция К (t). В (2.2) амплитуды возмущений eu, et много меньше единицы (см. [24]). Течение с функцией тока (2.3) рассматривалось в работе [6], где вводилась 5-образная структура возмущений.
Введена частота и внешних возмущений. В (2.2), (2.3) функции U(£), К(£), а(£) заданы далее как sin £ или cos £ с точностью до знака, где £ = ut. В данной работе, по сравнению с [4, 18], используются периодические вариации циркуляции, а не только скорости набегающего потока. Это сделано для общности рассмотрения. Кроме этого, нельзя исключать, что возможны изменения циркуляции, например, для океанских вихрей из-за бароклинных эффектов. Преимущество течения (2.3) по сравнению с (2.2) в более физически обоснованном убывании на бесконечности вихревой компоненты из-за экспоненты, т. е. в быстром переходе к однородному внешнему потоку. Течение (2.2), в отличие от (2.3), потенциально во всем пространстве, кроме точки р = 0. Течение (2.1), (2.2) при ао = 0 называют точечным вихрем, оно ему и соответствует, а динамику течения (2.1), (2.3) будем называть «гауссовым» вихрем, согласно терминологии статьи [25].
Для невозмущенного течения (2.2) (eu = et = 0) при 2аоUo < Ко имеется критическая
Kq 1 + \/1 — {2oqUo/Ко)2 о
точка с полярными координатами р = ро = 77-------------------------, 0 = тг7Г, в которой
Цо 2 2
скорость равна 0. В дальнейшем для (2.2) рассматривается точечный вихрь, для которого
а,о = 0. Для невозмущенного течения (2.3) тоже существует при Ао > \/ехр(1)/2 ~ 1.16 критическая точка 0 = -7Г на некотором зависящем от Ао расстоянии от начала координат.
Вводим для точечного вихря (2.1), (2.2) переменную ф = sin в, безразмерное время T =
= t ■ частоту возмущений Q = (при этом QT = uit) и радиус г = Из (2.1), (2.2)
Ро Uo Р0
имеем уравнения
фё'£ = ^ М
Из первого уравнения в (2.4) видно, что при малых возмущениях, когда U(O.T) > 0, изменение r монотонно между положениями частицы ф = ±1. Тогда, перейдя в (2.4) от T к независимой переменной г, получим
dtp = _<£_ _ к [ОГ(г)] ^ = sgn (cos в) 1 (9 м
dr Г U[QT(r)} г2’ dr .^1 - и[ПТ(г)}'
Для «гауссова» вихря (2.1), (2.3) имеем уравнения (ф = sin в, r = р), тоже в итоге основывающиеся на монотонности поведения г, когда sgn (cos в) фиксирован:
dr de ~ . _r2
— = cos в, г— = —2A(ujt)re — sm 0, dt dt
^ + |+2ЛИ(г)]в-'’=0, * <2'6)
dr 7 dr _ цр.
Общая картина невозмущенного движения частиц (линий тока) для точечного вихря (2.1), (2.2) представлена на рисунке 1a. Течение (2.1), (2.3) выглядит аналогично, но центр вихря rc не совпадает с началом координат O, и сами линии тока очень быстро выпрямляются по мере удаления от центра (см. рис. 1b).
При eu = 0, et = 0, ea = 0 в течениях имеются сепаратрисы, отделяющие движение в вихре от движения из бесконечности в бесконечность вокруг вихря. Сепаратриса в каждом
Рис. 1. Общий вид течения точечного вихря (2.1), (2.2) (а) и «гауссова» вихря (2.1), (2.3) при Ао = = 2 (Ь).
из течений замыкается на критическую точку Гнсг = 1> днсг = Для (2.2) и Гнсг ~ 1.27707,
3
днсг = 9^ Для (2.3) при Ао = 2. В случае течения (2.3) величина Гнсг определяется из второ-
£4
го уравнения в (2.6) как решение уравнения Е(Тнсг) = 2АоТисге-Гнсг = 1. Далее определим для каждого течения в области замкнутых линий тока в невозмущенном случае вблизи сепаратрисы значения величины «нижней» точки Тн, где линия тока пересекает ось симметрии при ф = —1, и величины «верхней» точки ТВ, где линия тока пересекает ось симметрии
при ф = +1. Для течения (2.2) тогда можно записать Тн = 1 +1п Тн + 5, Тн 1 — л/26 + |<5,
ТВ + 1 + 1п ТВ + 5 = 0, ТВ/тВсг ~ 1 — 5/(1 + ТВсг)+0.552/(1 + ТВст)2, ТВсг 0.2785. Величина 5 мала и положительна, она характеризует отклонение значения функции тока от ее значения — 1 на сепаратрисе. Для течения (2.3) можно записать Е(тн) 1 + £ (2Тнсг — 1/тнсг),
е = Тнсг — тн > 0, ТВсг & 0.22913.
Для дальнейшего понадобится вычисленный в виде рядов период т движения частицы в области собственно вихря около сепаратрисы. Вычисления подобного рода характеристик проводились, например, в [15, 18, 19]. Метод вычисления приведен в [21]. Результаты для течений (2.2) и (2.3) соответственно:
, 4 ч N
— ~ ZJ^l + 1п ( — | , ZJ^l = 1п (1 — Г Всг) + --(1 — ГВсгУ\
2 \е/ п=1 П (2.7)
е = 2л/26 >0, N 23, гм ^ -0.81141940;
2гнс
2Г%сг - 1
ZN - ІП (ГЕст ■ є)
N
а
ZN — ІП (Г2Нсг - г'всг) (гЯсг “ гВсг)
п=1
2 \п гв ст)
(2.8)
Є = ГЕст - ГЕ > о, N - 1500, ZN & +0.85927181.
Все ряды в этих формулах асимптотические, коэффициенты ап, Сп вычислялись по специальной программе. На примере течения (2.2) покажем, как была получена формула (2.7). Значение полупериода т/2 представлялось суммой интегралов:
т* тн
I = I В{г)) ■ (1,Г) + J 0(7]) ■ (1,Г), £>(т?)
тв
1
1 + 5 + 1п Г] V
-1/2
Первый соответствует участку линии тока, удаленному от критической точки, второй интеграл относится к окрестности критической точки. Выберем произвольную величину Т* таким образом, что 1 1 — г* л/§. Тогда второй интеграл легко вычисляется аналити-
чески, и в итоге получим:
1-тв г 2"| -1/2
! / + •^ + 1п2 + 1п(1-г*)-1п(|)
1-т* <-
Оставшийся интеграл, ввиду того, что |1 — тв| < 1, может быть представлен в следующем виде:
1 тв
1 т*
N
1
N
1 + 2-^Си”-2
п=3
п=2
-1/2
% £ ’
где С\ =0; Сп = ^ ^ + ... + п = 2,..а при п —> +оо имеется асимптотика Сп = С +
+ 1пп — 1, С = 0.5772157. В итоге для этого интеграла получим, разлагая далее по формуле Тейлора:
1-тв N
[ (1 +0-1 • С + • • • ■ ^Г- ~ 1п(1 - гв) - 1п(1 - г*) + Е — • (1 - Гв)
гі* * п=1 п
При подстановке этого выражения в формулу для т/2 член 1п(1 — т*) сокращается.
Результаты сравнения численных расчетов по (2.2) и (2.3) полупериода движения частицы в стационарных течениях и по формулам (2.7), (2.8) в интересующем далее диапазоне е приводятся на рис. 2, где для уравнений (2.2) (линии 1а и 2а) и уравнений (2.3) (линии 1Ь и 2Ь) показаны полупериоды т/2 в зависимости от е. Линии 2а, 2Ь получены с помощью расчетов по аппроксимациям (2.7), (2.8) соответственно. Из рис. 2 видно, что аппроксимации воспроизводят числовой расчет в окрестности сепаратрисы при малых е, что особенно важно для дальнейшего.
2
п
4 5|------------1------1-----1 г 1-----1 I |—|----------------------------------------------------------------------------1-]-?-1-1-III]-1-1-1-1-1-1—I—г
т/2
Рис. 2. Расчеты полупериода обращения частиц по замкнутым линиям тока в области вихря численно (линии 1а, 1Ь) и аналитически по (2.7) (линия 2а) для невозмущенного точечного вихря (2.1), (2.2) и по (2.8) (линия 2Ь) для невозмущенного «гауссова» вихря (2.1), (2.3).
3. «Невышедшие» траектории
Рассмотрим решение уравнений (2.5), (2.6) с начальными данными г = в = т. е.
из точки на оси симметрии вблизи сепаратрисы у верхней части вихря (рис. 1). В случаях обхода центра вихря против часовой стрелки на половину оборота sgn (cos в) = — 1.
Начальные условия для (2.5), (2.6): ф(^) = +1, T(п) = 0, t(r) = 0.
Представление решения уравнений (2.5) в виде ряда по возмущениям до второго порядка eu, et имеет вид
Г
ф) = _1 + 5 + lnr _ I . f
r r I (3.1)
^ _ U[QT(r)} , _ K[QT(r)} , _2 U2 [QT(r)] _ _ U [QT(r)} К [QT(r)}
j{r) — —£u-----------h Sk----------r £u---------------£u£k----------------}
r r r r
где параметр 5 определяется начальной точкой движения r = ri, ф = 1, а входящая в (3.1) величина T(r) — время перемещения частицы из положения ri в r — определяется вторым уравнением в (2.5). Для второго рассматриваемого течения («гауссов» вихрь) (2.3), (2.6) вычисления аналогичны. Решение уравнения для ф(г) в (2.6) имеет вид
Г
ф(г)г = Аое-г2 — ea j f (n) dn — ф, f (r) = 2re-r2a [ut(r)], (3.2)
ri
где постоянная ф определяется начальным положением частицы r = ri, ф = 1. Рассмотрим более подробно течение (2.1), (2.2), поскольку оно окажется сложнее «гауссова» вихря (2.3).
Пусть начальная точка Т1 находится в области замкнутых линий тока (собственно вих-
3
ря), которая ограничена сепаратрисой при еи = 0, ьд. = 0. Если при в = ^7Г (ф = —1) траектория возмущенного движения пересекает ось симметрии в точке т = г, так что ф(г) = — 1, то далее частица останется в вихре на данном обороте. Это можно показать из анализа скоростей в декартовой системе координат для течения (2.1), (2.3).
Для течения (2.1), (2.2) пересечение траекторией нижней части оси симметрии в = 3п
= — не означает, что частица далее монотонно по углу достигнет какой-то малой окрестности начальной точки Тi. На рис. 3 показана часть траектории частицы вблизи в = = Щ- (у критической ТОЧКИ ГЦсг = 1)) полученной при счете уравнений точечного вихря (2.1), (2.2) для частоты возмущений П = л/2 при еи = ед. = е8 = 0.08, и (ПТ) = = +008 (ПТ), К (ПТ) = + йп (ПТ), ег = твсг — ^ = 5.38 • 10-3. При движении слева направо траектория из-за возмущений делает петлю, горизонтальный размер которой порядка 0.96 • (4.74623 — 4.67965) ~ 0.06 (см. рис. 3). Цифры у траектории показывают: значения
угла в в радианах вблизи в = Щ- ~ 4.71, расстояние г до начала координат и безразмерное
время Т движения от точки г = Гг, в = Пунктиром показана сходящаяся в Гнсг = 1 часть
невозмущенной сепаратрисы ^в целом весь рисунок — это сильное увеличение окрестности
точки г = 1, в = . Можно показать, что для «гауссова» вихря (2.1), (2.3) такие петли
невозможны.
Гпгг - 1 г =0.9788
Г= 5.4514
Рис. 3. Пример поведения частицы для точечного вихря (2.1), (2.2) в окрестности гиперболической особой точки. Безразмерная частота П = \/2, функции возмущения и (ПТ) = со в(ПТ), к(ПТ) = = вш(ПТ), амплитуды возмущения еи = ек = 0.08. Начальные координаты частицы в^ = п/2,
Тг = ТВсг — ег, ег = 5.38 • 10-3.
Если решение уравнения ф(г) = —1 не существует, то частица не пересекает ось сим-3
метрии при 9 = т^7г и выходит из области циркуляционного движения.
Решение г уравнения ф(г) = —1, которое для (3.1) записывается в виде
Г
1 + 1пг + 5 + | / (п) йп — г = 0, (3.3)
Гг
можно искать как отклонение от величины «нижней» точки гн вида г = гн (1 + ГГ).
Можно из (3.3) получить квадратное относительно а уравнение (условимся обозначать это как первое приближение)
С(а) = о? Ан + аВн + Сн = 0, Ан = 1/2,
ГН
Вн = — 1+Гн — Гн/(гн), Сн = —Бн, Бн = J /(г) йг,
(3.4)
которое использовано в работах [18, 19, 21] для получения критериев выхода частиц из вихря и вычисления ширины области захвата частиц вихрем. При вычислении /(г) по (3.1) следует для (3.4) использовать только члены первого порядка по еи, £и, а входящие в /(г) значения времени Т(г) брать на невозмущенной траектории. Таким образом, входящая в /(гн) величина Т(гн) равна полупериоду движения т/2 частицы по линии тока в невозмущенном случае. Конечно, на самом деле это недостаточно точно, поскольку при возмущениях. особенно значительных, траектория частицы изменяется, а значение радиуса гн частица будет иметь при ином времени, чем т/2. Другими словами, Т(гн) может заметно отличаться
от т/2. В решении уравнения (3.4) а = —Вн + sgn(Вн) ^— 2Сн член Вц в подкоренном выражении имеет второй порядок по возмущениям еи, ей в отличие от слагаемого первого порядка 2Сн, которое и определяет величину а ~ -\/|ьд.|, т. е. величину отклонения г
от гн. При этом отрицательный знак дискриминанта В^—2Сн < 0 свидетельствует для этого приближения о невозможности существования решения Г уравнения (3.3). Выбор знака перед радикалом в выражении (3.4) для а соответствует отсутствию выхода частицы из вихря при нулевых возмущения. Этот анализ показывает, что в (3.4) квадратичным слагаемым по возмущению а2Ан нельзя пренебрегать при аппроксимации решения уравнения (3.3).
Рассмотрим более точную, чем (3.4), аппроксимацию решения уравнения (3.3). Теперь уже целесообразно учитывать и квадратичные по возмущениям члены для /(г) в представлении (3.1). Но это вовсе не главный вклад в поправку. Поскольку радиус Г пересечения траекторией оси симметрии (в случае его существования) находится вблизи критической
ТОЧКИ Г Нет = 1 при 0 = ^7Г, то во втором уравнении (2.5) имеется особенность \/1 — ф2(г)
в знаменателе, она дает быстрый рост функции Т(г) при г ^ Г. В малой окрестности решения Г, т. е. при г = Г(1 + ц), ц << 1, используя представление для <^(г) по (3.1), можно
2
получить с точностью до квадратичных членов по /л :
1 — ^>2(г) & к1ц + к2л2, Л = г/г — 1, (3.5)
где
к1 ■ г2 = +2Г(Г — 1) — 2Г2/(Г),
к2 ■ Г2 = +Г2 + (Гг — 1) ■ (1 — 4Г) — 2/(Г)Г + 4/(Г)Г2 — /2(Г)Г2.
При этом вблизи г ~ 1 имеем К2 ~ 1 и положительно, а |к 11 достаточно мало (см. ниже рисунок 8с).
Производная по г функции / (г) в (3.1) имеет вид
//(г) = _/М+_^гіГп = _Лг)+ д/ П 1
= Фгед(г) +
г д {ПТ) гіг г д (ПТ) и [ПТ] л/1 - ф2{г)
Ф(г)1 (г)
л/КіЦ + К.2^2 '
где
д/
(3.7)
Т = Т(г), Ф Гед(г) = -£Р-, Ф(Г) = +-^^П, 1(г) = 1/и[ПТ(г)}. (3.8)
подразумевается дифференцирование выражения /(г) в (3.1)
Под обозначением
а (ПТ)
по всему аргументу ПТ. Отсюда приближенно имеем:
г /л
/(г) - /(г) = J /'(п)йп ~ Фгеа(г) (г - ?) + Ф(?)1 (?)? J
йи
\/«1V + /4:2
(3.9)
К )2
Последний интеграл в (3.9) легко вычисляется заменой г/ = —^ ^— и равен
+
^п(^)
\/«2
1п№)|-1п|£(0)|), ^(^() = кі^(0) = кь (3.10)
Используем результаты (3.9), (3.10) для приближенного представления интеграла в уравнении (3.3) от п до Г, используя в качестве основы интеграл от функции возмущения в пределах от п до гн. Последний интеграл снова берется на невозмущенной траектории. Это дает ошибки при вычислении интеграла (см. ниже рис. 9).
Введем цн = гн/Г — 1 = —г/ (1 + сг). В итоге:
Г Гн г Гн
! /(г) йг = ! /(г) йг + У /(г) йг = ! /(г) йг + /(?) (? - гн) - Фгед(Г)
(г - гн)2
— sgn(Z )Ф(?)
гн
I (г) г
(3.11)
н
л/«2
(1 + й)2
- цн 1п(0)| +
ГЩ-Г(О)
4к2
Здесь
іп\г([і)\ = ДГІГ^, вд = +я(^) 1п\г{^)\-1 + 1{ц) іп\г([і)\ +1
Выделяя в (3.3) часть квадратного уравнения (3.4), используя (3.11), перепишем уравнение (3.3) в виде
й3 г2
С(5) - у + [-/(?) + /(гн)] г на + &і„д(г)-£ +
+ е^)Щ^гк (1 + а) [+з ь 2(0) +(! + «)г(№) - г(0)
\/«2
4К2
(3.12)
0.
2
2
Отметим, что в (3.12) был выделен еще и кубический по й член в разложении 1п (1 + й), как это было сделано в работах [18, 19] для несколько другой задачи.
В (3.12) слагаемое С(й) определяется согласно (3.4), но функция / (гн), в отличие от (3.4), содержит члены также и второй степени малости по еи, е^. Эти же члены используются в (3.12) и для /(й), /(гн) по формуле (3.1), и при вычислении Фгед(й), Ф(й) согласно (3.8), а также для величин к\ и К2 согласно (3.6). Для коэффициента К2 допустимо использовать только члены первого порядка по еи, ей. Необходимо подчеркнуть, что параметр й достаточно мал, однако в (3.12) в последнем слагаемом нельзя проводить разложение
в ряд функции Г(— а _), поскольку первая производная функции Г(р) в нуле сократится 1 + й
со слагаемым от 1п \Z(0)|, а вторая производная функции Г(р) в нуле имеет особенность.
Важнейшим элементом уравнения (3.12) является неявное присутствие времен Т(гн), Т(й). Время Т(й) вычислялось на возмущенной траектории. Для этого использовалась вышеуказанная методика представления 1 — ^>2(г) в малой окрестности решения й, так что
Поскольку при в = 3п/2 величина г достигает экстремума (максимума), то г < й, следовательно, р = г/й — 1 < 0. Если в расчетах получится, что < 0, и в выражении (3.10) для Z(р) взять знак минус перед радикалом, то, так как К2 > 0, выражение для Z(р) будет знакоопределенное ^п^) = —1). По мере изменения величины р в нем от 0 при г = = й до произвольного отрицательного значения при г < й функция Z (р) будет монотонно увеличиваться по модулю, так как р увеличивается по модулю. Если же в формуле (3.13) \Z(рн)\ > \Z(р)\, то, с учетом, что sgn(Z) = —1, получим Т(г) > Т(гн). Как и должно быть, если жидкая частица движется от положения гн до г. Так как рн = гн/й — 1, можно заключить, что \Z(рн)\ > \Z(р)\ при рн < р < 0. Значит, величина гн < г. Таким образом, если К1 < 0, и взять знак минус перед радикалом для Z(р) в формуле (3.10), то г > гн, Т(г) > Т(гн), как и должно быть.
Для вычисления Т(й) по формуле (3.13) (т. е. при р = 0) требуются К1, К2, а они, согласно (3.6), имеют в своих выражениях величину функции возмущения /(й), для вычисления которой по формуле (3.1) необходимо знать время Т(й). Поэтому вычисления по формуле (3.13) проводились итерационно: подставлялось отправное значение времени в выражения (3.6), которое выбиралось, чтобы была сходимость, затем вычислялось новое значение времени по (3.13). Далее процесс повторялся до установления.
Для второго рассматриваемого течения («гауссов» вихрь) (2.3), (2.6) вычисления аналогичны. С возмущением в (3.2) уравнение <^(й) = —1 имеет вид
Используя представление г = гн (1 + си), из (3.14) получим квадратное уравнение первого приближения:
Г
с г]Т й г
Г(г)=Г(гн) + J ^^»Г(гн)-8ёп(адг)—[+1п|^(^н)|-1п|аду. (3.13)
/
г
(3.14)
С((й) — йй2Ан + аВн + Сн — 0, Ан — Аог“н (2гн — 1) е Гн
(3.15)
Все комментарии, относящиеся к квадратному уравнению (3.4) для течения (2.1), (2.2), справедливы и здесь. Но следует учесть, что, согласно (3.2), в этом течении нет членов второго и высшего порядков малости по амплитуде возмущения.
В следующем приближении вычисления аналогичны (3.5)—(3.13), и вместо (3.12), (3.13) имеем уравнения
G(a) + а3 ( 1 — | г2н) ГцЕ(гн) +ea[-f(r) + /(гн)} г на + £аа2ФГед(г)Чг+
H
Ф(г)г
+ sgn(Z)ea——=— (1 + а) +аIn |^(0)| + (1 + а) \/«2
Г(рн) - г(0)
Фгед(г) = 2 (1 — 2r2) e r a [ut(r)], Ф(г) = 2re r a' [ut(r)] ■ u, a'(£) =
= 0, (3.16)
MO
/dt ~ г
— dr) = t(rH) -sgn(Z)-^= +\n\Z{nH)\ — In \ Z(p)\
Г H
Аналог формулы (3.5) имеет для течения (2.3), (2.6) вид
1 — p2(r) W тр + К2р2, р = r/r — 1,
(3.17)
(3.18)
где
К1 = +2 ■ [1 — E(r)] — 2f (r),
К2 = —3 ■ [1 — E(r)] + E(r) ■ [2Г2 — E(r) + ea2f (r) [2 — E(r)] — ef2(r).
(3.19)
Имеются те же оценки для коэффициентов при й ~ г нет, т. е. при Е(й) ~ 1: К2 ~ 2й2 — 1 > 0, \к1\ мал.
4. «Незахваченные» траектории
Рассмотрим теперь случай, когда частица, стартовав из точки на оси симметрии в = = ^ вблизи сепаратрисы вне ограниченной ею области собственно вихря, т. е. при п > Гвсг
(см. рисунок 1а, Ь), не пересекает ось симметрии при в = Щ-, т. е. не захватывается вихрем. В невозмущенном движении для такой частицы на линии тока имеется точка с радиусом
dp
dp dO
rm, где -f = cos в—
drdr
dp dr
= 0 (рис. 1). При наличии возмущений такая же точка, где
= 0, тоже имеется, но с радиусом г = rm (1 + а).
Для течения точечного вихря (2.1), (2.2) в отсутствие возмущений уравнение ^
= 0 дает простейшее решение при 5 < 0: rm = exp (—5) > 1, pm = p(rm) = —1/rm > —1. В нестационарном течении уравнение ^ =0 эквивалентно
Г
ln г + 5 + J f (r) dr — f (r)r = 0,
(4.1)
2
r=r
m
что требует снова анализировать интеграл от функции возмущения. В отличие от задачи (3.3), допустима аппроксимация подынтегральной функции в виде /(г) / (гт) +
+ /'(гт)(г — гт), поскольку фт = 1 и входящая в /'(гт) производная от времени ограничена. Первое приближение из (4.1) приводит к квадратному уравнению:
а
Ат + аВт + Ст — 0, Ат — — — (l + r%nf'(rm)')
Г m
Bm = 1 rm f (rm), Cm = rmfm + Sm, Sm = J f (r) dr.
(4.2)
В отличие от уравнения (3.4), в (4.2) наличие квадратичного члена не является обязательным, он лишь слегка поправляет линейную часть, поскольку Bm ~ 1. Решение уравнения (4.2) при малых возмущениях существует всегда. Знак перед радикалом в решении этого квадратного уравнения определяется однозначно по факту обязательного пролета при отсутствии возмущений. Однако одного лишь решения уравнения (4.2) еще недостаточно для критерия пролета. Необходимо также, чтобы при а, определяемом по значению а из (4.2), величина |р(а)| = | sin O\ < 1. В противном случае найденный а = rm (1 + а) не отвечает какой-либо точке на возмущенной траектории частицы.
В итоге указанное неравенство и уравнение (4.2) дают в рассматриваемом приближении следующий критерий пролета частицы:
(rm — 1) — Sm + O [ — f (rm)rm + Sm] + <°2
+
1 - f'(rm)rm
2
+ f (rm)rm Sm
> О. (4.3)
В численных расчетах в этом приближении были учтены квадратичные поправки по єи, є к в выражениях для f (гт) и f'(гт), хотя это не обязательно. В самом грубом упрощении этот критерий можно записать как 1 — гт + Бт < 0.
Можно также получить для уравнения (4.1) приближение, аналогичное выписанному выше (3.12), (3.13). Останавливаться на этом не будем.
5. Результаты расчетов
Численные расчеты, представленные далее, проводились для приближений первого (уравнения (3.4), (3.15)) и второго ((3.12)—(3.13), (3.16)—(3.17)) порядков. Результаты для них сравнивались с прямым численным счетом уравнений (2.1)—(2.3), который для конечного времени порядка полупериода т/2 может быть достаточно точным (метод Рунге-Кутты-Мерсона пятого порядка с переменным шагом). Начальное значение шага было 10-5, допустимый предельно малый шаг 10-9, погрешность 2 • 10-5.
Для «гауссова» вихря (2.1), (2.3) были проведены расчеты для пяти частот возмущения и = 0.2, \/2, 1.8, 2.0, 2.5 в диапазоне 0 < еа < 0.3 и начальных положений частицы 10-4 < ег = гвсг — гг < 5 • 10-2. Вся рассматриваемая область параметров еа, ег (см. рисунок 4а, Ь) разбивается на область малых возмущений еа при относительно больших ег (т. е. относительно больших удалений начальной точки гг от сепаратрисы у ее верхнего максимума гвст) ниже линии «В» на рисунке 4а, Ь, когда частица не выходит из вихря, и область относительно больших возмущений еа при малых ег — выше линии «В» (в этом случае частица выходит из вихря).
Рис. 4. (а) Результаты расчетов в параметрах ег, еа для «гауссова» вихря (2.1), (2.3) Ао = 2.0 траекторий частиц, первоначально находившихся в вихревой области на оси ординат, численно по (2.6) и аналитически по критерию (3.15). Частота возмущения и; = \/2, функция возмущения А(ш1) = = — 8ш(иЬ). В области «ЕЕ» численный и аналитический результаты совпадают и предсказывают выход частицы из вихря. В области отсутствия выхода относительная погрешность аеч для 1 от 0% до 0.5%, 2 — 0.5% ^ 0.75%, 3 — 0.75% ^ 1.25%, 4 — 1.25% ^ 1.75%, 5 — 1.75% ^ 2.0%, 6 — 2.0% ^ 3.0%, 7 — 3.0% + 4.0%, 8 — 4.0% + 5.0%, 9 — от 5.0%. (Ь) Результаты расчетов для «гауссова» вихря (2.1), (2.3) Ао = 2.0 по уравнению (3.16) с вычисленным по (3.17) значением Ь(й). Сравнение точности с уравнением (3.15) по параметру ■ для 1 величина 7^ меняется в пределах от 0% до 10%, для 2 — 10% + 20%, для 3 — 20% + 30%, для 4 — 30% + 80%, для 5 — 80% +100%, для 6 — более 100%.
На рис. 4а показаны результаты численного счета уравнений (2.1), (2.3) при ш = \/2, А (иЬ) = — 8ш(шЬ), Ао = 2 и приближения квадратного уравнения (3.15). В расчете было задействовано 128 х 128 точек на плоскости параметров еа, ег. В ней выделяются две главные области, где численный счет по уравнениям (2.6) (то же, что и (2.1), (2.3)) и критерий (3.15) показывают выход частицы из вихря (область «ЕЕ») и отсутствие выхода. В об-
ласти, где численный расчет по (2.6) и критерий (3.15) на рис. 4а показывают отсутствие выхода из вихря, отмечены области относительной погрешности aeq = \ (Fe — Fq) /Fe| = = | (Fe — Fq) / (Fe — Гн) | определения F = Fq по уравнению (3.15) по отношению к эталонному Fe из (2.6). На рисунке 4а область 1 соответствует погрешности ueq от 0% до 0.5%, 2 — от 0.5% до 0.75%, 3 — от 0.75% до 1.25%, 4 - от 1.25% до 1.75%, 5 - от 1.75% до 2.0%, 6 — от 2.0% до 3.0%, 7 — от 3.0% до 4.0%, 8 — от 4.0% до 5.0%, 9 — от 5.0%. В области, указанной на рисунке 4а символом «I», численный счет (2.6) давал отсутствие выхода из вихря (ниже линии «В»), а критерий (3.15) показал выход частицы, т. е. в этой области критерий квадратного уравнения (3.15) оказался несостоятельным.
На рисунке 4b показаны результаты расчетов F = Ft по второму приближению (3.16) с вычислением времени t(F) по формуле (3.17). Рисунок сравнивает эталонное решение Fe с Ft по отношению к разности |Fe — Fq| (Fq найдено из квадратного уравнения (3.15)) с помощью введения характеристики Yetq = |Fe — Ft|/|Fe — Fq|. Области 1-5 на рисунке 4b соответствуют лучшей точности по уравнению (3.16) по отношнению к квадратному уравнению
(3.15): Yetq меняется в области 1 в пределах от 0% до 10%, в 2 — от 10% до 20%, в 3 — 20% до 30%, в 4 — 30% до 80%, в 5 — 80% до 100%. В области 6 (где %tq > 100%) квадратное уравнение (3.15) работает лучше. В области «I» численный расчет давал выход частиц из вихря, а приближение (3.16), (3.17) — отсутствие выхода. В области «II» результаты по уравнению (3.16) и численному счету по (2.6) показывают отсутствие выхода частиц из вихря, в отличие от уравнения (3.15). В области «ЕЕ» численный счет и уравнения (3.16). (3.17) давали выход частиц.
При увеличении частоты и область выхода «ЕЕ» постепенно «наступает» на область отсутствия выхода, устанавливая неизменную границу. Все анализы до и = 1.8 показали хорошее совпадение предложенной теории — уравнения (3.15), (3.16), и непосредственным численным расчетом. Однако расхождение увеличивается по мере увеличения частоты возмущения. При и = 2.0 оно существенно, а при и = 2.5 становится относительно большим. Отметим, что вычисления по уравнению (3.16) различаются в зависимости от того, какое время t(F) подставлялось в уравнение: либо из численного решения, либо вычисляемое по формуле (3.17), либо вовсе из невозмущенного решения. Расчеты со временем из невозмущенного решения оказались заметно хуже двух других вариантов, как в определении положения линии «В», так и в точности вычисления Ft. Оказалось, что положения граничной линии «В» лучше всего определяется по уравнению (3.16) с расчетом времени по (3.17), а точность определения Fпо сравнению с (3.15) улучшается заметнее всего, если взять время из численного решения.
Для точечного вихря (2.1), (2.2) результаты зависят не только от безразмерной частоты возмущений Q, но и от того, сколько параметров из £u, £k в (3.1) представляют возмущения. Проведены расчеты для частот Q = 0.2, 0.8, 1.0, 1.2, \/2, 1.5, \/3, 1.8. Диапазон параметров еи, £k был от 0.0 до 0.1. При возмущении es = £u, £k = 0 функция U (QT) = +cos (QT). При es = еи = £k функции U (QT) = + cos(QT), K (QT) = + sin(QT). При es = £k, £u = 0
функция K (QT) = — cos (QT).
На рисунке 5а показаны результаты численного счета уравнений точечного вихря (2.1),
(2.2) (Q = 1.2, £s = £и, £k = 0, 10_4 < £r < 10_2) и расчетов по уравнению первого приближения (3.4). Разбиение области то же, что и для «гауссова» вихря (рисунок 4а): область 1 соответствует погрешности aeq от 0% до 1.5%, 2 — от 1.5% до 2.25%, 3 — от 2.25% до 3.25%, 4 — от 3.25% до 5.0%, 5 — от 5.0%. Показана также область «I», в которой результаты по численному счету (нет выхода) и по первому приближению (3.4) (выход) не совпадают.
Рис. 5. (а) Результаты расчетов в параметрах ег, е8 для точечного вихря (2.1), (2.2) траекторий частиц, первоначально находившихся в вихревой области на оси ординат, численно по (2.2) и аналитически по критерию (3.4). Частота возмущения П = 1.2, функция возмущения и (ПТ) = ссв(ПТ). В области «ЕЕ» численный и аналитический результаты совпадают и предсказывают выход частицы из вихря. В области отсутствия выхода относительная погрешность аед для 1 от 0% до 1.5%, 2 — 1.5% + 2.25%, 3 — 2.25% + 3.25%, 4 — 3.25% + 5.0%, 5 — от 5.0%. (Ь) Результаты расчетов для точечного вихря (2.1), (2.2) по уравнению (3.12) с вычисленным по (3.13) значением Т(Г). Частота возмущения П = 0.2, функция возмущения и (ПТ) = ссв(ПТ). Сравнение точности с уравнением (3.4) по параметру ^е1(1■ для 1 величина 7^ меняется от 0% до 10%, 2 — 10% + 20%, 3 — 20% + 30%, 4 — 30% ^ 100%, 5 — более 100%.
На рисунке 5Ь для частоты П = 0.2 и параметров возмущений е3 = ек, еи = 0 показаны результаты расчетов величины по уравнению (3.12) с временем Т(г) из формулы (3.13): область 1 — величина ^ел^ меняется от 0% до 10%, 2 — от 10% до 20%, 3 — от 20% до 30%, 4 — от 30% до 100%. В области 5 квадратное уравнение (3.4) работает лучше. Области несовпадения результатов обозначены так же, как и на рис. 4Ь: в области «I» численный счет давал отсутствие выхода, а уравнения (3.12), (3.13) — выход, в области «II» результаты (3.12),
(3.13) и численный счет совпадают и показывают отсутствие выхода, в отличие от уравнения первого приближения (3.4). В области «ЕЕ» численный счет и уравнения (3.12), (3.13) давали выход частиц.
При увеличении частоты область выхода из вихря (большие є3, малые єг) постепенно «отступает» от области отсутствия выхода (большие єг, малые є3), в отличие от «гауссова» вихря.
В работе [21] для точечного вихря (2.1), (2.2) были получены результаты при частотах возмущений О = 0.2, 1.0, УЗ при ь3 = єи = Єкі и (ПТ) = +соз(ОТ), К (ПТ) = + 8Іп(ПТ). Однако фактическая область расчета была меньше: 10_4 < єг < 6.5 • 10_3, 0 < є3 < 0.06, а сами результаты представлялись не посредством координаты єг = Гвсг — г і, а через отклонение г і от гн, т. е. через ширину стохастического слоя, которой соответствовала зависимость л/єі. Отметим, что эти координаты использовались и в работе [20].
Рассмотренные до сих пор результаты использовали небольшие величины безразмерной частоты О и один из параметров возмущения єи или єк, который был меньше 0.1. Здесь можно говорить об удовлетворительном совпадении результатов по квадратному уравнению (3.4) и по уравнению (3.12) (с вычислением времени Т(?) по (3.13)) со значениями из численного счета (2.1), (2.2). Однако при увеличении частоты возмущений О (О > 1.2) и наличии двух параметров возмущений є3 = єи = єк результаты в значительной степени ухудшаются, хотя предсказания обоими приближениями о выходе частицы из вихря или отсутствии выхода в значительной части области параметров є3, єг совпадают с полученным при численном счете. Заметим также, что для частот О > 0.8 вычисления по уравнениям (3.12), (3.13) не улучшают точности по сравнению с квадратным уравнением (3.4).
На рисунке 6а, Ь в области «I» между ломаными линиями теоретические расчеты по квадратному уравнению (3.4) (рис. 6а) и по уравнению (3.12) с вычислением времени Т(?) по (3.13) (рис. 6Ь) не совпадают с численным решением системы (2.1), (2.2), хотя вне этих областей («ЕЕ», «пЕпЕ») имеется совпадение, кроме узкой области «II» вблизи оси є3 на рисунке 6Ь. Для рисунка 6а, Ь О = 1.5, но имелся только один параметр є3 = єк, который менялся от 0 до 0.1. Если имеется и другой параметр возмущения, т. е. є3 = єи = єк, то области несовпадения «I» по квадратному уравнению (3.4) и по уравнениям (3.12), (3.13) с численным счетом значительно расширяются. Заметим, что при расчетах уравнений (3.12), (3.13) в некоторых областях наблюдалась перемежаемость небольших участков, где результаты совпадали, и рядом с ними участков, где были несовпадения. В остальной части диаграммы («ЕЕ», «пЕпЕ») результаты совпадали.
Для результатов, показанных на рисунке 4а, Ь, можно привести другую кумулятивную форму представления. Для области (например, на рисунке 4Ь ниже линии «В»), в которой численные расчеты и вычисления по уравнениям (3.15), (3.16), (3.17) (на сетке 16 х 16) показывают отсутствие выхода частиц из вихря, проводится нумерация точек расчета єа. єг так, как показано на вставке к рисунку 7а. В основной части рисунка по оси абсцисс отложены значения ае численного решения системы (2.1), (2.3) для каждой рассматриваемой точки (всего точек 67). Эти точки отмечены вверху рисунка 7а на вставке к диаграмме. По оси ординат показаны величины ае/ая (окружности) и (?е/а (крестики), т. е. отношения численных значений а к вычисленным по квадратному уравнению (3.15) (окружности) и по уравнениям (3.16), (3.17) (крестики). Из рисунка видно, что погрешность увеличивается с ростом значения ае, при этом для ае < 0.1 вычисления по уравнениям (3.16), (3.17) лучше соответствуют численному счету (отношение ае/а ближе к 1). При ае > 0.1, как видно из рисунка 7а, наблюдается систематическое уменьшение значения по квадратному уравнению ад по сравнению с численным ае, и увеличение а по сравнению с ае (окружности уходят вверх, а крестики вниз).
Рис. 6. Результаты расчетов при еи =0, е8 = ей для точечного вихря (2.1), (2.2) численно (а), (Ь) и по уравнению (3.4) (а), по уравнению (3.12) со временем Т(г) из уравнения (3.13) (Ь). Частота возмущения П = 1.5, функция возмущения К (ПТ) = — ссв(ПТ).
Рисунок 8а показывает для тех же точек, что и на рисунке 7а, относительную роль различных слагаемых в уравнении (3.16). По оси абсцисс отложены номера точек, показанных на вставке к рисунку 7а. Здесь 1 — функция О(а) с учетом возмущений второго порядка, 2 — слагаемое в (3.16), содержащее произведение первой степени а и разницы ( —/(г) + /(гн)), 3 — последние слагаемые левых частей, происходящие от Ф(г), 4 — кубический по а член, 5 — квадратичный по а член, имеющий Фгед(г). Сумма всех этих членов в каждой точке равна нулю, так как а — решение уравнения (3.16). Из графика видно, что члены 2 и 3 в (3.16) имеют наибольшую амплитуду, но они практически компенсируют друг друга. Члены 1 и 4 — квадратичная функция О(а) и кубическая по а поправка, в сумме также почти дают нуль, т. е. они могут адекватно определять решение а, что, учитывая его
Рис. 7. (а) Отклонения значений решений ач уравнения (3.15) (окружности), значений решений 5г уравнения (3.16) (крестики) от значений ае, полученных в численном счете. Все данные соответствуют рис. 4. Вставка — нумерация расчетных точек в области отсутствия выхода по рис. 4Ь. (Ь) Увеличенная часть рис. 7а при ае < 0.1.
1 1 0 1 1
— - -
- 1 1 1
.э! 1. 1 I------------------------------------1 I I------------------------------------
О 10 20 30 40 50 60 70
Рис. 8. (а) Относительная роль различных слагаемых в уравнении (3.16). Все данные по рис. 4. (Б) Сравнение времени £(г) по формуле (3.17), из численного решения и из невозмущенной системы. (о) Величины «4 и «2, рассчитанные по формулам (3.19).
малость, объясняет неплохую возможность простейшей (квадратное уравнение с поправками, как в [18, 19]) аппроксимации уравнения (3.14). Это также объясняет точность этого приближения, что и демонстрируют рисунки 4а, 5а.
Рисунки 8Ь, о, а также 9а, Ь имеют технический характер и демонстрируют оправданность предположений, взятых за основу при выводе приведенных выше аппроксимаций.
Рисунок 8Ь показывает вычисление времени Ь(г) по формуле (3.17) («гауссов» вихрь) (кривая 1), из численного решения (кривая 2) и из невозмущенной системы (кривая 3). Видно, что формула (3.17) хорошо работает в большинстве рассмотренных точек плоскости еа, ег.
Рисунок 8с показывает величины к\ (кривая 1) и К2 (кривая 2), найденных по соответствующим формулам (3.19).
Следующие рисунки иллюстрируют точность вычисления функции возмущения /(г) для «гауссова» вихря из (3.2) и ее интеграла 3(г) = /(п)^п для фиксированной точки
на плоскости еа, ег. На сетке 16 х 16 была выбрана точка еа = 8 • 10_2, ег = 3.66933 • 10_2 (номер 12 по ег, номер 5 по еа, номер 52 в ряду рассмотренных на рис. 7а), в которой точность применения (3.16) с вычислением по (3.17) времени Ь(г) была лучше по сравнению с квадратным уравнением (3.15) (\те — г\/\Г —гг\ ~ 16%). Результаты расчетов ряда величин для этой точки приводятся в приложении.
На рис. 9а показан график /(г) из численного решения (2.6) (кривая 1) по формуле, аналогичной (3.9), приведенной в приложении как (А1), (кривая 2). Линия 3 показывает постоянную /(гн), которая используется в уравнении (3.15). Все три графика строились в малой окрестности точек соответственно ге для кривой 1, Г* для кривой 2, тя для кривой 3. Видно, что кривые 1 и 2 вблизи правого конца указывают на бесконечность производной функции /(г) в точках ге, Г соответственно. Рисунок 9Ь показывает расчет 3(г) в тех же окрестностях численно по (2.6) (кривая 1), по формуле, аналогичной (3.11), приведенной
приложении как (А2), (кривая 2). Для уравнения (3.15) получится прямая линия, поскольку выполняется интегрирование постоянной (кривая 3). Из рисунков видна удовлетворительная аппроксимация кривых 1 линиями 2.
-0.4
■0.405
-0.410
-0.415
-0.420
-0.425
-0.430
-0.435
-0.440
-0.445
1.1
|. 1 . (Ь)
2>\
. . ! ■ ■ -
01
1.02
1.03
1.04
1.05
1.06
1.07
1.08
1.0?
1.10 1.11
Рис. 9. (а) Сравнение графиков /(г) из численного решения (2.6), по формуле, аналогичной (3.9) ((А1) в приложении), и постоянной /(гн), которая используется в уравнении (3.15). (Ь) Сравнение расчетов интеграла 3(г) из численного решения (2.6), по формуле, аналогичной (3.11) ((А2) в приложении), и для уравнения (3.15).
Для точечного вихря (2.1), (2.2) при старте частицы выше точки гвсг невозмущенной сепаратрисы на оси 9 = ^, так что ег = — Гвсг > 0 (см. рис. 1а), были проведены
для П = 1.2, е3 = еи, ек = 0, и (ПТ) = +008 (ПТ), 10-4 < ег < 10-2 расчеты численно и с использованием уравнений (4.2) и (4.3). Оказалось, что имеется адекватное определение границы захвата, но точность определения довольно низкая.
6. Заключение
В статье решается задача о траектории жидкой частицы, которая перемещается вблизи разрушенной сепаратрисы в гамильтоновой системе, подвергнутой малым периодическим нестационарным возмущениям. В отсутствие возмущений траектория представлялась бы замкнутой линией тока. Какова судьба жидкой частицы при возмущениях? В результате получены два приближения, совершенно различных по своей технологии, сложности, объему записи и вычислениям. Оба они могут ответить на поставленный вопрос, конкурируя друг с другом. Отметим, что приведенный в статье метод может быть применен не только для рассматриваемых течений (2.2), (2.3), но и для многих других конфигураций.
Сложное приближение (3.12), (3.16) для рассматриваемых двух течений в целом дает лучшую точность вычисления отклонения от критической точки при пересечении нижней части оси симметрии при в = 3п/2. Об этом свидетельствуют рисунки 4Ь, 5Ь, где проведено сравнение с результатом, полученным из простых приближений квадратного уравнения (3.4) и (3.15) (области 1-5 на рисунке 4Ь, 1-4 на рисунке 5Ь).
Однако при относительно больших величинах возмущений еа и е3 или начальных удалений ег квадратные уравнения (3.4), (3.15) могут давать лучшую точность вычисления указанного расстояния. Особенно заметно это проявляется вблизи линии раздела «В» в параметрах еа (или е3), ег, отделяющей области выхода и отсутствия выхода частицы. Это можно объяснить тем, что именно вблизи этой границы, как показывают кривые 2 и 3 на рисунке 8а, выделяющие часть членов в приближенном уравнении (3.16), амплитуда этих членов возрастает (у номеров точек 27-28, 47-48, 61-62 — см. вставку к рисунку 7а). Это, учитывая противоположность знаков слагаемых 2 и 3, приводит к ошибкам вычисления отклонения траектории частицы вблизи критической точки. Заметим, что для квадратичных членов в (3.16) (плюс кубическая поправка) такого возрастания амплитуд вблизи кривой «В» не наблюдается. Также возрастают ошибки вычисления времени Ь(Г) по аналитической формуле (3.17) (кривая 1 на рисунке 8Ь) в сравнении с численными значениями (линия 2 на рисунке 8Ь). А именно значение Ь(Г) (или Т(Г)) определяет точность вычисления интеграла 3(г) в уравнениях (3.3) или (3.14) (рисунок 9Ь). Собственно, недостаточная точность вычисления времени Ь(г) (Т(Г)) по невозмущенному движению (кривая 3 на рисунке 8Ь) и приводит к ухудшению предсказания по квадратным уравнениям (3.4) и (3.15) и к необходимости строить улучшеннные приближения (3.12), (3.13) или (3.16), (3.17).
Тем не менее, в качестве первого приближения простейшие уравнения (3.4), (3.15) могут работать, и это вызывает удивление. Вероятно, это связано с тем, что эти приближения связаны с критерием Мельникова разрушения сепаратрисы. Чтобы продемонстрировать это, возьмем пример системы, для которой мельниковское расстояние между устойчивым и неустойчивым многообразиями легко вычисляется (для рассматриваемых уравнений (2.1),
(2.2), (2.3) его пришлось бы рассчитывать численно). Это пример нелинейного осциллятора с трением и внешней накачкой (см. [26]):
= р, = — х + х2 — ур + /Зсов(Ш). (6.1)
аЬ аЬ
р2 х2
Невозмущенная система имеет гамильтониан Н(х,р) = — + Общий вид изолиний
гамильтониана соответствует рассмотренным выше конфигурациям, так как в системе имеется одна гиперболическая особая точка Хнсг = 1, р = 0. Сепаратриса Нсг = 1/6 проходит через нее и точку Хвсг = —1/2, р = 0. Уравнения (6.1) не описывают какое-либо физическое течение жидкости, хотя невозмущенная система бездивергентна.
Расстояние между указанными сепаратрисами имеет вид интеграла
+ГС)
д = у (+ая [тт]/да<)1РМ|„ + 01,(тт«) , (м)
где /(Х,Р,Ь) = 0, д(Х,Р,Ь) = 7Р + в008 (П^), а X(Ь), Р(Ь) — решение невозмущенных уравнений, отвечающих движению по петле сепаратрисы из седла в седло:
= совЬ (* + 0) - 2 = 3 эшЬ (£ + 0)
С08Ь(( + 0) + 1' ^совЬ(; + в) + 12
где 0 — произвольная постоянная, она отвечает за начальное положение частицы при Ь = 0 на сепаратрисе. Если начальное положение X(0 ) = Хвсг, Р(0 ) = 0, то 0 = 0.
Интеграл может быть вычислен (см. [26]):
Д(0) = —^J Р2(і) ' йі + /З І Р(і) ■ Ш8 (Ш) ■ йі, (6.3)
или, в конечном виде.
Д(Є) = -|7+/}.6іг!}2 8Іп(П0)
5 81пИ (Пп)
Согласно критерию Мельникова, при знакопеременности величины Д(0), когда второе слагаемое больше первого, имеет место возникновение хаоса, т. е. пересечение в бесконечном числе точек устойчивого и неустойчивого многообразий (или устойчивой и неустойчивой возмущенных сепаратрис, согласно [26]). Видно, что трение ^ > 0 и внешнее возмущение противоборствуют. Чем больше 7, тем больше необходимо задать величину в, чтобы состоялся хаос.
С другой стороны, имеется возможность столь же легко провести решение задачи (6.1) в первом и втором приближениях, аналогичных вышеприведенным. Выпишем только первое приближение:
С(а) = а2 Ан + аБн + Сн = 0, Ан = хН (2хн — 1) ~ 1 — 4е,
Вн = +хН (хн — 1) + вхн 008 [ПЬ(хн)] — 7ХнРо(хн) ~ £ + в008 [ПЬ(хн)] — 7Р0(хн),
хн хы (6-4)
Сн = +2
З сов [Ш(х)] ■ йх — у р0(х) ■ йх
где
Ро{х) = у +2Я - х2 + ^ • ж3, ро{хн) = О,
хн ~ 1 —л/25 — — • 8, хв ~ хвсг + — -8 + — • 82, 8 = Нсг— Н > 0, хвсг = —1/2, Нсг = 1/6.
Как и выше в (3.4), (3.15), под і(х) подразумевается невозмущенное время движения частицы из точки Хі = хв, задаваемой отклонением от точки хвсг на сепаратрисе с помощью разницы 8 значений гамильтониана, в точку х.
Видно, что свободный член Сн квадратного уравнения можно представить в виде
^-= 0 J Р(1) ■ сое (П • I) ■ <И — 7 J Р2{1) ■ <И & А(©)/2, о о
поскольку йх = Р(£) • й£, Р(£) = р0(Х(£)).
, ~ - Вн ± л/Вя _ 4А# С#
гешение уравнения (6.4], имеющее вид о; =-------------------------, может не существо-
2Ан
вать при Сн > 0 (поскольку величина Б‘Н имеет второй порядок малости по коэффициентам 7, в, е; это же наблюдается для квадратных уравнений (3.4) и (3.15)), что совпадает с критерием Мельникова — величина Д(0) > 0 при 0 = 0. Отсутствие решения при Сн > 0 означает, что частица, находившаяся вначале в вихре, выходит из него и не пересекает ось симметрии.
На рис. 10 показаны расчеты в параметрах 10_4 < ех = хг — хвсг < 5 • 10_2, 0 < в < 0.2 системы уравнений (6.1) численно и по формуле (6.4). Частота возмущений П = 0.5, параметр трения 7 = 0.1. Показана область «I» несовпадения расчетов (численно — выход, аналитически — нет выхода). Линия «В» определяет по численному решению границу выхода и отсутствия выхода.
Рис. 10. Результаты расчетов в параметрах ех, в системы уравнений (6.1) траекторий частиц, первоначально находившихся в области замкнутых траекторий на оси х, численно и аналитически по (6.4). Частота возмущения П = 0.5, параметр трения 7 = 0.1. В области «EE» численный и аналитический результаты показывают выход частицы. В области «nEnE» численный и аналитический результаты показывают отсутствие выхода. В области «I» численный счет дает выход, а теория по (6.4) — отсутствие выхода.
Автор благодарен О. Г. Чхетиани за интерес к работе по данной тематике и обсуждение результатов.
7. Приложение
Для «гауссова» вихря (2.1), (2.3) можно получить аналоги формул (3.9) и (3.11) соответственно:
/(г) = /(г) + Фге9(г) (г-г)+ sgn(Z)^^ (Ь\г(ц)\ - 1п |/н|); (А1)
f(r])d,r] = f(r)rH а + {1 + а)ц + Фгед{г)-^- -а2 + (1 + а)2 ц
22
rH
- sgn(Z)
4«2\/«2
(1 + а) Г(уН) - Г(у)
- sgn(Z)
Ф{¥)г% л/«2
(1 + а) +а + (1 +а) у ln |Z(0)|.
(A2)
Результаты расчетов величин для точки с координатами єа = В • 10 2, єг = 3.66933 x
х 10 2 для рис. 8, 9: тНсг = 1.2770446, тБсг = 0.2291321; тн = 1.0019324, тб = 0.1924388; дискриминант В2Н - 4АНСн = 8.94 • 10“2; ае = 8.88452671 • 10“2, ад = 8.68436798 • 10“2, а* = 8.85629457 • 10“2; те = 1.090949, гд = 1.088944, Т* = 1.090667; £(?) численно 2.06964, £(?) с помощью формулы (3.17) равно 2.14323, по невозмущенной траектории численным методом 1.66198; к\ = -0.643892, К2 = 2.371433. Численные величины, соответствующие кривым на рис. 8а для этой точки: (1) -5.11754 • 10_4, (2) -3.18715 • 10_3, (3) 3.33033 • 10_3,
(4) 3.39391 • 10“4, (5) 2.91844 • 10“5. Сумма их всех равна -1.91124 • 10“13.
r
2
Список литературы
[1] Пуанкаре А. Избранные труды: Т. 1,2. М.: Наука, 1971, 1972. 773 с.; 1001с.
[2] Мельников В. К. О поведении траекторий системы близкой к автономной гамильтоновой системе // Докл. АН СССР, 1962, т. 142, №3, с.542-545.
[3] Мельников В. К. Качественное описание сильного резонанса в нелинейной системе // Докл. АН СССР, 1963, т. 148, №6, с. 1257-1260.
[4] Kadtke J. B., Novikov T. F. Chaotic capture of vortices by a moving body: I. The single point vortex case // Chaos, 1993, vol.3, no. 4, pp. 543-553.
[5] Luithardt H.H., Kadtke J.B., Pedrizzetti G. Chaotic capture of vortices by a moving body: II. Bound pair model // Chaos, 1994, vol. 4, no. 4, pp. 681-691.
[6] Weiss J. B. Hamiltonian maps and transport in structured fluids // Phys. D, 1994, vol. 76, pp. 230-
238.
[7] Wiggins S. Chaotic transport in dynamical systems. New York: Springer, 1992. 301 p.
[8] Rom-Kedar V., Leonard A., Wiggins S. An analytical study of transport, mixing, and chaos in
an unsteady vortical flow // J. Fluid Mech., 1990, vol. 214, pp. 347-394.
[9] Кошель К. В., Пранц С. В. Хаотическая адвекция в океане. М.-Ижевск: НИЦ «Регулярная и хаотическая динамика», 2008. 360 с.
[10] Lozier M.S., Pratt L.J., Rogerson A.M., Miller P. D. Exchange geometry revealed by float trajectories in the Gulf Stream // J. Phys. Oceanogr., 1997, vol. 27, pp. 2327-2341.
[11] Rogerson A. M., Miller P. D., Pratt L. J., Jones C. K. R. T. Lagrangian motion and fluid exchange in
a barotropic meandering jet // J. Phys. Oceanogr., 1999, vol. 29, pp. 2635-2655.
[12] Coulliette C., Wiggins S. Intergyre transport in a wind-driven, quasigeostrophic double gyre: An application of lobe dynamics // Nonlinear Process. Geophys., 2001, vol. 8, pp. 69-94.
[13] Poje A. C., Haller G. Geometry of cross-stream mixing in a double-gyre ocean model // J. Phys.
Oceanogr., 1999, vol. 29, pp. 1649-1665.
[14] Del-Castillo-Negrete D., Morrison P. J. Chaotic transport by Rossby waves in shear flows // Phys. Fluids A, 1993, vol. 5, no. 4, pp. 948-965.
[15] Del-Castillo-Negrete D. Asymmetric transport and non-Gaussian statistics of passive scalars in vortices in shear // Phys. Fluids, 1998, vol. 10, no. 3, pp. 576-594.
[16] Козлов В.Ф., Кошель К. В. Об одной модели хаотического переноса в баротропном фоновом течении // Изв. РАН, ФАО, 2000, т. 36, №1, с. 119-128.
[17] Козлов В. Ф., Кошель К. В. Некоторые особенности хаотизации пульсирующего баротропного потока над осесимметричной подводной возвышенностью // Изв. РАН, ФАО, 2001, т. 37, №3, с. 378-389.
[18] Козлов В.Ф., Кошель К. В., Степанов Д. В. О влиянии границы на хаотическую адвекцию в баротропных квазигеострофических моделях фоновых течений // Дальневост. матем. журн., 2001, т. 2, №2, с. 89-98.
[19] Козлов В. Ф., Кошель К. В., Степанов Д. В. Влияние границы на хаотическую адвекцию в простейшей модели топографического вихря // Изв. РАН, ФАО, 2005, т. 41, №2, с. 242-252.
[20] Ryzhov E., Koshel K., Stepanov D. Background current concept and chaotic advection in an oceanic vortex flow // Theoret. Comput. Fluid Dyn., 2010, vol. 24, pp. 59-64.
[21] Гледзер А. Е. Захват и высвобождение массы в вихревых структурах океана // Изв. РАН, ФАО, 1999, т. 35, №6, с. 838-845.
[22] Гледзер А. Е. Динамика конечномерных геофизических течений и влияние различных механизмов вязкости: Дис. на соискание уч. ст. канд. ф.-м. н. Москва: ИФА РАН им. А. М. Обухова, 2004.
[23] Гледзер А. Е. Динамика лагранжевых частиц около колеблющегося вихря над плоской поверхностью // Междунар. конф. «Потоки и структуры в жидкостях: Физика геосфер» (Москва, ИПМех РАН, 23-27.06.2009): Сб. тезисов: Ч. 1. Москва, 2009. С.67-68.
[24] Милн-Томсон Л.М. Теоретическая гидродинамика. М.: Мир, 1964. 660с.
[25] Максименко Н. А., Орлов О. И. Интегральные характеристики ядра квазистационарного «гауссова» вихря в однородном или сдвиговом потоках // Океанология, 1991, т. 31, вып. 1, с. 34-41.
[26] Кузнецов С. П. Динамический хаос. М.: Физматлит, 2001. 296 с.
On the Lagrangian transport near oscillating vortex in running flow
Alexey E. Gledzer
A. M. Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences Pyzhevsky 3, Moscow, 119017 Russia aegledzer@gmail.com
Passive particles advection is considered in the vicinity of hyperbolic stationary point of the separatrix destroyed by insteady perturbations. For different frequencies of the disturbancies the trajectories of advected particles are investigated analytically and numerically. The approximate criteria of capture and release of particles are obtained. The results are linked with known law for the stochastic layer width near separatrix. The obtained criteria are connected with analytical Melnikov’s integral.
MSC 2010: 37D45, 76U05, 76B47, 76E20
Keywords: chaotic dynamics, vortex structures, stochastic layer
Received December 2, 2010, accepted December 22, 2010
Citation: Rus. J. Nonlin. Dyn., 2011, vol. 7, no. 1, pp. 75-100 (Russian)
