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

Краткое содержание работы

Введение

Глава I. Физическая и математическая постановка задачи

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

 для плазмы и фона рассматривалась сплошная односкоростная плазменная среда;

 при плазма и фон разделены не только заданием различных термодинамических характеристик, но и степенью ионизации , где и – концентрации электронов и тяжелых частиц;

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

С учётом сделанных предположений система МГД уравнений, описывающая развитие возмущённой области, начиная со стадии интенсивного торможения плазмы, имеет вид:

здесь: , – ускорение силы тяжести, – единичный тензор. Диффузия поля учитывалась в квазиоднородном приближении:

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

Выражение для коэффициента диффузии приведено в [4]. Остальные обозначения в (1.1)-(1.6) общепринятые.

Начальные параметры плазмы задавались исходя из расчётов, выполненных в [4], однако сам момент варьировался в процессе численного эксперимента. Степень ионизации фоновой среды предполагалась пространственно однородной, примерно соответствующей средней степени ионизации воздуха под действием жёсткого излучения плазмы, выходящего из неё в самые первые моменты времени.

Для изотермической однокомпонентной атмосферы имеем распределение Больцмана: . Формула отличается от стандартной, так как учтена нелинейность гравитационного потенциала и, соответственно, изменение g с высотой. Соответственно для n-компонентной изотермической атмосферы аналогичным образом получается:

где - плотность j-й компоненты “на бесконечности” – находится из соотношения по известным экспериментальным значениям (, T), - молярная масса j-го вещества, g₀ - значение g на расстоянии r₀, r₀ - расстояние до центра Земли от “начального уровня” высоты, начиная с которой мы считаем атмосферу изотермической, T - температура атмосферы, R - универсальная газовая постоянная.

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



компоненты поля равны: , . Здесь 810²⁵ Гссм³ – магнитный дипольный момент Земли, – расстояние от центра Земли, – магнитная широта, причём , где – угол между местом наблюдения и магнитной осью Земли.

Глава II. Построение высокоточных монотонных разностных схем для гиперболических систем уравнений

В данной главе с использованием простейшего уравнения переноса

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

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

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



Которая может быть представлена в следующем виде:

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

Аналоги разностных схем для уравнения переноса (2.1), сохраняющие все их свойства (порядок аппроксимации, устойчивость, монотонность и т.д.), имеют вид:



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

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

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

представлялась в виде суммы трех одномерных операторов (заключены в скобки):

Каждое из слагаемых в скобках уравнения (2.7)

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

Глава III. Обобщение и тестирование полученных высокоточных численных методов на случай уравнений газовой динамики

Здесь –плотность, – три компоненты вектора скорости движения газа, – удельная внутренняя энергия, – полная энергия единицы объёма и – давление. Замыкает систему уравнений (3.1) уравнение состояния (УРС) идеального совершенного газа

где – показатель адиабаты. Матрица Якоби системы уравнений (3.1) может быть легко найдена при помощи тождества . Таким образом, семейство схем в виде (2.5) может быть легко построено для системы уравнений (3.1). Предиктор в схеме (2.5) в случае присутствия разрывов большой интенсивности (на два и более порядка) должен быть взят из схемы Годунова, в случае отсутствия сильных разрывов в решении поставленной задачи, это может быть консервативный вариант сеточно-характеристического метода.

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

Нестационарная система уравнений Эйлера газовой динамики в трехмерной постановке в дивергентной форме имеет следующий вид:



Консервативные уравнения (3.4) описывают законы сохранения массы газа, его импульса и энергии. Замыкает систему уравнений (3.4) – уравнение состояния идеального совершенного газа (3.2). Расщепление по пространственным переменным осуществлялось в соответствии с построениями (2.7-2.8).

