Автомодельное решение нелинейного уравнения диффузии для спектральной плотности энергии турбулентности Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 19, № 1, 2014

Автомодельное решение нелинейного уравнения диффузии для спектральной плотности энергии турбулентности*

В.Н. Гребенёв1, С. В. Назаренко2, И. В. Шваб1,4, Ю. А. Чиркунов1,

Г. Г. Лазарева3'4, О. В. Штырина1'4, С. Б. Медведев1 1 Институт вычислительных технологий СО РАН, Новосибирск, Россия 2Институт математики университета г. Варвик, Англия

3Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия 4Новосибирский государственный университет, Россия e-mail: [email protected]

Изучается нелинейное вырождающееся уравнение диффузии для феноменологического описания турбулентности, ассоциированное с моделью Лейта. Дано аналитическое обоснование существования автомодельного режима для спектральной плотности энергии турбулентности Е(к, Ь) в пространстве волновых чисел к и проведено численное исследование поведения траекторий в фазовом пространстве соответствующей динамической системы. Рассматриваемое уравнение построено таким образом, что возникают два стационарных решения: колмогоровский спектр, соответствующий каскадному состоянию, и термодинамическое распределение, которое устанавливается в системе. Стационарное состояние в данной модели состоит из "нелинейной смеси" постоянного потока и термодинамической компоненты. Автомодельный режим модели реализуется как автомодельное решение второго рода, для которого формирование спектра на больших волновых числах происходит за конечное время, что было показано ранее в результате численных экспериментов в случае исчезающей вязкости и отсутствия внешних сил воздействия.

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

Введение

Уравнение, предложенное в [1], записывается в виде, близком к модели Лейта [2]:

где £ — время, к — абсолютное значение волнового числа, f — внешняя сила, V — коэффициент кинематической вязкости. Спектр энергии турбулентности Е(к,£) нормирован

*Работа выполнена при финансовой поддержке РФФИ (гранты № 12-01-00234-а и 12-01-00648-а) и ФЦП "Научные и научно-педагогические кадры инновационной России на 2012-2013 годы" (контракт 8504).

(1)

таким образом, что плотность кинетической энергии равна у Ес1к. Отметим, что уравнение (1) является частным случаем семейства уравнений, изучаемых в работе [3], и при отсутствии внешних сил и диссипативных членов имеет в качестве стационарного решения в дополнение к колмогоровскому спектру Е(к) = р2/3к-5/3 термодинамический спектр Е(к) = Q2/3k2, где Р и Q постоянны. Общее стационарное решение уравнения является нелинейной комбинацией этих спектров. Как показано в [1], стационарный спектр возникает либо при наличии внешних сил, либо в случае, когда начальные условия сосредоточены в окрестности достаточно малых значений волновых чисел к. При этом происходит распространение фронта плотности кинетической энергии турбулентности в направлении больших к ^ то и за фронтом устанавливается степенной колмогоровский спектр.

Введённое в рассмотрение автомодельное решение (см. [1]) уравнения (1) в случае исчезающей вязкости и отсутствия внешний сил есть следствие допустимости уравнением двухпараметрической группы растяжений, особенностью которой является невозможность вычисления её параметров на основе законов сохранения и анализа размерности задачи. Данный тип автомодельных решений относится к понятию автомодельно-сти второго рода, введённому в пионерных работах Л.Д. Ландау и К.П. Станюковича, а также Г. Гудерлея [4-6]. Отметим, что уравнение (1) близко связано с уравнением нелинейной диффузии в неоднородной среде

,дм 8 1 адп\

1 у 1 '« = ау( 1 Уц)' (2)

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

Сохранение рассматриваемой моделью симметрийных свойств уравнений Навье — Стокса делают её полезной при исследовании формирования промежуточного асимптотического режима, найденного численно в [1] для так называемого явления "тёплого" каскада, возникающего в турбулентности, который в рамках уравнения (1) описывается изучаемым автомодельным решением.

Существование автомодельного решения

Рассмотрим уравнение (1) в случае исчезающей вязкости и отсутствия внешних сил воздействия на неоднородную среду

