Описание расчетных моделей и используемых численных методов

Для определения карты поля по магнитным измерениям необходимо решить уравнение относительно токов i_j(Γ). Для измерения внешнего потока используются 24 петель, расположенных на ВК как показано на Рис. 10 зелеными точками.

Алгоритм нахождения токов i_j(Γ) по заданным измерениям потока (в программном пакете он реализован как Режим 1) состоит в следующем.

В полоидальной плоскости r, z строится замкнутый контур, состоящий из отрезков прямых, соединяющих точки пересечения соседних петель плоскостью r, z.

Каждый отрезок делится на n=20*ll/L+1 равных частей (численный коэффициент 20 может быть и другим).

ll - длина отрезка (расстояние между соседними петлями),

L - расстояние от петли до геометрического центра контура.

Получившиеся отрезки (i) принимаются за сечения плоскостью rz токовых элементов (ТЭ), которые имеют форму боковой поверхности усеченного конуса с постоянной по сечению плотностью поверхностного тока, текущего нормально к плоскости r, z.

Т. к. замкнутый контур (см. Рис. 10) имеет внутри ТЭ (пассивные витки, обозначаемые как ps), то теория п. 2 в этом случае неприменима. Поэтому поток от ТЭ ps (координаты которого известны, а ток измеряется) вычитается из измеренных магнитными петлями потоков. Значение потока Ps на отрезках (i) вычисляется с помощью линейной интерполяции получившегося потока на петлях по расстоянию между отрезками, т. е.:

где lm и lp расстояния ТЭ "i" от петель l и l+1 соответственно, между которыми расположен ТЭ "i".

Обозначим через G(a,b) коэффициент взаимной индукции ТЭ " a " и " b ". Тогда для нахождения токов отрезков (k) необходимо решить систему линейных уравнений с симметричной, положительно определенной матрицей G(i,k):

(3.2.2)

Здесь индексы i, k пробегают все ТЭ, построение которых было описано выше.

В программном пакете для вычисления коэффициентов матрицы G(i,k) (коэффициентов взаимной индукции ТЭ " a " и " b ") используется следующая формула:

(3.2.3)

где интегрирование проводится по объему V_a ТЭ " a " и по площади S_b сечения ТЭ " b " какой-либо полоидальной плоскостью rz.

В формуле (3.2.3) использованы следующие обозначения:

f_ab - угол между направлением тока в ТЭ " a " в точке интегрирования и нормалью к плоскости rz;

dS_b - элемент площади сечения ТЭ " b " плоскостью rz (для ТЭ с поверхностным током элемент площади вырождается в элемент отрезка);

r_b - его расстояние от оси системы;

dV_a - элемент объема ТЭ " a " (для ТЭ с поверхностным током элемент объема вырождается в элемент поверхности);

j(a), j(b) - плотности токов в точках интегрирования в элементах " a " и " b ";

r_ab - расстояние между элементами интегрирования dV_a и dS_b;

I_a, I_b - полные токи в ТЭ " a " и " b ".

Коэффициент взаимной индукции для далеко расположенных друг от друга ТЭ находится разложением по размерам элементов, а для близко расположенных или совпадающих (коэффициент самоиндукции) ТЭ находится аналитически. Там, где не удается получить удовлетворительный аналитический результат, используется численное интегрирование. Этот прием позволяет избежать сингулярности в подынтегральном выражении. Численное интегрирование проводится рекурсивно применяемым методом Симпсона (более подробно см. файлы F.C (функция _L) и H.C).

При решении системы линейных уравнений (3.2.1) используется метод разложения матрицы G в произведение двух треугольных матриц (нижней L и верхней L', G=L'L, где L' - транспонированная матрица) с положительными диагональными элементами в силу положительной определенности матрицы G. Решение системы уравнений реализуется функцией sol (файл MK.C).

В результате решения системы (3.2.2) находятся токи отрезков j(k) выбранного контура. Искомый поток в точке с координатами r,z внутри выбранного контура находится по формуле:

(3.2.4)

Здесь индекс суммирования k пробегает все ТЭ c найденными из системы линейных уравнений (3.2.1) токами. Т.к. при нахождении токов j(k) не учитывалось влияние тока пассивных витков внутри контура (элемент ps), то поле от этого тока необходимо добавить для получения суммарного поля. Поэтому в формуле (3.2.3) добавлен член G(r,z|ps)I(ps), который определяет поле в точке с координатами r,z от тока I(ps) пассивных витков (оба коэффициента G(r,z|k) и G(r,z|ps) вычисляются согласно (3.2.3)).

Компоненты вектора магнитной индукции получается с помощью аналитического дифференцирования этого потока по координатам r, z.

Поток и магнитная индукция ТЭ, расположенного далеко от точки наблюдения, находится численным интегрированием методом Гаусса с числом точек, зависящим от расстояния ТЭ от точки наблюдения. Для близко расположенного ТЭ - аналитическим интегрированием (более подробно см. файлы F.C (функция PS) и H.C)

Для определения карты магнитного поля по заданным токам обмоток I(l) (в программном пакете - Режим 2) поток в точке с координатами r, z находится по формуле:

(3.2.5)

где индекс l пробегает все ТЭ с заданными токами в них, а величины G(r,z|l) вычисляется согласно (3.2.2). Компоненты вектора магнитной индукции получается с помощью дифференцирования этого потока по координатам r, z.

Рисунок 10 - Положение полоидальных обмоток и петель магнитной диагностики токамака КТМ.

Проведение расчетов

4.1 Расчет карты магнитных полей плазмы с помощью PlasmalessTokScen

Для расчета магнитных полей плазмы на токамаке используется ряд специалных программ. Здесь расчет велся с помощью PlasmalessTokScen. Первым делом с помощью Plasmaless открываем сценарий пробоя плазмы на токамаке. В начале расчета необходимо ввести в программу следующий набор входных данных:

1) параметры датчиков ДНО

Параметры данных датчиков напряжения обхода, вводятся во вкладке «loops» диалогового окна «Edit window» (рисунок 11).

В данном окне вводятся следующие параметры:

Рисунок 11 - Окно редактирования данных вкладка параметров датчиков напряжения обхода

2) параметры датчики зондов

Параметры датчики данных зондов, вводятся во вкладке «probs» диалогового окна «Edit window» (рисунок 12).

В данном окне вводятся следующие параметры:

Рисунок 12 - Окно редактирования данных вкладка параметров датчиков зондов

После редактирования данных программы можна начать расчет. Расчет будет веcти до временного шага t = 0.0051c. После расчета берем карту магнитного поля и данные с зондов в окне «additional functions» - «show eddy currents table». Данные вводим в таблицу.

Рисунок 13 – карта магнитного поля с Plasmaless во временном шаге t=0.0051c

Таблица 3 – данные компонентов поля (Т) с зондов

В меню «additional functions» реализованы команды для просмотра таблицы индуктивностей для датчиков, наведенных на них, таблиц магнитных полей по координатам, заданным пользователем и др. И из окна «show eddy currents table» берем данные во временном шаге t=0.0051c по магнитным потокам на датчиках ДНО. Эти данные используются в расчетном коде для восстановления карты магнитного поля плазмы.

Таблица 4 – магнитные потоки на датчиках ДНО