Численное моделирование динамики вихревого клубка в сверхтекучем гелии

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

Том 15, № 2, 2010

Численное моделирование динамики вихревого клубка

в сверхтекучем гелии*

Л.П. Кондаурова, В. А. Андрющенко, С. К. Немировский Институт теплофизики им,. С. С. Кутателадзе СО РАН, Новосибирск, Россия

e-mail: theory@itp.nsc.ru

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

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

Введение

Динамика вихревых структур во многом определяет гидродинамические и термодинамические свойства сверхтекучего гелия [1, 2]. Численное изучение совокупности вихревых нитей в сверхтекучем гелии на основе уравнения динамики вихревой нити (уравнение Био — Савара) при различных начальных и граничных условиях проводилось многими авторами [3-7]. Изначально простая вихревая структура в своей динамике претерпевает существенные изменения, превращаясь в сильно запутанный вихревой клубок. При своем движении вихревые нити могут пересечься, в результате чего происходит либо слияние двух петель в одну. либо дробление одной петли Hci две (процессы реконнекции). Эти процессы оказывают существенное влияние на эволюцию вихревого клубка, так как изменяют топологию вихревой структуры.

В вышеупомянутых численных исследованиях процессы реконнекции, следуя Шварцу [4], осуществлялись при достижении между отрезками вихревой нити/нитей некоторого расстояния А, которое находилось из условия равенства нелокального и локального вкладов в определении скорости точки, расположенной на вихревой нити.

*Работа выполнена при финансовой поддержке РФФИ (гранты № 10-08-00369 и 10-02-00514).

© ИВТ СО РАН, 2010.

Согласно оценке^ СДбЛсШНОИ Шварцем, это расстояние определяется выражением А = 2Я/[с 1п(Я/а0)], иначе говоря, если расстояние между элементами нитей (нити) становится меньше Д, то процесс реконнекций осуществляется. Здесь Я — радиус кривизны в рассматриваемой точке, с — постоянная (~ 1) а0 — радиус ядра вихря. Однако малое расстояние между элементами вихревой нити/нитей не гарантирует того, что при своем дальнейшем движении произойдет их пересечение. В настоящей статье предлагается новый критерий осуществления реконнекций, основанный на анализе динамики вихреВЫХ НИТ6И. Авторами проведено численное исследование динамики вихревого клубка в безграничном пространстве, которое моделировалось с помощью периодических граничных условий.

1. Математическая модель

1.1. Уравнение движения вихревой нити

Гидродинамика сверхтекучего гелия построена в терминах двухжидкостной модели [8], суть которой состоит в том, что в сверхтекучем гелии может существовать одновременно два движения: нормальное и сверхтекучее. Каждое из этих движений связано со своей эффективной массой. Сумма этих масс равна полной истинной массе жидкости. Обозначим скорости нормального и сверхтекучего движения соответственно как vn и vs. Плотность потока массы жидкости определяется в виде j = pnvn + psvs, где pn, ps — нормальная и сверхтекучая плотность жидкости. Их сумма равна истинной плотности гелия в сверхтекучем состоянии: р = pn + ps.

Уравнение движения вихревых нитей можно получить следующим образом. Рассмотрим уравнение неразрывности

ft +divj = 0-

При постоянной температуре T плотности нормальной и сверхтекучей компонент остаЮТСЯ ПОСТОЯННЫМИ^ Т. 6. Pn = consti, ps = const2. Тогда для несжимаемой жидкости (см. (1)) pnVvn + psVvs = 0, или Vvn = 0 Vvs = 0. Поскольку квантованные вихревые нити могут возникать только в сверхтекучей компоненте, то в дальнейшем нас интересует только соотношение

Vvs = 0.

Уравнение (2) означает, что поле соленоидальное и скорость можно представить в виде ротора некоторого векторного потенциала

vs = го t A.

Определим ротор скорости, т.е. завихренность: ш = rotvs = rot rot A (см. (3)). Используя соотношение из векторного анализа, получим: rot rot A = graddivA — ДА. Векторный потенциал определяется с точностью до градиента некоторой произвольной

A=0

Таким образом, для поля завихренности получим уравнение Пуассона

AA = - ш.

Решение этого уравнения имеет следующий вид:

=¡л / тТ—^<* <в>

