Содержание к диссертации
Введение
Глава 1. Решение трехмерных нелинейных задач магнитостатики с использованием двух скалярных потенциалов 14
1.1. Основная схема решения трехмерных задач магнитостатики с использованием двух скалярных потеницалов 14
1.1.1. Математическая модель 14
1.1.2. Вариационная постановка 16
1.1.3. Вычисление скачка потенциалов 18
1.1.4. Вычисление значений полного и неполного потенциалов 25
1.2. Модификация метода скалярных потенциалов, основанная на выделении главной части поля 27
1.2.1. Основной принцип математических моделей с выделением главной части поля 27
1.2.2. Математическая модель для метода скалярных потенциалов 28
1.2.3. Вариационная постановка и конечноэлементная дискретизация... 30
1.3. Решение задач магнитостатики в анизотропных средах 36
1.3.1. Пересчет коэффициента магнитной проницаемости для шихтованных материалов 36
1.3.2. Модификация метода скалярных потенциалов для моделирования полей в анизотропных средах 41
Глава 2. Решение нелинейных систем уравнений, получаемых при моделировании трехмерных магнитостатических полей 43
2.1. Основные итерационные методы решения нелинейных систем уравнений 43
2.1.1. Метод простой итерации 43
2.1.2. Метод Ньютона з
2.1.3. Линеаризация конечноэлементных нелинейных систем уравнений для метода Ньютона 45
2.2. Учет нелинейности в коэффициенте магнитной проницаемости при решении задачи с использованием двух скалярных потенциалов 47
2.3. Учет нелинейности в коэффициенте магнитной проницаемости при решении задачи с выделением главной части поля 52
2.4. Учет нелинейности в коэффициенте магнитной проницаемости при решении задачи в шихтованных средах
2.4.1. Линеаризация конечноэлементной системы уравнений для решения задачи без выделения главной части поля 56
2.4.2. Линеаризация конечноэлементной системы для решения задачи с выделением главной части поля 60
2.5. Вычисление значений магнитной проницаемости и ее производных по заданной таблице значений 67
Глава 3. Примеры решения задач 69
3.1. Моделирование поля в вигглере 69
3.1.1. Описание расчетной области 69
3.1.2. Математическая модель для решения двумерной задачи магнитостатики в декартовых координатах 72
3.1.3. Результаты решения двумерной задачи 74
3.1.4. Результаты решения трехмерной задачи 76
3.2. Моделирование поля в магнитной системе циклотрона 80
3.2.1. Описание конструкции 80
3.2.2. Математическая модель для решения двумерной осесимметричной задачи магнитостатики 84
3.2.3. Результаты решения двумерной задачи 91
3.2.4. Результаты решения трехмерной задачи 92
3.3. Моделирование поля в магните с шихтованным железом 96
3.3.1. Описание конструкции 96
3.3.2. Результаты решения задачи 100
3.4. Моделирование поля в С-образном диполе 102
3.4.1. Описание конструкции 102
3.4.2. Результаты решения задачи 104
Заключение 107
Список использованных источников
- Вычисление скачка потенциалов
- Линеаризация конечноэлементных нелинейных систем уравнений для метода Ньютона
- Математическая модель для решения двумерной задачи магнитостатики в декартовых координатах
- Моделирование поля в магните с шихтованным железом
Введение к работе
Актуальность темы. Основным инструментом современной экспериментальной физики атомного ядра и элементарных частиц более полувека служат ускорители заряженных частиц. Существующие ускорительные установки постоянно модернизируются, сооружаются новые ускорители. В циклическом ускорителе магнитная система является основным узлом, обеспечивающим устойчивость движения пучка.
В связи с тем, что ускорители являются достаточно дорогими установками, одним из основных методов, используемых при проектировании и создании магнитных систем, является математическое моделирование. Расчет магнитных систем ускорителей - достаточно сложная задача математической физики, требующая математических исследований при разработке численных методов, эффективных программных реализаций, а также больших ресурсов вычислительных машин. Математическое моделирование дает возможность резко уменьшить время анализа поля в магните выбранной конфигурации, повысить точность, сократить стоимость и такого анализа, и самого магнита, так как непосредственное измерение самого магнитного поля является трудоемкой и дорогостоящей проблемой. Наряду с этим, математическое моделирование позволяет исследовать и те части конструкции магнита, измерения в которых крайне затруднительны или даже невозможны (например, распределение индукции в магнитопроводе традиционных магнитов), но распределение поля в этих частях сказывается существенным образом на характеристиках и работе магнита.
В конечном счете, только математическое моделирование магнитной системы позволяет сделать выбор оптимальной конструкции магнита в каждом конкретном случае. Основной (и очень важной) характеристикой решения, получаемого в результате моделирования, является его точность. Как правило, получение решения с высокой точностью требует очень больших вычислительных ресурсов. Так, например, при конечноэлементном моделировании самым простым и наиболее часто используемым способом повышения точности является дробление сетки, на которой решается задача. Однако, при решении задач, где точность расчета трехмерного магнитного должна достигать сотых долей процента (в частности, в задачах проектирования установок с электронным охлаждением), такой способ неэффективен. Это связано с тем, что использование чрезмерно подробной сетки приводит к существенному увеличению необходимой памяти и, что еще более важно, времени, требуемых для решения задачи и последующей выдачи.
Таким образом, несмотря на то, что к настоящему времени разработан достаточно большой объем программно-математического обеспечения, реализованного в таких широко распространенных программных комплексах, как ANSYS, FLUX3D, OPERA3D, NASTRAN и т.д., остается множество важных практических задач, для которых получение результатов требуемой точности весьма затруднительно.
Важно отметить, что подавляющее большинство современных магнитных систем выполнено из шихтованных материалов, обладающих анизотропными свойствами. Производство таких систем значительно дешевле, чем производство систем из цельных ферромагнитных деталей. Еще одним достоинством систем из шихтованных материалов является отсутствие вихревых токов при включении магнитов. К тому же, благодаря возможности соединения блоков из шихтованных материалов в единые фрагменты непосредственно при сборке конструкции магнита, для этих систем существенно упрощается их транспортировка. Таким образом, моделирование магнитных систем из шихтованных материалов является актуальной задачей, которая, тем не менее, на сегодняшний день не достаточно глубоко исследована, а потому большинство распространенных программных комплексов не позволяют получать результаты необходимого качества при расчетах магнитных полей в таких конструкциях.
В данной работе предложены вычислительные схемы моделирования магнитных полей в сложных магнитных системах, позволяющие получать более точное решение без значительного увеличения затрачиваемых ресурсов компьютера. Эти схемы позволяют проводить высокоточное моделирование магнитных полей в системах со сложной геометрией, в том числе и в системах, выполненных из анизотропных материалов, на персональных компьютерах.
При отсутствии поверхностных токов и токов, протекающих по ферромагнетику, задача нахождения распределения магнитного поля, созданного стационарными токами в проводниках, сводится к нахождению векторных функций напряженности магнитного поля и магнитной индукции из системы уравнений Максвелла для стационарного магнитного поля. При нелинейной зависимости магнитной проницаемости от напряженности магнитного поля и сложной геометрии областей решение этой системы возможно только численными методами, базирующимися на различных постановках, включающих интегральные, дифференциальные и комбинированные. Для решения трёхмерных задач магнитостатики наиболее распространенными являются постановки относительно скалярного и векторного потенциалов.
Векторная постановка приводит, как правило, к гораздо более трудоёмким вычислительным схемам, так как в этом случае при решении задачи либо нужно искать три неизвестные функции, удовлетворяющие системе из трёх скалярных уравнений, либо решать задачу как векторную с использованием специальной технологии представления векторного потенциала. И в том, и в другом случае вычислительные затраты гораздо выше, чем при решении одного скалярного уравнения для достижения той же точности вычисления характеристик магнитного поля. Дополнительной трудностью является неоднозначность определения векторного потенциала. Однако векторный потенциал широко используется при расчете магнитных полей в двумерной постановке.
Гораздо более эффективным для решения трёхмерных задач магнитостатики является использование скалярного магнитного потенциала. В
этой постановке вектор напряженности магнитного поля определяется как разность напряженности магнитного поля, создаваемого токами в вакууме, (определить ее можно, например, по закону Био-Савара) и градиента скалярной функции, называемой неполным потенциалом. Такая постановка позволяет получать приемлемые результаты при вычислении магнитного поля вне ферромагнитных материалов, однако внутри этих материалов значение поля вычисляется с очень низкой точностью. Это связано с тем, что в материалах с большой относительной магнитной проницаемостью напряженность магнитного поля может быть много меньше соответствующих значений поля, создаваемого теми же токами в вакууме. Таким образом, напряженность магнитного поля в ферромагнетиках получается как разность двух близких по значению величин, и даже если погрешность расчета этих величин составляет десятые доли процента, то в результате их вычитания она может стать неприемлемой.
В результате, использование только одного (неполного) скалярного потенциала при решении нелинейных задач магнитостатики неэффективно, так как требует очень высокой точности вычисления самого потенциала и его градиента для того, чтобы обеспечить необходимую точность вычисления коэффициента магнитной проницаемости. Эти трудности могут быть преодолены введением так называемого полного потенциала с тем, чтобы в ферромагнетиках напряженность магнитного поля представлялась непосредственно как градиент этой функции.
В данной работе подробно рассматривается метод решения трехмерных нелинейных задач магнитостатики с использованием двух скалярных потенциалов в дифференциальной постановке, являющийся одним из наиболее эффективных методов численного решения трехмерных нелинейных задач магнитостатики. В то же время, относительно полного скалярного потенциала возможно получить интегральную постановку, а в тех случаях, когда коэффициент магнитной проницаемости не зависит от поля, применима также и граничная постановка.
Отметим, что этот метод возможно сочетать с другими методами моделирования магнитных полей (например, с методом граничных элементов), используя конечные элементы со скалярным потенциалом лишь в части области, содержащей нелинейные ферромагнитные материалы, что позволяет избежать построения сетки в вакууме и тем самым значительно сократить размерность получаемой системы уравнений. При этом все предложенные в работе схемы для метода скалярных потенциалов остаются корректными.
Цель исследования состоит в разработке эффективных вычислительных схем для решения трехмерных нелинейных задач магнитостатики, учитывающих шихтованность материалов, и их программной реализации.
Научная новизна
1. Разработаны и реализованы вычислительные схемы решения трехмерных нелинейных задач магнитостатики с учетом анизотропных свойств материалов, в том числе и схема с выделением главной части поля.
Предложена методика учета шихтованных материалов при решении трехмерных нелинейных задач магнитостатики.
Проведены исследования эффективности предложенных вычислительных схем при моделировании магнитных систем ускорителей заряженных частиц.
Проведены исследования влияния шихтовки на характеристики С-образного дипольного магнита при постоянном токе.
Личный вклад автора работы заключается в разработке вычислительных схем моделирования трехмерных нелинейных магнитостатических полей с выделением нормального поля и с учетом анизотропии в коэффициенте магнитной проницаемости среды, программной реализации всех описываемых методов в программном комплексе MASTAC и проведении исследований эффективности предложенных методов.
На защиту выносятся:
Вычислительная схема решения трехмерных нелинейных задач магнитостатики с выделением главной части поля, основанная на использовании двух скалярных потенциалов, позволяющая учитывать анизотропные свойства среды;
Методика учета шихтованных материалов при решении трехмерных нелинейных задач магнитостатики, позволяющая с достаточной точностью моделировать шихтовку материала в произвольном направлении;
Алгоритмы решения систем нелинейных уравнений, получаемых в результате применения разработанных вычислительных схем, на основе метода Ньютона.
Результаты исследования влияния шихтовки на характеристики С-образного дипольного магнита.
Практическая ценность работы и реализация результатов Разработанные автором вычислительные схемы реализованы в программном комплексе MASTAC и широко применяются для решения многих сложных практических задач. В диссертационной работе приводятся несколько примеров решения таких задач:
высокоточное моделирование магнитной системы вигглера с использованием поэтапного выделения поля;
высокоточное моделирование магнитной системы циклотрона с использованием поэтапного выделения поля;
моделирование поля в магните из шихтованного железа;
моделирование поля с С-образном диполе.
Результаты диссертационной работы использовались при выполнении научно-исследовательских хоздоговорных работ НГТУ:
«Конечноэлементные исследования трёхмерных магнитных полей дипольных магнитов» (2006 г., НИУ ИЯФ СО РАН).
«Конечноэлементные исследования магнитных полей дипольных магнитов НЕВТ и МЕВТ с учетом шихтованности и сложной геометрии» (2007 г., НИУ ИЯФ СО РАН);
Разработанное программное обеспечение активно используется при выполнении расчетов в ИЯФ СО РАН им.Г.И.Будкера.
Достоверность результатов подтверждается как решением модельных задач, так и сравнением результатов численного моделирования с экспериментальными данными.
Апробация работы
Основные результаты работы были представлены и докладывались на: восьмой и десятой международной научно-технической конференции «Актуальные проблемы электронного приборостроения» АПЭП-2006, 2010 (Новосибирск, 2006г. и 2010 г.); V Всероссийской межвузовской конференции молодых ученых (Санкт-Петербург, 2008 г.); а также на семинарах ИЯФ СО РАН (г. Новосибирск), ОИЯИ (г. Дубна).
Результаты исследований изложены автором в 9 печатных работах, из которых 2 опубликованы в изданиях, включенных в перечень ВАК ведущих рецензируемых научных журналов и изданий, выпускаемых в Российской Федерации, 3 - в сборниках научных трудов, 4 - в сборниках трудов конференций.
Структура и объем работы
Вычисление скачка потенциалов
Для того чтобы уравнение (1.13) имело решение, соответствующее физическому смыслу задачи, то есть, чтобы результирующее поле удовлетворяло на границе Su между областями ОУ и Qp условиям сопряжения, на этой границе необходимо задать соотношение между значениями полного и неполного потенциалов. Это можно сделать, например, в следующем виде: ц/ = р + и. (1.14) Функция и в соотношении (1.14) называется скачком (или разрывом) потенциалов. Определим значения этой функции на поверхности S" из условия непрерывности тангенциальной составляющей напряженности магнитного поля Н на этой поверхности, а именно, из равенства
Это уравнение определяет скачок потенциалов и на границе Su таким образом, что если полный и неполный потенциал связаны соотношением (1.14), то напряженность магнитного поля Н имеет непрерывную тангенциальную составляющую на этой поверхности.
Как следует из уравнения (1.18), для определения скачка потенциалов необходимо знать величину напряженности магнитного поля токовых обмоток Нс, которая при решении реальных задач, как правило, не имеет аналитического представления, но может быть вычислена в произвольной точке пространства, исходя из закона Био-Савара [56]: где V — радиус-вектор точки расчета, V — объем токовой обмотки, г — радиус-вектор элемента объема dV, а у (г) - плотность тока этого элемента объема. Таким образом, процедура вычисления значений скачка потенциалов и должна учитывать тот факт, что величина Щ может быть вычислена лишь в отдельных точках поверхности Su.
Учитывая, что вся область Q. разбита на трехмерные конечные элементы, в дальнейшем будем рассматривать - поверхность Su как объединение некоторого числа непересекающихся многоугольников S%, являющихся гранями исходных элементов. В этом случае, с учетом того, что функции у/ и р ищутся в пространстве кусочно-линейных функций, уравнение (1.18) можно переписать в виде -—-НІ, (1.19) діі где { Ц } - ребра элементов S% поверхности Su, а Hf. \ — проекции вектора Нс на эти ребра. От решения системы уравнений вида (1.19) удобно перейти к решению вариационной задачи /W=Zj\w+Hh - min (L2) минимизирующей сумму квадратов отклонений производной функции и по ребрам элементов S% от проекций вектора Нс на эти ребра. Полученную задачу также удобно решать с помощью МКЭ, а именно, представить искомую функцию и в виде линейной комбинации базисных функций, определенных на ребрах: Подставляя данное равенство в функционал (1.20), получим l(u) = l{q) =
Чтобы найти точку минимума квадратичной формы (1.22), продифференцируем ее по компонентам вектора q и приравняем полученные производные к нулю. В результате получим СЛАУ следующего вида: решая которую с помощью итерационных методов получаем веса q, и, соответственно, аппроксимацию функции и на поверхности Su. Отметим, что в подавляющем большинстве задач достаточная точность вычисления функции и обеспечивается выбором в качестве базиса линейных полиномов. Тем не менее, для решения этой задачи можно использовать и квадратичный базис, причем в случае использования иерархического базиса квадратичного пространства, увеличение вычислительных затрат при этом не существенно, так как матрица, соответствующая СЛАУ (1.23), остается такой же, как и при использовании линейных базисных функций. Кроме того, портрет этой матрицы совпадает с реберным портретом треугольной конечноэлементной сетки, на которой решается поставленная задача, что позволяет применить стандартный алгоритм его построения.
Основной особенностью реализации расчета скачка потенциала является
тот факт, что поверхность Su, являющаяся областью расчета, хранится как двумерная конечноэлементная сетка, тогда как задача решается на конечных элементах, являющихся ребрами этих двумерных элементов.
Линеаризация конечноэлементных нелинейных систем уравнений для метода Ньютона
Программная реализация этого метода чуть более сложна, чем реализация метода простой итерации, но зато скорость сходимости метода Ньютона как правило значительно выше. Важно отметить, что сходимость этого метода также не гарантирована. Иногда хорошего результата можно достичь путем последовательного применения этих методов. Например, можно решать нелинейную систему методом Ньютона до тех пор, пока невязка решения существенно меняется, затем сделать несколько итераций методом простой итерации, после чего вернуться обратно к использованию метода Ньютона.
Важно отметить, что для улучшения сходимости можно ввести дополнительный параметр релаксации со є(0,і]. В этом случае после решения СЛАУ решение системы ищут в виде qu}=COq+(l-CD)q\ (2.4) подбирая со таким образом, чтобы невязка решения q03 для системы была минимальной, и при переходе на следующую итерацию полагают q := qa.
Рассмотрим решение конечноэлементной системы нелинейных уравнений с помощью метода Ньютона. На каждом шаге итерационного процесса метода Ньютона линеаризуем все уравнения полученной системы. При этом будем обозначать q решение, полученное на предыдущей итерации (в случае первой итерации, q — начальное приближение). Учитывая тот факт, что матрица А и вектор правых частей F конечноэлементной системы являются сумой локальных матриц Ат и локальных векторов правых частей Fm, которые вычисляются на элементах Qw, а затем дополняются до размерности глобальной матрицы, перенесем процесс линеаризации на этап генерации локальных матриц и векторов.
В дальнейшем изложении будем пользоваться локальной нумерацией неизвестных, полагая, для определенности, что мы работает с узловыми базисными функциями такими, что на каждом конечном элементе определено N базисных функций. Для перехода к глобальной нумерации достаточно каждому локальному индексу конечного элемента С1т поставить в соответствие глобальный номер соответствующего узла этого элемента.
Линеаризуя у-ое слагаемое z -ro уравнения через разложение его в ряд Здесь индексы і и j, как и индекс к, принимают значения от 1 до 4, что определяется числом неизвестных на одном конечном элементе. При этом частные производные, стоящие в выражении (2.5), имеют следующий вид:
Собирая линеаризованные локальные матрицы и векторы в глобальную матрицу, получим глобальную линеаризованную систему уравнений, которая используется в методе Ньютона.
Дальнейшие выкладки, связанные с линеаризацией системы уравнений,
уже непосредственно зависят от функций A"? (g) и F (g), поэтому они будут рассмотрены для конкретных задач в следующем разделе этой главы.
Учет нелинейности в коэффициенте магнитной проницаемости при решении задачи с использованием двух скалярных потенциалов Вернемся к решению нелинейной системы (1.45). Правая часть этой системы не зависит от вектора неизвестных q, поэтому чтобы решить эту систему с помощью метода Ньютона, нам остается линеаризовать только левую ее часть A{jq)q.
Как описано в 2.1.3, перенесем линеаризацию системы на уровень локальных матриц. Для того чтобы вычислить необходимые для линеаризации частные производные, стоящие в формулах (2.6), вернемся к выражению (1.36). На каждом конечном элементе Qm, лежащем в области полного потенциала Q , напряженность магнитного поля Н определяется соотношением N Н = -gradi// = -qkgrad ., k=l где рь, k = l,N - локальные базисные функции, определенные на конечном элементе Qm. Если эти функции являются линейными, то их градиенты являются константными векторами, а значит и напряженность магнитного поля Н на каждом конечном элементе является постоянной величиной, квадрат модуля которой может быть вычислен следующим образом:
В случае выбора базиса более высокого порядка, несмотря на то, что градиенты базисных функций уже не будут постоянными на элементе, мы также можем воспользоваться этим подходом, взяв в качестве некоторой аппроксимации реального значения Н на элементе его среднеинтегральное значение. Для этого проинтегрируем равенство (2.8) по конечному элементу С1т и разделим результирующее равенство на объем этого конечного элемента
Математическая модель для решения двумерной задачи магнитостатики в декартовых координатах
Вигглер - (от англ. wiggle — вихлять, изгибаться, ёрзать) устройство для генерации синхротронного излучения в электронном накопителе-синхротроне. Вигглер представляет собой магнит, создающий сильное поперечное (как правило, вертикальное) знакопеременное магнитное поле. Его можно представить себе как последовательность коротких дипольных магнитов, полярность каждого следующего из которых противоположна предыдущему.
Рассмотрим задачу моделирования постоянного магнитного поля в вигглере. Прежде всего, дадим описание расчетной области. Отметим, что конструкция вигглера симметрична, поэтому в качестве расчетной области возьмем одну ее четверть, при этом на границах зададим соответствующие краевые условия, которые опишем ниже.
На рис. 1 представлена часть расчетной области, обладающая ферромагнитными свойствами. Длина рассматриваемого фрагмента — 66.8 см, высота задней стенки — 16 см. Ширина дальнего полюса — 6 см, ширина остальных полюсов — 7.4 см, ширина боковой стенки - 2 см, расстояние между полюсами — 4.8 см. Высота основания конструкции — 5.1 см, высота полюсов — 9.3 см. Отметим, что в расчетную область попала только половина шестого полюса, так как одна из плоскостей симметрии прошла через его середину. Вторая плоскость симметрии проходит параллельно основанию рассматриваемой конструкции через верхнюю плоскость задней стенки. График зависимости /л(В) для описанной ферромагнитной области приведен на рис. 2 .
График зависимости ju(B) Таким образом, в полной конструкции было 22 полюса, на каждом из которых была расположена токовая обмотка. Поэтому, даже при решении задачи в четверти полной конструкции, необходимо задать все 22 обмотки. Плотность тока во всех обмотках одинакова и составляет 40.2299 А/см2, сечение обмоток на крайних четырех полюсах (расположенным возле боковых стенок) — 2х 4.3см, сечение остальных 18 обмоток — 2x8.7 см, при этом направление токов в обмотках меняется при переходе к следующему полюсу как в верхней, так и в нижней половинах, в результате чего распределение потенциала имеет вид, представленный на рис. 1.
Для решения этой задачи необходимо также дополнить расчетную область, поместив рассматриваемую конструкцию достаточно большим баком 2x2x1 м. Ввиду симметричности полной конструкции на верхней границе расчетной области следует задать первое однородное, а на плоскости, проходящей через середину шестого полюса - второе однородное краевые условия. На остальных границах задается первое однородное краевое условие ввиду удаленности этих границ.
Искомой величиной в рассматриваемой задаче является индукция магнитного поля в той части горизонтальной плоскости симметрии, которая расположена непосредственно между полюсами.
В этой области задача имеет хорошее двумерное приближение, поэтому для ее решения можно очень эффективно применить технологию выделения нормального поля, а именно, мы можем очень точно решить двумерную задачу в плоскости, проходящей через середины полюсов (параллельно задней стенке), а затем перейти к решению трехмерной задачи на аномалию. Следует отметить, однако, что нам придется решить две аномальные задачи, так как решить трехмерную задачу в исходной конструкции от ее двумерного приближения в данной постановке невозможно: токи двумерной задачи будут проходить через заднюю стенку, обладающую ферромагнитными свойствами. В связи с этим необходимо решить последовательность из двух аномальных задач: 1) от двумерного приближения без задней стенки; 2) от предыдущей задачи (трехмерной без стенки, посчитанной как аномальная) уже с задней стенкой. Отметим, что решения этих двух задач в исследуемой области практически не отличаются, что будет показано ниже. Математическая модель для решения двумерной задачи магнитостатики в декартовых координатах
В рассматриваемой двумерной задаче вектор плотности тока J и векторный потенциал А имеют только одну ненулевую компоненту, т.е. J = (0,0,Jz) и A = (0,0,AZ). Это позволяет переписать уравнение (3.2) в скалярном виде: На внешней границе Г расчетной области положим Az равным нулю, то есть = JZ. (3-3) Используя найденное из уравнения (3.3) с краевыми условиями (3.4) распределение Az, компоненты вектора магнитной индукции вычисляются по следующим формулам: Для решения поставленной задачи используем метод конечных элементов, предварительно перейдя к эквивалентной вариационной формулировке. Умножим обе части уравнения (3.3) на пробную функцию у/, которая принимает нулевое значение на границах расчетной области Q, и проинтегрируем полученное равенство по Q: - \ div — grad Az у/dQ. = \ Jzy/dQ. (3.6)
Применяя к первому слагаемому в левой части уравнения (3.6) формулу Грина (интегрирования по частям) и учитывая, что у/ на границах Q обращается в ноль, получаем вариационное уравнение Г—grad Azgrad y/dQ. = Г Jzy/dQ.. (3.7) Для дискретизации расчетной области будем использовать треугольные конечные элементы. Введем базисные функции ц/:, номера которых совпадают с номерами узлов сетки. Каждая базисная функция равна единице в узле с соответствующим номером, нулю в остальных узлах и является линейной функцией координат х и у на треугольниках.
Моделирование поля в магните с шихтованным железом
Искомой величиной при решении данной задачи являлась магнитная индукция поля в центральной плоскости между полюсами магнита, причем необходимая точность была порядка 5 гаусс. Получение такой точности при решении задачи прямым методом невозможно, так как даже для решения задачи в одной четверти конструкции (учитывая ее симметрию) необходима очень подробная сетка, работа с которой потребовала бы очень больших вычислительных ресурсов, недоступных пользователям персональных компьютеров. В связи с этим возникла необходимость в решении задачи с выделением главной части поля.
Учтем сразу все плоскости симметрии, чтобы уменьшить расчетную область и получить возможность использования более подробной сетки. Прежде всего, конструкция симметрична относительно плоскости z = 0, расположенной посередине между полюсами магнита. Более того, конструкция симметрична относительно плоскости у = 0, делящей ярмо магнита на две равные Ш-образные части. Таким образом, от решения задачи в полной конструкции можно сразу перейти к решению задачи в ее четверти, поэтому в дальнейшем будем рассматривать часть конструкции, расположенную в области у О, z О.
Также следует заметить, что если прямоугольную часть ярма магнита дополнить таким образом, чтобы оно стало осесимметричным, то задачу можно решать в 90 -ом секторе, что позволяет задать в расчетной области еще более подробную сетку. Полученное в результате решения такой задачи поле можно использовать как нормальное при решении задачи в четверти исходной конструкции.
Задачу в 90-ом секторе также можно решать как аномальную относительно двумерного поля, полученного при решении соответствующей осесимметричной задачи, не учитывающей спиральные шиммы. Для этого в двумерной задаче достаточно все пространство, расположенное между шиммами, положить также выполненным из ферромагнетика.
Стационарное магнитное поле описывается уравнением (3.1). Вводя векторный потенциал А такой, что В = rot А, преобразуем это уравнение к виду (3.2). В осесимметричных конструкциях плотность тока и векторный потенциал А имеют только одну ненулевую ср -компоненту, зависящую от г и z. Это позволяет переписать уравнение (3.2) в скалярном виде: -div dz Для решения поставленной задачи используем метод конечных элементов, предварительно перейдя к эквивалентной вариационной формулировке. Умножим обе части уравнения (3.12) на пробную функцию ц/, которая принимает нулевое значение на границах расчетной области Г2, и проинтегрируем полученное равенство по Q.:
Применяя к первому слагаемому в левой части уравнения (3.15) формулу Грина (интегрирования по частям) и учитывая, что у/ на границах Q обращается в ноль, получаем вариационное уравнение
Для дискретизации расчетной области будем использовать треугольные конечные элементы. Введем базисные функции у/.-, номера которых совпадают с номерами узлов сетки. Каждая базисная функция равна единице в узле с соответствующим номером, нулю в остальных узлах и является линейной функцией координат г и z на треугольниках. Будем считать, что на каждом конечном элементе JU = JU(B) является константой и определяется по среднему значению В на элементе. Рассмотрим третье слагаемое левой части вариационного уравнения (3.16). Оно представляет собой объемный интеграл, содержащий поверхностную 8-функцию (т.к. на границах некоторых конечных элементов м(В) терпит разрыв, то есть дг А ). производная от разрывной функции). Этот объемный интеграл может быть вычислен как интеграл по границам, где /л терпит разрыв. В результате уравнение (3.16) преобразуется к виду
Важно заметить, что полученные уравнения нелинейны относительно неизвестных qi, поскольку компоненты матрицы G неявно зависят от qt (через зависимость /л от величины В).
Для решения нелинейной конечноэлементной системы (3.20) будем использовать метод Ньютона, проводя линеаризацию отдельных слагаемых уравнений системы на уровне локальных матриц. На каждом конечном элементе коэффициент /л = /л(В) считается постоянной величиной и определяется средним значением модуля индукции на элементе. На практике удобнее использовать В — среднее значение квадрата модуля магнитной индукции. Получим выражения для вычисления значения Л 2 / Л 2 В на конечном элементе Q . Представим В в виде:
