Построение численной модели тонкой магнитной трубки с помощью гамильтонова формализма Текст научной статьи по специальности «Физика»

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

Том 13, № 3, 2008

Построение численной модели тонкой магнитной трубки с помощью гамильтонова формализма

Д. В. РОМАНОВ

Иркутский государственный университет путей сообщения, Россия e-mail: [email protected]

К.В. Романов

Красноярский институт железнодорожного транспорта, Россия

A novel approach to the construction of conservative schemes for a model case of the evolution of a thin magnetic tube is applied. It has been shown how to introduce an approximating discrete system which is used to derive the corresponding Lagrangian. As a result, robust up-to-date numerical schemes can be applied to such a system to provide a very accurate description of long-term evolution. The potential pitfalls and solutions are highlighted.

Введение

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

Пока это единственная из существующих моделей, дающая возможность исследовать подъем магнитного поля от основания конвективной зоны в атмосферу Солнца (в настоящее время уже появились работы, где распределенное поле описывается с помощью упрощенных МГД-моделей, но поддающийся моделированию диапазон глубин очень мал [5, 6]). Из-за того что совместное действие гравитации и конвекции действительно структурирует магнитное поле в тонкие (диаметром до 300 км на поверхности Солнца) жгуты c интенсивностью поля внутри на 3-4 порядка выше фоновой [4], разработка более точных моделей тонких магнитных трубок по-прежнему актуальна.

До настоящего времени потеря полем устойчивости и подъем в атмосферу были изучены с акцентами на временные характеристики начального этапа данного процесса (определение инкрементов неустойчивостей, времен насыщения, времени нелинейного подъема) и на определение связи между исходными параметрами и измеримыми на уровне поверхности Солнца величинами (структура поля в ведущем и ведомом

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

пятнах, угол наклона биполярной активной области и т.п.). Использовались конечно-разностные схемы, при разработке которых внимание уделялось исключительно устойчивости и экономичности [3, обзор лит-ры]. Тестирование проводилось на линейных задачах с небольшими амплитудами, и законам сохранения не уделялось внимания — целью было воспроизведение инкрементов неустойчивостей для определения времени удержания поля, интерпретируемого как ограничение сверху на время работы динамо-процесса [3].

Перед численным исследованием циклической активности, долговременной эволюции или генерации магнитного поля за счет взаимодействия с конвективными потоками плазмы следует пересмотреть используемые схемы с целью избежать усиления поля за счет численных эффектов при многократных прохождениях области в присутствии значительного (больше 10 порядков) перепада параметров внутри Солнца. Известно, что консервативные схемы (соблюдающие точный баланс сохраняющихся физических величин) можно построить, например, записав систему уравнений в дивергентном виде: это позволяет выписать потоки величин через границы ячеек сетки и получить разностные выражения для соответствующих компонент энергии. При использовании этого подхода конкретные выражения для различных видов энергии выбираются специально, ради соблюдения баланса в разностном виде [7]. Другим подходом стало использование переопределенных совместных систем уравнений и двухэтапный шаг "предиктор— корректор" с использованием дополнительных законов сохранения на этапе коррекции (метод Годунова). Подобные методы были успешно использованы в численной газодинамике, при решении уравнений МГД и Максвелла [8].

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

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

Достоинствами такого подхода являются физическая чистота постановки, полностью совместные определения всех видов энергии через гамильтониан, существование большого набора высокоточных методов интегрирования гамильтоновых систем, разработанных в смежных областях [9] (задачи небесной механики [10], вычислительной химии [11], численного решения волновых уравнений, уравнений Клейна—Гордона и многое другое). Также открывается возможность улучшенного описания поперечного распределения плазмы трубки без сложных выкладок при усреднении по сечению — дополнительные параметры сразу же перейдут в переменные гамильтониана. Еще заметим, что в литературе описаны продуктивные методы анализа численных схем, когда

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

Структура статьи следующая: в разд. 1 описана дискретная система, приведены гамильтониан системы и выражения для обобщенных импульсов. В разд. 2 найдена (с учетом всех голономных связей) обобщенная сила и получены уравнения движения. В разд. 3 выписаны схемы второго и четвертого порядков для интегрирования уравнений движения как системы обыкновенных дифференциальных уравнений.

Вторая часть работы посвящена исследованию полученных конечно-разностных схем. Для компактного представления последовательности шагов, используемых обозначений и критерия устойчивости системы аналитически исследована устойчивость тонкой магнитной трубки в стратифицированной плазме (разд. 4). Затем аналогичным образом выполнен анализ устойчивости схемы с перешагиванием (разд. 5) и симплекти-ческой схемы четвертого порядка (разд. 6). После численного моделирования нелинейных осцилляций и анализа результатов дано заключение, суммирующее все полученные результаты.

Обоснование используемого приближения и независимый вывод системы уравнений тонкой магнитной трубки в частных производных приведены в работах [1, 3, 12].

1. Гамильтониан дискретной системы

Для иллюстрации метода построим дискретную систему, близкую к использованным ранее [1,2], — замкнутую в кольцо цепь из набора сегментов цилиндрической формы с параметрами плазмы, однородными внутри каждого цилиндра (рис. 1). Лагранжевыми переменными данной системы являются координаты г и скорости vi концов сегментов. Параметры среды (такие как напряженность поля Н+1/2, давление р^+1/2, плотность Р-1+1/2 и площадь сечения трубки а^+1/2) определяются из условий адиабатичности, баланса давлений в полуцелых узлах и законов сохранения магнитного потока и массы [1]. Для простоты изложения концы цилиндров будут индексироваться как целые (а их середины — как полуцелые) узлы некоторой сетки.

Лагранжиан состоит из пяти групп слагаемых: кинетической энергии, внутренней энергии плазмы, энергии магнитного поля, потенциальной энергии в поле тяжести и

Рис. 1. Тонкая магнитная трубка и аппроксимирующий ее набор однородных цилиндров. Кружки отмечают концы цилиндров, координаты и скорости которых являются лагранжевыми переменными системы. Параметры плазмы внутри цилиндров задаются условиями сохранения и условием баланса давлений в полуцелых точках (помеченных крестиками)

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

Обезразмерив магнитное поле на р0/8п, чтобы избавиться от коэффициента 1/8п в дальнейших выкладках, получим лагранжиан системы в безразмерных переменных:

ь = ^ к

г,г+1 — ¿¿,¿+1

иг,г+1 = ^¿+1/2 Гг+1

Кг,г+1 =

^З^2 + #2+1/2 + Ре^(Гг+1/2)

+ Мг+1/2ф(Гг+1/2),

^(г)|2

Рг+1/2^г+1/2 -^- (1г+1/2 ' ¿г),

¿¿+1/2 = |Г^+1 - Гг|, I = (г,+1 - г,)/ь,

(1) (2)

(3)

(4)

где Ь — лагранжиан системы, К^+1 — кинетическая энергия сегмента, £/¿,¿+1 — потенциальная энергия сегмента, ^¿+1/2 — объем сегмента, М,+1/2 — масса сегмента, ф(г) — гравитационный потенциал, ^¿+1/2 — единичный вектор вдоль оси, ¿¿+1/2 — длина сегмента.

Индекс г меняется от 0 до N — 1, все величины периодичны по г с периодом N. Скорость у(г) меняется линейно, следовательно,

к

м,

¿,¿+1

¿+1/2

2

(avi + (1 — = Mi+l/2

V

2 + ^+1 + vi+l) 6 '

(5)

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

#¿+1/2 ^¿+1/2 = Ф = СОПЭ^ ^¿+1/2 + #¿+1/2 = (г,+1/2) ,

1/2 ¿¿+1/2 0Ч+1/2 = ^¿+1/2 , р,+1/2/р7+1/2 = ^¿+1/2 •

(6) (7)

Для перехода к гамильтонову формализму найдем обобщенный импульс р, = Ь и гамильтониан Н = ^, р,г, — Ь:

д^

Ем

v2+3 + к-,

3+1/2"

м,

¿+1/2

2vi + Vi+l 6

+ М,_1/2

2vi + Vi_l

6

(8)

н

Е

ь = ^ Ки+1 + /¿,¿+1]

(9)

2. Обобщенная сила

Наличие связей при получении уравнений движения учитывается многими способами (например, с помощью множителей Лагранжа [13]), но в нашем случае возможно значительно упростить выкладки, так как требуется найти исключительно градиент потенциальной энергии УГг / при наличии налагаемых связями ограничений.

1

Р,

6

Потенциальная энергия (2) является функцией от параметров плазмы сегмента и координат концов, а параметры вещества в свою очередь удовлетворяют системе (6)-(7). Опустим для краткости индекс г + 1/2 и рассмотрим уравнения (6)-(7) подробнее. Если обозначить г¿+1/2 как И,, вся система подстановками

Н = Ф/а, р=М/(аЬ), р = С (М/(аЬ))7 сводится к нелинейному уравнению на площадь сечения

п( м у ф2 с аь) + 02=

(10)

(11)

Это уравнение всегда имеет решение, и только одно, так как справа стоит положительное число, а левая часть монотонно убывает от бесконечности до нуля. С учетом (10) получаем, что параметры плазмы сегмента определяются двумя скалярными функциями от координат концов, Ь и рет = рет4(И). Именно этот факт позволяет сократить объем выкладок.

Предположим, что все параметры в (2) являются известными решениями (6)-(7), выписанными как функции Ь и рет4(И). Тогда несложно вычислить частные производные от параметров плазмы трубки по координатам. Для этого берется дифференциал от (11):

- г с( МУ

а V аЬ) а а2

йа - фет - ьС (

йЬ = 0.

Отсюда находятся да/дгретЛ и да/дЬ, а с помощью (10) и остальные частные производные:

да

а

да

а

1'Р

дретл ур + 2Н2 ' дЬ Ьур + 2Н2 '

дН Н дН Н ур

дретЛ 7р + 2Н2' дЬ Ьур + 2Н2 '

др Р др р 2Н2

дрет1 ~/р + 2Н2' дЬ Ьур + 2Н2 '

др _ тр др 2Н2 Тр

дретл IV + 2Н2' Из (4) дополнительно следует дЬг+1/2

дЬ

Ь ур +2Н2'

(12)

(13)

(14)

(15)

дГ;

-I

дЬ;

¿+1/2'

¿+1/2

дг

¿+1

'¿+1/2'

дрел*(Гг+1/2^ 1

д г

^ (гг+1/2) '

дреж4(Гг+1/2) _ 1

дг

¿+1

2 Vpexí(г¿+1/2)•

(16) (17)

Теперь вернемся в нахождению обобщенной силы дг. Н = дг. и, состоящей из произ-

водных от четырех вкладов в потенциальную энергию:

_д_ дг,-

¿,¿+1

дг,

£

^+1/2 / Р¿+1/2

_ Р¿+1/2 \7 - 1

+ рет^^^Н ^+1/2 + М+1/2ф(гТ+1/2)

(18)

Учитывая, что

дф(гт/2) _ дф(гт/2)

2 vф

- 1 ё (гг+1/2 ),

дгг дг^+1 2 ' г '.¿+1/2 2

последнее слагаемое уравнения (18) можно переписать как

д 1

^Мт/2ф(гг+1/2) = - 2 [М,+1/2ё(г,+1/2) + М, —1/2ё(г, —1/2^ •

дг

(19)

Остальные члены (18) дифференцируются последовательно как функции от двух аргументов рех^(гг + гг+1 )/2) и Ь(гг, гг+1). При вычислении производных от произведения потребуются уже полученные выражения (16), (17) для нахождения производных множителей (индекс г + 1/2 опускается для краткости):

д_ дгг

"М" М 1

_ р. р р

др дрехг + др дь д'РехЬ д г г дь д г г

м

1

р 7Р + 2Н2

2 ^Рвхг

2Н2

¿+1/2

ь

4+1/2

(20)

д_ дгг

Р

7 - 1

+ Рех* (г+1/2) + Н2

2Н2 7 - 2

"■4+1/2

ТР

1

+ "Vpexí I ; 2 .¿+1/2

1+

1

ь

7Р

+

7 - 1 тР + 2Н2 2Н 2

7 - 1 7Р + 2Н2 7Р + 2Н2 Производная от всего произведения имеет более компактный вид:

(21)

А

д гг

Мг+1/2 / Рг+1/2

Рг+1/2 \ 7 - 1

+ Рвх4 (гг+1/2) + Н2+1/2

мг

г+1/2

рг+1/2

1 I. ,

2 .¿+1/2

2Нг+1/2

ь

г+1/2

-г+1/2

причем с помощью (6) последнее слагаемое можно упростить:

■ 1г+1/2 = 2 Ф Нг+1 /2 1г+1/2

Мг+1/2 2Нг2+1/2

рг+1/2 ьг+1/2

После подстановок и суммирования уравнения движения примут законченный вид:

У,

(22)

ф, _ м (2У, + У,-1) + м (2У, + УУ,+1) = М,-1/2-6--+ М,+1/2-

<и

6

2ф ГН , Н , 1 + М,-1/2ё(г,-1/2) + М,+1/2ё(г,+1/2)

М,

,-1/2

2

I М,+1/2 .

^ '"-1/2 2р,+1/2 ^ г .7+1/2

(23)

_ 2р, —1/2

Первое слагаемое обобщенной силы описывает силу натяжения магнитного поля, второе — силу тяжести, последнее — действие внешнего давления. Полученная система обыкновенных дифференциальных уравнений может быть решена аналитически или численно, как описано ниже. Заметим, что обобщенный импульс р, является линейной комбинацией скоростей узлов, поэтому возможно обращение левой части уравнения (23) циклической прогонкой с получением альтернативного уравнения на У,.

.

+

3. Интегрирование уравнений движения

В предыдущем разделе получены аналитические уравнения движения (22), (23) рассматриваемой гамильтоновой системы, образующие систему обыкновенных дифференциальных уравнений. Для интегрирования подобных систем существует богатый набор методов [14], и выбор схемы во многом определяется решаемой задачей — в нашем случае, временным интервалом, на котором требуется получить решение. Как правило, область применимости классических методов (метод Эйлера, Рунге—Кутты и др. [14]) ограничена окрестностью начальной точки, и даже схемы высокого порядка, не гарантирующие сохранения нужных инвариантов, на больших временах дают качественно неверные решения.

В это же время для исследования эволюции гамильтоновых систем на больших временах разработаны симплектические схемы, сохраняющие дополнительный геометрический инвариант и, как следствие, структуру течения в фазовом пространстве системы [9, 15]. Такие схемы продемонстрировали хороший порядок аппроксимации и высокие консервативные свойства на асимптотически больших временах, вследствие чего сейчас активно используются во многих областях [9-11].

Для тестирования используемого метода построения схем система (22), (23) из 2Ы уравнений была проинтегрирована схемой с перешагиванием и симплектической схемой четвертого порядка. Обе схемы известны тем, что на периодических решениях не дают искусственной диссипации/генерации энергии, обратимы и являются явными.

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

ги+1 ги

ri - Гг = «+1/2

т '

и+1/2 и-1/2 (24)

Р — Р

1 i_1 ^ -р (ги „и „и \

= Г пЧ-1' Ii ' 1i+1)'

где ЕДгП^' ги, „и+1) = н обозначает обобщенную силу.

При построении симплектической схемы [15] используется разделимость гамильтониана, который можно представить в виде суммы н(р' г) = К (р) + и (г), вследствие чего силы и скорости зависят только от координат и обобщенных импульсов:

¿и(г) . ¿К(р)

В этой схеме четвертого порядка аппроксимации координаты и импульсы определены на одной сетке, но временной шаг теперь выполняется в четыре этапа:

р(к) = р(к-1) + тЬк¥ (г(к-1))' г(к) = г(к-1) + таку(р(к))

Ъг = 0' Ьз = 1/(1 - 22/3)' «1 = аА = (2 + 21/3 + 2-1/3)/6'

Ь2 = Ь4 = 1/(2 - 21/3)' «2 = аз = (1 - 21/3 - 2-1/3)/6. (25)

Этапы пронумерованы переменной к = 1'..' 4, а начальные и конечные значения определяются следующим образом:

г(0) = ги' ги+1 = г(4)' (26)

р(0) = ри' ри+1 = р(4). (27)

4. Линейный анализ устойчивости

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

Порядок аппроксимации находится разложением точного решения в ряд и подстановкой его в разностную схему. Для равномерных сеток (М^+\/2 = М) обе схемы имеют второй порядок аппроксимации по массовой переменной; порядок аппроксимации по времени для схемы с перешагиванием (24) — второй, а для симплектической схемы (25) — четвертый.

Анализ устойчивости будет проведен методом фон Неймана, т. е. исследованием устойчивости отдельных фурье-гармоник. В качестве равновесного положения выберем горизонтально расположенную трубку с однородными параметрами.

Нам потребуются некоторые предварительные сведения об условии физической устойчивости тонкой магнитной трубки [12, 16], кратко изложенные ниже. В пределе бесконечно малого размера цилиндров уравнения (6), (7), (22), (23) переходят в систему уравнений тонкой магнитной трубки в частных производных, с лагранжевой массовой

г

координатой, отложенной вдоль трубки, в(г) = / ра I ¿г [1]: д г д V д 1

й = V та = 2Н*дга(Н1) +8(г) - рур«(г)'

На = я°а"' £ = р0 • (28)

дг

р + Н2 = РезЛ( г), I = ар—, (I,1) = 1.

дв

Пусть трубка горизонтальна и покоится. Среда находится в состоянии гидростатического равновесия УрехДг) = рех4(г^(г), поэтому система уравнений (28) упрощается [1]:

дг дV т д . тт/,.

« = ^ ж = 2Фд^в(Н1) +

1 _ РвхЛ(г)

g(г),

На = Ф, = Р0, (29)

Р7 Ро

дг

р + Н2 = рехь (г), I = ар—, (I, 1) = 1.

дв

В используемых безразмерных единицах квадрат альфвеновской скорости равен V2 = 2Н2/р, квадрат скорости звука а2 = тр/р, квадрат касповой скорости Б2 = V2а2/(V2 + а2). Направим ось X вдоль оси трубки, а ось У — вертикально вверх и найдем невозмущенную плотность из условия равновесия:

р |4=0= рех1(Уо)-

р

Линеаризованная система уравнений (29) записывается как

¿v, ^ = 2Ф#{8HI + - g^¿y + g^ (30)

dt ' dt ds p p

¿p + 2H5H = p'ext8y = pg5y, ¿1 = (— + — J еж + pa^,

\ p a J ds

(31)

+ — = 0, -р = а2 -р, (-£, 1) = 0. (32)

Н а

Представим возмущения параметров через их комплексные амплитуды, волновой вектор к и частоту ш, т.е. в виде фурье-моды. Например, для возмущения а(в,1)

5а(э,г) = 5агкз/р0,Т0-шЬ.

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

—гш- г = -V, (33)

—гш-ч = гкУ2( -НI + — g ^ -у + g5р, (34)

VН ) р р

г тт г

-+ — = 0, -р = а2-р, -р + 2Н-Н = рд-у, (35)

Н а

-I = ^+ ех + гк-г, (36)

(-1,1) = 0. (37)

Так как -1х = 0 из (37), х-компонента уравнения (36) дает уравнение непрерывности:

+ — + гк-х = 0. (38)

ра

Выразив возмущения -р, -р и -а через -Н с помощью (35) и подставив их в уравнение (38), получим -Н и -р как функции малого смещения:

-Н а2 г., , д г и -р У2 ., г д

H a2 + V2

a2

ik¿x + — ¿y , — = — ik¿x + ¿y.

p a2 + V 2 a2 + V 2

Подстановка -I = гк-у еу и этих выражений в правую часть (34) позволяет переписать (33)-(34) как систему волновых уравнений:

—ш2 -х = —к2Б2-х + гк одУ2 „-у, (39)

а2 + У2

—ш2 -у = —к2У2-у — П2-у — гк^— -х, П2 = + д^. (40)

а2 + У2 а2 + У2 р

Отметим, что П2 является квадратом частоты Брента—Вяйсяля и его отрицательность сигнализирует о конвективной неустойчивости среды.

Ненулевое решение системы (39)-(40) существует, если ш и к удовлетворяют биквадратному дисперсионному уравнению, имеющему два корня:

= ^^ + к2 Б + П ± У (к2 V2 - к2^2 + П2)2 + 4к^

(41)

Показано [12, 16], что ш+ отвечает изгибной (альфвеновской), а ш- — медленной волне в присутствии гравитации. Квадрат частоты медленной волны ш- всегда меньше, чем квадрат частоты альфвеновской волны, поэтому ш- первой становится отрицательной: медленная волна первой теряет устойчивость при изменении параметров системы. Подкоренное выражение всегда положительно, следовательно, колебательных неустойчивостей нет.

Условием устойчивости системы, таким образом, является условие устойчивости медленной волны ш- ^ 0. Подставив выражение для ш- и использовав положительность дискриминанта, получим условие устойчивости

к ^ тах

П2

92

V2 + Б2' а2(а2 + V2)

П!

V2

(42)

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

5 1 г к» к» «-» _

. Линейный анализ устойчивости, схема с перешагиванием

Решим задачу о малых колебаниях системы, описанной в предыдущем разделе, численно, используя схему с перешагиванием на однородной сетке (шаг по массовой переменной равен М). Система (29) заменится дискретными уравнениями (6), (7), (24) в

виде

п+1/2 . Л п+1/2 . п+1/2 га-1/2 л п-1/2 п-1/2 ' + + - - - V,-

'¿+1

4-1

'¿+1

6т

1 -

Рвх4(ГП+1/2 )

рп р ¿+1/2

е(гп+х/2)

2

+

+

1

реж4(гП-1/2) рп

р г—1/2

ё(Гп-1/2) + 2ф Гнп лп _ нп 1

2 + М 1-Лг+1/21г+1/2 н г-1/21г-1/2 ]

(43)

Гга+1 _ гп

п+1/2 п п

V , Нг+1/2аг+1/2 = ф

РП+1/2 / (РП+1/2)Т = ССП81, рП+1/2 + (НП+1/2)2 = Ре*^^ ,

п

1г+1/2

г+1

п

Ьг+1/2

Гп = |Гп Гп| Ьг+1/2 = |Гг+1 - Гг

Напомним, что г = 0,...,Ж - 1 и все переменные периодичны по г с периодом N. Линеаризованные уравнения баланса аналогично (35), (38) сведутся к

г 77 г

"н + — = 0, -р = а2"р, -р + 2И-И = ре^Г) 9(У )-У, Н а

-р -а -Ь — + — + -г = 0,

р а Ь

(44)

2

где индекс г + 1/2 опущен для краткости, У = (уг+уг+1)/2, выражение (45) — уравнение непрерывности. Как и в аналитическом случае, (44)-(45) достаточно для определения возмущений плотности и напряженности поля из возмущений координат:

5Н = а2 = а2 + V2

т + а2гу

5р р

V2 ьь

+

9

а2 + V2 Ь а2 + V2

5У.

Для получения волновых уравнений потребуется линеаризовать отдельные слагаемые из правой части уравнения движения (43) и получить вспомогательные выражения:

51

г+1/2

=5

гг+1 - г1

ь

5гг+1 - 5гг - (I, 5гг+1 - 5гг)

ь

5Нг

г+1/2

ех Н

5ь

г+1/2

ё(уг+1/2)

1

а2 + V2 Ь

рех4 (У+1/2)

+

+ V2

5Уг

г+1/2

+ еу Н

5Уг+1 - 5Уг Ь

-еуП25Уг+1/2 - еу9

V2 5Ьг+1/2

рг+1/2

5Ьг+1/2 = 5хг+1 - 5-г, 5Уг+1/2 =

а2 + V2 Ь

5уг+1 + 5уг 2 •

(46)

(47)

(48)

Как и в предыдущем разделе, подстановка этих выражений и переход к представлению возмущений в виде дискретных фурье-мод

5гп = 5г7гавгл', 5уп+1/2 = 5у7 гаегв7

завершает получение из линеаризованного уравнения (43) системы волновых уравнений:

2 + сов(в) 3

х5х = 252С08(в) - 1• г 9^ ^(Я

ь2

-5- + г

а2 + V2 Ь

5У,

(49)

2 + ео8(в) х5у = 2^ ео8(в) - 15у - ^сов(в) + 1 5у г 9V2 в1п(в)

3

Ь2

2

15У - г

а2 + V2 Ь

5-,

1 / 7 - 1

X = I • (50)

Параметр х аппроксимирует -(см. (39) - (40)). Условие существования нетривиального решения этой системы сводится к дисперсионному уравнению:

X2 + ЬХ + с = 0,

где

X =

2 + сов(в),

"X,

Ь = -2(V2 + 52)со8(в2 - 1 +П2 сов(в) + 1

Ь2

2

с = —

9V2

а2 + V2

22

ап2 (в) Ь2

+ 4V2 52

(сов(в) - 1)^ _ 52П2 Сов2 (в) - 1

Ь4

Ь2

(51)

(52)

(53)

(54)

2

а

9

2

а

5

2

Уравнение (51) снова имеет два вещественных корня (изгибная и медленная волны):

3(2 + сс8(в))х± = 2(У2 + £2)С°8(в2 1 - П2

Ь2

ссв(в) + 1 2

±

± ,, , 2(^2 - £2)С°5М-1 - П2с^вШ)2 + 4!^ ((55)

Ь2

2

Ь2

а2 + V2

Найдем интервал устойчивых т, для которых |7| ^ 1 при любых в• Из (50) следует

1 , Хг , •

7 = ± г

1

1 +

Хт 2

22 = 1 + хт! ± 2

1+

Хт 2

- 1,

(56)

откуда очевидно, что условием устойчивости является выполнение неравенства

^2

1+

хт

2

^ 1.

(57)

Действительно, при его выполнении существует такое ф £ Л, что 1 + хт2/2 = ссв(ф), а 7 = е±г^ • При нарушении условия (57) одно из решений (56) по модулю больше 1, а другое меньше, отвечая экспоненциально затухающему и экспоненциально растущему решениям.

Из (57) следует, что знак х±, как и знак ранее (41), определяется физической устойчивостью трубки. Для устойчивых систем х± ^ 0 и (57) играет роль критерия Куранта—Леви явной схемы с перешагиванием, ограничивая т сверху.

Покажем, что переход к дискретной системе не внес новых неустойчивостей. Для этого выпишем решение квадратного уравнения (51) в виде х = (—Ь ± VЬ2 — 4с)/2 (дискриминант всегда неотрицателен, см. (55)). Система физически устойчива (х± ^ 0), если Ь ^ 0 и с ^ 0 для любого в. Необходимое условие устойчивости системы, таким образом, принимает вид

—2(У 2 + £ 2)сс8(в 2 — 1 + П С°«(в ) + 1

Ь2

2

> 0,

^ 2

а2 + V2

ап2 (в) Ь2

+ 4V 2£2

(ссв(в) — 1)2 , «1П2 (в)

Ь4

+ £ 2П2

Ь2

> 0.

(58)

(59)

Первое условие аппроксимирует условие к2 ^ — П2/^2 + £2) и всегда выполняется для конвективно устойчивой среды. Для второго условия при поиске интервалов знакопостоянства удобно избавиться от тригонометрических функций заменой

ссв(в) = г, вш2(в) = 1 — г2, г £ [—1,1].

(60)

Сделав замену, сократив на 1 — г и исключив корень г =1 (в = 0) из дальнейшего рассмотрения, приходим к неравенству

£ 2П2-

^2

+ V2

(1 + г) + 4V2 £2

1 — г

> 0, г £ [—1,1),

или

2 + ь2^) ^ 4V2 — ь2^, ^ =

£2 V2

а2 (а2 + V2)

— П2, г £ [—1,1).

2

2

2

2

а

Если Г ^ 0, то система неравенств верна для всех Ь (схема безусловно устойчива). Если же Г > 0, то итоговое неравенство принимает вид

4V2 - Ь2Г 4V2 + Ь2Г •

(62)

Важным моментом при дальнейшем рассмотрении этого неравенства является то, что фактор в дискретен. Максимально возможное значение Ь, лежащее в интервале [-1,1), равно соъ(2п/И). Подставив это значение, получим необходимое условие устойчивости:

, , , 4V2 - Ь2Г , Л

^^ ^ 4^2 + Ь2Г • (63)

Это аппроксимация второго условия физической устойчивости системы (42). Действительно, для больших N (малых Ь) разложение в ряд дает

1 к2Ь2 4V2 - Ь2Г 1 ГЬ2 Ь = сов(в) — 1--1—, -тт^-— 1 -

2

приводит (63) к виду

к2 >

4V2 + Ь2Г

П2

2У 2

а2(а2 + V2) V2'

Таким образом получаем, что для физически устойчивой системы схема с перешагиванием устойчива при выполнении условия Куранта (57). На практике вместо (57) удобнее использовать мажоранту. Так как = (-Ь - VЬ2 - 4с)/2 и Ь, с > 0, имеем оценку снизу х- ^ -Ь. Из нее следует (57) более жесткое ограничение сверху на т2:

т

т

3(2 + ео8(в))

-1

-2(V2 + 52)С08(в2 - 1 + Псо«(в) + 1

Ь2

2

> - 2.

(64)

После использования замены (60) и перегруппировки оно сводится к неравенству

т

V2 + 52 Ь - 1 П2 Ь +1

— 2--Ч---

Ь2 Ь + 2 2 Ь + 2

о

^ 8' Ь е [-1' 1].

(65)

Экстремум по Ь достигается на краях интервала, откуда следует удобная мажоранта:

(66)

2 V2 + 52 . 2

т

Ь2 « 3■ т2п2 « 8

Оба условия очень прозрачны, так как от явной схемы требуется разрешать частоту как гравитационных колебаний, так и магнитоакустических волн, бегущих вдоль трубки.

2

д

6 1 г «_» *.» *.» _

. Линеиныи анализ устойчивости, симплектическая схема

В симплектической схеме (25) для определения обобщенной силы используется тот же самый гамильтониан, что и в схеме с перешагиванием, поэтому линеаризация и переход от малых возмущений к их фурье-модам совпадает с проделанным в предыдущем

разделе. В дальнейшем нам потребуется только матрица Е для определения возмущения силы как функции возмущения координат, выписанная покомпонентно в правой части (49):

6//х

6 и у

2 СОв(в) - 1

дУ2 яп(в)

—г

а2 + У2 Ь

. дУ2 эш(в) 1 а2 + У2 Ь 2У2 СОй(в) - 1 ^2_СОй(в) + 1

6х 6у

Е 6т.

Ь2 2

Здесь V = р/М, V = V. Коэффициенты схемы (25) удовлетворяют тождествам а1 + а2 + аз + а4 = 2(а1 + а2) = 1, Ь\ + 62 + 63 + 64 = 262 + 63 = 1,

(67)

следовательно, матрицу полного шага по времени можно выписать, используя только два параметра а1 и Ь2, что позволяет получить следующее уравнение на фурье-амплитуды:

7

6т

¿V

2Ы Р(62)2(1/2 - а1 )Р(1 - 262)2(1/2 - ах)Т(кЬ2)^{а1)

6т 6v

А

6т

¿V

(68)

Здесь 2 и Р — матрицы подшага по т и V, действующие на вектор (6гх, 6гу, 6их, 6иу). В блочном виде их можно представить (25) как

2(а)

I тац1 0 I

Р (6)

I 0 т6Е I

I

10 01

3

2 + еов(в)'

Множитель п для перевода ¿V в ¿V получен фурье-преобразованием уравнения (8).

Матрица системы (68) имеет размер 4 х 4 и содержит много степеней т. Подобно матрицам 2 и Р, ее компактнее записать как матрицу 2 х 2 из четырех блоков:

Ац = 1 + ПГгЕ + ^62(2а1 - 1)(262 - 2а1 - 1)Е2 - ^а^^ - 1)2(262 - 1)Е

з

П2т 3

п3 т 5

А12 = пт + -4-(1 - 262 + 8а162 - 8а162)е - а162(2а1 - 1)(2а^2 - 262 + 1)е2-

4

4 7

2

4

а162(2а1 - 1)2(262 - 1)Е3

12 _з

т12т 5

А21 = те + П^62(2а1 - 1)(62 - 1)е2 - 62(2а1 - 1)2(262 - 1)е3 А22 = А

11.

Исследование устойчивости схемы начнем с того, что (68) является уравнением на собственные значения матрицы А, но прямое нахождение ее собственных значений нецелесообразно: необходимо найти не столько корни, сколько интервал значений т, на котором |7| ^ 1 для всех в. При поиске самих собственных значений объем работы уже слишком велик1, не говоря о дополнительном исследовании полученных выражений на максимум.

коэффициенты характеристического полинома являются полиномами 6-й и 12-й степеней т и формулы Кардана будут исключительно громоздкими.

Ряд свойств симплектических матриц позволяет заметно упростить задачу и перейти непосредственно к получению условий на т. Матрица А удовлетворяет условию

Ат ЗА = 3, 3

0 —I

1 0

которое проверяется представлением А в виде произведения (68) и использованием тождеств От 3 <2 = Рт 3 Р = 3. Это значит, что матрица А является симплектической, а условие устойчивости схемы |71 ^ 1 совпадает с определением устойчивости симплети-ческого преобразования А. Показано [17]2, что характеристический полином симплек-тической матрицы возвратный, т. е. имеет вид

'А4 + аА3 + ЪА2 + аА +1 = 0,

а = —А1 — А2 — А3 — А4, (69)

Ъ = А1А2 + А1А3 + А1А4 + А2А3 + А2А4 + А3А4.

Подставив матрицу е в общем виде и использовав е21 = е12 (черта обозначает комплексное сопряжение), получим

а = —4 — Г2 {ец + е22} — г4 Т {е21 + е22 — 2 е?2} — Г6 Ш + Е32 — 3 Е22 (Еп + е22)}, Ъ = 6 + 2 Г2 {е11 + е22} + Г4 {2 Т (е21 + е22 — 2 е22) + е22 + е22}+ + Г6 {Т (е11 + е22) (е11 е22 + е?2) + 2 Ш (е31 + е22 — 3е22 (е11 + е22))}+ + Г8 {Ш (е?1 + е22 — 2 е22) (е11 е22 + е22) + Т2 (е21 е22 + 2 е11 е22 + е12)}+ + Г10 ТШ (е11 + е22) {е21 е22 + 2 е11 е22 е22 + е42}+

+ г12 ш 2 {е31 е32 + + з е22 + з Е22 е?2},

где

22 г = Пт ,

Т = Ъ2(2а1 — 1)(22Ъ2 — 2а1 — 1) * 8.3333333333 ■ 10-02,

Ш = а1Ъ2(2а1 — 1)2(2Ъ2 — 1) * 1.2950839901 ■ 10"01.

2

Из возвратности характеристического полинома следует, что для любого корня А величина 1/А тоже является корнем. Вещественность коэффициентов а и Ъ означает, что все корни образуют комплексно-сопряженные пары. Многочлен четвертой степени имеет четыре корня, поэтому вследствие указанных симметрий любой набор собственных чисел принадлежит к одной из следующих групп (А, В,ф,ф — вещественные числа):

(70)

(71)

(72)

А1,2 = А±1е±гф, Аз,4 = Ат1в±гф. (73)

А1,2 = е±гф, А3,4 = е±гф,

А1,2 = А±1, А3,4 = е±г'ф,

А1,2 = А±1, А3,4 = В±1,

= А±1 е±гф, А3,4 = Ат1е±гф

2В указанной работе, в § 42, введены определения сильной устойчивости и описаны симметрии и

свойства спектра симплектических матриц, используемые далее в разделе.

Покажем, что при т ^ 0 матрица А является сильно устойчивой, т. е. имеется четыре различных корня (69), расположенных на комплексной окружности (решение (70) с ф = 0). Очевидно, что Л = 0 не является собственным значением, поэтому можно ввести и использовать две вспомогательные величины:

С+ = (Л1 + Л2)/2 = (Л1 + 1/Л1)/2, С- = (Лз + Л4)/2 = (Лз + 1/Лз)/2.

Подставив Л2 = 1 /Л1, Л4 = 1/Лз в (69), получим систему уравнений

2((+ + С-) = -а, 2 + 4С+С- = Ь,

решением которой являются

а

4

а

4

2 1 Ь

С± = -т±\/7 + 7 ^. (74)

24

Подставим коэффициенты а и Ь и разложим (74) в ряд по г2:

С± - 1 - ^ (-Е11 - Е22 ± \/(ел - ^22)2 - .

Если оба решения лежат на единичной окружности (70), то = оов(ф), = оо8(0). Так как е11, е22, ^ 0 и Г2 > 0, оба решения £± смогут попасть в интервал [-1,1] исключительно, если т достаточно мало и

- е22 ^

^л - е22)2 - 4е22

11 - Е22 ^ у (Е11 - Е22) - 4^12-

Подставив (67) — выражения для компонент матрицы Е — в это неравенство, получим уже рассмотренное условие физической устойчивости дискретной системы (59). Это значит, что для физически устойчивой системы при малых т матрица А сильно устойчива.

Выполнение условия сильной устойчивости означает, что при непрерывном изменении параметров системы (в нашем случае т) матрица остается сильно устойчивой, пока не произойдет потеря устойчивости, причем только при встрече решений (70) на единичной окружности [17].

Если подставить полные выражения для коэффициентов а и Ь и выделить общие множители, то £± можно представить в компактном виде: г2 г4 Г6

(± = 1 + Т(Ец + Е22) + 4Т(Е?1 + Е22 - 2Е22) - тж(Е31 + Е22 - 3Е22 (Е11 + Е22)) ± ±^ \/{(Е11 - Е22)2 - 4Е22>{1 + Г2Т(Е11 + Е22) - (Е21 + Е22 + ЕЦЕ22 - Е22)}2- (75)

Как было упомянуто [17], условие столкновения корней — необходимое, но не достаточное для схода решений (70) с единичной окружности. В нашем случае подкоренное выражение в (75) всегда неотрицательно, поэтому его нули кратные и смены знака в них не происходит. Вследствие этого £± после встречи остаются вещественными и независимые решения (70) проходят друг сквозь друга3 без потери устойчивости.

3Это подтверждает интуитивно ожидаемую потерю устойчивости сначала одного из решений £±

(при нарушении условия Куранта только высокочастотной волной), а потом оставшегося. Аналогично

рассмотренному в предыдущей секции случаю, пара решений е±гф = е±гшт при этом перейдет в А±1 =

„±Тт

Таким образом, решение становится неустойчивым только при столкновении корня со своим комплексным сопряжением на единичной окружности, т.е. при £ = ±1. Минимальное значение т, при котором |£+| или |£-| станет больше единицы, ограничивает сверху область устойчивости схемы.

В дальнейших выкладках подкоренное выражение будет часто приводиться к тем или иным полным квадратам, т. е. представляться в виде VX2. Знак ± стоит исключительно перед корнем, благодаря чему обе ветви решения можно представить аналитическими функциями без введения модулей, если переписать ±V/X2 как ±X. Это эквивалентно исключительно перестановке ^ в области X < 0 и упрощает работу без потери общности. С учетом вышесказанного вынесем полный квадрат из-под корня и перепишем (75) как

|2 ~4 ~6

I ™ Ч I ^ ~ Ч 1

Z± = 1 + 4 (f11 + f22) + - T (fÎ! + f22 - 2 f22) - - W (f3i + f22 - 3 f22 (f11 + f22))±

±^ [1 +12 T (f11 + f22) -14 W (f21 + F22 + f11f22 - F22)^(f11 - f22)2 - 4 f22.

Так как размерность элементов матрицы f равна квадрату частоты, удобно ввести безразмерные параметры:

X = -I2f11, Y = -I2f22, Z = -I4f22,

и переписать Z± как

1 T W

Z± = 1 - 4 (X + Y) + - (X2 + Y2 + 2 Z) + —(X3 + Y3 + 3 Z (X + Y))±

± 1 [1 -т(X + у)- w(X2 + У2 + Xу + z)] \fiX-У^Т^.

Это позволит найти область устойчивости в пространстве параметров X, У, Z и затем определить максимально допустимое значение т проверкой на выход X, У, Z из нее. Область допустимых параметров определяется системой неравенств

-1 ^ (± ^ 1, X ^ 0, У ^ 0, Z ^ 0, X + У ^ у/(X - У)2 + 4 Z. (76)

Покажем, что область допустимых параметров имеет вид

0 ^ X ^ А, 0 ^ У ^ А, 0 ^ Z ^ шт(XУ, (А - X)(Д - У)), (77)

где А будет определено ниже. На нижней границе ^ = 0) выражение для £± значительно упрощается после выноса полного квадрата (X - У )2 из-под корня

1 _, Т,_0 W

Z± Iz=0 = 1 - 4 (X + Y) + 4 (X2 + Y2) + — (X3 + Y3)± 1

4

что дает

[1 - T (X + Y) - W (X2 + Y2 + XY))] (X - Y),

I X T X2 W X3 I Y T Y2 W Y3

C- |z=o= 1 - y + 2 1 2 , |z=o= 1 - y + — + — • (78)

Входящий в (78) полином обладает следующими свойствами:

WX (X - А) (X + Х)

Р (X) А X

X ТХ2 WX3

1--+-+-

2 2 2

УТ2Т4"Ж - Т

= 1 + 1821/3(2 - 21/3)

2

2 W (1 + 21/3)2(1 - 21/3 + 22/3)

т +

~ 2.47559368819,

2 W

21/3А ~ 3.11905259874,

¡Т 2 I з т

тт Р(X) = Р ( --+-— 1 ~ 0.5593, тах Р(X) = 1.

хе[о,Д] \ 3 W ' Xе[о,Д] 1 ;

(79)

(80) (81) (82)

Из (82) и (78) следует, что для X и У из (77) на нижней поверхности ^ = 0) удовлетворены все неравенства (76). Верхняя граница разрешенной области определяется уравнениями = 1. Решением уравнения £+ = 1 является Z = XУ, что легко

доказывается подстановкой

4.

1 Т о W

= XV ) = 1 - 4 (X + у ) + — (X + у )2 + — (X + у )3 ±

1

± - [1 - Т (X + У) - W (X + У(X + У )2,

4

что с учетом положительности X + У дает

1 Т о W

С+ = 1, С- = 1 - ¿(X + У) + 2 (X + У)2 + — (X + У)3 = Р(X + У).

Выражение для принимает вид исследованного полинома (79), откуда следует, что при 0 С X + У С А область устойчивости ограничена сверху поверхностью Z = XV. Решением уравнения = 1 является Z = (А - X)(А - У). На этой поверхности

С- = 1+ -А + Т А + А3 = Р (А) = 1,

1 Т о W

С+ = 1 + 2 (А - X - У) + 2 (А - X - У)2 - — (А - X - У)3 = Р(X + У - А).

При А С X + У С 2А величина X + У - А € [0, А], откуда (82) следует, что эта поверхность и есть верхняя граница разрешенной области для А С X + У С 2А. Таким образом, область устойчивости схемы определяется условиями

X = 6

С А,

У=6

т2Б21 - сое (в) 2 + сов(в) т2У2 1 - соэ(в) 3 2П2 1 + соэ(в) Ь2 2 + сов (в) + 2 Т 2 + сов(в)

С А,

Z = 9 Т,

г4

ЬЬ2

9V2

а2 + V2

1 - соэ2(в) (2 + со8(в ))2

С min(XУ, (А - X)(А - У)).

4Заметим, что условие физической устойчивости системы -¥ц - ¥22 ^ л/С^ГТ-1^записывается как Z С ХУ.

2

Так как V2 >82, второе неравенство этой системы сильнее первого. Его исследование полностью аналогично уже проведенному (64) и в итоге дает

Т 2 V 2

т2П2 ^ А' 12 —— ^ Д.

Ь2

Неравенство Z ^ ХУ уже исследовано в разд. 5 (см. (59)-(63)). Оно сводится к

д2 V2 4V2 1 - Ь

д < П2 + ^^' Ь, = соз(2п/И).

а2 (а2 + V2) " Ь2 1 + и

Неравенство Z ^ (А - X)(Д - У) не дает удобных мажорант и сильно ограничивает сверху параметр Z при X ~ А' У ~ А. Вследствие этого практичнее ограничиться областью X + У ^ А, что дает условие устойчивости в виде

т2П2 , А. , А (83)

Это примерно 56 % от предельно допустимого временного шага (66) для схемы с перешагиванием.

7. Проверка консервативных свойств схем

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

Т = Т0

вт(аг)

д

до-

сов(аг)

Р = рТ = Ро

вт(аг) в1п(аг0)

30 81п(аг0) аТ0 сов(аго )

-о

вт(аг0) сов(аг0)

Ускорение свободного падения направлено к оси цилиндра, профили температуры, давления, плотности и д(г) показаны на рис. 2. В положении равновесия трубка замкнута в кольцо радиуса 3 и для проверки консервативных свойств элементам плазмы

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

Рис. 3. Исходная форма трубки и максимальные отклонения (отвечают пикам на рис. 4). Единицы времени и расстояния на всех рисунках одинаковы

сообщается начальная скорость (ф) = 0.7соя(3ф), т.е. возбуждается в основном третья мода быстроосциллирующей изгибной волны. На рис. 3 показана форма трубки в нескольких крайних положениях5.

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

Амплитуда колебаний полной энергии для схемы с перешагиванием (¿ьр) и симплектической схемы (дзут)

N = 50 N = 150

т дьр двут т дьр двут

0.300 0.150 0.075 0.010 4.660 ■ 10-3 1.436 ■ 10-3 5.459 ■ 10-4 6.397 ■ 10-5 1.707 ■ 10-4 (т = 0.29) 1.180 ■ 10-5 7.313 ■ 10-7 2.301 ■ 10-10 0.100 0.050 0.025 0.003333 5.317 ■ 10-4 1.673 ■ 10-4 6.409 ■ 10-5 7.177 ■ 10-6 2.164 ■ 10-6 (т = 0.095) 1.655 ■ 10-7 1.034 ■ 10-8 3.197 ■ 10-12

Примечание. Единицы абсолютные, в этих единицах интервал изменения гравитационной энергии имеет ширину 0.515, тепловой энергии — 0.379, кинетической — 0.232, магнитного поля — 0.162 и работы на вытеснение внешней среды — 0.292 (см. рис. 4). Для симплектической схемы пришлось уменьшить максимальное значение т до показанных в скобках из-за условия устойчивости Куранта (83). Выделенная ячейка отвечает данным рис. 4.

Рис. 5. Погрешность закона сохранения энергии на больших временах: кривые 1 и 2 отвечают схеме с перешагиванием для N = 60, т = 0.05 (1) и N = 600, т = 0.005 (2); кривые 3 и 4 отвечают симплектической схеме с N = 60, т = 0.05 (3) и N = 600, т = 0.005 (4)

Рис. 4. Распределение энергии по степеням свободы для выделенного в таблице численного эксперимента: 1 — тепловая энергия, 3 — гравитационная, 4 — энергия магнитного поля, 5 — кинетическая энергия (отложены по левой оси); 2 — отклонение полной энергии от начальной (отложена по правой оси); форма трубки для нескольких пиков кривой 2 показана на рис. 3

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

временах (десятки тысяч периодов колебаний) с целью продемонстрировать отсутствие ярко выраженных долговременных трендов для обеих схем в течение большого интервала времени.

8. Заключение

В данной работе сформулирован новый метод построения конечно-разностных схем для исследования эволюции тонкого магнитного волокна в незамагниченной плазме. Преимущества метода — физическая чистота и прямой переход от модельного приближения к схеме: в большинстве работ сначала предполагается некоторое модельное приближение (однородность параметров плазмы поперек оси [1] или фиксированный поперечный профиль, характеризуемый несколькими параметрами [3, обзор лит-ры]), потом производится усреднение уравнений МГД поперек сечения трубки с последующими попытками получить систему одномерных уравнений. После получения система еще раз аппроксимируется для решения численными методами.

При использовании разработанного формализма конечно-разностная схема выписывается непосредственно как уравнения движения дискретной системы, при этом гарантируется совместность определения всех вкладов в полную энергию и другие инварианты системы. Метод позволяет использовать более сложные приближения для описания формы трубки или ее поперечной структуры — все дополнительные параметры войдут как лагранжевы переменные на равных правах с г и V. Для примера можно аппроксимировать ось трубки безье-кривой или сплайном, либо предположить гауссов профиль напряженности магнитного поля поперек оси, либо учесть линейный и квадратичный вклады в распределение внешних параметров по высоте — построение схемы в любом случае сведется к определению соответствующих вкладов в (1) и нахождению обобщенной силы без промежуточных выкладок по усреднению МГД-уравнений в поперечном сечении трубки, часто нуждающихся в дополнительном обосновании.

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

Отметим, что выбор схемы определяется решаемой задачей и часто слепые усилия, направленные на достижение точного сохранения энергии (или любого другого выбранного наобум инварианта), приводят к полному искажению качественного поведения системы на больших временах [9, примеры из обзора лит-ры]. В частности, несмотря на то, что симплектические схемы не могут одновременно сохранять с разностной точностью симплектический инвариант и полную энергию системы [18], возможно сохранение энергии с точностью примерно 10-10 % на больших временах без чрезмерных численных затрат (рис. 5, кривая 4); при этом сохранение симплектического инварианта не возмущает структуру течения в фазовой плоскости и его сепаратрис [15].

Обобщение используемой дискретной физической модели для случая адаптивных сеток очевидно, хотя при этом симплектические схемы теряют ряд своих преимуществ (разработка симплектических схем с адаптивными шагами — все еще открытая задача [9]).

В работе были рассмотрены явные схемы только для экономии ресурсов при численном решении нелинейной системы (6)-(7)6. Вынужденное многократное решение этой системы при обращении неявной системы уравнений каким-либо итерационным методом могло свести на нет преимущества большего временного шага; тем не менее это не исключает использования неявных схем для решения системы (22)-(23), если это будет выгодно для конкретных задач.

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

[1] Романов Д.В., Романов К.В. Численное моделирование развития неустойчивости медленной волны тонкой магнитной трубки в конвективной зоне Солнца // Вычисл. технологии. 2001. Т. 6, № 6. С. 81-92.

[2] Fan Y., Fisher G.H., DeLuca E.E. The origin of morphological asymmetries in bipolar active regions // Astrophys. J. 1993. Vol. 405. P. 390-401.

[3] Романов Д.В. Математическое моделирование влияния многомерности на эволюцию магнитных полей и структуру аномального прогрева солнечной атмосферы: Дис. ... канд. физ.-мат. наук. Красноярск, 2003. 129 с.

[4] Parker E.N. Stellar fibril magnetic system. I. Reduced energy state // Astrophys. J. 1984. Vol. 283. P. 343-348.

[5] Abbett W.P., Fisher G.H. The three-dimensional evolution of rising, twisted magnetic flux tubes in a gravitationally stratified model convection zone // Astrophys. J. 2000. Vol. 540. P. 548-562.

[6] Abbett W.P., Fisher G.H. The effects of rotation on the evolution of rising omega loops in a stratified model convection zone // Astrophys. J. 2001. Vol. 546. P. 1194-1203.

[7] Самарский А.А., Попов Ю.П. Разностные схемы газовой динамики. М.: Наука. 1975. 351 c.

[8] Годунов С.К. Разностные аппроксимации переопределенных гиперболических уравнений // XIII Междунар. конф. по методам аэрофиз. исследований. Новосибирск. 2007.

[9] Stuchi T.J. Symplectic Integrators Revisited // Braz. J. Phys. 2002. Vol. 32, N 4. P. 958-979.

[10] Yoshida H. Recent progress in the theory and application of symplectic integrators // Celest. Mech. 1993. Vol. 56, May. P. 27-43.

[11] Okumura H., Itoh S.G., Okamoto Y. Explicit symplectic integrators of molecular dynamics algorithms for rigid-body molecules in the canonical, isobaric-isothermal, and related ensembles // J. of Chem. Phys. 2007. Vol. 126. P. 084103.

[12] Alekseenko S.V., Dudnikova G.I., Romanov V.A. et al. Magnetic field instabilities in the solar convective zone // Rus. J. of Eng. Thermophys. 2000. Vol. 10. P. 243-262.

[13] Ландау Л.Д., Лифшиц Е.М. Механика. М.: Наука, 1988. 215 с.

[14] Numerical recipes in C: the art of scientific computing / W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery. Cambridge: Cambridge University Press, 2002.

[15] Candy J., Rozmus W. A symplectic integration algorithm for separable hamiltonian functions // J. of Comp. Phys. 1991. Vol. 92. P. 230-256.

6,

Одна из самых ресурсоемких задач, на решение которой уходит до 80 % ресурсов.

[16] Романов Д.В., Романов К.В., Ничковл Н.М., Шалагина Е.В. Влияние углового вращения на подъем тонкой магнитной трубки из зоны Динамо // Тр. IX Междунар. байкальской молодежной научной школы по фундаментальной физике. Иркутск, 2006. С. 219-221.

[17] Арнольд В.И. Математические методы классической механики. М.: Наука, 1989. 472 с.

[18] Zhong Ge, Marsden J.E. Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators // Phys. Letters A. 1988. Vol. 133, N 3. P. 134-139.

Поступила в редакцию 15 октября 2007 г., в переработанном виде — 1 марта 2008 г.