Физика Земли, 2020, № 2, стр. 148-160

1. ВВЕДЕНИЕ

Физические условия на поверхности Земли существенно зависят от поведения вещества в пространстве, окружающем планету. Солнечный ветер – поток ионизированных частиц, выброшенных из Солнечной короны, – способен оказывать негативное воздействие как на технические устройства (спутники и средства связи), так и непосредственно на самочувствие людей. Поэтому изучение течения плазмы и электромагнитных процессов в космическом пространстве – одна из практически важных задач физики Земли. Эти процессы удовлетворяют системе уравнений идеальной магнитогидродинамики, состоящей из уравнения Эйлера, описывающего динамику течения жидкости и включающего в себя силу Лоренца, характеризующую воздействие магнитного поля на жидкость, и уравнения магнитной индукции, описывающего эволюцию во времени магнитного поля. До настоящего времени важнейшими алгоритмами для численного решения задач идеальной гидро- и магнитогидродинамики остаются численные схемы фиксированного относительно малого порядка (такие, например, как схемы Рунге–Кутта).

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

До сих пор неизвестно, имеет ли уравнение Эйлера (и даже уравнение Навье–Стокса, описывающее течение вязкой жидкости) глобальное (т.е. существующее при сколь угодно больших временах) решение, и единственны ли решения при данных начальных данных. Сложность этих фундаментальных вопросов связана с тем, что неизвестно, не теряют ли они с течением времени свою гладкость, или у решений возможно развитие коллапса (т.е. возникновение сингулярности) за конечное время. Для простоты далее будем рассматривать только решения периодические в пространстве. Из известного к настоящему времени, значительного массива результатов о решениях бессиловых уравнений Эйлера (см., например, [Majda, Bertozzi, 2002]) отметим следующие: если в начальный момент течение – аналитическое (по пространству) векторное поле, то оно остается аналитическим по пространственным переменным и времени, как минимум на конечном промежутке времени строго положительной длины [Benachour, 1976], а если начальное течение двумерно, то эта аналитичность гарантирована в любой момент времени [Bardos et al., 1976]. Гладкость решения влечет его единственность.

Если на идеальную жидкость не действуют никакие силы, то аналитичность по времени траекторий жидких частиц можно гарантировать на конечном интервале, даже когда начальное течение ${{{\mathbf{V}}}^{{({\text{in}})}}}$ имеет только минимальную гладкость, например, если ряд Фурье завихренности $\nabla \times {{{\mathbf{V}}}^{{({\text{in}})}}}$ сходится абсолютно [Zheligovsky, Frisch, 2014]. Уравнение Эйлера в лагранжевых переменных можно тогда один раз проинтегрировать по времени; так устанавливается факт независимости от времени т.н. инвариантов Коши [Cauchy, 1827]. Используя этот результат, в лагранжевых переменных можно построить разложение поля смещений частиц в ряд Тейлора по времени, имеющий заведомо ненулевой радиус сходимости. Это дало возможность разработать новый метод [Podvigina et al., 2016] для численного интегрирования бессилового уравнения Эйлера, в котором шаг по времени выполняется посредством суммирования ряда Тейлора.

Принадлежность решений уравнений идеальной магнитогидродинамики различным классам Жевре по пространственным эйлеровым переменным рассмотрена в работе [Cheng, Xu, 2019]; из принадлежности функции классу Жевре индекса $\alpha = 1$ следует его пространственная аналитичность. В следующем разделе мы даем простое оригинальное доказательство пространственной аналитичности решений уравнений идеальной магнитогидродинамики при условии аналитичности начальных данных. Оно основано на применении метода, предложенного в статье [Zheligovsky, 2011] и использующего идеи работ [Foias, Temam, 1989; Levermore, Oliver, 1997].

