рядок и повысить разрешение. Так, для задач планетарного динамо не исключена возможность, что на малых масштабах l < в основном объеме жидкого ядра энергия полей сравнима с энергией на больших масштабах [Kutzner, Cristensen, 2002; Simitev, 2004] (см. также обсуждение в [Решетняк, 2006]). Поскольку мелкомасштабные поля невидимы для наблюдателя на поверхности Земли и их характеристики слабо зависят как от геометрии области, так и от граничных условий [Sarson, Jones, 1997; Reshetnyak, Steffen, 2005], анализ задачи в ящике может быть весьма полезным. Ниже мы рассмотрим как используя достаточно простую модель во вращающемся прямоугольном ящике, имитирующем поведение магнитогидро-динамической турбулентности вдали от твердых границ, удается воспроизвести ряд интересных свойств системы планетарного динамо.
2. УРАВНЕНИЯ ДИНАМО
Рассмотрим уравнения динамо для несжимаемой жидкости (V ■ V = 0) во вращающемся с угловой скоростью О относительно вертикальной оси z бесконечном слое 0 < z < 1. Введя следующие единицы измерения для скорости V, времени t, давления P и магнитного поля B: к/L, L2/k, pK2/L2 и 2Оркц0, где L - единица длины, к - коэффициент молекулярной теплопроводности, р - плотность вещества, - магнитная постоянная, запишем систему уравнений динамо в Декартовой системе координат (x, y, z) в виде (см. для сферической геометрии [Jones, 2000a]):
d B
= rot( V x B)
d t
1 AB,
EPr
-i
dV
dV_ Vx(V x V)
(1)
= rotB x B - V P - 1z x V + Rar 1z + EAV, дТ + (V -V)(T + To) = AT. Безразмерные числа Прандтля, Экмана, Рэлея и
Робертса заданы a go 8 TL
виде: Pr = -, E =
к'
2 ü L
;, Ra =
к
2 йк
и q = n, где v - коэффициент кинема-
тической вязкости, а - коэффициент объемного расширения, £0 - ускорение свободного падения, 5Т - единица возмущения температуры Т относительно "диффузионного" распределения температуры Т0 = 1 - г, п - коэффициент магнитной диффузии. Число Россби введем как Ro = ЕРг-1.
Поскольку далее нас будут интересовать мелкомасштабные решения, то в горизонтальной
плоскости (х, у) мы сводим задачу в бесконечном слое к задаче в ящике с периодическими относительно координат х и у граничными условиями [Сайапео, 2003]. Для границ г = 0, 1 принимаем нулевые возмущения температуры Т = 0, что с учетом выбранного профиля Т0, эквивалентно заданию температур на границах: ^ = Т + Т0 = 1, 0. Для поля скорости принимаем условие непроникновения и равенство нулю градиентов тангенциаль-
ЭУх ЭУу
ных компонент на г = 0, 1: Vг = ---— = = 0. Для
д г дг
магнитного поля потребуем непрерывность нормальной компоненты поля Вг на границе и равенство нулю тангенциальных компонент поля: Вх = д Вг
= Ву = = 0. Такая постановка граничных устоев
вий гарантирует равенство нулю тангенциальных компонент тензора вязких напряжений и тороидальной компоненты магнитного поля на границах.
3. ЧИСЛЕННЫЕ МЕТОДЫ 3.1. Псевдо-спектральный метод
Идея псевдо-спектрального метода состоит в том [Orszag, 1971], чтобы проводить численное интегрирование уравнений в частных производных вида (1) в волновом пространстве. В этом случае операции дифференцирования-интегрирования сводятся к умножению-делению на соответствующее волновое число, что дает машинную точность вычислений. Для вычисления же нелинейных членов, делается переход в физическое пространство, поля перемножаются и далее осуществляется обратный переход в волновое пространство. При таком подходе наибольшее число операций приходится на переход из одного пространства в другое, и такие переходы должны быть "быстрыми". Для рассматриваемой декартовой системы координат мы используем быстрое преобразование Фурье, адаптированное для параллельных вычислений.
В волновом пространстве система уравнений (1) может быть записана в виде [Buffett, 2003]:
d B -Ь2„' — + q к B dt
= [V x (V x B)]
k>
n -i dV ,2,/ Pr -т— + к V
dt
= k ®к + Fk,
+ к2 T
t
= -[(V -V)(T + To)]k,
k
k
k
где
Перепишем его в виде
^ —^--г, ^ = kpkp, в = 1,..., 3,
^ (3)
Ек = [ БРг-1 V х (Ух V) + ЯаТ - 12 х V + (В • У)В ]к,
а T, V и В - трехмерные образы Фурье исходных физических полей, к - волновой вектор. Пусть в физическом пространстве поле f задано на сетке G = (1, ..., N.., 1, ..., N, 1, ..., Nz). Тогда физическое и волновое представление в ящике высотой 1 и квадратным основанием длиной X связаны соотношением
dAe
kyt
д t
= Ue
kyt
(7)
и далее используем метод АВ2 для новых переменных А = Аек<*, и = иек <*.
Ниже нам будут полезны некоторые величины: разложение векторного соленоидального поля на полоидальную и тороидальную компоненты:
V = VP + VT, B = BP + BT,
(8)
где
f( x, У, z) =
(4)
XXX f(nx, ny, nx)e'{k'X+VV(kzZ),
где
K = 2-nnx, K = 2jn, K = ^, (5)
X
X '
dA + k2A = U. д t
(6)
Vp = VxVx( 1,f), Vt = Vx( I, e), Bp = VxVx( I, h), Bt = V x (1zg),
(9)
ае, к, g - скалярные потенциалы, зависящие от радиус-вектора г и времени I. Тогда в волновом пространстве имеем:
e=
i (k x V )z
i( k x B )z
|k|2
|k| ^
V,
B,
(10)
и ф = sin или cos (или их комбинации) в зависимости от вида граничных условий. Далее для удобства мы опускаем символ .
Пояснения в (2-5) требует уравнение Навье-Стокса. Во-первых, для снижения численных не-устойчивостей, связанных с ошибками типа "aliasing" мы записали его в (1) в вихревой, консервативной форме [Canuto et al., 1988]. Далее, исключили давление, умножив уравнение движения скалярно на k и использовав условие несжимаемости k • V = 0. Отметим, что а) все поля действительны, и можно использовать более экономные модификации преобразований Фурье для половины периода; б) поля обладают дополнительной симметрией по направлению z и могут быть разложены либо по синусам, либо по косинусам, что также сокращает объем вычислений (см. подробней в [Press, Teukolsky, Flannery, 2002]). В качестве дополнительных граничных условий в ^-пространстве мы полагаем T, V, B = 0 для k = 0.
Для интегрирования по времени использована явная схема 2-го порядка Адамса-Башфора (AB2) для всех членов кроме диффузионного. Для вычисления же диссипативных членов использован известный "аналитический" прием [Canuto et al., 1988]. Рассмотрим уравнение
f = ik2' " ik2'
В силу ортогональности полоидальных и тороидальных компонент
V2 р т В2 р т
Ек = ~2 = Ек + Ек, Ем = 2Я- = Ем + Ем, (11)
удобно ввести энергии покомпонентно:
р Vр р Вр т Ут т Вт Ек = т, Ем = т¡7::, Ек = т, Ем = тгт:- (12)
2 '
-м'
2Ro'
2Ro
Число Пекле, Рейнольдса, и магнитное число Рей-нольдса зададим в виде:
Pe = J2EK, Re = PePr1, Rm = Pe q. (13)
Введем гидродинамическую спиральность ЖЖ (см. подробнее в [Решетняк, 2006]) токовую магнитную Жл и кросс-спиральности Ж4 [Brandenburg, Subramanian, 2005]:
0Ж
= < V • rot V>, Ж* = Ro-1 < B • rotB>,
= Ro-1 < A • B>, Ж% = Ro-12 < V • rotB>,
(14)
где А - вектор-потенциал (В = У х А). С учетом используемых граничных условий для V и В все спиральности (14) равны нулю на границах г = 0, 1.
nx = -NX/2 n = -N /2 nz = -Nz/2
Также введем трехмерный F(k) = |/(k)|2 и энергетический спектры E(k) = k2F(k)1.
3.2. Параллелизация и тесты
Следующим важным моментом является адаптация кода для многопроцессорной техники. Нами использован наиболее гибкий в настоящее время метод MPI, позволяющий проводить вычисления на компьютерах с различной архитектурой: как на кластерах, так и на многопроцессорных суперкомпьютерах с общей шиной (в частности, суперкомпьютерах IBM Regatta p69+). Суть параллелизации состоит в том, что волновое и физическое пространства разбиваются между процессорами в одном, или более, направлениях, и далее вычисления производятся на каждом процессоре независимо с периодическим обменом данными между процессорами. Поскольку мы ориентируемся на машины и с небольшим числом процессоров (несколько десятков), то выбрана одномерная разбивка в z-направлении для k-пространства и y-разбивка для конфигурационного пространства (см. подробнее о двумерной разбивке, ориентированной для Крэя в [Cattaneo, Emonet, Weis, 2003]). Нетривиальным является адаптация трехмерного быстрого преобразования, являющегося последовательностью трех одномерных быстрых преобразований Фурье (БПФ). Существуют две стратегии: сделать два одномерных БПФ (без обмена данными между процессорами) и одно многопроцессорное БПФ (с обменом данными), либо применить два одномерных (т.е. сделать двумерное преобразование по переменным x, y) без обмена, транспонировать матрицы, например, в плоскостях y = const, произведя обмен данными между процессорами, и далее сделать заключительное БПФ по направлению z. В этом случае само БПФ уже не потребует обмена данными. В работе использован второй вариант. Один дополнительный процессор отведен для синхронизации коллективных операций и операций ввода-вывода.
Описанный выше алгоритм реализован на языке Fortran-95 и протестирован на многочисленных аналитических примерах. Общий подход в проводимых тестах состоял в следующем: рассмотрим произвольное параболическое уравнение вида:
ff = Й ( X ), (15)
где Й - известный дифференциальный оператор, а вектор X = (T, V, B). Пусть X0 = X(t0). Вычислим
h = Й (X0(x, y, z)) аналитически. Добавим в правую
1 Строго множитель k2 справедлив для однородного, изотропного поля.
часть (15) источник g = -h. Интегрируя полученное уравнение численно необходимо убедиться, что для произвольного начального X0, удовлетворяющего постановке задачи, имеем = 0. Отметим, что для подобных тестов удобно использовать пакеты символьной алгебры, например, Maple или Mathematica, позволяющие преобразовывать полученные аналитические выражения для g в тексты программы на языках Fortran или C. Используемый подход существенно облегчает поиск ошибок, позволяя находить те точки в пространстве, где невязка с аналитическим решением больше точности, определяемой схемой интегрирования по времени.
4. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
4.1. Тепловая конвекция без магнитного поля.
Случай без вращения
Режиму без вращения как без магнитного поля (см. обзор [Гетлинг, 1999]), так и полной задаче динамо (см. ссылки в [Jones, Roberts, 2000b; Cattaneo, Emonet, Weis, 2003; Brandenburg, Subramanian, 2005]) посвящ
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.