Содержание к диссертации
Введение
1 Предварительный анализ 39
1.1 Одномерные расчетные схемы 39
1.1.1 Модифицированная схема Лакса-Фридрихса . 43
1.1.2 Схема коррекции потоков на основе решения линеаризованной задачи Римана 46
1.1.3 Схема Курганова-Тадмора 49
1.1.4 Сравнение результатов расчетов тестовых задач . 52
1.1.4.1 Тестовая задача 1 — «магнитный поршень» 53
1.1.4.2 Тестовая задача 2 — распад МГД-разрыва 54
1.1.4.3 Тестовая задача 3 — распространение ударной волны по холодному фону 55
2 Разработка МГД-солвера 61
2.1 Система уравнений двухтемпературной МГД 61
2.2 Алгоритм построения двойственных сеток для сложных двумерных областей 67
2.2.1 Формальная постановка задачи 69
2.2.2 Алгоритм построения сетки контрольных объемов 71
2.2.3 Основные достоинства предлагаемого метода построения двойственной сетки 77
2.3 Расчет диффузии магнитного поля 78
2.3.1 Вывод уравнений диффузии магнитного поля в случае изотропной проводимости СГц = <7 L = <7 . 78
2.3.2 Вывод уравнений диффузии магнитного поля в случае анизотропной проводимости сгц ф а± . 80
2.3.3 Расчетная схема для плоскосимметричного случая 86
2.3.3.1 Аппроксимация Е* в серединах ребер двойственной сетки 86
2.3.3.2 Аппроксимация rot Е* в узле треугольной сетки 94
2.3.3.3 Сборка матрицы для расчета диффузии магнитного поля 98
2.3.3.4 Решение системы линейных уравнений 105
2.3.3.5 Аппроксимация вектора напряженности электрического поля в узлах 113
2.3.4 Расчетная схема для случая осевой симметрии (ци
линдрическая геометрия) 114
2.4 Расчет идеальной двухтемпературной МГД 116
2.4.1 Решаемая система уравнений в плоскосимметричном случае 116
2.4.2 Расчетная схема для плоскосимметричного случая 119 2.4:2.1 Немонотонная реконструкция на этапе
«предиктор» 121
2.4.2.2 Квазимонотонная реконструкция на этапе «корректор» 123
2.4.2.3 Этап «предиктор» 126
2.4.2.4 Этап «корректор» 132
2.4.2.5 Вычисление потоков Fy 137
2.4.2.6 Построение явной вычислительной схемы 141
2.4.3 Расчетная схема для случая осевой симметрии (цилиндрическая геометрия) 144
2.5 Расчет теплопроводности и ионно-электронного обмена . 145
2.6 Расчет джоулева нагрева и излучения 147
2.6.1 Расчет джоулева нагрева 147
2.6.2 Расчет излучения 148
2.7 Расчет в околовакуумных (силыюразреженных) областях 150
2.7.1 Проблема расчета в околовакуумных областях . 150
2.7.2 Обеспечение сквозного расчета в околовакуумных областях 152
2.7.2.1 Общий вид корректирующих функций . 154
2.7.2.2 Корректировка проводимости а 155
2.7.2.3 Корректировка джоулева нагрева и излучения 156
2.7.2.4 Корректировка переноса (истечения) из околовакуумных ячеек 158
2.7.2.5 Корректировка уравнения импульса . 161
2.7.3 Оценка требуемого значения <хвакуум 163
2.8 Сборка общего МГД-солвера 168
2.8.1 Программная система единиц (ПСЕ) 168
2.8.2 Итерация общего МГД-солвера 170
2.8.3 О выборе временного шага At 177
3 Расчеты задач мод 180
3.1 Разгон и сжатие плазменной оболочки (лайнера) 180
3.1.1 Разгон невозмущенного плазменного лайнера . 186
3.1.2 Разгон плазменного лайнера с возмущениями . 197
3.2 Сжатие магнитного потока (MFC) 212
3.2.1 Постановки с одним лайнером (компрессором) . 216
3.2.2 Постановки с двумя лайнерами (компрессором и нагрузкой) 229
Заключение 243
Список литературы
- Схема коррекции потоков на основе решения линеаризованной задачи Римана
- Вывод уравнений диффузии магнитного поля в случае анизотропной проводимости сгц ф а± .
- Квазимонотонная реконструкция на этапе «корректор»
- Разгон невозмущенного плазменного лайнера
Введение к работе
Страница 6
ВВЕДЕНИЕ
Актуальность проблемы
Современные фундаментальные и прикладные исследования в области физики плотной высокотемпературной плазмы характеризуются большим и постоянно возрастающим объемом вычислительных экспериментов на компьютерах [1], [2]. Бурное развитие вычислительной физики плазмы оказалось возможным благодаря широкому внедрению новых ЭВМ, высокопроизводительных аппаратных средств для выполнения параллельных вычислений и сетевых технологий. С другой стороны, это развитие обусловлено и возрастающими потребностями анализа нелинейных и неравновесных плазменных процессов, т.е. усложнением собственно математических моделей плазмы различных видов.
Для проведения численных исследований различных задач создаются специализированные программные комплексы, которые могут включать инструментальные средства препроцессинга (предварительной подготовки данных перед расчетом), средства проведения расчетов (солверы), а также средства постпроцессинга (анализа результатов расчета). Прогресс в области компьютерной техники (повышение вычислительных мощностей, развитие мультимедийных средств) приводит к возможности и даже необходимости модернизации программного обеспечения. Например, широкое распространение параллельных компьютерных систем стимулирует создание и развитие программных
ВВЕДЕНИЕ
Страница 7
средств, реализующих параллельные алгоритмы расчетов и визуализации результатов. Развитие графических возможностей компьютеров позволяет поднять на новый уровень инструментальные средства пре-процессинга и постпроцессинга.
Очевидно, что прогресс в развитии научно-технического программного обеспечения обусловлен не только усовершенствованием вычислительной техники, но и развитием самого объекта численного исследования. Под этим мы подразумеваем совершенствование известных и создание новых физико-математических моделей, уточнение известных данных по уравнениям состояний веществ, появление новых экспериментальных результатов, для объяснения которых необходимо изучение новых постановок задач. Большую роль играет и создание новых крупных плазменных установок (типа ITER или NIF). Все вышеуказанные факторы приводят к необходимости значительной модификации существующего программного обеспечения (ПО) или создания новых прикладных программ.
Более того, существует еще один немаловажный фактор стимулирования, не связанный непосредственно ни с прогрессом вычислительной техники, ни с развитием объектов численного анализа. Этим фактором является развитие численных методов: усовершенствование существующих и создание новых. Например, в настоящее время весьма широкое распространение получили различные алгоритмы на неструктурированных расчетных сетках: тетраэдральных, призматических, треугольных и т.д. [27]. Ранее, как известно, в расчетных программах доминировали структурированные (часто даже просто кубические или прямоугольные) сетки. Использование новых классов расчетных сеток привело к пересмотру вычислительных алгоритмов и, соответственно, к модернизации соответствующего ПО.
Для исследований течений плотной плазмы применяются модели га-
ВВЕДЕНИЕ
Страница 8
зовой динамики (ГД) и магнитной газовой динамики (МГД). В рамках этих моделей плазма представляется континуальной (непрерывной) квазинейтральной средой, состоящей из ионов различных сортов, электронов и, в общем случае, нейтральных частиц. Уравнения МГД выводятся в кинетической теории в предположении локального термодинамического равновесия компонент плазмы.
Характерные временной масштаб Tq и пространственный масштаб
Lo, при которых применимы ГД- и МГД-модели, имеют [58] следующие
1
ограничения: по времени То ^ > гДе ^плазма — плазменная часто-
^плазма
та, а по пространству — Lq ^$> Яцсбая, где Ддебая — радиус Дебая. Указанные ограничения не позволяют использовать ГД- и МГД-модели для исследования высокочастотных движений плазмы. Мощность экспериментальной аппаратуры в современных исследованиях высокотемпературной плазмы является величиной порядка нескольких тераватт, причем нередко она выделяется за весьма короткое время порядка 102-f-103 наносекунд, а размеры исследуемых объектов варьируются от Ю-2 см до нескольких десятков сантиметров. При этом плазмодинамические ГД- и МГД-модели достаточно широко и в целом успешно применяются для аналитического и численного анализа и прогнозирования результатов многих физических экспериментов.
Задачи магнитной гидродинамики плотной плазмы можно условно разделить на два семейства, используя в качестве критерия существенность процессов взаимодействия с собственным магнитным полем плазмы по сравнению с процессами взаимодействия с внешним магнитным полем. К первому классу будем относить задачи, в которых зарождение и развитие плазменных образований осуществляются в условиях весьма существенного влияния собственного магнитного поля плазмы, по которой течет электрический ток. К этому семейству относятся задачи изучения плазменных объектов, возникающих в электроразрядных
ВВЕДЕНИЕ
Страница 9
и электродинамических устройствах. Имеются в виду такие плазменные объекты, как различного рода самосжатые разряды, Z-пинчи, плазменные лайнеры и т.п. К этой же группе относятся различные задачи генерации в импульсной плазме спонтанных магнитных полей.
Ко второму классу можно отнести физические задачи, в которых основную роль играют процессы взаимодействия плотной импульсной плазмы с внешним магнитным полем, а процессы, связанные с собственным магнитным полем плазмы, не играют ключевой роли. Внешнее магнитное поле может быть создано некоторым обособленным устройством, таким как генератор токовых импульсов и т.п. Например, ко второму семейству можно отнести многие задачи магнитного удержания плазмы, задачи формирования разрядных каналов с требуемыми свойствами для различных прикладных целей, задачи инжекции термоядерного топлива в системы с магнитным удержанием плазмы и т.д.
Очевидно, что деление физических задач на указанные выше классы достаточно условно. Более того, в каждом отдельном случае ситуация может динамически меняться с течением времени. Например, роль процессов, связанных с собственным магнитным полем плазмы, может постепенно расти, и, начиная с некоторого момента, задача станет относиться не ко второму классу, а к первому.
Постановки задач численного моделирования плазмы, рассматриваемых в первую очередь в настоящей диссертации, возникли из прикладных исследований физического метода электродинамического ускорения и сжатия плазмы по схеме «Z-пинч». Задачи динамики внешней токонесущей оболочки пинча (лайнера) относятся к первому семейству, в то время как задачи моделирования процесса формирования и распространения электромагнитной волны, порожденной сильноточным ко-роткоимпульсным генератором, следует отнести ко второму классу.
Электродинамические разгон и сжатие плазменных оболочек были
ВВЕДЕНИЕ
Страница 10
впервые предложены в 60-х годах XX века [59]. Этот метод предназначался, во-первых, для нагрева дейтериево-тритиевой смеси до термоядерных температур и, во-вторых, для создания мощных импульсов рентгеновского излучения за счет преобразования кинетической энергии ускоренной плазменной оболочки в энергию излучения [60]. Создание подобных источников рентгеновского излучения требуется как для различных физических исследований, так и для развития новых технологических процессов (например, рентгеновской фотолитографии). В России и за рубежом уже имеются генераторы электромагнитных импульсов длительностью порядка 102 наносекунд при мощности порядка 10 и более тераватт. Это известные установки PBFA-Z [61], Saturn [62], Ангара-5-1 [63], [81] и другие. В последние годы активно исследуется схема нагрузки сильноточного генератора, называемая быстрым Z-пинчем (или плазменным лайнером).
Одним из наиболее впечатляющих результатов, полученных в последнее десятилетие в области высоких импульсных мощностей (Pulsed-Power), явилась успешная разработка источника мягкого рентгеновского излучения в Sandia National Laboratories (США). Плазма, инициированная взрывом цилиндрической сборки из 240 вольфрамовых проводников диаметром 7.5 мкм, оказалась способной конвертировать магнитную энергию электрофизической установки в излучение с полным КПД порядка 10% и обеспечить импульс излучения в мягком рентгеновском диапазоне с энергией ~ 2 МДж и мощностью более 200 ТВт. Электротехническая установка, потребовавшаяся для создания такого излучения, обеспечивала амплитуду тока в нагрузке ~ 20 МА со временем его нарастания 100 не. Этот результат открыл новый реалистичный подход к исследованиям в области физики высоких энергий и инерционного управляемого термоядерного синтеза (УТС) [31].
Исходя из полученного результата, стало возможным предположить,
ВВЕДЕНИЕ
Страница 11
что при использовании вышеупомянутой нагрузки на уровне токов до 60 МА, можно было бы получить условия, необходимые для реализации УТС посредством непрямого обжатия термоядерной мишени. Характерно, однако, что строительство такого электрофизического генератора не было немедленно начато. Дело в том, что технология получения импульсных токов, использованная в установке лабораторией Sandia (увеличение мощности посредством стадийного накопления энергии электрического поля), имеет внутренние ограничения, связанные с максимальной электрической прочностью существующих диэлектриков. Альтернативой этому подходу, который развивается уже более 50 лет, могли бы послужить различные схемы, появившиеся относительно недавно и использующие стадийное накопление не электрической, а магнитной составляющей энергии электромагнитного поля. В этом случае магнитная энергия накапливается в вакуумных полостях, и максимальная плотность мощности уже не ограничена электрическим полем пробоя материалов, а зависит от механической прочности проводящих стенок полости.
Одной из схем, реализующих накопление магнитной энергии с последующей ее концентрацией в пространстве и времени, является схема сжатия магнитного потока (magnetic flux compression, MFC). Этот метод генерации больших магнитных полей и больших токов был предложен уже давно, ещё в конце 50-х годов прошлого столетия независимо в СССР [65] и в США [66].
Концепция MFC состоит в сжатии полости с проводящими стенками настолько быстро, чтобы запертый в ней магнитный поток не успевал продиффундировать наружу через проводящую поверхность. В этом случае магнитное поле в полости растет, как следствие сохранения внутреннего магнитного потока, обратно пропорционально её уменьшающемуся сечению. В цилиндрической геометрии возможно сжатие либо
ВВЕДЕНИЕ
Страница 12
аксиального магнитного поля, либо азимутального поля, захваченного между центральным электродом («статором») и сжимающейся проводящей стенкой («ротором»). Первая конфигурация, первоначально предложенная академиком А. Д. Сахаровым [65], позволяет получать сверхвысокие магнитные поля (до ~ 30 МГс), тогда как вторая может быть использована для получения сверхвысоких токов в полезной нагрузке (до ~ 100 МА).
Для сжатия полости с магнитным полем долгое время использовалась химическая энергия взрывчатых веществ (поэтому такие генераторы получили название взрывомагнитных — ВМГ) [67], что ограничивало максимальную скорость ускоряемой внешней проводящей стенки на уровне ~ 5 км/с. При первоначальном внешнем радиусе магнитной полости порядка 5 см (ограничен сверху максимальной допустимой амплитудой рэлей-тейлоровской неустойчивости, развивающейся во время сжатия) взрывомагнитные генераторы способны, таким образом, создавать импульсные магнитные поля и токи микросекундного диапазона.
С другой стороны, решение задачи УТС с источником излучения на основе магнитной имплозии взрываемых многопроволочных сборок требует на порядок меньших значений времени нарастания тока. При этом скорость сжатия полости (т.е. скорость внешней стенки магнитной полости) должна достигать ~ 100 км/с. На взрывомагнитных генераторах подобные скорости кажутся недостижимыми.
В то же время такие параметры являются типичными для плазменных оболочек малой массы, ускоряемых магнитным полем. В этом случае максимальная скорость тонкой плазменной стенки магнитной полости уже не ограничена скоростью разлета продуктов взрыва, а скорее может быть оценена как альфвеновская скорость в плазме оболочки при данном значении ускоряющего поля. Магнитный поток сжимается в области между неподвижным «статором» и внешней оболочкой
ВВЕДЕНИЕ
Страница 13
(«ротором»), ускоряемой импульсным магнитным полем, создаваемым электрофизическим генератором. Ожидаемое время сжатия оболочки настолько мало, что запертое между «статором» (электрод) и «ротором» (плазма) магнитное поле не успевает значительно продиффунди-ровать в плазму. Данный подход (его можно условно назвать «плазменным MFC») был предложен сравнительно недавно [40]. Поэтому схема плазменного MFC ещё слабо изучена экспериментально.
Одним из практических приложений схемы плазменного сжатия магнитного потока может стать создание новых установок для исследования проблемы УТС в соответствии с [31]. Другим приложением рассматриваемой схемы является создание мощных источников рентгеновского излучения. Наконец, предлагается использовать MFC-схемы вместо традиционных взрывомагнитных генераторов для создания импульсных давлений магнитного поля с целью экспериментального изучения уравнений состояния веществ методами безударного и ударного сжатия [68].
Несмотря на различия возможных экспериментальных установок MFC, обусловленные различными сферами их применения, в них можно выделить общие элементы (по крайней мере, на некотором уровне абстрагирования). Во-первых, это электрофизический генератор тока и основной плазменный лайнер-компрессор. Внешняя поверхность лайнера-компрессора и цепь генератора тока образуют первичный электрический контур. Лайнер-компрессор разгоняется и сжимается к оси системы магнитным полем, порождаемым кратковременным импульсом тока генератора. Вторичный электрический контур образуется электродами установки, внутренней поверхностью плазменного компрессора и нагрузкой. Магнитный поток, запертый во вторичном контуре, сжимается за счет имплозии компрессора, что ведет к росту давления запертого магнитного поля.
ВВЕДЕНИЕ
Страница 14
В качестве нагрузки может, например, использоваться второй плазменный лайнер (в задачах по исследованию УТС или по созданию радиационного источника) или вакуумная индуктивность (в задачах по изучению уравнений состояния материалов). В первом случае сжатие магнитного потока во вторичном контуре приводит к имплозии лайнера-нагрузки к оси системы. Во втором случае — к резкому нарастанию давления магнитного поля в нагрузке (вакуумной области); передача этого давления на исследуемый образец позволяет реализовать требуемый режим его сжатия.
Сжатие лайнера-компрессора используется для усиления мощности токового импульса электротехнического генератора, в первую очередь, за счет существенного (в несколько раз) уменьшения длительности импульса. Как было сказано выше, подобная схема с накоплением энергии магнитного поля представляется в настоящий момент более перспективной, чем классическая схема накопления энергии электрического поля, например, в конденсаторных батареях.
Для корректного согласования параметров плазменных пинчей (лайнеров) с параметрами генераторов импульсов тока, для профилирования формы электродов разрядных камер и т.п. необходим численный расчет электродинамического разгона и сжатия пинчей. Численное моделирование позволяет получить обоснованные количественные оценки энергетического баланса установки, устойчивости плазменных объектов, характеристики излучения и т.д., а также осуществлять оптимизацию параметров установок.
Необходимость развития старых и создания новых численных кодов для расчета задач плазмодипамики и, в частности, для расчета задач разгона и сжатия плазменных лайнеров (в том числе задач плазменного MFC) вызвана разными причинами.
Во-первых, это естественное требование увеличения точности чис-
ВВЕДЕНИЕ
Страница 15
ленных решений. Здесь имеются в виду как распределения расчетных величин по пространству, так и интегральные характеристики решений.
Во-вторых, это требование учета последних достижений физической теории, физического и численного эксперимента. Например, это может быть появление уточненных или расширенных таблиц уравнений состояния вещества, получение улучшенных аналитических оценок параметров плазмы и т.п.
Третья причина — это изменение постановок самих физических задач. Поясним это на примере интересующих нас в первую очередь задач разгона и сжатия лайнеров. Раньше вполне удовлетворительными считались расчеты, в которых расчетная область представляла собой простейшую фигуру — например, прямоугольник или призму (имеется в виду сечение плоскостью RZ цилиндрической системы координат; предполагается осевая симметрия). При дискретизации областей с такой простой геометрией использовались простейшие индексные, т.е. структурированные, расчетные сетки [27]. В современных постановках лайнерных задач расчетная область имеет достаточно сложную геометрию. Например, это может быть система из двух или нескольких разрядных камер, соединенных узкими каналами сложной формы. Характерные размеры элементов таких расчетных областей могут меняться в десятки, сотни раз и даже более. Используются криволинейные профили (поверхности) электродов. Использование простейших индексных сеток для дискретизации таких областей не всегда приемлемо. Если эти сетки и удается построить, то в них, вследствие их структуры, содержится физически неоправданно большое число расчетных узлов, возникают проблемы с границами сложной формы и т.п. Это приводит к тому, что старые коды оказываются непригодными для решения задач в областях сложной формы. Требуется создание кодов, использующих современные сеточные технологии. Например, в настоящей диссертации рассматривается
ВВЕДЕНИЕ
Страница 16
код MARPLE, в котором реализована поддержка неструктурированных треугольных сеток.
Четвертой и немаловажной причиной является необходимость учета современных требований к прикладному программному обеспечению (ПО). В частности, требуется обеспечение надежности ПО (за счет использования современных систем управления качеством разработки [64]), повышение производительности ПО (что связано с развитием аппаратного и программного обеспечения компьютеров), уменьшение затрат на модификацию ПО (за счет использования современных подходов к разработке) и т.д. Заметим, что обеспечение высокой производительности ПО на современных аппаратно-программных платформах особенно важно при проведении серийных расчетов (например, при оптимизации профиля электрода) и при наличии жестких требований относительно допустимого времени расчета.
В последние годы разработано достаточно много численных схем, позволяющих достичь высокой точности расчета на нерегулярных сетках [3], [19], [12], [72] (см. также краткий обзор разностных схем в разделе 1.1 настоящей работы). Это открывает возможности широкого внедрения в вычислительной физике и, в частности, в вычислительной илаз-модинамике новых сеточных технологий (блочные сетки, неструктурированные сетки и т.д.).
Разработка подобных технологий (построение вычислительных схем и алгоритмов, их программная реализация на современном уровне и интеграция в составе пакета прикладного ПО), а также их введение в практику вычислительных экспериментов в области физики сильноточных разрядов является некоторой идеальной «глобальной» целью настоящей диссертационной работы.
ВВЕДЕНИЕ
Страница 17
Цели и задачи
Целью работы являются: создание моделей и вычислительных методов для моделирования МГД-процессов в диссипативной двухтемпе-ратурной плазме на неструктурированных сетках, разработка и реализация соответствующего программного обеспечения, расчет ряда задач сжатия плазменных лайнеров и плазменного сжатия магнитного потока (MFC) в областях нетривиальной геометрии.
Настоящая диссертация посвящена следующим вопросам:
Разработка и анализ разностных схем суммарной аппроксимации МГД-уравнений на треугольных неструктурированных сетках и алгоритмов решения сеточных уравнений с целью моделирования процессов в плотной электродинамически ускоряемой плазме.
Программная реализация разработанных алгоритмов в коде MARPLE. Разработанное программное обеспечение ориентировано в первую очередь на задачи моделирования разгона и сжатия плазменных лайнеров, в том числе задачи MFC, однако оно также применимо и к другим задачам расчета плотной плазмы. При разработке программного обеспечения использован современный объектно-ориентированный подход; языком реализации является C++. Использование объектно-ориентированного подхода приводит к повышению удобочитаемости кода, существенному уменьшению затрат на поддержку и модификацию кода, к возможности повторного использования кода.
Использование созданного программного обеспечения для проведения вычислительных экспериментов, посвященных исследованию задач плазмодинамики:
1. Задачи о разгоне кратковременным (102 + 103 не) токовым им-
ВВЕДЕНИЕ
Страница 18
пульсом плазменных лайнеров, как однородных, так и с внесенными начальными возмущениями плотности. Эти задачи используются главным образом для оценки качества работы созданного кода и сравнения результатов с другими источниками (валидация кода).
Задачи о сжатии ускоряемым плазменным лайнером магнитного потока в экспериментальных установках сложной геометрии. Использование лайнера позволяет провести усиление токового импульса генератора за счет существенного уменьшения его длительности. Проведение численных экспериментов позволяет еще на этапе проектирования выявить потенциальные проблемы в работе и оптимизировать параметры устройств (профили электродов, формы токовых импульсов, параметры плазменных лайнеров и т.д.).
Задачи о сжатии магнитного потока в установках сложной геометрии в постановках с двумя плазменными лайнерами. Этот класс задач отличается от предыдущего тем, что в качестве нагрузки используется дополнительный подвижный лайнер. Основной лайнер-компрессор, ускоряемый токовым импульсом от генератора, сжимает в режиме усиления мощности магнитный поток, который, в свою очередь, инициирует движение лайнера-нагрузки. Целью является магнитная имплозия лайнера-нагрузки, в то время как лайнер-компрессор служит только для усиления токового импульса генератора за счет уменьшения его длительности. Численные эксперименты предназначены для тех же целей, что и в задачах предыдущего класса, — диагностика и оптимизация параметров установки на этапе проектирования.
ВВЕДЕНИЕ Страница 19
Научная новизна
В данном разделе сформулированы научно-технические задачи, которые были решены лично автором настоящей диссертации в рамках проекта создания кода MARPLE (Magnetically Accelerated Radiative PLasma Explorer). Разработка кода ведется в ИММ РАН.
1. Разработан ряд сеточных и вычислительных алгоритмов, предназначенных для работы на нерегулярных неструктурированных сетках:
Предложен оригинальный алгоритм построения двойственных сеток (сеток контрольных объемов), необходимых для расчетов методами конечных объемов.
Построены различные варианты явной двухэтапной разностной схемы расчета идеальной МГД с использованием немонотонных и квазимонотонных реконструкций сеточных величин на разных этапах. Предложенная схема является экономичной с вычислительной точки зрения и при этом обеспечивает достаточно хорошую точность расчета.
С помощью интегро-интерполяционного метода построена неявная разностная схема расчета анизотропной диффузии магнитного поля на треугольных сетках. Использование современных программных средств для обращения матриц позволяет получить высокоточные результаты за приемлемое время.
Предложена двухэтапная схема общего МГД-солвера, реализующего принцип суммарной аппроксимации для рассматриваемой системы уравнений двухтемпературной диссипативной МГД. Предложены специальные околовакуумные корректировки, позволяющие проводить сквозной расчет в областях
ВВЕДЕНИЕ
Страница 20
очень низкой плотности за счет определенной корректировки базовых уравнений МГД.
Разработанный программный код MARPLE сочетает развитую и расширяемую физико-математическую модель, использование современных технологий неструктурированных сеток и вычислительных методов, использование объектно-ориентированного подхода (ООП). Все это открывает широкие возможности использования кода для научных и практических приложений, а также перспективы его дальнейшего развития.
С помощью созданного программного обеспечения проведен ряд оригинальных вычислительных экспериментов. В частности, расчеты рассмотренных нами задач плазменного сжатия магнитного потока (MFC) в областях нетривиальной геометрии и с привлечением развитых физико-математических моделей ранее не проводились. Сделанные ранее расчеты либо выполнялись в областях самой тривиальной геометрии, либо не учитывали многих важных МГД-процессов. На конкретных примерах демонстрируется недостаточность расчетов задач MFC по упрощенным моделям и необходимость проведения полноценных двумерных или трехмерных расчетов.
Научная и практическая значимость
Научная и практическая значимость полученных результатов обусловлена следующими причинами:
1. В разработанном программном обеспечении (код MARPLE) сочетаются:
ВВЕДЕНИЕ
Страница 21
Достаточно развитая и при этом легко расширяемая физико-математическая модель. В модели учтены анизотропная диффузия магнитного поля, джоулев нагрев, ионная и электронная теплопроводности и ионно-электронный обмен, излучение, двухтемпературность.
Использование современных численных методов и сеточных технологий (так называемые неструктурированные сетки), позволяющих вести расчет в областях практически любой геометрии, в том числе с меняющимися на несколько порядков характерными размерами элементов. В настоящей работе под неструктурированными сетками понимаются нерегулярные сетки произвольной структуры из правильно примыкающих треугольных элементов. Возможно использование сеток в несколько сотен тысяч элементов (треугольников) в расчете на современном серийном ПК.
Использование современного объектно-ориентированного подхода (ООП), что повышает удобочитаемость исходного кода, уменьшает затраты на его модификации, способствует повторному использованию кода и т.д.
Таким образом, разработанный код приспособлен для решения современных практически интересных физических задач в областях очень сложной геометрии на сетках с большим количеством элементов. При этом за счет использования объектно-ориентированного подхода существенно облегчена модификация кода, что положительно сказывается на перспективах дальнейшего развития.
2. Рассчитанные с помощью разработанного кода MARPLE задачи плазменного сжатия магнитного потока (MFC), представляя самостоятельный практический интерес, являются также прототипа-
ВВЕДЕНИЕ
Страница 22
ми более сложных и еще более интересных в практическом плане задач этого класса. Схема плазменного сжатия магнитного потока может быть использована во многих приложениях: определение уравнений состояния методами безударного и ударного сжатия, получение рентгеновских источников, исследования по управляемому термоядерному синтезу. Указанные приложения, наряду с развитием новых численных методов магнитной гидродинамики, определяют практическую значимость и актуальность настоящей работы. Сравнение наших двумерных расчетов задач MFC с расчетами по упрощенным физико-математическим моделям убедительно доказало необходимость проведения полноценных двумерных и трехмерных расчетов, особенно для конфигураций со сложной геометрией. При этом такие полноценные расчеты для рассмотренных классов задач MFC ранее еще не были проведены. Таким образом, нет сомнений и в необходимости проведения расчетов задач MFC с помощью кода MARPLE. Как было сказано выше, это особенно актуально для расчетных областей с нетривиальной геометрией.
3. Многие разработанные для MARPLE вычислительные алгоритмы допускают обобщение на случай трехмерных неструктурированных сеток (например, тетраэдральных). Такая работа ведется уже в настоящее время, в частности, в рамках проекта ISTC-2830 (Международного научно-технического центра). Создание трехмерной версии кода существенно облегчается за счет использования в MARPLE объектно-ориентированного подхода (что способствует повторному использованию кода и т.д.). Разработанные модели, методы и непосредственно программное обеспечение ускоряют и облегчают ведущуюся в настоящее время разработку трехмерной версии кода, что имеет большое практическое значение.
ВВЕДЕНИЕ Страница 23
В силу вышеизложенного можно заключить, что совокупность полученных результатов (модели, вычислительные методы, программное обеспечение, расчеты задач сжатия лайнеров и плазменного сжатия магнитного потока) имеет как научную, так и практическую значимость. При этом открываются широкие перспективы дальнейшего развития моделей, методов, программного обеспечения, а также проведения расчетов других интересных задач плазмодинамики.
Основные положения, выносимые на защиту
Разработано программное обеспечение MARPLE для моделирования МГД-процессов в плотной электродинамически ускоряемой плазме, в котором использована технология неструктурированных сеток и современный объектно-ориентированный подход к созданию кода.
Разработаны и программно реализованы численные алгоритмы расчета идеальных и диссипативных МГД-процессов на неструктурированных сетках, разработан и реализован общий МГД-солвер, использующий принцип суммарной аппроксимации.
Предложена и реализована методика обеспечения сквозного расчета газодинамических величин и электромагнитного поля в околовакуумных (с низкой плотностью вещества) областях, практически всегда встречающихся в задачах разгона и имплозии лайнеров и задачах плазменного сжатия магнитного потока.
Проведены расчеты задач разгона и имплозии плазменного лайнера, как однородного, так и с гармоническими возмущениями плотности. Параметры расчетной установки соответствовали известному эксперименту Сандийской Национальной лаборатории (США)
ВВЕДЕНИЕ
Страница 24
на установке «Z». Данный расчет можно рассматривать как вали-дацию разработанного программного обеспечения. Сравнение результатов с другими источниками (экспериментальные данные, численные эксперименты на ЭВМ) показало совпадение с практически приемлемой точностью. Это позволяет использовать код MARPLE для прогноза результатов экспериментов, выполняемых на современных сильноточных устройствах.
Проведены расчеты задач плазменного сжатия магнитного потока (MFC) в различных постановках. В качестве нагрузки схемы MFC, обостряющей токовый импульс, были использованы как вакуумная индуктивность, так и дополнительный плазменный лайнер. Осуществлено сравнение с расчетами по упрощенной физической модели (упрощенная геометрия, пренебрежение диссипацией). Показано, что в задачах сжатия магнитного потока возможно увеличение мощности токового импульса более чем на порядок. Продемонстрировано отрицательное влияние на работу схемы MFC эффекта проникновения плазмы лайнера-компрессора в область нагрузки. С помощью численного эксперимента проверена эффективность способа предотвращения указанного эффекта, заключающегося в увеличении затравочного тока во вторичном (нагрузочном) контуре.
Проведенные расчеты доказывают пригодность разработанного программного обеспечения для расчета задач разгона и имплозии лайнеров и плазменного сжатия магнитного потока в областях сложной геометрии. Использование неоднородных неструктурированных расчетных сеток позволяет моделировать физические и инженерные плазменные установки высокой степени геометрической сложности, содержащие существенно разномасштабные элементы.
ВВЕДЕНИЕ
Страница 25
Достоверность результатов
Достоверность полученных результатов подтверждается, во-первых, сравнением их с физическими экспериментами и другими расчетами. В частности, результаты наших расчетов разгона и сжатия лайнера (см. раздел 3.1) достаточно хорошо согласуются с известными результатами, полученными с помощью двумерного лагранжево-эилерова кода на структурированных сетках (РАЗРЯД).
Результаты расчета более интересных в практическом плане задач сжатия магнитного потока MFC (см. раздел 3.2) сравнивались только с сильно упрощенными физическими оценками. Сравнение с кодом РАЗРЯД для этих задач оказалось невозможным по той причине, что данный код малопригоден для их решения. Это вызвано в первую очередь использованием в нем устаревших сеточных технологий (индексные моноблочные сетки). Кроме того, модификация кода РАЗРЯД затруднена процедурным подходом, использованным при его создании.
Таким образом, сравнение наших расчетов задач MFC с другими полноценными двумерными или трехмерными расчетами невозможно за отсутствием таковых. Однако сравнение даже с весьма упрощенными физическими оценками демонстрирует достоверность и непротиворечивость получаемых результатов. При этом наши расчеты позволили исследовать важные процессы, которыми упрощенные физические модели полностью пренебрегают. Естественно, такого рода отличия от примитивных оценок никак не могут снизить достоверность наших результатов.
Во-вторых, следует заметить, что каждый из программных модулей MARPLE подвергался автономному тестированию. В качестве тестов использовались задачи, для которых известны аналитические или высокоточные численные решения. Таким образом, блок расчета идеальной
ВВЕДЕНИЕ
Страница 26
МГД тестировался на ряде задач идеальной МГД, блок расчета диффузии магнитного поля тестировался на ряде задач диффузии поля и т.д.
В-третьих, достоверность проведенных расчетов косвенно подтверждается теоретическими оценками использованных разностных схем (аппроксимация, устойчивость, сходимость) и алгоритмов.
Вышесказанное позволяет сделать заключение о достоверности полученных в работе результатов.
Апробация результатов
Основные результаты были представлены на следующих конференциях и семинарах:
XXX Звенигородская конференция по физике плазмы и УТС, г. Звенигород, 24-28 февраля 2003 года.
13th International Symposium on High Current Electronics, Tomsk, Russia, 25-30 July 2004.
Семинар лаборатории LPTP (физика и технология плазмы) Ecole Polytechnique под руководством проф. J.M. Rax, 91128 Palaiseau Cedex, France, июль 2005 г.
The 25th International Symposium on Shock Waves - ISSW25, Bangalore, India, 17th - 22nd July, 2005.
International conference on parallel computational fluid dynamics, PARCFD, 2005.
International Conference on Parallel Computing, PARCO-2005, Malaga, Spain, 13-16 September, 2005.
ВВЕДЕНИЕ Страница 27
Третий международный научный семинар LPM : Математические модели и моделирование в лазерно-плазменных процессах, МОСТУ, Москва, Россия, 31 января - 4 февраля 2006 года.
Совместный семинар сотрудников ИММ РАН, ИТЭФ РАН, ФИ РАН и «Нейрок Техсофт» под руководством к.ф.-м.н. Диянко-ва О.В., г. Троицк, 31 января 200G года.
XXXIII Звенигородская конференция по физике плазмы и УТС, г. Звенигород, 13-17 февраля 2006 года.
Совместный семинар сотрудников ИММ РАН и ТРИНИТИ под руководством к.т.н. Грабовского Е.В., Троицк, февраль 2006 года.
Семинар ИММ РАН под руководством проф., д.ф.-м.н. Левано-ва Е.И., Москва, март 2006 года.
Личный вклад автора
Работа автора является частью проекта создания кода MARPLE, выполняемого коллективом сотрудников восьмого отдела ИММ РАН (Га-силов В.А., Болдарев А.С, Дьяченко СВ., Карташева Е.Л., Ольховская О.Г.). Автором настоящей диссертации лично решены следующие задачи:
Разработка специального алгоритма построения двойственной сетки (см. раздел 2.2). Программная реализация соответствующих классов.
Разработка общего МГД-солвера (имеется в виду основной управляющий модуль, см. раздел 2.8). Программная реализация соответствующих классов.
ВВЕДЕНИЕ
Страница 28
Разработка блока расчета анизотропной диффузии магнитного поля (см. раздел 2.3). Программная реализация соответствующих классов.
Разработка блока расчета идеальной двухтемпературной МГД (см. раздел 2.4). Программная реализация соответствующих классов. К данному пункту можно отнести также раздел 1.1, посвященный предварительному анализу ряда существующих одномерных разностных схем расчета идеальной МГД, и раздел 2.7, в котором предлагаются специальные корректирующие функции, предназначенные для обеспечения сквозного расчета в областях с низкой плотностью вещества (например, в низкоплотной плазменной короне или в остаточном газе).
Проведение вычислительных экспериментов (см. главу 3). Анализ полученных численных результатов. Постановка физических задач и анализ результатов расчетов проводились, как правило, совместно с Чуватиным А.С. (лаборатория LPTP физики и технологии плазмы Политехнической школы, Франция). Сравнение с известными источниками. В частности, решались задачи разгона и сжатия плазменных лайнеров, а также задачи сжатия магнитного потока (MFC). Рассмотренные постановки задач сжатия магнитного потока имеют практические приложения: определение уравнений состояния веществ методами безударного и ударного сжатия, создание мощных источников рентгеновского излучения, исследования в области управляемого термоядерного синтеза.
Разработка ряда вспомогательных модулей (например, модуль поддержки работы с разреженными матрицами, модуль ввода—вывода результатов и т.п.). Программная реализация соответствующих классов.
ВВЕДЕНИЕ
Страница 29
Публикации
По результатам диссертационной работы имеется 13 публикаций, список которых приведен в конце настоящего документа (на стр. 246).
Содержание работы
В настоящем введении описаны актуальность темы диссертации, сформулированы цели исследований, рассмотрены научная новизна и практическая значимость полученных результатов, перечислены основные положения, выносимые на защиту. Рассмотрены также вопросы достоверности и апробации результатов, описан личный вклад автора диссертации, приведен перечень научных публикаций по теме диссертации. Приведено развернутое описание содержания работы.
Первая глава посвящена предварительному сравнительному анализу ряда одномерных схем МГД, в частности, модифицированной схемы Лакса-Фридрихса [3],[14], схемы Курганова-Тадмора [16], а также схемы коррекции потоков на основе линеаризованной задачи Рішана [17]. Первая из указанных схем описана в разделе 1.1.1 рассматриваемой главы, вторая — в разделе 1.1.2, а третья — в разделе 1.1.3.
В разделе 1.1.4 осуществлено сравнение указанных схем. Его целью являлся выбор некоторой известной базовой одномерной МГД-схемы для дальнейшей ее модификации и обобщения на случай использования неструктурированных треугольных сеток. Для проведения сравнения схем использовался вычислительный эксперимент, заключавшийся в расчете ряда известных тестовых задач магнитной гидродинамики: «магнитный поршень» (подраздел 1.1.4.1), распад МГД-разрыва (подраздел 1.1.4.2), распространение ударной волны по холодному фону (подраздел 1.1.4.3) и другие задачи. Кроме того, при сравнении
ВВЕДЕНИЕ
Страница 30
схем учитывались известные достоинства и недостатки соответствующих им алгоритмов применительно к реализации на ЭВМ.
На основании проведенного сравнения сделан вывод о целесообразности использования одномерной модифицированной схемы Лакса-Фридрихса в качестве базы для обобщения на случай неструктурированных сеток. Кроме того, при разработке двумерной МГД-схемы использовались определенные свойства схемы Курганова-Тадмора, например, построение немонотонных и квазимонотонных реконструкций сеточных величин. Предложенная в настоящей работе разностная схема обеспечивает порядок аппроксимации 0(г2 + /г2) на однородных треугольных сетках, О (г2 + h) — на произвольных треугольных сетках и является устойчивой при выполнении условия Куранта. Обобщенная схема предназначена для использования в модуле расчета недиссипа-тивной двухтемпературной МГД, являющемся основным компонентом разрабатываемого МГД-солвера (MARPLE).
В главе 2 описаны вычислительные алгоритмы, реализованные в отдельных компонентах общего МГД-солвера, а также способ их взаимодействия друг с другом. Общий МГД-солвер MARPLE реализует известный принцип суммарной аппроксимации как на методологическом уровне, так и на уровне программной реализации.
Реализация принципа суммарной аппроксимации на методологическом уровне осуществляется разбиением членов уравнений системы дис-сипативной двухтемпературной МГД на группы, исходя из соображений физического и вычислительного характера. Отдельные группы затем аппроксимируются различными способами в различных вычислительных блоках. Затем результаты работы отдельных блоков объединяются в единое целое определенным образом, обеспечивая суммарную аппроксимацию полной системы уравнений.
Соблюдение принципа суммарной аппроксимации на уровне про-
ВВЕДЕНИЕ
Страница 31
граммной реализации основано на применении в процессе разработки современного объектно-ориентированного подхода. В соответствии с этим подходом осуществляется построение иерархической системы слабо взаимодействующих классов, инкапсуляция данных и методов внутри классов и т.д.
В разделе 2.1 приведена решаемая система уравнений двухтемпе-ратурной диссипативной МГД.
В разделе 2.2 предложен оригинальный алгоритм построения двойственных сеток (сеток контрольных объемов), необходимых, наряду с основными треугольными сетками, для расчета. Алгоритм основан на осуществлении коррекций известных диаграмм Вороного (Дирихле) с целью выполнения определенных условий. В разделе подробно описан алгоритм, показаны отличия получаемых с его помощью сеток контрольных объемов от известных аналогов (диаграммы Вороного, барицентрические контрольные объемы). Предложенный алгоритм программно реализован в MARPLE и успешно использовался при решении различных задач плазмодинамики.
Схема коррекции потоков на основе решения линеаризованной задачи Римана
Рассмотрим систему квазилинейных уравнений переноса: дї ut + їх = 0; f = f(u); A = —, (1.1.12) au. где u = (u\u2,...,un), f = f(u) - (Л/2,...,/"), (x,t) є R1 x {t Є [0, oo)}. Пусть на момент времени t = to заданы начальные условия и(х, t = to) = щ(х) и требуется найти u(x,t) при t to. Предположим также, что функция f (и) непрерывно дифференцируема, а матрица А размерности (п х п) имеет п действительных собственных чисел. Обозначим через Л диагональную матрицу собственных чисел матрицы А, через ft — матрицу левых собственных векторов (строк) матрицы А: А = Г2_1ЛП5. Наряду с системой (1.1.12) будем рассматривать также эквивалентную систему
В области решения введем равномерную сетку по переменным х и t: и = ujh х ит = {fa,tk), ХІ = ih, tk - to + jr, і є I, j = 0,1,2,...}. Конструкцию разностной схемы рассмотрим вначале при условии постоянства матрицы А. Умножив систему (1.1.13) на матрицу Г2, перейдем к системе уравнений для инвариантов Римана г = ft и:
Исходя из монотонной схемы первого порядка точности для уравнения (1.1.14) рассматривавшейся ранее в [18] и других работах, в [17] построена схема повышенной точности, обладающая свойствами монотонности и консервативности:
Функции ak± = a(Rk±) служат для ограничения антидиффузионных потоков [19], [3], Rk± — анализаторы гладкости решения. Подробный анализ свойств этих функций и требований, предъявляемых к их гладкости, проведен, например, в работах [21], [22].
Обобщение схемы (1.1.16) на нелинейные уравнения достигается переходом от инвариантов Римана г к исходным переменным и по формуле г = ft и и умножением уравнения (1.1.16) слева на матрицу S7-1. Один из возможных вариантов схемы для нелинейных уравнений, рассматривавшийся в [17] и использованный для численных экспериментов в данной работе, имеет вид
Для нелинейной системы уравнений магнитной газовой динамики используются анализаторы гладкости решения в виде, предложенном в
Дополнительные ограничения антидиффузионных потоков с целью предотвращения возможных нефизичных изменений энтропии учитывались на основе рекомендаций [22]. На сравнительных графиках, приведенных в разделе 1.1.4, видно, что данная схема дает самые качественные результаты из трех рассмотренных нами одномерных схем. Однако при этом она сравнительно сложна в реализации, требует весьма существенных затрат вычислительных ресурсов, и ее обобщение на случай треугольных сеток не является очевидным. Схемы, основанные на решении линейной задачи Римана, будем называть схемами Куранта-Изаксона-Риса (КИР) в соответствии с терминологией, принятой в [3].
Наряду с приведенными модификациями схем ЛФ и КИР в вычислительных экспериментах с одномерными схемами применялась также центрально-разностная схема КТ [16]. Вычисления в схеме КТ организованы по принципу «предиктор-корректор» и начинаются с построения сглаженной интерполяции разностного решения, построенной по узловым значениям решения гЦ, известным на текущий момент времени расчета t = tn: .22) где Xj(x) характеристическая функция сеточного интервала If
В кусочно-линейной интерполяции, выражаемой формулой (1.1.22), использованы сглаженные коэффициенты наклонов интерполирующих прямых u j. Они вычисляются по приращениям интерполируемой величины посредством функции MINMOD: u j = ММв {uj-uuj j+x) = ММ (в Дг +і/2, - (A -i/2 + А +1/2) , в Aiij_i/2),
Здесь в — безразмерный параметр: 9 Є [1,2], введено обозначение AWJ+I/2 = Uj+\ — Uj. В восполнении разностного решения для всей области изменения переменной х, построенном с помощью (1.1.22), на границах интервалов Ij, вообще говоря, могут быть разрывы. В схеме КТ вычисления на шаге tn — tn+1 производятся с учетом информации о движении этих возможных разрывов и, соответственно, об изменении участков гладкости восполненного разностного решения.
Расчет поля течения на момент времени t = tn+1 — tn + г состоит из следующих этапов. Вначале рассчитываются величины в узлах сетки х = Xj на промежуточном слое по времени t — tn+1/2:
В уравнении (1.1.27) vJ/j — u(xj,tn) = Uj, величины /j вычисляются с помощью A = (df/du) по формуле f- — A- u j, а величины и ,, в свою очередь, определяются выражением (1.1.24). Штрих ( ) обозначает, как и в предыдущих формулах, дифференцирование по пространству —. Найденные из (1.1.27) промежуточные величины позволяют со вторым порядком аппроксимации по пространству рассчитать те же величины на время t = tn+l в полуцелых точках сетки х = ж;-+і/2
Завершающая стадия вычислений на шаге tn — tn+l состоит в расчете значений искомых величин в целых узлах сетки х = Xj. С этой целью снова строится восполнение разностного решения по его значениям г +і/2 в полу целых точках х — 2 +1/2, вычисленным по формуле (1.1.28):
Вывод уравнений диффузии магнитного поля в случае анизотропной проводимости сгц ф а± .
В данном разделе осуществляется вывод уравнений диффузии поля при анизотропной проводимости плазмы. Под анизотропностью проводимости а подразумевается, что значение проводимости в направлении вдоль вектора напряженности магнитного поля И. ( сгц) отличается от значения проводимости в перпендикулярном И направлении ( 7j_). Рассматривается общий трехмерный случай, т.е. векторы плотности тока, напряженности электрического поля, напряженности магнитного поля имеют три компоненты. Напомним, что частный случай изотропной проводимости описан в разделе 2.3.1 данной работы.
Пусть Е — вектор напряженности электрического поля в лабора —) торной (покоящейся) системе координат (далее — Л.С.К.), Е — вектор напряженности электрического поля в системе координат, движущейся вместе с плазмой (далее — Д.С.К.), Н — вектор напряженности магнитного поля в Л.С.К., j — вектор плотности тока в Л.С.К. Далее в данном разделе подразумевается, что величины с индексом « » относятся к движущейся вместе с плазмой Д.С.К., а величины без « » — к покоящейся Л.С.К.
Рассмотрим в трехмерном пространстве плоскость, образованную векторами Е и п. Поскольку вектор Е лежит в данной плоскости, его можно разложить в базисе, состоящем из двух векторов, лежащих в этой плоскости (Е , Н): вектора, параллельного Я, и вектора, перпендикулярного Й. Пусть для определенности направление второго из указанных векторов (v±) задается следующим образом:
Используя известное соотношение между векторами напряженности электрического поля в Л.С.К. (Е) И Д.С.К. (Е ) (ги - скорость плазмы в рассматриваемой точке пространства, членами порядка ги2/с2 и далее пренебрегаем) и выведенное выше выражение (2.3.21), получим:
Уравнения (2.3.17) и (2.3.23) определяют (при известном профиле скорости плазмы w = w(r,t)) эволюцию магнитного и электрического полей в пространстве с течением времени в случае с анизотропной проводимостью а = а. Из уравнений (2.3.17), (2.3.23) можно исключить электрическое поле, в этом случае диффузия магнитного поля описывается одним уравнением:
Для вычисления диффузии магнитного поля в плазме можно решать либо одно дифференциальное уравнение второго порядка (2.3.24), либо систему из двух дифференциальных уравнений первого порядка (2.3.17), (2.3.23). В данной работе мы остановились на первом способе — решении (2.3.24). Однако следует отметить, что выведенные уравнения (2.3.21) для Е и (2.3.23) для Е все равно используются в расче те. Это связано с тем, что в процессе расчета различных задач возникает необходимость в вычислении электрических полей в плазме. Однако векторы напряженности электрического поля Е в Д.С.К. и Е в Л.С.К. не фигурируют ни в каких уравнениях в качестве независимых переменных, относительно которых разрешаются уравнения. Т.е. указанные векторы можно рассматривать просто как некоторые известные дифференциальные операторы вида Е = Е (н,а±,(тА,
Е = Е I И, ох, сц j, значения которых вычисляются (точнее — аппроксимируются на расчетной сетке) по ходу расчета в случае необходимости.
Следует отметить еще и следующее. В случае идеальной МГД (а± = сгц -» со) уравнения (2.3.21), (2.3.23), (2.3.24) принимают следующий вид:
В этом идеальном случае диффузия магнитного поля, как таковая, отсутствует, а магнитное поле «вморожено» в плазму и перемещается в пространстве и времени в соответствии с уравнением (2.3.27). Легко заметить, что это уравнение является составной частью системы уравнений идеальной магнитной гидродинамики, см. раздел 2.4. Соответствующий (2.3.27) член обобщенного уравнения (2.3.24), таким образом, отвечает за перемещение магнитного поля по «вмо-роженности» в плазму, а не за диссипацию.
Квазимонотонная реконструкция на этапе «корректор»
Таким образом, в предлагаемом алгоритме расчета диффузии магнитного поля вместо одного обращения полной диффузионной матрицы М размерности 3m х Зт осуществляется многократное обращение трех ее матриц-подблоков Мхх, Муу, Mzz размерности тхт. Многократные итерации требуются для учета вклада в решение недиагональных подблоков MXY, Mxz, Мух, Myz, Mzx, Mzy. Способ осуществления итераций похож на метод Зейделя, где в качестве неизвестных фигурируют векторы размерности т, а в качестве коэффициентов при них — матрицы размерности тхт.
Экономия вычислительных ресурсов осуществляется за счет отказа от обращения полной диффузионной матрицы размерности 3m х Зт. Поскольку вычислительные затраты на обращение матрицы нелинейно возрастают при увеличении ее размерности, то даже многократное обращение матриц размерности тхт требует меньших затрат, чем обращение полной диффузионной матрицы. Заметим также, что высокая точность обращения диффузионной матрицы совсем не требуется. В частности, при решении различных описанных в данной работе задач магнитной гидродинамики установленная требуемая точность Єо обращения диффузионной матрицы (оцениваемая по бесконечной норме невязки L( о )) составляла 5 10 5.
Еще одно преимущество работы не с полной диффузионной матрицей, а с ее подблоками, заключается в следующем. Во многих практически важных задачах структура магнитного поля имеет заранее известные особенности. Например, в задаче о разгоне плазменного лайнера азимутальным (полоидальным) магнитным полем вектор Н имеет только одну ненулевую компоненту — Нр, в то время как Нц и Hz равны нулю (рассматриваемая задача решается в цилиндрической геометрии (R,(p,Z)). Причем эти условия выполняются в течение всего расчета, т.е. при t 0 (это следует из анализа решаемых уравнений, начальных и краевых условий задачи). Соответственно, в задаче о разгоне лайнера в используемой нами постановке достаточно работать только с одной матрицей-подблоком Mw — операции с полной диффузионной матрицей в этом случае не имеют смысла. Итак, декомпозиция полной диффузионной матрицы М на подблоки и, соответственно, замена операций с этой матрицей на операции с отдельными подблоками позволяют лучше учесть особенности структуры магнитного поля некоторых решаемых задач. В частности, это помогает избежать выполнения лишних вычислительных операций.
Для осуществления шагов IV, V, VI рассмотренного выше алгоритма требуется обратить разреженные матрицы размерности тхт. В данной работе использован один из вариантов итерационного метода сопряженных градиентов (preconditioned bi-conjugate gradient stabilized method). Для предварительной обработки (preconditioning) матрицы применяется неполное LU-разложение. Использована программная реализация метода сопряженных градиентов, входящая в состав свободно распространяемого пакета Iterative Methods Library [29] / Sparse Lib [30], созданного в Национальном институте стандартов и технологий (США).
В данном разделе выполнено построение алгоритма, осуществляющего один шаг (tn — tn+i = tn + At) неявной разностной схемы расчета диффузии магнитного поля. Предлагаемая неявная схема обеспечивает порядок аппроксимации 0(h2) на однородных треугольных сетках, 0(h) — на произвольных треугольных сетках. Заметим, что модуль расчета диффузии магнитного поля используется совместно с модуля ми расчета идеальной магнитной гидродинамики, теплопроводности и иошю-электронного обмена, джоулева нагрева и излучения. Поэтому шаги построенной схемы специальным образом чередуются с шагами вычислительных схем из других модулей с целью обеспечения суммарной аппроксимации полной системы решаемых уравнений. Подробнее о построении общей вычислительной схемы см. в разделе 2.8.
Выведенное ранее выражение (2.3.65) позволяет аппроксимировать величину Е в серединах ребер двойственной сетки. В некоторых случаях (например, при определенных способах вычисления джоулева нагрева) требуется аппроксимировать эту величину в узлах треугольной сетки. Пусть нужно найти Е\ — аппроксимацию Е в узле і.
Разгон невозмущенного плазменного лайнера
В приведенных выше уравнениях CV(i) — множество всех контрольных объемов, смежных по ребру (грани) с i-ым контрольным объемом (ij — соответствующее ребро (грань) двойственной сетки); Vi = Si L — объем г -ой призматической ячейки двойственной сетки (трехмерного контрольного объема); Si — площадь г -ой плоской ячейки двойственной сетки (двумерного контрольного объема). О вычислении потоков Fij написано ниже в подразделе 2.4.2.5.
На рассматриваемом этапе «предиктор» для восполнения сеточного решения (и других сеточных функций) на временном слое t = tn используется немонотонная реконструкция (2.4.20), описанная в подразделе 2.4.2.1. пе «предиктор»; г-ая ячейка двойственной сетки (контрольный объем) выделена жирными линиями. Треугольники, па которые она разбивается, пронумерованы. Каоїсдая область постоянства интерполяционных коэффициентов (в пределах рассматриваемого контрольного объема) выделена своей штриховкой.
Наличие однозначно реконструированного (на текущий момент времени t = tn) по своим узловым значениям решения позволяет однозначно аппроксимировать интегралы недивергентных членов (w VPe) и (w VPj) из уравнений (2.4.47) и (2.4.48) соответственно. Напомним, что эти члены относятся к n-му слою по времени, хотя верхний индекс у них опущен. Пусть требуется аппроксимировать выражение / (w VP&) dV, Vi где к 6 {є, і} позволяет выбрать электронную или ионную компоненту давления, і — номер рассматриваемого узла основной сетки (ему соответствует контрольный объем Vi). Разобьем контрольный объем на треугольники в соответствии с рис. 2.4.3.
Для всех внутренних точек каждого из треугольников (ЛТт), выделенных и пронумерованных на рис. 2.4.3, в соответствии с используемой реконструкцией (2.4.20) можно записать: !$ = -? (х- xnode{r)) +t(y- ynode{r)) + , (2.4.49) Рк = аргк (х - xnode{r)) + 6ffc (у - ynode(r)) + cffc, (2.4.50) где г = г(АТт) — однозначно определяемая треугольная ячейка основной сетки, подобластью которой является рассматриваемый треугольник (ЛТт С Лг), а коэффициенты ar, br, сг определяются в соответствии с формулой (2.4.22). Очевидно, что указанные коэффициенты не меняются в пределах рассматриваемого ЛТт. В соответствии с (2.4.49), (2.4.50), получаем: (2.4.51)
Здесь г = г(АТт) — однозначно определяемая треугольная ячейка основной сетки, подобластью которой является рассматриваемый треугольник АТт; V(ATm) = S(ATm) L — объем трехмерной призматической ячейки, основанием которой является АТт (S(ATm) — площадь этого треугольника, L — высота призматических ячеек); (я?центр(дтт) 2/цеіітр(лтт)) — координаты геометрического центра треугольника АТт, определяемые по формуле
На данном этапе осуществляется переход от временного слоя t = tn к временному слою t = tn+\ = tn + At, соответственно во всех вычислениях используется полный шаг по времени At. Величины (без верхнего индекса), фигурирующие в правых частях приведенных ниже разностных уравнений, берутся с промежуточного (п + 1/2)-го временного слоя t = tn + At/2, рассчитанного на этапе «предиктор». Верхний индекс у величин соответствует временному слою, а нижний — узлу треугольной сетки. Разностные уравнения на данном этапе имеют следующий вид:
В приведенных выше уравнениях CV(i) — множество всех контрольных объемов, смежных по ребру (грани) с г-ым контрольным объемом (ij — соответствующее ребро (грань) двойственной сетки); Vi = Si L — объем г -ой призматической ячейки двойственной сетки (трехмерного контрольного объема); Si — площадь г -ой плоской ячейки двойственной сетки (двумерного контрольного объема). О вычислении потоков Fij написано ниже в подразделе 2.4.2.5.
На рассматриваемом этапе «корректор» для восполнения сеточного решения (и других сеточных функций) на промежуточном временном слое t = tn+i/2 = tn + At/2 используется квазимонотонная реконструкция (2.4.27), описанная в подразделе 2.4.2.2.
Наличие однозначно реконструированного (на промежуточный момент времени t = п+і/г) по своим узловым значениям решения позволяет однозначно аппроксимировать интегралы недивергентных членов (wVPe) и (wVPi) из уравнений (2.4.61) и (2.4.62) соответственно.
