автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Численное моделирование нелинейной трехмерной МГД эволюции тороидальной плазмы

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ измени М.В. ЛОМОНОСОВА Факультет вычислительной математики и кибернетики кафедра автоматизации научных исследований

Кандидатская диссертация УДК 519.6:533.9

Лю Юсцян

Численное моделирование нелинейной трехмерной МГД эволюции тороидальной плазмы

Специальность 05.13.16 - применение вычислительной техники, математического моделирования и математических методов

в научных исследованиях

Научный руководитель доктор физико-математических наук профессор А.М.Попов

Москва - 1999

Содержание

Введение 5

Глава 1 Трехмерные МГД модели для описания нелинейной эволюции тороидальной плазмы. Постановка математической задачи , 25

§ 1. Трехмерные МГД модели для описания нелинейной эволюции плазмы в тороидальной геометрии............ 25

п. 1. Модель устойчивости стационарного состояния плазмы с быстрым тороидальным вращением....... 29

п.2. Модель нелинейной устойчивости плазмы при наличии внешнего винтового магнитного поля...... 30

п.З. Учет влияния неоклассических эффектов на устойчивость тиринг-моды. Роль бутстрэп-тока...... 31

п.4. Самосогласованная модель для описания взаимодействия транспортных процессов и диссипативной

МГД неустойчивости плазмы ............. 34

§ 2. Масштабы МГД модели плазмы в токамаке и оценки безразмерные параметры ...................... 36

§ 3. Исходное равновесие и потоковая система координат. Запись дифференциальных операторов МГД уравнений в произвольной криволинейной системе координат ........ 37

§ 4. Постановка задачи об устойчивости и нелинейной эволюции плазмы в тороидальной геометрии............ 41

Глава 2 Разработка и исследование численных алгоритмов решения задачи о МГД устойчивости плазмы в произвольной тороидальной геометрии. Структура и тестирование нелинейного трехмерного кода МГТС 45

§ 1. Полуспектральный метод решения системы МГД уравнений в криволинейной системе координат........... 45

§ 2. Запись и формирование сверток................ 47

§ 3. Потоковая разностная аппроксимация по радиальному направлению ............................ 49

§ 4. Разностная аппроксимация по времени. Блочный метод Гаусса-Зейделя-Ньютона для решения нелинейных уравнений ................................. 50

§ 5. Метод обратной итерации для нахождения собственных значений краевой задачи для полной системы МГД уравнений ............................ ... 53

§ 6. Вычисление комплексных собственных значений в задачах

устойчивости вращающейся плазмы.............. 54

§ 7. Разработка и структура кода КИГС ............. 56

§ 8. Исследование численных свойств кода ЫРТС при решении

нелинейной МГД задачи !................... 58

§ 9. Тестирование численных результатов исследования устойчивости с использованием кода №ТС ............ 69

§ 10. Тестирование алгоритмов расчета осесимметричного равновесия для вращающейся плазмы с точными аналитическими решениями уравнения Машке-Перрина........ 73

Глава 3 Численные исследования МГД устойчивости и нелинейной эволюции плазмы в реальной геометрии тока-мака 78

§ 1. Исследование влияния тороидального вращения плазмы на

МГД устойчивость диссипативной плазмы.......... 78

п. 1. Расчет стационарного состояния вращающейся плазмы ............................. 78

п.2. Постановка задачи о линейной устойчивости плазмы с тороидальным вращением............ 81

п.З. Исследование влияния вращения плазмы на устойчивость тиринг-мод................... 83

§ 2. Моделирование эффектов внешнего винтового магнитного

поля на нелинейную МГД эволюцию вращающейся плазмы 90 § 3. Самосогласованное моделирование нелинейной МГД эволюции плазмы в разрядах с обращенным центральным магнитным широм..........................102

п.1. Анализ экспериментальных данных на токамаке

БШ-Б...........................104

п. 2. Линейная устойчивость исходного равновесия .... 115

п.З. Расчет нелинейной МГД устойчивости плазмы без

учета влияния источников...............119

п. 4. Нелинейная эволюция плазмы в присутствии источников. Моделирование МГД вспышек и процесса срыва в разрядах с обращенным центральным ши-ром. Сравнение результатов вычислений и экспериментальными данными на токамаке DIII-D .....121

п. 5. Выводы.......................... 124

§ 4. Численное моделирование нелинейной неустойчивости неоклассической тиринг-моды...................130

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

п.2. Самосогласованное моделирование нелинейной неустойчивости неоклассической тиринг-моды.....142

§ 5. Моделирование эффектов стабилизации неоклассических тиринг-мод с помощью радиально локализованного тока

ВЧ-увлечения...........................150

Заключение 169

Литература 170

Введение

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

Целью диссертации является разработка специфических математических моделей для описания нелинейной трехмерной МГД эволюции тороидальной плазмы; построение численных алгоритмов для решения нелинейных задач МГД устойчивости плазмы в реальной тороидальной геометрии; создание и тестирование нелинейного трехмерного численного кода NFTC (Nonlinear Full Toroidal Code); проведение численного моделирования ряда нелинейных эффектов, обнаруженных в экспериментах на современных установках токамак.