Из этого результата следует аналитичность траекторий жидких частиц и магнитного поля по времени. Как и при отсутствии магнитного поля в рассматриваемой физической системе, используя пространственные лагранжевые переменные, траектории жидких частиц можно разложить в ряд Тейлора по времени (магнитное поле оказывается подчиненным неизвестным полем, т.к. оно обладает свойством вмороженности [Моффат, 1980]). Из аналитичности решения по времени следует сходимость этого ряда. Таким образом, оказывается возможным сформулировать алгоритм для решения магнитогидродинамической задачи, аналогичный алгоритму для чисто гидродинамического случая. Отметим, что пространственная аналитичность начальных данных теперь оказывается существенным условием, т.к. в отличие от случая решений бессилового уравнения Эйлера, из-за наличия в уравнении, описывающем динамику течения жидкости, силы Лоренца, каждый последующий коэффициент в разложении решения в ряд Тейлора по времени зависит теперь от вторых пространственных производных предыдущих коэффициентов. Это явление “потери производных” оказывается существенным усложнением задачи – в чисто гидродинамическом случае его не было.

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

2. АНАЛИТИЧНОСТЬ РЕШЕНИЙ УРАВНЕНИЙ ИДЕАЛЬНОЙ МАГНИТОГИДРОДИНАМИКИ

В этом разделе мы выводим оценки для норм Жевре–Соболева решений уравнений идеальной магнитогидродинамики в эйлеровых переменных на некотором конечном интервале времени $[0,t*)$ в предположении, что начальные данные – аналитические функции пространственных переменных.

Течение жидкости, скорость которой ${\mathbf{V}}({\mathbf{x}},t),$ в присутствии магнитного поля ${\mathbf{B}}({\mathbf{x}},t)$ описывается в эйлеровых пространственных координатах уравнениями:

(считаем, что на жидкость действует только магнитная сила Лоренца), где $P$ – давление. Жидкость предполагается несжимаемой, а магнитное поле соленоидально:

В начальный момент (при $t = 0$) скорость жидкости ${{{\mathbf{V}}}^{{({\text{in}})}}}$ и магнитное поле ${{{\mathbf{B}}}^{{({\text{in}})}}}$ заданы. Считаем, что поля периодичны в пространстве с кубом периодичности ${{\mathbb{T}}^{3}} = {{[0,2\pi ]}^{3}}$ (наши результаты легко распространить на случай различных периодов по разным декартовым переменным) и имеют нулевые средние. Будет использовано их разложение в ряды Фурье:

Поскольку поля действительны, (черта сверху означает комплексное сопряжение). Условия соленоидальности (1.3)–(1.4) принимают вид Норма Жевре–Соболева пространственно-периодического векторного или скалярного поля ${\mathbf{f}}$ с коэффициентами Фурье $\mathop {\widehat {\mathbf{f}}}\nolimits_n $ определяется для $\sigma > 0$ равенством где в случае векторного поля $\left| {\, \cdot \,} \right|$ обозначает евклидову норму вектора. Поле ${\mathbf{f}},$ норма ${{\left| {\, \cdot \,} \right|}_{{{\sigma },q}}}$ которого конечна для некоторого $\sigma > 0,$ аналитично (см. [Levermore, Oliver, 1997]), а $\sigma $ тогда является оценкой снизу для ширины области аналитичности вокруг действительного подпространства в комплексном пространстве ${{\mathbb{C}}^{3}}.$

При $\sigma = 0,$ (5) определяет норму в пространстве Соболева ${{H}_{q}}({{\mathbb{T}}^{3}}),$ которую обозначим ${{\left\| \cdot \right\|}_{q}}.$ Пусть ${{\left| \cdot \right|}_{q}}$ обозначает норму в пространстве Лебега ${{L}_{q}}({{\mathbb{T}}^{3}}).$ По теореме вложения для соболевских пространств [Берг, Oliver, 1980; Трибель, 1980], для любого положительного $q < 3{\text{/}}2$ существует константа ${{C}_{q}}$ такая, что для любой функции $f \in {{H}_{q}}({{\mathbb{T}}^{3}})$ выполнено неравенство

Обозначим ${{\mathcal{P}}_{{\mathbf{n}}}}$ линейный оператор проекции трехмерного вектора на плоскость, нормальную к ${\mathbf{n}} \ne 0{\text{:}}$

Эволюция коэффициентов Фурье скорости жидкости и магнитного поля описывается уравнениями

полученными подстановкой рядов Фурье (2) в уравнения (1).

Будем считать, что при некотором $\sigma > 0$ и $s > 3{\text{/}}2$ начальные данные имеют конечные нормы ${{\left| {{{{\mathbf{V}}}^{{({\text{in}})}}}} \right|}_{{{\sigma },s}}}$ и ${{\left| {{{{\mathbf{B}}}^{{({\text{in}})}}}} \right|}_{{{\sigma },s}}}.$ Следуя [Zheligovsky, 2011], рассмотрим замену переменных

где а g(p) – монотонно строго убывающая функция, положительная при любом $\rho > 0,$ причем

Легко показать, используя монотонность $g,$ что

Наша цель – доказать, что ${\mathbf{v}}$ и ${\mathbf{b}}$ ограничены в ${{H}_{s}}({{\mathbb{T}}^{3}})$ на некотором временном отрезке $[0,t*),$ где $t* > 0.$ Подстановка (8) приводит систему (7) к виду:

Скалярно умножая (9.1) и (9.2) на ${{\left| {\mathbf{n}} \right|}^{{2s}}}{{{\mathbf{\hat {v}}}}_{{ - {\mathbf{n}}}}}$ и ${{\left| {\mathbf{n}} \right|}^{{2s}}}{{{\mathbf{\hat {b}}}}_{{ - {\mathbf{n}}}}}$ соответственно и суммируя оба полученные равенства по всем ${\mathbf{n}} \ne 0,$ находим:

Если преобразовать последнее слагаемое, используя свойства векторного произведения и ортогональность (4) и делая замену индексов суммирования ${\mathbf{n}} \to - {\mathbf{k}},$ ${\mathbf{k}} \to {\mathbf{n}} - {\mathbf{k}}{\text{:}}$

то легко показать, что при $s = \sigma = 0$ сумма по n и k в (10) равна 0; полученное равенство выражает закон сохранения энергии.

По неравенству треугольника, показатель экспоненты в сумме по n и k в (10) неположителен, и она не превосходит 1. Оценим эту сумму:

Каждое из четырех слагаемых оценим аналогично первому из них. Очевидно,

где – скалярные функции. Согласно неравенству Коши–Буняковского–Шварца, для любой функции $f = \sum\nolimits_{\mathbf{n}} {{{{\hat {f}}}_{{\mathbf{n}}}}{{{\text{e}}}^{{{\text{i}}n \cdot {\mathbf{x}}}}}} $ с нулевым средним ${{\hat {f}}_{0}} \ne 0$ и любого $s > 3{\text{/}}2,$

Используя неравенство Гельдера, теорему вложения Соболева (6) и равенство Парсеваля, продолжим при $s > 3{\text{/}}2$ оценку (12):

Таким образом, получаем из (10):

откуда,

(Строго говоря, это неравенство выведено в предположении конечности норм ${{\left\| {{\mathbf{v}}({\mathbf{x}},t)} \right\|}_{{s + 1{\text{/}}2}}}$ и ${{\left\| {{\mathbf{b}}({\mathbf{x}},t)} \right\|}_{{s + 1{\text{/}}2}}}.$ Точное доказательство состоит в установлении неравенства (14) для приближений Фурье–Галеркина к решению – для них, очевидно, указанные нормы конечны, – а затем применении методов компактности.) Решая дифференциальное неравенство (14), находим

Это неравенство устанавливает контроль над нормами решения $p(t)$ на конечных временах $t < t{\text{*}}.$ Для оценки $t*,$ рассмотрим (15) для функций ${{g}_{{\gamma }}}(p) = {{({{p(0)} \mathord{\left/ {\vphantom {{p(0)} p}} \right. \kern-0em} p})}^{{\gamma }}}.$ В силу (8.3), функция $p(t)$ при $t > 0$ также зависит от $\gamma ,$ что будем обозначать нижним индексом $\gamma .$ Интегрируя (15) по времени от 0 до $t,$ получаем

при

Отсюда следует

Таким образом, решение ${\mathbf{V}}(x,t),$ ${\mathbf{B}}(x,t)$ аналитично в слое ширины

при условии (16). Поскольку эти оценки выполнены при любом $\gamma > 0$ и правые части (16) и (17) растут при увеличении $\gamma ,$ получаем в пределе $\gamma \to \infty ,$ что имеет место