где ¿V — элемент объема, б — радиус-вектор, направленный к некоторой точке вихря, г — радиус-вектор, направленный к некоторой точке пространства. Взяв ротор от выражения (5), получим уравнение для скорости сверхтекучей компоненты:

/ ч 1 С - х (г — б) 1тг ..

Для нахождения этого интеграла необходимо знать поле завихренности. Согласно Онза-геру [9] и Фейнману [10], циркуляция Г скорости сверхтекучей компоненты квантована:

Г = V. • ¿1 = пк.

с

Интегрирование идет по замкнутой (без самопересечений) кривой С в сверхтекучей жидкости. Величиной к обозначен квант циркуляции для поля скорости, равный к/ш4, где к — постоянная Планка, ш4 — масса атома гелия, п — целое число. Согласно тео-

С

открытую поверхность Б\

Г = j • ¿1 = Ц

с

где е — единичный вектор, перпендикулярный поверхности ¿Б. Согласно теореме Кельвина, циркуляция скорости по замкнутому контуру в идеальной жидкости при предположении непрерывности скорости и гладкости деформации контура является инвариантом движения: ¿Г/йЬ = 0.

В сверхтекучем гелии образуются вихри с радиусом ядра а0 порядка ангстрема. Иначе говоря, распределение завихренности отлично от нуля только вдоль некоторой пространственной кривой (линии). Такое распределение завихренности называют бесконечно тонкой вихревой нитью. Далее под термином "вихревая нить" подразумевается бесконечно тонкая вихревая нить. Рассмотрим замкнутую вихревую нить. Вектор завихренности будет параллелен элементу нити ¿е. Тогда для циркуляции скорости из уравнений (7) и (8) получим следующее выражение:

Г = -5Б = пк

(здесь и далее полагаем п =1, так как существование таких вихрей является энергетически выгодным). Вернемся к уравнению для скорости сверхтекучей компоненты. В уравнении (6) элемент объема можно представить в виде произведения длины элемента нити и площади ее поперечного сечения: ¿V = 6Бds. Далее, используя соотношение (9), получим уравнение, представляющее собой интеграл Био — Савара:

¿бв 1 у! ш х (г — s)¿Бds к Г (б — г) х ¿б

L L

О

Рис. 1. Схема векторов для фрагмента вихревой нити; О — начало координат, 1 — элемент вихревой линии, 2 — некоторая точка пространства

Это выражение определяет скорость жидкости в некоторой точке г (рис. 1). Для нахождения скорости вихревой нити величина г должна принимать значения всех элементов нити:

Здесь и далее б(£) — радиус-вектор, направленный к некоторой точке вихря, параметр £ — длина дуги, которая пробегает значения от 0 до Ь (полная длина вихревой нити), б' — производная радиус-вектора по параметру, или касательный вектор. Интеграл (11) расходится в точке = б(£,1). Для его нахождения вводится параметр обрезания 6 [11-13]. Величина 6 определяется из выражения

Лучшее согласие с экспериментом [14] наблюдается при предположении, что в сверхтекучем гелии вихрь имеет ядро круглого сечения, в котором завихренность ш постоянна, что соответствует модели вихря Рэнкина; тогда g = 1/4 Подставив это значение g в уравнение (12), получим 6 = 0.64 • a0, т.е. величина параметра обрезания имеет тот же порядок, что и значение радиуса ядра.

Полная скорость движения элементов вихревой нити vsi (без учета взаимодействия с нормальной компонентой жидкости) является суммой скорости движения сверхтекучей компоненты жидкости vs0 и скорости vB, индуцированной всеми присутствующими в системе вихревыми нитями. Влияние границ объема приводит к возникновению некоторой коррекции полной скорости на величину Vs,b:

Для описания динамики вихревого клубка необходимо также учитывать силу трения возникающую при взаимодействии вихрей с нормальной компонентой жидкости. В общем случае скорость элемента вихревой нити (у^) отлична от скорости сверхтекучей жидкости в месте расположения вихревой нити (у8г), что приводит к появлению силы Магнуса. Учтем влияние сил трения и Магнуса. Величина силы трения зависит от разности скоростей нормальной компоненты жидкости уп и вихря у^. Согласно [2],

ln(26/ao) = g.

Vsl = Vß + Vso + Vs,b.

ftr = -Tos' x [s' X (Vn - VL)] + yOs' X (Vn - VL),