дЕ 1 д ( 11/2 1 д . , 2,

■ж = 8 дк к /2Е2 дк <Е/к2>

Уравнение (3) допускает следующее представление решения в автомодельном виде:

E = (t* - t)aF(n), n = k/k*, k* = c(t* - t)b, (4)

где a, b, c — константы такие, что a, c — положительные, а b — отрицательная величины, при этом a и b удовлетворяют соотношению

a = -2 - 3b. (5)

Тогда уравнение для функции F запишем как

dF C3/2 d / nn i d , „ %

(3b + 2)F + ^ = "Г dn (n'1/2F 1 ¡Sí (n-2F)

и будем искать F в классе неотрицательных функций, где она ведёт себя как F0n-x при b2

П ^ 0 и как — (1 - n)2 при n ^ 1, где x = -a/b, а F0 — произвольная положительная константа. Подходящим масштабным преобразованием константа c и коэффициент C3/2/8 могут быть приравнены к единице, и тогда (6) примет следующий вид, где a и b выражены в терминах параметра x:

х^з (nf + xF) = dn (,11/2f2J (n-2F)) . (7)

Для уравнения (7) получаем нелинейную краевую задачу на собственные значения, в которой необходимо найти значение параметра x вместе с функцией F, удовлетворяющей вышеприведённым краевым условиям. Уравнение (7) может быть сведено к автономному уравнению с помощью подстановки

1 о 0 dF 3 . F=25п щ (8)

где f (s) и g (s) — функции переменной s = ln n, что позволяет получить на них систему дифференциальных уравнений первого порядка

df 3 7 = 2(f+»>•

f dS = 3 í5f2 + 6/9 - 9»2 + T-J*3» + xf)

Далее преобразованием

d 1 d

dS = /d? р(т) =f(s), a(T) =g(s)

система (9) сводится к виду, пригодному для её последующего изучения в рамках теории

динамических систем:

dp 3

Тт = 2Р +а),

da 1 / 2 2 10 / ч . / ч

— = - 5р2 + 6pa - 9a2 +--(3a + xp) |. (10)

ят 3 \ x — 3

Состояниями равновесия системы (10) на фазовой плоскости (р, а) являются следующие точки фазовой плоскости: Р1 = (0,0) , Р2 = (0,10/3(х — 3)) и Р3 = (1, —1). Сразу отметим, что искомая траектория системы расположена в квадранте р > 0, а < 0. Переменной п ^ 0 соответствует т ^ —то, и так как Р ~ п-х, СР/Сп ~ —хп-х-1 при П ^ 0 получаем f ~ 5п(3-х)/2 и д ~ — (5/3)хп(3-х)/2. Отметим также, что значение искомого параметра х принадлежит интервалу (0, 3). Граница фронта автомодельного решения (6) соответствует значению п =1 или т ^ то. Фактически границей фронта автомодельного решения является точка равновесия Р2, где f ~ —5/2(6п3/2(1 — п)) и д ~ 56/3 (см. [8]). Точка равновесия Р3 соответствует точному решению уравнения (7), равному Р = (1/25)п-3. Дифференциал правой части системы (10) имеет следующий вид:

(

Д(р,а)

3р + 2 а

1 10х

3 10р + 6а + х-^^

Тогда

Д(0, 0)

0

10х

р

^ 6р — 18а +

0 10

\

30

х3

11)

3(х — 3) х — 3

и точка Р1 = (0, 0) является вырожденной особой точкой линеаризованной системы с собственными значениями А1 = 10/(х — 3), Л2 = 0 и собственными векторами

б1 = (0,1)

е2

1 х

Следовательно, устойчивое многообразие динамической системы (10) лежит на оси а, в то время как центральное многообразие, которое мы обозначим через их, выходит из начала координат полуплоскости р > 0 с коэффициентом угла наклона, равным —х/3. В точке Р2 = (0,10/3(х — 3)) имеем

