Содержание к диссертации
Введение
ГЛАВА 1. Расчет динамики высокотемпературной плазмы в приближении ЭМГ 22
Физическая модель 22
Математическая модель 23
Обезразмеривание и основные масштабы величин 25
Численный метод расчета 27
ГЛАВА II. Описание численного метода и результаты тестовых расчетов 39
Результаты тестирования 42
Тестирование модулей комплекса 42
Описание тестов 44
ГЛАВА III Результаты вычислительных экспериментов 54
Исследование сжатия тяжелого Z-пинча 67
Заключение 71
Список использованных источников
- Математическая модель
- Обезразмеривание и основные масштабы величин
- Тестирование модулей комплекса
- Исследование сжатия тяжелого Z-пинча
Введение к работе
Актуальность темы. Математическое моделирование как научное направление еще достаточно молодо, и его активное развитие продолжается. Основные черты современного математического моделирования связаны с тем, что в последнее десятилетие оно быстро теряет «академические» черты чисто научного и узкопрофессионального направления. Наряду с классическими задачами в настоящее время решаются многочисленные проблемы, возникающие при практическом использовании методов и результатов моделирования.
Методология математического моделирования в кратком виде выражена знаменитой триадой «модель — алгоритм — программа», сформулированной академиком А. А, Самарским [1], одним из «отцов-основателей» отечественного математического моделирования. Эта методология получила свое развитие в виде технологии «вычислительного эксперимента», разработанной школой А. А, Самарского, — одной из информационных технологий, предназначенной для изучения явлений окружающего мира, когда натурный эксперимент оказывается слишком дорогим и сложным.
Хотя традиции применения математического аппарата и математического моделирования насчитывают сотни лет, первые серьезные работы в области моделирования процессов в плазме относятся к 60-ым годам XX века.
Дата появления первых серьезных результатов вычислительного
эксперимента в СССР зафиксирована вполне официально — 1968 год, когда
Госкомитет СССР по делам открытий и изобретений засвидетельствовал
открытие явления, которого на самом деле никто не наблюдал. Это было
открытие, так называемого, эффекта Т-слоя (температурного токового слоя в
плазме, которая образуется в МГД-генераторах) [I], Свидетельство на это
открытие было выдано академикам А. Н. Тихонову и А. А. Самарскому,
члену-корреспонденту АН СССР С. П. Курдюмову, докторам физико-
математических наук П. П. Волосевичу, Л. М. Дегтяреву,
Л. А. Заклязьминскому, Ю. П. Попову, В. С. Соколову и А. П. Фаворскому
И-
В данном случае вычислительный эксперимент предшествовал
натурному. Натурные эксперименты «заказывались» по результатам
математического моделирования. Через несколько лет в трех физических
лабораториях на разных экспериментальных установках практически
одновременно был надежно зарегистрирован Т-слой, после чего технологам
и инженерам стал окончательно ясен принцип работы МГД-генератора с Т-
слоем.
В настоящее время плазма с ее нелинейными свойствами стала одним
из важнейших объектов математического моделирования и вычислительного
эксперимента. В качестве исследуемых объектов, представляющих интерес,
как с научной, так и с практической точки зрения, можно отметить процессы,
протекающие в плазме в электромагнитном поле при сверхвысоких
температурах. Экспериментальные исследования проводятся в нескольких
лабораториях различных стран, в том числе на крупнейших установках
России - Ангара-5-1 и С-300, в Великобритании на установке MAGPIE и в
США на крупнейшей в мире установке «Z».
Сейчас установки по исследованию динамики высокотемпературной плазмы представляют собой сложные системы. В качестве частей таких установок применяются плазменные прерыватели тока (ППТ) и плазменные размыкатели. Физические процессы, происходящие в ППТ и размыкателях, представляют самостоятельный интерес.
Для комплексного исследования динамики плазмы на подобных установках не достаточно только натурного эксперимента — он не дает полной и подробной картины происходящего. Конструктивные особенности диагностических частей позволяют измерять характеристики плазмы лишь в отдельные дискретные моменты времени. Таким образом, для получения полной информации о протекающих процессах, актуальным является подход, сочетающий проведение натурных и численных экспериментов.
В последние годы мощность современных компьютеров и кластерных систем значительно увеличилась, вследствие чего, стало возможным применять более ресурсоемкие с вычислительной точки зрения алгоритмы. В основе таких алгоритмов лежат сложные математические модели, включающие описания тонких эффектов, в частности, эффекта Холла, аномального сопротивления и т.д.
Цели и задачи диссертации Основными целями и задачами диссертационной работы являются:
Построение эффективного численного метода для моделирования плазмы в приближении электронной магнитной гидродинамики (ЭМГ) на основе неявной схемы.
Разработка соответствующих численных алгоритмов.
Разработка комплекса программ, реализующих численный метод.
4. Проведение исследования динамики поведения плазмы при
параметрах модели, соответствующих экспериментальной установке
«Стенд-300»; исследование динамики плазмы в ПИТ и Z-пинчах.
Научная новизна
Для расчета динамики высокотемпературной плазмы в приближении ЭМГ применен неявный вариационный метод.
Исследована динамика плазменной перетяжки в ПИТ, в частности, численно исследован эффект проникновения магнитного поля в материал плазмы с образованием пузырей с магнитным полем.
Численно исследован характер неоднородностей, образующихся в плазме при сжатии медного Z-пинча в рамках приближения ЭМГ.
Научная и практическая ценность Разработан эффективный численный метод для моделирования высокотемпературной плазмы в приближении ЭМГ. Данный метод был реализован в виде комплекса программ ЭВМ. Полученные результаты численного моделирования расширяют представления о физических явлениях, происходящих в высокотемпературной плазме на различных экспериментальных установках» в частности С-300. Они позволяют численно воспроизвести ряд явлений, предсказанных теоретически (в частности, образование плазменных пузырей) и картину их возникновения.
Разработанный комплекс программ также может применяться для моделирования широкого круга научно-практических задач, связанных моделированием высокотемпературной плазмы.
Апробация работы Материалы, отражающие содержание диссертационной работы докладывались на двух международных конференциях (11-я Международная конференция «Математика, компьютер, образование» — Дубна, 2004 г, 15-я Крымская осенняя математическая школа —2004 г), на 2-ой Курчатовской молодежной научной школе — Москва, 2004 и на XLVII научной конференции Московского физико-технического института— Долгопрудный, 2004 г. Результаты также докладывались на научных семинарах кафедр информатики и вычислительной математики МФТИ и отделения прикладной физики Института ядерного синтеза РНЦ «Курчатовский институт».
Публикации По теме диссертации опубликовано 6 работ [3 — 8].
Благодарности, гранты Проведенная работа выполнена при поддержке РФФИ (проект №04-01-00774). Автор выражает огромную благодарность А. С. Кингсепу, К- В, Чукбару, Ю- Г- Калинину за постановку задачи и полезные обсуждения, А. В. Богатову за техническое содействие- Особую
благодарность автор выражает А. В. Зырянову за предоставленную возможность использования вычислительного оборудования.
Математическая модель
В приведенных выше формулах индекс є относится к электронам, а /—- к ионам. Большинство обозначений традиционные Р без индекса означает полное давление — сумму парциальных давлений электронов и ионов, функция J{zeff) учитывает потери энергии на ионизацию, zejf — эффективный заряд иона. В — напряженность магнитного поля, j — плотность тока, ve — электронная токовая скорость, а — проводимость, пе — концентрация электронов, Те — электронная температура, Q — джоулев нагрев. В дальнейшем считаем, что магнитное поле имеет лишь азимутальную компоненту, ЭМГ-эффекты вносят дополнительную нелинейность в систему уравнений. Приведем соотношения, необходимые для замыкания системы. Уравнения состояния: е »Р А Р. А» 1 00 =АЛ (8) РІ=ЛРЩ (9 Ъ=А Тп (10) здесь и далее А — число нуклонов в ядре, Ае = 14,38; Ар = 9,68 10 2,
Зависимость эффективного заряда иона от электронной температуры дается выражением zeJf =9T/3,zef Z, иначе zeJf=Z, где 2— максимально возможный заряд иона для данного элемента (материала). Для рассматриваемого далее материала — углерода — Z принимает значение 6. Для обмена энергиями между электронной и ионной компонентами плазмы имеем: A f73 є Для описания потерь на ионизацию применяется скейлинг гг ч 0Л31 0.85 + 0.15Z/3 г 3 г г г 2g/ Zeff l Zeff К6 2 6 ; Граничные условия На правой и левой границах области ставились условия жесткой стенки (условия непротекания).
На нижней и верхней границе ставились условия границы с вакуумом. При этом азимутальное магнитное поле на нижней границе связано с полным током во внешней электрической цепи известным соотношением, записанным в безразмерных переменных [73] 0.2/(0 „х г . «ti —Ч /0) = /0 sin — где /— ток во внешней цепи, R — текущий радиус соответствующей точки поверхности плазмы. В случае учета радиационных потерь на этой границе ставилось также условие отсутствия падающего извне излучения.
Обезразмеривание и основные масштабы величин
Для анализа получаемых уравнений и организации устойчивых численных расчетов необходимо перейти к безразмерным переменым. В таблице 1Л приводятся значения констант, входящих в уравнения при их записи в безразмерных единицах и в единицах системы СГС [74] Для численного решения нелинейной системы применяется консервативная разностная схема на криволинейной подвижной сетке [75] с использованием расщепления по физическим процессам — стандартным методом для решения сложных задач МСС [76]. На первом этапе расчета (модуль МГД) используются уравнения модели движения двухтемпературнои плазмы в магнитном поле в приближении МГД аналогичные [71]. Расчет эволюции магнитного поля проводится на «замороженной» на данном шаге по времени сетке с применением ЭМГ системы [52], в качестве начальных данных используются результаты модуля МГД.
На следующем этапе расчета учитываются изменения электронной и ионной температур за счет теплопроводности. На заключительном этапе расщепления проводится коррекция значений температуры, учитывается обмен энергиями за счет электрон-ионных столкновений [71]. Далее процесс расчета повторяется на следующем шаге по времени. Для расчета МГД течений характеризующихся большими сдвиговыми деформациями чаще всего применяются разностные методы, основанные на явной аппроксимации уравнений на подвижной (лагранжевой) сетке [75], При расчете по явным схемам приходится использовать малый шаг по времени т — это следует из условия устойчивости явных схем с% к=— 1, где h — характерный размер ячейки сетки. h Неявные схемы формально свободны от ограничений такого рода.
Для построения неявного метода (модуль МГД) использовались: система уравнений двухтемпературной магнитной гидродинамики (МГД) [52,71], уравнения состояния и условие «вмороженности» магнитного поля, В рассматриваемой области, занятой плазмой, была введена подвижная лагранжева сетка. В начальный момент времени граница сетки прямоугольна и совпадает с границей области, занятой плазмой. Разностная схема во внутренних точках сетки имеет второй порядок аппроксимации по пространственным переменным. Для уменьшения осцилляции, возникающих при расчете разрывных решений (ударных волн, контактных разрывов или тангенциальных разрывов), к давлению добавляется искусственная вязкость. В расчетах применяется комбинация линейной и квадратичной вязкости [77], Разностная схема строилась следующим образом.
Обезразмеривание и основные масштабы величин
Система уравнений математической модели (1 —6) является нелинейной. Непосредстве ное численное решение такой системы с помощью неявных разностных схем при существующем уровне развития вычислительной техники не представляется возможным. Для решения подобных задач обычно применяются методы операторного расщепления (суммарной аппроксимации, расщепления по физическим процессам). Основные идеи методов расщепления были предложены Н. Н. Янеико [80], для задач вычислительной гидродинамики методы расщепления развивались научной школой О. М. Белоцерковского [76].
Кроме того, опыт решения задач двухтемпературной гидродинамики в условиях, соответствующих экспериментальным для установки С-300. показывает, что ионное давление является более консервативной величиной и на малых временах меняется значительно слабее, чем электронное [52].
Таким образом, физическая природа задачи диктует выбор слеующей схемы расщепления. На каждом шаге по времени проводится операторное расщепление, решение системы уравнений сводится к последовательности более простых задач. В рамках решения сеточных уравнений, аппроксимирующих систему уравнений МГД, проводится расщепление по давлениям. Сначала решается система нелинейных разностных уравнений для определения ионного давления на новом временном слое, затем с учетом вычисленного значения ионного давления решается система уравнений для определения электронного давления и значения магнитной индукции с учетом вмороженности поля на данном этапе расчета.
Остановимся на описании последовательности расчетов подробнее. Расчет одного шага по времени в модуле МГД осуществляется следующим образом: 1. С помощью метода Ньютона рассчитывается поля давлений и значение вектора магнитной индукции; Обновляются значения на границах и в фиктивных ячейках; 2. По новым давлениям обновляются поля скоростей, координаты узлов, площади и объемы ячеек, плотности, температур; 3. Вычисляются энергии, степень ионизации, потери на ионизацию.
В программе расчета заложены ограничения на минимальный размер ячейки разностной сетки, В том случае, когда пространственный размер какой-либо ячейки становится меньше допустимого, проводится консервативный пересчет всех величин на новую сетку. При этом пересчитываются не только «традиционные» величины — температуры, количество движения, масса плазмы, но и средний заряд иона в ячейке, потери на ионизацию.
На следующих этапах решаются уравнения электронной и ионной теплопроводности, а также диффузии магнитного поля. Для этого применяются схемы потоковой прогонки [81],
На заключительном этапе учитываются обмены энергиями между электронами и ионами и потери на излучение. Жесткие системы ОДУ, возникающие на данном этапе расщепления, для каждой ячейки решаются методом трапеций. На рис. 4 представлена блок-схема программного комплекса для исследования динамики плазмы.
Применяемая неявная схема, основанная на вариационной методике, позволила снять существенные ограничения на шаги по времени, диктуемые условиями устойчивости. Похожие неявные методы ранее реализовывались для задач однотемпературной МГД [57,78]. Характерной особенностью предлагаемой расчетной методики является расщепление по ионному и электронному давлениям при решении системы уравнений для давления на следующем временном слое. Модуль расчетов динамики плазмы на основе неявной схемы реализован в виде нескольких файлов на языке FORTRAN77 для совместного использования с имеющимися на кафедре программами.
Результаты тестирования
Для проверки реализации модуля МГД применялась система тестов. Процедура тестирования вычислительного комплекса проводилась в следующем порядке; сначала тестировались отдельные блоки системы, а именно отдельные функции, далее проводилось тестирование составных блоков, то есть модулей системы, после этого проводилось тестирование всей системы в целом. Тестирование модулей комплекса Модуль расчета магнитной гидродинамики плазмы (МГД), основанный на неявной схеме.
В данном модуле для численного решения нелинейной системы применяется консервативная разностная схема на криволинейной подвижной сетке [75] с использованием расщепления по физическим процессам [82].
Ранее использовался модуль МГД, основанный на явной разностной схеме. Использование этого модуля было лимитировано соответствующими ограничениями на параметры расчета для устранения численной неустойчивости.
Тестирование модулей комплекса
Численное моделирование динамики плазмы на Стенде - 300 проводилось при различных значениях параметров плазмы. Приводимые ниже результаты расчетов получены для следующих вариантов модели: радиусы электродов Rj = 0,2 см, R2 - 0,5 см, начальная толщина оболочки z2 0,075 см, начальная плотность вещества плазмы 0,3 мг/см , материалом плазмы считался углерод. Длительность импульса электрического тока Гітр = 100 не, амплитуда импульса тока /0 = —1 МА, По предварительным качественным оценкам параметр замагниченности в условиях, соответствующих установке С - 300, находятся в пределах 0,01—1,0.
В работе также учитывается неравномерное начальное распределение температуры плазмы Te=Tt= Т(г). Считается, что в результате прохождения предварительного импульса плазма нагревалась. Так как плотность тока выше в окрестности катода, то начальная температура в этой области выше. Будем считать, что начальное распределение температуры можно описать распределением: T(r) = T0+ATexV{-(r Rl)2/L2), где Лі — радиус внутреннего электрода, L — ширина зоны прогрева. Представленные результаты получены при значениях Г0 = 2эВ, AT— 8,5 эВ, L = (R2-Ri)/2.
По данным численных расчетов [52], динамика плазменной перемычки существенно зависит от начального прогрева плазмы, то есть от начального распределения температур.
Так как в расчетах использована подвижная перестраеваемая сетка, по динамике сетки можно восстановить эволюцию плазменной перемычки. Вначале плазма разлеталась под действием газокинетического давления- По мере роста тока генератора начиналось подтормаживание нижней границы, изменялось направление движения перемычки. Плазма приходила в движение при времени 17-20 не. По мере возрастания тока и нагрева в приэлектродной области по материалу перемычки начинала распространяться ударная волна. На рисунках 12а и 12Ь показано изменениеДля таких режимов сжатия характерен эффект «снежного плуга», когда вещество концентрируется в достаточно узкой области [84], Эффект «снежного плуга» виден на рис. 13а и 13Ь, где приведены распределения плотности (концентрации) в соответствующие моменты времени.
На рис. 14а и 14Ь для того лее варианта расчета приведены распределения электронной и ионной температур. В области «снежного плуга» они практически равны, наблюдается незначительное преобладание электронной температуры. Это связано с джоулевым нагревом.
Вблизи внутреннего электрода (анода) плотность тока максимальна, соответственно, там же наблюдается и максимальное тепловыделение за счет джоулева нагрева. Но в силу большого значения магнитной индукции вблизи анода формируется сильная ударная волна. Из-за наличия в разностной схеме искусственной вязкости [77], фронт ударной волны размазан на несколько ячеек разностной сетки. Формирование ударной волны можно проследить по превышению значений ионной температуры над электронной (рис. На и 14Ь). Так, для 40 не. Максимальное значение электронной температуры в области за ударной волной составляет 18,5 эВ, в то время как для ионов температура достигает 27,5 эВ. К 50 не эта разность увеличивается, и температура электронов достигает 36,4 эВ, а ионов — 86 эВ.
При времени около 5 0 не существенно изменяется характер проникновения магнитного поля в материал плазмы- Если до этого времени проникновение поля носило диффузионный характер, то в связи с повышением электронной температуры вблизи анода магнитное поле оказывается вытесненным из плазмы (рис, 15), В силу нелинейных эффектов на границе ударной волны вблизи границы с вакуумом формируются токовые петли. После формирования таких петель (а в наших расчетах такие петли являются устойчивыми) магнитное поле начинает распространяться от анода вдоль нижней границы. Видно, что такое распространение уже не носит диффузионный характер, а представляет собой нелинейную волну. Пока остается открытым вопрос о немонотонности магнитного поля вблизи фронта этой волны, возможно, она носит численный характер.
Оценки характерных размеров неоднородности в этом случае указывают на то, что может наблюдаться режим пространственного разделения зарядов, и в окрестности фронта нелинейной волны сформулированная выше математическая модель уже не описывает динамику плазмы.
Исследование сжатия тяжелого Z-пинча
Эволюция плазменной перемычки показана на рис. 16. Мелкомасштабные эффекты при проникновении магнитного поля в материал плазмы и джоулев нагрев в этом случае приводят к формированию пузырей малой плотности с магнитным полем внутри. Теоретически такая возможность проникновения магнитного поля в материал была впервые рассмотрена Л. И. Рудаковым с соавторами [86,87]. В последнее время появились экспериментальные указания на то, что такой характер проникновения возможен, по крайней мере, для Z-пинчей [88],
Отметим, что в случае формирования пузыря расчетная область становится неодносвязной, что влечет за собой невозможность продолжения счета по используемым разностным схемам на подвижной сетке.
Отметим, что проникновение в материал магнитного поля с образованием малоплотных пузырей приводит к формированию более слабых ударных волн в перемычке. В частности, это видно на рис. 18, где показаны распределения электронной и ионной температур. В то время как характерные значения электронной температуры остались примерно теми же, что и в предыдущем варианте, значение ионной температуры уменьшилось в полтора раза (58 эВ).
Как и в предыдущем расчете, вещество скапливалось на нижней границе вдали от анода из-за эффекта снежного плуга (рис. 17). Для вариантов с образованием пузырей так же характерна смена типа проникновения магнитного поля в материал плазмы с диффузионного на нелинейную волну. Рассмотрение направления проникновения показывает, что проникновение магнитного поля соответствует направлению распространения волны КМЧ [18].
При дальнейшем увеличении значения параметра замагниченности характер проникновения магнитного поля в материал плазмы вновь менялся, малоплотный пузырь с магнитным полем не формировался. Отличия от вариантов с малой замагниченностыо заключались в изменении амплитуд ударных волн и появлению мелкомасштабных образований в материале плазмы.
В качестве приложений данной программы рассмотрена задача о сжатии тяжелого медного Z-пинча. Основной целью исследования являлось изучение влияния мелкомасштабных приэлектродных ЭМГ-эффектов на однородность сжатия. Такая оценка необходима для аккуратного исследования его излучательных свойств.
Приводимые ниже результаты расчетов получены для следующих начальных данных: внутренний и внешний радиусы Д/ = 0,05 см, R2= 1,6 см, длина пинча Z2 = 2 см, полная масса 740 микрограмм (погонная плотность 340мкг/см), начальная плотность вещества распределена по следующему закону р(г) = р0 + Дрехр(-С- )2/і2), где плотность р0 10" ; масса пинча сосредоточена в окрестности К= 1,2 см; /, = 0,3.
Материалом Z-пинча считалась медь. Время наростання электрического тока Т\тр — 100 не, амплитуда импульса тока /0 = -6 МА. Начальные электронная и ионная температуры распределены равномерно Те — Т,\— 1 эВ.
По распределению плотности рис. 20 виден фронт ударной волны, идущий по заполнению. Магнитное поле экранируется перегретой областью за фронтом ударной волны. В прианодной области видно проникновение магнитного поля с образованием узкого языка. Такой характер проникновения объясняется вмороженностью обобщенного импульса в электронное течение и является специфическим эффектом ЭМГ. Повышение плотности тока в окрестности языка обычно приводит к резкому росту электронной температуры в заполнении и эффектам анодного взрыва.
Разница между ионной и электронной температурами рис. 22 на порядок за фронтом распространяющейся ударной волны характерна для таких течений. Распространение электронной температуры за счет струи направленной по заполнению к оси пинча приводит к образованию горячих точек в момент пинчевания и существенному нарушению однородности сжатия.
