В последнее время изучение процессов нелинейной эволюции плазмы в токамаке в магнитогидродинамическом (МГД) приближении с помощью численного моделирования становится все более актуальным [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 — характерная величина тороидального магнитного поля в центре вакуумной камеры, тд = а^/уд — аль-фвеновское время, т&bdquo- = р (}а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 вычитается член LY = t2LqV, причем член с левой стороны аппроксимируется неявно, что обеспечивает абсолютную устойчивость разностной схемы в целом по отношению к быстрому альфвеновскому движению, и тем самым существенно поднимает временной шаг численного интегрирования. По существу полунеявный метод эквивалентен добавлению к исходному уравнению движения (1) регуляризирующего оператора эллиптического типа L = г2Lo-' р (1 — T2Lo) d*V/dt = F, где F — правая часть уравнения (1). Когда временной шаг т стремится к нулю, мы переходим к исходному уравнению. В полунеявном методе относительно простая структура оператор-регуляризатора LoV с одной стороны не сильно усложняет процедуру решения полученных разностных уравнений, и это является преимуществом данного метода, с другой стороны, из-за анизотропии магнитного поля В нельзя точно вычесть члены, связанные с быстрыми альфвеновскими колебаниями. Более того, в нелинейной МГД задаче с изменением магнитного поля вычитание становится все более не точным. На развитой нелинейной стадии, когда нелинейные члены типа (V-V)V или (V х В) х В становятся большими, полунеявный метод уже не обеспечивает безусловную устойчивость схемы. Другая полунеявная схема для системы МГД уравнений впервые применялась в нелинейном коде FAR [21][23][24][22], разработанном в ОкРидже (США). Этот код решает полную МГД задачу в тороидальной геометрии с некруглым по. перечным сечением плазменного шнура. Линейная часть МГД уравнений в коде FAR интегрируется неявным образом. Абсолютная устойчивость неявной схемы в линейных операторах позволяет коду выбрать любой временной шаг т в линейных вычислениях и тем самым выбрать его в соответствии с интересующим нас движением. При разумном выборе г данная схема, как и в методе обратной итерации, позволяет очень быстро найти собственные значения, причем не только максимальное собственное значение. Нелинейные операторы в коде FAR учитываются явным образом, это связано с тем, что в нелинейном расчете неявная схема приводит к необходимости обращения больших матриц со сложной структурой на каждом шаге по времени. Кроме того, на каждом временном шаге для решения нелинейных уравнений приходится применять итерационные методы, например, метод Ньютона. Явное интегрирование нелинейных операторов в коде FAR сильно ограничивает шаг по времени. Типичные расчеты кода FAR на машине CRAY-II для линейной задачи занимает десятки минут, а для нелинейной задачи — десятки часов.
В системе МГД уравнений существуют квадратичные члены типа (В, V) B. Наличие такой квадратичной нелинейности часто обусловливает переход энергии Фурье-гармоник решения от длинноволновых гармоник в коротковолновые в численных расчетах. В результате нелинейность способствует накоплению большого коротковолнового шума, который может полностью разрушить решение. Для фильтрации таких коротковолновых гармоник в нелинейной задаче необходимо использовать процедуру сглаживания. Такая процедура была применена в работе Шнака [25] [26], где после каждой временной итерации сеточные функции явным образом модифицируются путем добавления разностного аналога оператора диффузии по радиальному направлению. Другая процедура сглаживания была предложена в нелинейном МГД коде FAR. В этом коде при аппроксимации уравнения 1-ого порядка относительно дифференцирования переменной функции, а по радиальной переменной р интегрируется МГД уравнение относительно, а + Daвместо вместо а, где Da в нелинейных вычислениях служит параметрам сглаживания. В нелинейном коде XTOR был-использован другой подход сглаживания. В этом коде к уравнению движения для скорости плазмы V был добавлен эллиптический оператор сД V, который доминирует квадратичные члены от магнитных возмущений. Константа с составляет 10 ~ 10″ 4. Этот диффузионный член добавлен неявно.
С точки зрения корректности более адекватной, чем нелинейная идеальная модель, является постановка резистивной нелинейной МГД задачи. При численном решении такой задачи принципиальным является четкое разделение диссипатывных эффектов за счет резистивного члена в уравнениях, за счет процедуры сглаживания и за счет дискретности разностной схемы. Очень важно построить такую разностную схему, позволяющую правильно описывать физическую неустойчивость самой системы МГД уравнений. Для этого нужно добиваться, чтобы физический диссипативный эффект всегда превалировал над численным. В то же время численная диссипация должна быть достаточно большой, чтобы ликвидировать коротковолновой шум, возникший из-за квадратичной нелинейности. В связи с этим следует напомнить работу С. И. Мухина.
27], в которой автор предложил класс разностных схем с искусственной дисперсией для одномерных уравнений магнитной гидродинамики в ла-гранжевых координатах. Существенным моментом в этих схемах является учет баланса между численной дисперсией и диффузией. Применение таких схем в трехмерной задаче со сложной геометрией оказалось не так просто, хотя попытка рассмотрения данного метода в диссертации была бы целесообразной.
Нелинейная многомерная МГД задача характерна наличием множества локальных минимумов энергии системы около исходного аксиально-симметричного равновесного состояния. Каждый такой минимум соответствует локальному винтовому равновесию плазмы. Цель нелинейного моделирования МГД эволюции плазмы по сути дела состоит именно в отыскании ближайшего локального винтового равновесия около заданного аксиально-симметричного начального состояния. Но часто бывает и другая ситуация, когда начальное состояние уже находится на точке локального минимума. Другими словами, начальное состояние системы линейно устойчиво, при этом существует эффект жесткого возбуждения системы, т. е. система за счет конечного возмущения выходит из одного локального минимума в другой. Такой пороговый эффект наблюдается в экспериментах при наличии внешних резонансных винтовых магнитных полей, при дестабилизации тиринг-моды из-за неоклассических эффектов. В этих случаях выбор начальных условий существенно влияет на дальнейшую эволюцию системы. Поэтому очень важно выбрать адекватные начальные условия. В случае, когда отсутствует явный критерий выбора таких условий, приходится найти адекватные условия на основе серийных расчетов с различной начальной амплитудой.
Следующая нелинейная проблема связана с возможностью перезамыкания магнитных поверхностей вследствие конечной проводимости плазмы, что не позволяет при вычислениях использовать лагранжевы координаты, часто применяемые для решения уравнений газовой динамики.
28], из-за постоянного изменения топологии магнитных поверхностей.
В процессе нелинейного развития МГД неустойчивости, в результате, перезамыкания магнитных силовых линий в плазме возникают магнитные острова. Эти магнитные острова образуются в результате разрушения рациональных поверхностей и имеют различные винтовые структуры. При своем неустойчивом развитии они могут перекрываться друг с другом, в результате чего происходит процесс стохастизации магнитных силовых линий в области перекрытия [29]. Этот процесс стохастизации сильно изменяет процессы теплопереноса в плазме на нелинейной стадии развития. Поэтому в нелинейном расчете необходимо учитывать изменение коэффициентов переноса в процессе эволюции плазмы.
Многомерность задачи представляет существенную трудность в численном решении МГД задачи. Равновесная конфигурация плазмы имеет аксиально-симметричную двумерную структуру. За счет МГД активности в плазме возбуждается трехмерная квазивинтовая структура, близкая к исходному равновесному состоянию в силу малости амплитуды винтовой составляющей. В реальной геометрии современных токамаков винтовые гармоники различной структуры тесно связаны между собой. В развитой стадии нелинейной эволюции в результате взаимодействия этих гармоник с различной винтовой1 структурой в плазме нарушается винтовая симметрия, поэтому в такой стадии учет трехмерности неустойчивости плазмы является существенным. Характер поведения системы в трехмерном пространстве может быть совсем иным по сравнению с двумерным случаем. Имеется очень мало теоретических результатов и оценок относительно нелинейной трехмерной МГД эволюции, что создает главную трудность при тестировании численных расчетов. С точки зрения вычислительного аспекта трехмерная задача требует на порядки больших машинных ресурсов. Полная система МГД устойчивости сжимаемой плазмы в трехмерной тороидальной геометрии состоит из 7 уравнений. Все эти уравнения связаны друг с другом. Каждая физическая функция из этих уравнений зависит от трех пространственных и одной временной переменных. Поэтому здесь возникает необходимо, сть экономить машинную память. Один из способов решения такой проблемы — использование так называемого полу спектрального метода для численного решения МГД уравнений в тороидальной системе координат. Полунеявный метод используется в современных МГД кодах CTD, FAR, XTOR и MH3D [30][31][32] в Принстоне (США), а также в коде NFTC, разработанном в диссертации. В данном методе искомые функции разлагаются в ряды Фурье по пространственным угловым переменным тора в и <�р: и (р, е,.
)],.
V т где через U обозначена совокупность функций V, В и Р. По радиальному направлению р пишется конечно-разностная схема. Использование Фурье представления решения по угловым переменным тора создает сложные свертки в нелинейных вычислениях. Например, для вычисления произведения типа w (p, в, ф) = и (р, в, ф) • г>(р, 0, (р) через свертки относительно полоидальных гармоник требует 0(М2) операций при фиксированном рис/?, где М — число гармоник по в. В кодах Шна-ка и СТБ предлагают так называемый псевдоспектральный алгоритм для вычисления значения т (р, 0,<�р). Суть этого алгоритма заключается в следующем. Сначала получают значения функций и и V в точках (рг-, в], фи), используя дискретное Фурье-представление, затем вычисляют гу (рг, 9 $, (рк) = и (рг, 9<�рк) 9^ (рк) как простое перемножение двух чисел, после чего вычисляют Фурье-гармоники гитп (р?) методом обратного преобразования Фурье. На первом и последнем шагах применяется метод быстрого дискретного преобразования Фурье (БДПФ), который требует 0(Мк^2М) операций умножения в отличие от 0(М2) для метода свертки. Недостаток псевдоспектрального алгоритма состоит в генерации коротковолнового шума из-за квадратичной нелинейности. Для устранения этого шума приходится брать массив длины ~ |М для представления Фурье-гармоник, что требует большей машинной памяти. В наиболее крупном проекте ШМЫОБ [33], разрабатываемом в настоящее время в США, используется метод конечных элементов для аппроксимации функций на поперечном сечении тора с тщательно разработанной сеткой блочного вида. Этот код должен решать одно-жидкостные и дву-жидкостные МГД уравнения в реальной тороидальной геометрии тока-мака.
Равновесная плазменная конфигурация, получаемая в современных разрядах, имеет поперечное сечение с большой вытянутостью, треуголь-ностью и сложной формой магнитных поверхностей. Форма поперечного сечения плазменного шнура в экспериментах часто оптимизирована так, чтобы все крупномасштабные идеальные МГД неустойчивости были стабилизированы в плазме. Оставшиеся резистивные неустойчивости имеют очень малый инкремент 7 ^ Ю-41, что создает определенную трудность при численном вычислении 7. Кроме того, форма магнитных поверхностей в реальной равновесной конфигурации создает сильное сцепление гармоник. Для адекватного описания такой формы требуется минимум 40 полоидальных гармоник в Фурье представлении равновесного решения и соответствующих метрических элементов. Для сравнения в торе с круглым поперечным сечением плазмы 3 гармоники может быть, достаточно для сравнительно точного представления метрических элементов.
В диссертации система МГД уравнений, описывающая процесс нелинейной эволюции диссипативной плазмы, содержит несколько малых параметров, в том числе коэффициенты Де-1, Ле~} и х±-, связанные с вязкостью, сопротивлением и поперечной теплопроводностью плазмы. Эти малые параметры играют существенную роль в диссипативном слое около рациональных магнитных поверхностей. В отсутствии сопротивления и вязкости, в «идеальном» случае уравнения имеют сингулярные решения. Наличие сопротивления и вязкости ликвидирует сингулярность решений идеальных МГД уравнений. Рассмотрим влияние малого параметра Ре" 1. В современных разрядах основные МГД эффекты связаны с неустойчивостью резистивных мод, зависящих от сопротивления плазмы. Моделирование МГД эволюции высокотемпературной плазмы с учетом конечной проводимости представляет собой задачу с малым параметром .Re" 1 ~ Ю-8 при старшей производной по пространству АВ. В плазме токамака на рациональной (резонансной) магнитной поверхности конвективный член (Beq, V) V обращается в нуль. На равновесной резонансной поверхности у правой части уравнения магнитного поля главный член, связанный с оператором (Beq, V), исчезает, остается только резистивный член с малым параметром Ре" 1. Поэтому на резонансной поверхности даже очень малое сопротивление важно. При численном исследовании плохая аппроксимация конвективного члена (Beq, V) V может создать большую численную диффузию, которая полностью подавляет резистивный эффект на резонансной поверхности плазмы. Необходимо иметь высокий порядок аппроксимации конвективного члена, входящего в уравнение для магнитного поля. В коде FAR эта проблема решается путем использования так называемой натуральной системы координат [34], которая связана с магнитными силовыми линиями равновесного поля Beq. При Фурье разложении решения по угловым переменным в такой системе координат, конвективный оператор дифференцирования вдоль поля (Beq, V) становится алгебраическим и ошибки в описании резистивных диссипативных слоев при численном моделировании наименьшие. Однако, натуральная система координат в общем случае не ортогональна. Появляются сложные метрические коэффициенты, что осложняет дальнейшие вычисления. Недостаток натуральной системы координат также состоит в том, что в нелинейных расчетах естественная на начальной стадии система координат затем становится неестественной.
При разработке численной схемы для резистивной плазмы возникает еще одна проблема. Дело в том, что с ростом проводимости ширина рези-стивного слоя быстро убывает, что создает принципиальные трудности при создании дискретной математической модели. В нелинейном расчете ширина резистивиого слоя все время меняется. Для обеспечения адекватного моделирования поведения системы внутри резистивного слоя, нужно зарезервировать достаточное количество точек разностной сети внутри данного слоя в процессе вычислений. Другим выходом из данной ситуации является создание адаптирующейся сетки. Для полной, многомерной МГД задачи очень трудно найти хороший алгоритм адаптации сетки. В работах [35] [36] был развит метод подгоночных коэффициентов. Данный метод хорош тем, что он равномерно сходится относительно малого параметра Ae" 1 при шаге сетки h, стремящемся к нулю, и позволяет исследовать диссипативную неустойчивость плазмы для значения Rem до 1010. Обобщение данного метода на случай нелинейной многомерной МГД системы сталкивается с определенными трудностями, связанными с необходимостью нахождения двумерной асимптотики решения около диссипативного слоя.
Наличие диссипативных членов в МГД уравнениях усложняет задачу нахождения собственных значений линеаризованной системы. Эта задача имеет важное значение для дальнейшего моделирования нелинейной эволюции плазмы. Линейная резистивная МГД задача является несамосопряженной, поэтому традиционные методы, основанные на вариационном принципе и используемые в спектральных кодах PEST [37] и ERATO [38], не работают. В работе [39] решается непосредственно спектральная задача для линеаризованных уравнений, содержащих искомую частоту сложным образом, поэтому для ее нахождения используется метод стрельбы. Решение уравнений проводится специальным вариантом метода исключения Гаусса, который учитывает специфику и разряженность матрицы. Практические затруднения при реализации данного метода возникают при увеличении числа Remr когда уменьшается расстояние между точками дискретного спектра (так как при Rem —" сю спектр становится непрерывным).
В работе японских исследователей [34] аналогично [37] непосредственно решается исходная временная задача методом сквозного счета по пространству на каждом шаге по времени. Метод не учитывает наличия внутренних диссипативных слоев, поэтому для достижения достаточной точности требуется большое количество точек в пространственной сети. Однако, стандартность матричной задачи позволяет проводить вычисления в любой геометрии и общая производительность повышается по срав-. нению с [37] [40], ориентированными на специфику задачи в той или иной геометрии.
В линейном коде NOVA [41] предлагается подход не вариационного типа с использованием метода конечных элементов (метода Галеркина) для исследования линейной идеальной МГД устойчивости. Код NOVA позволяет точнее найти собственные значения с меньшими машинными затратами по сравнению с вариационными кодами PEST и ERATO. Ре-зистивный линейный код NOVA-R [42] использует подход, подобный подходу в коде NOVA. Разница заключается в том, что вместо одного окончательного обыкновенного дифференциального уравнения П-ого порядка получается три таких уравнения. Код MOVA-R может вычислить инкремент резистивной моды при реальной проводимости плазмы в экспериментах. Недостатками серии кодов NOVA являются: во-первых, сложная аналитическая процедура сведения исходной задачи МГД устойчивости к алгебраической задаче нахождения собственных значений матрицы, при малом изменении постановки исходной задачи приходится заново проделать всю аналитическую процедуру, что ограничивает универсальность кода NOVAво-вторых, сложность модификации для исследования нелинейной МГД эволюции.
В современных разрядах на установках токамак многие МГД эффекты тесно связаны с наличием давления и сжимаемости плазмы. В разрядах с обращенным центральным магнитным широм изменение давления в центральной зоне плазмы играет принципиальную роль в эффектах типа пилообразных колебаний температуры плазмы. В современных экспериментах наблюдается неустойчивость так называемой неоклассической тиринг-моды, связанной с наличием бутстрэп-тока (неоклассического МГД тока), который создается в плазме за счет градиента профиля давления. Учет сжимаемости плазмы также создает определенную трудность для дискретного моделирования МГД эволюции. Дело в том, что в сжимаемой плазме с высоким давлением возникает новые типы колебаний. Для численного расчета таких колебаний в явной схеме требуется малый временной шаг, что часто делает нелинейные вычисления невозможными из-за слишком больших затрат машинного времени. Поэтому в диссертации используется полностью неявная схема в нелинейных вычислениях.
Наконец, наличие аксиально-симметричного вращения плазмы также порождает определенную проблему в численном моделировании. Тороидальное вращение плазмы возникает за счет инжекции пучков нейтралов в плазму. В экспериментах на установке DIII-D было замечено, что тороидальное вращение экранирует плазму от паразитного резонансного внешнего магнитного поля и стабилизирует МГД неустойчивость. Торможение же вращения приводит к вспышке неустойчивых МГД мод. Аксиально-симметричное вращение меняет исходное равновесие плазменной конфигурации. Плазма переходит в некоторое стационарное состояние, у которого может радикально изменяться характер устойчивости.
Разработанный в настоящей работе код NFTC является нелинейным развитием линейного полуспектрального кода FTC [43]. Код NFTC решает полные сжимаемые МГД уравнения в общей тороидальной геометрии с использованием так называемой потоковой системы координат. Учитываются разные диссипативные процессы (сопротивление и вязкость). По идеологии разработки код NFTC близок к нелинейному коду FAR. Но в отличие от последнего в коде NFTC используется полностью неявная безусловно устойчивая схема. Хотя в реальном нелинейном расчете шаг по времени конечно всегда должен соответствовать характерному масштабу развития исследуемого процесса. Другое различие кода NFTC от FAR состоит в том, что код FAR решает МГД задачу в специальных переменных-потенциалах. Код NFTC интегрирует МГД уравнения непосредственно относительно (V, В, Р).
Отличительными чертами кода NFTC, развитого в настоящей диссертации, по сравнению с имеющимися кодами являются: 1) применение полностью неявной схемы для решения нелинейной трехмерной МГД задачи- 2) включение произвольного ширового вращения плазмы- 3) наличие внешнего винтового магнитного поля- 4) включение неоклассических эффектов- 5) возможность исследования несимметричной вдоль главной оси тора равновесной диверторной конфигурации- 6) гибкость кода: в коде линейные, квазилинейные и нелинейные операторы реализованы в отдельных процедурах. Задавая соответствующие параметры, можно легко исключать любые из этих операторов, что позволяет исследовать линейную, квазилинейную, нелинейную задачи отдельно и вместе. Кроме того, в коде тороидальные числа п можно задавать любые и в любом количестве. Полоидальные числа т для фиксированного п задаются на диапазоне [Мт[п (п), Мтах (п)]} что позволяет в ряде случаев специально исследовать МГД устойчивость интересующих нас полоидальных гармоник- 7) вычисление комплексных собственных значений [44]. Наличие вращения плазмы требует вычислять комплексные собственные значения, для чего были использованы метод обратной итерации для комплексных собственных значений и метод, использующий значения решений в нескольких шагах по времени- 8) неравномерная сетка по радиальному направлению, которая дает возможность очень точно описать диссипативный слой методом сгущения сетки.
Диссертация состоит из взедения, трех глав, заключения и литературы.
В первой главе разработаны математические модели для описания нелинейной МГД устойчивости плазмы токамака в тороидальной геометрии. Поставлена математическая задача нахождения нелинейной трехмерной МГД эволюции тороидальной плазмы. Первый параграф посвящен разработке МГД моделей нелинейной неустойчивости сжимаемой одно-жидкостной плазмы. В токамаке равновесное состояние обладает аксиальной симметрией. Задача нелинейной МГД эволюции плазмы заключается в прослеживании временного развития совокупности неустойчивых винтовых возмущений плазмы относительно заданного равновесия. В данном параграфе сформулирована исходная система МГД уравнений для винтовых возмущений и соответствующей квазилинейной аксиально-симметричной поправки к равновесию. В полученной системе уравнений учитываются разные физические эффекты, для их численного моделирования разработаны соответствующие нелинейные МГД модели. В п. 1. рассмотрена модель МГД устойчивости стационарного состояния плазмы с шировым тороидальным вращением. В п. 2. построена модель нелинейной МГД устойчивости плазмы при наличии внешнего резонансного винтового магнитного поля. В п.З. разработана модель порогового возбуждения неустойчивости тиринг-моды с включением неоклассических эффектов, в том числе бут^трэп-тока и конечной продольной теплопроводности плазмы вдоль магнитных силовых линий. Записан оператор продольной теплопроводности через возмущенное магнитное поле. В п. 4. рассматривается самосогласованная модель для описания взаимодействия транспортных процессов и диссипативной МГД неустойчивости плазмы. В данном пункте построены модели разного рода источников: источника тока, возбуждаемого внешним источником * ВЧ-поля, источник тепла за счет омического нагрева фон и дополнительного нагрева ЯтхВ системе МГД уравнений выделяются следующие безразмерные параметры: коэффициент /?0, число Рейнольдса Ре, магнитное число Рей-нольдса Рет, число Маха М, коэффициенты поперечной и продольной хц теплопроводности. Во втором параграфе проведена оценка значений этих коэффициентов для установки токамак ОIII-Б (Дженерал Ато-микс, США). Также вычислены безразмерные значения альфвеновской скорости, скорости звука и скорости тороидального вращения плазмы в экспериментальных условиях. В третьем параграфе первой главы описа-, на потоковая система координат, в которой силовые линии равновесного магнитного поля являются прямыми. Для данной системы координат вычислены соответствующие метрические элементы. В данном параграфе также показано, как записать дифференциальные операторы МГД уравнений в потоковой системе координат для контравариантных компонент векторов решения. В последнем параграфе главы 1 определена структура численного исследования МГД устойчивости и нелинейной эволюции тороидальной плазмы. Процесс численного моделирования МГД устойчивости делится на три этапа: первый этап посвящен численному нахождению осесимметричного стационарного состояния плазмы с шировым тороидальным вращениемвторой этап связан с предварительным исследованием линейной МГД устойчивости плазмы относительно выделенного стационарного состоянияна третьем этапе непосредственно проводится численное моделирование нелинейной МГД эволюции плазмы.
Вторая глава диссертации посвящена построению численных алгоритмов для решения задачи нелинейной МГД устойчивости плазмы в трехмерном пространстве, а также созданию трехмерного, нелинейного численного кода МГТС. В первом параграфе описано полуспектральное представление решения исходной МГД задачи. Искомые функции разлагаются в ряд Фурье по полоидальному и тороидальному углам. Задача МГД устойчивости плазмы решается относительно Фурье гармоник, зависящих от радиальной переменной и времени. Во втором параграфе проведена запись и формирование сверток Фурье гармоник искомых функций и метрических элементов. Классифицировано три типа сверток. Проведена оценка количества арифметических операций для сверток разной размерности. В параграфе § 3 построена потоковая разностная схема по радиальному направлению р. В данной схеме радиальные компоненты векторов решения определяются в точках, сдвинутых на полшага относительно тангенциальных компонент, что обеспечивает центральную разностную аппроксимацию, по радиальной переменной р с точностью до второго порядка по радиусу без расширения шаблона аппроксимации. Центральным местом в численном решении нелинейных МГД уравнений занимается проблема-интегрирования—по времени. Построение разностной схемы по времени является содержанием четвертого параграфа. Здесь разработана полностью неявная схема для системы нелинейных МГД уравнений относительно аксиально-симметричных и винтовых возмущений с различными тороидальными числами. Неявная аппроксимация по времени порождает систему нелинейных уравнений относительно искомых функций на верхнем временном слое. Для решения нелинейных уравнений разработан блочный метод Гаусса-Зейделя-Ньютона с учетом специфики системы МГД уравнений. Полученная после этого линейная система содержит порядка 104 уравнений. Соответствующая матрица имеет блочно-трехдйагональную структуру и для ее обращения применяется метод матричной прогонки. В четвертом параграфе разработан алгоритм сглаживания решения по радиальной переменной р на каждом временном слое для устранения эффекта филаментации профилей решения при численном интегрировании нелинейных уравнений. Пятый параграф посвящен описанию метода обратной итерации нахождения собственных значений краевой задачи для полной системы МГД уравнений. Модификация данного метода для поиска комплексных собственных значений в задачах устойчивости вращающейся плазмы предложена в шестом параграфе. Здесь также описан метод Уилкинсона с использованием значений решения в трех последовательных временных шагах для вычисления комплексных собственных значений в задачах МГД устойчивости. В § 7 перечислены проблемы и их решения при реализации численного трехмерного кода NFTC. Полученный код представляет собой сложную систему с взаимодействием разных программных компонент. В данном параграфе также построена схема вычислений с помощью кода NFTC. Параграфы § 8 — § 10 посвящены тестированию созданного кода. В восьмом параграфе подробно исследована сходимость вычислений собственных значений, А в зависимости от численных параметров разностной схемы в задачах линейной МГД устойчивости, исследована точность вычисления, А при различных значениях физических параметров задачи, исследована устойчивость и сходимость нелинейных расчетов в зависимости от численных параметров. Выяснен профиль затрат машинных ресурсов кодом NFTC на современных вычислительных машинах. В § 9 проведено сравнение дисперсионных кривых, вычисленных кодом NFTC, с результатами спектрального кода ERATO. Рассмотрен случай идеальной МГД устойчивости плазмы для известного аналитического найденного равновесия. В последнем параграфе главы 2 показано сравнение численных расчетов с теоретическими решениями уравнения Машке-Перрина для стационарных состояний плазмы с постоянным радиальным профилем тороидального вращения.
В третьей главе диссертации с помощью кода NFTC проведено численное моделирование ряда нелинейных физических эффектов, обнаруженных в современных экспериментах на установках токамак. Первый параграф посвящен исследованию влияния тороидального вращения плазмы на МГД устойчивость диссипативной плазмы. В п. 1. данного параграфа найдено стационарное состояние вращающейся плазмы при разных значениях числа Маха. Обнаружена немонотонная зависимость величины стационарной коррекции магнитного поля от коэффициентаRe" 1. Исследовано влияние профиля тороидальной скорости вращения плазмы на характер полученного стационарного состояния. В п. 2. дана постановка задачи о линейной МГД устойчивости диссипативной плазмы с произвольным тороидальным вращением. В следующем пункте проведено подробное исследование влияния ширового тороидального вращения плазмы на МГД устойчивость тиринг-моды в токамаке. Вычислены кривые зависимости инкремента тиринг-моды от частоты вращения плазмы при различных значениях коэффициента вязкости. Для разного значения тороидального вращения также исследована зависимость инкремента тиринг-моды от сопротивления плазмы, коэффициента Д) и коэффициента запаса q = pB! p/(RBe) (где р, R — малый и большой радиусы тора, B^^Bq — тороидальная и полоидальная компоненты равновесного магнитного поля Beq). Во втором параграфе изучена проблема нелинейной МГД неустойчивости вращающейся плазмы при наличии внешнего резонансного винтового магнитного поля. Р’асчеты проводились для реального равновесия, полученного в разряде токамака DIII-D. Обнаружено, что без тороидального вращения плазма усиливает резонансное винтовое внешнее магнитное поле примерно в 6 раз. Вращение экранирует плазму от внешнего поля и тем самым, стабилизирует МГД устойчивость плазмы. На основе численных расчетов найдена критическая кривая зависимости амплитуды внешнего поля от скорости тороидального вращения, выше которой проявляется неустойчивость. Проведено сравнение полученной кривой с упрощенной теоретической оценкой и экспериментальными данными. Выяснена роль диссипативных эффектов в нелинейном процессе жесткого возбуждения МГД неустойчивости плазмы. Третий параграф посвящен самосогласованному моделированию нелинейной МГД эволюции плазмы в разрядах с обращенным центральным магнитным широм (Negative Central Shear, NCS). В первом пункте данного параграфа проведен анализ экспериментально наблюдаемых эффектов МГД вспышек на токамаке DIII-D в NCS разрядах. Предварительный анализ позволяет определить характеристики временного поведения амплитуды и радиального профиля винтового возмущения и вызванного им квазилинейного возмущения внутреннего полоидального магнитного поля. Эти величины можно сравнивать с экспериментальными измерениями. Исследование линейной МГД устойчивости для равновесия с обращенным центральным широм проведено в п. 2. Обнаружено две основные неустойчивые моды в таких разрядах: резистивная перестановочная мода, локализованная около внутренней рациональной поверхности q = 3, и двойная тиринг-мода, локализованная около внутренней и внешней поверхностей q — 3. Исследована количественная зависимость инкремента неустойчивости от сопротивления и тороидального вращения плазмы. В третьем пункте приводятся результаты моделирования нелинейной МГД эволюции плазмы в NCS разрядах без учета влияния источников. В этом случае численные расчеты показывают наличие двух стадий в нелинейной эволюции, качественно соответствующих экспериментальным наблюдениям: первых МГД вспышек и окончательной неустойчивости срыва в конце нелинейного развития. Для количественного сопоставления результатов численных расчетов с экспериментальными данными потребовалось включение источников токов и тепла. Исследованию нелинейной самосогласованной эволюции плазмы с учетом действия источников посвящен п. 4. данного параграфа. Основные выводы исследования МГД неустой-чивостей, наблюдаемых в разрядах с обращенным центральным широм сделаны в последнем пункте. Четвертый параграф посвящен численному моделированию нелинейной неустойчивости неоклассической тиринг-моды, наблюдаемой в современных экспериментах при наличии большой доли бутстрэп-тока в плазме. В п. 1. проведены численные исследования эффекта наличия большой, но конечной продольной теплопроводности с коэффициентом Х\ вдоль магнитных силовых линий. Расчеты показали, что конечная величина х\ существенно влияет на радиальный профиль возмущенного давления внутри магнитных островов, тем самым влияет на величину бутстрэп-тока. В первом пункте подробно сопоставлены результаты численных расчетов с теоретическими оценками. В следующем пункте проведено самосогласованное моделирование нелинейной неустойчивости неоклассической тиринг-моды при наличии бутстрэп-тока и конечной продольной теплопроводности. Обнаружен пороговый эффект нелинейного возбуждения неоклассической тиринг неустойчивости, связанный с конечной величиной начального возмущения. Дана количественная оценка критической ширины магнитных островов, больше которой возбуждается неустойчивость. В последнем параграфе третьей главы исследованы эффекты стабилизации неоклассической тиринг-моды с помощью радиального локализованного тока jccj, созданного внешним источником ВЧ-поля.
В заключении приведены основные результаты диссертации.
Основные результаты опубликованы в работах [97] [98] [99] [100] [101] [102] [103] [104] [105] [106] [107] [108]. Автор выражает глубочайшую благодарность своему научному руководителю, профессору A.M. Попову за постановку задач, руководство и большую помощь в работе. Благодарю также доцента A.B. Педоренко за предоставление программы расчета метрических элементов и многочисленные полезные обсуждения проблем МГД устойчивости. Считаю своим приятным долгом выразить признательность доценту H.H. Поповой за бесценные советы в программной реализации кода NFTC и предоставление программ для визуализации результатов численных расчетов. Хочу также поблагодарить профессора Д. П. Костомарова и профессора Ю. Н. Днестровского, регулярно принимавших участие в обсуждении полученных результатов и сделавших ряд ценных замечаний.
Заключение
.
В диссертации разработаны и исследованы МГД модели нелинейной эволюции неустойчивой плазмы в реальной тороидальной геометрии установок токамак.
1) Разработаны математические модели для описания сжимаемой одно-жидкостной магнитной гидродинамики тороидальной плазмы.
2) Разработаны полуспектральные численные алгоритмы для решения нелинейной трехмерной задачи МГД устойчивости плазмы в реальной тороидальной геометрии. Исследованы свойства блочного метода Гаусса-Зейделя-Ньютона при решении нелинейных многомерных МГД уравнений.
3) Реализован и тестирован трехмерный численный код №ТС для моделирования нелинейной МГД эволюции тороидальной плазмы на установках токамак.
4) Проведено трехмерное численное моделирование процессов в плазме токамака. Впервые в рамках трехмерной модели исследованы стабилизирующее влияние тороидального ширового вращения плазмы на МГД устойчивость тиринг-мод. Изучен эффект жесткого возбуждения МГД неустойчивости плазмы при наличии внешнего резонансного винтового магнитного поля. Проведено моделирование МГД неустойчивостей в перспективных разрядах с обращенным центральным магнитным широм. Исследованы процесс пороговой нелинейной дестабилизации неоклассической тиринг-моды, приводящей к ограничению максимального давления плазмы.
