CC BY
14
3. D'yakonov E.G. Operator problems in strengthened Sobolev spaces and numerical methods for them // Lect. Notes Сотр. Science. 1997. 1196. P. 161-169.
4. Дьяконов Е.Г. Оценки iV-поперечников в смысле Колмогорова для некоторых компактов в усиленных пространствах Соболева // Изв. вузов. Математика. 1997. № 4 (119). С. 32-50.
5. Дьяконов Е.Г. Энергетические пространства и их применения. М.: Издательский отдел ф-та ВМиК МГУ, 2001.
6. Дьяконов Е.Г. Некоторые типы усиленных пространств Соболева в случае многомерных областей с нерегулярной границей // Докл. РАН. 2003. 390. № 1. С. 19-23.
7. D'yakonov E.G. Asymptotically optimal algorithms for stationary problems in energy spaces // HERMIS-дтг. 2002. 3. P. 25-50.
8. Maz'ya V.G., Poborchii S.V. Differentiable functions on bad domains. Singapore: World Scientific Publishers, 1997.
9. Дьяконов Е.Г. Усиленные пространства Соболева для областей с нерегулярной границей // Тр. МИАН. 2003. 243. С. 213-219.
Ю.Дьяконов Е.Г. Метод разрезов для некоторых многомерных стационарных задач // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 1999. № 2. С. 9-16.
11. Волевич Л.Р.,Панеях Б.П. Некоторые пространства обобщенных функций и теоремы вложения // УМН. 1965. XX. Вып. 1 (121). С. 3-74.
Поступила в редакцию 07.06.04
УДК 517.958:533.9
В. А. Вознесенский, Д. Ю. Сычугов
МОДЕРНИЗИРОВАННЫЙ ЧИСЛЕННЫЙ МЕТОД РАСЧЕТА МГД-РАВНОВЕСИЯ ПЛАЗМЕННОГО ШНУРА В ТОКАМАКЕ
(кафедра автоматизации научных исследований факультета ВМиК, Российский научный центр "Курчатовский институт")
1. Введение. Задача расчета МГД-равновесия является классической и включается во всякий пакет программ, моделирующих процессы в плазме токамака. Первые работы в этом направлении появились свыше 20 лет назад [1, 2]. По мере развития эксперимента требования к таким программам ужесточались. Примером современного высокоразвитого кода может служить ЕР1Т [3], а также ТОКАМЕС] [4], применяемый при проектировании токамака Т-15М [5, 6].
В данной работе проведена модернизация классической схемы расчета МГД-равновесия. Создан код ТОРОЬ, не требующий подбора параметров для сходимости итерационных процедур. Результаты методических и физических расчетов показывают высокую эффективность нового метода.
2. Постановка задачи. Рассмотрим уравнения МГД-равновесия в токамаке (расстояние измеряется в метрах, ток в мегаамперах (MA), магнитное поле в теслах (Тл), магнитный поток в вольт-секундах (B.c.)) [7]:
±Д*ф ЕЕ #-г r or
1 дгр
г дг
+ = -0,8тг2<|
' jv{r, ф-фр), ф > фр, (в плазме),
JV
Ф{ 0,г) = 0,
lim ф(г, z) = 0.
Е hS(r -rk,z- zk), ф < фр, k= 1
r.z^-oс
Существуют следующие варианты постановки задачи (1) (в вариантах 1, 2 плазма и внешние токи считаются симметричными по г, в вариантах 3-6 — асимметричными).
1. Прямая задача для лимитерной конфигурации (плазма ограничена элементом конструкции). Величина фр задается соотношением фр = ф(гцт, гцт), где гцт, гцт — координаты касания плазменного шнура лимитером.
2. Прямая задача диверторного типа. Величина фр определяется соотношением фр = ф(г3,г3), где Гц, — координаты Х-точки сепаратрисы.
3. 4. Квазипрямая задача (лимитерного и диверторного типов). Все внешние токи полагаются заданными. Ищется такое положение плазменного шнура, при котором он уравновешен внешними токами.
5, 6. Квазиобратная задача. Ищется равновесие с наперед заданным положением центра шнура по г и по г, горизонтальным диаметром шнура, эллиптичностью (и, возможно, треугольностью). Величины внешних токов подбираются так, чтобы удовлетворить этим характеристикам.
3. Описание численного метода
3.1. Переход к задаче в ограниченной области. Выберем счетную область предположительно содержащую в себе плазменный шнур и не содержащую внешних токов. Пусть ф(г,г) — решение задачи (1). Значение функции ф на границе может быть вычислено по формуле объемного потенциала. Таким образом, задача (1) сводится к задаче в области с самосогласованным граничным условием (процедура Лакнера [1]).
3.2. Итерационный процесс. При решении задач 5 и 6 итерационно проводятся 4 вида вычислений.
1. Расчет граничного условия от плазмы.
2. Расчет тока в плазме.
3. Расчет потока, создаваемого током в плазме.
4. Вычисление геометрических параметров плазмы, подбор внешних токов.
Для создания экономичного кода самые трудоемкие расчеты должны проводиться минимальное число раз на внешнем цикле, а прочие — на внутреннем.
Код ТОКАМЕС] реализован как два вложенных цикла: на внешнем выполняется процесс 1, на внутреннем — 2, 3 и 4. Такой выбор обусловлен ограничениями оперативной памяти, существовавшими на момент разработки кода.
В новой схеме во внешнем цикле осуществляются процессы 1, 2 и 3, а во внутреннем — 4. Это представляется оптимальным для ЭВМ с большим ОЗУ.
3.3. Процедура Лакнера и расчет магнитного потока. В начале каждой внешней итерации по I для нахождения на потока фр1, создаваемого плазмой, применяется формула, следующая из (1):
1(Г' ап = // С(Г' Гр' Ф' ~ ^ ^ (2)
Здесь функция jlíp(r, г) — плотность тороидального тока в плазме, бг(г, г, г к) — функция источника кольцевого тока, расположенного в точке
С(г, г, гр1 гр) ее 0,8тг[(2 - к2) К (к) - 2 Е(к)} , (3)
где
2 4г-гр
(г + гр) + (г - гр)
К (к) и Е(к) — полные эллиптические интегралы первого и второго рода [8]. Поток от внешних проводников находится с помощью функции Грина (3).
Полный поток ф1'п для точек (г, г) 6 ^ в начале внутренней итерации по п ищется в виде суммы фр1 и суммарного потока от внешних проводников ф1е
ф1'"(г1г) = ф1^(г1г) + ф1р1(г1г). (5)
Для нахождения потока внутри счетной области решается уравнение Пуассона методом Ьи-разло-жения [9]. Значения сеточной функции источника и Ьи-разложение матрицы вычисляются единожды и хранятся в файлах.
3.4. Геометрические параметры плазмы и управление ими. У плазменного шнура ищутся самые верхняя и нижняя точки (г\'"р, а также концы максимального горизонтального
отрезка, лежащего внутри области плазмы г1^ахП), {г1£дМ, г1^ахВ).
По этим точкам вычисляются следующие управляемые геометрические параметры.
Координаты центра (гс,гс): г1с'п = 'е/< 2 г'яН*, г'с,п = <ор 2 Ыт .
Горизонтальный диаметр с11,п = — При явном задании величины горизонтального
диаметра ф1^п подбирается так, чтобы с11,п = б?треб-
Эллиптичность е: е =
ztop ~zbt 1
rright—rieft
Для достижения заданных геометрических параметров поток от внешних токов разделен на четыре части: фиксированный (в частности, от индуктора), управляющий положением по г (обычно от токов на внешнем обводе тора), управляющий положением по z и управляющий эллиптичностью (от токов сверху и снизу от плазмы). Таким образом,
ф1£(г, z) = ф}1Х{г, z) + т1;пфг{г, z) + т1^фг{г, z) + т1е'"фе(г, z). (6)
Множители mr, mz, те (фактически величины внешних токов) меняются в ходе итераций с помощью обратных связей.
3.5. Обратные связи. При выборе обратных связей следует учитывать, что производные зависимостей геометрических параметров шнура от mr, mz и те могут терпеть разрыв (например, при переходе от лимитерной конфигурации к сепаратрисной).
Алгоритм работает в двух режимах.
1. Режим накопления точек. На первой итерации случайным образом варьируется первый управляющий параметр, на второй — второй и т.д. Когда количество точек становится больше количества параметров и можно вычислять разностные производные, алгоритм переходит к режиму 2.
2. Режим секущих (хорд). В этом режиме осуществляется замещение накопленных точек т^ с
вычисленными в них значениями функций 1^(т^) новыми тк и 1^(т^). Вытесняются, во-первых,
точки, у которых Yi(mУ,(т> 0 для любого i, и, во-вторых, имеющие максимальное отклонение от предписанных значений.
Вычисляются аппроксимации первых производных управляемых параметров по управляющим, после чего делается шаг по методу Ньютона с использованием этих аппроксимаций.
Во избежание "уплощения" множества точек каждые 8 итераций значения функций в точках "забываются" и алгоритм переключается на режим 1.
3.6. Расчет профиля тока. В конце внешней итерации происходит пересчет профиля тока jl (г, z) по вычисленному потоку ф1(г,г) в области (г, z) £ Qlp по формуле
R2
jl = \1и(г,Ф> - ФР) = \l{ßr {ф1 - фРУ' + (1 - ß)— {ф1 - фРУ2 \ , (7)
где /3 — "бета токовая"; 71, 72 — показатели распределения плотности тока;
Л г (ф1 — фр)11 (1г(1г
- 9 а R2 ЕЕ
Я 7 W ~ ФрУ'2 drdz
Slp
среднее по шнуру значение г2; А' вычисляется исходя из заданной величины полного тока I в плазме:
А1 [[ и(г,ф(г,2)-фр)с1гс12 = 1. (8)
4. Методические вычисления
4.1. Сравнительная верификация кода. Задача расчета МГД-равновесия плазмы нелинейна и в общем случае может быть решена только численно. Это значит, что обязательным тестом является сравнение результатов вычислений с другими кодами.
1,0 1,5 2,0 2,5
R
Рис. 1. Пример линий уровня ф = const и аппроксимации лимитера в токамаке Т-15М
На рис. 1 представлены линии уровня ф = const для момента одного из сценариев разряда токамака Т-15М г = 2,06 с. Расчеты проводились по кодам T0KAMEQ и T0P0L на сетке 40 X 60 при одинаковых входных величинах.
Визуально линии уровня, полученные с помощью обоих кодов, отличить невозможно, относительная норма разности в ¿2 составляет 0,000328. Аналогичные вычисления были проделаны для моментов т = 0,03 си г = 0,36 с.
Приведем результаты сопоставления кодов для момента г = 2,06 с в зависимости от числа узлов сетки. В таблице представлена относительная норма в ¿2 разности потоков, расчитанных по двум кодам в зависимости от п.
Размер сетки 30 X 45 40 X 60 50 X 75
W^iopol ^ïofcameg || 0,0211 0,000328 0,000232
IIФtovol II
4.2. Сравнение скорости вычислений. Важнейшей характеристикой метода является скорость вычислений. В таблице приведено время расчетов для стадий начала разряда, подъема тока и стационарного этапа по кодам T0KAMEQ и T0P0L для различного числа узлов на сетке. Расчеты проводились на ПЭВМ с процессором Intel Pentium III с тактовой частотой 800 Mhz под управлением ОС Linux. Видно, что новый алгоритм работает быстрее старого.
Сетка Момент 0,03 с Момент 0,36 с Момент 2,06 с
30 X 45 5,51 с/0,78 с 6,47 с/1,06 с 4,78 с/0,79 с
40 X 60 11,07 с/1,52 с 11,82 с/1,86 с 8,53 с/1,41 с
50 X 75 17,29 с/2,52 с 19,75 с/3,07 с 13,18 с/2,39 с
4.3. Сходимость метода хорд. На рис. 2 представлен пример хода подстройки гс при фиксированном I для момента г = 0,36 с для старого и модернизированного методов. Итерации 1, 2, 9 и 10 проходят в режиме 1. Видно, что в режиме 2 сходимость постепенно ухудшается.
10"1 г
10"3
10"5
Ю-7
Ю"9 -
Ю-111—_,_,_,_,_,
0 5 10 15 20 25
Рис. 2. Зависимости отклонений гс от предписанного значения в ходе итераций по п
5. Выводы. На основе модернизации стандартного алгоритма расчета асимметричного равновесия разработан код, обладающий следующими двумя преимуществами по сравнению с кодом Т0КАМЕС].
1. Не требует методических вычислений по подбору параметров, управляющих итерационным процессом.
2. Работает существенно быстрее.
СПИСОК ЛИТЕРАТУРЫ
1. Lackner К. Computation of ideal MHD equilibria // Computer Physics Communications. 1976. 12. N 1. P. 33-44.
2. Попов A.M., Сычугов Д.Ю. Равновесие и устойчивость тороидального плазменного шнура с дивертором // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 1980. № 3. С. 32-37.
3. Lao L.L.,John H.St.,Stambaugh R. D., К el 1 m an A.G.,Pfeiffer W. Reconstruction of Current Profile Parameters and Plasma Shapes in Tokamaks // Nuclear Fusion. 1985. 25. N 11.
4. Вознесенский В.А., Гасилов H.A., Днестровский Ю.Н., Кузнецов А.Б., Сычугов Д.Ю., Цаун С. В. TOKAMEQ - код для расчета равновесия плазмы в токамаке. Препринт РНЦ "Курчатовский институт". М., 2001.
5. Бондарчук Э.Н., Вознесенский В.А., Днестровский Ю.Н., Леонов В.М., Максимова И. И., Сычугов Д.Ю., Цаун С. В. Вертикальная МГД-устойчивость плазмы в Т-15М // Вопросы атомной науки и техники. Сер. Термоядерный синтез. 2003. Вып. 2. С. 55-60.
6. Bondarchuk E.N.,Dnestrovskij Yu. N.,Leonov V.M.,Maksimova I.I.,Sychugov D. Yu., Tsaun S. V., Voznesensky V.A. Vertical MHD stability of the T-15M tokamak plasma // Plasma Devices and Operations. 2003. 11. N 4. P. 219-227.
7. Шафранов В.Д. Равновесие плазмы в магнитном поле // Вопросы теории плазмы. Т. 2. М.:
Госатомиздат, 1963.
8. Thacher Н.С., Jr. Complete Elliptic Integrals // Communications of the ACM. 1963. Apr. 6. P. 163.
9. Stewart D.E., Leyk Z. Meschach: matrix computations in С // CMA Proceedings. 32. Canberra: Australian National University, Australia, 1993.
Поступила в редакцию 01.10.04
УДК 517.958:535.14
А. Г. Волков, В. А. Трофимов
О РОЛИ СПЕКТРАЛЬНОГО ИНВАРИАНТА
ПРИ КОМПЬЮТЕРНОМ МОДЕЛИРОВАНИИ НЕЛИНЕЙНОГО
РАСПРОСТРАНЕНИЯ ФЕМТОСЕКУНДНЫХ ИМПУЛЬСОВ1
(кафедра вычислительных методов факультета ВМиК)
1. Введение. Как известно [1, 2], анализ распространения фемтосекундного лазерного импульса основан на так называемом комбинированном (обобщенном) нелинейном уравнении Шредингера (КНУШ), которое отличается от традиционного нелинейного уравнения Шредингера (НУШ) наличием производной по времени от нелинейного отклика среды. Ее присутствие предъявляет новые требования к построению разностных схем для данного класса уравнений, которые обусловлены новыми (по сравнению с традиционными НУШ) свойствами КНУШ. Так, производная от нелинейного отклика приводит к зависимости скорости распространения волны от ее интенсивности. Результатом этого является возможность формирования оптических ударных волн [1, 2].
Другие свойства КНУШ удалось выявить, используя подход [3, 4], развиваемый авторами для подобных задач. Так, построенные инварианты показали присутствие "особой" спектральной гармоники в начальном условии для уравнения КНУШ: ее амплитуда изменяется по определенному закону, т.е. в процессе распространения светового импульса имеет место спектральный инвариант. Очевидно, данное свойство необходимо учитывать при построении разностных схем для уравнения КНУШ. Именно влияние точности сохранения спектрального инварианта на получаемое при компьютерном моделировании решение исследуется в настоящей работе. При этом рассмотрение проводится на примере нескольких схем, в том числе и консервативной схемы.
2. Основные уравнения. Процесс распространения фемтосекундного импульса в среде с кубичной нелинейностью с учетом производной по времени от нелинейного отклика среды (дисперсии нелинейности) при отсутствии влияния дифракции оптического излучения (что реализуется либо при его распространении в оптическом волокне, либо при длине среды много меньшей дифракционной длины светового пучка) и с учетом дисперсии второго порядка описывается следующим безразмерным комбинированным нелинейным уравнением Шредингера (КНУШ):
дА д2 А д
— + 10— + 1а\А\2А + а1—(\А\2А) = 01 ¿>0, 0 < I < 1и (1)
£ = 0,5(А ехр(шЬ — 1кг)) + к. с. с начальными и граничными условиями
А|г=0 = А0(;), (2)
Здесь £ — напряженность электрического поля; и и к соответственно безразмерная частота и волновое число светового импульса; А(г,£) — нормированная на максимальное значение на входе в нелинейную
1 Работа выполнена при частичной финансовой поддержке РФФИ (грант № 02-01-727).