(

Д (0,10/3(х — 3))

5

х — 3 10(2 + х)

V 3(х — 3)

\

10

х

—3

поэтому она является вырожденной особой точкой линеаризованной системы с соответствующими собственными значениями А1 = 5/(х — 3), Л2 = —10/(х — 3) и собственными векторами

е1 = 1,

2х + 4 9"

е2 = (0,1)

Неустойчивое многообразие лежит на оси а, и устойчивое многообразие, которое обо-

2х + 4

значим через £х, имеет коэффициент угла наклона с осью р, равный-в точке Р2.

Отметим, что многообразия (орбиты) их и Бх зависят от параметра х непрерывно [9].

3

0

В точке P3 = (1, -1) дифференциал равен

( 3

3 \

А(1,-1)= 4 10x)

\3+ 3(x - 3)

+

2

10

Собственные значения

являются комплексными для x Е (1.02,2.29) и вещественными для x Е (0,1.02) и x Е (2.29, 3), соответственно ReAi,2 > 0 при x < 1.95 и ReAi,2 < 0 при x > 1.95. Бифуркации Хопфа существуют при x = 1.95. Таким образом, точка (1, -1) является неустойчивым узлом при x Е (0,1.02), и неустойчивым фокусом при x Е (1.02,1.95) и соответственно устойчивым узлом при x Е (2.28, 3) и устойчивым фокусом при x Е (1.95, 2.28). Чтобы показать разрешимость рассматриваемой нелинейной задачи на собственные значения, необходимо доказать существование значения x* параметра x такого, что многообразия Ux* и Sx* совпадают, т. е. установить существование гетероклинической траектории, соединяющей точки P1 = (0, 0) и P2 = (0,10/3(x — 3)) на фазовой плоскости (р, а). Сначала приведём несколько лемм о свойствах траекторий (орбит) динамической системы (10).

Лемма 1. Существует решение системы (10) в виде а(т) = — (5/9)р(т) = —(x/3)x р(т) с параметром x = 5/3.

Доказательство этой леммы основано на прямых вычислениях. Заметим, что чисто колмогоровский спектр E = const • k-5/3 соответствует x = 5/3 (см. [1]) и соответствующее неустойчивое многообразие (орбита) U5/3 "уходит" на бесконечность при т ^ то.

Лемма 2. Не существует периодических траекторий системы (10), расположенных целиком в областях р + а > 0 и р + а < 0.

Из первого уравнения системы (10) легко показать, что выражение ^р/^т меняет знак, когда траектория пересекает линию р = —а, что и доказывает утверждение леммы.

Существование точного решения системы (10) позволяет установить нижнюю оценку значения x* .

Лемма 3. Неустойчивое многообразие (орбита) Ux не пересекает прямую а(т) =

— (5/9)р(т) для x < 5/3.

Предположим, что существует точка (ро, а0), в которой Ux пересекает прямую а(т) =

— (5/9)р(т). Можно считать, что многообразие (орбита) Ux задано как функция ax = а(р, x) в окрестности (р0, а0). Для достаточно малых значений р и а, т. е. вблизи начала координат (0, 0), имеем Ux > а(т) для x < 5/3 и Ux > а(т) на интервале (0,т0), где т0

определяет (р0,а0). Таким образом, в точке (р0,а0) производная меньше или

равна ^а/^р. С другой стороны, из системы (10) имеем

dax/dp = | 5р2 + 6pa — 9а2 +--+ xp) J j (2р(р + a)).

_2 , , „„м / л,^ , (12)

Рассмотрим разность da/dp и dax/dp в точке (po, ao), используя формулу (12). В результате получим

da/dp — dax/dp = x — g(3ao + xpo) ) / (2po(po + ao))

или

