Применение технологии выделения поля при конечноэлементном моделировании квадрупольной линзы Текст научной статьи по специальности «Математика»

ПРИМЕНЕНИЕ ТЕХНОЛОГИИ ВЫДЕЛЕНИЯ ПОЛЯ ПРИ КОНЕЧНОЭЛЕМЕНТНОМ МОДЕЛИРОВАНИИ КВАДРУПОЛЬНОЙ ЛИНЗЫ

М.М. Корсун, А.Н. Игнатьев (Новосибирский государственный технический университет)

Научный руководитель - к. т.н., доцент М.Э. Рояк (Новосибирский государственный технический университет)

Приводится вариационная и конечноэлементная постановки трехмерных задач магнитостатики с выделением нормального поля и использованием двух скалярных потенциалов. Применение рассматриваемого подхода показано на примере решения задачи численного моделирования квадрупольной линзы.

Введение

Довольно часто при моделировании трехмерных физических процессов, в том числе и магнитостатических, искомое поле имеет достаточно хорошее приближение, получаемое как решение другой, возможно, двумерной задачи, которое можно получить с более высокой точностью, чем требуется от решения исходной задачи. В том случае, если разница решений рассматриваемых задач составляет не более 10-15%, можно построить более эффективные как в плане вычислительных затрат, так и в плане точности расчетные схемы, учитывающие это обстоятельство. В этих схемах, основанных на выделении основной части поля, ставится задача на нахождение разницы полей, являющихся решением двух рассматриваемых задач, причем требования к точности решения такой задачи существенно ослабляется в силу того, что разница является достаточно малой по сравнению с искомым решением [1, 2].

Рассмотрим вычислительную схему, являющуюся модификацией схемы решения трехмерных задач магнитостатики с использованием двух потенциалов, с выделением главной части поля.

Математическая модель для метода скалярных потенциалов

Будем считать, что рассматриваемое магнитное поле имеет довольно хорошее приближение в виде решения более простой задачи, определяющей напряженность

магнитного поля Н0 , которая удовлетворяет системе уравнений Максвелла

rot Н 0 = J0, div (д0H 0 ) = 0,

(1)

где J0 - плотность стороннего тока, д0 - относительная магнитная проницаемость

среды. Заметим, что J0 и д0 могут, как совпадать с соответствующими величинами полной задачи, так и отличаться от них. Решение исходной задачи rot H = J,

' div(дН) = 0 будем искать в виде

(2)

Н = Н0 + Н +,

(3)

считая вектор-функцию Н известной и удовлетворяющей системе (1), при этом будем

называть поле Н0 нормальным, а поле Н + - аномальным (или добавочным).

Учитывая соотношения (1), (3), систему (2) можно переписать в следующем виде:

rot (Й + ) = J - J0

div (|Н ++(ц-|д° )) 0 ) = 0. ()

Важно отметить, что в случае равенства токов полной и нормальной задач первое уравнение системы (4) будет иметь нулевую правую часть, и в этом случае решение задачи нахождения аномального поля Н + является достаточно простым за счет возможности представления его в виде градиента некоторой функции. В дальнейшем будем

рассматривать общий случай, учитывающий, что поле Н + не является безвихревым.

Аналогично тому, как это сделано в [3], выделим в расчетной области О полной

задачи две подобласти: область Ор , содержащую токи J - J0, с относительной магнитной проницаемостью | = 1, и область Оу, содержащую ферромагнитные материалы. Обратим особое внимание на тот факт, что токи аномальной задачи представляют собой разность токов исходной и нормальной задач, при этом по-прежнему должно выполняться условие, что никакой контур, лежащий в области Оу, не должен охватывать

ненулевой ток.

В области Оу добавочное поле Н + является безвихревым, а значит, его можно определить как градиент некоторой функции у (полного потенциала):

Н + = Я+=- grad у. (5)

В области Ор добавочное поле Н + имеет ненулевой ротор, потому определим

его как сумму двух полей: поля, создаваемого разностью токов J - J0 в однородном

пространстве, и градиентом функции р (неполного потенциала). Таким образом, в обО р

Н + = Н+ + Н + = Н+ - grad р, (6)

ласти Q p

где rot (й+) = J - J0, div (й+) = 0.

Вариационная постановка и конечноэлементная дискретизация

При определении добавочного поля H + , согласно формулам (5), (6), первое уравнение системы (4) автоматически выполнилось, поэтому потенциалы у и p будем искать как решение второго дифференциального уравнения этой системы в области Q, при этом выберем область Q настолько большой, чтобы на ее границах можно было

положить H = H ++ H0 = 0. Для решения этого уравнения перейдем к эквивалентной вариационной формулировке, умножив это уравнение на пробную функцию и и проинтегрировав по Q:

J div+ + (|-|0)0 )udQ = 0. (7)

Представив область Q в виде объединения областей полного и неполного потенциалов Qy и Q p , перейдем от уравнения (7) к уравнению

"у -"p :

J div(|H ++й-|0))0)uJQ+ J div(|H ++(|-|0)H0)uJQ = 0.

Qy Qp

Далее, подставляя в полученное равенство соотношения (5), (6) и учитывая, что в

области О р д = 1, а поле Н+ является соленоидальным, получим

| (дЙ,

Ои

и

+ +(д-д0 ))0 )иёО + | ё1у (Й+ +(1 -д0) Я0 )иёО = 0.

Ор

(8)

Применяя к каждому слагаемому в левой части равенства (8) формулу Грина интегрирования по частям и учитывая, что на внешней границе Я области О

Н0 + Н + = 0, перейдем к следующему соотношению:

| д0Н - | Н+ §гаё иёО - | (1 - д0) Н 0 §гаё иёО

Ор

Ор

> +

| Н+иёЯ + | (1 -д0)Й0иёЯ- | дН+ БгаёиёО- | (д-д0)Н^гаёиёО

0\ £Т0

Я,

Я,

I дЙ+ийЯ + I (д-д0)Й0иёЯ = 0.

(9)

+ и

В данном равенстве - граничная поверхность между областями О+ и Ор с нормалью п+ , внешней для области О+, а поверхность - та же поверхность, но с нормалью Пр, внешней для области Ор . Таким образом, нормали п+ и пр являются

противоположно направленными, а значит, определив нормаль единым образом, например, п = п+, и обозначив соответствующую поверхность , поверхностные интегралы в соотношении (9) можно свести в один. Используя в объемном интеграле единое обозначение д для магнитной проницаемости во всей области О, получим .0 п 0

^0 -* 0 —* г ~~* + р -»• +

д Н иёЯ - | Нр §гаёиёО- | дН+ §гаёиёО

Я

Ор

О+

-1 (д- д0)Н^гаёиёО + | (дЙ+ -Н

О

+ + (д-д+)Н+ -(1 -дР)Нр = 0, (10)

где Й+ и Нр - нормальное поле, а д+ и д°р - относительные магнитные проницаемости в прилегающих к поверхности областях полного и неполного потенциалов, соответственно.

Рассмотрим более подробно выражение, стоящее внутри полученного поверхностного интеграла. Заметим, что на поверхности Би нормальные составляющие вектора

магнитной индукции как нормального, так и суммарного поля должны быть непрерывны. Это условие можно записать в виде следующей системы:

д(( И + Н, +

„0 н 0

и

Й°Р + Й+ + ЙС

1,0 Й0 \хрг1 р

(11)

Из пе

рвого уравнения системы (11) получаем

ЯС

д(( ++Н и)- Н р- Н+

и

Таким образом, - Й р

дЙИ -Й + +(д- дИИ)) -(1 -д"р))0

Я

Н + № Н0 +„° Н0

Нс -„ун у+„ Vм I

н+

(12)

'с 1 г р V

С учетом равенства (12), справедливого на поверхности Яи, уравнение (10) может быть

записано в более простом виде: .0 гт о

^0 ~* 0 -* р ~* + р I

„ Н иёЯ - ] Нр БгаёиёО- ] „Ну §гаёиёО

я Ор Оу

-|(„-„0)Н^гаёиёО + | Н+иёЯ = 0.

О ^

Перенося в правую часть все слагаемые, не содержащие неизвестного поля Н + в явном виде, и подставляя выражения (5), (6) для добавочного поля через градиенты скалярных потенциалов, получим окончательное уравнение: | §гаё р ёгаё и ёО + | „ §гаё у §гаё иёО =

= {„0Н0иёЯ +1(„-„0 ))^гаё иёО- { Н+иёЯ. (13)

Я О Би

Это уравнение обеспечивает выполнение равенства

Ау („Н ) = 0 в области О, а

также необходимую для этого непрерывность нормальной составляющей вектора „Н на поверхности Би = О, I Ор .

Для того чтобы на границе Яи между областями О, и Ор выполнялись требуемые для решения исходного уравнения условия сопряжения, зададим соотношение между значениями полного и неполного потенциалов в следующем виде:

у = р + и, (14)

где функция и - скачок (или разрыв) потенциалов аномальной задачи. Значения этой функции на поверхности Яи определяются из условия непрерывности тангенциальной

составляющей напряженности суммарного магнитного поля Н на этой поверхности. Далее для определенности представим область О в виде тетраэдров с линейными базисными функциями, а в качестве пробной функции и в уравнении (13) будем поочередно выбирать глобальные базисные функции, получаемые сшивкой локальных.

Представим функции у и р в виде линейных комбинаций базисных функций фг-:

у= X дуф/, р = X Чл Ф/ . (15)

/е/(оу) /е/(Ор)

Здесь / (оу) - множество индексов узлов сетки, принадлежащих области Оу, включая ее границы, а / (ор) - аналогичное множество для области Ор. Обозначим множество индексов узлов, принадлежащих поверхности как / (и), т.е.

/ (и ) = / (оу)п / (о р)

. Тогда из равенства (14) следует, что

ду = д[ + и/, / е / (Г), (16)

где П} - значение скачка потенциалов и в / -ом узле заданной сетки. Подставим выражения (15) в уравнение (13). В результате получим следующее уравнение:

Ор

X Чр Бгаё ф/

/е/(О р )

Л Г

§гаё и ёО + | д

О+

X 4+ Бгаё ф/

/е/(О+)

§гаё иё О =

(17)

= |д0Я0иёЯ + |(д-д0))^гаёиёО- | ЙС+иёЯ. Я О

Неизвестными в данном уравнении являются , / е /(ои) и qf , / е /(ор ), т.е. узлам с номерами / е / (Яи ) соответствуют по две неизвестные. Чтобы ввести единый

вектор неизвестных q, размерность которого совпадает с числом узлов конечноэле-ментной сетки (а, следовательно, и с числом базисных функций), используем соотношение (14). В этом случае компоненты вектора q можно определить следующим образом:

И, / е/(О+

q ,

qР, / е /(0р )\/(и).

qi =

Вводя дополнительный вектор q той же размерности, что и вектор q, компоненты которого определяются соотношениями

щ, / е / ((

0, / е (/ (О+)и / (Ор)) \ / (),

Л

§гаё иё О =

qU =

перепишем уравнение (17) в виде

Г N Л Г N

| д Xqi Бгаёф/ §гаёиёО- | X§гаёф/

О V/=1 У ОД / =1

= |д0Я0иёЯ + | (д - д0)^гаё иёО- | Й+иёЯ, Я О Яи

где N - полное число узлов в конечноэлементной сетке.

Выбирая поочередно в качестве пробной функции и базисные функции, получим систему уравнений, которая в матричном виде может быть записана как Aq = Ё , где

А = | д §гаё ф/ §гаё ф] ёО, /, у = 1, N О

Г N Л

=1 X qui §гаё ф/

О р ч/=1

У

.0 ^0,

О

§гаё ф у ё О + I (д-д0 )Й 0§гаё ф ¡ё О-

- | Й+Пфу ёЯ + | д0Й0фу ёЯ , у = 1, N.

Полученная система уравнений, очевидно, является в общем случае нелинейной, так как коэффициент относительной магнитной проницаемости д в большинстве реальных задач является функцией от напряженности магнитного поля, аномальная часть которого в области О+ определяется через полный потенциал + .

Результаты численного моделирования

Рассмотрим задачу моделирования постоянного магнитного поля в квадруполь-ной линзе, предложенную ИЯФ СО РАН им. Г.И. Будкера. Прежде всего, дадим описание расчетной области. Конструкция магнита симметрична относительно плоскости г = 0, что дает возможность решать задачу в одной, например, верхней, его половине, при задании соответствующих краевых условий на плоскости симметрии. Также плоскостями симметрии являются плоскости х = 0 и у = 0, а значит, мы можем перейти к решению задачи в области х > 0, у > 0, г > 0, что позволяет задавать в этой области более подробную сетку, чем это было бы возможно при решении задачи во всей конструкции. Требуемая симметрия обеспечивается краевыми условиями первого рода на границах х = 0, у = 0 и г = 0.

На рис. 1 представлена часть расчетной области, обладающая ферромагнитными свойствами. Высота рассматриваемого фрагмента - 68.8 см, внешний радиус - 42 см, толщина боковых стенок - 16.5 см. Полюс магнита в плоскости z = const описывается гиперболой, радиус гиперболы 23 см, высота полюса - 56 см.

Таким образом, в полной конструкции было четыре полюса, на каждом из которых была расположена токовая обмотка. Плотность тока во всех четырех обмотках одинакова и составляет 1780 кА/см2, сечение обмоток - 0.5 х 11.26 см, при этом направление токов в обмотках меняется при переходе к следующему полюсу.

Базовое сечение расчетной области, с учетом удаленных границ большого объема для задания первых краевых условий, имеет вид, представленный на рис. 2, а; на рис. 2, б, представлена дискретизация базового сечения конечными элементами.

X

Рис. 1. Трехмерная конструкция

а) б)

Рис. 2. Базовое сечение расчетной области

Квадрупольная линза имеет хорошее двумерное приближение, поэтому для ее моделирования можно очень эффективно применить технологию выделения нормального поля, поскольку мы можем очень точно решить двумерную задачу, а затем перейти к решению трехмерной задачи на аномалию. Для решения двумерной задачи воспользуемся постановкой с вектор-потенциалом. Стационарное магнитное поле описывается уравнением

rot — = J , (18)

где J - плотность тока, B - вектор магнитной индукции, | - коэффициент магнитной

проницаемости среды. Вводя векторный потенциал A такой, что B = rot A, преобразуем уравнение (18) к виду

rot — rot A = J. (19) I

В рас сматриваемой двумерной задаче вектор плотности тока J и векторный потенциал A имеют только одну ненулевую компоненту, т.е. J = (0,0, Jz) и

A = (0,0, Az ) . Это позволяет переписать уравнение (19) в скалярном виде:

(1 ^

-div

I

-grad Az

= Jz

(20)

На внешней границе Г расчетной области положим Аг равным нулю:

А |Г= о.

(21)

Используя найденное из уравнения (20) с краевыми условиями (21), распределение Аг , компоненты вектора магнитной индукции вычисляем по формулам

Вх = А , Бу =. ду у дх

Так как решение рассматриваемой двумерной задачи используется как нормальное поле для решения другой задачи, его необходимо получить с очень высокой точностью, поэтому конечноэлементная сетка для решения этой задачи была взята очень подробной. Число узлов при конечноэлементной дискретизации составило 688168. Ис-

пользование столь подробной сетки для решения задачи в трехмерной постановке очевидно нереально. На рис. 3 представлено распределение вектор-потенциала А.

Рис. 3. Решение двумерной задачи

Искомыми характеристиками, которые применяются на практике для оценки качества квадрупольной линзы, являются гармоники магнитного поля B. Гармоники представляют собой коэффициенты ряда Фурье:

An = — í Iв sin In —I rd ф, Bn = — í ¡b cos In —I rd ф, nr 0 ^ r J nr 0 ^ r J

где n - номер гармоники, r - радиус, на котором проводятся измерения, а величина

в r

¡B рассчитывается по формуле ¡в = J B—(r, —)dz . В силу симметрии конструкции все

гармоники An = 0, а ненулевыми являются только гармоники Bn с номерами n = 4i - 2,

i = 1, да .

Оценим эффективность применения технологии выделения поля, вычислив гармоники с номерами 2, 6 и 10 при r = 19 см, а = 0 см, в = 56 см, в сравнении с результатами, полученными непосредственно при решении трехмерной исходной задачи (2). Для этого решим задачи (2) и (4) на вложенных конечноэлементных сетках, при этом за грубую примем сетку с числом узлов 28368 (дискретизация базового сечения такой сетки изображена на рис. 2, б). Следовательно, удвоенная и учетверенная трехмерные сетки будут содержать 217469 и 1708689 узлов соответственно. При этом оценку точности конечноэлементной аппроксимации построенной трехмерной сетки для решения задачи (2) напрямую дополнительно проверим следующим образом: зададим ферромагнитную часть расчетной области и обмотки как бесконечно длинные. В этом случае мы фактически будем решать двумерную задачу, т.е. в любом сечении z = const полу-

а

чим двумерное решение, точность которого будем оценивать путем сравнения с решением двумерной задачи на подробной сетке (используемым как нормальное поле).

В2 Вб В10

Двумерное решение (нормальное поле) 0.50528427 0.00004262 0.00000254

Решение исходной трехмерной задачи как двумерной на грубой сетке 0.50477600 0.00029760 0.00001469

Решение исходной трехмерной задачи как двумерной на удвоенной сетке 0.50510539 0.00028654 0.00005079

Решение исходной трехмерной задачи как двумерной на учетверенной сетке 0.50521215 0.00010346 0.00001654

Таблица 1. Результаты вычисления гармоник на основе двумерных решений

В2 В6 В10

Решение исходной трехмерной задачи на грубой сетке 0.49887337 -0.00011270 0.00052052

Решение исходной трехмерной задачи на удвоенной сетке 0.49909997 0.00000941 0.00057038

Решение исходной трехмерной задачи на учетверенной сетке 0.49910367 -0.00011407 0.00053633

Решение задачи с выделением аномальной части на грубой сетке 0.50067227 -0.00060202 0.00039387

Решение задачи с выделением аномальной части на удвоенной сетке 0.50094538 -0.00057088 0.00039105

Решение задачи с выделением аномальной части на учетверенной сетке 0.501045418 -0.00056978 0.00041062

Таблица 2. Результаты вычисления гармоник на основе трехмерных решений

Из полученных результатов (см. табл. 1) решения двумерной задачи (нормальное поле) и трехмерной задачи как двумерной видим, что трехмерность задачи вносит существенную погрешность в решение. С дроблением конечноэлементной сетки, конечно, трехмерное решение постепенно сходится к точному решению, однако вычислительные затраты, например, на учетверенной сетке являются уже критическими, и дальнейшие дробления практически невозможны.

Рассмотрим результаты решения трехмерных задач, приведенные в табл. 2. Здесь следует отметить, что аномалия поля достаточно низка, т.е. на 2-3 порядка меньше нормального поля. Вследствие этого можем утверждать, что технология выделения поля уже на самой грубой сетке дает результаты гораздо точнее, чем решение трехмерной задачи (2) напрямую на самой подробной сетке.

Заключение

Разработанная вычислительная схема реализована в конечноэлементном программном комплексе ТБЬМЛ и довольно эффективно применяется для решения сложных практических задач.

Построенная вычислительная схема решения трехмерных задач магнитостатики с использованием скалярных потенциалов, основанная на выделении главной части напряженности магнитного поля, позволяет получить более высокую точность решения, чем исходная схема, при проведении расчетов на одной и той же конечноэлементной сетке. Чтобы получить такую же точность базовым методом, необходима гораздо более подробная сетка, использование которой приводит к значительному увеличению потребностей в ресурсах ЭВМ и временных затрат на решение поставленной задачи.

Высокая точность и эффективность такой вычислительной схемы, основанной на том, что главная часть поля определяется из решения более простых (двумерных) задач с высокой точностью, продемонстрирована на примере решения практической задачи моделирования магнитостатического поля в квадрупольной линзе.

Литература

1. Игнатьев А.Н., Рояк М.Э. Выделение основной части поля при решении трехмерных нелинейных задач магнитостатики // Материалы 8 международной конференции Актуальные проблемы электронного приборостроения. - Новосибирск, 2006. - № 6. - С. 37-44.

2. Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач: учеб. пособие. - Новосибирск: Изд-во НГТУ, 2007. - 869 с.

3. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Решение трехмерных нелинейных маг-нитостатических задач с использованием двух потенциалов. - Новосибирск, 1996. -28 с. (Препринт / РАН. Сиб. отд-ние. ВЦ; № 1070).