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

2

МОДЕЛИРОВАНИЕ

УДК 517.946:518.12:538.3:538.5

ПРИМЕНЕНИЕ ТЕХНОЛОГИИ ВЫДЕЛЕНИЯ ПОЛЯ ПРИ КОНЕЧНОЭЛЕМЕНТНОМ МОДЕЛИРОВАНИИ КВАДРУПОЛЬНОЙ ЛИНЗЫ М.М. Корсун, А.Н. Игнатьев, М.Э. Рояк

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

Введение

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

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

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

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

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

rot Н 0 = J0,

div (д0Н 0 ) = 0, (1)

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

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

rot H = J,

1 i (2) div (дН ) = 0

будем искать в виде

Н = Н 0 + Н +, (3)

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

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

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

rot (H+ ) = J - J0,

(4)

div (|H + +(-|i0)) 0 ) = 0.

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

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

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

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

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

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

Н += НУ=- вгаё у. (5)

Q p

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

7 70

его как сумму двух полей - создаваемого разностью токов J - J в однородном пространстве и градиентом функции p (неполного потенциала). Таким образом, в области

Q p

H + = H++ H+= H+ - grad p, (6)

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

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

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

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

{ div (|H ++(|i-|0 )H 0 )uJQ = 0. (7)

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

| 61У(|ДЯ++(|Д-|Д°)Н° )иёО = 0.

| 61У(|Я ++(|-|°))0)иёО +

О, Ор

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

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

| 61У(|Я, +(|-|°))0 )иёО +

| 61У(Н+ +(1 )Н0 )иёО = 0.

(8)

О, Ор у ' ; 1

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

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

Н°иё£- | Н + §гаёиёО- | (1 )Н^гаёиёО

Ор

Ор

о \

! +

I Н+иёЯ + | (1 )Н°иё£- | |Н+ §гаёиёО- | )Н^гаёиёО

О,

О,

о\ л-о

I дН+иёЯр I (I-1°))0иёЯ = 0.

(9)

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

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

-||°Н°иё£- | Н+ §гаёиёО- | НБгаёиёО-

Ор

О,

))^гаёиёО+ | (|Н+-Н

О

рр(|-|,)Н,-(1 -д°р)Н° )иё£ = 0, (10)

¿7° ¿7° 0 0

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

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

, + Н+

№ н о

нр+Н++Н+

№ Н о

р

(11)

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

£

£

и:

+ Иу+)-Ир -И+

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

Vйш- и +- (1- ^р))р

и+ - мШИШ +

^ ¥

Р ¿>Р

р"р

и:

С учетом равенства (12), справедливого на поверхности 5и, уравнение (1Р) может быть записано в более простом виде:

^Р —► р -* р —►+ р —» +

V И - ] Ир §гаёиёО- ] дИу ёгаёиёО-

Ор

-{(ц-ДР)РБгаёиёО + | И+иё5 = Р. О 5„

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

Ор оу

= {|дРиРиё5 + )^гаё иёО- { И+иё5 . (13)

5 О Би

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

также необходимую для этого непрерывность нормальной составляющей вектора дИ на поверхности 5и = Оу I Ор. Чтобы на границе 5и между областями Оу и Ор выполнялись требуемые для решения исходного уравнения условия сопряжения, зададим соотношение между значениями полного и неполного потенциалов в следующем виде: Ш = р + и , (14)

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

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

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

ш= X ^ф/, р = X Л. (15)

/е/(ош) /е/(о р )

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

чая ее границы, а / (ор) - аналогичное множество для области Ор . Обозначим множество индексов узлов, принадлежащих поверхности 5и, как / (и / (5и ) = / (Ош) 1 / (Ор). Тогда из равенства (14) следует, что

т. е.

5

дШ = др + и, / е / (),

(16)

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

Л Г

§гаё и ёО + | V

О, ,

I

Ор

X др ф/

/е/(О р )

X дш вгаа ф/

/е/(О,)

§гаё иё О =

(17)

: | ИРиё5 + | ( - |дР))Р Бгаё иёО - | И+иё5.

5 О Би

Неизвестными в данном уравнении являются дш, /е / (Ош) и др, /е / (Ор ), т.е. узлам с номерами / е / (5и ) соответствуют по две неизвестные. Чтобы ввести единый век-

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

V, / е /(ош), р, /е/(ор)\/(и).

д/,

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

и/, /е / ((

Р, / е (/(Ош)и/(Ор))\/(

д/ =

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

ГN Л Г н

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

О \'=1 У ОД/=1

§гаё иё О =

||РИРиё5 + |(ц-|Р))^гаё иёО- | И+иё5,

5 О 5и

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

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

А = 11 §гаё ф/ §гаё фу ёО, /, у = 1, N, О

Г N Л

¥ = I X ди §гаё ф/

О р ^=1

О

§гаёф у ёО + I (|Д-|Р ))Р§гаёф ёО-

- | И+Яфуё5 +1 |РИРфуё5 , у = 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. Трехмерная конструкция

B r

rot— = J . д

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

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

rot — rot A - J. (19)

д

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

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

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

(

-div

1 д

grad Az

Л

= Jz.

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

(20)

(21)

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

B =dAz

Bx --

dy

By

y dx

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

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

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

1 2nr i \ i 2nr I \

An = — [ Ib sin In —I rd ф, Bn = — [ Ib cos In — I rd ф, nr 0 v r J nr 0 v r J

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

в r

Ib рассчитывается по формуле Ib = 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 получим двумерное решение, точность которого будем оценивать путем сравнения с решением двумерной задачи на подробной сетке (используемым как нормальное поле).

B2 B6 ^0

Двумерное решение (нормальное поле) 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. Результаты вычисления гармоник на основе двумерных решений

B2 B6 ^0

Решение исходной трехмерной задачи на грубой сетке 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).

Корсун Мария Михайловна Игнатьев Александр Николаевич Рояк Михаил Эммануилович

— Новосибирский государственный технический университет, аспирант, maria.korsun@gmail.com

— Новосибирский государственный технический университет, аспирант, ignat@hotmail.ru

Новосибирский государственный технический университет, кандидат технических наук, доцент, royak@fmp.ami.nstu.ru

УДК 519.6

ТРЕХМЕРНОЕ МОДЕЛИРОВАНИЕ САМОГРАВИТИРУЮЩЕГО ГАЗА

И.М. Куликов, В.А. Вшивков

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

Ключевые слова: газовая динамика, математическое моделирование, уравнение Пуассона, тесты ударной трубы, равновесные конфигурации

Введение

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

В настоящее время из всего широкого диапазона численных методов используются следующие: лагранжев метод сглаженных частиц SPH (Smoothed Particle Hydrodynamics), в основе которого лежит интерполяция расчетных ячеек в области сглаживания [3], и эйлеровы методы на адаптивных сетках AMR (Adaptive Mesh Refinement), базирующиеся в основном на кусочно-параболическом методе PPM (piece-parabolic method) [4], который является конечно-разностным методом высокого порядка точности типа метода Годунова. Метод сглаженных частиц SPH, разработанный в 1977 г. [5, 6], имеет большие воз-