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

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

Александр Солженицын, «В круге первом».

Актуальность работы.

В последние десятилетия проведено большое количество экспериментов с выбросом в верхнюю атмосферу и магнитосферу плазменных и слабоионизованных сгустков и струй сравнительно высокой массы [12−14]. К наиболее интересным из подобного рода геофизических экспериментов относятся эксперименты с использованием мощных возмущений [15,16]. Научное и прикладное значение подобных экспериментов зависит, прежде всего, от того насколько удаётся понять и численно промоделировать их результаты [17−20]. Как известно, эволюция крупномасштабных возмущений в верхней атмосфере и магнитосфере сопровождается весьма сложными и многообразными физическими процессами, по некоторым из которых пока нет достаточно полного понимания. Кроме того, эволюция образующейся возмущённой области, вообще говоря, имеет пространственный трёхмерный характер, и лишь на отдельных временных стадиях своего развития может приближенно моделироваться геометрией меньшей размерности. Эти два обстоятельства существенно осложняют процесс численного моделирования эволюции плазменных сгустков в ионосфере и магнитосфере, имеющих.

23 6 высокую начальную энергию и массу (Е < 10 эрг, М < 10 г). Проделать это на основе имевшихся ЗО-кодов не представляется возможным, так как они обычно применимы лишь для слабоградиентных течений идеальных сплошных сред [21−24], и не учитывали определяющих внутренних процессов, связанных с неравновесной ионизацией и диффузией магнитного поля.

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

С определённой степенью условности весь процесс динамики плазмы можно подразделить на ряд стадий — в зависимости от содержания и важности физических процессов. В соответствии с этим различаются и физико-математические постановки задач на различных стадиях.

В [4], где с использованием достаточно полной физической модели в одномерной постановке со сферической геометрией рассчитывалась начальная стадия взрыва в верхней атмосфере. Показано (высоты 100−1000 км), что характер поведения параметров плазмы на ранней стадии разлёта и развитие возмущенной области в целом на поздней стадии её торможения и конвекции существенно зависят, прежде всего, от удельного начального энерговыделения. Инерционный разлёт плазмы и начальная стадия её торможения в атмосфере за счёт взаимодействия с сильно разреженным воздухом и геомагнитным полем сопровождается разнообразными кинетическими процессами, определяющими ионизационный состав плазмы и её взаимодействие с геомагнитным полем. Достаточно подробное решение задачи неравновесной ионизационной кинетики на ранней стадии динамики плазмы по существу определяют всё последующее её развитие. Однако корректное решение задачи с учётом ионизационной и температурной неравновесности требует достаточно большого расчетного времени, так как основано на неявных или многоточечных алгоритмах. Это в свою очередь накладывает определенные ограничения на пространственную размерность задачи: использование ЗЭ и даже 20 алгоритмов с учетом всего комплекса физических процессов делает численный алгоритм чрезвычайно трудоемким. Однако, как показали расчёты, выполненные в [4], на ранней стадии разлёта плазмы вполне оправдано использование одномерного лагранжева численного кода, учитывающего всю совокупность физических процессов. По мере торможения плазмы течение приобретает пространственный характер, становится существенным азимутальное перетекание массы, значительную долю в которой начинает составлять воздух.

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

Цели настоящей работы.

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

В данной работе были поставлены следующие задачи:

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

• всестороннее тестирование программного комплекса и его апробация на прикладных задачах;

• разработка удобного интерфейса для представления и визуализации результатов расчетов.

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

Кроме того, большое внимание должно быть уделено объединению различных элементов в работоспособный программный комплекс, допускающий его дальнейшее развитие.

Научная новизна.

Для трехмерных уравнений магнитной газодинамики (МГД — модели с учетом диффузии магнитного поля) разработаны новые модификации монотонных консервативных вариантов сеточно-характеристического метода и метода типа Годунова 2−3-го порядка аппроксимации.

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

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

Научная и практическая ценность.

Построен и реализован в виде комплекса программ эффективный монотонный численный алгоритм решения трёхмерных задач плазмодинамики, обладающий высоким порядком аппроксимации. Алгоритм применён для решения крупномасштабных МГД — задач о поведении плазмы в околоземном космическом пространстве при существенно различных условиях — на высотах от 150 до 1000 км.

Определены основные особенности и закономерности в поведении плазменных областей. В частности, показано, что для 150 км магнитное поле начинает заметно влиять на поведение плазменной области уже при / > 1 с, и с этого же времени характер её поведения становится трёхмерным.

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

Для умеренных энерговыделений реализуется восходящее течение в виде гигантской струи. Для больших широт струя будет направлена вертикально по полю. Для средних широт картина в целом остаётся той же. Формируется крупномасштабная плазменная струя под углом к силовым линиям поля, которая затем разворачивается и происходит широтное растекание плазмы по полю. Таким образом, при умеренных энерговыделениях, хотя плазменная струя и может подниматься на большие высоты, но в целом её движение направляется геомагнитным полем.

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

Апробация работы.

Результаты, полученные в работе, докладывались на конференциях:

1. 53rd Annual Meeting of the Division of Fluid Dynamics, November 19−21 2000, Washington, DC, USA.

2. First MIT Conference on Computational Fluid and Solid Mechanics, June 12 -15 2001, the Massachusetts Institute of Technology, Cambridge, MA, USA.

3. Международный российско-японский симпозиум — Актуальные проблемы вычислительной механики, 5−10 августа 2002, г. С-Петербург, Россия.

4. Second MIT Conference on Computational Fluid and Solid Mechanics, June 1720 2003, the Massachusetts Institute of Technology, Cambridge, MA, USA.

5. 30th EPS Conference on Controlled Fusion and Plasma Physics, July 7−11 2003, St Petersburg, Russia.

6. 18th International Conference on Numerical Simulation of Plasmas, September 7−10 2003, at the Sea Crest Oceanfront Resort and Conference Center in Falmouth, Massachusetts, USA.

7. 31st EPS Conference on Plasma Physics, 28 June — 2 July 2004, Imperial College, London, Great Britain.

8. Third MIT Conference on Computational Fluid and Solid Mechanics, June 1417 2005, the Massachusetts Institute of Technology, Cambridge, MA, USA.

Публикации.

Основные результаты сделанной работы опубликованы в следующих изданиях:

1. О. В. Воробьев, Я. А Холодов. Об одном методе численного интегрирования одномерных задач газовой динамики // Математическое моделирование — 1996 — Т. 8, № 1 — С. 77−92.

2. Y.A. Kholodov, A.S. Kholodov, E.L. Stupitsky, A.Y. Repin. Numerical Simulations of the Plasma Cloud Expansion in the Rarefied Ionosphere with Magnetic Field by Means of 3D MHD Equations // 30th EPS Conference on Controlled Fusion and Plasma Physics — 2003 — St Petersburg — Vol. — 27A — P. 2.210.1−2.210.4.

3. E.L. Stupitsky, A.S. Kholodov, A.Y. Repin, Y.A. Kholodov. Numerical Modeling Plasmadynamics Progress in Large Scale Ionosphere Experiments // 30th EPS Conference on Controlled Fusion and Plasma Physics — 2003 — St Petersburg — V. 27A — P. 2.209.1−2.209.4.

4. E.JI. Ступицкий, А. Ю. Репин, A.C. Холодов, Я. А. Холодов. Численное исследование поведения высокоэнергетичного плазменного сгустка в верхней ионосфере. Часть 1. Начальная стадия разлёта и торможения плазменного сгустка // Математическое моделирование — 2004 — Т. 16, № 7 -С. 43−58.

5. А. С. Холодов, Я. А. Холодов, Е. Л. Ступицкий, А. Ю. Репин. Численное исследование поведения высокоэнергетичного плазменного сгустка в верхней ионосфере. Часть 2. Разработка трёхмерной модели // Математическое моделирование — 2004 — Т. 16, № 8, — С. 3−23.

6. Y.A. Kholodov, A.S. Kholodov, E.L. Stupitsky, A.Y. Repin. Numerical simulation of the convective plasma dynamics stage at the ionosphere motion by means of 3D MHD equations // Computer Physics Communications — 2004 -Vol. 164,№ 1−3,-P. 91−97.

7. E.L. Stupitsky, A.S. Kholodov, A.Y. Repin, Y.A. Kholodov. Numerical modeling of behavior of high-energy plasmoid in upper ionosphere // Computer Physics Communications — 2004 — Vol. 164, № 1−3, — P. 258−261.

8. A.S. Kholodov, Y.A. Kholodov, E.L. Stupitsky, A.Y. Repin. Numerical modeling of behavior of high-energy plasmoid in upper ionosphere and geomagnetic field // 31 st EPS Conference on Plasma Physics — 2004 — LondonVol. 28G-P. 1.069.1−1.069.4.

9. M.O. Vasiliev, A.S. Kholodov, Y.A. Kholodov, E.L. Stupitsky, A.Y. Repin. Numerical researches of the plasma jet stream formation in the large-scale geophysical experiment // 31st EPS Conference on Plasma Physics — 2004 -London — Vol. 28G — P. 1.070.1−1.070.4.

10. A.C. Холодов, Я. А. Холодов, E.JI. Ступицкий, А. Ю. Репин. Численные исследования поведения плазменного облака в верхней ионосфере // Математическое моделирование — 2005 — Т. 17, № 11 — С. 43−62.

11. М. О. Васильев, А. С. Холодов, Я. А. Холодов, E. J1. Ступицкий, А. Ю. Репин. Формирование крупномасштабного струйного течения в результате развития желобковой неустойчивости // Математическое моделирование -2005;Т. 18, № 1-С. 17−28.

Обзор существующих работ.

Из огромного количества существующих книг, рассматривающих вопросы численного решения гиперболических систем уравнений, следует упомянуть две из них. Это книги: «Сеточно-характеристические численные методы» К. М. Магомедова и А. С. Холодова [25] и «Математические вопросы численного решения гиперболических систем уравнений» А. Г. Куликовского, Н. В. Погорелова и А. Ю. Семенова [26]. И если первая из них посвящена в основном самому сеточно-характеристические методу и использованию его применительно к различным практическим задачам, то вторая имеет гораздо более общий прикладной характер. В своей книге авторы А. Г. Куликовский, Н. В. Погорелов и А. Ю. Семенов смогли не только дать наиболее подробное на сегодняшний день описание различных математических аспектов численного решения гиперболических систем уравнений, но и изложить материал в тесной взаимосвязи с такими важными областями применения этих систем, как теория мелкой воды, газовая динамика, магнитная гидродинамика, динамика твердого деформируемого тела и ряд неклассических областей механики сплошной среды. Считаю необходимым заметить, что большинство вопросов, рассматриваемых в данной диссертационной работе, напрямую отталкиваются от материала изложенного в этих двух книгах.

Противопоточные и симметричные TVD схемы в последние время пользуются большой популярностью при расчете сложных газодинамических течений с многочисленными ударными волнами. Это связанно с их численной устойчивостью при расчете сильных ударных волн. Перенос этих методов на уравнения магнитной газовой динамики не очевидно. Точное решение задачи о распаде произвольного МГД разрыва (задачи Римана) [27] слишком многовариантно для использования при проведении численных расчетов. По этой причине на практике обычно используются различные приближенные решения задачи Римана. Важным моментом является то, что вместо нахождения полного решения, достаточно получить только численный поток, как аппроксимацию задачи Римана. При этом, правда, сохраняется проблема соответствия полученного решения точному, в том случае если решение задачи Римана не является единственным.

В работе [28] предложен метод решения МГД уравнений, основанный на разновидности схемы Ошера, предложенной Bell [29]. В дальнейшем метод был усовершенствован [30]. Тесты показали хорошее разрешение ударных волн, однако оказалось, что в окрестности волн разрежения могут появиться малые возмущения.

В работе [31 ] был обобщен на МГД задачи кусочно-параболический метод (РРМ), который первоначально был разработан для уравнений газовой динамики [32]. Как уже говорилось, точное решение задачи Римана слишком сложно для включения его в алгоритм. По этой причине в РРМ используется приближенное решение нелинейной задачи. Упрощение достигается путем аппроксимации волн разрежения ударными волнами разрежения, удовлетворяющими условиям на разрывах и приводящими к уменьшению энтропии. Это позволяет организовать итерационный процесс для нахождения скоростей волн и значений параметров в областях постоянного решения. Так как описанный метод является итерационным, могут возникнуть трудности со сходимостью итераций. В частности в работе [33] имелась ошибка в результатах расчетов, в дальнейшем исправленная в работе [31]. Разработчики метода предлагают оставаться в таких случаях в рамках решения линеаризованной задачи, которая имеет инварианты вдоль соответствующих характеристик и берется в качестве начального приближения для итерационного процесса. Сильной стороной этого метода является возможность точного удовлетворения соотношений на всех разрывах, а также постоянное присутствие альфвеновского разрыва в возможной конфигурации решения.

Другие методы линеаризации исходной задачи Римана предложены в работах [34−40]. В этих работах замороженная якобиева матрица не является функцией единого усредненного вектора переменных, как это имеет место в обычной газовой динамике, но сложным образом зависит от значений переменных на правой и левой сторонах поверхности разрыва.

В работах [41−43] было показано, что такая процедура усреднения не может быть единственной. Полученное в работах многопараметрическое семейство решений линеаризованной задачи Римана позволяет точно выполнять законы сохранения на разрывах, также показано, что для получения физически допустимых решений исходной задачи важным является правильный выбор параметров вектора усреднения.

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

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

Другой метод [45] использует для уничтожения магнитного заряда искусственный скалярный потенциал, удовлетворяющий уравнению Пуассона для напряженности магнитного поля, при этом возникает необходимость в численном интегрировании этого уравнения. На его решение, как показано в [46], часто расходуется до 30% полного расчетного времени. Отметим также, что граничные условия для уравнения Пуассона могут быть довольно сложными. Несмотря на наличие этих недостатков, метод искусственного скалярного потенциала, в настоящее время широко используется при моделировании астрофизических задач [47−49].

Принципиально другой подход используется в следующих работах [50−52]. Основная идея этого подхода состоит в определении консервативной формы МГД уравнений, которая являлась бы следствием её исходной формы в пренебрежении условием бездивергентности магнитного поля. В результате получаемая таким образом система уравнений будет не строго консервативной, но это как раз будет предпочтительно для устранения магнитного заряда, и её якобиева матрица будет иметь в этом случае размерность 8×8. То есть, можно решать одномерную задачу Римана с непостоянной нормальной компонентой магнитного поля и использовать её для получения решения многомерной задачи. Следует отметить, что такая процедура приводит к появлению дополнительных нефизических волн в решении задачи. И этот подход, по-видимому, хорошо работает для задач с открытыми внешними границами [5354]. Однако сомнительно, что он также может быть использован в нестационарных задачах. Это сомнение вызвано тем фактом, что методы, основанные на расширенной МГД системе, удовлетворяют условию соленоидальности магнитного поля только в бесконечном пределе по времени, когда устанавливается стационарное решение задачи.

Другой важный момент заключается в том, что для точной реализации условия соленоидальности магнитного поля его надо знать скорее на границах расчетных сеток, чем в их центрах. Это стимулирует применение смещенных расчетных сеток [55−58], в которых магнитное поле задается на гранях ячеек, а остальные величины — в их центрах. У этого подхода есть только один недостаток, обусловленный тем, что для нахождения величин в центре ячеек, при численном интегрировании, также понадобятся значения магнитного поля в центре этих ячеек. Это можно сделать с помощью интерполяции, но при этом точность вычислений уменьшится, особенно если магнитное поле будет разрывным.

ЗАКЛЮЧЕНИЕ

.

Для 30 уравнений магнитной газодинамики (МГД — модели с учетом диффузии магнитного поля) разработаны новые модификации монотонных консервативных вариантов сеточно-характеристического метода и метода типа Годунова 2−3-го порядка аппроксимации.

Общий алгоритм перехода от известного состояния при t = tn к искомому состоянию в момент времени / = /" +1 включал расщепление по пространственным переменным, а, при наличии разрывов большой интенсивности, также по физическим процессам («газодинамический» этап, «магнитный» этап, этап расчета диффузии магнитного поля).

В результате двух и трёхмерных расчётов на основе разработанных численных методов и теоретического анализа выяснена картина поведения плазменного течения взрывного типа в околоземном космическом пространстве для времен порядка 100 секунд. С определённой степенью условности в развитии течения можно выделить несколько стадий, продолжительность и структура которых зависит от трёх основных параметров взрыва (Е0,И,(р):

• начальная стадия разлёта плазмы, когда формируется её динамические и ионизационные характеристики;

• стадия торможения плазмы и формирование тепловой возмущённой области и волнового МГД-возмущения;

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

• глобальное дрейфово-диффузионное движение плазмы и образование искусственного радиационного пояса.

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