Теорема. Пусть поля ${{{\mathbf{V}}}^{{({\text{in}})}}}$ и ${{{\mathbf{B}}}^{{({\text{in}})}}}$ аналитичны в слое ширины $\sigma > 0$ и нормы ${{\left| {{{{\mathbf{V}}}^{{({\text{in}})}}}} \right|}_{{{\sigma },s}}}$ и ${{\left| {{{{\mathbf{B}}}^{{({\text{in}})}}}} \right|}_{{{\sigma },s}}}$ конечны при некотором $s > 3{\text{/}}2.$ Тогда решение системы (1) с такими начальными данными аналитично в слое ширины $\sigma (t) = ({{1 - t} \mathord{\left/ {\vphantom {{1 - t} {t{\text{*}}}}} \right. \kern-0em} {t{\text{*}}}})\sigma $ при любом $t < t*,$ где

Решения допускают оценку

Отметим, что построенные оценки времени и размера области аналитичности грубы, т.к. в них не учтено существенное сокращение, имеющее место, например, при вычислении векторного произведения в (1.2).

Коэффициенты разложения в конечный ряд Фурье приближения Фурье–Галеркина к решениям исходной системы (1) удовлетворяют автономной системе обыкновенных дифференциальных уравнений конечной размерности, правые части которых аналитически зависят от неизвестных коэффициентов Фурье. По теореме Коши [Петровский, 1970] эти коэффициенты являются аналитическими функциями времени. Гармоники Фурье аналитичны по пространственной переменной, поэтому каждое приближение аналитично по паре переменных $({\mathbf{x}},t)$. Поскольку (как уже указано выше) для этих приближений верны равномерные по номеру приближения оценки, аналогичные выведенным, их можно аналитически продолжить во всю область

Опираясь на выведенные априорные оценки, легко показать, что приближения Фурье–Галеркина сходятся к решению исходной системы (1) равномерно в любой компактной подобласти области аналитичности $\mathfrak{W},$ поэтому по теореме Вейерштрасса [Шабат, 1976] и их предельные поля ${\mathbf{V}}({\mathbf{x}},t)$ и ${\mathbf{B}}({\mathbf{x}},t)$ – аналитические функции пары переменных в $\mathfrak{W}.$

Лагранжевы координаты ${\mathbf{a}}$ частицы жидкости – это ее декартовы координаты в начальный момент времени. Траектории частиц ${\mathbf{x}}({\mathbf{a}},t)$ удовлетворяют обыкновенному дифференциальному уравнению:

поэтому по теореме Коши при $t < t{\text{*}}$ траектории и решения ${\mathbf{V}}({\mathbf{x}}({\mathbf{a}},t),t),$ ${\mathbf{B}}({\mathbf{x}}({\mathbf{a}},t),t)$ системы (1) в лагранжевых координатах также аналитичны по паре переменных $({\mathbf{a}},t).$

3. РЕКУРРЕНТНЫЕ СООТНОШЕНИЯ ДЛЯ КОЭФФИЦИЕНТОВ РЯДОВ ТЕЙЛОРА ПО ВРЕМЕНИ

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

Наиболее просто они выглядят для решений в эйлеровых координатах. Именно, для коэффициентов рядов Тейлора по времени

имеем для $n > 0{\text{:}}$ и для $n \geqslant 0$ выполнены условия соленоидальности: (коэффициенты ${{P}_{n}}$ для давления определяются из условия (20.3)).

В лагранжевых координатах уравнения (1.1) и (1.2) принимают соответственно вид:

Точка обозначает т.н. материальную производную, т.е. производную по времени при фиксированной лагранжевой координате: а градиент $\nabla ,$ как и ранее, производные по эйлеровым пространственным переменным.

Уравнение (21.2) имеет следующее решение в лагранжевых переменных (см. [Моффат, 1980]):

где ${\mathbf{x}}({\mathbf{a}},t)$ – решение обыкновенного дифференциального уравнения (18), описывающее траекторию частицы с лагранжевой координатой ${\mathbf{a}}$. Подставив его в (21.1), воспользуемся для выведения рекуррентных соотношений процедурой, использованной в работе [Zheligovsky, Frisch, 2014]. Умножив (21.1) скалярно на ${{\partial {\mathbf{x}}({\mathbf{a}},t)} \mathord{\left/ {\vphantom {{\partial {\mathbf{x}}({\mathbf{a}},t)} {\partial {{a}_{m}}}}} \right. \kern-0em} {\partial {{a}_{m}}}},$ получаем:

После взятия ротора по лагранжевым переменным от обеих сторон этого уравнения градиент справа исчезает, и уравнение принимает вид:

(символ ${{\nabla }^{{\text{L}}}}$ обозначает градиент по лагранжевым переменным a). В случае бессилового уравнения Эйлера, т.е. при отсутствии магнитного поля (${{{\mathbf{B}}}^{{({\text{in}})}}} = 0$), полученное уравнение можно один раз проинтегрировать по $t;$ так устанавливается независимость от времени выражений, известных как инварианты Коши ([Cauchy, 1827], см. также [Zheligovsky, Frisch, 2014]).

Рассмотрим ряд Тейлора по времени для траекторий частиц жидкости:

Здесь, очевидно,

Подстановка ряда (24) в (23) дает для $s \geqslant 0$ равенство при степени ${{t}^{s}}{\text{:}}$

позволяющее найти ${{\nabla }^{L}} \times {{{\mathbf{x}}}^{{(s + 2)}}},$ если известны все ${{{\mathbf{x}}}^{{(n)}}}$ с меньшими индексами. При s = 0 оно существенно упрощается:

Условие несжимаемости (1.3) эквивалентно равенству $det\left\| {{{\partial {{x}_{i}}} \mathord{\left/ {\vphantom {{\partial {{x}_{i}}} {\partial {{a}_{j}}}}} \right. \kern-0em} {\partial {{a}_{j}}}}} \right\| = 1,$ из которого было получено (см. [Zheligovsky, Frisch, 2014]) соотношение:

Здесь индексы суммирования положительны: $n \geqslant 1,$ $m \geqslant 1,$ $l \geqslant 1,$ а ${{\varepsilon }_{{ijk}}}$ – единичный антисимметричный тензор.

Итак, для любого индекса $s \geqslant 2$ можно вычислить ротор (25) и градиент (26) коэффициента Тейлора с этим индексом, если известны все коэффициенты с меньшими индексами. Соответственно, можно вычислить

и, обращая лапласиан ${{\Delta }^{L}},$ сам коэффициент.

Таким образом, можно сформулировать два метода интегрирования системы уравнений (1) по времени, использующие аналитичность решений по времени. Первый из них (который естественно назвать эйлеровым) предполагает построение рядов Тейлора (19) для решений в эйлеровых переменных, используя соотношения (20) и суммируя на каждом шагу по времени эти ряды Тейлора. Второй (который назовем лагранжевым) отличается тем, что опирается на аналитичность решений в лагранжевых переменных. В этом алгоритме приведенные выше рекуррентные соотношения позволяют последовательно вычислять коэффициенты ряда Тейлора по времени (24) для перемещения частиц жидкости ${\mathbf{x}}({\mathbf{a}},t),$ и далее следовать плану [Podvigina et al., 2016] решения бессилового уравнения Эйлера: суммируя ряд Тейлора (24), определяем перемещение частиц, в начальный момент времени находящихся в узлах регулярной пространственной сетки; в конечных точках перемещения находим скорости частиц (суммируя ряд Тейлора (24), продифференцированный по времени) и магнитное поле (22), и, используя каскадную интерполяцию ([Purser, Leslie, 1991; Podvigina et al., 2016]), находим решение ${\mathbf{V}}({\mathbf{x}},t)$ и ${\mathbf{B}}({\mathbf{x}},t)$ в новый момент времени на регулярной сетке. Эти шаги можно повторять, пока не будет достигнуто время, на которое желательно проинтегрировать систему (1), и пока решение не продемонстрирует сингулярное поведение, утеряв пространственную аналитичность (сохраняя или не сохраняя при этом пространственную гладкость).

4. ПРИМЕР РАСЧЕТА ЗАДАЧИ ИДЕАЛЬНОЙ МАГНИТОГИДРОДИНАМИКИ ЛАГРАНЖЕВЫМ АЛГОРИТМОМ

Мы провели вычисления для случайным образом выбранных начальных гладких трехмерных полей:

Из двух рассмотренных алгоритмов был выбран лагранжев, поскольку тестовые расчеты, выполненные в работе [Podvigina et al., 2016], показали, что в применении для чисто гидродинамической задачи этот алгоритм более эффективен (позволяя провести расчет решения на заданном интервале времени за меньшее время процессора). Вычисление коэффициентов временного ряда Тейлора с использованием рекуррентных соотношений (25) и (26) проводили стандартным псевдоспектральным методом с разрешением ${{256}^{3}}$ гармоник Фурье. На каждом шагу по времени суммировали 6 или 12 членов ряда Тейлора. Временной шаг выбирался автоматически из условия, что среднеквадратичная величина последнего суммируемого слагаемого (векторного поля, зависящего от лагранжевых переменных) в обрезанном ряду Тейлора равна заданной погрешности ${{10}^{{ - 12}}}.$

Расчеты проводили до времени $t = 2.$ На рис. 1 и рис. 2 показаны изоповерхности кинетической и магнитной энергий, соответственно, с шагом $ \sim 0.5$ единиц времени. (Рисунки построены для моментов окончания шагов интегрирования, в которых были превышены целевые временные пороги, кратные 0.5, поэтому эти времена не кратны 0.5 точно.) На рис. 1 видно, что с течением времени в поле скорости образуются двумерные структуры и концентрированные сгустки высокой кинетической энергии течения. В магнитном поле (рис. 2) в процессе эволюции возникают одномерные структуры (жгуты). Это согласуется с существующими представлениями о развитии в гидродинамических и магнитогидродинамических системах структур меньшей размерности.

В табл. 1 приведены данные об изменении с течением времени максимальной плотности кинетической, $E_{{max}}^{{\text{K}}},$ и магнитной, $E_{{max}}^{{\text{M}}},$ энергий и пространственные средние этих плотностей, ${{\bar {E}}^{{\text{K}}}}$ и ${{\bar {E}}^{{\text{M}}}}$ соответственно. Максимальная плотность магнитной энергии ко времени $t = 2$ увеличивается в 1.9 раз, а кинетической энергии – в 5.4 раза (этим объясняется относительная малость объема, ограниченного приведенными на рис. 1 изоповерхностями при $t \geqslant 1.5$). Рост максимумов энергий на фоне мало меняющихся их пространственных средних также свидетельствует о формировании контрастных пространственных структур. Проверено выполнение сохранения полной энергии ${{({{{\left| {\mathbf{V}} \right|}}^{2}} + {{{\left| {\mathbf{B}} \right|}}^{2}})} \mathord{\left/ {\vphantom {{({{{\left| {\mathbf{V}} \right|}}^{2}} + {{{\left| {\mathbf{B}} \right|}}^{2}})} 2}} \right. \kern-0em} 2}{\text{:}}$ с течением времени эта величина несколько уменьшается (в связи с тем, что использованное в расчете пространственное разрешение ограничено и поэтому не может правильно описать развитые процессы прямого каскада энергии), однако это уменьшение не превышает 0.006 при $t = 2.$

Из рекуррентных соотношений (20) и (25) следует, что (в отличие от (26) и случая лагранжевых рекуррентных соотношений (25) для бессилового уравнения Эйлера) каждый коэффициент зависит от пространственных производных коэффициентов меньших индексов первого (20) или второго (25) порядка, и, таким образом, при вычислении последовательности коэффициентов неявно используются производные коэффициентов меньших индексов все более высоких порядков. Поскольку вычисление производных – численно неустойчивая процедура (неустойчивость которой может усугубляться ошибками, возникающими при интерполировании полей после завершения шага по времени с равномерной сетки лагранжевых переменных на равномерную сетку эйлеровых переменных), естественно предположить, что уже при вычислении коэффициентов с достаточно малыми индексами возможна катастрофическая потеря точности. Чтобы проверить, насколько этот эффект существенен, расчеты проведены 4 раза: была использована стандартная так называемая “двойная” точность представления чисел с плавающей запятой (данные типа real*8 Фортрана) и “четверная” (данные типа real*16), а также суммировались 6 или 12 членов ряда Тейлора по времени. При этом существенного различия в результатах расчетов не обнаружено (несмотря на то, что векторные поля – коэффициенты ряда Тейлора – оказываются велики по норме, что качественно соответствует ожидаемому проявлению влияния многократного повторного дифференцирования, возникающего в рекуррентных соотношениях при наличии магнитного поля).