где коэффициенты трения 70 и 70 зависят от температуры и находятся экспериментально. Сила Магнуса определяется следующим выражением [11]:

41 = РзК&' х

Поскольку масса ядра вихря очень мала, то ее можно принять равной нулю. Тогда сумма сил, действующих на вихрь, должна быть равна нулю, т. е. + ^ = 0, или, подставляя их значения из (14), (15), получим:

б' х [р»£

Это уравнение выполняется в случае, если вектор, равный сумме векторов в квадратных скобках, параллелен б' или равен нулю. Движение вдоль вихря (вдоль б ^ не влияет на значение скорости вихря. Отсюда следует, что выра^кение в квадратных скобках уравнения (16) равно нулю. Из последнего можно найти скорость вихря:

Уь = V» + а

где а и а' выражаются через 7 и следующим образом:

1окр8 ' То + т0 (т0 - ^ а =- а =-

То + (Т0 - ' 70 + (Т0 - кр»)2 '

Численные исследования проводились при Т = 1.6 К. При этой температуре, согласно [12], а = 0.098, а' = 0.012. Для нахождения скорости движения элементов вихревых нитей (см. уравнения (17), (13), (11)) необходимо определить значение индуцированной скорости ув. В проведенных расчетах использовалось локальное приближение, согласно которому движение точки нити определяется только ближайшими к ней отрезками (участками вихревой нити). Рассмотрим уравнение (11). Разложим в(£') в ряд Тейлора в точке в(£): б(£') - б(£) = £б' + £2б''/2 + ..., и продифференцируем это выражение по £: ¿Б(£')/ё£ = б' + £ б'' + ... Подставив полученное выражение в (11), получим

У =

В 4п/ |£б' + £2б''/2 + ...|3 1 '

ь

Пренебрегая в (18) производными высших порядков, вычислим скорость в локальном (самоиндуцированном) приближении:

/ ч, „"Л.

v», = — б x б

8п

ь

Уравнение (19) можно записать в виде

V», = в б' х б''.

В работе [7] показано, что лучшая аппроксимация формулы (11) выражением (20) достигается при в = с (к/4п)\п(К/а0) и константе с = 1.1, что несколько отличается от выражения для в = (к/4п) • 1п(сЯ/а0), где с имеет величину порядка 1. Послед-

в

[7]. Расчеты проводились для периодических во всех направлениях граничных условий: У»(х) = У»(х + I), У»(у) = У»(у + I), У3(г) = У3(г + /) (период I = Хшах - Хшш = Ушах Ушт — ^тах ^Шт)  Если вихревои сегмент нити выходит за границы расчетной области с одной стороны, то входит в нее с противоположной стороны и при этом сохраняет направление б'.

1.2. Процессы реконнекции

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

т

Из условия совместности системы уравнений

+ У,гК + (¿г+1 — ¿г)в1 = + Угл К + (¿л+1 — )в2, 0 < 31 < 1, 0 < $2 < 1

определяем, возможно ли пересечение отрезков нитей. Здесь (хг,уг, гг,хг+1,уг+1, гг+1) — координаты первой и (xj,yj, Zj,Xj+1,yj+1, Zj+1) второй тар точек, а Ух,г, Уу, г, Уг,г, Ухл, Уу л, Угл — проекции скоростей точек на координатные оси. Решение этой системы имеет следующий вид:

К = — (—(хг+1 — хг)(^+1 — ^ )Уг + (хг+1 — хг)(^+1 — ^ )УЛ + (хг+1 — хг)(УЛ+1 — Уj )хг + (^+1 —

Ъ )(Уг+1 — Уг)хг — (^+1 — Ъ )(Уг+1 — Уi)xj — (хг+1 — ¿г)(хЛ+1 — xj )УЛ — (Уг+1 — Уi)(xj+1 — xj )хг) — (хг+1 — хг)(Ул+1 — Уj)^ + (Уг+1 — Уi)(xj+l — xj)^ + (Ул+1 — Уj )(хг+1 — zi)xj — (Ул+1 — Уj )(хг+1 — хг)хг + (¿г+1 — — Xj )Уг)/(—(хг+1 — хг)^+1 — Уj )УХ л + (хг+1 — хг)^+1 —

Уj)у ,г + (хг+1 — хг)(^+1 — ^ )Уул — (хг+1 — zi)(xj+l — xj )Уу л + (Уг+1 — Уi)(xj+l — xj)у л + (Уj+l — Уj )(хг+1 — ¿г)ух л — (Уj+l — Уj )(хг+1 — ¿г)ух, г — (Уг+1 — Уi)(xj+l — xj)у ,г — (хг+1 — хг)(^+1 — Ъ )уу, г + (^+1 — Ъ )(Уг+1 —Уг)ух, г — (^+1 — Ъ )(Уг+1 — Уг)ух л + (хг+1 — Хг )(хл+1 — хл )Уу, г)-

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