(30po / ao x — X—3 l P" + 3 II / (2Po(Po + ao)) . (13)

Здесь ao/po = —5/9, 2po(po + ao) > 0. Так как x < 5/3, то получим, что правая часть в (13) является отрицательной, что ведёт к противоречию. Таким образом, орбита Ux не может пересечь линию a(r) = — (5/9)р(т).

Следствие 1. Параметр x* > 5/3.

Лемма 4. Неустойчивое многообразие Ux лежит "ниже" линии a(r) = —(5/9)р(т), если x > 5/3.

Доказательство этой леммы аналогично доказательству леммы 3. Рассмотрим устойчивое многообразие Sx в области p + a < 0. Докажем свойство монотонности Sx при изменении параметра x.

Лемма 5. Если xi < x2, то Sxi лежит "выше" Sx2 в области p + a < 0.

Очевидно, что Sxi лежит "выше" Sx2 при достаточно малых p и xi < x2. Предположим, как и при доказательстве леммы 3, что существует точка (po,ao), в которой Sxi пересекает Sx2. Пусть aj = a(p, xj) определяет Sxi, i = 1, 2, локально и ai = a2 при p = po. Тогда dai/dp < da2/dp в точке пересечения Sxiи Sx2. Используя, как и в лемме 3, уравнения системы (10), получим следующую формулу для разницы dai/dp — da2/dp в точке po:

dai/dp — da2/dp = 10pЛ —op° + i--op° + ) / (2po(po + ao)) , (14)

\ xi — 3 x2 — 3 //

где a(po,xj) = ao в точке пересечения aj. Анализируя правую часть формулы (14) и принимая во внимание условие po + ao < 0, приходим к неравенству dai/dp > da2/dp при p = po. Последнее означает, что Sxi не может пересечь Sx2. Такое же утверждение справедливо и для неустойчивого многообразия Ux в области p + a > 0.

Лемма 6. Если xi < x2, то Uxi лежит "выше" Ux2 в области p + a > 0.

Доказательство этого утверждения проводится аналогично доказательству леммы 5. Для завершения доказательства монотонности £х и их вплоть до прямой р + а = 0 предположим сначала, что $Х1 и $х2 пересекают прямую р + а = 0 в точке (р0,а0), ро = —а0. В окрестности (р0, а0) можно представить БХ1 в виде р^ = р(а, ж^). Тогда из леммы 5 следует, что р1 < р2 вблизи (р0,а0). С другой стороны, в окрестности (р0,а0), р0 = —а0 имеем следующее разложение в ряд Тейлора:

/ ч 9оо 2

рг(а) = р0 +20(а0 + а0)(р + а) +

9( 9а02 — 24а0 + 30/(ж — 3)\ _ , Л

+ 2а 100(а0 + а0)2 ^ +а)3 + ...■ (15)

из которого следует, что р1 (а) > р2(а) вблизи (р0,а0) для р + а < 0. Полученное противоречие р1 < р2 доказывает, что БХн не пересекает линию р + а = 0. Аналогичные вычисления верны для траектории [/х..

При доказательстве свойства монотонности для Бх и Цх, неявно предполагалось, что £х и Цх пересекают прямую р + а = 0. Сначала установим, что для каждого фиксированного ж орбита £х лежит в ограниченной области квадранта р > 0, а < 0. Обозначим производную ^а/^р через С(р, а) и рассмотрим прямоугольный треугольник Г, ограниченный прямой р = (3/2)а — 5/(ж — 3) (которую обозначим через /) и осями р и а. Заметим, что коэффициент угла наклона прямой I мажорирует коэффициент угла наклона орбиты £х с осью а для значений ж > 5/3. Докажем, что векторное поле, определяемое С(р,а), направлено внутрь треугольника Г. Другими словами, докажем следующее утверждение.

Лемма 7. Для а Е (0,10/3(ж — 3)) имеет место неравенство

С(/,а) < 2/3. (16)

Рассмотрим векторное поле С(р,а) и подставим р = /. Тогда (16) переписывается в виде

45/4а2 — 15Аа + Аа2 + 2Аж/

—-< 1,

45/4а2 — 12Аа + 3А2 < '

что эквивалентно неравенству -2А2 + 3Аа — 2Аж/ > 0, где А = 5/(ж — 3). Неравенство имеет место, если А(ж — 1)(2А — 3а) > 0. Так как А< 0 и ж> 5/3, то для доказательства достаточно показать, что 2А — 3а < 0, или 10/3(ж — 3) < а. Последнее неравенство с очевидностью верно, что и завершает доказательство леммы.

Следствие 2. Орбита Бх всегда пересекают прямую р + а = 0 при ж > 5/3.

Из леммы 7 следует, что орбита £х, если её рассматривать как неустойчивое многообразие, выходящее из точки равновесия Р2, не может уйти в области р + а < 0 на бесконечность с изменением времени и не имеет предельного цикла (см. лемму 2). Если предположить, что орбита Бх не пересекает прямую р + а = 0, то она должна достигать точки равновесия Р1, приближаясь к ней ниже прямой р + а = 0, что ведёт к неравенству ж/3 > 1 и противоречит условию ж < 3. Следовательно, орбита Бх должна пересекать прямую р + а = 0.

Рассмотрим орбиты для значений x > 5/3.

Лемма 8. Существует x** такое, что орбита пересекает прямую р + а = 0 для всех x > x**

Предположим противное, т.е. орбита при всех значениях x > 5/3 не пересекает прямую р+а = 0. Тогда орбита полностью расположена в угле со сторонами р+а = 0 и а = (—x/3)p. В силу монотонной и непрерывной зависимости от параметра x, координаты р(т) и а(т) орбиты при каждом фиксированном т сходятся к конечным значениям при x ^ 3. Интегрируя второе уравнение системы (8), получим

а(т) - а(то) = I ( 5р2 - 6ра - 9а2 + x—3 (3а + xp) ) . (17)

Перейдём к пределу в (17), устремляя x ^ 3. Предел в левой части существует и конечен. Для существования предела в правой части необходимо положить р(т) = —а(т) при x ^ 3. Таким образом, последовательность Ux(t) сходится к р(т) + а(т) = 0 при x ^ 3 для всех т G (—го, го). Рассмотрим формулы (8) для F(п) и dF(n)/dn. Так как р(т) = —а(т), то получим, что f (п) = — g(n), П ^ (0,1). Дифференцируя первое соотношение из формул (8) и сравнивая полученное выражение со вторым соотношением, приходим к равенствам f (п) = const, g(n) = —const. Следовательно, р(т) = const, а(т) = —const для т G (—го, го). Последнее противоречит тому, что р(т) + а(т) = 0 определяет прямую на фазовой плоскости (р, а).

Полученные результаты позволяют доказать следующее утверждение.

Теорема. Существует гетероклиническая траектория, соединяющая точки равновесия Pi и P2.

Рассмотрим последовательность точек sm, m ^ го, на фазовой плоскости (р,а), которая образуется в результате пересечения орбит SXm с прямой р + а = 0, где xm — возрастающая последовательность значений параметра x. Обозначим через (рт, ат) координату точки sm, где рт = —ат. В силу монотонной и непрерывной зависимости от параметра x значения рт монотонно и непрерывно возрастают при m ^ го, соответственно ат монотонно и непрерывно убывают при изменении m. Другими словами, точки sm движутся вниз по биссектрисе угла квадранта (р, а), р > 0, а = —р. Принимая во внимание нижнюю оценку значения искомого параметра x* > 5/3, выберем начальное значение параметра x0 последовательности xm равным x0 = 5/3. Последовательность точек sm является неограниченным множеством на прямой р + а = 0 при m ^ го или x ^ 3. Предположим противное, тогда существует ат = — рт = s < го.

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

7 2 2 10 ,

а(т) — ат = J ( 5р — 6ра — 9а +--- (3а + xp)

xm 3

где SXm (тт) определяет точку sm на прямой р + а = 0. Выбирая m достаточно большим (или xm близким к 3), получим, что правая часть интегрального соотношения является сколь угодно большой величиной. В силу конечности значений а(т) и ат приходим

х = 1.8! >091

/

л-1-1-1-1-1-1-

"О 0.2 0.4 0.6 0.8 1 1.2 1.4

Р

Рис. 1. Гетероклиничекая траектория, соединяющая точки равновесия Р\ и Рз при х* = 1.85091

к противоречию с предположением существования конечного предела последовательности 5т. Что касается точек пересечения орбиты их с прямой р + а = 0 при увеличении параметра х, то ввиду установленных свойств эти точки движутся вдоль прямой по направлению к Р2. Следовательно, существует единственная точка 5* (и соответственно значение параметра х*) последовательности эт, в которой £х и их пересекаются на прямой р + а = 0. Таким образом, £х* и их * порождает гетероклиническуюю траекториюю, соединяющую точки Р\ и Р2 (рис. 1).

Л. Визуализация орбит динамической системы

Рисунки 2, а-г иллюстрируют поведение траекторий динамической системы (10) при различных значениях х (штриховые линии соответствуют гетероклинической траектории).

Приведённый анализ фазового портрета системы (10) показывает существование решения изучаемой краевой задачи на собственные значения, которое описывается гетеро-клинической траекторией динамической системы (10). Численное значение параметра х для гетероклинической траектории х* = 1.85091.

В. Граничные условия и область значений для т

Граничные условия для функции Р(п) имеют вид

Ь2

Р ~ Р0п х при п ^ 0, Р (1 - п)2 при п ^ 1, (18)

где Р0 — произвольная константа. После перехода к функции f и независимой переменной 5 = 1п п граничные условия приобретают вид

1 3-х 5Ь 3

f ~ 5Р02е5 при 5 ^—то, f ~—— е23(1 - е3) при 5 ^ 0. (19)

4 0 6 0 8 1 12 1.4

Р

в

1 1.2 14

Рис. 2. Траектории, выходящие из точки Р1 (а, б) и точки Р2 (в, г) для параметров в диапазонах х € [х*, 1.86] (а), х € [1.84, ж*] (б), х € [х*, 1.86] (в), х € [1.84, ж*] (г)

а

г

Переход к независимой переменной т определяется заменой

Г

/ (я')-

80

Константу в0 необходимо определять из граничных условий (19) на функцию /(я). Поскольку она имеет нули в начальной (в = —го) и фронтальной (в = 0) точках, то интеграл имеет особенности и необходимо выбирать точку в0 € (—го, 0) так, что /(в0) > 0. Интегрируя в граничных точках, получим

[ ^в' [ ^в 12 (х-з)„

-в 2 ^ —го при в ^ —го,

I /(в') 3 5^в3-"8 5Р022 х — 3

Г ^в' Г ^в 2

— (— 1п(—в)) ^ +го при в ^ —0.

J /и у 5ь зв(1 8) 5Ь

80 — — в 2 8(1 — в8)

Таким образом, граничные условия (18) задают область значений для независимой переменной т.

Список литературы

[1] CONNAUGHTON D., Nazarenko S. Warm cascade and anomalous scaling in a diffusion model of turbulence // Phys. Rev. Lett. 2004. Vol. 92(4).

[2] Leith С. Diffusion approximation to inertial energy transfer in isotropic turbulence // Phys. Fluids. 1967. Vol. 10. P. 1409-1416.

[3] Besnard D., Harlow F., Rauenzahn R., and Zemaoh C. // Theor. Comput. Fluid Dyn. 1996. Vol. 8(1).

[4] Станюкевич К.П. Неустановившееся движение сплошной среды. Изд. 2-е. М.: Наука, 1971.

[5] von Guderley G. Starke kugelige und zylindrische Verdichtungsstosse in der Nahe des Kugelmittelpunktes bzw. Der Zylinderachse // Luftfahrtforschung. 1942. Vol. 19.9. P. 302-312.

[6] Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. Изд. 2-е, доп. М.: Наука, 1966. 688 с.

[7] Grundy R.E. Large time solution of an inhomogeneous nonlinear diffusion equation // Proc. R. Soc. London. 2004. Vol. A386. P. 347-372.

[8] Connaughton D., Nazarenko S. A Model Equation for Turbulence. arXiv:physics/0304044. 2004.

[9] Marsden J.E., MoCraoken M. The Hopf Bifuraction and its Application. New York: Springer-Verlag, 1976.

Поступила в 'редакцию 11 сентября 2013 г., с доработки — 12 ноября 2013 г.