На рис. 1 и рис. 2 видно, что при $t = 1.5073$ в зонах большой кривизны изоповерхностей энергии магнитного поля на них возникает “бахрома”; в меньшей степени этим поражены изоповерхности кинетической энергии. При $t = 2.0001$ это явление имеет уже весьма развитый характер – “мелкой рябью” покрыты многие изоповерхности и вдалеке от их “кромок”. Это характерное проявление недостаточности пространственного разрешения, использованного в вычислениях. Именно недостаточное разрешение, а не ограниченная точность вычислений оказывается главным препятствием для продолжения численного интегрирования уравнений идеальной магнитогидродинамики. Вместе с тем, эти два явления взаимосвязаны: численное дифференцирование рядов Фурье сводится к умножению коэффициентов на соответствующие волновые числа, а при использовании в вычислениях большего числа гармоник Фурье для увеличения пространственного разрешения растет величина этих сомножителей, что предъявляет повышенные требования к точности.

5. ЗАКЛЮЧЕНИЕ

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

Используя это соотношение и основываясь на аналитичности решения, мы построили разложение поля скорости течения в лагранжевых переменных в ряд Тейлора по времени. Для коэффициентов этого ряда для скорости, представленной в лагранжевых координатах, выведены рекуррентные соотношения (25), обобщающие рекуррентные соотношения для решений бессилового уравнения Эйлера [Zheligovsky, Frisch, 2014; Podvigina et al., 2016]. Для решений, представленных в эйлеровых переменных, коэффициенты ряда Тейлора удовлетворяют равенствам (20). Сформулированы лагранжев и эйлеров алгоритмы для численного интегрирования уравнений идеальной магнитогидродинамики, основанные на использовании этих рекуррентных соотношений. В принципе, они точнее стандартных алгоритмов, т.к. имеют произвольный порядок точности; какие порядки реально достижимы в вычислениях, зависит от имеющихся в наличии компьютерных ресурсов (прежде всего, оперативной памяти).

Лагранжев алгоритм опробован в расчетах. Результаты показывают, что, в согласии с существующими представлениями, в процессе временной эволюции в течении и магнитном поле образуются структуры меньших размерностей, например, жгуты магнитного поля. С течением времени энергия передается в структуры малого пространственного масштаба (прямой каскад), в решении растут величины пространственных градиентов. Для корректного описания этих процессов необходимо проводить расчеты с достаточно большим пространственным разрешением. Ограниченность используемого разрешения не позволила нам продвинуться далее $t = 2.$ В связи с этим представляется целесообразным в последующей работе исследовать вопрос, какой из двух описанных алгоритмов менее чувствителен к ограниченному пространственному разрешению на больших временных интервалах.

Успех применения лагранжева алгоритма для численного решения бессилового уравнения Эйлера основан на том, что при переходе к лагранжевым координатам в уравнении исчезает нелинейное адвекционное слагаемое $({\mathbf{V}} \cdot \nabla ){\mathbf{V}}.$ Сила Лоренца $({\mathbf{B}} \cdot \nabla ){\mathbf{B}},$ имеющая аналогичную структуру, в уравнении остается, что имеет ряд негативных с вычислительной точки зрения последствий, главное из которых – эффект потери производных. Как известно, в переменных Эльзаcсера [Elsasser, 1950] ${{{\mathbf{S}}}^{ \pm }} = {\mathbf{V}} \pm {\mathbf{B}}$ уравнения идеальной магнитогидродинамики (1) принимают вид

т.е. каждое из этих векторных полей “переносится потоком”, скорость которого – векторное поле, являющееся другой переменной Эльзассера. Таким образом, представляла бы большой интерес разработка алгоритма, в котором использованы две “лагранжевы” системы координат, отвечающие двум потокам ${{{\mathbf{S}}}^{ \pm }},$ при переходе в которые соответствующие переносные члены в уравнениях (27) пропадают, и две эти “системы координат Лагранжа–Эльзассера” синхронизируюся с помощью каскадной интерполяции при завершении временного шага.