хг + Ух, гК + (хг+1 — хг)в1 = хл + Ух л К + (хл+1 — хл )$2, Уг + Уу, гК + (Уг+1 — Уг)в1 = Ул + Уу л К + (ул+1 — Ул )в2,

1'

Рис. 2. Пересекающиеся участки вихревых нитей

2. Вычислительная схема

2.1. Численная схема решения уравнения динамики квантованных вихревых нитей

Для нахождения скорости V», (см. уравнение (20)) необходимо определить производные б', б''. Они могут быть вычислены с помощью вспомогательных векторов 1+ = б,+1 - б, и 1_ = б, - б^. Разложим век торы б^+1, б_ в ряд Тейлора в точке б,; тогда выражения для 1+ 1_ имеют вид

, 8+Б» 8+8™ 8+3'"'

1+ = 8+Б' + -^ + -1Т + -+Г + -'

, в! б'' 83 б''' 8_ Б'"'

1_ = 8_Б' - ~_Г + — - -_2Г + -

Здесь 8+ — длина дуги между б^+1 и б,, а в_ — между б, и б,_1. Из этих уравнений можно построить аппроксимацию производных и получить выражение для их векторного произведения:

Б'х Б'' = +1_ |Г+2+1_|'+х 1_ +е

которое определялось с точностью

е = б' х б''' + ...

Для нахождения скорости точек вихревых петель необходимо подставить вычисленное значение V», в уравнение (17).

2.2. Регулирование количества точек на вихрях

Точность аппроксимации уравнения (22) определяется выражением (23). Из (23) следует, что при уменьшении разности 8+ и в_ точность аппроксимации увеличивается. С целью уменьшения е в вычислительной схеме используется следующая процедура: если расстояние между точками становится в два раза меньше изначально выбранного

пространственного шага вдоль вихревои нити Д£0 (шаг выбирался из условия согласования численного и аналитического значений скорости движения вихревого кольца), то промежуточная точка убирается, и наоборот, если расстояние между точками б, и б,+1 становится в два раза больше Д£0, то между ними вставляется еще одна точка с радиус-вектором б^. Точка г' вставляется на середине дуги окружности радиусом Я,' = |б"| = |(б,' + б,'+1)/2|. Координата ее определяется выражением б= (б, + б,+1)/2 + Я2 - /+/4 - Я)Я'б". Петли из трех точек в процессе расчета убирались. Уменьшение длины нитей за счет удаления малых петель составляло величину примерно на три порядка меньше по сравнению с уменьшением ее за счет силы трения.

Поэтому можно считать, что эта процедура не является источником ошибок.

2.3. Отбор отрезков — кандидатов для реконнекции

При отборе пар соседних точек, являющихся концами отрезков, которые могут пересечься в течение временного шага, применялась следующая процедура. Все пространство, занимаемое вихревым клубком, разбивалось на ячейки, размеры которых определялись как Лх = Кшахт, Ау = Уушахт, Лг = Ушахт. Здесь У?™, Ушах — максимальные проекции скоростей точек вихревых нитей на соответствующие оси координат.

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

3. Результаты расчетов

Вычисления проводились в кубическом объеме, на границах которого задавались периодические граничные условия по всем координатам. Расчеты были проведены для различных скоростей противотока Упз = Уп — а имени о: Упз = 6, 8 и 12 см/с. Начальная вихревая конфигурация представляет собой двадцать четыре вихревых кольца (рис. 3), расположенных в кубе со стороной I = 50 мкм таким образом, чтобы полный начальный импульс системы был равен нулю. Для интегрирования уравнения движения по времени была использована схема Рунге Куттьт четвертого порядка точности. Максимальный радиус кольца Я был равен 5 • 10"4 см? пространственный шаг вдоль линии Д£о = 2пЯ/100. Шаг по времени выбирался из условия устойчивости схемы: т = 2\/2а2/(4в) (где а — пространственный шаг) для а = Д^0/2.

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

а б

Рис. 3. Развитие вихревого клубка при скорости противотока Упз = 12 см/с; I = 0

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

Ь

ционарному состоянию, от скорости противотока Упз, было получено выражение

Ь =

с 7 ~ 280 (рис. 5). Такая же зависимость Ь от Упз наблюдалась в эксперименте [15]. Заметим, что в работе [4] равновесное состояние достигалось при использовании специальной процедуры: поворота вихрей на 90° относительно направления скорости противотока.

о 0.01 0.02 0.03

с

Рис. 4. Плотность вихревых нитей Ь как функция от времени при скорости противотока Упз, равной 6

Рис. о. Плотность вихревых нитей как функция от скорости противотока в логарифмических координатах; точки — полученные численные результаты

Рис. 6. Зависимость количества перезамыканий N с образованием и делением петель от времени; верхние линии соответствуют процессу деления петель, нижние — процессу образования

I, см

Рис. 7. Зависимость скорости рекоттпекций от плотности вихревого клубка в логарифмических координатах; точки — полученные численные результаты

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

Из полученных результатов можно определить скорость реконнекций N' = dN (число перезамыканий в единицу времени) для каждого стационарного состояния

N' « СкЬ5/2

с С ~ 2.47 (рис. 7). Эта зависимость ранее была предсказана в расчетных [3, 16] и аналитической [17] работах.

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

вихревых нитей L от скорости противотока Vns и скорости реконнекций N' от L. Полученные результаты хорошо согласуются с данными экспериментальных, расчетных и аналитических работ.

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

[1] Nemirovskii S.K., Fiszdon W. Chaotic quantized vortices and hydrodynamic processess superfluid helium // Rev. Mod. Phys. 1995. Vol. 67, No. 1. P. 37-84.

Donnelly R.J. Quantized Vortices in Helium II. Cambridge Univ. Press, 1991.

Tsubota Makoto, Araki Tsunehiko, Nemirovskii S.K. Dynamics of vortex tangle without mutual friction in superfluid 4He // Phys Rev. B. 2000. Vol. 62, No. 17. P. 11751-11762.

4

boundary interactions // Ibid. 1985. Vol. 31, No. 9. P. 5782-5804.

Buttke T.F. A numerical study of superfluid turbulence in the self-induction approximation // J. Сотр. Phys. 1988. Vol. 76. P. 301-326.

Samuels D.C. Velocity matching and poiseuille pipe flow of superfluid helium // Phys. Rev. B. 1992. Vol. 46, No. 18. P. 11714-11724.

Aarts R.G.M., de Waele A.T.A.M. Numerical investigation of the flow properties of Hell // Ibid. 1994. Vol. 50, No. 14. P. 10069-10079.

Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 736 с.

Onsager L. Unpublished Remark at Low Temperature Physics Conference at Shelter. Island, 1948.

Feynman R.P. Application of quantum mechanics to liquid helium / Ed. C.J. Gorter / / Progress in Low Temherature Physics. North-Holland, 1955. Vol. 1. P. 17-53.

Moore D.W., Saffman P.G. The motion of a vortex filament with axial flow // Philos Trans. R. Soc. London. A. 1972. Vol. 272. P. 403-429.

Saffman P.G. The velocity of viscous vortex rings // Stud. Appl. Math. 1970. Vol. 49. P. 371380.

Saffman P.G., Baker G.R. Vortex interactions // Annual Rev. Fluids Mech. 1979. Vol. 11. P. 95-122.

Rayfield G.W., Reif F. Quantized vortex rings in superfluid helium // Phys. Rev. 1964. Vol. 136. P. A 119! A 1208.

Hall H.E., Vinen W.F. The rotation of liuid helium II: I. Experiments on the propagation of second sound in uniformly rotation helium II // Proc. R. Soc. London. A. 1956. Vol. 238. P. 204-214.

Barenghi C.F., Samuels D.C. Scaling laws of vortex reconnections //J. Low Temper. Phys. 2004. Vol. 136. P. 281-293.

Nemirovskii S.K. Evolution of a network of vortex loops in He-II: Exact solution of the rate equation // Phys. Rev. Lett. 2006. Vol. 96. P. 015301-015304.

Поступила в редакцию 29 января 2009 г., с доработки — 20 ноября 2009 г.