Для тестирования предложенного численного алгоритма была выбрана задача о сильном взрыве. Моделировалось распространение сферической ударной волны большой мощности в декартовой системе координат, возникшей в результате сильного взрыва, то есть мгновенного выделения в некотором небольшом объёме большого количества энергии. Газ, в котором волна распространяется, считался идеальным. Будем рассматривать волну на расстояниях, не слишком удаленных от источника, а именно там где волна ещё обладает большой интенсивностью. В то же время, эти расстояния полагаются большими по сравнению с размерами источника. Это дает возможность считать, что выделение энергии происходит в одной точке (начале координат). Будем считать, что давление за ударной волной намного больше давлением перед ней.

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

Численные расчеты были проведены для следующих начальных данных: радиус заряда ; значения скорости, плотности и внутренней энергии заряда ; для невозмущенного газа ; уравнение состояния – идеальный газ с , ; число Куранта . Количество сеточных узлов расчетной области постоянно – . Размеры расчетной области изменялись по мере распространения ударной волны от куба с ребрами до куба с ребрами . Заряд помещался в центр куба, который располагался в начале координат.

На рисунке 3.3 показаны зависимость внутренней энергии сразу за фронтом ударной волны от времени в сравнении с аналитическим решением (3.33) для трех различных значений порядка аппроксимации – первого, второго и второго-третьего.



Рис. 3.3 Зависимости внутренней энергии за фронтом ударной волны от времени в сравнении с аналитическим решением (3.5) (сплошная линия) для трех различных значений порядка аппроксимации – первого (бледно серые точки), второго (серые точки) и второго-третьего (черные точки).

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

Глава IV. Обобщение и тестирование полученных высокоточных численных методов на случай уравнений магнитной гидродинамики

где , , .

Для уравнений МГД с сильными разрывами, в принципе, можно построить аналог схемы С.К.Годунова и использовать его в качестве предиктора в рассматриваемых монотонных схемах высокого порядка аппроксимации (2.5) или в качестве самостоятельной схемы, однако, технически это достаточно сложная задача, т.к. необходимо предусмотреть несколько десятков возможных комбинаций из ударных волн, контактных разрыв, волн разрежения, альфвеновских волн, быстрых и медленных магнитозвуковых волн и ряда вырожденных ситуаций (в отличие от уравнений газовой динамики, в которых таких ситуаций существенно меньше – различные комбинации из ударных волн, контактных разрывов и волн разрежения). Поэтому не случайно, что в литературе по МГД отсутствует описание подобных реализаций.

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

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

Итак, при наличии сильных газодинамических разрывов и умеренных магнитных полей, в одномерной нелинейной системе уравнений гиперболического типа в дивергентной форме (4.1) представим вектор магнитогазодинамического потока в виде суммы векторов обычного газодинамического потока и его дополнения до полного вектора , а саму систему (4.1) – в виде суммы двух (гиперболических, в чем можно непосредственно убедиться) операторов

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

Обозначая , каждое из слагаемых в скобках уравнения (4.3)



как и ранее, можно аппроксимировать с использованием монотонных схем 2-3 порядка аппроксимации (2.5), т.к. эквивалентное (4.1) представление (4.2) не изменяет гиперболичности (4.1).

Для 1-го (газодинамического, с сильными разрывами) этапа на каждом шаге по времени в качестве предиктора будем использовать схему Годунова, а для 2-го («магнитного», с умеренными разрывами магнитных характеристик) – также консервативную сеточно-характеристическую схему, в которой не используется точное решение задачи о распаде произвольного разрыва и которая, в этом смысле, более универсальна. В принципе, и для второго этапа в качестве предиктора можно построить аналог схемы С.К. Годунова, однако, без необходимости это делать нецелесообразно, учитывая довольно значительные затраты на итерационное решение задачи о распаде разрыва в этом методе.

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

Наряду с расщеплением (4.2) можно использовать последовательную аппроксимацию одномерных уравнений (4.3) при (аддитивные схемы). При этом шаг интегрирования по времени становится следующим:

Для тестирования предлагаемого метода и реализующего его программных модулей было проведено решение ряда модельных МГД задач с разрывами в начальных данных, с использованием в расчетах на каждом из этапов монотонной разностной схемы 2-3-го порядка аппроксимации (2.5). Некоторые результаты этих тестовых расчетов приведены ниже для плоской одномерной МГД задачи о распаде разрыва – рис. 4.1.

В качестве модельного примера рассматривалась плоская задача о распаде разрыва со следующими начальными данными: , , В данном расчете число узлов равномерной по разностной сетки – число Куранта – . На рис. 4.1 результаты расчетов по схеме (2.5) (штриховые линии) для плотности , компонент вектора скорости давления и компонент вектора напряженности магнитного поля в момент времени сравниваются с результатами известных расчетов (точки) и демонстрируют вполне приемлемое согласие.

Переход к расчетам для 3D уравнений магнитной газодинамики (1.1-1.4) осуществлялся расщеплением по пространственным переменным в соответствии с построениями (2.7-2.8). На заключительном этапе каждого временного шага производился окончательный расчет параметров, в котором учитывалась неоднородность атмосферы Земли и диффузия магнитного поля. При этом использовалось стандартная неявная схема второго порядка аппроксимации для аппроксимации членов с диффузией магнитного поля и метод прогонки для решения соответствующих трехточечных разностных уравнений.

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

Возникающее под действием разлетающейся плазмы МГД-возмущение распространяется на тысячи километров, поэтому существенное влияние на его динамику оказывает не только ближняя от центра энерговыделения область (~100300 км), ионизованная жёстким излучением, но и ионизация в дальней зоне, соответствующая естественным условиям. Поэтому в расчётах в качестве внешних условий задавались естественные значения электронной температуры и степени ионизации , через которые определялся коэффициент диффузии магнитного поля. Отметим, что на км степень ионизации в естественной атмосфере близка к нулю и на эти высоты приходят лишь возмущения геомагнитного поля, среда же остаётся покоящейся. Для правильного расчёта этих возмущений необходимо знать проводимость грунта, и требуется специальная доработка самого расчётного алгоритма. В данной работе вопросы, связанные с поведением возмущённого геомагнитного поля вблизи поверхности Земли не рассматривались.

В расчетах использовалась, декартова система координат с началом в центре Земли так, что (x,z)–плоскость совпадает с плоскостью меридиана, а z–координата с началом на поверхности Земли направлена по вертикали к поверхности Земли и составляет угол с плоскостью экватора. Центр взрыва помещался на положительной части оси на высоте над поверхностью Земли (в южном полушарии). В этой же точке размещался центр начальной области интегрирования (параллелепипед со сторонами ). В последующем, по мере необходимости, линейный размер области интегрирования удваивался (в расчетах на времена в сотни секунд до 14 раз). Область интегрирования покрывалась равномерной сеткой с узлами. В целом, расчеты показали принципиальную возможность решения данной задачи с использованием разработанных методов и программного комплекса на доступных в настоящее время вычислительных системах.

В представленных ниже результатах расчетов приводятся безразмерные значения параметров. За основные характерные параметры принимались: линейный размер =1 км, плотность г/см³ и удельная внутренняя энергия (км/с)² невозмущенной атмосферы в центре области энерговыделения (где К – температура невозмущенной атмосферы в этой же точке). Остальные параметры обезразмерены соответственно: км/с=, дин/см²=, , (с)=.

Линейный размер начальной области интегрирования для 3D вариантов – куб со стороной равной примерно учетверенному радиусу км начальной (сферической) области возмущенной взрывом плазмы, размер которого последовательно удваивался при приближении фронтов возмущений к границам области интегрирования. Внутри начальной возмущенной взрывом сферы с линейным распределением по радиусу от центра сферы задавалось направленное вдоль радиуса начальное распределение скорости , а также постоянные всюду внутри сферы начальные значения плотности г/см³, температуры К и нулевое магнитное поле. Кроме указанных параметров взрыва, задавались высота центра взрыва над поверхностью Земли км и магнитная широта центра взрыва (отрицательная для южного полушария).