В диссертации для описания нелинейной эволюции плазмы в токамаке используется следующая МГД система уравнений, которая в безразмерной форме имеет вид: 8V

р^г = -р (V, V) V - VP + {(V х В) х В] + Re AV, (1)

\J Ь

f = Vx[V.xB]-Re-1? х [г, (V х В - jbs)], (2)

дР = V • (PV) - (Г- 1)PV .V + Vj. • (xjVXP) +

dt

d_ /р;

dt V/9r,

V|,P=(!f,v)p, v,p = vp-v„p, ! = | + (V,V)

+ V„ • (XfflV„P) + Qoh, (3)

= o, (4)

Первое уравнение описывает движение плазмы. Здесь через р, V, Р обозначаются соответственно плотность, скорость движения и газокинетическое давление плазмы. Уравнение (2) для магнитного поля В отражает закон Фарадея с учетом закона Ома. Третье уравнение описывает закон изменения давления сжимаемой плазмы, показатель адиабаты обозначается через Г. Последнее уравнение (4") описывает уравнение состояния плазмы в виде адиабатического закона. В системе уравнений (1)-(4) г) — функция профиля удельного сопротивления плазмы, ^ — плотность бутстрэп-тока, связанного с градиентом давления Р, фон — тепловой источник омического нагрева.

В уравнениях присутствуют следующие безразмерные параметры: число Рейнольдса Re = роУлан/щ = т„/та, магнитное число Рей-нольдса Rem = ¿¿оп^ан/щ = тц/тд, коэффициент продольной теплопроводности плазмы х\\ = {уао>н/х\\о)-1 — (тхц/та)-1 вдоль направления магнитного поля В, коэффициент поперечной теплопроводности х± = (^ао-н/х±о)~1 = гДе и х±о — характерные вели-

чины плотности, вязкости, продольной и поперечной теплопроводности

плазмы, /¿о = 47г • 10-7(kg • m/k2) — вакуумная магнитная проницаемость, 77о = 1/(^00")} а — электропроводность плазмы, а# — горизонтальный радиус тороидального плазменного шнура, va = А)/л/ЦоЩ — альфвеновская скорость, Bq — характерная величина тороидального магнитного поля в центре вакуумной камеры, тд = а^/уд — аль-фвеновское время, т„ = р(}а2н/щ, тд = Цоа^/щ — резистивное время, тх = а2н/х\\о, тх± = а2н/х±о- В задаче также определяются безразмерный коэффициент (3q = 2Р0/(-бо//хо), характеризующий отношение газокинетического давления Ро плазмы на магнитной оси системы к магнитному давлению, и число Маха М = Vq/Cs, Vo — характерная скорость вращения плазмы, Cs = уТРо/ро — скорость звука. Для современных экспериментов на установках токамак Re и Rem достигают значения 108, число Маха М ~ 0.09-0.6, коэффициент (3q достигает 12%, х± ~ 10-4-10"6,X|| - 102-103.

В диссертации предполагается рассматривать главным образом два аспекта вопросов, связанных с МГД моделированием плазмы в токама-ках. Во-первых, мы считаем плазму электропроводящей сжимаемой жидкостью, находящейся в сильно анизотропном магнитном поле. В приближении "токамака" продольная составляющая Вц магнитного поля вдоль большого обвода тора по величине обычно на порядок больше поперечной составляющей В±. Анизотропия магнитного поля и сжимаемость плазмы приводят к появлению нескольких характерных колебаний плазмы с существенно различными временами развития. Самым быстрым из этих колебаний является поперечная альфвеновская волна, связанная с продольным магнитным полем Рц и распространяющаяся со скоростью va-Такие колебания в токамаке устойчивы и нас не интересует. В диссертации исследуются неустойчивые длиноволновые и долгоживущие относительно альфвеновского времени та процессы. Во-вторых, в плазме существуют магнитные поверхности двух типов: рациональные, где силовые линии магнитного поля замыкаются сами на себя, совершая ровно п оборотов по большому обводу тора и га оборотов по малому обводу тора (т и п — натуральные числа) и иррациональные поверхности, где силовая линия не замыкаются на себя и плотно заполняет магнитную поверхность. Интересующие дас длиноволновые колебания развиваются в окрестности нескольких рациональных магнитных поверхностей с малыми числами тип. Эти поверхности разделяют плазму на несколько областей. Между указанными рациональными поверхностями плазму можно считать идеально проводящей и малыми параметрами Ре-1 и Ре"1 можно пренебречь. Здесь в основном соблюдается баланс между силой Ампера и

силой газокинетического давления плазмы. Каждая такая рациональная поверхность окружена тонким внутренним диссипативным слоем с шириной 6 ~ Ре"1/3. В этих внутренних слоях важны сопротивление, вязкость и динамика плазмы. Именно эти внутренние слои существенно определяют развитие неустойчивых диссипативных процессов с характерным временем тд в плазме. В диссертации в основном исследуется МГД эволюция плазмы, связанная с наличием внутренних диссипативных слоев.

Задача решается в так называемой квазитороидальной геометрии (р, где (р — тороидальный угол вокруг главной оси тора. Переменные (р, в) образуют полярную систему координат в поперечном сечении тора. Задача нелинейной МГД эволюции плазмы заключается в прослеживании временного развития совокупности трехмерных неустойчивых винтовых (д/д<р ф 0) возмущений плазмы относительно исходного аксиально-симметричного (д/дср = 0) равновесия {Уеч = 0, Веф Peq}, удовлетворяющего уравнениям:

УРеа = [Лед х ВеС1], V X Веч = Леф V • Вед = 0. (5)

Учитывая аксиальную симметрию решения система (5) сводится к скалярному уравнению Трэда-Шафранова [4], которое в диссертации численно решается методом обращения переменных [5], разработанным П.Н. Вабшцевичим, Л.М. Дегтяревым и А.П. Фаворским. Этот метод обеспечивает необходимую точность вычисления равновесных профилей для задачи МГД устойчивости плазмы.

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

Основная проблема численного решения связана с жесткостью нелинейной системы МГД уравнений. Характерные временные масштабы процессов, описываемых этими уравнениями, различаются на несколько порядков. Это продольные акустические волны, поперечные альфвенов-ские волны и быстрые магнйтозвуковые волны. Во многих задачах, связанных с моделированием МГД-процессов в высокотемпературной плазме, эта трудность обходится аналитическим подходом, предложенным Б.Б. Кадомцевым и О.П. Погуце [6]. При таком подходе выводится так называемая редуцированная система уравнений, в которой удалены члены, ответственные за быстрые магнйтозвуковые волны, связанные с сильным продольным полем. Редуцированная МГД модель, несмотря на су-

щественные упрощения, хорошо описывает устойчивость винтовых дис-сипативных мод и внешней винтовой моды со свободной поверхностью. Подход Кадомцева-Погуце позволяет найти решение для многих МГД проблем. В работах [7] [8] [9] [10] [11] [12] [13] были развиты линейные и нелинейные коды для редуцированной системы МГД уравнений. Были применены явная и полунеявные разностные схемы для численного решения системы уравнений. Особенно интересно отметить предложенную в работе [13] идею исследования нелинейной эволюции плазмы методом проектирования на неустойчивое многообразие, в которой приближенное решение строится в виде функциональных рядов по амплитудам неустойчивых в линейном приближении гармоник. Недостаток редуцированной модели состоит в том, что данная модель работает только в цилиндрическом приближении, не учитывает сжимаемости среды и давления плазмы. В современных токамаках с большим давлением, малым аспектным отношением и некруглым поперечным сечением, все колебания сцеплены друг с другом, их трудно разделить; кроме того, за счет большой тороидальности токамака МГД моды сильно связываются и основными могут оказаться сразу несколько мод. Поэтому необходимо рассмотреть полную, нередуцированную МГД модель, в которой с учетом сжимаемости и давления присутствуют одновременно перечисленные выше типы колебаний. В условиях современных экспериментов время тд развития ре-зистивной неустойчивости на несколько порядков больше, чем характерное время та быстрых альфвеновских колебаний. Явные схемы связаны с сильным ограничением на временной шаг по условию Куранта и ориентируются на, устойчивые быстрые альфвеновские колебания, не представляющие интереса. Эти схемы не позволяют исследовать нелинейные диссипативные МГД процессы плазмы даже на самой современной машине. Явная схема была использована в первом трехмерном нелинейном коде ТОР, созданном в работе [14][15]. В коде использована схема Лакса-Вендроффа класса предиктор-корректор на эйлеровой сети. С помощью кода ТОР авторы исследовали взаимодействие различных мод в полной МГД модели. Жесткое ограничение на шаг по времени практически не позволяет исследовать нелинейные резистивные процессы, связанные с высокой проводимостью плазмы. Для того, чтобы преодолеть ограничение на шаг по времени при численном решении нелинейной жесткой МГД системы, необходимо построить такую разностную схему, которая исключает быстрые альфвеновские волны. Во второй половине 80-х годов был разработан полунеявный метод [16], который был применен в нелинейном МГД коде СТБ [17] [18], созданном в Техасском университе-

те (США) и в коде XTOR [19] [20], разработанный во Франции. Основная идея полунеявного метода для МГД задачи заключается в следующем. Сначала построим такой оператор Lo, который с одной стороны имеет сравнительно простую структуру, с другой стороны, существенно отражает быстрые альфвеновские колебания в МГД системе. Оказывается, в качестве такого оператора нужно выбрать LqV = [V х V х [V х Со] х С0], где вектор Со должен быть направлен параллельно основному магнитному полю Во- Далее, исходная система (1)-(3) численно интегрируется по двушаговой схеме предиктор-корректор, при этом на этапе корректор из левой и правой частей разностного уравнения для скорости V выч