CC BY
38
С.Н. Андрианов, С. А. Артамонов, А. Н. Дубровин, В. А. Елисеев
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЗБ-МАГНИТНОГО ПОЛЯ ИЗОХРОННОГО ЦИКЛОТРОНА В ПИЯФ
1. Введение. В Петербургском институте ядериой физики им. Б. П. Константинова Российской Академии наук (ПИЯФ РАН) ведется строительство Гатчинского Изохронного Циклотрона (ГИЦ) по ускорению Н--ионов до энергии 75-80 МэВ с током выведенного пучка около 100 мкА [1]. Одной из важнейших проблем при создании циклотрона является задача формирования требуемого изохронного поля только с помощью железных масс, т. е. без кольцевых подстроенных катушек. Широко применяемая ранее процедура последовательного приближения к такой изохронной зависимости за счет неоднократного изготовления шиммов постепенно уходит в прошлое. Дело заключается в том, что большая стоимость проектируемых физических установок выдвигает в последнее время на первый план процедуру их численного моделирования. Прогресс в развитии вычислительной техники дает возможность существенно удешевить решение этой задачи, а сама процедура численного моделирования позволяет заранее предсказать основные свойства и характеристики проектируемого объекта. Кроме того, такой подход приводит к возможности изготовления необходимого «железа» лишь на последнем этапе строительства, что существенно сокращает стоимость установки, а также время на ее ввод в эксплуатацию. Поэтому задача компьютерного моделирования трехмерного поля ГИЦ представляется своевременной и актуальной.
2. Постановка задачи, основные уравнения и методы. Сложность поставленной задачи в первую очередь определяется существенной нелинейностью проблемы, которая может быть преодолена при помощи созданной в Новосибирском институте ядерной физики программы MERMAID [2], позволяющей вычислять трехмерные статические магнитные поля при задании тех или иных граничных условий.
2.1. Конечно-элементные уравнения, прямые треугольные призмы. Для дискретизации уравнений Максвелла стационарной магнитостатики
в программе используется эквивалентная вариационная постановка задачи. Однако в трехмерном случае формулировка задачи в терминах векторного потенциала (как в двумерном случае) слишком громоздка для решения. В то же время формулировка в терминах скалярного потенциала требует выделения роторной части поля. Представим поле в виде Н = Н' I ’Уу\ где [V х Н'] = Подставляя введенные обозначения в (1), получим следующее уравнение для скалярного потенциала: = У^Н'. Для
численной реализации решения сформулированной выше задачи требуется осуществить дискретизацию последнего уравнения. Как известно, приведенные уравнения могут быть получены минимизацией интеграла:
[V х Н] = j, VB = 0.
(1)
© С.Н. Андрианов, С. А. Артамонов, А.Н. Дубровин, В. А. Елисеев, 2008
где р = У/хН'. Необходимым условием минимума, как известно, является равенство нулю первой вариации этого интеграла:
,)1 = I (ГОВ + (Ш(Ш - 2(хр) (IV = 0.
Разобьем магнит на тетраэдры. Обычно предполагается, что поле В внутри тетраэдра постоянно. В этом случае потенциал будет изменяться линейно, тогда можно написать Н = (р1<тг, где гг' = к'/61р. н к' вектор, перпендикулярный г-й грани тетраэдра и с длиной, равной площади соответствующей грани. Нетрудно видеть, что имеет место равенство НВ = (здесь мы воспользовались правилом суммирования
Эйнштейна).
Прежде чем вычислять вариацию интеграла, обсудим форму конечных элементов. Наиболее общим является разбиение исследуемой области на произвольные тетраэдры, но, в отличие от двумерного случая, в трехмерном моделировании форма элементов существенным образом влияет на эффективность решения системы уравнений. Для того чтобы сохранить возможность пользоваться высокой эффективностью циклического метода сопряженных градиентов [3], необходимо реализовать разбиение с использованием прямых призм. При этом возникает связь с переменными только вдоль оси призмы (для определенности пусть это ось 2). Если провести разбиение переменных по линиям, например, по оси У, то в плоскости ХЕ получим картину связей, аналогичную двумерному случаю с прямоугольной сеткой.
В нашем случае мы можем использовать так называемое красно-черное разбиение переменных - т. е. такое разбиение переменных на две группы, при котором связи внутри группы разрешаются явно. В случае треугольников подобное разделение переменных достигается разбиением по линиям: четные и нечетные линии сетки образуют соответственно красные и черные группы переменных, связи же внутри линий легко разрешаются явно, поскольку представляют систему уравнений с трехдиагональной матрицей. Следовательно, можно реализовать описанное красно-черное разбиение, меняя четные и нечетные переменные в каждой строке по
Покажем, каким образом возникает связь переменных только вдоль оси Для этого возьмем в качестве поля в призме его среднее значение на образующих призму тетраэдрах. В этом случае в вариационных уравнениях возникают перекрестные связи, потому для формулировки вариационной задачи достаточно знать только выражение квадрата поля через потенциал. Запишем его в виде среднего квадрата поля по образующим тетраэдрам:
1
У' Н2
2 / у тетр'
тетр
Для каждого отдельного тетраэдра вариационные уравнения дают связи вида
Мг] = гг'гг'.
следовательно, из-за ортогональности соответствующих векторов гг' связи «не вдоль оси» Е будут равны нулю. Отсутствие связей не вдоль оси, кроме уменьшения числа арифметических действий в расчете на один узел, снижает также требования к памяти компьютера, что дает возможность увеличить число узлов и таким образом частично скомпенсировать неточность ступенчатой аппроксимации по оси. Несмотря на то, что можно получить явное выражение для матрицы метода Ньютона, в эту матрицу входят
перекрестные связи, что не позволяет использовать эффективный метод Ньютона одновременно с циклическим методом сопряженных градиентов. Необходимо заметить, что наличие перекрестных связей приводит также к усложнению распараллеливания описываемого алгоритма.
2.2. Метод бисекции. Как было сказано выше, необходимо произвести расщепление поля на градиентную и роторную части. Простым и эффективным способом конструирования роторной части поля является метод бисекции [4]. Наиболее просто его формулировка выглядит в случае топологически регулярной сетки. Для ограничивающего область расчета куба всегда можно задать значения циркуляции вдоль каждого его ребра так, чтобы для каждой грани куба выполнялось условие
где I - суммарный ток через грань. Рассечем теперь куб одной из плоскостей неиз-На получившихся новых ребрах можно задать значения циркуляции, чтобы выполнялось то же соотношение (2) для новых появившихся граней. Таким образом, просто перераспределяем циркуляцию на каждом из рассеченных ребер с сохранением суммы. Требуемые значения циркуляции получаются теперь автоматически. Можно показать, что соотношение циркуляции выполняется одновременно для всех граней в случае консервативности тока, т. е. суммарный ток через любой куб должен быть равен нулю. Получаемые значения Н', вообще говоря, не малы в областях с железом, что плохо для нелинейных итераций, поскольку на начальных итерациях железо будет перенасыщено. Однако поле Н' можно сделать малым в областях с железом следующим образом. Сначала делаем бисекцию, как описано выше. Затем рассчитываем этот магнит, заменяя «нелинейное железо» на железо с постоянным, но таким большим /л, что его можно считать бесконечным. Тогда, очевидно, получившееся поле Н' в железе будет нулевым, поскольку поле В конечно. Однако данный путь закрывает возможность решения таких задач, в которых поле Н' стремится к бесконечности при неограниченном росте ц,. Однако, к счастью, это достаточно узкий класс задач, в которых существует путь вокруг тока, полностью проходящий по железу (задачи типа расчета трансформаторов).
2.3. Расчет распределения токов. Строгое соблюдение консервативности токов может оказаться непростой задачей при вводе в программу их распределения, особенно для катушек, выполненных в виде одного толстого витка сложной формы. Поэтому представляется вполне естественным сделать ввод только геометрии катушек, а распределение токов вычислять. Фактически необходимо рассчитать распределение электрического потенциала внутри катушки с фиксированным (постоянным) потенциалом на концах и отпущенным на остальной части поверхности, поскольку компонента тока, перпендикулярная поверхности, равна нулю. Итак, представим плотность тока в виде j = оЧц), где а - проводимость (будем считать ее независимой от плотности тока). Тогда уравнения для потенциала можно записать в виде
(2)
((ТпрР]^^) УПр = 0.
ПР
Или, используя соотношение j = ГТу ,(т‘ и то, что гг'Г ос к', получим
£/=().
ПР
Здесь I - токи через те грани малых призм, которые находятся внутри больших призм. Тем самым найдено консервативное решение, но на смещенной сетке. Предстоит «переместить» решение на исходную сетку, сохранив консервативность и добавив требование на равенство нулю тока через поверхность. Выполнение последнего требования дает важное свойство - сохраняется полный ток через любое сечение катушки.
Рассмотрим двенадцать маленьких призм, окружающих данный узел. Найдем такие значения токов через их грани, чтобы суммарный ток через поверхность каждой маленькой призмы был нулевым. Тогда, очевидно, будет нулевым и суммарный ток через поверхность каждой большой призмы, что и нужно. Сначала покажем, что точное решение существует. Для этой цели реализуем следующие последовательные шаги:
1. Зададим произвольные значения токов для пяти граней, лежащих в плоскости.
2. Зададим произвольный ток для верхней вертикальной грани на границе с оставшейся неопределенной грани в плоскости.
3. Из сохранения суммарного потока вычислим по порядку все токи верхних вертикальных граней.
4. Определим ток для оставшейся неопределенной грани в плоскости.
5. Зададим произвольный ток для нижней вертикальной грани.
6. Вычислим по порядку все токи нижних вертикальных граней.
Остается показать сохранение тока для последней маленькой призмы. Ясно, что существует такой суммарный ток через лежащие внутри большой призмы грани, при котором будет сохраняться ток для данной маленькой призмы. Следовательно, ток будет сохраняться для всей фигуры, состоящей из малых призм, в том числе будет равен нулю суммарный ток через поверхность этой фигуры. Но таким же свойством обладает решение на исходной сетке. Значит, суммарный ток через лежащие внутри большой призмы грани совпадает с тем, который нужен для сохранения суммарного тока малой призмы.
Аналогично можно доказать существование точного решения, когда отсутствует ряд призм сверху и/или снизу, но так, чтобы отсутствующие призмы шли подряд и верх с низом были связаны хотя бы через одну грань. На токи на границе с отсутствующими призмами накладывается строгое требование равенства их нулю. В то же время можно показать, что в общем случае нет точного решения для конфигурации, когда существует участок нулевой толщины внутри тока. На таких участках нужно ввести либо ненулевую толщину, либо явный разрыв.
Итак, при соблюдении некоторых, весьма естественных ограничений, есть точное решение задачи, более того, существует некоторый произвол в его выборе. Естественно выбрать то точное решение, которое наиболее близко к нашему приближенному решению. Поскольку в качестве такого приближенного решения нужно взять некоторую комбинацию из токов двух призм, имеющих общую данную грань, то удобно взять комбинацию токов с весами, обратными проводимости:
_ (1/0-11] + 1/02/2) _ (о2/х + а1 Ь)
1сред ~ (1/(7! + 1/(72) “ (01+02) '
Тогда для границы тока, если возьмем для призмы вне тока 1\ и ох, получим 1сред = О, что и требовалось. В то же время для призмы с током, лежащей за границей задачи, можно положить произвольный ток и очень большое значение проводимости, тогда
/сред = Ь- Таким образом, получаем совместные линейные уравнения для токов каж-
дой из маленьких призм вида
VI = - V
/ ^1внутр / ^1внешн •
V (I -1° ) = - V (1° - Ґ
/ и 1/внутр 1внешн/ / ^ ^внутр 21
внешн /
Нужно найти ближайшее к нулю решение этого уравнения. Как известно, оно может быть получено с использованием метода наименьших квадратов: для линейного уравнения МХ = В ближайшее к нулю решение ищется как решение линейной системы уравнений
(М*М + А) X = М*В,
где А - диагональная матрица с положительными, стремящимися к нулю элементами Агг а1 > 0. Заметим, что таким образом минимизируется сумма ^ ((Д/''.г; Ь') 2 | + агх2). Обсудим численные аспекты выбора элементов матрицы А. Если выбрать элементы а1 слишком большими, то исказится исходное уравнение, и решение будет не консервативно. При слишком малых элементах а1 решение будет далеко от оптимального. Поскольку все элементы а1 порядка единицы, разумным будет выбор как две трети разрядного диапазона между 1 и минимальным числом. Теперь ясно, что, выбрав а1 к, 1 для граней на границе тока, получим нулевые значения токов с достаточной точностью.
2.4■ Метод нелинейных итераций с релаксацией. Итак, необходимо решить систему уравнений вида ,4Х = Г с нелинейным оператором А. Поскольку мы отказались от метода Ньютона, для его решения воспользуемся общим итерационным методом, который запишем в виде [5]
В^к+1 - Yk)/т + ЛYk = ¥.
Скорость сходимости итераций определяется близостью оператора В к А', т. е. матрице метода Ньютона. Вводя операторную норму А = ||Д||, для наших целей можно записать это утверждение следующим образом.
Утверждение. Если 71В < А' < 72В, то при т = 2/(71 + 72) итерации сходятся, и в некоторой операторной норме ошибка за одну итерацию убывает со скоростью геометрической прогрессии с показателем р = (72 — 71) / (72 +71)-
Выпишем в явном (матричном) виде вклад в оператор А' от одной призмы:
А'Ч = ц ^ гг'гг' + 2ц' ^ Нгг' ^ Нгг'.
тетр тетр тетр
Тогда А'^угу] = цЬ2 + 2/л' ^ НЬ ^ НЬ. Отсюда ц + 2//'Н2 < А' < ц при ц' < 0 и ц <
тетр тетр
А' < /л + 2//'Н2 при р! > 0. Матрица В оператора В выбирается таким образом, чтобы не вызвать заполнения ненулевыми элементами, т. е. сохраняем структуру исходной матрицы. Для того чтобы параметр г не зависел от конкретной призмы, выберем матрицу В с пропорциональными р + р'Ш2 элементами, тогда 71+72 = р/ (р + р'Н2) + {р + р'Ш2) = 2 (это справедливо при любом знаке р) и г = 1.
Итак, полное выражение для матрицы
в1;> = (р + ^'Н2) V ^2
тетр
Оценим скорость сходимости нелинейных итераций: р к, 1 + 2 {р + р'Н2). Здесь мы использовали малость фактора Л = (р + 2//Н2) //i, равного логарифмической производной поля В = |В| по Н = |В|: Л = (dB/dH) / (В/Н). Практически фактор Л мал при поле В от нескольких килогаусс до индукции насыщения железа. Для низкоуглеродистого железа в минимуме Л фактор составляет 0.1-0.15, но для ваннадиевого пермендюра достигает 0.01.
Обладая важным свойством глобальной сходимости, изложенный выше метод все же заметно уступает по скорости методу Ньютона, имеющего, как известно, квадратичную скорость сходимости вблизи решения. Тем не менее данный метод вполне подходит для начальной стадии итерационного процесса. Он имеет важное достоинство хорошей скорости сходимости итераций для сильно насыщенных областей железа, в то время как точное значение проницаемости на ненасыщенных участках, где он сходится хуже, не важно для решения вне них самих.
Еще одним достоинством изложенного метода является его вычислительная простота, поскольку время вычисления нелинейной матрицы становится сравнимым со временем решения линейных уравнений с этой матрицей. Заметим также, что степень обусловленности матрицы В не более чем вдвое превосходит степень обусловленности исходной матрицы А. И еще одно замечание: метод Ньютона менее эффективен для скалярного потенциала, чем для векторного, поскольку степень обусловленности матрицы метода больше как минимум в 1/Л раз. Мы имеем хорошее начальное приближение к решению, полученное благодаря решению линейной задачи и последующим уточнением с помощью некоторого числа итераций метода релаксации. Поэтому в качестве следующего шага необходимо разработать эффективный метод, который необязательно обладает свойством глобальной сходимости. За основу возьмем локально оптимальный метод сопряженных градиентов. Вообще говоря, он применяется для решения линейных задач, но считаем, что найденное приближение находится в окрестности точного решения, где проблему можно считать в некотором смысле линейной; более точно, вариационный интеграл можно считать квадратичным и, следовательно, А' рз const. Фактически, будем решать систему уравнений метода Ньютона
A' (yk+i ~ У к) = F - Ayfc
с помощью итерационного процесса
Ук+i = У* + (1 - Ctk) (у*_1 - у к) - акТкМ^Тк,
в котором г к = А' (у к — Ук-i) + Ay*._i — F. Здесь матрица В взята той же, что и в методе релаксации. В известном смысле применяется процедура ускорения этого метода. Наложив условие минимума погрешности по норме А', получим выражение для ак и Тк через скалярные произведения а*. = (w*.,r*.), 6*. = (w*.,r*._i), с*. = (r*.,y*._i), dk = (гк \-Ук ~ Ук-i), е* = (A'wfe, wfe):
ак = ((ак ~ bk) Ьк - dkek) / ((ск ~ dk) ек - (ак - Ь*)2) ,
тк = ак/ек + (1 - ак)/(акЬкек) ■
Заметим, что можно избавиться от г*_1, полагая r*_i = г* + А' (у*-1 — у*), тогда можно записать: bk = ак - (A'wfe,yfe - yfe_i) и dk = ек - (А' (yfe - yk-i) ,Ук ~ Ук-i)-Расчетная процедура выглядит так:
1) решается методом сопряженных градиентов система уравнений Bw*. = г*.;
2) вычисляются скалярные произведения а* и с*;
3) вычисляются скалярные произведения типа (A'w*.,y* — yk-i) и тут же r*+i;
4) определяются оставшиеся параметры и вычисляются у*+1.
Для скалярного произведения (А'а, Ь) получим
(Л'a. I)) = I'^T ab I 2//' ^ На ^ НЬ. а = {</,......)\ Ь = {/>,...)\
тетр тетр тетр
Можно показать, что требование минимума погрешности через одну итерацию приводит к методу, совпадающему с методом сопряженных градиентов при некотором выборе начального приближения. В то же время традиционный метод сопряженных градиентов основан на выборе последовательности взаимно перпендикулярных векторов. Однако из-за нелинейности матрицы перпендикулярность быстро нарушится, и сходимость замедлится или даже может появиться расходимость. Локально же оптимальный метод подстраивает свою последовательность векторов под изменяющуюся матрицу и более устойчив. Практика показывает, что комбинация методов релаксации и сопряженных градиентов позволяет решить большинство нелинейных задач за 25-30 итераций до точности 10”“4 10”“5.
Заметим, что хорошую оценку достигнутой точности дает разница решения между двумя нелинейными итерациями. Точный критерий остановки итераций существенен, поскольку данный метод решения не обладает, вообще говоря, свойством квадратичной сходимости.
2.5. Особенности магнитной структуры ГИЦ. Другой, и не менее важной, проблемой является чрезвычайная сложность магнитной структуры ГИЦ, которая обусловлена большой спиральностью секторов.
На рис. 1 видно, что каждый сектор ГИЦ простирается на три квадранта, что приводит к невозможности использования широко применяемого стандартного приема при расчете трехмерных полей. Прием заключается в том, что для расчета подготавливается, как правило, 1/8 часть магнитной структуры с заданием соответствующих граничных условий. Более того, если у секторов отсутствует спиральность, т. е. сектора прямые, как, например, в циклотроне У-400 ЛЯР ОИЯИ, то для расчета достаточно использовать 1/16 часть магнитной структуры. Это обстоятельство позволяет существенно сократить число узлов трехмерной сетки и, тем самым, время вычислений того или иного варианта проектируемого трехмерного магнитного поля.
Следует также отметить, что все программы для вычисления трехмерных магнитных полей ориентированы на использование указанных симметрий задачи, поэтому в них для расчета необходимо задавать незамкнутую токовую катушку.
Сложная магнитная структура ГИЦ не позволяет воспользоваться указанными преимуществами, что приводит к существенным трудностям и дополнительным требованиям к программе трехмерных вычислений. Стандартный вариант программы MERMAID разрешает задавать в медианной плоскости (х, у) сетку размером не более чем (223 х 223) узлов. Поскольку магнитная структура ГИЦ имеет только полюс диаметром 2.05 м, то задача достижения точностей, необходимых для формирования изохронного расчетного поля, становится более чем проблематичной. Здесь уместно напомнить, что магнит ГИЦ простирается в длину примерно на 5.7 м, а в ширину - приблизительно на 2.5 м. Кроме того, для правильного расчета магнитного поля необходимо задавать вокруг магнита довольно большое воздушное пространство. Это позволяет, установив
Рис. 1. Секторы ГИЦ (проекция Рис. 2. Оптимальное расчетное трех-
на медианную плоскость). мерное поле ГИЦ.
на его границах соответствующие граничные условия, добиться практически полного отсутствия влияния внешних границ на поведение магнитного поля хотя бы внутри рабочей области.
Указанные трудности не закрывают путь к компьютерному моделированию трехмерного магнитного поля ГИЦ. Они всего лишь предупреждают о «подводных камнях» сформулированной проблемы, намечают пути к их преодолению и требуют грамотной интерпретации полученных результатов вычислений. Сама же процедура представляется своевременной и современной.
3. Основные результаты компьютерного моделирования. Для моделирования трехмерного магнитного поля ГИЦ была применена как стандартная, так и нестандартная версия MERMAID. Это дает возможность задавать в (ж, 2/)-плоскости максимальное число узлов: (400 х 400). Данное обстоятельство привело к существенному повышению точности расчетов. Поскольку каждый сектор магнитной структуры ГИЦ простирается на три квадранта, то минимальным объектом, дающим возможность использовать симметрию задачи, является половина магнита. При задании геометрии магнитной структуры требуется все описывающие ее линии спроектировать на медианную плоскость, т. е. предусмотреть заранее.
Поэтому первым этапом работы являются вопросы выбора шага сетки, отдаленности границ области для выставления соответствующих граничных условий, разбиение основной области на подобласти, где наиболее часто будут вноситься изменения в геометрию задачи, и тому подобные операции. Следующий этап заключается в разработке пространственной модели магнитной системы ГИЦ, максимально приближенной к реальной геометрии. Эта задача решается за счет детального описания геометрической формы элементов конструкции магнитопровода, катушек, границ раздела различных сред и нелинейных магнитных свойств применяемых электротехнических
сталей. И только после всех операций, используя тс или иные заложенные и структуру модели и необходимые кривіле намагничивания, можно провести расчет и получить трехмерное магнитное поле циклотрона (рис. 2). Стандартная версия MERMAID позволила использовать 8.5 млн призм для описания пространственной модели магнитной структуры ГИЦ. В то время как нестандартная версия программы повысила их число до 20 млн. что существенно увеличило точность вычислений и их предсказательную силу.
Расчетное ЗБ-поле позволяет при дальнейшем анализе установить тс или иные свойства изучаемой магнитной структуры, а также наметить пути, улучшающие их или устраняющие недостатки. Разработанная и оптимизированная трехмерная модель широко использовалась для всестороннего анализа пространственного поля и его интегральных характеристик. В частности, сравнивались с экспериментом среднее поле В(г). флаттер F. показатель роста магнитного поля к, абсолютные Ам (г) и относительные A'N(r) амплитуды вариации магнитного поля и другие параметры рассчитанного поля.
Радиальное распределение среднего по азимуту магнитного поля В (г) в медианной плоскости определяется, как известно, выражением
■2тт
B(r) = ± J Bz(r,<p,U)dip. (3)
о
Оно является первым и достаточно хорошим приближением к необходимому изохронному ПОЛЮ.
R, см R, см
Рис. 3. Сравнение вычисленного (сплошная линия) и экспериментального (пунктирная) среднего поля В (г).
Рис. 4■ Относительные амплитуды А]у(г) вариации магнитного поля ГИЦ.
Пунктирная линия - эксперимент, сплошная - расчет.
На рис. 3 видно, что всюду, начиная с радиуса I? ~ 20 см. в расчете достигнуто хорошее согласие с опытом по среднему полю. Отклонение от эксперимента не превышает 10-20Гс. В центре же магнита погрешность несколько выше и достигает примерно 50 Гс. Скорее всего, это различие связано с не вполне точным заданием геометрии исследуемой структуры в данной области. На вертикальную устойчивость ускоряемых в циклотроне частиц существенное влияние оказывает глубина вариации поля, так называемый флаттер і7, который вычисляется следующим образом:
В2 {г)
2п
где В2 (г) = [ В:(г.р.0)с1.(р.
-7Г >о
На рис. 4 видно, что в широком диапазоне радиусов до /? ~ 65 см совпадение с опытом можно считать полным.
Далее погрешность нарастает от близких к нулю значений до максимального, которое находится приблизительно в районе выводного радиуса. Погрешность при /? ~ 90 см составляет АР ~ 0.001, что соответствует примерно 4% от его расчетного значения.
Такое расхождение представляется вполне удовлетворительным для количественного описания экспериментальных данных и предсказания определенных свойств описываемой магнитной структуры. Кроме того, на вертикальную устойчивость оказывает не менее важное влияние и показатель роста магнитного ноля к(г) = (?’/Л) / (йВ/йг). Абсолютные Ар}(г) и относительные А н(г) = -4лг( г)/В(г) амплитуды вариации расчетного и экспериментального магнитных нолей извлекаются, как обычно, из стандартного разложения ноля В .(г, ср. 0) в ряд Фурье при каждом фиксированном значении г:
Вг (г, ер, 0) = В (г) + ^ [а„(г) соя пер + Ьп(г) зш тр\
71=1
В (г) ( 1 + ^ -4?г((г) соя [тир + Ф'(г)]
\ и= I
Здесь -4„(г) = \/а.2 (г) + Ь2 (г) - амплитуды вариации магнитного ноля и введено обозначение tg Ф(г) = Ъп(г)/ап(г).
Результаты, полученные для показателя роста магнитного ноля в этом варианте оптимизации, представляются также вполне хорошими. Степень согласия расчета с экспериментом продемонстрирована на рис. 5.
На рис. 6 видно, что расчетные и экспериментальные абсолютные амплитуды вариации магнитного ноля в широком диапазоне радиусов практически совпадают. Приблизительно вблизи радиуса II ~ 70 см начинается небольшое отличие в амплитудах.
Разница между расчетными и опытными значениями амплитуд монотонно нарастает примерно до 70 Гс в районе II ~ 90 см. Такое максимальное отклонение соответствует погрешности около 2%. Можно с уверенностью утверждать, что полученная степень согласия является вполне достаточной для того, чтобы иметь предсказательную силу но конструированию в последующих расчетах поправок к среднему нолю для перехода к изохронному. Из рис. 7 видно, что степень согласия расчета с экспериментом имеет тот же уровень, что и для абсолютных амплитуд вариации магнитного ноля ГИЦ. В районе II ~ 90 см погрешность снова не превышает 2%.
Я, см
Рис. 5. Сравнение расчетного (сплошная линия) и экспериментального (пунктирная) показателя роста магнитного поля к.
R, см
R, см
Рис. 6. Абсолютные амплитуды *4 = А к (Л) вариации магнитного
поля ГИЦ, N = 4.
Пунктирная линия - эксперимент, сплошная - расчет.
Рис. 7. Относительные амплитуды *4(1?)4(г) вариации магнитного поля ГИЦ.
Пунктирная линия - эксперимент, сплошная - расчет.
4. Заключение. Настоящая работа посвящена проблеме полномасштабного моделирования трехмерного поля магнитной структуры ГИЦ. Сложность поставленной задачи в первую очередь определяется существенной нелинейностью проблемы, а во вторую - большой спиральностью секторов ГИЦ, она была успешно преодолена при помощи программы MERMAID, позволяющей рассчитывать трехмерные статические магнитные поля сложной конфигурации. При этом были получены следующие основные результаты. Процедура построения оптимальной модели исследуемой магнитной структуры позволила достигнуть хорошего согласия с опытом по среднему полю. Отклонение от эксперимента не превышает 10-20 Гс. В районе выводного радиуса при I? = 90 см девиация составляет всего 32 Гс. Напомним, что эти величины следует сравнивать с уровнем магнитного поля в центре магнита 13 520Гс и с 14 670Гс на выводном радиусе. Также хорошо воспроизводится в расчетах глубина вариации магнитного поля ГИЦ, так называемого флаттера F. Можно считать, что в широком диапазоне радиусов совпадение с опытом является практически полным. Максимальная погрешность находится приблизительно в районе выводного радиуса при I? ~ 90 см и составляет
0.001, что соответствует примерно 4% от его расчетного значения. Такое расхождение представляется вполне удовлетворительным для количественного описания экспериментальных данных и предсказания определенных свойств описываемой магнитной структуры. Результаты, полученные для показателя роста магнитного поля к, демонстрируют почти полное согласие расчета с экспериментом и представляются также вполне хорошими. Логарифмическая производная по среднему полю за выводным радиусом и до него описывается достаточно правильно. Расчетные и экспериментальные абсолютные и относительные амплитуды вариации магнитного поля в широком диапазоне радиусов практически совпадают. Небольшое отличие в расчетных и опытных амплитудах начинается приблизительно вблизи радиуса с Л ~ 70 см и понемногу нарастает до приблизительно 70 Гс в районе I? ~ 90 см. Такое отклонение соответствует максимальной погрешности около 2% от расчетного значения. Полученная степень согласия является вполне достаточной для того, чтобы иметь предсказательную силу по конструированию поправок к среднему полю для перехода к изохронному полю. Сформулированный подход позволил с хорошей точностью описать реальный эксперимент, соответствующий определенной геометрии магнитной структуры ГИЦ. Тем самым открывается путь для проведения целого спектра численных экспериментов по созданию
необходимого изохронного поля ГИЦ. При этом не потребуются затраты материальных ресурсов на изготовление реальных шиммов. Основная работа будет выполняться на компьютере, опираясь на процедуру построения оптимальной модели исследуемой магнитной структуры.
Summary
Andrianov S. N., Artamonov V. A., Dubrovin A. N., Eliseev V. A. Computer 3D-modeling of magnetic field of isochronic cyclotron for GIZ.
The paper is devoted to the problem of magnetic structure modeling for GIZ. The complexity this problem is defined by its essential nonlinearity and large helicity of GIZ sectors. For this purpose the authors used the special programming package MERMAID, assigned for static magnetic fields. The well coincidence with experiments was achieved.
Литература
1. Абросимов H.K., Артамонов С. А., Елисеев В. А., Рябов Г. А. Разработка магнита и магнитной структуры циклотрона для ускорения Н--ионов до энергии 80 МэВ. Ч. 2. Исследования на моделях: Препринт Петерб. ин-та ядерной физики, № 2285. Гатчина, 1998. С. 3-32.
2. Dubrovin N. A. User’s guide MERMAID: Magnet design in two and three dimensions. SIM Limited, Novosibirsk department. Russia, 1994. P. 3-60.
3. Хейгеман Л., Янг Д. Прикладные итерационные методы / Пер. с англ. А. Ю. Еремина, И. Е. Капорина; Под ред. Ю. А. Кузнецова. М.: Мир, 1986. 448 с.
4. Suttman G., Steffen В. High-order compact solvers for the three-dimensional Poisson equation // J. of Comput. and Appl. Math. 2006. Vol. 187, Iss. 2. P. 142-170.
5. Самарский А. А., Николаев E.G. Методы решения сеточных уравнений. М.: Наука, 1978. 591 с.
Статья рекомендована к печати проф. Д. А. Овсянниковым.
Статья принята к печати 21 февраля 2008 г.