По этой же программе проводились расчеты в некоторой квазитрехмерной постановке с цилиндрической геометрией на различных сетках с числом узлов по каждому из пространственных направлений до: M_k=600 (по х), K_k=600 (по z), L_k=10 (по y), что косвенно позволяло провести оценки влияния сеточных параметров на точность трехмерных расчетов (в которых не имелось возможности существенно варьировать сетки), а также правомерность расчетов в такой квазитрехмерной постановке на большие времена (до сотен секунд), где возникают те же проблемы недостаточности имеющихся вычислительных ресурсов. Такие двумерные расчеты соответствуют срединному меридиональному сечению реальных трехмерных расчетов. Хотя при этом ряд количественных характеристик (времена затухания возмущений, окончательные размеры возмущенной области и др.) могут заметно различаться, сопоставление таких двух- и трехмерных расчетов позволяет, в определенной мере, получить представление для реальных трехмерных ситуаций.

	|

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

Расчёты полностью подтверждают принципиальные изменения в характере движения плазмы при высоком начальном энерговыделении. На рисунке 5.3 показаны результаты 3D расчёта для высокого начального энерговыделения (=8,410²³ эрг) на высоте 400 км. На раннее время характер течения плазмы в меридиональной и широтной плоскостях подобны (=0,65 с), то есть геополе не оказывает существенного влияния на течение плазмы, и с течением времени наблюдается ламинарный прорыв плазмы через геомагнитное поле. Последующее образование струйного течения если и происходит, то на очень поздних временах, имеет существенно меньшую протяженность и невысокую интенсивность, так как давление внутри области из-за сильного расширения уже очень мало. Магнитозвуковая волна существенно интенсивнее на всём протяжении её фронта, который вблизи поверхности Земли расположен практически по нормали к ней. На малых высотах плотная атмосфера и практически нулевая проводимость на 50 км "не пропускают" возмущения к поверхности Земли и происходит их распространение преимущественно вдоль этой границы атмосферы, которая повышается над полярными областями. Сопоставление 3D и 2D расчётов на ранние моменты времени для данного варианта показали, что в обоих случаях процессы формирования и эволюции возмущений, оставаясь одинаковыми качественно, количественно отличаются достаточно сильно (например, с запаздыванием для 3D расчётов в 1,52 раза по времени выхода на одинаковые размеры возмущённой области). Как показано ранее с ростом высоты требуется меньшая энергия плазмы для прорыва геомагнитного поля. Расчёты в целом подтверждают этот вывод.

	|

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

Основные результаты работы

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

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

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

Список публикаций по теме диссертации

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. Е.Л. Ступицкий, А.Ю. Репин, А.С. Холодов, Я.А. Холодов. Численное исследование поведения высокоэнергетичного плазменного сгустка в верхней ионосфере. Часть 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 // 31st EPS Conference on Plasma Physics – 2004 – London – Vol. 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. А.С. Холодов, Я.А. Холодов, Е.Л. Ступицкий, А.Ю. Репин. Численные исследования поведения плазменного облака в верхней ионосфере // Математическое моделирование – 2005 – Т. 17, № 11 – C. 43–62.

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

Холодов Ярослав Александрович

РАЗРАБОТКА ВЫСОКОТОЧНЫХ ЧИСЛЕННЫХ МЕТОДОВ РАСЧЕТА ПРОСТРАНСТВЕННОГО ПОВЕДЕНИЯ ПЛАЗМЫ ПОД ДЕЙСТВИЕМ СИЛЬНЫХ МАГНИТОГИДРОДИНАМИЧЕСКИХ ВОЗМУЩЕНИЙ

Подписано в печать 26.10.06

Формат 60x84 1/16. Усл. Печ. Л. 1,1

Тираж 80 экз. Заказ № 467

Московский физико-технический институт

(государственный университет)

НИЧ МФТИ

141700, Московская обл., г. Долгопрудный, Институтский пер. 9

жүктеу/скачать 271.4 Kb.

Достарыңызбен бөлісу: